The TOP500 list comprehensively represents global HPC infrastructure.

The energy consumption structure of HPC centers in each country follows the overall energy structure of that country.

HPC facilities consist primarily of CPU, memory and storage.

There will not be a sudden decrease in the demand for HPC caused by unpredictable factors, such as economic downturns.

This paper collected data on worldwide HPC centers’ P_max, N_c, and R_max from the TOP500 HPC list from 2014 to 2024, enabling the calculation of average power consumption (6, 7).

After arranging the HPCs in ascending order based on N_c, this study identified missing P_max values of some HPC centers. To address this issue, linear regression was implemented to model the relationship between the existing N_c and P_max values. Missing P_max values were estimated by substituting N_c into the derived function.

To address this issue, we implemented linear regression to estimate these values due to the strong linear relationships indicated by the Pearson correlation coefficients, as shown in Figure 3A.

The resulting regression function is as follows:

The R² of this regression function is as high as 0.813, which means that by using this function, we can fill in all missing P_max effectively.

To calculate the average utilization, our team refers to a formal essay that records the η_cpu, η_mem, η_stg of 50 HPCs in China in 2016 (9). The relevant equations, adopted from another research, apply to our calculation of η_avg. The calculations of P_cpu and P_mem were based on the models of another study (10).

The final equation utilized was:

The calculated η_avg for 50 Chinese HPCs in 2016 was 67.6%.

This study used a rating provided by a previous study to calculate the average utilization rate of the United States, Germany, and France in 2016 based on the utilization of China HPCs (11). To estimate the average utilization rate for HPCs in various countries other than these four countries, and in other years other than 2016, we use these four countries’ ratio of market size to the number of HPCs and the average utilization rate to plot a scatter graph (shown in Figure 3B) (12, 13).

We observed the relationship between the ratio and the rate. Because there was an upper and a lower limit of the utilization rate of HPCs, we use logistic fit instead of linear regression:

The logistic model helps us estimate the average use rate of different countries at any year using the ratio of Market size to HPC number. France, Germany, America, and China have distinct utilization rates because of available market size data. The market size of the rest of the countries was recorded as a whole, so we categorized them as “other” and calculated their market size per HPC by subtracting the HPCs in the four major countries from 500. The use rate of other was then calculated. Table 1 displays the rating and ratio used and the resulting utilization rate.

For each country group, energy consumption at the average utilization rate and at full capacity was computed. This study calculated and compared an HPC’s energy consumption at the average utilization rate by multiplying the EC_max of that facility by the average utilization rate in the country where that HPC facility was located. Global totals were derived by summing values for all TOP500 HPCs. Thus, 19 representative countries that possess HPC had been selected. This paper calculated the annual energy consumption of HPC facilities in each country in the last decade.

The changes of energy consumption of HPC facilities in each country during the last decade were shown in the following graphs, with separate graphs for the USA, Japan, and China (Figures 4A–F).

Global annual energy consumption ranges were:

Both EC_max and EC_avg were highest in 2017 and lowest in 2014.

Country-specific distributions (Figures 4G,H) showed:

When divided by countries, both EC_max and EC_avg were the highest in America and the lowest in Poland.

Our analysis employed EC_avg as the primary metric for representing the real-world energy usage patterns. Drawing upon energy structure data from 20 representative HPC-intensive countries (12, 13). This study categorized national energy portfolios into 5 distinct types: coal, natural gas, petroleum, nuclear energy, and renewable energy (14).

The calculation framework utilized the following notation:

Notably, α_r carried an emission factor of zero as renewable energy generation produces no direct carbon emissions. Thus, the studies’ calculations focused on the four non-renewable energy categories.

The total carbon emissions from HPC countries in 2016, 2019, and 2022 were as shown in Figure 5A. China and the United States exhibited carbon emissions significantly higher than those of other countries (China: 0.7681 billion tons in 2016, 0.651 billion tons in 2019, 0.4647 billion tons in 2022; United States: 0.8256 billion tons in 2016, 0.4884 billion tons in 2019, 0.4938 billion tons in 2022), while Japan’s carbon emissions were of the same order of magnitude (174.5 million tons in 2016, 93.55 million tons in 2019, 184.1 million tons in 2022). Generally, the total carbon emissions from HPC decreased in countries with larger economic scales and increased in countries with lower economic scales (e.g., Brazil: 0.626 million tons in 2016, 3.212 million tons in 2019, 11.61 million tons in 2022).

The Environmental Impact Index (EII) incorporated energy-specific sensitivity coefficients (Sensitivity _i ) that quantified the relative environmental impact of each energy type. This study derived the EII through a weighted summation of carbon emissions multiplied by their respective sensitivity coefficients, normalized by total emissions. This index provided a standardized measure of HPC’s environmental impact across different national contexts.

Figure 5B displayed the results of the EII calculation, which quantified the environmental impact of HPC in these countries for the three selected years. As the figure indicated, the EII of the United States and China was far higher than that of other countries, while Japan, Germany, France, and the United Kingdom followed closely.

This study quantified the economic consequences of HPC-related emissions using established climate-economy relationships. The model applied a conversion factor whereby each trillion tons of CO₂ emissions corresponded to 0.5% reduction in global GDP (15).

The total world GDP in 2016, 2019, and 2022 was 76.59 trillion dollars, 87.752 trillion dollars, and 99.67 trillion dollars, respectively. The economic loss caused by HPC in each HPC country is then calculated and shown in Figure 5C. The total economic loss due to carbon emission from HPC in 2016, 2019, and 2022 reached as high as 2,183,499 dollars.

To quantify the impact of the evolution of HPC over time on energy consumption, the peak theoretical performance of supercomputers (R_peak) was used as a key metric for HPC development. Data for R_peak from the TOP500 list during the years 2014–2023 were collected.

A quadratic regression model was adopted due to its strong fit and better alignment with recent trends. The general form of the fitted model was:

Subsequently, the predicted R_peak values were obtained for the years 2014 to 2030 using the fitted model. These values were used as input to the energy consumption projection model established in Section 3.2.3.

The regression analysis yielded the following quadratic equation for HPC performance:

This model achieved a coefficient of determination of R² = 0.9519, indicating a strong fit with the observed data and reflecting the decelerated growth of computing power in recent years.

In qualifying the extent to which increasing energy demand in other sectors would impact future energy consumption of HPC, the energy price was considered (14). In this respect, data of the global energy commodity price index, which was expressed as PI, from 2014 to 2023 were correlated and analyzed for their correlation with energy consumption EC_other in other sectors, which was derived via subtracting EC of HPC from the global energy consumption (16–19). Since energy price was highly interrelated with various social factors such as current affairs, data points with large fluctuations were intentionally removed. Based on the remaining data, a logistic fitting model was constructed to examine how PI would change with increasing EC_other.

To evaluate the relationship between energy demand in other sectors and global energy prices, a logistic regression model was applied. The resulting equation was:

With an R² of 0.7868, this model was capable of explaining about 79% of the increase in energy prices due to the growth of energy demand in other sectors. Assuming that there were no other variables playing roles in the system, a reasonable prediction of the future energy prices could be made.

Building upon the previous analysis regarding HPC development and rising energy demand in other sectors, these two factors were integrated into the energy structure model to improve its representational accuracy. A multivariate polynomial regression was used to fit the relationship between average energy consumption EC_avg2 and the two variables: peak computing power R_peak and energy price index PI. The specific degrees of two independent variables and interaction terms were determined based on the R-squared value and consistency with reality. The regression model was in the following form:

A similar model was constructed for estimating the maximum energy consumption EC_max2, using the same independent variables and polynomial form.

Subsequently, the regression functions were brought back into the carbon emission function from Section 3.2:

Where

By altering the percentage of various energy sources, this approach could provide insights into the change in carbon emissions.

Indeed, building upon these inputs, the following polynomial regression equations were established for average and maximum energy consumption:

With R² = 0.8425. In this model, the variable PI represented the normalized energy price index, where the original price values were standardized using a mean of 75.09 and a standard deviation of 18.6. The variable R_peak represented the normalized peak computing power of HPC systems, standardized using a mean of 3.528 × 10⁹ and a standard deviation of 3.341 × 10⁹. Moreover, for maximum energy consumption:

With R² = 0.9983. In this model, the variable PI and R_peak were normalized based on the same standard above.

Here, Figures 6A–D provided visual representations of the polynomial fitting surfaces for EC_avg and EC_max, respectively. In Figures 6A,B, when energy prices were extremely low, average energy consumption decreased as R_peak increases, reflecting the energy-saving benefits of improved computational efficiency. When computing power was low, a rise in energy prices also led to reduced consumption, likely due to suppressed computational demand. However, at higher levels of R_peak, the opposite trend emerged—higher prices coincided with increased energy use. This may be explained by the behavior of large-scale HPC centers, where organizations concentrated intensive workloads under rising energy costs to optimize budget efficiency, resulting in greater total consumption. Figures 6C,D revealed similar patterns for EC_max, though the trends appeared more pronounced. In particular, maximum energy consumption increased rapidly when both computing power and energy price were high. Overall, despite potential uncertainties, both models demonstrated strong predictive performance, with high goodness-of-fit and consistency with observed real-world behavior.

To make a prediction of the future EC_avg and EC_max2, the growth models of R_peak and PI ought to be constructed. The value of R_peak in 2030 could be readily obtained, while for PI, future data of EC_other was required. A linear regression model was constructed to make this projection, which was given by:

By substituting the value of EC_other into the previous function of PI, the change in PI in the following years was calculated. Then, the scope of EC in 2030 can be acquired.

To further estimate the scope of carbon emissions in the 2030s, projections of national energy structures were required. In this part of the study, the Grey Prediction Model GM(1,1) was primarily employed.

To construct the input for the GM model, data from the years 2013, 2016, 2019, and 2022 were collected for each country because of availability, detailing the share of four major energy types—nuclear, coal, oil, and natural gas—in total national energy consumption (20–23). Thus, for each energy type and country, a separate GM(1,1) model was built to project its share in 2031, representing the carbon emission scenario in the 2030s. All grey model forecasts followed the standard GM(1,1) procedure: the original data series x⁽⁰⁾ was first accumulated to generate the series x⁽¹⁾; the background values z⁽¹⁾ were then computed as the average of consecutive elements in x⁽¹⁾; finally, the parameters a and b were estimated via least squares using the equation x ( 0 ) ( k ) + a z ( 1 ) ( k ) = b . The initial value x ( 0 ) ( 1 ) , corresponding to the first observation year (2013), was used as the initial condition to solve the model analytically:

Although the original years were 2013, 2016, 2019, 2022, for modeling convenience, the time series was set as k = 1, 2, 3, 4, assuming uniform step intervals.

The validity and accuracy of these models were assessed through level ratio tests and relative error evaluations, which provided guidance on both the applicability of the GM model and the reliability of its forecasts.

In cases where the GM model failed the diagnostic tests or yielded low-accuracy results, alternative curve-fitting methods were employed. These included linear regression, logarithmic regression, inverse regression, and S-curve (logistic) fitting. For each problematic data series, multiple models were tested and the one with the highest R² value and best visual fit to empirical patterns was selected.

Data limitations were also addressed appropriately. For China, India, Japan, and South Korea, where nuclear energy data for 2019 were unavailable, the missing value was replaced by the average of the 2016 and 2022 figures. Negative values produced by any forecasting method were adjusted to zero, and the share of renewable energy was subsequently calculated as the residual from unity (i.e., 1—sum of other shares). If this residual yielded a negative value, it was likewise adjusted to zero.

Subsequently, using the energy consumption predictions and the carbon emission function introduced in Problem 2, the final estimates of national and global carbon emissions in 2030 were calculated. To obtain the energy consumption in 2031 in each country, GM(1,1) was again utilized, with data collected from the HPC TOP500 List in 2013, 2016, 2019, and 2022. After classification and summation, the researcher calculated the power ratio R_P by dividing the total average power P_{avg, i} of each country by the total average power of the world P_avg. All the data was in original percentage shares without normalization. Due to some data deficiency, missing values were filled up using average values, and abnormal values were processed as well. Then, the projected power ratio in 2031 was calculated through GM(1,1), with level ratio tests and relative error evaluations determining whether the model was suitable. Again, for cases where the GM model did not fit, linear regression, logarithm regression, or inverse regression were incorporated as alternatives.

Since

Where E_{avg, i} refers to the average energy consumption of each country, and E_avg2 refers to the average energy consumption in the world, which could be obtained through the equation in Section 3.3.3. Therefore, the energy consumption EC_i in each country was given by:

Through Equation 13, the final carbon emission value could be obtained.

In order to estimate the energy price index P in 2030, the following linear regression was used to forecast EC_other:

With an R² value of 0.9092, indicating a strong linear growth trend. Substituting the projected value of EC_other into Equation 12, future values of PI were obtained. These, together with the extrapolated value of R_peak, enabled the calculation of both EC_avg and EC_max for 2030, which were 1.127180873 and 1.623270845, respectively. In 2031, EC_avg and EC_max were 1.35256363 and 1.967819606.

To further quantify the resulting carbon emissions, national energy structures were forecasted for each country. In evaluating the accuracy of the national energy mix projections, the level ratio test and relative error evaluation demonstrated that the majority of the GM(1,1) models achieved acceptable predictive accuracy and fit. In cases where the GM model showed instability or failed the tests, alternative modeling approaches were adopted. Notably, Sweden’s nuclear energy share was forecasted using linear regression due to the instability of GM results, while France’s natural gas data required the removal of outliers and subsequent linear fitting. Germany, Norway, Poland, and Sweden also employed logarithmic regression for natural gas projections after outlier removal. For Switzerland, the natural gas share was forecasted using an S-curve (logistic) model, which captured the gradual growth pattern more effectively than alternatives. The average power share for each country was calculated using consistent methods; detailed procedures were not included here due to the extensive volume of data. Missing values were imputed by taking the arithmetic mean of adjacent years, while abnormal entries—caused by limitations in the original data source—were corrected using linear interpolation. This approach was justified by the fact that HPC systems were large-scale infrastructure that were rarely dismantled or abruptly shut down, and thus their energy consumption rates were expected to remain relatively stable over time. The results were shown in the table below. Through calculation, the total carbon emission generated by maximal power operation was 1.071 × 10²⁰ kg.

Detailed information is found in Supplementary Table 2.

Finally, a visual representation of the resulting emission bounds under different structural assumptions was provided in Figure 7, which illustrated the potential variability in national and global carbon emissions based on the projected changes in energy consumption and structure.

To quantify the impact of input uncertainties on the model output, a global sensitivity analysis was conducted using the Sobol method, a variance-based technique capable of decomposing the output variance into contributions attributable to individual input variables and their interactions. The objective was to identify the relative influence of the energy price index and peak computing power on the predicted average energy consumption per unit of computing power in 2030.

The Sobol method was grounded in the functional ANOVA decomposition of a model function f ( X ) , where f ( x ) = ( X 1 , X 2 ) represents the input parameter vector. The total variance of the model output Y = f ( X ) was decomposed as:

Where V i = VarX i ( E [ Y ∣ X i ] ) denoted the variance contribution of input X_i, V_ij represented the second-order interaction between X_i and X_j, and the final term accounts for higher-order interactions. From this decomposition, three key sensitivity indices were defined: the first-order index S i = V i Var ( Y ) the total-order index S Ti = 1 − V ~ i Var ( Y ) (where V ~ i denoted the variance excluding X_i), and the second-order index S ij = V ij Var ( Y ) , which quantified pairwise interactions.

In this study, the model inputs included the original energy price index P I original and the original peak computing power R_{peak original}), both defined within the range (10⁻⁶, 10¹⁰). To ensure adequate coverage of the input space and accurate estimation of sensitivity indices, input samples were generated using the Saltelli sampling scheme, which was a quasi-Monte Carlo method designed specifically for Sobol analysis. A base sample size of 1,024 was used, resulting in a total of 6,144 evaluations, consistent with the sampling requirement of N(k + 2), where N was the base sample size and k the number of input variables.

Each sampled input pair was evaluated through a standardized model function. The original variables were first standardized to improve numerical stability and interpretability. The standardized variables were then substituted into the following predictive equation:

Where PI = P I original – 75.09 18.6 and R peak = R peak original − 3.528 × 10 9 3.341 × 10 9

The Sobol sensitivity analysis was then performed using the SALib Python library. The output variance was decomposed to compute the first-order sensitivity indices (S₁), representing the direct contribution of each variable to output variance, and the total-order indices (S_T), accounting for both direct and interaction effects. Additionally, second-order indices (S₂) were calculated to evaluate the joint interaction effect between the two inputs.

This approach provided a comprehensive assessment of how variability in the input space propagates through the model and affects predictions of future energy efficiency. The results were visualized through bar plots for S₁ and S_T, and a heatmap for S₂, enabling a detailed interpretation of individual and interactive parameter influences.

For the first-order Sobol sensitivity indices, the parameter R_{peak original} showed a high first-order index of 0.7214, indicating it contributed the most to the output variance when considered individually. In contrast, PI_original had a much lower first-order index of 0.0459, suggesting a relatively minor direct effect. The corresponding confidence intervals were ±0.0878 for R_{peak original} and ±0.0459 for PI_original.

Regarding the total-order Sobol indices, R_{peak original} again exhibited a dominant influence with a total-order index of 0.9578, compared to 0.2815 for PI_original. This indicates that R_{peak original} not only has a strong direct effect but also plays a substantial role through interactions with other parameters. The confidence intervals were ±0.0862 and ±0.0372 for R_{peak original} and PI_original, respectively.

As for the second-order interaction effects, the interaction between PI_original and R_{peak original} yielded a Sobol index of 0.2363, with a confidence interval of ±0.0810. This suggested that there was a notable synergistic effect between the two parameters that contributes appreciably to the overall model uncertainty.

The results indicated that R_{peak original} was the dominant influencing factor, contributing the mostly to both the first order and total variance. A strong interaction effect was also observed between the two input variables.

The followings were the basic model of Fuzzy Decision, the detailed model was in Supplementary Table 3.

Subsequently, the AHP was utilized to determine the weight of each evaluation criterion. The consistency index was calculated to ensure that the results met the required threshold (see Table 2).

The table presented the AHP results, which systematically derived the relative weights of evaluation criteria through pairwise comparisons. The weight distribution demonstrated that risk assessment carried the highest priority (23.65%), followed by environmental impact (22.355%) and sustainability (18.65%), while feasibility (12.737%), economic investment (12.712%), and policy support (9.897%) constituted relatively lower weights in the decision framework.

The procedure began with the construction of a judgment matrix A, where each element a_ij represented the relative importance of factor a_i compared to a_j, thus forming a subjective evaluation matrix.

Then, when determining the weights of each indicator, this study adopted the square root method to calculate the eigenvector.

The specific steps were as follows: First, this paper computed the geometric mean of each row in the judgment matrix to derive the initial weights. Next, these values were normalized to obtain the relative weights (eigenvector), followed by a consistency check to ensure the reliability of the results.

This approach ensured an objective and systematic determination of the weights for each factor in the AHP model.

Then, normalized the weight vector.

The final weight vector:

Next, a consistency check was conducted. The first step involved calculating the maximum eigenvalue (λ_max) of the judgment matrix.

Among them, (AW) _i was the i-th element of the product of the matrix A and the weight vector W.

Then the consistency index CI and the consistency ratio CR were calculated.

Here, RI was the random consistency index, obtained by looking up the table. If CR < 0.1, the matrix passed the consistency test; otherwise, the judgment matrix needed to be adjusted.

The consistency validation confirmed the reliability of these judgments, with the maximum eigenvalue (λ_max) calculated as 6.406. Using the standard random consistency index (RI) of 1.25 for the matrix dimension, the consistency ratio (CR = CI/RI = 0.065) satisfied the acceptable threshold (CR < 0.1), thereby validating the consistency of the pairwise comparison matrix. This mathematical verification ensured the logical coherence of the expert judgments in the AHP implementation.

Finally, weights for the scores of t technologies and policy were evaluated to obtain the final total plan score. The formula for the final total scores, considering the political and technology policy score, was:

Since this study considered both policies to be equally important, we simply added their scores together and divided by two.

Finally, we calculated both Double-Weight Method Scores and Total Score Method Scores for each policy and method. First, the Double-Weight Method Scores for the eight policies—including the use of low-power processors and GPUs, establishment of energy efficiency standards for HPC centers, installation of solar panels or wind turbines, encouragement of HPC centers to invest in renewable energy infrastructure, design of waste heat recovery systems, regulated recovery and utilization of waste heat, and establishment of an environmental monitoring system to track energy consumption and emissions in real time—were 0.172, 0.172, 0.162, 0.164, 0.161, 0.173, 0.162, and 0.159, respectively. The corresponding Total Score Method Scores were 58.553, 58.512, 54.973, 55.829, 54.659, 58.819, 55.122, and 54.079.

The final Double-Weight Method Scores for the methods—including Energy Efficiency Optimization, Renewable Energy Integration, Waste Heat Recovery System, and Environmental Monitoring and Reporting—after incorporating both political and technological policies were 0.1722, 0.1630, 0.1670, and 0.1606, respectively. The final Total Score Method Scores for these methods were 58.5329, 55.4015, 56.7392, and 54.6009.

In conclusion, based on both the data and the graph, the Energy Efficiency Optimization scheme achieved higher scores in both evaluation systems. Therefore, Energy Efficiency Optimization was highly recommended.

To incorporate energy efficiency schemes into our model and optimize its performance, this study first collected energy efficiency data from 2016 to 2023 based on the Top500 HPC rankings (6, 7). It was assumed that the energy efficiency (EE) of HPC systems improved over time at a natural rate, influenced by the current state of efficiency and technological constraints. To evaluate the effectiveness of incentive mechanisms in improving energy efficiency and reducing energy consumption, this study first constructed a dynamic model based on differential equations to describe the temporal evolution of energy efficiency (EE). A logistic growth model was employed to represent the natural progression of technology:

Here, EE_old denoted the energy efficiency (i.e., computing power per unit of energy) at a given time t in the baseline model. The parameter b represented the intrinsic growth rate, reflecting the pace of technological advancement, while a denoted the theoretical upper bound of energy efficiency. The initial condition was set as EE(0) = EE₀. This model captured the typical logistic growth behavior: rapid initial improvement, gradually slowing progress due to increasing technological difficulty, and eventual saturation as the efficiency approaches the upper limit.

Building on the policy incentive mechanism for high-efficiency HPC systems discussed in the previous section, the model was extended as follows:

In this formulation, EE_new represents the energy efficiency under policy intervention, and δ > 0 captures the additional growth rate attributed to policy incentives. This modification allowed the efficiency to improve more rapidly and approach the upper limit an earlier.

Analytical solutions to both differential equations were derived, and the parameters were subsequently determined via least squares fitting. Based on this, a logarithmic regression model was established between EC_max and EE, in order to predict how changes in energy efficiency influence the maximum energy consumption.

Here, a denoted the maximum achievable energy efficiency, b represented the natural growth rate, and δ reflected the additional growth rate induced by policy incentives.

By solving the above differential equations analytically, this study derived the following logistic functions for EE under the two scenarios:

In these equations, c indicated the inflection point, corresponding to the time at which EE reached half of its maximum value and experienced the most rapid growth.

Based on historical data fitting, the study determined the model parameters as follows: the theoretical upper bound of energy efficiency b was 56.7134, the intrinsic growth rate b was 0.3675, and the policy-induced additional growth rate δ was 0.2325. Furthermore, the predicted saturation years for energy efficiency under the baseline and policy-intervention scenarios were c_old = 2026.4 and c_new = 2024.7, respectively. These results indicated that policy incentived accelerate the improvement of energy efficiency, leading to an earlier convergence toward the maximum achievable efficiency.

The baseline model exhibited a high goodness-of-fit with R² = 0.96286, while the incentivized model achieved a slightly lower but acceptable R² = 0.84245, still capturing the accelerated growth trend effectively. The value of δ = 0.2325 was selected through multiple simulations, balancing empirical fit with plausible policy impact.

Building upon the EE projection, this study further established a regression model to investigate the relationship between EE and the maximum potential energy consumption per unit time ( E C max ). The following logarithmic regression function was obtained:

This model yielded a goodness-of-fit of R² = 0.7715 and presented a concave-up, decreasing trend, indicating that improvements in energy efficiency were associated with reductions in energy consumption.

To quantify the impact of the incentive mechanism, this study forecasted EE and E C max values for the year 2025 under both baseline and incentivized scenarios. The results are presented in Table 3.

As shown in Table 3, under the incentivized scenario in 2025, energy efficiency increased significantly (from 21.22 to 30.90), while the maximum energy consumption decreased (from 0.3449 to 0.3278). These findings suggested that enhancing energy efficiency through policy incentives played a critical role in mitigating energy consumption.