- Email: [email protected]
Expansion planning model of multi-energy system with the integration of active distribution network
Applied Energy 253 (2019) 113517
Contents lists available at ScienceDirect
Applied Energy journal homepage: www.elsevier.com/locate/apenergy
Expansion planning model of multi-energy system with the integration of active distribution network Jueying Wang, Zhijian Hu , Shiwei Xie ⁎
T
⁎
School of Electrical Engineering and Automation, Wuhan University, Wuhan City, Hubei Province, China
HIGHLIGHTS
for multi-energy system integrating active distribution network is modeled. • Planning scenario generation method is proposed considering the robustness. • AA probabilistic modified piecewise linearization method is developed. • Mixed integer second order cone programming is employed to solve the model. • ARTICLE INFO
ABSTRACT
Keywords: Multi-energy system expansion planning Active distribution network Active network management Probabilistic scenario generation Second-order cone programming
As an effective pattern to promote efficient use of energy, multi-energy system has become the main focus of development. On the other hand, active distribution network is now also the promoting concept for developing power distribution system, since the active network management technologies can make the network more flexible and controllable. To facilitate and utilize the advantages of both systems, this paper proposes an expansion planning model for multi-energy system integrating active distribution network, natural gas network and energy hub; and the positive impact of active network managements on expanding multi-energy system is originally investigated. With the aim of minimizing the total cost over the planning horizon, this model centers on the optimal determination of the type, location and size of all the infrastructures, where the active network managements are modeled and incorporated. In addition, in order to improve the robustness of the planning results, a probabilistic scenario generation approach is proposed to test the model. In order to solve the proposed mixed-integer nonlinear model, the second order cone programming as well as a modified piecewise linearization approach is applied to convert the original model to a mixed integer second order cone programming model. A multi-energy system (including a modified IEEE 33-node distribution system, a new 23-node gas system and 8 energy hubs) are employed to verify the effectiveness of the prosed model and methods. The simulation results exhibit the superiority of the joint expansion planning model for multi-energy system and the beneficial impact of considering active network management.
1. Introduction Global energy crisis and its relevant problems have led to a considerable movement into efficient utilization of renewable energy [1]. Therefore, making progress in supplying a sustainable, clean, secure and affordable energy is one of the basic challenges for industry in the 21st century [2]. In order to seek ways to deal with these challenges, the concept of multi-energy system (MES) [3] is introduced as an effective pattern to improve the energy efficiency and reduce energy loss by incorporating the multiple forms of energy. However, utilizing the benefits of MES requires an integrated modeling framework. Currently,
⁎
a new window has been opened in the field of MES optimization by presenting the promising concept of energy hub (EH), which is defined as a framework where production, consumption, conversion and storage of different forms of energy carriers can be done [4]. Based on the EH concept, energy systems can be investigated and controlled as a whole rather than as separate energy systems like distribution network, natural gas network (NGN) and so on. In addition, emerging new technologies such as combined heat and power (CHP) and power to gas (P2G) have made it possible to convert different forms of energy into each other, and made the EH concept practical for MES [5]. At present, the research on multi-energy system expansion planning
Corresponding authors. E-mail addresses: [email protected] (Z. Hu), [email protected] (S. Xie).
https://doi.org/10.1016/j.apenergy.2019.113517 Received 31 January 2019; Received in revised form 13 June 2019; Accepted 9 July 2019 0306-2619/ © 2019 Elsevier Ltd. All rights reserved.
Applied Energy 253 (2019) 113517
J. Wang, et al.
Nomenclature
units/thermal storage units/electric storage systems c M _lori/ ckM _lup/ ckM _lnew maintenance cost coefficients of original/replacing/feeder ckM _wt / ckM _pv / ckM _cb/ maintenance cost coefficients of WTG/PVG ckM _svc / ckM _vr /CB/SVC/VR gas price coefficient 1/ 0 / 1/ 0 electricity price coefficient m/ wt / pv power factor limits for WTG/PVG control c c ratios of existent/adding compressors ij / ij, k transmission coefficients through pipelines ij
Indices
i, j ij, jk k /t /s
indices for nodes indices for branches index for types of each component index for time periods/stages/scenarios
Sets set of existent nodes/pipelines/source nodes in the gas network set of existent nodes connected with energy hub HUB E _B / Lori / E _Sup set of existing fixed buses/feeders/electricity source buses G _B / E _B set of newly-added gas/distribution network nodes WTG / PVG / CB / SVC / VR set of candidate buses for WTG/PVG/ CB/ SVC/VR Lnew / Lup set of candidate branches for addition/replacement of feeders Pipe / C set of candidate branches for adding pipelines/compressors GS / HS / ESS set of candidate nodes for gas storage units/thermal storage units/electric storage systems set of candidate nodes for thermal storage units HS G _B / Pipe set of all nodes/branches in the gas network, G _B = { G _B , G _B } , Pipe = { Pipe , Pipe } E _B / Line / L set of all buses/branches/load nodes in the distribution network, E _B = { E _B , E _B } , Line = { Lori , Lup , Lnew } KPipe/ K C set of alternative types of adding pipelines/compressors K GS /KHS /KESS set of alternative types of adding gas storage units/ thermal storage units/electric storage systems KLup/ KLnew set of alternative types for replacing/adding feeders KWTG / KPVG/ K CB /K SVC / KVR set of alternative types for WTG/PVG/ CB/SVC/VR JK (j )/ IJ (j) set of feeders whose parent/child is j G _B / Pipe / G _Sup
energy conversion efficiency of P2G/CHP/GF storage/release efficiency of gas storage units sh / rh heat storage/release efficiency of thermal storage units se / re charge/discharge efficiency of ESS DG DG max / min maximum/minimum allowable DG permeability Cur _wt Cur _pv Cur _L / max / max maximum curtailment rate of wind/solar max power/electrical load kijvr incremental ratio between VR taps Rijvr maximum gear of VR j,max / j,min upper/lower limits for gas nodal PRESSURE C FijPipe ,max / Fij,max flow limits of existing pipes/compressors C GS FijPipe / F , k ,max ij, k ,max / Fk,max flow limits of adding pipes/compressors/gas storage units Sup F jSup ,max / F j,min upper/lower flow limit for gas source node j Sup Sup Pgas max / Pgas max upper/lower power limit for the sum of gas supply Sup Sup Sup Sup Pelec max / Pelec min, Qelec max /Qelec min upper/lower power limit for the sum of electricity supply Lori Lori Lori PijLori ,max / Qij,max , Pij,min / Qij,min maximum and minimum permissible active/reactive power flow through existent fix feeders Lup Lup Lup PijLup , k,max / Qij, k ,max , Pij, k ,min / Qij, k,min maximum and minimum permissible active/reactive power flow through replacing feeders Lnew Lnew Lnew PijLnew , k,max / Qij, k ,max , Pij, k ,min / Qij, k,min maximum and minimum permissible active/reactive power flow through adding feeders QkCB / QkSVC reactive power compensated per CB unit/SVC P jL,max / P jL,min maximum/minimum electrical load of node j ES PkHS ,max / Pk,max power flow limits of adding thermal storage units/ESS HS ES SkGS / S ,max k ,max / Sk,max storage capacity limits of adding gas storage units/thermal storage units/ESS Iij,max maximum allowable current amplitude Vmax /Vmin upper/lower limits of nodal voltage amplitude O _CB O _VR Nmax /Nmax maximum switching time for CB/SVC pv N jwt ,max /Nj,max maximum installation number of WTG/PVG NS number of scenarios NT number of time stages in a scenario Ny number of years during the whole planning horizon M very large number p2g chp gf eg / gh / gh sg / rg gas
Parameters duration time of time staget annual interest rate (s ) probability of scenarios lijpnew / lijlup/ lijlnew length of adding pipeline/replacing feeder/adding feeder rij /x ij resistance/reactance of feederij ckI _pipe/ ckI _c investment cost coefficients of pipelines/compressors ckI _gs/ckI _hs/ckI _ess investment cost coefficients of gas storage units/ thermal storage units/electric storage systems ckI _lup/ ckI _lnew investment cost coefficients of replacing/adding feeders ckI _wt /ckI _pv /ckI _cb/ ckI _svc / ckI _vr investment cost coefficients of WTG/ PVG/ESS/CB/SVC/VR c O _gs /c O _hs/ c O _ess operation cost coefficients of gas storage units/ thermal storage units/electric storage systems c O _p2g /c O _chp/ cO _gf operation cost coefficients of P2G/CHP/GF c O _loss operation cost coefficients of electrical loss c O _wt / c O _pv / operation cost coefficients of WTG/PVG /CB/SVC/ c O _cb/ c O _svc / cO _vr ESS/VR c Cur _wt / c Cur _pv / c Cur _L operation cost coefficients of curtailing wind power/solar power/electrical load ckM _pipe/ ckM _c maintenance cost coefficients of pipelines/compressors/gas storage units ckM _gs/ ckM _hs / ckM _ess Maintenance cost coefficients of gas storage
t
d
Constants fp
flow-power conversion coefficient
Variables gs c x ijpipe , k / x ij, k / x j, k binary investment variables for pipelines/compressors/gas storage units x jgs, k / x jhs, k / x jess , k binary investment variables for gas storage units/ thermal storage units/electric storage systems lnew x ijlup , k / x ij, k binary investment variables for replacing feeders/adding feeders/SVC/VR vr x jcb, k /x jsvc binary investment variables for CB/SVC/VR , k /x ij pv cb njwt , k / nj, k /nj, k integer investment variables for WTG/PVG/CB yjcb, t , s operation number of CB units yijvr, t , s operation tap position of VR rijvr, t , s operation ratios of VR
2
Applied Energy 253 (2019) 113517
J. Wang, et al.
nij+, tvr, s /nij, tvr, s , n+j, tcb, s / nj, tcb, s positive/negative regulation times of CB/SVC ujsg, t , s / urg j, t , s operation state of storage/release gas for gas storage units ujsh, t , s / urh j, t , s operation state of storage/release heat for thermal storage units ujse, t , s / ure j, t , s operation state of charging/discharging for ESSs gas nodal pressure j, t , s f jSup gas supply from source nodes ,t,s fij, t , s / f jk, t , s flows through pipelines Sup pjSup , t , s / qj , t , s active/reactive power supply from source node j
pjp, t2,gs
power flow into P2G from distribution network
pjLh ,t,s
thermal load power
WTG/PVG/load pj,Lt , s /qj,Lt , s active/reactive load after curtailment f jsg, t , s /f jrg, t , s gas flow into/out of gas storage units
rh pjsh , t , s / pj, t , s power flow into/out of thermal storage units pjse, t , s / pjre, t , s charging/discharging power of ESS
HS ESS S GS j, t , s / Sj, t , s /Sj, t , s storage capacity of gas storage units/ thermal storage units/ESS Iij, t , s current flow amplitude Vj, t , s nodal voltage amplitude Vm, t , s voltage amplitude of dummy bus in VR CTotal total cost CInv total cost of investment COpe /CMa/CSup annual cost of operation/maintenance/energy purchase CInv _Gas/CInv _Elec /CInv _Hub cost of investment for NGN/ADN/EH _Gas _Hub _Elec CsOpe /CsOpe /CsOpe cost of NGN/EH/ADN operation at time ,t ,t ,t stage t under scenario s CMa _Gas /CMa _Hub/CMa _Elec Annual NGN/EH/ADN maintenance cost _Gas _Elec CsSup /CsSup cost of gas/electricity purchase at time stage t ,t ,t under scenario s
gf f jchp , t , s / f j, t , s flow into CHP/GF from gas network pij, t , s / pjk, t , s , qij, t , s / qjk, t , s active and reactive power flows through feeders qjcb, t , s / qjsvc , t , s reactive power compensated by CB/SVC
f jL, t , s
gas load flow
pjL, t , s / qjL, t , s active/reactive load before curtailment pv wt pv pjwt , t , s /pj, t , s , qj, t , s / qj, t , s active and reactive power of WTG/PVG generation before curtailment _wt _pv _wt pjcur / pjcur , qjcur , t,s ,t,s ,t ,s / curtailment of active and reactive power for _pv cur _L _L qjcur /qjcur , t , s , pj, t ,t
management (DM), on load tap changers (OLTCs), and static var compensations (SVCs) was fully examined. Authors in [15] proposed a hierarchical robust expansion planning model integrating active management elements that includes OLTCs, SVCs, capacitor banks (CBs) and electric storage systems (ESSs). In [16], a flexible expansion planning model was proposed for ADN. In this work, DGs and ESSs were considered as capable expansion options, and in order to enhance the system’s reliability, island mode operation on the network was also incorporated in ANM schemes. In [17], the power electric vehicle and the solar parking lots can be applied as the sources of active and reactive power. Thus, a planning problem combing the economic aspects and technical factors of the system is proposed to achieve realistic results. Authors in [18] discussed the DG planning in ADN with the consideration of voltage regulators (VRs) and demand management. In [19], authors have developed an ADN planning model that focused on efficient renewable energy harvesting. In the ADN planning modeling of [20], authors mainly considered the dispatch of centralized and distributed ESSs. In [21], authors proposed an ADN planning with consideration of DM as well as the control of reactive compensators, OLTCs, ESSs and DGs. In the abovementioned literature, different parts of ANM schemes were modeled and examined selectively, mainly including adjustment of voltage regulators (VRs), reactive compensation of CBs and SVCs, dispatch of ESSs and DGs, demand management. However, there still remains a research gap in expansion planning model with full consideration of ANM schemes aforementioned. Another crucial challenge for planning is the method to address the uncertain source injection. The variability and uncertainty issues of renewable energy sources (RESs) have been well addressed in some research. In [22], a stochastic framework is adopted to solve the distribution system reconfiguration problem while adaptively and dynamically dealing with the incidental charging pattern of power electric vehicles and the variable and uncertain power of RESs. The reconfiguration of electrical distribution model in [23] is the first work to investigate the dynamic and adaptive reconfiguration problem, where the uncertainty of RESs at every time step (hour) of the day for the updated optimization time horizon is considered. In the same way, a model in [24] employed predictive control technique to minimize the energy loss with consideration of the uncertain RESs. Although these advanced dynamic techniques for RESs have been applied to model the reconfiguration problems in many studies, a long-term expansion planning issue requires a lager sample than a 24-h period data to cover
(MESEP) has attracted more attention than ever before. This issue includes the determination of the options for investment of renewable energy based distributed generation (DG) and the new infrastructures to be added to the existing system, and the size, location and type of them to be installed to meet the increasing energy demand over the planning horizon [6]. However, the expansion planning of MES is still confronted with a lot of huge challenges in planning and modeling approaches. Research in [7] presented a two-stage mixed-integer linear programming (MILP) model for MES planning to optimize the MES configuration from scratch considering distributed renewable energy source based on the EH concept. The proposed approach enabled planners to optimize the equipment investment and the MES configuration in a coordinated manner. In [8], simultaneous expansion planning of electricity and heat sources along with the transmission network was considered together. By using linear models for energy transmission network, the expansion problem was converted to a MILP problem to be solved efficiently. Authors in [9] proposed a new expansion planning formulation, which sought out the optimum usage of energy hubs for minimizing the utility costs, and at the same time maximizing the system’s reliability. These existing research works have made great contributions to planning the energy system based on the concept of energy hub; meanwhile, all the aforementioned works investigated the energy system expansion planning problems focusing on the part of the traditional power system. The electrical distribution planning problem has been well modeled in some studies like [10] and [11]. In [10], distributed generation placement planning problem has been modeled considering reliability level and power loss. In [11], the optimal placements of DGs and capacitors were investigated simultaneously in distribution network. However, due to renewable energy policy-making and its substantial advantages to accommodate the widespread renewable energy sources, the traditional distribution network is undergoing the transformation from a passive pattern to another, which is called active distribution network (ADN) [12]. The emerging active network management (ANM) technologies integrating into ADN make the network more flexible and controllable [13]. This transition brings huge challenges to distribution planners and operators. To model this issue, there are many research works concentrating on ADN expansion planning. In [14], a novel expansion planning model of ADN was proposed. In this work, the coordinated optimization of ANMs including the controls of DGs, demand 3
Applied Energy 253 (2019) 113517
J. Wang, et al.
more possible operating conditions. To the best of our knowledge, there still remains research gaps from the following aspects:
contains not only normal but also extreme conditions, which is more suitable for long term planning of larger systems. Due to the fact that mixed-integer nonlinear nature of the proposed expansion planning model makes the model particularly difficult to be solved, several equivalent reformulation methods are employed to linearize the model to be convex in order to ensure the global optimization results. Our linearization approach relies on the application of (i) second order cone relaxation (SOCR) method, which converts the original model into a mixed integer second order cone programming (MISOCP) model [27]; (ii) a modified piecewise linearization approach based on the work in [28], which can address the nonlinear polynomial expression efficiently; (iii) big M method, which can cope with the bilinear terms particularly; and (iv) auxiliary variable approach, which can replace the nonlinear product terms like quadratic terms and absolute value expressions. Finally, a modified IEEE 33-node distribution system and a new 23-node gas system coupled by 8 energy hub units are employed to verify the effectiveness of the proposed model. The main contributions of this paper can be summarized as follows:
(1) The comprehensive expansion planning model of multi-energy system where active distribution network is incorporated as an electricity part has never been reported in literatures. (2) The impact of ANM schemes on multi-energy system expansion planning has not been explored so far, and has not been applied to model the MES planning problem. To fill the knowledge gap of the state-of-the-art research, a new expansion planning model of multi-energy system is originally proposed based on the energy hub concept. In this model, the considered MES is composed of active distribution network, gas network, and energy hub unit. As a comprehensive expansion planning model, this paper centers on the optimal determination of the type, location and size of the feeders, DGs, CBs, SVCs, VRs and ESSs in ADN, natural gas pipelines, compressors and gas storage units (GSU) in NGN, and thermal storage units (TSUs) in EH. With the aim of minimizing the total cost over the planning horizon, the final investment decision scheme and the optimal operation strategy (ANM schemes, EH management, gas flow dispatch, etc.) for the expansion of MES are obtained. It is worth noting that, unlike the existing literatures, aforementioned ANMs including adjustment of VRs, reactive compensation of CBs and SVCs, dispatch of ESSs and DGs, as well as DM are fully modeled and tested in this paper. In addition, in order to consider the robustness [25] of the planning results, the stochastic framework is adopted and a robust probabilistic scenario generation approach extended from [26] based on the corresponding annual historical data is innovatively proposed to address the uncertainty of RESs, as well as the uncertain demands (electric, gas and heat demands). The resulting scenario set
(1) An expansion planning model of MES with the consideration of integrating ADN is originally proposed, in which the full ANM schemes are modeled and incorporated. (2) A novel scenario generation approach considering the robustness of planning results is designed for testing the proposed model. (3) A modified piecewise linearization as well as the SOCR approach is applied to convert the original nonlinear model into a mixed-integer SOCP model, and hence the model can be efficiently solved. The remainder of this paper is structured as follows: the MES is modeled in Section 2. In Section 3, a probabilistic scenario generation approach is proposed for testing the proposed model. In Section 4, the MES expansion planning problems integrated with ADN, NGN and EH
Substation Active Distribution Network
Natural Gas Network
Thermal Storage Unit
Gas Storage Unit
Gas Furnace
Combined Heat and Power
Power to Gas
Distribute Generation
Energy Hub
Gas Load
Electrical Energy Storage
Thermal Load
Energy Demand Fig. 1. Multi-energy system framework with energy hub. 4
Electrical Load
Applied Energy 253 (2019) 113517
J. Wang, et al.
power and solar power) throughout a year. Each of them is divided by the corresponding maximum value to obtain the per-unit value, and hence total 8760 hourly scenarios with normalized data is prepared. Step 2: Sort out the largest and smallest data for each hour through the year separately. Meanwhile, the mean value for each hour is computed. Therefore, three 24-h curves (i.e. the maximum, minimum and mean value curves) for each of them can be obtained. Step 3: Generate the scenario set based on three curves. The maximum value curves of wind power and solar power, as well as the minimum value curves of all demands are put together to constitute extreme scenario A . Meanwhile, extreme scenario B is generated by combining the maximum value curves of all demands and the minimum value curves of wind and solar power. Finally, scenario C is the combination of all mean values, expressed as vt , t = {1, 2, , 24} . Step 4: Calculate the probability for each scenario. Due to the data have been normalized in step 1, each value of the data is between 0 and 1. By introducing a fluctuation coefficient , data can be classified into ¯ t 1} , low level three intervals: the high level { ¯ t : vt + { : 0 vt } , and middle level { t : vt vt + } . In t t t detail, for a specific hour, if the demand is in the low level and energy power is in the high level, this hour will belong to scenario A . The reverse situation will belong to scenario B , and the rest belong to scenario C . Therefore, the probability of each scenario situation can be obtained by counting their frequency throughout the data (i.e. (s ) is calculated by dividing the total number of the data belonging to scenario s by the total number of the data in the year, s { A , B , C} and s (s ) = 1). In this paper, the proposed generation procedure of the extreme scenarios is implemented based on the historical data (download from http://www.elia.be/en/grid-data) of generation and demand for 2016. Accordingly, the maximum value curve and the minimum value curve of wind power, photovoltaic power, electricity demand, gas demand and heat demand can be obtained respectively, as shown in Figs. 2 and 3. In addition, the average values curves are shown and denoted by dotted line in each figure. For clarity, we give three points A, B and C at time t in the figures as examples. Each point includes five values (e.g. A = [A1, A2, A3, A4, A5]) representing the wind power, photovoltaic power, electricity, gas and heat demands in one operating condition (A) respectively. If the operating condition at time t is A, this scenario will be classified into scenario A since the generation power output A1 and A2 are in the high level while demands A3, A4 and A5 belong to lower level. In the same way, B and C will be classified into scenarios B and C , respectively.
units are modeled and formulated. The linearization methods, mainly including modified piecewise linearization and SOCR approach are described in Section 5. The flowchart of the overall proposed modeling methodology is given in Section 6. Numerical results and analysis based on a multi-energy system including a modified IEEE 33-node distribution system, a new 23-node gas system and 8 energy hub units are presented in Section 7. Finally, conclusions are drawn in Section 8. 2. Multi-energy system In this section, the framework of MES is modeled, as shown in Fig. 1. In MES, converters, storage and energy transfer equipment are used to combine and share the energy carriers under the so-called EH units. In this way, the MES is represented as an integrated system which is mainly composed of natural gas network, active distribution network and energy hubs that connect two networks. Furthermore, a EH unit is constituted by power to gas, combined heat and power, thermal storage unit and gas furnace (GF) to make its function possible. In the output port of EH, energy demand is provided by sharing and converting electricity and natural gas. In particular, the input port of EH is linked to energy networks such as ADN and NGN and the input electricity and gas are purchased from the upper distribution network or transformer station and natural gas distribution network, respectively; the output port of EH unit is connected to the loads or the energy networks. In this paper, thermal load is considered to be located in the output port of EH. In this regard, the conversions between natural gas and electricity as well as the supply of the thermal demand can be achieved and satisfied in the proposed approach. 3. Probabilistic scenario generation The key of probabilistic programming is to formulate the decisionmaking scheme under the combination of potential circumstances [29]. In order to obtain the robustness of investments, this paper tries to find the optimal planning scheme based on probabilistic scenarios which includes both extreme and typical situations. Given the fact that wind speed, solar radiation intensity and load demands all change constantly with time, it is likely to generate the scenarios based on the historical hourly data. In this paper, 24-h data through a year are used as the research objects of scenarios, and thus main characteristics of time sequence of wind power, photovoltaic power generation and load can also be completely reflected. Based on this point, the probabilistic scenario set is generated by the following steps. Step 1: Firstly, normalize the historical data of demands (including electricity, gas and heat demands) and renewable energy (i.e., wind
Fig. 2. Renewable energy power curves of 24-h under 3 scenarios (including the maximum value and the minimum value curves, as well as the mean value curves of the wind and photovoltaic power). 5
Applied Energy 253 (2019) 113517
1 B3 C3 A3
Gap2
0.5
Min value curve Max value curve Mean value curve
0
6
1
Gas flow (pu)
Gap1
v3 t +Δ v3t v3 t - Δ
B4 C4
A4
0.5
0
12
Time (h) (a) Electricity demand
6
12
18 Thermal power (pu)
Demand power (pu)
J. Wang, et al.
v4 t + Δ
v4t v4 t - Δ
18
1
B5 v + Δ 5t
C5 v5t A5 v5 t - Δ
0.5
0
24
Time (h) (b) Natual gas demand
24
6
12
Time (h) (c) Heat demand
18
24
Fig. 3. Demand power curves of 24-h under 3 scenarios (including the maximum value, the minimum value and the mean value curves for electrical, gas and thermal demands).
4. Modeling of the multi-energy system expansion planning
leverage the investment and operation costs.
Based on the framework of MES and the approach of probabilistic scenario generation, a model of expansion planning for MES can be mathematically formulated using the approach of scenario-based deterministic equivalent, as described as follows. To increase the readability of this passage, some brief descriptions about key symbols which are easy to be confused or used repeatedly are given here: t refers to the time stage; s is the index of scenario; i and ij denote the identifier of node and pipe/feeder respectively; k is the index of type of investment equipment. Omitting the superscripts and subscripts of the notations for simplification, c represents the unit cost of investment or operation; Γ, Ω and Λ are the collections of the existing, candidate and all locations of the investment equipment respectively, and = ; K is the collection of alternative types; x denotes the corresponding binary investment variable; n is the integer investment variable; y is the integer operating variable of infrastructures; parameter l denotes the length of pipelines or feeders; f / F denotes the variable/parameter related to the gas flow; represents the nodal pressure in the gas network; p / P and q/Q are the variables/parameters representing the corresponding active and reactive power. N denotes the number of corresponding things (such as years and facilities). In addition, other notations in this paper are fully explained in Nomenclature part.
4.1.1. Investment cost In this paper, the expected investment cost is firstly formulated as: It can be seen that the investments in the proposed model are constituted by three items: (i) the cost in gas network with consideration of installing pipelines, compressors and gas storage units (Eq. (3)); (ii) the cost in energy hubs where the installation of thermal storage unit is involved (Eq. (4)); and (iii) the cost in active distribution network considering the upgrade and expansion of feeders, as well as the installations of wind turbine generations (WTGs), photovoltaic generations (PVGs), ESSs, CBs, SVCs and VRs (Eq. (5)).
CInv _Gas =
=1
Pipe k KPipe
j
GS k KGS
ckI _pipe lijpnew x ijpipe ,k +
(COpe + CMa + CSup) 1
ckI _c x ijc, k + (3)
Lup k KLup
ckI _lup lijlup xijlup ,k +
j
+ PVG k KPVG
WTG k KWTG
ckI _pv njpv ,k +
+ SVC k K SVC
ij
Lnew k KLnew
ckI _lnew lijlnew x ijlnew ,k
ckI _wt njwt ,k
+
j
(4)
HS k KHS
ij
1 (1 + d )
C k KC
ckI _gs x jgs, k
CInv _Elec =
j
ij
ckI _hs x jhs, k j
The objective (1) is to minimize the expected total cost in the planning horizon, which consists of 4 terms: (i) investment cost (Eqs. (2)–(5)), (ii) operation cost (Eqs. (6)–(9)), (iii) maintenance cost (Eqs. (10)–(13)), and (iv) energy purchase cost (Eqs. (14)–(16)) in the multienergy system. In other words, based on the investment and construction, the corresponding operation cost, maintenance cost and energy purchase cost are taken into account. Ny
ij
CInv _Hub =
4.1. Objective function
min CTotal = CInv +
(2)
CInv = CInv _Gas + CInv _Hub + CInv _Elec
j
ESS k KESS
ckI _ess x jess ,k +
j
CB k KCB
ckI _cb njcb, k
ckI _svc x jsvc ,k + ij
VR k KVR
ckI _vr x jvr, k (5)
Eqs. (3)–(5) is written in terms of the variables that represent the investment decisions. For example, the cost of the newly-added pipelines depends on the binary investment variable of newly-added pipe ij I _pipe ( x ijpipe ) and the length , k ), the cost of the k -type pipe per unit length (ck pnew of each newly-added pipeline (lij ); the investment cost of
(1)
where Ny represents number of years during the whole planning horizon and d represents annual discount rate. Herein, the economic factor (i.e., discount rate d [30]) is considered in the planning problem to 6
Applied Energy 253 (2019) 113517
J. Wang, et al.
compressors is calculated according to the binary variable of whether to use k -type compressor ( x ijc, k ) and its price ckI _c ; the investment cost of DG pv is calculated by the number of the k -type WTG/PVG (njwt , k /nj, k ) in-
Eq. (13) represents the additional maintenance costs brought by upgraded lines and the maintenance costs corresponding to newly-added lines and newly-installed WTGs, PVGs, ESSs, CBs, SVCs and VRs in active distribution network.
stallation and its price (ckI _wt /ckI _pv ). The other investment costs are similar.
CMa _Gas =
4.1.2. Operation cost The operation cost is given as: NT
COpe = 365
_Gas _Hub _Elec (CsOpe + CsOpe + CsOpe ) ,t ,t ,t
(s ) s
(6)
t=1
c O _gs (f jsg, t , s + f jrg, t , s )
j
fp
(cO _p2gpjp, t2,gs + c O _chpf jchp ,t ,s
j
fp
+ cO _gf f jgf, t , s
ckM _hs x jhs, k
ij
_Elec CsOpe = ,t
ij
j
c O _vr |yijvr, t , s
WTG k KWTG
j
+
ij
1, s |
c O _cb |yjcb, t , s
+ j
yjcb, t
CB k K CB
ESS k KESS
ckM _cb njcb, k +
c O _svcqjsvc ,t,s +
j
NT
CSup = 365
(s ) s
c O _ess (pjse, t , s + pjre, t , s ) +
j
j
PVG
j
L
_L c Cur _Lpjcur + ,t, s
+
j
c O _wtpjwt ,t,s + WTG
j
1, s |
_wt c Cur _wtpjcur , t,s WTG
cO _pvpjpv ,t ,s PVG
(9) 4.1.3. Maintenance cost
CMa = CMa _Gas + CMa _Hub + CMa _Elec
PVG k KPVG
ckM _pv njpv ,k
ckM _svc x jsvc ,k +
ij
VR k KVR
ckM _vr x ijvr, k
_Gas _Elec (CsSup + CsSup ) ,t ,t
(14)
t,
s
(15)
t,
s
(16)
where the annual energy purchase cost (CSup ) is calculated by the similar way to the annual operation cost (COpe ) in terms of the method of _gas _elec probabilistic scenario generation. ptsup and ptsup indicate the total ,s ,s gas and electricity supply power at time stage t of scenario s, (i.e., sup _elec _gas = j E _Sup pjSup ptsup = j G _Sup f jSup , t , s respectively). ,s , t , s fp and pt , s In addition, the real time pricing options will gradually become a trend in energy trade in the future. In order to optimize energy utilizing from the perspective of the MES, the real-time change of the supply elec prices of gas and electricity (ctgas , s and ct , s ) are considered in the model, as done in [31]. That is, the MES is allowed to access to the real time price data and then it can wisely manage its energy consumption according to them. By this means, the energy converters (i.e., CHP and P2G) of EHs is controlled and operated according to weighing the cost of using different forms of energy. If the gas price is low while electricity price is in the high level, CHP will generate electricity using natural gas to support the supply of electricity. On the contrary, if the gas price is high while electricity price is in the low level, P2G will work and convert the excess electrical energy into the natural gas form. It is assumed that the gas and electricity price are determined by the _gas _elec demands ( ptsup and ptsup ) in real time. ,s ,s
_pv c Cur _pvpjcur ,t, s
+
j
ckM _ess x jes, k
SVC k K SVC
t=1
_Gas sup _gas CsSup = ctgas t, ,t , s pt , s
CB
ESS
j
ckM _lnew lijlnew x ijlnew ,k
4.1.4. Energy purchase cost The energy purchase cost is formulated as Eq. (14) including gas purchase and electricity purchase and the energy costs associated with them are formulated as (15) and (16) respectively.
Line
SVC
ckM _wt njwt ,k +
c M _lori ) lijlup x ijlup ,k
(13)
c O _lossIij2, t , s rij
+ t +
yijvr, t
(ckM _lup
Lnew k KLnew
+
(8)
VR
(12)
Lup k KLup
+
fp ) +
c O _hs (pjsh, t , s + pjrh, t , s )
ckM _c x ijc, k +
HS k KHS
ij
j
C k KC
(11)
+
Hub
Hs
ij
ckM _gs x jgs, k
_Elec sup _elec CsSup = ctelec t, ,t , s pt , s
= t
j
GS k KGS
ckM _pipe lijpnew x ijpipe + ,k
CMa _Elec =
_Hub CsOpe ,t
j
j
j
(7)
GS
Pipe k KPipe
CMa _Hub =
It is worth to be noted that the operational cost in (6) is computed based on the approach of probabilistic scenario generation. By this means, there are 365· (s ) days belonging to scenario s and NT time stages in each scenario. Thus, the annual operation cost (COpe ) is calculated by the sum of the total operation cost at each time stage t under Ope each scenario s (i.e., CsOpe , t ) through the year and Cs, t is equal to the sum of operation costs in the gas network, distribution network and energy Ope _Gas _Hub _Elec + CsOpe + CsOpe ). hubs (i.e., CsOpe , t = Cs, t ,t ,t Ope _Gas C In detail, s, t is characterized in (7) where the operation of gas storage units is included; Eq. (8) models the operation costs related to _Hub energy hub (CsOpe ) consisting of P2Gs, CHPs, GFs and thermal sto,t rage units; In (9), not only power losses and operations of VRs, CBs, SVCs, EESs, WTGs and PVGs but also active management measures such as curtailments of wind power, solar power and electrical load are taken into account. In general, the operational cost for each component _Gas _Elec _Hub (CsOpe /CsOpe /CsOpe ) can be calculated by the production of the ,t ,t ,t corresponding power (such as f jsg, t , s fp / pjp, t2,gs ), the unit cost (such as c O _gs /c O _p2g ) and the duration ( t ). Considering the artificial adjustment expenses (for VRs or CBs), the regulation times in the duration (such as |yijvr, t , s yijvr, t 1, s |), multiplied by the relevant unit cost (such as c O _vr ) is included in the operation cost of ADN. _Gas CsOpe = t ,t
ij
(1) Real-time gas price
(10)
In a similar way to the formulation of the investment cost, the maintenance cost is formulated as Eq. (11)–(13). The maintenance costs of newly-added pipelines, compressors and gas storage units in the gas network are modeled in (11); in (12), the maintenance costs of newlyinstalled thermal storage units are taken into account in energy hubs;
It is assumed in (17) that the relationship between supply price and demand of natural gas can be described by a monotonic incremental linear function.
ctgas ,s = 7
sup _gas 1 pt , s
+
0
(17)
Applied Energy 253 (2019) 113517
J. Wang, et al.
(2) Real-time electricity price
x are established in (21)–(30). This set of constraints indicates that a maximum of one replacement, addition, or installation is permitted for each component and location as well. Expressly, the connectivity and the radiation of the expanded integrated network are ensured by constraints (21) and (22). Particularly for a newly-added bus j ( j E _B , and E _B denotes the collection of the newly-added buses), it is restricted to connect to one and only one existent bus (the index i is to represent the candidate existent buses). Meanwhile, a maximum pv number of installation for WTG/PVG (N jwt ,max /Nj,max ) at each candidate location is imposed in (31)–(32). (33) limits the maximum permitted installation number of capacitors in each CB unit (N CB j,max ).
The electricity cost of upper level substation can be expressed as follows: _elec ce (ptsup )= ,s
sup _elec m ) m (pt , s
_elec + ...+ 1 ptsup + ,s
(18a)
0
In (18b), it is supposed that the upper-level substation sets the supply price according to the maximization of electricity sale income. NT
max s
t=1
_elec elec (ptsup ct , s ,s
_elec ce (ptsup )) ,s
(18b)
The solution of the optimization problem can be obtained as follows:
4.3. Operational constraints
sup _elec ctelec ) , s = ce (pt , s
4.3.1. Constraints related to the natural gas network
(18c)
To simplify the problem, set m = 2: _elec ce (ptsup ) ,s
ctelec ,s
=2
=
sup _elec 2 ) 2 (pt , s
sup _elec 2 pt , s
+
sup _elec 1 pt , s
+
(1) Gas supply constraints
+
(18d)
0
0
(18e)
1
sup _gas 0 ) pt , s
_Elec _elec CsSup = (2 2 ptsup + ,t ,s
t,
sup _elec 1) pt , s
t, t,
s t,
j,min
(20)
where and
j
G _B ,
lnew x ijlnew , k = 1, x ij, k = \{ 0, 1\} , i
j
E _B , i
k KLnew
x ijc, k
1, x ijc, k = \{ 0, 1\} ,
ij
ij
Pipe
x jgs, k
1,
x jgs, k
= \{ 0, 1\} ,
j
1, x jhs, k = \{ 0, 1\} ,
GS
j
1, x ijlup , k = \{ 0, 1\} ,
HS
x jess ,k
1, x jess , k = \{ 0, 1\} ,
x jcb, k
1, x jcb, k = \{ 0, 1\} ,
ij
ESS
CB
j
SVC
1, xijvr = \{ 0, 1\} ,ij njwt ,k
wt N jwt ,max , nj, k
Z+,
j
WTG
njpv ,k
pv Njpv ,max , nj, k
Z+,
j
PVG
j
CB ,
(35)
s
ij IJ (j)
f jL, t , s + f jSup + f jrg, t , s ,t ,s
fij, t , s
f jsg, t , s
f jgf, t , s
f jchp , t,s
lg Lg j, t , s P j
j t,
s
K CB
sgn(
2 i, t , s
2 2 j, t , s ) fij, t , s
= =
ij (
2 i, t , s
2 j, t , s ),
ij (
2 i, t , s
2 j, t , s )
ij
Pipe /
C,
x ijpipe ,k ,
ij
t,
s Pipe ,
(37)
t,
s
ij (
2 i, t , s
2 j, t , s )
ij (
2 i, t , s
2 j, t , s )
t,
s
sgn(
2 i, t , s
2 2 j, t , s ) fij, t , s
+M
sgn(
2 i, t , s
2 2 j, t , s ) fij, t , s
M
k KC k KC
x ijc, k xijc, k
,
ij (39)
by the maximum gas demand of node j (P jLg ) and the gas load rate at
time t of scenario s ( j, tlg, s ) and the rate can be obtained by the approach of probabilistic scenario generation. Additionally, the Weymouth steady state power flow constraints [32] are formulated as expression (37)–(39) where ij represents the transmission coefficient through pipeline ij and sgn( i2, t , s - j2, t , s ) is a sign function (for sgn(a) , if a > 0 , the value of sgn(a) is 1; if a = 0 , the value of sgn(a) is 0 and otherwise, the value of sgn(a) is −1). It should be noted that the flow through pipeline ij ( fij, t , s ) is imposed
(32)
k
sgn(
2 2 j, t , s ) fij, t , s
Gas flow balance is formulated by Eq. (36), which indicates that the sum of gas inflows is equal to the sum of gas outflows in any node at any time stage and any scenario. Note that the gas load ( f jL, t , s ) is computed
(31)
k KPV
Z +,
f jk, t , s =
2 i, t , s
C,
(30)
VR
k KWT
cb cb N CB j,max x j, k , nj, k
t,
(38)
(29)
k K SVC
x ijvr
G _B ,
k KPipe
(28)
k KCB
1, x jsvc , k = \{ 0, 1\} ,
j
represents nodal pressure at time stage t of scenario s; are the lower and upper pressure limits.
G _B ,
(27)
k KESS
x jsvc ,k
j, t , s
j,max
(22)
(26)
k KLup
j
j,max ,
f jL, t , s =
(25) Lup
j
(34)
+ f jp, t2,gs
(21)
(24)
k KHS
xijlup ,k
Lnew
jk JK (j)
(23)
k KGS
x jhs, k
ij
s
(36)
C
k KC
njcb, k
i
k KPipe
t,
(3) Gas flow constraints
Constraints related to the investment and utilization decisions of all facilities in this planning are given as: i
G _Sup ,
denotes the gas supplied by the gas source node j at time
j, t , s
j,min
4.2. Investment and utilization constraints
pipe x ijpipe , k = 1, x ij, k = \{ 0, 1\} ,
where
f jSup ,t,s
j
(2) Pressure inequality constraints
(19)
s
F jSup ,max ,
stage t of scenario s and F jSup ,max refers to the corresponding maximum value.
Therefore, Eqs. (15) and (16) can be transformed into (19) and (20), respectively. _Gas _gas CsSup = ( 1 ptsup + ,t ,s
f jSup ,t ,s
to be 0 if no pipeline is added (i.e. k KPipe x ijpipe , k = 0 ) in (38). Constraint (39) guarantees that compressors can be added in pipelines that have been previously built. Given that M is a very large number, (39) will be
(33)
The relevant constraints in terms of the binary investment decision 8
Applied Energy 253 (2019) 113517
J. Wang, et al.
equivalent to (37) if not to add the compressor ( restricted otherwise.
k KC
for each gas storage unit, its initial capacity (S GS j,0, s ) is imposed to be 0.
x ijc, k = 0 ), un-
4.3.2. Constraints relevant to energy hubs
(4) Pipeline operation constraints
FijPipe ,max
FijPipe ,max ,
fij, t , s
pipe FijPipe , k ,max x ij, k
ij
Pipe /
C,
t,
s
pipe FijPipe , k,max x ij, k ,
fij, t , s
k KPipe
(1) Thermal power balance constraints
(40)
ij
Pipe ,
t,
sh phgf j, t , s + phchpj, t , s = pjLh , t , s + pj, t , s
s
pjLh ,t ,s =
k KPipe
(41)
FijPipe ,max
x ijc, k
M
FijPipe ,max + M
fij, t , s
k KC
x ijc, k ,
ij
C,
t,
k KC
j, t , s
FijC,max for
s:
c ij i, t , s
=
(43)
FijC,max
fij, t , s
ij
C,
c ij, k i, t , s
k
t,
x ijc, k )
M (1
FijC, k,max
(44)
KC,
xijc, k )
M (1
c ij, k i, t , s
x ijc, k )
+ M (1
FijC, k,max + M (1
fij, t , s
x ijc, k )
pjgf, t , s = f jgf, t , s
(45) (46)
GS ,
ujsg, t , s + urg j, t , s
f jsg, t , s
0
f jrg, t , s
S GS j, t , s =
0
t
gs rg FkGS ,max x j, k uj, t , s ,
t
(48b)
k K GS
0+
S GS j, t 1, s
+
t
f jsg, t , s sg
t
gs SkGS ,max x j, k ,
S GS j, t , s k
GS
t
f jrg, t , s / rg
t
t=1 t
j
HUB ,
t,
s
(51)
pechpj, t , s = f jchp ,t ,s
chp fp ge ,
f jchp , t,s
j
HUB ,
t,
s
chp Fmax
(52)
t
2
pjp, t2,gs
p2g eg / fp p2g Pmax
,
j
HUB ,
t,
s
(53)
Note that in (53) the gas output of P2G ( f jp, t2,gs ) is calculated in terms of the electric power injecting into P2G from distribution network ( pjp, t2,gs ), multiplied by the corresponding efficiency ( egp2g ) and divided by the flow-power conversion coefficient ( fp ). In (53) P2G is also imposed to transform from electric power into gas flow by the corresponding p2g non-negative constraint and Pmax is the upper power input limit for P2G.
(48a)
k K GS
f jrg, t , s / rg
,
chp fp gh
0
(47)
f jsg, t , s sg
gf Fmax
f jp, t2,gs = pjp, t2,gs
t
gs sg FkGS ,max x j, k uj, t , s ,
(50)
(4) Operational constraints for power to gas
ujrg,1, s = 0 0
gf fp gh
phchpj, t , s = f jchp ,t ,s 0
s:
1,
s
(3) Operational constraints for combined heat and power
(6) Gas storage unit operation constraints
j
f jgf, t , s
0
The operation of existing compressors ( ij C ) is modeled by (43)–(44) where ijc / ijc, k represents the booster coefficient of existent and newly-added compressors. On the other hand, the operation of those potential to be installed ( ij C ) are subject to (45) and (46), which are still associated with the investment decision on them ( x ijc, k ).
for
t,
(2) Operational constraints for gas furnace
s:
j, t , s
HUB ,
gf chp / gh ) and the flow-power ergy conversion efficiency of GF/CHP ( gh conversion coefficient fp . Constraints (51) and (52) also give the fact that the gas flow over a given range can input into GF and CHP. In addition, P2G is utilized in the energy hubs to convert electric power into natural gas.
(5) Compressor operation constraints
t,
j
computed according to the maximum load (P jLh ) and the corresponding
c also on the investment decision x ijpipe , k / x ij, k .
C,
,
rate ( j, tlg, s ). It also can be noted that the construction of thermal storage unit provides the possibility of storing heat that could not be consumed or generated in time. Based on the model of EH, the thermal power could be supplied by GF and CHP. In (51)/(52), the thermal power supplied by GF/CHP is calculated according to the gas flow injecting into GF/CHP from gas network ( f jgf, t , s / f jchp , t , s ), multiplied by relevant en-
Expression (40) sets the bounds on the gas flow in pipelines which has been built previously. In (41) and (42), it can be seen that the gas flow through the newly-added pipelines or candidate pipelines for Pipe compressors depends not only on the upper flow limit FijPipe , k ,max /Fij,max but
ij
lh Lh j, t , s P j
For energy hub nodes, the sum of thermal power outputs and the sum of thermal power inputs is equal, which is stated in Eq. (50). In the way that is similar to computing the gas load, the thermal load ( pjLh , t , s ) is
s (42)
for
pjrh, t , s
(5) Operation constraints for thermal storage units
(49a)
Exactly as gas storage units, thermal storage unit constraints can be modeled as follows:
(49b)
for
It is not allowed to storage and release gas simultaneously for each gas storage unit and it is also not allowed to release gas at the first stage (t = 1), which can be enforced by (47). Constraint (48) represents that each gas storage unit storages or releases gas only when the conditions sg and rg in conx gs = 1 and or urg j, t , s = 1 are both satisfied. k K GS j, k straint (49) denote the efficiencies of gas storage and release respectively. S GS j, t , s represents the capacity at time stage t of scenario s and it is associated with the capacity in the former time stage (S GS j, t 1, s ), the flow of gas storage/release in the current time stage ( f jsg, t , s / f jrg, t , s ), storage/ release efficiency of gas storage units ( sg / rg ) and the duration t . And
j
HS ,
ujsh, t , s + urh j, t , s ujrh,0, s 0
s:
1,
t
=0
pjsh ,t, s
(54) hs sh PkHS ,max x j, k uj, t , s ,
t
hs rh PkHS ,max x j, k uj, t , s ,
t
k KHS
0
pjrh, t , s k KHS
9
(55a) (55b)
Applied Energy 253 (2019) 113517
J. Wang, et al.
0 + pjsh, t , s
SjHS ,t ,s =
SjHS , t 1, s
sh
+
SjHS ,t, s
0
pjrh, t , s /
t
sh pjsh ,t, s
t
(4) Feeder operation constraints
t=1
pjrh, t , s / rh
t
hs SkHS ,max x j, k , k
rh
t
t
2
for
(56a)
t
(56b)
HS
for
4.3.3. Constraints related to the active distribution network
j
E _Sup ,
t,
Lori ,
t,
s:
PijLori ,max
pij, t , s
PijLori ,max
QijLori ,max
qij, t , s
QijLori ,max
ij
Lup ,
PijLori ,max (1
(1) Electricity supply constraints
for
ij
t,
k KLup
s: x ijlup ,k )
PijLori ,max (1
pjSup ,t, s
P jSup ,max
(57)
0
qjSup ,t, s
QjSup ,max
(58)
k KLup
QijLori ,max (1
k KLup
x ijlup ,k )
k KLup
for
(2) Security constraints
0
Iij, t , s
Iij,max ,
j
E _B ,
ij
t,
Line ,
k KLnew
(60)
s
k KLnew
where Vmin and Vmax are the lower and upper voltage level limits. Vj, t , s represents the nodal voltage amplitude and Iij,max denotes the current amplitude for feeder ij.
jk JK (j)
ij IJ (j)
pij, t , s
ij IJ (j)
pv + pj,wt t , s + pj, t , s
+ pjre, t , s jk JK (j)
qjk, t , s =
pjse, t , s + pechpj, t , s
ij IJ (j)
qij, t , s
ij IJ (j)
pjp, t2,gs qj,Lt , s + qjSup ,t, s
Iij2, t , s xij
pij2, t , s + qij2, t , s = Vi2, t , s Iij2, t , s , V j2, t , s
=
Vi2, t , s t,
V j2, t , s
ij
Line ,
2(pij, t , s rij + qij, t , s x ij ) +
(61)
s
t, Iij2, t , s (rij2
s
for
(62)
+
x ij2 ),
ij
Line /
Vi2, t , s
k KLnew
2(pij, t , s rij + qij, t , s x ij ) +
M (1 Lnew ,
k KLnew
t,
s
j
Lnew ,
0
lnew QijLnew , k ,max x ij, k
qij, t , s
k KLnew k KLnew
lnew PijLnew , k ,max x ij, k lnew QijLnew , k ,max x ij, k
ESS ,
(67)
s:
1,
t (68)
pjse, t , s
ess se PkESS ,max x j, k uj, t , s ,
t
ess re PkESS ,max x j, k uj, t , s ,
t
(69a)
k KES
0
x ijlnew ,k )
x ijlnew ,k
pij, t , s
ujre,0, s = 0
2(pij, t , s rij + qij, t , s x ij ) + Iij2, t , s (rij2 + x ij2)
Vi2, t , s
s:
lnew PijLnew , k ,max x ij, k
ujse, t , s + ure j, t , s
(63)
s
+ M (1 V j2, t , s
t,
(66)
In a similar way to GSU and TSU modeling, the EES can be modeled as (68)–(70).
+ qjcb, t , s + qjsvc ,t,s E _B ,
lup QijLup , k,max x ij, k
(5) Operation constraints for electrical energy storage units
pv + qj,wt t , s + qj, t , s
j
x ijlup ,k )
feeder ij is non-upgraded ( k KLup xijlup , k = 0 ), the boundary of power flow will be forced to keep the same with the original one; otherwise, it will be changed into the upgraded one. Analogously, the value of pij, t , s /qij, t , s for candidate newly-added feeder ij depends on whether feeder ij is added. If added, pij, t , s /qij, t , s is limited to the ranges of the new one, being 0 otherwise.
pj,Lt , s + pjSup ,t, s
Iij2, t , s rij
t,
k KLup
qij, t , s
The power limits of original fixed feeders ( ij Lori ) are formulated in (65) whereas that of replacing and adding ones ( ij Lup and ij Lnew , respectively) are stated in (66) and (67). It can be seen that x ijlup , k plays a critical role in restricting the boundary of pij, t , s /qij, t , s . If
(3) Power flow balance constraints
pjk, t , s =
Lnew ,
(59)
s
t,
ij
x ijlup ,k )
lup PijLup , k,max x ij, k
QijLori ,max (1 +
pij, t , s
lup QijLup , k ,max x ij, k
k KLup
Sup former station at time stage t of scenario s and P jSup ,max / Q j,max denotes the corresponding upper active/reactive power limit.
Vmax ,
k KLup
+
Sup where pjSup , t , s /qj , t , s is the active/reactive power supplied by the trans-
Vj, t , s
lup PijLup , k,max x ij, k
k KLup
s:
0
Vmin
(65)
pjre, t , s
(69b)
k KES
Iij2, t , s (rij2
+
x ij2)
,
ij SjESS ,t ,s =
) (64)
0
0+
pjse, t , s se
SjESS , t 1, s
SjESS ,t, s
t
pjse, t , s se
t
ess SkESS ,max x j, k , k
The balance of the power flow for the ADN can be ensured by (61)–(64). pv wt L q pv It should be noted that pj,Lt , s , pj,wt t , s , pj, t , s and qj, t , s , qj, t , s , j, t , s in (61) represent active and reactive power of electrical load, WTG generation, PVG generation after curtailment, respectively. And constraint (64) indicates that if k KLnew xijlnew , k = 0 , which means that candidate feeder ij will not be added into the network, there will be no voltage connection between nodes i and j .
+
pjre, t , s / re
t
pjre, t , s / re
t=1 t
t
t
ES
2
(70a) (70b)
(6) Operation constraints for capacitor banks
for 0
j
CB ,
yjcb, t , s
njcb, k , k K CB
10
s: t
(71)
Applied Energy 253 (2019) 113517
J. Wang, et al.
qjcb, t , s =
QkCB yjcb, t , s x jcb, k ,
t
NT s t=1
|yjcb, t , s
yjcb, t
(9) Operation constraints for distributed generations
(72)
k K CB
DG DG and min represent DG penetration is restricted in (80), where max the maximum and minimum permissible of DG penetration rates.
O _CB Nmax
1, s |
yjcb,0, s = 0
DG min
(73)
pjwt ,t, s =
j
SVC ,
t,
s
k K SVC
pjpv ,t ,s =
(74)
Vm2, t , s = Vi2, t , s
t,
Vmax / kijvr,min
Vm, t , s
V j, t , s
Vm, t , s
Mx ijvr
(77)
kijvr, t , s = 1 + ( kijvr % × yijvr, t , s )
(79)
Rijvr
yijvr, t , s
Rijvr , yijvr, t , s
Z
ij NT
s t=1
|yijvr, t , s
1, s |
yijvr,0, s = 0
t,
s
(83)
pv PkPV ,max nj, k
j
PVG ,
t,
s
(84)
p pv ( pjwt , t , s / j, t , s )
is written in
Apart from power curtailment of DG, load-shedding is also considered in this paper as one of the important ANM schemes. The relevant constraint is given in expressions (85) where the definition of L L Cur _L L L j,max , pj, t , s /qj, t , s and pj, t , s /qj, t , s is similar to (83) and (84). Additionally,
(80)
the active and reactive loads ( pjL, t , s and qjL, t , s ) can be obtained in a similar way to the aforementioned gas/thermal load ( f jL, t , s / pjLh , t , s ).
0
pjL, t , s =
le Le j, t , s P j
qjL, t , s
le Le j, t , s Q j
=
_L pjcur ,t,s
pj,Lt , s
=
pjL, t , s
Cur _L L j,max pj, t , s
j
L,
t,
s
_L pjcur ,t,s
qj,Lt , s = qjL, t , s × pj,Lt , s / pjL, t , s
Ui
VR :
yijvr, t
k KPVG
WTG ,
(10) Operation constraints relevant to load shedding
The voltage of the virtual node Vm, t , s can be computed by (75). The upper and lower voltages of virtual node m are limited by (76) to ensure that the voltage of node j will not exceed the limits (Vmin and Vmax ). In (77), the value of Vj, t , s is associated with the binary variable x ijvr and the value ofVm, t , s , where for x ijvr = 0 (i.e., the VR is not to be installed), the values of Vj, t , s and Vm, t , s will be imposed to be equal. And for x ijvr = 1 (i.e., the VR is to be installed), Vj, t , s will be equal to kijvr, t Vm, t , s by (77) and (78). Moreover, the operational VR ratio kijvr, t , s is reflected in (79) and (80), where an integer variable, namely yijvr, t , s denotes the tap position of VR at time stage t of scenario s. It also can be shown in (80) that the tap position variable can only vary from Rijvr to Rijvr , corresponding to the VR’s minimum and the maximum ratios. In addition, analogously to constraint (73), changing times of VR tap are restricted by (81) and meanwhile, the initial value of yijvr, t , s is imposed to be 0.
for
pv j, t , s
j
wt pv pv refers to the power factor of DG; pj,wt t , s / pj, t , s and qj, t , s /qj, t , s represent the active and reactive power of DG after curtailment. It also can be obtained that the certain triangular region in Fig. 5 is actually corresponding to the DG operation control according to (83) or (84), which indicates that the controllable range of reactive power in DG relies on the active power output after curtailment.
(76)
(78)
wt PkWT ,max nj, k
The original DG (WTG/PVG) generation
(75)
Vj, t , s = kijvr, t , s Vm, t , s
(82)
PV terms of the upper output limit for DG using type k (PkWT ,max /Pk ,max ), wt multiplied by the installation number at node j (nj, k ) and the output coefficient at time stage t of scenario s ( j, t ,wt s ). In addition, it is assumed that DGs can compensate or absorb reactive power to a certain extent, as shown in Fig. 5. Thus, the curtailment of wind and solar power is subject to constraints (83) and (84) respectively. In this set of conCur _pv _wt wt straints, jCur / pv ,max / j,max denotes the maximum curtailment rate;
s:
2(pij, t , s rij + qij, t , s xij ) + Iij2, t , s (rij2 + x ij2 )
Vmin / kijvr,max
Mx ijvr
VR ,
P jL,max
_pv Cur _pv pv pjcur ,t ,s j,max pj, t , s pv pv _pv pj, t , s = pj, t , s pjcur ,t,s pv pv qj, t , s max = pj, t , s tan(cos 1 pv ) qj,pv qj,pv qj,pv t , s max t ,s t , s max
For the purpose of improving voltage quality in the network, the investment and utilization of VR is considered in this planning as well. It is worth emphasizing that, the VR is modeled as a branch connecting buses i and j with an ideal and discretely adjustable transformer, as shown in Fig. 4. On this basis, the modeling of VR can be described by expressions (75)–(80).
m (ij), ij
pv PkPV ,max nj, k
0
(8) Operation constraints for voltage regulators
m
k KWTG
PVG k KPVG
L
_wt Cur _wt wt pjcur ,t, s j,max pj, t , s wt wt _wt pj, t , s = pj, t , s pjcur ,t,s wt 1 wt ) qj,wt t , s max = pj, t , s tan(cos qj,wt qj,wt qj,wt t , s max t ,s t , s max
where the maximum reactive power imported by SVC at node j (qjsvc , t , s ) is associated with the type of SVC and the investment decision ( x jsvc , k ).
for
wt j, t , s
j
0
Unlike CB, SVC can be continuously adjusted in a certain range, which can be stated as:
0
WTG k KWTG
j
(7) Operation constraints for static var compensations
QkSVC x jsvc ,k ,
wt PkWT ,max nj, k + j
L
DG max
where yjcb, t , s in (71) and qjcb, t , s in (72) represent the number and the reactive power output of the CB units in operation at time stage t of scenario s respectively. In order to extend life expectancy of CB switches, switching times are constrained by (73) during the planning period.
qjsvc ,t, s
P jL,min j
rij
(85)
xij
Um
1: kijvr
O _VR Nmax
Fig. 4. Voltage regulator model.
(81) 11
Uj
Applied Energy 253 (2019) 113517
J. Wang, et al.
B. Modified piecewise linearization approach
p j wt,t ( p j ,pvt )
Nonlinear polynomials aforementioned in (19), (20) and (37)–(39) can be handled by applying the method of modified piecewise linearization. _gas 2 ) . In As can be seen in (19), there exist nonlinear terms (ptsup ,s order to solve this problem, equivalent transformation is firstly made by (95) and (96).
Q
q
wt j ,t ,max
( q
pv j ,t ,max
) q
wt j ,t ,max
(q
pv j ,t ,max
Q
)
Fig. 5. Distributed generator inverter power control curve.
Sup Pgas
n + 1, t , s
In (86)–(88), auxiliary variables are introduced to replace the original square term.
t,
0
(86)
s
V j, t , s = V j2, t , s,
j
E _B ,
t,
s
(87)
Iij, t , s = Iij2, t , s ,
ij
Line ,
t,
s
(88)
j, t , s
j, t , s c 2 ij i, t , s ,
=
c 2 ij, k i, t , s
M (1
2 j,max ,
j
ij
C,
G _B ,
t,
xijc, k )
KC ,
t,
+ M (1
ij
0
2 Vmax ,
V j, t , s
2 Imax ,
Iij, t , s for:
m
ij
t,
t,
s
VR ,
Vm, t , s V j, t , s
E _B ,
Lori ,
m (ij ), ij
(Vmin/ kijvr,max ) 2 Mx ijvr
j
t,
s
Mx ijvr
s
(96)
s
(97)
Sup Pgas
max ,
t,
s
(98)
2
+
0
Sup Sup [Pgas min, Pgas max ]
,
mn, t , s , mn, t , s
n, t , s ,
(99a)
n = 1, 2, ...,Nplg
1,
1, mn, t , s = {0, 1}, n = 1, 2, ...,Nplg ,
n, t , s
t,
t,
s
(99b) (99c)
s
Nplg 0
+
(
n
n 1) n, t , s ,
t,
(Cg ( n )
Cg (
s
(100a)
n=1
Nplg n 1)) n, t , s ,
t,
n=1
s
(100b)
t,
s
(101)
sup _elec , 1) pt , s
t,
s
(102)
Eq. (103) indicates the total electricity supply power at time stage t of scenario s. On the basis of (103) and constraint (57), a preliminary _elec boundary of ptsup can be obtained as (104) ,s
C,
_elec ptsup = ,s
I j2, t , s
pjSup ,t,s , j
Sup pelec
(92)
min
t,
s
(103)
E _Sup
_elec ptsup ,s
Sup pelec
max ,
t,
s
(104)
The approach of modified piecewise linearization is also applied in (105)–(106).
(93)
Ce (p) = 2 2 p2 +
s
(Vmax /kijvr,min ) 2
Vm, t , s
t,
_elec _elec E (ptsup ) = (2 2 ptsup + ,s ,s
Analogously, the nonlinear and can be eliminated by substituting V j, t , s and Iij, t , s into (61)–(64) and (75). Expressions (59), (60) and (76) can be reformulated by (92), (93) and (94), respectively. 2 Vmin
fp,
_Elec _elec CsSup = E (ptsup ) t, ,t ,s
(91)
s
terms V j2, t , s
x ijc, k ),
t,
Similarly, the reformulation on (20) is presented in (101) and (102) as well.
(90) c 2 ij, k i, t , s
f jSup ,t,s
_gas G (ptsup ) = Cg ( 0 ) + ,s
(89)
s
s
j, t , s
k
t,
1
_gas ptsup = ,s
Consequently, substituting j, t , s into expressions, (35), (43) and (45) can be rewritten as linear constraints (89)–(91), respectively. In a similar way, i2, t , s can be removed from (37)–(39) as well. 2 j,min
+
(95)
G _Sup
_gas ptsup ,s
min
Cg ( ) =
A. Auxiliary variable substitution
G _B ,
=(
sup _gas 1 pt , s
s
By applying the modified piecewise linearization, the boundary Sup Sup [Pgas min, Pgas max ] is equally divided into Npl g intervals. In order to ensure the equivalent transformation, two auxiliary variables are introduced: (i) a binary variable mn, t , s , which implies that whether the nth interval is selected, (ii) a continuous variable n, t , s which is the optimal incremental value in each interval. Their relationships are character_gas ized by (99), and the gas supply power ptsup as well as the gas pur,s _gas ) ) can be chase cost at each time stage of the scenario (i.e., G (ptsup ,s represented by the mentioned auxiliary variables respectively, as shown in (100a) and (100b). Based on this, the nonlinear polynomial function (96) can be rewritten as (100). Thus, Eq. (19) is finally replaced by (95) and (99)–(100).
In what follows, we will derive equivalent linear expressions for each of them.
j
sup _gas , 0 ) pt , s
j
a. Quadratic terms such as i2, t , s , Iij2, t , s , V j2, t , s b. The nonlinear expressions like Eq. (19), Eq. (20) and the constraints (37)–(39) c. The nonlinear power flow constraint (62) hs sh hs rh d. Bilinear terms consisting of x jgs, k ujsg, t , s , x jgs, k urg j, t , s , x j, k u j, t , s , x j, k u j, t , s , ess se ess re vr cb cb x j, k uj, t , s , x j, k uj, t , s , yj, t , s x j, k and kij, t , s Vm, t , s e. Items with absolute value, including the switching times limits for CBs and VRs and objective expression (9).
=
_gas G (ptsup ) ,s
_gas ptsup = ,s
The programming model established in this paper is a complicated mixed integer nonlinear programming (MINLP) with the following nonlinear terms:
j, t , s
t,
Given the flow-power coefficient fp , the total gas supply power at each time stage of the scenario can be computed by (97) and in (98), the _gas boundary of ptsup can be obtained which is based on (34) and (97). ,s
5. Equivalent reformulation—linearization
2 j, t , s ,
_Gas _gas CsSup = G (ptsup ) t, ,t ,s
n + 1, t , s
(94)
0 12
n, t , s
1 p,
mn, t , s , mn, t , s
p
Sup [pelec n,
min ,
Sup pelec
n = 1, 2, ...,Nple
1, n = 1, 2, ...,Nple ,
t,
s
(105a)
max ]
1,
t,
s
(105b) (105c)
Applied Energy 253 (2019) 113517
J. Wang, et al.
C. Second order cone relaxation
Nple
_elec ptsup = p0 + ,s
(pn
pn 1 )
n, t , s ,
t,
s
(106a)
n=1
Although Iij2, t , s and V j2, t , s in constraint (62) have been substituted by Iij, t , s and V j, t , s in part A, there remain some nonlinear items including pij2, t , s , qij2, t , s and Vi, t , s Iij, t , s . In order to cope with them, SOCR is applied to transform non-convex constraint (62) into (113), as done in [14].
Npl _elec E (ptsup ) = Ce (p0 ) + ,s
(Ce (pn )
Ce (pn 1 ))
n, t , s,
t,
s
(106b)
n=1
Consequently, Eq. (20) can be rewritten as (101) and (105)–(106). Even though constraints (37)–(39) have been addressed in the part A, there remain some issues to be solved due to the complicated nonlinear term sgn( i, t , s, j, t , s ) fij2, t , s . In order to solve this problem, equivalent transformation is firstly made in (107) where sgn( i, t , s, j, t , s ) fij2, t , s can be seen as fij, t , s |fij, t , s | and the sign of fij, t , s is up to the comparison of the pressure at both ends for the facts that nodal pressure is nonnegative and gas always flow from higher pressure node to lower pressure node.
sgn(
2 j, t , s ) fij, t , s
i, t , s
= fij, t , s |fij, t , s |
2pij, t , s 2qij, t , s Iij, t , s
Pipe Pipe [ Fmax , Fmax ]
Pipe Fmax = max{FijPipe ,max }
, ij
t,
s
0
(108)
mn, t , s , mn, t , s
n + 1, t , s
0
n, t , s
n, t , s ,
0
(109a)
Pipe Pipe [ Fmax , Fmax ]
n = 1, 2, ...,Npl
1, n = 1, 2, ...,Npl,
t,
1,
t,
(109c)
s
0
Based on above, the constraints (37)–(39) are rewritten as (110)–(112) in combination with M (a very large number).
fij, t , s = f0 + ij ( i, t , s Pipe ,
j, t , s )
t,
n= 1
= Y (f0 ) +
(fn Npl n= 1
j, t , s )
Npl n=1
Y (f0 ) +
(fn Npl n=1
j, t , s )
Y (f0 ) +
Npl n=1
k KPipe
Pipe ,
fij, t , s fij, t , s
f0 + f0 +
n=1 Npl n=1
i, t
Y (f0 ) +
ij ( i, t
j, t )
Y (f0 ) +
t,
s
t,
x ijpipe ,k
k KHS
j
s
GS ,
t,
s (114b)
hs PkHS ,max x j, k
j
Mujsh, t , s
k KHS
t,
(114a)
gs FkGS ,max x j, k
Murg j, t , s
pjrh, t , s
HS ,
t,
s (115a)
hs FkHS ,max x j , k
j
Murh j, t , s
,
k K GS
pjse, t , s
pjre, t , s 0
n, t , s
)
pjse, t , s 0
0
Y (fn 1 ))
k KGS
GS ,
HS ,
t,
s (115b)
n, t , s
n, t , s
Y (fn 1 ))
j
Mujsg, t , s
pjsh , t,s
pjrh, t , s
gs FkGS ,max x j, k
ess PkESS ,max x j , k
j
Mujse, t , s
k KES
pjre, t , s
ESS ,
t,
s (116a)
ess FkESS ,max x j, k
j
Mure j, t , s
ESS ,
t,
s (116b)
In constraint (117), M is used again to eliminate the nonlinear item yjcb, t , s x jcb, k in (72).
ij
for
n, t , s
j
QkCB yjcb, t , s
x ijpipe ,k )
CB ,
k
K CB ,
x jcb, k )
M (1 M
(111)
s
k KCB
x jcb, k
t,
s:
qjcb, t , s
QkCB yjcb, t , s + M (1
qjcb, t , s
M k K CB
x jcb, k )
x jcb, k
(117)
In addition, there still exists a non-convex item related to in VR operational constraint, which is a crucial problem to be solved. In order to solve it, a set of binary variables aijvr, n, t , s is introduced in (118) to formulate the VR tap position firstly.
kijvr, t , s2 Vm, t , s
Npl
j, t )
C,
Y (fn 1 ))
(Y (fn )
(Y (fn )
pjsh , t,s
0
ij
0
k KPipe
M (1
ij (
0
(Y (fn )
s (113)
k K GS
f jrg, t , s
n, t , s
,
fn 1 )
f jrg, t , s
0
(110)
+ M (1 ij ( i, t , s
fn 1 )
s
fij, t , s = f0 + ij ( i, t , s
Npl
t,
2
f jsg, t , s
0
(109b)
s
f jsg, t , s 0
Therefore, the method of modified piecewise linearization is applied as follows:
Y (f ) = f |f |, f
Line ,
x jgs, k ujsg, t , s and x jgs, k urg j, t , s in (48) are further resolved in (114). It still can be seen that the value of gas flow into/out of the gas storage unit ( f jsg, t , s / f jrg, t , s ) is limited by the boundary [0, FkGS ,max ] only when the k -type gas storage unit is installed at node j (x jgs, k = 1) and in the operating state of gas storage/release (/urg j, t , s = 1); otherwise, it is imposed to be 0 (only if x jgs, k = 0 or ujsg, t , s = 0 /urg j, t , s = 0 ). By the same way, thermal storage constraint (55) and ESS constraint (69) can be reformulated respectively as (115)–(116).
(107)
Pipe ,
ij
D. Bilinear term relaxation
Then, the range of fij, t , s is determined by (40), and hence gas flow through each pipeline is also subjected to the boundary as (108). It Pipe depends on the largest flow of all pipelines (Fmax ).
fij, t , s
Iij, t , s + Vi, t , s,
Vi, t , s
(fn
fn 1 )
n, t , s
(fn
fn 1 )
n, t , s
Npl n=1 Npl n=1
+M k KC
M k KC
(Y (fn )
Y (fn 1 ))
n, t , s
(Y (fn )
Y (fn 1 ))
n, t , s
x ijc, k x ijc, k , +M k KC
M k KC
ij
yijvr, t , s =
x ijc, k
Rijvr +
2Rijvr
aijvr, n, t , s n,
ij
VR ,
t,
s
n=0
2Rijvr
x ijc, k
aijvr, n, t , s = 1,
ij
VR ,
t,
s
n=0
(112)
Then, the square of the VR ratio can be replaced by (119). 13
(118a)
(118b)
Applied Energy 253 (2019) 113517
J. Wang, et al. 2Rijvr
kijvr, t , s 2 =
1+
Rijvr
n
2Rijvr
n=0
Xijvr = kijvr % × 2Rijvr ,
2
Xijvr
ij
aijvr, n, t , s , ij
VR ,
t,
s
V j, t , s =
(119a)
M (1
V j, t , s =
1+ n= 0
j
VR ,
Rijvr
n
m
2Rijvr t,
s
aijvr, n, t , s Vm, t , s,
m
Rijvr
Xijvr ) 2uijvr, n, t , s
2Rijvr
Vm, t , s
uijvr, n, t , s
m (ij), ij
M (1
aijvr, n, t , s )
Maijvr, n, t , s VR ,
t,
s
(121)
E. Absolute value processing
2
Xijvr
n
uijvr, n, t , s
Maijvr, n, t , s
After that, the constraint (78) can be rewritten as (120). 2Rijvr
(1 +
n=0
aijvr, n, t , s )
(119b)
VR
2Rijvr
m (ij ), i
In (122), integer variables n+j, tcb, s and nj, tcb, s are introduced for linearizing the absolute value |yjcb, t , s yjcb, t 1, s | in (73).
(120)
for
Furthermore, the bilinear term aijvr, n, t , s Vm, t , s in (120) is also an important problem, which can be solved by introducing uijvr, n, t , s and actually the values of uijvr, n, t , s and aijvr, n, t Vm, t , s are equal. In order to ensure that, the equivalent transformation can be finally shown as:
j
CB ,
n+j, tcb, s
0, nj, tcb, s
n+j, tcb, s
nj, tcb, s =
s:
0, yjcb, t , s yjcb, t , s
(122a)
t 0 yjcb, t 1, s
Fig. 6. Flowchart of the proposed modeling methodology. 14
t=1 t
2
(122b)
Applied Energy 253 (2019) 113517
J. Wang, et al. NT s
n+j, tcb, s + nj, tcb, s
min
O _CB Nmax
pipe gs gs ess lup lnew cb svc vr wt {xij, k , xijc, k , x j, k , x j, k , x hs j, k , x j, k , xij, k , xij, k , x j, k , x j, k , xij , nj, k,
(122c)
t=1
pv +vr cb vr vr nj, k , ncb j, k, y j, t , s , yij, t , s , rij, t , s, nij, t , s, j, t , s, V j, t , s,Iij, t , s, Vm, t , s, sg rg + cb cb se re vr rh nij, t , s, nj, t , s, n j, t , s, u j, t , s, uj, t , s, ush j, t , s, u j, t , s, uj, t , s, u j, t , s }
Analogously, absolute value in (81) can be substituted by (123).
for
ij
nij+, tvr, s
0,
nij+, tvr, s
subject to : 1.linear constraints: {(1)
s:
nij, tvr, s
0,
yijvr, t , s
nij+, tvr, s + nij, tvr, s
+
VR
yijvr, t
t
1, s
2
(123b)
O _VR Nmax
+ j
SVC
c O _svcqjsvc , t, s +
Line
j
PVG
j
L
+
(75), (82) (95), (97), (117), (121) (124)}
In order to give a clear description of the proposed modeling methodology, a flowchart is presented in Fig. 6. The modeling procedure listed in block diagram can be summarized as follows: Step 1: Uncertainty modeling. In order to describe the operating conditions, a probabilistic scenario generation method is developed with the consideration of extreme and normal situations. The details which have been illustrated in Section 3 are briefly presented in the right side of Fig. 6 as well. Step 2: Multi energy system modeling. To establish the comprehensive planning model, the NGN, ADS models are investigated separately and then EH is modeled for coupling these two networks. By incorporating NGN, ADN and EH, an expansion planning model for the MES can be formulated as an MINLP problem. The specific modeling has been clearly given in Section 4. Step 3: Equivalent reformulation. Since the proposed MINLP model is difficult to be solved, the mathematical techniques are developed to transform the original version into a MISOCP. The relevant processing methods and procedures have been illustrated in Section 5. Step 4: Solve the model by calling CPLEX.
j
CB
cb cb cO _cb (n+ j, t + n j, t )
cO _lossIij2, t , s rij
j
ESS _pv c Cur _pvpjcur ,t ,s
+
(42), (44),
6. Framework
(123c)
vr c O _vr (nij+, vr t + nij, t ) +
t ij
(34), (36), (40)
t=1
According to the analysis (122) and (123), Eq. (9) can be rewritten as (124). ij
(14), (21)
(125)
0
t=1
_Elec CsOpe = ,t
(8), (10)
(46) (47), (49) (54), (56) (58), (61), (63) (68), (70) (71), (74) (99) (101), (103), (105) (106), (109) (112), (114) 2.conic constraints: \{ (113)\}
(123a)
t
yijvr, t , s
nij, tvr, s =
NT s
VR ,
CTotal
_L c Cur _Lpjcur , t, s +
j
cO _ess (pjse, t , s + pjre, t , s ) +
cO _wtpjwt , t, s + WTG
j
WTG
_wt c Cur _wtpjcur ,t ,s
cO _pvpjpv ,t ,s
j
(124)
PVG
Thus, the complex nonconvex MINLP has been transformed as a MISOCP model, which is formulated as:
19 1
18 17
16
20 2
37 36
3
33
23
32
31
15
Feeders to be upgraded Candidate WTG
Candidate PVG
22
4
5
24 30
13
14
Substation
21
35
7
6
8
26
25 29
34
28
12
9
27
11
10
Slack bus
Load bus
Newly added bus
Existing feeders
Hub Unit (No.1)
Candidate VR
Candidate ESS
Candidate CB
Candidate SVC
Fig. 7. Schematic of the modified IEEE 33-node system (adding 4 new load bus, 8 candidate feeders and other potential candidate units). 15
Applied Energy 253 (2019) 113517
J. Wang, et al.
Fig. 8. Schematic of the 23-node gas network system (with 3 new load bus, 7 candidate pipelines, 2 potential candidate compressors and 2 candidate gas storage units). Table 1 Parameters for the candidate distribution feeders.
Table 4 Parameters for available alternatives of CB and SVC units.
No.
From bus
To bus
Length (km)
Resistance (Ω)
Reactance (Ω)
CB
33 34 35 36 37 38 39 40
7 8 21 22 31 32 18 33
34 34 34 35 35 35 36 37
1.00 0.75 0.63 1.13 0.33 0.75 0.50 0.88
0.40 0.30 0.25 0.45 0.15 0.30 0.20 0.35
0.19 0.23 0.19 0.38 0.11 0.23 0.15 0.26
Capacity (MVar)
Cost (103 $/unit)
Capacity (MVar)
Cost (103 $/unit)
0.05 0.12
200 500
0.40 0.80
8000 14,000
I II
Table 5 Parameters for available alternatives of ESS units.
Table 2 Parameters for available alternatives of feeders. Feeders to be upgraded
Ⅰ II
I II
Feeders to be added
Capacity (MW)
Cost (103 $/km)
Capacity (MW)
Cost (103 $/km)
1.4 3.0
12,000 14,000
1.4 4.0
15,000 22,000
SVC
Capacity (MW)
Power limit (MW)
Efficiency
Cost (103 $)
10 30
3 9
98% 98%
6000 10,000
system is employed, as shown in Fig. 7 and Fig. 8 respectively. Simulations have been implemented by using solver CPLEX 12.6.0. 7.1. Overview of the test system
Table 3 Parameters for available alternatives of DG units. WTG
I II
In the original system, there are 33 buses and 32 feeders totally, and corresponding system parameters can be found in [33]. As can be seen from Fig. 7, compared with the IEEE 33-node network, four load buses (with the same power 220 kW and 110 kVar) and 8 candidate feeders are added to the modified distribution system. Meanwhile, several potential locations (i.e., buses or feeders) for the investment of candidate active equipment are marked in the system. The related parameters for all alternative components are listed in Tables 1–5 respectively. According to Fig. 8, the presented gas system is composed of 2 gas source nodes, 2 gas compressors, 26 gas pipelines, 20 existing gas load nodes and 3 newly added nodes. However, in order to meet the growing demand for gas, expansion investment is still needed. Hence, 7 candidate gas pipelines (marked with black dashed line), 2 candidate gas compressors and 2 candidate gas storage units are selected to be
PVG
Rated power (MVA)
Cost (103 $)
Rated power (MVA)
Cost (103 $)
2 4
22,000 31,000
3 5
24,000 29,000
7. Case study To demonstrate the effectiveness of the proposed planning model, a modified IEEE 33-node distribution system coupled with a 23-node gas 16
Applied Energy 253 (2019) 113517
J. Wang, et al.
Table 6 Data of all gas pipelines. Existing gas pipelines
Candidate gas pipelined
No.
Pipeline (i–j)
Length (km)
ϕij (1012)
No.
Pipeline (i–j)
Length (km)
ϕij (1012)
No.
Pipeline (i–j)
Length (km)
ϕij (1012)
1 2 3 4 5 6 7 8 9 10 11 12 13
1–2 1–2 2–3 2–3 3–4 5–6 5–8 6–7 4–7 4–14 8–9 8–9 9–10
0.4 0.4 0.6 0.6 2.6 4.3 2.6 2.9 1.9 5.5 0.5 0.5 0.2
0.0052 0.0052 0.0026 0.0052 0.0035 0.0017 0.0035 0.0024 0.0035 0.0017 0.0035 0.0017 0.0035
14 15 16 17 18 19 20 21 22 23 24 25 26
9–10 10–11 10–11 11–12 12–13 13–14 14–15 15–16 11–17 11–17 17–18 18–19 19–20
2 2.5 2.5 4.2 4 0.5 1.0 2.5 1.05 1.05 2.6 9.8 0.6
0.0069 0.0069 0.0069 0.0035 0.0017 0.0122 0.0035 0.0017 0.0026 0.0026 0.0026 0.0035 0.0035
27 28 29 30 31 32 33
18–21 19–21 20–21 12–22 13–22 15–23 16–23
0.6 1.05 1.05 2.6 2.5 0.5 0.5
0.0035 0.0026 0.0026 0.0026 0.0069 0.0122 0.0035
Table 7 Parameters for available alternatives of gas pipelines and gas compressors. Candidate gas compressor
I II
Table 10 Parameters for available alternatives of thermal storage units.
Candidate gas pipeline
Flow limit (106 m3/h)
Ratio
Cost (103 $)
Flow limit (106 m3/h)
Cost (103 $)
10 10
1.03 1.08
10,000 15,000
1.2 6
10,000 15,000
I II
Thermal power limit (MW)
Efficiency
Cost (103 $)
10 30
3 9
95% 95%
4200 8000
7.2. Scenario set Based on the generated scenarios proposed in Section 3, the scenario set can be constituted as follows:
Table 8 Parameters for available alternatives of gas storage units.
I II
Capacity (MWh)
Capacity (m6)
Flow limit (106 m3/h)
Efficiency
Cost (103 $)
0.06 0.12
0.005 0.015
95% 95%
10,000 20,000
(1) Scenario ΦA: The combination of the maximum value curves in Fig. 2 and the minimum value curves in Fig. 3. (2) Scenario ΦB: The combination of the minimum value curves in Fig. 2 and the maximum value curves in Fig. 3. (3) Scenario ΦC: The combination of the mean value curves in all graph.
available investment options, as shown in Tables 6–8. In addition, ADN and GN systems are coupled by totally 8 energy hub units. The related data of hub unit is given in Table 9. It is worth mentioning that thermal storage units are only optional investments embraced in hub units, whose parameters are shown in Table 10. Table 11 depicts other supplement parameters in the original system. Meanwhile, there are some reasonable assumptions as follows: (1) Annual maintenance costs for all units are set to be to 5% of investment cost. (2) Unit production cost for WTGs and PVGs are assumed to be 0. (3) Unit costs of energy losses of TSs and feeders are all 200 $/MWh. (4) Unit active curtailment costs of wind, photovoltaic powers are both 3000 $/MWh and unit load-shedding cost is 5000 $/MWh. (5) All existing feeders are considered to be upgraded, and all load buses are schedulable except the slack bus. (6) The planning scheme is implemented in 10 years and the interest rate (d) is 10%. (7) Electricity, gas and thermal demands are doubled in next 10 years. (8) Big constant M is 10000000.
By applying the Step 4 in Section 3, the corresponding probabilities of the 3 scenarios are 0.0465, 0.1079 and 0.8457, respectively. 7.3. Results of expansion planning model In this section, the aforementioned probabilistic scenario set is employed to test the model. In order to analyze the performance of the proposed model, four cases are compared and presented as follows: Case 1: model (125) with probabilistic scenarios: ΦA, ΦB and ΦC. Case 2: model (125) with mean value scenario ΦC only. Case 3: model (125) with three seasonal average scenarios: daily average scenarios for summer, intermediate (spring and autumn) and winter. The corresponding probabilities are 0.25, 0.5 and 0.25, respectively. Case 4: model (125) with five scenarios: two extreme scenarios (ΦA and ΦB) and three average seasonal scenarios with probabilities 0.0465, 0.1079, 0.2281, 0.3926 and 0.2250, respectively.
Table 9 Data of energy hub units. No.
1
2
3
4
5
6
7
8
Bus number in ADN Node number in GN Heat demand (MW)
10 14 1.07
14 10 2.15
2 12 0.54
28 16 2.68
32 6 1.61
4 15 0.54
21 2 1.07
23 11 1.13
Conversion efficiency
P2G: 65%, GF: 92%, CHP: 1) From gas to heat 36%; 2) From gas to electricity 18%
17
Applied Energy 253 (2019) 113517
J. Wang, et al.
Table 11 Other parameter setting in the original system. Gas system
Active distribution system
Gas nodal pressure limits Gas flow limits (in existing pipelines) Gas source flow limits GHV (ζfp) Ratio of existing compressor ( ijc )
50–70 bar 0–5 m6/h 0–20 m6/h 0.0096 MW/m3 1.1
Voltage magnitude range Maximum current magnitude Maximum power in existing feeders Maximum switch times of CB and VR Maximum investment numbers of SVCs, CBs, WTs and PVs
Gas price coefficients
γ1 = 0.5 γ0 = 2
DG power factor (
0.94–1.06 pu 3 pu 1.6 MW/1.2Mvar 24
SVC CB Nmax = 1, Nmax = 16 ,
wt pv Nmax = 5, Nmax =5 0.8
wt / pv )
Maximum curtailment rates of WTs, PVs and demand Electricity price coefficients
20% 2 = 0.3 1 = 1.2 0 = 0.1
Table 12 Investment decision of all components (including type, location and installation number).
Active distribution network
Infrastructures
Type
Case 1
Upgraded feeder
I II I II I II I II I II I II I I II
23(1), 6(1) 34(1), 37(1), – 13(1), – 6(1) – – 33(5) 3(1) – 1(1) 7(1) –
Added feeder WTG PVG CB SVC VR ESS
24(1) 35(1), 40(1) 36(1)
Case 2
Case 3
Case 4
22(1), 23(1) – 34(1), 35(1), 37(1), 40(1) – 13(1) – – 8(1) – 33(8) 3(1) – 1(1) 7(1) –
14(1), 23(1) 6(1) 34(1), 35(1), 37(1) 39(1) 13(1) – 6(1) – – 33(6) 3(1) – 1(1) 7(1) –
23(1) 8(1), 9(1) 34(1), 35(1), 38(1), 40(1) – 13(1), 36(1) – 8(1) – – 33(6) 3(1) – 1(1) 7(1) –
Hub unit
Thermal storage
I II
3(1), 8(1) –
8(1) –
8(1) –
7(1), 8(1)
Natural gas network
Added gas pipeline
I
1(1), 5(1), 7(1) – – 8(1) 9(1) –
2(1), 4(1), 6(1) – 8(1) – 9(1) –
1(1), 5(1), 6(1) – 8(1) –
1(1), 5(1), 7(1) – 8(1) – – –
Gas storage unit Compressor
II I II I II
Note: number within () is the amount of investment; number outside () is the location.
thermal storage and gas storage units are installed to satisfy the operation constraints. Compared with Case 1, two feeders are required to be upgraded and there is less development in renewable DGs. Since investment cost (see Table 13) in active management equipment is far less than Case 1, Case 2 has a lower total planning cost; however, its annual energy purchase cost has a 13.1% increase compared with Case 1, which means that Case 2, with less investment decisions, requires more energy consumption than Case 1 during the planning horizon. In this respect, the probability model with extreme scenario set has a positive impact on reducing the energy purchase cost, but it increases the cost of investment and makes the 10-year planning cost slightly increase by $ 5.61 × 106. It should be noted that there is usually a trade-off between the investment costs and the operation performance (including energy consumption during planning period). Although less investment can be obtained under the average scenario ΦC (i.e., Case 2), such an investment scheme may lead to a significant increase in comprehensive operation cost in some extreme cases, and even can not meet the requirements of system operation under some circumstances. On the contrary, the probabilistic model in Case 1 synthesizes the extreme scenarios and considers the robustness of the investment decision, which makes the investment scheme competent for the occurrence of
Table 13 Planning cost component comparisons with the consideration of different scenario sets ($ 108).
Investment cost Annual comprehensive operation cost Annual maintenance cost Annual energy purchase cost Annual DG curtailment cost Annual demand curtailment cost Total planning cost
Case 1
Case 2
Case 3
Case 4
1.5059 0.1624 0.0753 0.0496 0.0004 0.0002 2.6035
1.3668 0.1747 0.0683 0.0561 0.0005 0.0002 2.5474
1.3731 0.1650 0.0686 0.0553 0.0002 0.0002 2.5220
1.5209 0.1584 0.0765 0.0483 0.0001 0.0001 2.6018
Note: the annual comprehensive operation cost represents the sum of operation (including curtailment), maintenance and energy purchase costs.
Tables 12 and 13 present the simulation results including all the investment decisions and planning costs under Cases 1–4. In Table 12, the investment results of Case 1 indicate that during the planning horizon, the system has been developed to response the demand growth by upgrading three feeders, adding four new feeders, three gas pipelines and three DGs. In view of the objective that is to minimize the total costs, CB, SVC, VR, compressor and energy storage units including ESS, 18
Applied Energy 253 (2019) 113517
Input power (MW)
Input gas (103m3/h)
Input gas (10 3m3/h)
J. Wang, et al.
Scenario ΦA
0.2 0.1
0
Scenario
0.2
ΦB
0.1
0
12
24
0.4
0
0.1
0
0.5
12
(a) Input gas to GF
24
0
12
24
0
12
24
0
12
24
0.2
0
12
24
0
0.4
0.4
0.2
0.2
0
0 0.4
0.2 0
Scenario ΦC
0.2
0
12
24
0
0
12
(b) Input gas to CHP
24
0 0.4 0.2
0
12
24
(c) Input electricity to P2G
0
Thermal power Thermal residual capacity
0.05
0
0
5
10
15
20
Thermal power Thermal residual capacity
0.2
0.1
0
-0.1
0
5
10
15
20
Thermal power/capacity (pu)
Scenario ΦB
Scenario ΦA
0.1
Thermal power/capacity (pu)
Thermal power/capacity (pu)
Fig. 9. Energy input to the EH (No.1) in each time stage of the 3 different scenarios.
Scenario ΦC
0.1
Thermal power Thermal residual capacity
0.05
0
0
5
10
15
20
Fig. 10. Storing/Releasing thermal power status and capacity for thermal storage unit in the EH (No. 8) at each time stages in 3 different scenarios.
any uncertain scenarios. For further investigating the impact of seasonal scenarios on the planning results, three scenarios (summer, intermediate and winter) are calculated and constituted by averaging the hourly power values throughout all days in each season. Note that spring and autumn are together considered as the intermediate scenario. In Case 4, the probabilities of summer, intermediate and winter scenarios are obtained by counting their frequency of occurrences except for the extreme hours in ΦA and ΦB. Without the consideration of extreme scenarios, the resulting schemes in Case 3 has one less DG investment than that in both Case 1 and Case 4, thereby allowing a fewer investment cost. Moreover, it makes no great difference to each resulting cost between Case 2 and Case 3, even if Case 3 uses a combination of seasonal scenarios. It can also be observed that the results in Case 4 are relatively similar to that in Case 1 (the variations in investment, comprehensive operation, energy purchase and total planning costs are no more than 1%, 2.44%, 2.63% and 0.066%, respectively), indicating that the mean value scenario ΦC (total 24 h/time stages) has similar effect to the three seasonal average scenarios (total 3 × 24 h/time stages). Since the reduction of scenario number is significant for a complicated long-term planning issue, using the combination of ΦA, ΦB and ΦC is better than the
combination of ΦA, ΦB and three seasonal scenarios via fewer operating conditions. According to the above discussion, the different results on Cases 1–4 show the dominant status of extreme scenarios ΦA and ΦB over other scenarios. To make further analysis of the energy conversion in EH, the simulations of Case 1 are presented in Figs. 9, 10 and 11. Among them, Energy input changes to EH in each scenario are shown in Fig. 9; the operation of thermal storage unit is shown in Fig. 10; and the variation curves of electricity and gas price are shown in Fig. 11. It can be observed that: (1) In the case of scenario ΦA, due to the abundant wind and solar resources and the low electrical demand (especially from 12 h to 20 h), as well as the low electricity price level (see Fig. 11), the additional power generation is converted into natural gas through P2G, which together supplies natural gas demand; in the mean time, thermal demand is mainly provided by GF. During the rest of the period, thermal demand is together provided by natural gas through GF and CHP (see the left side of Fig. 9). In addition, the thermal storage unit (see Fig. 10) stores the excess heat energy from 1 h to 5 h, and then releases them from 18 h to 24 h to respond the heat demand peak at night. 19
Applied Energy 253 (2019) 113517
Electrical price $/(103MWh)
Gas price $/(103m3)
J. Wang, et al.
1000
Scenario ΦC
Scenario ΦB
Scenario ΦA
500 0
12
120
24
12
time (h) (a) Gas price Scenario ΦB
Scenario ΦA
80
24
12
24
Scenario ΦC
40 0
12
24
12
24
time (h) (b) Electricity price
12
24
Fig. 11. Energy price variation curves with time stages under three scenarios.
meet the higher demand: storing energy from 1 h to 9 h, and releasing them from 10 h to 24 h. (3) In the normal Scenario ΦC, although it has a lower requirement of energy supplies than Scenario ΦB, because of the low price of natural gas and the the overall high level electricity price in this scenario (see Fig. 11), natrual gas plays a leading role in the energy conversion and supply. Most of the gas injected into CHP to supply electrical and heat demand, while relatively little gas inputs GF to secondarily supply heat demand from moring to night (see the right side of Fig. 9). Simarly to other scenarios, heat release occurs mainly at night (see Fig. 10).
Table 14 Planning cost components for the separate models (ADN and NGN) ($ 108).
Investment cost Annual maintenance cost Annual energy purchase cost Annual DG curtailment cost Annual demand curtailment cost Total planning cost
The single ADN model
The single NGN model
Sum
1.5213 0.0761 0.0496 0.0126 0.0022
0.5440 0.0272 0.0430 – –
2.0653 0.1033 0.0926 0.0126 0.0022
3.2394
1.0186
4.2580
7.4. Comparison of the comprehensive model and separate models
Table 15 Optimal planning costs of different cases with full/partial consideration of ANMs (unit $ 106).
VR adjustment ESS scheduling CB and SVC control DG control Investment cost Annual maintenance cost Annual comprehensive operation cost Annual energy purchase cost Annual active curtailment cost Total planning cost
Case 1
Case 5
Case 6
✓ ✓ ✓ ✓ 150.5914 7.5296 16.2391
✓ ✓
✓ ✓ ✓
4.9636 0.062 260.3517
6.3366 0.0276 292.2024
✓ 148.9606 7.4480 21.3329
To verify the superiority of the comprehensive planning model where ADN and NGN are optimized coordinately (i.e., Case 1), the separate planning models for single ADN and NGN are conducted without the EH energy conversion, respectively. The cooresponding results on ADN and NGN are given in Table 14. It is obvious that the separate models present less economic than the joint expansiong planning model in Case 1. Specifically, compared with Case 1, there are a 37.15% increase of investment cost and a 86.69% increase of energy purchase cost; thus, the sum of planning cost is increased by 63.55%. Apparently this results show the importance and necessity of considering the interdependency of ADN and NGN in the MES planing problem.
Case 7
138.5041 6.9252 20.4271
✓ ✓ 148.0040 7.4002 19.2761
6.1791 0.0338 273.4036
6.2141 0.0163 276.0944
7.5. Analysis of the active network management
Note: the annual comprehensive operation cost represents the sum of operation (including curtailment), maintenance and energy purchase costs.
In order to verify the advantages of ANMs, the proposed model with full consideration of ANMs is compared with the model with partial consideration of ANMs. The different cases and results conducted based on the probabilistic scenario set are listed in Table 15. Since the full ANM schemes are considered in the proposed model, the investment cost of Case 1 is relatively higher than other cases in Table 15, however the comprehensive cost of annual operation cost as well as the total planning cost is lower. It is obvious that the model with full consideration of ANMs is more economical than that with partial consideration of ANMs, since ANMs make the multi-energy system more controllable and adjustable in controlling active and reactive power. By comparing Case 1 with Case 5, we can observe that the
(2) In the extreme Scenario ΦB, since the electrical demand outstrips the renewable energy supply, ADN needs more natural gas to inject into CHP (see the second graph of (b) in Fig. 9) for converting the gas into electricity. Due to thermal power along with CHP operating, the heat requirements are supplied by natrual gas primarily through CHP and secondarily through GF. It can also be seen that P2G has no power input in the meantime. Also, energy prices are higher than other scenarios in most of the time (espeacially 0–10 h and 16–24 h) according to Fig. 11, indicating a shortage of energy supplies due to the high-level load demand. Moreover, the thermal storage uint (see Fig. 10) works more frequently through the day to 20
Applied Energy 253 (2019) 113517
ΦC
Electric Demand
0
12
24
(a)
24
WTG
PVG
12
24 WTG(MW)
1.1
12
1 0
12
24
1
12
(c)
24
12
0
12
24
12
(e)
24
12
5
0
12
24
12
(g)
24
12
24
0
Scenario
Curtailment
12
ΦB
Scenario ΦC
Without curtailment
24
2
12
(b)
With curtailment
24
12
24
1 Curtailment
12
Without curtailment
24
12
(d)
With curtailment
24
12
24
2 1 0
24
Scenario ΦA
0.2
0
24
0.5
0.4
Curtailment
12
Without curtailment
24
0.2
12
(f )
With curtailment
24
12
24 1
0
-0.2
0 ESS power
12
time(h)
ESS capacity
24
12
(h)
24
12
24
-1
ESS capacity(pu)
Scenario
Demand(MW)
ΦB
PVG(MW)
SVC(MVar)
Scenario
0.5
0.9
Number of CBs
Scenario ΦA
1
ESS power (pu)
VR ratio
Active power(pu)
J. Wang, et al.
time(h)
Fig. 12. The 24-h operation schemes for each of the 3 scenarios considering active distribution management. (a) electric demand, wind and photovoltaic power generation curves; (b) load-shedding for bus 34; (c) ratio changes of VR between buses 1 and 2; (d) active power curtailment for WTG at bus 13; (e) reactive power output of SVC at bus 3; (f) active power curtailment for PVG at bus 6; (g) reactive power output of CBs at bus 33; (h) charging/discharging power and capacity of ESS at bus 7.
model without the reactive compensation of CBs and SVCs increases the comprehensive operation cost by 31.37%, and total planning cost by 12.09%. Furthermore, 27.66% more energy consumptions are required in Case 5 during the planning horizon. In Case 6, without considering the control of DG, the results also show an increase in the operation, energy purchase, and total planning costs by 25.79%, 24.45%, 5.01%, respectively. It indicates that the flexible control of DG has a positive impact on reducing the operation and planning cost for MES. Since the investment costs for DG are higher than that of other infrastructures, the total investment cost of Case 6 without DG is lower than that of other cases. Similarly, compared with the full consideration of ESS and VR, the more expensive cost results in Case 7 show that the model with the consideration of ESS and VR can also enjoy more economic value. In order to more intuitively reflect the influence of ANM schemes, Fig. 12 gives the results of the ANM operation strategies in the 24-h period under each of the three scenarios in Case 1. The power variation curves of electric demand, wind and photovoltaic power curves are presented in Fig. 12(a). The corresponding curtailment results are given in Fig. 12(b), (d) and (f), respectively. From Fig. 12(b), we can observe that the load-shedding is only adopted in extreme scenario ΦB when the electric demand attains peaks (16–20 h); while the load demand in scenarios ΦA and ΦC is relatively lower and no load-shedding is required in this circumstance. In addition, it can be noted from Fig. 12(d) as well as Fig. 12(f) that DGs are only curtailed in extreme scenario ΦA during the period of high-level DG output (16–20 h for WTG and 9–15 h for PVG). In comparison, since the power output is low in scenarios ΦB and ΦC, there is no curtailment and the renewable energy can be utilized fully. Fig. 12(c), (e) and (g) show the changes of VR tap position, SVC output and CB operating status with the variation of time stages respectively. Since the VR is installed closed to substation, the voltage can be controlled to adapt the changes of demand. If the total DG power output is lower than the load demand, especially in scenario ΦB, the overall VR tap is modulated in a lower position than that in ΦA and ΦB in order to limit the voltage drops, and restrict the energy loss over feeders. Furthermore, the reactive power output of SVC unit is kept
basic stability except for 20–24 h in ΦA since the high proportion DG power injection. Similarly, the operating CB units are the most during scenario ΦB, while the CB compensation during scenario ΦA is in the lowest level. Finally, according to Fig. 12(h), it can be briefly concluded that EES is charged during off-peak electric demand periods or on-peak DG output periods and discharged during on-peak demand periods, which can effectively smooth the power fluctuation caused by the injection of fluctuating loads or DG power. In detail, (i) in scenario ΦA, EES is continuously charged from 8 h to 16 h (during the on-peak DG output period) and discharged from 17 h to 24 h (during the on-peak power consumption period); (ii) in scenario ΦB, ESS is charged during the period from 1 h to 6 h (in the night when power consumption is low) and discharged from17h to 20 h (i.e., when power consumption attains peaks); however, due to the low DG output and high load in scenario ΦB, charging and discharging of ESS are not as obvious as scenarios ΦA and ΦC; (iii) in scenario ΦC, EES is charged during the period 10–15 h (along with large PVG output) and discharged during the period 17–24 h (along with the peak power consumption). 8. Conclusion This work proposes a new expansion planning model for multi-energy system where active distribution network, natural gas network and energy hub are incorporated. In this paper, we have shown how this comprehensive expansion planning problem can be modeled and formulated. This paper addresses the complex mixed integer nonlinear programming problem by applying second order cone relaxation as well as a modified piecewise linearization approach. With the objective of minimizing the total cost over the planning horizon, the final planning solution identifies the optimal determination of the type, location and size of all infrastructures, and hence their optimal operation strategy can be obtained simultaneously. Numerical computational results and analysis are carried out based on a multi-energy system with the integration of a modified IEEE33-node system, a 23-node gas system and 8 energy hub units. 21
Applied Energy 253 (2019) 113517
J. Wang, et al.
In order to conduct the simulations, a probabilistic scenario generation approach is proposed with the consideration of extreme cases. The proposed method has been verified to provide an efficient and robust scenario set for the expansion planning model. Compared with usual scenarios, the probabilistic scenario set increases the investment opportunities and reduces the energy consumption costs greatly. Additionally, the results based on different scenarios show the dominant impact of extreme scenarios on the planning decisions. Since there is usually a trade-off between the investment costs and the operation performance, several methods of scenario generation in this paper allow the planners to make a reasonable choice according to the actual situation. Moreover, the comparison of the comprehensive model and separate models has also been conducted in the case studies, and the results show the great superiority of the proposed joint expansion planning model over the separate models in the economic aspect. In addition, the impact of active network management on expansion planning of multi-energy system is investigated originally. The numerical results illustrate the superiority of considering the impact of active network managements on multi-energy system expansion planning. It has been proven that the full consideration of active network management can reduce the annual operation cost as well as the total planning cost. Given the fact that the proposed model with active network management enjoys more economic value, the incorporation of active distribution network in the joint multi-energy system expansion planning is indispensable. In summary, we believe the multi-energy system with consideration of active network managements in line with the emerging trend of active distribution network will be a promising research direction in the future.
[10] Rahmani-Andebili M. Distributed generation placement planning modeling feeder's failure rate and customer’s load type. IEEE Trans Ind Electron 2016;63(3):1598–606. [11] Rahmani-Andebili M. Simultaneous placement of DG and capacitor in distribution network. Electr Power Syst Res 2016;131:1–10. [12] Mokryani G, Hu YF, Papadopoulos P, Niknam T, Aghaei J. Deterministic approach for active distribution networks planning with high penetration of wind and solar power. Renew Energy 2017;113:942–51. [13] Koutsoukis NC, Georgilakis PS, Hatziargyriou ND. Multistage coordinated planning of active distribution networks. IEEE Trans Power Syst 2017;33(1):32–44. [14] Xie S, Hu Z, Zhou D, Li Y, Kong S, Lin W, et al. Multi-objective active distribution networks expansion planning by scenario-based stochastic programming considering uncertain and random weight of network. Appl Energy 2018;219:207–25. [15] Wu M, Kou L, Hou X, Ji Y, Xu B, Gao H. A bi-level robust planning model for active distribution networks considering uncertainties of renewable energies. Int J Electr Power Energy Syst 2019;105:814–22. [16] Moradijoz M, Moghaddam MP, Haghifam MR. A flexible active distribution system expansion planning model: a risk-based approach. Energy. 2018;145:442–57. [17] Rahmani-Andebili M. Optimal power factor for optimally located and sized solar parking lots applying quantum annealing. IET Gener Transm Distrib 2016;10(10):2538–47. [18] Zhang S, Cheng H, Wang D, Zhang L, Li F, Yao L. Distributed generation planning in active distribution network considering demand side management and network reconfiguration. Appl Energy 2018;228:1921–36. [19] Zeng B, Wen J, Shi J, Zhang J, Zhang Y. A multi-level approach to active distribution system planning for efficient renewable energy harvesting in a deregulated environment. Energy. 2016;96:614–24. [20] Shen X, Shahidehpour M, Han Y, Zhu S, Zheng J. Expansion planning of active distribution networks with centralized and distributed energy storage systems. IEEE Trans Sustain Energy 2017;8(1):126–34. [21] Xing H, Cheng H, Zhang Y, Zeng P. Active distribution network expansion planning integrating dispersed energy storage systems. IET Gener Transm Distrib 2016;10(3):638–44. [22] Andebili MR, Fotuhi-Firuzabad M. A robust approach for PEVs charging management and reconfiguration of electrical distribution system penetrated by renewables. IEEE Trans Ind Inf 2017;PP(99):1. [23] Rahmani-Andebili M. Dynamic and adaptive reconfiguration of electrical distribution system including renewables applying stochastic model predictive control. IET Gener Transm Distrib 2017;11(16):3912–21. [24] Rahmani-Andebili M. Stochastic, adaptive, and dynamic control of energy storage systems integrated with renewable energy sources for power loss minimization. Renew Energy 2017;113:1462–71. [25] Dehghan S, Amjady N, Kazemi A. Two-stage robust generation expansion planning: a mixed integer linear programming model. IEEE Trans Power Syst 2014;29(2):584–97. [26] Muñoz-Delgado G, Contreras J, Arroyo JM. Multistage generation and network expansion planning in distribution systems considering uncertainty and reliability. IEEE Trans Power Syst 2016;31(5):3715–28. [27] Bai L, Jiang T, Li F, Chen H, Li X. Distributed energy storage planning in soft open point based active distribution networks incorporating network reconfiguration and DG reactive power capability. Appl Energy 2017;210. [28] Markowitz HM, Manne AS. On the solution of discrete programming problems. Econometrica 1957;25(1):84–110. [29] Álvarez-Miranda E, Campos-Valdés C, Rahmann C. Two-stage robust UC including a novel scenario-based uncertainty model for wind power applications. Energy Convers Manage 2015;101:94–105. [30] Rahmani-Andebili M. Reliability and economic-driven switchable capacitor placement in distribution network. Gen Transm Distrib Iet 2015;9(13):1572–9. [31] Sheikhi A, Bahrami S, Ranjbar AM. An autonomous demand response program for electricity and natural gas networks in smart energy hubs. Energy 2015;89:490–9. [32] Wang C, Wei W, Wang J, Bi T. Convex optimization based adjustable robust dispatch for integrated electric-gas systems considering gas delivery priority. Appl Energy 2019;239:70–82. [33] Baran ME, Wu FF. Network reconfiguration in distribution systems for loss reduction and load balancing. IEEE Trans Power Delivery 1989;4(2):1401–7.
References [1] Vahid-Pakdel MJ, Nojavan S, Mohammadi-ivatloo B, Zare K. Stochastic optimization of energy hub operation with consideration of thermal energy market and demand response. Energy Convers Manage 2017;145:117–28. [2] Noorollahi Y, Yousefi H, Mohammadi M. Multi-criteria decision support system for wind farm site selection using GIS. Sustain Energy Technol Assess 2016;13:38–50. [3] Wang J, Zhong H, Ma Z, Xia Q, Kang C. Review and prospect of integrated demand response in the multi-energy system. Appl Energy 2017;202. [4] Mohammadi M, Noorollahi Y, Mohammadi-ivatloo B, Hosseinzadeh M, Yousefi H, Khorasani ST. Optimal management of energy hubs and smart energy hubs – a review. Renew Sustain Energy Rev 2018;89:33–50. [5] Beigvand SD, Abdi H, La Scala M. A general model for energy hub economic dispatch. Appl Energy 2017;190:1090–111. [6] Dagoumas AS, Koltsaklis NE. Review of models for integrating renewable energy in the generation expansion planning. Appl Energy 2019;242:1573–87. [7] Huang W, Zhang N, Yang J, Wang Y, Kang C. Optimal configuration planning of multi-energy systems considering distributed renewable energy. IEEE Trans Smart Grid 2018:1. [8] Fathtabar H, Barforoushi T, Shahabi M. Dynamic long-term expansion planning of generation resources and electric transmission network in multi-carrier energy systems. Int J Electr Power Energy Syst 2018;102:97–109. [9] Nazar MS, Haghifam MR. Multiobjective electric distribution system expansion planning using hybrid energy hub concept. Electr Power Syst Res 2009;79(6):899–911.
22
Contents lists available at ScienceDirect
Applied Energy journal homepage: www.elsevier.com/locate/apenergy
Expansion planning model of multi-energy system with the integration of active distribution network Jueying Wang, Zhijian Hu , Shiwei Xie ⁎
T
⁎
School of Electrical Engineering and Automation, Wuhan University, Wuhan City, Hubei Province, China
HIGHLIGHTS
for multi-energy system integrating active distribution network is modeled. • Planning scenario generation method is proposed considering the robustness. • AA probabilistic modified piecewise linearization method is developed. • Mixed integer second order cone programming is employed to solve the model. • ARTICLE INFO
ABSTRACT
Keywords: Multi-energy system expansion planning Active distribution network Active network management Probabilistic scenario generation Second-order cone programming
As an effective pattern to promote efficient use of energy, multi-energy system has become the main focus of development. On the other hand, active distribution network is now also the promoting concept for developing power distribution system, since the active network management technologies can make the network more flexible and controllable. To facilitate and utilize the advantages of both systems, this paper proposes an expansion planning model for multi-energy system integrating active distribution network, natural gas network and energy hub; and the positive impact of active network managements on expanding multi-energy system is originally investigated. With the aim of minimizing the total cost over the planning horizon, this model centers on the optimal determination of the type, location and size of all the infrastructures, where the active network managements are modeled and incorporated. In addition, in order to improve the robustness of the planning results, a probabilistic scenario generation approach is proposed to test the model. In order to solve the proposed mixed-integer nonlinear model, the second order cone programming as well as a modified piecewise linearization approach is applied to convert the original model to a mixed integer second order cone programming model. A multi-energy system (including a modified IEEE 33-node distribution system, a new 23-node gas system and 8 energy hubs) are employed to verify the effectiveness of the prosed model and methods. The simulation results exhibit the superiority of the joint expansion planning model for multi-energy system and the beneficial impact of considering active network management.
1. Introduction Global energy crisis and its relevant problems have led to a considerable movement into efficient utilization of renewable energy [1]. Therefore, making progress in supplying a sustainable, clean, secure and affordable energy is one of the basic challenges for industry in the 21st century [2]. In order to seek ways to deal with these challenges, the concept of multi-energy system (MES) [3] is introduced as an effective pattern to improve the energy efficiency and reduce energy loss by incorporating the multiple forms of energy. However, utilizing the benefits of MES requires an integrated modeling framework. Currently,
⁎
a new window has been opened in the field of MES optimization by presenting the promising concept of energy hub (EH), which is defined as a framework where production, consumption, conversion and storage of different forms of energy carriers can be done [4]. Based on the EH concept, energy systems can be investigated and controlled as a whole rather than as separate energy systems like distribution network, natural gas network (NGN) and so on. In addition, emerging new technologies such as combined heat and power (CHP) and power to gas (P2G) have made it possible to convert different forms of energy into each other, and made the EH concept practical for MES [5]. At present, the research on multi-energy system expansion planning
Corresponding authors. E-mail addresses: [email protected] (Z. Hu), [email protected] (S. Xie).
https://doi.org/10.1016/j.apenergy.2019.113517 Received 31 January 2019; Received in revised form 13 June 2019; Accepted 9 July 2019 0306-2619/ © 2019 Elsevier Ltd. All rights reserved.
Applied Energy 253 (2019) 113517
J. Wang, et al.
Nomenclature
units/thermal storage units/electric storage systems c M _lori/ ckM _lup/ ckM _lnew maintenance cost coefficients of original/replacing/feeder ckM _wt / ckM _pv / ckM _cb/ maintenance cost coefficients of WTG/PVG ckM _svc / ckM _vr /CB/SVC/VR gas price coefficient 1/ 0 / 1/ 0 electricity price coefficient m/ wt / pv power factor limits for WTG/PVG control c c ratios of existent/adding compressors ij / ij, k transmission coefficients through pipelines ij
Indices
i, j ij, jk k /t /s
indices for nodes indices for branches index for types of each component index for time periods/stages/scenarios
Sets set of existent nodes/pipelines/source nodes in the gas network set of existent nodes connected with energy hub HUB E _B / Lori / E _Sup set of existing fixed buses/feeders/electricity source buses G _B / E _B set of newly-added gas/distribution network nodes WTG / PVG / CB / SVC / VR set of candidate buses for WTG/PVG/ CB/ SVC/VR Lnew / Lup set of candidate branches for addition/replacement of feeders Pipe / C set of candidate branches for adding pipelines/compressors GS / HS / ESS set of candidate nodes for gas storage units/thermal storage units/electric storage systems set of candidate nodes for thermal storage units HS G _B / Pipe set of all nodes/branches in the gas network, G _B = { G _B , G _B } , Pipe = { Pipe , Pipe } E _B / Line / L set of all buses/branches/load nodes in the distribution network, E _B = { E _B , E _B } , Line = { Lori , Lup , Lnew } KPipe/ K C set of alternative types of adding pipelines/compressors K GS /KHS /KESS set of alternative types of adding gas storage units/ thermal storage units/electric storage systems KLup/ KLnew set of alternative types for replacing/adding feeders KWTG / KPVG/ K CB /K SVC / KVR set of alternative types for WTG/PVG/ CB/SVC/VR JK (j )/ IJ (j) set of feeders whose parent/child is j G _B / Pipe / G _Sup
energy conversion efficiency of P2G/CHP/GF storage/release efficiency of gas storage units sh / rh heat storage/release efficiency of thermal storage units se / re charge/discharge efficiency of ESS DG DG max / min maximum/minimum allowable DG permeability Cur _wt Cur _pv Cur _L / max / max maximum curtailment rate of wind/solar max power/electrical load kijvr incremental ratio between VR taps Rijvr maximum gear of VR j,max / j,min upper/lower limits for gas nodal PRESSURE C FijPipe ,max / Fij,max flow limits of existing pipes/compressors C GS FijPipe / F , k ,max ij, k ,max / Fk,max flow limits of adding pipes/compressors/gas storage units Sup F jSup ,max / F j,min upper/lower flow limit for gas source node j Sup Sup Pgas max / Pgas max upper/lower power limit for the sum of gas supply Sup Sup Sup Sup Pelec max / Pelec min, Qelec max /Qelec min upper/lower power limit for the sum of electricity supply Lori Lori Lori PijLori ,max / Qij,max , Pij,min / Qij,min maximum and minimum permissible active/reactive power flow through existent fix feeders Lup Lup Lup PijLup , k,max / Qij, k ,max , Pij, k ,min / Qij, k,min maximum and minimum permissible active/reactive power flow through replacing feeders Lnew Lnew Lnew PijLnew , k,max / Qij, k ,max , Pij, k ,min / Qij, k,min maximum and minimum permissible active/reactive power flow through adding feeders QkCB / QkSVC reactive power compensated per CB unit/SVC P jL,max / P jL,min maximum/minimum electrical load of node j ES PkHS ,max / Pk,max power flow limits of adding thermal storage units/ESS HS ES SkGS / S ,max k ,max / Sk,max storage capacity limits of adding gas storage units/thermal storage units/ESS Iij,max maximum allowable current amplitude Vmax /Vmin upper/lower limits of nodal voltage amplitude O _CB O _VR Nmax /Nmax maximum switching time for CB/SVC pv N jwt ,max /Nj,max maximum installation number of WTG/PVG NS number of scenarios NT number of time stages in a scenario Ny number of years during the whole planning horizon M very large number p2g chp gf eg / gh / gh sg / rg gas
Parameters duration time of time staget annual interest rate (s ) probability of scenarios lijpnew / lijlup/ lijlnew length of adding pipeline/replacing feeder/adding feeder rij /x ij resistance/reactance of feederij ckI _pipe/ ckI _c investment cost coefficients of pipelines/compressors ckI _gs/ckI _hs/ckI _ess investment cost coefficients of gas storage units/ thermal storage units/electric storage systems ckI _lup/ ckI _lnew investment cost coefficients of replacing/adding feeders ckI _wt /ckI _pv /ckI _cb/ ckI _svc / ckI _vr investment cost coefficients of WTG/ PVG/ESS/CB/SVC/VR c O _gs /c O _hs/ c O _ess operation cost coefficients of gas storage units/ thermal storage units/electric storage systems c O _p2g /c O _chp/ cO _gf operation cost coefficients of P2G/CHP/GF c O _loss operation cost coefficients of electrical loss c O _wt / c O _pv / operation cost coefficients of WTG/PVG /CB/SVC/ c O _cb/ c O _svc / cO _vr ESS/VR c Cur _wt / c Cur _pv / c Cur _L operation cost coefficients of curtailing wind power/solar power/electrical load ckM _pipe/ ckM _c maintenance cost coefficients of pipelines/compressors/gas storage units ckM _gs/ ckM _hs / ckM _ess Maintenance cost coefficients of gas storage
t
d
Constants fp
flow-power conversion coefficient
Variables gs c x ijpipe , k / x ij, k / x j, k binary investment variables for pipelines/compressors/gas storage units x jgs, k / x jhs, k / x jess , k binary investment variables for gas storage units/ thermal storage units/electric storage systems lnew x ijlup , k / x ij, k binary investment variables for replacing feeders/adding feeders/SVC/VR vr x jcb, k /x jsvc binary investment variables for CB/SVC/VR , k /x ij pv cb njwt , k / nj, k /nj, k integer investment variables for WTG/PVG/CB yjcb, t , s operation number of CB units yijvr, t , s operation tap position of VR rijvr, t , s operation ratios of VR
2
Applied Energy 253 (2019) 113517
J. Wang, et al.
nij+, tvr, s /nij, tvr, s , n+j, tcb, s / nj, tcb, s positive/negative regulation times of CB/SVC ujsg, t , s / urg j, t , s operation state of storage/release gas for gas storage units ujsh, t , s / urh j, t , s operation state of storage/release heat for thermal storage units ujse, t , s / ure j, t , s operation state of charging/discharging for ESSs gas nodal pressure j, t , s f jSup gas supply from source nodes ,t,s fij, t , s / f jk, t , s flows through pipelines Sup pjSup , t , s / qj , t , s active/reactive power supply from source node j
pjp, t2,gs
power flow into P2G from distribution network
pjLh ,t,s
thermal load power
WTG/PVG/load pj,Lt , s /qj,Lt , s active/reactive load after curtailment f jsg, t , s /f jrg, t , s gas flow into/out of gas storage units
rh pjsh , t , s / pj, t , s power flow into/out of thermal storage units pjse, t , s / pjre, t , s charging/discharging power of ESS
HS ESS S GS j, t , s / Sj, t , s /Sj, t , s storage capacity of gas storage units/ thermal storage units/ESS Iij, t , s current flow amplitude Vj, t , s nodal voltage amplitude Vm, t , s voltage amplitude of dummy bus in VR CTotal total cost CInv total cost of investment COpe /CMa/CSup annual cost of operation/maintenance/energy purchase CInv _Gas/CInv _Elec /CInv _Hub cost of investment for NGN/ADN/EH _Gas _Hub _Elec CsOpe /CsOpe /CsOpe cost of NGN/EH/ADN operation at time ,t ,t ,t stage t under scenario s CMa _Gas /CMa _Hub/CMa _Elec Annual NGN/EH/ADN maintenance cost _Gas _Elec CsSup /CsSup cost of gas/electricity purchase at time stage t ,t ,t under scenario s
gf f jchp , t , s / f j, t , s flow into CHP/GF from gas network pij, t , s / pjk, t , s , qij, t , s / qjk, t , s active and reactive power flows through feeders qjcb, t , s / qjsvc , t , s reactive power compensated by CB/SVC
f jL, t , s
gas load flow
pjL, t , s / qjL, t , s active/reactive load before curtailment pv wt pv pjwt , t , s /pj, t , s , qj, t , s / qj, t , s active and reactive power of WTG/PVG generation before curtailment _wt _pv _wt pjcur / pjcur , qjcur , t,s ,t,s ,t ,s / curtailment of active and reactive power for _pv cur _L _L qjcur /qjcur , t , s , pj, t ,t
management (DM), on load tap changers (OLTCs), and static var compensations (SVCs) was fully examined. Authors in [15] proposed a hierarchical robust expansion planning model integrating active management elements that includes OLTCs, SVCs, capacitor banks (CBs) and electric storage systems (ESSs). In [16], a flexible expansion planning model was proposed for ADN. In this work, DGs and ESSs were considered as capable expansion options, and in order to enhance the system’s reliability, island mode operation on the network was also incorporated in ANM schemes. In [17], the power electric vehicle and the solar parking lots can be applied as the sources of active and reactive power. Thus, a planning problem combing the economic aspects and technical factors of the system is proposed to achieve realistic results. Authors in [18] discussed the DG planning in ADN with the consideration of voltage regulators (VRs) and demand management. In [19], authors have developed an ADN planning model that focused on efficient renewable energy harvesting. In the ADN planning modeling of [20], authors mainly considered the dispatch of centralized and distributed ESSs. In [21], authors proposed an ADN planning with consideration of DM as well as the control of reactive compensators, OLTCs, ESSs and DGs. In the abovementioned literature, different parts of ANM schemes were modeled and examined selectively, mainly including adjustment of voltage regulators (VRs), reactive compensation of CBs and SVCs, dispatch of ESSs and DGs, demand management. However, there still remains a research gap in expansion planning model with full consideration of ANM schemes aforementioned. Another crucial challenge for planning is the method to address the uncertain source injection. The variability and uncertainty issues of renewable energy sources (RESs) have been well addressed in some research. In [22], a stochastic framework is adopted to solve the distribution system reconfiguration problem while adaptively and dynamically dealing with the incidental charging pattern of power electric vehicles and the variable and uncertain power of RESs. The reconfiguration of electrical distribution model in [23] is the first work to investigate the dynamic and adaptive reconfiguration problem, where the uncertainty of RESs at every time step (hour) of the day for the updated optimization time horizon is considered. In the same way, a model in [24] employed predictive control technique to minimize the energy loss with consideration of the uncertain RESs. Although these advanced dynamic techniques for RESs have been applied to model the reconfiguration problems in many studies, a long-term expansion planning issue requires a lager sample than a 24-h period data to cover
(MESEP) has attracted more attention than ever before. This issue includes the determination of the options for investment of renewable energy based distributed generation (DG) and the new infrastructures to be added to the existing system, and the size, location and type of them to be installed to meet the increasing energy demand over the planning horizon [6]. However, the expansion planning of MES is still confronted with a lot of huge challenges in planning and modeling approaches. Research in [7] presented a two-stage mixed-integer linear programming (MILP) model for MES planning to optimize the MES configuration from scratch considering distributed renewable energy source based on the EH concept. The proposed approach enabled planners to optimize the equipment investment and the MES configuration in a coordinated manner. In [8], simultaneous expansion planning of electricity and heat sources along with the transmission network was considered together. By using linear models for energy transmission network, the expansion problem was converted to a MILP problem to be solved efficiently. Authors in [9] proposed a new expansion planning formulation, which sought out the optimum usage of energy hubs for minimizing the utility costs, and at the same time maximizing the system’s reliability. These existing research works have made great contributions to planning the energy system based on the concept of energy hub; meanwhile, all the aforementioned works investigated the energy system expansion planning problems focusing on the part of the traditional power system. The electrical distribution planning problem has been well modeled in some studies like [10] and [11]. In [10], distributed generation placement planning problem has been modeled considering reliability level and power loss. In [11], the optimal placements of DGs and capacitors were investigated simultaneously in distribution network. However, due to renewable energy policy-making and its substantial advantages to accommodate the widespread renewable energy sources, the traditional distribution network is undergoing the transformation from a passive pattern to another, which is called active distribution network (ADN) [12]. The emerging active network management (ANM) technologies integrating into ADN make the network more flexible and controllable [13]. This transition brings huge challenges to distribution planners and operators. To model this issue, there are many research works concentrating on ADN expansion planning. In [14], a novel expansion planning model of ADN was proposed. In this work, the coordinated optimization of ANMs including the controls of DGs, demand 3
Applied Energy 253 (2019) 113517
J. Wang, et al.
more possible operating conditions. To the best of our knowledge, there still remains research gaps from the following aspects:
contains not only normal but also extreme conditions, which is more suitable for long term planning of larger systems. Due to the fact that mixed-integer nonlinear nature of the proposed expansion planning model makes the model particularly difficult to be solved, several equivalent reformulation methods are employed to linearize the model to be convex in order to ensure the global optimization results. Our linearization approach relies on the application of (i) second order cone relaxation (SOCR) method, which converts the original model into a mixed integer second order cone programming (MISOCP) model [27]; (ii) a modified piecewise linearization approach based on the work in [28], which can address the nonlinear polynomial expression efficiently; (iii) big M method, which can cope with the bilinear terms particularly; and (iv) auxiliary variable approach, which can replace the nonlinear product terms like quadratic terms and absolute value expressions. Finally, a modified IEEE 33-node distribution system and a new 23-node gas system coupled by 8 energy hub units are employed to verify the effectiveness of the proposed model. The main contributions of this paper can be summarized as follows:
(1) The comprehensive expansion planning model of multi-energy system where active distribution network is incorporated as an electricity part has never been reported in literatures. (2) The impact of ANM schemes on multi-energy system expansion planning has not been explored so far, and has not been applied to model the MES planning problem. To fill the knowledge gap of the state-of-the-art research, a new expansion planning model of multi-energy system is originally proposed based on the energy hub concept. In this model, the considered MES is composed of active distribution network, gas network, and energy hub unit. As a comprehensive expansion planning model, this paper centers on the optimal determination of the type, location and size of the feeders, DGs, CBs, SVCs, VRs and ESSs in ADN, natural gas pipelines, compressors and gas storage units (GSU) in NGN, and thermal storage units (TSUs) in EH. With the aim of minimizing the total cost over the planning horizon, the final investment decision scheme and the optimal operation strategy (ANM schemes, EH management, gas flow dispatch, etc.) for the expansion of MES are obtained. It is worth noting that, unlike the existing literatures, aforementioned ANMs including adjustment of VRs, reactive compensation of CBs and SVCs, dispatch of ESSs and DGs, as well as DM are fully modeled and tested in this paper. In addition, in order to consider the robustness [25] of the planning results, the stochastic framework is adopted and a robust probabilistic scenario generation approach extended from [26] based on the corresponding annual historical data is innovatively proposed to address the uncertainty of RESs, as well as the uncertain demands (electric, gas and heat demands). The resulting scenario set
(1) An expansion planning model of MES with the consideration of integrating ADN is originally proposed, in which the full ANM schemes are modeled and incorporated. (2) A novel scenario generation approach considering the robustness of planning results is designed for testing the proposed model. (3) A modified piecewise linearization as well as the SOCR approach is applied to convert the original nonlinear model into a mixed-integer SOCP model, and hence the model can be efficiently solved. The remainder of this paper is structured as follows: the MES is modeled in Section 2. In Section 3, a probabilistic scenario generation approach is proposed for testing the proposed model. In Section 4, the MES expansion planning problems integrated with ADN, NGN and EH
Substation Active Distribution Network
Natural Gas Network
Thermal Storage Unit
Gas Storage Unit
Gas Furnace
Combined Heat and Power
Power to Gas
Distribute Generation
Energy Hub
Gas Load
Electrical Energy Storage
Thermal Load
Energy Demand Fig. 1. Multi-energy system framework with energy hub. 4
Electrical Load
Applied Energy 253 (2019) 113517
J. Wang, et al.
power and solar power) throughout a year. Each of them is divided by the corresponding maximum value to obtain the per-unit value, and hence total 8760 hourly scenarios with normalized data is prepared. Step 2: Sort out the largest and smallest data for each hour through the year separately. Meanwhile, the mean value for each hour is computed. Therefore, three 24-h curves (i.e. the maximum, minimum and mean value curves) for each of them can be obtained. Step 3: Generate the scenario set based on three curves. The maximum value curves of wind power and solar power, as well as the minimum value curves of all demands are put together to constitute extreme scenario A . Meanwhile, extreme scenario B is generated by combining the maximum value curves of all demands and the minimum value curves of wind and solar power. Finally, scenario C is the combination of all mean values, expressed as vt , t = {1, 2, , 24} . Step 4: Calculate the probability for each scenario. Due to the data have been normalized in step 1, each value of the data is between 0 and 1. By introducing a fluctuation coefficient , data can be classified into ¯ t 1} , low level three intervals: the high level { ¯ t : vt + { : 0 vt } , and middle level { t : vt vt + } . In t t t detail, for a specific hour, if the demand is in the low level and energy power is in the high level, this hour will belong to scenario A . The reverse situation will belong to scenario B , and the rest belong to scenario C . Therefore, the probability of each scenario situation can be obtained by counting their frequency throughout the data (i.e. (s ) is calculated by dividing the total number of the data belonging to scenario s by the total number of the data in the year, s { A , B , C} and s (s ) = 1). In this paper, the proposed generation procedure of the extreme scenarios is implemented based on the historical data (download from http://www.elia.be/en/grid-data) of generation and demand for 2016. Accordingly, the maximum value curve and the minimum value curve of wind power, photovoltaic power, electricity demand, gas demand and heat demand can be obtained respectively, as shown in Figs. 2 and 3. In addition, the average values curves are shown and denoted by dotted line in each figure. For clarity, we give three points A, B and C at time t in the figures as examples. Each point includes five values (e.g. A = [A1, A2, A3, A4, A5]) representing the wind power, photovoltaic power, electricity, gas and heat demands in one operating condition (A) respectively. If the operating condition at time t is A, this scenario will be classified into scenario A since the generation power output A1 and A2 are in the high level while demands A3, A4 and A5 belong to lower level. In the same way, B and C will be classified into scenarios B and C , respectively.
units are modeled and formulated. The linearization methods, mainly including modified piecewise linearization and SOCR approach are described in Section 5. The flowchart of the overall proposed modeling methodology is given in Section 6. Numerical results and analysis based on a multi-energy system including a modified IEEE 33-node distribution system, a new 23-node gas system and 8 energy hub units are presented in Section 7. Finally, conclusions are drawn in Section 8. 2. Multi-energy system In this section, the framework of MES is modeled, as shown in Fig. 1. In MES, converters, storage and energy transfer equipment are used to combine and share the energy carriers under the so-called EH units. In this way, the MES is represented as an integrated system which is mainly composed of natural gas network, active distribution network and energy hubs that connect two networks. Furthermore, a EH unit is constituted by power to gas, combined heat and power, thermal storage unit and gas furnace (GF) to make its function possible. In the output port of EH, energy demand is provided by sharing and converting electricity and natural gas. In particular, the input port of EH is linked to energy networks such as ADN and NGN and the input electricity and gas are purchased from the upper distribution network or transformer station and natural gas distribution network, respectively; the output port of EH unit is connected to the loads or the energy networks. In this paper, thermal load is considered to be located in the output port of EH. In this regard, the conversions between natural gas and electricity as well as the supply of the thermal demand can be achieved and satisfied in the proposed approach. 3. Probabilistic scenario generation The key of probabilistic programming is to formulate the decisionmaking scheme under the combination of potential circumstances [29]. In order to obtain the robustness of investments, this paper tries to find the optimal planning scheme based on probabilistic scenarios which includes both extreme and typical situations. Given the fact that wind speed, solar radiation intensity and load demands all change constantly with time, it is likely to generate the scenarios based on the historical hourly data. In this paper, 24-h data through a year are used as the research objects of scenarios, and thus main characteristics of time sequence of wind power, photovoltaic power generation and load can also be completely reflected. Based on this point, the probabilistic scenario set is generated by the following steps. Step 1: Firstly, normalize the historical data of demands (including electricity, gas and heat demands) and renewable energy (i.e., wind
Fig. 2. Renewable energy power curves of 24-h under 3 scenarios (including the maximum value and the minimum value curves, as well as the mean value curves of the wind and photovoltaic power). 5
Applied Energy 253 (2019) 113517
1 B3 C3 A3
Gap2
0.5
Min value curve Max value curve Mean value curve
0
6
1
Gas flow (pu)
Gap1
v3 t +Δ v3t v3 t - Δ
B4 C4
A4
0.5
0
12
Time (h) (a) Electricity demand
6
12
18 Thermal power (pu)
Demand power (pu)
J. Wang, et al.
v4 t + Δ
v4t v4 t - Δ
18
1
B5 v + Δ 5t
C5 v5t A5 v5 t - Δ
0.5
0
24
Time (h) (b) Natual gas demand
24
6
12
Time (h) (c) Heat demand
18
24
Fig. 3. Demand power curves of 24-h under 3 scenarios (including the maximum value, the minimum value and the mean value curves for electrical, gas and thermal demands).
4. Modeling of the multi-energy system expansion planning
leverage the investment and operation costs.
Based on the framework of MES and the approach of probabilistic scenario generation, a model of expansion planning for MES can be mathematically formulated using the approach of scenario-based deterministic equivalent, as described as follows. To increase the readability of this passage, some brief descriptions about key symbols which are easy to be confused or used repeatedly are given here: t refers to the time stage; s is the index of scenario; i and ij denote the identifier of node and pipe/feeder respectively; k is the index of type of investment equipment. Omitting the superscripts and subscripts of the notations for simplification, c represents the unit cost of investment or operation; Γ, Ω and Λ are the collections of the existing, candidate and all locations of the investment equipment respectively, and = ; K is the collection of alternative types; x denotes the corresponding binary investment variable; n is the integer investment variable; y is the integer operating variable of infrastructures; parameter l denotes the length of pipelines or feeders; f / F denotes the variable/parameter related to the gas flow; represents the nodal pressure in the gas network; p / P and q/Q are the variables/parameters representing the corresponding active and reactive power. N denotes the number of corresponding things (such as years and facilities). In addition, other notations in this paper are fully explained in Nomenclature part.
4.1.1. Investment cost In this paper, the expected investment cost is firstly formulated as: It can be seen that the investments in the proposed model are constituted by three items: (i) the cost in gas network with consideration of installing pipelines, compressors and gas storage units (Eq. (3)); (ii) the cost in energy hubs where the installation of thermal storage unit is involved (Eq. (4)); and (iii) the cost in active distribution network considering the upgrade and expansion of feeders, as well as the installations of wind turbine generations (WTGs), photovoltaic generations (PVGs), ESSs, CBs, SVCs and VRs (Eq. (5)).
CInv _Gas =
=1
Pipe k KPipe
j
GS k KGS
ckI _pipe lijpnew x ijpipe ,k +
(COpe + CMa + CSup) 1
ckI _c x ijc, k + (3)
Lup k KLup
ckI _lup lijlup xijlup ,k +
j
+ PVG k KPVG
WTG k KWTG
ckI _pv njpv ,k +
+ SVC k K SVC
ij
Lnew k KLnew
ckI _lnew lijlnew x ijlnew ,k
ckI _wt njwt ,k
+
j
(4)
HS k KHS
ij
1 (1 + d )
C k KC
ckI _gs x jgs, k
CInv _Elec =
j
ij
ckI _hs x jhs, k j
The objective (1) is to minimize the expected total cost in the planning horizon, which consists of 4 terms: (i) investment cost (Eqs. (2)–(5)), (ii) operation cost (Eqs. (6)–(9)), (iii) maintenance cost (Eqs. (10)–(13)), and (iv) energy purchase cost (Eqs. (14)–(16)) in the multienergy system. In other words, based on the investment and construction, the corresponding operation cost, maintenance cost and energy purchase cost are taken into account. Ny
ij
CInv _Hub =
4.1. Objective function
min CTotal = CInv +
(2)
CInv = CInv _Gas + CInv _Hub + CInv _Elec
j
ESS k KESS
ckI _ess x jess ,k +
j
CB k KCB
ckI _cb njcb, k
ckI _svc x jsvc ,k + ij
VR k KVR
ckI _vr x jvr, k (5)
Eqs. (3)–(5) is written in terms of the variables that represent the investment decisions. For example, the cost of the newly-added pipelines depends on the binary investment variable of newly-added pipe ij I _pipe ( x ijpipe ) and the length , k ), the cost of the k -type pipe per unit length (ck pnew of each newly-added pipeline (lij ); the investment cost of
(1)
where Ny represents number of years during the whole planning horizon and d represents annual discount rate. Herein, the economic factor (i.e., discount rate d [30]) is considered in the planning problem to 6
Applied Energy 253 (2019) 113517
J. Wang, et al.
compressors is calculated according to the binary variable of whether to use k -type compressor ( x ijc, k ) and its price ckI _c ; the investment cost of DG pv is calculated by the number of the k -type WTG/PVG (njwt , k /nj, k ) in-
Eq. (13) represents the additional maintenance costs brought by upgraded lines and the maintenance costs corresponding to newly-added lines and newly-installed WTGs, PVGs, ESSs, CBs, SVCs and VRs in active distribution network.
stallation and its price (ckI _wt /ckI _pv ). The other investment costs are similar.
CMa _Gas =
4.1.2. Operation cost The operation cost is given as: NT
COpe = 365
_Gas _Hub _Elec (CsOpe + CsOpe + CsOpe ) ,t ,t ,t
(s ) s
(6)
t=1
c O _gs (f jsg, t , s + f jrg, t , s )
j
fp
(cO _p2gpjp, t2,gs + c O _chpf jchp ,t ,s
j
fp
+ cO _gf f jgf, t , s
ckM _hs x jhs, k
ij
_Elec CsOpe = ,t
ij
j
c O _vr |yijvr, t , s
WTG k KWTG
j
+
ij
1, s |
c O _cb |yjcb, t , s
+ j
yjcb, t
CB k K CB
ESS k KESS
ckM _cb njcb, k +
c O _svcqjsvc ,t,s +
j
NT
CSup = 365
(s ) s
c O _ess (pjse, t , s + pjre, t , s ) +
j
j
PVG
j
L
_L c Cur _Lpjcur + ,t, s
+
j
c O _wtpjwt ,t,s + WTG
j
1, s |
_wt c Cur _wtpjcur , t,s WTG
cO _pvpjpv ,t ,s PVG
(9) 4.1.3. Maintenance cost
CMa = CMa _Gas + CMa _Hub + CMa _Elec
PVG k KPVG
ckM _pv njpv ,k
ckM _svc x jsvc ,k +
ij
VR k KVR
ckM _vr x ijvr, k
_Gas _Elec (CsSup + CsSup ) ,t ,t
(14)
t,
s
(15)
t,
s
(16)
where the annual energy purchase cost (CSup ) is calculated by the similar way to the annual operation cost (COpe ) in terms of the method of _gas _elec probabilistic scenario generation. ptsup and ptsup indicate the total ,s ,s gas and electricity supply power at time stage t of scenario s, (i.e., sup _elec _gas = j E _Sup pjSup ptsup = j G _Sup f jSup , t , s respectively). ,s , t , s fp and pt , s In addition, the real time pricing options will gradually become a trend in energy trade in the future. In order to optimize energy utilizing from the perspective of the MES, the real-time change of the supply elec prices of gas and electricity (ctgas , s and ct , s ) are considered in the model, as done in [31]. That is, the MES is allowed to access to the real time price data and then it can wisely manage its energy consumption according to them. By this means, the energy converters (i.e., CHP and P2G) of EHs is controlled and operated according to weighing the cost of using different forms of energy. If the gas price is low while electricity price is in the high level, CHP will generate electricity using natural gas to support the supply of electricity. On the contrary, if the gas price is high while electricity price is in the low level, P2G will work and convert the excess electrical energy into the natural gas form. It is assumed that the gas and electricity price are determined by the _gas _elec demands ( ptsup and ptsup ) in real time. ,s ,s
_pv c Cur _pvpjcur ,t, s
+
j
ckM _ess x jes, k
SVC k K SVC
t=1
_Gas sup _gas CsSup = ctgas t, ,t , s pt , s
CB
ESS
j
ckM _lnew lijlnew x ijlnew ,k
4.1.4. Energy purchase cost The energy purchase cost is formulated as Eq. (14) including gas purchase and electricity purchase and the energy costs associated with them are formulated as (15) and (16) respectively.
Line
SVC
ckM _wt njwt ,k +
c M _lori ) lijlup x ijlup ,k
(13)
c O _lossIij2, t , s rij
+ t +
yijvr, t
(ckM _lup
Lnew k KLnew
+
(8)
VR
(12)
Lup k KLup
+
fp ) +
c O _hs (pjsh, t , s + pjrh, t , s )
ckM _c x ijc, k +
HS k KHS
ij
j
C k KC
(11)
+
Hub
Hs
ij
ckM _gs x jgs, k
_Elec sup _elec CsSup = ctelec t, ,t , s pt , s
= t
j
GS k KGS
ckM _pipe lijpnew x ijpipe + ,k
CMa _Elec =
_Hub CsOpe ,t
j
j
j
(7)
GS
Pipe k KPipe
CMa _Hub =
It is worth to be noted that the operational cost in (6) is computed based on the approach of probabilistic scenario generation. By this means, there are 365· (s ) days belonging to scenario s and NT time stages in each scenario. Thus, the annual operation cost (COpe ) is calculated by the sum of the total operation cost at each time stage t under Ope each scenario s (i.e., CsOpe , t ) through the year and Cs, t is equal to the sum of operation costs in the gas network, distribution network and energy Ope _Gas _Hub _Elec + CsOpe + CsOpe ). hubs (i.e., CsOpe , t = Cs, t ,t ,t Ope _Gas C In detail, s, t is characterized in (7) where the operation of gas storage units is included; Eq. (8) models the operation costs related to _Hub energy hub (CsOpe ) consisting of P2Gs, CHPs, GFs and thermal sto,t rage units; In (9), not only power losses and operations of VRs, CBs, SVCs, EESs, WTGs and PVGs but also active management measures such as curtailments of wind power, solar power and electrical load are taken into account. In general, the operational cost for each component _Gas _Elec _Hub (CsOpe /CsOpe /CsOpe ) can be calculated by the production of the ,t ,t ,t corresponding power (such as f jsg, t , s fp / pjp, t2,gs ), the unit cost (such as c O _gs /c O _p2g ) and the duration ( t ). Considering the artificial adjustment expenses (for VRs or CBs), the regulation times in the duration (such as |yijvr, t , s yijvr, t 1, s |), multiplied by the relevant unit cost (such as c O _vr ) is included in the operation cost of ADN. _Gas CsOpe = t ,t
ij
(1) Real-time gas price
(10)
In a similar way to the formulation of the investment cost, the maintenance cost is formulated as Eq. (11)–(13). The maintenance costs of newly-added pipelines, compressors and gas storage units in the gas network are modeled in (11); in (12), the maintenance costs of newlyinstalled thermal storage units are taken into account in energy hubs;
It is assumed in (17) that the relationship between supply price and demand of natural gas can be described by a monotonic incremental linear function.
ctgas ,s = 7
sup _gas 1 pt , s
+
0
(17)
Applied Energy 253 (2019) 113517
J. Wang, et al.
(2) Real-time electricity price
x are established in (21)–(30). This set of constraints indicates that a maximum of one replacement, addition, or installation is permitted for each component and location as well. Expressly, the connectivity and the radiation of the expanded integrated network are ensured by constraints (21) and (22). Particularly for a newly-added bus j ( j E _B , and E _B denotes the collection of the newly-added buses), it is restricted to connect to one and only one existent bus (the index i is to represent the candidate existent buses). Meanwhile, a maximum pv number of installation for WTG/PVG (N jwt ,max /Nj,max ) at each candidate location is imposed in (31)–(32). (33) limits the maximum permitted installation number of capacitors in each CB unit (N CB j,max ).
The electricity cost of upper level substation can be expressed as follows: _elec ce (ptsup )= ,s
sup _elec m ) m (pt , s
_elec + ...+ 1 ptsup + ,s
(18a)
0
In (18b), it is supposed that the upper-level substation sets the supply price according to the maximization of electricity sale income. NT
max s
t=1
_elec elec (ptsup ct , s ,s
_elec ce (ptsup )) ,s
(18b)
The solution of the optimization problem can be obtained as follows:
4.3. Operational constraints
sup _elec ctelec ) , s = ce (pt , s
4.3.1. Constraints related to the natural gas network
(18c)
To simplify the problem, set m = 2: _elec ce (ptsup ) ,s
ctelec ,s
=2
=
sup _elec 2 ) 2 (pt , s
sup _elec 2 pt , s
+
sup _elec 1 pt , s
+
(1) Gas supply constraints
+
(18d)
0
0
(18e)
1
sup _gas 0 ) pt , s
_Elec _elec CsSup = (2 2 ptsup + ,t ,s
t,
sup _elec 1) pt , s
t, t,
s t,
j,min
(20)
where and
j
G _B ,
lnew x ijlnew , k = 1, x ij, k = \{ 0, 1\} , i
j
E _B , i
k KLnew
x ijc, k
1, x ijc, k = \{ 0, 1\} ,
ij
ij
Pipe
x jgs, k
1,
x jgs, k
= \{ 0, 1\} ,
j
1, x jhs, k = \{ 0, 1\} ,
GS
j
1, x ijlup , k = \{ 0, 1\} ,
HS
x jess ,k
1, x jess , k = \{ 0, 1\} ,
x jcb, k
1, x jcb, k = \{ 0, 1\} ,
ij
ESS
CB
j
SVC
1, xijvr = \{ 0, 1\} ,ij njwt ,k
wt N jwt ,max , nj, k
Z+,
j
WTG
njpv ,k
pv Njpv ,max , nj, k
Z+,
j
PVG
j
CB ,
(35)
s
ij IJ (j)
f jL, t , s + f jSup + f jrg, t , s ,t ,s
fij, t , s
f jsg, t , s
f jgf, t , s
f jchp , t,s
lg Lg j, t , s P j
j t,
s
K CB
sgn(
2 i, t , s
2 2 j, t , s ) fij, t , s
= =
ij (
2 i, t , s
2 j, t , s ),
ij (
2 i, t , s
2 j, t , s )
ij
Pipe /
C,
x ijpipe ,k ,
ij
t,
s Pipe ,
(37)
t,
s
ij (
2 i, t , s
2 j, t , s )
ij (
2 i, t , s
2 j, t , s )
t,
s
sgn(
2 i, t , s
2 2 j, t , s ) fij, t , s
+M
sgn(
2 i, t , s
2 2 j, t , s ) fij, t , s
M
k KC k KC
x ijc, k xijc, k
,
ij (39)
by the maximum gas demand of node j (P jLg ) and the gas load rate at
time t of scenario s ( j, tlg, s ) and the rate can be obtained by the approach of probabilistic scenario generation. Additionally, the Weymouth steady state power flow constraints [32] are formulated as expression (37)–(39) where ij represents the transmission coefficient through pipeline ij and sgn( i2, t , s - j2, t , s ) is a sign function (for sgn(a) , if a > 0 , the value of sgn(a) is 1; if a = 0 , the value of sgn(a) is 0 and otherwise, the value of sgn(a) is −1). It should be noted that the flow through pipeline ij ( fij, t , s ) is imposed
(32)
k
sgn(
2 2 j, t , s ) fij, t , s
Gas flow balance is formulated by Eq. (36), which indicates that the sum of gas inflows is equal to the sum of gas outflows in any node at any time stage and any scenario. Note that the gas load ( f jL, t , s ) is computed
(31)
k KPV
Z +,
f jk, t , s =
2 i, t , s
C,
(30)
VR
k KWT
cb cb N CB j,max x j, k , nj, k
t,
(38)
(29)
k K SVC
x ijvr
G _B ,
k KPipe
(28)
k KCB
1, x jsvc , k = \{ 0, 1\} ,
j
represents nodal pressure at time stage t of scenario s; are the lower and upper pressure limits.
G _B ,
(27)
k KESS
x jsvc ,k
j, t , s
j,max
(22)
(26)
k KLup
j
j,max ,
f jL, t , s =
(25) Lup
j
(34)
+ f jp, t2,gs
(21)
(24)
k KHS
xijlup ,k
Lnew
jk JK (j)
(23)
k KGS
x jhs, k
ij
s
(36)
C
k KC
njcb, k
i
k KPipe
t,
(3) Gas flow constraints
Constraints related to the investment and utilization decisions of all facilities in this planning are given as: i
G _Sup ,
denotes the gas supplied by the gas source node j at time
j, t , s
j,min
4.2. Investment and utilization constraints
pipe x ijpipe , k = 1, x ij, k = \{ 0, 1\} ,
where
f jSup ,t,s
j
(2) Pressure inequality constraints
(19)
s
F jSup ,max ,
stage t of scenario s and F jSup ,max refers to the corresponding maximum value.
Therefore, Eqs. (15) and (16) can be transformed into (19) and (20), respectively. _Gas _gas CsSup = ( 1 ptsup + ,t ,s
f jSup ,t ,s
to be 0 if no pipeline is added (i.e. k KPipe x ijpipe , k = 0 ) in (38). Constraint (39) guarantees that compressors can be added in pipelines that have been previously built. Given that M is a very large number, (39) will be
(33)
The relevant constraints in terms of the binary investment decision 8
Applied Energy 253 (2019) 113517
J. Wang, et al.
equivalent to (37) if not to add the compressor ( restricted otherwise.
k KC
for each gas storage unit, its initial capacity (S GS j,0, s ) is imposed to be 0.
x ijc, k = 0 ), un-
4.3.2. Constraints relevant to energy hubs
(4) Pipeline operation constraints
FijPipe ,max
FijPipe ,max ,
fij, t , s
pipe FijPipe , k ,max x ij, k
ij
Pipe /
C,
t,
s
pipe FijPipe , k,max x ij, k ,
fij, t , s
k KPipe
(1) Thermal power balance constraints
(40)
ij
Pipe ,
t,
sh phgf j, t , s + phchpj, t , s = pjLh , t , s + pj, t , s
s
pjLh ,t ,s =
k KPipe
(41)
FijPipe ,max
x ijc, k
M
FijPipe ,max + M
fij, t , s
k KC
x ijc, k ,
ij
C,
t,
k KC
j, t , s
FijC,max for
s:
c ij i, t , s
=
(43)
FijC,max
fij, t , s
ij
C,
c ij, k i, t , s
k
t,
x ijc, k )
M (1
FijC, k,max
(44)
KC,
xijc, k )
M (1
c ij, k i, t , s
x ijc, k )
+ M (1
FijC, k,max + M (1
fij, t , s
x ijc, k )
pjgf, t , s = f jgf, t , s
(45) (46)
GS ,
ujsg, t , s + urg j, t , s
f jsg, t , s
0
f jrg, t , s
S GS j, t , s =
0
t
gs rg FkGS ,max x j, k uj, t , s ,
t
(48b)
k K GS
0+
S GS j, t 1, s
+
t
f jsg, t , s sg
t
gs SkGS ,max x j, k ,
S GS j, t , s k
GS
t
f jrg, t , s / rg
t
t=1 t
j
HUB ,
t,
s
(51)
pechpj, t , s = f jchp ,t ,s
chp fp ge ,
f jchp , t,s
j
HUB ,
t,
s
chp Fmax
(52)
t
2
pjp, t2,gs
p2g eg / fp p2g Pmax
,
j
HUB ,
t,
s
(53)
Note that in (53) the gas output of P2G ( f jp, t2,gs ) is calculated in terms of the electric power injecting into P2G from distribution network ( pjp, t2,gs ), multiplied by the corresponding efficiency ( egp2g ) and divided by the flow-power conversion coefficient ( fp ). In (53) P2G is also imposed to transform from electric power into gas flow by the corresponding p2g non-negative constraint and Pmax is the upper power input limit for P2G.
(48a)
k K GS
f jrg, t , s / rg
,
chp fp gh
0
(47)
f jsg, t , s sg
gf Fmax
f jp, t2,gs = pjp, t2,gs
t
gs sg FkGS ,max x j, k uj, t , s ,
(50)
(4) Operational constraints for power to gas
ujrg,1, s = 0 0
gf fp gh
phchpj, t , s = f jchp ,t ,s 0
s:
1,
s
(3) Operational constraints for combined heat and power
(6) Gas storage unit operation constraints
j
f jgf, t , s
0
The operation of existing compressors ( ij C ) is modeled by (43)–(44) where ijc / ijc, k represents the booster coefficient of existent and newly-added compressors. On the other hand, the operation of those potential to be installed ( ij C ) are subject to (45) and (46), which are still associated with the investment decision on them ( x ijc, k ).
for
t,
(2) Operational constraints for gas furnace
s:
j, t , s
HUB ,
gf chp / gh ) and the flow-power ergy conversion efficiency of GF/CHP ( gh conversion coefficient fp . Constraints (51) and (52) also give the fact that the gas flow over a given range can input into GF and CHP. In addition, P2G is utilized in the energy hubs to convert electric power into natural gas.
(5) Compressor operation constraints
t,
j
computed according to the maximum load (P jLh ) and the corresponding
c also on the investment decision x ijpipe , k / x ij, k .
C,
,
rate ( j, tlg, s ). It also can be noted that the construction of thermal storage unit provides the possibility of storing heat that could not be consumed or generated in time. Based on the model of EH, the thermal power could be supplied by GF and CHP. In (51)/(52), the thermal power supplied by GF/CHP is calculated according to the gas flow injecting into GF/CHP from gas network ( f jgf, t , s / f jchp , t , s ), multiplied by relevant en-
Expression (40) sets the bounds on the gas flow in pipelines which has been built previously. In (41) and (42), it can be seen that the gas flow through the newly-added pipelines or candidate pipelines for Pipe compressors depends not only on the upper flow limit FijPipe , k ,max /Fij,max but
ij
lh Lh j, t , s P j
For energy hub nodes, the sum of thermal power outputs and the sum of thermal power inputs is equal, which is stated in Eq. (50). In the way that is similar to computing the gas load, the thermal load ( pjLh , t , s ) is
s (42)
for
pjrh, t , s
(5) Operation constraints for thermal storage units
(49a)
Exactly as gas storage units, thermal storage unit constraints can be modeled as follows:
(49b)
for
It is not allowed to storage and release gas simultaneously for each gas storage unit and it is also not allowed to release gas at the first stage (t = 1), which can be enforced by (47). Constraint (48) represents that each gas storage unit storages or releases gas only when the conditions sg and rg in conx gs = 1 and or urg j, t , s = 1 are both satisfied. k K GS j, k straint (49) denote the efficiencies of gas storage and release respectively. S GS j, t , s represents the capacity at time stage t of scenario s and it is associated with the capacity in the former time stage (S GS j, t 1, s ), the flow of gas storage/release in the current time stage ( f jsg, t , s / f jrg, t , s ), storage/ release efficiency of gas storage units ( sg / rg ) and the duration t . And
j
HS ,
ujsh, t , s + urh j, t , s ujrh,0, s 0
s:
1,
t
=0
pjsh ,t, s
(54) hs sh PkHS ,max x j, k uj, t , s ,
t
hs rh PkHS ,max x j, k uj, t , s ,
t
k KHS
0
pjrh, t , s k KHS
9
(55a) (55b)
Applied Energy 253 (2019) 113517
J. Wang, et al.
0 + pjsh, t , s
SjHS ,t ,s =
SjHS , t 1, s
sh
+
SjHS ,t, s
0
pjrh, t , s /
t
sh pjsh ,t, s
t
(4) Feeder operation constraints
t=1
pjrh, t , s / rh
t
hs SkHS ,max x j, k , k
rh
t
t
2
for
(56a)
t
(56b)
HS
for
4.3.3. Constraints related to the active distribution network
j
E _Sup ,
t,
Lori ,
t,
s:
PijLori ,max
pij, t , s
PijLori ,max
QijLori ,max
qij, t , s
QijLori ,max
ij
Lup ,
PijLori ,max (1
(1) Electricity supply constraints
for
ij
t,
k KLup
s: x ijlup ,k )
PijLori ,max (1
pjSup ,t, s
P jSup ,max
(57)
0
qjSup ,t, s
QjSup ,max
(58)
k KLup
QijLori ,max (1
k KLup
x ijlup ,k )
k KLup
for
(2) Security constraints
0
Iij, t , s
Iij,max ,
j
E _B ,
ij
t,
Line ,
k KLnew
(60)
s
k KLnew
where Vmin and Vmax are the lower and upper voltage level limits. Vj, t , s represents the nodal voltage amplitude and Iij,max denotes the current amplitude for feeder ij.
jk JK (j)
ij IJ (j)
pij, t , s
ij IJ (j)
pv + pj,wt t , s + pj, t , s
+ pjre, t , s jk JK (j)
qjk, t , s =
pjse, t , s + pechpj, t , s
ij IJ (j)
qij, t , s
ij IJ (j)
pjp, t2,gs qj,Lt , s + qjSup ,t, s
Iij2, t , s xij
pij2, t , s + qij2, t , s = Vi2, t , s Iij2, t , s , V j2, t , s
=
Vi2, t , s t,
V j2, t , s
ij
Line ,
2(pij, t , s rij + qij, t , s x ij ) +
(61)
s
t, Iij2, t , s (rij2
s
for
(62)
+
x ij2 ),
ij
Line /
Vi2, t , s
k KLnew
2(pij, t , s rij + qij, t , s x ij ) +
M (1 Lnew ,
k KLnew
t,
s
j
Lnew ,
0
lnew QijLnew , k ,max x ij, k
qij, t , s
k KLnew k KLnew
lnew PijLnew , k ,max x ij, k lnew QijLnew , k ,max x ij, k
ESS ,
(67)
s:
1,
t (68)
pjse, t , s
ess se PkESS ,max x j, k uj, t , s ,
t
ess re PkESS ,max x j, k uj, t , s ,
t
(69a)
k KES
0
x ijlnew ,k )
x ijlnew ,k
pij, t , s
ujre,0, s = 0
2(pij, t , s rij + qij, t , s x ij ) + Iij2, t , s (rij2 + x ij2)
Vi2, t , s
s:
lnew PijLnew , k ,max x ij, k
ujse, t , s + ure j, t , s
(63)
s
+ M (1 V j2, t , s
t,
(66)
In a similar way to GSU and TSU modeling, the EES can be modeled as (68)–(70).
+ qjcb, t , s + qjsvc ,t,s E _B ,
lup QijLup , k,max x ij, k
(5) Operation constraints for electrical energy storage units
pv + qj,wt t , s + qj, t , s
j
x ijlup ,k )
feeder ij is non-upgraded ( k KLup xijlup , k = 0 ), the boundary of power flow will be forced to keep the same with the original one; otherwise, it will be changed into the upgraded one. Analogously, the value of pij, t , s /qij, t , s for candidate newly-added feeder ij depends on whether feeder ij is added. If added, pij, t , s /qij, t , s is limited to the ranges of the new one, being 0 otherwise.
pj,Lt , s + pjSup ,t, s
Iij2, t , s rij
t,
k KLup
qij, t , s
The power limits of original fixed feeders ( ij Lori ) are formulated in (65) whereas that of replacing and adding ones ( ij Lup and ij Lnew , respectively) are stated in (66) and (67). It can be seen that x ijlup , k plays a critical role in restricting the boundary of pij, t , s /qij, t , s . If
(3) Power flow balance constraints
pjk, t , s =
Lnew ,
(59)
s
t,
ij
x ijlup ,k )
lup PijLup , k,max x ij, k
QijLori ,max (1 +
pij, t , s
lup QijLup , k ,max x ij, k
k KLup
Sup former station at time stage t of scenario s and P jSup ,max / Q j,max denotes the corresponding upper active/reactive power limit.
Vmax ,
k KLup
+
Sup where pjSup , t , s /qj , t , s is the active/reactive power supplied by the trans-
Vj, t , s
lup PijLup , k,max x ij, k
k KLup
s:
0
Vmin
(65)
pjre, t , s
(69b)
k KES
Iij2, t , s (rij2
+
x ij2)
,
ij SjESS ,t ,s =
) (64)
0
0+
pjse, t , s se
SjESS , t 1, s
SjESS ,t, s
t
pjse, t , s se
t
ess SkESS ,max x j, k , k
The balance of the power flow for the ADN can be ensured by (61)–(64). pv wt L q pv It should be noted that pj,Lt , s , pj,wt t , s , pj, t , s and qj, t , s , qj, t , s , j, t , s in (61) represent active and reactive power of electrical load, WTG generation, PVG generation after curtailment, respectively. And constraint (64) indicates that if k KLnew xijlnew , k = 0 , which means that candidate feeder ij will not be added into the network, there will be no voltage connection between nodes i and j .
+
pjre, t , s / re
t
pjre, t , s / re
t=1 t
t
t
ES
2
(70a) (70b)
(6) Operation constraints for capacitor banks
for 0
j
CB ,
yjcb, t , s
njcb, k , k K CB
10
s: t
(71)
Applied Energy 253 (2019) 113517
J. Wang, et al.
qjcb, t , s =
QkCB yjcb, t , s x jcb, k ,
t
NT s t=1
|yjcb, t , s
yjcb, t
(9) Operation constraints for distributed generations
(72)
k K CB
DG DG and min represent DG penetration is restricted in (80), where max the maximum and minimum permissible of DG penetration rates.
O _CB Nmax
1, s |
yjcb,0, s = 0
DG min
(73)
pjwt ,t, s =
j
SVC ,
t,
s
k K SVC
pjpv ,t ,s =
(74)
Vm2, t , s = Vi2, t , s
t,
Vmax / kijvr,min
Vm, t , s
V j, t , s
Vm, t , s
Mx ijvr
(77)
kijvr, t , s = 1 + ( kijvr % × yijvr, t , s )
(79)
Rijvr
yijvr, t , s
Rijvr , yijvr, t , s
Z
ij NT
s t=1
|yijvr, t , s
1, s |
yijvr,0, s = 0
t,
s
(83)
pv PkPV ,max nj, k
j
PVG ,
t,
s
(84)
p pv ( pjwt , t , s / j, t , s )
is written in
Apart from power curtailment of DG, load-shedding is also considered in this paper as one of the important ANM schemes. The relevant constraint is given in expressions (85) where the definition of L L Cur _L L L j,max , pj, t , s /qj, t , s and pj, t , s /qj, t , s is similar to (83) and (84). Additionally,
(80)
the active and reactive loads ( pjL, t , s and qjL, t , s ) can be obtained in a similar way to the aforementioned gas/thermal load ( f jL, t , s / pjLh , t , s ).
0
pjL, t , s =
le Le j, t , s P j
qjL, t , s
le Le j, t , s Q j
=
_L pjcur ,t,s
pj,Lt , s
=
pjL, t , s
Cur _L L j,max pj, t , s
j
L,
t,
s
_L pjcur ,t,s
qj,Lt , s = qjL, t , s × pj,Lt , s / pjL, t , s
Ui
VR :
yijvr, t
k KPVG
WTG ,
(10) Operation constraints relevant to load shedding
The voltage of the virtual node Vm, t , s can be computed by (75). The upper and lower voltages of virtual node m are limited by (76) to ensure that the voltage of node j will not exceed the limits (Vmin and Vmax ). In (77), the value of Vj, t , s is associated with the binary variable x ijvr and the value ofVm, t , s , where for x ijvr = 0 (i.e., the VR is not to be installed), the values of Vj, t , s and Vm, t , s will be imposed to be equal. And for x ijvr = 1 (i.e., the VR is to be installed), Vj, t , s will be equal to kijvr, t Vm, t , s by (77) and (78). Moreover, the operational VR ratio kijvr, t , s is reflected in (79) and (80), where an integer variable, namely yijvr, t , s denotes the tap position of VR at time stage t of scenario s. It also can be shown in (80) that the tap position variable can only vary from Rijvr to Rijvr , corresponding to the VR’s minimum and the maximum ratios. In addition, analogously to constraint (73), changing times of VR tap are restricted by (81) and meanwhile, the initial value of yijvr, t , s is imposed to be 0.
for
pv j, t , s
j
wt pv pv refers to the power factor of DG; pj,wt t , s / pj, t , s and qj, t , s /qj, t , s represent the active and reactive power of DG after curtailment. It also can be obtained that the certain triangular region in Fig. 5 is actually corresponding to the DG operation control according to (83) or (84), which indicates that the controllable range of reactive power in DG relies on the active power output after curtailment.
(76)
(78)
wt PkWT ,max nj, k
The original DG (WTG/PVG) generation
(75)
Vj, t , s = kijvr, t , s Vm, t , s
(82)
PV terms of the upper output limit for DG using type k (PkWT ,max /Pk ,max ), wt multiplied by the installation number at node j (nj, k ) and the output coefficient at time stage t of scenario s ( j, t ,wt s ). In addition, it is assumed that DGs can compensate or absorb reactive power to a certain extent, as shown in Fig. 5. Thus, the curtailment of wind and solar power is subject to constraints (83) and (84) respectively. In this set of conCur _pv _wt wt straints, jCur / pv ,max / j,max denotes the maximum curtailment rate;
s:
2(pij, t , s rij + qij, t , s xij ) + Iij2, t , s (rij2 + x ij2 )
Vmin / kijvr,max
Mx ijvr
VR ,
P jL,max
_pv Cur _pv pv pjcur ,t ,s j,max pj, t , s pv pv _pv pj, t , s = pj, t , s pjcur ,t,s pv pv qj, t , s max = pj, t , s tan(cos 1 pv ) qj,pv qj,pv qj,pv t , s max t ,s t , s max
For the purpose of improving voltage quality in the network, the investment and utilization of VR is considered in this planning as well. It is worth emphasizing that, the VR is modeled as a branch connecting buses i and j with an ideal and discretely adjustable transformer, as shown in Fig. 4. On this basis, the modeling of VR can be described by expressions (75)–(80).
m (ij), ij
pv PkPV ,max nj, k
0
(8) Operation constraints for voltage regulators
m
k KWTG
PVG k KPVG
L
_wt Cur _wt wt pjcur ,t, s j,max pj, t , s wt wt _wt pj, t , s = pj, t , s pjcur ,t,s wt 1 wt ) qj,wt t , s max = pj, t , s tan(cos qj,wt qj,wt qj,wt t , s max t ,s t , s max
where the maximum reactive power imported by SVC at node j (qjsvc , t , s ) is associated with the type of SVC and the investment decision ( x jsvc , k ).
for
wt j, t , s
j
0
Unlike CB, SVC can be continuously adjusted in a certain range, which can be stated as:
0
WTG k KWTG
j
(7) Operation constraints for static var compensations
QkSVC x jsvc ,k ,
wt PkWT ,max nj, k + j
L
DG max
where yjcb, t , s in (71) and qjcb, t , s in (72) represent the number and the reactive power output of the CB units in operation at time stage t of scenario s respectively. In order to extend life expectancy of CB switches, switching times are constrained by (73) during the planning period.
qjsvc ,t, s
P jL,min j
rij
(85)
xij
Um
1: kijvr
O _VR Nmax
Fig. 4. Voltage regulator model.
(81) 11
Uj
Applied Energy 253 (2019) 113517
J. Wang, et al.
B. Modified piecewise linearization approach
p j wt,t ( p j ,pvt )
Nonlinear polynomials aforementioned in (19), (20) and (37)–(39) can be handled by applying the method of modified piecewise linearization. _gas 2 ) . In As can be seen in (19), there exist nonlinear terms (ptsup ,s order to solve this problem, equivalent transformation is firstly made by (95) and (96).
Q
q
wt j ,t ,max
( q
pv j ,t ,max
) q
wt j ,t ,max
(q
pv j ,t ,max
Q
)
Fig. 5. Distributed generator inverter power control curve.
Sup Pgas
n + 1, t , s
In (86)–(88), auxiliary variables are introduced to replace the original square term.
t,
0
(86)
s
V j, t , s = V j2, t , s,
j
E _B ,
t,
s
(87)
Iij, t , s = Iij2, t , s ,
ij
Line ,
t,
s
(88)
j, t , s
j, t , s c 2 ij i, t , s ,
=
c 2 ij, k i, t , s
M (1
2 j,max ,
j
ij
C,
G _B ,
t,
xijc, k )
KC ,
t,
+ M (1
ij
0
2 Vmax ,
V j, t , s
2 Imax ,
Iij, t , s for:
m
ij
t,
t,
s
VR ,
Vm, t , s V j, t , s
E _B ,
Lori ,
m (ij ), ij
(Vmin/ kijvr,max ) 2 Mx ijvr
j
t,
s
Mx ijvr
s
(96)
s
(97)
Sup Pgas
max ,
t,
s
(98)
2
+
0
Sup Sup [Pgas min, Pgas max ]
,
mn, t , s , mn, t , s
n, t , s ,
(99a)
n = 1, 2, ...,Nplg
1,
1, mn, t , s = {0, 1}, n = 1, 2, ...,Nplg ,
n, t , s
t,
t,
s
(99b) (99c)
s
Nplg 0
+
(
n
n 1) n, t , s ,
t,
(Cg ( n )
Cg (
s
(100a)
n=1
Nplg n 1)) n, t , s ,
t,
n=1
s
(100b)
t,
s
(101)
sup _elec , 1) pt , s
t,
s
(102)
Eq. (103) indicates the total electricity supply power at time stage t of scenario s. On the basis of (103) and constraint (57), a preliminary _elec boundary of ptsup can be obtained as (104) ,s
C,
_elec ptsup = ,s
I j2, t , s
pjSup ,t,s , j
Sup pelec
(92)
min
t,
s
(103)
E _Sup
_elec ptsup ,s
Sup pelec
max ,
t,
s
(104)
The approach of modified piecewise linearization is also applied in (105)–(106).
(93)
Ce (p) = 2 2 p2 +
s
(Vmax /kijvr,min ) 2
Vm, t , s
t,
_elec _elec E (ptsup ) = (2 2 ptsup + ,s ,s
Analogously, the nonlinear and can be eliminated by substituting V j, t , s and Iij, t , s into (61)–(64) and (75). Expressions (59), (60) and (76) can be reformulated by (92), (93) and (94), respectively. 2 Vmin
fp,
_Elec _elec CsSup = E (ptsup ) t, ,t ,s
(91)
s
terms V j2, t , s
x ijc, k ),
t,
Similarly, the reformulation on (20) is presented in (101) and (102) as well.
(90) c 2 ij, k i, t , s
f jSup ,t,s
_gas G (ptsup ) = Cg ( 0 ) + ,s
(89)
s
s
j, t , s
k
t,
1
_gas ptsup = ,s
Consequently, substituting j, t , s into expressions, (35), (43) and (45) can be rewritten as linear constraints (89)–(91), respectively. In a similar way, i2, t , s can be removed from (37)–(39) as well. 2 j,min
+
(95)
G _Sup
_gas ptsup ,s
min
Cg ( ) =
A. Auxiliary variable substitution
G _B ,
=(
sup _gas 1 pt , s
s
By applying the modified piecewise linearization, the boundary Sup Sup [Pgas min, Pgas max ] is equally divided into Npl g intervals. In order to ensure the equivalent transformation, two auxiliary variables are introduced: (i) a binary variable mn, t , s , which implies that whether the nth interval is selected, (ii) a continuous variable n, t , s which is the optimal incremental value in each interval. Their relationships are character_gas ized by (99), and the gas supply power ptsup as well as the gas pur,s _gas ) ) can be chase cost at each time stage of the scenario (i.e., G (ptsup ,s represented by the mentioned auxiliary variables respectively, as shown in (100a) and (100b). Based on this, the nonlinear polynomial function (96) can be rewritten as (100). Thus, Eq. (19) is finally replaced by (95) and (99)–(100).
In what follows, we will derive equivalent linear expressions for each of them.
j
sup _gas , 0 ) pt , s
j
a. Quadratic terms such as i2, t , s , Iij2, t , s , V j2, t , s b. The nonlinear expressions like Eq. (19), Eq. (20) and the constraints (37)–(39) c. The nonlinear power flow constraint (62) hs sh hs rh d. Bilinear terms consisting of x jgs, k ujsg, t , s , x jgs, k urg j, t , s , x j, k u j, t , s , x j, k u j, t , s , ess se ess re vr cb cb x j, k uj, t , s , x j, k uj, t , s , yj, t , s x j, k and kij, t , s Vm, t , s e. Items with absolute value, including the switching times limits for CBs and VRs and objective expression (9).
=
_gas G (ptsup ) ,s
_gas ptsup = ,s
The programming model established in this paper is a complicated mixed integer nonlinear programming (MINLP) with the following nonlinear terms:
j, t , s
t,
Given the flow-power coefficient fp , the total gas supply power at each time stage of the scenario can be computed by (97) and in (98), the _gas boundary of ptsup can be obtained which is based on (34) and (97). ,s
5. Equivalent reformulation—linearization
2 j, t , s ,
_Gas _gas CsSup = G (ptsup ) t, ,t ,s
n + 1, t , s
(94)
0 12
n, t , s
1 p,
mn, t , s , mn, t , s
p
Sup [pelec n,
min ,
Sup pelec
n = 1, 2, ...,Nple
1, n = 1, 2, ...,Nple ,
t,
s
(105a)
max ]
1,
t,
s
(105b) (105c)
Applied Energy 253 (2019) 113517
J. Wang, et al.
C. Second order cone relaxation
Nple
_elec ptsup = p0 + ,s
(pn
pn 1 )
n, t , s ,
t,
s
(106a)
n=1
Although Iij2, t , s and V j2, t , s in constraint (62) have been substituted by Iij, t , s and V j, t , s in part A, there remain some nonlinear items including pij2, t , s , qij2, t , s and Vi, t , s Iij, t , s . In order to cope with them, SOCR is applied to transform non-convex constraint (62) into (113), as done in [14].
Npl _elec E (ptsup ) = Ce (p0 ) + ,s
(Ce (pn )
Ce (pn 1 ))
n, t , s,
t,
s
(106b)
n=1
Consequently, Eq. (20) can be rewritten as (101) and (105)–(106). Even though constraints (37)–(39) have been addressed in the part A, there remain some issues to be solved due to the complicated nonlinear term sgn( i, t , s, j, t , s ) fij2, t , s . In order to solve this problem, equivalent transformation is firstly made in (107) where sgn( i, t , s, j, t , s ) fij2, t , s can be seen as fij, t , s |fij, t , s | and the sign of fij, t , s is up to the comparison of the pressure at both ends for the facts that nodal pressure is nonnegative and gas always flow from higher pressure node to lower pressure node.
sgn(
2 j, t , s ) fij, t , s
i, t , s
= fij, t , s |fij, t , s |
2pij, t , s 2qij, t , s Iij, t , s
Pipe Pipe [ Fmax , Fmax ]
Pipe Fmax = max{FijPipe ,max }
, ij
t,
s
0
(108)
mn, t , s , mn, t , s
n + 1, t , s
0
n, t , s
n, t , s ,
0
(109a)
Pipe Pipe [ Fmax , Fmax ]
n = 1, 2, ...,Npl
1, n = 1, 2, ...,Npl,
t,
1,
t,
(109c)
s
0
Based on above, the constraints (37)–(39) are rewritten as (110)–(112) in combination with M (a very large number).
fij, t , s = f0 + ij ( i, t , s Pipe ,
j, t , s )
t,
n= 1
= Y (f0 ) +
(fn Npl n= 1
j, t , s )
Npl n=1
Y (f0 ) +
(fn Npl n=1
j, t , s )
Y (f0 ) +
Npl n=1
k KPipe
Pipe ,
fij, t , s fij, t , s
f0 + f0 +
n=1 Npl n=1
i, t
Y (f0 ) +
ij ( i, t
j, t )
Y (f0 ) +
t,
s
t,
x ijpipe ,k
k KHS
j
s
GS ,
t,
s (114b)
hs PkHS ,max x j, k
j
Mujsh, t , s
k KHS
t,
(114a)
gs FkGS ,max x j, k
Murg j, t , s
pjrh, t , s
HS ,
t,
s (115a)
hs FkHS ,max x j , k
j
Murh j, t , s
,
k K GS
pjse, t , s
pjre, t , s 0
n, t , s
)
pjse, t , s 0
0
Y (fn 1 ))
k KGS
GS ,
HS ,
t,
s (115b)
n, t , s
n, t , s
Y (fn 1 ))
j
Mujsg, t , s
pjsh , t,s
pjrh, t , s
gs FkGS ,max x j, k
ess PkESS ,max x j , k
j
Mujse, t , s
k KES
pjre, t , s
ESS ,
t,
s (116a)
ess FkESS ,max x j, k
j
Mure j, t , s
ESS ,
t,
s (116b)
In constraint (117), M is used again to eliminate the nonlinear item yjcb, t , s x jcb, k in (72).
ij
for
n, t , s
j
QkCB yjcb, t , s
x ijpipe ,k )
CB ,
k
K CB ,
x jcb, k )
M (1 M
(111)
s
k KCB
x jcb, k
t,
s:
qjcb, t , s
QkCB yjcb, t , s + M (1
qjcb, t , s
M k K CB
x jcb, k )
x jcb, k
(117)
In addition, there still exists a non-convex item related to in VR operational constraint, which is a crucial problem to be solved. In order to solve it, a set of binary variables aijvr, n, t , s is introduced in (118) to formulate the VR tap position firstly.
kijvr, t , s2 Vm, t , s
Npl
j, t )
C,
Y (fn 1 ))
(Y (fn )
(Y (fn )
pjsh , t,s
0
ij
0
k KPipe
M (1
ij (
0
(Y (fn )
s (113)
k K GS
f jrg, t , s
n, t , s
,
fn 1 )
f jrg, t , s
0
(110)
+ M (1 ij ( i, t , s
fn 1 )
s
fij, t , s = f0 + ij ( i, t , s
Npl
t,
2
f jsg, t , s
0
(109b)
s
f jsg, t , s 0
Therefore, the method of modified piecewise linearization is applied as follows:
Y (f ) = f |f |, f
Line ,
x jgs, k ujsg, t , s and x jgs, k urg j, t , s in (48) are further resolved in (114). It still can be seen that the value of gas flow into/out of the gas storage unit ( f jsg, t , s / f jrg, t , s ) is limited by the boundary [0, FkGS ,max ] only when the k -type gas storage unit is installed at node j (x jgs, k = 1) and in the operating state of gas storage/release (/urg j, t , s = 1); otherwise, it is imposed to be 0 (only if x jgs, k = 0 or ujsg, t , s = 0 /urg j, t , s = 0 ). By the same way, thermal storage constraint (55) and ESS constraint (69) can be reformulated respectively as (115)–(116).
(107)
Pipe ,
ij
D. Bilinear term relaxation
Then, the range of fij, t , s is determined by (40), and hence gas flow through each pipeline is also subjected to the boundary as (108). It Pipe depends on the largest flow of all pipelines (Fmax ).
fij, t , s
Iij, t , s + Vi, t , s,
Vi, t , s
(fn
fn 1 )
n, t , s
(fn
fn 1 )
n, t , s
Npl n=1 Npl n=1
+M k KC
M k KC
(Y (fn )
Y (fn 1 ))
n, t , s
(Y (fn )
Y (fn 1 ))
n, t , s
x ijc, k x ijc, k , +M k KC
M k KC
ij
yijvr, t , s =
x ijc, k
Rijvr +
2Rijvr
aijvr, n, t , s n,
ij
VR ,
t,
s
n=0
2Rijvr
x ijc, k
aijvr, n, t , s = 1,
ij
VR ,
t,
s
n=0
(112)
Then, the square of the VR ratio can be replaced by (119). 13
(118a)
(118b)
Applied Energy 253 (2019) 113517
J. Wang, et al. 2Rijvr
kijvr, t , s 2 =
1+
Rijvr
n
2Rijvr
n=0
Xijvr = kijvr % × 2Rijvr ,
2
Xijvr
ij
aijvr, n, t , s , ij
VR ,
t,
s
V j, t , s =
(119a)
M (1
V j, t , s =
1+ n= 0
j
VR ,
Rijvr
n
m
2Rijvr t,
s
aijvr, n, t , s Vm, t , s,
m
Rijvr
Xijvr ) 2uijvr, n, t , s
2Rijvr
Vm, t , s
uijvr, n, t , s
m (ij), ij
M (1
aijvr, n, t , s )
Maijvr, n, t , s VR ,
t,
s
(121)
E. Absolute value processing
2
Xijvr
n
uijvr, n, t , s
Maijvr, n, t , s
After that, the constraint (78) can be rewritten as (120). 2Rijvr
(1 +
n=0
aijvr, n, t , s )
(119b)
VR
2Rijvr
m (ij ), i
In (122), integer variables n+j, tcb, s and nj, tcb, s are introduced for linearizing the absolute value |yjcb, t , s yjcb, t 1, s | in (73).
(120)
for
Furthermore, the bilinear term aijvr, n, t , s Vm, t , s in (120) is also an important problem, which can be solved by introducing uijvr, n, t , s and actually the values of uijvr, n, t , s and aijvr, n, t Vm, t , s are equal. In order to ensure that, the equivalent transformation can be finally shown as:
j
CB ,
n+j, tcb, s
0, nj, tcb, s
n+j, tcb, s
nj, tcb, s =
s:
0, yjcb, t , s yjcb, t , s
(122a)
t 0 yjcb, t 1, s
Fig. 6. Flowchart of the proposed modeling methodology. 14
t=1 t
2
(122b)
Applied Energy 253 (2019) 113517
J. Wang, et al. NT s
n+j, tcb, s + nj, tcb, s
min
O _CB Nmax
pipe gs gs ess lup lnew cb svc vr wt {xij, k , xijc, k , x j, k , x j, k , x hs j, k , x j, k , xij, k , xij, k , x j, k , x j, k , xij , nj, k,
(122c)
t=1
pv +vr cb vr vr nj, k , ncb j, k, y j, t , s , yij, t , s , rij, t , s, nij, t , s, j, t , s, V j, t , s,Iij, t , s, Vm, t , s, sg rg + cb cb se re vr rh nij, t , s, nj, t , s, n j, t , s, u j, t , s, uj, t , s, ush j, t , s, u j, t , s, uj, t , s, u j, t , s }
Analogously, absolute value in (81) can be substituted by (123).
for
ij
nij+, tvr, s
0,
nij+, tvr, s
subject to : 1.linear constraints: {(1)
s:
nij, tvr, s
0,
yijvr, t , s
nij+, tvr, s + nij, tvr, s
+
VR
yijvr, t
t
1, s
2
(123b)
O _VR Nmax
+ j
SVC
c O _svcqjsvc , t, s +
Line
j
PVG
j
L
+
(75), (82) (95), (97), (117), (121) (124)}
In order to give a clear description of the proposed modeling methodology, a flowchart is presented in Fig. 6. The modeling procedure listed in block diagram can be summarized as follows: Step 1: Uncertainty modeling. In order to describe the operating conditions, a probabilistic scenario generation method is developed with the consideration of extreme and normal situations. The details which have been illustrated in Section 3 are briefly presented in the right side of Fig. 6 as well. Step 2: Multi energy system modeling. To establish the comprehensive planning model, the NGN, ADS models are investigated separately and then EH is modeled for coupling these two networks. By incorporating NGN, ADN and EH, an expansion planning model for the MES can be formulated as an MINLP problem. The specific modeling has been clearly given in Section 4. Step 3: Equivalent reformulation. Since the proposed MINLP model is difficult to be solved, the mathematical techniques are developed to transform the original version into a MISOCP. The relevant processing methods and procedures have been illustrated in Section 5. Step 4: Solve the model by calling CPLEX.
j
CB
cb cb cO _cb (n+ j, t + n j, t )
cO _lossIij2, t , s rij
j
ESS _pv c Cur _pvpjcur ,t ,s
+
(42), (44),
6. Framework
(123c)
vr c O _vr (nij+, vr t + nij, t ) +
t ij
(34), (36), (40)
t=1
According to the analysis (122) and (123), Eq. (9) can be rewritten as (124). ij
(14), (21)
(125)
0
t=1
_Elec CsOpe = ,t
(8), (10)
(46) (47), (49) (54), (56) (58), (61), (63) (68), (70) (71), (74) (99) (101), (103), (105) (106), (109) (112), (114) 2.conic constraints: \{ (113)\}
(123a)
t
yijvr, t , s
nij, tvr, s =
NT s
VR ,
CTotal
_L c Cur _Lpjcur , t, s +
j
cO _ess (pjse, t , s + pjre, t , s ) +
cO _wtpjwt , t, s + WTG
j
WTG
_wt c Cur _wtpjcur ,t ,s
cO _pvpjpv ,t ,s
j
(124)
PVG
Thus, the complex nonconvex MINLP has been transformed as a MISOCP model, which is formulated as:
19 1
18 17
16
20 2
37 36
3
33
23
32
31
15
Feeders to be upgraded Candidate WTG
Candidate PVG
22
4
5
24 30
13
14
Substation
21
35
7
6
8
26
25 29
34
28
12
9
27
11
10
Slack bus
Load bus
Newly added bus
Existing feeders
Hub Unit (No.1)
Candidate VR
Candidate ESS
Candidate CB
Candidate SVC
Fig. 7. Schematic of the modified IEEE 33-node system (adding 4 new load bus, 8 candidate feeders and other potential candidate units). 15
Applied Energy 253 (2019) 113517
J. Wang, et al.
Fig. 8. Schematic of the 23-node gas network system (with 3 new load bus, 7 candidate pipelines, 2 potential candidate compressors and 2 candidate gas storage units). Table 1 Parameters for the candidate distribution feeders.
Table 4 Parameters for available alternatives of CB and SVC units.
No.
From bus
To bus
Length (km)
Resistance (Ω)
Reactance (Ω)
CB
33 34 35 36 37 38 39 40
7 8 21 22 31 32 18 33
34 34 34 35 35 35 36 37
1.00 0.75 0.63 1.13 0.33 0.75 0.50 0.88
0.40 0.30 0.25 0.45 0.15 0.30 0.20 0.35
0.19 0.23 0.19 0.38 0.11 0.23 0.15 0.26
Capacity (MVar)
Cost (103 $/unit)
Capacity (MVar)
Cost (103 $/unit)
0.05 0.12
200 500
0.40 0.80
8000 14,000
I II
Table 5 Parameters for available alternatives of ESS units.
Table 2 Parameters for available alternatives of feeders. Feeders to be upgraded
Ⅰ II
I II
Feeders to be added
Capacity (MW)
Cost (103 $/km)
Capacity (MW)
Cost (103 $/km)
1.4 3.0
12,000 14,000
1.4 4.0
15,000 22,000
SVC
Capacity (MW)
Power limit (MW)
Efficiency
Cost (103 $)
10 30
3 9
98% 98%
6000 10,000
system is employed, as shown in Fig. 7 and Fig. 8 respectively. Simulations have been implemented by using solver CPLEX 12.6.0. 7.1. Overview of the test system
Table 3 Parameters for available alternatives of DG units. WTG
I II
In the original system, there are 33 buses and 32 feeders totally, and corresponding system parameters can be found in [33]. As can be seen from Fig. 7, compared with the IEEE 33-node network, four load buses (with the same power 220 kW and 110 kVar) and 8 candidate feeders are added to the modified distribution system. Meanwhile, several potential locations (i.e., buses or feeders) for the investment of candidate active equipment are marked in the system. The related parameters for all alternative components are listed in Tables 1–5 respectively. According to Fig. 8, the presented gas system is composed of 2 gas source nodes, 2 gas compressors, 26 gas pipelines, 20 existing gas load nodes and 3 newly added nodes. However, in order to meet the growing demand for gas, expansion investment is still needed. Hence, 7 candidate gas pipelines (marked with black dashed line), 2 candidate gas compressors and 2 candidate gas storage units are selected to be
PVG
Rated power (MVA)
Cost (103 $)
Rated power (MVA)
Cost (103 $)
2 4
22,000 31,000
3 5
24,000 29,000
7. Case study To demonstrate the effectiveness of the proposed planning model, a modified IEEE 33-node distribution system coupled with a 23-node gas 16
Applied Energy 253 (2019) 113517
J. Wang, et al.
Table 6 Data of all gas pipelines. Existing gas pipelines
Candidate gas pipelined
No.
Pipeline (i–j)
Length (km)
ϕij (1012)
No.
Pipeline (i–j)
Length (km)
ϕij (1012)
No.
Pipeline (i–j)
Length (km)
ϕij (1012)
1 2 3 4 5 6 7 8 9 10 11 12 13
1–2 1–2 2–3 2–3 3–4 5–6 5–8 6–7 4–7 4–14 8–9 8–9 9–10
0.4 0.4 0.6 0.6 2.6 4.3 2.6 2.9 1.9 5.5 0.5 0.5 0.2
0.0052 0.0052 0.0026 0.0052 0.0035 0.0017 0.0035 0.0024 0.0035 0.0017 0.0035 0.0017 0.0035
14 15 16 17 18 19 20 21 22 23 24 25 26
9–10 10–11 10–11 11–12 12–13 13–14 14–15 15–16 11–17 11–17 17–18 18–19 19–20
2 2.5 2.5 4.2 4 0.5 1.0 2.5 1.05 1.05 2.6 9.8 0.6
0.0069 0.0069 0.0069 0.0035 0.0017 0.0122 0.0035 0.0017 0.0026 0.0026 0.0026 0.0035 0.0035
27 28 29 30 31 32 33
18–21 19–21 20–21 12–22 13–22 15–23 16–23
0.6 1.05 1.05 2.6 2.5 0.5 0.5
0.0035 0.0026 0.0026 0.0026 0.0069 0.0122 0.0035
Table 7 Parameters for available alternatives of gas pipelines and gas compressors. Candidate gas compressor
I II
Table 10 Parameters for available alternatives of thermal storage units.
Candidate gas pipeline
Flow limit (106 m3/h)
Ratio
Cost (103 $)
Flow limit (106 m3/h)
Cost (103 $)
10 10
1.03 1.08
10,000 15,000
1.2 6
10,000 15,000
I II
Thermal power limit (MW)
Efficiency
Cost (103 $)
10 30
3 9
95% 95%
4200 8000
7.2. Scenario set Based on the generated scenarios proposed in Section 3, the scenario set can be constituted as follows:
Table 8 Parameters for available alternatives of gas storage units.
I II
Capacity (MWh)
Capacity (m6)
Flow limit (106 m3/h)
Efficiency
Cost (103 $)
0.06 0.12
0.005 0.015
95% 95%
10,000 20,000
(1) Scenario ΦA: The combination of the maximum value curves in Fig. 2 and the minimum value curves in Fig. 3. (2) Scenario ΦB: The combination of the minimum value curves in Fig. 2 and the maximum value curves in Fig. 3. (3) Scenario ΦC: The combination of the mean value curves in all graph.
available investment options, as shown in Tables 6–8. In addition, ADN and GN systems are coupled by totally 8 energy hub units. The related data of hub unit is given in Table 9. It is worth mentioning that thermal storage units are only optional investments embraced in hub units, whose parameters are shown in Table 10. Table 11 depicts other supplement parameters in the original system. Meanwhile, there are some reasonable assumptions as follows: (1) Annual maintenance costs for all units are set to be to 5% of investment cost. (2) Unit production cost for WTGs and PVGs are assumed to be 0. (3) Unit costs of energy losses of TSs and feeders are all 200 $/MWh. (4) Unit active curtailment costs of wind, photovoltaic powers are both 3000 $/MWh and unit load-shedding cost is 5000 $/MWh. (5) All existing feeders are considered to be upgraded, and all load buses are schedulable except the slack bus. (6) The planning scheme is implemented in 10 years and the interest rate (d) is 10%. (7) Electricity, gas and thermal demands are doubled in next 10 years. (8) Big constant M is 10000000.
By applying the Step 4 in Section 3, the corresponding probabilities of the 3 scenarios are 0.0465, 0.1079 and 0.8457, respectively. 7.3. Results of expansion planning model In this section, the aforementioned probabilistic scenario set is employed to test the model. In order to analyze the performance of the proposed model, four cases are compared and presented as follows: Case 1: model (125) with probabilistic scenarios: ΦA, ΦB and ΦC. Case 2: model (125) with mean value scenario ΦC only. Case 3: model (125) with three seasonal average scenarios: daily average scenarios for summer, intermediate (spring and autumn) and winter. The corresponding probabilities are 0.25, 0.5 and 0.25, respectively. Case 4: model (125) with five scenarios: two extreme scenarios (ΦA and ΦB) and three average seasonal scenarios with probabilities 0.0465, 0.1079, 0.2281, 0.3926 and 0.2250, respectively.
Table 9 Data of energy hub units. No.
1
2
3
4
5
6
7
8
Bus number in ADN Node number in GN Heat demand (MW)
10 14 1.07
14 10 2.15
2 12 0.54
28 16 2.68
32 6 1.61
4 15 0.54
21 2 1.07
23 11 1.13
Conversion efficiency
P2G: 65%, GF: 92%, CHP: 1) From gas to heat 36%; 2) From gas to electricity 18%
17
Applied Energy 253 (2019) 113517
J. Wang, et al.
Table 11 Other parameter setting in the original system. Gas system
Active distribution system
Gas nodal pressure limits Gas flow limits (in existing pipelines) Gas source flow limits GHV (ζfp) Ratio of existing compressor ( ijc )
50–70 bar 0–5 m6/h 0–20 m6/h 0.0096 MW/m3 1.1
Voltage magnitude range Maximum current magnitude Maximum power in existing feeders Maximum switch times of CB and VR Maximum investment numbers of SVCs, CBs, WTs and PVs
Gas price coefficients
γ1 = 0.5 γ0 = 2
DG power factor (
0.94–1.06 pu 3 pu 1.6 MW/1.2Mvar 24
SVC CB Nmax = 1, Nmax = 16 ,
wt pv Nmax = 5, Nmax =5 0.8
wt / pv )
Maximum curtailment rates of WTs, PVs and demand Electricity price coefficients
20% 2 = 0.3 1 = 1.2 0 = 0.1
Table 12 Investment decision of all components (including type, location and installation number).
Active distribution network
Infrastructures
Type
Case 1
Upgraded feeder
I II I II I II I II I II I II I I II
23(1), 6(1) 34(1), 37(1), – 13(1), – 6(1) – – 33(5) 3(1) – 1(1) 7(1) –
Added feeder WTG PVG CB SVC VR ESS
24(1) 35(1), 40(1) 36(1)
Case 2
Case 3
Case 4
22(1), 23(1) – 34(1), 35(1), 37(1), 40(1) – 13(1) – – 8(1) – 33(8) 3(1) – 1(1) 7(1) –
14(1), 23(1) 6(1) 34(1), 35(1), 37(1) 39(1) 13(1) – 6(1) – – 33(6) 3(1) – 1(1) 7(1) –
23(1) 8(1), 9(1) 34(1), 35(1), 38(1), 40(1) – 13(1), 36(1) – 8(1) – – 33(6) 3(1) – 1(1) 7(1) –
Hub unit
Thermal storage
I II
3(1), 8(1) –
8(1) –
8(1) –
7(1), 8(1)
Natural gas network
Added gas pipeline
I
1(1), 5(1), 7(1) – – 8(1) 9(1) –
2(1), 4(1), 6(1) – 8(1) – 9(1) –
1(1), 5(1), 6(1) – 8(1) –
1(1), 5(1), 7(1) – 8(1) – – –
Gas storage unit Compressor
II I II I II
Note: number within () is the amount of investment; number outside () is the location.
thermal storage and gas storage units are installed to satisfy the operation constraints. Compared with Case 1, two feeders are required to be upgraded and there is less development in renewable DGs. Since investment cost (see Table 13) in active management equipment is far less than Case 1, Case 2 has a lower total planning cost; however, its annual energy purchase cost has a 13.1% increase compared with Case 1, which means that Case 2, with less investment decisions, requires more energy consumption than Case 1 during the planning horizon. In this respect, the probability model with extreme scenario set has a positive impact on reducing the energy purchase cost, but it increases the cost of investment and makes the 10-year planning cost slightly increase by $ 5.61 × 106. It should be noted that there is usually a trade-off between the investment costs and the operation performance (including energy consumption during planning period). Although less investment can be obtained under the average scenario ΦC (i.e., Case 2), such an investment scheme may lead to a significant increase in comprehensive operation cost in some extreme cases, and even can not meet the requirements of system operation under some circumstances. On the contrary, the probabilistic model in Case 1 synthesizes the extreme scenarios and considers the robustness of the investment decision, which makes the investment scheme competent for the occurrence of
Table 13 Planning cost component comparisons with the consideration of different scenario sets ($ 108).
Investment cost Annual comprehensive operation cost Annual maintenance cost Annual energy purchase cost Annual DG curtailment cost Annual demand curtailment cost Total planning cost
Case 1
Case 2
Case 3
Case 4
1.5059 0.1624 0.0753 0.0496 0.0004 0.0002 2.6035
1.3668 0.1747 0.0683 0.0561 0.0005 0.0002 2.5474
1.3731 0.1650 0.0686 0.0553 0.0002 0.0002 2.5220
1.5209 0.1584 0.0765 0.0483 0.0001 0.0001 2.6018
Note: the annual comprehensive operation cost represents the sum of operation (including curtailment), maintenance and energy purchase costs.
Tables 12 and 13 present the simulation results including all the investment decisions and planning costs under Cases 1–4. In Table 12, the investment results of Case 1 indicate that during the planning horizon, the system has been developed to response the demand growth by upgrading three feeders, adding four new feeders, three gas pipelines and three DGs. In view of the objective that is to minimize the total costs, CB, SVC, VR, compressor and energy storage units including ESS, 18
Applied Energy 253 (2019) 113517
Input power (MW)
Input gas (103m3/h)
Input gas (10 3m3/h)
J. Wang, et al.
Scenario ΦA
0.2 0.1
0
Scenario
0.2
ΦB
0.1
0
12
24
0.4
0
0.1
0
0.5
12
(a) Input gas to GF
24
0
12
24
0
12
24
0
12
24
0.2
0
12
24
0
0.4
0.4
0.2
0.2
0
0 0.4
0.2 0
Scenario ΦC
0.2
0
12
24
0
0
12
(b) Input gas to CHP
24
0 0.4 0.2
0
12
24
(c) Input electricity to P2G
0
Thermal power Thermal residual capacity
0.05
0
0
5
10
15
20
Thermal power Thermal residual capacity
0.2
0.1
0
-0.1
0
5
10
15
20
Thermal power/capacity (pu)
Scenario ΦB
Scenario ΦA
0.1
Thermal power/capacity (pu)
Thermal power/capacity (pu)
Fig. 9. Energy input to the EH (No.1) in each time stage of the 3 different scenarios.
Scenario ΦC
0.1
Thermal power Thermal residual capacity
0.05
0
0
5
10
15
20
Fig. 10. Storing/Releasing thermal power status and capacity for thermal storage unit in the EH (No. 8) at each time stages in 3 different scenarios.
any uncertain scenarios. For further investigating the impact of seasonal scenarios on the planning results, three scenarios (summer, intermediate and winter) are calculated and constituted by averaging the hourly power values throughout all days in each season. Note that spring and autumn are together considered as the intermediate scenario. In Case 4, the probabilities of summer, intermediate and winter scenarios are obtained by counting their frequency of occurrences except for the extreme hours in ΦA and ΦB. Without the consideration of extreme scenarios, the resulting schemes in Case 3 has one less DG investment than that in both Case 1 and Case 4, thereby allowing a fewer investment cost. Moreover, it makes no great difference to each resulting cost between Case 2 and Case 3, even if Case 3 uses a combination of seasonal scenarios. It can also be observed that the results in Case 4 are relatively similar to that in Case 1 (the variations in investment, comprehensive operation, energy purchase and total planning costs are no more than 1%, 2.44%, 2.63% and 0.066%, respectively), indicating that the mean value scenario ΦC (total 24 h/time stages) has similar effect to the three seasonal average scenarios (total 3 × 24 h/time stages). Since the reduction of scenario number is significant for a complicated long-term planning issue, using the combination of ΦA, ΦB and ΦC is better than the
combination of ΦA, ΦB and three seasonal scenarios via fewer operating conditions. According to the above discussion, the different results on Cases 1–4 show the dominant status of extreme scenarios ΦA and ΦB over other scenarios. To make further analysis of the energy conversion in EH, the simulations of Case 1 are presented in Figs. 9, 10 and 11. Among them, Energy input changes to EH in each scenario are shown in Fig. 9; the operation of thermal storage unit is shown in Fig. 10; and the variation curves of electricity and gas price are shown in Fig. 11. It can be observed that: (1) In the case of scenario ΦA, due to the abundant wind and solar resources and the low electrical demand (especially from 12 h to 20 h), as well as the low electricity price level (see Fig. 11), the additional power generation is converted into natural gas through P2G, which together supplies natural gas demand; in the mean time, thermal demand is mainly provided by GF. During the rest of the period, thermal demand is together provided by natural gas through GF and CHP (see the left side of Fig. 9). In addition, the thermal storage unit (see Fig. 10) stores the excess heat energy from 1 h to 5 h, and then releases them from 18 h to 24 h to respond the heat demand peak at night. 19
Applied Energy 253 (2019) 113517
Electrical price $/(103MWh)
Gas price $/(103m3)
J. Wang, et al.
1000
Scenario ΦC
Scenario ΦB
Scenario ΦA
500 0
12
120
24
12
time (h) (a) Gas price Scenario ΦB
Scenario ΦA
80
24
12
24
Scenario ΦC
40 0
12
24
12
24
time (h) (b) Electricity price
12
24
Fig. 11. Energy price variation curves with time stages under three scenarios.
meet the higher demand: storing energy from 1 h to 9 h, and releasing them from 10 h to 24 h. (3) In the normal Scenario ΦC, although it has a lower requirement of energy supplies than Scenario ΦB, because of the low price of natural gas and the the overall high level electricity price in this scenario (see Fig. 11), natrual gas plays a leading role in the energy conversion and supply. Most of the gas injected into CHP to supply electrical and heat demand, while relatively little gas inputs GF to secondarily supply heat demand from moring to night (see the right side of Fig. 9). Simarly to other scenarios, heat release occurs mainly at night (see Fig. 10).
Table 14 Planning cost components for the separate models (ADN and NGN) ($ 108).
Investment cost Annual maintenance cost Annual energy purchase cost Annual DG curtailment cost Annual demand curtailment cost Total planning cost
The single ADN model
The single NGN model
Sum
1.5213 0.0761 0.0496 0.0126 0.0022
0.5440 0.0272 0.0430 – –
2.0653 0.1033 0.0926 0.0126 0.0022
3.2394
1.0186
4.2580
7.4. Comparison of the comprehensive model and separate models
Table 15 Optimal planning costs of different cases with full/partial consideration of ANMs (unit $ 106).
VR adjustment ESS scheduling CB and SVC control DG control Investment cost Annual maintenance cost Annual comprehensive operation cost Annual energy purchase cost Annual active curtailment cost Total planning cost
Case 1
Case 5
Case 6
✓ ✓ ✓ ✓ 150.5914 7.5296 16.2391
✓ ✓
✓ ✓ ✓
4.9636 0.062 260.3517
6.3366 0.0276 292.2024
✓ 148.9606 7.4480 21.3329
To verify the superiority of the comprehensive planning model where ADN and NGN are optimized coordinately (i.e., Case 1), the separate planning models for single ADN and NGN are conducted without the EH energy conversion, respectively. The cooresponding results on ADN and NGN are given in Table 14. It is obvious that the separate models present less economic than the joint expansiong planning model in Case 1. Specifically, compared with Case 1, there are a 37.15% increase of investment cost and a 86.69% increase of energy purchase cost; thus, the sum of planning cost is increased by 63.55%. Apparently this results show the importance and necessity of considering the interdependency of ADN and NGN in the MES planing problem.
Case 7
138.5041 6.9252 20.4271
✓ ✓ 148.0040 7.4002 19.2761
6.1791 0.0338 273.4036
6.2141 0.0163 276.0944
7.5. Analysis of the active network management
Note: the annual comprehensive operation cost represents the sum of operation (including curtailment), maintenance and energy purchase costs.
In order to verify the advantages of ANMs, the proposed model with full consideration of ANMs is compared with the model with partial consideration of ANMs. The different cases and results conducted based on the probabilistic scenario set are listed in Table 15. Since the full ANM schemes are considered in the proposed model, the investment cost of Case 1 is relatively higher than other cases in Table 15, however the comprehensive cost of annual operation cost as well as the total planning cost is lower. It is obvious that the model with full consideration of ANMs is more economical than that with partial consideration of ANMs, since ANMs make the multi-energy system more controllable and adjustable in controlling active and reactive power. By comparing Case 1 with Case 5, we can observe that the
(2) In the extreme Scenario ΦB, since the electrical demand outstrips the renewable energy supply, ADN needs more natural gas to inject into CHP (see the second graph of (b) in Fig. 9) for converting the gas into electricity. Due to thermal power along with CHP operating, the heat requirements are supplied by natrual gas primarily through CHP and secondarily through GF. It can also be seen that P2G has no power input in the meantime. Also, energy prices are higher than other scenarios in most of the time (espeacially 0–10 h and 16–24 h) according to Fig. 11, indicating a shortage of energy supplies due to the high-level load demand. Moreover, the thermal storage uint (see Fig. 10) works more frequently through the day to 20
Applied Energy 253 (2019) 113517
ΦC
Electric Demand
0
12
24
(a)
24
WTG
PVG
12
24 WTG(MW)
1.1
12
1 0
12
24
1
12
(c)
24
12
0
12
24
12
(e)
24
12
5
0
12
24
12
(g)
24
12
24
0
Scenario
Curtailment
12
ΦB
Scenario ΦC
Without curtailment
24
2
12
(b)
With curtailment
24
12
24
1 Curtailment
12
Without curtailment
24
12
(d)
With curtailment
24
12
24
2 1 0
24
Scenario ΦA
0.2
0
24
0.5
0.4
Curtailment
12
Without curtailment
24
0.2
12
(f )
With curtailment
24
12
24 1
0
-0.2
0 ESS power
12
time(h)
ESS capacity
24
12
(h)
24
12
24
-1
ESS capacity(pu)
Scenario
Demand(MW)
ΦB
PVG(MW)
SVC(MVar)
Scenario
0.5
0.9
Number of CBs
Scenario ΦA
1
ESS power (pu)
VR ratio
Active power(pu)
J. Wang, et al.
time(h)
Fig. 12. The 24-h operation schemes for each of the 3 scenarios considering active distribution management. (a) electric demand, wind and photovoltaic power generation curves; (b) load-shedding for bus 34; (c) ratio changes of VR between buses 1 and 2; (d) active power curtailment for WTG at bus 13; (e) reactive power output of SVC at bus 3; (f) active power curtailment for PVG at bus 6; (g) reactive power output of CBs at bus 33; (h) charging/discharging power and capacity of ESS at bus 7.
model without the reactive compensation of CBs and SVCs increases the comprehensive operation cost by 31.37%, and total planning cost by 12.09%. Furthermore, 27.66% more energy consumptions are required in Case 5 during the planning horizon. In Case 6, without considering the control of DG, the results also show an increase in the operation, energy purchase, and total planning costs by 25.79%, 24.45%, 5.01%, respectively. It indicates that the flexible control of DG has a positive impact on reducing the operation and planning cost for MES. Since the investment costs for DG are higher than that of other infrastructures, the total investment cost of Case 6 without DG is lower than that of other cases. Similarly, compared with the full consideration of ESS and VR, the more expensive cost results in Case 7 show that the model with the consideration of ESS and VR can also enjoy more economic value. In order to more intuitively reflect the influence of ANM schemes, Fig. 12 gives the results of the ANM operation strategies in the 24-h period under each of the three scenarios in Case 1. The power variation curves of electric demand, wind and photovoltaic power curves are presented in Fig. 12(a). The corresponding curtailment results are given in Fig. 12(b), (d) and (f), respectively. From Fig. 12(b), we can observe that the load-shedding is only adopted in extreme scenario ΦB when the electric demand attains peaks (16–20 h); while the load demand in scenarios ΦA and ΦC is relatively lower and no load-shedding is required in this circumstance. In addition, it can be noted from Fig. 12(d) as well as Fig. 12(f) that DGs are only curtailed in extreme scenario ΦA during the period of high-level DG output (16–20 h for WTG and 9–15 h for PVG). In comparison, since the power output is low in scenarios ΦB and ΦC, there is no curtailment and the renewable energy can be utilized fully. Fig. 12(c), (e) and (g) show the changes of VR tap position, SVC output and CB operating status with the variation of time stages respectively. Since the VR is installed closed to substation, the voltage can be controlled to adapt the changes of demand. If the total DG power output is lower than the load demand, especially in scenario ΦB, the overall VR tap is modulated in a lower position than that in ΦA and ΦB in order to limit the voltage drops, and restrict the energy loss over feeders. Furthermore, the reactive power output of SVC unit is kept
basic stability except for 20–24 h in ΦA since the high proportion DG power injection. Similarly, the operating CB units are the most during scenario ΦB, while the CB compensation during scenario ΦA is in the lowest level. Finally, according to Fig. 12(h), it can be briefly concluded that EES is charged during off-peak electric demand periods or on-peak DG output periods and discharged during on-peak demand periods, which can effectively smooth the power fluctuation caused by the injection of fluctuating loads or DG power. In detail, (i) in scenario ΦA, EES is continuously charged from 8 h to 16 h (during the on-peak DG output period) and discharged from 17 h to 24 h (during the on-peak power consumption period); (ii) in scenario ΦB, ESS is charged during the period from 1 h to 6 h (in the night when power consumption is low) and discharged from17h to 20 h (i.e., when power consumption attains peaks); however, due to the low DG output and high load in scenario ΦB, charging and discharging of ESS are not as obvious as scenarios ΦA and ΦC; (iii) in scenario ΦC, EES is charged during the period 10–15 h (along with large PVG output) and discharged during the period 17–24 h (along with the peak power consumption). 8. Conclusion This work proposes a new expansion planning model for multi-energy system where active distribution network, natural gas network and energy hub are incorporated. In this paper, we have shown how this comprehensive expansion planning problem can be modeled and formulated. This paper addresses the complex mixed integer nonlinear programming problem by applying second order cone relaxation as well as a modified piecewise linearization approach. With the objective of minimizing the total cost over the planning horizon, the final planning solution identifies the optimal determination of the type, location and size of all infrastructures, and hence their optimal operation strategy can be obtained simultaneously. Numerical computational results and analysis are carried out based on a multi-energy system with the integration of a modified IEEE33-node system, a 23-node gas system and 8 energy hub units. 21
Applied Energy 253 (2019) 113517
J. Wang, et al.
In order to conduct the simulations, a probabilistic scenario generation approach is proposed with the consideration of extreme cases. The proposed method has been verified to provide an efficient and robust scenario set for the expansion planning model. Compared with usual scenarios, the probabilistic scenario set increases the investment opportunities and reduces the energy consumption costs greatly. Additionally, the results based on different scenarios show the dominant impact of extreme scenarios on the planning decisions. Since there is usually a trade-off between the investment costs and the operation performance, several methods of scenario generation in this paper allow the planners to make a reasonable choice according to the actual situation. Moreover, the comparison of the comprehensive model and separate models has also been conducted in the case studies, and the results show the great superiority of the proposed joint expansion planning model over the separate models in the economic aspect. In addition, the impact of active network management on expansion planning of multi-energy system is investigated originally. The numerical results illustrate the superiority of considering the impact of active network managements on multi-energy system expansion planning. It has been proven that the full consideration of active network management can reduce the annual operation cost as well as the total planning cost. Given the fact that the proposed model with active network management enjoys more economic value, the incorporation of active distribution network in the joint multi-energy system expansion planning is indispensable. In summary, we believe the multi-energy system with consideration of active network managements in line with the emerging trend of active distribution network will be a promising research direction in the future.
[10] Rahmani-Andebili M. Distributed generation placement planning modeling feeder's failure rate and customer’s load type. IEEE Trans Ind Electron 2016;63(3):1598–606. [11] Rahmani-Andebili M. Simultaneous placement of DG and capacitor in distribution network. Electr Power Syst Res 2016;131:1–10. [12] Mokryani G, Hu YF, Papadopoulos P, Niknam T, Aghaei J. Deterministic approach for active distribution networks planning with high penetration of wind and solar power. Renew Energy 2017;113:942–51. [13] Koutsoukis NC, Georgilakis PS, Hatziargyriou ND. Multistage coordinated planning of active distribution networks. IEEE Trans Power Syst 2017;33(1):32–44. [14] Xie S, Hu Z, Zhou D, Li Y, Kong S, Lin W, et al. Multi-objective active distribution networks expansion planning by scenario-based stochastic programming considering uncertain and random weight of network. Appl Energy 2018;219:207–25. [15] Wu M, Kou L, Hou X, Ji Y, Xu B, Gao H. A bi-level robust planning model for active distribution networks considering uncertainties of renewable energies. Int J Electr Power Energy Syst 2019;105:814–22. [16] Moradijoz M, Moghaddam MP, Haghifam MR. A flexible active distribution system expansion planning model: a risk-based approach. Energy. 2018;145:442–57. [17] Rahmani-Andebili M. Optimal power factor for optimally located and sized solar parking lots applying quantum annealing. IET Gener Transm Distrib 2016;10(10):2538–47. [18] Zhang S, Cheng H, Wang D, Zhang L, Li F, Yao L. Distributed generation planning in active distribution network considering demand side management and network reconfiguration. Appl Energy 2018;228:1921–36. [19] Zeng B, Wen J, Shi J, Zhang J, Zhang Y. A multi-level approach to active distribution system planning for efficient renewable energy harvesting in a deregulated environment. Energy. 2016;96:614–24. [20] Shen X, Shahidehpour M, Han Y, Zhu S, Zheng J. Expansion planning of active distribution networks with centralized and distributed energy storage systems. IEEE Trans Sustain Energy 2017;8(1):126–34. [21] Xing H, Cheng H, Zhang Y, Zeng P. Active distribution network expansion planning integrating dispersed energy storage systems. IET Gener Transm Distrib 2016;10(3):638–44. [22] Andebili MR, Fotuhi-Firuzabad M. A robust approach for PEVs charging management and reconfiguration of electrical distribution system penetrated by renewables. IEEE Trans Ind Inf 2017;PP(99):1. [23] Rahmani-Andebili M. Dynamic and adaptive reconfiguration of electrical distribution system including renewables applying stochastic model predictive control. IET Gener Transm Distrib 2017;11(16):3912–21. [24] Rahmani-Andebili M. Stochastic, adaptive, and dynamic control of energy storage systems integrated with renewable energy sources for power loss minimization. Renew Energy 2017;113:1462–71. [25] Dehghan S, Amjady N, Kazemi A. Two-stage robust generation expansion planning: a mixed integer linear programming model. IEEE Trans Power Syst 2014;29(2):584–97. [26] Muñoz-Delgado G, Contreras J, Arroyo JM. Multistage generation and network expansion planning in distribution systems considering uncertainty and reliability. IEEE Trans Power Syst 2016;31(5):3715–28. [27] Bai L, Jiang T, Li F, Chen H, Li X. Distributed energy storage planning in soft open point based active distribution networks incorporating network reconfiguration and DG reactive power capability. Appl Energy 2017;210. [28] Markowitz HM, Manne AS. On the solution of discrete programming problems. Econometrica 1957;25(1):84–110. [29] Álvarez-Miranda E, Campos-Valdés C, Rahmann C. Two-stage robust UC including a novel scenario-based uncertainty model for wind power applications. Energy Convers Manage 2015;101:94–105. [30] Rahmani-Andebili M. Reliability and economic-driven switchable capacitor placement in distribution network. Gen Transm Distrib Iet 2015;9(13):1572–9. [31] Sheikhi A, Bahrami S, Ranjbar AM. An autonomous demand response program for electricity and natural gas networks in smart energy hubs. Energy 2015;89:490–9. [32] Wang C, Wei W, Wang J, Bi T. Convex optimization based adjustable robust dispatch for integrated electric-gas systems considering gas delivery priority. Appl Energy 2019;239:70–82. [33] Baran ME, Wu FF. Network reconfiguration in distribution systems for loss reduction and load balancing. IEEE Trans Power Delivery 1989;4(2):1401–7.
References [1] Vahid-Pakdel MJ, Nojavan S, Mohammadi-ivatloo B, Zare K. Stochastic optimization of energy hub operation with consideration of thermal energy market and demand response. Energy Convers Manage 2017;145:117–28. [2] Noorollahi Y, Yousefi H, Mohammadi M. Multi-criteria decision support system for wind farm site selection using GIS. Sustain Energy Technol Assess 2016;13:38–50. [3] Wang J, Zhong H, Ma Z, Xia Q, Kang C. Review and prospect of integrated demand response in the multi-energy system. Appl Energy 2017;202. [4] Mohammadi M, Noorollahi Y, Mohammadi-ivatloo B, Hosseinzadeh M, Yousefi H, Khorasani ST. Optimal management of energy hubs and smart energy hubs – a review. Renew Sustain Energy Rev 2018;89:33–50. [5] Beigvand SD, Abdi H, La Scala M. A general model for energy hub economic dispatch. Appl Energy 2017;190:1090–111. [6] Dagoumas AS, Koltsaklis NE. Review of models for integrating renewable energy in the generation expansion planning. Appl Energy 2019;242:1573–87. [7] Huang W, Zhang N, Yang J, Wang Y, Kang C. Optimal configuration planning of multi-energy systems considering distributed renewable energy. IEEE Trans Smart Grid 2018:1. [8] Fathtabar H, Barforoushi T, Shahabi M. Dynamic long-term expansion planning of generation resources and electric transmission network in multi-carrier energy systems. Int J Electr Power Energy Syst 2018;102:97–109. [9] Nazar MS, Haghifam MR. Multiobjective electric distribution system expansion planning using hybrid energy hub concept. Electr Power Syst Res 2009;79(6):899–911.
22