Geomorphology-Guided Optimization of Hybrid Solar-Wind Systems: A Model with Comparative Evidence in Colombia

1. Introduction

The decarbonization of the electricity matrix requires systems that combine different renewable sources in proportions adapted to the characteristics of each territory. Compared to single-technology installations, hybrid solar-wind matrices offer more stable generation throughout the year and a better utilization of the available area. Nevertheless, the optimal size depends on a resource that could vary according to the geomorphology of the terrain [1,2].

The influence of the relief on the wind resource is widely documented: wind speed at hub height is modulated by terrain roughness, slope, and exposure, which has motivated the development of statistical and machine learning-based downscaling methods that correct satellite reanalysis products in complex terrain [3,4,5]. Recent studies couple computational fluid dynamics (CFD) simulation with machine learning models to predict wind production in wind farms located on complex orography [6], while others have quantified the bias of the fifth-generation European Centre for Medium-Range Weather Forecasts atmospheric reanalysis (ERA5) dataset products in regions with steep topography, finding systematic errors that are only corrected by incorporating terrain variables [7,8,9]. The solar resource exhibits an analogous dependence, although of opposite sign: horizon shading and orographic cloudiness reduce available irradiance with elevation, an effect that literature on photovoltaic potential estimation from digital elevation models has characterized in detail [10]. Taken together, this evidence suggests that a hybrid sizing model that ignores the geomorphological heterogeneity of the territory risks assigning an optimal configuration for the municipality, being not optimal for any of its individual zones.

To solve the sizing problem, literature has increasingly turned to mixed-integer linear programming (MILP), which guarantees a globally optimal solution within a discrete set of candidate configurations. It has recently been used to size hybrid photovoltaic systems with hydrogen storage [11], to optimize wind farm layout by coupling a wake model with integer programming [12], for the optimal design of offshore wind farms under resource uncertainty [13], for the stochastic planning of multiple microgrids [14,15], and for optimizing the collector system topology of large-scale photovoltaic plants [16]. In parallel, another body of work has combined geographic information systems (GIS) with machine learning for site selection, integrating multiple criteria and variable reduction through spatial clustering [17], or fuzzy multi-criteria analysis for wind farm siting [18,19,20]. These two research problems, namely, the mathematical optimization of sizing and the spatial selection of the site have largely developed separately, and few studies integrate both within a single methodological workflow that is applied and empirically validated for different territories.

A practical limitation of integrating detailed physical simulations within an optimization process is their high computational cost. Evaluating the configuration of each candidate through simulations at every point in the territory can become infeasible when the analysis is conducted at the municipal scale, owing to the large number of possible combinations and the required spatial resolution. Surrogate models address this limitation by learning the relationship between the input variables of the system and its performance, thereby replacing repeated simulation with near-instantaneous prediction. Machine learning-based surrogate models have been used to accelerate the optimization of energy systems at the national scale, with reported computation time reductions of up to 64% compared to the conventional approach [21]; to quantify uncertainty reduction in wind measurement campaigns using linear and machine learning comparison models [22]; and to estimate, through polynomial response surfaces, the uncertainty of CFD-based wind site assessment [23,24]. This body of evidence supports the use of a surrogate model as a bridge between high-fidelity physical simulation and optimization. However, to date, no references were identified that combine an explicit geomorphological stratification of the territory with a formal hybridization constraint subject to MILP optimization.

The present work addresses this gap through a four-stage model: geomorphological stratification of the territory by elevation ranges, solar and wind resource characterization specific to each stratum, physical simulation of generation coupled with a decision tree surrogate model that reproduces this simulation with high fidelity, and optimization through mixed-integer linear programming subject to a hybridization constraint that prevents a single technology from occupying more than 80% of the installed area. The model was applied and contrasted in two Colombian territories with opposite resource profiles: Villavicencio (Meta state), a municipality in the Llanos foothills, interconnected to the national electrical system, and Manaure (La Guajira state), a coastal municipality with high wind potential. The contrast between the two territories constitutes the main verification of the proposed approach. In one of the cases, the model determined that the wind resource was not sufficient to justify a hybrid configuration and selected an exclusively solar solution. In the other, it identified favorable conditions for a hybrid matrix, with a wind share of 73%. These results show that the same methodology, applied autonomously and without external intervention, was able to adapt the technological combination to the local characteristics of each territory.

The main contribution of this work is the development of an integrated and reproducible workflow that couples MILP optimization with a surrogate model, in which geomorphological stratification is hybridized using globally available satellite-derived data to support renewable energy source selection across a territory. A complementary contribution is the demonstration, through explicitly contrasting case studies, that the proposed framework maintains consistent performance under opposing resource conditions, thereby supporting its generalizability without territory-specific parameterization.

This paper is organized as follows Section 2 details the materials, data, and methodological approach; Section 3 presents the results obtained in two geographically different locations; Section 4 discusses these results in the context of the literature; and Section 5 presents the conclusions and limitations of the work.

2. Materials and Methods

This section presents the model, formulates and explains how it will be applied. The proposed model is developed through a sequence of four stages. First, a geomorphological stratification is performed based on a digital elevation model (DEM), classifying the territory into units with flat, undulating, or mountainous relief. Next, the energy resource characterization is carried out, using the National Solar Radiation Database (NSRDB) for solar and the System Advisor Model (SAM) for wind data, extrapolated to the operating height of the wind turbines. In the third stage, simulations and surrogate modeling are performed, using PVWatts, developed by the National Renewable Energy Laboratory (NREL), for photovoltaic generation and PyWake, a Python-based wind-farm simulation framework developed by DTU Wind Energy, for wind generation. An objective was to build simplified, high-precision models (R² > 0.99). Finally, optimization is applied through MILP, implemented with Pyomo and solved with the GNU Linear Programming Kit (GLPK) solver, to select the optimal technological combination. As a result, the model assigns to each zone the most suitable technology, required area, and energy production, configuring the optimal solution according to the local conditions of each territory.

The analysis was carried out in Python 3.11.11, integrated within the ArcGIS Pro 3.5 notebook environment, and is applied independently to two localities with significant differences in geomorphology and climate. The first case corresponds to the municipality of Villavicencio (Meta state in Colombia), for each of its four geomorphological strata, and, as a contrasting case, to the municipality of Manaure (La Guajira state in Colombia). Figure 1 summarizes the methodological workflow of the model, whose four stages will be described in Section 2.a through 2.d.

2.1. Study Area and Digital Elevation Model

Villavicencio is a municipality in the Llanos foothills, interconnected to the national electrical system, whose elevation ranges from the low savanna to the Eastern Cordillera; it is also humid and transitional. Manaure, in contrast, is a coastal municipality in La Guajira with predominantly flat topography, exposed to the Caribbean trade winds, and characterized by dry, sunny, and windy conditions. This contrast in relief and wind regime motivated its selection as a contrasting case within the same methodological framework.

The DEM of Villavicencio was constructed from official contour line cartography using the Topo to Raster interpolation tool based on the Australian National University Digital Elevation Model (ANUDEM) algorithm, which generates hydrologically coherent surfaces, with the resulting DEM projected in the EPSG:32618 coordinate reference system (WGS 84 / UTM zone 18N) with a spatial resolution of 12.5 m. The resulting model covered an elevation range between 208 m and 3,612 m. Slope, calculated in degrees, and the Terrain Ruggedness Index (TRI), which quantifies local elevation variability, were derived from the DEM [25]; both attributes describe the morphological complexity of each zone and constitute the basis of the stratification.

2.2. Geomorphological Stratification and Resource Characterization

The territory of Villavicencio was segmented into four geomorphological strata defined by elevation ranges, reflecting the transition from the plain to the mountain range: low savanna (0 m–500 m), low foothills (500 m–800 m), high foothills (800 m–1500 m), and mountain (above 1500 m) [26]. Each stratum was represented by an actual point from the National Solar Radiation Database (NSRDB) [27], located within its boundaries, which provides hourly time series of solar radiation, temperature, wind speed, and pressure over the course of one year. The Manaure locality was characterized in the same way, using its own representative NSRDB point, without subdivision into strata, given that its relief does not present the altitudinal variability found in Villavicencio. Representing each zone by a single point constitutes an arbitrary simplification that does not capture the spatial variability within each stratum; this limitation is revisited in Section 5.

The NSRDB time series were converted to the SAM format, the standard required by the simulation models used. Descriptive resource statistics were calculated from each file: mean global horizontal irradiance (GHI), mean temperature, mean wind speed at 10 m above ground level, and mean pressure. Since wind turbines operate at a hub height of 110 m rather than at the 10 m measurement height, wind speed was extrapolated using the power law according to Equation (1)

$$v\left(h\right)=vo{\left({\displaystyle \frac{h}{ho}}\right)}^{\alpha },$$

(1)

where $v\left(h\right)$ is the wind speed at height $h$, $vo$ is the wind speed is measured at the reference height $ho$ (10 m), and $\alpha$ is the shear exponent, set at 0.143, a characteristic value for open terrain. Machine learning-based extrapolation methods have shown improvements over the classical power law under variable atmospheric stability conditions [28], an option identified as a direction for future work.

2.3. Generation Simulation

Photovoltaic generation is simulated using the PVWatts model, through its PySAM interface [29]. PVWatts estimates the annual energy produced by a photovoltaic system based on solar resource and a set of technical system parameters. The optimal orientation of the panels, tilt and azimuth is not assumed: the model is run for each combination of tilt, within a range of ±15° around the latitude of the site, and azimuth, between 90° and 270°. The combination that maximizes annual production is selected. The simulation also incorporates the horizon shading effect, calculated from the surrounding terrain elevation profile derived from the DEM.

Wind generation is simulated using an engine that combines three components: the Jensen wake model, also known as NOJ (N. O. Jensen), implemented through the open-source library PyWake [30], which estimates the energy loss due to interaction between neighboring turbines within the same wind farm; the Hellmann shear profile, described in Equation (1), which extrapolates the measured wind speed to hub height; and an orographic correction derived from the DEM, which estimates the wind acceleration or deceleration effect caused by local topography. The coupling of a wind farm wake model, a vertical shear profile, and a topographic correction seek a higher-fidelity estimate of wind production than would be provided by an isolated power curve. Recent reviews of wake models document equivalent coupling approaches [31]. Before optimization, the viability of the wind resource is evaluated using hub-height wind speed, comparing it against two thresholds: 3.0 m/s for conventional horizontal-axis turbines and 1.5 m/s for vortex-induced vibration (VIV) technology. When hub-height wind speed does not reach the conventional threshold, the model restricts candidate configurations to zero wind capacity, rather than forcing an infeasible hybrid solution.

2.4. Surrogate Model: Decision Tree

The physical simulation described above is accurate but computationally expensive. Evaluating all candidate configurations in each zone through direct simulation is impractical. To address this, a surrogate model, a decision tree is introduced, which learns the relationship between conditions of the system and its energy efficiency, to predict that efficiency for new configurations in milliseconds, without repeating the simulation, which is in line with the growing use of surrogate models to accelerate the optimization of energy systems [32].

The tree was trained using the actual resource data of each zone, generating 500 synthetic sites through sampling of statistical distributions fitted to that resource. Irradiance was sampled from a gamma distribution, wind speed from a Weibull distribution, both standard choices for these variables and temperature and pressure from normal distributions. Each of these 500 sites was crossed with twelve solar and wind capacity configurations, generating 6000 cases. For each case, the full physical simulation was run, and the target variable was calculated, defined as the total energy produced divided by the total area occupied, i.e., the energy density, in kWh per square kilometer. The eight input variables of the tree were solar capacity, panel efficiency, wind capacity, hub height, solar radiation, wind speed, temperature, and pressure.

The tree structure was tuned through an exhaustive hyperparameter search using GridSearchCV from the scikit-learn library in Python. A three-fold cross-validation scheme was employed to evaluate combinations of maximum tree depth (4, 6, 8, and 10), minimum number of samples required to split an internal node (5, 10, and 20), minimum number of samples required at a leaf node (5 and 10), and the cost-complexity pruning parameter (0, 0.001, and 0.005), with the aim of preventing overfitting. The data were split into 80% for training and 20% for testing. Across all zones evaluated, the coefficient of determination (R²) on the test set exceeded 0.99, indicating that the surrogate reproduces the physical simulation with high fidelity, consistent with the accuracy reported for decision tree-based surrogate models applied to energy systems [21,33].

The choice of a decision tree over more complex surrogate modeling approaches, such as neural networks and ensemble methods, was based on three criteria. First, its interpretability: the tree structure and the variable importance analysis using the Gini index allow clear identification of which factors determine efficiency, something less evident in low-transparency models. Second, its parsimony: a single tree, regularized through pruning, reduces the risk of overfitting on a synthetic dataset without requiring the design and tuning of a complex architecture. Finally, its predictive sufficiency: by achieving a coefficient of determination above 0.99, the surrogate model reproduces the physical simulation with minimal error, so a more complex alternative would be unlikely to provide a significant improvement.

2.5. Selection of the Optimization Method

Error! Reference source not found.summarizes the mathematical properties of each family across six dimensions relevant to renewable energy project design: computational complexity, convergence behavior, global optimality guarantee, constraint handling capacity, treatment of non-linearity, and solution traceability. These dimensions correspond directly to the requirements that the proposed framework places on the optimization component: the method must handle mixed discrete-continuous variables, enforce hard spatial and hybridization constraints, and produce solutions whose optimality can be verified and documented for engineering and regulatory purposes.

Figure 2 provides a graphical synthesis of these properties as a radar chart comparing all six dimensions simultaneously across the four method families. Scores represent a qualitative assessment of the properties summarized in Table 1, on a 1 (low) to 5 (high) scale for visualization purposes.

Table 2 and Figure 2 converge on the same conclusion: the surrogate-assisted MILP occupies the most favorable position across all six dimensions, achieving global optimality and full constraint enforcement while reducing computational cost relative to exact MILP by eliminating repeated simulation calls through the surrogate. HOMER and RETScreen were excluded from the workflow because full enumeration over the joint solar-wind capacity space does not certify optimality beyond the pre-defined configuration grid and has no native mechanism for the spatial hybridization constraint introduced in the following subsection. Metaheuristics offer the best native handling of non-linearity, but cannot certify optimality or enforce hard constraints algebraically, a disqualifying limitation for applications where solution traceability is required for engineering documentation or regulatory submission [1,34]. Exact MILP without a surrogate remains a valid alternative in principle, but requires linearizing the coupled PVWatts–PyWake simulation beforehand, which the decision-tree surrogate achieves empirically with R² > 0.99, as detailed in the following subsections.

2.6. Formulation of the Optimization Model

The selection of the optimal configuration for each zone is formulated as a MILP problem, implemented using the Pyomo library [35], and solved with the GLPK 4.65 solver. The problem starts from a set of candidate configurations, each defined by a combination of solar and wind capacity in discrete increments (0, 1, 2, and 3 MW). Each configuration $i$ is assigned to a binary variable $x_{\mathrm{i}}$, which takes the value 1 if the configuration is selected and 0 otherwise. The objective function, Equation (2), maximizes the efficiency predicted by the surrogate model, $eff\mathrm{i}$, summed over the selected configurations:

$$\max {\sum }_i{eff}_i x_i ,$$

(2)

subject to the constraints of Equations (3) to (5). Equation (3) prevents the total installed area from exceeding the available area $A$:

$${\sum }_i\left(a_i^s+ a_i^w\right)x_i \le A,$$

(3)

Where $a_i^s$ and $a_i^w$ are the solar and wind areas required by configuration$i$, respectively, and $A$ is the available area.

Equation (4) ensures that exactly one configuration is selected:

$${\sum }_i{x}_i=1,$$

(4)

And Equation (5), which is activated only when the wind resource is viable, requires that at least one configuration belonging to the set$H$ of admissible hybrid configurations be selected:

$${\sum }_{i\in H}{x}_i \ge 1,$$

(5)

The set $H$ of admissible hybrid configurations is formally defined according to Equation (6):

$$H = \{i \in C: a_{\mathrm{i}}^{\mathrm{s}} \le 4·a_{\mathrm{i}}^{\mathrm{w}} \wedge a_{\mathrm{i}}^{\mathrm{w}} \le 4·a_{\mathrm{i}}^{\mathrm{s}}\},$$

(6)

Where $C$ is the set of all candidate configurations, a condition equivalent to no technology occupying more than 80% of the total installed area. The combination of Equations (4) and (5) guarantees that, when the wind resource is viable, the only active configuration belongs to $H$. Thus, Equation (4) establishes that the number of selected configurations is exactly one, and Equation (5) requires that at least one configuration from $H$ be active, such that this single configuration must belong to$H$. When the wind resource is not viable, the set of candidate configurations is restricted beforehand to those with zero wind capacity, so that the model selects the best solar-only option without applying Equation (5) and formally documents this limitation. This branching does not modify the structure of the model, but rather the set over which it operates, and is handled as preprocessing, prior to constructing the decision variables. The model was solved repeatedly over a sequence of available-area scenarios, ranging from 0.10 to 2.00 km², to observe how the optimal configuration evolves as more area becomes available.

The MILP approach for sizing hybrid systems, with analogous area and technology-mix constraints, has precedents that address related sizing problems with storage [11] and joint layout-and-technology optimization [12,16]; however, unlike these, the formulation used here couples the objective function directly to a validated surrogate model, rather than to a closed-form analytical model, which allows the nonlinearity of the physical simulation to be transferred into the MILP without sacrificing the global optimality guarantee characteristic of integer linear programming.

2.7. Algorithms

The three core computational stages of the workflow, coupled resource simulation, surrogate model training, and MILP optimization, are formalized in Algorithm 1, Algorithm 2 and Algorithm 3. Algorithm 1 governs the overall procedure for a given zone, from data loading through selection of the optimal configuration, and invokes Algorithm 2 and Algorithm 3 as subroutines in Steps 6 and 7, respectively. The complete code, executed on the territories of Villavicencio and Manaure, is documented in the Supplementary Materials.

Algorithm 1. Characterization and resource assessment per zone
Input	DEM of the zone NSRDB→SAM resource file for the representative point.
Output	Optimal renewable configuration, allocated areas, and estimated production.
Step: 1	Validate the computational environment and the integrity of the inputs (raster processing libraries, GLPK solver, DEM, and resource file).
Step: 2	Load the DEM; verify its reference system (EPSG:32618) and resolution (12.5 m); derive the slope and the Terrain Ruggedness Index (TRI).
Step: 3	Segment the territory into four strata according to elevation ranges. (Only in the stratified run).
Step: 4	Read the resource data from the representative point; calculate the resource statistics (GHI, wind speed at 10 m, temperature, pressure); extrapolate the wind speed to hub height using the power law (Equation (1)).
Step: 5	Simulate generation: solar with PVWatts, optimizing tilt and azimuth and incorporating horizon shading and wind with PyWake—Jensen wake model; evaluate wind viability against the thresholds of 3.0 m/s (conventional) and 1.5 m/s (VIV).
Step: 6	Train the surrogate model (Algorithm 2)
Step: 7	Solve the optimization for each available-area scenario (Algorithm 3); consolidate the optimal configuration, areas, and production.

Algorithm 2 details the procedure invoked in Step 6 of Algorithm 1, in which a synthetic dataset spanning the resource space of the zone is generated to train the decision-tree surrogate that subsequently replaces the coupled PVWatts–PyWake simulation during optimization.

Algorithm 2. Surrogate model training
Input	Resource statistics of the zone $\mathrm{Set} \mathrm{of} \mathrm{capacity} \mathrm{configurations} C$.
Output	Trained decision tree.
Step: 1	Generate 500 synthetic sites by sampling distributions fitted to the actual resource: irradiance ∼ Gamma, wind ∼ Weibull, temperature and pressure ∼ Normal (fixed seed); bound to physical ranges.
Step: 2	Cross$\mathrm{the} 500 \mathrm{sites} \mathrm{with} \mathrm{the} 12 \mathrm{configurations} \mathrm{in} C$, generating 6000 cases.
Step: 3	For each case: simulate solar (PVWatts) and wind (PyWake) energy; calculate total energy and total area; calculate the efficiency (Equation (7)).
	$efficiency = {\displaystyle \frac{total energy}{total area}}$,	(7)
Step: 4	Build the feature matrix X (8 variables) and the target y (efficiency).
Step: 5	Split the data into 80% training and 20% testing.
Step: 6	Tune the hyperparameters through exhaustive search (GridSearchCV) with 3-fold cross-validation, over maximum depth, minimum samples per split and per leaf, and pruning parameter; select the best tree.
Step: 7	Evaluate the coefficient of determination (R²) on the test set and calculate feature importance (Gini).

With the surrogate model trained, Algorithm 3 details the procedure invoked in Step 7 of Algorithm 1. For each available-area scenario, candidate solar–wind configurations are evaluated using the surrogate and the optimal configuration is selected subject to the hybridization constraint.

Algorithm 3. MILP optimization per area scenario
Input	Trained surrogate model Resource conditions $\mathrm{Available} \mathrm{area} A$ Wind viability indicator.
Output	Selected configuration, areas, and energy.
Step: 1	Construct the candidate configurations (solar capacity, wind capacity) from the discrete set {0, 1, 2, 3} MW.
Step: 2	If the wind resource is not viable, restrict the candidates to zero wind capacity.
Step: 3	$\mathrm{For} \mathrm{each} \mathrm{candidate} i$$: \mathbf{calculate} \mathrm{the} \mathrm{required} \mathrm{areas} (a_{\mathrm{i}}^{\mathrm{s}}$$\mathrm{and} a_{\mathrm{i}}^{\mathrm{w}}$$) \mathrm{and} \mathrm{predict} \mathrm{their} \mathrm{efficiency} eff$ᵢ using the surrogate model.
Step: 4	$\mathrm{If} \mathrm{the} \mathrm{resource} \mathrm{is} \mathrm{viable}, \mathbf{identify} \mathrm{the} \mathrm{set} H$$\mathrm{of} \mathrm{admissible} \mathrm{hybrid} \mathrm{configurations} (a_{\mathrm{s}} \le 4·a_{\mathrm{w}}$$\mathrm{and} a_{\mathrm{w}} \le 4·a_{\mathrm{s}},$$\mathrm{that} \mathrm{fit} \mathrm{within} A$).
Step: 5	Define$\mathrm{the} \mathrm{binary} \mathrm{variables} x_{\mathrm{i}}.$
Step: 6	$\mathrm{Objective} \mathrm{function}: \mathbf{maximize} \Sigma {eff}_{\mathrm{i}}·x_{\mathrm{i}}$(Equation (2)).
Step: 7	Constraint$: total area \le A$(Equation (3)); exactly one configuration (Equation (4)); if viable, at least one hybrid (Equation (5)).
Step: 8	Solve using GLPK through Pyomo and return the selected configuration, its areas, and its energy.

This complete procedure was executed independently for Villavicencio considered as a whole, for each of its four geomorphological strata, and for Manaure as a contrasting case. For each zone, the model delivered the highest-efficiency configuration, its estimated daily production, the allocated solar and wind areas, and, when applicable, documentation of wind resource infeasibility. The results of this procedure are presented in Section 3.

3. Results

3.1. Villavicencio: Municipality as a Whole

To perform generalization of the model, two sites in Colombia are analyzed, Figure 3, where Colombian Andes mountains are pictured. The first location is Villavicencio, for which resource characterization yielded a mean global horizontal irradiance of 4.65 kWh/m²/day and a wind speed of 0.87 m/s at 10 m. When extrapolated to hub height using Equation (1), reached 1.22 m/s. This speed did not reach either of the two operating thresholds evaluated, 3.0 m/s for conventional technology or 1.5 m/s for vortex-induced vibration technology. So, the wind resource was classified as non-viable, and the set of candidate configurations was restricted to zero wind capacity, formally documenting the limitation rather than forcing a hybrid solution.

The surrogate model reproduced the physical simulation with a coefficient of determination of 0.9995. The feature importance analysis, using the Gini index, indicated that wind capacity was the most determinant variable for efficiency (71%), followed by solar radiation (15%) and solar capacity (14%), refer to Figure 4; this hierarchy, even with the wind resource classified as non-viable, reflects that wind capacity continues to govern the variance of the synthetic training set, on which the non-viability condition was not imposed a priori.

The optimization was solved over the sequence of available-area scenarios, ranging from 0.10 to 2.00 km², and all of them selected exclusively solar configurations. In the first scenario (0.10 km²), 2 MW were proposed for installation, with a production of 6,883 kWh/day; from 0.31 km² onward, capacity saturated at 3 MW, with a constant production of 10,325 kWh/day, revealing that, beyond a certain surface area, the active constraint ceases to be the available area and becomes the maximum installable capacity considered instead (Figure 5). The optimal configuration, with a tilt of 8° and an azimuth of 240°, reached an annual production of 3,768,628 kWh, equivalent to an energy density close to 25.12 GWh per year per square kilometer.

3.2. Results by Geomorphological Stratum

To overcome the limitation of a single representative point, Villavicencio was stratified into four elevation-based strata, each characterized by its own NSRDB representative point (Figure 6). The low savanna stratum (E1, 0–500 m) dominates the municipal area (1009 km², 77%), while the low foothills (E2, 500–800 m), high foothills (E3, 800–1500 m), and mountain (E4, 1500–3612 m) strata together account for the remaining 23% of the territory, concentrated in its northwestern portion. Each representative point was selected from NSRDB grid coordinates falling within its corresponding stratum. As noted in the technical limitations, this approach characterizes each climate of the stratum using data from a single point, rather than an average across its full spatial extent.

The model was run independently for each of the four geomorphological strata of Villavicencio, using their respective representative points. Across all four strata, the surrogate model reproduced the physical simulation with a coefficient of determination above 0.996, and the feature importance analysis maintained the same pattern observed for the municipality as a whole: wind capacity was dominant (between 72% and 75%), followed by solar capacity and solar radiation. The optimal panel tilt decreased from 10° in the low savanna to 4° in the mountain stratum, resulting from the orientation optimization at each representative point (Figure 7).

The solar resource exhibited a decreasing gradient with altitude: irradiance decreased from 5.10 kWh/m²/day in the low savanna to 3.75 kWh/m²/day in the mountain stratum, a difference of 1.35 kWh/m²/day. The wind resource, in contrast, showed an increasing gradient: hub-height wind speed rose from 0.51 m/s in the low savanna to 1.58 m/s in the mountain stratum. This inverse behavior between the two resources, with the solar resource decreasing and the wind resource increasing with altitude, is consistent with the climatic transition from the plain to the Eastern Cordillera.

Despite this increase, the wind resource did not reach the conventional operating threshold of 3.0 m/s in any stratum, so all four strata were classified as non-viable for wind generation, and the model selected exclusively solar configurations in all of them. The optimal production accordingly followed the solar resource: the low savanna led with 11,129 kWh/day and the highest energy density (27.08 GWh/year/km²), while the mountain stratum recorded the lowest values (8,451 kWh/day and 20.56 GWh/year/km²). Stratification thus fulfilled its methodological purpose: it differentiated the resource across zones and revealed an altitudinal gradient that an aggregated analysis of the municipality would have concealed, even though the operational conclusion generalized wind non-viability, coincided with that obtained for the municipality.

3.3. Contrasting Case: Manaure

To examine the behavior of the model under an abundant wind resource, the same procedure was applied to Manaure, in the state of La Guajira. The resource of Manaure differed substantially from that of Villavicencio: mean irradiance was higher (5.97 kWh/m²/day), but the decisive difference lay in the wind. Hub-height wind speed reached 7.17 m/s, well above the conventional threshold of 3.0 m/s, so that only 0.6% of the evaluated sites fell below that threshold; the wind resource was accordingly classified as viable for conventional technology.

This regime shift was reflected in the structure of the surrogate model: whereas in Villavicencio efficiency depended mainly on wind capacity, in Manaure the feature importance shifted toward the wind resource: wind speed (67.5%) and wind capacity (19.4%) dominated, relegating solar radiation and solar capacity to a secondary role. Gini importance is presented in Figure 8, which constitutes the statistical signature of a wind-dominant regime.

The optimization showed behavior dependent on the available area. In the smallest scenario (0.10 km²), the model selected an exclusively solar configuration of 2 MW (9227 kWh/day): with such a reduced surface area, a single turbine did not satisfy the 80% hybridization constraint, so the pure solar solution was the best feasible option. From 0.31 km² onward, in contrast, the model selected a genuine hybrid configuration of 1 MW solar and 1 MW wind, with a production of 17,098 kWh/day and a split of 73% wind and 27% solar (Figure 9). The optimal configuration occupied 0.170 km² (0.050 km² solar and 0.120 km² wind) and reached an annual production of 6,240,594 kWh, equivalent to an energy density close to 36.71 GWh per year per square kilometer, higher than that of the best stratum of Villavicencio. The 80% hybridization constraint operated as expected: above the minimum required area, it forced an effective mix of both technologies rather than concentrating all capacity in one.

3.4. Contrast and Synthesis

The results for Villavicencio and Manaure, presented in Table 1, obtained using the same model and the same procedure, allow the behavior of the methodology to be contrasted under opposing resource conditions. Manaure outperformed Villavicencio in solar resource (by 1.32 kWh/m²/day), but the decisive difference was in wind: 5.95 m/s more in hub-height speed. That gap, rather than the solar one, defined the optimal configuration in each territory: where the wind was insufficient, the model selected a pure solar solution and documented the infeasibility; where the wind exceeded the threshold, it built a genuine hybrid matrix subject to the hybridization constraint. The difference in production between the two optimal solutions was 6,773 kWh/day.

This contrast constitutes the central validation of the work: applied to two territories with opposite resources, the model responded consistently in both, adapting the renewable composition to local conditions without external intervention or manual reparameterization. It did not impose a predetermined configuration, but rather let each zone's resource determine the solution, which empirically demonstrates the generalization of the approach proposed in Section 1 beyond a single calibration territory.

4. Discussion

The results obtained can be interpreted considering the reviewed literature, which confirms their coherence while also clarifying their scope. The inverse altitudinal gradient identified in Villavicencio, with irradiance decreasing and wind speed increasing with elevation, is consistent with the control exerted by relief on both resources, as documented for statistical and machine learning-based downscaling of the wind resource in complex terrain [36] and for topography-conditioned photovoltaic potential estimation, also including LiDAR technology [37,38,39]. This agreement supports the relevance of having assessed the resource separately by stratum, rather than assuming a single value for the municipality, even though, in this case, stratification did not alter the final operational conclusion regarding wind viability. This stratification approach represents a baseline methodology. As such, each stratum was characterized using a single representative point, and the terrain attributes used for the wind component were derived from the municipal area rather than from the specific boundaries of each stratum. Previous studies have shown that greater spatial detail can improve representation in regions with heterogeneous renewable energy potential, whereas aggregated approaches remain suitable for preliminary planning and prefeasibility analyses [40]. Therefore, the adopted stratification is considered sufficiently representative for the purposes of this study while preserving the possibility of further refinement in subsequent stages of analysis.

The wind non-viability obtained across the four strata of Villavicencio must be interpreted within the framework of the data source used. The literature repeatedly warns that reanalysis products underestimate or incompletely represent wind in complex terrain, and that extrapolation to hub height introduces additional uncertainties [41,42,43]. The result is, therefore, robust within the framework of the satellite data used, but does not equate to a direct measurement, which reinforces the recommendation for local instrumentation prior to any definitive wind design in the region. For low-wind contexts such as the foothills, moreover, bladeless wind turbine configurations based on vortex-induced vibration have recently been explored, which operate efficiently at speeds lower than those required by conventional horizontal-axis technology [44]; this alternative, identified as a secondary viability threshold in Section 2.c but not reached in any stratum of Villavicencio, defines a concrete direction for future work in the territory.

The behavior of the model in Manaure, where it autonomously selected a hybrid matrix with a 73% wind share, is consistent with optimization frameworks that integrate spatial selection and technology allocation in territories with a dominant wind resource [45]. Due to the fact that the model derived this composition without imposing it confirms that the hybridization constraint operates as an effective diversification mechanism rather than a rigid rule: below the minimum area required by a single turbine, the model maintained the pure solar solution, and only from that area onward did it activate the mix, a non-trivial behavior that supports the formulation of Equation (5) over fixed-proportion alternatives. This hybrid configuration also achieved a higher energy density than the best solar solution of Villavicencio, which is consistent with the literature on renewable resource complementarity, documenting how combining sources with complementary generation patterns enables more efficient use of the available area compared to single-technology installations [46]. For the Colombian context specifically, the Caribbean Coast, where Manaure is located, has been reported to exhibit particularly favorable complementarity potential between solar and wind resources [47], which is consistent with the result obtained by the model.

The fidelity of the surrogate model, with coefficients of determination above 0.99 across all zones evaluated, is in line with the computation time reductions and accuracy levels reported for machine learning-based surrogate models applied to energy system optimization [21], and for response surfaces used in uncertainty quantification of CFD-based wind site assessments [48]. This agreement supports the methodological decision to delegate the nonlinearity of the resource to the surrogate while preserving the accuracy of the optimization in the MILP formulation, rather than also approximating the optimization stage itself, as other approaches in the literature do by replacing the solver itself with a learned model [21].

The proposed model is limited to maximizing production and energy efficiency over a discrete set of capacities, without incorporating local demand, storage, or economic costs. This scope delimitation is consistent with the literature on hybrid system sizing, which identifies these three dimensions as determinants of a project's ultimate feasibility, beyond resource availability alone [49]. Incorporating them would allow the model to move from a resource-optimal configuration toward one that is also optimal in terms of demand-matching capability and levelized cost of energy, which constitutes a natural extension of the model toward a subsequent detailed design stage.

Before detailed system design becomes relevant, the model fulfills a screening function by reducing the set of candidate technologies to those supported by local resource conditions. This role is consistent with early-stage renewable energy decision-support frameworks, whose objective is not detailed engineering design but the efficient allocation of analytical effort toward the most promising alternatives [50]. In this context, the exclusion of a technology on resource grounds constitutes an informative outcome, as it narrows the feasible solution space and reduces the need for further evaluation of non-competitive configurations.

Finally, the results highlight dimensions that fall beyond the scope of the model but are identified in literature as relevant to territorial planning. Energy density, used here as an efficiency criterion, coincides with concerns about efficient land use in the deployment of distributed generation [51], and the implementation of large-scale photovoltaic plants entails measurable environmental effects on the surrounding ground surface temperature, documented even in arid climates comparable to foothills [52], which should be evaluated in an eventual implementation phase. In this direction, digital twin frameworks and AI-assisted spatial intelligence offer a pathway to integrate these environmental considerations, together with continuous resource updating, into energy planning decision-making [53,54,55,56,57], which opens a natural avenue for extending the proposed model.

5. Conclusions

The proposed four-stage model geomorphological stratification, resource characterization, simulation coupled with a surrogate model, and MILP optimization subject to a hybridization constraint demonstrated that explicitly integrating the geomorphology of the territory with a formal technology-mix constraint allows a single methodological workflow to autonomously adapt the solar-wind composition to local conditions, without manual reparameterization between territories. Faced with opposing resource profiles, the model yielded an exclusively solar solution where wind did not reach the operating threshold, and a wind-dominated hybrid matrix where it did, evidencing that the response of the model is governed by the local resource rather than by a predetermined configuration.

The main contribution of this work is to propose geomorphological stratification as an organizing criterion that articulates, within a single reproducible workflow, resource characterization, a validated surrogate model and optimization; built from globally available satellite data and therefore extensible to other territories without depending on prior local instrumentation. The model not only delivers the optimal configuration for each zone, but also explicitly documents when local resource conditions do not support a given technology, a property that makes it valuable both for guiding investment decisions and for ruling out options that lack a foundation in the local resource base.

The framework developed in this study provides a basis for several avenues of future research. These include refining resource characterization through higher-resolution spatial representations within each stratum, validating wind resource estimates with direct field measurements, and expanding the optimization model to account for demand dynamics, storage systems, and economic performance metrics. Integrating these elements would allow the transition from a baseline planning tool toward a more detailed design and decision-support framework.