China’s Inter-Provincial Energy Security Resilience Assessment over Space and Time: An Improved Gray Relational Projection Model

Abstract

In recent years, with the increasing impact of extreme weather events on energy security, energy vulnerability has increased significantly, and more and more international institutions and departments have begun to incorporate resilience governance into energy security. This paper focuses on China’s inter-provincial energy security assessment. Compared with existing relevant research, the significant features of our work are (i) introducing the concept of energy resilience and presenting its evolution mechanism and evaluation criteria, (ii) developing a gray relational projection model by using the level difference maximization and optimization theory, (iii) measuring the energy resilience of 30 Chinese provinces over space and time. Our results show that the spatial–temporal patterns of energy resilience in China changed significantly from 2005 to 2018. High energy resilience moved from provinces with abundant nonrenewable energy before 2010 to provinces with high energy diversity. Energy endowment is a primary condition to ensure a region’s energy resilience. Renewable energy development, energy investment, economic development, and policy coordination play vital roles in ensuring regional energy resilience. Energy investment and economic development can effectively improve the energy resilience of resource-poor areas. This study’s results will serve as a reference for China and contribute to expanding knowledge in this field.

1. Introduction

Energy resilience is one of the major standards to measure energy security. Facing the increasingly complex background of COVID-19 and globalization, especially with the conflict between Russia and Ukraine and the natural gas crisis in Europe, energy security is attracting more and more attention from the international community. With the rapid development of China’s economy, energy resilience has become particularly important. As the major emerging economy and the largest energy-consuming country globally [1], the gap between China’s energy supply and demand is increasing [2]. In 2020, China’s dependence on foreign oil and natural gas reached 72% and 43%, respectively [3]. It is estimated that by 2030, China’s energy imports will surpass Europe’s and become the world’s largest energy importer [4], thus placing great pressure on China’s energy security. More importantly, China is also accelerating the process of carbon peaking and carbon neutralization. The threats facing the energy system in the future will be more difficult to predict, such as technology and power system flexibility and so on. The Chinese government attaches great importance to energy security; energy supply risk aversion has been included in the national energy security strategy since 1996. For the first time, energy security has been included in China’s 14th five-year plan for national economic and social development. However, the traditional concept of energy security emphasizes the security of supply and the availability of abundant fossil fuels, such as oil, natural gas, and coal [5]. Therefore, it is of great value to explore China’s energy resilience to mitigate the new risks and challenges of energy security in the process of achieving the dual-carbon goal.

The word “resilience” has been in the English language for hundreds of years; its original meaning is the ability to recover from adversity [6]. Since the 1980s, ongoing global changes have generated considerable interest in research on resilience [7]. “Resilience thinking” is a new scientific and political paradigm [8] widely used in disaster reduction, engineering, ecology, and other fields [9]. Energy resilience research started relatively late and was accompanied by research on energy security. In recent years, many scholars have paid attention to energy resilience, for example, Yergin pointed out that resilience is a “margin of security” in energy supply systems, which provides a buffer to withstand shocks and promote recovery after interruptions [10]. Further, Thomas and Kerner characterized the need for energy resilience metrics [11]. McLellan et al. incorporated the concept of energy elasticity into sustainability; the authors foresaw a triangular relationship between resilience, sustainability, and risk management in response to disasters [12]. The resilience of a system is not solely dependent on physical disruptions but also on dynamic factors, such as societal and geo-political influences. Roege adopted a matrix format and proposed a framework for measuring energy resilience from four aspects: physical, information, cognitive, and social [13]. Additionally, Sandia National Laboratories proposed a seven-step resilience framework for local and national energy infrastructure [14]. Ding et al. pointed out that due to the low resilience of the natural gas importation network to import disruptions, China should increase the natural gas storage for sudden demand shortages [15]. Bento et al. discussed the concept of resilience in the oil and gas industry [16]. Abdin et al. reported that considering resilience during the planning process can significantly increase energy systems’ resilience [17]. Durán-Romero et al. proposed a framework for action to governments, businesses, and society based on the contribution of the Circular Economy (CE) towards sustainability [18]. Some scholars believe that energy diversity is the key to improving energy resilience [19]; when there are redundant and diverse parameters, the energy system’s resilience will be enhanced, and the supply risk will be reduced [20].

Most of the literature mentioned above are largely qualitative studies, which focused on the energy resilience of a specific region or community, or discussed the energy resilience of different social sectors. It is worth mentioning that, as early as 2009, Kruyt et al. pointed out that assessments are fundamental to providing adequate resilient actions [21]. Jansen also noted that integrating methods to assess society’s resilience could meet their need for energy services in a longer time frame [22]. Multi-Criteria Decision Analysis (MCDA) is most commonly used for addressing multiple conflicting objectives [23]. However, quantitative research on energy resilience is still lacking [24,25], a unified evaluation standard for resilience has not been formed, and scholars have not reached a consensus on the definition of energy resilience. Moreover, research on the mechanism of energy resilience has not yet been reported, quantitative research methods for energy system resilience are lacking, and scholars have not studied China’s energy resilience at the provincial level. Based on the scholars’ research mentioned above, this article proposed a definition of energy resilience as a vital attribute of energy security. It is the shock absorber of a regional energy system that can alleviate the system’s complexity resulting from political, economic, technological, and environmental issues. It can also predict the negative consequences of unexpected interference of factors to ensure the operation of the energy system or minimize negative consequences. Hence, this paper attempts to explore the spatial and temporal evolution characteristics of energy resilience in 30 provinces of China from 2005 to 2018. Due to the limitation of data availability, four provinces of Hong Kong, Macau, Taiwan, and Tibet are excluded in this paper. The remarkable contributions of this paper can be clearly illustrated as follows: First, this paper proposes the energy resilience evolution mechanism and constructs an energy resilience index for 30 Chinese provinces. Second, this paper constructs an improved gray relational projection evaluation model. This model overcomes the limitations of the single evaluation and gray correlation methods. Third, this paper furnishes key information on energy resilience in 30 provinces in China.

The remainder of this paper is organized as follows: Section 2 provides the evolution mechanism of energy resilience, the selection process of the energy resilience index, and the evaluation method in detail. Section 3 lays out the ranking, evolution trend, and spatial evolution law of energy resilience of each province. Section 4 discusses the results of this paper considering relevant research. Section 5 concludes the paper.

2. Methodology

2.1. The Evolution Mechanism of Energy Resilience

The basis of, and key to, investigating energy resilience is the evolution mechanism of energy resilience. A region should consider the following capabilities during the energy-planning process to improve energy resilience: first, improving energy efficiency to conserve resources and minimize energy supply costs; second, increasing the diversity of the energy supply to reduce energy supply risks; third, maintaining spare capacity or redundancy to manage unexpected surges in demand or interference from sources of uncertainty. These three methods increase the energy system’s ability to absorb interference, improve the energy system’s flexibility, and reduce the risk of sudden disturbances to the energy supply. How are these capabilities measured? What should the core required capabilities of energy resilience building be for a region? There is not yet a consensus on these issues [26]. The concept of energy resilience has not been formed in China and most other countries. In order to comprehensively evaluate China’s energy resilience, this article proposes a four-stage energy resilience evolution mechanism: “preparation, absorption, mitigation, and adaptation” (see Figure 1). The closest reference to this article is the work of Chen et at [27]. The basic ideas and methods of our evaluation mechanism are as follows:

• Preparation $0<\mathrm{t}<{\mathrm{t}}_1$. At this stage, the energy system is in a state of resilience. Set energy resilience in this period as S. Predict and prepare for disruption using various planning and design measures (by identifying and improving critical thresholds) designed to avoid and withstand potential disruptions and keep energy services available and assets functioning during emergency disturbances to achieve reliable energy services. Specifically, this includes strengthening energy resource reserves, maintaining energy infrastructures, and ensuring energy investments. Sufficient investment is a priority for resilient systems [28]. Additionally, it also includes minimizing the impact of the energy system on the environment [29].
• Absorption ${\mathrm{t}}_1\le \mathrm{t}<{\mathrm{t}}_2$. The main purpose of this stage is to examine the affordability of the energy system. Regardless of how well the energy system is prepared to withstand disturbances, potential shocks may exceed the system’s resistance threshold. When the system’s disaster-bearing capacity is insufficient to absorb the impact from ${\mathrm{t}}_1$, a disturbance occurs, and the energy system’s performance begins to decline. The energy system’s self-sufficiency rate, energy investment, the utilization efficiency of energy facilities, and other indicators are related to the regional energy system’s affordability. Efficient use of existing facilities can also reduce demand for new facilities and improve the energy system’s economic and environmental resilience [30].
• Mitigation ${\mathrm{t}}_2\le \mathrm{t}<{\mathrm{t}}_3$. This stage means that the energy system recovers to the preparation stage after the disturbance. In this stage, we should establish a risk-management method to rapidly restore the availability of all system operations and services to achieve efficiency in advance. Ideally, planning for the recovery process should begin before a disruptive event occurs. If the planning and absorption activities are appropriately implemented, the recovery process can be accelerated. Human intervention measures to reduce greenhouse gas and pollutant emissions decrease energy use, improve the diversity of energy production, help shorten the mitigation time after disturbances, and improve the energy system [31].
• Adaptation $\mathrm{T}\ge {\mathrm{t}}_3$. In this stage, adaptation is related to the degree of disturbance, policy support from the government, and recovery time. After an energy crisis, the state and affected departments should introduce more policies, build more effective crisis-response mechanisms, and support the resilient development of the energy system by learning and absorbing the experience and lessons from a disaster. Specifically, it includes increasing local energy reserves and improving the level of energy diversification [32]. The performance of the energy system is improved from the lessons learned from interference to achieve higher energy resilience (S⁺).

2.2. Energy Resilience Evaluation Index Construction

The evaluation index system of energy resilience should include four stages of the evolution mechanism of energy resilience. The four-stage energy resilience evolution mechanisms are mutual influence and promotion, a closed-loop process [33,34]. Therefore, in the design process of an energy resilience and evaluation index system, we should fully consider the four stages in Figure 1. Energy resilience is also an interdisciplinary concept [35,36]. In the current context of rapidly increasing complexities and deep uncertainty, the concepts of resilience and sustainability are poorly or too narrowly defined in the urban context and can be used interchangeably [37]. To the best of our knowledge, so far, no studies have clearly clarified the boundaries of resilience and sustainability.

Moreover, many energy security indicators also apply to resilient environments [38,39,40]. The main difference between energy resilience and energy security indicators is that energy security indicators include all energy resilience indicators, but resilience indicators are only one type of energy security indicator. Resilience indicators focus on the energy system’s resistance to destructive events and indicators to improve the system’s ability to rebuild. Therefore, based on this analysis and guided by the resilience requirements in the international sustainable development goals, referring to the national academies [41,42], the interactions of energy resilience with society, economy, and governance should be considered comprehensively according to the actual conditions in each province in China and the data availability. China’s energy resilience index (CERI) is measured by five dimensions: availability, diversity, economic resilience, environmental resilience, and technological resilience. Figure 2 provides the visual relationship between the evolution mechanism of energy resilience and its performance indicators. The detailed calculation process of each indicator is shown in

Table 1. The calculation process of each indicator and sources for CERI.

Dimensions	Indicators	No.	Attribute	Equations	Variable Description	Indicator Source
Availability	Reserve and production ratio	I1	Positive	${\displaystyle {\displaystyle \sum }_{\mathrm{i}}}{\mathrm{p}}_{\mathrm{i}}{\mathrm{s}}_{\mathrm{i}}$	${\mathrm{p}}_{\mathrm{i}}$: the ratio of energy reserve to production in each province. ${\mathrm{s}}_{\mathrm{i}}$: the proportion of the$\mathrm{ith}$, energy in the province’s total energy production.	[43,44,45]
	Energy self-sufficiency	I2	Positive	$\raisebox{1ex}{${\mathrm{e}}_{\mathrm{pro}}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{e}}_{\mathrm{con}}$}\right.$	${\mathrm{e}}_{\mathrm{pro}}$: total energy production in each province. ${\mathrm{e}}_{\mathrm{con}}$: total energy consumption in each province.	[46,47,48]
Diversity	Production diversity index	I3	Positive	$\sqrt{{\displaystyle {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}}{\mathrm{s}}_{\mathrm{i}}^2}$	${\mathrm{s}}_{\mathrm{i}}$: the proportion of the$\mathrm{ith}$ energy in the province’s total energy production.	[46,49,50]
	Consumption diversity index	I4	Positive	$\mathrm{SWI}=-{\displaystyle {\displaystyle \sum }_{\mathrm{I}}}{\mathrm{C}}_{\mathrm{i}}\ln \left({\mathrm{C}}_{\mathrm{i}}\right)$	${\mathrm{c}}_{\mathrm{i}}$: the proportion of the$\mathrm{ith}$ energy in the province’s total energy consumption.	[51,52,53]
Economic resilience	Energy investment	I5	Positive	------	Each province’s investment in fixed assets of the energy industry.	[54,55,56]
Environmental resilience	GDP energy intensity	I6	Negative	$\raisebox{1ex}{${\mathrm{e}}_{\mathrm{con}}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{GDP}$}\right.$	${\mathrm{e}}_{\mathrm{con}}$: total energy consumption in each province. GDP: gross domestic product in each province.	[48,50,53,57]
Technical resilience	Utilization time of power generation equipment	I7	Positive	$\raisebox{1ex}{${\mathrm{t}}_{\mathrm{pow}}$}\!\left/ \!\raisebox{-1ex}{${\mathrm{h}}_{\mathrm{tot}}$}\right.$	${\mathrm{t}}_{\mathrm{pow}}$: annual operating hours of power plants in each province. ${\mathrm{h}}_{\mathrm{tot}}$: total annual hours in each province.	[47,48,49,58]
	Distribution loss of power system	I8	Negative	------	------	[48,53,58]

.

2.3. Energy Resilience Evaluation Model Construction

In this subsection, our main objective is to develop a dynamic evaluation model for the evolution of energy resilience. Due to grey relational analysis (GRA) being an evaluation model based on the data of the influencing factors of the research object, this method uses mathematical methods to study the geometric correspondence between the factors to determine the order of each research project [59,60]. One of the main advantages of GRA is that it uses a relatively small amount of data or elements with considerable variability to produce satisfactory results. It has been widely used in agriculture, industry, and energy and has achieved significant results [61,62,63,64]. However, this method only obtains the energy resilience ranking of each province but cannot accurately assess each province’s energy resilience strength. Traditional comprehensive evaluation methods based on a combination of indicators and weights that were calculated by principal component analysis (PCA) [65] and the entropy method [66,67] can obtain accurate evaluation values of energy resilience. However, when solving complex high-latitude problems, they are easily affected by outliers and lack robustness. In this paper, we thoroughly considered the advantages and disadvantages of GRA and comprehensively evaluated techniques [68], based on the fusion of five methods: PCA, gray correlation, entropy, mean square error, and projection methods [69]; thus, to a certain extent, this paper can be regarded as the development of the GRA model by improving a gray relational projection model that maximizes the level differences. The flowchart of the proposed improved gray relational model is shown in Figure 3. The key calculation steps of the improved gray relational projection model are as follows.

If the index value corresponding to the $\mathrm{j}\left(\mathrm{j}=1, 2,\dots \mathrm{n}\right)$ evaluation indicator of a $i(i=1, 2,\dots \mathrm{m})$ evaluation object at time $\mathrm{t}$ is ${\mathrm{x}}_{\mathrm{ij} }\left(\mathrm{t}\right)$, then according to the indicator data’s characteristics, as referenced by [70], in this paper, we use the following formula to normalize each indicator:

$${\mathrm{b}}_{\mathrm{ij}}^{}=\frac{{\mathrm{x}}_{\mathrm{ij} -}{ \overline{\mathrm{x}}}_{\mathrm{j}}}{{\mathrm{s}}_{\mathrm{j}}}$$

(1)

where ${ \overline{\mathrm{x}}}_{\mathrm{j}}{ \mathrm{and} \mathrm{s}}_{\mathrm{j}}$ (j = 1, 2, …, m) are the sample averages and mean square deviation of the observation value of the jth indicator, respectively.

Let the decision matrix after normalization be B(t), then

$$\mathrm{B}\left(\mathrm{t}\right)=\left[\begin{array}{cc}\begin{array}{ccc}{\mathrm{b}}_{11}\left(\mathrm{t}\right)& \dots & {\mathrm{b}}_{1\mathrm{j}}\left(\mathrm{t}\right)\\ \cdots & \cdots & \cdots \\ {\mathrm{b}}_{\mathrm{i}1}\left(\mathrm{t}\right)& \cdots & {\mathrm{b}}_{\mathrm{ij}}\left(\mathrm{t}\right)\end{array}& \begin{array}{c}\begin{array}{cc}\dots & {\mathrm{b}}_{1\mathrm{n}}\left(\mathrm{t}\right)\end{array}\\ \begin{array}{cc}\cdots & \cdots \end{array}\\ \begin{array}{cc}\cdots & {\mathrm{b}}_{\mathrm{in}}\left(\mathrm{t}\right)\end{array}\end{array}\\ \begin{array}{ccc}\cdots & \cdots & \cdots \\ {\mathrm{b}}_{\mathrm{m}1}\left(\mathrm{t}\right)& \cdots & {\mathrm{b}}_{\mathrm{mj}}\left(\mathrm{t}\right)\end{array}& \begin{array}{c}\begin{array}{cc}\cdots & \dots \end{array}\\ \begin{array}{cc}\cdots & {\mathrm{b}}_{\mathrm{mn}}\left(\mathrm{t}\right)\end{array}\end{array}\end{array}\right]$$

(2)

Among the eight indicators used in this paper, there are six positive and two negative indicators. At time t, we assume that the indicator values of the positive and negative ideal objects are ${\mathrm{b}}^{+} \left(\mathrm{t}\right){ \mathrm{and} \mathrm{b}}^{-} \left(\mathrm{t}\right)$, then

$${\mathrm{b}}^{+}\left(\mathrm{t}\right)=\left[{\mathrm{b}}_{01}^{+}\left(\mathrm{t}\right)\dots {\mathrm{b}}_{0\mathrm{j}}^{+}\left(\mathrm{t}\right)\dots {\mathrm{b}}_{0\mathrm{n}}^{+}\left(\mathrm{t}\right)\right]$$

(3)

$${\mathrm{b}}^{-}\left(\mathrm{t}\right)=\left[{\mathrm{b}}_{01}^{-}\left(\mathrm{t}\right)\dots {\mathrm{b}}_{0\mathrm{j}}^{-}\left(\mathrm{t}\right)\dots {\mathrm{b}}_{0\mathrm{n}}^{-}\left(\mathrm{t}\right)\right]$$

(4)

According to the gray relational analysis, the correlation coefficient ${\mathrm{h}}_{\mathrm{ij} }\left(\mathrm{t}\right)$ is

$${\mathrm{h}}_{\mathrm{ij}}\left(\mathrm{t}\right)=\frac{{}_{1\le \mathrm{i}\le \mathrm{m}1\le \mathrm{j}\le \mathrm{m}}{}^{\min \min } \left\{{\Delta }_{\mathrm{ij}}\left(\mathrm{t}\right)\right\}+\mathsf{\rho }{}_{1\le \mathrm{i}\le \mathrm{m}1\le \mathrm{j}\le \mathrm{m}}{}^{\max \max } \left\{{\Delta }_{\mathrm{ij}}\left(\mathrm{t}\right)\right\}}{{\Delta }_{\mathrm{ij}}\left(\mathrm{t}\right)+\mathsf{\rho }{}_{1\le \mathrm{i}\le \mathrm{m}1\le \mathrm{j}\le \mathrm{m}}{}^{\max \max } \left\{{\Delta }_{\mathrm{ij}}\left(\mathrm{t}\right)\right\}}$$

(5)

where ${\Delta }_{\mathrm{ij}}\mathrm{t}=\left|{\mathrm{b}}_{\mathrm{ij}}\left(\mathrm{t}\right)-{\mathrm{b}}_{\mathrm{ij}}^{*}\left(\mathrm{t}\right)\right| \mathrm{and} {}_{1\le \mathrm{i}\le \mathrm{m}1\le \mathrm{j}\le \mathrm{m}}{}^{\min \min } \left\{{\Delta }_{\mathrm{ij}}\left(\mathrm{t}\right)\right\}$ and ${}_{1\le \mathrm{i}\le \mathrm{m}1\le \mathrm{j}\le \mathrm{m}}{}^{\max \max } \left\{{\Delta }_{\mathrm{ij}}\left(\mathrm{t}\right)\right\}$ are the smallest secondary and largest secondary errors, respectively, and$\mathsf{\rho }\in \left[0, 1\right]$ is the resolution coefficient.

When ${\mathrm{b}}_{\mathrm{ij}}^{*}\left(\mathrm{t}\right)$ is ${\mathrm{b}}_{\mathrm{ij}}^{+}\left(\mathrm{t}\right)$, ${\mathrm{h}}_{\mathrm{ij}}^{}\left(\mathrm{t}\right)={\mathrm{h}}_{\mathrm{ij}}^{+}\left(\mathrm{t}\right)$, and, at time t, the positive ideal gray incidence coefficient matrix is ${\mathrm{H}}^{+}\left(\mathrm{t}\right)$, then

$${\mathrm{H}}^{+}\left(\mathrm{t}\right)=\left[\begin{array}{cc}\begin{array}{ccc}{\mathrm{h}}_{11}^{+}\left(\mathrm{t}\right)& \dots & {\mathrm{h}}_{1\mathrm{j}}^{+}\left(\mathrm{t}\right)\\ \dots & \dots & \dots \\ {\mathrm{h}}_{\mathrm{i}1}^{+}\left(\mathrm{t}\right)& \dots & {\mathrm{h}}_{\mathrm{ij}}^{+}\left(\mathrm{t}\right)\end{array}& \begin{array}{c}\begin{array}{cc}\dots & {\mathrm{h}}_{1\mathrm{n}}^{+}\left(\mathrm{t}\right)\end{array}\\ \begin{array}{cc}\cdots & \cdots \end{array}\\ \begin{array}{cc}\cdots & {\mathrm{h}}_{\mathrm{in}}^{+}\left(\mathrm{t}\right)\end{array}\end{array}\\ \begin{array}{ccc}\cdots & \cdots & \cdots \\ {\mathrm{h}}_{\mathrm{m}1}^{+}\left(\mathrm{t}\right)& \cdots & {\mathrm{h}}_{\mathrm{mj}}^{+}\left(\mathrm{t}\right)\end{array}& \begin{array}{c}\begin{array}{cc}\cdots & \dots \end{array}\\ \begin{array}{cc}\cdots & {\mathrm{h}}_{\mathrm{mn}}^{+}\left(\mathrm{t}\right)\end{array}\end{array}\end{array}\right]$$

(6)

When ${\mathrm{b}}_{\mathrm{ij}}^{*}\left(\mathrm{t}\right)$ is ${\mathrm{b}}_{\mathrm{ij}}^{-}\left(\mathrm{t}\right)$, ${\mathrm{h}}_{\mathrm{ij}}^{}\left(\mathrm{t}\right)={\mathrm{h}}_{\mathrm{ij}}^{-}\left(\mathrm{t}\right)$, and, at time t, the negative ideal gray incidence coefficient matrix is ${\mathrm{H}}^{-}\left(\mathrm{t}\right)$, then

$${\mathrm{H}}^{-}\left(\mathrm{t}\right)=\left[\begin{array}{cc}\begin{array}{ccc}{\mathrm{h}}_{11}^{-}\left(\mathrm{t}\right)& \dots & {\mathrm{h}}_{1\mathrm{j}}^{-}\left(\mathrm{t}\right)\\ \cdots & \cdots & \cdots \\ {\mathrm{h}}_{\mathrm{i}1}^{-}\left(\mathrm{t}\right)& \cdots & {\mathrm{h}}_{\mathrm{ij}}^{-}\left(\mathrm{t}\right)\end{array}& \begin{array}{c}\begin{array}{cc}\dots & {\mathrm{h}}_{1\mathrm{n}}^{-}\left(\mathrm{t}\right)\end{array}\\ \begin{array}{cc}\cdots & \cdots \end{array}\\ \begin{array}{cc}\cdots & {\mathrm{h}}_{\mathrm{in}}^{-}\left(\mathrm{t}\right)\end{array}\end{array}\\ \begin{array}{ccc}\cdots & \cdots & \cdots \\ {\mathrm{h}}_{\mathrm{m}1}^{-}\left(\mathrm{t}\right)& \cdots & {\mathrm{h}}_{\mathrm{mj}}^{-}\left(\mathrm{t}\right)\end{array}& \begin{array}{c}\begin{array}{cc}\cdots & \dots \end{array}\\ \begin{array}{cc}\cdots & {\mathrm{h}}_{\mathrm{mn}}^{-}\left(\mathrm{t}\right)\end{array}\end{array}\end{array}\right]$$

(7)

We use PCA, entropy, and mean square error methods to calculate the target weight of each indicator, then the weight matrix is $\mathrm{U}\left(\mathrm{t}\right)$ is given by

$$\mathrm{U}\left(\mathrm{t}\right)=\left[\begin{array}{cc}\begin{array}{cc}{\mathrm{u}}_{11}\left(\mathrm{t}\right)& \dots \\ \dots & \dots \\ {\mathrm{u}}_{\mathrm{g}1}\left(\mathrm{t}\right)& \dots \end{array}& \begin{array}{cc}{\mathrm{u}}_{1\mathrm{j}}\left(\mathrm{t}\right)& \begin{array}{cc}\dots & {\mathrm{u}}_{1\mathrm{n}}\left(\mathrm{t}\right)\end{array}\\ \dots & \begin{array}{cc}\dots & \dots \end{array}\\ {\mathrm{u}}_{\mathrm{gj}}\left(\mathrm{t}\right)& \begin{array}{cc}\dots & {\mathrm{u}}_{\mathrm{gn}}\left(\mathrm{t}\right)\end{array}\end{array}\\ \begin{array}{cc}\dots & \dots \\ {\mathrm{u}}_{\mathrm{k}1}\left(\mathrm{t}\right)& \dots \end{array}& \begin{array}{cc}\dots & \begin{array}{cc}\dots & \dots \end{array}\\ {\mathrm{u}}_{\mathrm{kj}}\left(\mathrm{t}\right)& \begin{array}{cc}\dots & {\mathrm{u}}_{\mathrm{kn}}\left(\mathrm{t}\right)\end{array}\end{array}\end{array}\right]$$

(8)

where ${\mathrm{u}}_{\mathrm{gj} }\left(\mathrm{t}\right)$ is the weight of the jth index under the gth weighting method at time t, g = 1, 2, …, 3, and j = 1, 2, …, n.

If the combined weight to be solved at time t is ${\mathrm{w}}_{\left(\mathrm{t}\right)}$, then the reasonable interval of the combined weight of the jth indicator is expressed as

$${\mathrm{w}}_{\mathrm{j}}\left(\mathrm{t}\right)\in \left[{\mathrm{u}}_{\mathrm{j}}^{-1}\left(\mathrm{t}\right),{\mathrm{u}}_{\mathrm{j}}^{+1}\left(\mathrm{t}\right)\right]$$

(9)

$$\mathrm{If}{ \mathrm{Z}}_{\left(\mathrm{t}\right)}=\mathrm{w}\left(\mathrm{t}\right)\mathrm{X}\left(\mathrm{t}\right)=\left[\mathrm{w}\left(\mathrm{t}\right){\mathrm{x}}_1\left(\mathrm{t}\right)\dots \mathrm{w}\left(\mathrm{t}\right){\mathrm{x}}_{\mathrm{i}}\left(\mathrm{t}\right)\dots \mathrm{w}\left(\mathrm{t}\right){\mathrm{x}}_{\mathrm{m}}\left(\mathrm{t}\right)\right]$$

(10)

Let ${\mathrm{x}}_0\left(\mathrm{t}\right)=\frac{1}{\mathrm{m}}\left[{\mathrm{x}}_1\left(\mathrm{t}\right)+{\mathrm{x}}_2\left(\mathrm{t}\right)+,\dots ,+{\mathrm{x}}_{\mathrm{m}}\left(\mathrm{t}\right)\right]$ and the mean value of $\mathrm{Z}\left(\mathrm{t}\right)$is $\overline{{\mathrm{Z}}_{\left(\mathrm{t}\right)}}$, then

$$\begin{array}{ll}{\mathrm{Z}}_{\left(\mathrm{t}\right)}& =\frac{1}{\mathrm{m}}\left[\mathrm{w}\left(\mathrm{t}\right){\mathrm{x}}_1\left(\mathrm{t}\right)+\mathrm{w}\left(\mathrm{t}\right){\mathrm{x}}_2\left(\mathrm{t}\right)+,\dots ,+\mathrm{w}\left(\mathrm{t}\right){\mathrm{x}}_{\mathrm{m}}\left(\mathrm{t}\right)\right]\\ & =\frac{1}{\mathrm{m}}\mathrm{w}\left(\mathrm{t}\right)\left[{\mathrm{x}}_1\left(\mathrm{t}\right)+{\mathrm{x}}_2\left(\mathrm{t}\right)+, \dots ,+{\mathrm{x}}_{\mathrm{m}}\left(\mathrm{t}\right)\right]\\ & =\mathrm{w}\left(\mathrm{t}\right){\mathrm{x}}_0\left(\mathrm{t}\right)\end{array}$$

(11)

Suppose ${\mathrm{x}}_{\mathrm{i}}^{*}\left(\mathrm{t}\right)={\mathrm{x}}_{\mathrm{i}}\left(\mathrm{t}\right)-{\mathrm{x}}_0\left(\mathrm{t}\right)$ and the variance of ${\mathrm{Z}}_{\left(\mathrm{t}\right)}$ is ${\left[\mathrm{S}\left(\mathrm{t}\right)\right]}^2$, then

$$\begin{array}{ll}{\left[\mathrm{S}\left(\mathrm{t}\right)\right]}^2& =\frac{1}{\mathrm{m}-1}{\displaystyle {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{m}}}{\left[\mathrm{w}\left(\mathrm{t}\right){\mathrm{x}}_{\mathrm{i}}\left(\mathrm{t}\right)-\mathrm{w}\left(\mathrm{t}\right){\mathrm{x}}_0\left(\mathrm{t}\right)\right]}^2\\ & =\frac{1}{\mathrm{m}-1}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{m}}{\left[\mathrm{w}\left(\mathrm{t}\right){\mathrm{x}}_{\mathrm{i}}^{*}\left(\mathrm{t}\right)\right]}^2\\ & =\frac{1}{\mathrm{m}-1}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{m}}\mathrm{w}\left(\mathrm{t}\right){\mathrm{x}}_{\mathrm{i}}^{*}\left(\mathrm{t}\right){\left[\mathrm{w}\left(\mathrm{t}\right){\mathrm{x}}_{\mathrm{i}}^{*}\left(\mathrm{t}\right)\right]}^{\mathrm{T}}\\ & =\frac{1}{\mathrm{m}-1}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{m}}\mathrm{w}\left(\mathrm{t}\right)\left\{{\mathrm{x}}_{\mathrm{i}}^{*}\left(\mathrm{t}\right){\mathrm{x}}_{\mathrm{i}}^{*}{\left(\mathrm{t}\right)}^{\mathrm{T}}\right\}{\left[\mathrm{w}\left(\mathrm{t}\right)\right]}^{\mathrm{T}}\end{array}$$

(12)

Using ${\left[\mathrm{S}\left(\mathrm{t}\right)\right]}^2$ max as the objective function, the reasonable weight interval determined by the sum of the combined weights is 1, the maximum level difference model is constructed, and the optimal combination weight of each indicator is obtained. Then, we have the following control system:

$$\max \frac{1}{\mathrm{m}-1}{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{m}}\mathrm{w}\left(\mathrm{t}\right)\left\{{\mathrm{x}}_{\mathrm{i}}^{*}\left(\mathrm{t}\right){\mathrm{x}}_{\mathrm{i}}^{*}{\left(\mathrm{t}\right)}^{\mathrm{T}}\right\}{\left[\mathrm{w}\left(\mathrm{t}\right)\right]}^{\mathrm{T}}$$

(13)

$$\mathrm{s}.\mathrm{t}.\{\begin{array}{c}{\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{w}}_{\mathrm{j}}\left(\mathrm{t}\right)=1\\ {\mathrm{u}}_{\mathrm{j}}^{-1}\left(\mathrm{t}\right)\le {\mathrm{w}}_{\mathrm{j}}\left(\mathrm{t}\right)\le {\mathrm{u}}_{\mathrm{j}}^{+1}\left(\mathrm{t}\right)\end{array}$$

Accordingly, at time t, one can obtain the following positive and negative ideal weighted gray correlation coefficient matrix ${\mathrm{F}}^{+}\left(\mathrm{t}\right){ \mathrm{and} \mathrm{F}}^{-}\left(\mathrm{t}\right)$:

$${\mathrm{F}}^{+}\left(\mathrm{t}\right)=\left[\begin{array}{cc}\begin{array}{ccc}{{\mathrm{h}}_{11}^{+}\mathrm{w}}_1\left(\mathrm{t}\right)& \dots & {{\mathrm{h}}_{1\mathrm{j}}^{+}\mathrm{w}}_{\mathrm{j}}\left(\mathrm{t}\right)\\ \cdots & \cdots & \cdots \\ {{\mathrm{h}}_{\mathrm{i}1}^{+}\mathrm{w}}_1\left(\mathrm{t}\right)& \cdots & {{\mathrm{h}}_{\mathrm{i}\mathrm{j}}^{+}\mathrm{w}}_{\mathrm{j}}\left(\mathrm{t}\right)\end{array}& \begin{array}{c}\begin{array}{cc}\dots & {{\mathrm{h}}_{1\mathrm{n}}^{+}\mathrm{w}}_{\mathrm{n}}\left(\mathrm{t}\right)\end{array}\\ \begin{array}{cc}\cdots & \cdots \end{array}\\ \begin{array}{cc}\cdots & {{\mathrm{h}}_{\mathrm{i}\mathrm{n}}^{+}\mathrm{w}}_{\mathrm{n}}\left(\mathrm{t}\right)\end{array}\end{array}\\ \begin{array}{ccc}\cdots & \cdots & \cdots \\ {{\mathrm{h}}_{\mathrm{m}1}^{+}\mathrm{w}}_1\left(\mathrm{t}\right)& \cdots & {{\mathrm{h}}_{\mathrm{m}\mathrm{j}}^{+}\mathrm{w}}_{\mathrm{j}}\left(\mathrm{t}\right)\end{array}& \begin{array}{c}\begin{array}{cc}\cdots & \dots \end{array}\\ \begin{array}{cc}\cdots & {{\mathrm{h}}_{\mathrm{m}\mathrm{n}}^{+}\mathrm{w}}_{\mathrm{n}}\left(\mathrm{t}\right)\end{array}\end{array}\end{array}\right]$$

(14)

$${\mathrm{F}}^{-}\left(\mathrm{t}\right)=\left[\begin{array}{cc}\begin{array}{ccc}{{\mathrm{h}}_{11}^{-}\mathrm{w}}_1\left(\mathrm{t}\right)& \dots & {{\mathrm{h}}_{1\mathrm{j}}^{-}\mathrm{w}}_{\mathrm{j}}\left(\mathrm{t}\right)\\ \cdots & \cdots & \cdots \\ {{\mathrm{h}}_{\mathrm{i}1}^{-}\mathrm{w}}_1\left(\mathrm{t}\right)& \cdots & {{\mathrm{h}}_{\mathrm{i}\mathrm{j}}^{-}\mathrm{w}}_{\mathrm{j}}\left(\mathrm{t}\right)\end{array}& \begin{array}{c}\begin{array}{cc}\dots & {{\mathrm{h}}_{1\mathrm{n}}^{-}\mathrm{w}}_{\mathrm{n}}\left(\mathrm{t}\right)\end{array}\\ \begin{array}{cc}\cdots & \cdots \end{array}\\ \begin{array}{cc}\cdots & {{\mathrm{h}}_{\mathrm{i}\mathrm{n}}^{-}\mathrm{w}}_{\mathrm{n}}\left(\mathrm{t}\right)\end{array}\end{array}\\ \begin{array}{ccc}\cdots & \cdots & \cdots \\ {{\mathrm{h}}_{\mathrm{m}1}^{-}\mathrm{w}}_1\left(\mathrm{t}\right)& \cdots & {{\mathrm{h}}_{\mathrm{m}\mathrm{j}}^{-}\mathrm{w}}_{\mathrm{j}}\left(\mathrm{t}\right)\end{array}& \begin{array}{c}\begin{array}{cc}\cdots & \dots \end{array}\\ \begin{array}{cc}\cdots & {{\mathrm{h}}_{\mathrm{m}\mathrm{n}}^{-}\mathrm{w}}_{\mathrm{n}}\left(\mathrm{t}\right)\end{array}\end{array}\end{array}\right]$$

(15)

Each province is regarded as a row vector, and the energy resilience of the ith province in year t is expressed as

$${\mathrm{F}}_{\mathrm{i}}\left(\mathrm{t}\right)=\left[{\mathrm{h}}_{\mathrm{i}1}\left(\mathrm{t}\right){\mathrm{w}}_1\left(\mathrm{t}\right)\dots {\mathrm{h}}_{\mathrm{ij}}\left(\mathrm{t}\right){\mathrm{w}}_{\mathrm{j}}\left(\mathrm{t}\right)\dots {\mathrm{h}}_{\mathrm{in}}\left(\mathrm{t}\right){\mathrm{w}}_{\mathrm{n}}\left(\mathrm{t}\right)\right]$$

(16)

Using the least squares criterion, we establish the following objective function:

$$\mathrm{minF}{({\mathrm{y}}_{\mathrm{i}}\left(\mathrm{t}\right)={\left[{\mathrm{y}}_{\mathrm{i}}\left(\mathrm{t}\right)\right]}^2.{\mathrm{D}}_{\mathrm{i}}^{+}\left(\mathrm{t}\right)-{\mathrm{D}}_{\mathrm{i}}^{+}\left(\mathrm{t}\right)]}^2+{\left[1-{\mathrm{y}}_{\mathrm{i}}\left(\mathrm{t}\right).{\mathrm{D}}_{\mathrm{i}}^{-}\left(\mathrm{t}\right)-{\mathrm{D}}_{\mathrm{i}}^{-}\left(\mathrm{t}\right)\right]}^2$$

(17)

Differentiating the above expression $\mathrm{F}({\mathrm{y}}_{\mathrm{i}}\left(\mathrm{t}\right)$) with respect to ${\mathrm{y}}_{\mathrm{i}}\left(\mathrm{t}\right)$, one can obtain

$${\mathrm{y}}_{\mathrm{i}}\left(\mathrm{t}\right)=\frac{{\left[{\mathrm{D}}_{\mathrm{i}}^{+}\left(\mathrm{t}\right)\right]}^2}{{\left[{\mathrm{D}}_{\mathrm{i}}^{+}\left(\mathrm{t}\right)\right]}^2+{\left[{\mathrm{D}}_{\mathrm{i}}^{-}\left(\mathrm{t}\right)\right]}^2}$$

(18)

where $0<{\mathrm{y}}_{\mathrm{i}}\left(\mathrm{t}\right)<1$. The closer ${\mathrm{y}}_{\mathrm{i}}\left(\mathrm{t}\right)$ is to 1, the higher the energy resilience.

Finally, the dynamic change rate of energy resilience is used to measure the relative change degree of energy resilience in various provinces in different years. Let the dynamic change rate (%) of energy resilience of each province be${ \mathrm{d}}_{\mathrm{i}}\left(\mathrm{t}\right)$,

$${\mathrm{d}}_{\mathrm{i}}\left(\mathrm{t}\right)=\frac{{\mathrm{y}}_{\mathrm{i}+1}\left(\mathrm{t}\right)-{\mathrm{y}}_{\mathrm{i}}\left(\mathrm{t}\right)}{{\mathrm{y}}_{\mathrm{i}}\left(\mathrm{t}\right)}$$

(19)

2.4. Evaluation Criterion and Data

The CERI represents the ability of a region’s energy supply system to prevent interruption in the event of uncertainty or sudden interference, and its value range is [0, 1]. Combining with the evolution mechanism of energy resilience, referring to Zhao et al. [71,72] and following the step length of 0.2, the corresponding level of energy resilience is divided into five grades, as shown in

Table 2. Energy resilience level table.

Number	Security Grade	Score Range	Basic Characteristics
1	Ⅰ	0.8–1	When it is disturbed by uncertainty, the energy supply in the area is in a safe state
2	Ⅱ	0.6–0.8	When it is disturbed by uncertainty, the energy supply in the area is basically safe
3	Ⅲ	0.4–0.6	When disturbed by uncertainty, individual energy sources with high external dependence may be slightly short of supply during a specific period
4	Ⅳ	0.2–0.4	When disturbed by uncertainty, individual energy sources with high external dependence may be in short supply during a specific period
5	Ⅴ	0–0.2	When disturbed by uncertainty, the energy system in the region is very tight

Data: Associated data from 2005 to 2018 were collected from the 30 provinces examined in this study: Shaanxi, Ningxia, Xinjiang, Qinghai, Gansu, Inner Mongolia, Shanxi, Beijing, Tianjin, Hebei, Jilin, Liaoning, Heilongjiang, Henan, Hunan, Hubei, Guangdong, Guangxi, Hainan, Sichuan, Yunnan, Chongqing, Guizhou, Shandong, Fujian, Jiangsu, Jiangxi, Shanghai, Zhejiang, and Anhui. There are 31 provinces in China, but data for Xizang province were not available. The relevant operation data of indicators I1, I2, I3, I4, and I5 were obtained from the respective Statistical Yearbooks from the 30 provinces (2000–2020) and the China National Bureau of Statistics Database [3]. The relevant operation data of indicator I6 were obtained from the China Statistics Yearbook (2000–2020), and the power generation equipment utilization times and distribution losses data were acquired from the Wind Database.

.

3. Results

3.1. Optimal Weight of Each Indicator

From Section 2.3, we can obtain the optimal combination weight of each indicator, as shown in Figure 4 (${\mathrm{u}}_1$, ${ \mathrm{u}}_2$, and ${ \mathrm{u}}_3$ represent the indicator weights obtained using the mean square error, PCA, and the entropy weight methods, respectively, and ${ \mathrm{u}}_{\mathrm{t}}$ represents the optimal combination weight of each indicator). Simultaneously, we also use the vertical line chart to clearly show the dynamic change of the optimal weight of each indicator in Figure 5.

Figure 4 shows that the weights of each indicator calculated using the three methods were different, but the optimal combination weight is better distinguished. The optimal combination weights of the six positive indicators I1, I2, I3, I4, I5, and I7 were closer to the maximum value of their single weights. Conversely, the optimal combination weights of the two negative indicators I6 and I8 were closer to the minimum value of their single weights. Therefore, the optimal combination weights of the eight indicators in this paper fully reflected the importance of each indicator on China’s overall energy resilience.

Figure 5 clearly shows that, from 2005 to 2018, the optimal combination weight of each indicator had an unstable dynamic state of change. Before 2011, the weights of each indicator varied significantly in different years. After 2012, this characteristic changed significantly, and the average weight of the eight indicators was 0.13 in 2012 and 2013. After 2013, the changing trend of each indicator’s weight became obvious. The weights of the energy self-sufficiency, production diversity index, GDP energy intensity, and energy investment were higher than the other indicators. In particular, the weights of the energy investment indicators increased rapidly after 2014 and were far higher than the other indicators. The utilization time of power generation equipment and distribution loss decreased significantly compared with the prior years. On the one hand, this decrease reflected the improvement of power plant efficiency in China; on the other hand, it also reflected that these two indicators have reached a fairly high level and will not threaten the power interruption in China.

3.2. Suitability of the Proposed Improved Gray Relational Projection Model

This paper used the improved grey model to calculate the energy resilience ranking of 30 provinces in China, and compared it with the traditional grey model, which is seen in

Table 3. Energy resilience ranking of each province in China by comparison of two gray relational analyses.

Year	2005		2006		2007		2008		2009		2010		2011		2012		2013		2014		2015		2016		2017		2018
Province	R1	R2	R1	R2	R1	R2	R1	R2	R1	R2	R1	R2	R1	R2	R1	R2	R1	R2	R1	R2	R1	R2	R1	R2	R1	R2	R1	R2
Yunnan	9	6	8	11	14	4	5	4	9	12	7	9	3	1	8	2	2	4	4	29	3	24	9	27	7	20	1	11
Heilongjiang	10	4	7	7	5	12	6	20	10	14	12	27	6	20	9	9	10	2	12	2	10	8	11	8	8	11	10	19
Shaanxi	1	11	2	5	3	9	2	17	3	10	1	3	2	11	1	28	4	29	1	25	1	30	2	30	1	25	3	30
Shanghai	29	19	11	21	30	20	25	12	29	11	28	12	29	28	25	19	29	21	28	24	28	21	18	15	26	29	25	25
Chongqing	11	5	21	14	13	23	17	16	13	25	11	30	13	26	15	18	12	15	11	7	9	5	15	10	14	30	12	18
Qinghai	5	15	12	19	11	19	8	29	7	18	6	29	9	30	2	27	6	27	8	23	7	16	7	9	6	6	4	9
Inner Mongolia	6	9	5	2	1	5	10	9	5	3	3	5	5	27	6	13	1	17	7	27	8	4	4	4	5	8	8	5
Guizhou	12	8	13	10	15	6	11	6	15	1	10	1	11	6	28	7	11	14	10	10	12	25	14	25	9	3	11	16
Jilin	13	27	25	27	16	27	12	26	14	22	13	21	17	17	12	22	13	20	15	19	14	18	10	29	10	19	13	14
Zhejiang	30	25	27	30	8	14	30	7	24	5	29	6	28	3	27	6	30	3	29	21	30	14	30	19	28	14	30	22
Liaoning	14	12	4	8	6	22	13	15	16	24	16	24	12	14	11	11	14	10	16	18	13	26	12	26	11	26	18	28
Shandong	15	16	23	17	17	29	19	24	19	28	17	26	14	23	13	24	16	24	13	16	15	11	13	18	15	24	20	23
Henan	16	13	6	12	18	8	18	27	11	27	14	11	15	29	14	30	17	30	18	28	16	29	17	20	17	13	19	8
Shanxi	2	3	1	9	2	7	1	8	4	30	2	19	1	10	3	1	5	9	3	5	4	12	1	7	3	1	5	6
Hunan	27	28	19	29	24	21	24	13	23	17	27	14	27	21	29	17	28	26	26	22	27	23	16	21	29	17	26	12
Fujian	17	24	14	24	21	24	14	23	17	15	15	15	16	9	16	20	15	6	14	12	18	13	23	23	16	21	22	17
Guangdong	18	30	30	13	9	13	21	21	12	16	18	20	19	22	17	21	19	22	17	17	19	20	24	12	23	7	14	1
Hubei	28	26	26	26	29	30	29	25	25	9	22	25	30	19	26	16	26	16	27	8	29	9	19	28	25	28	29	27
Jiangsu	19	17	20	20	19	15	20	11	18	23	24	22	21	18	19	4	18	11	19	6	17	19	25	14	20	15	15	13
Sichuan	3	7	22	1	7	2	3	3	6	6	5	8	7	16	5	3	7	5	5	3	6	1	3	1	12	2	7	4
Guangxi	20	14	24	15	20	18	22	22	26	19	23	16	18	8	20	23	21	25	22	11	20	10	26	3	19	5	21	7
Ningxia	4	2	10	4	10	1	7	5	2	7	4	4	10	13	10	14	3	13	6	14	2	22	5	17	2	16	6	10
Xinjiang	7	1	17	6	12	3	4	2	1	4	9	2	4	4	4	26	8	23	2	26	5	28	6	24	4	22	2	21
Beijing	21	21	16	23	25	26	15	30	20	20	19	28	20	12	21	8	22	18	20	20	22	7	22	5	18	10	16	29
Gansu	26	18	15	28	28	17	28	28	28	26	26	18	23	24	30	25	24	19	30	15	26	17	27	16	27	9	27	3
Tianjin	22	22	9	22	22	25	16	10	21	2	20	7	22	5	18	12	20	7	21	4	24	2	21	6	22	12	17	15
Anhui	8	10	3	3	4	11	9	19	8	29	8	23	8	25	7	29	9	28	9	30	11	27	8	22	13	27	9	26
Hebei	23	23	29	25	27	10	23	1	27	8	25	10	26	2	23	5	25	8	24	9	23	3	20	2	21	4	23	2
Jiangxi	24	20	18	16	23	16	27	14	22	21	21	17	24	15	24	15	27	12	23	13	21	15	29	13	24	23	24	20
Hainan	25	29	28	18	26	28	26	18	30	13	30	13	25	7	22	10	23	1	25	1	25	6	28	11	30	18	28	24

(R1—energy resilience ranking of each province calculated by gray relational analysis, R2—energy resilience ranking of each province calculated by improved gray relational analysis). The results showed that the energy resilience of the five provinces of Shaanxi, Shanxi, Inner Mongolia, Ningxia, and Xinjiang ranked in the top five in most years, and Shaanxi and Shanxi ranked in the top two in most years. In contrast, in most years, Zhejiang, Hunan, Hubei, Gansu, and Hainan Provinces ranked in the bottom five. Except for Sichuan Province, the ranking results of the gray relational analysis and the improved gray relational analysis are the same. This result supports the suitability of the improved gray relational analysis constructed in this study, which is also in line with the actual situation.

Sichuan is a province with a high degree of energy diversification. Hydropower accounts for more than 75% of total energy consumption. It is the largest hydropower development base in China. In recent years, many natural gas resources have been discovered, which effectively guarantees the energy supply in the region and large amounts of energy are exported to other provinces each year. Moreover, in R2 ranking, Shanghai also ranked slightly higher than R1, which is in line with the actual situation. Although Shanghai lacks energy, it has a large energy investment and low energy consumption per unit GDP. In the past, when Shanghai suffered from extreme weather events, there was no energy shortage. Therefore, the improved grey relational analysis results are closer to the reality of each province in China, making them suitable for this article.

3.2.1. Energy Resilience of the 30 Chinese Provinces

As shown in Figure 6, from 2005 to 2018, the energy resilience of 30 provinces in China showed significant dynamics and differences. The energy resilience in most provinces is greater than 0.6, and the resilience level is above level II. The results indicate that if it is disturbed by uncertain emergencies, China’s overall energy system will not be interrupted and will be in a safe supply state. Specifically, the five provinces of Sichuan, Xinjiang, Shanxi, Inner Mongolia, and Shaanxi, show significantly higher resilience figures than the other provinces; the energy resilience of these provinces is above 0.8 in most years, indicating that the energy supply of these five provinces is sufficient and not affected by sudden disturbances. Except for Sichuan, Xinjiang, Shanxi, Beijing, Tianjin, Hebei, Inner Mongolia, and Shaanxi, the energy resilience of most provinces is unstable.

Energy resilience in Beijing, Tianjin, and Hebei has been maintained at around 0.6 for many years, reflecting the lack of significant improvement in energy resilience in these regions. The trend of energy resilience changes in the three provinces of Chongqing, Guizhou, and Yunnan, which showed a high degree of consistency before 2013, with a significant “W”-shaped trend, and after 2013, a significant “V”-shaped trend. Although after 2013 the volatility of these three provinces was high, in comparison, Yunnan Province has gradually improved since 2016. The ten provinces in eastern and southern China and Gansu Province in western China have significantly higher fluctuations than other provinces. The energy resilience of the nine provinces of Yunnan in the southwest, Xinjiang and Qinghai in the northwest, Inner Mongolia in the north, Heilongjiang and Jilin in the northeast, Anhui and Fujian in the east, and Henan in the central region experienced fluctuations before 2015, and they have seen an upward trend in the past three years. The energy resilience of Chongqing and Guizhou in the southwest, Gansu in the northwest, and Hainan and Shaanxi in the south have also declined slightly in recent years. The five provinces of Liaoning, Fujian, Jiangxi, Hunan, and Hubei showed a downward trend. If disturbed by uncertain factors, the energy systems of these regions will show significant vulnerabilities. In eight provinces in the eastern region and three provinces in the central region, energy resilience has been concentrated in the range of 0.4–0.6 for most years. If disturbed by sudden uncertain factors, these regions are prone to a short supply of energy with high external dependence for some time. In particular, in the five provinces of Chongqing, Gansu, Hainan, Hunan, and Zhejiang, the energy resilience in the sample interval was sometimes lower than 0.4, indicating that they are more susceptible to energy supply interruptions. These results require the attention of the government.

3.2.2. Spatial Evolution of the 30 Chinese Provinces

From the research results in Figure 6, we explored the spatial–temporal changes in the CERI scores across 30 Chinese provinces in 2005, 2010, 2015, and 2018 as shown in Figure 7, Figure 8, Figure 9 and Figure 10. The figures indicated that Jilin, Liaoning, Hebei, Henan, Hubei, Hunan, and Guangxi formed a clear dividing line from 2005 to 2010. Northwest of this dividing line were provinces with better energy resilience, which are also China’s energy-rich areas, except for Ningxia and Qinghai. In 2010, the energy resilience in provinces northwest of this dividing line did not change significantly, except for Qinghai and Ningxia. However, southeast of the dividing line, the number of provinces with energy resilience higher than 0.6 decreased, and the number of provinces with energy resilience lower than 0.4 increased. For example, the energy resilience in Guangdong and Hainan in 2010 improved compared with 2005, and the energy resilience of Shandong, Jiangsu, Henan, and the other provinces decreased to below 0.6 in 2010.

In 2015, the spatial–temporal patterns of China’s energy resilience changed significantly. In provinces that mainly rely on nonrenewable energy, such as Inner Mongolia and Qinghai, the energy resilience declined. Conversely, in provinces that focused on diversified energy development, such as Yunnan, Guangxi, and Anhui, the energy resilience has greatly improved. In 2018, the number of provinces with energy resilience 0.8–1 remained unchanged, while the number of provinces with energy resilience lower than 0.4 increased by one. The energy resilience of more than half of China, such as Shaanxi and Shanxi Provinces, is lower than in 2015. Among them, Guangxi and Chongqing have changed the most. The energy resilience of these two provinces was reduced from level II in 2015 to level III in 2018. From 2015 to 2018, China’s energy resilience has slowly declined, and its spatial distribution has shown significant imbalance. Regions with well-developed energy diversification have shown higher energy resilience.

4. Discussion

4.1. Key Factors Affecting Energy Resilience

This paper developed five dimensions of reliability, diversity, economic resilience, environmental resilience, and technological resilience to evaluate energy resilience; these five dimensions are intertwined. Most of the indicators of the five dimensions are helpful to improve energy resilience and can also help to improve the level of energy supply security [73]. Energy self-sufficiency is the key indicator of energy reliability and the basic element of energy resilience. Exner emphasizes the significance of energy self-sufficiency for energy reliability. In this paper, Shaanxi, Shanxi, Inner Mongolia, Xinjiang, and other provinces with abundant energy resources have a high energy supply security level and high energy resilience [74]. In contrast, Hainan Province and Gansu Province lack energy, and their energy resilience is low. Energy diversity is the key factor in increasing regional energy resilience and improving the adaptability and resistance of the energy system. In the context of climate change, improving regional energy diversity is of great significance for indirectly improving regional environmental resilience. In recent years, energy investment has been significant for energy resilience.

Weak energy infrastructure is the key factor influencing energy resilience, especially in underdeveloped areas in Hunan, Guizhou, Guangxi, etc. Therefore, it is vulnerable to extreme weather events, and the year of extreme weather events showed poor energy resilience. On the one hand, it reflects the lack of redundancy in China’s energy transformation policy, which leads to the lack of reliability of energy systems and seasonal supply interruption. On the other hand, it reflects that China’s energy transformation is faster than expected. Ultimately, this will be of great help to China’s achievement of the carbon-neutral goals for 2030 and 2060. Energy self-sufficiency, production diversity index, GDP energy intensity, and energy investment are the most significant factors affecting energy resilience in a region, and they are also interrelated, reflecting the complexity of energy resilience as a synergistic concept. It also shows the importance and necessity of using comprehensive evaluation methods to evaluate energy resilience.

4.2. Energy Resilience and Sustainable Development

The concept of resilience is still under development, and it is not static. Especially in the energy sector, the concept of energy resilience is undergoing continuous development and may involve significant energy policy issues, such as energy vulnerability, security, poverty, and justice. The concept of energy resilience proposed in this paper is consistent with the viewpoints of the Editorial of Nature Sustainability [75,76], and other institutions and scholars who believe that energy resilience cannot exist independently. Energy resilience is closely related to sustainable development factors such as economic development and social ecosystems. Compared with the energy system’s resilience, the concept of sustainable energy development mainly emphasizes the intergenerational equity of energy development and its impact on the environment, rather than the impact of external shocks on the energy system. Energy resilience refers to the adaptability and resistance of the energy system with full consideration of social, political, economic, and other uncertain factors. Energy resilience emphasizes the stability of the energy system during the impact.

4.3. Economic Development and Energy Resilience

Economic development and energy resilience complement each other. For the areas that lack energy resources, the rapid development of the economy will also impact the energy resilience of the region at the initial stage; a typical case is Guangdong Province. From 2005 to 2007, the energy resilience of most provinces in China was above 0.8, and those with energy resilience below 0.4 only appeared twice in Guangdong Province, where the economy was developing rapidly. This paper further collates the energy resilience trends of the top and bottom three provinces, as shown in Figure 11 and Figure 12, to confirm this result. The study results showed that from 2005 to 2018, the energy resilience of the top three cities in China, Beijing, Shanghai, and Tianjin, was approximately 0.6. All Shanghai’s energy resources were supplied from other places, but this city’s energy resilience is not the worst. However, it was evident that the energy resilience of the three provinces of Gansu, Yunnan, and Guizhou showed greater fluctuations, indicating that energy systems in economically underdeveloped regions are more susceptible to the impact of uncertainty. Energy resource endowment is the basic condition to ensure energy resilience; economic stability is dependent on the effective functioning and resilience of energy systems, and at the national level, investment scale played the dominant role [77]. The advantages of economic and social green development in the eastern region also help to improve energy resilience [78,79,80]. This is particularly noticeable in Germany, as their energy transition implies deep transformations in the economy, technology, and politics, and in the paramount fields of transition economies and rentier states [81]. Conversely, countries with poor economic development are more likely to cause people to worry about energy resilience [82].

Energy resilience is a complex scientific issue. Due to the difficulties of research methods and limitation of data, the results of this study have certain limitations. In the future, we will also continue to study it to explore more scientific decision-making methods for energy resilience.

5. Conclusions and Policy Implications

China is a typical representative developing country that uses coal as its energy source. Under the constraints of carbon neutrality goals, China is in a critical period of energy transition, and new problems may arise during the transition process, bringing new challenges to energy security. It is very important to investigate China’s energy resilience for China’s energy security and the energy security of coal-based developing countries such as China. This paper introduced the concept of regional energy resilience, presented the evolution mechanism and evaluation criteria of energy resilience, and constructed an improved gray relational projection model to measure the energy resilience of 30 provinces in China. We also analyzed the spatial evolution of energy resilience in various provinces and the relationship between energy resilience and regional economic development. The results indicated that the spatial–temporal patterns of energy resilience in China changed significantly from 2005 to 2018. High-energy resilience moved from provinces with abundant nonrenewable energy before 2010 to provinces with high energy diversity. Energy endowment is a primary condition to ensure a region’s energy resilience. Energy investment and economic development can effectively improve the energy resilience of resource-poor areas. Due to the limitation of data availability, the indicators selected in this paper cannot represent all the indicators for evaluating energy resilience. However, the results of this study are consistent with the actual conditions of various provinces in China, indicating the reliability of the research in this paper. The research results of this paper may serve as a reference for the Chinese government and other countries.

Although energy resilience is different from sustainable development, in recent years, the rapid decline in energy resilience in individual provinces in China is mainly due to energy transformation and extreme weather events. Climate issues have become the primary factor affecting energy resilience in China’s provinces. Therefore, China’s future provincial energy resilience construction plans should be combined with sustainable development goals and carbon-neutral goals. Factors such as energy self-sufficiency, production diversity index, and energy investment are the main factors affecting energy resilience. In the context of increasing uncertainty and interference, economic development and policy coordination are important guarantees for stabilizing regional energy resilience. Provincial disparities in economic development are significant. In some economically developed areas, despite the lack of energy resources, their policy is highly adaptable. When encountering extreme interference, they can maintain the stability of the energy system in the region. Therefore, while considering the source of energy supply, each province should strengthen policy coordination. Policy coordination among provinces and between provinces and energy companies should be strengthened. Provinces with developed economies but with a lack of energy should provide financial support to provinces with rich energy resources but poor infrastructure or energy. This can not only improve the energy resilience of these areas with poor infrastructure, but also ensure the energy supply of economically developed cities. At the same time, it is crucial to set up appropriate province-specific emission peak targets and raise province-specific emission reduction policies by considering the local realities.

With the development of the energy system towards complete de-carbonization, natural gas may play an important role in ensuring energy toughness with relatively low carbon characteristics and enough supply. Additionally, greater focus should be placed on the value of gas storage. Technological progress is an important method of improving energy resilience. According to the data analysis method, the Chinese government should establish an early warning mechanism for energy supplies, and collect and monitor data on energy supplies, energy consumption, and environmental emissions to optimize the frequency of energy data sharing. The Chinese government should also invest in exascale computing technology and advanced artificial intelligence methods for weather and climate forecasting so that the energy sector can better prepare for emergencies and extreme weather events. Moreover, the Chinese government should continue to develop carbon capture and storage technology, and general technology such as bioenergy to better conserve energy, protect the environment, and ameliorate the effects of extreme weather events and other events through the integration of multiple technologies to improve China’s energy resilience.