An inexact two-stage dynamic stochastic model for regional electricity and heat supply management with pollutants mitigation control

Abstract

Background

Energy system management is an important tool for regional energy and environmental development, and many parameters and their interrelationships in energy-environmental management model appear complexity and uncertain. How to deal with these uncertainties and make a reasonable decision schemes are desired for managers.

Results

In this study, an inexact two-stage dynamic programming model is developed for regional electricity and heat supply management under considering the complexities and uncertainties in regional energy system. The model can reflect not only uncertainties expressed as probability distribution but also those being available as intervals. The developed model is applied to a case of planning regional electricity and heat supply as well as pollution emission reduction considered.

Conclusions

A number of scenarios corresponding to different pollutants emission reduction levels are examined; the results indicated that reasonable solutions have been generated under different pollutants reduction levels. They can be used for generating plans for energy resource/electricity/heat allocation and capacity expansion and help decision makers identify desired regional electricity and heat supply which need minimum cost under various standards of pollutants emission reduction control.

Background

Currently, with the speedup of urbanization and the further economic development, the limitation of electricity resources to regional development becomes more and more obvious. On the one hand, globally, energy consumption grows more rapidly than the economy, meaning the energy intensity of economic activity has rosed. On the other hand, urban air quality is one of the hot topics, and some air pollutant emissions are from burning fossil fuels such as coal, natural gas which mainly generated in electricity and heat supply, especially, during long winter heating period which is the most serious air quality damage season. The goal of energy conservation and pollutant emission reduction in the regional electricity and heat supply is in accordance with the stable and orderly development of the electricity and heat supply enterprises. Therefore, effective electricity and heat supply system planning method with pollutants emission reduction is desired urgently.

Previously, some deterministic models for energy management were developed. For example, Kwaczek et al. ([1996]) put forward an optimization model for comprehending economic impacts of various emission-reduction strategies on energy activities in Saskatchewan. Sailor ([1997]) conducted a comprehensive assessment of climate change influence on renewable energy demand and supply technologies in many places. Zhang et al. ([2001]) reviewed the relationship between global warming and energy system structural shift in power generation sector in south China. Heinrich et al. ([2007]) presented the South African ESI using a partial equilibrium E3 model approach, and extend the approach to include multiple objectives under selected future uncertainties. Klaassen and Riahi ([2007]) utilized the long-term MESSAGE (Model for Energy Supply Strategy Alternatives and their General Environmental Impact) to analyze energy planning management and climate change response. Chung et al. ([2009]) conducted a hybrid E-IO (Energy top-down approach) table which has higher classification sector resolutions to determine the strength of optimization model for making clear of economic impacts from various emission-reduction strategies on energy activities in Saskatchewan, Canada.

However, the above models emphasized on the planning of either electricity supply or heat supply, however, could hardly achieve energy efficiency. Though, some studies about combined heat and electricity generation have been done (Motevasel et al. [2011]; Dong et al. [2012]; Mehdi et al. [2013]). However, there are defects: the uncertainties are not addressed in the previous studies which are critical for analyzing system reliability; on the other hand, some studies only pursue the maximum benefits in the operation of enterprises ignoring the strategic thinking of pollutants emission reduction control (Gustavsson and Madlener [2003]; Cai et al. [2009]).

Thus, in this study, an inexact two-stage dynamic programming (ITSDP) model is developed for regional electricity and heat supply management under uncertainties, aimed at the co-win between the enterprises profit and social benefits. Furthermore, this model performs satisfactory role on both cost reduction and pollutants emission reduction, with the ability to supply necessary information for decision makers, which could be applied to similar planning.

Methods

A general inexact two-stage stochastic programming can be formulated as follows:

$$\mathrm{Min}{\mathit{f}}^{\pm }={\displaystyle \sum _{\mathit{j}=1}^{{\mathit{n}}_1}{\mathit{c}}_{\mathit{j}}^{\pm }{\mathit{x}}_{\mathit{j}}^{\pm }}+{\displaystyle \sum _{\mathit{j}=1}^{{\displaystyle {\mathit{n}}_2}}{\displaystyle \sum _{\mathit{h}=1}^{\mathit{v}}{\mathit{p}}_{\mathit{jh}}{\mathit{d}}_{\mathit{j}}^{\pm }{\mathit{y}}_{\mathit{jh}}^{\pm }}}$$

(1)

subject to

$${\mathit{a}}_{\mathit{r}}^{\pm }{\mathit{x}}^{\pm }\le {\mathit{b}}_{\mathit{r}}^{\pm },\mathit{r}=1,2,\dots ,{\mathit{m}}_1$$

(2)

$${\mathit{a}}_{\mathit{t}}^{\pm }{\mathit{x}}^{\pm }+{\mathit{a}}_{\mathit{t}}^{\text{'}\pm }{\mathit{y}}^{\pm }\ge {\mathit{\omega }}_{\mathit{h}}^{\pm },\mathit{t}=1,2,\dots ,{\mathit{m}}_2;\mathit{h}=1,2,\dots ,\mathit{v}$$

(3)

$${\mathit{x}}_{\mathit{j}}^{\pm }\ge 0,\mathit{j}=1,2,\dots ,{\mathit{n}}_1$$

(4)

$${\mathit{y}}_{\mathit{jh}}^{\pm }\ge 0,\mathit{j}=1,2,\dots ,{\mathit{n}}_2;\mathit{h}=1,2,\dots ,\mathit{v}.$$

(5)

where ${\mathit{a}}_{\mathit{r}}^{\pm }\in {\left\{{\mathit{R}}^{\pm }\right\}}^{{\mathit{m}}_1\times {\mathit{n}}_1}$, ${\mathit{a}}_{\mathit{t}}^{\pm }\in {\left\{{\mathit{R}}^{\pm }\right\}}^{{\mathit{m}}_2\times {\mathit{n}}_2}$, ${\mathit{b}}_{\mathit{r}}^{\pm }\in {\left\{{\mathit{R}}^{\pm }\right\}}^{{\mathit{m}}_1\times 1}$, ${\mathit{c}}_{\mathit{j}}^{\pm }\in {\left\{{\mathit{R}}^{\pm }\right\}}^{1\times {\mathit{n}}_1}$, ${\mathit{d}}_{\mathit{j}}^{\pm }\in {\left\{{\mathit{R}}^{\pm }\right\}}^{1\times {\mathit{n}}_2}$, ${\mathit{x}}^{\pm }\in {\left\{{\mathit{R}}^{\pm }\right\}}^{{\mathit{n}}_1\times 1}$, ${\mathit{y}}^{\pm }\in {\left\{{\mathit{R}}^{\pm }\right\}}^{{\mathit{n}}_2\times 1}$ and {R ^±} denote a set of interval parameters and/or variables. According to Huang et al. [2001], model (1) can be transformed into two deterministic submodels that correspond to the lower and upper bounds of the desired objective function. The objective function value corresponding to f ⁻ is desired first because the objective is to minimize net system costs, and it can be formulated as follows (assume that ${\mathit{b}}_{\mathit{r}}^{\pm }\ge 0$):

$$\mathrm{Min}{\mathit{f}}^{-}={\displaystyle \sum _{\mathit{j}=1}^{{\mathit{k}}_1}{\mathit{c}}_{\mathit{j}}^{-}{\mathit{x}}_{\mathit{j}}^{-}}+{\displaystyle \sum _{\mathit{j}={\mathit{k}}_1+1}^{{\displaystyle {\mathit{n}}_1}}{\mathit{c}}_{\mathit{j}}^{-}{\mathit{x}}_{\mathit{j}}^{+}}+{\displaystyle \sum _{\mathit{j}=1}^{{\displaystyle {\mathit{k}}_2}}{\displaystyle \sum _{\mathit{h}=1}^{\mathit{v}}{\mathit{p}}_{\mathit{jh}}{\mathit{d}}_{\mathit{j}}^{-}{\mathit{y}}_{\mathit{jh}}^{-}}}+{\displaystyle \sum _{\mathit{j}={\mathit{k}}_2+1}^{{\mathit{n}}_2}{\displaystyle \sum _{\mathit{h}=1}^{\mathit{v}}{\mathit{p}}_{\mathit{jh}}{\mathit{d}}_{\mathit{j}}^{-}{\mathit{y}}_{\mathit{jh}}^{+}}}$$

(6)

subject to:

$${\displaystyle \sum _{\mathit{j}=1}^{{\mathit{k}}_1}{\left|{\mathit{a}}_{\mathit{rj}}\right|}^{+}\mathrm{sign}\left({\mathit{a}}_{\mathit{rj}}^{+}\right){\mathit{x}}_{\mathit{j}}^{-}}+{\displaystyle \sum _{\mathit{j}={\mathit{k}}_1+1}^{{\displaystyle {\mathit{n}}_1}}{\left|{\mathit{a}}_{\mathit{rj}}\right|}^{-}\mathrm{sign}\left({\mathit{a}}_{\mathit{rj}}^{-}\right){\mathit{x}}_{\mathit{j}}^{+}}\le {\mathit{b}}_{\mathit{r}}^{+},\forall \mathit{r}$$

(7)

$$\begin{array}{l}{\displaystyle \sum _{\mathit{j}=1}^{{\mathit{k}}_1}{\left|{\mathit{a}}_{\mathit{tj}}\right|}^{+}\mathrm{sign}\left({\mathit{a}}_{\mathit{tj}}^{+}\right){\mathit{x}}_{\mathit{j}}^{-}}+{\displaystyle \sum _{\mathit{j}={\mathit{k}}_1+1}^{{\displaystyle {\mathit{n}}_1}}{\left|{\mathit{a}}_{\mathit{tj}}\right|}^{-}\mathrm{sign}\left({\mathit{a}}_{\mathit{tj}}^{-}\right){\mathit{x}}_{\mathit{j}}^{+}}+{\displaystyle \sum _{\mathit{j}=1}^{{\mathit{k}}_2}{\left|{\mathit{a}}_{\mathit{tj}}^{\text{'}}\right|}^{+}\mathrm{sign}\left({\mathit{a}}_{\mathit{tj}}^{\text{'}+}\right){\mathit{y}}_{\mathit{jh}}^{-}}\\ +{\displaystyle \sum _{\mathit{j}={\mathit{k}}_2+1}^{{\displaystyle \mathit{n}2}}{\left|{\mathit{a}}_{\mathit{tj}}^{\text{'}}\right|}^{-}\mathrm{sign}\left({\mathit{a}}_{\mathit{tj}}^{\text{'}-}\right){\mathit{y}}_{\mathit{jh}}^{+}}\ge {\mathit{\omega }}_{\mathit{h}}^{-},\forall \mathit{t},\mathit{h}\end{array}$$

(8)

$${\mathit{x}}_{\mathit{j}}^{-}\ge 0,\mathit{j}=1,2,\dots ,{\mathit{k}}_1$$

(9)

$${\mathit{x}}_{\mathit{j}}^{+}\ge 0,\mathit{j}={\mathit{k}}_1+1,{\mathit{k}}_1+2,\dots ,{\mathit{n}}_1$$

(10)

$${\mathit{y}}_{\mathit{jh}}^{-}\ge 0,\forall \mathit{h};\mathit{j}=1,2,\dots ,{\mathit{k}}_2$$

(11)

$${\mathit{y}}_{\mathit{jh}}^{+}\ge 0,\forall \mathit{h};\mathit{j}={\mathit{k}}_2+1,{\mathit{k}}_2+2,\dots ,{\mathit{n}}_2$$

(12)

where ${\mathit{x}}_{\mathit{j}}^{\pm }$, j = 1, 2,…, k ₁, are positive coefficients; ${\mathit{x}}_{\mathit{j}}^{\pm }$, j = k ₁ + 1, k ₁ + 2,…, n ₁ are negative variables; ${\mathit{y}}_{\mathit{jh}}^{\pm }$, j = 1, 2,…, k ₂ and h = 1, 2,…, v, are random variables with positive coefficients; ${\mathit{y}}_{\mathit{jh}}^{\pm }$,j = k ₂ + 1, k ₂ + 2,…, n ₂ and h = 1, 2,…, v, are random variables with negative coefficients. Solutions of ${\mathit{x}}_{\mathit{j}\mathrm{opt}}^{-}$ (j = 1, 2,…, k ₁), ${\mathit{x}}_{\mathit{j}\mathrm{opt}}^{+}$(j = k ₁ + 1, k ₁ + 2,…, n ₁),${\mathit{y}}_{\mathit{jh}\mathrm{opt}}^{-}$ (j = 1, 2,…, k ₂), and ${\mathit{y}}_{\mathit{jh}\mathrm{opt}}^{+}$ (j = k ₂ + 1, k ₂ + 2,…, n ₂) can be obtained through submodel (2). Based on the above solutions, the submodel f ⁺ can be formulated as follows:

$$\mathrm{Min}{\mathit{f}}^{+}={\displaystyle \sum _{\mathit{j}=1}^{{\mathit{k}}_1}{\mathit{c}}_{\mathit{j}}^{+}{\mathit{x}}_{\mathit{j}}^{+}}+{\displaystyle \sum _{\mathit{j}={\mathit{k}}_1+1}^{{\displaystyle {\mathit{n}}_1}}{\mathit{c}}_{\mathit{j}}^{+}{\mathit{x}}_{\mathit{j}}^{-}}+{\displaystyle \sum _{\mathit{j}=1}^{{\displaystyle {\mathit{k}}_2}}{\displaystyle \sum _{\mathit{h}=1}^{\mathit{v}}{\mathit{p}}_{\mathit{jh}}{\mathit{d}}_{\mathit{j}}^{+}{\mathit{y}}_{\mathit{jh}}^{+}}}+{\displaystyle \sum _{\mathit{j}={\mathit{k}}_2+1}^{{\mathit{n}}_2}{\displaystyle \sum _{\mathit{h}=1}^{\mathit{v}}{\mathit{p}}_{\mathit{jh}}{\mathit{d}}_{\mathit{j}}^{+}{\mathit{y}}_{\mathit{jh}}^{-}}}$$

(13)

subject to:

$${\displaystyle \sum _{\mathit{j}=1}^{{\mathit{k}}_1}{\left|{\mathit{a}}_{\mathit{rj}}\right|}^{-}\mathrm{sign}\left({\mathit{a}}_{\mathit{rj}}^{-}\right){\mathit{x}}_{\mathit{j}}^{+}}+{\displaystyle \sum _{\mathit{j}={\mathit{k}}_1+1}^{{\displaystyle {\mathit{n}}_1}}{\left|{\mathit{a}}_{\mathit{rj}}\right|}^{+}\mathrm{sign}\left({\mathit{a}}_{\mathit{rj}}^{+}\right){\mathit{x}}_{\mathit{j}}^{-}}\le {\mathit{b}}_{\mathit{r}}^{-},\forall \mathit{r}$$

(14)

$$\begin{array}{l}{\displaystyle \sum _{\mathit{j}=1}^{{\mathit{k}}_1}{\left|{\mathit{a}}_{\mathit{tj}}\right|}^{-}\mathrm{sign}\left({\mathit{a}}_{\mathit{tj}}^{-}\right){\mathit{x}}_{\mathit{j}}^{+}}+{\displaystyle \sum _{\mathit{j}={\mathit{k}}_1+1}^{{\displaystyle {\mathit{n}}_1}}{\left|{\mathit{a}}_{\mathit{tj}}\right|}^{+}\mathrm{sign}\left({\mathit{a}}_{\mathit{tj}}^{+}\right){\mathit{x}}_{\mathit{j}}^{-}}+{\displaystyle \sum _{\mathit{j}=1}^{{\mathit{k}}_2}{\left|{\mathit{a}}_{\mathit{tj}}^{\text{'}}\right|}^{-}\mathrm{sign}\left({\mathit{a}}_{\mathit{tj}}^{\text{'}-}\right){\mathit{y}}_{\mathit{jh}}^{+}}\\ +{\displaystyle \sum _{\mathit{j}={\mathit{k}}_2+1}^{{\displaystyle {\mathit{n}}_2}}{\left|{\mathit{a}}_{\mathit{tj}}^{\text{'}}\right|}^{+}\mathrm{sign}\left({\mathit{a}}_{\mathit{tj}}^{\text{'}+}\right){\mathit{y}}_{\mathit{jh}}^{-}}\ge {\mathit{\omega }}_{\mathit{h}}^{+},\forall \mathit{t},\mathit{h}\end{array}$$

(15)

$${\mathit{x}}_{\mathit{j}}^{+}\ge {\mathit{x}}_{\mathit{j}\mathrm{opt}}^{-},\mathit{j}=1,2,\dots ,{\mathit{k}}_1$$

(16)

$$0\le {\mathit{x}}_{\mathit{j}}^{-}\le {\mathit{x}}_{\mathit{j}\mathrm{opt}}^{+},\mathit{j}={\mathit{k}}_1+1,{\mathit{k}}_1+2,\dots ,{\mathit{n}}_1$$

(17)

$${\mathit{y}}_{\mathit{jh}}^{+}\ge {\mathit{y}}_{\mathit{jh}\mathrm{opt}}^{-},\forall \mathit{h};\mathit{j}=1,2,\dots ,{\mathit{k}}_2$$

(18)

$$0\le {\mathit{y}}_{\mathit{jh}}^{-}\le {\mathit{y}}_{\mathit{jh}\mathrm{opt}}^{+},\forall \mathit{h};\mathit{j}={\mathit{k}}_2+1,{\mathit{k}}_2+2,\dots ,{\mathit{n}}_2$$

(19)

Solutions of ${\mathit{x}}_{\mathit{j}\mathrm{opt}}^{+}$ (j = 1, 2,…, k ₁), ${\mathit{x}}_{\mathit{j}\mathrm{opt}}^{-}$ (j = k ₁ + 1, k ₁ + 2,…, n ₁), ${\mathit{y}}_{\mathit{jh}\mathrm{opt}}^{+}$ (j = 1, 2,…, k ₂), and ${\mathit{y}}_{\mathit{jh}\mathrm{opt}}^{-}$ (j = k ₂ + 1, k ₂ + 2,…, n ₂) can be obtained through submodel (3). Through integrating solutions of submodels (2) and (3), the solution for model (1) can be obtained.

Model formulation

Typically, electricity and heat supply system often contains some components such as energy supply/demand, processing and transformation technologies, and electricity and heat generation (Heinrich et al. [2007]). These components generally involve an array of economic activities, energy consumption and pollutants discharge. Electricity and heat supply options are usually classified as fossil energy and renewable resources. Each of those has its own industry representing the characteristics of the technologies.

Consider a case wherein a regional electricity and heat supply manager is responsible to allocate electricity and heat flows from the enterprise to users all the year (12 months). The manager can formulate the problem as minimizing the expected value of net system cost in the region during one year. Based on the local electricity and heat supply policies, a promised allowable supply quantity is defined. If this level for each user is not reached, it will result in the higher cost to the system, the system will then be subject to penalties of the system failure. In the mean time, the manager always seek a project which can assure the emissions of some environmental pollutants (e.g., nitrogen oxides (NO_x), sulfur dioxide (SO₂), particulate matter (PM)), and greenhouse gas (GHG) could meet the regional environmental standard with sulfur dioxide emission fees deduced. As waste-generation and energy requirement amounts from the region are uncertain at the time when the planning decisions must be made, the problem under consideration can be formulated as a TSDP model as follows:

$$min{\mathit{f}}^{\pm }=\left(1\right)+\left(2\right)+\left(3\right)+\left(4\right)+\left(5\right)+\left(6\right)+\left(7\right)+\left(8\right)+\left(9\right)$$

(20)

$$\left(1\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{i}=1}^{\mathit{I}}\mathit{E}{\mathit{S}}_{\mathit{it}}^{\pm }}}\ast \mathit{E}{\mathit{P}}_{\mathit{it}}^{\pm }$$

(21)

$$\left(2\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=1}^2\mathit{T}{\mathit{C}}_{\mathit{kt}}^{\pm }\ast {\mathit{W}}_{\mathit{kt}}^{\pm }+{\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=1}^3{\displaystyle \sum _{\mathit{h}=1}^{\mathit{H}}{\mathit{P}}_{\mathit{ht}}\ast \left(\mathit{T}{\mathit{C}}_{\mathit{kt}}^{\pm }+\mathit{P}{\mathit{P}}_{\mathit{kt}}^{\pm }\right)\ast \mathit{E}{\mathit{Q}}_{\mathit{kth}}^{\pm }}}}}}+\mathit{T}{\mathit{C}}_{3\mathit{t}}^{\pm }\ast \left({\mathit{W}}_{3\mathit{t}}^{\pm }+\mathit{SA}{\mathit{E}}_{\mathit{t}}^{\pm }\right)$$

(22)

$$\left(3\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}\mathit{TP}{\mathit{H}}_{\mathit{jt}}^{\pm }\ast {\mathit{Y}}_{\mathit{jt}}^{\pm }+{\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}{\displaystyle \sum _{\mathit{L}=1}^{\mathit{l}}\mathit{PMlt}\ast \left(\mathit{TP}{\mathit{H}}_{\mathit{jt}}^{\pm }+\mathit{PC}{\mathit{H}}_{\mathit{jt}}^{\pm }\right)\ast \mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{\pm }}}}}}$$

(23)

$$\left(4\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=\mathit{i}}^{\mathit{K}}{\displaystyle \sum _{\mathit{h}=1}^{\mathit{H}}{\displaystyle \sum _{\mathit{r}=1}^{\mathit{R}}({\mathit{W}}_{\mathit{kt}}^{\pm }}}}}+{\mathit{P}}_{\mathit{ht}}\ast \mathit{E}{\mathit{Q}}_{\mathit{kth}}^{\pm })\ast \mathit{P}{\mathit{E}}_{\mathit{rkt}}^{\pm }\ast \mathit{C}{\mathit{E}}_{\mathit{rkt}}^{\pm }$$

(24)

$$\left(5\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}{\displaystyle \sum _{\mathit{l}=1}^{\mathit{L}}{\displaystyle \sum _{\mathit{q}=1}^{\mathit{Q}}\left({\mathit{Y}}_{\mathit{jt}}^{\pm }+\mathit{P}{\mathit{M}}_{\mathit{lt}}\ast \mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{\pm }\right)\ast \mathit{P}{\mathit{H}}_{\mathit{qjt}}^{\pm }\ast \mathit{C}{\mathit{H}}_{\mathit{qjt}}^{\pm }}}}}$$

(25)

$$\left(6\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=1}^{\mathit{K}}{\displaystyle \sum _{\mathit{m}=1}^{\mathit{M}}{\displaystyle \sum _{\mathit{h}=1}^{\mathit{H}}{\mathit{P}}_{\mathit{ht}}}}}}\ast \mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{\pm }\ast \mathit{DI}{\mathit{V}}_{\mathit{kmt}}\ast \mathit{DC}{\mathit{O}}_{\mathit{kmt}}$$

(26)

$$\left(7\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}{\displaystyle \sum _{\mathit{n}=1}^{\mathit{N}}{\displaystyle \sum _{\mathit{l}=1}^{\mathit{L}}\mathit{P}{\mathit{M}}_{\mathit{lt}}\ast \mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{\pm }}}}}\ast \mathit{DV}{\mathit{H}}_{\mathit{jnt}}\ast \mathit{DC}{\mathit{H}}_{\mathit{jnt}}$$

(27)

$$\left(8\right)={\displaystyle \sum _{\mathit{t}=1}^{12}\mathit{SA}{\mathit{E}}_{\mathit{t}}^{\pm }\cdot \mathit{CE}{\mathit{S}}_{\mathit{t}}^{\pm }}$$

(28)

$$\left(9\right)={\displaystyle \sum _{\mathrm{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=1}^{\mathit{K}}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}\left[\mathit{P}{\mathit{E}}_{1\mathit{kt}}^{\pm }\ast \left(1-{\mathit{\eta }}_{1\mathit{kt}}^{\pm }\right)+\mathit{P}{\mathit{H}}_{1\mathit{jt}}^{\pm }\ast \left(1-{\mathit{\mu }}_{1\mathit{jt}}^{\pm }\right)\right]\ast \mathit{PD}{\mathit{F}}_{\mathit{t}}^{\pm }}}}$$

(29)

where: i denotes the energy sources, i = 1 for coal, i = 2 for natural gas; k denotes the power generation technologies, k = 1 for coal-fired power, k = 2 for natural gas-fired power, k = 3 for wind power; t denotes the planning periods; r is the air pollutants, r = 1 for sulfur dioxide, r = 2 for nitrogen oxides, r = 3 for particulate matter; h is the demand level. j denotes the heat supply technologies, j = 1 for coal-fired heat, j = 2 for natural gas-fired heat, j = 3 for electric heat; l is capacity expansion size option for heat generation; f ^± = the net expected system total cost (million dollar); $\mathit{E}{\mathit{P}}_{\mathit{it}}^{\pm }$ = the supply of energy resources in month t (PJ); $\mathit{E}{\mathit{S}}_{\mathit{it}}^{\pm }$ = the supply cost of energy sources i in month t($million/PJ); $\mathit{T}{\mathit{C}}_{\mathit{kt}}^{\pm }$ = the variable cost for electricity generated by technology k in month t ($million/GWh); ${\mathit{W}}_{\mathit{kt}}^{\pm }$ = allowable power generation by technology k during month t (GWh); P _ht = probability of occurrence for scenario h during month t; PP _kt = penalty cost of excess electricity generated by technology k in month t ($million/GWh); $\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{\pm }$ = the additional power generation because of the shortage of electricity generated by technology k in scenario h during month t (GWh); $\mathit{TP}{\mathit{H}}_{\mathit{jt}}^{\pm }$ = the variable cost for heat generated by technology j in month t ($million/GWh); ${\mathit{Y}}_{\mathit{jt}}^{\pm }$ = allowable heat generation by technology j during month t (PJ); PM _lt = probability of occurrence for scenario l during month t; PCH _jt = penalty cost of excess heat generated by technology j in month t ($million/PJ); $\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{\pm }$ = the additional heat generation because of the shortage of heat generated by technology k in scenario l during month t (PJ); $\mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{\pm }$ = binary variable for technology k with expansion option m in scenario h; DIV _kmt = capacity expansion size option m for power generation technology k in month t (GW); DCO _kmt = capacity cost of capacity expansion size m for power generation technology k in month t ($Million/GW); JUH _jml = binary variable for technology j with expansion option n in scenario l during month t; DVH _jnt = capacity expansion size option n for power generation technology j in month t (PJ); DCH _jnt = capacity cost of capacity expansion size n for heat generation technology j in month t ($Million/PJ); $\mathit{P}{\mathit{E}}_{\mathit{krt}}^{\pm }$ = the emission intensity of pollutant r from power generation technology k in month t (kiloton/GWh); $\mathit{C}{\mathit{E}}_{\mathit{krt}}^{\pm }$ = the removal cost of pollutant r from power generation technology k in month t (dollar/kiloton); $\mathit{P}{\mathit{H}}_{\mathit{jqt}}^{\pm }$ = the emission intensity of pollutant q from power generation technology j in month t (kiloton/PJ);

$\mathit{C}{\mathit{H}}_{\mathit{jrt}}^{\pm }$ = the removal cost of pollutant r from power generation technology j in month t (dollar/kiloton); $\mathit{CE}{\mathit{S}}_{\mathit{t}}^{\pm }$ = the cost of per unit wind power storage in month t (dollar/kWh);

$\mathit{SA}{\mathit{E}}_{\mathit{t}}^{\pm }$ = the stored wind power in month t (GW); $\mathit{PD}{\mathit{F}}_{\mathit{t}}^{\pm }$ = pollutants discharge fee for sulfur dioxide (dollar/tonnes); ${\mathit{\eta }}_{\mathit{krt}}^{\pm }$ = the removal efficiency of pollutant r from power generation technology k in month t; ${\mathit{\mu }}_{\mathit{lqt}}^{\pm }$ = the removal efficiency of pollutant q from heat generation technology j in month t.

Mass balance constraints

The mass balance constraints describe the balance of resource and energy flows in the system. They can be classified into three groups:

• (1) balance for energy resource (4k) and (4l);

• (2) balance for electricity generation (4m), (4n), (4o), (4p) and (4q);

• (3) balance for heat generation (4r) and (4s).

These constraints are established to ensure that the input energy is greater than the output one.

$$\left({\mathit{W}}_{\mathit{kt}}^{\pm }+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{\pm }\right)\ast \mathit{EF}{\mathit{E}}_{\mathit{kt}}^{\pm }+\left({\mathit{Y}}_{\mathit{jt}}^{\pm }+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{\pm }\right)\ast \mathit{EF}{\mathit{H}}_{\mathit{jt}}^{\pm }+{\displaystyle \sum _{\mathit{s}=1}^{\mathit{S}}\mathit{E}{\mathit{D}}_{\mathit{its}}^{\pm }}\le {\mathit{T}}_{\mathit{it}}^{\pm },\forall \mathit{t},\mathit{k},\mathit{j},\mathit{l},\mathit{h},\mathit{i}$$

(30)

$$\left({\mathit{W}}_{3\mathit{t}}^{\pm }+\mathit{E}{\mathit{Q}}_{3\mathit{th}}^{\pm }+\mathit{SA}{\mathit{E}}_{\mathit{t}}^{\pm }\right)\ast \mathit{EF}{\mathit{E}}_{3\mathit{t}}^{\pm }\le {\mathit{T}}_{\mathit{it}}^{\pm },\forall \mathit{t},\mathit{k},\mathit{h}$$

(31)

$${\displaystyle \sum _{\mathit{k}=1}^3\left({\mathit{W}}_{\mathit{kt}}^{\pm }+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{\pm }\right)}+\mathit{SA}{\mathit{E}}_{\mathit{t}-1}^{\pm }\ge \mathit{DT}{\mathit{E}}_{\mathit{th}}^{\pm },\forall \mathit{t},\mathit{h}$$

(32)

$${\displaystyle \sum _{\mathit{k}=1}^3{\mathit{W}}_{\mathit{kt}}^{\pm }}-\mathit{DT}{\mathit{E}}_{\mathit{th}}^{\pm }+\mathit{SA}{\mathit{E}}_{\mathit{t}-1}^{\pm }=\mathit{SA}{\mathit{E}}_{\mathit{t}}^{\pm },\forall \mathit{t},\mathit{h}$$

(33)

$${\mathit{W}}_{\mathit{kt}}^{\pm }+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{\pm }\le \mathit{E}{\mathit{T}}_{\mathit{kt}}^{\pm }\ast \left(\mathit{R}{\mathit{C}}_{\mathit{k}}+\mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{\pm }\ast \mathit{DI}{\mathit{V}}_{\mathit{kmt}}\right),\forall \mathit{t},\mathit{k},\mathit{h},\mathit{m}$$

(34)

$${\mathit{W}}_{3\mathit{t}}^{\pm }-\left[\mathit{DT}{\mathit{E}}_{\mathit{th}}^{\pm }-\left({\mathit{W}}_{1\mathit{t}}^{\pm }+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{\pm }\right)-\left({\mathit{W}}_{2\mathit{t}}^{\pm }+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{\pm }\right)\right]+\mathit{SA}{\mathit{E}}_{\mathit{t}-1}^{\pm }=\mathit{SA}{\mathit{E}}_{\mathit{t}}^{\pm },\forall \mathit{t}\ge 2,\mathit{k},\mathit{h}$$

(35)

$${\mathit{W}}_{3\mathit{t}}^{\pm }+\mathit{E}{\mathit{Q}}_{3\mathit{th}}^{\pm }-\left[\mathit{DT}{\mathit{E}}_{\mathit{th}}^{\pm }-\left({\mathit{W}}_{1\mathit{t}}^{\pm }+\mathit{E}{\mathit{Q}}_{1\mathit{th}}^{\pm }\right)-\left({\mathit{W}}_{2\mathit{t}}^{\pm }+\mathit{E}{\mathit{Q}}_{2\mathit{th}}^{\pm }\right)\right]+\mathit{SAE}{\mathrm{O}}^{\pm }=\mathit{SA}{\mathit{E}}_{\mathit{t}}^{\pm },\mathit{t}=1,\forall \mathit{h}$$

(36)

$${\displaystyle \sum _{\mathit{j}=1}^3\left({\mathit{Y}}_{\mathit{jt}}^{\pm }+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{\pm }\right)}\ge \mathit{TD}{\mathit{H}}_{\mathit{tl}}^{\pm },\forall \mathit{t},\mathit{l}$$

(37)

$${\mathit{Y}}_{\mathit{jt}}^{\pm }+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{\pm }\le \mathit{H}{\mathit{T}}_{\mathit{jt}}^{\pm }\ast \left(\mathit{RC}{\mathit{H}}_{\mathit{j}}+\mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{\pm }\ast \mathit{DV}{\mathit{H}}_{\mathit{jnt}}\right),\forall \mathit{t},\mathit{j},\mathit{n},\mathit{l}$$

(38)

Capacity constraints of technologies

For a single technology, it is assumed that its output or production should be less than the amount that total installed capacity. If this requirement is not satisfied, cost will be increased for additional capacities.

$$\mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{\pm }=\left\{\begin{array}{l}1\\ 0\end{array}\right.,\forall \mathit{k},\mathit{t},\mathit{m},\mathit{h}$$

(39)

$$\mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{\pm }=\left\{\begin{array}{l}1\\ 0\end{array}\right.,\forall \mathit{j},\mathit{t},\mathit{n},\mathit{l}$$

(40)

$${\displaystyle \sum _{\mathit{m}=1}^{\mathit{M}}\mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{\pm }}\le 1,\forall \mathit{k},\mathit{t},\mathit{h}$$

(41)

$${\displaystyle \sum _{\mathrm{n}=1}^{\mathit{N}}\mathit{JU}{\mathrm{H}}_{\mathit{jtnl}}^{\pm }\le 1},\forall \mathit{j},\mathit{t},\mathit{l}$$

(42)

Environmental constraints

For an electricity and heat supply system planning, it is assumed that environmental requirement should be considered as an important constraint. Equation (4x) is the constraint of pollutants emission of power generation technologies; Equation (4y) is the constraint of pollutants emission of heat generation technologies. Equation (4z) is the constraint of GHG (Greenhouse Gas) emission of the electricity and heat supply system.

$${\displaystyle \sum _{\mathit{k}=1}^3\left({\mathit{W}}_{\mathit{kt}}^{\pm }+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{\pm }\right)}\ast \mathit{P}{\mathit{E}}_{\mathit{krt}}^{\pm }\ast \left(1-{\mathit{\eta }}_{\mathit{krt}}^{\pm }\right)\le \mathit{T}{\mathit{P}}_{\mathit{rt}}^{\pm },\forall \mathit{t},\mathit{h},\mathit{r}$$

(43)

$${\displaystyle \sum _{\mathit{j}=1}^3\left({\mathit{Y}}_{\mathit{jt}}^{\pm }+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{\pm }\right)}\ast \mathit{P}{\mathit{H}}_{\mathit{jqt}}^{\pm }\ast \left(1-{\mathit{\mu }}_{\mathit{jqt}}^{\pm }\right)\le \mathit{TH}{\mathit{P}}_{\mathit{qt}}^{\pm },\forall \mathit{t},\mathit{l},\mathit{q}$$

(44)

$${\displaystyle \sum _{\mathit{s}=1}^{\mathit{S}}{\displaystyle \sum _{\mathit{i}=1}^2\mathit{E}{\mathit{D}}_{\mathit{its}}^{\pm }}}\ast \mathit{IN}{\mathit{T}}_{\mathit{s}}^{\pm }+{\displaystyle \sum _{\mathit{k}=1}^3({\mathit{W}}_{\mathit{kt}}^{\pm }}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{\pm })\ast \mathit{CO}{\mathit{E}}_{\mathit{kt}}^{\pm }+{\displaystyle \sum _{\mathit{j}=1}^3({\mathit{Y}}_{\mathit{jt}}^{\pm }+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{\pm }})\ast \mathit{CO}{\mathit{H}}_{\mathit{jt}}^{\pm }\le \mathit{T}{\mathit{C}}_{\mathit{t}}^{\pm },\forall \mathit{t},\mathit{h},\mathit{l}$$

(45)

where:

$\mathit{EF}{\mathit{E}}_{\mathit{kt}}^{\pm }$ = the conversion efficiency of power generation technology k in month t (PJ/GW);

$\mathit{EF}{\mathit{H}}_{\mathit{jt}}^{\pm }$ = the conversion efficiency of heat generation technology j in month t (PJ/PJ);

$\mathit{E}{\mathit{D}}_{\mathit{its}}^{\pm }$ = the demand of energy resource i in sector s during month t;

${\mathit{T}}_{\mathit{it}}^{\pm }$ = the supply of energy i in month t;

$\mathit{TD}{\mathit{E}}_{\mathit{th}}^{\pm }$ = electricity demand in scenario h during month t (GWh);

$\mathit{TD}{\mathit{H}}_{\mathit{tl}}^{\pm }$ = heat demand in scenario l during month t (PJ);

$\mathit{E}{\mathit{T}}_{\mathit{kt}}^{\pm }$ = the working hours of power generation technology k in month t (hour);

$\mathit{H}{\mathit{T}}_{\mathit{jt}}^{\pm }$ = the working hours of heat generation technology j in month t (hour);

$\mathit{R}{\mathit{C}}_{\mathit{k}}^{\pm }$ = residual capacity of power generation technology k (GW);

${\mathit{\eta }}_{\mathit{krt}}^{\pm }$ = residual capacity of heat generation technology j (PJ);

$\mathit{T}{\mathit{P}}_{\mathit{rt}}^{\pm }$ = the total allowable emissions of pollutant r in month t (kiloton);

$\mathit{TH}{\mathit{P}}_{\mathit{qt}}^{\pm }$ = the total allowable emissions of pollutant q in month t (kiloton);

$\mathit{IN}{\mathit{T}}_{\mathit{s}}^{\pm }$ = CO₂ emission intensity of sector s (kiloton/PJ);

$\mathit{CO}{\mathit{E}}_{\mathit{kt}}^{\pm }$ = CO₂ emission intensity of power generation technology k in month t (kiloton/GWh);

$\mathit{CO}{\mathit{H}}_{\mathit{jt}}^{\pm }$ = CO₂ emission intensity of power generation technology j in month t (kiloton/PJ);

$\mathit{T}{\mathit{C}}_{\mathit{t}}^{\pm }$ = the total allowable CO₂ emissions in scenario h during month t (kiloton).

In the ITSDP model, when the allowable amount of power generation ${\mathit{W}}_{\mathit{kt}}^{\pm }$ and heat generation ${\mathit{Y}}_{\mathit{jt}}^{\pm }$ are known, the above model can be transformed into two branches of deterministic sub-models, which related to the upper and lower bounds of the conceived objective-function value. The transformation is on the base of an interactive algorithm, which is different from traditional interval analysis as well as best/worst case analysis, and the previous methods for solving inexact linear programming problems cannot be used for granted directly (Li et al. [2006a], [b]; Klaassen and Riahi [2007]). In this study, an optimized set would be made and correspond to minimize the system cost under the uncertain electricity and heat demands and supplies. Accordingly, let ${\mathit{W}}_{\mathit{kt}}^{\pm }={\mathit{W}}_{\mathit{kt}}^{-}+\mathrm{\Delta }{\mathit{W}}_{\mathit{kt}}\cdot {\mathit{u}}_{\mathit{kt}}$, ${\mathit{Y}}_{\mathit{jt}}^{\pm }={\mathit{Y}}_{\mathit{jt}}^{-}+\mathrm{\Delta }{\mathit{Y}}_{\mathit{jt}}\cdot {\mathit{v}}_{\mathit{jt}}$, where $\mathrm{\Delta }{\mathit{W}}_{\mathit{kt}}={\mathit{W}}_{\mathit{kt}}^{+}-{\mathit{W}}_{\mathit{kt}}^{-}$, $\mathrm{\Delta }{\mathit{Y}}_{\mathit{kt}}={\mathit{Y}}_{\mathit{kt}}^{+}-{\mathit{Y}}_{\mathit{kt}}^{-}$ u _kt  ∈ [0, 1], v _jt  ∈ [0, 1], u _kt and v _jt are decision variables that can identify an optimized set of target values ${\mathit{W}}_{\mathit{kt}}^{\pm }$ and ${\mathit{Y}}_{\mathit{jt}}^{\pm }$ so that the related policy analyses could be supported. For example, when ${\mathit{W}}_{\mathit{kt}}^{\pm }$ and ${\mathit{Y}}_{\mathit{jt}}^{\pm }$ reach their upper bounds (i.e., when u _kt = 1 and v _jt = 1), a relatively low cost would be procured when the electricity and heat demands are satisfied; a lower penalty might have to be paid if the promised electricity or heat is delivered. Thus, if ${\mathit{W}}_{\mathit{kt}}^{\pm }$ and ${\mathit{Y}}_{\mathit{jt}}^{\pm }$ approach their lower bounds (i.e., when u _kt = 0 and v _jt = 0), there might exist a higher cost as well as more risk of breaking the promise. Consequently, according to Huang and Loucks ([2000]) and Li et al. ([2009]), two transformed deterministic sub-models based on an interactive algorithm through introducing decision variables u _kt and v _jt could obtain the net system cost. Since the first objective is to minimize this cost, the sub-model relevant to lower-bound objective function value f ⁻ is preferential. Therefore, we have model B:

$$min{\mathit{f}}^{-}=\left(1\right)+\left(2\right)+\left(3\right)+\left(4\right)+\left(5\right)+\left(6\right)+\left(7\right)+\left(8\right)+\left(9\right)$$

(46)

$$\left(1\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{i}=1}^{\mathit{I}}\mathit{E}{\mathit{S}}_{\mathit{it}}^{-}}}\ast \mathit{E}{\mathit{P}}_{\mathit{it}}^{-}$$

(47)

$$\left(2\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=1}^2\mathit{T}{\mathit{C}}_{\mathit{kt}}^{-}\ast {\mathit{W}}_{\mathit{kt}}+{\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=1}^3{\displaystyle \sum _{\mathit{h}=1}^{\mathit{H}}{\mathit{P}}_{\mathit{ht}}\ast \left(\mathit{T}{\mathit{C}}_{\mathit{kt}}^{-}+\mathit{P}{\mathit{P}}_{\mathit{kt}}^{-}\right)\ast \mathit{E}{\mathit{Q}}_{\mathit{kth}}^{-}}}}}}+\mathit{T}{\mathit{C}}_{3\mathit{t}}^{-}\ast \left({\mathit{W}}_{3\mathit{t}}+\mathit{SA}{\mathit{E}}_{\mathit{t}}^{-}\right)$$

(48)

$$\left(3\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}\mathit{TP}{\mathit{H}}_{\mathit{jt}}^{-}\ast {\mathit{Y}}_{\mathit{jt}}+{\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}{\displaystyle \sum _{\mathit{L}=1}^{\mathit{l}}{\mathit{PM}}_{\mathit{lt}}\ast \left(\mathit{TP}{\mathit{H}}_{\mathit{jt}}^{-}+\mathit{PC}{\mathit{H}}_{\mathit{jt}}^{-}\right)\ast \mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{-}}}}}}$$

(49)

$$\left(4\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=\mathit{i}}^{\mathit{K}}{\displaystyle \sum _{\mathit{h}=1}^{\mathit{H}}{\displaystyle \sum _{\mathit{r}=1}^{\mathit{R}}({\mathit{W}}_{\mathit{kt}}}}}}+{\mathit{P}}_{\mathit{ht}}\ast \mathit{E}{\mathit{Q}}_{\mathit{kth}}^{-})\ast \mathit{P}{\mathit{E}}_{\mathit{rkt}}^{-}\ast \mathit{C}{\mathit{E}}_{\mathit{rkt}}^{-}$$

(50)

$$\left(5\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}{\displaystyle \sum _{\mathit{l}=1}^{\mathit{L}}{\displaystyle \sum _{\mathit{q}=1}^{\mathit{Q}}\left({\mathit{Y}}_{\mathit{jt}}+\mathit{P}{\mathit{M}}_{\mathit{lt}}\ast \mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{-}\right)\ast \mathit{P}{\mathit{H}}_{\mathit{qjt}}^{-}\ast \mathit{C}{\mathit{H}}_{\mathit{qjt}}^{-}}}}}$$

(51)

$$\left(6\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=1}^{\mathit{K}}{\displaystyle \sum _{\mathit{m}=1}^{\mathit{M}}{\displaystyle \sum _{\mathit{h}=1}^{\mathit{H}}{\mathit{P}}_{\mathit{ht}}}}}}\ast \mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{-}\ast \mathit{DI}{\mathit{V}}_{\mathit{kmt}}\ast \mathit{DC}{\mathit{O}}_{\mathit{kmt}}$$

(52)

$$\left(7\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}{\displaystyle \sum _{\mathit{n}=1}^{\mathit{N}}{\displaystyle \sum _{\mathit{l}=1}^{\mathit{L}}\mathit{P}{\mathit{M}}_{\mathit{lt}}\ast \mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{-}}}}}\ast \mathit{DV}{\mathit{H}}_{\mathit{jnt}}\ast \mathit{DC}{\mathit{H}}_{\mathit{jnt}}$$

(53)

$$\left(8\right)={\displaystyle \sum _{\mathit{t}=1}^{12}\mathit{SA}{\mathit{E}}_{\mathit{t}}^{-}\cdot \mathit{CE}{\mathit{S}}_{\mathit{t}}^{-}}$$

(54)

$$\left(9\right)={\displaystyle \sum _{\mathrm{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=1}^{\mathit{K}}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}\left[\mathit{P}{\mathit{E}}_{1\mathit{kt}}^{-}\ast \left(1-{\mathit{\eta }}_{1\mathit{kt}}^{+}\right)+\mathit{P}{\mathit{H}}_{1\mathit{jt}}^{-}\ast \left(1-{\mathit{\mu }}_{1\mathit{jt}}^{+}\right)\right]\ast \mathit{PD}{\mathit{F}}_{\mathit{t}}^{-}}}}$$

(55)

Subject to:

$${\mathit{W}}_{\mathit{kt}}={\mathit{W}}_{\mathit{kt}}^{-}+\left({\mathit{W}}_{\mathit{kt}}^{+}-{\mathit{W}}_{\mathit{kt}}^{-}\right)\cdot {\mathit{u}}_{\mathit{kt}},\forall \mathit{k},\mathit{t}$$

(56)

$$0\le {\mathit{u}}_{\mathit{kt}}\le 1,\forall \mathit{k},\mathit{t}$$

(57)

$${\mathit{Y}}_{\mathit{jt}}={\mathit{Y}}_{\mathit{jt}}^{-}+\left({\mathit{Y}}_{\mathit{jt}}^{+}-{\mathit{Y}}_{\mathit{jt}}^{-}\right)\cdot {\mathit{u}}_{\mathit{jt}},\forall \mathit{j},\mathit{t}$$

(58)

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

(59)

$$\left({\mathit{W}}_{\mathit{kt}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{-}\right)\ast \mathit{EF}{\mathit{E}}_{\mathit{kt}}^{+}+\left({\mathit{Y}}_{\mathit{jt}}+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{-}\right)\ast \mathit{EF}{\mathit{H}}_{\mathit{jt}}^{+}+{\displaystyle \sum _{\mathit{s}=1}^{\mathit{S}}\mathit{E}{\mathit{D}}_{\mathit{its}}^{-}}\le {\mathit{T}}_{\mathit{it}}^{-},\forall \mathit{t},\mathit{k},\mathit{j},\mathit{l},\mathit{h},\mathit{i}$$

(60)

$$\left({\mathit{W}}_{3\mathit{t}}+\mathit{E}{\mathit{Q}}_{3\mathit{th}}^{-}+\mathit{SA}{\mathit{E}}_{\mathit{t}}^{-}\right)\ast \mathit{EF}{\mathit{E}}_{3\mathit{t}}^{+}\le {\mathit{T}}_{\mathit{it}}^{-},\forall \mathit{t},\mathit{h},\mathit{i}$$

(61)

$${\displaystyle \sum _{\mathit{k}=1}^3\left({\mathit{W}}_{\mathit{kt}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{-}\right)}+\mathit{SA}{\mathit{E}}_{\mathit{t}-1}^{-}\ge \mathit{DT}{\mathit{E}}_{\mathit{th}}^{-},\forall \mathit{t},\mathit{h}$$

(62)

$${\displaystyle \sum _{\mathit{k}=1}^3{\mathit{W}}_{\mathit{kt}}}-\mathit{DT}{\mathit{E}}_{\mathit{th}}^{+}+\mathit{SA}{\mathit{E}}_{\mathit{t}-1}^{-}=\mathit{SA}{\mathit{E}}_{\mathit{t}}^{-},\forall \mathit{t},\mathit{k},\mathit{h}$$

(63)

$${\mathit{W}}_{\mathit{kt}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{-}\le \mathit{E}{\mathit{T}}_{\mathit{kt}}^{-}\ast \left(\mathit{R}{\mathit{C}}_{\mathit{k}}+\mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{-}\ast \mathit{DI}{\mathit{V}}_{\mathit{kmt}}\right),\forall \mathit{t},\mathit{k},\mathit{h},\mathit{m}$$

(64)

$${\mathit{W}}_{3\mathit{t}}-\left[\mathit{DT}{\mathit{E}}_{\mathit{th}}^{+}-\left({\mathit{W}}_{1\mathit{t}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{-}\right)-\left({\mathit{W}}_{2\mathit{t}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{-}\right)\right]+\mathit{SA}{\mathit{E}}_{\mathit{t}-1}^{-}=\mathit{SA}{\mathit{E}}_{\mathit{t}}^{-},\forall \mathit{t}\ge 2,\mathit{k},\mathit{h}$$

(65)

$${\mathit{W}}_{3\mathit{t}}+\mathit{E}{\mathit{Q}}_{3\mathit{th}}^{-}-\left[\mathit{DT}{\mathit{E}}_{\mathit{th}}^{+}-\left({\mathit{W}}_{1\mathit{t}}+\mathit{E}{\mathit{Q}}_{1\mathit{th}}^{-}\right)-\left({\mathit{W}}_{2\mathit{t}}+\mathit{E}{\mathit{Q}}_{2\mathit{th}}^{-}\right)\right]+\mathit{SAE}{\mathrm{O}}^{-}=\mathit{SA}{\mathit{E}}_{\mathit{t}}^{-},\mathit{t}=1,\forall \mathit{h}$$

(66)

$${\displaystyle \sum _{\mathit{j}=1}^3\left({\mathit{Y}}_{\mathit{jt}}+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{-}\right)}\ge \mathit{TD}{\mathit{H}}_{\mathit{tl}}^{-},\forall \mathit{t},\mathit{l}$$

(67)

$${\mathit{Y}}_{\mathit{jt}}+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{-}\le \mathit{H}{\mathit{T}}_{\mathit{jt}}^{-}\ast \left(\mathit{RC}{\mathit{H}}_{\mathit{j}}+\mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{-}\ast \mathit{DV}{\mathit{H}}_{\mathit{jnt}}\right),\forall \mathit{t},\mathit{j},\mathit{n},\mathit{l}$$

(68)

$$\mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{-}=\left\{\begin{array}{l}\le 1\\ \ge 0\end{array}\right.,\forall \mathit{k},\mathit{t},\mathit{m},\mathit{h}$$

(69)

$${\displaystyle \sum _{\mathit{m}=1}^3\mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{-}}\le 1,\forall \mathit{k},\mathit{t},\mathit{h}$$

(70)

$$\mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{-}=\left\{\begin{array}{l}\le 1\\ \ge 0\end{array}\right.,\forall \mathit{j},\mathit{t},\mathit{n},\mathit{l}$$

(71)

$${\displaystyle \sum _{\mathit{n}=1}^3\mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{-}}\le 1,\forall \mathit{j},\mathit{t},\mathit{l}$$

(72)

$${\displaystyle \sum _{\mathit{k}=1}^3\left({\mathit{W}}_{\mathit{kt}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{-}\right)}\ast \mathit{P}{\mathit{E}}_{\mathit{krt}}^{-}\ast \left(1-{\mathit{\eta }}_{\mathit{krt}}^{-}\right)\le \mathit{T}{\mathit{P}}_{\mathit{rt}}^{-},\forall \mathit{t},\mathit{h},\mathit{r}$$

(73)

$${\displaystyle \sum _{\mathit{j}=1}^3\left({\mathit{Y}}_{\mathit{jt}}+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{-}\right)}\ast \mathit{P}{\mathit{H}}_{\mathit{jqt}}^{-}\ast \left(1-{\mathit{\mu }}_{\mathit{jqt}}^{-}\right)\le \mathit{TH}{\mathit{P}}_{\mathit{qt}}^{-},\forall \mathit{t},\mathit{l},\mathit{q}$$

(74)

$${\displaystyle \sum _{\mathit{s}=1}^{\mathit{S}}{\displaystyle \sum _{\mathit{i}=1}^2\mathit{E}{\mathit{D}}_{\mathit{its}}^{-}}}\ast \mathit{IN}{\mathit{T}}_{\mathit{s}}^{+}+{\displaystyle \sum _{\mathit{k}=1}^3({\mathit{W}}_{\mathit{kt}}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{-})\ast \mathit{CO}{\mathit{E}}_{\mathit{kt}}^{+}+{\displaystyle \sum _{\mathit{j}=1}^3({\mathit{Y}}_{\mathit{jt}}+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{-}})\ast \mathit{CO}{\mathit{H}}_{\mathit{jt}}^{+}\le \mathit{T}{\mathit{C}}_{\mathit{t}}^{-},\forall \mathit{t},\mathit{h},\mathit{l}$$

(75)

Where ${\mathit{T}}_{\mathit{it}}^{-}$, $\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{-}$, $\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{-}$, u _kt and v _jt are continuous decision variables, $\mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{\pm }$ and $\mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{\pm }$ are binary ones. Solution for f ⁻ reach a extreme lower bound of system cost under uncertainties. Then, the optimized electricity and heat targets would be ${\mathit{W}}_{\mathit{ktopt}}={\mathit{W}}_{\mathit{kt}}^{-}+\mathrm{\Delta }{\mathit{W}}_{\mathit{kt}}\cdot {\mathit{u}}_{\mathit{ktopt}}$, ${\mathit{Y}}_{\mathit{jtopt}}={\mathit{Y}}_{\mathit{jt}}^{-}+\mathit{\Delta }{\mathit{Y}}_{\mathit{jt}}\cdot {\mathit{v}}_{\mathit{jtopt}}$. Conversely, the upper bound of the objective function value f ⁺ is:

$$min{\mathit{f}}^{+}=\left(1\right)+\left(2\right)+\left(3\right)+\left(4\right)+\left(5\right)+\left(6\right)+\left(7\right)+\left(8\right)+\left(9\right)$$

(76)

$$\left(1\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{i}=1}^2\mathit{E}{\mathit{S}}_{\mathit{it}}^{+}}}\ast \mathit{E}{\mathit{P}}_{\mathit{it}}^{+}$$

(77)

$$\left(2\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=1}^3\mathit{T}{\mathit{C}}_{\mathit{kt}}^{+}\ast {\mathit{W}}_{\mathit{kt}}+{\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=1}^3{\displaystyle \sum _{\mathit{h}=1}^{\mathit{H}}{\mathit{P}}_{\mathit{ht}}\ast \left(\mathit{T}{\mathit{C}}_{\mathit{kt}}^{+}+\mathit{P}{\mathit{P}}_{\mathit{kt}}^{+}\right)\ast \mathit{E}{\mathit{Q}}_{\mathit{kth}}^{+}}}}}}+\mathit{T}{\mathit{C}}_{3\mathit{t}}^{+}\ast \left({\mathit{W}}_{3\mathit{t}}+\mathit{SA}{\mathit{E}}_{\mathit{t}}^{+}\right)$$

(78)

$$\left(3\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}\mathit{TP}{\mathit{H}}_{\mathit{jt}}^{+}\ast {\mathit{Y}}_{\mathit{jt}}+{\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}{\displaystyle \sum _{\mathit{L}=1}^{\mathit{l}}{\mathit{PM}}_{\mathit{lt}}\ast \left(\mathit{TP}{\mathit{H}}_{\mathit{jt}}^{+}+\mathit{PC}{\mathit{H}}_{\mathit{jt}}^{+}\right)\ast \mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{+}}}}}}$$

(79)

$$\left(4\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=\mathit{i}}^{\mathit{K}}{\displaystyle \sum _{\mathit{h}=1}^{\mathit{H}}{\displaystyle \sum _{\mathit{r}=1}^{\mathit{R}}({\mathit{W}}_{\mathit{kt}}}}}}+{\mathit{P}}_{\mathit{ht}}\ast \mathit{E}{\mathit{Q}}_{\mathit{kth}}^{+})\ast \mathit{P}{\mathit{E}}_{\mathit{rkt}}^{+}\ast \mathit{C}{\mathit{E}}_{\mathit{rkt}}^{+}$$

(80)

$$\left(5\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}{\displaystyle \sum _{\mathit{l}=1}^{\mathit{L}}{\displaystyle \sum _{\mathit{q}=1}^{\mathit{Q}}\left({\mathit{Y}}_{\mathit{jt}}+\mathit{P}{\mathit{M}}_{\mathit{lt}}\ast \mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{+}\right)\ast \mathit{P}{\mathit{H}}_{\mathit{qjt}}^{+}\ast \mathit{C}{\mathit{H}}_{\mathit{qjt}}^{+}}}}}$$

(81)

$$\left(6\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=1}^{\mathit{K}}{\displaystyle \sum _{\mathit{m}=1}^{\mathit{M}}{\displaystyle \sum _{\mathit{h}=1}^{\mathit{H}}{\mathit{P}}_{\mathit{ht}}}}}}\ast \mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{+}\ast \mathit{DI}{\mathit{V}}_{\mathit{kmt}}\ast \mathit{DC}{\mathit{O}}_{\mathit{kmt}}$$

(82)

$$\left(7\right)={\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}{\displaystyle \sum _{\mathit{n}=1}^{\mathit{N}}{\displaystyle \sum _{\mathit{l}=1}^{\mathit{L}}\mathit{P}{\mathit{M}}_{\mathit{lt}}\ast \mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{+}}}}}\ast \mathit{DV}{\mathit{H}}_{\mathit{jnt}}\ast \mathit{DC}{\mathit{H}}_{\mathit{jnt}}$$

(83)

$$\left(8\right)={\displaystyle \sum _{\mathit{t}=1}^{12}\mathit{SA}{\mathit{E}}_{\mathit{t}}^{+}\cdot \mathit{CE}{\mathit{S}}_{\mathit{t}}^{+}}$$

(84)

$$\left(9\right)={\displaystyle \sum _{\mathrm{t}=1}^{12}{\displaystyle \sum _{\mathit{k}=1}^{\mathit{K}}{\displaystyle \sum _{\mathit{j}=1}^{\mathit{J}}\left[\mathit{P}{\mathit{E}}_{1\mathit{kt}}^{+}\ast \left(1-{\mathit{\eta }}_{1\mathit{kt}}^{-}\right)+\mathit{P}{\mathit{H}}_{1\mathit{jt}}^{+}\ast \left(1-{\mathit{\mu }}_{1\mathit{jt}}^{-}\right)\right]\ast \mathit{PD}{\mathit{F}}_{\mathit{t}}^{+}}}}$$

(85)

subject to:

$$\left({\mathit{W}}_{\mathit{ktopt}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{+}\right)\ast \mathit{EF}{\mathit{E}}_{\mathit{kt}}^{-}+\left({\mathit{Y}}_{\mathit{jtopt}}+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{+}\right)\ast \mathit{EF}{\mathit{H}}_{\mathit{jt}}^{-}+{\displaystyle \sum _{\mathit{s}=1}^{\mathit{S}}\mathit{E}{\mathit{D}}_{\mathit{its}}^{+}}\le {\mathit{T}}_{\mathit{it}}^{+},\forall \mathit{t},\mathit{k},\mathit{j},\mathit{l},\mathit{h},\mathit{i}$$

(86)

$$\left({\mathit{W}}_{3\mathit{topt}}+\mathit{E}{\mathit{Q}}_{3\mathit{th}}^{+}+\mathit{SA}{\mathit{E}}_{\mathit{t}}^{+}\right)\ast \mathit{EF}{\mathit{E}}_{3\mathit{t}}^{-}\le {\mathit{T}}_{\mathit{it}}^{+},\forall \mathit{t},\mathit{h},\mathit{i}$$

(87)

$${\displaystyle \sum _{\mathit{k}=1}^3\left({\mathit{W}}_{\mathit{ktopt}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{+}\right)}+\mathit{SA}{\mathit{E}}_{\mathit{t}-1}^{+}\ge \mathit{DT}{\mathit{E}}_{\mathit{th}}^{+},\forall \mathit{t},\mathit{k},\mathit{h}$$

(88)

$${\displaystyle \sum _{\mathit{k}=1}^3{\mathit{W}}_{\mathit{ktopt}}}-\mathit{DT}{\mathit{E}}_{\mathit{th}}^{+}+\mathit{SA}{\mathit{E}}_{\mathit{t}-1}^{+}=\mathit{SA}{\mathit{E}}_{\mathit{t}}^{+},\forall \mathit{t},\mathit{k},\mathit{h}$$

(89)

$${\mathit{W}}_{\mathit{ktopt}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{+}\le \mathit{E}{\mathit{T}}_{\mathit{kt}}^{+}\ast \left(\mathit{R}{\mathit{C}}_{\mathit{k}}+{\displaystyle \sum _{\mathit{t}=1}^{12}{\displaystyle \sum _{\mathit{m}=1}^{\mathit{M}}\mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{+}\ast \mathit{DI}{\mathit{V}}_{\mathit{kmt}}}}\right),\forall \mathit{t},\mathit{k},\mathit{h},\mathit{m}$$

(90)

$${\mathit{W}}_{3\mathit{topt}}-\left[\mathit{DT}{\mathit{E}}_{\mathit{th}}^{-}-\left({\mathit{W}}_{1\mathit{topt}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{+}\right)-\left({\mathit{W}}_{2\mathit{topt}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{+}\right)\right]+\mathit{SA}{\mathit{E}}_{\mathit{t}-1}^{+}=\mathit{SA}{\mathit{E}}_{\mathit{t}}^{+},\forall \mathit{t}\ge 2,\mathit{k},\mathit{h}$$

(91)

$${\mathit{W}}_{3\mathit{topt}}+\mathit{E}{\mathit{Q}}_{3\mathit{th}}^{+}-\left[\mathit{DT}{\mathit{E}}_{\mathit{th}}^{-}-\left({\mathit{W}}_{1\mathit{topt}}+\mathit{E}{\mathit{Q}}_{1\mathit{th}}^{+}\right)-\left({\mathit{W}}_{2\mathit{topt}}+\mathit{E}{\mathit{Q}}_{2\mathit{th}}^{+}\right)\right]+\mathit{SAE}{\mathrm{O}}^{+}=\mathit{SA}{\mathit{E}}_{\mathit{t}}^{+},\mathit{t}=1,\forall \mathit{h}$$

(92)

$${\displaystyle \sum _{\mathit{j}=1}^3\left({\mathit{Y}}_{\mathit{jtopt}}+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{+}\right)}\ge \mathit{TD}{\mathit{H}}_{\mathit{tl}}^{+},\forall \mathit{t},\mathit{l}$$

(93)

$${\mathit{Y}}_{\mathit{jtopt}}+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{+}\le \mathit{H}{\mathit{T}}_{\mathit{jt}}^{+}\ast \left(\mathit{RC}{\mathit{H}}_{\mathit{j}}+\mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{+}\ast \mathit{DV}{\mathit{H}}_{\mathit{jnt}}\right),\forall \mathit{t},\mathit{j},\mathit{n},\mathit{l}$$

(94)

$$\mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{+}=\left\{\begin{array}{l}\le 1\\ \ge 0\end{array}\right.,\forall \mathit{k},\mathit{t},\mathit{m},\mathit{h}$$

(95)

$${\displaystyle \sum _{\mathit{m}=1}^3\mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{+}}\le 1,\forall \mathit{k},\mathit{t},\mathit{h}$$

(96)

$$\mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{+}=\left\{\begin{array}{l}\le 1\\ \ge 0\end{array}\right.,\forall \mathit{j},\mathit{t},\mathit{n},\mathit{l}$$

(97)

$${\displaystyle \sum _{\mathit{n}=1}^3\mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{+}}\le 1,\forall \mathit{j},\mathit{t},\mathit{l}$$

(98)

$${\displaystyle \sum _{\mathit{k}=1}^3\left({\mathit{W}}_{\mathit{ktopt}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{+}\right)}\ast \mathit{P}{\mathit{E}}_{\mathit{krt}}^{+}\ast \left(1-{\mathit{\eta }}_{\mathit{krt}}^{-}\right)\le \mathit{T}{\mathit{P}}_{\mathit{rt}}^{+},\forall \mathit{t},\mathit{h},\mathit{r}$$

(99)

$${\displaystyle \sum _{\mathit{j}=1}^3\left({\mathit{Y}}_{\mathit{jtopt}}+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{+}\right)}\ast \mathit{P}{\mathit{H}}_{\mathit{jqt}}^{+}\ast \left(1-{\mathit{\mu }}_{\mathit{jqt}}^{-}\right)\le \mathit{TH}{\mathit{P}}_{\mathit{qt}}^{+},\forall \mathit{t},\mathit{l},\mathit{q}$$

(100)

$${\displaystyle \sum _{\mathit{s}=1}^{\mathit{S}}{\displaystyle \sum _{\mathit{i}=1}^2\mathit{E}{\mathit{D}}_{\mathit{its}}^{+}}}\ast \mathit{IN}{\mathit{T}}_{\mathit{s}}^{-}+{\displaystyle \sum _{\mathit{k}=1}^3({\mathit{W}}_{\mathit{ktopt}}}+\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{+})\ast \mathit{CO}{\mathit{E}}_{\mathit{kt}}^{-}+{\displaystyle \sum _{\mathit{j}=1}^3({\mathit{Y}}_{\mathit{jtopt}}+\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{+}})\ast \mathit{CO}{\mathit{H}}_{\mathit{jt}}^{-}\le \mathit{T}{\mathit{C}}_{\mathit{t}}^{+},\forall \mathit{t},\mathit{h},\mathit{l}$$

(101)

$${\mathit{T}}_{\mathit{it}}^{+}\ge {\mathit{T}}_{\mathit{itopt}}^{-},\forall \mathit{i},\mathit{t}$$

(102)

$$\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{+}\ge \mathit{E}{\mathit{Q}}_{\mathit{kthopt}}^{-},\forall \mathit{k},\mathit{t},\mathit{h}$$

(103)

$$\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{+}\ge \mathit{H}{\mathit{Q}}_{\mathit{jtlopt}}^{-},\forall \mathit{j},\mathit{t},\mathit{l}$$

(104)

$$\mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{\pm }\ge \mathit{JU}{\mathit{E}}_{\mathit{ktmhopt}}^{-},\forall \mathit{k},\mathit{t},\mathit{m},\mathit{h}$$

(105)

$$\mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{\pm }\ge \mathit{JU}{\mathit{H}}_{\mathit{jtnlopt}}^{-},\forall \mathit{j},\mathit{t},\mathit{n},\mathit{l}$$

(106)

where ${\mathit{T}}_{\mathit{it}}^{+}$, $\mathit{E}{\mathit{Q}}_{\mathit{kth}}^{+}$ and $\mathit{H}{\mathit{Q}}_{\mathit{jtl}}^{+}$ are continuous decision variables, and $\mathit{JU}{\mathit{E}}_{\mathit{ktmh}}^{+}$ and $\mathit{JU}{\mathit{H}}_{\mathit{jtnl}}^{+}$ are binary ones; ${\mathit{T}}_{\mathit{itopt}}^{-}$, $\mathit{E}{\mathit{Q}}_{\mathit{kthopt}}^{-}$, $\mathit{H}{\mathit{Q}}_{\mathit{jtlopt}}^{-}$, $\mathit{JU}{\mathit{E}}_{\mathit{ktmhopt}}^{-}$ and $\mathit{JU}{\mathit{H}}_{\mathit{jtnlopt}}^{-}$ are solutions of the model B. As a result, the solutions for electricity and heat supply model under the optimized targets could be obtained by incorporating the solutions of the above two sub-models.

Case study

The following regional electricity and heat supply management problem is used to demonstrate the applicability of the developed ITSDP model. In this study system, a decision maker is responsible for electricity and heat to multiple users through several technologies on an annual cycle based on different pollutants emission permissible levels. Since different environmental restrictions, the decision maker would take expected electricity and heat demand into account to make fossil and renewable resources utilization optimal, manage electricity and heat generation and plan the equipment expansion. Mostly, increasing electricity and heat demand could lead to facility expansion, fuel exploitation and energy imports. However, it is unpractical to achieve sustainable development with the ever-increasing economic and environmental costs, unlimited facility expansion and energy exploitation (Liu et al. [2000]). Thus, the problem is how to incorporate different pollutants emission reduction control levels into regional electricity and heat planning. In this study system, planning period is on an annual cycle, with heating seasons and no heating seasons. Multiple energy resources/technologies need to be allocated to multiple end users. Traditional energy resources (e.g., coal, natural gas) with limited availabilities and renewable resources (e.g., wind) are utilized to meet the electricity and heat demand. More specifically, coal and natural gas are both used for electricity generation; wind is a supplement for electricity generation; the cogeneration system of electricity and heat is another energy-efficient way. The electricity and heat demands during one year are affected by many uncertainties, such as the growing population, energy-transformation rate, electricity and heat prices and changing weather; all of those factors would produce many complicated uncertainties (Li et al. Li et al. [2006a], [b]; Liu [2007]; Li et al. [2008]). In addition, these uncertainties are intricate because of a series of imprecise information (e.g., social, economic, environmental, seasonal and geographic conditions, energy carrier characteristics). They could only be expressed as distribution information or intervals when most data can hardly be available determinately. Once these uncertainties are determined, system costs, efficiencies and capacities of each technology could be defined.

Under different environmental quality standards, decision maker is responsible for these followings (Lehtila and Pirila [1996]):

• (1) Assigning electricity load to four convention technologies (e.g., coal-fired power, natural-gas-fired power, wind power and cogeneration of power and heat) and heat load to three convention technologies (e.g., coal boiler heating generation, natural-gas boiler technology and cogeneration of power and heat);

• (2) Planning the facility expansion of electricity and heat generation;

• (3) Managing the fossil fuels purchases (including coal and natural gas).

If actual supply could not meet users’ demands, decision maker would invest more funds on capacity expansion or turn to purchase shortfall from other enterprises with higher costs; both of the programs would make the region economic interests punished.

In this study, sulfur dioxide, nitrogen oxide and particulate matters (PM10) are considered as main gaseous emission generated from fuel combustion in electricity and heat conversion; based on the gaseous pollutants emission intensity and specific demand of fuel in power and heat supply, the total amount of sulfur dioxide, nitrogen oxide and particulate matters (PM10) emission could be calculated. The fuel demand has been predicted in each month, whether or not consider pollutants emission reduction, fuel supply would be determined to meet users demand. Therefore, under the scenarios of pollutants emission reduction, three electricity generation processes are forced to decarbonize in energy system. At the same time, pollutants emission intensity are affected by many uncertain factors (e.g., pollutants emission inventory, reduction control measures, related costs, weather situation), which can be presented as interval numbers without distribution information. The availabilities of electricity and heat demand are directly influenced by natural fluctuations, which can be expressed as probability distributions. Most of other parameters (energy demand, technological convention efficiency, and utilization factors) are indicated as intervals. Tables 1 and 2 show the available electricity demands and heat demands under different probability distributions; Table 3 presents some corresponding economic data; gaseous pollutants emission intensity is listed in Table 4. Besides, coal-fired power has a residual capacity of 0.9 GW, natural-gas-fired power has a residual capacity of 0.6 GW, CHP (Combined Heat and Power) has a residual capacity of 0.75 GW, wind power has a residual capacity of 0.27 GW; coal boiler heating generation has a residual capacity of 1.8 TJ, natural-gas boiler technology has a residual capacity of 1 TJ and cogeneration of heat and power has a residual of capacity of 1.5 TJ. The technical data and representative costs are investigated in government report and other related documents (USDE and USEPA [2000]; Iniyan and Sumath [2000]; Cormio et al. [2003]; Kristoffersen [2007]; Fleten and Kristoffersen [2008]; Cai et al. [2008], [2009]; Carla and Carlos [2011]).

Table 1 Electricity demands under different probability distributions

Table 2 Heat demands under different probability distributions

Table 3 Economic data

Table 4 Parameters of contamination

Results and discussion

The objective of the ITSDP model is to minimize the costs of the overall system under different gaseous pollutants reduction levels during one year. Solutions provide an effective relationship between the preexisting set of environmental standardize policies and the associated economic interests, for example, improper policies would cause losses and penalties. The results contain a combination of some deterministic, interval and distributional information, it means, therefore, different forms of uncertainties could be reflected out with less error (Li et al., Li et al. [2006a],[b]). The interval results provide managers multiple decision alternatives and form the basis for further study of trade-offs between electricity and heat management cost and so-caused gaseous pollutants emission reduction in electricity and heat generation; the binary variable solutions stand for the decisions of facility expansion, so that several alternative schemes are generated; the continuous variable solutions are interrelated to electricity and heat generation and energy resources supply.

Solutions without pollutants emission reduction considered

This case is raised as a reference to show a pattern of energy production and system development without constrains on gaseous pollutants emission reduction (that is 0% pollutants-emission reduction).

Figure 1 shows the results of electricity schemes under this case, from which we can see thermoelectricity would be the largest energy resource among all electricity generation options during the whole year. Because combined heat and power (CHP) is a key energy-saving project as well as technology popularization system advocated by government. The next most likely methods are, in order, coal-fired power and natural gas-fired power; as the traditional energy resources, they are still an essential part of power supply and share the load. Wind power would provide considerable amounts of electricity as a kind of clean energy; however, the weak wind power system and long distance from load center limit the development of wind farm, which would have a direct impact on wind installed capacity. As seen from Figure 1, electricity production would fluctuate smoothly in one year; overall, the main way of electricity production influenced by some factors (e.g., economy growth, colder/warmer season) would seasonally vary. The pre-regulated coal-fired power generation would be increased from [91.27, 124.30] GWh in February to [157.60, 207.42] GWh in May; for the natural gas-fired power, its plan would be [74.02, 95.36] GWh in March and up to [146.59, 158.32] GWh in July; for the CHP, its pre-regulated targets would be decreased from [247.96, 283.61] GWh in January to [181.97, 226.59] GWh in August. The pre-regulated wind power targets would remain 54 GWh over the year. The solution of planned heat generation is shown in Figure 2. If the planned electricity and heat cannot meet the random demand, the insufficient electricity has to be produced under different demand levels.

Figure 1

Optimized electricity generation in different periods.

Figure 2

Optimized heat generation in different periods.

Gaseous pollutants (e.g., SO₂, NO_X, PM10) emission associated with energy exploitation mostly from the burning of fuels can be categorized into generation in power supply and generation in heat supply. Figure 3 and Figure 4 present the detailed solutions of these three pollutants emission during a year. Different energy resources would be supplied to these electricity and heat generation technologies to meet the demands; the amount of pollutants emission would be related to energy activities. From Figure 3, the greatest output pollutant would be NO_X in electricity and SO₂ in heat from Figure 4; PM10 would be a small value. SO₂ and NO_X are still the main pollutants; in addition, CO₂ would be under control as a major kind of greenhouse gases. The coal-fired power generation and coal boiler heat generation technologies would be the largest pollutants and greenhouse gas emission source.

Figure 3

Quantity of pollutants emission in electricity generation.

Figure 4

Quantity of pollutants emission in heat generation.

Solutions under pollutants emission reduction

In this study, two scenarios of pollutants-emission reduction are considered (i.e. 10%, and 20% of total pollutants -emission reduction). The results indicate that increased substantive capacity expansion investment for clean energy (to reduce pollutants emissions) could lead to an increased system cost. The results indicate that the increased capacity expansion investment for clean energy (proven to be effective for pollutants emission reduction) could lead to a higher system cost. Based on the results of the scenario of 0% pollutants emission reduction, coal and natural gas would be supplied to guarantee the electricity and heat demand. Next, Figure 5 shows clean energy resources storage. Compared with the result that no considering pollutants emission reduction, the amount of wind power would be largely decreased, especially, excess wind power would be stored to meet electricity demand in the future. The utilization of wind power under scenarios of 10% and 20% pollutants emission reduction would increase, compared to the consumption under 0% pollutants emission reduction condition. That is to say, under the scenarios of 10% and 20% pollutants reduction, the fuel-fired power would decrease with pollutants reduction increasing. For example, in April, electricity generated from coal-fired power conversion technology would be [157.60, 207.42] GWh under 0% pollutants emission reduction, [146.30, 159.60] GWh under 10% pollutants emission reduction and [108.85, 145.83] GWh under 20% pollutants emission reduction. Compared the supplies of two energy resources under pollutants emission condition, natural gas supplies would be higher than coal supplies. Therefore, it recommends that natural gas would be more popular than coal in a condition of considering the pollutants emission reduction. This is because the amount of pollutants emission would be confined in a certain level and coal-fired electricity conversion technology corresponds to a higher pollutants emission rate than natural gas-fired conversion technology. In addition, the increased clean energy (wind power) electricity conversion technology and power conversion technology with higher energy utilization efficiency (CHP) also share the excess electricity load caused by the decreased coal supplies to meet different pollutants emission-cutting policies.

Figure 5

Storage of wind power under 20% pollutants emission reduction.

Discussion

Compared the contribution of various electricity generation technologies to the medium power demand, it presents that different electricity conversion technologies have various generation quantities under changed pollutants emission- cutting scenarios. As the above analysis, in the respect of energy type, coal-fired electricity would be the mainly power supply source under 0% pollutants emission reduction. Natural gas-fired electricity conversion technology would play a key role in the power generation activities, coal-fired power would be in a second place and wind power would be the supplement under 10% and 20% pollutants emissions reduction. This is because coal-fired electricity conversion technology has comparatively low operating and penalty costs and relatively low capital cost for capacity expansion, and the related cost of natural gas-fired electricity conversion technology is a little higher than coal-fired power and the pollutants emission from natural gas-fired power generation process is smaller relatively. The maximum optimized wind power generation would be 0.5 GWh in April, July and October, due to the relatively higher operating cost and capital cost for its storage, which limits the development of wind power. The dominant role of coal-fired electricity would be displaced by some other conversion technologies with an increased demand for pollutants emission reduction. For example, under 20% pollutants emission reduction, the optimized coal-fired power target would decrease to [39.94, 46.20] GWh in January; although fuel-fired power would decrease, the optimized targets of wind power would have a growth and play an important role to meet electricity demand. It indicates that environment-friendly electricity conversion technologies should better serve interests of the pollutants emission reduction and be chosen for power generation to meet the ever-increasing electricity demands and gradually enhancing pollutants emission reduction requirements. These are optimized conclusions in terms of pollutants emission-cutting; from view point of energy saving, CHP would be properly adopted in electricity and heat generation progress.

As shown in Figures 6, 7, 8, the system cost would rise up along with the growing intensity of pollutants emission reduction. Without pollutants emission reduction, the system cost would be $[172.20, 269.96] × 10⁶, and the system cost would become $[187.64, 297.45] × 10⁶, $[242.33, 471.41] × 10⁶ under 10% and 20% pollutants emission reduction respectively. One of the major reason is that the traditional power and heat generation technologies (coal-fired and natural gas-fired) would gradually be replaced by lean energy (wind power)and more effective technologies (CHP) when the restrictions on pollutants emission are considered. Besides, the growing electricity and heat demands also lead to some degree of electricity and heat generating facilities to be expanded, resulting in a high capital cost.

Figure 6

System cost under 20% pollutants emission reduction.

Figure 7

System cost under 10% pollutants emission reduction.

Figure 8

System cost under 0% pollutants emission reduction.

Without ILP, the pollutants emission management and programming problem can also be solved by fixed-mix stochastic programming approach though adopting their mid-point values instead of the interval parameters. Undoubtedly further sensitivity analysis could be undertaken, but the model still cannot efficiently reflect the interactions among these uncertainties because each solution can only provide a single result corresponding to variations of the uncertain inputs. Likewise, if best/worst sub-models are solved, we can only obtain two solutions under extreme scenarios (best condition and worst condition). They are serviceable to judge the possibility that desired goal realizes but not necessarily build a series of stable intervals for decision variables.

Thus, the best/worst case analysis is not efficiently useful to produce decision alternatives. Actually, it is a special type of sensitivity analysis for extreme conditions (Huang et al. [2001]).

From the above analysis, we can conclude that the solutions worked out from the ITSDP are value for making decisions of energy resources allocation, capacity expansion of power and heat generation as well as pollutants emission management. In the progress of decision alternatives generation, the interval solutions are effective to represent various options which reflect environmental-economic trades off. Cost-effective options can be obtained with a least-cost strategy though planning pollutants emission management in regional electricity and heat supply. However, if pollutants emission reduction is considered, the pre-regulated targets of electricity and heat supply from various technologies have to be reallocated, the pre-arranged capacity expansion options of power and heat generation technologies would be reselected.

Conclusions

An inexact two-stage dynamic programming (ITSDP) model has been developed for planning pollutants emission management and electricity and heat supply systems under uncertainties. This method based on interval-parameter programming and two-stage stochastic programming so that it allows uncertainties expressed as both interval values and probability distributions to be incorporated within a common optimization framework. Furthermore, ITSDP also addresses various dynamics such as capacity expansion, storage of wind power and pollutants emission reduction scenarios related to different levels of economic implications. Probability distributions of power and heat demand can be integrated into a optimization process under a train of fixed levels by introducing fix-stochastic programming, which has significance advantages in reflecting uncertainties in large-scale problems. Accordingly, the developed model has been applied to a case of electricity and heat supply management planning. The results of this case study suggest that the model is appropriate for reflecting complexities of regional electricity and heat supply management systems and incorporating pollutants emission reduction issue during one year.

The prospect method could help decision-makers identify required management policies under changeable environmental and economic considerations. Nevertheless, there still exists much space for improvement of this model. As a powerful supplement to TSP, ITSDP can reflect the dynamic variations such as sequential structure and simplify large quantities design scenarios in case the problem of “dimension disaster” happen. This study is attempted to integrate IPP and TSP methods into a general framework, and apply the ITSDP for pollutants emission reduction management under uncertainty. In addition, the optimization algorithm is of value on other environmental issues or energy program problems within complex uncertainties. There are also other programming techniques can be integrated with ITSDP for handling more complicated cases, fuzzy programming for instance.