Multistage Expansion Co-Planning of Integrated Natural Gas and Electricity Distribution Systems

This paper focuses on expansion co-planning studies of natural gas and electricity distribution systems. The aim is to develop a mixed-integer linear programming (MILP) model for such problems to guarantee the finite convergence to optimality. To this end, at first the interconnection of electricity and natural gas networks at demand nodes is modelled by the concept of energy hub (EH). Then, mathematical model of expansion studies associated with the natural gas, electricity and EHs are extracted. The optimization models of these three expansion studies incorporate investment and operation costs. Based on these separate planning problems, which are all in the form of mixed-integer nonlinear programming (MINLP), joint expansion model of multi-carrier energy distribution system is attained and linearized to form a MILP optimization formulation. The presented optimization framework is illustratively applied to an energy distribution network and the results are discussed.


Introduction
In the era of modern energy systems, electricity and natural gas networks are becoming more and more interdependent [1][2][3][4]. It is therefore essential to run integrated studies of electricity and gas systems to improve the overall efficiency of the energy system. Consequently, planning and operation studies of integrated gas and electricity systems have been of great interest in new publications [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16]. In Reference [2], an optimization model is proposed for multi-time period combined gas and electricity networks aimed at minimizing costs associated with gas supply, storage and line pack management as well as electricity generation and load shedding. Stochastic models based on the Monte Carlo simulation are developed in Reference [3,4] to perform day-ahead scheduling and obtain energy flow solution for integrated natural gas and electricity power systems, respectively. In Reference [5], a two-stage stochastic look-ahead dispatch is proposed for integrated electricity power and gas systems which can adjust the improper day-ahead dispatch plan (Stage 1) in the second stage. A datamining technique called decision tree (DT) is utilized in Reference [6] to develop a systemic security dispatch method for integrated natural gas and power systems. In Reference [7], a multi-area, multi-stage model is proposed to integrate the long-term expansion planning studies of generation and transmission of natural gas and electricity systems. An optimization model is developed in Reference [8] to reach optimal expansion plans for gas power plants, electricity transmission lines and gas pipelines in terms • Developing a MILP model for co-expansion planning of EHs, natural gas and electricity distribution grids • Modelling the impacts of autonomous DG units on natural gas and electricity distribution networks using the EH concept • Deriving a novel mixed-integer linear formulation for natural gas distribution system expansion planning • Taking into account the effect of recovered heat from gas-fired DG units on heat demand reduction.

General Structure of the Proposed Framework
The general structure of the proposed framework is illustrated in Figure 1. According to this figure, the problem modelling starts with the formulation of three general parts of multi-carrier energy distribution systems, that is, energy hub, natural gas distribution system and electricity distribution network expansion planning problems. Subsequently, the co-expansion planning model is extracted from the combination of these three sets of equations. As will be shown, this model is a MINLP optimization problem for which obtaining a global optimal solution could not be guaranteed. Hence, in the next step, some accurate linear approximations are employed to handle nonlinearities and derive a MILP model. Finally, this MILP problem can be solved by standard mathematical programming optimization algorithms to acquire the global optimal expansion plan. • Taking into account the effect of recovered heat from gas-fired DG units on heat demand reduction.

General Structure of the Proposed Framework
The general structure of the proposed framework is illustrated in Figure 1. According to this figure, the problem modelling starts with the formulation of three general parts of multi-carrier energy distribution systems, that is, energy hub, natural gas distribution system and electricity distribution network expansion planning problems. Subsequently, the co-expansion planning model is extracted from the combination of these three sets of equations. As will be shown, this model is a MINLP optimization problem for which obtaining a global optimal solution could not be guaranteed. Hence, in the next step, some accurate linear approximations are employed to handle nonlinearities and derive a MILP model. Finally, this MILP problem can be solved by standard mathematical programming optimization algorithms to acquire the global optimal expansion plan.

Expansion Planning of Energy Hubs
As stated previously, in this paper, demand of energy carriers in different nodes of electricity and natural gas distribution networks are modelled via energy hubs. This allows for accurate modelling of the interaction between natural gas and electricity infrastructures at the load points [17,18]. General structure of the proposed energy hub is depicted in Figure 2.

Expansion Planning of Energy Hubs
As stated previously, in this paper, demand of energy carriers in different nodes of electricity and natural gas distribution networks are modelled via energy hubs. This allows for accurate modelling of the interaction between natural gas and electricity infrastructures at the load points [17,18]. General structure of the proposed energy hub is depicted in Figure 2. • Taking into account the effect of recovered heat from gas-fired DG units on heat demand reduction.

General Structure of the Proposed Framework
The general structure of the proposed framework is illustrated in Figure 1. According to this figure, the problem modelling starts with the formulation of three general parts of multi-carrier energy distribution systems, that is, energy hub, natural gas distribution system and electricity distribution network expansion planning problems. Subsequently, the co-expansion planning model is extracted from the combination of these three sets of equations. As will be shown, this model is a MINLP optimization problem for which obtaining a global optimal solution could not be guaranteed. Hence, in the next step, some accurate linear approximations are employed to handle nonlinearities and derive a MILP model. Finally, this MILP problem can be solved by standard mathematical programming optimization algorithms to acquire the global optimal expansion plan.

Expansion Planning of Energy Hubs
As stated previously, in this paper, demand of energy carriers in different nodes of electricity and natural gas distribution networks are modelled via energy hubs. This allows for accurate modelling of the interaction between natural gas and electricity infrastructures at the load points [17,18]. General structure of the proposed energy hub is depicted in Figure 2.  Input energy carriers of the energy hub are considered as grid electricity and the natural gas. Received natural gas from the network is dispatched between CHP unit and furnace based on a dispatch factor (υ n,t,ll ) [19]. Then, the CHP unit would supply a portion of electricity and heat demand [19].
In EH planning, the objective is to find optimal size of EH components, that is, transformer, CHP unit and furnace to efficiently meet the electricity and heat demand at the output layer. In this respect, the EH cost function which is comprised of the present value of the total investment and operating costs of EHs is given by Equation (1.a), where, N h and N L are number of EHs and load levels, respectively.
Du t,ll P Elec n,t,ll Pr Elec t,ll + P Gas n,t,ll Pr Gas t,ll + η Tra P Elec n,t,ll OC Tra +η CHP ge υ n,t,ll P Gas n,t,ll OC CHP +η Fur (1 − υ n,t,ll )P Gas n,t,ll OC Fur (1.c) η Tra P Elec n,t,ll + η CHP ge υ n,t,ll P Gas n,t,ll = D Elec n,t,ll η Fur (1 − υ n,t,ll )P Gas n,t,ll + η CHP gh υ n,t,ll P Gas n,t,ll = D Heat n,t,ll η CHP ge υ n,t,ll P Gas n,t,ll ≤ Ca CHP n,t  η Tra P Elec n,t,ll ≤ Ca Tra Ca CHP n,t = Ca CHP n,t−1 + NCa CHP n,t Ca Tra n,t = Ca Tra n,t−1 + NCa Tra In (1.a), the investment cost at stage t, represents the cost of newly installed equipment, that is, CHP units, furnaces and transformers. Moreover, the operating cost includes the cost of purchasing electricity and natural gas from the grids and also operating cost of EHs components.
This model is also subject to some constraints including energy balance in the EHs and capacity limitations of the EHs equipment [20] which are expressed as Equations (2)- (6). Constraints (2), (3) represent electricity and natural gas power balance in each EH. Constraints (4.a)-(4.c) set the bounds on the output power of EHs' components to their capacity and (4.d) is to limit penetration of CHP units within the electricity distribution network [21][22][23][24]. Capacity of EHs' components in each time stage is determined based on their capacity in previous time stage and newly added capacity using (5.a)-(5.c). Furthermore, the value of dispatch factors are limited by (6).
Finally, it is worth noting that this model is nonlinear due to the product of two decision variables, that is, dispatch factor (υ n,t,ll ) and natural gas power from grid (P Gas n,t,ll ), in (1.b)-(4.b).

Electricity Distribution Network Expansion Planning
Traditionally, this problem is solved to determine the optimal sets of feeders and substations which should be constructed or reinforced in anticipation of new electricity demand [24][25][26][27][28]. Mathematical model of this problem is presented in Equation (7.a). According to (7.b), the investment cost includes reinforcement cost of existing feeders and substations as well as construction cost of new feeders and substations. Various alternatives are considered for reinforcement and construction of each candidate (.),k,t ) is associated with each network component, that is, a feeder or substation, which indicates whether it is in-service in stage t (by taking value 1) or not. Subsequently, network operating cost is calculated as the sum of operating cost of all in-service components using (7.c).
Energies 2019, 12, 1020 6 of 16 This optimization problem is also subject to a variety of technical constraints which are formulated by (8.a)-(17.h). The first set of constraints are related to the availability of network facilities, that is, feeders and substations, which are expressed as (8.a)-(9.c). Constraints (8.a), (9.a) guarantee that by investing in any reinforcement alternatives of a component, utilization of its initial state is not allowable anymore [28]. Constraints (8.b), (8.c) and (9.b), (9.c) ensure that a candidate alternative can be utilized in network operation only after its associated investment is made [24,29].
As presented in Reference [24,[27][28][29], constraints (10.a)-(11.b) express that a maximum of one investment is allowed to be performed on each component during the planning horizon.
Kirchhoff's current law (KCL) at demand and substation nodes are imposed by (12.a) and (12.b), respectively. It is worth noting that V l p,t,ll and I b,t,ll are complex voltage and current phasors, respectively.
Constraints (13.a)-(13.d) represent the Kirchhoff's voltage law (KVL) [28]. Note that according to these equations, in case a branch is in-service (i.e., its binary utilization variable ϕ is 1), its relevant voltage drop constraint becomes active, otherwise the constraint would be relaxed.
Flow limits of feeders and substations are modelled by (14.a)-(15.c). As can be inferred, flow limit of an in-service component is the maximum capacity of its utilized alternative, being zero otherwise. Moreover, constraint (16) sets the upper and lower bounds of voltages at different nodes.
At last, radial operation of distribution networks should also be taken into account in the model. In fact, although electricity distribution networks are designed as meshed grids, they are still operated radially [24][25][26][27][28]. In this respect, radiality constraint of the network is modelled by (17.a)-(17.c). These equations are based on a graph theory which states that the number of nodes of a forest (herein, whole Energies 2019, 12, 1020 7 of 16 distribution network) is equal to the number of branches (in-service feeders) plus the number of trees (number of separate parts of the network which is equal to the number of substations) that are within the forest [30]. Accordingly, this set of equations simply indicates that at a given stage t, the number of in-service buses must be equal to the sum of in-service feeders and substations [31]. In this respect, the mode of network nodes are determined by (17.a) and (17.b). These equations force LPM lp,t to be one when node lp is in-service and set it to zero otherwise. It is worth mentioning that a bus is considered in-service as long as it is connected to at least one in-service feeder.
Although, this set of equations is sufficient condition for radiality of passive distribution networks, it is only necessary condition in case of active distribution networks. Hence, in order to ensure the radiality of an active distribution network, set of constraints (17.d)-(17.h) representing a fictitious current flow should also be considered in the model, as described in Reference [32].
Finally, it should be noted that all the equations of the electricity distribution planning model are mixed-integer linear, except for the ones associated with Kirchhoff's laws.

Natural Gas Distribution Network Expansion Planning
This problem intends to determine the optimal investment plan for expansion of natural gas network, that is, installing new city gates and pipelines or reinforcing the existing ones, to serve the predicted heat loads [14]. Hence, by considering the city gates and pipelines of the gas network equivalent to substations and feeders of electricity network, the general form of this problem is almost similar to the electricity distribution network planning. Therefore, the mathematical model of gas distribution expansion planning can be formulated as (18.a). Similar to electricity network expansion planning problem, the binary decision variables associated with investment and utilization (i.e., σ  ∑ b∈ l p χ Gas l p,b f b,t,ll = αGD l p,t,ll ; ∀ lp ∈ Ξ D (19.a) ∑ b∈ l p χ Gas l p,b f b,t,ll = α GD l p,t,ll − ∑ c∈O ξ Gas Energies 2019, 12, 1020 8 of 16 Moreover, constraints regarding the nodal gas balance as well as gas flow through pipelines [8]

Co-Expansion Formulation
The overall expansion planning model is extracted from the combination of the three aforementioned models by using the following equations: GD l p,t,ll = P Gas n,t,ll ED l p,t,ll = P Elec n,t,ll p f n,t,ll ∠ cos −1 (p f n,t,ll ) (24.c) As previously mentioned, there is no guarantee to acquire the global optimal solution for this optimization problem due to nonlinearities. Hence, in the next section, some modifications will be applied to linearize the model.

Linearization of EHs Planning Model
The source of nonlinearities in EH planning model is the product of dispatch factor (υ n,t,ll ) and natural gas power from grid (P Gas n,t,ll ). By introducing a new variable, NuP Gas n,t,ll and substituting all the υ n,t,ll P Gas n,t,ll terms with it and also replacing Equation (6) with Equation (25)

Linearization of Electricity Distribution Planning Model
As mentioned before, this model is nonlinear due to constraints addressed in (12.a)-(13.d). In order to linearize these constraints, the approximate model which was initially proposed in Reference [28] and then successfully used in Reference [24,27] is employed. According to this approximation, the complex nodal voltages, line currents, apparent powers and line impedances can be replaced by their absolute values [28]. This allows for elimination of nonlinear complex algebraic manipulations and results in linearization of the model. Moreover, owing to the slight variation of nodal voltages, the nominal bus Moreover, constraint (24.c) becomes: ED l p,t,ll = P Elec n,t,ll p f n,t,ll

Linearization of Natural Gas Distribution Network Expansion Planning
This model is nonlinear due to the gas flow equations through pipelines, that is, (20.a)-(20.d). Since the nonlinear terms of these equations are quite similar, in the following, the linearization process of (20.a) is explained. In this way, at first, in order to eliminate the square value of nodal pressures a new set of variables are introduced, as follows: π l p,t,ll = λ 2 l p,t,ll Then, by introducing some new positive variables and auxiliary binary variables to the model, the absolute value calculation underneath the radical sign can be omitted: note that since ∆π b,t,ll is a positive variable, only one of the constraints (31.c), (31.d) can be feasible for a given in-service pipeline b, depending on the sign of the result of ∑ l p∈Ξ χ Gas l p,b π l p,t,ll . For instance, for an in-service branch b (i.e., ϕ Pi,Fi b,t = 1), constraints (31.e) together with (31.f) express that one of the auxiliary binary variablesφ Pi,Fi b,t ,φ Pi,Fi b,t must be 1. Hence, if the outcome of ∑ l p∈Ξ χ Gas l p,b π l p,t,ll is positive, for example,φ Pi,Fi b,t is forced to 1 which indicates that the flow is in the predetermined direction and (31.a) becomes active so that the value of the branch flow would be assigned tof b,t,ll . In contrast, when ∑ l p∈Ξ χ Gas l p,b π l p,t,ll is negative,φ Pi,Fi b,t takes value 1. Hence, (31.b) would be active andf b,t,ll represents the flow of branch b.
Consequently, Equation (21.a) must be replaced by the following equations: Finally, a piecewise linear approximation can be employed to accurately calculate the square root of the pressure drop through pipelines. By doing so, Equations (31.a,b) are rewritten as follows: where, N pln is the number of blocks of piecewise linear pipeline flow function. Now, we have reached a linear model for co-planning study of a multicarrier energy system which can be efficiently solved using available commercial optimization software packages.

Case Study
In this section, the proposed co-expansion planning model is applied to a typical multicarrier energy distribution system. As shown in Figure 3, this energy system comprises electricity and natural gas networks, each consisting of 18 nodes: 16 load nodes and 2 supplying nodes (i.e., substation node for electricity network and city gate node for natural gas network). Moreover, two and three alternatives are considered for reinforcement and construction of candidate branches, respectively. The gas flow equations are also approximated by a piecewise linear function comprising four segments. As in Reference [28], the planning horizon is four years which is divided to three stages: duration of the first two is one year, while the third stage is two-year long. Furthermore, three load levels are considered as peak, medium and low load periods, with durations of 1000, 5760 and 2000 h, respectively. The whole model was implemented in General Algebraic Modelling System (GAMS) environment and solved by CPLEX solver, version 12.6.3 on a Fujitsu Celsius W530 POWER with a Quad 3.30 GHz Intel Xenon E3-1230 processor and 32 GB of RAM. The optimal gap tolerance was also set to 1%. It is worth mentioning that all the input data are available from [33].
Running the proposed MILP model of co-expansion studies, the optimal network topologies depicted in Figure 4, are obtained. In this figure, indices of the selected alternatives are placed alongside the candidate branches. For instance, label of the branch between nodes 13 and 17, at the first stage of electricity network, implies that the first alternative for construction of this feeder should be chosen. Accordingly, four feeders must be added at the first stage and four new feeders are required at the second stage. Moreover, the natural gas network needs five and four new pipelines at the first and second stages, respectively. Note that the dashed lines between nodes 17 and 1 at the second stage and nodes 1 and 5 at the third stage of gas network represent the unused pipelines to satisfy constraints  Running the proposed MILP model of co-expansion studies, the optimal network topologies depicted in Figure 4, are obtained. In this figure, indices of the selected alternatives are placed alongside the candidate branches. For instance, label of the branch between nodes 13 and 17, at the first stage of electricity network, implies that the first alternative for construction of this feeder should be chosen. Accordingly, four feeders must be added at the first stage and four new feeders are required at the second stage. Moreover, the natural gas network needs five and four new pipelines at the first and second stages, respectively. Note that the dashed lines between nodes 17 and 1 at the second stage and nodes 1 and 5 at the third stage of gas network represent the unused pipelines to satisfy constraints regarding the radial operation of natural gas distribution network. The grey-shaded nodes indicate the load points at which CHP units are installed. Finally, it is worth noting that no pipeline or feeder reinforcements are required.     Running the proposed MILP model of co-expansion studies, the optimal network topologies depicted in Figure 4, are obtained. In this figure, indices of the selected alternatives are placed alongside the candidate branches. For instance, label of the branch between nodes 13 and 17, at the first stage of electricity network, implies that the first alternative for construction of this feeder should be chosen. Accordingly, four feeders must be added at the first stage and four new feeders are required at the second stage. Moreover, the natural gas network needs five and four new pipelines at the first and second stages, respectively. Note that the dashed lines between nodes 17 and 1 at the second stage and nodes 1 and 5 at the third stage of gas network represent the unused pipelines to satisfy constraints regarding the radial operation of natural gas distribution network. The grey-shaded nodes indicate the load points at which CHP units are installed. Finally, it is worth noting that no pipeline or feeder reinforcements are required.  The investment and operating costs associated with optimal co-expansion planning studies of Case I and separate expansion planning of Case II are also illustrated in Table 1. According to this table, Case II in which the EH, electricity and gas distribution expansion planning studies are carried out separately, results in a lower EH cost than Case I. However, in case of co-expansion planning, this increment is dominated with the networks expansion costs reduction so that a lower cost solution is obtained, despite the fact that the values of peak demand for both networks have not changed. It is Energies 2019, 12, 1020 12 of 16 worth noting that in both cases, due to the cost-effectiveness, total capacity of installed CHP units is at the maximum value determined by (4.d).
Nevertheless, as can be seen from the table, the selected load points for installation of CHP units are considerably different. Moreover, as illustrated in Table 1, in collaborative expansion planning mode (Case I), the energy production of CHP units as well as total energy demand of gas network are reduced. In fact, in separate planning regime (Case II), EH owner tends to install CHP units at higher number of load points in order to locally serve heat demands at lower costs. However, it causes the injection of electricity power generated by CHP units at more diverse nodes, which calls for providing more electricity network capacity to handle it. At the same time, in comparison with Case I, more investment should be made in natural gas network to be able to serve increased gas demand (owing to the natural gas consumption of CHP units) at such nodes. Therefore, optimizing the EHs objective, without considering its impact on the distribution infrastructures, results in the increased investment costs of both electricity and natural gas networks. However, comparing the obtained results of the investigated cases, it can be inferred that this is not efficient from the viewpoint of social welfare.
Finally, as can be expected, solving the integrated planning model (Case I) needs considerably more running time as seen in Table 1. The full test system database and detailed results are also available from [33].

Conclusions
In this paper, a new framework has been proposed for expansion co-planning studies of natural gas and electricity distribution networks. In this regard, the concept of energy hub is employed to effectively model the interactions of different energy carriers. Extracting the mathematical model of co-planning problem, an effective algorithm has been introduced to reach a MILP model for this problem. The proposed framework was applied to a multi-carrier energy test system and the achieved results were illustrated and compared with traditional separate planning studies. Lowering the required investment for expansion of distribution networks and optimizing the future vision of EHs based on operating conditions of energy networks are superiorities of the co-planning studies compared to separate planning studies.     Element of node-branch incidence matrix for electricity and gas networks which is −1 or +1 if branch b is connected to load point lp and the predetermined current or flow direction is toward or away from node lp, respectively and is 0 otherwise. δ t,Inv , δ t,Op Present value factors for investment and operating costs.