Robust Planning of Energy and Environment Systems through Introducing Traffic Sector with Cost Minimization and Emissions Abatement under Multiple Uncertainties

Abstract

Motor vehicles have been identified as a growing contributor to air pollution, such that analyzing the traffic policies on energy and environment systems (EES) has become a main concern for governments. This study developed a dual robust stochastic fuzzy optimization—energy and environmental systems (DRSFO-EES) model for sustainable planning EES, while considering the traffic sector through integrating two-stage stochastic programming, robust two-stage stochastic optimization, fuzzy possibilistic programming, and robust fuzzy possibilistic programming methods into a framework, which can be used to effectively tackle fuzzy and stochastic uncertainties as well as their combinations, capture the associated risks from fuzzy and stochastic uncertainties, and thoroughly analyze the trade-offs between system costs and reliability. The proposed model can: (i) generate robust optimized solutions for energy allocation, coking processing, oil refining, heat processing, electricity generation, electricity power expansion, electricity importation, energy production, as well as emission mitigation under multiple uncertainties; (ii) explore the impacts of different vehicle policies on vehicular emission mitigation; (iii) identify the study of regional atmospheric pollution contributions of different energy activities. The proposed DRSFO-EES model was applied to the EES of the Beijing-Tianjin-Hebei (BTH) region in China. Results generated from the proposed model disclose that: (i) limitation of the number of light-duty passenger vehicles and heavy-duty trucks can effectively reduce vehicular emissions; (ii) an electric cars’ policy is enhanced by increasing the ratio of its power generated from renewable sources; and (iii) the air-pollutant emissions in the BTH region are expected to peak around 2030, because the energy mix of the study region would be transformed from one dominated by coal to one with a cleaner pattern. The DRSFO-EES model can not only provide scientific support for the sustainable managing of EES by cost-effective ways, but also analyze the desired policies for mitigating pollutant emissions impacts with a risk adverse attitude under multiple uncertainties.

1. Introduction

1.1. Background

China has experienced severe environmental pollution in recent years, which pose a critical threat to public health and sustainable development [1,2]. Energy-related activities are the dominant sources of air pollution [3], with the amounts of carbon dioxide (CO₂) and air pollutant emitted from electricity generation plants accounting for approximately 40% and 30% of the total CO₂ and air pollutant emissions, respectively [2,4]. Further, motor vehicles have been identified as growing contributors to air pollution due to the rapid growth of transportation, accounting for approximately 20–67% of carbon monoxide (CO) emissions, 12–36% of oxynitride (NOx) emissions, and 12–39% of hydrocarbon compound (HC) emissions [5,6]. The scale of emissions of most of China’s regions has exceeded the capacity of self-purification and air-pollutants’ diffusion from the atmosphere. There is currently a severe conflict between increasing energy demand, excessive vehicle population, and “high coal” energy mix on the one hand, and the imperative of mitigating air pollution on the other hand [7]. To tackle the above-mentioned problems, several policies and measures have been implemented with regard to the development of renewable energy resources: adjustment of the energy structure; encouragement of the use of electric cars (EVs); improvement of energy conversion efficiencies; and enhancement of vehicular emission standards. However, it remains unclear how much pollution reduction can be achieved by these control measures and policies. This situation has forced local managers to propose ambitious schemes for planning energy and environment systems (EES), and to deeply analyze the impacts of different emission mitigation policies and measures on these EES [8,9]. However, EES are complicated by many systemic uncertainties regarding the relevant environmental, economic, energy, and social factors. For instance, electricity demands are often shown as stochastic uncertainties that vary over time based on extant policies and highly variable conditions [10]. Moreover, many economic data and energy demands often exhibit ambiguity [11,12]. Failure to consider these uncertainties may result in less robust decision support [8,13]. Therefore, EES planning as well as considering uncertainty information are required to help confront such problems of EES, and to ensure sustainable economic development and environmental protection [14,15].

1.2. Literature Review

Numerous non-deterministic programming approaches have been used to handle uncertainties in EES.

Table 1. Previous studies related to the subject.

Ref. No.	Non-Deterministic Programming				Research Area		Considering Traffic Sector
	TSP	RTSO	FPP	RFPP	Energy Systems	Others	Yes	No
[5]	✓				✓			✓
[17]	✓				✓			✓
[18]	✓				✓		✓
[19]	✓				✓		✓
[20]	✓				✓			✓
[21]		✓			✓
[22]	✓	✓			✓			✓
[23]			✓		✓			✓
[24]	✓	✓				✓		✓
[25]			✓	✓	✓			✓
[26]		✓				✓		✓
[27]			✓	✓		✓		✓
[28]	✓		✓			✓		✓
[29]			✓			✓		✓
[30]			✓			✓		✓
[31]			✓		✓			✓
[32]			✓	✓		✓		✓

Note: TSP, two-stage stochastic programming; RTSO. robust two-stage stochastic optimization; FPP, fuzzy possibilistic programming; and RFPP, robust fuzzy possibilistic programming.

lists some previous studies on non-deterministic programming problems.

Among them, two-stage stochastic programming (TSP) has been widely used to tackle uncertainties expressed as a probability distribution [3,16,17,18]. For instance, Gong et al. [16] proposed a two-stage programming method to optimize electric power systems considering air pollutant emissions and CO₂ mitigation. Mavromatidis et al. [19] developed a two-stage integer linear program model to optimize distributed energy systems, enable cost-optimized design decisions regarding technology selection and sizing before the determination of uncertain parameters. Mohan et al. [20] presented a two-stage stochastic method for managing the energy reserve of a microgrid system, with an emphasis on the different levels and sources of uncertainties. However, the TSP is unable to regard the variability of stochastic recourse values because it is based on the assumption that the manager adopts a risk-neutral attitude. Thus, TSP may become infeasible when managers are risk-averse under the conditions of high variability [21].

The robust two-stage stochastic optimization (RTSO) method is an attractive method for tackling the above shortcomings of TSP. It is specifically used to penalize the costs of the second-stage that are greater than the expected values and capture the associated risk of stochastic uncertainties [22,23,24]. In the last few decades, the RTSO method has been extensively employed in many research areas, such as supply chain systems, electric power systems, solid waste management, and water resource allocation [3,23,25]. For example, Govindan and Cheng [26] developed a stochastic robust programming method for improving retail supply chain planning through supply chain coordination, risk reduction, vendor selection, and sustainability assessment. Xu et al. [24] proposed a robust TSP method for tackling water resource allocation problems, enabling the handling of uncertainties expressed as stochastic, and analysis of policy scenarios regarding economic penalties for the violation of predefined policies.

However, in the real world, many economic parameters and energy demands often exhibit ambiguities, which can be shown as fuzzy sets [27,28]. Fuzzy possibilistic programming (FPP) theory can effectively address the fuzzy uncertainties of goals and constraints [29]. For instance, Vahdani et al. [30] employed an FPP method for closed-loop recycling collection networks, in which some uncertainty information (e.g., distance, capacity, demand, costs, as well as returned products quantity) were tackled by FPP. Lu et al. [31] proposed an interval FPP method for managing China’s energy systems with CO₂ emissions constraints, which could address the uncertainties presented in terms of fuzzy-boundary intervals in both the objective and constraints.

However, an FPP algorithm is unable to ensure the minimization of objective function under all conditions because minimizing of the expected objective value is used as the objective function. This can result in significant deviations of the optimized decision schemes, even in the event of system optimization failure. Robust fuzzy possibilistic programming (RFPP) was developed by Pishvaee et al. [32] to overcome the drawbacks of FPP methods and involves the extension of robust optimization from stochastic algorithms to fuzzy algorithms. RFPP considers three sections in objective function: (i) the minimization of the weight sum of the expected objective values; (ii) the difference between two possible extreme objective values; and (iii) the penalty for constraint violation as the objective function [27,33]. It has, however, been limitedly applied to EES planning.

Generally, although previous research works can effectively deal with EES issues under multiple uncertainties, several gaps still need to be remedied. Firstly, few of the previous studies considered the traffic sector, which has been identified as a growing contributor to air pollution. Currently, a series of traffic policies have been adopted to alleviate the air pollution caused by motor vehicles. However, it remains unclear how much pollution reduction can be achieved by these control measures and policies. Secondly, most of these studies are incapable of considering the system risks from the stochastic and fuzzy uncertainties during the optimization process, which may lead to significant deviations in the optimized decision schemes, even in the event of system optimization failure. Thirdly, the RFPP method is commonly used to plan water resource allocation, and solid management, and is scarcely applied to plan EES systems.

1.3. Objective

The objective of this study was the development of a dual robust stochastic fuzzy optimization—energy and environmental systems (DRSFO-EES) model for planning EES while considering the traffic sector. This study is the first attempt at planning an EES while considering the traffic sector by integrating TSP, FPP, RTSO, and RFPP into a single framework. The proposed model can effectively tackle stochastic and fuzzy uncertainties as well as their combinations, capture associated risks from fuzzy and stochastic uncertainties, and thoroughly analyze trade-offs between system costs and reliability. The proposed DRSFO-EES model was applied to the Beijing-Tianjin-Hebei (BTH) region in China, which experiences severe smog and haze associated with high concentrations of air pollutants. The following is a detailed enumeration of the capabilities of the proposed model: (i) exploration of the impacts of different vehicle policies (such as regarding EVs usage, EV’s power source, and vehicular emission standards) on vehicular emission mitigation via scenario analysis; (ii) generation of robust optimized solutions for energy allocation, coking processing, oil refining, heat processing, electricity generation, electricity power expansion, electricity importation, energy production, as well as emission mitigation under multiple uncertainties; (iii) identification of the study regional atmospheric pollution contributions of different energy activities such as coke processing, heat processing, oil refining, electricity generation, and motor vehicle operation.

2. Methodology

The robust two-stage stochastic optimization (RTSO) method brings risk aversion into stochastic programming methods, and finds robust schemes for system management [33,34]. The RTSO method can deal with the stochastic uncertainties of real-world management problems, analyze the policy scenarios associated with economic penalties when the predefined policies of the first-stage are violated, capture the variability of the second-stage costs that are greater than the expected values, and evaluate trade-offs between system economy and risk [23,24]. An RTSO method is formulated as follows:

$$\mathrm{Min}{f}^{}=C_{T_1}^{}{X}^{}+{\displaystyle \sum _{h=1}^s{p}_h{D}_{T_2}^{}}{Y}_h+{\rho }_1{\displaystyle \sum _{h=1}^s{p}_h V^{}}$$

(1)

subject to

$$A_{r}X\le B_r,r\in M; M=1,2,\dots , m_1$$

(2)

$$A_{r}X+{A^{\prime }}_{rh}{Y}_h\ge {w^{\prime }}_{ih}, i\in M; i= 1,2, \dots , m_2; h=1, 2, \dots , s$$

(3)

$$x_j\ge 0, x_j\in X;j=1, 2, \dots , n_1$$

(4)

$$y_{jh}\ge 0, y_{jh}\in Y;j=1, 2, \dots , n_2; h=1, 2,\dots , v$$

(5)

$$V^{}=D_{T_2}^{}{Y}_h-{\displaystyle \sum _{h=1}^s{p}_h{D}_{T_2}^{}{Y}_h}+2{\theta }_h, h=1,2, \dots ,s$$

(6)

$$v_{jh}\ge 0,v_{jh}\in V; j=1,2, \dots ,n_2; h=1,2, \dots ,s$$

(7)

where, $X$ and $Y_h$ denote the decision variables of the first-stage and the second-stage, respectively; $p_h$ are occurrence probability of scenario h, ${\displaystyle \sum _{h=1}^H{p}_h}=1$; $C$ represent coefficients of X and $D_h$ are coefficients of $Y_h$; ${w^{\prime }}_{ih}$ are random variables with probability levels $p_h$; $A_r$ is the fixed coefficient of X; $A_{rh}{}^{\prime }$ are coefficients of $Y_h$. $B_r$ represent the boundary vectors of the right-hand side of constraints. ${\theta }_h$ denote slack variables used for achieving looser constraints. ρ₁ denotes a goal programming weight of stochastic uncertainties; the managers can regulate the variability of the stochastic recourse cost through adjusting the ρ₁ level [24,35]. When ρ₁= 0, the RTSO model becomes a conventional TSP, which indicates that the managers adopt risk-neutral attitudes and the variability of the stochastic uncertain recourse costs is not considered. However, when ρ₁ = 1, the managers adopt risk- aversive attitudes and the variability of the second-stage cost is considered. $V$ denotes the deviation of an expected value from the given scenario’s cost [3]. Besides, constraints (6) and (7) can define the positive variability of the recourse costs.

However, RTSO is inefficient in addressing the uncertainties expressed by fuzzy sets. Thus, fuzzy possibilistic programming (FPP) is joined to RTSO as a hybrid robust stochastic-fuzzy optimization (RSFO) model as follows:

$$\mathrm{Min} {\tilde{f}}^{}={\tilde{C}}_{T_1}^{}{X}^{}+{\displaystyle \sum _{h=1}^s{p}_h{\tilde{D}}_{T_2}^{}}{Y}_h+{\rho }_1{\displaystyle \sum _{h=1}^s{p}_h V^{}}$$

(8)

subject to

$$Cr\{\tilde{B}|A_{r}X\le {\tilde{B}}_r\}\ge \alpha ,r\in M; M=1,2,\dots , m_1$$

(9)

$$A_{r}X+{A^{\prime }}_{rh}{Y}_h\ge {w^{\prime }}_{ih}, i\in M; i= 1,2, \dots , m_2; h=1, 2, \dots , s$$

(10)

$$x_j\ge 0, x_j\in X;j=1, 2, \dots , n_1$$

(11)

$$y_{jh}\ge 0, y_{jh}\in Y;j=1, 2, \dots , n_2; h=1, 2,\dots , v$$

(12)

$$V^{}=D_{T_2}^{}{Y}_h-{\displaystyle \sum _{h=1}^s{p}_h{D}_{T_2}^{}{Y}_h}+2{\theta }_h, h=1,2, \dots ,s$$

(13)

$$v_{jh}\ge 0,v_{jh}\in V; j=1,2, \dots ,n_2; h=1,2, \dots ,s$$

(14)

where $Cr\{·\}$ represents the credibility measure of a fuzzy event in $\{·\}$; ${\tilde{C}}_{T_1}^{}$ and ${\tilde{D}}_{T_2}^{}$ are cost coefficients expressed as a triangular fuzzy number; ${\tilde{B}}_r(B_r{}^1,B_r{}^2,B_r{}^3)$ is the boundary vectors of the right-hand side of constraints, which expresses as triangular fuzzy sets with its membership functions $\mu ({\tilde{B}}_r)$. α denotes the predetermined confidence-level. $Cr\{\tilde{B}|A_{r}X\le {\tilde{B}}_r\}\ge \alpha$ denotes the credibility of satisfying $A_{r}X\le {\tilde{B}}_r$ is higher than or equal to confidence-level α. According to reference [29], the detailed solution procedures of model (2) can be summarized as: firstly, transforming objective function (8) into its expected value form; secondly, converting constraints (9) into their crisp equivalents. A series of solutions can be obtained under difference confidence-levels [25].

However, an FPP algorithm cannot ensure the minimization of objective function under all conditions because the expected objective value is used as the objective function. This may result in significant deviations of the optimized decision solutions, even in the event of system optimization failure. Robust fuzzy possibilistic programming (RFPP) was developed by Pishvaee et al. [32] to overcome the above-mentioned drawbacks of FPP methods. Additionally, the RFPP method is scarcely applied to EES systems, which are commonly used to plan water resource allocation, and solid management [25,27].

Therefore, this study developed a dual robust stochastic fuzzy optimization (DRSFO) method through integrating TSP, RTSP, FPP, and RFPP methods in a single framework, which can be used to effectively tackle fuzzy and stochastic uncertainties as well as their combinations, capture the associated risks from fuzzy and stochastic uncertainties, as well as thoroughly analyze the trade-offs between system costs and reliability. In detail, it is formulated as follows:

$$\begin{array}{cc}\hfill \mathrm{Min} {\tilde{f}}^{}& =\overline{f}+{\lambda }_1·(f_{\max }-\overline{f})+{\lambda }_2·(\overline{f}-f_{\min })+{\rho }_2·[B^R(\alpha )-B_r{}^1]\hfill \\ & =0.25·(C_{T_1}^1+2·C_{T_1}^2+C_{T_1}^3)·X+{\lambda }_1·[C_{T_1}^3-0.25·(C_{T_1}^1+2·C_{T_1}^2+C_{T_1}^3)]·X\hfill \\ & +{\lambda }_2·[0.25·(C_{T_1}^1+2·C_{T_1}^2+C_{T_1}^3)-C_{T_1}^1]·X\hfill \\ & +0.25·(D_{T_2}^1+2·D_{T_2}^2+D_{T_2}^3)·{\displaystyle \sum _{h=1}^s{p}_h}·Y_h\hfill \\ & +{\lambda }_1·[D_{T_2}^3-0.25·(D_{T_2}^1+2·D_{T_2}^2+D_{T_2}^3)]·{\displaystyle \sum _{h=1}^s{p}_h}·Y_h\hfill \\ & +{\lambda }_2·[0.25·(D_{T_2}^1+2·D_{T_2}^2+D_{T_2}^3)-D_{T_2}^1]·{\displaystyle \sum _{h=1}^s{p}_h}·Y_h\hfill \\ & +{\rho }_1·{\displaystyle \sum _{h=1}^s{p}_h V+} {\rho }_2·[B^R(\alpha )-B_r{}^1]\hfill \end{array}$$

(15)

subject to:

$$Cr\{\tilde{B}|A_{r}X\le {\tilde{B}}_r\}\ge \alpha ,r\in M; M=1,2,\dots , m_1$$

(16)

$$A_{r}X+{A^{\prime }}_{rh}{Y}_h\ge {w^{\prime }}_{ih}, i\in M; i= 1,2, \dots , m_2; h=1, 2, \dots , s$$

(17)

$$x_j\ge 0, x_j\in X;j=1, 2, \dots , n_1$$

(18)

$$y_{jh}\ge 0, y_{jh}\in Y;j=1, 2, \dots , n_2; h=1, 2,\dots , v$$

(19)

$$V^{}=D_{T_2}^{}{Y}_h-{\displaystyle \sum _{h=1}^s{p}_h{D}_{T_2}^{}{Y}^{}}+2{\theta }_h, h=1,2, \dots ,s$$

(20)

$$v_{jh}\ge 0,v_{jh}\in V; j=1,2, \dots ,n_2; h=1,2, \dots ,s$$

(21)

where $\overline{f}$ denotes the expected value of objective function (8); $f_{\max }$ and $f_{\min }$ are the maximum value and minimum value of objective function (8); $B^R(\alpha )$ is the maximum values of all potential values, i.e., $B^R(\alpha )=\sup \{B|B={\mu }^{-1}(\alpha )\}$, ${\mu }^{-1}$ is the inverse of $\mu$. The first section in the objective function (i.e., $\overline{f}$) is used to denote the system cost, while minimizing the weight sum among the expected objective values; the second section [i.e., ${\lambda }_1· (f_{\max }-\overline{f}) + {\lambda }_2· (\overline{f}-f_{\min })$] represents the differences between two extreme possible values, which can enhance the robustness of optimization solutions, where ${\lambda }_1$ and ${\lambda }_2$ represent their respective weight coefficients; and the last section (i.e., ${\rho }_2·[B^R(\alpha )-B_r{}^1]$) denotes the difference between the extreme values of fuzzy uncertainty information, which can control the feasibility robustness of the generated solution. ρ₂ is a goal programming weight for fuzzy uncertainties [25]. The evaluation of the trade-offs between system cost and risk can be obtained under different coefficients ${\lambda }_1$, ${\lambda }_2$, ρ₁, ρ₂, and α. The constraint (16) can be tackled by a conventional FPP algorithm.

3. Applications

3.1. Statement of Problem

The BTH region is regarded as China’s economic, political, cultural, and technological innovation center, representing the most dynamic urban cluster in the country, including 11 prefecture-level cities of Hebei Province, as well as two megacities (Beijing and Tianjin) [36,37]. The region has experienced rapid economic development, with a high annual growth rate of 6.64% for gross domestic product (GDP) over the last five years [1,2,3]. The BTH region’s energy demands have dramatically increased along with its rapid urbanization, industrialization, and economic development. According to the National Bureau of Statistics of the People’s Republic of China, the region’s electricity demand reached 50.92 million kWh in 2016, representing an increase of 13.31% compared with five years earlier. The domestic energy structure is heavily dependent on coal resources, with coal-based power accounting for 85.04% of the region’s electricity power generation, while renewable energy contributes only 8%.

Fossil fuels in general are the major contributor to air pollution. For example, approximately 50%, 70%, 90%, and 65% of the total CO₂ emissions, nitrogen oxide (NO_x), sulfur dioxide (SO₂), and particulate matter (PM) are contributed by coal resources [4,38]. In addition, the BTH region is one of the high-traffic regions of China, with a vehicle population of 20.67 million in 2016, representing an increase of 63% from 2011. Motor vehicles are recognized as a major contributor of regulated pollutants such as CO, NO_x, HC, and PM. Based on several references, vehicular emissions contributed approximately 20–67% of CO emission, 12–36% of NO_x emission, and 12–39% of HC emissions [5].

The BTH region experiences severe air pollution due to increasing energy consumption and vehicle population, as well as a “high coal” energy mix characteristic. Cities in the region generally occupy more than half of the top-ten spots for the most polluted cities in China [39]. The scale of emissions in the BTH region has exceeded the capacity of self-purification and air-pollutants’ diffusion from the atmosphere [25,33]. This situation has forced local managers to propose ambitious schemes for planning the EES of the BTH region.

To tackle the above-mentioned problems, several policies and measures have been implemented with regard to the development of renewable energy resources: adjustment of the energy structure; encouragement of the use of EVs; improvement of energy conversion efficiencies; and enhancement of vehicular emission standards. For example, several strategies have been proposed for reducing vehicular emissions, such as the deployment of EVs and the implementation of the China VI vehicular emission standard. However, it remains unclear how much pollution reduction can be achieved by these control measures and policies. As noted earlier, this paper proposed the DRSFO-EES model for the assessment of the impacts of different emission mitigation policies and measures on EES; development of robust optimization solutions for EES planning; determination of the regional atmospheric pollution contributions of different energy sectors such as energy processing, electricity generation, and vehicular traffic.

3.2. Schematic Overview of This Study

The proposed DRSFO-EES model was applied to EES planning for the BTH region. As shown in Figure 1, the considered EES included several energy-related activities such as heat processing, coke processing, oil refining, electricity conversion, EVs expansion, and energy import. The electricity demands of the BTH region are supplied by local power plants and import from adjacent power grids. The major local electricity conversion technologies include natural gas-based, coal-based, wind and solar power. The vehicle population consists of heavy-duty passenger vehicles (HDVs), light-duty passenger vehicles (LDVs), light duty trucks (LDTs), heavy duty trucks (HDTs), EVs, and “others”. The DRSFO-EES model was used in this study to determine the air pollutants and greenhouse gases (NOx, SO₂, PM, CO, HC, and CO₂) emitted from the different energy related-activities in the region. The EES considered a number of uncertainties (such as electricity demand, and many economic parameters expressed as random variables and fuzzy sets) that would affect the optimization scheme.

3.3. Development of DRSFO-EES Model

Minimizing the system cost is the objective of DRSFO-EES, which includes costs for coke processing, heat processing, input for coal-based power, oil refining, input for natural gas-based power, imported electricity, first-stage of electricity generation, second-stage of electricity generation, electricity expansion, energy production, subsidy for solar power generation, subsidy for wind power generation, EVs charging piles, EVs charging stations, pollutants treatment, and risk recourse for the stochastic and fuzzy uncertainties. The proposed DRSFO-EES can be solved by Lingo 11.0 software. The schematic diagram of DRAOM-EWNS model is presented in Figure 2.

In the DRSFO-EES, t denotes the planning periods, period 1 (2020), period 2 (2025), and period 3 (2030); i denotes the primary energy production type, i = 1 is coal, i = 2 is natural gas, and i = 3 is crude oil; k denotes the electric conversion technology, k = 1 is coal-based power, k = 2 is natural gas-based power, k = 3 is solar power, and k = 4 is wind power; h expresses the electricity demand-level, h = 1 (low demand level), h = 2 (medium demand level), h = 3 (high demand level), ${\displaystyle \sum _{h=1}^3{p}_h}=1$; g are the vehicles types, g = 1 is HDVs, g = 2 is LDVs, g = 3 is TDTs, g = 4 is HDTs, and g = 5 is others. The proposed DRSFO-EES is formulated as follows:

Objective:

$$\begin{array}{cc}\hfill \mathrm{Min} \tilde{f}& =\overline{f}+{\lambda }_1(f_{\max }-\overline{f})+{\lambda }_2(\overline{f}-f_{\min })+V{F}_t\hfill \\ & =f_1 + f_2 + f_3 + f_4+ f_5 + f_6 + f_7\hfill \\ & + f_8+ f_9 + f_{10} - f_{11}^{}-f_{12}+f_{13} + f_{14} \hfill \\ & + f_{15}+ f_{16} \hfill \end{array}$$

(22)

where $\tilde{f}$ represents the system cost; ${\tilde{f}}_1$ denotes the heat processing cost; ${\tilde{f}}_2$ denotes oil refining cost; ${\tilde{f}}_3$ denotes coke processing cost; ${\tilde{f}}_4$ denotes the cost of coal inputs for coal-based power; ${\tilde{f}}_5$ denotes the cost of natural gas inputs for natural gas-based power; ${\tilde{f}}_6$ denotes importing electricity cost; ${\tilde{f}}_7$ denotes the cost of first-stage electricity generation; ${\tilde{f}}_8$ denotes the cost of second-stage electricity generation; ${\tilde{f}}_9$ denotes electricity expansion cost; ${\tilde{f}}_{10}$ denotes the energy production cost; ${\tilde{f}}_{11}$ denotes the subsidy for solar power; ${\tilde{f}}_{12}$ denotes the subsidy for wind power; ${\tilde{f}}_{13}$ denotes the cost of charging pile; ${\tilde{f}}_{14}$ denotes the cost of changing station; ${\tilde{f}}_{15}$ denotes air-pollutants removal cost; ${\tilde{f}}_{16}$ denotes the risk recourse costs of stochastic and fuzzy uncertainties. $V{F}_t^{}$ represents the positive deviation between maximum values and worst value of fuzzy parameters. ${\lambda }_1$ and ${\lambda }_2$ represent their respective weight coefficients, where it is assumed ${\lambda }_1$ and ${\lambda }_2$ are fixed (i.e., ${\lambda }_1$ = ${\lambda }_2$ = 1).

(1) Cost of heat processing. This cost is used for heat processing, and it is calculated based on the unit-price and amount of heat processing.

$$\begin{array}{l}{\tilde{f}}_1={\displaystyle \sum _{t=1}^3[0.25·(HP{P}_t{}^1+2HP{P}_t{}^2+HP{P}_t{}^3) }·HG{A}_t\\ + {\lambda }_1·{\displaystyle \sum _{t=1}^3[HP{P}_t{}^3-0.25·(HP{P}_t{}^1+2·HP{P}_t{}^2+HP{P}_t{}^3)]·} HG{A}_t\\ + {\lambda }_2·{\displaystyle \sum _{t=1}^3[0.25·(HP{P}_t{}^1+2·HP{P}_t{}^2+HP{P}_t{}^3)-HP{P}_t{}^1]·} HG{A}_t\end{array}$$

(23)

where $H\tilde{P}{P}_t(HP{P}_{it}{}^1,HP{P}_{it}{}^2,HP{P}_{it}{}^3)$ is the cost for unit of heat processing, which is expressed as a triangular fuzzy number; $HG{A}_t$ is the heat processing amount.

(2) Cost for oil refining. This cost is calculated in terms of the unit-price and the amount of oil refining.

$$\begin{array}{l}{\tilde{f}}_2={\displaystyle \sum _{t=1}^3[0.25·(PV{O}_t{}^1+2·PV{O}_t{}^2+PV{O}_t{}^3) }·OFOI{L}_t\\ + {\lambda }_1·{\displaystyle \sum _{t=1}^3[PV{O}_t{}^3-0.25·(PV{O}_t{}^1+2·PV{O}_t{}^2+PV{O}_t{}^3)]·}OFOI{L}_t\\ + {\lambda }_2·{\displaystyle \sum _{t=1}^3[0.25·(PV{O}_t{}^1+2·PV{O}_t{}^2+PV{O}_t{}^3) -PV{O}_t{}^1]·} OFOI{L}_t\end{array}$$

(24)

where $P\tilde{V}{O}_t (PV{O}_t{}^1,PV{O}_t{}^2,PV{O}_t{}^3)$ is the cost for unit of oil refining; $OFOI{L}_t$ is the crude oil consumption of oil refining. (3) Cost for coke processing. This cost is calculated in terms of the unit-price and the amount of coke processing.

$$\begin{array}{l}{\tilde{f}}_3={\displaystyle \sum _{t=1}^3[0.25·(PW{O}_t{}^1+2·PW{O}_t{}^2+PW{O}_t{}^3) }·CKP{A}_t\\ + {\displaystyle \sum _{t=1}^3[PW{O}_t{}^3-0.25·(PW{O}_t{}^1+2·PW{O}_t{}^2+PW{O}_t{}^3)] ·}CKP{A}_t\\ + {\displaystyle \sum _{t=1}^3[0.25·(PW{O}_t{}^1+2·PW{O}_t{}^2+PW{O}_t{}^3) -PW{O}_t{}^1]·} CKP{A}_t\end{array}$$

(25)

where $P\tilde{W}{O}_t(PW{O}_t{}^1,PW{O}_t{}^2,PW{O}_t{}^3)$ is the cost for unit of coke processing; and $CKP{A}_t$ is the coke processing amount.

(4) Cost for purchasing coal resources. This cost is used for purchasing coal resources, and it is calculated based on the unit-price and the amount of coal resource.

$$\begin{array}{l}{\tilde{f}}_4={\displaystyle \sum _{t=1}^3[0.25·(NS{C}_t{}^1+2·NS{C}_t{}^2+NS{C}_t{}^3)] }·ECOAL{M}_t\\ + {\lambda }_1·{\displaystyle \sum _{t=1}^3[NS{C}_t{}^3-0.25·(NS{C}_t{}^1+2·NS{C}_t{}^2+NS{C}_t{}^3)]·}ECOAL{M}_t\\ + {\lambda }_2·{\displaystyle \sum _{t=1}^3[0.25·(NS{C}_t{}^1+2·NS{C}_t{}^2+NS{C}_t{}^3) -NS{C}_t{}^1]·} ECOAL{M}_t\end{array}$$

(26)

where $N\tilde{S}{C}_t(NS{C}_t{}^1,NS{C}_t{}^2,NS{C}_t{}^3)$ is the price for unit of coal resource; and $ECOAL{M}_t$ is the coal consumption of coal-based power.

(5) Cost for purchasing natural gas resources. This cost is used for purchasing natural gas resources, and it is calculated based on the unit-price and the amount of natural gas resource.

$$\begin{array}{l}{\tilde{f}}_5={\displaystyle \sum _{t=1}^3[0.25·(NS{N}_t{}^1+2·NS{N}_t{}^2+NS{N}_t{}^3)] }·ENG{M}_t\\ + {\lambda }_1·{\displaystyle \sum _{t=1}^3[NS{N}_t{}^3-0.25·(NS{N}_t{}^1+2·NS{N}_t{}^2+NS{N}_t{}^3)]·}ENG{M}_t\\ + {\lambda }_2·{\displaystyle \sum _{t=1}^3[0.25·(NS{N}_t{}^1+2·NS{N}_t{}^2+NS{N}_t{}^3) -NS{N}_t{}^1]·} ENG{M}_t\end{array}$$

(27)

where $N\tilde{S}{N}_t(NS{N}_t{}^1, NS{N}_t{}^2, NS{N}_t{}^3)$ is the price for unit of natural gas resource; and $ENG{M}_t$ is the natural gas consumption of natural gas-based power.

(6) Cost for importing electricity. This cost is calculated in terms of the unit-price and the amount of importing electricity.

$$\begin{array}{l}{\tilde{f}}_6={\displaystyle \sum _{t=1}^3{\displaystyle \sum _{h=1}^3[}0.25·(N{E}_t{}^1+2·N{E}_t{}^2+N{E}_t{}^3) }·E{D}_{th}\\ + {\lambda }_1·{\displaystyle \sum _{t=1}^3{\displaystyle \sum _{h=1}^3[}N{E}_t{}^3-0.25·(N{E}_t{}^1+2·N{E}_t{}^2+N{E}_t{}^3)]·}E{D}_{th}\\ + {\lambda }_2·{\displaystyle \sum _{t=1}^3{\displaystyle \sum _{h=1}^3[}0.25·(N{E}_t{}^1+2·N{E}_t{}^2+N{E}_t{}^3) -N{E}_t{}^1]·} E{D}_{th}\end{array}$$

(28)

where $N{\tilde{E}}_t(N{E}_t{}^1, N{E}_t{}^2, N{E}_t{}^3)$ is the price for unit of imported electricity; $E{D}_{th}$ is the amount of imported electricity.

(7) Operation cost for first-stage electricity generation. The cost represents the operation cost of first-stage electricity generation facilities (i.e., coal-based power, gas-based power, wind power and solar power) during the planning periods. It is calculated in terms of the operation costs and the amount of electricity generation for each of the electricity generation facilities. Furthermore, the first-stage is given by $W_{kt}^{\pm }=W_{kt}^{-}+r{r}_{kt}·\Delta W$, where $r{r}_{kt}$ denotes the decision variables, $\Delta W=W_{kt}^{+}-W_{kt}^{-}$.

$$\begin{array}{l}{\tilde{f}}_7={\displaystyle \sum _{k=1}^4{\displaystyle \sum _{t=1}^3[0.25·(P{V}_{kt}{}^1}} + 2·P{V}_{kt}{}^2+P{V}_{kt}{}^3) ]· (W_{kt}{}^{-}+r{r}_{kt}·\Delta W_{kt})\\ +{\lambda }_1{\displaystyle \sum _{k=1}^4{\displaystyle \sum _{t=1}^3[P{V}_{it}{}^3-0.25·(P{V}_{kt}{}^1}} + 2·P{V}_{kt}{}^2+P{V}_{kt}{}^3) ]· (W_{kt}{}^{-}+r{r}_{kt}·\Delta W_{kt})\\ +{\lambda }_2{\displaystyle \sum _{k=1}^4{\displaystyle \sum _{t=1}^3[0.25·(P{V}_{kt}{}^1}} + 2·P{V}_{kt}{}^2+P{V}_{kt}{}^3) -P{V}_{kt}{}^1]· (W_{kt}{}^{-}+r{r}_{kt}·\Delta W_{kt})\end{array}$$

(29)

where $P{\tilde{V}}_{kt}(P{V}_{kt}{}^1, P{V}_{kt}{}^2, P{V}_{kt}{}^3)$ denotes the operation cost of different power conversion technologies; $W_{kt}{}^{\pm }=W_{kt}{}^{-}+r{r}_{kt}·\Delta W_{kt}$ is the first-stage electricity generation amount, which is determined by decision variables $r{r}_{kt}$ ($\Delta W=W_{kt}^{+}-W_{kt}^{-}$ and $r{r}_{kt}\in [0, 1]$).

(8) Penalty cost for second-stage electricity generation amounts (i.e., shortage electricity amount of first-stage). It is calculated in terms of the penalty costs and the amount of second-stage electricity generation by each electricity generation technology.

$$\begin{array}{l}{\tilde{f}}_8={\displaystyle \sum _{k=1}^4{\displaystyle \sum _{t=1}^3{\displaystyle \sum _{h=1}^3[}0.25·(P{P}_{kt}{}^1+}}2·P{P}_{kt}{}^2+ P{P}_{kt}{}^3)]·p_{th}·Y_{kth}\\ + {\lambda }_1·{\displaystyle \sum _{k=1}^4{\displaystyle \sum _{t=1}^3{\displaystyle \sum _{h=1}^3[}P{P}_{kt}{}^3-0.25·(P{P}_{kt}{}^1+}}2·P{P}_{kt}{}^2+ P{P}_{kt}{}^3)]· p_{th}·Y_{kth}\\ +{\lambda }_2·{\displaystyle \sum _{k=1}^4{\displaystyle \sum _{t=1}^3{\displaystyle \sum _{h=1}^3[}0.25·(P{P}_{kt}{}^1+}}2·P{P}_{kt}{}^2+ P{P}_{kt}{}^3)-P{P}_{kt}{}^1]· p_{th}·Y_{kth}\end{array}$$

(30)

where $P{\tilde{P}}_{kt}(P{P}_{kt}{}^1, P{P}_{kt}{}^2, P{P}_{kt}{}^3)$ is the operating cost for the second-stage electricity generation amount ($Y_{kth}^{}$). $Y_{kth}^{}$ is the second-stage electricity generation amount.

(9) Fixed and variable costs for electric capacity expansion. These costs include the fixed and variable electric capacity expansion costs of four electricity generation technologies in the planning horizon.

$$\begin{array}{l}{\tilde{f}}_9={\displaystyle \sum _{i=1}^4{\displaystyle \sum _{t=1}^3{\displaystyle \sum _{h=1}^3{p}_{th}·}[0.25·(A_{kt}{}^1+}}2·A_{kt}{}^2+ A_{kt}{}^3)]·Q_{kth}\\ +{\lambda }_1·{\displaystyle \sum _{i=1}^4{\displaystyle \sum _{t=1}^3{\displaystyle \sum _{h=1}^3{p}_{th}·}[A_{kt}{}^3-0.25·(A_{kt}{}^1+}}2·A_{kt}{}^2+ A_{kt}{}^3)]·Q_{kth}\\ +{\lambda }_2·{\displaystyle \sum _{i=1}^4{\displaystyle \sum _{t=1}^3{\displaystyle \sum _{h=1}^3{p}_{th}·}[0.25·(A_{kt}{}^1+}}2·A_{kt}{}^2+ A_{kt}{}^3)-A_{kt}{}^1]·Q_{kth}\\ +{\displaystyle \sum _{i=1}^4{\displaystyle \sum _{t=1}^3{\displaystyle \sum _{h=1}^3{p}_{th}·}[0.25·(B_{kt}{}^1+}}2·B_{kt}{}^2+ B_{kt}{}^3)]·Z_{kth}\\ +{\lambda }_1·{\displaystyle \sum _{i=1}^4{\displaystyle \sum _{t=1}^3{\displaystyle \sum _{h=1}^3{p}_{th}·}[B_{kt}{}^3-0.25·(B_{kt}{}^1+}}2·B_{kt}{}^2+ B_{kt}{}^3)]·Z_{kth}\\ +{\lambda }_2·{\displaystyle \sum _{i=1}^4{\displaystyle \sum _{t=1}^3{\displaystyle \sum _{h=1}^3{p}_{th}·}[0.25·(B_{kt}{}^1+}}2·B_{kt}{}^2+ B_{kt}{}^3)-B_{kt}{}^1]·Z_{kth}\end{array}$$

(31)

where ${\tilde{A}}_{kt}(A_{kt}{}^1,A_{kt}{}^2,A_{kt}{}^3)$ denote fixed-charge cost for different electric capacity expansion. $Q_{kth}^{}$ are binary variables for identifying whether or not the capacity expansion needs to be undertaken by power conversion technology k. ${\tilde{B}}_{kt}(B_{kt}{}^1,B_{kt}{}^2,B_{kt}{}^3)$ denote the variable cost of capacity expansion. $Z_{kth}$ denotes the continuous variable of the capacity expansion amount.

(10) Generation costs for primary energy resources. This cost is calculated in terms of the unit-price of primary energy (i.e., coal, natural gas, and crude oil) and the amount of primary generation.

$$\begin{array}{l}{\tilde{f}}_{10}={\displaystyle \sum _{i=1}^3{\displaystyle \sum _{t=1}^{3}0.25·(EP{P}_{it}{}^1}} +2·EP{P}_{it}{}^2+EP{P}_{it}{}^3)· EP{A}_{it}\\ +{\lambda }_1·{\displaystyle \sum _{i=1}^3{\displaystyle \sum _{t=1}^3[}}EP{P}_{it}{}^3-0.25·(EP{P}_{it}{}^1+2·EP{P}_{it}{}^2+EP{P}_{it}{}^3)]· EP{A}_{it}\\ +{\lambda }_2·{\displaystyle \sum _{i=1}^3{\displaystyle \sum _{t=1}^3[}}0.25·(EP{P}_{it}{}^1+2·EP{P}_{it}{}^2+EP{P}_{it}{}^3)-EP{P}_{it}{}^1]· EP{A}_{it}\end{array}$$

(32)

where $E\tilde{P}{P}_{it}(EP{P}_{it}{}^1,EP{P}_{it}{}^2,EP{P}_{it}{}^3)$ is the production cost per unit of energy resource; and $EP{A}_{it}$ is the production amount.

(11) Subsidy for solar power. The government provides a subsidy for units of electricity generated by solar power, and it is calculated in terms of the unit-subsidy of solar power and the amount of electricity generation amount of solar power.

$$\begin{array}{l}{\tilde{f}}_{11}={\displaystyle \sum _{t=1}^3[0.25·(S{P}_t{}^1+2·S{P}_t{}^2+S{P}_t{}^3) }·(W_{3t}{}^{-}+r{r}_{3t}·\Delta W_{3t})\\ + {\displaystyle \sum _{t=1}^3[S{P}_t{}^3-0.25·(S{P}_t{}^1+2·S{P}_t{}^2+S{P}_t{}^3)] ·}(W_{3t}{}^{-}+r{r}_{3t}·\Delta W_{3t})\\ + {\displaystyle \sum _{t=1}^3[0.25·(S{P}_t{}^1+2·S{P}_t{}^2+S{P}_t{}^3) -S{P}_t{}^1]·} (W_{3t}{}^{-}+r{r}_{3t}·\Delta W_{3t})\end{array}$$

(33)

where $S{\tilde{P}}_t(S{P}_t{}^1,S{P}_t{}^2,S{P}_t{}^3)$ is the solar subsidy provided by government for unit of electricity generated by solar power.

(12) Subsidy for wind power. The government provides a subsidy for units of electricity generated by wind power, and it is calculated in terms of the unit-subsidy of wind power and the amount of electricity generation amount of wind power.

$$\begin{array}{l}{\tilde{f}}_{12}={\displaystyle \sum _{t=1}^3[0.25·(W{P}_t{}^1+2·W{P}_t{}^2+W{P}_t{}^3) }·(W_{4t}{}^{-}+r{r}_{4t}·\Delta W_{4t})\\ + {\displaystyle \sum _{t=1}^3[W{P}_t{}^3-0.25·(W{P}_t{}^1+2·W{P}_t{}^2+W{P}_t{}^3)] ·}(W_{4t}{}^{-}+r{r}_{4t}·\Delta W_{4t})\\ + {\displaystyle \sum _{t=1}^3[0.25·(W{P}_t{}^1+2·W{P}_t{}^2+W{P}_t{}^3) -W{P}_t{}^1]·} (W_{4t}{}^{-}+r{r}_{4t}·\Delta W_{4t})\end{array}$$

(34)

where $W{\tilde{P}}_t(W{P}_t{}^1,W{P}_t{}^2,W{P}_t{}^3)$ is the wind subsidy provided by government for unit of electricity generated by wind power.

(13) Costs for battery charging piles. It is calculated in terms of the average charging pile amount for unit of electric vehicle, the EVs population, as well as the investment cost per EV charging pile.

$$f_{13}={\displaystyle \sum _{t=1}^{3}EV{A}_t·RE{C}_t·FC{D}_t}$$

(35)

where $EV{A}_t$ is the EVs population; $RE{C}_t$ is the average charging pile amount for unit of electric vehicle; $FC{D}_t$ is the investment cost per EV charging pile.

(14) Costs for battery changing station. It is calculated in terms of the average changing station amount for unit of electric vehicle, the EVs population, as well as the investment cost per EV changing station.

$$f_{14}={\displaystyle \sum _{t=1}^{3}EV{A}_t·RE{V}_t·FV{D}_t}$$

(36)

where $RE{V}_t$ is the average changing station amount for unit of electric vehicles; $FV{D}_t$ is the investment cost per EV changing station.

(15) Costs for pollutant reduction. These costs include reduction costs for three pollutants, i.e., SO₂, NO_X, and PM emissions.

$${\tilde{f}}_{15}= {\tilde{f}}_{15s}+{\tilde{f}}_{15N}+{\tilde{f}}_{15PM}$$

(37)

where ${\tilde{f}}_{15s}$ is the cost of desulfurization; $SO2K_t$ is the SO₂ emissions from coking processing; $SO2C{E}_t$ is the SO₂ emissions from coal-based power; $SO2N{E}_t$ is the SO₂ emissions from natural gas-based power; $SO2H_t$ is the SO₂ emissions from heating processing; $SO2O_t$ is the SO₂ emissions from oil refining; $T\tilde{P}{S}_t(TP{S}_t{}^1,TP{S}_t{}^2,TP{S}_t{}^3)$ is the desulfurization cost of per unit of SO₂ emission; ${\eta }_s$ is the desulfurization efficiency.

(15a) Costs for SO₂ emissions reduction. The cost is calculated in terms of the total SO₂ emissions, the desulfurization cost per unit of SO₂ emission, and the desulfurization efficiency.

$$\begin{array}{l}{\tilde{f}}_{15S}= \{{\displaystyle \sum _{t=1}^3(}SO2K_t+SO2C{E}_t+ SO2N{E}_t+SO2H_t+ SO2O_t)\\ ·0.25· (TP{S}_t{}^1+TP{S}_t{}^2+TP{S}_t{}^3 )· {\eta }_s\}\\ +{\lambda }_1·\{{\displaystyle \sum _{t=1}^3(}SO2K_t+SO2C{E}_t+ SO2N{E}_t+SO2H_t+ SO2O_t)\\ ·[TP{S}_t{}^3- 0.25· (TP{S}_t{}^1+TP{S}_t{}^2+TP{S}_t{}^3 )]· {\eta }_s\}\\ +{\lambda }_2·\{{\displaystyle \sum _{t=1}^3(}SO2K_t+SO2C{E}_t+ SO2N{E}_t+SO2H_t+ SO2O_t)\\ ·[0.25· (TP{S}_t{}^1+TP{S}_t{}^2+TP{S}_t{}^3 )-TP{S}_t{}^1]· {\eta }_s\}\end{array}$$

(38)

where ${\tilde{f}}_{15s}$ is the cost of desulfurization; $SO2K_t$ is the SO₂ emissions from coking processing; $SO2C{E}_t$ is the SO₂ emissions from coal-based power; $SO2N{E}_t$ is the SO₂ emissions from natural gas-based power; $SO2H_t$ is the SO₂ emissions from heating processing; $SO2O_t$ is the SO₂ emissions from oil refining; $T\tilde{P}{S}_t(TP{S}_t{}^1,TP{S}_t{}^2,TP{S}_t{}^3)$ is the desulfurization cost per unit of SO₂ emission; ${\eta }_s$ is the desulfurization efficiency.

(15b) Costs for NO_X emissions reduction. The cost is calculated in terms of the total NOx emissions, the desulfurization cost per unit of NOx emission, and the desulfurization efficiency.

$$\begin{array}{l}{\tilde{f}}_{15N}= \{{\displaystyle \sum _{t=1}^3(}NOX{K}_t+NOXC{E}_t+ NOXN{E}_t+NOX{H}_t+ NOX{O}_t)\\ ·0.25· (TP{N}_t{}^1+TP{N}_t{}^2+TP{N}_t{}^3 )· {\eta }_N\}\\ +{\lambda }_1·\{{\displaystyle \sum _{t=1}^3(}NOX{K}_t+NOXC{E}_t+ NOXN{E}_t+NOX{H}_t+ NOX{O}_t)\\ ·[TP{N}_t{}^3- 0.25· (TP{N}_t{}^1+TP{N}_t{}^2+TP{N}_t{}^3 )]· {\eta }_N\}\\ +{\lambda }_2·\{{\displaystyle \sum _{t=1}^3(}NOX{K}_t+NOXC{E}_t+NOXN{E}_t+NOX{H}_t+NOX{O}_t)\\ ·[0.25· (TP{N}_t{}^1+TP{N}_t{}^2+TP{N}_t{}^3 )-TP{N}_t{}^1]· {\eta }_N\}\end{array}$$

(39)

where ${\tilde{f}}_{15N}$ is the cost of denitration; $NOX{K}_t$ is the NO_x emissions from coking processing; $NOXC{E}_t$ is the NO_x emissions from coal-based power; $NOXN{E}_t$ is the NO_x emissions from natural gas-based power; $NOX{H}_t$ is the NO_x emissions from heat processing; $NOX{O}_t$ is the NO_x emissions from oil refining; $T\tilde{P}{N}_t(TP{N}_t{}^1,TP{N}_t{}^2,TP{N}_t{}^3)$ is the denitration cost per unit of NO_x emission; ${\eta }_N$ is the denitration efficiency.

(15c) Costs for PM emissions reduction. The cost is calculated in terms of the total PM emissions, the desulfurization cost per unit of PM emission, and the desulfurization efficiency.

$$\begin{array}{l}{\tilde{f}}_{15PM}= \{{\displaystyle \sum _{t=1}^3(}PM{K}_t+PMC{E}_t+ PMN{E}_t+PM{H}_t+ PM{O}_t)\\ ·0.25· (TPP{M}_t{}^1+TPP{M}_t{}^2+TPP{M}_t{}^3 )· {\eta }_{PM}\}\\ +{\lambda }_1·\{{\displaystyle \sum _{t=1}^3(}PM{K}_t+PMC{E}_t+ PMN{E}_t+PM{H}_t+PM{O}_t)\\ ·[TPP{M}_t{}^3- 0.25· (TPP{M}_t{}^1+TPP{M}_t{}^2+TPP{M}_t{}^3 )]· {\eta }_{PM}\}\\ +{\lambda }_2·\{{\displaystyle \sum _{t=1}^3(}PM{K}_t+PMC{E}_t+PMN{E}_t+PM{H}_t+PM{O}_t)\\ ·[0.25· (TPP{M}_t{}^1+TPP{M}_t{}^2+TPP{M}_t{}^3 )-TPP{M}_t{}^1]· {\eta }_{PM}\}\end{array}$$

(40)

where ${\tilde{f}}_{15PM}$ is the PM removal cost; $PM{K}_t$ is the PM emissions from coking processing; $PMC{E}_t$ is the PM emissions from coal-based power; $PMN{E}_t$ is the PM emissions from natural gas-based power; $PM{H}_t$ is the PM emissions from heating processing; $PM{O}_t$ is the NO_X emissions from oil refining; $T\tilde{P}P{M}_t(TPP{M}_t{}^1,TPP{M}_t{}^2,TPP{M}_t{}^3)$ is the treatment cost per unit of PM emission; ${\eta }_{PM}$ is the PM removal efficiency.

(16) Costs for risk recourse of stochastic and fuzzy uncertainties.

$${\tilde{f}}_{16}=V\tilde{T}S+V\tilde{T}F={\rho }_1·{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{h=1}^H{p}_{th}}}}(V_{ith}+V{C}_{ith}]+{\rho }_2·{\displaystyle \sum _{t=1}^{3}V{F}_t}$$

(41)

where $V\tilde{T}S$ denotes the risk recourse cost for stochastic uncertainties; $V\tilde{T}F$ denotes the risk recourse cost of fuzzy uncertainties, with capturing the difference between the extreme values of fuzzy parameters. The nonnegative factors ρ₁ and ρ₂ represent robust levels (its range from 0 to 1), which can help managers make trade-offs between system economy and reliability.

Constraints:

The constraints of the proposed DRSFO-EES model include traffic sector, heat processing, coke processing, oil refining, electricity generation, energy production, air-pollutants treatment, and risk recourse cost for stochastic uncertainties and risk recourse for fuzzy uncertainties.

(1) Constraints for traffic sector.

This constraint represents that the optimized vehicles population must be not less than the lower bounds of vehicle population. $Cr\{M\tilde{A}{L}_{gt}|M{A}_{gt}\ge M\tilde{A}{L}_{gt}\}\ge \alpha$ means the credibility of the $M{A}_{gt}\ge M\tilde{A}{L}_{gt}$ is higher than or equal to confidence-level α.

$$Cr\{M\tilde{A}{L}_{gt}|M{A}_{gt}\ge M\tilde{A}{L}_{gt}\}\ge \alpha$$

(42)

where $M{A}_{gt}$ is the optimized solutions of vehicles population; $M\tilde{A}{L}_{gt}(MA{L}_{gt}{}^1, MA{L}_{gt}{}^2, MA{L}_{gt}{}^3)$ is the lower bounds of vehicle population, which is expressed as a triangular fuzzy number.

Constraints for vehicle-emissions, including CO, NOx, HC, and PM emissions. These constraints represent that the vehicle-emissions are calculated in terms of the annual average mileage, the proportion of electric vehicles and the emission factors of CO, NOx, HC and PM emissions.

$$TRC{O}_{gt|g=1,3,4,5}^{}=M{A}_{gt|g=1,3,4,5}^{}·M{Q}_{gt|g=1,3,4,5}^{}·TRCO{E}_{gt|g=1,3,4,5}$$

(43)

$$TRC{O}_{2t}^{}=(M{A}_{gt}-M{A}_{gt}·RE{V}_t)·M{Q}_{gt}^{}·TRCO{E}_{gt}$$

(44)

$$TRNO{X}_{gt|g=1,3,4,5}^{}=M{A}_{gt|g=1,3,4,5}^{}·M{Q}_{gt|g=1,3,4,5}^{}·TRNO{E}_{gt|g=1,3,4,5}$$

(45)

$$TRNO{X}_{2t}^{}=(M{A}_{gt}-M{A}_{gt}·RE{V}_t)·M{Q}_{gt}^{}·TRNO{E}_{gt}$$

(46)

$$TRH{C}_{gt|g=1,3,4,5}^{}=M{A}_{gt|g=1,3,4,5}^{}·M{Q}_{gt|g=1,3,4,5}^{}·TRH{C}_{gt|g=1,3,4,5}$$

(47)

$$TRH{C}_{2t}^{}=(M{A}_{gt}-M{A}_{gt}·RE{V}_t)·M{Q}_{gt}^{}·TRH{C}_{gt}$$

(48)

$$TRP{M}_{gt|g=1,3,4,5}^{}=M{A}_{gt|g=1,3,4,5}^{}·M{Q}_{gt|g=1,3,4,5}^{}·TRP{M}_{gt|g=1,3,4,5}$$

(49)

$$TRP{M}_{2t}^{}=(M{A}_{gt}-M{A}_{gt}·RE{V}_t)·M{Q}_{gt}^{}·TRP{M}_{gt}$$

(50)

where $TRC{O}_{gt}$ is CO emissions from vehicle g in period t; $M{Q}_{gt}^{}$ is the annual average mileage; $RE{V}_t$ is the proportion of electric vehicles; $TRCO{E}_{gt}$ the emission factor of CO; $TRNO{E}_{gt}$ the emission factor of NO_X; $TRH{C}_{gt}$ the emission factor of HC; $TRP{M}_{gt}$ the emission factor of PM.

Constraints of the electricity demands of EVs. This constraint is calculated in terms of the electricity consumption amounts per hundred kilometers of EVs and the amount of EVs.

$$EVE{C}_t=M{A}_{2t}·RE{V}_t·M{Q}_{2t}·ECPE{V}_t$$

(51)

where $EVE{C}_t$ is the electricity consumption amounts of EVs; $ECPE{V}_t$ is the electricity consumption amounts per hundred kilometers of EVs.

(2) Constraints of heat processing. Constraint (52) depicts that the heat generation amounts must not be less than the heat demands, and $Cr\{D\tilde{M}{H}_t|HG{A}_t\ge D\tilde{M}{H}_t\}\ge \alpha$ means the credibility of the $HG{A}_t\ge D\tilde{M}{H}_t$ is higher than or equal to confidence-level α. Constraints (53) and (54) represent the consumption amounts of natural gas and coal for heat processing.

$$Cr\{D\tilde{M}{H}_t|HG{A}_t\ge D\tilde{M}{H}_t\}\ge \alpha$$

(52)

$$CFHG{C}_t=CF{H}_t·HG{A}_t·HC{P}_t$$

(53)

$$CFHG{N}_t=CF{H}_t·HG{A}_t·HN{P}_t$$

(54)

where $D\tilde{M}{H}_t(DM{H}_t{}^1,DM{H}_t{}^2,DM{H}_t{}^3)$ is the heat demand; $CF{H}_t$ denotes the energy consumption for the unit of kerosene processing. $CFHG{C}_t$ denotes the coal consumption of heat processing; $CFHG{N}_t$ is the natural gas consumption of heat processing; $HC{P}_t$ denotes the ratio of coal consumption of heat processing; $HN{P}_t$ denotes the ratio of natural gas consumption of heat processing.

(3) Constraints of oil refining. The constraint depicts the input amounts of crude oil for oil refining.

$$\begin{array}{l}OFOI{L}_t=ZC{Y}_t·ZCY{E}_t+ZQ{Y}_t·ZQY{E}_t+ZM{Y}_t·ZMY{E}_t\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+ ZRL{Y}_t · ZRLY{E}_t + ZSN{Y}_t · ZSNY{E}_t\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+ ZSJ{Y}_t · ZSJY{E}_t + ZLP{G}_t · ZLPG{E}_t + ZQ{T}_t\end{array}$$

(55)

where $ZC{Y}_t$ is the diesel processing amount; $ZCY{E}_t$ denotes the crude oil consumption for unit of diesel processing; $ZQ{Y}_t$ is the processing amount of gasoline; $ZQY{E}_t$ denote the crude oil consumption for unit of gasoline processing; $ZM{Y}_t$ is the processing amount of kerosene; $ZMY{E}_t$ denotes the crude oil consumption for unit of kerosene processing; $ZRL{Y}_t$ is the fuel oil processing amount; $ZRLY{E}_t$ is the crude consumption for unit of fuel oil processing; $ZSN{Y}_t$ is the processing amount for naphtha; $ZSNY{E}_t$ denotes the crude oil consumption for unit of naphtha processing; $ZSJ{Y}_t$ is the processing amount of asphaltic pyrobitumen; $ZSJY{E}_t$ denotes the crude oil consumption for unit of asphaltic pyrobitumen processing; $ZLP{G}_t$ is the processing amount of LPG; $ZLPG{E}_t$ denotes the crude consumption for unit of LPG processing; $ZQ{T}_t$ denotes the crude oil consumption of other oil products.

(4) Constraints of coke processing. Constraint (56) means that the amount of coke processing is not less than the coke demands, and $Cr\{DM\tilde{C}{J}_t|CKP{A}_t\ge DM\tilde{C}{J}_t\}\ge \alpha$ means the credibility of the $CKP{A}_t\ge DM\tilde{C}{J}_t$ is higher than or equal to confidence-level α. Constraint (57) depicts the coal consumption amount for coke processing.

$$Cr\{DM\tilde{C}{J}_t|CKP{A}_t\ge DM\tilde{C}{J}_t\}\ge \alpha$$

(56)

$$CFC{J}_t = CKP{A}_t/CTCK{E}_t$$

(57)

where $DM\tilde{C}{J}_t(DMC{J}_t{}^1,DMC{J}_t{}^2,DMC{J}_t{}^3)$ is the coke demand expressed as a triangular fuzzy number. $CTCK{E}_t$ denotes the coal consumption for unit of coke processing; $CFC{J}_t$ is the coal consumption of coke processing.

(5) Constraints of electricity generation.

Constraints for mass balance of coal and natural gas resources. These constraints are established to calculate the consumption amounts of coal and natural gas for electricity generation.

$$[(W_{1t}^{-}+r{r}_{1t}^{}·W{C}_{1t}^{})+Y_{1th}+EVE{C}_t]·F{E}_{1t}=ECOAL{M}_t$$

(58)

$$[(W_{2t}^{-}+r{r}_{2t}^{}·W{C}_{2t}^{})+Y_{2th}]·F{E}_{2t}=ENG{M}_t$$

(59)

where $F{E}_{it}$ is the coal (i = 1) and natural gas (i = 2) consumption for a unit of electricity generation.

Constraints of electricity demand and supply balance. Constraint (60) is established to ensure the electricity demand be satisfied by domestic electricity generation and importation. And $Cr\{{\tilde{D}}_{th}|{\displaystyle \sum _{i=1}^4{\displaystyle \sum _{t=1}^3[(W_{it}^{-}+r{r}_{it}·W{C}_{it})+Y_{ith}}}+EVE{C}_t]+E{d}_{th}\ge {\tilde{D}}_{th}\}\ge \alpha$ means the credibility of the ${\displaystyle \sum _{i=1}^4{\displaystyle \sum _{t=1}^3[(W_{it}^{-}+r{r}_{it}·W{C}_{it})+Y_{ith}}}+EVE{C}_t]+E{d}_{th}\ge {\tilde{D}}_{th}$ is higher than or equal to confidence-level α. For constraint (61), the optimized amount of electricity generated in the first-stage is given by $W_{kt}^{\pm }=W_{kt}^{-}+r{r}_{kt}·\Delta W$, where $r{r}_{kt}$ denotes the decision variables, and $r{r}_{kt}\in [0, 1]$.

$$Cr\{{\tilde{D}}_{th}|{\displaystyle \sum _{i=1}^4{\displaystyle \sum _{t=1}^3[(W_{it}^{-}+r{r}_{it}·W{C}_{it})+Y_{ith}}}+EVE{C}_t]+E{d}_{th}\ge {\tilde{D}}_{th}\}\ge \alpha$$

(60)

$$0\le r{r}_{it}\le 1$$

(61)

where ${\tilde{D}}_{th}$ is the electricity demand under various electricity demand levels, which show the characteristics of stochastic and fuzzy sets. $r{r}_{it}$ denotes the decision variable of first-stage electricity generation.

Constraints for electricity capacities. These constraints mainly depict that the amount of generated electricity must not exceed its existing and expanded capacities, which are established to ensure that the available electricity-generation capacity is greater than the generated electricity.

$${\displaystyle \sum _{i=1}^4(W_{it}^{-}}+r{r}_{it}^{}·W{C}_{it}^{})+Y_{ith}^{}+EVE{C}_t\le {\displaystyle \sum _{i=1}^2[C{F}_{it}^{}}·(R{C}_{it}^{}+Z_{ith)}^{}]$$

(62)

$${\displaystyle \sum _{i=1}^4(R{C}_i}+Z_{ith})\ge U_t;\forall t, h$$

(63)

$${\displaystyle \sum _{i=1}^4(R{C}_i}+Z_{ith})\le U{S}_t;\forall t, h$$

(64)

$${\displaystyle \sum _{i=1}^4{Y}_{ith}^{}}\le {\displaystyle \sum _{i=1}^{4}C{F}_{it}^{}}·Z_{ith}^{}$$

(65)

$$0\le Y_{ith}\le (W_{it}^{-}+r{r}_{it}·W{C}_{it});\forall i, t, h$$

(66)

$$Q_{ith}=\{\begin{array}{l}=1, \mathrm{if} \mathrm{capacity} \mathrm{expansion} \mathrm{of} \mathrm{is} \mathrm{undertaken} \forall i, t, h\\ =0, \mathrm{if} \mathrm{otherwise}\end{array}$$

(67)

where $C{F}_{kt}$ is the operating hours of electricity conversion technology k; $R{C}_{kt}$ is the original capacity; $U_t^{}$ is the lower bound of load demands; $U{S}_t$ is the upper bound of load demands.

(6) Constraints of energy production. The constraint is established to ensure the primary energy production amount must be less than the lower bound and higher than the upper bound.

$$EPAL\le EP{A}_{it} \le EPA{U}_{it}$$

(68)

where $EPA{U}_{it}$ is the upper bound of amount of primary energy production. $EPA{L}_{it}$ is the lower bound of amount of primary energy production.

(7) Constraints for air-pollutants management.

SO₂ emissions. Constraints (69) to (74) represent the SO₂ emissions from energy activities (i.e., coke processing, heat processing, electricity generation, and oil refining). Constraints (74) are used for ensuring that the SO₂-emissions are satisfied by the pollutant-emission permits.

$$SO2K_t=CFC{J}_t·C{S}_t$$

(69)

$$SO2C{E}_t=ECOAL{M}_t· C{S}_t$$

(70)

$$SO2N{E}_t=NEG{M}_t·C{S}_t$$

(71)

$$SO2H_t=CFHG{C}_t·C{S}_t+CFHG{N}_t·N{S}_t$$

(72)

$$SO2O_t=OFOI{L}_t·O{S}_t$$

(73)

$$(SO2K_t+SO2C{E}_t+SO2N{E}_t+SO2H_t+SO2O_t) · (1-{\eta }_s) \le T{S}_t$$

(74)

where $C{S}_t$ is the SO₂ emission factor of coal; $N{S}_t$ is the SO₂ emission factor of natural gas; $O{S}_t$ is the SO₂ emission factor of crude oil; ${\eta }_s$ is the desulfurization efficiency; $T{S}_t$ is the upper bounds of SO₂ emission.

NO_X emissions. Constraints (75) to (79) represent the NOx emissions from energy activities (i.e., coke processing, heat processing, electricity generation, and oil refining). Constraints (80) are used for ensuring that the NOx-emissions are satisfied by the pollutant-emission permits.

$$NOX{K}_t=CFC{J}_t·C{N}_t$$

(75)

$$NOXC{E}_t=ECOAL{M}_t· C{N}_t$$

(76)

$$NOXN{E}_t=NEG{M}_t·N{N}_t$$

(77)

$$NOX{H}_t=CFHG{C}_t·C{N}_t+CFHG{N}_t·N{N}_t$$

(78)

$$NOX{O}_t=OFOI{L}_t·O{N}_t$$

(79)

$$(NOX{K}_t+NOXC{E}_t+NOXN{E}_t+NOX{H}_t+NOX{O}_t) · (1-{\eta }_N) \le T{N}_t$$

(80)

where $C{N}_t$ is the NO_x emission factor of coal; $N{N}_t$ is the NO_x emission factor of natural gas; $O{N}_t$ is the NO_x emission factor of crude oil; ${\eta }_N$ is the denitration efficiency; $T{N}_t$ is the upper bounds of NO_x emission.

PM emissions. Constraints (81) to (85) represent the PM emissions from energy activities (i.e., coke processing, heat processing, electricity generation, and oil refining). Constraints (86) are used for ensuring that the PM-emissions are satisfied by the pollutant-emission permits.

$$PM{K}_t=CFC{J}_t·CP{M}_t$$

(81)

$$PMC{E}_t=ECOAL{M}_t· CP{M}_t$$

(82)

$$PMN{E}_t=NEG{M}_t·NP{M}_t$$

(83)

$$PM{H}_t=CFHG{C}_t·CP{M}_t+CFHG{N}_t·NP{M}_t$$

(84)

$$PM{O}_t=OFOI{L}_t·OP{M}_t$$

(85)

$$(PM{K}_t+PMC{E}_t+PMN{E}_t+PM{H}_t+PM{O}_t) · (1-{\eta }_{PM}) \le TP{M}_t$$

(86)

where $CP{M}_t$ is the PM emission factor of coal; $NP{M}_t$ is the PM emission factor of natural gas; $OP{M}_t$ is the PM emission factor of crude oil; ${\eta }_{PM}$ is the PM removal efficiency; $TP{M}_t$ is the upper bound of PM emission.

CO₂ emissions. Constraints (81) to (85) represent the CO₂- emissions from energy activities (i.e., coke processing, heat processing, electricity generation, and oil refining).

$$CO2K_t=CFC{J}_t·CCO{2}_t$$

(87)

$$CO2C{E}_t=ECOAL{M}_t· CCO{2}_t$$

(88)

$$CO2N{E}_t=NEG{M}_t·NCO{2}_t$$

(89)

$$CO2H_t=CFHG{C}_t·CCO{2}_t+CFHG{N}_t·NCO{2}_t$$

(90)

$$CO2O_t=OFOI{L}_t·OCO{2}_t$$

(91)

where $CCO{2}_t$ is the CO₂ emission factor of coal; $NCO{2}_t$ is the CO₂ emission factor of natural gas; $OCO{2}_t$ is the CO₂ emission factor of crude oil.

(8) The expected deviations for stochastic uncertainties. These constraints are used for capturing the risk from stochastic uncertainties.

$$V_{ith}^{}=P{\tilde{P}}_{it}^{}·Y_{ith}^{\pm }-{\displaystyle \sum _{h=1}^H{p}_{th}P{\tilde{P}}_{it}^{}·Y_{ith}^{\pm }}+2{\theta }_h;\forall i,t,h$$

(92)

$$V{C}_{ith}^{}=({\tilde{A}}_{it}^{}{Q}_{ith}^{}+{\tilde{B}}_{it}^{}{Z}_{ith}^{})-{\displaystyle \sum _{h=1}^H{p}_{th}({\tilde{A}}_{it}^{}{Q}_{ith}^{}+{\tilde{B}}_{it}^{}{Z}_{ith}^{})}+2{\theta }_h;\forall i,t,h$$

(93)

$$V_{ith}^{\pm }\ge 0, V{C}_{i{j}_{c}th}^{\pm }\ge 0;$$

(94)

where ${\theta }_h\ge 0$ are slack variables that can achieve looser constraints; $V_{ith}^{}$ and $V_{i{j}_{c}th}^{}$ are the weighted values of the expected deviations from stochastic uncertainties.

(9) The expected deviations for fuzzy uncertainties. These constraints are used for capturing the risk from fuzzy uncertainties.

$$\begin{array}{l}V{F}_t=(DMC{J}_t{}^R(\alpha )-DMC{J}_t{}^1)+(DM{H}_t{}^R(\alpha )-DM{H}_t{}^1)\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}+{\displaystyle \sum _{g=1}^5(MA{L}_{gt}{}^R(\alpha )}-MA{L}_{gt}{}^1)+{\displaystyle \sum _{h=1}^3[}{D}_{th}{}^R(\alpha )-D_{th}{}^1];\forall g,t,h\end{array}$$

(95)

$$V{F}_t^{}\ge 0$$

(96)

where $V{F}_t^{}$ represents the positive deviation between maximum values and worst value of fuzzy parameters of right-side of constraints. $DMC{J}_t{}^R(\alpha )$, $DM{H}_t{}^R(\alpha )$, $MA{L}_t{}^R(\alpha )$, and $D_{th}{}^R(\alpha )$ are the maximum values of all potential values, i.e., $DMC{J}_t{}^R(\alpha )=\sup \{B|B={\mu }^{-1}(\alpha )\}$, ${\mu }^{-1}$ is the inverse of $\mu$. $DMC{J}_t{}^1(\alpha )$, $DM{H}_t{}^1(\alpha )$, $MA{L}_t{}^1(\alpha )$, and $D_{th}{}^1(\alpha )$ represent the possible worst values of fuzzy numbers $DM\tilde{C}{J}_t(\alpha )$, $D\tilde{M}{H}_t(\alpha )$, $M\tilde{A}{L}_t(\alpha )$, and ${\tilde{D}}_{th}(\alpha )$, respectively.

(10) Non-negative constraints. This constraint assures that only positive variables are considered in the solutions, eliminating infeasibility while calculating the solution.

$$\begin{array}{l}M{A}_{gt},EP{A}_{it}, HG{A}_t, CFHG{C}_t, CFHG{N}_t, OFOI{L}_t,\\ CKP{A}_t,ECOAL{M}_t,ENG{M}_t,E{d}_{th},r{r}_{it},EVE{C}_t,\\ Y_{ith},Q_{ith},Z_{ith}\ge 0\end{array}$$

(97)

3.4. Data Acquirement

The data sources of the DRSFO-EES are related references, government reports, and statistical yearbooks. The data regarding the vehicle population of the BTH region in particular was obtained from reference [38]. The annual average vehicle kilometers travelled (VKT) data were obtained from related reference [40,41,42].

Table 2. Vehicular emission factors with different emissions standards (mg/km•unit).

		CO	NOX	HC	PM
HDV	China-Ⅴ	300	4610	35	100
	China-Ⅵ	300	4610	35	100
LDV	China-Ⅴ	1400	60	230	5
	China-Ⅵ	700	35	115	5
LDT	China-Ⅴ	5800	60	1200	5
	China-Ⅵ	2900	35	600	5
HDT	China-Ⅴ	200	3530	35	100
	China-Ⅵ	200	3530	35	100
Other	China-Ⅴ	2750	150	855	20
	China-Ⅵ	2750	150	855	20

details the vehicular emission factors of China-V and China-VI [40,41,43,44,45,46,47,48].

Table 3. Electricity demand and technological data.

	Period
	Period 1	Period 2	Period 3
Electricity demand (109 kWh)
Low demand level	[517.71, 537.71, 557.71]	[553.82, 573.82, 593.82]	[594.67, 614.67, 634.67]
Medium demand level	[547.25, 567.25, 587.25]	[595.73, 615.73, 635.73]	[636.18, 656.18, 676.18]
High demand level	[559.07, 579.07, 599.07]	[646.05, 666.05, 686.05]	[686.87, 706.87, 726.87]
Electricity generation target (109 kWh)
Coal-fired power	[221.07, 279.81]	[213.27, 268.01]	208.93, 269.50]
Gas-fired power	[58.41, 63.41]	[76.91, 80.91]	[91.88, 95.88]
Wind	[40.19, 49.55]	[45.79, 56.48]	[51.01, 61.94]
Solar power	[1.20, 1.70]	[3.20, 3.60]	[5.99, 6.22]
Energy consumption amounts per unit of electricity production
Coal (ton of SCE/103 kWh)	30.50	30.50	30.50
Natural gas (m3/103 kWh)	142.80	142.80	142.80
Coal consumption amounts per unit of coke processing (ton of SCE/ton)
	1.35	1.35	1.35
Energy consumption amounts for unit of heat processing (ton of SCE/109 kJ)
	36.00	36.00	36.00

gives the electricity demands and representative technical data [1,23,49,50,51,52]. The electricity demands were expressed as stochastic uncertain with three probability levels (20%, 60%, and 20% corresponding to low, medium, and high levels of electricity demand, respectively). The residual capacities and operation times of the different electricity generation technologies were obtained from [49]. The CO₂ and air-pollutant emission factors were acquired from related references and Intergovernmental Panel on Climate Change (IPCC) reports [2,3,53,54].

4. Results analysis

4.1. Analysis of Vehicular Emissions of BTH Region

Several policies on vehicular emission mitigation, such as the development of EVs, EV power sources and vehicular emission standards, have been proposed by the government of China. In this study, five scenarios labeled S1–S5 were designed to analyze the potentials of different emission mitigation strategies and policies for reducing vehicular emissions (NOx, HC, CO, and PM) in the BTH region.

Table 4. Parameter settings of scenarios.

	Vehicular Emissions Standards	The Proportion of EVs	Power Sources for EVs
S1	China Ⅴ	0%	—
S2	China Ⅴ	1.50%	100% coal-fired power based
S3	China Ⅴ	1.50%	50% coal-fired power based, 50% renewable energy based
S4	China Ⅴ	1.50%	100% renewable energy based
S5	China Ⅵ	1.50%	100% renewable energy based

gives the parameter settings of the different scenarios. Following is a description of each scenario:

• S1: Without consideration of EVs, and with the local governments implementing the China V vehicular emission standard.
• S2: With consideration of EVs, which account for 1.5% of the LDVs population. EV power sources are 100% based on coal-based power, and the China V vehicular emission standard is implemented.
• S3: With consideration of EVs, which account for 1.5% of the LDVs population. EV power sources are based 50% on coal-based and 50% on renewable power. It is assumed that coal-based power is local coal-based power, and renewable power included wind, solar power and imported electricity. The China V vehicular emission standard is implemented.
• S4: With consideration of EVs, which account for 1.5% of the LDVs population and with power sources based 100% on renewable energy. The China V vehicular emission standard is implemented.
• S5: With consideration of EVs, which account for 1.5% of the LDVs population and with power sources based 100% on renewable energy. The China VI vehicular emission standard (gasoline standard) is implemented.

4.1.1. Contributions of Different Vehicle Categories to Vehicular Emissions

Numerous results were obtained by the DRSFO-EES model, and the situation of α = 0.5, ρ = 0.2, and scenario S2 is used here as an example to illustrate the optimized solutions for the EES. The rapid development of traffic systems has led to large amounts of vehicular emissions such as CO, HC, NO_x, and PM. Figure 3 presents the contributions of the different vehicle categories to vehicular emissions in the BTH region over the planning periods. As indicated in Figure 3a–d, due to the sharp growth of LDVs, the total CO, NO_x, HC, and PM emissions are expected to increase by 23.28%, 17.07%, 17.19%, and 22.22% over the planning periods, respectively. The results indicated that LDVs are the major contributors of CO and HC emissions, to which they contribute 65.46–73.72% and 60.15–69.15%, respectively, over the planning horizon. It is therefore necessary for local governments to limit the number of LDVs through appropriate measures and policies such as improvement of the public transit system and encouragement of the use of EVs.

Further, HDTs are the major contributor of NO_x and PM, to which they are expected to contribute 67.07–63.08% and 58.58%–54.98%, respectively, between 2020 and 2030. Limiting the development of HDTs would therefore be effective for reducing vehicular NO_x and PM emissions. The results also indicated that the traffic sector is set to be one of the major emitters for NO_x in the BTH region. For example, the NO_x emission from vehicles will account for 58.65%, 60.13%, and 61.90% of the total NO_x emission (the total NO_x emission includes contributions from oil refining, coke processing, heat processing, natural gas-based power, and coal-based power) in 2020, 2025, and 2030, respectively.

4.1.2. Analysis of the Impact of Traffic Policies on Vehicle Pollutants

Figure 4 shows the vehicular emissions for different scenarios and can be used to analyze the pollutant emission mitigation potential of different vehicle policies. As per the results, implementation of the stringent China VI emission standard is effective for mitigating vehicular emissions, especially CO and HC emissions. The standard is specifically projected to afford CO, HC, NO_x, and PM emission reductions of 49.08%, 49.00%, 4.88%, and 0.00%, respectively, in 2020. As also indicated in Figure 4, the adoption of an EVs policy could reduce CO, NO_x, HC, and PM emissions from vehicles by 25.01×10³, 1.07×10³, 4.11×10³, and 0.08×10³ ton, respectively, by 2020. However, approximately 2.6×10⁹, 3.2×10⁹, and 3.7×10⁹ kWh of electricity would have to be added by 2020, 2025, and 2030, respectively, to meet the power requirements of EVs. This implies the emission of additional air pollutants (SO₂, NO_x, PM, and CO₂) through the increased electricity generation. Further analysis is thus required to explore the comprehensive impact of the adoption of EVs on pollutant emissions.

The sources of electricity for EVs generally include coal-based power and renewable energy. In this study, it is assumed that electricity from coal-based power was locally generated and renewable power means solar, wind power and imported electricity. Based on different power sources of EVs, there would be different impacts on the environment. Figure 5 shows the additional CO₂ and air pollutant (PM, SO₂, and NO_x) emissions that result from the additional electricity generation required by the adoption of EVs. As can be seen, an EVs policy promises to effectively mitigate CO and HC emissions in all the scenarios. However, it increases SO₂ and CO₂ emissions. Moreover, no reduction in NO_x and PM emissions would be achieved by the EVs policy if the power required by the EVs were entirely generated by coal-based power. Actually, in such a case, the total NO_x emission would increase by 0.72×10³ ton by 2020, 0.87×10³ ton by 2025, and 1.32×10³ ton by 2030. Conversely, if the power for EVs was generated 50% by coal-based and 50% from renewable sources, as in scenario S3, the EVs policy would reduce NO_x emission by 0.17×10³ ton by 2020, 0.21×10³ ton by 2025, and 0.24×10³ ton by 2030. If the EVs power was entirely generated from renewable sources, as in scenario S4, the policy would remarkably reduce NO_x emission by 1.07×10³ ton by 2020, 1.29×10³ ton by 2025, and 1.49×10³ by 2030. These results indicate that the vehicle-emissions could be reduced directly; but through the extra electricity generated from local region, leading to additional air-pollutants. Thus, an EVs policy should be enhanced by increasing the ratio of the needed power generated from renewable sources.

4.2. Optimized Robust Solutions For Energy and Environment Systems

4.2.1. Optimized Schemes of Energy Allocation

Figure 6 shows the optimized energy consumption of energy processing (heat processing, coke processing, and oil refining) and electricity generation (natural gas-based power and coal-based power) between period 1 and period 3. The use of natural gas is expected to substantially increase toward achieving a sustainably developed society. For example, the natural gas inputs to heat processing and natural gas-based power generation are projected to respectively increase by 47.16% (from 3.15 × 10⁹ m³ in 2020 to 4.63 × 10⁹ m³ in 2030) and 53.86% (from peak m 9.51 × 10⁹ m³ to 14.63 × 10⁹ m³) between 2020 and 2030. Crude oil input to oil refining, coal input to coke processing, and coal input to heating processing are expected gradually increase by 10.88%, 10.95%, and 3.86%, respectively, in periods 1, 2 and 3. Conversely, the coal consumption of coal-based power generation is expected to decrease by 2.69% between 2020 and 2030, specifically from 79.83 × 10⁶ ton to 77.68 × 10⁶ ton between period 1 and period 3. An appropriate energy mix can thus be effectively used to reduce coal consumption, and hence mitigate pollutant emissions.

4.2.2. Optimized schemes of electricity supply

Figure 7a shows the optimized electricity generation scheme between 2020 and 2030. The major electricity power conversion technologies in the BTH region include gas-based, coal-based, wind, and solar power. Optimized electricity generation is defined by the following equation: $OE{G}_{kth opt} = W_{kt opt}^{\pm } + Y_{kth opt}$, where $W_{kt\mathrm{opt}}^{\pm }$ and $Y_{kt\mathrm{opt}}$ are the amounts of electricity generated in the first and second stages. Coal-based power is expected to play the dominant role, contributing what would be 259.06 × 10⁹ kWh, 265.41 × 10⁹ kWh, and 250.97 × 10⁹ kWh in periods 1, 2 and 3 (h = 3), respectively. With the implementation of a series of energy mix policies for pollutant emission mitigation, such as the “Paris Agreement” and “Chinese Action Plan of Air Pollution Prevention and Control”, the ratio of coal-based power generation would decrease from 68.79% in 2020 to 55.48% in 2030. The utilization of natural gas-based and renewable energy would rapidly develop from 2020 to 2030. Natural gas-based power, wind power, solar power and generation are projected to increase by 53.86%, 87.96%, and 266.34% over the planning horizon, contributing 22.65%, 20.44%, and 1.44%, respectively to the total power generation by 2030.

Further, the amount of locally generated electricity is insufficient for the needs of the BTH region, with approximately 36% of the consumed electricity imported from other regions. Figure 7b describes imported electricity under different electricity demand levels between 2020 and 2030. Imported electricity is expected to increase from 202.47 × 10⁹ to 254.47 × 10⁹ kWh over the planning horizon under a high electricity demand level. This indicates that imported electricity is expected to play an increasingly important role from period 1 to period 3, especially at a high demand level.

4.2.3. Optimized schemes of energy processing

Figure 8 shows the amounts of processed secondary energy (i.e., coke, heat, diesel, gasoline, fuel oil, liquefied petroleum gas, kerosene, naphtha, and tar) during periods 1 to 3. With the rapid development of society over the entire considered time, the energy processing amounts for heat, gasoline, coke, kerosene, naphtha, fuel oil, tar, and liquefied petroleum gas (LPG) are expected to respectively increase by 10.37%, 29.98%, 10.59%, 19.97%, 1.11%, 7.48%, 6.90%, and 8.84%, between 2020 and 2030.

4.2.4. Pollutant Emissions from Energy Activities

Figure 9 shows the CO₂ and air pollutant (NO_x, SO₂, and PM) emissions from natural gas-based power generation, coal-based power generation, heat processing, coke processing, and oil refining) over the planning horizon. As can be observed, CO₂, NO_x, SO₂, and PM emissions are expected to gradually increase by 9.37%, 5.25%, 4.43%, and, 3.59%, respectively, between 2020 and 2030. The air-pollutants emissions would approach their peak values around 2030. For example, the SO₂ emission would increase by 4.14% between 2020 and 2025, and by 0.27% between 2025 and 2030, indicative of peaking around 2030. Coal-based power generation would be the major contributor of NO_x, SO₂, and PM emissions, respectively accounting for 36.66–33.90%, 41.89–39.03%, and 47.33–44.46% of these emissions between 2020 and 2030, with the specific contributions progressively decreasing. Coke processing would contribute 30.86–32.53%, 35.26%–37.46%, and 39.84%–42.67% of NO_x, SO₂, and PM emissions, respectively.

As indicated in Figure 9d, oil refining, coal-based power generation, and coke processing would be the major sources of CO₂ emissions between 2020 and 2030 in BTH region, accounting for 33.83–34.30%, 28.08–24.99%, and 23.64–23.98%, respectively. In conclusion, coal-based power generation and coke processing would be the major enablers of SO₂, NO_x, and PM emissions reduction, while oil refining, coal-based power generation, and coke processing would be the major contributors of CO₂ emission.

4.2.5. Analysis of system cost

The system costs are composed of 16 main components, namely, the costs of heat processing, coke processing, oil refining, coal input to coal-based power generation, natural gas input to natural gas-based power generation, imported electricity, first-stage of electricity generation, second-stage of electricity generation, electricity expansion, subsidy for solar power generation, subsidy for wind power generation, energy production, EVs charging piles, EVs charging stations, pollutants treatment, and risk recourse for the stochastic and fuzzy uncertainties. As shown in Figure 10, the costs of electricity generation, imported electricity, energy production, and coal input for coal-based power generation account for 30.51%, 16.90%, 16.67%, and 11.98% of the total system costs, respectively. Solar and wind power would be subsidized by the local government because their current price and technology limitations make them uncompetitive in the market. Government subsidies for wind and solar power would be RMB ¥ 40.99 × 10⁹, and RMB ¥ 6.52 × 10⁹, respectively.

5. Discussion

5.1. Analysis of Atochastic Uncertainties

In EES in the real-world, many parameters are expressed as random distribution (i.e., electricity demand, electricity generation, and electricity expansion). In this study, the TSP method was used to tackle the stochastic uncertainties. The decision variables of the electricity generation were divided into two subsets, namely, $W_{kt\mathrm{opt}}^{\pm }$ and $Y_{kt\mathrm{opt}}$, which are respectively the amounts of electricity generated in the first-stage (which must be predetermined) and the second-stage (obtained after the determination of the random variables) [55]. Generally, shortages may occur if the electricity demand-levels are continuously high, and the second-stage electricity generation amount $Y_{kt\mathrm{opt}}$ would be undertaken to avoid insufficient electricity supply.

In the TSP method, the final optimized electricity generation schemes were equal to $W_{kt}^{\pm }+Y_{kth}^{\pm }$. In detail, the optimized amount of electricity generated in the first-stage is given by $W_{kt}^{\pm }=W_{kt}^{-}+r{r}_{kt}·\Delta W$, where $r{r}_{kt}$ denotes the decision variables, $\Delta W=W_{kt}^{+}-W_{kt}^{-}$, and $r{r}_{kt}\in [0, 1]$. Take period 1 for example, ${\lambda }_{11opt}$ = 0.55, ${\lambda }_{21opt}$ = 1, ${\lambda }_{31opt}$ = 1.00, and ${\lambda }_{41opt}$ = 1.00, indicating that the optimized electricity generation amounts for coal-based power, natural gas-based power, solar power and wind power are supposed to be (221.08 + 0.52 × 58.74) × 10⁹, (58.41 + 1.00 ×5.00) × 10⁹ kWh, (1.78 + 1.00 × 0.50) × 10⁹ kWh, and (49.54 + 1.00 × 9.36) × 10⁹ kWh, respectively. Generally, the variation of $r{r}_{kt}$ represents diverse policies of electricity generation under stochastic electricity demands. When $r{r}_{kt}$ = 0, the cost would be relatively low, although a higher penalty may have to be paid when the generated electricity does not meet the demand. On the contrary, when $r{r}_{it}$ = 1, the cost would be higher, but accompanied by a lower risk of violating the target, and hence of incurring a lower penalty.

5.2. Analysis of Fuzzy Uncertainties

In the real world, energy price and energy demands often exhibit vagueness and ambiguity because of the subjectivity of human judgment [25]. According to FPP theory, the minimum, medium, and maximum values of these parameters are sufficient for expressing a triangular fuzzy parameter. And the proposed DRSFO-EES model can be used to effectively address the uncertainties expressed as triangular fuzzy parameter of the objectives and constraints. For FPP, the confidence level α is an indication of the manager’s violation risk attitude towards imprecise information [56,57]. In the present application of the DRSFO-EES model, four confidence levels (α = 0.5, 0.6, 0.8, and 0.9) were used to examine the impacts of different confidence levels on the EES. Generally, a higher confidence level implies a higher likelihood of satisfying the fuzzy confidence constraints, resulting in less uncertainty about the imprecise constraints. For instance, a confidence level of 0.8 indicates that the credibility of the constraint (e.g., $Cr\{M\tilde{A}{L}_{gt}|M{A}_{gt}\ge M\tilde{A}{L}_{gt}\}\ge \alpha$) is greater than or equal to 0.8. However, a higher confidence level increases the system costs. Contrarily, a lower confidence level implies a more aggressive attitude of the decision maker regarding the expected total system costs, and increases the uncertainty of the fuzzy constraints, resulting in a higher risk of violating the energy demand. Table 5 gives the optimized solution under different confidence levels α during the planning periods. As can be seen, a higher α increases the coefficients of the right-side constraints, further necessitating electricity import and increasing pollutant emission and the system costs. Vehicle ownership, heat processing, and coke processing also increase with α, all accompanied by pollutant emissions.

5.3. Risk Analysis

A robust optimization method can be effectively used to determine the associated risk from stochastic and fuzzy uncertainties. In this study, two risk recourse actions were adopted to make the model robust, which were used to capture the risks from stochastic and fuzzy uncertainties, respectively. Figure 11 presents the weighted values of the expected deviations from the stochastic uncertainties ($V\tilde{T}S=\rho ·{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{h=1}^H{p}_{th}}}}(V_{ith}+V{C}_{i{j}_{c}th}]$) and fuzzy uncertainties ($V\tilde{T}F=\rho ·{\displaystyle \sum _{t=1}^{3}V{F}_t}$) under different robustness levels of 0.2, 0.6, 0.8, and 1.0, respectively. ρ is a goal programming weight, through varying the ρ level, the decision makers can then control the variability of the recourse cost. Generally, a lower ρ corresponds to a lower weight value of the expected deviations and system costs, indicating an aggressive attitude of the manager regarding the system costs. However, this might be associated with a higher risk level because of the expected deviations from the uncertainties. On the contrary, a plan with a higherρwould better resist a deviation from the uncertainties of the EES. A decision with a higher robust level would thus correspond to a lower risk of system failure and higher system reliability. As results, the weighted values of the expected deviations from the stochastic and fuzzy uncertainties increased with increasing robustness level ρ. For instance, $V\tilde{T}S$ would be RMB 5.17 × 10⁹ when ρ = 0.2, RMB 15.54 × 10⁹ when ρ = 0.6, RMB 20.71 × 10⁹ when ρ = 0.8, and 25.89 × 10⁹ when ρ = 1. This analysis demonstrates the trade-offs between system costs and reliability. The results enable the manager to plan with a reasonable consideration of both system costs and reliability.

6. Conclusions

In this study, a DRSFO-EES model was developed for planning an EES while considering the traffic sector, which integrates the TSP, FPP, RTSO, and RFPP methods in a single framework for the effective handling of EES uncertainties expressed as fuzzy sets and stochastic uncertainties as well as their combinations, capturing of the associated risks from the stochastic and fuzzy uncertainties, and analyzing the trade-offs between system costs and reliability. Four confidence levels (α = 0.5, α = 0.6, α = 0.8, and α = 0.9) and four robust levels (ρ = 0.2, ρ = 0.6, ρ = 0.8, and ρ = 1.0) were used for examining the impacts of uncertainties on the objective function, constraints and optimized solutions of the DRSFO-EES model.

The proposed model was applied to the EES of the BTH region in China. Following is a summary of the findings and the identified policy implications:

(1)

Limiting the numbers of LDVs and HDTs could effectively reduce vehicular emissions. LDVs are expected to be the major contributors of CO and HC emissions, and HDTs are expected to be the major contributors of NO_x and PM emissions.

(2)

A EVs policy would be enhanced by increasing the ratio of power generated for EVs from renewable sources. The emission reduction effect of an EVs policy would thus be limited, especially with regard to NO_x and PM emissions, if the EVs power source was entirely coal-based.

(3)

Optimizing the energy mix and developing the renewable energy can effectively reduce air-pollutant and CO₂ emissions. Air-pollutant amounts of NO_x, SO₂, and PM emissions in the BTH region are expected to peak around 2030, because the energy mix of the study region would be transformed from one dominated by coal to one with a cleaner pattern, with vigorous development of the utilization of natural gas and renewable energy.

(4)

Enhancement of the energy utilization efficiencies of coal-based power generation, oil refining, and coke processing would effectively reduce CO₂ and air-pollutant emissions. Coal-based power generation and coke processing are expected to be the major contributors of air-pollutant missions, while oil refining, coal-based power generation and coke processing would be the chief sources of CO₂ emissions.

Although the DRSFO method was the first attempt for planning an EES while considering the traffic sector of the BTH region, results indicated that DRSFO-EES could: (i) explore the impacts of different vehicle policies (i.e., EVs deployment, EVs power source for EVs, and vehicular emission standards) on vehicular emissions; (ii) generate robust optimized solutions for energy allocation, oil refining, coking processing, heat processing, electricity generation and expansion, electricity importation, as well as emission mitigation under multiple uncertainties; (iii) identify the atmospheric pollution contributions of different energy activities such as coke processing, electricity generation, heat processing, oil refining, and motor vehicle operation. The proposed model could help to balance the contradiction between increasing energy demands and an increasing vehicle population, “high coal” energy systems, and the pressures of emission mitigation. Moreover, the proposed model could be applied at both city and regional scales, which would support policymakers adjusting current energy and environmental strategies in sustainable and robust ways.

However, the DRSFO-EES also has potential limitations and extensions should be addressed in future study. Firstly, the developed model, based on historical data of annual electricity demand to predict future electricity demands, does not consider the specific parameters such as hourly or seasonal electricity load curves, which may result in significant deviations from the optimized decision schemes, even in the event of electricity shortage. Thus, further study is required for considering the hourly and seasonal electricity load curves; secondly, the TSP, FPP, RTSO, and RFPP methods were combined into a single framework to formulate the DRSFO-EES model, leading to relatively high computational requirements. As a result, simplifying the calculation procedure could be required in the further study work; thirdly, the DRSFO-EES model mainly focused on economic objectives, whilst scarcely considering the trade-off between economic and environmental objectives. Therefore, further study should make improvements in the handling of multi-objective problems and better balance the tension between energy and environmental systems.