An Innovative Double-Frontier Approach to Measure Sustainability Efficiency Based on an Energy Use and Operations Management Perspective

Abstract

China’s economic development has achieved great success in recent years, but the problems of energy scarcity and environmental pollution have become increasingly serious. To enhance the reliability and efficiency between energy, the environment and the economy, sustainable development is an inevitable choice. In the context of measuring sustainability efficiency, a network data envelopment analysis model is proposed to formulate the two-stage process of energy use and operations management. A double frontier is derived to optimize the available energy for sustainable development. Due to nonlinearity, previous linear methods are not directly applicable to identify the double frontier and calculate stage efficiencies for inefficient decision-making units. To address this problem, this study develops the primal-dual relationship between multiplicative and envelopment network models based on the Lagrange duality principle of parametric linear programming. The newly developed approach is used to evaluate the sustainability efficiency of 30 administrative regions in China. The results show that insufficient sustainability efficiency is a systemic problem. Different regions should take different measures to conserve energy and reduce pollutant emissions for sustainable development. To increase sustainability efficiency, regions should support energy-saving and emission-reducing technologies in production processes and strengthen their capacity for technological innovation. Compared with energy use efficiency, operations management efficiency in China has a wider range of changes. During the operations management stage, there is not much difference between the capacity and quantity of each region. Based on benchmark regions at the efficiency frontier, there is an opportunity to improve operations management in the near future. Blockchain technology can effectively improve energy allocation efficiency.

1. Introduction

Although sustainable economic development in China has achieved remarkable results, the overall efficiency of energy use is still not high. Compared with 2015, China’s primary resource production rate increased by about 26% in 2020, energy consumption per unit of GDP continued to decline significantly, water consumption per unit of GDP decreased by 28%, and the overall recovery rate of large solid waste increased 56% [1]. It has been shown that the efficiency of energy use has been significantly improved. However, the efficiency of energy use is still at a low level. China’s energy consumption per unit of GDP is 1.5 times the global average, and total energy consumption is above the acceptable levels in many regions [2,3]. The main reason is that limited energy limits economic and social development [4]. In other words, the energy allocation mode does not meet the requirements of sustainability. Therefore, optimizing energy allocation is particularly important for China’s regional sustainable development.

Blockchain technology effectively improves sustainability efficiency through its unique decentralization, transparency, and immutability in the process of energy use [5]. First, blockchain mitigates the costs incurred by third-party institutions by allowing different merchants to conduct transactions directly, thereby significantly reducing transaction costs [6]. For example, in the fields of construction project management and fund payment, blockchain technology can enable decentralized data exchange and verification, reduce the dependence on third-party institutions in traditional models, and reduce operating costs. Second, through technologies such as smart contracts, blockchain can create effective incentive measures, encourage people to contribute their idle energy, achieve interaction and control in the energy allocation process, and significantly improve sustainability efficiency [7]. Third, blockchain is based on distributed storage, enabling real-time information sharing between multiple parties and improving sustainability efficiency [8]. Blockchain technology can ensure the lossless flow of data between the upstream and downstream of the supply chain, effectively avoid information distortion and distortion, and improve energy integration capabilities.

Data Envelopment Analysis (DEA), a data-driven benchmarking tool, is used to identify best practices or efficient decision-making units (DMUs) that apply multiple inputs to produce multiple outputs. This tool has been used since the work of Charnes et al. [9] and has been extensively applied to various practical problems. However, the standard DEA technique ignores the internal structure of DMUs and leads to misleading evaluation results [10,11]. From the perspective of opening the “black box” of DMUs, network DEA has been an important area in efficiency measurement [12]. By considering multi-stage processes with inputs, outputs, and intermediate measures, the efficiency relationship can be determined and the weaknesses between the internal components of DMUs can be identified. As a non-parametric technique for analyzing internal structures of DMUs, network DEA offers many advantages in assessing sustainability efficiency and energy allocation [13,14].

Studies on sustainability efficiency have been conducted in the literature of previous DEA networks. The assessment of supply chains, regions, transport, power distribution systems, and energy structures are the focus of relevant research. Kahi et al. [15] used a dynamic network DEA approach to assess the sustainability of supply chains. In the same way, Badiezadeh et al. [16] assessed the sustainability of supply chains using a big data methodology. An inverse network dynamic DEA model was used by Kalantary and Saen [17] to assess the sustainability of supply chains. Samavati et al. [18] went one step further and developed a dynamic network DEA model to assess the sustainability of supply chains. Azadi et al. [19] used the DEA network as an assessment framework to examine the resilience and sustainability of healthcare supply chains. A distinct investigation by Sarkhosh-Sara et al. [20] examined the sustainability of countries with different income levels using a network DEA model that considered null data and undesirable outputs. The sustainability of Chinese domestic transport was the focus of Stefaniec et al. [21], which used a network DEA approach based on the triple bottom line. The sustainability of Iran’s electricity distribution system was discussed in Tavassoli et al. [22] using a network DEA model. Geng et al. [23] discussed the global static and dynamic energy structure for energy use optimization using the total factor productivity approach based on slacks-based measurement. In their study, Lu et al. [24] used a modified dynamic parallel three-stage network DEA model to quantify the sustainability efficiency of climate change and global disasters based on greenhouse gas emissions from different parallel production sectors.

Energy allocation is one of the important aspects of sustainable development. Energy allocation should pay attention to the optimization and reuse of energy, which can not only reduce waste and energy loss, but also protect the environment and achieve a win–win situation for energy use and operations management [25]. In addition to sustainability efficiency, DEA also provides best-practice frontiers for energy allocation. DEA-based strategies for centralized energy allocation were initially put forth by Lozano and Villa [26]. Chen et al. [27] attempted to improve the input/output and developed an envelopment model to determine the efficient frontier for the basic two-stage network structure. Bi and others. [28] offered a goal-setting-focused DEA-based method for allocating energy. Chen et al. [29] showed that the envelopment network model should be used to determine the frontiers, while the multiplier network model should be used to determine the division efficiency. Wu et al. [30] used a DEA-based approach to allocate energy while taking the environment and economy into consideration. Kao [31] introduced frontiers identification for a relational model. Lim and Zhu [32] demonstrated that a simple two-stage network structure can be used to compute theoretical frontiers for the multiplicative and additive models. A lexicographic goal programming technique was created by Lozano and Contreras [33] for the Spanish public university system’s centralized energy allocation

Although such a study is crucial for sustainable development, these approaches cannot be directly applied to consider joint optimization of energy allocation for both energy use and operations management. During the 13th Five-Year Plan period, the total contribution of sustainable development to reducing carbon dioxide emissions in China exceeded 25% [34]. From 2012 to 2019, China’s average annual energy consumption growth of 2.8% supported the annual average energy use growth of 7% [35]. To increase the effectiveness of energy use and operations management working together for sustainable development, it is crucial to allocate energy as efficiently as possible.

Significantly, the above works focus on linear models. These linear techniques currently used to optimize energy allocation in standard DEA may not always be applicable to general network DEA models [36,37]. In this study, the theoretical framework of sustainability efficiency measurement is expanded to include the double frontiers of the two stages of energy use and operations management. Due to nonlinearity, the double frontier and stage efficiencies for inefficient DMUs cannot be determined directly using previous linear methods. The primal-dual relationship between multiplicative and envelopment network models based on the Lagrange duality principle of parametric linear programming is formulated. For the two-stage network process of sustainable development, the production possibility set is established under the multiplicative distance measures. The envelopment model can then be used to derive the frontiers and stage efficiency scores. The new method is illustrated using the sustainability efficiency assessment of thirty administrative regions in China.

The contribution of this study is as follows: (1) The derived envelopment network model is proven to be equivalent to the multiplier network model under the primal-dual relationship of parametric linear programming. (2) The impractical issues from the previous research can be avoided in the new approach. (3) It can optimize China’s energy allocation by merging the two double frontiers of operations management and energy use in the envelopment model.

The remainder of the paper is organized as follows. Section 2 proposes a double-frontier approach. Section 3 optimizes energy allocation for China’s regional sustainable development. Finally, the policy implications are presented in Section 4.

2. Methodology

2.1. The Framework of Sustainability Efficiency Measurement

In order to truly present the characteristics of low consumption, low emissions, and high efficiency of sustainability, sustainable development is divided into two stages: energy use and operations management. Figure 1 illustrates the general two-stage network process in 30 regions in China. In this setting, it is assumed each DMU_j (j = 1, 2, …, 30) has m inputs $x_{ij}$, (i = 1, 2, …, m), G outputs $y_{gj}^1$, (g = 1, 2, …, G), and D outputs $z_{dj}$, (d = 1, 2, …, D) in the stage of energy use. These D outputs referred as intermediate measures become the inputs in the stage of operations management, which uses some external inputs $x_{hj}^2$, (h = 1, 2, …, H) in addition to the intermediate measures. The outputs from this stage are $y_{rj}$, (r = 1, 2,…, s).

It is noted that some inputs or outputs in the above process, such as urban sewage discharge, exhaust gas discharge, and PM2.5, are all undesirable outputs. In particular, intermediate measures play a dual role, as they are not only undesirable outputs in the first stage but also inputs in the second stage. When constructing the DEA model, several methods were developed to deal with undesirable outputs or inputs, see Song et al. [38]. The undesirable factors are usually converted into negative numbers in network models [39]. However, this method leads to contradictions when optimizing the two stages of the intermediate measures. In addition, this method cannot guarantee that the multipliers of intermediate measures in network models have the same requirements in both stages [40,41]. In this study, the reciprocal undesirable output/input method can not only overcome the contradiction between the two stages, but also satisfy the assumption that the network model is the same for the intermediate measures in the two stages.

The following input-oriented general two-stage multiplicative network model is developed for DMUo in Figure 1.

$$\begin{array}{l}{e}_{\mathrm{o}}=\max {\displaystyle \frac{{\displaystyle \sum _{d=1}^D{\eta }_d{z}_{do}}+{\displaystyle \sum _{g=1}^G{\omega }_g{y}_{go}^1}}{{\displaystyle \sum _{i=1}^m{\nu }_i{x}_{io}}}}\ast {\displaystyle \frac{{\displaystyle \sum _{r=1}^s{u}_r{y}_{ro}}}{{\displaystyle \sum _{d=1}^D{\eta }_d{z}_{do}}+{\displaystyle \sum _{h=1}^H{Q}_h{x}_{ho}^2}}}\\ s.t. {\displaystyle \frac{{\displaystyle \sum _{d=1}^D{\eta }_d{z}_{dj}}+{\displaystyle \sum _{g=1}^G{\omega }_g{y}_{gj}^1}}{{\displaystyle \sum _{i=1}^m{\nu }_i{x}_{ij}}}}\le 1 \forall j \left(1\mathrm{a}\right)\\ {\displaystyle \frac{{\displaystyle \sum _{r=1}^s{u}_r{y}_{rj}}}{{\displaystyle \sum _{d=1}^D{\eta }_d{z}_{dj}}+{\displaystyle \sum _{h=1}^H{Q}_h{x}_{hj}^2}}}\le 1 \forall j \left(1\mathrm{b}\right)\\ v_i, {\eta }_d, {\omega }_g, Q_h, u_r\ge 0, \forall i, d, g, h, r \end{array}$$

(1)

In model (1), constraints (1a) and (1b) are used to model the first and second stages, respectively. $v_i, {\eta }_d, {\omega }_g, Q_h, u_r$ are weights. Note that model (1) cannot be solved using linear methods. Chen et al. [27] first proposed the solution of an additional equivalent envelopment model to determine the frontiers of model (1) for simple two-stage network process. Chen et al. [29] showed that multiplier and envelopment network models are two different methods for measuring the efficiency of two-stage network structures in Figure 1. Lim and Zhu [42] showed that the duality of standard DEA naturally transitions into model (1) for simple two-stage network process. They further point out that it is still unclear whether the same result can also be derived for general two-stage network processes.

Below, it is shown that the primal-dual correspondence can be extended to the general two-stage network structure. An envelopment network model is developed for model (1) in parametric linear form. The multiplier and envelopment network models are equivalent and avoid infeasible problems in Chen et al. [29]. By creating the production possibility set, the envelopment network model allows us to identify frontiers and stage efficiencies.

2.2. Envelopment Network Model and Duality

Model (1) is a nonlinear fractional model and cannot be converted to a linear model. Using the technique of Zhang and Chen [43], model (1) can be transformed into the following parametric linear model:

$$\begin{array}{l}{e}_{\mathrm{o}}=\max {\theta }_o^1{\displaystyle \sum _{r=1}^s{u}_r{y}_{ro}}\\ s.t. {\displaystyle \sum _{d=1}^D{\eta }_d{z}_{dj}}+{\displaystyle \sum _{g=1}^G{\omega }_g{y}_{gj}^1}-{\displaystyle \sum _{i=1}^m{\nu }_i{x}_{ij}}\le 0 \forall j \left(2\mathrm{a}\right)\\ {\displaystyle \sum _{r=1}^s{u}_r{y}_{rj}}-{\displaystyle \sum _{d=1}^D{\eta }_d{z}_{dj}}-{\displaystyle \sum _{h=1}^H{Q}_h{x}_{hj}^2}\le 0 \forall j \left(2\mathrm{b}\right)\\ {\displaystyle \sum _{d=1}^D{\eta }_d{z}_{do}}+{\displaystyle \sum _{g=1}^G{\omega }_g{y}_{go}^1}-{\theta }_o^1{\displaystyle \sum _{i=1}^m{\nu }_i{x}_{io}}=0 \left(2\mathrm{c}\right)\\ {\displaystyle \sum _{d=1}^D{\eta }_d{z}_{do}}+{\displaystyle \sum _{h=1}^H{Q}_h{x}_{ho}^2}-1=0 \left(2\mathrm{d}\right)\\ v_i, {\eta }_d, {\omega }_g, Q_h, u_r\ge 0, \forall i, d, g, h, r \end{array}$$

(2)

where ${\theta }_o^1$ is the efficiency of stage 1 and can be treated as a parameter in $[{\theta }_o^{1\min }, {\theta }_o^{1\max }]$. By using leader–follower models in Guo et al. [36], it can be determined the lower and upper bounds of ${\theta }_o^1$ as $[{\theta }_o^{1\min }, {\theta }_o^{1\max }]$. ${\theta }_o^1$ can be treated as a parameter in $[{\theta }_o^{1\min }, {\theta }_o^{1\max }]$. This makes model (2) a parametric linear model with parameter ${\theta }_o^1$. Model (2) is a linear model when ${\theta }_o^1$ is set as a constant in very linear search.

Based on the Lagrange duality principle of linear programming, the following dual model is derived for very fixed ${\theta }_o^1$:

$$\begin{array}{l}{e}_{\mathrm{o}}=\min t_2\\ s.t. {\displaystyle \sum _{j=1}^n{{\lambda }^{\prime }}_j}{x}_{ij}\le t_1{\theta }_o^1{x}_{io} i=1, \dots , m (3\mathrm{a})\\ {\displaystyle \sum _{j=1}^n({{\lambda }^{\prime }}_j}-{{\mu }^{\prime }}_j)z_{dj}\ge (t_1-t_2)z_{do} d=1, \dots , D (3\mathrm{b})\\ {\displaystyle \sum _{j=1}^n{{\lambda }^{\prime }}_j}{y}_{gj}^1\ge t_1{y}_{go}^1 g=1, \dots , G (3\mathrm{c})\\ {\displaystyle \sum _{j=1}^n{{\mu }^{\prime }}_j}{x}_{hj}^2\le t_2{x}_{ho}^2 h=1, \dots , H (3\mathrm{d})\\ {\displaystyle \sum _{j=1}^n{{\mu }^{\prime }}_j}{y}_{rj}\ge {\theta }_o^1{y}_{r_{2}o} r=1, \dots , s (3\mathrm{e})\\ {{\lambda }^{\prime }}_j, {{\mu }^{\prime }}_j\ge 0 \forall j; t_1, t_2 \mathrm{are} \mathrm{free} \mathrm{in} \mathrm{sign}\end{array}$$

(3)

Note that model (3) is also a parametric linear model about parameter ${\theta }_o^1$ and can be solved using the parametric linear method. Since ${\theta }_o^1>0$ and $t_1>0$, set ${\lambda }_j={{\lambda }^{\prime }}_j/t_1$ and ${\mu }_j={{\mu }^{\prime }}_j/{\theta }_o^1$. The following model is derived:

$$\begin{array}{l}{e}_{\mathrm{o}}=\min t_2\\ s.t. {\displaystyle \sum _{j=1}^n{\lambda }_j}{x}_{ij}\le {\theta }_o^1{x}_{io} i=1, \dots , m \left(4\mathrm{a}\right)\\ {\displaystyle \sum _{j=1}^n(\frac{{\lambda }_j}{{\theta }_o^1}}-{\displaystyle \frac{{\mu }_j}{t_1}})z_{dj}\ge ({\displaystyle \frac{1}{{\theta }_o^1}}-{\displaystyle \frac{t_2}{t_1\ast {\theta }_o^1}})z_{do} d=1, \dots , D \left(4\mathrm{b}\right)\\ {\displaystyle \sum _{j=1}^n{\lambda }_j}{y}_{gj}^1\ge y_{go}^1 g=1, \dots , G \left(4\mathrm{c}\right)\\ {\displaystyle \sum _{j=1}^n{\mu }_j}{x}_{hj}^2\le {\displaystyle \frac{t_2}{{\theta }_o^1}}{x}_{ho}^2 h=1, \dots , H \left(4\mathrm{d}\right)\\ {\displaystyle \sum _{j=1}^n{\mu }_j}{y}_{rj}\ge y_{ro} r=1, \dots , s \left(4\mathrm{e}\right)\\ {\lambda }_j, {\mu }_j\ge 0 \forall j\end{array}$$

(4)

Let ${\theta }_o^2={\displaystyle \frac{t_2}{{\theta }_o^1}}$ and $f={\displaystyle \frac{{\theta }_o^1}{t_1}}$, and model (4) becomes

$$\begin{array}{l}{e}_{\mathrm{o}}=\min {\theta }_o^1\ast {\theta }_o^2\\ s.t. {\displaystyle \sum _{j=1}^n{\lambda }_j}{x}_{ij}\le {\theta }_o^1{x}_{io} i_1=1, \dots , m \left(5\mathrm{a}\right)\\ {\displaystyle \sum _{j=1}^n{\lambda }_j}{z}_{dj}-z_{do}\ge f\ast ({\displaystyle \sum _{j=1}^n{\mu }_j}{z}_{dj}-{\theta }_o^2{z}_{do}) d=1, \dots , D \left(5\mathrm{b}\right)\\ {\displaystyle \sum _{j=1}^n{\lambda }_j}{y}_{gj}^1\ge y_{go}^1 g=1, \dots , G \left(5\mathrm{c}\right)\\ {\displaystyle \sum _{j=1}^n{\mu }_j}{x}_{hj}^2\le {\theta }_o^2{x}_{ho}^2 h=1, \dots , H \left(5\mathrm{d}\right)\\ {\displaystyle \sum _{j=1}^n{\mu }_j}{y}_{rj}\ge y_{ro} r=1, \dots , s \left(5\mathrm{e}\right)\\ {\lambda }_j, {\mu }_j\ge 0 \forall j; f>0\end{array}$$

(5)

In model (5), the constraints (5a) and (5d) are used to model the inputs $x_{ij}$ and $x_{hj}^2$, respectively. The constraints (5c) and (5e) are used to model the outputs $y_{ij}$ and $y_{gj}^1$, respectively. The constraint (5b) is used to model the intermediate measures $z_{dj}$ between the two stages. Based on the strong duality theorem, ${\theta }_o^{1*} \mathrm{and} {\theta }_o^{2*}$ are two-stage efficiencies, and $e_{\mathrm{o}}={\theta }_o^{1*}\ast {\theta }_o^{2*}$ is the system efficiency score. The inputs $x_{io}$ to stage 1 should be proportionally reduced to ${\theta }_o^{1\ast }{x}_{io}$ when the slack variables are zero in constraint (5a). Similarly, the inputs $z_{do}$ and $x_{ho}^2$ to stage 2 should be proportionally reduced to ${\theta }_o^{2*}{z}_{do}$ and ${\theta }_o^{2*}{x}_{ho}^2$ when the slack variables are zero in constraints (5b) and (5d). The new variable $f$ reflects the relationship of the intermediate measures between stage 1 and stage 2. Chen et al. (2013) [29] assumed that the intermediate measures would be treated as a “fixed link”, and ${\displaystyle \sum _{j=1}^n{\lambda }_j}{z}_{dj}={\displaystyle \sum _{j=1}^n{\mu }_j}{z}_{dj}$. The stage efficiencies between multiplier and envelopment models become unfeasible as a result of their approach. It is evident that this approach can avoid the impractical issues since models (1) and (5) are equal.

In standard DEA, an envelopment model can be obtained by defining its production possibility set [44]. Assume that $z_{do}$ are freely disposed of and define the following production possibility set for the general two-stage network structure:

$$\begin{array}{l}PPS=\{[(x_{io}, z_{do}^1, y_{go}^1), (z_{do}^2, x_{ho}^2, y_{ro})]\\ |{\displaystyle \sum _{j=1}^n{\lambda }_j{x}_{ij}\le x_{io}}, i=1, 2, \dots , m;{\displaystyle \sum _{j=1}^n{\lambda }_j{z}_{dj}^1}\ge z_{do}^1, d=1, 2, \dots , D; \left(6\mathrm{a}\right)\\ {\displaystyle \sum _{j=1}^n{\lambda }_j}{y}_{gj}^1\ge y_{go}^1 g=1, \dots , G; {\displaystyle \sum _{j=1}^n{\mu }_j{z}_{dj}^2}\le z_{do}^2, d=1, 2, \dots , D; \left(6\mathrm{b}\right)\\ {\displaystyle \sum _{j=1}^n{\mu }_j{x}_{hj}^2}\le x_{ho}^2, h=1, 2, \dots , H;{\displaystyle \sum _{j=1}^n{\mu }_j{y}_{rj}}\ge y_{ro}, r=1, 2, \dots , s; \left(6\mathrm{c}\right)\\ z_{do}^1\ge z_{do}^2; {\lambda }_j\ge 0; {\mu }_j\ge 0\}\end{array}$$

(6)

Normally, it holds $z_{do}^1=z_{do}^2$ for a normal observed data set for intermediate measures. In definition (6), $z_{do}^1\ge z_{do}^2$. With regard to intermediate measures, there is a distinction between the two stages. This is because the amount $z_{do}^1$ is firstly produced in stage 1 and then the amount $z_{do}^2$ (This is a portion of $z_{do}^1$) is consumed in stage 2. $z_{do}^1$ produced in stage 1 are disposed of before $z_{do}^2$ become a portion of inputs in stage 2.

Let $[({\theta }_o^1{x}_{io},z_{do}^1+f{\tilde{z}}_{do},y_{go}^1),({\theta }_o^2{z}_{do}^2+{\tilde{z}}_{do},{\theta }_o^2{x}_{ho}^2,y_{ro})]\in PPS$ and $z_{do}^1+f{\tilde{z}}_{do}\ge {\theta }_o^2{z}_{do}^2+{\tilde{z}}_{do}$. ${\tilde{z}}_{do}$ is free in sign. $f$ is a decision variable and a positive number.

${\theta }_o^{1*}$ is computed via the amount of the contraction of inputs $x_{io}$ to stage 1.

$$\begin{array}{l}{\theta }_o^{1*}=\min {\theta }_o^1\\ s.t.{\displaystyle \sum _{j=1}^n{\lambda }_j}{x}_{ij}\le {\theta }_o^1{x}_{io} i=1, \dots , m \left(7\mathrm{a}\right)\\ {\displaystyle \sum _{j=1}^n{\lambda }_j}{z}_{dj}^1\ge z_{do}^1+f{\tilde{z}}_{do} d=1, \dots , D \left(7\mathrm{b}\right)\\ {\displaystyle \sum _{j=1}^n{\lambda }_j}{y}_{gj}^1\ge y_{go}^1 g=1, \dots , G \left(7\mathrm{c}\right)\\ {\lambda }_j\ge 0 \forall j\end{array}$$

(7)

${\theta }_o^{2*}$ is computed via the amount of the contraction of inputs $z_{do}^2 \mathrm{and} x_{ho}^2$ to stage 2.

$$\begin{array}{l}{\theta }_o^{1*}=\min {\theta }_o^2\\ s.t.{\displaystyle \sum _{j=1}^n{\mu }_j}{z}_{dj}^2\le {\theta }_o^2{z}_{do}^2+{\tilde{z}}_{do} d=1, \dots , D\left(8\mathrm{a}\right)\\ {\displaystyle \sum _{j=1}^n{\mu }_j}{x}_{hj}^2\le {\theta }_o^2{x}_{ho}^2 h=1, \dots , H \left(8\mathrm{b}\right)\\ {\displaystyle \sum _{j=1}^n{\mu }_j}{y}_{rj}\ge y_{ro} r=1, \dots , s \left(8\mathrm{c}\right)\\ {\mu }_j\ge 0 \forall j \end{array}$$

(8)

Under the multiplicative efficiency aggregation method, the system efficiency $e_o$ is defined as the multiplicativity of the two-stage efficiencies ${\theta }_o^1\ast {\theta }_o^2$. The system efficiency can be calculated by the following model:

$$\begin{array}{l}\min {\theta }_o^1\ast {\theta }_o^2\\ s.t. \\ \mathrm{constraints} \mathrm{in} \mathrm{model} \left(7\right):\\ \{\begin{array}{l}{\displaystyle \sum _{j=1}^n{\lambda }_j}{x}_{ij}\le {\theta }_o^1{x}_{io} i=1, \dots , m \left(7\mathrm{a}\right)\\ {\displaystyle \sum _{j=1}^n{\lambda }_j}{z}_{dj}^1\ge z_{do}^1+f{\tilde{z}}_{do} d=1, \dots , D \left(7\mathrm{b}\right)\\ {\displaystyle \sum _{j=1}^n{\lambda }_j}{y}_{gj}^1\ge y_{go}^1 g=1, \dots , G \left(7\mathrm{c}\right)\\ {\lambda }_j\ge 0 \forall j\end{array}\\ \mathrm{constraints} \mathrm{in} \mathrm{model} \left(8\right):\\ \{\begin{array}{l}{\displaystyle \sum _{j=1}^n{\mu }_j}{z}_{dj}^2\le {\theta }_o^2{z}_{do}^2+{\tilde{z}}_{do} d=1, \dots , D \left(8\mathrm{a}\right)\\ {\displaystyle \sum _{j=1}^n{\mu }_j}{x}_{hj}^2\le {\theta }_o^2{x}_{ho}^2 h=1, \dots , H \left(8\mathrm{b}\right)\\ {\displaystyle \sum _{j=1}^n{\mu }_j}{y}_{rj}\ge y_{ro} r=1, \dots , s \left(8\mathrm{c}\right)\\ {\mu }_j\ge 0 \forall j \end{array}\end{array}$$

(9)

By removing ${\tilde{z}}_{do}$ from model (9), model (5) is derived when $z_{do}^1=z_{do}^2$.

2.3. Double Frontier and Stage Efficiency

From model (5), the frontier for DMUo can be determined as follows:

$$\begin{array}{ll}{\widehat{x}}_{io}={\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}{x}_{ij}={\theta }_o^{1\ast }{x}_{io} i=1, \dots , m& \left(10\mathrm{a}\right)\\ {\widehat{z}}_{do}^1={\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}{z}_{dj}=z_{do} d=1, \dots , D& \left(10\mathrm{b}\right)\\ {\widehat{y}}_{go}^1={\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}{y}_{gj}^1=y_{go}^1 g=1, \dots , G& \left(10\mathrm{c}\right)\\ {\widehat{z}}_{do}^2={\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{z}_{dj}={\theta }_o^{2\ast }{z}_{do} d=1, \dots , D& \left(10\mathrm{d}\right)\\ {\widehat{x}}_{ho}^2={\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{x}_{i_{2}j}^2={\theta }_o^{2\ast }{x}_{ho}^2 h=1, \dots , H& \left(10\mathrm{e}\right)\\ {\widehat{y}}_{ro}={\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{y}_{rj}=y_{ro} r=1, \dots , s& \left(10\mathrm{f}\right)\end{array}$$

(10)

(${\widehat{x}}_{io}, {\widehat{z}}_{do}^1, {\widehat{y}}_{go}^1$) in the Equations (10a)–(10c) forms the frontier of stage 1. Similarly, (${\widehat{x}}_{ho}^2, {\widehat{z}}_{do}^2, {\widehat{y}}_{ro}$) in the equations (10d)–(10f) forms the frontier of stage 2. Lim and Zhu [42] showed that the frontier of the intermediate measure $z_{do}$ would be any choice such that ${\widehat{z}}_{do}^1\ge {\widehat{z}}_{do}\ge {\widehat{z}}_{do}^2$ in Figure 1. Similarly, there is a possible difference for $z_{do}$ in two stages: ${\widehat{z}}_{do}^1\ge {\widehat{z}}_{do}^2$ in Formula (10). $z_{do}$ produced in stage 1 are disposed of before they become a portion of the inputs to stage 2.

As demonstrated below, Formula (10) is a frontier.

First, it will be proven that [(${\widehat{x}}_{io}, {\widehat{z}}_{do}^1, {\widehat{y}}_{go}^1$), (${\widehat{x}}_{ho}^2, {\widehat{z}}_{do}^2, {\widehat{y}}_{ro}$)] has a system efficiency score of unity. The efficiency of [(${\widehat{x}}_{io}, {\widehat{z}}_{do}^1, {\widehat{y}}_{go}^1$), (${\widehat{x}}_{ho}^2, {\widehat{z}}_{do}^2, {\widehat{y}}_{ro}$)] is evaluated using model (11):

$$\begin{array}{ll}\max {\theta }_o^1{\displaystyle \sum _{r=1}^s{u}_r{\widehat{y}}_{ro}}& \\ s.t. {\displaystyle \sum _{d=1}^D{\eta }_d{z}_{dj}}+{\displaystyle \sum _{g=1}^G{\omega }_g{y}_{gj}^1}-{\displaystyle \sum _{i=1}^m{\nu }_i{x}_{ij}}\le 0 \forall j& \left(11\mathrm{a}\right)\\ {\displaystyle \sum _{r=1}^s{u}_r{y}_{rj}}-{\displaystyle \sum _{d=1}^D{\eta }_d{z}_{dj}}-{\displaystyle \sum _{h=1}^H{Q}_h{x}_{hj}^2}\le 0 \forall j& \left(11\mathrm{b}\right)\\ {\displaystyle \sum _{d=1}^D{\eta }_d{\widehat{z}}_{do}^1}+{\displaystyle \sum _{g=1}^G{\omega }_g{\widehat{y}}_{go}^1}-{\displaystyle \sum _{i=1}^m{\nu }_i{\widehat{x}}_{io}}\le 0& \left(11\mathrm{c}\right)\\ {\displaystyle \sum _{r=1}^s{u}_r{\widehat{y}}_{ro}}-{\displaystyle \sum _{d=1}^D{\eta }_d{\widehat{z}}_{do}^2}-{\displaystyle \sum _{h=1}^H{Q}_h{\widehat{x}}_{ho}^2}\le 0& \left(11\mathrm{d}\right)\\ {\displaystyle \sum _{d=1}^D{\eta }_d{\widehat{z}}_{do}^1}+{\displaystyle \sum _{g=1}^G{\omega }_g{\widehat{y}}_{go}^1}-{\theta }_o^1{\displaystyle \sum _{i=1}^m{\nu }_i{\widehat{x}}_{io}}=0& \left(11\mathrm{e}\right)\\ {\displaystyle \sum _{d=1}^D{\eta }_d{\widehat{z}}_{do}^2}+{\displaystyle \sum _{h=1}^H{Q}_h{\widehat{x}}_{ho}^2}-1=0 \left(11\mathrm{f}\right)& \left(11\mathrm{f}\right)\\ v_i, {\eta }_d, {\omega }_g, Q_h, u_r\ge 0, \forall i, d, g, h, r & \end{array}$$

(11)

Consider a solution $({v^{\prime }}_i, {{\eta }^{\prime }}_d{{\omega }^{\prime }}_g, {Q^{\prime }}_h, {u^{\prime }}_r)=({\displaystyle \frac{v_i^{\ast }}{\gamma }}, {\displaystyle \frac{{\eta }_d^{\ast }}{\gamma }}, {\displaystyle \frac{{\omega }_g^{\ast }}{\gamma }}, {\displaystyle \frac{Q_h^{\ast }}{\gamma }}, {\displaystyle \frac{u_r^{\ast }}{\gamma }})$, where $(v_i^{\ast }, {\eta }_d^{\ast }, {\omega }_g^{\ast }, Q_h^{\ast }, u_r^{\ast })$ is an optimal solution to model (2) and $\gamma ={\displaystyle \sum _{j=1}^n{\displaystyle \sum _{d=1}^D{\eta }_d^{\ast }{z}_{dj}{\mu }_j^{\ast }}}+{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{h=1}^H{Q}_h^{\ast }{x}_{hj}^2{\mu }_j^{\ast }}}>0$. Plugging this solution into constraints (11a) and (11b), then

$$\begin{array}{ll}{\displaystyle \sum _{d=1}^D{{\eta }^{\prime }}_d{z}_{dj}+{\displaystyle \sum _{g=1}^G{{\omega }^{\prime }}_g{y}_{gj}^1}-{\displaystyle \sum _{i=1}^m{{\nu }^{\prime }}_i{x}_{ij}}}={\displaystyle \frac{1}{\gamma }}({\displaystyle \sum _{d=1}^D{\eta }_d^{\ast }{z}_{dj}+{\displaystyle \sum _{g=1}^G{\omega }_g^{\ast }{y}_{gj}^1}-{\displaystyle \sum _{i=1}^m{\nu }_i^{\ast }{x}_{ij}}})\le 0& \forall j\\ {\displaystyle \sum _{r=1}^s{u^{\prime }}_r{y}_{rj}}-{\displaystyle \sum _{d=1}^D{{\eta }^{\prime }}_d{z}_{dj}}-{\displaystyle \sum _{h=1}^H{Q^{\prime }}_h{x}_{hj}^2}={\displaystyle \frac{1}{\gamma }}({\displaystyle \sum _{r=1}^s{u}_r^{\ast }{y}_{rj}}-{\displaystyle \sum _{d=1}^D{\eta }_d^{\ast }{z}_{dj}}-{\displaystyle \sum _{h=1}^H{Q}_h^{\ast }{x}_{hj}^2})\le 0& \forall j\end{array}$$

Plugging $({v^{\prime }}_i, {{\eta }^{\prime }}_d, {{\omega }^{\prime }}_g, {Q^{\prime }}_h, {u^{\prime }}_r)$ into constraints (11c), (11d) and (11e), the following inequalities are derived:

$$\begin{array}{l} {\displaystyle \sum _{d=1}^D{{\eta }^{\prime }}_d{\widehat{z}}_{do}^1+{\displaystyle \sum _{g=1}^G{{\omega }^{\prime }}_g{\widehat{y}}_{go}^1}-{\displaystyle \sum _{i=1}^m{{\nu }^{\prime }}_i{\widehat{x}}_{io}}}\\ ={\displaystyle \frac{1}{\gamma }}[{\displaystyle \sum _{d=1}^D{\eta }_d^{\ast }({\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}{z}_{dj})}+{\displaystyle \sum _{g=1}^G{{\omega }^{\prime }}_g({\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}{y}_{gj}^1)}-{\displaystyle \sum _{i=1}^m{v}_i^{\ast }({\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}{x}_{ij})]}\\ ={\displaystyle \frac{1}{\gamma }}({\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }})({\displaystyle \sum _{d=1}^D{\eta }_d^{\ast }{z}_{dj}}+{\displaystyle \sum _{g=1}^G{\omega }_g^{\ast }{y}_{gj}^1}-{\displaystyle \sum _{i=1}^m{v}_i^{\ast }{x}_{ij})}\\ \le 0\end{array}$$

$$\begin{array}{l} {\displaystyle \sum _{r=1}^s{u^{\prime }}_r{\widehat{y}}_{ro}}-{\displaystyle \sum _{d=1}^D{{\eta }^{\prime }}_d{\widehat{z}}_{do}^2}-{\displaystyle \sum _{h=1}^H{Q^{\prime }}_h{\widehat{x}}_{ho}^2} \\ ={\displaystyle \frac{1}{\gamma }}[{\displaystyle \sum _{r=1}^s{u}_r^{\ast }({\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{y}_{rj})}-{\displaystyle \sum _{d=1}^D{\eta }_d^{\ast }({\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{z}_{dj})}-{\displaystyle \sum _{h=1}^H{Q}_h^{\ast }({\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{x}_{ij}^2)}]\\ ={\displaystyle \frac{1}{\gamma }}({\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }})({\displaystyle \sum _{r=1}^s{u}_r^{\ast }{y}_{rj}}-{\displaystyle \sum _{d=1}^D{\eta }_d^{\ast }{z}_{dj}}-{\displaystyle \sum _{h=1}^H{Q}_h^{\ast }{x}_{ij}^2)}\\ \le 0\\ {\displaystyle \sum _{d=1}^D{{\eta }^{\prime }}_d{\widehat{z}}_{do}^1}+{\displaystyle \sum _{g=1}^G{{\omega }^{\prime }}_g{\widehat{y}}_{go}^1}-{\theta }_o^1{\displaystyle \sum _{i=1}^m{{\nu }^{\prime }}_i{\widehat{x}}_{io}}\\ ={\displaystyle \frac{1}{\gamma }}[{\displaystyle \sum _{d=1}^D{\eta }_d^{\ast }({\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}{z}_{dj})}+{\displaystyle \sum _{g=1}^G{\omega }_g^{\ast }({\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}{y}_{gj}^1)}-{\theta }_o^{1\ast }{\displaystyle \sum _{i=1}^m{v}_i^{\ast }({\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}{x}_{ij})]}\\ ={\displaystyle \frac{1}{\gamma }}({\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }})({\displaystyle \sum _{d=1}^D{\eta }_d^{\ast }{z}_{dj}}+{\displaystyle \sum _{g=1}^G{\omega }_g^{\ast }{y}_{gj}^1}-{\theta }_o^{1\ast }{\displaystyle \sum _{i=1}^m{v}_i^{\ast }{x}_{ij})}\\ =0\end{array}$$

In addition, the equality for constraint (11f) is derived as follows:

$$\begin{array}{l}{\displaystyle \sum _{d=1}^D{{\eta }^{\prime }}_d{\widehat{z}}_{do}^2}+{\displaystyle \sum _{h=1}^H{Q^{\prime }}_h{\widehat{x}}_{ho}^2}-1={\displaystyle \frac{1}{\gamma }}[{\displaystyle \sum _{d=1}^D{\eta }_d^{\ast }({\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{z}_{dj})}+{\displaystyle \sum _{h=1}^H{Q}_h^{\ast }({\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{x}_{ij}^2)}]-1\\ ={\displaystyle \frac{1}{\gamma }}({\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }})({\displaystyle \sum _{d=1}^D{\eta }_d^{\ast }{z}_{dj}}+{\displaystyle \sum _{h=1}^H{Q}_h^{\ast }{x}_{ij}^2)}-1\\ =0\end{array}$$

The following equalities holds true for the objection function based on the complementary slackness condition:

$$\begin{array}{l}{\theta }_o^1{\displaystyle \sum _{r=1}^s{u^{\prime }}_r{\widehat{y}}_{ro}}={\displaystyle \frac{{\displaystyle \sum _{d=1}^D{{\eta }^{\prime }}_d{\widehat{z}}_{do}^1}}{{\displaystyle \sum _{i=1}^m{{\nu }^{\prime }}_i{\widehat{x}}_{io}}}}\ast {\displaystyle \sum _{r=1}^s{u^{\prime }}_r{\widehat{y}}_{ro}}={\displaystyle \frac{\frac{1}{\gamma }{\displaystyle \sum _{d=1}^D{\eta }_d^{\ast }({\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{z}_{dj})}}{\frac{1}{\gamma }{\displaystyle \sum _{i=1}^m{v}_i^{\ast }({\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}{x}_{ij})}}}\ast {\displaystyle \frac{1}{\gamma }}{\displaystyle \sum _{r=1}^s{u}_r^{\ast }({\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{y}_{rj})}\\ ={\displaystyle \frac{1}{\gamma }}{\displaystyle \sum _{r=1}^s{u}_r^{\ast }({\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{y}_{rj})}={\displaystyle \frac{1}{\gamma }}({\displaystyle \sum _{j=1}^n{\displaystyle \sum _{d=1}^D{\eta }_d^{\ast }{z}_{dj}{\mu }_j^{\ast }}}+{\displaystyle \sum _{j=1}^n{\displaystyle \sum _{h=1}^H{Q}_h^{\ast }{x}_{hj}^2{\mu }_j^{\ast }}})\\ =1\end{array}$$

Therefore, the solution $({v^{\prime }}_i, {{\eta }^{\prime }}_d, {{\omega }^{\prime }}_g, {Q^{\prime }}_h, {u^{\prime }}_r)$ and the multipliers [$({\widehat{x}}_{io}, {\widehat{z}}_{do}^1, {\widehat{y}}_{go}^1)$, (${\widehat{x}}_{ho}^2, {\widehat{z}}_{do}^2, {\widehat{y}}_{ro}$)] have a system efficiency score of unity in model (11).

It is demonstrated in the following that adding [$({\widehat{x}}_{io}, {\widehat{z}}_{do}^1, {\widehat{y}}_{go}^1)$, (${\widehat{x}}_{ho}^2, {\widehat{z}}_{do}^2, {\widehat{y}}_{ro}$)] to the DMU data set does not change the efficiency scores of the original DMUs.

The following model can be used to calculate the efficiency of an original DMU_k:

$$\begin{array}{ll}\max {\theta }_o^1{\displaystyle \sum _{r=1}^s{u}_r{y}_{rk}}& \\ s.t.{\displaystyle \sum _{d=1}^D{\eta }_d{z}_{dj}+{\displaystyle \sum _{g=1}^G{\omega }_g{y}_{gj}^1}-{\displaystyle \sum _{i=1}^m{\nu }_i{x}_{ij}}}\le 0 \forall j& \left(12\mathrm{a}\right)\\ {\displaystyle \sum _{r=1}^s{u}_r{y}_{rj}}-{\displaystyle \sum _{d=1}^D{\eta }_d{z}_{dj}}-{\displaystyle \sum _{h=1}^H{Q}_h{x}_{hj}^2}\le 0 \forall j& \left(12\mathrm{b}\right)\\ {\displaystyle \sum _{d=1}^D{\eta }_d{\widehat{z}}_{do}^1+{\displaystyle \sum _{g=1}^G{\omega }_g{\widehat{y}}_{go}^1}-{\displaystyle \sum _{i=1}^m{\nu }_i{\widehat{x}}_{io}}}\le 0& \left(12\mathrm{c}\right)\\ {\displaystyle \sum _{r=1}^s{u}_r{\widehat{y}}_{ro}}-{\displaystyle \sum _{d=1}^D{\eta }_d{\widehat{z}}_{do}^2}-{\displaystyle \sum _{h=1}^H{Q}_h{\widehat{x}}_{ho}^2}\le 0& \left(12\mathrm{d}\right)\\ {\displaystyle \sum _{d=1}^D{\eta }_d{z}_{dk}+{\displaystyle \sum _{g=1}^G{\omega }_g{y}_{gk}^1}-{\theta }_o^1{\displaystyle \sum _{i=1}^m{\nu }_i{x}_{ik}}}=0& \left(12\mathrm{e}\right)\\ {\displaystyle \sum _{d=1}^D{\eta }_d{z}_{dk}}+{\displaystyle \sum _{h=1}^H{Q}_h{x}_{hk}^2}-1=0& \left(12\mathrm{f}\right)\\ v_i, {\eta }_d, {\omega }_g, Q_h, u_r\ge 0, \forall i, d, g, h, r\end{array}$$

(12)

Assume that $(v_{ik}^{\ast }, {\eta }_{dk}^{\ast }, {\omega }_{gk}^{\ast }, Q_{hk}^{\ast }, u_{rk}^{\ast })$ is an optimal solution to model (2) for DMU_k. It suffices to demonstrate that $(v_{ik}^{\ast }, {\eta }_{dk}^{\ast }, {\omega }_{gk}^{\ast }, Q_{hk}^{\ast }, u_{rk}^{\ast })$ is an optimal solution to model (12). All that is needed to establish feasibility is that this solution meets the constraints (12c) and (12d). These constraints allow us to incorporate this solution, giving us the following:

$$\begin{array}{l}{\displaystyle \sum _{d=1}^D{\eta }_{dk}^{\ast }{\widehat{z}}_{do}^1+{\displaystyle \sum _{g=1}^G{\omega }_{gk}^{\ast }{\widehat{y}}_{go}^1}-{\displaystyle \sum _{i=1}^m{\nu }_{ik}^{\ast }{\widehat{x}}_{io}}}={\displaystyle \sum _{d=1}^D{\eta }_{dk}^{\ast }({\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}{z}_{dj})+{\displaystyle \sum _{g=1}^G{\omega }_{gk}^{\ast }({\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}{y}_{gj}^1)}-}{\displaystyle \sum _{i=1}^m{v}_{ik}^{\ast }({\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}{x}_{ij}})\\ ={\displaystyle \sum _{j=1}^n{\lambda }_j^{\ast }}({\displaystyle \sum _{d=1}^D{\eta }_{dk}^{\ast }{z}_{dj}}+{\displaystyle \sum _{g=1}^G{\omega }_{gk}^{\ast }{y}_{gj}^1}-{\displaystyle \sum _{i=1}^m{v}_{ik}^{\ast }{x}_{ij}})\le 0 \end{array}$$

$$\begin{array}{l}{\displaystyle \sum _{r=1}^s{u}_{rk}^{\ast }{\widehat{y}}_{ro}}-{\displaystyle \sum _{d=1}^D{\eta }_{dk}^{\ast }{\widehat{z}}_{do}^2}-{\displaystyle \sum _{h=1}^H{Q}_{hk}^{\ast }{\widehat{x}}_{ho}^2}={\displaystyle \sum _{r=1}^s{u}_{rk}^{\ast }({\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{y}_{rj})}-{\displaystyle \sum _{d=1}^D{\eta }_{dk}^{\ast }({\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{z}_{dj})}-{\displaystyle \sum _{h=1}^H{Q}_{hk}^{\ast }({\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}{x}_{hj}^2)}\\ ={\displaystyle \sum _{j=1}^n{\mu }_j^{\ast }}( {\displaystyle \sum _{r=1}^s{u}_{rk}^{\ast }{y}_{rj}}-{\displaystyle \sum _{d=1}^D{\eta }_{dk}^{\ast }{z}_{dj}}-{\displaystyle \sum _{h=1}^H{Q}_{hk}^{\ast }{x}_{hj}^2})\\ \le 0\end{array}$$

Given that the optimal objective value of model (12) for DMU_k is not greater than that of model (2), $(v_{ik}^{\ast }, {\eta }_{dk}^{\ast }, {\omega }_{gk}^{\ast }, Q_{hk}^{\ast }, u_{rk}^{\ast })$ is an optimal solution to model (12). The proof with these inspections is finished.

3. Empirical Study

3.1. Variables

In order to measure China’s regional sustainability efficiency in the new double-frontier approach, it is important to select appropriate inputs and outputs to accurately reflect the characteristics of China’s energy use efficiency and operations management efficiency.

(1) Inputs/outputs in energy use

Inputs: The inputs of this stage are capital, labor, and energy consumption. Coal, petroleum, and natural gas represent the main structure of energy consumption in China, but a dominant role in the consumption structure is played by coal, which is the main source of undesirable output such as sulfur dioxide. In this study, coal consumption is chosen as a substitution variable for energy consumption.

Outputs: Each region’s GDP is the desirable output. Undesirable outputs serve to comprehensively measure the level of pollution in a region and reflect the reduction principle of sustainable development. Reducing COD, SO₂, and nitrogen oxides are the main goals for energy conservation and emission reduction in China. Industrial waste gas is the sum of SO₂, nitrogen oxides, and dust emissions. COD, SO₂, and nitrogen oxides are the main pollutants in wastewater. In this study, wastewater and waste gas are therefore substitution variables for undesirable outputs. Given the availability and editability of the data, municipal sewage is the most important wastewater variable in the China Statistical Yearbook.

(2) Inputs/outputs in operations management

Inputs: The undesirable outputs become inputs of the operations management stage. The inputs of this stage also include wastewater and waste gas. After the destruction of the ecological environment, ecological restoration requires large investments in humans and energy. The inputs of this stage also include investments in air environment management and wastewater disposal. The investment in completing operations management is selected as a capital investment in the China Statistical Yearbook.

Outputs: In combination with the characteristics of the operations management stage, the indices are mainly selected from the two dimensions of ecological and social output. As an output index of ecological dimension, the number of days with good air quality is the desirable output index. This index indicates the positive direction that is beneficial for humans. In addition, the sewage treatment and PM_2.5 are two outputs of the social dimension. The sewage treatment is the desirable output, and PM_2.5 is the undesirable output. The two indices reflect to some extent the level of operations management.

3.2. Data Sources

This study focuses on China’s regional data from 2013 to 2021. During this period, China implemented reforms on both the supply and demand sides, and the economy transitioned from a period of rapid growth to a period of high-quality development. The data in this period can ensure timeliness. The input–output data for the sustainability efficiency assessment come from the China Statistical Yearbook, China Environmental Yearbook, Regional Statistical Yearbook, and Statistical Bulletin. The number of days with good air quality in the desirable outputs refers to the number of days on which the air quality reaches Grade II or better. Due to the incomplete air quality data of some cities, the number of days can be utilized when the air quality in the state capital is not worse than Grade II to replace the number of days with good air quality in these regions. There is a lack of some input–output data to evaluate the efficiency of operations management in Tibet, Hong Kong, Macau, and Taiwan, which are not considered in this study.

Table 1. Descriptive statistics.

Variables	Gross Fixed Capital Formation	Number of Employees	Coal Consumption	Urban Sewage Discharge	Exhaust Emissions	Completed Investment in Operations Management	GDP	Sewage Treatment	PM2.5	Number of Days with Good Air Quality
Number of DMUs	150	150	150	150	150	150	150	150	150	150
Minimum	2223.000	314.210	1720.330	16,496.000	9.530	21.100	2122.060	10,476.000	20.000	49.000
Maximum	116,628.000	6767.000	38,722.800	712,678.000	450.010	952.500	89,705.230	673,323.000	154.000	361.000
Median	10,709.550	2073.000	12,123.170	121,706.500	133.240	251.350	18,536.225	106,165.000	53.500	260.000
Mean	13,824.790	2726.329	14,566.256	154,947.927	153.366	302.081	24,422.983	141,643.693	57.126	253.255
Standard deviation	11,753.476	1783.981	8425.239	132,299.741	101.854	206.992	18,528.356	124,332.902	20.255	58.799

presents descriptive statistical data for 30 provinces in China in 2013, 2015, 2017, 2019, and 2021. It can be seen that there are significant differences between the maximum and minimum values of most indicators. For example, the minimum value of gross fixed capital formation is only 2223, while the maximum value exceeds 110,000; the median is only 10,709.55, and the mean and standard deviation do not exceed 15,000. The difference between urban sewage discharge and sewage treatment is also significant, with the maximum value exceeding 600,000, the minimum value not exceeding 20,000, and the median not exceeding 130,000, with the mean and standard deviation not exceeding 160,000. This suggests that there is an imbalance in energy use and energy allocation between regions in China, with a significant difference between coastal areas such as Jiangsu and Guangdong and regions such as Qinghai and Ningxia. Compared to other indicators, indicators such as coal consumption, exhaust emissions, and the number of days with good air quality perform relatively well in terms of data. Although there is still a significant difference between the maximum and minimum values, the mean and median remain at a moderate level, with a smaller standard deviation compared to the maximum value. This suggests that most regions in China have achieved effective control of carbon dioxide and other exhaust emissions, with a generally higher number of air quality compliance days, which can help protect the climate and quality of life in cities.

3.3. Results

Using the newly developed approach to evaluate China’s sustainability efficiency, the results are shown in

Table 2. Regional sustainability efficiency of China by using new method.

Regions	e					e1					e2
	2013	2015	2017	2019	2021	2013	2015	2017	2019	2021	2013	2015	2017	2019	2021
Beijing	0.34	0.20	0.24	0.20	0.24	1.00	1.00	1.00	1.00	1.00	0.34	0.20	0.24	0.20	0.24
Tianjin	0.42	0.30	0.48	0.97	0.80	0.95	1.00	0.98	0.97	0.95	0.44	0.30	0.49	1.00	0.84
Hebei	0.32	0.31	0.45	0.52	0.50	0.61	0.46	0.61	0.65	0.66	0.53	0.67	0.73	0.80	0.75
Shanxi	0.40	0.37	0.19	0.13	0.28	0.51	0.39	0.55	0.63	0.90	0.79	0.95	0.35	0.21	0.31
Inner Mongolia	0.63	0.65	0.63	0.52	0.51	0.63	0.65	0.63	0.58	0.57	1.00	1.00	1.00	0.90	0.90
Liaoning	0.40	0.48	0.61	0.84	0.80	0.60	0.61	0.84	0.92	0.90	0.67	0.78	0.73	0.91	0.89
Jilin	0.49	0.43	0.33	0.55	0.52	0.53	0.52	0.51	0.55	0.55	0.93	0.82	0.65	1.00	0.94
Heilongjiang	0.31	0.36	0.36	0.48	0.48	0.57	0.50	0.57	0.65	0.61	0.55	0.72	0.63	0.74	0.79
Shanghai	1.00	0.74	0.63	0.82	0.99	1.00	1.00	0.99	1.00	0.99	1.00	0.74	0.64	0.82	1.00
Jiangsu	0.59	0.61	0.65	0.80	0.70	0.80	0.78	0.83	0.91	0.87	0.74	0.78	0.78	0.88	0.80
Zhejiang	0.53	0.39	0.38	0.33	0.49	0.86	0.79	0.83	0.86	0.88	0.62	0.49	0.46	0.38	0.56
Anhui	0.22	0.29	0.21	0.25	0.49	0.73	0.65	0.70	0.75	0.74	0.30	0.45	0.30	0.33	0.66
Fujian	0.71	0.67	0.60	0.70	0.65	0.71	0.67	0.65	0.70	0.68	1.00	1.00	0.92	1.00	0.96
Jiangxi	0.45	0.47	0.39	0.37	0.34	0.80	0.74	0.80	0.77	0.78	0.56	0.64	0.49	0.48	0.44
Shandong	0.61	0.60	0.62	0.76	0.78	0.69	0.62	0.67	0.76	0.78	0.88	0.96	0.92	1.00	1.00
Henan	0.32	0.33	0.29	0.22	0.22	0.51	0.49	0.48	0.52	0.54	0.62	0.67	0.60	0.43	0.40
Hubei	0.36	0.22	0.27	0.23	0.26	0.67	0.61	0.62	0.67	0.64	0.53	0.36	0.43	0.34	0.40
Hunan	0.40	0.35	0.30	0.48	0.58	0.67	0.62	0.68	0.72	0.73	0.59	0.56	0.44	0.66	0.79
Guangdong	0.93	0.84	0.88	0.92	0.87	0.93	0.84	0.88	0.92	0.87	1.00	1.00	1.00	1.00	1.00
Guangxi	0.47	0.38	0.32	0.43	0.62	0.57	0.56	0.54	0.58	0.77	0.82	0.68	0.59	0.74	0.80
Hainan	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
Chongqing	0.50	0.38	0.53	0.47	0.41	0.73	0.69	0.70	0.76	0.74	0.69	0.55	0.76	0.62	0.55
Sichuan	0.39	0.34	0.43	0.40	0.42	0.73	0.65	0.74	0.79	0.77	0.54	0.53	0.58	0.51	0.55
Guizhou	0.59	0.25	0.43	0.46	0.32	0.59	0.54	0.54	0.56	0.60	1.00	0.46	0.79	0.83	0.54
Yunnan	0.48	0.44	0.42	0.44	0.43	0.48	0.44	0.42	0.44	0.43	1.00	1.00	1.00	1.00	1.00
Shaanxi	0.23	0.30	0.20	0.17	0.19	0.57	0.54	0.48	0.60	0.58	0.41	0.56	0.41	0.29	0.33
Gansu	0.39	0.31	0.39	0.41	0.76	0.69	0.67	0.55	0.56	0.81	0.57	0.47	0.70	0.73	0.94
Qinghai	0.59	0.59	0.69	0.26	0.77	1.00	1.00	1.00	1.00	1.00	0.59	0.59	0.69	0.26	0.77
Ningxia	0.71	0.53	0.45	0.39	0.29	0.71	0.71	0.72	0.71	0.83	1.00	0.75	0.63	0.55	0.35
Xinjiang	0.16	0.35	0.28	0.15	0.22	0.47	0.44	0.40	0.47	0.38	0.34	0.79	0.69	0.32	0.58
The average	0.50	0.46	0.46	0.48	0.53	0.71	0.67	0.70	0.73	0.75	0.70	0.68	0.65	0.66	0.70

. Shanghai, Guangdong, and Hainan have the highest sustainability efficiency. Hainan’s efficiency value is 1, which is always at the highest level. Shanghai’s efficiency value is 1 in 2013 and 2021. Guangdong’s efficiency value is relatively stable, at 0.93, 0.84, 0.88, 0.92, and 0.87. Although Guangdong’s efficiency value is not 1, it is still ahead of other regions. Guangdong, like Hainan and Shanghai, has demonstrated harmonious energy use and operations management in the past five years and is the most stable province in China. For Shanxi, Henan, Hubei, Shaanxi, Anhui, and Xinjiang, the efficiency values are relatively low.

In particular, Xinjiang’s efficiency score is the lowest in 2019. These regions are far from the efficient frontier in terms of energy use and operations management. They have no advantage for sustainable development. They are all distributed in the central and western areas. Shanxi, Henan, and Hubei are located in the central area and account for half of the total area of the central area. Other regions with low efficiency values show increasing trends. For example, Heilongjiang rose from 0.31 in 2013 to 0.48 in 2021, while Hebei rose from 0.32 in 2013 to 0.50 in 2021. The efficiency of sustainable development in Tianjin, Liaoning, and Gansu has changed greatly, Tianjin increased from 0.42 in 2013 to 0.97 in 2019, Liaoning increased from 0.40 in 2013 to 0.84 in 2019, and Gansu increased from 0.31 in 2015 to 0.76 in 2021.

Beijing, Tianjin, Shanghai, Jiangsu, Zhejiang, Guangdong, Hainan, and Qinghai performed well in energy use efficiency. Beijing, Hainan, and Qinghai are at the frontier of energy use, which means that these three regions make full use of energy, labor, and fixed assets. Although Tianjin, Shanghai, Jiangsu, Zhejiang, and Guangdong are not at the efficient frontier, they are stable and have good energy use. The efficiency values for Jilin, Henan, Yunnan, Shaanxi, and Xinjiang are relatively low at the energy use level. The industrial structures of these regions mainly comprise heavy industry and agriculture, and the technical level is obviously worse than that of the eastern areas.

The central area is responsible for extensive industrial relocation to the eastern area. Energy utilization tends to be extensive and pollution levels are high, making these regions inefficient at the stage of energy use. The efficiency of Shanxi, Liaoning, and Guangxi has changed greatly during this period. In particular, Shanxi’s efficiency value is 0.9 in 2021. This shows that these regions have adjusted their industrial structure and the supply-side structural reform has achieved remarkable results. For example, Shanxi is interested in financial investment and tax relief. In 2021, various types of taxes were reduced by CNY 31.115 billion, an increase of 43.57% year-on-year, including CNY 1.346 billion for high-tech incentives and CNY 2.034 billion for energy conservation and environmental protection. Liaoning considers the development of new generation information technology, high-end equipment manufacturing, biomedicine, and other strategic emerging industries as top priorities to accelerate supply-side structural reform.

Compared to the efficiency of sustainable development and energy use, the efficiency of operations management varies greatly between different regions, which are vulnerable to external factors. The efficiency of Fujian, Guangdong, Hainan, Yunnan, and Inner Mongolia is relatively high, while that of Beijing, Zhejiang, Hubei, Sichuan, and Shaanxi is relatively low. The efficiency of operations management in Beijing is very low. This is because Beijing overinvests in pollution removal, more than 4 times as much as Shanghai, which has a similar energy use efficiency, and 1.82 times as much as Guangdong, which has the highest energy consumption. In some regions, the level of operations management has changed. For example, Tianjin’s efficiency changed from 0.30 in 2015 to 1 in 2019.

Finally, the two-stage additive network DEA model in Guo et al. [36] is used to measure China’s regional sustainability efficiency. The results are shown in

Table 3. Regional sustainability efficiency of China by using Guo et al.’s [36] method.

Regions	e					e1					e2
	2013	2015	2017	2019	2021	2013	2015	2017	2019	2021	2013	2015	2017	2019	2021
Beijing	0.68	0.61	0.63	0.61	0.64	1	1	1	1	1	0.35	0.21	0.26	0.22	0.27
Tianjin	0.69	0.64	0.74	0.97	0.89	0.95	1	0.97	0.95	0.94	0.43	0.28	0.5	0.99	0.83
Hebei	0.57	0.57	0.67	0.73	0.71	0.61	0.46	0.61	0.65	0.66	0.53	0.67	0.73	0.8	0.75
Shanxi	0.65	0.67	0.45	0.42	0.61	0.51	0.39	0.55	0.63	0.9	0.79	0.95	0.35	0.21	0.31
Inner Mongolia	0.82	0.83	0.82	0.74	0.74	0.63	0.65	0.63	0.58	0.57	1	1	1	0.9	0.9
Liaoning	0.64	0.70	0.79	0.92	0.90	0.6	0.61	0.84	0.92	0.9	0.67	0.78	0.73	0.91	0.89
Jilin	0.73	0.67	0.58	0.78	0.75	0.53	0.52	0.51	0.55	0.55	0.93	0.82	0.65	1	0.94
Heilongjiang	0.56	0.62	0.60	0.70	0.70	0.57	0.51	0.57	0.65	0.61	0.55	0.73	0.63	0.74	0.79
Shanghai	1.00	0.87	0.82	0.91	1.00	1	1	0.99	1	0.99	1	0.74	0.64	0.82	1
Jiangsu	0.77	0.78	0.81	0.90	0.84	0.8	0.78	0.83	0.91	0.87	0.74	0.78	0.78	0.88	0.8
Zhejiang	0.74	0.64	0.65	0.62	0.72	0.86	0.79	0.83	0.86	0.88	0.62	0.49	0.46	0.38	0.56
Anhui	0.54	0.55	0.52	0.54	0.70	0.75	0.65	0.71	0.75	0.74	0.32	0.45	0.32	0.33	0.66
Fujian	0.86	0.84	0.79	0.85	0.82	0.71	0.67	0.65	0.7	0.68	1	1	0.92	1	0.96
Jiangxi	0.70	0.69	0.63	0.63	0.61	0.82	0.74	0.78	0.77	0.78	0.58	0.64	0.48	0.48	0.44
Shandong	0.79	0.79	0.80	0.88	0.89	0.69	0.62	0.67	0.76	0.78	0.88	0.96	0.92	1	1
Henan	0.57	0.58	0.54	0.48	0.47	0.51	0.49	0.48	0.52	0.54	0.62	0.67	0.6	0.43	0.4
Hubei	0.60	0.49	0.53	0.51	0.52	0.67	0.61	0.62	0.67	0.64	0.53	0.36	0.43	0.34	0.4
Hunan	0.63	0.59	0.56	0.69	0.76	0.67	0.62	0.68	0.72	0.73	0.59	0.56	0.44	0.66	0.79
Guangdong	0.97	0.92	0.94	0.96	0.94	0.93	0.84	0.88	0.92	0.87	1	1	1	1	1
Guangxi	0.70	0.62	0.57	0.66	0.79	0.57	0.56	0.54	0.58	0.77	0.82	0.68	0.59	0.74	0.8
Hainan	1.00	1.00	1.00	1.00	1.00	1	1	1	1	1	1	1	1	1	1
Chongqing	0.71	0.62	0.75	0.69	0.65	0.73	0.69	0.72	0.76	0.74	0.69	0.55	0.77	0.62	0.55
Sichuan	0.64	0.59	0.66	0.65	0.66	0.73	0.65	0.74	0.79	0.77	0.54	0.53	0.58	0.51	0.55
Guizhou	0.80	0.50	0.67	0.70	0.57	0.59	0.54	0.54	0.56	0.6	1	0.46	0.79	0.83	0.54
Yunnan	0.74	0.72	0.71	0.72	0.72	0.48	0.44	0.42	0.44	0.43	1	1	1	1	1
Shaanxi	0.49	0.55	0.45	0.47	0.46	0.57	0.54	0.48	0.62	0.58	0.41	0.56	0.41	0.31	0.33
Gansu	0.63	0.57	0.63	0.65	0.88	0.69	0.67	0.55	0.56	0.81	0.57	0.47	0.7	0.73	0.94
Qinghai	0.80	0.80	0.85	0.63	0.89	1	1	1	1	1	0.59	0.59	0.69	0.26	0.77
Ningxia	0.86	0.73	0.68	0.63	0.59	0.71	0.71	0.72	0.71	0.83	1	0.75	0.63	0.55	0.35
Xinjiang	0.41	0.62	0.55	0.40	0.48	0.47	0.44	0.4	0.47	0.38	0.34	0.79	0.69	0.32	0.58
The average	0.73	0.68	0.70	0.70	0.75	0.71	0.67	0.72	0.73	0.75	0.74	0.68	0.67	0.66	0.74

. The efficiency values of regional sustainability and the two stages for energy use and operations management are different from Guo et al. [36]. This is because the mechanisms for aggregating regional sustainability efficiency are different in the multiplicative network model and additive network model. However, it can be seen that the trend in efficiency value changes consistently over the years.

It is worth noting that Guo et al. [36] developed an improved golden section method to reduce the computational burden of line search. However, the derived optimal solution in their algorithm is still sensitive to computational precision in linear search, like the heuristic line search method. More importantly, it is difficult to see the changes in sustainability efficiency values and the two-stage efficiency values of energy use and operations management in Guo et al.’s [36] line search. It is obvious that the efficiency of energy use ${\theta }_o^1$ is the search parameter in this study, and the new method can reflect the varying sustainability efficiency and the two-stage efficiencies.

3.4. Findings

U-shaped change trends in efficiency. The efficiency averages for sustainable development, energy use, and operations management from 2013 to 2021 are 0.74, 0.76, and 0.73, respectively. Neither the sustainable development efficiency nor the two-stage efficiencies are high. There is still a lot of room for improvement. The efficiency of sustainable development, energy use, and operations management have similar changes in their characteristics and have U-shaped change trends. The efficiency of sustainable development reached its lowest point in 2017 and 2021.

After reaching its lowest point in 2015, the efficiency of energy use begins to increase every year. The efficiency of operations management lags behind the efficiency of energy use by a year, reaching its lowest point only in 2017 and then increasing from year to year. China’s economy entered a new normal in 2015, moving from high-speed growth to medium-high growth. The economic structure has now been constantly optimized and developed. Some industries have reduced production in the new and old shift periods, while pollutant emissions have not changed significantly; thus, the efficiency of energy use reached the lowest level in 2015. Improving energy use efficiency led to the lowest operations management efficiency in 2017. The change from the lowest efficiency of energy use and operations management led to the lowest sustainability efficiency in 2015 and 2017. From 2013 to 2021, the efficiency of operations management was higher than that of energy use only in 2015.

The efficiency of energy use is higher than that of operations management. In these periods, the efficiency of operations management is higher than that of energy use only in 2015, while that of sustainable development is somewhere in between. The efficiency of energy use and operations management has always been at the efficiency frontier. The efficiency of energy use in Beijing, Jiangsu, Zhejiang, Anhui, Jiangxi, Hubei, Sichuan, and Qinghai has always been higher than that of operations management. Beijing has the largest differences in the efficiency of energy use and operations management, both of which are above 0.65. In Tianjin, Shanghai, Chongqing, Hunan, and Shaanxi, the efficiency of energy use is lower than that of operations management in just one year, and the opposite is true in the remaining years. In Shanxi, Liaoning, and Ningxia, the efficiency of energy use is lower than that of operations management in 2013 and 2015 and higher than that of operations management after 2017. In Henan, the efficiency of energy use is higher than that of operations management after 2019. In Inner Mongolia, Jilin, Fujian, Shandong, Guangdong, Guangxi, and Yunnan, the efficiency of energy use has always been lower than that of operations management. Except for 2013 in Hebei and Heilongjiang, the efficiency of energy use was lower than that of operations management. Since 2015, the efficiency of energy use in Gansu has been lower than that of operations management. The efficiency of energy use and operations management has changed alternately in Guizhou and Xinjiang.

The prospects for sustainable development are still bleak, and the number of efficient DMUs for energy use and operations management is low. There are only 19 efficient DMUs for energy use, accounting for 12.7%. There are only three regions that were efficient at this time and accounted for 10%, namely Beijing, Hainan, and Qinghai. At this point, there are always 25 inefficient regions, which is 83.3%. The seven most affected regions are Jilin, Heilongjiang, Henan, Yunnan, Shaanxi, Guizhou, and Xinjiang. The efficiency of energy use in Shanxi and Guangxi was relatively low before 2019, but improved in 2019. There are only 29 efficient DMUs for operations management, accounting for 19.3%. At this time, only three regions were efficient, namely Guangdong, Hainan, and Yunnan. At this stage, 19 regions are always inefficient, which corresponds to 63.3%. The most serious regions are Beijing, Anhui, Hubei, and Shaanxi.

The eastern area is the best, followed by the western area, and the middle area is the worst. From a regional perspective, the efficiencies of sustainable development in eastern China are significantly higher than those in the central and western areas. The eastern area is superior to the central and western areas in terms of energy use and operations management. The efficiency of sustainable development, energy use and operations management is highest in the eastern area, followed by the western area, and the efficiency in the central area is lowest. The efficiency of energy use varies alternately in the central and western areas. The average efficiency of energy use in the eastern area varies between 0.80 and 0.88, while the average efficiency in the central and western areas varies between 0.56 and 0.69 and 0.61 and 0.68, respectively. In particular, the efficiency of energy use in the central area is lower than in the western area. This finding refutes the previous research finding that the performance of the central area is better than that of the western area.

3.5. Optimized Energy Allocation

The new approach is used to calculate the input–output improvements based on the data in 2021. The results are shown in

Table 4. Improvements of inputs in sustainable development of China in 2021.

Regions	Energy Use			Operations Management
	Gross Fixed Capital Formation	Number of Employees	Coal Consumption	Urban Sewage Discharge	Exhaust Emissions	Completed Investment in Operations Management
Beijing	0.0	0.0	0.0	553,750.6	57.7	503.76
Tianjin	535.0	47.2	412.7	19,163.9	5.1	11.48
Hebei	6437.5	1421.3	10,276.1	55,775.6	82.8	152.41
Shanxi	676.1	193.1	1698.9	162,241.9	335.6	191.17
Inner Mongolia	4430.7	607.5	8490.9	7271.5	17.5	41.66
Liaoning	969.4	229.8	2061.7	31,303.7	19.0	23.94
Gilling	4479.2	665.8	3585.2	6330.2	3.7	5.16
Heilongjiang	3761.9	783.9	4886.0	28,700.2	28.5	26.97
Shanghai	78.7	9.4	81.1	0.0	0.0	0.00
Jiangsu	4615.0	602.9	3983.1	107,640.1	43.0	143.84
Zhejiang	2689.2	466.9	2586.9	242,906.7	62.0	201.23
Anhui	3563.1	1149.6	3427.3	78,265.3	51.5	170.90
Fujian	5719.3	909.8	4179.7	5450.4	2.6	9.52
Jiangxi	2173.0	590.3	2007.2	114,830.4	106.2	175.17
Shandong	7804.0	1475.3	8698.9	0.0	0.0	0.00
Henan	13,915.2	3095.9	10,497.0	287,291.6	176.3	385.14
Hubei	7447.2	1305.9	6204.1	329,188.0	116.4	259.55
Hunan	4603.6	1024.1	4338.4	50,857.4	21.1	46.47
Guangdong	5168.6	853.7	4264.7	0.0	0.0	0.00
Guangxi	2107.6	663.0	2439.8	33,851.8	18.2	36.70
Hainan	0.0	0.0	0.0	0.0	0.0	0.00
Chongqing	2597.7	449.6	2215.3	90,013.0	43.4	98.92
Sichuan	3993.7	1100.0	3851.8	174,094.5	87.2	138.30
Guizhou	3626.1	807.5	3904.6	50,700.6	107.6	100.51
Yunnan	8478.9	1711.5	6343.1	0.0	0.0	0.00
Shaanxi	5917.5	866.9	5245.0	200,646.0	175.4	211.35
Gansu	658.9	287.8	1396.4	2649.9	4.4	5.69
Qinghai	0.0	0.0	0.0	4809.8	8.6	9.28
Ningxia	636.7	62.4	1077.0	45,603.6	102.9	54.75
Xinjiang	6665.2	814.9	11,466.6	49,262.3	94.9	162.13

. It can be seen that the improvements in Beijing, Hainan, and Qinghai in energy consumption are 0, because these three regions are at the efficiency frontier. Improvements are also very small in Shanghai, as well as in Tianjin, Shanxi, Liaoning, Ningxia, and Gansu. Shandong, Henan, and Yunnan occupy the top three places in terms of gross fixed capital formation and number of employees improvements, and the largest margin improvements are 7804.0, 13,915.2, 8478.9 and 1475.3, 3095.9, 1711.5, respectively. Xinjiang, Henan, and Hebei show the largest improvements in coal consumption. Coal consumption in Xinjiang can be reduced by 11,566.6. Likewise, the values of Henan and Hebei can be reduced by 10,497.0 and 10,276.1, respectively, closely followed by Inner Mongolia, Shandong, and Yunnan. This is closely related to the high energy consumption of heavy industry in these regions.

Investments in operations management in Beijing show a clear surplus. Urban sewage (553,750.6) and waste gas (57.7) should be treated, and investment in operations management (503.8) should be reduced in Beijing, followed by Henan and Hubei in terms of urban sewage discharge and investment in the stage of operations management. Shanxi has the largest need for improvement in exhaust emissions and should increase its disposal capacity by 335.6, followed by Henan and Shaanxi. Shanghai, Shandong, Guangdong, and Yunnan are at the efficiency frontier. Region improvements are 0.

The input improvements in the two stages are divided into five layers. Efficiency ratings are based on the ratio of actual improvement to maximum improvement. Improvement 0 is considered small. Then, the stage efficiencies are as follows: (0~0.20) is lower, (0.21~0.40) is normal, (0.41~0.80) is high, and (0.81~1) is higher. The regions with large improvements will make adjustments to this input to improve efficiency.

Table 5. Classifications of input improvements in inefficient regions.

Inputs\Levels	Low	Lower	Normal	Higher	High
Gross fixed capital formation	3	10	9	7	1
Number of employees	3	11	11	4	1
Coal consumption	3	8	10	6	3
Urban sewage discharge	5	17	4	3	1
Exhaust emissions	5	15	7	2	1
Completed investment in operations management	5	13	8	3	1

shows the number of regions in each layer for six inputs.

In the energy use stage, there are 17 regions with gross fixed capital formation at “normal”, “higher” and “high” levels, accounting for 63.0% of the inefficient regions. There are 15 regions with the number of employees at “normal”, “higher” and “high” levels, accounting for half of the inefficient regions. There are 19 regions with coal consumption at “normal”, “higher” and “high” levels, accounting for 73.0% of the inefficient regions. This shows that most inefficient regions at this stage have a lot of room for improvement in the three inputs, which is consistent with the result that the efficiency gap for energy use is large compared to the efficient frontier. The level of energy utilization needs to be improved at this stage.

The improvements in the operations management stage are slightly different from those in the energy use stage. Most regions are classified as “low”, “lower” and “normal”. Particularly in the inefficient regions, the proportion of completed investments in operations management is the lowest at two-thirds. The highest proportion of urban sewage discharge is in the inefficient regions at 81.5%. This suggests that there is not much difference in the amount and capacity of pollutant disposal at this stage. It is likely that the regions’ operations management efficiency will be improved in the short term according to the benchmark regions.

Blockchain technology can effectively improve energy allocation efficiency by reducing costs, improving operational efficiency, information sharing and transparency, and applying energy supply and distribution, energy trading and payment, and energy decentralization.

4. Conclusions and Policy Implications

As China’s economy transitions from high-speed growth to high-quality development, regions are vigorously pushing forward the transformation of the economic growth pattern. Due to limited energy, it is hoped that the larger the desired outputs, the smaller the undesirable outputs. The current study assesses China’s regional sustainability efficiency and optimizes energy allocation to improve its efficiency based on a double-frontier approach. From the perspective of input and output, sustainable development can be decomposed into the two stages of energy use and operations management. By establishing the dual equivalence between the multiplier and envelopment network models for the two-stage network structure, stage efficiencies and frontiers can be identified for inefficient DMUs. Then, the new approach is used to evaluate the sustainability efficiency and optimize energy allocation of 30 administrative regions of China during 2013–2021.

This study has several implications for sustainability development. Beyond energy use and operations management, the newly developed model can be extended to possible cases of economic development and environmental management, industrial production and waste disposal, technology development and economic application in the fields of supply chains, regional development, transport, power distribution systems, and energy structures. This newly developed method is a static sustainability efficiency measurement that simply ignores fluctuations in different time periods. The development of a dynamic multi-stage network model in the multi-period context is also an important research direction of the sustainability measurement efficiency method. It would also be interesting to explore the method of measuring sustainability efficiency in a hierarchical environment.

The policy implications of this study are as follows.

With increasing pressure on natural resources and the environment resulting from rapid urbanization and industrialization, the sustainability of China’s economic and social development has become a key concern. In response, the Chinese government has formulated a series of sustainability policies and set ambitious targets to ensure long-term prosperity without compromising environmental integrity.

The government’s aim is to significantly improve resource efficiency, reduce waste, and promote circular economy practices. One of the goals is to improve the efficiency of energy use, and the key is to improve the efficiency of operations management. By improving the efficiency of energy use and operations management, sustainability efficiency can be improved. The main reason for low sustainability efficiency is the combination of high inputs, undesirable outputs, and low desirable outputs. This suggests that low sustainability efficiency is still a structural problem. In order to improve sustainability efficiency, technological innovation capability, promote energy saving, and emission reduction technologies should be improved in production processes. The change range of operations management efficiency in China is larger than that of the energy use efficiency. The difference in the quantity and capacity of each region is not very large in the operations management stage. In the short term, there is an opportunity to improve operations management based on the benchmark regions at the efficiency frontier.

The government has implemented green industrial policies to promote green and low-carbon industries. These policies promote the adoption of clean technologies, energy efficiency improvements, and waste reduction measures. Different policies are promoted in different regions. The eastern region should play a pioneering role in sustainable development. Based on the advantages of the Yangtze River Delta, the Great Bay Area, and the coordinated development of Beijing, Tianjin, and Hebei, this area can enhance the ability of technological innovation and accelerate the development of a series of profitable, high-yield, and low-pollution industries and reduce pollution as much as possible. The central area is connected to urban metropolises such as the Yangtze Delta, the Pearl River Delta, and the urban agglomeration in the middle reaches of the Yangtze River. The area should make full use of the innovation and industrial spillover advantages of coastal areas, vigorously eliminate energy intensive industries, rehabilitate enterprises with significant environmental impact, improve the industrial level of central cities, and train talents in industrial clusters. The western area should promote the development of characteristic industries and promote the development of low-pollution and high-income industries that can meet local production conditions combined with local geographical environment, energy, and cultural customs.

The real-time updates of blockchain technology help achieve real-time information exchange between centralized and distributed energy sources, avoid the duplication of multiple energy sources, and reduce waste in energy supply systems [45] through facilitating a reasonable operation in energy supply systems [46,47]. In energy allocation, automatic execution and comprehensive sharing based on blockchain will greatly improve the rationality of energy allocation and solve the measurement problem of cross-conversion and spatio-temporal coupling between multiple energy sources in the energy system. Blockchain can help increase the level of decentralization in the energy industry. By enabling small providers to accept online payments quickly and easily, blockchain technology can increase the competitiveness and innovation of the energy industry, reduce energy prices, and provide consumers with more choice and cheaper prices.

There is a need to improve the statistical evaluation system of sustainable development. Many researchers hold the same view on the characteristics, connotations, and principles of sustainable development, but hold different views on the sustainable development evaluation system, resulting in complex and lengthy evaluation indices. This study presents the indices of undesirable outputs in the two stages of energy use and operations management based on network DEA, reflecting the principles of reduction, reuse, and recycling in sustainable development. However, how to reflect the evaluation system of the ten categories of renewable energy is a direction for future research, i.e., steel waste, non-ferrous metal waste, plastic waste, old tires, waste paper, old electrical and electronic equipment, old vehicles, old textiles, old glass, and old batteries.

China’s sustainability policies and goals are evidence of its commitment to achieving long-term prosperity while protecting the environment. The implementation of these policies and the pursuit of these goals are essential for the country’s sustainability development and its contribution to global efforts to combat climate change and environmental degradation.