A Stochastic Inexact Robust Model for Regional Energy System Management and Emission Reduction Potential Analysis—A Case Study of Zibo City, China

Abstract

In this study, in order to improve regional energy system adjustment, a multistage stochastic inexact robust programming (MSIRP) is proposed for electric-power generation planning and structure adjustment management under uncertainty. Scenario-based inexact multistage stochastic programming and stochastic robust optimization were integrated into general programming to reflect uncertainties that were expressed as interval values and probability distributions in the objective function and constraints. An MSIRP-based energy system optimization model is proposed for electric-power structure management of Zibo City in Shandong Province, China. Three power demand scenarios associated with electric-power structure adjustment, imported electricity, and emission reduction were designed to obtain multiple decision schemes for supporting regional sustainable energy system development. The power generation schemes, imported electricity, and emissions of CO₂ and air pollutants were analyzed. The results indicated that the model can effectively not only provide a more stable energy supply strategies and electric-power structure adjustment schemes, but also improve the balanced development between conventional and new clear power generation technologies under uncertainty.

1. Introduction

Rapid power consumption increment, increasing deterioration of environmental quality, and imperfect energy system management have led to unsustainable energy resources exploitation and utilization, unreasonable electric-power structure, and serious environmental issues [1,2,3]. In order to search effective and suitable energy development strategies for regional condition, energy system management and planning has become a priority for many countries and regions. However, multiple forms of uncertain information are involved in energy system management and the related social-economic factors and/or technical-economic parameters, causing a variety of complexities in decision support and policy analysis for regional energy planning [4]. In addition, such complexities would pose great challenges in formulating more scientific and reasonable development strategies for decision-makers, and have serious impact on the effectiveness of energy supply schemes. Therefore, it is desirable to develop effective uncertain optimization models/techniques for energy system management and the related decision analysis.

Previously, a great number of inexact programming approaches were proposed for helping energy system planning and management in different regional scales [5,6,7,8,9,10,11,12,13,14]. For example, Cai et al. (2009) advanced an interval parameter interactive decision support system for energy system management under reflecting uncertainties as interval values [15]. Li et al. (2010) proposed an inexact fuzzy multistage stochastic energy system management model for supporting regional electric-power generation and capacity planning, where interval parameter programming, mixed integer linear programming, multistage stochastic programming, and fuzzy linear programming were incorporated into a general optimization framework [16]. Li et al. (2011) proposed a fuzzy stochastic energy system optimization model associated with renewable energy development and greenhouse gas mitigation, where the uncertainties in the objective and constraints were expressed as fuzzy interval functions, interval values, and discrete probability distributions [17]. Huang et al. (2017) developed an inexact fuzzy stochastic chance-constrained programming for evacuation management of nuclear power plant, where interval parameter programming and fuzzy stochastic chance-constrained programming were integrated into a general framework for dealing with uncertainties [18]. Sheikhahmadi et al. (2018) proposed a risk-based two-stage stochastic programming for microgrid system operation management, where two-stage stochastic programming was to reflect uncertainties of renewable energy, and conditional value at risk index was used to avoid the system risk [19].

Among these methods, scenario-based interval multistage stochastic programming, as a hybrid method of interval parameter programming and scenario-based multistage stochastic programming, could deal with uncertainties presented as interval numbers and random distributions, and have been widely applied in energy system management [20,21,22]. For example, Xie et al. (2010) advanced an inexact fixed-mix multistage stochastic programming for long-term greenhouse gas emission reduction management in a regional scale energy system, where the fixed probability multistage stochastic programming and interval-parameter programming were integrated for expressing uncertainties in energy system management problems [23]. Wu et al. (2015) proposed an integrated method with interval-parameter programming, chance-constraint programming, and multistage stochastic programming for the coupled biomass–municipal solid waste power system operation management, which could reflect uncertainties as interval information and random distributions over a multistage context [24]. Golari et al. (2016) presented a production-inventory planning model in a multi-plant manufacturing system powered with onsite and grid renewable energy, where multistage stochastic programming was used to reflect system dynamic and uncertainties [25]. Fu et al. (2017) advanced an interval multistage fuzzy-stochastic programming for regional electric-power system management under considering environmental quality constraints, where interval-parameter programming, multistage stochastic programming, and fuzzy probability distribution was integrated to reflect the uncertain information and dynamic variation in the energy system [26]. Wang et al. (2018) developed multistage joint-probabilistic left-hand-side chance-constrained fractional programming for electric-power system planning considering climate change mitigation [27].

Although scenario-based inexact multistage stochastic programming had been successfully applied in many fields, it could not directly and effectively avoid the risk of random events, and the limitations would pose threats to system stability. Based on this point, stochastic robust optimization (SRO) is proposed for solving the problems through introducing the risk-aversion attitude into optimization models and obtaining robust solutions for stochastic system management [28,29,30]. The methods that coupled with the scenario-based inexact multistage stochastic programming and SRO have been used in solving many energy and environmental problems, such as water resources allocation, electric-power generation, and water/air quality management. For example, Chen et al. (2013) developed an interval robust-optimization model for CO₂ emission reduction management in energy systems, where the robustness measures were introduced to examine whether the second-stage cost variability could meet the expected levels or not [31]. Xie et al. (2014) proposed an inexact stochastic risk-aversion model for electric-power structure adjustment and pollutant emission management, where interval-parameter programming, stochastic robust optimization, and multistage stochastic programming were integrated to address system uncertainties [32].

Therefore, the aim of this study is to formulate a multistage stochastic inexact robust programming (MSIRP) model to support regional electric-power system management coupled with pollutant mitigation constraints and power structure adjustment requirements in Zibo City, China. The method could not only reflect multiple uncertainties expressed as interval values and probability distribution, but also make a tradeoff between system risk and cost according to the decision-makers’ attitudes. The modeling results will be helpful for local decision-makers to choose cost–risk electric-power generation schemes, and obtain reasonable electric-power structure adjustment strategies. The rest structure organization of this paper is provided as follows. The development process and solution algorithm of multistage stochastic inexact robust programming (MSIRP) is introduced in Section 2. The overview of the energy system of Zibo City are described, and a MSIRP-based energy structure adjustment model is proposed in Section 3. The obtained results and deep discovery of the case study are analyzed and discussed in Section 4. The main conclusions are presented in Section 5.

2. Methodology

In regional energy systems, dynamic characters, discrete probability distributions, intervals information, and policy implications were addressed through scenario-based inexact multistage stochastic programming, and SRO could effectively handle the system risk. The modeling framework of the MSIRP could obtain applicable and reasonable solutions under different random scenarios corresponding to power generation targets for decision-makers in order to support the energy system development in the future.

2.1. Inexact Scenario-Based Multistage Stochastic Programming

In the scenario-based multistage stochastic programming, the probabilities $p_{tk}$($t=1,2,\dots ,T;$ $k=1,2,\dots ,K_t,$) of the stochastic event have predefined values, and the parameters without probability can be reflected as interval values. Thus, the scenario-based inexact multistage stochastic programming can be expressed as follows [33,34]:

$$Min\hspace{0.17em}{f}^{\pm }={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{j=1}^{n_1}{c}_{jt}^{\pm }{x}_{jt}^{\pm }+{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{j=1}^{n_1}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{d}_{jt}^{\pm }{y}_{jtk}^{\pm }}}}}},$$

(1)

subject to

$${\displaystyle \sum _{j=1}^{n_1}{a}_{rjt}^{\pm }{x}_{jt}^{\pm }\le b_{rt}^{\pm }\hspace{0.17em},\forall r,t}$$

(2)

$${\displaystyle \sum _{j=1}^{n_1}{a}_{ijt}^{\pm }{x}_{jt}^{\pm }+{\displaystyle \sum _{j=1}^{n_1}{e}_{ijt}^{\pm }}}{y}_{jkt}^{\pm }\le {\tilde{w}}_{itk}^{\pm }\hspace{0.17em},\forall i,t,k$$

(3)

$$x_{jt}^{\pm }\ge 0,\forall t,j=1,2,\dots ,n_1$$

(4)

$$y_{jkt}^{\pm }\ge 0,\forall t,k,j=1,2,\dots ,n_1$$

(5)

where $p_{tk}$ is the probability for scenario k in period t; for each period t, the total number of scenarios is denoted as $K_t$, and ${\displaystyle \sum _{k=1}^{K_t}{p}_{tk}=1}$; and ${\tilde{w}}_{ikt}^{\pm }$ represents the random parameter in the model associated with the occurrence probability $p_{tk}$ in period t. $x_{jt}^{\pm }$ denotes the first-stage variables that have to be determined before the random event occurrence; and $y_{jkt}^{\pm }$ are the second-stage variables that have to be decided for making a recourse actions to fulfil validity of the decision-making after the random event occurrence.

2.2. Inexact Multistage Stochastic Robust Programming

The proposed inexact scenario-based multistage stochastic programming can effectively reflect stochastic information, interval values, and dynamic feature by means of discrete random variables in long-term planning problems. However, Model (1) could not effectively reflect the system risk introduced by random information, that directly affect the feasibility and reliability of the proposed model. SRO is an effective choice for solving such problems, and it can be introduced into Model (1), that leads to a multistage stochastic inexact robust programming (MSIRP) as follows [32]:

$$\begin{array}{cc}\hfill Min\hspace{0.17em}{f}^{\pm }& ={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{j=1}^{n_1}{c}_{jt}^{\pm }{x}_{jt}^{\pm }+{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{j=1}^{n_2}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{d}_{jt}^{\pm }{y}_{jtk}^{\pm }}}}}}\hfill \\ & +\omega {\displaystyle \sum _{t=1}^T{\displaystyle \sum _{j=1}^{n_2}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}\left|d_{jt}^{\pm }{y}_{jtk}^{\pm }-{\displaystyle \sum _{j=1}^{n_2}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{d}_{jt}^{\pm }{y}_{jtk}^{\pm }}}\right|}}}\end{array} ,$$

(6)

subject to

$${\displaystyle \sum _{j=1}^{n_1}{a}_{rjt}^{\pm }{x}_{jt}^{\pm }\le b_{rt}^{\pm }\hspace{0.17em},\forall r,t}$$

(7)

$${\displaystyle \sum _{j=1}^{n_1}{a}_{ijt}^{\pm }{x}_{jt}^{\pm }+{\displaystyle \sum _{j=1}^{n_1}{e}_{ijt}^{\pm }}}{y}_{jkt}^{\pm }\le {\tilde{w}}_{itk}^{\pm }\hspace{0.17em},\forall i,t,k$$

(8)

$$x_{jt}^{\pm }\ge 0,\forall t,j=1,2,\dots ,n_1$$

(9)

$$y_{jkt}^{\pm }\ge 0,\forall t,k,j=1,2,\dots ,n_1$$

(10)

where the non-negative factor $\omega$ denotes a trade-off weight coefficient; and $\left|d_{jt}^{\pm }{y}_{jtk}^{\pm }-{\displaystyle \sum _{j=1}^{n_2}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{d}_{jt}^{\pm }{y}_{jtk}^{\pm }}}\right|$ is a variability measure for reflecting the multistage recourse costs. The objective of Model (6) is a nonlinear function, and according to [35,36], the model can be converted into a linear programming model as follows:

$$\begin{array}{cc}\hfill Min\hspace{0.17em}{f}^{\pm }& ={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{j=1}^{n_1}{c}_{jt}^{\pm }{x}_{jt}^{\pm }+{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{j=1}^{n_1}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{d}_{jt}^{\pm }{y}_{jtk}^{\pm }}}}}}\hfill \\ & +\omega {\displaystyle \sum _{t=1}^T{\displaystyle \sum _{j=1}^{n_1}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}(d_{jt}^{\pm }{y}_{jtk}^{\pm }-{\displaystyle \sum _{j=1}^{n_1}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{d}_{jt}^{\pm }{y}_{jtk}^{\pm }+2{\theta }_{jkt}^{\pm }}})}}}\end{array},$$

(11)

subject to

$$d_{jt}^{\pm }{y}_{jtk}^{\pm }-{\displaystyle \sum _{j=1}^{n1}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{d}_{jt}^{\pm }{y}_{jtk}^{\pm }+{\theta }_{jkt}^{\pm }\ge 0,\forall k,j=1,2,\dots ,n_1}}$$

(12)

$${\displaystyle \sum _{j=1}^{n_1}{a}_{rjt}^{\pm }{x}_{jt}^{\pm }\le b_{rt}^{\pm }\hspace{0.17em},\forall r,t}$$

(13)

$${\displaystyle \sum _{j=1}^{n_1}{a}_{ijt}^{\pm }{x}_{jt}^{\pm }+{\displaystyle \sum _{j=1}^{n_1}{e}_{ijt}^{\pm }}}{y}_{jkt}^{\pm }\le {\tilde{w}}_{itk}^{\pm }\hspace{0.17em},\forall i,t,k$$

(14)

$$x_{jt}^{\pm }\ge 0,\forall t,j=1,2,\dots ,n_1$$

(15)

$$y_{jkt}^{\pm }\ge 0,\forall t,k,j=1,2,\dots ,n_1$$

(16)

where, through introducing the slack variable ${\theta }_{jkt}^{\pm }$, the objective can be transferred into a linear function as well as generate a specific control constraint (12). For Model (11), the first-stage variables $x_{jt}^{\pm }$ are considered/inputted as interval values with the lower and upper bound, and this cannot be directly solved using the existing methods. In this study, let ${\mu }_{jt}$ be a decision variable, ${\mu }_{jt}\in [0,1];$ $\Delta x_{jt}=x_{jt}^{+}-x_{jt}^{-}$, the first-stage variable $x_{jt}=x_{jt}^{-}+{\mu }_{jt}\Delta x_{jt}$, and ${\mu }_{jt}$ are intermediate decision variables for obtaining an optimized target values of the first-stage to support the related policy analyses [32]. According to [37], the MSIRP model can be transformed into two linear submodels, and the submodel corresponding to $f^{-}$ can be firstly transformed as follows (assume that $c_{jt}^{\pm }\ge 0,$ ${\widehat{w}}_{itk}^{+}>0,$ $b_{rt}^{\pm }>0,$ and $f^{\pm }>0$):

$$\begin{array}{cc}\hfill Min\hspace{0.17em}{f}^{-}& ={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{j=1}^{n_1}{c}_{jt}^{-}(x_{jt}^{-}+{\mu }_{jt}\Delta x_{jt})}}\hspace{0.17em}+{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}({\displaystyle \sum _{j=1}^{j_1}{d}_{jt}^{-}{y}_{jtk}^{-}+{\displaystyle \sum _{j=j_1+1}^{n_1}{d}_{jt}^{-}{y}_{jtk}^{+}}})}}\hfill \\ & +\omega {\displaystyle \sum _{t=1}^T{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}[{\displaystyle \sum _{j=1}^{j_1}(d_{jt}^{-}{y}_{jtk}^{-}-{\displaystyle \sum _{j=1}^{j_1}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{d}_{jt}^{-}{y}_{jtk}^{-}+2{\theta }_{jkt}^{-}}})}}}]\hfill \\ & +\omega {\displaystyle \sum _{t=1}^T{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}[{\displaystyle \sum _{j=j_1+1}^{j_1}(d_{jt}^{-}{y}_{jtk}^{+}-{\displaystyle \sum _{j=j_1+1}^{n_1}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{d}_{jt}^{-}{y}_{jtk}^{+}+2{\theta }_{jkt}^{-}}})}}}]\hfill \end{array},$$

(17)

subject to

$$\begin{array}{l}{\displaystyle \sum _{j=1}^{j_1}{\left|d_{jtk}\right|}^{-}Sign(d_{jtk}^{-})y_{jkt}^{-}+{\displaystyle \sum _{j=j_1+1}^{n_1}{\left|d_{jtk}\right|}^{-}Sign(d_{jtk}^{-})y_{jkt}^{+}}}\\ -{\displaystyle \sum _{j=1}^{j_1}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{\left|d_{jtk}\right|}^{-}Sign(d_{jtk}^{-})y_{jkt}^{-}}}-{\displaystyle \sum _{j=1}^{n_1}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{\left|d_{jtk}\right|}^{-}Sign(d_{jtk}^{-})y_{jkt}^{+}+{\theta }_{jtk}^{-}\ge 0,\forall i,j}}\end{array}$$

(18)

$${\displaystyle \sum _{j=1}^{n_1}{\left|a_{rjt}\right|}^{+}Sign(a_{rjt}^{+})(x_{jt}^{-}+{\mu }_{jt}\Delta x_{jt})\le b_{rt}^{-},\forall r,t}$$

(19)

$$\begin{array}{l}{\displaystyle \sum _{j=1}^{n_1}{\left|a_{ijt}\right|}^{+}Sign(a_{rjt}^{+})(x_{jt}^{-}+{\mu }_{jt}\Delta x_{jt})}+{\displaystyle \sum _{j=1}^{j_1}{\left|e_{ijt}\right|}^{+}Sign(e_{ijt}^{+})}{y}_{jkt}^{-}\\ +{\displaystyle \sum _{j=j_1+1}^{n_1}{\left|e_{ijt}\right|}^{-}Sign(e_{ijt}^{-})y_{jkt}^{+}\le {\tilde{w}}_{itk}^{-}},\forall i,t,k\end{array}$$

(20)

$$x_{jt}^{-}+{\mu }_{jt}\Delta x_{jt}\ge 0,\forall j,t$$

(21)

$$0\le {\mu }_{jt}\le 1,\forall j,t$$

(22)

$$y_{jkt}^{-}\ge 0,\forall t,k,j=1,2,\dots ,j_1$$

(23)

$$y_{jkt}^{+}\ge 0,\forall t,k,j=j_1+1,j_1+2,\dots ,n_1$$

(24)

where ${\mu }_{jt}$, $y_{jkt}^{-}$($j=1,2,\dots ,j_1$) and $y_{jkt}^{+}$($j=j_1+1,j_1+2,\dots ,n_1$) are the decision variables of model (17); $y_{jkt}^{-}$($j=1,2,\dots ,j_1$) and $y_{jkt}^{+}$($j=j_1+1,j_1+2,\dots ,n_1$) are the second-stage decision variables with positive and negative coefficients in the objective function; and the optimized solution of the first-stage variables are $x_{jtopt}=x_{jt}^{-}+{\mu }_{jtopt}\Delta x_{jt}$. Then, the submodel corresponding to $f^{+}$ can be expressed as follows:

$$\begin{array}{cc}\hfill Min\hspace{0.17em}{f}^{+}& ={\displaystyle \sum _{t=1}^T{\displaystyle \sum _{j=1}^{n_1}{c}_{jt}^{+}{x}_{jtopt}}}+{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}({\displaystyle \sum _{j=1}^{j_1}{d}_{jt}^{+}{y}_{jtk}^{+}+{\displaystyle \sum _{j=j_1+1}^{n_1}{d}_{jt}^{+}{y}_{jtk}^{-}}})}}\hfill \\ & +\omega {\displaystyle \sum _{t=1}^T{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}[{\displaystyle \sum _{j=1}^{j_1}(d_{jt}^{+}{y}_{jtk}^{+}-{\displaystyle \sum _{j=1}^{n_2}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{d}_{jt}^{+}{y}_{jtk}^{+}+2{\theta }_{jkt}^{+}}})}}}]\hfill \\ & +\omega {\displaystyle \sum _{t=1}^T{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}[{\displaystyle \sum _{j=j_1+1}^{n_1}(d_{jt}^{+}{y}_{jtk}^{-}-{\displaystyle \sum _{j=j_1+1}^{n_1}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{d}_{jt}^{+}{y}_{jtk}^{-}+2{\theta }_{jkt}^{+}}})}}}]\hfill \end{array},$$

(25)

subject to

$$\begin{array}{l}{\displaystyle \sum _{j=1}^{j_1}{\left|d_{jtk}\right|}^{+}Sign(d_{jtk}^{+})y_{jkt}^{+}+{\displaystyle \sum _{j=j_1+1}^{n_1}{\left|d_{jtk}\right|}^{+}Sign(d_{jtk}^{+})y_{jkt}^{-}}}\\ -{\displaystyle \sum _{j=1}^{j_1}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{\left|d_{jtk}\right|}^{+}Sign(d_{jtk}^{+})y_{jkt}^{+}}}-{\displaystyle \sum _{j=j_1+1}^{n_1}{\displaystyle \sum _{k=1}^{K_t}{p}_{tk}{\left|d_{jtk}\right|}^{+}Sign(d_{jtk}^{+})y_{jkt}^{-}+{\theta }_{jtk}^{+}\ge 0,\forall i,j}}\end{array}$$

(26)

$${\displaystyle \sum _{j=1}^{n_1}{\left|a_{rjt}\right|}^{-}Sign(a_{rjt}^{-})x_{jtopt}\le b_{rt}^{+},\forall r,t}$$

(27)

$${\displaystyle \sum _{j=1}^{n_1}{\left|a_{ijt}\right|}^{-}Sign(a_{rjt}^{-})\Delta x_{jtopt}}+{\displaystyle \sum _{j=1}^{j_1}{\left|e_{ijt}\right|}^{-}Sign(e_{ijt}^{-})}{y}_{jkt}^{+}+{\displaystyle \sum _{j=j_1+1}^{n_1}{\left|e_{ijt}\right|}^{+}Sign(e_{ijt}^{+})y_{jkt}^{-}\le {\tilde{w}}_{itk}^{+}},\forall i,t,k$$

(28)

$$y_{jkt}^{+}\ge y_{jktopt}^{-},\forall t,k,j=1,2,\dots ,j_1$$

(29)

$$y_{jktopt}^{+}\ge y_{jkt}^{-}\ge 0,\forall t,k,j=j_1+1,j_1+2,\dots ,n_1$$

(30)

where $y_{jkt}^{+}$($j=1,2,\dots ,j_1$) and $y_{jkt}^{-}$($j=j_1+1,j_1+2,\dots ,n_1$) are decision variables that can be obtained through solving Submodel (25). Thus, the optimal solutions of Model (11) can be expressed as follows:

$$x_{jtopt}=x_{jt}^{-}+{\mu }_{jtopt}\Delta x_{jt}$$

(31)

$$y_{jktopt}^{\pm }=[y_{jktopt}^{-},y_{jktopt}^{+}]$$

(32)

$$f_{opt}^{\pm }=[f_{opt}^{-},f_{opt}^{+}].$$

(33)

3. Case Study

3.1. Overview of Energy System in Zibo City

Zibo City (35°55′20″~37°17′14″ N, 117°32′15″~118°31′00″ E), as shown in Figure 1, is located in the middle of Shandong province, China. Zibo City governs Zhangdian district, Zichuan district, Boshan district, Zhoucun district, Linzi district, Huantai country, Gaoqing country, and Yiyuan country, with a total area of 5938 km² and a total population of 4.61 million in 2014 [38]. In Zibo City, the manufacturing industry plays a significant role in supporting regional economic development; especially the ceramics manufacturing industry is famous around the world. For example, in 2014, the income of ceramic industry reached 112.8 billion yuan. In addition, high-new-technology industries (e.g., new materials, fine chemicals, and biological medicines) and other traditional industries (e.g., petrochemical industry, pharmaceuticals, metallurgy, and machinery and textiles) are developing rapidly in recent years. Moreover, in 2014, gross agricultural product reached to RMB 25.22 billion yuan, and the tertiary industry increased by RMB 163.45 billion yuan compared with 2013. In general, the rapid social-economic development is closely related with a higher power consumption. According to regional energy system statistic data in recent years, local electric-power generation is far from satisfying increasing regional demands.

Generally, the main electricity generation in Zibo City mainly relies on coal-fired power. The cogeneration power plays a large proportion in all electricity generation in Zibo, which could not only meet the demand of the district heating, but also greatly improve efficiency of coal resource utilization. In order to meet environmental requirements, there would have to be total consumption control on coal resources, according to the regional development plan from Zibo Municipal Development and Reform Committee. In addition, Zibo is abundant in renewable energy resources, such as solar, biomass, and wind, that have been considered as the primary options for addressing the crisis of electric-power shortage, and air pollutant and greenhouse gas mitigation. For instance, the average annual sunshine time reaches up to 2542.6 h with a greater potential and space for solar power and heat utilization. Moreover, throughout the windy corridor, in the surrounding of Boshan District and the southern mountain areas of Zichuan District, Zibo possesses the excellent conditions to build wind farms. According to regional energy development strategy of Zibo City (2010–2020), a greater number of renewable energy development plans have been promoted for adjusting the existing electric-power system structure, including 114 MW, 50 MW, and 244.5 MW of biomass and garbage power, solar, and wind power generation capacity by 2015, respectively. As a result, it will be helpful for alleviating the contradiction between energy supply capacity and consumption demands, and reducing atmospheric pollutants and carbon emission.

Although renewable energy has achieved development, and the government has also made great efforts to change regional electric-power structure, it still faces many challenges in electric-power system management. As a result of regional economic development, urbanization advance, and population growth, electric-power consumption and environmental quality requirement would be increasingly prominent, leading to an urgent need for regional electric-power structure adjustment. In this study, an inexact regional electric-power system optimization model is developed through multistage stochastic inexact robust programming for solving the following questions: (1) how to develop electric-power generation schemes for different power conversion technologies under air pollution and carbon mitigation requirements; (2) how to plan the overall development of renewable power conversion technologies and the proportion of imported electricity; (3) how to formulate more reasonable decision alternatives for decision-makers under different trade-offs between system cost and risk.

3.2. Electric-Power System Optimization Model Formulation

The developed multistage stochastic inexact robust programming is considered for regional electric-power system management in Zibo City. The objective is to achieve the optimal plans of electric-power supply with minimized system costs. The renewable power generation development, capacity expansion, and air pollutant and carbon emission reduction were also considered. Thus, the optimized model can be developed as follows:

$$Min f^{\pm }=f_1^{\pm }+f_2^{\pm }+f_3^{\pm }+f_4^{\pm }+f_5^{\pm }+f_6^{\pm }-f_7^{\pm }+f_8^{\pm }$$

(34)

[Costs for energy resources consumption]

$$f_1^{\pm }={\displaystyle \sum _{p=1}^P{\displaystyle \sum _{t=1}^{T}PE{C}_{pt}^{\pm }}}\cdot (A{E}_{pt}^{\pm }+p_{th}\cdot D{E}_{pth}^{\pm })\cdot E{F}_{pt}^{\pm }$$

(35)

[Costs for power generation]

$$f_2^{\pm }= {\displaystyle \sum _{p=1}^P{\displaystyle \sum _{t=1}^{T}P{V}_{pt}^{\pm }\cdot A{E}_{pt}^{\pm }}}+{\displaystyle \sum _{p=1}^P{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{h=1}^H{p}_{th}\cdot \left(P{V}_{pt}^{\pm }+P{P}_{pt}^{\pm }\right)\cdot D{E}_{pth}^{\pm }}}}$$

(36)

[Cost for the district heating]

$$f_3^{\pm }={\displaystyle \sum _{p=1}^P{\displaystyle \sum _{t=1}^{T}C{V}_{pt}^{\pm }\cdot (A{H}_{pt}^{\pm }+D{H}_{pt}^{\pm })}}$$

(37)

[Costs for the expansion of installed capacity]

$$f_4^{\pm }={\displaystyle \sum _{p=1}^P{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{h=1}^H{p}_{th}\cdot \left(YE{H}_{pth}^{\pm }\cdot A_{pt}^{\pm }+XE{H}_{pth}^{\pm }\cdot B_{pt}^{\pm }\right)}}}$$

(38)

[Costs for atmospheric pollutants treatment]

$$f_5^{\pm }={\displaystyle \sum _{i=1}^I{\displaystyle \sum _{p=1}^P{\displaystyle \sum _{t=1}^{T}A{E}_{pt}^{\pm }\cdot {\xi }_{ipt}^{\pm }\cdot \left(1-{\eta }_{ipt}^{\pm }\right)\cdot CP{C}_{it}^{\pm }}}}+{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{p=1}^P{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{h=1}^H{p}_{th}\cdot D{E}_{pth}^{\pm }\cdot {\xi }_{ipt}^{\pm }\cdot \left(1-{\eta }_{ipt}^{\pm }\right)\cdot DP{C}_{it}^{\pm }}}}}$$

(39)

[Costs for imported electric power]

$$f_6^{\pm }={\displaystyle \sum _{t=1}^T{p}_{th}\cdot I{E}_{th}^{\pm }\cdot IP{E}_t^{\pm }}$$

(40)

[Subsidies for renewable energy generation]

$$f_7^{\pm }={\displaystyle \sum _{p=3}^P{\displaystyle \sum _{t=1}^T\left(A{E}_{pt}^{\pm }+p_{th}\cdot D{E}_{pth}^{\pm }\right)\cdot S{U}_{pt}^{\pm }}}$$

(41)

[Robust function]

$$f_8^{\pm }=\lambda {\displaystyle \sum _{p=1}^P{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{h=1}^H{p}_{th}[{\epsilon }_{pth}^{\pm }-{\displaystyle \sum _{p=1}^P{\displaystyle \sum _{h=1}^H{p}_{th}\cdot {\epsilon }_{pth}^{\pm }+2{\theta }_{pth}^{\pm }}}]}}}$$

(42)

where,

$$\begin{array}{cc}\hfill {\epsilon }_{ pth}^{\pm }& =PE{C}_{pt}^{\pm }\cdot D{E}_{pt}^{\pm }\cdot E{F}_{pt}^{\pm }+(P{V}_{pt}^{\pm }+P{P}_{pt}^{\pm })\cdot D{E}_{pth}^{\pm }\hfill \\ & +I{E}_t^{\pm }\cdot IP{E}_t^{\pm }+(YE{H}_{pth}^{\pm }\cdot A_{pt}^{\pm }+XE{H}_{pth}^{\pm }\cdot B_{pt}^{\pm })\hfill \\ & +{\displaystyle \sum _{i=1}^{I}D{E}_{pth}^{\pm }\cdot {\xi }_{ipt}^{\pm }\cdot \left(1-{\eta }_{ipt}^{\pm }\right)\cdot DP{C}_{it}^{\pm }}-D{E}_{pth}^{\pm }\cdot S{U}_{pt}^{\pm }\hfill \end{array}$$

(43)

where $f^{\pm }$ is the objective of the proposed model (million yuan ¥); $p$ is the power conversion technologies, $p$ = 1, 2, 3, 4, and 5 for combined heat and power (CHP), hydroelectric power, solar photovoltaic power, wind power, and garbage power and biomass power, respectively; $i$ denotes different atmospheric pollutants, $i$ = 1, 2, 3, 4 for CO₂, SO₂, NOx, and particulate matter, respectively; t is the planning period; h denotes the electric-power demand level, $h$ = 1 for low level, h = 2 for medium level, and h = 3 for high level, respectively. $Z_{pt}^{\pm }$ is the amount of energy resource consumption for power conversion technology $p$ (PJ); $PE{C}_{pt}^{\pm }$ represents the energy price for technology $p$ (million ¥/PJ); $P{V}_{pt}^{\pm }$ and $P{P}_{pt}^{\pm }$ are the variable cost for power generation and the penalty cost of excess power generation of technology $p$ (million ¥/GWh); $A{E}_{pt}^{\pm }$ denotes the pre-regular electric-power generation by technology $p$ (GWh); $D{E}_{pth}^{\pm }$ is the excess power generation by technology $p$ under different electric-power deficiency levels $h$ (GWh); $C{V}_{pt}^{\pm }$ represents the variable cost for heat generation by technology $p$ (million ¥/PJ); $A{H}_{pt}^{\pm }$ is the amount of district heat supply by technology $p$ (PJ); $D{H}_{pt}^{\pm }$ denotes the amount of district heat supply by expanded capacity $XE{H}_{pth}^{\pm }$ (PJ); $A_{pt}^{\pm }$ and $B_{pt}^{\pm }$ are the fixed-charge cost and variable cost for capacity expansion of technology $p$ (million ¥); $S{U}_{pt}^{\pm }$ is the subsidy for new renewable energy generation $p$ (million ¥/GW); $YE{H}_{pth}^{\pm }$ represents the binary variable for determining the capacity choice of technology $p$ expansion (0 denotes no expansion; 1 represents expansion); $XE{H}_{pth}^{\pm }$ is the capacity expansion amount for technology $p$ under different electric-power deficiency levels $h$ (GW); $I{E}_{th}^{\pm }$ denotes imported power amount (GWh); $IP{E}_t^{\pm }$ is the cost of imported power (million ¥/GWh); $CP{C}_{it}^{\pm }$ and $DP{C}_{it}^{\pm }$ are the removal cost of pollutant $i$ treatment and the penalty cost of excess pollutant $i$ treatment for technology $p$ (million¥/ton); ${\xi }_{ipt}^{\pm }$ is the generation rate of pollutant from technology (ton/GWh).

Constraint:

[Constraints for electric-power supply and demand balance]

$${\displaystyle \sum _{p=1}^P\left(A{E}_{pt}^{\pm }+D{E}_{pth}^{\pm }\right)}+I{E}_{th}^{\pm }\ge AD{E}_{th}^{\pm },\hspace{0.17em}\forall t,h$$

(44)

$$(A{E}_{pt}^{\pm }+D{E}_{pth}^{\pm })\le S{T}_{pt}^{\pm }\cdot I{C}_{pt}^{\pm },\hspace{0.17em}\forall p,t,h$$

(45)

$$A{E}_{pt}^{\pm }\ge D{E}_{pth}^{\pm }\ge 0,\forall p,t,h$$

(46)

$$I{E}_{th}^{\pm }\le 40\text{\%}AD{E}_{th}^{\pm },\hspace{0.17em}\forall t,h$$

(47)

[Constraints for the district heating supply and demand balance]

$${\displaystyle \sum _{p=1}^P(A{H}_{pt}^{\pm }+D{H}_{pt}^{\pm })}\ge T{H}_{th}^{\pm },\hspace{0.17em}\forall t,h$$

(48)

$$A{H}_{pt}^{\pm }\ge D{H}_{pth}^{\pm }\ge 0,\forall p,t,h$$

(49)

[Constraint for combined heat and power generation balance]

$$Q^{\pm }{m}_{1t}^{\pm }=B_Q(\frac{A{E}_{1t}^{\pm }}{1-E{S}^{\pm }}C{E}^{\pm }+\frac{A{H}_{1t}^{\pm }}{1-H{S}^{\pm }}C{H}^{\pm })\hspace{0.17em}\forall t;$$

(50)

[Constraints for the heat-to-electric ratio of cogeneration plant]

$$A{H}_{1t}^{\pm }+D{H}_{1t}^{\pm }=(A{E}_{1t}^{\pm }+D{E}_{1th}^{\pm })\cdot {\kappa }^{\pm },\forall t$$

(51)

$$XE{H}_{1th}^{\pm }\cdot S{T}_{1t}^{\pm }=D{H}_{1t}^{\pm }\cdot {\kappa }^{\pm },\forall t$$

(52)

[Constraint for the total thermal efficiency of thermal power plant from national policy]

$$(A{E}_{1t}^{\pm }+A{H}_{1t}^{\pm })\ge 45\%\cdot Q^{\pm }{m}_t^{\pm },\forall t$$

(53)

[Constraints for environment capacity (CO2, PM, SO2, and NOx emission)]

$${\displaystyle \sum _{p=1}^P\left(A{E}_{pt}^{\pm }+D{E}_{pth}^{\pm }\right)\cdot {\xi }_{ipt}^{\pm }\cdot \left(1-{\eta }_{ipt}^{\pm }\right)\le MAG{E}_{it}},\hspace{0.17em}\forall i,t,\hspace{0.17em}h$$

(54)

[Constraints for installed capacity]

$$I{C}_{pt}^{\pm }=IC{P}_p+YE{H}_{pth}^{\pm }\cdot XE{H}_{pth}^{\pm }-CI{C}_{pt}^{\pm },t=1,\forall p,h$$

(55)

$$I{C}_{pt}^{\pm }=I{C}_{p(t-1)}^{\pm }+YE{H}_{pth}^{\pm }\cdot XE{H}_{pth}^{\pm }-CI{C}_{pt}^{\pm },t>1,\forall p,h$$

(56)

[Constraints for capacity expansion]

$$YE{H}_{pth}^{\pm }\{\begin{array}{l}=1,\hspace{0.17em}\mathrm{if} \mathrm{capacity} \mathrm{expansion} \mathrm{is} \mathrm{undertaken}\\ =0,\hspace{0.17em}\mathrm{otherwise}\end{array}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\forall \hspace{0.17em}p,\hspace{0.17em}t,\hspace{0.17em}h$$

(57)

$$0\le XE{H}_{pth}^{\pm }\le M_{pt}\cdot YE{H}_{pth}^{\pm },\forall p,t,h$$

(58)

[Constraints for generation proportion of different technologies]

$$A{E}_{1t}^{\pm }+D{E}_{1th}^{\pm }\le {\gamma }_t^{\pm }\cdot AD{E}_{dth}^{\pm },\hspace{0.17em}\forall \hspace{0.17em}t,h$$

(59)

$${\displaystyle \sum _{p=3}^5\left(A{E}_{pt}^{\pm }+D{E}_{pth}^{\pm }\right)}\ge \delta \cdot AD{E}_{dth}^{\pm }\text{}, \forall t, h$$

(60)

[Constraints for availabilities of energy resources]

$$\left(A{E}_{pt}^{\pm }+D{E}_{pth}^{\pm }\right)·E{E}_{pt}^{\pm }\le Z_{pt}^{\pm }, \forall p,t,h$$

(61)

[Robust constraints]

$${\epsilon }_{pth}^{\pm }-{\displaystyle \sum _{p=1}^P{\displaystyle \sum _{h=1}^H{p}_{th}\cdot {\epsilon }_{pth}^{\pm }+{\theta }_{pth}^{\pm }}}\ge 0,\hspace{0.17em}\forall p,t,h$$

(62)

where $AD{E}_{th}^{\pm }$ denotes the electricity demand under different electric-power deficiency levels $h$ during period t (GWh); $T{H}_{th}^{\pm }$ is the district heat demand under different deficiency levels $h$ during period $t$ (PJ); $Q^{\pm }$ is the heating value of coal (PJ/ton); $m_{1t}^{\pm }$ represents the coal quantity fed to combined heat and power (CHP) (ton); $B_Q$ denotes the calorific value of coal (PJ/ton); $E{S}^{\pm }$ is the electricity consumption rate of thermal power plant; $C{E}^{\pm }$ is the standard coal consumption of power generation of thermal power plant (ton/PJ); $H{S}^{\pm }$ represents the heat loss of the facilities; $C{H}^{\pm }$ is the standard coal consumption of heat supply of thermal power plant (ton/PJ); ${\kappa }^{\pm }$ denotes heat-to-electric ratio; ${\gamma }_t^{\pm }$ denotes the proportion of thermal power; ${\eta }_{ipt}^{\pm }$ is the removal efficiency of pollutant $i$ from technology $p$; ${\xi }_{ipt}^{\pm }$ denotes the emission intensity of pollutant $i$ from technology $p$ (10³ ton/GWh); $MAG{E}_{it}$ is the total allowable amount of pollutant i emission (10³ ton); $M_{pt}^{\pm }$ and $N_{pt}^{\pm }$ are the constraints for the upper and lower capacity expansion bound of technology $p$ (GW); $S{T}_{pt}^{\pm }$ is the operation hours of technology P in period t (h); $\delta$ denotes the percentage of power generation amount by renewable energy resources; $IC{P}_p$ is the initial installed capacity of power conversion technology $p$ (GW); $I{C}_{pt}^{\pm }$ represents the total installed capacity of technology p (GW); $CI{C}_{pt}^{\pm }$ denotes the closed installed capacity of “developing large units and suppressing small ones” in period t (GW); $E{F}_{pt}^{\pm }$ is the resources conversion efficiency of technology $p$ (PJ/GWh).

The planning horizon is considered as being from 2016 to 2021, and divided into two periods with a 3-year interval for each period. The related technical-economic information was obtained through analyzing many representative energy-related governmental reports and plans.

Table 1. Regional electricity and heat demand during the planning period.

Energy Demand	Demand Level	Probability (%)	T = 1	T = 2
Electricity demand (103 GWh)	Low	20	[97.11, 98.53]	[97.73, 99.11]
	Medium	60	[98.53, 100.14]	[99.21, 100.60]
	High	20	[99.90, 100.60]	[101.00, 102.60]
District heat quantity (PJ)	Low	20	[253.59, 259.59]	[255.00, 263.00]
	Medium	60	[278.68, 288.68]	[285.00, 293.00]
	High	20	[288.87, 297.87]	[295.00, 302.00]

presents power demands and the occurrence probabilities of each demand level (25%, 55%, and 20%). According to Zibo Statistics Bureau (from 1990 to 2014), and the forecasting information of electric-power demand by the government, three electricity generation targets are selected. Table 1 also shows the district heating demands during the planning horizon. To achieve the targets of renewable power generation and emission reduction, in the electric-power system, some scenarios are designed, which corresponds to environmental constraints and renewable power development constraints (i.e., renewable energy generation in period 1 and 2 accounts for 5% and 10% of the total regional power consumption, respectively).

4. Result Analysis and Discussion

4.1. Electricity-Generation Plan

Table 2. The optimized power generation schemes under λ = 0.

Technology	Level	Probability (%)	Optimized Generation Target (GWh)	Optimized Shortage Quantity (GWh)	Optimized Generation Quantity (GWh)
CHP	L	25	56,989.23	1240.12	58,229.35
	M	55	56,989.23	1240.12	58,229.35
	H	20	56,989.23	1240.12	58,229.35
	L-L	6.25	59,978.03	[8302.54, 8737.05]	[68,280.57, 68,715.08]
	L-M	13.75	59,978.03	[8302.54, 8737.05]	[68,280.57, 68,715.08]
	L-H	5	59,978.03	[8302.54, 8737.05]	[68,280.57, 68,715.08]
	M-L	13.75	59,978.03	[8302.54, 8737.05]	[68,280.57, 68,715.08]
	M-M	30.25	59,978.03	[8302.54, 8737.05]	[68,280.57, 68,715.08]
	M-H	11	59,978.03	[8302.54, 8737.05]	[68,280.57, 68,715.08]
	H-L	5	59,978.03	[8302.54, 8737.05]	[68,280.57, 68,715.08]
	H-M	11	59,978.03	[8302.54, 8737.05]	[68,280.57, 68,715.08]
	H-H	4	59,978.03	[8302.54, 8737.05]	[68,280.57, 68,715.08]
Hydropower	L	25	29.01	29.01	58.02
	M	55	29.01	29.01	58.02
	H	20	29.01	29.01	58.02
	L-L	6.25	30.85	30.85	61.70
	L-M	13.75	30.85	30.85	61.70
	L-H	5	30.85	30.85	61.70
	M-L	13.75	30.85	30.85	61.70
	M-M	30.25	30.85	30.85	61.70
	M-H	11	30.85	30.85	61.70
	H-L	5	30.85	30.85	61.70
	H-M	11	30.85	30.85	61.70
	H-H	4	30.85	30.85	61.70
Solar power	L	25	303.74	303.74	607.48
	M	55	303.74	303.74	607.48
	H	20	303.74	303.74	607.48
	L-L	6.25	406.53	[279.34, 315.44]	[685.87, 721.97]
	L-M	13.75	406.53	[279.34, 315.44]	[685.87, 721.97]
	L-H	5	406.53	[279.34, 315.44]	[685.87, 721.97]
	M-L	13.75	406.53	[279.34, 315.44]	[685.87, 721.97]
	M-M	30.25	406.53	[279.34, 315.44]	[685.87, 721.97]
	M-H	11	406.53	[279.34, 315.44]	[685.87, 721.97]
	H-L	5	406.53	[279.34, 315.44]	[685.87, 721.97]
	H-M	11	406.53	[279.34, 315.44]	[685.87, 721.97]
	H-H	4	406.53	[279.34, 315.44]	[685.87, 721.97]
Wind power	L	25	1691.76	1691.76	3383.52
	M	55	1691.76	1691.76	3383.52
	H	20	1691.76	1691.76	3383.52
	L-L	6.25	2077.83	2077.83	4155.66
	L-M	13.75	2077.83	2077.83	4155.66
	L-H	5	2077.83	2077.83	4155.66
	M-L	13.75	2077.83	2077.83	4155.66
	M-M	30.25	2077.83	2077.83	4155.66
	M-H	11	2077.83	2077.83	4155.66
	H-L	5	2077.83	2077.83	4155.66
	H-M	11	2077.83	2077.83	4155.66
	H-H	4	2077.83	2077.83	4155.66
Biomass and garbage power	L	25	877.53	877.53	1755.06
	M	55	877.53	877.53	1755.06
	H	20	877.53	877.53	1755.06
	L-L	6.25	2522.63	2522.63	5045.26
	L-M	13.75	2522.63	2522.63	5045.26
	L-H	5	2522.63	2522.63	5045.26
	M-L	13.75	2522.63	2522.63	5045.26
	M-M	30.25	2522.63	2522.63	5045.26
	M-H	11	2522.63	2522.63	5045.26
	H-L	5	2522.63	2522.63	5045.26
	H-M	11	2522.63	2522.63	5045.26
	H-H	4	2522.63	2522.63	5045.26

,

Table 3. The optimized power generation schemes under λ = 1.

Technology	Level	Probability (%)	Optimized Generation Target (GWh)	Optimized Shortage Quantity (GWh)	Optimized Generation Quantity (GWh)
CHP	L	25	56,989.23	0	56,989.23
	M	55	56,989.23	0	56,989.23
	H	20	56,989.23	0	56,989.23
	L-L	6.25	59,978.03	1536.75	61,514.78
	L-M	13.75	59,978.03	1536.75	61,514.78
	L-H	5	59,978.03	1536.75	61,514.78
	M-L	13.75	59,978.03	1536.75	61,514.78
	M-M	30.25	59,978.03	1536.75	61,514.78
	M-H	11	59,978.03	1536.75	61,514.78
	H-L	5	59,978.03	1536.75	61,514.78
	H-M	11	59,978.03	1536.75	61,514.78
	H-H	4	59,978.03	1536.75	61,514.78
Hydropower	L	25	29.01	29.01	58.02
	M	55	29.01	29.01	58.02
	H	20	29.01	29.01	58.02
	L-L	6.25	30.85	30.85	61.70
	L-M	13.75	30.85	30.85	61.70
	L-H	5	30.85	30.85	61.70
	M-L	13.75	30.85	30.85	61.70
	M-M	30.25	30.85	30.85	61.70
	M-H	11	30.85	30.85	61.70
	H-L	5	30.85	30.85	61.70
	H-M	11	30.85	30.85	61.70
	H-H	4	30.85	30.85	61.70
Solar power	L	25	303.74	303.74	607.48
	M	55	303.74	303.74	607.48
	H	20	303.74	303.74	607.48
	L-L	6.25	406.53	[255.63, 290.48]	[662.16, 697.01]
	L-M	13.75	406.53	[255.63, 290.48]	[662.16, 697.01]
	L-H	5	406.53	[255.63, 290.48]	[662.16, 697.01]
	M-L	13.75	406.53	[255.63, 290.48]	[662.16, 697.01]
	M-M	30.25	406.53	[255.63, 290.48]	[662.16, 697.01]
	M-H	11	406.53	[255.63, 290.48]	[662.16, 697.01]
	H-L	5	406.53	[255.63, 290.48]	[662.16, 697.01]
	H-M	11	406.53	[255.63, 290.48]	[662.16, 697.01]
	H-H	4	406.53	[255.63, 290.48]	[662.16, 697.01]
Wind power	L	25	1691.76	1691.76	3383.52
	M	55	1691.76	1691.76	3383.52
	H	20	1691.76	1691.76	3383.52
	L-L	6.25	2077.83	2077.83	4155.66
	L-M	13.75	2077.83	2077.83	4155.66
	L-H	5	2077.83	2077.83	4155.66
	M-L	13.75	2077.83	2077.83	4155.66
	M-M	30.25	2077.83	2077.83	4155.66
	M-H	11	2077.83	2077.83	4155.66
	H-L	5	2077.83	2077.83	4155.66
	H-M	11	2077.83	2077.83	4155.66
	H-H	4	2077.83	2077.83	4155.66
Biomass and garbage power	L	25	877.53	0	877.53
	M	55	877.53	[0, 80.6]	[877.53, 958.13]
	H	20	877.53	[68.45, 103.45]	[945.98, 980.98]
	L-L	6.25	2534.49	2212.79	4747.28
	L-M	13.75	2534.49	2212.79	4747.28
	L-H	5	2534.49	2212.79	4747.28
	M-L	13.75	2534.49	2357.88	4892.37
	M-M	30.25	2534.49	2357.88	4892.37
	M-H	11	2534.49	2357.88	4892.37
	H-L	5	2534.49	2534.49	5068.98
	H-M	11	2534.49	2534.49	5068.98
	H-H	4	2534.49	2534.49	5068.98

,

Table 4. The optimized power generation schemes under λ = 5.

Technology	Level	Probability (%)	Optimized Generation Target (GWh)	Optimized Shortage Quantity (GWh)	Optimized Generation Quantity (GWh)
CHP	L	25	56,989.23	226.47	57,215.70
	M	55	56,989.23	226.47	57,215.70
	H	20	56,989.23	226.47	57,215.70
	L-L	6.25	59,978.03	0	59,978.03
	L-M	13.75	59,978.03	0	59,978.03
	L-H	5	59,978.03	0	59,978.03
	M-L	13.75	59,978.03	0	59,978.03
	M-M	30.25	59,978.03	0	59,978.03
	M-H	11	59,978.03	0	59,978.03
	H-L	5	59,978.03	0	59,978.03
	H-M	11	59,978.03	0	59,978.03
	H-H	4	59,978.03	0	59,978.03
Hydropower	L	25	29.01	29.01	58.02
	M	55	29.01	29.01	58.02
	H	20	29.01	29.01	58.02
	L-L	6.25	30.85	30.85	61.70
	L-M	13.75	30.85	30.85	61.70
	L-H	5	30.85	30.85	61.70
	M-L	13.75	30.85	30.85	61.70
	M-M	30.25	30.85	30.85	61.70
	M-H	11	30.85	30.85	61.70
	H-L	5	30.85	30.85	61.70
	H-M	11	30.85	30.85	61.70
	H-H	4	30.85	30.85	61.70
Solar power	L	25	303.74	303.74	607.48
	M	55	303.74	303.74	607.48
	H	20	303.74	303.74	607.48
	L-L	6.25	406.53	[255.63, 290.48]	[662.16, 697.01]
	L-M	13.75	406.53	[255.63, 290.48]	[662.16, 697.01]
	L-H	5	406.53	[255.63, 290.48]	[662.16, 697.01]
	M-L	13.75	406.53	[255.63, 290.48]	[662.16, 697.01]
	M-M	30.25	406.53	[255.63, 290.48]	[662.16, 697.01]
	M-H	11	406.53	[255.63, 290.48]	[662.16, 697.01]
	H-L	5	406.53	[255.63, 290.48]	[662.16, 697.01]
	H-M	11	406.53	[255.63, 290.48]	[662.16, 697.01]
	H-H	4	406.53	[255.63, 290.48]	[662.16, 697.01]
Wind power	L	25	1691.76	1691.76	3383.52
	M	55	1691.76	1691.76	3383.52
	H	20	1691.76	1691.76	3383.52
	L-L	6.25	2077.83	2077.83	4155.66
	L-M	13.75	2077.83	2077.83	4155.66
	L-H	5	2077.83	2077.83	4155.66
	M-L	13.75	2077.83	2077.83	4155.66
	M-M	30.25	2077.83	2077.83	4155.66
	M-H	11	2077.83	2077.83	4155.66
	H-L	5	2077.83	2077.83	4155.66
	H-M	11	2077.83	2077.83	4155.66
	H-H	4	2077.83	2077.83	4155.66
Biomass and garbage power	L	25	945.98	0	945.98
	M	55	945.98	[0, 12.15]	[945.98, 958.13]
	H	20	945.98	[0, 35]	[945.98, 980.98]
	L-L	6.25	2787.65	1959.62	4747.27
	L-M	13.75	2787.65	1959.62	4747.27
	L-H	5	2787.65	1959.62	4747.27
	M-L	13.75	2787.65	2104.72	4892.37
	M-M	30.25	2787.65	2104.72	4892.37
	M-H	11	2787.65	2104.72	4892.37
	H-L	5	2787.65	2281.33	5068.98
	H-M	11	2787.65	2281.33	5068.98
	H-H	4	2787.65	2281.33	5068.98

and

Table 5. The optimized power generation schemes under λ = 50.

Technology	Level	Probability (%)	Optimized Generation Target (GWh)	Optimized Shortage Quantity (GWh)	Optimized Generation Quantity (GWh)
CHP	L	25	56,989.23	0	56,989.23
	M	55	56,989.23	0	56,989.23
	H	20	56,989.23	0	56,989.23
	L-L	6.25	59,978.03	1244.54	61,222.57
	L-M	13.75	59,978.03	1244.54	61,222.57
	L-H	5	59,978.03	1244.54	61,222.57
	M-L	13.75	59,978.03	1244.54	61,222.57
	M-M	30.25	59,978.03	1244.54	61,222.57
	M-H	11	59,978.03	1244.54	61,222.57
	H-L	5	59,978.03	1244.54	61,222.57
	H-M	11	59,978.03	1244.54	61,222.57
	H-H	4	59,978.03	1244.54	61,222.57
Hydropower	L	25	29.01	29.01	58.02
	M	55	29.01	29.01	58.02
	H	20	29.01	29.01	58.02
	L-L	6.25	30.85	30.85	61.70
	L-M	13.75	30.85	30.85	61.70
	L-H	5	30.85	30.85	61.70
	M-L	13.75	30.85	30.85	61.70
	M-M	30.25	30.85	30.85	61.70
	M-H	11	30.85	30.85	61.70
	H-L	5	30.85	30.85	61.70
	H-M	11	30.85	30.85	61.70
	H-H	4	30.85	30.85	61.70
Solar power	L	25	303.74	303.74	607.48
	M	55	303.74	303.74	607.48
	H	20	303.74	303.74	607.48
	L-L	6.25	406.53	[255.63, 272.34]	[662.16, 678.87]
	L-M	13.75	406.53	[255.63, 272.34]	[662.16, 678.87]
	L-H	5	406.53	[255.63, 272.34]	[662.16, 678.87]
	M-L	13.75	406.53	[255.63, 272.12]	[662.16, 678.65]
	M-M	30.25	406.53	[255.63, 272.12]	[662.16, 678.65]
	M-H	11	406.53	[255.63, 272.12]	[662.16, 678.65]
	H-L	5	406.53	[255.63, 290.11]	[662.16, 696.64]
	H-M	11	406.53	[255.63, 290.11]	[662.16, 696.64]
	H-H	4	406.53	[255.63, 290.11]	[662.16, 696.64]
Wind power	L	25	1691.76	1691.76	3383.52
	M	55	1691.76	1691.76	3383.52
	H	20	1691.76	1691.76	3383.52
	L-L	6.25	2077.83	2077.83	4155.66
	L-M	13.75	2077.83	2077.83	4155.66
	L-H	5	2077.83	2077.83	4155.66
	M-L	13.75	2077.83	2077.83	4155.66
	M-M	30.25	2077.83	2077.83	4155.66
	M-H	11	2077.83	2077.83	4155.66
	H-L	5	2077.83	2077.83	4155.66
	H-M	11	2077.83	2077.83	4155.66
	H-H	4	2077.83	2077.83	4155.66
Biomass and garbage power	L	25	945.98	0	945.98
	M	55	945.98	[0, 12.15]	[945.978, 958.13]
	H	20	945.98	[0, 35]	[945.978, 980.98]
	L-L	6.25	2787.65	1959.62	4747.27
	L-M	13.75	2787.65	1959.62	4747.27
	L-H	5	2787.65	1959.62	4747.27
	M-L	13.75	2787.65	2104.72	4892.37
	M-M	30.25	2787.65	2104.72	4892.37
	M-H	11	2787.65	2104.72	4892.37
	H-L	5	2787.65	2281.33	5068.98
	H-M	11	2787.65	2281.33	5068.98
	H-H	4	2787.65	2281.33	5068.98

present the optimal solutions of electric-power generation schemes of different technologies with different λ values, under different demand levels, during the whole planning horizon. The optimal combined electricity and heat generation targets in period 2 would be greater than that in period 1. In period 1, the generation amount of combined electricity and heat would be 56.99 × 10³ GWh in period 1, and 59.98 × 10³ GWh in period 2 under different λ values. Furthermore, power generation amount of CHP would increase. For example, in period 1, under medium demand level, power generation amount by CHP would be 58.23 × 10³ GWh, 56.99 × 10³ GWh, 57.22 × 10³ GWh, and 56.99 × 10³ GWh, as λ is fixed with the values of 0, 1, 5 and 50, respectively; under medium–medium level (with the probability of 30.25%) in period 2, power generation amount would be (68.28, 68.72) × 10³ GWh, 61.51 × 10³ GWh, 59.98 × 10³ GWh, and 61.22 × 10³ GWh, respectively. It indicated that the CHP is a more economical and stable way for power supply with the demand of electricity increasing, and along with regional electric-power structure optimization, the combined heat and power would still be the main choice for supporting regional electric-power supply.

Among these renewable power generation technologies, clean electricity would mainly come from solar power, wind power, and biomass and garbage power (BGP). The optimized electricity generation for wind power would be 3.38 × 10³ GWh and 4.16 × 10³ GWh in periods 1 and 2, respectively. The wind power would play a significant role in renewable power development during the planning horizon. For example, in period 1, wind power generation would occupy about 3% of total electricity consumption, and 60% of total renewable power generation under different demand level; in period 2, the proportion would increase from about 3% to 4% of total electricity consumption, and be 40% of total renewable power generation. Since wind power possesses the characteristic of cleanliness and the condition of convenience in this region, wind power would be developed as a priority. In addition, BGP power generation would increase significantly during the whole planning horizon. For example, in period 1 under medium level (with the probability of 55%), power generation amount of BGP would be 1.76 × 10³ GWh, (877.53, 958.13) GWh, (945.98, 958.13) GWh, and (945.98, 958.13) GWh under λ with the values of 0, 1, 5, and 50, respectively; in period 2 under medium–medium level, power generation amount of BGP would be 5.05 × 10³ GWh, 5.07 × 10³ GWh, 5.07 × 10³ GWh, and 5.07 × 10³ GWh with λ fixed as 0, 1, 5 and 50, respectively. The proportion of biomass and garbage power generation would rise from about 1% in period 1, to 5% in period 2 of total electricity consumption, and 20% in period 1 to 50% in period 2 of total renewable power generation. In Zibo city, the hydropower would have a smaller scale under water resource and geography limitation. In general, renewable power generation amount would change as λ values vary, and the stability of the regional electric-power supply would be enhanced as the total renewable power generation amount increases.

Figure 2, Figure 3, Figure 4 and Figure 5 show the optimized solutions for electric-power generation schemes under different λ values. Electric-power generation amount of CHP would be decreased as λ increases. For example, in period 2 under medium–medium level, electric-power amount generated by CHP would be (68.28, 68.72) × 10³ GWh, 61.51 × 10³ GWh, 59.98 × 10³ GWh, and 61.22 × 10³ GWh under λ fixed as 0, 1, 5, and 50, respectively. It indicated that the risk of system failure, which means higher CO₂ and pollutants discharged from cogeneration exceeding the regulated limitation, would decrease as λ increases. In general, relatively lower power generation of CHP would promote emissions reduction and evade the risk of regional energy system.

As shown in Figure 3, solar power generation amount would be decreased as λ value increases. For example, in period 2 under medium–medium level, the electricity generated by solar power would be (685.87, 721.97) GWh, (662.16, 697.01) GWh, (662.16, 697.01) GWh, and (662.16, 678.65) GWh as λ is fixed with the values of 0, 1, 5, and 50, respectively. Since the regional power supply of solar power has the characteristic of instability and higher cost, the stability and security of system power supply would increase as λ increases. Electric power generated by BGP would decrease as λ increases (Figure 5). For instance, under λ fixed with the values of 0, 1, 5, and 50, power generation amount of BGP would be 1.76 × 10³ GWh, (877.53, 958.13) GWh, (945.98, 958.13) GWh, and (945.98, 958.13) GWh in period 1 under medium level, respectively; the generation amount would be 5.05 × 10³ GWh, 5.07 × 10³ GWh, 5.07 × 10³ GWh, and 5.07 × 10³ GWh under medium–medium level in period 2, respectively. A higher power generation of BGP would lead to a higher pollutants and CO₂ emission, which would violate environmental constraints of the system. As λ increases, the power generation of BGP would be reduced. In summary, the total renewable power generation amount would increase as λ values increase. Thus, as λ values increases, the system failure risk would be lessened; meanwhile, the security and stability would be enhanced.

4.2. Imported Electricity Scheme

Figure 6 presents the imported electric power amount during the planning horizon. It would decrease from period 1 to 2 under different power demand levels. For example, with λ fixed as 0,1,5, and 50, the imported power amount would be (34.5, 35.41) × 10³ GWh, (36.62, 38.15) × 10³ GWh, (36.32, 37.92) × 10³ GWh, and (36.55, 38.15) × 10³ GW h under medium level in period 1, respectively; (20.98, 21.9) × 10³ GWh, (28.1, 29.46) × 10³ GWh, (29.46, 30.82) × 10³ GWh, and (28.21, 29.59) × 10³ GWh under medium–medium demand level in period 2, respectively. It indicated that the imported power amount would be decreased with regional power generation and power structure adjustment increasing, and the amount of imported electricity would rise as λ increases. For instance, under medium level in period 1, the amount would be (34.5, 35.41) × 10³ GWh under λ = 0 and (36.55, 38.15) × 10³ GWh under λ = 50. As a result, it would lead to a smaller system risk and enhanced system feasibility, which could also promote the energy conservation and emissions reduction to some degree.

4.3. CO₂ and Air Pollution Control

Table 6. The amount of CO₂ and air pollution emissions under different λ values.

Gaseous Emission	λ Level	T = 1	T = 2
CO2 (106 ton)	λ= 0	[61.88, 63.11]	[63.07, 64.10]
	λ = 1	[60.3, 60.77]	[56.82, 57.40]
	Λ = 5	[60.54, 61.01]	[55.60, 56.17]
	λ = 50	[60.31, 60.77]	[56.72, 57.30]
SO2 (103 ton)	λ = 0	[70.35, 82.75]	[51.58, 76.43]
	λ = 1	[68.76, 79.91]	[46.37, 68.28]
	λ = 5	[69.03, 80.23]	[45.35, 66.78]
	λ = 50	[68.76, 79.91]	[46.28, 68.16]
NOx (103 ton)	λ = 0	[60.36, 86.29]	[40.51, 81.00]
	λ = 1	[59.04, 83.38]	[36.41, 72.34]
	λ = 5	[59.27, 83.71]	[35.61, 70.75]
	λ = 50	[59.03, 83.38]	[36.35, 72.21]
PM (103 ton)	λ = 0	[9.74, 12.82]	[5.57, 8.04]
	λ = 1	[9.53, 12.38]	[5.00, 7.18]
	λ = 5	[9.56, 12.43]	[4.89, 7.02]
	λ = 50	[9.52, 12.38]	[5.00, 7.17]

shows the solutions of optimized air pollutants and CO₂ emission. The air pollutants and CO₂ emission amount would decrease. For example, under λ = 5, the amount of CO₂ emissions would decrease from (60.54, 61.01) × 10⁶ ton in period 1 to (55.6, 56.17) × 10⁶ ton in period 2; the amount of SO₂ emissions would decrease from (69.03, 80.23)× 10³ ton in period 1 to (45.35, 66.78) × 10³ ton in period 2; the amount of NO_x emissions would be (59.27, 83.71) × 10³ ton and (35.61, 70.75) × 10³ ton in period 1 and 2; the amount of PM₁₀ emissions would decrease from (9.56, 12.43) × 10³ ton in period 1 to (4.89, 7.02) × 10³ ton in period 2, respectively. The reasons for decreasing emission are that firstly, the technology and facilities would be updated to reduce the average emissions level; secondly, due to the power structure optimization in the first period, renewable energy has been developed to some degree, which could make contributions to energy conservation and emissions reduction. In addition, the effect of emission mitigation could be better under considering system risk aversion. It indicates that the results would lead to a lower system risk and more robust regional energy system, which is important for achieving sustainable development and better environment quality.

4.4. System Cost

Figure 7 shows the total system costs under different scenarios during the planning periods. The energy system cost in Zibo city would have a slight increase trend as λ levels increase. For instance, under the scenarios of λ with the values of 0, 1, 5, and 50, the system cost would be RMB¥ (490.63, 651.16) × 10⁹, RMB¥ (499.65, 659.71) × 10⁹, RMB¥ (502.62, 662.94) × 10⁹, and RMB¥ (502.82, 663.56) × 10⁹, respectively. As λ levels increasing, the system failure risk would be reduced, and the system cost would be increased. Conversely, a lower λ level would bring about a higher system risk and a lower system cost. It indicated that if the decision-makers aim to lower costs, a higher system risk may occur.

Based on the above analyses, these indicated that the optimized solutions are able to support regional energy system management for making integrated schemes of power generation, capacity expansion, air pollutant and CO₂ emission reduction under different renewable energy development targets and environmental quality requirements. The solutions with lower and upper bounds are helpful for generating decision alternatives representing various options. Cost–risk analysis can be obtained through integrated the stochastic robust optimization method into the multistage stochastic programming in regional energy system management.

5. Conclusions

In this study, a multistage stochastic inexact robust programming was proposed for supporting regional electric-power system structure optimization and management. The model covered the district heating supply, power generation, and air pollutant mitigation coupled with relevant technique constraints and governmental policies. Comparing the solutions optimized by the model with the strategies carried out in the real world, the former emerged with obvious advantages, which can lead to a more prosperous future. In addition, the developed method could be valuable for obtaining trade-off schemes between system economy and risks that are introduced by system’s uncertainties according to decision-makers’ willingness. An energy system structure management of Zibo City, China, is used as a case study for verifying the efficiency of the developed model. Optimized schemes of power generation, capacity expansion, air pollutant and CO₂ emission reduction, and system cost were analyzed. The results indicated that under different requirements of renewable energy development, and pollutant and CO₂ mitigation, traditional power generation technology would still be increased, attributing to its lower costs and traditional energy resources structure based on the thermal power generation. In addition, renewable energy would also play an important role in solving energy, resource, and environmental pressures; renewable power generation amount would be rising continuously, though it might develop slowly for a certain period of time.

However, a number of limitations also exist in the proposed model of this study. Firstly, in the optimization model, many energy industrial processes are not considered, and only generation processes and energy-related environmental problems are involved in this study. In order to obtain more comprehensive management schemes, more energy development and utilization patterns could be considered. Second, compared with other optimization methods, the model would be infeasible in addressing the high uncertainties in the model parameters; and through introducing different λ values in the model, regional energy-managers cannot directly obtain suitable management schemes. Therefore, further research can strengthen knowledge and mitigate these limitations in the future.