A Comprehensive Forecasting–Optimization Analysis Framework for Environmental-Oriented Power System Management—A Case Study of Harbin City, China

Abstract

In this study, a comprehensive research framework coupled with electric power demand forecasting, a regional electric system planning model, and post-optimization analysis is proposed for electric power system management. For dealing with multiple forms of uncertainties and dynamics concerning energy utilization, capacity expansions, and environmental protection, the inexact two-stage stochastic robust programming optimization model was developed. The novel programming method, which integrates interval parameter programming (IPP), stochastic robust optimization (SRO), and two-stage stochastic programming (TSP), was applied to electric power system planning and management in Harbin, China. Furthermore, the Gray-Markov approach was employed for effective electricity consumption prediction, and the forecasted results can be described as interval values with corresponding occurrence probability, aiming to produce viable input parameters of the optimization model. Ten scenarios were analyzed with different emissions reduction levels and electricity power structure adjustment modes, and the technique for order of preference by similarity to ideal solution (TOPSIS) was selected to identify the most influential factors of planning decisions by selecting the optimal scheme. The results indicate that a diversified power structure that dominates by thermal power and is mainly supplemented by biomass power should be formed to ensure regional sustainable development and electricity power supply security in Harbin. In addition, power structure adjustment is more effective than the pollutants emission control for electricity power system management. The results are insightful for supporting supply-side energy reform, generating an electricity generation scheme, adjusting energy structures, and formulating energy consumption of local policies.

1. Introduction

Top speed economic development, sustained industrial expansion, and soaring energy demands all place a higher pressure on decision-makers to establish energy-related plans that address, for example, energy utilization, structure adjustment, and pollutant emission reduction, especially for China. Confronted with a rigorous energy situation, a coal-dominating energy consumption structure constitutes a serious impediment to sustainable development [1,2]. In addition, resource depletion, total quantity control of atmospheric pollutants, and defective energy system management exacerbate this imbalance to various degrees. However, in the process of establishing a regional electric power management scheme, decision-makers are confronted with multiple forms of uncertain information relating to generation and consumption parameters, which are introduced from the features and purposes of the energy system itself, as well as to social–economic and technology factors. Such uncertainties directly weaken the viability of generating expected electric power system management strategies [3,4,5].

Therefore, it is desirable to develop a comprehensive and effective research framework with an optimization model as a core for electric power management under unexpected multiple uncertainties. Previously, focused on electric power planning and policy-making, a number of optimization models were proposed worldwide for tackling the aforementioned uncertainties and complexities at different temporal and spatial scales [6,7,8,9]. For example, Santos and Legey proposed an optimization formulation for long-term electricity system expansion planning, taking environmental and operation costs into consideration [10]. Sheikhahmadi et al. developed a risk-based two-stage stochastic optimization model for microgrid system operation management to address the uncertainties of renewable energy [11]. Liu et al. developed a hybrid optimization model based on the chance-constrained method for the planning of coupled coal and power management systems under China’s special coal pricing mechanism [12]. Wu et al. proposed a general optimization framework by integrating the inexact fuzzy programming method and the inexact stochastic programming method to the heat supply management of an actual wind power heating system, where uncertainties were presented in multiple forms [13].

In addition, many applications have been enforced at different scales through interval parameter, stochastic, fuzzy mathematical programming, and diverse coupled optimization methods. Among these methods, the two-stage stochastic programming (TSP) and interval parameter programming (IPP) coupled methods can tackle random interval information, and have been extensively applied to step-forward planning and management in energy systems. Huang et al. established a multi-region two-stage stochastic optimization model to address demand uncertainty, and the proposed model was applied to Taiwan’s electricity sector [14]. Ahn and Han proposed a two-stage stochastic model to form a mixed network, and optimized public utility supplement and carbon dioxide emission reduction strategies [15]. Ji et al. proposed a method that combined two-stage stochastic programming and interval parameter programming to reflect the uncertainties of management, technology, economy, and policy in the power system of Ningxia Hui Autonomous Region, China [16]. Zhang et al. proposed an interval two-stage stochastic programming model for grasping the direction of regional energy structure adjustment, controlling resources allocation patterns, and formulating local energy consumption policies in Heilongjiang Province, China [17].

In interval two-stage stochastic programming (ITSP), the first decision needs to be made before the random event occurs, and the second decision is then undertaken to minimize the “punishment” to avoid the occurrence of infeasibilities after random events. Although the ITSP method is valid in dealing with uncertainties in the form of probability in both sides of the model, it is infeasible in directly sheltering the risk of stochastic events, and this defect could threaten the stability of the whole system. Following this view, stochastic robust optimization (SRO) was introduced to overcome this shortcoming, relying on the application of a risk-aversion coefficient to optimization approaches and acquiring robust optimization schemes of energy management issues [18]. Xie et al. established a multi-stage stochastic robust model to study the direction of energy structure adjustment and the policy of pollutants emission prevention in Jining City, China [19]. Guo et al. established an optimization model that integrated ITSP and SRO for programming industrial systems and the ecosystem carbon market in Zhangjiakou, China [20]. Xie et al. presented a multi-level SRO model for outputting detailed electric power production plans and power structure adjustment reform strategies under uncertainty in Zibo City, China [21]. Jin et al. proposed an SRO model combined with interval fuzzy programming for addressing various uncertainties and for avoiding inherent risks in order to support the low-carbon transformation of power systems [22]. Considering that the SRO method is expert at balancing quantitative assessment between economy and stability, integrated interval two-stage stochastic programming and stochastic robust optimization (ITSP–SRO) within a general model is a comprehensive approach for energy system management.

Electric power consumption forecasting is an essential process in an energy system optimization model, and the accuracy of forecasting is driven by a number of complicated factors. A mathematical model can be an efficient means to predict and analyze regional electric power consumption, and previous research has made valuable attempts in the field of developing forecasting models, including artificial neural networks, adaptive networks, regression analysis, and multiple linear regression analysis. Ivanin and Direktor applied artificial neural networks to forecast the short-term electricity power consumption of independent users located outside a centralized power supply system [23]. Chahkoutahi and Khashei built several hybrid models through integrating adaptive networks and fuzzy inference for forecasting the electricity load of the energy market [24]. Cao and Wu developed a forecasting monthly electricity consumption model by coupled regression analysis and the fruit fly optimization algorithm [25]. Lasso et al. predicted the short-term electric power consumption in South Africa during the peak period from 2009 to 2012 by using the multiple linear regression analysis model [26]. Differing from these techniques, the Gray forecasting model incorporated with the Markov chain forecasting model can make accurate predictions even with limited historical data, and can deal with dynamic varying time series in the system [27]. In addition, the Gray-Markov model (GMM) is well suited for forecasting the development of the system on the basis of transition probabilities among different states that could reflect the influence of all stochastic and uncertainty factors, and it is useful for assisting decision-makers to precisely forecast electricity consumption/demand. Hence, introducing the GMM into an electric power planning system for determining the significant input parameters of the ITSP–SRO optimization model is a potential approach, although few researches have investigated this.

Furthermore, post-optimization analysis is an attractive technique that can capture and reflect the critical factors of electric power system management programs through the use of the multi-criteria decision-making (MCDM) approach [28]. A number of innovated MCDM methods have been extensively applied in practical decision problems, such as the simple additive weighted method [29], the weighted product method [30], and the analytic hierarchy process [31], as well as the interactions among them [32,33]. Among these methods, the technique for order preference by similarity to ideal solution (TOPSIS) is an available tool for solving real-world decision puzzles by employing understandable mathematical concepts to quantify the relative performance of decision schemes [34].

Although many studies have been performed based on approaches for coupled electricity load forecasting and energy system optimization with multiple uncertainties, a few challenges and issues concerning system application of prediction–optimization and post-optimization under multi-scenario analysis still exist, according to the aforementioned literature: (a) Few studies focused on the mathematical correlation between electricity prediction and optimization process. For the real-world energy system, electricity consumption in the future and the corresponding electricity generation both represent strong stochasticity. If the relationship of the two stochastic variables cannot be properly reflected and handled, the reliability of the optimized strategies would be reduced; (b) Most of the conventional stochastic optimization methods neglect the risk of random events in the optimization process, and this limitation poses threats to system stability; (c) Few studies applied reasonable methods to analyze post-optimization for multi-scenario schemes. The related energy policies and measures (e.g., energy structure adjustment and controlling of the amount of atmospheric pollutants) could create an incentive for sustainable development, and the effectiveness of these measures should be quantified and compared.

Aiming to resolve these problems, the purpose of this study is to develop a comprehensive forecasting-optimization framework, coupled with electric power consumption forecasting, the regional electric system optimization model, and the post-optimization analysis approach, as shown in Figure 1. By integrating the methods of IPP, TSP, and SRO, a novel ITSP–SRO model is developed for supporting regional electric power system administration, considering power structure adjustment and pollutant mitigation scenarios. Compared with typical hybrid inexact optimization tools, this method can not only reflect multiple uncertainties expressed as interval values and probability distribution, but can also make a trade-off between system risk and cost according to the decision-makers’ attitudes. Moreover, the Gray-Markov forecasting approach is selected for electricity consumption prediction, which can achieve the goal of dynamic forecasting, continuously modifying the consumption values. The forecasting results are presented in the form of interval parameters with corresponding random probability, and are highly consistent with the formal requirements for the input parameters of the ITSP–SRO model. Through scenario analysis, multiple decision-making plans can be obtained under different thermal power structure and atmospheric pollutant mitigation scenarios. Additionally, the TOPSIS method is employed to identify the most influential factors of planning decisions through selecting the optimal scheme. The TOPSIS method is an operational evaluation and decision support approach that is expert at addressing conflicting objectives in the decision-making process. Hence, a comprehensive forecasting–optimization framework is desired to address the multiple uncertainties and to support policy analysis in practical electricity system planning.

The major contributions of this study are as follows:

(a)

A comprehensive framework of forecasting, optimization, and post-optimization for electric power system management is developed for providing a new perspective to optimization-centred research;

(b)

The proposed ITSP–SRO model can effectively reflect uncertain variables, can evaluate the trade-off between system economy and stability, and can generate robust solutions for electric power system management;

(c)

The designed approach is applied to the electric power system in Harbin, China, for detailed planning and management, which can provide reasonable strategies for multi-scenarios from different conversion technologies under uncertainty;

(d)

By comparing the importance of energy structure adjustment and controlling the amount of atmospheric pollutants, the most influential factors of planning decisions are identified, and the optimum planning schemes are obtained.

In general, these solutions can provide meticulous and overall management schemes for decision-making departments from the perspective of comprehensive trade-offs among energy utilization, economic development, and environmental protection.

2. Methodology

2.1. Gray-Markov Forecasting Approach

Accurate forecasting for electricity consumption has profoundly influenced electricity planning and management. An electricity consumption system is a gray system that exhibits the features of uncertainty and dynamic; thus, the Gray-Markov forecasting approach is an effective method for electricity consumption prediction of the regional electric power system, which is based on the Gray forecasting model (GM) and the Markov chain forecasting model. The main feature of Gray forecasting is to obtain the approximate exponential law by processing the original data sequence, and the GM (1,1) model is a prominent one that can utilize differential equations to excavate the variation essence of the system [35]. Moreover, for the comprehensive prediction of an event, it is not only that all possible results should be exhibited, but also that the probability of each prediction result should be given. By introducing the transition probability matrix, the Markov chain can predict the fluctuating trend according to the current situation of events [36]. The integrated approach is summarized as follows (Equation (1)):

$$X^{(0)}=\left\{X^{(0)}(1),X^{(0)}(2),X^{(0)}(3),\dots ,X^{(0)}(n)\right\}.$$

(1)

Formulate an original series $X^{(0)}$ as Equation (1), where $X^{(0)}$ is the baseline data series and n denotes the sample size of $X^{(0)}$. Obtain a new sequence by a one-time accumulated generating operation (AGO) as Equation (2).

$$X^{(1)}=\left\{X^{(1)}(1),X^{(1)}(2),X^{(1)}(3),\dots ,X^{(1)}(n)\right\}.$$

(2)

where $X^{(1)}(k)={\displaystyle \sum _{i=1}^k{x}^{(0)}}(i),k=1,2,3,\dots ,n.$

$$\left[d{x}^{(1)}/dt\right]+az=b.$$

(3)

Establish a first-order differential equation that can be expressed as Equation (3), where $z^{(1)}(k)=\lambda x^{(1)}(k)+(1-\lambda )x^{(1)}(k+1)$, $k=1,2,3,\dots ,(n-1).$ $\lambda$ denotes a horizontal adjustment coefficient whose theoretical value range is (0,1), and the value of $\lambda$ is 0.5 in this study. $a$ and $b$ are undetermined parameters in differential equations, and they can be obtained by using the least squares method operation from Equations (4) and (5).

$${\widehat{x}}^{(1)}(k+1)=\left[x^{(0)}(1)-\frac{b}{a}\right]\cdot e^{-ak}+\frac{b}{a}.$$

(4)

$$\widehat{Y}(k)={\left\{x^{(0)}(2),x^{(0)}(3),x^{(0)}(4),\dots ,x^{(0)}(n)\right\}}^T.$$

(5)

where

$$\widehat{\theta }=\left[\begin{array}{c}a\\ b\end{array}\right]={(B^{T}B)}^{-1}{B}^{T}Y,B=\left[\begin{array}{cc}-z^{(1)}(2)& 1\\ -z^{(1)}(3)& 1\\ \dots & \dots \\ -z^{(1)}(n)& 1\end{array}\right]$$

As the forecasting sequence is built by AGO, inverse accumulated generation operational (IAGO) can be employed for reverse prediction in Equation (6):

$${\widehat{x}}_0^{(0)}(k)={\widehat{x}}_{}^{(1)}(k)-{\widehat{x}}_{}^{(1)}(k-1),k=2,3,4,\dots ,n.$$

(6)

Therefore, a trend prediction equation can be proposed by Equation (7):

$$\widehat{Y}(k)={\widehat{x}}_{}^{(0)}(k+1)=r{e}^{-ak}$$

(7)

where

$$r=\left[x^{(0)}(1)-\frac{b}{a}\right](1-e^{\lambda }),k=0,1,2,\dots ,n-1.$$

$\widehat{Y}(k)$ denotes the forecasting data. The relative error is utilized to examine the predication accuracy of the forecasting model as Equation (8):

$$\xi (k)=\frac{x_{}^{(0)}(k+1)-{\widehat{x}}_{}^{(0)}(k+1)}{x_{}^{(0)}(k+1)}\times 100\%,k=1,2,3,\dots n.$$

(8)

$X^{(0)}(k+1)$ is distributed in the range of $\widehat{Y}(k)$, which can be classified as several continuous intervals. $X^{(0)}(k+1)$ is treated as corresponding to state $E_i$ in a Markov chain when it drops into interval i of U. $E_i$ can be expressed as Equations (9)–(11):

$$E_i=\left[E_{1i},E_{2i}\right]$$

(9)

$$E_{1i}=\widehat{Y}(k)+A_i$$

(10)

$$E_{2i}=\widehat{Y}(k)+B_i$$

(11)

where i = 1, 2, 3…, U, U is the number of states, and $\widehat{Y}(k)$ is a time function. $A_i$ and $B_i$ denote the lower bound and upper bound compensate value of $\widehat{Y}(k)$ in state $E_i$, respectively.

$${\pi }_j=P_{ij}(m)=\left\{\frac{M_{ij}(m)}{M_i}\right\}i,j=1,2,3,\dots ,U,{\displaystyle \sum _{j=1}^J{\pi }_j=1}.$$

(12)

where $P_{ij}\left(m\right)$ denotes the transition probability from $E_j$ to $E_i$ through m steps; ${\displaystyle \prod =({\pi }_1},{\pi }_2,{\pi }_3,\dots ,{\pi }_j)$ represents the state transition probability vector of state $E_j$; $M_{ij}(m)$ is the number of state $E_j$; $M_i$ is the number of state $E_i$; and the values of $P_{ij}\left(m\right)$ denote a transition probability matrix W(m) as Equation (13).

$$W(m)=\left[\begin{array}{cccc}{p}_{11}(m)& p_{12}(m)& \dots & p_{1j}(m)\\ p_{21}(m)& p_{22}(m)& \dots & p_{2j}(m)\\ & & \dots & \\ p_{j1}(m)& p_{j2}(m)& \dots & p_{ij}(m)\end{array}\right]i,j=1,2,3,\dots ,U.$$

(13)

Furthermore, the one-step transition matrix W (1) should be observed. Suppose the target to be predicted is in state $E_K(1\le K\le U)$, row K in W (1) should be noticed. In the case of $\max P_{Kj}(1)=P_{KV}(1)$ $j=1,2,3,\dots ,U;1\le K\le U$, then the transition from state $E_K$ to $E_V$ is most likely to occur in the system in the next moment. Subsequently, calculate the forecasting data by Equation (14):

$$\widehat{Y}(k+1)=\frac{\widehat{Y}(k)}{1\pm \frac{1}{2}\left(A_i+B_i\right)}$$

(14)

When the Markov state reaches equilibrium, the valve of the transition probability vector does not change, no matter what further iteration happens, and the state transition probability vector is defined as the limiting transition probability vector of the Markov process (Equations (15) and (16)).

$${\displaystyle \prod =[{\pi }_1(k)},{\pi }_2(k),\dots ,{\pi }_j(k)]=[\underset{k\to \infty }{\lim }{\pi }_1(k),\underset{k\to \infty }{\lim }{\pi }_2(k),\dots ,\underset{k\to \infty }{\lim }{\pi }_j(k)]$$

(15)

$$j=1,2,\dots ,U;k=1,2,3,\dots n;{\displaystyle \sum _{j=1}^J{\pi }_j=1}.$$

(16)

2.2. Inexact Two-Stage Stochastic Robust Programming

The ITSP–SRO model was integrated with IPP programming, TSP programming, and SRO programming. These approaches make distinctive contributions to strengthen the functionality of the programming model, such as by reflecting uncertainties through discrete intervals, combining target and punishment in decision-making, and executing risk aversion of the system.

Generally, TSP is summarized as Equations (17)–(20) [37]:

$$z=\min C^{T}X+E_{\omega \in \mathsf{\Omega }}[Q(X,\omega )]$$

(17)

Subject to:

$$x\in X$$

(18)

$$Q(X,\omega )=\min f{(x)}^{T}y$$

(19)

Subject to:

$$D(\omega )y\ge h(\omega )+T(\omega )x$$

(20)

where $\omega$ is a random variable in probability space $\left(\mathsf{\Omega },F,P\right)$ with $\mathsf{\Omega }\to R^k$, $f:\mathsf{\Omega }\to R^{n_2}$, $D:\mathsf{\Omega }\to R^{m_2\times n_2}$, $T:\mathsf{\Omega }\to R^{m_2\times n_1}$, $X\subseteq R^{n_1}$, and $Y\subseteq R^{n_2}$. $C^T$, $D$, and $T$ represent the model input parameters, and $E_{\omega \in \mathsf{\Omega }}[Q(X,\omega )]={\displaystyle \sum _{h=1}^H{p}_h}Q(X,{\omega }_h)$. The first-stage decision is made of Equation (17) and parameter $x$, and it is decided before random variable $\omega$ is determined. The second-stage decision is made of Equation (20) and parameter $y$. Through taking discrete values ${\omega }_h$ with probability levels $P_h$ (h = 1, 2, …, v and ${\displaystyle \sum P_h}=1$), TSP can be transformed into a linear programming model.

In the actual decision-making process, the random disturbances, measurements and statistical errors, or incomplete data can directly affect the model’s stability and practicality. SRO is an effective method that can carry out risk aversion in the programming process and can capture robust optimization results with respect to a given disturbance range by satisfying the worst case within the range. It can also balance the trade-off between the variable random values and the expected recourse costs [38]. In addition, the parameters and information obtained in the actual planning process have a low degree of accuracy, most of which cannot be presented as probability distribution; therefore, introducing interval numbers by IPP is an efficacious approach. Thus, by integrating SRO and IPP into TSP, the model can be transformed by deriving Equations (21)–(26):

$$\mathrm{Min}{f}^{\pm }=C_{T_1}^{\pm }{X}^{\pm }+{\displaystyle \sum _{h=1}^v{P}_h{D}_{T_2}^{\pm }{Y}^{\pm }}+\omega {\displaystyle \sum _{h=1}^v{p}_h}\left(D_{T_2}^{\pm }{Y}^{\pm }-{\displaystyle \sum _{h=1}^v{p}_h{D}_{T_2}^{\pm }{Y}^{\pm }+2{\theta }_h^{\pm }}\right)$$

(21)

Subject to:

$$A_r^{\pm }{X}^{\pm }\le B_r^{\pm },r=1,2,\dots ,m_1$$

(22)

$$A_t^{\pm }{X}^{\pm }+A_t^{’\pm }{Y}^{\pm }\ge {\omega }_h^{\pm },t=1,2,\dots ,m_2;h=1,2,\dots ,v$$

(23)

$$x_j^{\pm }\ge 0,x_j^{\pm }\in X^{\pm },j=1,2,\dots ,n_1$$

(24)

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

(25)

$$\left(D_{T_2}^{\pm }{Y}^{\pm }-{\displaystyle \sum _{h=1}^v{p}_h{D}_{T_2}^{\pm }{Y}^{\pm }+2{\theta }_h^{\pm }}\right)\ge 0$$

(26)

where $A_r^{\pm }\in {\left\{R^{\pm }\right\}}^{m_1\times n_1}$, $A_t^{\pm }\in {\left\{R^{\pm }\right\}}^{m_2\times n_2}$, $B_r^{\pm }\in {\left\{R^{\pm }\right\}}^{m_1\times 1}$, $C_{T_1}^{\pm }\in {\left\{R^{\pm }\right\}}^{1\times n_1}$, $D_{T_2}^{\pm }\in {\left\{R^{\pm }\right\}}^{1\times n_2}$, $X^{\pm }\in {\left\{R^{\pm }\right\}}^{n_1\times 1}$, $Y^{\pm }\in {\left\{R^{\pm }\right\}}^{n_2\times 1}$, and $\left\{R^{\pm }\right\}$ represent a series of interval parameters; the mathematical symbol “±” denotes parameters with interval attribute values; the symbols “+” and “–“ denote the upper and lower bounds of an interval parameter, respectively; $\left(D_{T_2}^{\pm }{Y}^{\pm }-{\displaystyle \sum _{h=1}^v{p}_h{D}_{T_2}^{\pm }{Y}^{\pm }+2{\theta }_h^{\pm }}\right)$ represents the deviation between the expectation value of the second stage and the value under a specific scenario; ω represents the objective planning trade-off coefficient, which ranges from 0 to 1, reflecting the decision-maker’s attitude toward the economy and robustness of the system; and ${\theta }_h$ is a slack variable whose function is to ensure the stability and nonnegativity of the model solution. The model can be divided into two definite submodules with upper and lower bounds of the objective function, and the conversion process is completed by an interactive algorithm [39]. Therefore, $f^{-}$ is dismantled first by deriving Equations (27)–(34):

$$\begin{array}{l}\mathrm{Min}{f}^{-}={\displaystyle \sum _{j=1}^{{}_{k_1}}{c}_j^{-}{x}_j^{-}}+{\displaystyle \sum _{j=k_1+1}^{{}_{{{\displaystyle n}}_1}}{c}_j^{-}{x}_j^{+}}+{\displaystyle \sum _{j=1}^{{}_{{{\displaystyle k}}_2}}{\displaystyle \sum _{h=1}^v{p}_{jh}{d}_j^{-}{y}_{jh}^{-}}}+{\displaystyle \sum _{j=k_2+1}^{{}_{n_2}}{\displaystyle \sum _{h=1}^v{p}_{jh}{d}_j^{-}{y}_{jh}^{+}}}\\ +\omega {\displaystyle \sum _{h=1}^v{p}_h}\left[\left({\displaystyle \sum _{j=1}^{k_2}{d}_j^{-}{y}_{jh}^{-}+{\displaystyle \sum _{j=k_2}^{n_2}{d}_j^{-}{y}_{jh}^{+}}}\right)-{\displaystyle \sum _{h=1}^v{p}_h}\left({\displaystyle \sum _{j=1}^{k_2}{d}_j^{-}{y}_{jh}^{-}+{\displaystyle \sum _{j=k_2+1}^{n_2}{d}_j^{-}{y}_{jh}^{+}}}\right)+2{\theta }_h^{-}\right]\end{array}$$

(27)

Subject to:

$${\displaystyle \sum _{j=1}^{{}_{k_1}}{\left|a_{rj}\right|}^{+}\mathrm{sign}(a_{rj}^{+})x_j^{-}}+{\displaystyle \sum _{j=k_1+1}^{{}_{{{\displaystyle n}}_1}}{\left|a_{rj}\right|}^{-}\mathrm{sign}(a_{rj}^{-})x_j^{+}}\le b_r^{+},\forall r$$

(28)

$$\begin{array}{l}{\displaystyle \sum _{j=1}^{{}_{k_1}}{\left|a_{tj}\right|}^{+}\mathrm{sign}(a_{tj}^{+})x_j^{-}}+{\displaystyle \sum _{j=k_1+1}^{{}_{{{\displaystyle n}}_1}}{\left|a_{tj}\right|}^{-}\mathrm{sign}(a_{tj}^{-})x_j^{+}}+{\displaystyle \sum _{j=1}^{{}_{k_2}}{\left|a_{tj}^{’}\right|}^{+}\mathrm{sign}(a_{tj}^{’+})y_{jh}^{-}}\\ +{\displaystyle \sum _{j=k_2+1}^{{}_{{{\displaystyle n}}_2}}{\left|a_{tj}^{’}\right|}^{-}\mathrm{sign}(a_{tj}^{’-})y_{jh}^{+}}\ge {\omega }_h^{-},\forall t,h\end{array}$$

(29)

$$\left({\displaystyle \sum _{j=1}^{k_2}{d}_j^{-}{y}_{jh}^{-}+{\displaystyle \sum _{j=k_2}^{n_2}{d}_j^{-}{y}_{jh}^{+}}}\right)-{\displaystyle \sum _{h=1}^v{p}_h}\left({\displaystyle \sum _{j=1}^{k_2}{d}_j^{-}{y}_{jh}^{-}+{\displaystyle \sum _{j=k_2+1}^{n_2}{d}_j^{-}{y}_{jh}^{+}}}\right)+{\theta }_h^{-}\ge 0,\forall t,h$$

(30)

$$x_j^{-}\ge 0,j=1,2,\dots ,k_1$$

(31)

$$x_j^{+}\ge 0,j=k_1+1,k_1+2,\dots ,n_1$$

(32)

$$y_{jh}^{-}\ge 0,\forall h;j=1,2,\dots ,k_2$$

(33)

$$y_{jh}^{+}\ge 0,\forall h;j=k_2+1,k_2+2,\dots ,n_2$$

(34)

where $a_{rj}^{\pm }$, $c_j^{\pm }$, and $d_j^{\pm }$ are the parameters or the coefficients of the decision variables; $x_j^{\pm }$, j = 1, 2, …, k₁ are the interval positive variables; $x_j^{\pm }$, $j=k_1+1,k_1+2,\dots ,n_1$ are the interval negative variables; $y_{jh}^{\pm }$, $j=1,2,\dots ,k_2$ and $h=1,2,\dots ,v$ are the random positive variables; and $y_{jh}^{\pm }$, $j=k_2+1,k_2+2,\dots ,n_2$ and $h=1,2,\dots ,v$ are the random negative variables. Submodule (4) can produce solutions of $x_j^{-}{}_{opt}\left(j=1,2,\dots ,k_1\right)$, $x_j^{+}{}_{opt}\left(j=k_1+1,k_1+2,\dots ,n_1\right)$, $y_{jh}^{-}{}_{opt}\left(j=1,2,\dots ,k_2\right)$, and $y_{jh}^{+}{}_{opt}\left(j=k_1+1,k_1+2,\dots ,n_2\right)$. Thus, the second submodule $f^{-}$ is dismantled as follows by deriving Equations (35)–(42) [40]:

$$\begin{array}{l}\mathrm{Min}{f}^{+}={\displaystyle \sum _{j=1}^{{}_{k_1}}{c}_j^{+}{x}_j^{+}}+{\displaystyle \sum _{j=k_1+1}^{{}_{{{\displaystyle n}}_1}}{c}_j^{+}{x}_j^{-}}+{\displaystyle \sum _{j=1}^{{}_{{{\displaystyle k}}_2}}{\displaystyle \sum _{h=1}^v{p}_{jh}{d}_j^{+}{y}_{jh}^{+}}}+{\displaystyle \sum _{j=k_2+1}^{{}_{n_2}}{\displaystyle \sum _{h=1}^v{p}_{jh}{d}_j^{+}{y}_{jh}^{-}}}\\ +\omega {\displaystyle \sum _{h=1}^v{p}_h}\left[\left({\displaystyle \sum _{j=1}^{k_2}{d}_j^{+}{y}_{jh}^{+}+{\displaystyle \sum _{j=k_2}^{n_2}{d}_j^{+}{y}_{jh}^{-}}}\right)-{\displaystyle \sum _{h=1}^v{p}_h}\left({\displaystyle \sum _{j=1}^{k_2}{d}_j^{+}{y}_{jh}^{+}+{\displaystyle \sum _{j=k_2+1}^{n_2}{d}_j^{+}{y}_{jh}^{-}}}\right)+2{\theta }_h^{+}\right]\end{array}$$

(35)

Subject to:

$${\displaystyle \sum _{j=1}^{{}_{k_1}}{\left|a_{rj}\right|}^{-}\mathrm{sign}(a_{rj}^{-})x_j^{+}}+{\displaystyle \sum _{j=k_1+1}^{{}_{{{\displaystyle n}}_1}}{\left|a_{rj}\right|}^{+}\mathrm{sign}(a_{rj}^{+})x_j^{-}}\le b_r^{-},\forall r$$

(36)

$$\begin{array}{l}{\displaystyle \sum _{j=1}^{{}_{k_1}}{\left|a_{tj}\right|}^{-}\mathrm{sign}(a_{tj}^{-})x_j^{+}}+{\displaystyle \sum _{j=k_1+1}^{{}_{{{\displaystyle n}}_1}}{\left|a_{tj}\right|}^{+}\mathrm{sign}(a_{tj}^{+})x_j^{-}}+{\displaystyle \sum _{j=1}^{{}_{k_2}}{\left|a_{tj}^{’}\right|}^{-}\mathrm{sign}(a_{tj}^{’-})y_{jh}^{+}}\\ +{\displaystyle \sum _{j=k_2+1}^{{}_{{{\displaystyle n}}_2}}{\left|a_{tj}^{’}\right|}^{+}\mathrm{sign}(a_{tj}^{’+})y_{jh}^{-}}\ge {\omega }_h^{+},\forall t,h\end{array}$$

(37)

$$\left({\displaystyle \sum _{j=1}^{k_2}{d}_j^{+}{y}_{jh}^{+}+{\displaystyle \sum _{j=k_2}^{n_2}{d}_j^{+}{y}_{jh}^{-}}}\right)-{\displaystyle \sum _{h=1}^v{p}_h}\left({\displaystyle \sum _{j=1}^{k_2}{d}_j^{+}{y}_{jh}^{+}+{\displaystyle \sum _{j=k_2+1}^{n_2}{d}_j^{+}{y}_{jh}^{-}}}\right)+{\theta }_h^{+}\ge 0,\forall t,h$$

(38)

$$x_j^{+}\ge x_j^{-}{}_{opt},j=1,2,\dots ,k_1$$

(39)

$$0\le x_j^{-}\le x_j^{+}{}_{opt},j=k_1+1,k_1+2\dots n_1$$

(40)

$$y_{jh}^{+}\ge y_{jh}^{-}{}_{opt},\forall h;j=1,2,\dots ,k_2$$

(41)

$$y_{jh}^{+}{}_{opt}\ge y_{jh}^{-}\ge 0,\forall h;j=k_2+1,k_2+2,\dots ,n_2$$

(42)

The solutions of $x_j^{+}{}_{opt}\left(j=1,2,\dots ,k_1\right)$, $x_j^{-}{}_{opt}\left(j=k_1+1,k_1+2,\dots ,n_1\right)$, $y_{jh}^{+}{}_{opt}\left(j=1,2,\dots ,k_2\right)$, and $y_{jh}^{-}{}_{opt}\left(j=k_1+1,k_1+2,\dots ,n_2\right)$ can be produced by the second submodule. Thus, by reintegrating the solutions of the two submodules, an interval solution with upper and lower bound forms for ITSP–SRO can be obtained.

2.3. Technique for Order of Performance by Similarity to Ideal Solution

Each optimization scheme has multiple attribute characteristics, and post-optimization analysis can capture and reflect critical factors through the multi-criteria decision-making (MCDM) approach. TOPSIS can be applied to rank the evaluation of objects according to their distance from the idealized solution. Moreover, the method assists decision-makers in identifying key criteria, evaluating the optimization results, and sorting the alternative schemes. The procedure of the TOPSIS method is expressed as follows [41]:

$$R=\begin{array}{c}{A}^1\\ A^2\\ ⋮\\ A^m\end{array}\left[\begin{array}{cccc}{x}_{11}& x_{12}& \dots & x_{1n}\\ x_{21}& x_{22}& \dots & x_{2n}\\ ⋮& ⋮& ⋮& ⋮\\ x_{m1}& x_{m2}& \dots & x_{mn}\end{array}\right]$$

(43)

Formulate the decision matrix R, where $A^i$ denotes alternative schemes with $i=1,2,\dots ,m;$ $j=1,2,\dots ,n.$ $x_{ij}$ denotes the performance of attribute $X_j$.

$$n_{ij}=x_{ij}/{({\displaystyle \sum _{i=1}^m{x}_{ij}^2})}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$$

(44)

Normalize the matrix R, denoting T and standardizing the attributes of each dimension in the decision matrix, with $i=1,2,\dots ,m;j=1,2,\dots ,n.$

$$WT=\left[\begin{array}{cccc}{w}_1{n}_{11}& w_2{n}_{12}& \dots & w_n{n}_{1n}\\ w_1{n}_{21}& w_2{n}_{22}& \dots & w_n{n}_{2n}\\ ⋮& ⋮& ⋮& ⋮\\ w_1{n}_{m1}& w_2{n}_{m2}& \dots & w_n{n}_{mn}\end{array}\right]=\left[\begin{array}{cccc}{r}_{11}& r_{12}& \dots & r_{1n}\\ r_{21}& r_{22}& \dots & r_{2n}\\ ⋮& ⋮& ⋮& ⋮\\ r_{m1}& r_{m2}& \dots & r_{mn}\end{array}\right]$$

(45)

Based on Equations (43)–(45),calculate the weighted matrix of WT, where $w_j$ is the weight of the $j$ attribute and ${\displaystyle \sum _{j=1}^n{w}_j}=1$. Then, calculate the identification of the positive/negative idealized solutions.

The positive idealized solution can be determined as Equation (46):

$$r^{+}=\{(\max r_{ij}\left|j\in J),(\min r_{ij}\right|j\in J^{’})\}=\{r_1^{+},r_2^{+},\dots ,r_j^{+},\dots ,r_n^{+}\left|i=1,2,\dots ,m\right.\}$$

(46)

The negative idealized solution can be determined as Equation (47):

$$r^{-}=\{(\min r_{ij}\left|j\in J),(\max r_{ij}\right|j\in J^{’})\}=\{r_1^{-},r_2^{-},\dots ,r_j^{-},\dots ,r_n^{-}\left|i=1,2,\dots ,m\right.\}$$

(47)

where J and J’ denote the benefit criteria and the cost criteria, respectively.

The distance of each alternative scheme from the positive idealized solution can be determined as Equation (48):

$$d_i^{+}={[{\displaystyle \sum _{j=1}^n{(r_{ij}-r_j^{+})}^2}]}^{1/2}.i=1,2,\dots ,m;j=1,2,\dots ,n.$$

(48)

The distance of each alternative scheme from the negative idealized solution can be determined as Equation (49):

$$d_i^{-}={[{\displaystyle \sum _{j=1}^n{(r_{ij}-r_j^{-})}^2}]}^{1/2}.i=1,2,\dots ,m;j=1,2,\dots ,n.$$

(49)

The relative closeness $C_i^{*}$ is defined as:

$$C_i^{*}=d_i^{-}/(d_i^{+}+d_i^{-}).i=1,2,\dots ,m.$$

(50)

According to the relative closeness of each scheme calculated by Equation (50), the candidate schemes are ranked in descending order, and the optimal scheme is determined by comparing the main influence degree of each factor.

3. Case Study

3.1. Overview of the Study Area

Harbin City (125°42′–130°10′ E, 44°04′–46°40’ N), situated in the northeast of China, is the capital of Heilongjiang Province and a large central metropolis in the region of Northeast Asia. The city covers an area of 53,100 km², ranks first among the provincial-level cities in China, and had a permanent population of 10.85 million in 2018 (Figure 2). As the major industrial base of China, Harbin has a strong foundation, forming an industrial system composed of cars, food, medicine, heavy machinery, petrochemical equipment, power generation equipment, and electromechanical equipment. In 2018, the regional gross domestic product was RMB 630.05 billion yuan, and the secondary industry increased by RMB 168.93 billion yuan from the previous year. Referring to the Statistics Bureau of Harbin 2002 to 2018, the energy consumption of industrial enterprises above the designated size amounted to 17.96 million tons of coal equivalent in 2018, and this was more than double the comparable figure over 2002. Coal consumption accounts for an enormous proportion of total energy consumption, which has continuously exceeded 65% of the total local energy consumption over the past decade. Considering the high speed of economic growth, energy demand has also increased substantially, which has definitely impacted on the electric power system of the city. Currently, the capacity of domestic electricity production in Harbin cannot satisfy the city’s rising demands. In 2017, the total electricity consumption reached 22.18 × 10³ GWh, and in order to make up the shortage caused by the poor productivity, importing electricity was utilized, which amounted to 5.58 × 10³ GWh.

With the continuous implementation of supply-side energy structure adjustment, the energy end-use system of Harbin is moving toward a health mode. However, the energy consumption structure is relatively fixed, and electricity generation still relies heavily on coal. Although the proportion of coal in the total power production systems shows a slight declining trend, the actual consumption amount of coal is gradually increasing. By the end of 2017, the proportion of coal power generation dropped to 87.83% from 92.61% of the year 2013. In addition, the amount of coal consumption in the thermal power plant industry surged to 34.48 million tons, and led to a tremendous potential hazard of atmospheric environmental issues in the city. For instance, the total emission amounts of SO₂ and NO_x in 2017, respectively, increased to 103.60 thousand tons and 64.80 thousand tons, and Harbin became one of the heaviest-polluted cities in Northeast China. Faced with the above problems, a number of positive and effective measures aiming at reducing air pollutants emission have been adopted, such as boosting the utilization of clean energy, enlarging the percentage of imported electricity, and the enhancing environmental legal enforcement. In addition, Harbin has an immense potential in the development and utilization of renewable energy resources, such as wind energy, hydro energy, solar energy, and bio-energy. Harbin is qualified for the conditions of building wind farms, with the average wind speed distributed within 2.4–5.0 m/s in the mountainous region along Songhua River. The city is rich in hydro power resources: the amount of river with a watershed area greater than 4 km² is up to 186 in the area, including Ashihe River, Yunlianghe River, and Hulan River. Solar energy resources have a broad perspective of utilization and development, and the average annual solar radiation in Harbin is up to 5040 MJ/m². Moreover, the available amount of crop straw is more than 20,630 thousand tons per year. According to the energy development strategy of Harbin Municipal Development and Reform Commission, 24 biomass power generation projects will be built in the period of 2018–2020.

In general, decision-makers still face enormous challenges to optimize the structure of electric power systems, to ensure power supply, and to satisfy air pollution reduction requirements. Electric power system management tends to be dynamic, multi-objective, and complicated [42,43,44]. Moreover, it involves multiple uncertain factors that can hardly be expressed as deterministic values, such as electricity consumption, economic cost, and resource availability. Such inherent uncertainties can interfere with the decision schemes and optimization processes. Thus, the purpose of this study is to formulate a comprehensive research framework for realizing the whole management cycle. In detail, the approaches proposed in this study can effectively (a) forecast the electricity consumption based on limited historical data characterized by stochasticity, and the forecasting results are presented in the form of interval parameters with corresponding random probability; (b) make stable and reliable power generation programs for multiple conversion technologies under electric power structure adjustment plans and regional pollutants emission reduction policies; and (c) identify the most influential factors of planning decisions by comprehensively analyzing and generating optimum electricity generation schemes.

3.2. Electricity Consumption Predictions

In this study, two planning periods were considered, the first was 2019–2023, and the second was 2024–2028. Future electricity consumption has a wide range of uncertainties, and the way in which to transform it into stochastic parameters with probability distribution is key to improving the accuracy of the optimization results. The GMM can effectively tackle the uncertainties in the system, and can describe the forecasted results as interval values associated with probabilities. According to the procedure of the GMM, made explicit in the Methodology section, the future electricity consumption of Harbin based on data from 2009 to 2018 was predicted in this study (

Table 1. Time series variations of electricity consumption for Harbin (10³ GWh).

k	1	2	3	4	5	6	7	8	9	10
Year	2009	2010	2011	2012	2013	2014	2015	2016	2017	2018
X(0)(k)	161.6	173.8	185.6	191.4	197.4	200.3	204.8	215.3	221.8	238.7
X(1)(k)	161.6	335.5	521.1	712.6	910.1	1110.5	1315.3	1530.6	1752.4	1991.1

).

Based on Equations (1)–(8) and the historical data, the associated coefficients in GM (1,1) were obtained: a = −0.034887 and u = 167.36. Therefore, Harbin’s electricity consumption GM (1,1) forecast model can be expressed as:

$${\widehat{x}}^{(1)}(k+1)=\left[x^{(0)}(1)-\frac{b}{a}\right]\cdot e^{-ak}+\frac{b}{a}=4958.969e^{0.034887k}-4797.36,k=0,1,2,\dots ,n.$$

(51)

The prediction trend curves for electricity consumption produced by GM (1,1) are shown in Figure 3. Remarkably, the time series of total electricity consumption in Harbin increased irregularly around the trend curves, and the fluctuations interfered with the forecasting accuracy. Following Equation (8), the relative errors were distributed in [−2.349%, 2.498%]. Based on the specific state division method of the Gray–Markov forecasting approach, the ranges of the relative error series from GM (1,1) can be divided into contiguous intervals, [−3%, −1%], [−1%, 1%], and [1%, 3%], treated as identifiable intervals of state E(1), E(2), and E(3), respectively. The relative errors of electricity consumption prediction and the states of each year for 2009–2018 are represented in

Table 2. Time series variations of electricity consumption for Harbin (10³ GWh).

k	1	2	3	4	5	6	7	8	9	10
Year	2009	2010	2011	2012	2013	2014	2015	2016	2017	2018
ζ(k)(%)	0	−1.247	−1.821	1.409	1.016	−1.020	−2.349	−0.814	−1.334	2.498
E(i)	$E(2)$	$E(1)$	$E(3)$	$E(3)$	$E(3)$	$E(1)$	$E(1)$	$E(2)$	$E(1)$	$E(3)$

.

Subsequently, the transition probability matrix of each state can be obtained by Equations (12) and (13), and the transition matrix is as follows:

$$W=\left[\begin{array}{ccc}0.25& 0.25& 0.5\\ 1& 0& 0\\ 0.333& 0& 0.667\end{array}\right]$$

(52)

The future predictive state vectors can be generated by the m steps transition matrix, and the related probabilities can thus be obtained. Since each state is defined by relative errors expressed as interval parameters, and the valve of the prediction results corresponding to various states can be similarly determined as interval parameters. The interval values of electricity consumption have the advantage of reducing the risk by addressing the variability and uncertainties when they are regarded as input parameters of the optimization model. Thus, the forecasting results of electricity consumption during the planning horizon are presented in

Table 3. The forecasted values of electricity consumption and the associated probabilities in 2019–2028.

Gray–Markov Model (GMM) Forecasting Data
Year	Prediction Interval (103 GWh)	Probability	Year	Prediction Interval (103 GWh)	Probability
2019	[237.44, 239.80]	0.444	2024	[282.68, 285.50]	0.363
	[239.80, 242.21]	0.112		[285.50, 288.37]	0.092
	[242.21, 244.67]	0.444		[288.37, 291.29]	0.545
2020	[245.86, 248.31]	0.370	2025	[292.72, 295.63]	0.364
	[248.31, 250.81]	0.112		[295.63, 298.61]	0.091
	[250.81, 253.85]	0.518		[298.61, 301.64]	0.545
2021	[254.59, 257.13]	0.376	2026	[303.11, 306.13]	0.364
	[257.13, 259.71]	0.093		[306.13, 309.21]	0.091
	[259.71, 262.35]	0.531		[309.21, 312.35]	0.545
2022	[263.63, 266.26]	0.364	2027	[313.88, 317.01]	0.364
	[266.26, 268.93]	0.094		[317.01, 320.18]	0.091
	[268.93, 271.66]	0.542		[320.18, 323.44]	0.545
2023	[272.99, 275.71]	0.365	2028	[325.02, 328.25]	0.364
	[275.71, 278.48]	0.091		[328.25, 331.55]	0.091
	[278.48, 281.31]	0.544		[331.55, 334.92]	0.545

.

In addition, the prediction interval values of states E(1), E(2), and E(3) showed a growing trend during the planning period; thus, the three states could form three contiguous intervals, which correspond to the low, medium, and high electricity demand levels, respectively. Considering the existence of the limiting transition probability vector in the Markov process, the values of the vector can reflect the probabilities of each state occurrence when the transition reaches stability. Therefore, the limiting transition probability vector can be used to reflect the probability of occurrence for various demand levels within the forecasting period, and the vector can be calculated as follows:

$${\displaystyle \prod =[{\pi }_1(k)},{\pi }_2(k),{\pi }_3(k)]=[\underset{k\to \infty }{\lim }{\pi }_1(k),\underset{k\to \infty }{\lim }{\pi }_2(k),\underset{k\to \infty }{\lim }{\pi }_3(k)]=\left[0.364,0.091,0.545\right].$$

(53)

In general,

Table 4. Time series variations of electricity consumption for Harbin (10³ GWh).

Demand Level	Probability	Electricity Demand (103 GWh)
		Period 1 (2019–2023)	Period 2 (2024–2028)
Low	36.4%	[1274.54, 1287.22]	[1517.43, 1532.53]
Medium	9.1%	[1287.22, 1300.16]	[1532.53, 1547.94]
High	54.5%	[1300.16, 1313.36]	[1547.94, 1563.65]

shows the amount of regional electricity demand and the corresponding occurrence probabilities. The results clearly indicate that the GMM model has enormous potential in electricity consumption prediction and can generate feasible input parameters of the optimization model.

3.3. Optimization Model Formulation

In this study, the ITSRP model was applied to Harbin’s power system management to demonstrate its advantages in addressing multiple complexities and uncertainties. The objective function of the model is to minimize the total cost concerning different power generation technologies. In detail, the items of system cost and investment contain purchasing primary energy resources, importing power, capacity expansions, electricity conversion, and controlling emissions (such as $S{O}_2$, $N{O}_x$, and $P{M}_{}$). Thus, the model can be formulated as follows:

$$Min{f}^{\pm }=(1)+(2)+(3)+(4)+(5)$$

(1)

The cost of supplying resources and importing electric power

$${\displaystyle \sum _{i=1}^I{\displaystyle \sum _{t=1}^T(PE{C}_{it}^{\pm }Z{1}_{it}^{\pm }})}+{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{h=1}^H{p}_{th}PI{E}_t^{\pm }}}Z{2}_{th}^{\pm }$$

(54)

(2)

The cost of electricity conversion

$${\displaystyle \sum _{i=1}^I{\displaystyle \sum _{t=1}^T\left(P{V}_{it}^{\pm }{W}_{it}^{\pm }\right)}}+{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{h=1}^H{p}_{th}\left(P{V}_{it}^{\pm }+P{P}_{it}^{\pm }\right)Q_{ith}^{\pm }}}}$$

(55)

(3)

The cost of capacity expansion

$${\displaystyle \sum _{i=1}^I{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{h=1}^H{p}_{th}(A_{it}^{\pm }{Y}_{ith}^{\pm }+B_{it}^{\pm }{X}_{ith}^{\pm })}}}$$

(56)

(4)

The cost of atmospheric pollutants treatment

$${\displaystyle \sum _{i=1}^I{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{j=1}^J{W}_{it}^{\pm }\cdot V_{itj}^{\pm }}}}\cdot {\lambda }_{itj}^{\pm }\cdot P{C}_{itj}^{\pm }+{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{h=1}^H{\displaystyle \sum _{j=1}^J{p}_{th}\cdot Q_{ith}^{\pm }\cdot V_{itj}^{\pm }}}}}\cdot {\lambda }_{itj}^{\pm }\cdot P{D}_{itj}^{\pm }$$

(57)

(5)

Robust function

$$\omega {\displaystyle \sum _{i=1}^I{\displaystyle \sum _{t=1}^T{\displaystyle \sum _{h=1}^H{p}_{th}}}}\left({\xi }_{ith}^{\pm }-{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{h=1}^H{p}_{th}{\xi }_{ith}^{\pm }+2{\theta }_{ith}^{\pm }}}\right)$$

(58)

$${\xi }_{ith}^{\pm }=PI{E}_t^{\pm }Z{2}_{th}^{\pm }+\left(P{V}_{it}^{\pm }+P{P}_{it}^{\pm }\right)Q_{ith}^{\pm }+\left(A_{it}^{\pm }{Y}_{ith}^{\pm }+B_{it}^{\pm }{X}_{ith}^{\pm }\right)+{\displaystyle \sum _{j=1}^J{Q}_{ith}^{\pm }\cdot V_{itj}^{\pm }}\cdot {\lambda }_{itj}^{\pm }\cdot P{D}_{itj}^{\pm }$$

(59)

Subject to:

(1)

Constraints for primary energy and control of total coal consumption amount

$$(W_{1t}^{\pm }+Q_{1th}^{\pm })F{E}_{1t}^{\pm }\le Z{1}_t^{\pm }\le UP{C}_t^{\pm },\forall t,h;$$

(60)

(2)

Constraints for renewable energy

$$(W_{2t}^{\pm }+Q_{2th}^{\pm })F{E}_{2t}^{\pm }\le UP{W}_t^{\pm },\forall t,h;$$

(61)

$$(W_{3t}^{\pm }+Q_{3th}^{\pm })F{E}_{3t}^{\pm }\le UP{H}_t^{\pm },\forall t,h;$$

(62)

$$(W_{4t}^{\pm }+Q_{4th}^{\pm })F{E}_{4t}^{\pm }\le UP{S}_t^{\pm },\forall t,h;$$

(63)

$$(W_{5t}^{\pm }+Q_{5th}^{\pm })F{E}_{5t}^{\pm }\le UP{B}_t^{\pm },\forall t,h;$$

(64)

$$(W_{6t}^{\pm }+Q_{6th}^{\pm })F{E}_{6t}^{\pm }\le UP{G}_t^{\pm },\forall t,h;$$

(65)

(3)

Constraints for electricity demand and energy structure adjustment

$${\displaystyle \sum _{i=1}^I(W_{it}^{\pm }+Q_{ith}^{\pm }+Z{2}_{th}^{\pm })\ge D_{th}^{\pm }},\forall t,h;$$

(66)

$$(W_{1t}^{\pm }+Q_{1th}^{\pm })\le ({\mu }_t-\alpha )D_{th}^{\pm },\forall t,h;$$

(67)

(4)

Constraints for limited value of capacity expansion

$$Y_{ith}^{\pm }\left\{\begin{array}{cc}=1\hfill & ifcapacityexpansionundertaken\\ =0\hfill & ifotherwise\end{array}\right.,\forall i,t,h;$$

(68)

$$N_{it}\le X_{ith}^{\pm }\le M_{it}{Y}_{ith}^{\pm },\forall i,t,h;$$

(69)

(5)

Constraints for power generation technology

$$W_{it}^{\pm }+Q_{ith}^{\pm }\le S{T}_{it}^{\pm }\cdot \left\{R{C}_i+{\displaystyle \sum _{t^{’}=1}^t{X}_{ith}^{\pm }}\right\},\forall i,t,h;$$

(70)

$$W_{it}^{\pm }\ge Q_{ith}^{\pm }\ge 0,\forall i;t;h=1,2,\dots ,H_t$$

(71)

(6)

Constraints for atmospheric environmental capacity

$${\displaystyle \sum _{i=1}^I(W_{it}^{\pm }+Q_{ith}^{\pm })}\cdot (1-{\lambda }_{itj}^{\pm })\cdot V_{itj}^{\pm }\le (1-\beta )\cdot E{P}_{tj}^{\pm },\forall t,h,j;$$

(72)

(7)

Constraints for technology

$$Z{1}_{th}^{\pm },Z{2}_{th}^{\pm },W_{it}^{\pm },Q_{ith}^{\pm }\ge 0,\forall i,t,h;$$

(73)

(8)

Robust constraint

$${\xi }_{ith}^{\pm }-{\displaystyle \sum _{i=1}^I{\displaystyle \sum _{h=1}^H{p}_{th}{\xi }_{ith}^{\pm }+2{\theta }_{ith}^{\pm }}}\ge 0,\forall i,t,h;$$

(74)

In this study, the involved technical–economic–environmental data were obtained from Harbin governmental reports and municipal development plans (www.hlj.stats.gov.cn). The specific definitions for the variables and parameters are shown in Appendix

Table A1. Definition for variables and parameters of optimization model.

Variables/Parameters	Specific Physical Meaning
f	Total cost during the planning horizon
t	Planning period, t = 1 is period 1st, t = 2 is period 2nd
i	Different type of power conversion technology, i = 1 for thermal power, i = 2, for wind power, i = 3, for hydro power, i = 4, for solar power, i = 5, for biomass power, i = 6, for garbage power
h	Electricity demand-level, h = 1 for the low demand-level, h =2 for the medium demand-level, h =3 for the high demand-level
j	type of air pollutant, j = 1 for SO2, j = 2 for NOx, and j = 3 for particulate matter (PM)
PECit	Price of energy resource i in period t
Z1it	Consumption amount of energy resource under conversion technology i in period t
pth	Probability of demand-level h occurrence in period t
PIEt	Purchased electricity prices in period t
Z2th	Amount of electricity purchased under demand-level h occurrence in period t
PVit	Operating cost of power conversion technology i for pre-regulated electricity generation in period t
Wit	Pre-regulated electricity generation target of power conversion technology i which is promised to end-users during period t
Qith	Excess electricity generation of power conversion technology i by which electricity generation target is exceeded when electricity demand-level is h in period t
PPit	Operating cost of power conversion technology i for excess electricity generation in period t
Ait	Fixed-charge cost for capacity expansion of power conversion technology i in period t
Bit	Variable cost for capacity expansion of power conversion technology i in period t
Yith	Binary variables for identifying whether or not a capacity expansion action of power conversion technology i needs to be undertaken when electricity demand-level is h in period t
Xith	Continuous variables about the amount of capacity expansion of power conversion technology i when electricity demand-level is h in period t
Vitj	Discharge intensity of pollutant j from power generation technology i in period t
λitj	Average removal efficiency of pollutant j from power conversion technology i in period t
PCitj	Removal cost of pollutant j treatment for conversion technology i in period t
PDitj	Penalty cost of excess pollutant j treatment for conversion technology i in period t
ω	Weight coefficient corresponding to the risk-aversion levels, different levels would produce various output solutions
θith	The slack variable
ξith	The robust variable
FEit	Conversion efficiency of power conversion technology i in period t
UPCt	Amount of available coal resources in period t
UPWt	Amount of available wind resources in period t
UPHt	Amount of available hydro resources in period t
UPSt	Amount of available solar resources in period t
UPBt	Amount of available biomass resources in period t
UPGt	Amount of available garbage resources in period t
Dth	Total electricity demand under demand-level h during period t
μt	Planning target of thermal power proportion of total electricity generation during period t
Nit	Variable lower bound for capacity expansion of power conversion technology i in period t
Mit	Variable upper bound for capacity expansion of power conversion technology i in period t
STit	Unit operating time for power conversion technology i in period t
RCi	Residual capacity of conversion technology i
EPtj	Upper limit amount of pollutant j in period t
α	The thermal power scale coefficient
β	The environmental permission coefficient

, where $\omega$ denotes the objective planning trade-off weight coefficient, which ranges from 0 to 1, and different levels produce various output solutions; $\omega =1$ is considered in this study to mean that decision-makers hold an anti-risk attitude, and the economy can be guaranteed on the basis of the system’s robustness. To verify the energy structure adjustment and the atmospheric pollutants emission reduction target is critical to regional electric power system management planning; thus, a series of scenarios with different thermal power proportions and pollutants reduction levels were designed. Furthermore, $\alpha$ denotes the thermal power scale coefficient, and $\beta$ denotes the environmental permission coefficient. Ten representative scenarios were generated by combining the different values of $\alpha$ and $\beta$, and the specific descriptions of the fetching values are presented in

Table 5. Descriptions of the fetching values for scenarios 1–10.

Scenario	S1	S2	S3	S4	S5	S6	S7	S8	S9	S10
$\alpha$	0	1%	1%	1%	2%	2%	2%	3%	3%	3%
$\beta$	0	2%	4%	6%	2%	4%	6%	2%	4%	6%

Note: $\alpha$ denotes the thermal power scale coefficient; $\beta$ denotes the environmental permission coefficient. Scenarios 1–10 were established by using interactive combinations of the variable assignments of $\alpha$ and $\beta$.

. Scenario 1 is the base thermal power scale and pollutants emission scenario, in which the required thermal power proportion target and atmospheric environmental capacity limitation for the electric industry refer to those in the Program Outline for Medium—and Long-Term Development of Energy in Harbin City. The other scenarios were established by using interactive combinations of variable assignments of $\alpha$ and $\beta$, corresponding to $\alpha =${1%, 2%, 3%} and $\beta =${2%, 4%, 6%}, respectively. Therefore, an approach based on multi-scenario analysis is very effective, leading to a trade-off between power supply security and environmental quality improvement. In addition, it is also a preparatory work for choosing the optimal scheme by applying TOPSIS to make comparisons among the subsequent scenarios.

4. Results Analysis

4.1. Electricity Generation Schemes

Figure 4 describes the optimal power generation scheme without considering further energy structure adjustment and pollutants emission reduction scenarios (e.g., scenario 1) during the whole planning horizon. From period 1 to 2, the pre-regulated electricity produced by different generation modes grows, and the vast majority of electricity is supplied by thermal power units. For instance, the pre-regulated electricity generated by thermal power at the secondary demand level reaches 93.04 × 10³ GWh and 94.91 × 10³ GWh in the two periods, which occupy 77.76% and 75.61% of the total amount, respectively. Moreover, with intensified renewable power development, the amount of pre-regulated electricity produced by renewable sources assumes an obvious increasing trend. In detail, for wind power, hydro power, solar power, biomass power, and garbage power, the corresponding pre-regulated targets grow from 8.08 × 10³ GWh, 0.40 × 10³ GWh, 0.13 × 10³ GWh, 10.00 × 10³ GWh, and 8.00 × 10³ GWh in period 1 to 9.00 × 10³ GWh, 0.45 × 10³ GWh, 0.16 × 10³ GWh, 11.00 × 10³ GWh, and 10.00 × 10³ GWh in period 2, respectively. In general, based on regional massive coal reserves and low-cost running technology, thermal power would have a dominant position in the power system for a long time.

Optimized electricity power generation is a combination of pre-regulated generation and excess generation. Once the pre-regulated electricity targets can no longer meet the random demand, the action for excess generation is activated. With the increase of electricity demand, the excess generation changes to some extent, and contributes to the variability of optimized electricity power generation. Specifically, the optimized generation amount of thermal power is [93.04, 97.79] × 10³ GWh, [102.86, 103.87] × 10³ GWh, and [93.96, 98.73] × 10³ GWh under the three progressive demand levels in period 1. For biomass power in period 2, the optimized generation amount is [21.73, 22.00] × 10³ GWh under the low level, [15.04, 15.39] × 10³ GWh under the medium level, and 22.00 × 10³ GWh under the high level, respectively. Furthermore, the optimized generation of wind power, hydro power, solar power, and garbage power all have a slight growth benchmarked on their pre-regulated targets during the planning horizon. Conclusively, this reveals that in the case of an electricity shortage, thermal power would be crucial as the first choice to compensate the deficits, and the other generation modes would be supplements.

Table 6. The optimized power generation schemes of the different scenarios.

Planning Period	Demand Level	Conversion Technology	Scenario
			S1	S2	S3	S4	S5
t = 1	Low	thermal	[93.04,97.79]	[91.76,92.67]	[91.76,92.67]	[91.76,92.67]	[90.49,91.39]
		wind	8.07	8.07	8.07	8.07	8.07
		hydro	0.8	0.8	0.8	0.8	0.8
		solar	0.139	0.139	0.139	0.139	0.139
		biomass	[17.39,17.74]	[18.67,19.02]	[18.67,19.02]	[18.67,19.02]	[19.94,20]
		garbage	8	8	8	8	8
	Medium	thermal	[102.8,103.8]	[99.04,104.6]	[99.04,105.4]	[99.04,106.2]	[99.04,104.6]
		wind	8.07	8.07	8.07	8.07	8.07
		hydro	0.8	0.8	0.8	0.8	0.8
		solar	0.139	0.139	0.139	0.139	0.139
		biomass	18	18	18	18	18
		garbage	16	16	16	16	16
	High	thermal	[93.96,98.73]	[92.67,93.61]	[92.67,93.61]	[92.67,93.61]	[91.39,92.31]
		wind	8.07	8.07	8.07	8.07	8.07
		hydro	0.8	0.8	0.8	0.8	0.8
		solar	0.139	0.139	0.139	0.139	0.139
		biomass	[17.74,18.09]	[19.02,19.39]	[19.02,19.39]	[19.02,19.39]	20
		garbage	8	8	8	8	8
t = 2	Low	thermal	[111.1,116.9]	[110.3,111.4]	[110.3,111.4]	[110.3,111.4]	[108.8,109.9]
		wind	9	9	9	9	9
		hydro	0.85	0.85	0.85	0.85	0.85
		solar	0.163	0.163	0.163	0.163	0.163
		biomass	[21.73,22]	22	22	22	22
		garbage	10	[10.84,11.27]	[10.84,11.27]	[10.84,11.27]	[12.37,12.82]
	Medium	thermal	[94.91,99.77]	[93.61,94.56]	[93.61,94.56]	[93.61,94.56]	[92.31,93.24]
		wind	9	9	9	9	9
		hydro	0.9	0.9	0.9	0.9	0.9
		solar	0.163	0.163	0.163	0.163	0.163
		biomass	[15.04,15.39]	[16.34,16.71]	[16.34,16.71]	[16.34,16.71]	[17.64,18.02]
		garbage	10	10	10	10	10
	High	thermal	[111.1,116.9]	[110.3,111.4]	[110.3,111.4]	[110.3,111.4]	[108.8,109.9]
		wind	9	9	9	9	9
		hydro	0.9	0.9	0.9	0.9	0.9
		solar	0.163	0.163	0.163	0.163	0.163
		biomass	22	22	22	22	22
		garbage	10	10	10	10	10
Planning Period	Demand Level	Conversion Technology	Scenario
			S6	S7	S8	S9	S10
t = 1	Low	thermal	[90.49,91.39]	[90.49,91.39]	[89.2,90.10]	[89.2,90.10]	[89.2,90.10]
		wind	8.07	8.07	8.07	8.07	8.07
		hydro	0.8	0.8	0.8	0.8	0.8
		solar	0.139	0.139	0.139	0.139	0.139
		biomass	[19.94,20]	[19.44,20]	20	20	[19.55,20]
		garbage	8	8	8	8	8
	Medium	thermal	[99.04,105.4]	[99.04,106.2]	[99.04,104.6]	[99.04,105.4]	[97.39,99.04]
		wind	8.07	8.07	8.07	8.07	[8.07,16.15]
		hydro	0.8	0.8	0.8	0.8	0.8
		solar	0.139	0.139	0.139	0.139	0.139
		biomass	18	18	18	18	18
		garbage	16	16	16	16	16
	High	thermal	[91.39,92.31]	[91.39,92.31]	[90.10,91.01]	[90.10,91.01]	[90.10,91.01]
		wind	8.07	8.07	8.07	8.07	8.07
		hydro	0.8	0.8	0.8	0.8	0.8
		solar	0.139	0.139	0.139	0.139	0.139
		biomass	20	[19.35,20]	20	20	[19.47,20]
		garbage	8	8	8	8	8
t = 2	Low	thermal	[108.8,109.9]	[108.8,109.9]	[107.2,108.3]	[107.2,108.3]	[107.2,108.3]
		wind	9	9	9	9	9
		hydro	0.85	0.85	0.85	0.85	0.85
		solar	0.163	0.163	0.163	0.163	0.163
		biomass	22	22	22	22	22
		garbage	[12.37,12.82]	[12.04,12.82]	[13.91,14.37]	[13.25,14.37]	[12.20,14.37]
	Medium	thermal	[92.31,93.24]	[92.31,93.24]	[91.01,91.93]	[91.01,91.93]	[91.01,91.93]
		wind	9	9	9	9	9
		hydro	0.9	0.9	0.9	0.9	0.9
		solar	0.163	0.163	0.163	0.163	0.163
		biomass	[17.64,18.02]	[17.64,18.02]	[18.94,19.33]	[18.94,19.33]	[18.94,19.33]
		garbage	10	10	10	10	10
	High	thermal	[108.8,109.9]	[108.8,109.9]	[107.2,108.3]	[107.2,108.3]	[107.2,108.3]
		wind	9	9	9	9	9
		hydro	0.9	0.9	0.9	0.9	0.9
		solar	0.163	0.163	0.163	0.163	0.163
		biomass	22	22	22	22	22
		garbage	10	10	10	10	10

displays the optimized electricity power generation schemes under the ten typical scenarios for different generation modes during the whole period. Generally, the amount of traditional and renewable power generation varies as the energy structure and pollutants emission limitation change. As $\alpha$ increases, the generation amount of thermal power shows a declining trend. For example, the amount of thermal power generation is [91.76, 92.67] × 10³ GWh, [90.49, 91.39] × 10³ GWh, and [89.2, 90.10] × 10³ GWh at the low level during period 1 in scenarios 2, 5, and 8, respectively. On the contrary, under different demand levels, the generation amount of biomass power is raised as $\alpha$ increases. In period 2, at the medium level, the generation amount of biomass power is [16.34, 16.71] × 10³ GWh, [17.64, 18.02] × 10³ GWh, and [18.94, 19.33] × 10³ GWh in scenarios 3, 6, and 9, respectively. However, during the planning horizon, the optimized generation amounts of the other renewable powers (e.g., wind, hydro, solar, and garbage) are approximately constant under different demand levels. This indicates that biomass power would be a significant electricity supplement for Harbin when there is a decrease of the thermal power generation ratio. However, other conversion technologies have evolved slowly as resource constraints and technology limitations, and government managers facilitate new energy construction and promote renewable energy industry development by making policy changes and encouraging measures. In addition, the amount of power generation for each mode during the planning period has only a very slight variation as the $\beta$ value changes. Comparing the two policies, energy structure adjustment has a greater impact on regional optimized electricity power generation schemes than pollution reduction. Thus, implementing energy structure adjustment is the most efficient method to promote the healthy and stable development of the electricity system of Harbin in the long term.

4.2. Capacity Expansion Schemes

Figure 5 describes the schemes of capacity expansion for each power generation mode under different scenarios during the whole planning period. Generally, in the case that the current power production facilities cannot meet the continuously growing demand, a series of capacity expansion schemes can be adopted to alleviate power shortage crises. For different scenarios, thermal power, hydro power, and solar power have no need for expansion when considering energy structure adjustment and pollutants emission reduction. However, due to the advantages of resource costs, renewable capability and the increasing exploitation efforts of new energy, biomass power, garbage power, and wind power can be expanded adequately in the planning horizon. Under scenario 1, the total expansion capacity reaches [0.68, 0.69] $\mathrm{GW}$, [0.78, 0.99] $\mathrm{GW}$, and [0.69, 0.70] $\mathrm{GW}$ under different demand levels in period 1, respectively, corresponding to [0.86, 0.87] $\mathrm{GW}$, [0.64, 0.65] $\mathrm{GW}$, and 0.87 $\mathrm{GW}$ in period 2. As the $\alpha$ value increases, the capacity expansion of biomass power and garbage power shows a growth tendency, while the $\beta$ value only weakly influences the schemes of each scenario. For example, in period 2, under the medium demand level, biomass power capacity expansion is [0.44, 0.45] $\mathrm{GW}$, [0.48, 0.49] $\mathrm{GW}$, and [0.52, 0.54] $\mathrm{GW}$ in scenarios 2, 5, and 8, respectively; while under the low demand level, garbage power capacity expansion is [0.25, 0.27] $\mathrm{GW}$, [0.30, 0.31] $\mathrm{GW}$, and [0.32, 0.34] $\mathrm{GW}$ in scenarios 3, 6, and 9, respectively. In addition, with the large-scale construction of wind power plants in recent years, there would be nearly no expansion in the future, only 0.007 $\mathrm{GW}$ in period 2 under each scenario. On the whole, renewable energy would be a new frontier of power industry development, becoming the first choice of power capacity expansion, especially in terms of biomass power and garbage power. As the most famous grain production base in the world, Harbin City has rich biomass resources; thus, biomass utilization is one of the main aspects of exploiting renewable energy. In addition, with the expansion of the city, municipal waste production is continuously growing; therefore, garbage power is an effective measure to reduce rubbish, particularly under the harmless and resource recycling targets. At present, the rapid progress of biomass power and garbage power is also supported by efficient technologies and favorable policy guidance. If the bottlenecks of the other renewable energy (e.g., hydro and solar power) are technically broken in the future, the ratio of thermal power will decline dramatically, and the energy structure in Harbin must then undergo further adjustment.

4.3. Total System Cost

Figure 6 presents the total system cost under the different scenarios. Obviously, the total system cost grows in line with the increasing $\alpha$ and $\beta$ values. The cost reaches RMB ¥ [115.11, 127.87] × 10⁹ in the condition of the original scenario. However, the cost of system is RMB ¥ [115.72, 128.45] × 10⁹, RMB ¥ [116.23, 129.05] × 10⁹, and RMB ¥ [116.77, 129.68] × 10⁹ under scenarios 2, 5, and 8, corresponding to $\alpha$ with the value of 1%, 2%, and 3%, respectively; the cost of the system is RMB ¥ [116.23, 129.05] × 10⁹, RMB ¥ [116.36, 129.05] × 10⁹, and RMB ¥ [116.52, 129.05] × 10⁹ under scenarios 5, 6, and 7, corresponding to $\beta$ with the value of 2%, 4%, and 6%, respectively. Therefore, a conclusion can be drawn that the total system cost under different scenarios would be determined through combining power structure adjustment and pollutants emission reduction. For the sake of meeting stricter environmental requirements and up-scaling the proportion of renewable power generation, a large number of environmental protection facilities and renewable power projects would need to be expanded, which would lead to a high capital cost.

4.4. Post-Optimization Analysis by TOPSIS

For the sake of further exploring the mechanism that influences the selection of the “optimal” scheme of the electricity system in the different scenarios, this study applied TOPSIS to provide a systemic framework for ranking alternatives by comparing their distance to the ideal solution. This is a useful approach that is able to identify and reflect the critical factors of decision-making in electric power system management during the post-optimization period. From the viewpoints of electricity supply-demand balance, environmental pressures, and power structure adjustment in the region, five characteristic variables are specified as attribute indicators. According to expert consultation and questionnaires on the importance of evaluation, as shown in

Table 7. The weights of attribute indicators for TOPSIS.

Attribute Indicators	Total System Cost	Emission Reduction of Air Pollutants	Coal Consumption	Renewable Power Generation	Electricity Purchased
Weights	0.35	0.2	0.2	0.15	0.1

, the weights of the different attributes ($w_{}$ = 0.35, 0.2, 0.2, 0.15, and 0.1) were selected. Then, the values of the relevant decision attributes and the decision matrix under the ten scenarios were obtained. As all of the parameters in the study are expressed in the form of interval values, there exists the lower bound and upper bound matrices.

According to the preference selection of Harbin electricity power system management, positive/negative ideal solutions can be identified through comparing the associated attribute vectors. The positive/negative ideal solution in the lower bound is (0.10952, 0.06386, 0.06271, 0.05110, 0.02312) and (0.11212, 0.06263, 0.06405, 0.04458, 0.04127), respectively; and the corresponding upper bound is (0.10976, 0.06344, 0.06054, 0.04901, 0.02249) and (0.11131, 0.06275, 0.06418, 0.04429, 0.03958), respectively. After calculating the deviation of each scenario from its positive and negative ideal solutions, the relative closeness degree can be estimated for illustrating the distance of the attributes from the best ideal solution (as shown in

Table 8. The relative degree of closeness of the different scenarios.

Scenario	Parameter Setting		Relative Closeness
	$\mathit{\alpha }$	$\mathit{\beta }$	${\mathit{C}}_{\mathit{i}}^{*}$
S1	0%	0%	[0.7305,0.7407]
S2	1%	2%	[0.7367,0.7432]
S3	1%	4%	[0.7292,0.7321]
S4	1%	6%	[0.7212,0.7245]
S5	2%	2%	[0.4272,0.4517]
S6	2%	4%	[0.5149,0.6342]
S7	2%	6%	[0.4635,0.5121]
S8	3%	2%	[0.2226,0.2788]
S9	3%	4%	[0.2201,0.2208]
S10	3%	6%	[0.2590,0.2660]

). From the results, it can be observed that the values of the relative degree of closeness decline as the $\alpha$ and $\beta$ values increase, and the variation of $\alpha$ has a greater effect on the calculation. For example, the relative degree of closeness is (0.7367, 0.7432), (0.4272, 0.4517), and (0.2226, 0.2788) under scenarios 2, 5, and 8, corresponding to fixed $\beta$ as 2% and variable $\alpha$ with the value of 1%, 2%, and 3%, respectively. Furthermore, the relative degree of closeness is (0.7367, 0.7432), (0.7292, 0.7321), and (0.7212, 0.7245) under scenarios 2, 3, and 4, corresponding to fixed $\alpha$ as 1% and variable $\beta$ with the value of 2%, 4%, and 6%, respectively. This demonstrates that electric power system structure adjustment would have a more profound influence on optimal scheme selection than controlling the total amount of pollutants emission in Harbin. Furthermore, scenario 2 has the maximum value in terms of the relative degree of closeness, which is considered the optimal scheme among the ten scenarios. This illustrates that cutting 1% of the proportion of thermal power and reducing 2% of the proportion of the total pollutants emission target will bring the most optimized effect for Harbin electric power system management planning. This discovery also demonstrates that depending on the compression of the thermal power ratio and controlling the pollutants emission on a large scale cannot produce a desirable effect on the current energy environment system. Therefore, under the precondition of regional power system stability, relevant government sectors should focus on the task of fine-tuning the power structure in Harbin City.

5. Discussion and Conclusions

In this study, a comprehensive management framework for environmental-oriented power system management was proposed. The framework integrates electric power consumption forecasting, the regional electric system planning model, and post-optimization analysis. Through the prediction of total electricity consumption by utilizing the Gray-Markov forecasting approach, an inexact two-stage stochastic robust programming model was developed and ten scenarios of different thermal power proportion levels and pollutants mitigation levels were designed. The proposed model is efficient in addressing various types of uncertainties expressed as stochastic probabilities and interval values, and it is expert at the quantitative assessment of system stabilization. Furthermore, the developed model is capable of striking an equilibrium between pre-regulated electric power generation targets and realistic generation schemes. This study carried out a detailed analysis of the optimized electricity power generation schemes, the capacity expansion schemes, and the total system cost. In addition, in order to identify the most influential impact factors of optimal scheme selection, the TOPSIS method was applied during the post-optimization period.

However, as this study mainly focused on macroscope planning of the electricity power industry, the impacts of the power transmission system [45], the distributed energy [46], the scheduling process [47], and the emission taxes [48] were not involved. These factors will be prioritized in future study. Furthermore, in the future, improvements can be conducted to incorporate more system components, more advanced techniques, as well as more uncertain information into the modeling framework.

To sum up, the conclusions drawn from this study are valuable for (a) forecasting electricity consumption based on limited historical data characterized by stochasticity, and the forecasted results are presented in the form of interval parameters with corresponding random probability; (b) making stable and reliable power generation schemes in terms of different power generation modes based on power structure adjustment plans and regional pollutants emission reduction policies; and (c) identifying the most influential factors of planning decisions by comprehensively analyzing and generating optimum electricity generation schemes. In detail, the results indicate that a diversified and distinctive power structure (i.e., thermal power as the pillar industry and renewable energy as supplementary) should be formed to maintain a sustainable development and to steadily supplement electricity in the study area. Biomass power and garbage power both have increasingly important roles in the power system as the thermal power generation ratio decreases. Analysis of various scenarios indicated that fine-tuning the power structure would bring the most optimized effect in regional electric power system management planning. According to the results of this study, decision-makers could establish a satisfactory electricity system management scheme under multiple uncertainties in real life.