Metaheuristic-Enhanced Multi-Layer Perceptron for Early-Stage Cooling Load Prediction of Courtyard Buildings in Tropical Climates

The building sector is one of the largest consumers of energy worldwide, accounting for approximately 40% of total primary energy consumption in many developed and developing nations [1]. In tropical climates, this issue is particularly pressing because cooling systems operate year-round due to persistently high ambient temperatures and humidity [2]. The increasing urbanization of Southeast Asian cities, including those in Malaysia and Indonesia, has accelerated the construction of buildings, many of which incorporate traditional courtyard configurations as passive design strategies [3]. Courtyards are well established as architectural elements that moderate indoor thermal conditions by promoting natural ventilation, shading, and daylighting [4]. However, the relationship between courtyard morphology and building energy performance is complex and nonlinear, making accurate cooling-load prediction challenging at the early design stage.

Conventional methods for estimating cooling load rely on detailed energy simulation software such as EnergyPlus, DesignBuilder, and IES Virtual Environment (IESVE) [5]. These tools require complete specification of building geometry, envelope properties, and occupancy schedules, which are rarely available during the conceptual design phase [6]. As a result, energy performance is often evaluated only after major design decisions have been finalized, limiting opportunities to implement effective passive cooling strategies. The disconnection between energy assessment and early-stage design has been widely recognized as a critical barrier to building energy efficiency improvement [7].

Machine-learning models offer a useful surrogate for this early-stage problem because they can learn nonlinear relationships between a limited set of design variables and energy outcomes. Artificial neural networks (ANNs), support vector regression (SVR), random forests (RF), and decision trees (DT) have all been applied to building-energy prediction [8]. Among these, the Multi-Layer Perceptron (MLP) is widely applied due to its universal approximation capability and suitability for capturing nonlinear relationships between input features and energy outputs. However, MLP predictive performance is sensitive to hyperparameter selection, including the number of neurons per hidden layer, the regularization coefficient, and the learning rate. Suboptimal hyperparameter settings can lead to underfitting or overfitting, reducing the reliability of predictions [9].

Metaheuristic optimization provides one way to search for this hyperparameter space without relying on gradient information [10]. These algorithms are population-based or trajectory-based search methods that explore the solution space without requiring gradient information, making them suitable for optimizing non-differentiable and multimodal objective functions [11]. Several studies have integrated mathematical algorithms with ANN or MLP models for building energy prediction. Zheng et al. [12] applied a shuffled complex evolution algorithm to optimize an MLP for cooling load prediction and reported improvements in learning and prediction accuracy exceeding 19% compared to a standard MLP. Afzal et al. [13] compared multiple metaheuristic optimizers including Satin Bowerbird Optimizer, Ant Lion Optimizer, Slime Mold Optimizer, and Particle Swarm Optimization (PSO) for predicting heating and cooling loads using SVR and extreme gradient boosting. Mu et al. [14] demonstrated that metaheuristic-optimized MLP models trained with systematic data preprocessing achieved lower prediction errors compared to models without preprocessing, highlighting the importance of both algorithm selection and data quality. Despite these advances, the literature has not reported a comparative assessment of multiple algorithmically diverse methods applied to courtyard building typologies in tropical climates.

This study addresses this gap by evaluating five metaheuristic algorithms from the evolutionary and swarm-based optimization categories to optimize MLP hyperparameters for cooling-load prediction. This study investigates three courtyard building shapes: L-shape, O-shape, and U-shape. These configurations represent common courtyard typologies found in tropical regions and differ in enclosure degree, spatial geometry, and exposure to solar radiation. The dataset was generated from parametric energy simulations conducted in IESVE, covering 8,535 combinations of building morphology and envelope variables. The objective function used during optimization is the mean cross-validated RMSE on the training set, ensuring that the selected hyperparameters generalize well to unseen data. Each optimized model evaluated on an independent test set using seven statistical metrics, and the statistical significance of differences between algorithms was assessed using the Friedman test across multiple independent runs.

The main contributions of this study are as follows:

(1) This work provides a comparative assessment of five metaheuristic search strategies for tuning an MLP-based cooling-load predictor for tropical courtyard buildings under a common experimental framework.

(2) This work evaluates optimizer behavior using multiple performance measures, convergence information, and repeated stochastic runs rather than relying on a single optimization run.

(3) This examines the influence of courtyard form in a tropical context and identifies an important limitation in the fully enclosed O-shaped subset, thereby clarifying where the proposed surrogate is currently most dependable and where additional design variables are needed.

3.1 Dataset characteristics

Table 4 summarizes the dataset. Annual cooling-load intensity ranges from 43.4962 to 49.2145 kWh/m², with a mean of 46.0419 kWh/m² and a standard deviation of 1.2616 kWh/m². The relatively narrow target range reflects the constrained simulation space considered in this study. The courtyard forms also differ in geometric descriptors such as S_V_RATIO and FORM_FACTOR, which is expected because enclosure changes the relationship between exposed surface area and building volume.

Table 4. Descriptive statistics of the dataset

3.2 Convergence analysis

Figure 2 presents the convergence curves of the five metaheuristic algorithms, showing the best cross-validated RMSE achieved on the training set as a function of iteration number, averaged across five independent runs. All five algorithms show a general decreasing trend in the objective function value across iterations, confirming that each algorithm improves the MLP hyperparameter configuration over successive iterations.

2.png

Figure 2. Convergence curves of the five metaheuristic algorithms

GA converged fastest, reaching a near-stable plateau by approximately iteration 3–4. However, this rapid convergence came at a cost: GA achieved the smallest total reduction in objective value across the full search (0.00081), compared to 0.00198 for DE. This pattern suggests premature convergence, where the population loses diversity too early, and the search settles into a limited region of the hyperparameter space before adequately exploring alternatives. This behavior is consistent with the higher run-to-run variability observed for GA on the test set, discussed in Section 3.3.

In contrast, DE, PSO, and HHO continued to improve steadily well into the later iterations, with DE still showing measurable gains as late as iteration 18. Averaged across the five runs, DE achieved the lowest mean training RMSE (0.1472), followed closely by HHO (0.1472) and PSO (0.1473); GWO followed with 0.1474, while GA showed the highest mean training RMSE (0.1480). Table 5 separately reports the hyperparameter configuration of the single best-performing run per algorithm (selected by test R²), which does not necessarily correspond to the lowest training RMSE, since the two are optimized on different criteria, cross-validated training RMSE during search versus final test-set accuracy after refitting.

Table 5. Training CV-RMSE and hyperparameter configuration of the best solution per algorithm (mean ± std over 5 runs)

3.3 Prediction performance on test set

Table 6 presents the mean and standard deviation of all seven-evaluation metrics across five independent runs for each algorithm on the independent test set. The results indicate that GWO achieved the highest overall prediction accuracy, with mean R² of 0.5836 ± 0.0031, RMSE of 0.8106 ± 0.0030 kWh/m², MAE of 0.5821 ± 0.0068 kWh/m², MAPE of 1.2558 ± 0.0145%, NSE of 0.5836 ± 0.0031, and Willmott’s d of 0.8540 ± 0.0044. These values indicate that the GWO-optimized MLP predicts cooling load with the best overall accuracy among the algorithms evaluated, although with modest absolute explanatory power (R² ≈ 0.58), a limitation discussed further in Section 3.7.

Table 6. Test set performance metrics for each algorithm (mean ± std over 5 runs)

GWO and PSO consistently ranked among the top performers across most metrics. Both algorithms are swarm-based, meaning each candidate solution updates its position using information from the best solutions found elsewhere in the population. This shared-information mechanism likely contributed to their stronger, more consistent performance than the evolutionary algorithms (GA, DE) tested in this study.

GA showed the highest run-to-run variability in test R² (std = 0.0153), more than four times that of DE (std = 0.0024). DE was the most stable and reproducible algorithm across independent runs. Interestingly, GA and DE are both evolutionary algorithms, yet they sit at opposite ends of the stability range. Algorithm category alone does not explain this difference. A more likely cause is GA’s crossover and mutation operators, which tend to produce less consistent exploration across different random initializations than DE’s difference-vector mutation strategy.

Figure 3 shows the boxplot distribution of R² across five runs for each algorithm. DE had the narrowest interquartile range, again confirming its stability. GA had the widest spread, suggesting its final solution quality depends more heavily on the search’s random initialization.

3.png

Figure 3. Boxplot distribution of coefficient of determination (R²) across five independent runs for each algorithm

Figure 4 presents the scatter plots of predicted versus observed cooling load values for the best-performing algorithm (GWO) on the training and test sets. The predicted values closely follow the diagonal line of perfect agreement, with no systematic over-prediction or under-prediction evident across the range of observed values.

4a.png

4b.png

Figure 4. Predicted versus observed cooling load using the best-performing algorithm (GWO): (a) training set, (b) test set

Figure 5 presents the residual analysis for the GWO model, showing residuals plotted against predicted values and the corresponding residual distribution. The residuals are approximately symmetric around zero with no strong systematic trend across the range of predicted values, indicating the absence of major systematic bias in the model predictions.

5a.png

5b.png

Figure 5. Residual analysis for the best-performing algorithm (GWO): (a) residuals versus predicted values, (b) residual distribution

3.4 Feature importance analysis

To identify which design variables strongly influence cooling load predictions, permutation feature importance was computed for the best-performing model (GWO) on the test set. Table 7 and Figure 6 present the ranked importance of each input feature, measured as the mean decrease in R² when the feature values are randomly permuted.

Table 7. Permutation feature importance (ΔR²) for the best-performing algorithm (GWO), ranked by mean importance across 10 permutation repeats

6.png

Figure 6. Permutation feature importance (ΔR²) for the best-performing algorithm (GWO)

The results show that WINDOWS_RATIO overwhelmingly dominates prediction accuracy, with permutation importance more than an order of magnitude higher than any other feature. Because R² is not bounded below, permuting a dominant feature can drive model predictions below a naive mean-based baseline, producing a permuted R² lower than zero. For Window Ratio, the permuted R² is approximately -0.4330 (0.5836 minus 1.0166), confirming that a ΔR² exceeding the baseline R² is a mathematically expected outcome of the permutation procedure rather than a computational error. The dominance of window-to-wall ratio observed here is consistent with established building physics and prior simulation-based studies. WWR directly governs solar heat gain through glazing, which is the primary driver of cooling demand in tropical and hot-humid climates. Li et al. [18] demonstrated a near-linear increase in cooling-energy demand as WWR increased from 0% to 100% in a coupled microclimate-building energy simulation.

3.5 Statistical comparison

The Friedman test was applied to the RMSE values of all five algorithms across five independent runs to assess whether the observed performance differences are statistically significant. The test yielded a Friedman statistic of 13.44 and a p-value of 0.0093. Since the p-value is lower than the significance level of 0.05, the null hypothesis that all algorithms perform equally is rejected. This result indicates that the observed differences in prediction accuracy among the five algorithms are statistically significant, and that the ranking of algorithms reported in Table 6 reflects genuine performance differences rather than random variation. Post-hoc Wilcoxon signed-rank tests comparing GWO against each of the remaining algorithms showed a consistent directional advantage for GWO, which outperformed DE, HHO, and PSO in all five paired runs (p = 0.0625 for each comparison). However, these pairwise comparisons did not reach the conventional significance threshold individually, a limitation attributable to the maximum achievable precision of the exact Wilcoxon test at n = 5 rather than to a weak underlying effect.

3.6 Implications for early-stage building design

The results show that metaheuristic-optimized MLP models can reliably predict cooling load for open courtyard configurations at the early design stage. These configurations comprise L-shape and U-shape. Predictive reliability for the fully enclosed O-shape configuration remains constrained. The dataset limitations discussed below explain this constraint. The model requires only basic morphological and envelope parameters. This information is typically available before detailed design documentation is complete. The best-performing model was optimized using GWO. This model achieved an R² of 0.5836 on the independent test set. It achieved an RMSE of 0.8106 kWh/m² on the same set. The selected input features explain roughly 58% of the variance in cooling load.

These findings carry practical value for architectural design workflows. A trained MLP model can be embedded in early-stage design decision support tools. Designers can use these tools to quickly screen open courtyard configurations. Designers can identify configurations with lower cooling demand before committing them to a specific design direction. This approach is faster than running a full energy simulation for every design variant. It does not require specialist energy modeling expertise from the designer.

This study examined three courtyard shapes. These shapes are L-shape, O-shape, and U-shape. The shapes differ in enclosure degree and self-shading capacity. Both factors influence solar heat gain within the courtyard space. Both factors also influence natural ventilation within the courtyard space. O-shape is the fully enclosed configuration. O-shape recorded the highest mean cooling load at 46.22 kWh/m². U-shape is the more open configuration. U-shape recorded the lowest mean cooling load at 45.61 kWh/m². This result runs counter to a common assumption. The assumption holds that greater enclosure reduces cooling demand through increased self-shading. A plausible explanation involves ventilation. Full enclosure also restricts natural ventilation within the courtyard. This restriction traps heat inside the courtyard. The trapped heat would otherwise dissipate through cross-ventilation in more open configurations such as U-shape. The dataset is also imbalanced across shape. O-shape comprises 64% of all samples. This imbalance may itself contribute to the wider spread of cooling load values observed for O-shape. The standard deviation for O-shape reaches 1.30. The standard deviation for both L-shape and U-shape reaches only 1.11. This finding suggests an interaction between courtyard enclosure and ventilation performance. Solar exposure alone does not capture this interaction. Section 3.7 makes the case for incorporating airflow and ventilation-related variables in future extensions of this dataset. This finding reinforces that case.

Per-shape evaluation of the best-performing model on the test set clarifies this imbalance further. L-shape is predicted with high accuracy. L-shape reaches an R² of 0.9913 and an RMSE of 0.1042 kWh/m². U-shape is also predicted with high accuracy. U-shape reaches an R² of 0.9883 and an RMSE of 0.1213 kWh/m². O-shape is predicted with substantially lower accuracy. O-shape reaches an R² of only 0.3892 and an RMSE of 1.0048 kWh/m². This disparity traces to the underlying dataset rather than to model capacity or training imbalance. The O-shape subset contains 1,514 unique combinations of the seven geometric input features. Among these, 1,481 combinations correspond to multiple samples with differing cooling load values. The pooled within-combination variance accounts for approximately 79% of total O-shape target variance. This finding is consistent with the correlation pattern between window-to-wall ratio and cooling load. This correlation reaches 0.97 for both L-shape and U-shape. This correlation reaches only 0.55 for O-shape. These results indicate a missing geometric degree of freedom specific to the fully enclosed configuration. The internal courtyard aperture dimension most plausibly represents this missing factor. The present input features do not capture this factor. This gap imposes a data-driven ceiling on achievable predictive accuracy for O-shape. This ceiling is independent of model architecture or hyperparameter tuning. Predictive claims regarding this MLP-based approach should therefore be scoped to open courtyard configurations. These configurations comprise L-shape and U-shape. This scope should hold until an aperture-related predictor can be incorporated for the fully enclosed case.

The absolute range of cooling load in this dataset is also narrow. This range spans 43.50 to 49.21 kWh/m². The standard deviation reaches 1.26. This narrow range constrains the practical magnitude of difference an R² of 0.58 can represent. In absolute terms, however, the best-performing model achieves an RMSE of 0.81 kWh/m². This RMSE corresponds to approximately 14% of the observed range. This magnitude remains practically informative for comparative design screening. R² itself appears modest despite this. Differences in R² at the third decimal place between algorithms should therefore be interpreted carefully. GWO reaches 0.5836. PSO reaches 0.5772. These differences indicate relative algorithm consistency. These differences do not indicate meaningfully different real-world prediction accuracy. Section 3.5 reports the Friedman and Wilcoxon tests. These tests address this distinction directly. These tests examine rank consistency across runs. These tests do not examine the practical magnitude of difference between algorithms.

The best hyperparameter configurations identified for several algorithms also reached the boundary of the defined search space. GWO’s optimal number of neurons in the first hidden layer reached 128. GWO’s optimal learning rate reached 0.01. Both values correspond to the upper bounds reported in Table 2. This raises the possibility that better-performing configurations exist outside the explored range. The present search already incurred substantial computational cost. Section 2.6 describes this cost. The search space was not extended further in this study. This saturation is instead noted here as a limitation. Future work should verify whether relaxing these bounds yields additional accuracy gains.

3.7 Beyond window-to-wall ratio: Toward a more comprehensive cooling load model

In general, cooling loads in tropical regions are strongly influenced by design and non-design factors. Non-design factors such as occupancy level, occupancy behavior, and building function can only be determined effectively during operation. However, design factors such as window design, building facades, and material types can influence architects’ design decisions [19]. This decision is especially insightful when supported by early cooling-load estimates before construction starts. It also aligns with this study’s finding, as shown in Figure 6, that the most important feature influencing courtyard building cooling-load prediction is the WWR. Because design factors are so important, future studies should explore combinations of design factors, including WWR ratios, glazing types and specifications, building envelope properties, color, and shading, to generate lower cooling-load predictions.

In addition, extending the prediction with hourly historical climate and outdoor temperature data measured specifically in the Shah Alam microclimate will provide more accurate estimates and support solid design decisions. Extending the dataset in this direction would allow future metaheuristic-optimized models to unravel the relative contributions of morphological design decisions (the focus of the present study) from envelope material to climate-driven factors. This will provide architects and building designers in regions with comparable tropical contexts with a more complete, physically grounded basis for early-stage cooling load estimation. Furthermore, direct comparison between the courtyard-shape-driven design strategies and retrofits involving building’s envelope-level interventions (e.g., low-SHGC glazing, external shading) is not possible within the scope of the present dataset. Thus, metaheuristic-optimized MLP models can help support informed design decisions in the future.

3.8 Economic implications of metaheuristic-optimized early-stage prediction

Beyond its technical contribution, the practical value of a metaheuristic-optimized MLP model for cooling load prediction lies in its potential to reduce the cost and time associated with early-stage building energy assessment. Detailed dynamic energy simulation using tools such as IESVE typically requires substantial input preparation, computational runtime, and specialist expertise, incurring direct consultancy costs that are often difficult to justify for architects and developers during the conceptual design phase, when multiple design alternatives are still being compared under tight budget and schedule constraints. The trained MLP model developed in this study produces a cooling load estimate in a fraction of a second once optimized, compared with the substantially longer setup and computation time typically required for a single detailed simulation run in dynamic energy modeling software [20]. This order-of-magnitude reduction in evaluation time and required expertise lets developers and design teams screen far more courtyard configurations within the same design timeline and consultancy budget, potentially reducing the direct cost of the design exploration phase and enabling earlier identification of lower-cooling-load design options before capital-intensive detailed design work begins. Because cooling load is directly proportional to ongoing HVAC operating expenditure over a building’s service life, design decisions informed at this early, low-cost screening stage carry disproportionate downstream economic value relative to their marginal assessment cost, as design changes become progressively more expensive to implement as a project moves toward construction documentation.