This section develops the SA-MCGM in mathematical form, proceeding from the base grey model through the Markov-chain extension to the simulated-annealing optimization layer, with the underlying equations and reasoning presented at each stage.

The (GM(1,1)) has three main stages: Summation operator, Inverse addition operator and forecasting.

Consider or the initial time series (the real data) with n samples in Vector (1).

where, is the time series value at ith time point. Using the Accumulated Generating Operator (AGO), the regular series is converted to the following time series in a monotonically increasing manner.

At this stage, can be calculated using the inverse operator with Eqs. (3) and (4).

The grey model, GM(1,1), is formed using the First-degree differential equation as below:

This method does not have any mathematical proof and is called the whitening process at which i is a time point, a is improvement coefficient and b is the controller coefficient. The data of the GM(1,1) are shown in matrix form in Relationships (8)-(10).

(10)

Using the least squares method, the parameters of the model (a and b) are calculated as Eq. (11) [24].

(11)

After calculating a and b, grey forecasting equation such as Eq. (12) [25] is obtained.

(12)

So series is smoothened and will be the forecast series which calculated by Eq. (12).

While GM(1,1) captures general trends effectively, its accuracy deteriorates when the underlying data exhibit sharp, non-linear fluctuations, since the model has no mechanism for representing stochastic deviation from the deterministic trend it fits. The MCGM addresses this by treating the residual error between GM(1,1) predictions and observed values as a Markov process, allowing the model to adapt to the non-stationary behavior characteristic of energy consumption data.

The residual at each time point, quantifying the deviation between the GM(1,1) forecast and the observed value, is calculated as in Eq. (13).

(13)

The next stage is incorporating Markov chains for residual adjustment. The Markov chain component analyzes and adjusts the residuals, enhancing the accuracy of predictions by incorporating probabilistic transitions. Residuals e(i) are partitioned into r states represented as Sj where j=1, 2, …, r. For each state, lower (Lj) and upper (Uj) bounderies are defined representing distinct ranges of error magnitudes. These boundaries are determined based on the distribution of e(i), ensuring the states capture meaningful variations by Eqs. (15) and (16).

(14)

(15)

(16)

The probability of transitioning from one state Si to another Sj is represented in the transition matrix (P). The elements of transition matrix ( ) are calculated by Eq. (17) where m represents the probability of transitioning from state i to state j after m steps. In this context, denotes the number of data points associated with the transition from state i to state j after m steps. is the data which belongs to the ith state.

(17) Consequently, the matrix of transition probabilities among states is given by Eq. (18).

(18)

The Markov adjusted residual is added to the values obtained from Eq. (12) (GM(1,1) forecasted values) in order to improve the accuracy of the forecasting results which are calculated through Relation (19):

(19)

In which eM(i) is the Markov-adjusted residual for time i and is calculated as Eq. (20).

(20)

In Eq. (20), count(e(i) ∈ Sj) is the count of e(i) refers to the number of residual values e(i) that belong to the state Sj . Specifically, the count is the number of times the residual values fall within the state .

Fig. 4 presents total energy consumption in Asia from 2000 to 2022, as reported by the World Energy Agency. Consumption rose from 59,584.29 TWh in 2000 to 121,234.56 TWh in 2022, nearly doubling over the 22- year period, with a minor dip in 2020 coinciding with the COVID-19-related economic slowdown. The corresponding energy mix, reported by the International Energy Agency, reflects the region's gradual diversification across fossil fuels, nuclear power, and renewables such as wind, solar, and hydropower.

A sensitivity analysis was conducted on the three SA hyperparameters, the initial temperature ,the final temperature and the number of iterations (steps), evaluating MSE on a held-out test set across a predefined grid of values. Results are shown in Fig. 5. The sensitivity analysis was performed over a set of predefined values:

For each combination, the proposed optimizer model was executed to adjust the parameters of the SA. The MSE Performance metric was computed and shown in Fig. 5 on a held-out test set to assess forecasting

The experimental analysis in Fig. 5 reveals that optimal SA performance hinges on strategic parameter interplay. Higher values (e.g., 100,000) promote exploration by escaping local optima early, but require pairing with sufficiently low (e.g., 1×10-6) to enable precise exploitation in later stages. Moderate iteration counts (Steps ≈ 5,000) balance computational efficiency and convergence quality, as excessive iterations yield diminishing returns, while insufficient steps hinder refinement. For instance, configurations like = 25,000, = 1×10-6, and Steps = 10,000 achieved peak performance, whereas mismatched parameters (e.g., high with inadequate ) led to instability. The findings advocate a calibrated hierarchy: Prioritize exploration with elevated , enforce precision through minimized , and align Steps with problem complexity to optimize both solution quality and computational practicality.

The implementation of the SA-MCGM incorporated carefully calibrated parameters to optimize its forecasting performance. After conducting the optimization using SA, the model achieved optimized parameter values that resulted in significantly improved forecasting accuracy. Under the optimal configuration, the development coefficient a was found to be approximately -0.0146, the grey action quantity b was approximately 83889.3901, and the number of Markov states r was set to 4. These values were determined through iterative refinement, which ensured that the model parameters were precisely tuned to minimize forecasting error.

Table1 compares the results of the SA-MCGM with the traditional MCGM, highlighting the differences in their respective parameter values and the objective function MSE. The development coefficient a in the SA-MCGM is slightly improved over the MCGM value of −0.0397. The grey action quantity b in SA-MCGM is 83,889.39, a 39.37% increase compared to MCGM’s 60,194.88. A key difference between the models is the number of Markov states r. The SA-MCGM optimized r to 4 states, compared to 8 in MCGM, reflecting a decrease of 50%. This adjustment is expected to contribute to improved model performance by allowing for a more granular error correction process. The most significant difference lies in the MSE, where the SA-MCGM achieved an MSE of 3818842, a remarkable reduction of approximately 98.17% compared to the MCGM's MSE of 208622806.8. This substantial improvement in MSE confirms the effectiveness of the SA optimization process in refining the model's forecasting capabilities. These results demonstrate that the SA-MCGM, with its ability to more effectively explore the parameter space, offers a significant performance advantage over the traditional MCGM, especially in terms of reducing forecasting error.

To assess the performance of the proposed model, Table 3 compares the actual energy consumption values with those predicted by different methods across both training and testing datasets. The statistical method ARIMA was employed to forecast the actual energy consumption values. Following the analysis of autocorrelation and partial autocorrelation function plots, an ARIMA model of order one (AR(1)) was implemented on the data. The column GM(1,1) shows the predicted annual energy consumption values using the first-order grey model with the parameters listed in Table 1. After applying the MCGM method to the actual values, predictions were calculated using this method. For the MCGM approach, a trial-and-error method was utilized to determine the optimal value of r that minimizes the MAPE, ensuring the highest prediction accuracy. This iterative process highlights the importance of parameter tuning in achieving reliable and precise forecasting results. Utilizing the parameters from the SA-MCGM method in Table 1, the annual energy consumption values in Asia were computed using the proposed method.

Following the precise calculation of the Average Percentage Error (APE) for each method from 2000 to 2022 in two different periods, the percentage of absolute deviation is illustrated in Fig 6. The analysis highlights SA-MCGM as the most robust model for long-term forecasting. While it initially shows high errors (e.g., 42.56% in 2002), its adaptive mechanisms drive significant improvement over time, achieving near-perfect accuracy (0.007% in 2022) during testing. In contrast, ARIMA exhibits moderate but inconsistent performance (0.14–10.56%), while GM(1,1) and MCGM degrade sharply in later years (up to 21.06% and 18.31%, respectively). SA-MCGM’s ability to minimize errors over extended horizons, coupled with its superior test-phase stability, positions it as the optimal choice for scenarios requiring long-term precision, despite initial calibration challenges.

To better compare the methods, two Analysis of Variance (ANOVA) test were performed. The results of the ANOVA on the APEs indicate that the performance of all four methods differs at the significant level of 0.000 in two datasets (learning and testing data). The results of ANOVA as box plots of APE to compare the algorithms are shown in Fig. 7.

The comparative performance of the prediction methods is further illustrated in the side-by-side boxplots of the APE for the training and testing datasets. As evident from the plots, the proposed SA-MCGM algorithm demonstrates superior performance for both datasets compared to the other methods, ARIMA, GM(1,1), and MCGM. In the training data, SA-MCGM shows higher variability and an extreme outlier (APE of 42.56), but in the testing data, it exhibits consistently low APE values (ranging from 0.007 to 4.28, with a median of approximately 0.57), indicating remarkable stability and accuracy. The other performance metrics, including MAPE, Mean Absolute Error (MAE), and MSE, are presented in the table below for the two periods of training and testing data. The analysis results highlight a significant relative improvement of the SA-MCGM model compared to the ARIMA, GM(1,1), and MCGM prediction methods, particularly in the testing phase, where it outperforms the others with lower and more consistent errors.

The performance metrics presented in Table 3 offer a comprehensive comparison of the forecasting models GM(1,1), MCGM, and SA-MCGM, highlighting their efficacy across training and testing phases. GM(1,1) demonstrates a low training MAPE of 2.74%, yet its testing MAPE escalates to 13.53%, with MSE surging from 7,405,457 to 302,223,168, indicative of overfitting and poor generalization to unseen data. Similarly, MCGM, while maintaining a modest training MAPE of 3.22%, sees a testing MAPE of 10.63% and an MSE increase to 208,622,844, suggesting a moderate overfitting tendency that compromises its predictive reliability. In stark contrast, SA-MCGM exhibits a training MAPE of 13.48% and an elevated MSE of 189,841,068, which might initially suggest underfitting; however, its testing performance is exceptional, with MAPE plummeting to 1.04% and MSE reducing to 3,818,842. This remarkable improvement underscores SA-MCGM’s superior generalization capability, a critical attribute in forecasting applications. The success of SA-MCGM can be attributed to its hybrid methodology, which integrates SA to optimize the grey forecasting framework. This approach appears to mitigate the overfitting pitfalls observed in GM(1,1) and MCGM by prioritizing robustness over excessive fit to training data, as evidenced by its tightly clustered testing errors. Furthermore, the MAE metrics reinforce this conclusion: While GM(1,1) and MCGM show substantial deviations in testing (15,712 and -1,111, respectively), SA-MCGM’s testing MAE of 615 reflects minimal absolute error, affirming its precision. These results align with contemporary trends in forecasting literature, where hybrid models leveraging optimization techniques consistently outperform traditional approaches by balancing complexity and adaptability.

Although GM(1,1) and MCGM offer competitive training performance, their inability to maintain accuracy in testing renders them less reliable for practical applications. SA-MCGM, by contrast, excels where it matters most, delivering dependable forecasts on new data. This robustness, coupled with its methodological sophistication, positions SA-MCGM as the preeminent choice among the evaluated models, offering a compelling case for its adoption in scenarios demanding high predictive accuracy and stability.

The projected values for annual energy consumption in Asia from 2023 to 2035 have been forecasted using GM(1,1), MCGM, and SA-MCGM methods, utilizing the optimal parameter values calculated in previous sections. These projections are illustrated in Fig. 8. It is anticipated that the energy consumption will reach approximately 144780.04 units by the year 2035 using the proposed method.

The analysis of the forecasting models GM(1,1), MCGM, and SA-MCGM, spanning the period from 2023 to 2035, reveals statistically significant disparities in their predictive outcomes. A one-way ANOVA yielded an F-statistic of 71.19 with a p-value effectively approaching zero (p < 0.001), indicating robust evidence to reject the null hypothesis of equality among the group means. Examination of the temporal trends demonstrates that GM(1,1) exhibits the most pronounced growth, projecting a value exceeding 245,000 by 2035, while MCGM follows a comparable trajectory, albeit with a slightly tempered increase, culminating at approximately 242,000. In contrast, SA-MCGM displays a markedly conservative profile, stabilizing at around 144,000 with minimal incremental change over the study period. This divergence is further corroborated by the box plot analysis (Fig. 9), which illustrates greater variability and higher median values for GM(1,1) and MCGM, indicative of their expansive forecasting ranges, whereas SA-MCGM presents a narrower interquartile range and lower median, suggesting a more constrained predictive scope. These findings underscore significant methodological differences with potential implications for applied forecasting contexts, highlighting the necessity for careful model selection in research and policy applications where projected outcomes may substantially influence strategic planning and resource allocation.

This divergence has a direct practical implication: Because GM(1,1) and MCGM extrapolate the historical upward trend largely unchecked, they imply substantially higher future infrastructure and supply requirements than SA-MCGM, which instead reflects the flattening pattern already visible in the most recent testing-period data (2016–2022). Given SA-MCGM's demonstrated testing-period accuracy in Section 3.4, its more conservative long-range projection is the more credible basis for capacity and policy planning, though the magnitude of the gap between models underscores the sensitivity of long-horizon forecasts to model choice.