Evaluation of Water—Energy—Food—Economy Coupling Efficiency Based on Three-Dimensional Network Data Envelopment Analysis Model

Abstract

In the process of human survival and sustainable development, water security, energy security and food security have become the three most prominent issues, but they are interrelated and directly affect each other, that is, to form a Water–Energy–Food (WEF) nexus. Scientific understanding and correct response to the relationship between WEF is important to realize the sustainable development of natural resources. There are some deficiencies in the existing research on the input–output efficiency of WEF system. There are few articles that can study the efficiency relationship between internal and external factors (such as the economy and environment) of the WEF system at the same time, or the research is not perfect. In view of the shortcomings of the existing research, this paper establishes a three-dimensional network structure to describe the water–energy–food–economy (WEF-Eco) system and establishes the corresponding network Data Envelopment Analysis (DEA) model. We use the data of 19 provinces in Northeast, East, and central China to show the application results of this model.

1. Introduction

Food production and agricultural development are inseparable from water resources and energy. Energy development also depends on the development of water resources and the food industry to a great extent. At the same time, the development and utilization of water resources need the support of energy and food. Any one of these three resources is affected by the other two resources. For example, the low reserves of oil, natural gas, coal, and other major energy resources in the Yangtze River Economic Belt basin; the decline in water quality in several water sources; and the serious degradation of ecological functions in the basin have gradually become important factors restricting the agricultural and economic sustainability of the region. If any one of these resources is short, the fragile balance between the three will be destroyed, which will lead to serious consequences. Especially in the context of global COVID-19 intensification, population growth, environmental degradation, and intensified impacts of climate change, the problem of resource shortage has had a major impact on global development. It is important to study the internal relationship and interaction between water–energy–food (WEF) [1].

Research on WEF mainly focuses on two aspects: (1) qualitative research, including elaborating the relationship between the internal WEF system or the relationship between the WEF system and external factors such as economy and environment [2,3,4,5,6,7], determining the boundary and core issues of the WEF system [3,8], and WEF research and analysis from different perspectives, such as collaboration, security, risk, and optimization [9,10,11,12]; and (2) quantitative research, focused on the construction of the operation framework of the WEF system and the screening of key indicators of the WEF system. The research methods involved mainly include Life Cycle Assessment (LCA) [13,14,15], Multiregional Input–Output (MRIO) model [16,17,18], Index System Method (ISM) [19,20], Data Envelopment Analysis (DEA) [21,22], System Dynamics (SD) model [23,24,25,26,27], Coupling Coordination Degree Model (CCDM) [28,29,30,31], Exploratory Spatial Data Analysis (ESDA) model [32], and Geographically Weighted Regression (GWR) [33].

DEA is to use mathematically program a model to evaluate the relative effectiveness (DEA effectiveness) between “departments” or “units” (called decision-making units, abbreviated as DMUs) with multiple inputs and outputs. According to the observed data of each DMU, it is judged whether a DMU is DEA effective or not. Essentially, it is judged whether a DMU is on the “frontier” of the production possibility set. A production frontier is a generalization of production function to multi-output situations in economics. The structure of a production frontier can be determined by the DEA method and model. Since the evaluation method was put forward, a large number of DEA models have been derived. The traditional DEA model [34,35] cannot meet the needs of more complex production systems due to its limitations. For example, the traditional DEA model regards the whole system as a “Black box”, ignoring the existence and differences of various subsystems that determine the internal functions of the system and the input–output relationship within the DMU. In addition, the traditional DEA model belongs to the self-evaluation mode, that is, each DMU selects a group of weights that are most beneficial to calculate its own efficiency [36]. In the self-assessment mode, many DMUs are effective, and effective DMUs cannot be further distinguished [37]. Aiming at the “Black box” problem, Färe et al. [38] constructed the network DEA model. After that, Kao et al. [39] modified the traditional DEA model by considering the sequence relationship of the two sub-processes in the entire process, decomposing the efficiency of the entire process into the product of the efficiency of the two sub-processes and proposed a two-stage network DEA model. Tone et al. [40] built a slacks-based network DEA model, which can formally deal with intermediate products. This scalar measure deals with the input excesses and the output shortfalls of the DMU concerned. It is units-invariant and monotonally decreasing with respect to input excess and output shortfall. Furthermore, this measure is determined only by consulting the reference-set of the DMU and is not affected by statistics over the whole data set. The SBM is a non-radial method and is suitable for measuring efficiencies when inputs and outputs may change non-proportionally. This model can decompose the overall efficiency into divisional ones. Chen et al. [41] established an additive two-stage network DEA model. In order to solve the problems existing in the self-evaluation efficiency model, Sexton et al. [42] proposed the cross-efficiency evaluation method to solve the shortcomings of the self-evaluation model, that is, each DMU maximizes its own efficiency through the traditional DEA model, and at the same time, uses a set of optimal weights of its input and output indicators to evaluate the efficiency of all DMUs. Doyle et al. [43] introduced two-level objective models in a cross-efficiency evaluation method to solve the problem that cross efficiency evaluation results are not unique. Wang et al. [44] proposed a neutral two-level objective model. Since then, many scholars have put forward other improved evaluation methods [45,46,47].

Due to the unique advantages of the DEA method, this method has been widely used in many fields. The DEA method can deal with the problem of multi-input and multi-output, and there is no need to build a production function to estimate the parameters. This method is not affected by the input–output dimension, uses the comprehensive index to evaluate the efficiency, is suitable to describe the situation of total factor production efficiency, and can compare the efficiency between DMUs. The weight of the DEA method is not affected by human subjective factors, and the evaluation of DMUs is relatively fair.

Because of the unique advantages of the DEA method, many scholars have applied the DEA method to the research of WEF. One type of research is to take single or multiple resources in WEF as input to explore the efficiency between them and the external elements of WEF. For example, Li et al. [21] considered the efficiency of taking the WEF system as input and the economy and environment as outputs. Sun et al. [48] divided the WEF system into “Water resources subsystem”, “Energy subsystem”, and “Food subsystem”; took water resources, energy, and food as part of the input indicators; and used the corresponding economic output indicators to calculate the efficiency of each subsystem separately. Li et al. [49] built a three-stage, dual-boundary network DEA (TD-NDEA) model and decomposed the WEF nexus into three stages, “W-E”, “WE-F”, and “WEF-Eco”; the efficiency of each stage and the overall efficiency of the system can be calculated separately.

The proposed DEA research methods for studying WEF problems have the following deficiencies:

(a) Some describe the structure of the WEF system as a series structure, such as taking water resource variables as input variables to produce energy variables, then taking the produced energy variables as input variables to produce food variables, and so on to produce other variables. This one-way production structure ignores some reverse input–output relationships, such as the input of energy variables in the production process of water resource variables. Therefore, there are the problems of lacking the information to describe the internal relationship of the WEF system with the series structure.

(b) Some consider three independent subsystems (the water resources subsystem, the energy subsystem, and the food subsystem), but they can only obtain the efficiency values of three independent subsystems and not the overall efficiency value of the WEF system.

(c) Others will only take the three indicators of water, energy, and food as inputs and take economic indicators as inputs to investigate the impact of water energy and food system on external factors, which ignores the investigation of the internal relationship of WEF system.

In the present paper, we explore a special comprehensive network structure and strive to comprehensively and accurately describe the internal relationship of the WEF system and the influence of the system on external factors. The main contributions of this paper are:

(1) We set up a new DEA model, namely, the three-dimensional network DEA model, and use this model to evaluate the efficiency of the WEF system in relevant regions.

(2) In contrast to the traditional additive DEA model, this paper proposes to determine the weight of the corresponding stage by the proportion of the total output of each stage in the total output of the system.

The remaining sections are structured as follows. In Section 2, we propose the three-dimensional network structure and the construction process of DEA models in detail. In Section 3, we investigate the WEF system of China’s 19 provincial administrative regions to show the calculation process of the model and analyze the calculation results. In Section 4, some conclusions are given, and future directions are pointed out.

Relevant nouns and their abbreviations are shown in

Table A1. Full names and abbreviations of nouns.

Full Name	Abbreviation	Full Name	Abbreviation
Water–Energy–Food	WEF	Life Cycle Assessment	LCA
Multiregional Input–Output	MRIO	Index System Method	ISM
Data Envelopment Analysis	DEA	System Dynamics	SD
Coupling Coordination Degree Model	CCDM	Exploratory Spatial Data Analysis	ESDA
Geographically Weighted Regression	GWR	Decision Making Unit	DMU
Slacks-based Model	SBM	Three-stage, dual-boundary Network DEA	TD-NDEA

in the Appendix A.

2. Methods

2.1. Network DEA Structure

According to the WEF system relationship, we establish the three-dimensional network structure diagram of WEF as shown in Figure 1. We can divide the structure into two stages. The bottom of the pyramid is the first stage. This stage includes the water resources subsystem, the energy subsystem, and the food subsystem. These three subsystems are in a parallel relationship. The arrows in the network structure diagram represent the production process, and the nodes are input or output indicators. Specifically, the water resources subsystem takes water as the output index and energy, investment, resources, and employed population as the input indices. Due to the lack of food input indicators that directly reflect the production needs of the water production industry, the possible food input indicators are the food consumption indicators of water production employees, which belong to indirect input indicators. Therefore, we do not consider the food input indicators in the water system. The energy subsystem takes energy as the output index and relevant water, investment, resources, and employed population as the input index. We do not consider the food input indicators required by the energy production industry because the energy generated by food accounts for a very small proportion of the total energy production, and the relevant indicators are not easy to statistically query. The food subsystem takes food as the output index and relevant water, energy, investment, resources, and employed population as the input indices. In the second stage, from the bottom of the pyramid to the top, economic indicators are taken as the output, and the three outputs of the first stage and related investment, resources, and employed population are taken as the input indicators.

Suppose there are n DMUs, and all indicators and corresponding weights of DMU_k (k = 1, 2, …, n) are shown in

Table 1. The notations and their meaning.

Notations	Definitions
E1k,i	Efficiency value of ith subsystem in the first stage of the DMUk (i = 1, 2, 3)
Ejk	Efficiency of the jth stage of the DMUk (j = 1, 2)
Ekk	Overall efficiency of the DMUk (self-assessment)
ω1k,i	Weight of ith subsystem in the first stage of the DMUk (i = 1, 2, 3)
ωjk	Weight of the jth stage of the DMUk (j = 1, 2)
τik,1	Weight of investment index of the ith subsystem of the first stage of the DMUk (i = 1, 2, 3)
τik,2	Weight of resources index of the ith subsystem of the first stage of the DMUk (i = 1, 2, 3)
τik,3	Weight of labor force index of the ith subsystem of the first stage of the DMUk (i = 1, 2, 3)
v1k,1	Weight of energy input index in the water resources subsystem of the first stage of the DMUk
v2k,1	Weight of water resources input index in the energy subsystem of the first stage of the DMUk
v3k,1	Weight of water resources input index in the food subsystem of the first stage of the DMUk
v3k,2	Weight of energy input index in the food subsystem of the first stage of the DMUk
uik	Weight of output index of the ith subsystem of the first stage of the DMUk (i = 1, 2, 3)

.

2.2. DEA Model

We first calculate the efficiency of the three subsystems in the first stage. Different from the traditional definition of efficiency, we define the ratio of input to output as the efficiency value. The notations and their meanings are shown in

Table 2. Efficiency evaluation index system of the WEF-Eco system.

Stage	Category	Index	Unit	Notation	Weight
W-E-F	Water supply subsystem	Total electric power consumption of water production and supply industry	100 million kWh	x1k,1	v1k,1
		Investment in water production and supply	100 million yuan	z1k,1	τ1k,1
		Total regional water resources	100 million cu.m	z1k,2	τ1k,2
		Employment in water production and supply	10,000 persons	z1k,3	τ1k,3
		Total water supply	100 million cu.m	y1k	u1k
	Energy production subsystem	Water consumption of electric power and heat power production and supply industry	100 million cu.m	x2k,1	v2k,1
		Investment in the production and supply of electric power and heat power	100 million yuan	z2k,1	τ2k,1
		Employment in the production and supply of electric power and heat power	10,000 persons	z2k,2	τ2k,2
		Total energy production	10,000 tons of SCE	y2k	u2k
	Food production subsystem	Total agricultural water consumption	100 million cu.m	x3k,1	v3k,1
		Total Energy consumption of agriculture, forestry, animal production and hunting	10,000 tons of SCE	x3k,2	v3k,2
		Investment in primary industry	100 million yuan	z3k,1	τ3k,1
		Grain sowing area	1000 hectares	z3k,2	τ3k,2
		Number of employees in the primary industry	10,000 persons	z3k,3	τ3k,3
		Total grain production	10,000 tons	y3k	u3k
WEF-Eco	Economic development subsystem	Total water supply	100 million cu.m	x4k,1	v4k,1
		Total energy production	10,000 tons of SCE	x4k,2	v4k,2
		Total grain production	10,000 tons	x4k,3	v4k,3
		Total investment in fixed assets	100 million yuan	z4k,1	τ4k,1
		R&D number of Industrial Enterprises above Designated Size	10,000 persons	z4k,2	τ4k,2
		Number of employed population at the end of the year	10,000 persons	z4k,3	τ4k,3
		Per capita GDP	100 million yuan	y4k	u4k

.

The efficiency values of these three subsystems are:

$$E_{1k,1}=\frac{{\displaystyle \sum _{i=1}^3{\tau }_{1k,i}{z}_{1k,i}}+v_{1k,1}{x}_{1k,1}}{u_{1k}{y}_{1k}}$$

(1)

$$E_{1k,2}=\frac{{\displaystyle \sum _{i=1}^2{\tau }_{2k,i}{z}_{2k,i}}+v_{2k,1}{x}_{2k,1}}{u_{2k}{y}_{2k}}$$

(2)

$$E_{1k,3}=\frac{{\displaystyle \sum _{i=1}^3{\tau }_{3k,i}{z}_{3k,i}}+v_{3k,1}{x}_{3k,1}+v_{3k,2}{x}_{3k,2}}{u_{3k}{y}_{3k}}$$

(3)

Because we use the ratio of input to output to calculate the efficiency value, the calculated efficiency value is greater than or equal to 1. The closer the efficiency value is to 1, the higher the efficiency level of the DMU is.

We use the efficiency aggregation method proposed by Chen et al. [41], and the overall efficiency of the system as the weighted sum of the efficiency values of each subsystem in the first stage is defined as:

$$E_{1k}={\omega }_{1k,1}{E}_{1k,1}+{\omega }_{1k,2}{E}_{1k,2}+{\omega }_{1k,3}{E}_{1k,3}$$

(4)

We use the proportion of the total output of each system in the total output of all systems to define the weight of each system, and the sum of the weights of the three subsystems is equal to 1:

$${\omega }_{1k,1}=\frac{u_{1k}{y}_{1k}}{u_{1k}{y}_{1k}+u_{2k}{y}_{2k}+u_{3k}{y}_{3k}}$$

(5)

$${\omega }_{1k,2}=\frac{u_{2k}{y}_{2k}}{u_{1k}{y}_{1k}+u_{2k}{y}_{2k}+u_{3k}{y}_{3k}}$$

(6)

$${\omega }_{1k,3}=\frac{u_{3k}{y}_{3k}}{u_{1k}{y}_{1k}+u_{2k}{y}_{2k}+u_{3k}{y}_{3k}}$$

(7)

Taking Formulas (1)–(3), (5)–(7) into Formula (4), the overall efficiency definition formula of the first stage of DMU_k (k = 1, 2, …, n) can be changed into the following form:

$$E_{1k}=\frac{{\displaystyle \sum _{i=1}^3{\tau }_{1k,i}{z}_{1k,i}}+{\displaystyle \sum _{i=1}^2{\tau }_{2k,i}{z}_{2k,i}}+{\displaystyle \sum _{i=1}^3{\tau }_{3k,i}{z}_{3k,i}}+v_{1k,1}{x}_{1k,1}+v_{2k,1}{x}_{2k,1}+v_{3k,1}{x}_{3k,1}+v_{3k,2}{x}_{3k,2}}{u_{1k}{y}_{1k}+u_{2k}{y}_{2k}+u_{3k}{y}_{3k}}$$

(8)

The efficiency expression of the second stage is:

$$E_{2k}=\frac{{\displaystyle \sum _{i=1}^3{\tau }_{4k,i}{z}_{4k,i}}+v_{4k,1}{x}_{4k,1}+v_{4k,2}{x}_{4k,2}+v_{4k,3}{x}_{4k,3}}{u_{4k}{y}_{4k}}$$

(9)

We still define the overall efficiency of the system as the weighted sum of the efficiency of the two sub stages:

$$E_k={\omega }_{1k}{E}_{1k}+{\omega }_{2k}{E}_{2k}$$

(10)

The weight of each stage is defined as the proportion of the total output of each stage in the total output of the system:

$${\omega }_{1k}=\frac{u_{1k}{y}_{1k}+u_{2k}{y}_{2k}+u_{3k}{y}_{3k}}{u_{1k}{y}_{1k}+u_{2k}{y}_{2k}+u_{3k}{y}_{3k}+u_{4k}{y}_{4k}}$$

(11)

$${\omega }_{2k}=\frac{u_{4k}{y}_{4k}}{u_{1k}{y}_{1k}+u_{2k}{y}_{2k}+u_{3k}{y}_{3k}+u_{4k}{y}_{4k}}$$

(12)

Bringing (8), (9), (11), and (12) into (10), the overall efficiency expression of the system is

$$E_{kk}=\frac{{\displaystyle \sum _{i=1}^3{\tau }_{1k,i}{z}_{1k,i}}+{\displaystyle \sum _{i=1}^2{\tau }_{2k,i}{z}_{2k,i}}+{\displaystyle \sum _{i=1}^3{\tau }_{3k,i}{z}_{3k,i}}+{\displaystyle \sum _{i=1}^3{\tau }_{4k,i}{z}_{4k,i}}+{\displaystyle \sum _{s=1}^4{v}_{sk,1}{x}_{sk,1}}+{\displaystyle \sum _{t=3}^4{v}_{tk,2}{x}_{tk,2}}+v_{4k,3}{x}_{4k,3}}{u_{1k}{y}_{1k}+u_{2k}{y}_{2k}+u_{3k}{y}_{3k}+u_{4k}{y}_{4k}}$$

(13)

The overall efficiency of the system can be obtained by calculating the following linear programming:

$$\min {\displaystyle \sum _{i=1}^3{\tau }_{1k,i}{z}_{1k,i}}+{\displaystyle \sum _{i=1}^2{\tau }_{2k,i}{z}_{2k,i}}+{\displaystyle \sum _{i=1}^3{\tau }_{3k,i}{z}_{3k,i}}+{\displaystyle \sum _{i=1}^3{\tau }_{4k,i}{z}_{4k,i}}+{\displaystyle \sum _{s=1}^4{v}_{sk,1}{x}_{sk,1}}+{\displaystyle \sum _{t=3}^4{v}_{tk,2}{x}_{tk,2}}+v_{4k,3}{x}_{4k,3}$$

$$\mathrm{s}.\mathrm{t}.\left\{\begin{array}{c}{u}_{1k}{y}_{1k}+u_{2k}{y}_{2k}+u_{3k}{y}_{3k}+u_{4k}{y}_{4k}=1\hfill \\ {\displaystyle \sum _{i=1}^3}{\tau }_{1k,i}{z}_{1j,i}+v_{1k,1}{x}_{1j,1}-u_{1k}{y}_{1j}\ge 0,j=1,2,\dots ,n\hfill \\ {\displaystyle \sum _{i=1}^2}{\tau }_{2k,i}{z}_{2k,j}+v_{2k,1}{x}_{2j,1}-u_{2k}{y}_{2j}\ge 0,j=1,2,\dots ,n\hfill \\ {\displaystyle \sum _{i=1}^3}{\tau }_{3k,i}{z}_{3k,j}+v_{3k,1}{x}_{3j,1}+v_{3k,2}{x}_{3j,2}-u_{3k}{y}_{3j}\ge 0,j=1,2,\dots ,n\hfill \\ {\displaystyle \sum _{i=1}^3}{\tau }_{4k,i}{z}_{4k,j}+v_{4k,1}{x}_{4j,1}+v_{4k,2}{x}_{2j,2}+v_{4k,3}{x}_{4j,3}-u_{4k}{y}_{4j}\ge 0,j=1,2,\dots ,n\hfill \\ {\tau }_{1k,i},{\tau }_{2k,i},{\tau }_{3k,i},{\tau }_{4k,i},v_{1k,1},v_{2k,1},v_{3k,1},v_{3k,2},v_{4k,1},v_{4k,2},v_{4k,3},u_{1k},u_{2k},u_{3k},u_{4k}\ge \epsilon \hfill \end{array}\right.$$

(14)

In order to ensure that the calculated result is the effective value and to make each index play a role, we set the lower bound of the optimal solution to an infinitesimal quantity ε, which can be adjusted according to different examples. The calculation model (14) obtains a set of optimal solutions as

$$\begin{array}{l}({\tau }_{1k,1}^{\ast },{\tau }_{1k,2}^{\ast },{\tau }_{1k,3}^{\ast },v_{1k,1}^{\ast },{\tau }_{2k,1}^{\ast },{\tau }_{2k,2}^{\ast },{\tau }_{2k,3}^{\ast },v_{2k,1}^{\ast },{\tau }_{3k,1}^{\ast },{\tau }_{3k,2}^{\ast },{\tau }_{3k,3}^{\ast },\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{v}_{3k,1}^{\ast },v_{3k,2}^{\ast },{\tau }_{4k,1}^{\ast },{\tau }_{4k,2}^{\ast },{\tau }_{4k,3}^{\ast },v_{4k,1}^{\ast },v_{4k,2}^{\ast },v_{4k,3}^{\ast },u_{1k}^{\ast },u_{2k}^{\ast },u_{3k}^{\ast },u_{4k}^{\ast }{)}^{{\rm T}}\end{array}$$

Then, the overall efficiency expression and the efficiency value expression of the three subsystems are respectively:

$$E_{kk}^{\ast }=\frac{{\displaystyle \sum _{i=1}^3{\tau }_{1k,i}^{\ast }{z}_{1k,i}}+{\displaystyle \sum _{i=1}^2{\tau }_{2k,i}^{\ast }{z}_{2k,i}}+{\displaystyle \sum _{i=1}^3{\tau }_{3k,i}^{\ast }{z}_{3k,i}}+{\displaystyle \sum _{i=1}^3{\tau }_{4k,i}^{\ast }{z}_{4k,i}}+{\displaystyle \sum _{s=1}^4{v}_{sk,1}^{\ast }{x}_{sk,1}}+{\displaystyle \sum _{t=3}^4{v}_{tk,2}^{\ast }{x}_{tk,2}}+v_{4k,3}^{\ast }{x}_{4k,3}}{u_{1k}^{\ast }{y}_{1k}+u_{2k}^{\ast }{y}_{2k}+u_{3k}^{\ast }{y}_{3k}+u_{4k}^{\ast }{y}_{4k}}$$

(15)

$$E_{1k,1}^{\ast }=\frac{{\displaystyle \sum _{i=1}^3{\tau }_{1k,i}^{\ast }{z}_{1k,i}}+v_{1k,1}^{\ast }{x}_{1k,1}}{u_{1k}^{\ast }{y}_{1k}}$$

(16)

$$E_{1k,2}^{\ast }=\frac{{\displaystyle \sum _{i=1}^2{\tau }_{2k,i}^{\ast }{z}_{2k,i}}+v_{2k,1}^{\ast }{x}_{2k,1}}{u_{2k}^{\ast }{y}_{2k}}$$

(17)

$$E_{1k,3}^{\ast }=\frac{{\displaystyle \sum _{i=1}^3{\tau }_{3k,i}^{\ast }{z}_{3k,i}}+v_{3k,1}^{\ast }{x}_{3k,1}+v_{3k,2}^{\ast }{x}_{3k,2}}{u_{3k}^{\ast }{y}_{3k}}$$

(18)

$$E_{1k}^{\ast }=\frac{{\displaystyle \sum _{i=1}^3{\tau }_{1k,i}^{\ast }{z}_{1k,i}}+{\displaystyle \sum _{i=1}^2{\tau }_{2k,i}^{\ast }{z}_{2k,i}}+{\displaystyle \sum _{i=1}^3{\tau }_{3k,i}^{\ast }{z}_{3k,i}}+{\displaystyle \sum _{s=1}^3{v}_{sk,1}^{\ast }{x}_{sk,1}}+v_{3k,2}^{\ast }{x}_{3k,2}}{u_{1k}^{\ast }{y}_{1k}+u_{2k}^{\ast }{y}_{2k}+u_{3k}^{\ast }{y}_{3k}}$$

(19)

$$E_{2k}^{\ast }=\frac{{\displaystyle \sum _{i=1}^3{\tau }_{4k,i}^{\ast }{z}_{4k,i}}+v_{4k,1}^{\ast }{x}_{4k,1}+v_{4k,2}^{\ast }{x}_{4k,2}+v_{4k,3}^{\ast }{x}_{4k,3}}{u_{4k}^{\ast }{y}_{4k}}$$

(20)

In view of the situation that multiple solutions will appear in the calculation of this model, we further expand this model with reference to the traditional treatment method of cross-efficiency. For other DMU_j (j = 1, 2, …, n), the optimal solution can also be obtained by using linear programming (14). Based on this, the overall cross-efficiency of DMU_j using the optimal weight of DMU_k in model (14) is defined as:

$$E_{kj}^{\ast }=\frac{{\displaystyle \sum _{i=1}^3{\tau }_{1k,i}^{\ast }{z}_{1j,i}}+{\displaystyle \sum _{i=1}^2{\tau }_{2k,i}^{\ast }{z}_{2j,i}}+{\displaystyle \sum _{i=1}^3{\tau }_{3k,i}^{\ast }{z}_{3j,i}}+{\displaystyle \sum _{i=1}^3{\tau }_{4k,i}^{\ast }{z}_{4j,i}}+{\displaystyle \sum _{s=1}^4{v}_{sk,1}^{\ast }{x}_{sj,1}}+{\displaystyle \sum _{t=3}^4{v}_{tk,2}^{\ast }{x}_{tj,2}}+v_{4k,3}^{\ast }{x}_{4j,3}}{u_{1k}^{\ast }{y}_{1j}+u_{2k}^{\ast }{y}_{2j}+u_{3k}^{\ast }{y}_{3j}+u_{4k}^{\ast }{y}_{4j}}, k,j=1,2,\dots ,n$$

(21)

The above cross-efficiency value determination process is shown in the cross-efficiency matrix E in

Table 3. Cross efficiency matrix (CEM).

DMUk	DMUj
	1	2	3	⋯	n
1	E11	E12	E13	⋯	E1n
2	E21	E22	E23		E2n
3	E31	E32	E33	⋯	E3n
⋯	⋯	⋯	⋯	⋯	⋯
⋯	⋯	⋯	⋯	⋯	⋯
⋯	⋯	⋯	⋯	⋯	⋯
n	En1	En2	En3	⋯	Enn
Average value	${\overline{E}}_1$	${\overline{E}}_2$	${\overline{E}}_{\mathrm{e}}$	⋯	${\overline{E}}_n$

. The element Ekj in the matrix is the cross-efficiency value obtained by DMU_j (j = 1, 2, …, n) using the weight of DMUk, and the element on the diagonal represents the efficiency value during self-assessment.

For DMU_j (j = 1, 2, …, n), the average value of all E_kj (k = 1, 2, …, n), that is, $E_j={\displaystyle {\sum }_{k=1}^n{E}_{kj}}$ (j = 1, 2, …, n), represents the cross-efficiency value of DMU_j. Therefore, the overall cross-efficiency, the cross-efficiency of the two stages, and the cross-efficiency of the three subsystems in the first stage are as follows:

$${\overline{E}}_j=\frac{1}{n}{\displaystyle \sum _{k=1}^n{E}_{kj}^{\ast }}=\frac{1}{n}{\displaystyle \sum _{k=1}^n\frac{{\displaystyle \sum _{i=1}^3{\tau }_{1k,i}^{\ast }{z}_{1j,i}}+{\displaystyle \sum _{i=1}^2{\tau }_{2k,i}^{\ast }{z}_{2j,i}}+{\displaystyle \sum _{i=1}^3{\tau }_{3k,i}^{\ast }{z}_{3j,i}}+{\displaystyle \sum _{i=1}^3{\tau }_{4k,i}^{\ast }{z}_{4j,i}}+{\displaystyle \sum _{s=1}^4{v}_{sk,1}^{\ast }{x}_{sj,1}}+{\displaystyle \sum _{t=3}^4{\pi }_{tk,2}^{\ast }{x}_{tj,2}}+v_{4k,3}^{\ast }{x}_{4j,3}}{u_{1k}^{\ast }{y}_{1j}+u_{2k}^{\ast }{y}_{2j}+u_{3k}^{\ast }{y}_{3j}+u_{4k}^{\ast }{y}_{4j}}}$$

(22)

$${\overline{E}}_{j,1}=\frac{1}{n}{\displaystyle \sum _{k=1}^n{E}_{kj,1}^{\ast }}=\frac{1}{n}{\displaystyle \sum _{k=1}^n\frac{{\displaystyle \sum _{i=1}^3{\tau }_{1k,i}^{\ast }{z}_{1k,i}}+v_{1k,1}^{\ast }{x}_{1k,1}}{u_{1k}^{\ast }{y}_{1k}}}$$

(23)

$${\overline{E}}_{j,2}=\frac{1}{n}{\displaystyle \sum _{k=1}^n{E}_{kj,2}^{\ast }}=\frac{1}{n}{\displaystyle \sum _{k=1}^n\frac{{\displaystyle \sum _{i=1}^2{\tau }_{2k,i}^{\ast }{z}_{2k,i}}+v_{2k,1}^{\ast }{x}_{2k,1}}{u_{2k}^{\ast }{y}_{2k}}}$$

(24)

$${\overline{E}}_{j,3}=\frac{1}{n}{\displaystyle \sum _{k=1}^n{E}_{kj,3}^{\ast }}=\frac{1}{n}{\displaystyle \sum _{k=1}^n\frac{{\displaystyle \sum _{i=1}^3{\tau }_{3k,i}^{\ast }{z}_{3k,i}}+v_{3k,1}^{\ast }{x}_{3k,1}+v_{3k,2}^{\ast }{x}_{3k,2}}{u_{3k}^{\ast }{y}_{3k}}}$$

(25)

$${\overline{E}}_{1j}=\frac{1}{n}{\displaystyle \sum _{k=1}^n{E}_{1k}^{\ast }}=\frac{1}{n}{\displaystyle \sum _{k=1}^n\frac{{\displaystyle \sum _{i=1}^3{\tau }_{1k,i}^{\ast }{z}_{1k,i}}+{\displaystyle \sum _{i=1}^2{\tau }_{2k,i}^{\ast }{z}_{2k,i}}+{\displaystyle \sum _{i=1}^3{\tau }_{3k,i}^{\ast }{z}_{3k,i}}+{\displaystyle \sum _{s=1}^3{v}_{sk,1}^{\ast }{x}_{sk,1}}+v_{3k,2}^{\ast }{x}_{3k,2}}{u_{1k}^{\ast }{y}_{1k}+u_{2k}^{\ast }{y}_{2k}+u_{3k}^{\ast }{y}_{3k}}}$$

(26)

$${\overline{E}}_{2j}=\frac{1}{n}{\displaystyle \sum _{k=1}^n{E}_{2k}^{\ast }}=\frac{1}{n}{\displaystyle \sum _{k=1}^n\frac{{\displaystyle \sum _{i=1}^3{\tau }_{4k,i}^{\ast }{z}_{4k,i}}+v_{4k,1}^{\ast }{x}_{4k,1}+v_{4k,2}^{\ast }{x}_{4k,2}+v_{4k,3}^{\ast }{x}_{4k,3}}{u_{4k}^{\ast }{y}_{4k}}}$$

(27)

3. Case Analysis

3.1. Establishment of Index System

All data in this paper are from the China Statistical Yearbook, the China Energy Statistical Yearbook, the Statistical Yearbook of Chinese Investment Field, the China Population and Employment Statistics Yearbook, the China Water Resources Bulletin, the National Bureau of Statistics, and provincial statistical yearbooks. For some data not published in the statistical data, we consulted the Statistics Bureau or water resources department of relevant provinces. According to the availability of data, the data of China’s 19 provincial administrative regions (in eastern China, Northeast China, and central China) in 2019 were selected for efficiency calculation. Individual missing values were processed by interpolation method. The specific original data values are shown in

Table A2. Data of evaluation index system.

Type	x1k,1	z1k,1	z1k,2	z1k,3	y1k	x2k,1	z2k,1	z2k,2	y2k	x3k,1	x3k,2
Beijing	24.56	54.20	24.6	1.4852	41.7	0.659	153.36	6.9453	691.1	3.7	55.75
Tianjin	8.48	107.46	8.1	0.8051	28.4	0.5779	280.11	2.6935	5106.83	9.2	107
Hebei	23.82	208.48	113.5	3.0811	182.3	3.7712	1578.00	13.7431	6820	114.4	50.04
Shanxi	6.49	78.78	97.3	2.4604	76	4.14512	687.45	11.0051	69,313.12	43.8	308.12
Liaoning	26.02	75.26	256	2.8808	130.3	0.1	439.46	11.1454	4441.1	80.7	302.12
Jilin	17.71	120.44	506.1	1.7157	115.4	5.39	171.99	8.0834	2288.2	81.5	149.68
Heilongjiang	9.57	216.83	1511.4	2.2036	310.4	7.1	380.53	11.6212	9766.93	274.2	658.9
Shanghai	6.47	116.62	48.3	0.964	100.9	49.87	123.21	1.8695	42,674.352	16.9	59.57
Jiangsu	41.36	266.59	231.7	4.122	619.1	80.5	869.88	8.4166	3489.29	619.1	550.24
Zhejiang	30.56	182.23	1321.5	3.3197	165.8	0.7	625.19	6.7434	2937	72.4	421
Anhui	19.32	434.87	539.9	1.985	277.7	4.77541	555.78	7.419	8943.62	150.2	259.48
Fujian	8.86	466.95	1363.9	1.7629	177.5	1.04436	634.48	9.1325	4353.87	83.7	252.57
Jiangxi	12.21	155.24	2051.6	2.0534	253.3	18.09	408.14	6.7125	1320.3	162.5	152.44
Shandong	40.97	350.79	195.2	3.6952	225.3	7.410074	1568.71	22.1594	12,539.1	138.2	599.7
Henan	11.1	299.89	168.6	3.7418	237.8	4.1549	1574.82	17.7199	10,303.95	121.8	572
Hubei	10.11	437.47	613.7	2.5862	303.2	45	678.72	11.1438	4892.2	155.6	452
Hunan	14.27	602.98	2098.3	2.9797	333	38.9	693.39	11.5254	3000	191.7	523.88
Guangdong	72.07	831.13	2068.2	6.4027	412.3	35.7	1338.49	17.437	8377.25	208.5	618.21
Hainan	3.61	50.56	252.3	0.6518	46.4	38.66	161.98	1.3652	494.08	34.2	109.66
Type	z3k,1	z3k,2	z3k,3	y3k	x4k,1	x4k,2	x4k,3	z4k,1	z4k,2	z4k,3	y4k
Beijing	125.95	46.52	42.4	28.76	41.7	691.1	28.76	7868.74	6.5486	1273	164,563
Tianjin	270.06	339.27	58.28	223.25	28.4	5106.83	223.25	12,112.08	6.6307	896.56	90,058
Hebei	1595.11	6469.17	1331.15	3739.24	182.3	6820	3739.24	717,100.30	11.6854	4182.46	46,182
Shanxi	285.50	3126.17	666.7	1361.8	76	69,313.12	1361.8	7094.59	4.474	1902.5	45,549
Liaoning	121.98	3488.73	510	2429.95	130.3	4441.1	2429.95	6716.57	7.8128	2238.4	57,067
Jilin	442.85	5644.93	466.2	3877.93	115.4	2288.2	3877.93	11,285.40	2.0809	1456.43	43,475
Heilongjiang	990.52	14,338.1	564.1	7503.01	310.4	9766.93	7503.01	8837.84	2.3123	1776.9	36,001
Shanghai	3.51	117.37	40.8	95.89	100.9	42,674.352	95.89	8012.22	11.3425	1376.2	156,587
Jiangsu	383.94	5381.48	734.51	3706.2	619.1	3489.29	3706.2	58,766.89	69.3442	4745.2	122,398
Zhejiang	99.38	977.44	406.83	592.15	165.8	2937	592.15	36,702.88	57.4571	3875.1	107,814
Anhui	761.48	7287	1346.9	4054	277.7	8943.62	4054	35,631.90	18.358	4384	58,072
Fujian	1327.45	822.43	548.85	493.9	177.5	4353.87	493.9	31,164.02	18.0365	2781.26	106,966
Jiangxi	471.50	3665.14	700.8	2157.45	253.3	1320.3	2157.45	26,794.15	12.2207	2631.95	52,865
Shandong	996.96	8312.81	1445.9	5357	225.3	12539.1	5357	591,438.30	30.4172	5987.9	70,129
Henan	2460.42	10,734.54	2277	6695.36	237.8	10,303.95	6695.36	51,241.12	20.6775	6562	55,825
Hubei	1052.38	4608.6	1164	2724.98	303.2	4892.2	2724.98	39,128.68	17.5724	3548	76,712
Hunan	2114.48	4616.38	1409.24	2974.84	333	3000	2974.84	37,941.46	16.088	3666.48	57,746
Guangdong	321.51	2160.64	823	1240.8	412.3	8377.25	1240.8	46,093.28	83.8891	7150.25	94,448
Hainan	36.84	272.65	175.15	144.96	46.4	494.08	144.96	3277.63	0.2832	102.2	56,740

in the Appendix A.

For the indicators of the water resources subsystem, since the total energy consumption data of water production and supply industry in many provinces cannot be obtained, we use the power consumption of this industry to replace it. In the energy production subsystem, according to the availability of data, we use the relevant indicators of the power and heat production and supply industries to reflect the energy production situation. The water consumption of the power and heat production and supply industries in a few provinces is replaced by the water consumption data of fire (nuclear) power. We choose R&D practitioners as the resources needed for the economic development of each province.

3.2. Efficiency Values and Ranking

We use MATLAB 2018a software to calculate the input–output efficiency of the WEF-Eco systems of 19 provincial administrative regions in 2019. In the calculation process of this example, on the premise of ensuring that the values of expressions (22)–(27) are non-null values, after testing, we determine that the minimum lower bound ε of the solution of model (14) should be 10⁻⁴. Before bringing the data into the model calculation, we use the following data standardization methods to process the data:

$$x_i^{\ast }=\frac{x_i-\underset{1\le j\le n}{\min }\left\{x_j\right\}}{\underset{1\le j\le n}{\max }\left\{x_j\right\}-\underset{1\le j\le n}{\min }\left\{x_j\right\}}$$

(28)

The results are shown in

Table 4. WEF-Eco system efficiency values.

Type	${\overline{E}}_j$	Ranking	${\overline{E}}_{1j}$	Ranking	${\overline{E}}_{2j}$	Ranking
Beijing	24.06	16	211.45	17	1.00	1
Tianjin	1264.69	18	3187.09	19	2.36	3
Hebei	21.09	15	2.67	5	220.11	18
Shanxi	13.90	10	2.67	6	42.81	13
Liaoning	8.19	4	2.71	8	20.76	7
Jilin	19.92	14	3.01	9	57.63	15
Heilongjiang	9897.34	19	1.81	2	707,919,722.99	19
Shanghai	2.10	1	1.72	1	2.82	4
Jiangsu	8.27	5	3.70	11	21.15	8
Zhejiang	5.12	2	4.08	12	12.68	6
Anhui	12.99	9	2.08	3	43.70	14
Fujian	5.17	3	4.51	15	7.58	5
Jiangxi	15.22	12	4.96	16	39.08	11
Shandong	10.14	8	2.68	7	69.21	17
Henan	17.29	13	2.41	4	61.81	16
Hubei	8.32	6	3.19	10	21.40	9
Hunan	15.07	11	4.36	13	41.78	12
Guangdong	8.73	7	4.38	14	28.02	10
Hainan	931.98	17	2676.13	18	1.38	2

and

Table 5. Efficiency values of three subsystems in the first stage.

Type	${\overline{E}}_{j,1}$	Ranking	${\overline{E}}_{j,2}$	Ranking	${\overline{E}}_{j,3}$	Ranking
Beijing	17.12	18	111.73	18	1,766,195.52	19
Tianjin	13,677,286.78	19	2.83	3	2.95	14
Hebei	2.81	10	18.89	11	1.84	7
Shanxi	3.76	15	1.00	1	2.51	12
Liaoning	3.73	14	12.38	6	1.39	2
Jilin	3.32	11	20.67	12	1.09	1
Heilongjiang	1.57	3	6.78	5	1.48	3
Shanghai	1.20	2	2.60	2	1.96	9
Jiangsu	1.17	1	78.02	16	2.30	10
Zhejiang	4.68	17	17.86	10	4.05	16
Anhui	2.01	6	6.10	4	1.77	5
Fujian	3.53	12	13.76	9	6.98	18
Jiangxi	2.17	8	85.85	17	1.95	8
Shandong	3.56	13	12.96	7	1.62	4
Henan	2.09	7	13.66	8	1.82	6
Hubei	1.75	5	35.77	14	2.38	11
Hunan	2.71	9	58.31	15	2.76	13
Guangdong	4.29	16	24.12	13	4.00	15
Hainan	1.74	4	123,927,931.15	19	5.47	17

.

4. Discussion

In order to better reflect the regional characteristics of efficiency, we put the provinces belonging to Northeast China, Eastern China, and Central China together. According to the efficiency value ranking, the efficiency value ranking is expressed in a gradual color from yellow to red. The redder the color of the color block, the higher the efficiency level ranking. The yellower the color of the color block, the lower the efficiency level ranking. The results are shown in Figure 2.

From Figure 2, we can intuitively see the regional characteristics and the differences in efficiency in each province. In terms of overall efficiency, the efficiency level of the eastern and central regions is high in the middle and low in the north and south ends. The efficiency level of the northeast region is gradually decreasing from south to north. The provinces with high efficiency levels include Shanghai, Zhejiang, and Fujian. The provinces with low efficiency mainly include Hainan, Tianjin, and Heilongjiang. In terms of the efficiency of the first stage system (WEF system), the regional characteristics are not obvious. The efficiency levels in the three regions are quite different. The efficiency levels of Zhejiang, Heilongjiang, and Anhui are relatively high, while those of Beijing, Hainan, and Tianjin are relatively low. In terms of the system efficiency in the second stage, the regional characteristics are relatively obvious. The provinces with higher efficiency levels are mostly located in the eastern region, while the efficiency levels of other regions are relatively low. Similarly, we can find that the efficiency level of northern provinces is low, while that of southern provinces is high. These characteristics are roughly consistent with the level of economic development in real life. In terms of water resources subsystem, the regional characteristics are not obvious, and the efficiency levels in each region differ greatly. The provinces with high efficiency level mainly include Jiangsu, Shanghai, and Heilongjiang, while the provinces with low efficiency level mainly include Zhejiang, Beijing, and Tianjin. As far as the energy subsystem is concerned, the regional characteristics of efficiency level are not obvious. There are great differences in efficiency level among the three regions. The efficiency value of Shanxi is 1, and the efficiency level is the highest. In addition, the efficiency level of Shanghai and Tianjin is relatively high. The provinces with low efficiency level mainly include Jiangxi, Beijing, and Hainan. As far as the food subsystem is concerned, the efficiency level of Northeast China is the highest, and the efficiency level of central and eastern regions is low. The efficiency value of Beijing is very large, which indicates that the efficiency level is very low, and special attention should be paid to it.

Provinces with higher efficiency levels often have higher management levels, more advanced technology levels, etc., which can make full use of resources for production activities, effectively reduce resource waste, or obtain more output. Through Figure 2, we can clearly find out the reasons for the high and low efficiency levels of systems in different provinces. For example, the overall efficiency level of the system in Shanghai is relatively high. We find that the efficiency level of the first stage (WEF system) and the efficiency level of the second stage in Shanghai are relatively high. The high efficiency level of the first stage is due to the relatively high efficiency of the three subsystems in the WEF system. This kind of province with a high efficiency level in each part tends to develop more comprehensively and is ahead of other provinces in all aspects. However, there may also be situations where the overall efficiency level is high, but some local efficiency levels are low, such as in Zhejiang Province and Fujian Province. The efficiency levels of these two provinces in the first stage are relatively low, but their efficiency levels in the second stage are relatively high, which leads to higher overall efficiency. The main reason for this is that the efficiency level of the second stage has obvious advantages, making up for the impact of the first stage, which has a low efficiency level. The reasons for the efficiency level of other provinces can also be seen from Figure 2. On the contrary, we can also give an improvement method to improve the system efficiency according to the results in Figure 2. For example, the overall efficiency level of Beijing is low. We can find that the efficiency level of the first phase of Beijing is low. To be more precise, the efficiency level of the three subsystems in the first phase is relatively low. We can reduce the waste of resources or increase the output to improve the efficiency level of the three subsystems by introducing advanced management experience and improving the level of science and technology, so as to improve the overall efficiency level of Beijing. Other provinces with low overall efficiency can also accurately find the links to be improved according to the results in Figure 2.

5. Conclusions

This paper establishes a three-dimensional network structure to describe the WEF-Eco system and establishes the corresponding complex network DEA model taking the data of 19 regions in China in 2019 as an example to show the application effect of this model. The innovation of this article is to expand the traditional two-dimensional network used to describe the structure of the WEF system into a three-dimensional network structure which can more accurately describe the structure of the WEF system. We build weights and corresponding models according to the output indicators, and expand the types of DEA models.

In this paper, the structure of the WEF system is shown more clearly, and the calculation method of the overall and local efficiency of the system is proposed, so that we can find the reason for the low efficiency level more accurately and then propose effective improvement methods. In addition, the construction method of the three-dimensional network DEA model proposed in this paper can inspire us to deal with other problems with three-dimensional network structure and even build more complex structures and DEA models. It is also worth noting that in the process of building the additive DEA model, we should not be limited to the traditional weight building method, but should reasonably select input indicators or output indicators to build the weight of each system according to the actual situation, so that the evaluation results are more consistent with the actual situation.

There are still some limitations in this paper. First, in addition to indicators of water resources, energy, and food, we increase the investigation of investment, resources, and labor indicators and do not distinguish the preference relationship between indicators. Second, this paper does not consider the unexpected output, such as the emission of pollutants in the production process. Finally, this paper only considers the efficiency value of a single period, and the dynamic model of multiple periods is not given. All of these are challenges for the future research.