Sustainable Cities and Society 97 (2023) 104787

Contents lists available at ScienceDirect

Sustainable Cities and Society

journal homepage: [Link]/locate/scs

A Convex Micro-Grid-Based Optimization Model for Planning of Resilient

and Sustainable Distribution Systems Considering Feeders Routing and
Siting/Sizing of Substations and DG Units
Saeed Behzadi, Amir Bagheri *
Department of Electrical Engineering, University of Zanjan, Zanjan, Iran

A R T I C L E I N F O A B S T R A C T

Keywords: Integration of distributed generation (DG) and renewable energy resources (DERs) has affected planning and
Distribution system expansion planning operation of recent distribution networks (DNs). On the other side, resiliency and sustainability of DNs is the
Micro-grid formation other concern of distribution network operators. In this paper, a new formulation is developed for optimal
Resiliency
planning of resilient DNs based on optimal formation of resilient micro-grids (MGs). The formulation aims to
Feeder routing
DER
implement optimal siting and sizing of conventional and renewable-based DGs, optimal routing and type of
Emission cost feeders, as well as optimal sizing and placement of HV/MV substations to configure adequate and resilient MGs
Uncertainty against extreme events. The investment, operational, emissions, and resiliency costs have been included in the
objective function, and all the problem constraints as well as uncertainty of renewable DGs have been taken into
account. The new formulation uses line-flow-based (LFB) model of AC power flow equations, and reserve feeders
(tie-lines) are employed to enhance resiliency after severe outages via reconfiguration concept. All the relations
have been convexfied in order to compose a mixed-integer quadratically-constrained (MIQCP) model to be
solved with global optimum solvers in GAMS. Efficiency of the conducted methodology has been evaluated by
different experiments on the 24-node system, and the results are investigated.

2015, Gilani et al., 2022, Galvani, Bagheri, Farhadi-Kangarlu, and

1. Introduction Nikdel, 2022). In order to better exploit the advantages of DG units,
nowadays, the distribution systems are configured as different
Distribution systems play an important role in supplying the final micro-grids containing several loads and DER units. By forming optimal
customers from the transmission grids. Along with the increase of cus­ micro-grids, installation and operation costs of DNs are reduced. This is
tomers’ load demand, the facilities of distribution networks lose their due to supplying the loads via nearby sources instead of investing on
adequacy in appropriate and reliable feeding of load centers (Sallam and upward grid reinforcement. Management of micro-grids is more
Malik, 2018). Due to radial structure of DNs, the distribution network convenient than an integrated large network (Ebadat-Parast, Nazari, and
operators (DNOs) are always engaged with high voltage drop and power Hosseinian, 2022). Moreover, in the emergency conditions, the faulted
losses, and these problems are intensified when the load growth occurs MGs can be isolated form the remaining ones in order to reduce the fault
(Bayat and Bagheri, 2019). Therefore, it is required to enhance the effect and economic damages, and increase the electric grid’s sustain­
distribution system adequacy by doing expansion planning studies. ability. In addition, performance of distribution network under normal
Conventionally, the planning is fulfilled either by reinforcing the ca­ operating condition is also improved such that there would be lower
pacity of existing lines and transformers or by installing new ones to power loss and voltage drop within the system (Dada, 2014).
cover the load demand growth (Seifi and Sadegh Sepasian, 2011). By On the other side, in recent years, due to climate changes, natural
emerging DG units in today’s electric grids, the passive DNs have been disasters with low probability and high impact (LPHI) such as floods,
transformed to active ones to benefit from different aspects of DERs. In hurricanes, earthquakes, ice storms, dust storms, etc. have been inten­
this regard, the renewable-energy-based DGs such as wind and photo­ sified in many countries (Abessi, Jadid, and Salama, 2020, Rajabzadeh
voltaic (PV) units have attracted the DNOs because of their fuel-free and and Kalantar, 2022). This implies that today’s DNs should be resilient
non-polluting nature (Lopes et al., 2007, Bagheri, Monsef, and Lesani, against these events, and the resiliency and sustainability should be

* Corresponding Author.
E-mail address: amir_bagheri@[Link] (A. Bagheri).

[Link]
Received 24 March 2023; Received in revised form 16 June 2023; Accepted 7 July 2023
Available online 8 July 2023
2210-6707/© 2023 Elsevier Ltd. All rights reserved.
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

Nomenclature PF Generation power factor of synchronous DG units

DGMax
Num Maximum number of DG units that can be installed at each
Indices bus
i, j Buses DGMax
IN Maximum number of DG units that can be allocated in the
LL Load level whole network
m Micro-grids Max
WTNum Maximum number of WT units that can be installed at each
s Uncertainty scenario bus
w Outage state Max
WTIN Maximum number of WT units that can be allocated in the
k Conductor type whole network
y Planning horizon year CRs WT units capacity rate under uncertainty scenario
capacity
Sets PWT Maximum active power capacity of WT units
ΩN Set of buses GenMax Maximum fictitious active power that can be injected into
ΩT Set of time the network
ΩK Set of conductors type SUBγ,Max Maximum number of existing substations
ΩMG Set of micro-grids SUBβ,Max Maximum number of new substations
ΩW Set of contingency states VOLLLL Value of lost load
ΩS Set of different uncertainty scenarios TF Time duration that equipment is out of service due to
ΩY Set of year severe events

Parameters Variables
LDi,j Line ij length Pnet net
i,j,LL,s,w /Qi,j,LL,s,w Active/reactive power passing through the lines
net a net a
Pnet net
Min /PMax Minimum/maximum active power passing through the Pi,j,LL,k,s,w /Qi,j,LL,k,s,w Active/reactive power passing through special
lines type of conductor
Qnet
Min / Q net
Max Minimum/maximum reactive power passing through PSUB SUB
i,LL,s,w /Qi,LL,s,w Active/reactive power injected into the system from
the lines the upstream network
PLi / QLi Peak value of active/reactive load demand PDG DG
i,LL,s,w /Qi,LL,s,w Active/reactive power injected into the system by
L,contingency
PTotal Total active load demand during contingency states synchronous generators
SMax Maximum capacity of substations PWT
i,LL,s,w Active power injected into the system by wind generators
λLL Load level factor STotal
i,LL,m,s,w Total apparent power injected to the network from all
RLine
i,j,k Line resistance substation
Line
Xi,j,k Line reactance Sγi,LL,m,s,w Apparent power injected to the network from existing
substation
V Min /V Max Minimum/maximum voltage of buses
Sαi,LL,m,s,w Apparent power injected to the network from expanded
UMin / UMax Square of minimum/maximum voltage of buses substation
bigM A big number
Sβi,LL,m,s,w Apparent power injected to the network from new
ICTrans Installation cost of transformers
substation
ICSub Installation cost of substations Pγi,LL,m,s,w /Qγi,LL,m,s,w Active/reactive power injected to network from
ICDG Installation cost of gas turbine (synchronous DG) existing substation
ICWT Installation cost of wind turbine (WT) Pαi,LL,m,s,w /Qαi,LL,m,s,w Active/reactive power injected to network from
PW Present worth factor expansion existing substation
TLL Duration of load level LL (hr)
Pβi,LL,m,s,w /Qβi,LL,m,s,w Active/reactive power injected to network from
Infr Inflation rate (%)
new substation
Intr Interest rate (%)
PLsh
i,LL,s,w /Q Lsh
i,LL,s,w Active/reactive power of load shedding
Probs Probability of each uncertainty scenario
PrSUB Price of energy purchase from the upstream network PLsh,MG Lsh,MG
i,LL,m,s,w /Qi,LL,m,s,w Active/reactive power of load shedding in micro-
LL
($/MWhr) grid
PrDG Operating cost of synchronous DG ($/MWhr) PLsh,Total
w Active power of load shedding during contingency state
emc Emission cost ($/ton) RIw System resilience index for each contingency state
GE Emission of CO2 related to the received power from the Transαi Number of expanded transformers in existing substation
transmission grid (tons/MWhr) Transβi Number of installed transformers in new substation
DGE Emission of CO2 related to the generated power by DGs OF Main objective function
(tons/MWhr) ICSUBα Expansion cost of existing substations
Transγ Number of transformers in existing substations ICSUBβ Installation cost of new substations
Transcapacity Capacity of each transformer ICLine Installation and expansion cost of lines
TransMax Maximum number of transformers that can be paralleled in ICDG,WT Installation cost of synchronous DGs and wind turbines
each substation OCSUB Cost of purchasing energy from upstream network
ωEx
i,j,k Existing type of conductor between buses i and j OCDG Total operating cost of synchronous DGs and wind turbines
Buy
Lk /LSell
k Cost of buying/selling of the feeders EMC Total emission cost due to pollutant CO2
IkAmp Lines ampacity in Ampere ENSC Energy not supplied cost
Vi,LL,s,w Voltage of buses
flowMin /flowMax Minimum/maximum fictitious active power that
Ui,LL,s,w Square of bus voltage
can be passed through lines
Ji,j,LL,s,w Square of line current passing through the lines
SMax
DG Capacity of synchronous generators a
Ji,j,LL,k,s,w Square of line current passing through each type of

2
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

conductor DGNum
i Integer variable to show number of DG units installation at
Xi,j
R
Binary variable of connection/disconnection of all lines bus i
(including main and tie lines) WTiIN Binary variable to allocate new WTs
Xi,j,w Binary variable of connection/disconnection of lines WTiNum Integer variable to show number of WT units installation at
MG
Xi,j,m,w Binary variable of connection/disconnection of lines in bus i
micro-grid ωNew
i,j,k Binary variable to show conductor type in line i-j after
BNew
i,j,k Binary variable of constructing new line expansion
Transαi Integer variable to show number of new transformer ωai,j,k,w Binary variable to show conductor type in line i-j after
installed in existing substations expansion for each contingency state
Transβi Integer variable to show number of installed transformers Zi,j,k Binary variable to show either constructing new feeder or
in new substations changing conductor of existing feeder
SUBγi,m Binary variable to show presence of existing substations in ai,m,w Binary variable indicating the presence of bus i in micro-
bus i grid m
neg pos neg
SUBαi,m Binary variable to show expansion of existing substations Apos
i,j,k /Ai,j,k ri,j,m,w /ri,j,m,w Positive/negative binary variables used for
linearization of absolute function
SUBβi,m Binary variable to allocate new substations
Geni,LL,m,w Fictitious injected active power
βwi,m,w Binary variable of new substations in contingency states flowi,j,LL,m,w Fictitious active power passing through lines
γwi,m,w Binary variable of existing substations in contingency
PLsh,island
i,LL,m,w Fictitious active power of load shedding in micro-grids
scenarios
MGNum RWi Optimal micro-grid number
SUBTotal
i,m Binary variable to show all substations
FWi,w ZWi,j,w Binary variables for contingency states
DGIN
i Binary variable to allocate new synchronous DGs

accounted as a major criterion of modern distribution systems (Gritza­ equipment outage. The considered constraints include voltage and
lis, Theocharidou, and Stergiopoulos, 2019, Haghshenas, Hooshmand, current limitations, substations loading limit, and the network radiality
and Gholipour, 2022). The resiliency improvement of DNs can be ful­ constraint. For the optimization purpose, the artificial immune system
filled through strengthening substations and feeders, installing reserve (AIS) as a meta-heuristic algorithm has been employed and imple­
feeders (tie-lines) and sectionalizers, and installing DG units in appro­ mented in MATLAB. Although both of (Bagheri, Monsef, and Lesani,
priate locations of the network. Therefore, the resiliency should be 2015) and (Seta, Oliveira, and Oliveira, 2019) considered almost all the
regarded in the expansion planning stage of distribution systems and problem aspects, their model is in the form of MINLP which is
formation stage of micro-grids such that the designed network could non-convex, and there is no guarantee of finding global optimum solu­
satisfy the load growth and recover the maximum possible load after the tion using the employed meta-heuristic algorithms. Ref. (Elkadeem
occurrence of disasters. Subsequently, the recent researches around et al., 2019) implements planning of distribution system considering
expansion planning of active distribution networks and micro-grids as integration of PV and WT units. The power loss reduction, voltage
well as their resiliency improvement are reviewed. profile improvement, and yearly economic saving are considered as the
Resiliency of a multi-infrastructure energy system (MIES), including problem objectives. The model of this paper is also MINLP, and hence,
natural gas and electricity networks is studied in (Ravadanegh, Jamali, the optimality concern still exists. In addition, this paper only considers
and Vaniar, 2022) under the effect of multi-hazard disasters such as DG units as the planning tool, and the reliability or resiliency aspects are
hurricane and earthquakes. The impacts of hazard levels on component not under consideration. A bi-level optimization methodology has been
and system performance are analyzed, and node resilience criteria are presented in (Prathapaneni and Detroja, 2019) for optimal planning and
employed for various energy sub-systems to evaluate the resilience total operation of micro-grid composed of different types of distributed en­
system. The results show that the simultaneous occurrence of the ergy resources (DERs). The planning decisions include installation of DG
different hazards on interdependent infrastructures can boost and units which is handled by accelerated particle swarm optimization
overlap the impact of a single hazard on the MIES. The resiliency index is (APSO) in MATLAB. Also, the operation stage consists of unit commit­
defined as restored energy which is covered by different DG units based ment problem of DGs as well as scheduling of demand side management
on the system topology. Distribution network expansion planning has fulfilled within GAMS. This paper has not considered the power flow
been formulated as an integrated mixed-integer non-linear program­ constraints, and the micro-grid is considered as a single-bus network to
ming (MINLP) model in (Bagheri, Monsef, and Lesani, 2015) to consider which all the DER units and loads are connected. In (Sun et al., 2021), a
most of the planning alternatives. The planning determines the optimal long-term distribution system expansion planning based on approximate
reinforcement of existing MV lines and HV/MV substations, or con­ dynamic programming (ADP) is proposed using both DG units and
struction of new facilities to cover the load growth in the planning ho­ network equipment reinforcement. The objective function considers
rizon while meeting technical and operational constraints. The investment and operational costs, generation cost of DGs, cost of energy
renewable and non-renewable DG units are considered, and the uncer­ purchased from the upward grid, and un-supplied energy cost. All the
tainty of renewable DGs has been regarded. The routing and type of relations are in convexified form to construct a mixed-integer linear
tie-lines (reserve feeders) are also optimized to restore maximum programming (MILP) model. This paper has not regarded the resiliency
possible loads in the emergency conditions. Finally, the genetic algo­ issues and outage event of system equipment. A two-stage long-term
rithm, as a meta-heuristic algorithm, solves the proposed MINLP model. expansion planning approach is proposed in (Saravi, Kalantar, and
A comprehensive MINLP model has been presented in (Seta, Oliveira, Anvari-Moghaddam, 2022) aiming at simultaneous optimal expansion
and Oliveira, 2019) for distribution network expansion planning of distribution networks and resilience improvement against hurricanes.
(DNEP). The objective function includes investment and operations In the first stage, the total cost composed of investment, operation, and
costs as well as CO2 emission and reliability ones. The planning tools are resilience costs is minimized, where the optimal expansion strategy
the same for (Bagheri, Monsef, and Lesani, 2015), but the along with optimal operation scheme is obtained. In the second stage, a
non-renewable DGs are not considered. The “N-1” reliability criterion stochastic model under hurricane occurrence scenarios is modeled for
has been taken into account to restore the maximum possible loads after the resiliency index maximization. The objectives are attained through

3
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

Table 1
Summary of recent literature review around distribution system expansion planning.
Ref. No (Year) Feeders Substations DGs Operating Condition Model

Hardening Routing Tie- Fixed MG-Based Non- Renewable Normal Resiliency Convex Non-
lines Planning (Sizing Renewable Operation Consideration Convex
and Placement of
Substations)

(Bagheri, Monsef, and ✓ ✓ ✓ - ✓ ✓ ✓ ✓ ✓ × ✓

Lesani, 2015)
(Seta, Oliveira, and ✓ ✓ ✓ - ✓ × ✓ ✓ ✓ × ✓
Oliveira, 2019)
(Elkadeem et al., 2019) × × × ✓ × × ✓ ✓ × × ✓
(Prathapaneni and × × × × × ✓ ✓ ✓ × × ✓
Detroja, 2019)
(Sun et al., 2021) ✓ ✓ × - ✓ × ✓ ✓ × ✓ -
(Saravi, Kalantar, and × × ✓ ✓ × ✓ ✓ ✓ ✓ ✓ -
Anvari-Moghaddam,
2022)
(Zakernezhad, Nazar, ✓ × ✓ ✓ × ✓ ✓ ✓ ✓ × ✓
Shafie-khah, and
Catalão, 2021)
(Borghei and Ghassemi, × × ✓ ✓ × ✓ ✓ ✓ ✓ ✓ -
2021)
(Shahbazi et al., 2021) ✓ × ✓ ✓ × ✓ × ✓ ✓ ✓ -
(Zhou et al., 2022) ✓ ✓ ✓ - ✓ × ✓ ✓ ✓ × ✓
(Saberi et al., 2023) ✓ ✓ ✓ ✓ × ✓ ✓ ✓ ✓ ✓ -
(DC-
PF)
(Bayani and Manshadi, ✓ ✓ ✓ - ✓ ✓ ✓ ✓ ✓ × ✓
2023)
Proposed model ✓ ✓ ✓ - ✓ ✓ ✓ ✓ ✓ ✓ -

optimal placement of different DER units and lines hardening as well as allocated to improve network resiliency in outage cases. The relations of
installation of tie-lines, but the feeders routing has not been considered, this work compose an MINLP model which is a non-convex one. An
and just the existing routes are reinforced. In (Zakernezhad, Nazar, expansion planning model is presented in (Saberi et al., 2023) to limit
Shafie-khah, and Catalão, 2021), a multi-stage optimization framework the wildfire ignition risk of electricity networks and model the long-term
is developed for expansion planning of resilient distribution networks at impacts of public safety power shut-off events. A trade-off between
the presence of utility and non-utility owned distributed energy re­ shutting off power to customers and grid expansion or hardening mea­
sources (UDERs and NDERs). The main and tie feeders are considered, sures has been established. The results verify that DER installation is the
but there is no placement of substations. The proposed approach is best approach to meet demand growth and inhibit wildfire risk in the
modeled as an MINLP one which is optimized by hybrid utilization of long run. The model does not account for HV/MV substation, and the
PSO algorithm (in MATLAB) and DICOPT solver (in GAMS). The simu­ power flow is based on DC model (DC-PF) which is not accurate
lation results show total cost reduction using the proposed methodology. compared to the AC one. The most recent research which considers
For strengthening distribution network against severe natural disasters, almost all the planning options for resilient distribution network is
Ref. (Borghei and Ghassemi, 2021) has addressed the expansion plan­ (Bayani and Manshadi, 2023) in which the substations, feeders, and
ning of MGs. The optimal nodes for connection of micro-grids and the renewable/non-renewable DG sources are the alternatives. The pro­
capacity of DG units are determined by using a heuristic approach as posed planning is solved using a multi-objective genetic algorithm (GA)
well as a multi-objective mixed-integer linear programming (MO-MILP) as a meta-heuristic approach. The results show that considering resil­
while meeting operational, topological, and power flow constraints. This ience in distribution networks changes the network design compared to
work does not consider planning of feeders and substations, and just the traditional plans and the distribution companies will be able to deter­
tie-lines are optimized. Ref. (Shahbazi et al., 2021) has focused on mine an acceptable level of distribution network resilience considering
planning of backup DGs, hardening of lines, and installation of tie-lines costs and benefits of resiliency improvement. The model of this paper is
in distribution networks in order to enhance system resiliency in front of an MINLP one which is solved by the GA, and hence, the results may not
natural disasters by minimizing investment, operation and costs. The be global optimum.
Benders decomposition (BD) is employed to convert the MINLP model of After investigating the recent works, a review summary can be
the problem into an MILP one to be optimized in GAMS using CPLEX expressed as Table 1. In the current paper, optimal planning of resilient
solver. Also, a hybrid stochastic/robust optimization (HSRO) and a micro-grid based DNs has been addressed using a new formulation. The
scenario-based stochastic optimization are used to model the un­ planning decisions include sizing and placement of conventional and
certainties of load, energy price and availability of the network equip­ renewable-based DG units, hardening and routing of main and reserve
ment under the extreme weather conditions. This paper does not feeders, and optimal sizing and placement of HV/MV substations in
consider MGs formation and substations planning. In (Zhou et al., 2022), order to configure adequate and resilient MGs. The proposed model is
optimal distribution network expansion planning is proposed regarding able to reconfigure the network into appropriate micro-grids in order to
MG concept in the presence of intermittent renewable energy sources prepare optimal normal operation as well as to withstand against
(RESs) to improve the topological flexibility in modern distribution possible outages during natural disasters. The investment, operational,
networks. The best locations of substations, lines, and RES-based DGs emission, and load shedding costs are included in the objective function,
are coordinately taken into consideration, and multiple uncertainty and all the problem constraints have been taken into consideration. All
scenarios are modeled to acquire robustness of the network against load of the problem relations have been convexfied in order compose a
and RES uncertain behavior. The graph theory is used to guarantee the MIQCP model; the model is implemented in GAMS and solved by global
radial topology of network, and the reserve feeders are optimally optimum solvers.

4
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

Fig. 1. A schematic overview of the considered items and contributions of the current paper.

In summary, the contributions of the current paper can be outlined as cost of DG units. As seen, the operational costs are converted to their
follows: present worth using PW which is yield based on inflation and interest
rates. As the air pollution has harmful effects on human and societies’
• New convex formulation is proposed for optimal sizing and alloca­ health, it is necessary to consider the effect of emissions in the planning.
tion of new substations as well as for expansion of existing On this basis, the emission cost has been formulated in (8), which
substations; depend on the emission of DG units as well as that incurred by the up­
• New convex formulation is proposed to form micro-grids and allo­ ward grid. In Eq. (9), unsupplied energy cost (ENSC) is presented, which
cate tie-lines; can be accounted as an indication of demand response to the probable
• New formulations for feeders routing and existing ones’ optimal fault occurrence. This cost is calculated based on the average values of
hardening are developed; energy not supplied in different substation outage states. All the vari­
• Resiliency is considered through several outage events in different ables and parameters have been defined in the nomenclature for
load levels; convenience.
• Emission reduction as well as uncertainty of wind generations is also
∀i, j ∈ ΩN , ∀LL ∈ ΩT , ∀k ∈ ΩK , ∀s ∈ ΩS , ∀m ∈ ΩMG , ∀w ∈ ΩW , ∀y ∈ ΩY
taken into account.
OF = ICSUBα + ICSUBβ + ICLine + ICDG,WT + OCSUB + OCDG + EMC + ENSC
Fig. 1 schematically illustrates different aspects of distribution sys­
(1)
tem planning which are addressed in this paper.
∑
ICSUBα = Transαi ICTrans (2)
2. Problem formulations i∈ΩN

∑ ∑ ∑
2.1. Objective function ICSUBβ = Transβi ICTrans + SUBβi,m ICSub (3)
i∈ΩN i∈ΩN m∈ΩMG
The objective function considered in this paper is given in (1), which
( (( ) ( )))
is composed of all investment and operation costs, emission cost of CO2, ∑∑ ∑ ∑
and cost of energy not supplied due to occurrence of the worst-case ICLine = LDi,j BNew Buy
i,j,k Lk − Zi,j,k ωEx Sell
i,j,k Lk (4)
disaster events. In (2)-(5), installation cost of new substations and
i∈ΩN j∈ΩN k∈ΩK k∈ΩK

expansion cost of existing ones, hardening cost of existing lines and ∑( ) ∑( )

ICDG,WT = DGIN
i IC
DG
+ WTiIN ICWT (5)
construction cost of new ones, and installation cost of DG units are i∈ΩN i∈ΩN
presented, respectively. Also, operation costs formulation under uncer­
tainty is shown in relations (6)-(7). Relation (6) denotes cost of pur­ OCSUB =
∑
(PW)y
∑∑ ∑∑
Total
Probs TLL Si,LL,m,s,w PrSUB (6)
LL ; ∀w = w1
chasing energy from the upstream network, and (7) gives the operation y∈ΩY i∈ΩN LL∈ΩT m∈ΩMG s∈ΩS

5
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

Table 2 ( )2 ( )2
Relationship between variables of (24). Pnet
i,j,LL,s,w + Qnet
i,j,LL,s,w = Ji,j,LL,s,w Uj,LL,s,w (14)
State No. ωNew
i,j,k ωEx
i,j,k
Zi,j,k
∑
1 0 0 0
Pnet net
i,j,LL,s,w = − Pj,i,LL,s,w = Pnet a
i,j,LL,k,s,w (15)
2 0 1 1
k∈ΩK

3 1 0 1 ∑
4 1 1 0 Qnet net
i,j,LL,s,w = − Qj,i,LL,s,w = Qnet a
i,j,LL,k,s,w (16)
k∈ΩK

∑ ∑∑∑ ωai,j,k,w Pnet net a a net

Min ≤ Pi,j,LL,k,s,w ≤ ωi,j,k,w PMax (17)
OCDG = (PW)y Probs TLL PDG
i,LL,s,w Pr
DG
; ∀w = w1 (7)
y∈ΩY i∈ΩN LL∈ΩT s∈ΩS
ωai,j,k,w Qnet net a a net
Min ≤ Qi,j,LL,k,s,w ≤ ωi,j,k,w QMax (18)
⎛( ) ⎞
∑∑∑ ∑
⎜ Probs TLL PSUB
i,LL,s,w GE ⎟ Ji,j,LL,s,w = a
Ji,j,LL,k,s,w (19)
∑⎜ i∈ΩN LL∈ΩT s∈ΩS ⎟
y ⎜ ⎟ k∈ΩK
EMC = emc (PW) ⎜ ( ) ⎟; ∀w = w1
⎜ ∑ ∑ ∑ ⎟
y∈ΩY ⎝ ⎠ ( )2
+ Probs TLL PDG
i,LL,s,w DGE
i∈ΩN LL∈ΩT s∈ΩS
a
Ji,j,LL,k,s,w ≤ ωai,j,k,w IkAmp (20)
(8) ( )2
Ui,LL,s,w = Vi,LL,s,w (21)

⎛ ⎞
⎜ ⎟∑ ∑ ∑ ∑∑
⎜ T ⎟
ENSC = ⎜
⎜ ( F )⎟
⎟ (PW)y Probs TLL PLsh
i,LL,s,w VOLLLL (9)
⎝ ∑ ⎠y∈ΩY i∈ΩN LL∈ΩT s∈ΩS w∈ΩW
365 w− 1
w∈ΩW

1 + Infr V Min ≤ Vi,LL,s,w ≤ V Max (22)

PW = (10)
1 + Intr
Eqs. (11) and (12) denote active/reactive power balance constraint
at each bus of the network in different load levels, scenarios, and
2.2. Problem constraints
operating mode. Aℓi is the ℓi − th element of the bus-line incidence
matrix. It is equal to 1 if bus i is the sending bus of line i-j. If bus i is the
2.2.1. Power flow equations
receiving bus of line i-j, Aℓi equals to -1. Otherwise, it will be zero. Also,
In distribution systems, the AC power flow relations are non-linear
Bℓi is the same as Aℓi by replacing 1 with 0. In (13), the voltage drop
and non-convex. In order to develop a convex optimization model, a
along a feeder connected between nodes i and j has been defined, where
new formation of line flow based (LFB) model, which considers selection
Xi,j,w is a binary variable indicating that the line between buses i and j in
of lines conductors’ type as well as reconfiguration of network using
contingency state of w is connected or not. Eq. (14) is the nodal rela­
main and reserve feeders, is proposed as follows.
tionship between power, voltage, and current. In (15) and (16), the
∑
Aℓi Pnet SUB DG WT Lsh L relation between active and reactive powers of sending and receiving
i,j,LL,s,w = Pi,LL,s,w + Pi,LL,s,w + Pi,LL,s,w + Pi,LL,s,w − λLL Pi
j∈ΩN buses can be seen. It should be noted that Pnet i,j,LL,s,w and Qi,j,LL,s,w are active
net

∑
( (
∑
)) (11) and reactive powers sent by node j toward node i at the side of node i.
Bℓi LDi,j RLine That is, the power loss of line i-j has not been considered in (15) and
a
− i,j,k Ji,j,LL,k,s,w
j∈ΩN k∈ΩK
(16). Relations (17) and (18) guarantee that the power passing through
∑ the disconnected line is zero. In the relations, ωai,j,k,w guarantees that after
Aℓi Qnet SUB DG Lsh L
i,j,LL,s,w = Qi,LL,s,w + Qi,LL,s,w + Qi,LL,s,w − λLL Qi
j∈ΩN
the expansion, each existing or newly installed line must have only one
( ) conductor type by considering emergency states for the network. The
∑ ∑
− D
Bℓi Li,j Line a
Xi,j,k Ji,j,LL,k,s,w (12) value of ωai,j,k,w is zero for the lines that are not installed and were not
j∈ΩN k∈ΩK present in the network before expansion. By considering ωai,j,k,w as a
variable, the type of line is selected by regarding the critical conditions,

( )
( ) ∑ ∑
− 1 − Xi,j,w bigM ≤ Uj,LL,s,w + 2LDi,j RLine P net a
i,j,k i,j,LL,k,s,w + X Line net a
i,j,k Qi,j,LL,k,s,w − Ui,LL,s,w
k∈ΩK k∈ΩK
( ) (13)
( )2 ∑ ( )2 ∑( )2 ( )
D
+ Li,j RLine
i,j,k J a
i,j,LL,k,s,w + X Line
i,j,k J a
i,j,LL,k,s,w ≤ 1 − Xi,j,w bigM
k∈ΩK k∈ΩK

6
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

and also, the routing of feeders is accomplished in an optimal manner

such that the system has the best performance both under normal and ωNew Ex pos neg
i,j,k − ωi,j,k = Ai,j,k − Ai,j,k (26)
emergency conditions. Eq. (19) shows the current flowing through line i- In Eq. (27), BNew
i,j,k is binary variable showing that a new line with type
j with regards to its conductor type. Relation (20) ensures that when
k must be installed between buses i and j through line i-j. In (27), product
ωai,j,k,w = 0, the current flow must be zero and when ωai,j,k,w = 1, the
of two binary variables makes the model non-convex. By using the
current flow must be lower than the thermal capacity of the corre­ relation (28), Eq. (27) will be relaxed and convexified. In order to pre­
sponding feeder of line (route) i-j connected between nodes i and j. The vent replacing an existing conductor with a weaker one (with lower
voltage magnitude limitation is expressed by (21) and (22). Except for capacity), relation (29) has been employed. To select just one type
(14), all the power flow constraints in the LFB model are convex and among various conductors for the conductor of each line, relation (30)
linear. To convexify (14), the conic relaxation method can be used. must be satisfied denoting that the sum of line types must be equal to
Hence, it is rewritten as (23). Xi,j,w . The type of conductor is expressed using (31).
( )2 ( )2
Pnet
i,j,LL,s,w + Qnet
i,j,LL,s,w ≤ Ji,j,LL,s,w Uj,LL,s,w (23) Zi,j,k ωNew New
i,j,k = Bi,j,k (27)

Now, the above problem is a MIQCP optimization model that can be 2BNew New New
(28)
i,j,k ≤ Zi,j,k + ωi,j,k ≤ Bi,j,k + 1
solved with GUROBI, MOSEK, CPLEX, and other appropriate solvers
(Behzadi, Bagheri, and Rabiee, 2023). ∑ ∑
Zi,j,k ωEx Sell
i,j,k Lk ≤ BNew Buy
i,j,k Lk (29)
k∈ΩK k∈ΩK
2.2.2. Conductor type selection for existing and new lines
∑
One of the effective elements in optimal distribution system planning Xi,j,w = ωai,j,k,w (30)
is choosing optimal type of conductor for the feeders. The selection of k∈ΩK

conductors’ type should be done in such a way that, in addition to ∑

having appropriate capacity, they are economical in terms of cost. For ωai,j,k,w
this aim, the existing feeders’ conductors should be optimally re- (31)
w∈ΩW
∑ ≤ ωNew
i,j,k
w
selected, and optimal conductor type for the new routes should be w∈ΩW
determined. In this paper, as one of contributions, a new convex
formulation is presented to satisfy the mentioned cases as follows. Eq. 2.2.3. Micro-grid formation
(24) indicates changing the type of lines in expansion models using a Nowadays, scientists in the electric power industry have come to the
binary variable Zi,j,k . If Zi,j,k equals to 1, it shows that a new feeder with conclusion that it is very difficult to maintain the security and stability of
type k is needed to be constructed between buses i and j (line i-j) or the a large grid in the face of extreme events. For this reason, by dividing a
existing conductor needs to be changed to increase the line capacity. large network into several optimal micro-grids, in addition to increasing
Also, if Zi,j,k is equal to 0, it means that there is neither need to construct a the resilience and stability of the system, it becomes easier to manage
new feeder nor to change the existing feeder between buses i and j (line i- the crisis in supplying electricity to the consumers in the micro-grids. In
j). In (24), totally four states can be occurred which are evaluated in addition, these MGs should be connected to each other via the tie-lines
Table 2: to be resilient and sustainable in the emergency condition. To achieve
The variable Zi,j,k has been used for calculation of installation or this goal, a new convex formulation is proposed in this paper as one of
reinforcement costs of lines. The first state in Table 2 indicates that there the contributions. Based on this formation, each bus of the network must
is no existing line between buses i and j, and also, no line has been be in a specific island. Therefore, there is no isolated bus within the
installed after the expansion planning; hence, the related cost must be network. In this paper, the system has the ability to be divided into
zero. The second state expresses that there is a line with type k in path i-j separated islands in which only one substation exists. In the following,
before the expansion, but, after the expansion, there is no line between i- the relations required for the micro-grid formation are discussed. For
j having type k; therefore, this shows that this line has been reinforced this aim, the first issue that should be considered is that each bus and
(its conductor has been reinforced by replacing its type with a new type each line in the network must belong to only one island (micro-grid) in
having higher capacity); hence, Zi,j,k must be equal to 1 such that the each outage state. This issue is expressed by (32)-(35). Eq. (36) relates
replacement of this line is considered in cost calculation (sale of old the value of lines binary variables to that of MGs.
conductor). In the third state, a new type of conductor is installed in path ∑
i-j which it was not present before the expansion. This state can be for ai,m,w = 1 (32)
m∈ΩMG
the case that the conductor of path i-j has been reinforced or a new line
with type k has been installed. Therefore, installation cost of the line ⃒ ⃒
⃒ai,m,w − aj,m,w ⃒ ≤ 1 − X MG (33)
having type k must be calculated. Finally, state 4 indicates that there is
i,j,m,w

no change in the type of line i-j, and there is no need to reinforce the MG
Xi,j,m,w ≤ ai,m,w (34)
conductor, and hence, there is not any related cost. Now, according to
relations (27) and (4), in the sates that Zi,j,k = 1 and ωNew i,j,k = 1, cost of MG
Xi,j,m,w ≤ aj,m,w (35)
new installed line is calculated. When Zi,j,k = 1 and ωEx
i,j,k = 1, the income
∑
due to sale of old conductor is calculated. Xi,j,w = MG
Xi,j,m,w (36)
⃒ ⃒ m∈ΩMG
⃒ New ⃒
⃒ωi,j,k − ωEx
i,j,k ⃒ = Zi,j,k (24)
Due to the presence of absolute function in (33), this relation is non-
Presence of the absolute function in (24) makes the model non- linear and non-convex. To linearize (33), relations (37) and (38) are
pos neg
convex. In order to convexify Eq. (24), relations (25) and (26) have employed in this paper. It should be noticed that ri,j,m,w and ri,j,m,w cannot
neg
been developed. According to (25), Apos
i,j,k and Ai,j,k cannot be equal to 1 at be equal to 1 at the same time. In relation (39), the variable Xi,j R
is
the same time. employed only to show the status of installed main/reserve feeders, and
it has not been used either in the constraints or in the objective function.
Zi,j,k = Apos neg
i,j,k + Ai,j,k (25)
In each state of w, some of the system feeders are close, and the others
must be open to satisfy the radiality constraint. The variable Xi,j,w just

7
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

shows the close status of lines by a binary variable of 1, and the open maximum installable number. Based on (56)-(57), the number of
status by 0. Therefore, the user (network operator/designer) cannot transformers must be zero when there is no substation in bus i. The
easily identify the newly installed lines unless he/she examines different number of transformers must be at least 1 if there is either new con­
states of contingency. However, by defining a new variable as Xi,j
R
, all the structed substation or expanded one.
newly installed lines (either main lines or reserve lines) can be easily Total
Si,LL,m,s,w α
= Si,LL,m,s,w β
+ Si,LL,m,s,w γ
+ Si,LL,m,s,w (46)
provided for the network operator.
∑( )
pos
ai,m,w − aj,m,w = ri,j,m,w neg
− ri,j,m,w (37) PSUB Pαi,LL,m,s,w + Pβi,LL,m,s,w + Pγi,LL,m,s,w (47)
i,LL,s,w =
m∈ΩMG
pos
ri,j,m,w neg
+ ri,j,m,w MG
≤ 1 − Xi,j,m,w (38)
∑( )
∑ QSUB
i,LL,s,w = Qαi,LL,m,s,w + Qβi,LL,m,s,w + Qγi,LL,m,s,w (48)
Xi,j,w ≤ Xi,jR ≤ ωNew
i,j,k (39) m∈ΩMG

k∈ΩK
( )2 ( )2 ( γ )2 ( )2
Pγi,LL,m,s,w + Qγi,LL,m,s,w ≤ Si,LL,m,s,w ≤ SUBγi,m SMax (49)
2.2.4. Radiality constraint
( )2 ( )2 ( )2
All micro-grids within the distribution network are operated in radial ( )2
Pαi,LL,m,s,w + Qαi,LL,m,s,w ≤ Si,LL,m,s,w
α
≤ SUBαi,m SMax (50)
structure in all conditions (either normal or emergency). To obtain
radiality, the following two conditions must be satisfied: ( )2 ( )2 ( β )2 ( )2
Pβi,LL,m,s,w + Qβi,LL,m,s,w ≤ Si,LL,m,s,w ≤ SUBβi,m SMax (51)
(a) The number of network lines must be equal to the number of
nodes minus one;
γ
Si,LL,m,s,w ≤ Transγ Transcapacity (52)
(b) There must be a path from the slack bus in each MG to all buses
located in the same MG.
α
Si,LL,m,s,w ≤ Transαi Transcapacity (53)

The reason why condition (a) is required is that if the condition (b) is β
Si,LL,m,s,w ≤ Transβi Transcapacity (54)
met, no loop can be formed in the MGs (it is met in Eq. (40)). The reason
for the requirement of condition (b) is that no bus nor a set of buses is Transγ + Transαi ≤ TransMax (55)
isolated in the MGs; if this happens (isolation of buses), a loop will be
∑ ∑
created in the system, which contradicts the radiality condition (Lav­ SUBαi,m ≤ Transαi ≤ TransMax SUBαi,m (56)
orato, Franco, Rider, and Romero, 2011). In this paper, the authors have m∈ΩMG m∈ΩMG

improved the method proposed in (Lavorato, Franco, Rider, and ∑ ∑

Romero, 2011) to establish a radial structure within the micro-grids in SUBβi,m ≤ Transβi ≤ TransMax SUBβi,m (57)
normal as well as emergency condition in which the load shedding oc­ m∈ΩMG m∈ΩMG

curs. Relation (40) guarantees connection of MG buses to each other. By Relation (58) shows the number of formed MGs, in which different
the relations (41)-(45), using the idea of a fictitious power flow, the numbers of MGs are compared to select the optimal value. In each MG,
desired goal in condition (b) can be achieved. there must be one substation as (59). Relation (60) implies that each
∑∑ ∑ substation in the network must belong to only one micro-grid. Relation
MG
Xi,j,m,w = ai,m,w − 1 (40)
i∈ΩN j∈ΩN i∈ΩN (61) allocates the bus i to micro-grid m in which the substation has been
installed. According to (62), the number of newly installed substations
∑
Geni,LL,m,w + PLsh,island L
(41) must be lower than its maximum value. Also, (63) states that the
i,LL,m,w − λLL ai,m,w Pi = flowi,j,LL,m,w
j∈ΩN expansion is allowed for the existing substations. Finally, (64) gives the
total number of substations as the sum of existing and new ones.
flowi,j,LL,m,w + flowj,i,LL,m,w = 0 (42) ∑ ∑ ∑
MGNum = SUBTotal
i,m ≤ m (58)
(43)
MG i∈ΩN m∈ΩMG m∈ΩMG
Xi,j,m,w flowMin ≤ flowi,j,LL,m,w ≤ Xi,j,m,w
MG
flowMax
∑
⎧ SUBTotal
i,m = 1 (59)
⎨ PLsh,island = 0; ∀w =w1 i∈ΩN
(44)
i,LL,m,w
⎩ 0 ≤ PLsh,island L
i,LL,m,w ≤ 0.99ai,m,w λLL Pi ; ∀otherwise ∑
SUBTotal
i,m ≤ 1 (60)
m∈ΩMG
0 ≤ Geni,LL,m,w ≤ SUBTotal
i,m Gen
Max
(45)
SUBTotal
i,m ≤ ai,m,w (61)
2.2.5. Substation constraints
∑ ∑
Always due to the high cost required for development and con­ SUBβi,m ≤ SUBβ,Max (62)
struction of substations, preparing optimal capacity of existing sub­ i∈ΩN m∈ΩMG

stations as well as siting and sizing of new substations is one of the main
concerns. As another novelty of this paper, a new convex formulation is SUBαi,m ≤ SUBγi,m (63)
presented to solve the problem. According to Eqs. (46)-(48), the total
apparent power, and active and reactive powers injected to the network SUBTotal γ β
i,m = SUBi,m + SUBi,m (64)
through all substations are given, respectively. By changing the equality
sign to inequality, on the left side of relations (49)-(51), the relations 2.2.6. Outage events
become relaxed, and the right side makes a limit on power transmission In order to make sure of distribution network resiliency and sus­
through substations after expansion. Also in (52)-(54), maximum ca­ tainability against natural disasters, the expansion planning must be
pacity of each substation is limited by determining the number of implemented considering severe contingency events. In this paper,
transformers after expansion. According to relation (55), the total different outage states are considered to make sure of network resil
number of transformers in existing substations must be lower than the

8
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

γ
Si,LL,m,s,w ≤ γwi,m,w SMax (68)

α
Si,LL,m,s,w ≤ γwi,m,w SMax (69)

β
Si,LL,m,s,w ≤ βwi,m,w SMax (70)

γ wi,m,w ≤ SUBγi,m (71)

βwi,m,w ≤ SUBβi,m (72)

( )(( ) )
∑ ∑ ∑ ∑
γwi,m,w = SUBγi,m w − 1 (73)
m∈ΩMG w∈ΩW m∈ΩMG w∈ΩW

( )(( ) )
∑ ∑ ∑ ∑
βwi,m,w = SUBβi,m w − 1 (74)
m∈ΩMG w∈ΩW m∈ΩMG w∈ΩW

( )
∑ ∑ ∑ ∑
γwi,m,w + βwi,m,w = SUBγi,m + SUBβi,m − 1; ∀w ∕
= w1 (75)
i∈ΩN m∈ΩMG i∈ΩN m∈ΩMG

( ( ))
∑ ∑
1− βwi,m,w + γwi,m,w ai,m,w ≤ 1 (76)
i∈ΩN i∈ΩN

Fig. 2. Single-line diagram of the 24-bus test system. ∑

FWi,w = RWi + βwi,m,w + γ wi,m,w (77)
iency. Hence, as another novelty of this paper, the following new convex m∈ΩMG

formulations are presented to model the different outage states. In each

ZWi,j,w = FWi,w FWj,w (78)
outage event, one of the substations and all its outgoing lines are taken
out of services, and on this basis, optimal resilient routing and hardening
Ji,j,LL,s,w ≤ ZWi,j,w bigM (79)
of feeders, sizing and siting of substations and renewable and non-
renewable DGs, and optimal micro-grids formation are determined. 2ZWi,j,w ≤ FWi,w + FWj,w ≤ ZWi,j,w + 1 (80)
The formulation is in such a way that it guarantees that all the sub­
stations and their outgoing feeders are taken out of service once, while
2.2.7. Wind turbine relations
the other equipment are in service. It should be noted that the most
The limitations concerning wind turbine units installation and
severe outage event is simultaneous outage of HV/MV substation and its
operation have been given in (81)-(83). The related variables and pa­
outgoing feeders. The constraint (65) guarantees that if there exists a
rameters are defined in the nomenclature for the sake of simplicity.
substation at a bus, the voltage of this bus is equal to 1 per unit (p.u.),
otherwise, the voltage of bus must satisfy the defined limitation. In mode PWT Num
i,LL,s,w ≤ WTi WTiIN CRs Pcapacity
WT (81)
w1, all the network equipment are in service, and the system operates
∑
under its normal condition, then, the relations (66) and (67) warranty WTiIN ≤ WTIN
Max
(82)
that the HV/MV substations are present in the network. Relations (68)- i∈ΩN
(70) express that when a substation is taken out of service due to a
contingency, its output power becomes zero. γ wi,m,w and βwi,m,w in (71) and WTiIN ≤ WTiNum ≤ WTiIN WTNum
Max
(83)
(72) confines the contingency occurrence only to the buses where the
substations located. Relations (73)-(75) express that only one substation 2.2.8. Distributed generation units
must be out of service in each outage state; also, all the substations must The installation and operation of dispatchable DGs are constrained
be in service during normal operating state. Relation (76) states that by by (84)-(87).
the outage of a HV/MV substation, its outgoing lines must be discon­ ( ) ( ) DG
− tan cos− 1 (PF) PDG DG − 1
i,LL,s,w ≤ Qi,LL,s,w ≤ tan cos (PF) Pi,LL,s,w (84)
nected as well. In order to make zero the current of lines connected to
the substation taken out of service, some binary variables have been ( )2 ( )2 ( Max )2
utilized in the proposed formulation. The variable RWi in (77) would be PDG
i,LL,s,w + QDG
i,LL,s,w ≤ DGNum
i DGIN
i SDG (85)
0 only when βwi,m,w and γ wi,m,w are equal to 1. Therefore, in case a substation
∑
is out of service, FWi,w will be 0 in the corresponding bus. In this way, DGIN Max
i ≤ DGIN (86)
ZWi,j,w equals to 0 which prevents flowing of current through the sub­ i∈ΩN

station’s outgoing feeders (according to (78)-(79)). To linearize (78), the

relation (80) has been employed. DGIN Num
i ≤ DGi ≤ DGIN Max
i DGNum (87)
( ) ( )
1 − βwi,m,w − γ wi,m,w U Min ≤ Ui,LL,s,w − βwi,m,w − γwi,m,w ≤ 1 − βwi,m,w − γwi,m,w U Max 2.2.9. Load shedding equations
(65) The curtailed active and reactive loads in emergency states as well as
resiliency index are expressed as (88)-(95). The reactive power shedding
SUBβi,m = βwi,m,w ; ∀w = w1 (66) is fulfilled considering loads power factor as (88). Also, the maximum
load shedding is constrained as (89)-(90). The resiliency index is yield
SUBγi,m = γ wi,m,w ; ∀w = w1 (67) from (93)-(95) based on the outage duration, loading level, and uncer­
tainty scenarios.

9
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

Table 3 Table 7
Load data of 24-bus system. Value of parameters used in simulation.
Bus Load Bus Load Bus Load Bus Load Parameter Value
(kVA) (kVA) (kVA) (kVA) Sub
IC (k$) 5000
1 5962 6 1584 11 3080 16 1342 ICTrans (k$) 200
2 1331 7 4796 12 1419 17 2640 TF (days a year) 20
3 4378 8 1034 13 1485 18 2310 Infr (%) 10
4 539 9 1947 14 3476 19 1991 Intr (%) 12
5 517 10 2640 15 1782 20 4169 Horizon year 5
emc ($/ton) 45
GE (tons/MWhr) 0.632
DGE (tons/MWhr) 0.365
Table 4 DGMax 4
Num
Line length data of 24-bus system. WTMax 3
Num

From To Length From To Length From To Length PF 0.9

{
(km) (km) (km) Buses voltage limit 0.95 ≤ Vi,LL,s,w ≤ 1.05; ∀w =w1
0.9 ≤ Vi,LL,s,w ≤ 1.05; ∀w ∕
= w1
1 5 3.885 4 9 2.1 7 23 1.575 Resiliency in normal operation RIw = 100; ∀w = w1
1 9 2.1 4 15 2.8 8 22 3.5
1 14 2.1 4 16 2.275 10 16 1.4
1 21 3.85 5 6 4.2 10 23 2.275
2 3 3.5 5 24 1.225 11 23 2.8 Table 8
2 12 1.925 6 13 2.1 14 18 1.75 Data of load level and energy price.
2 21 2.975 6 17 3.85 15 17 2.1
3 10 1.925 6 22 4.55 15 19 2.8 Load Load value in p. Duration in year Energy price VOLLLL
3 16 2.1 7 8 3.5 17 22 2.625 level u. (λLL ) (TLL ) (PrSUB
LL )
3 23 2.1 7 11 1.925 18 24 2.625
1 1 1000 hrs 60 $/MWhr 120
4 7 4.55 7 19 2.8 20 24 1.575
$/MWhr
2 0.75 5260 hrs 40 $/MWhr 100
$/MWhr
3 0.4 2500 hrs 24 $/MWhr 80
Table 5
$/MWhr
Existing feeders’ conductor type.
From To Conductor type (ωEx
i,j,k )

1 21 Type 2 Table 9
2 3 Type 1 Uncertainty scenarios of wind power generation.
2 21 Type 2
5 6 Type 1 Scenario Percent of rated power (CRs ) Probability (Probs )
7 8 Type 1 1 0.85 0.4
6 22 Type 2 2 0.60 0.5
8 22 Type 2 3 0.25 0.1

Table 6 Table 10
Different conductor types used in feeders. Specification of HV/MV substations.
Types RLine
i,j,k (ohms/ XLine
i,j,k (ohms/ IkAmp
Buy
Lk (k LSell
k (k Substation Location Existing capacity Expandable/Installable capacity
km) km) (Ampere) $/km) $/km)
1 Bus 21 1 × 7.5 MVA 3 × 7.5 MVA
Type 0.3080 0.1596 114 30 15 2 Bus 22 1 × 7.5 MVA 3 × 7.5 MVA
1 3 Unknown 0 3 × 7.5 MVA
Type 0.1208 0.1442 208 54 27 4 Unknown 0 3 × 7.5 MVA
2
Type 0.0487 0.1217 400 125 62.5
3
Table 11
Specifications of candidate DG units.
QLi Lsh DG Capacity of each Installation cost (k Operation cost
QLsh
i,LL,s,w = P (88)
PLi i,LL,s,w technology unit (MVA) $/MVA) ($/MWhr)

Gas turbine 1 400 46

PLsh,MG L
i,LL,m,s,w ≤ ai,m,w λLL Pi (89) Wind turbine 1 800 0

QLsh,MG L
i,LL,m,s,w ≤ ai,m,w λLL Qi (90)
TF ∑ ∑
∑ PL,contingency = TLL λLL PLi (94)
PLsh
i,LL,s,w = PLsh,MG
i,LL,m,s,w (91) Total
365 i∈ΩN LL∈ΩT
m∈ΩMG

∑ PL,contingency − PLsh,Total
QLsh
i,LL,s,w = QLsh,MG (92) RIw = Total w
× 100 (95)
m∈ΩMG
i,LL,m,s,w
PL,contingency
Total

TF ∑ ∑ ∑ 2.3. Final model

PLsh,Total
w = Probs TLL PLsh
i,LL,s,w (93)
365 i∈ΩN LL∈ΩT s∈ΩS
The mathematical formulation including the objective function and
constraints was described as relations (1)-(95). All the equations are

10
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

Fig. 3. Network configuration after expansion in case 1 in state w1.

Fig. 4. Network configuration after expansion in case 2 in state w1.

Table 12 Table 13
Substation expansion results in the two cases. Results of installed DGs in case 2.
Substation @ bus Optimal capacity DG @ bus Number of installed units
Existing New
Gas turbine 7 1×1
Case 1 21 —- 3× 7.5 MVA 11 1×1
22 —- 2× 7.5 MVA 16 1×1
—- 4 3× 7.5 MVA 18 1×1
Case 2 21 —- 3× 7.5 MVA Wind turbine 10 3×1
22 —- 3× 7.5 MVA 20 3×1
24 3×1

present in the final optimization model except the relations (14), (24),
(27), (33), and (78). The mentioned relations have non-linear nature,
and they were convexified using proposed formulations.

11
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

Table 14
Comparison of cost components in two cases.
Case Case
1 2

ICSUBα (M$) 0.60 0.80

ICSUBβ (M$) 5.60 0
ICLine (M$) 3.56 2.63
ICDG,WT (M$) 0 8.80
Total installation cost (M$) 9.76 12.23
OCSUB (M$) 39.79 30.34
OCDG (M$) 0 3.23
Total operation cost (M$) 39.79 33.57
EMC (M$) 23.08 17.61
Energy not supplied cost at each outage state and its W1 0 0
average value (M$) W2 0.58 1.12
W3 0.54 0.73
W4 0.49 —-
ENSC 0.54 0.92
OF (M$) 73.17 64.33
RIw (%) W1 100 100
W2 87.72 75.93
W3 88.30 84.77
W4 89.73 —-

Fig. 6. Network configuration in outage state w3 in case 2.

21.11MVAr at the peak load hour. Table 3 reports the loads data, and
Table 4 gives the length of lines (Wu et al., 2018). The conductor type of
existing feeders is given in Table 5; the characteristics of each conductor
type are shown in Table 6. In addition, the required constant parameters
of the network equipment are given in Table 7 (Bagheri, Monsef, and
Lesani, 2015).
Table 8 shows the considered load levels, their duration in a year,
energy purchasing price, and the value of lost load (Bagheri, Monsef,
and Lesani, 2015). In addition, the uncertainty states of wind power
generation have been given in Table 9.
Table 10 presents the specifications of new and candidate HV/MV
substation. As shown in Fig. 2, there are two existing substations located
at buses 21 and 22 with the capacities given in Table 10. Also, two new
substations can be installed which their location and capacity are un­
known in advance and are to be determined by the optimization model.
Also, Table 11 presents the specifications of candidate DG units
(Bagheri, Monsef, and Lesani, 2015).

3.2. Simulation results

To evaluate the performance of the proposed approach, it has been

Fig. 5. Network configuration in outage state w2 in case 2. programmed and accomplished in GAMS (General Algebraic Modeling
System) as a powerful optimization package for different linear and non-
3. Numerical study linear problems having real, binary, and integer variables. This software
can handle optimization problems with large number of variables and
In this part, the proposed approach has been numerically tested in constraints. Regarding the convexity of the proposed model, which is a
order to evaluate its performance. MIQCP one, the “GUROBI” solver is employed to find the optimal so­
lution. Two simulation cases have been set up as follows.

3.1. Test system and the input parameters 3.2.1. Case 1

In this case, distribution system expansion planning is implemented
The proposed model in this paper has been applied to a 24-bus sys­ without the presence of DG units. That is, the network reinforcement is
tem as shown in Fig. 2 and used in Ref. (Wu et al., 2018). This system is a fulfilled by construction of new lines and substations or by hardening the
20kV one having 33 routes (lines) in total. Among these 33 lines, 7 lines existing ones. Therefore, the aim is to optimally determine the route and
are the existing ones, and the remaining 26 lines are candidate lines that conductor type of new lines, possible hardening of existing feeders,
are to be constructed in a way that each MG maintains its radial struc­ location and capacity of new substations, and new capacities of existing
ture, and the tie-lines connect the MGs to each other during emergency substations;
conditions if required. There are two existing HV/MV substations at
buses 21 and 22. The total load of this system is 43.58MW and

12
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

Fig. 7. Active power of WT units in case 2, (a): bus 10, (b): bus 20, (c): bus 24.

Fig. 8. Apparent power imported from upward grid into the distribution system in case 1.

13
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

Fig. 9. Apparent power imported from upward grid into the distribution system in case 2.

Fig. 10. Voltage of buses in load level 1 (LL1) and normal operating condition (w1) in case 2.

Fig. 11. Voltage of buses in load level 1 (LL1) and outage state 2 (w2) in case 2 (bus 21 is out of service).

Fig. 12. Voltage of buses in load level 1(LL1) and outage state 3 (w3) in case 2 (bus 22 is out of service).

3.2.2. Case 2 number of DG units must be determined. It should be noticed that in the
In the second case, as well as the options of the first case, the DG units 24-bus system, the considered disaster states include outage of each
from both renewable and non-renewable types can be also incorporated substation along with its outgoing feeders.
in the planning problem. Therefore, the optimum type, location, and Figs. 3 and 4 respectively illustrate the configured network after the

14
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

Fig. 13. Active load shedding in load level 1 (LL1) in case 2.

Fig. 14. Reactive load shedding in load level 1 (LL1) in case 2.

expansion planning for cases 1 and 2 in normal operating mode (w1). As the resiliency of the system in emergency conditions when some parts of
seen, the routing of existing and new lines is in a way that the obtained the system are affected by natural disasters.
networks have radial structure. The installed reserve feeders are shown Total operation cost in case 2 is 33.57M$ while it is 39.79M$ in case
by dashed lines having the cross mark (×). These lines are open in 1. This means that the operation cost in the presence of DG units has a
normal state and will be used after the fault occurrence for the network reduction of about 15.7% compared to the case when DG units are not
reconfiguration to manage the emergency conditions. The Figs. 3 and 4 present. In addition, the emission cost in case 2 has reduced by 24%
show also the service area of each HV/MV substation, that is, which showing efficient scheduling of DG units. The resiliency index (RI) for
buses are supplied by the considered substation. The specifications of each outage state (w) has been given in last row of the Table 14. It is seen
the substations after the planning in the two cases are given in Table 12. that in all outage states, the restored load is at least about 76% which
As it can be seen, in case 1, the existing substations @ buses 21 and 22 means that only 24% of the total load is curtailed. It can be observed that
have been reinforced up to the capacity of 3 × 7.5MVA and 2 × 7.5MVA, the total installation cost (for substation, line, and DG) in case 2 is
respectively. Also, a new substation has been sited and installed @ bus 4 12.23M$ while it is 9.76M$ in case 1. Although the installation cost in
with the maximum installable capacity of 3 × 7.5MVA. In case 2, there is case 2 is 25.3% more than case 1, the total cost is lower by 12.1%. This
no newly installed substation and both of the existing substations have cost reduction is mostly due to decreasing operational and emission
been upgraded to the capacity of 3 × 7.5MVA. That is, the presence of costs, while the resiliency index is acceptable and the expansion plan­
DG units in case 2 has prevented installation of new high-cost substation. ning is cleaner from the CO2 emission viewpoint.
It can be said that the total installed capacities in cases 1 and 2 are In case 1, for state w1, as all the loads must be supplied, the resiliency
45MVA and 30MVA, respectively. The characteristics of installed DG index (RI) is 100%. In states w2, w3, and w4, the substations no. 21, 22,
units in case 2 according to Fig. 4 are given in Table 13. For comparing and 4 are out of service, respectively. The same issue is true for case 2 in
the conductors used in the feeders, in Figs. 3 and 4, the feeders have been which the substation no. 21 and 22 are out of service in states w2 and w3,
highlighted with colors representing their conductor type based on respectively. As in case 1, there is three installed substations, the resil­
Table 6. The feeders nearer to HV/MV substations usually have higher iency index is somehow better than case 2 in which the number of
capacity, and the other feeders are from low-capacity types. Installation installed substations is only 2. As the number of installed energy sources
of DGs in case 2 has caused that the employed lines have lower capac­ (DGs) is higher in case 2, the network relies on more power supply
ities and costs in total. Figs. 3 and 4 show that the network configuration points, and this would be helpful when other outage states are
in cases 1 and 2 are different from the viewpoint of main and reserve considered.
feeders (tie-lines) routing as well as their conductor type. The lines Figs. 5 and 6 depict the network configuration in two outage states in
expansion and installations cost in case 2 is lower than that in case 1. case 2. In Fig 5, the outage is defined on substation 21 and its outgoing
This is due to the presence of DG units and releasing the capacity of feeders. Also, in Fig. 6, the outage of substation case 2 is considered. As
feeders and eliminating the need for installing high-cost lines. Routing of seen, the network is configured in a way that the loads are supplied from
appropriate tie-lines with optimal type can be also accounted as the the other substations and DG units by means of installed tie-lines. It can
reason of this cost reduction. The tie-lines (reserve feeders) can improve be seen that the radiality constraint is also preserved under emergency

15
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

conditions. In Fig. 7, the generated power of wind units installed at the acquire optimal normal operation as well as to withstand against severe
buses has been shown for each load levels, outage states, and uncertainty outages during natural disasters. The objective function includes in­
scenarios. The outage states, shown by w1, w2, and w3 are respectively vestment, operational, emission, and load shedding costs. Also, the
for normal state, outage of the first HV/MV substation, and outage of the operational and other constraints are taken into consideration in order
second HV/MV substation. The generation level of these units is to compose a convex MIQCP optimization model. The simulations using
maximum in uncertainty scenario of s1 and is minimum in scenario s3. global optimum solvers in GAMS verified the effectiveness of the con­
On the other side, they are scheduled with maximum capacity during ducted model in optimal planning of substations, main and tie feeders,
peak load level (that is LL1); also, they generate the required power in and sizing and siting of DGs. The obtained results illustrated reduction of
the lower load levels. total cost, reduction of emission cost, and improvement of system
In Figs. 8 and 9, the imported power from the upward grid into the resiliency and sustainability. The proposed formulation is general and
distribution system though all HV/MV substations have been depicted not dedicated to the used case study, and it can be applied to any other
for different load levels and emergency states in cases 1 and 2. As seen in system having different buses and configuration.
Fig. 8, the imported power is high in peak load hours, and as the load
decreases (from LL1 to LL3), the imported power is reduced. On the other
side, in emergency conditions (w2, w3, and w4), regarding the outage of Declaration of Competing Interest
one of HV/MV substations, the imported power of existing substations is
more than that in normal operating condition (w1). All the described The authors declare that they have no known competing financial
issues are also seen in Fig. 9. Furthermore, in Fig. 9, the imported power interests or personal relationships that could have appeared to influence
is shown for different uncertainty scenarios. In scenario s1, the wind the work reported in this paper.
units have their maximum power, and they generate the minimum
power in scenario s3. Therefore, in Fig. 9, the imported power is Data availability
increased form s1 to s3.
Figs. 10-12 illustrate voltages of network buses for different load No data was used for the research described in the article.
levels, outage states, and uncertainty scenarios. In Fig. 10, voltages of
buses 21 and 22 at which the HV/MV substations have been installed are
1 per unit (pu), and the voltage of other buses are within the permissible References
range of [0.95-1.05] pu. In Fig. 11, due to the outage of substation 21 in
state w2, the voltage of bus 21 is not present, and the voltage of other Abessi, Ahad, Jadid, Shahram, & Salama, M. M. A. (2020). A new model for a resilient
HV/MV substation (bus 22) is 1 pu. The similar issue is true in Fig. 12 distribution system after natural disasters using microgrid formation and considering
ICE cars. IEEE Access, 9, 4616–4629.
about voltages of buses 21 and 22, as the HV/MV substations, when the Bagheri, A., Monsef, H., & Lesani, H. (2015). Integrated distribution network expansion
substation 22 is out of service. The permitted voltage range in states w2 planning incorporating distributed generation considering uncertainties, reliability,
and w3 is [0.9-1.05]; this is while all the voltages are more than 0.96pu, and operational conditions. International Journal of Electrical Power & Energy Systems,
73, 56–70.
which shows satisfaction of voltage constraints under the emergency Bagheri, A., Monsef, H., & Lesani, H. (2015). Renewable power generation employed in
conditions. In Figs. 13 and 14, the curtailed active and reactive powers an integrated dynamic distribution network expansion planning. Electric Power
of the network in states w3 and w2 have been shown for different buses Systems Research, 127, 280–296.
Bayani, Reza, & Manshadi, Saeed D. (2023). Resilient expansion planning of electricity
and uncertainty scenarios. The maximum load shedding corresponds to grid under prolonged wildfire risk. IEEE Transactions on Smart Grid.
bus 1 in state w2. Bayat, A., & Bagheri, A. (2019). Optimal active and reactive power allocation in
distribution networks using a novel heuristic approach. Applied Energy, 233, 71–85.
Behzadi, Saeed, Bagheri, Amir, & Rabiee, Abbas (2023). Optimal operation of
4. Discussion reconfigurable active distribution networks aiming at resiliency Improvement. In
2023 8th International Conference on Technology and Energy Management (ICTEM) (pp.
The proposed model comprehensively considered all the planning 1–7). IEEE.
Borghei, Moein, & Ghassemi, Mona (2021). Optimal planning of microgrids for resilient
options of resilient distribution systems including hardening of existing
distribution networks. International Journal of Electrical Power & Energy Systems, 128,
lines, routing and sizing of new main/reserve feeders, MG-based siting Article 106682.
and sizing of HV/MV substations, siting and sizing of dispatchable and Dada, Joseph O. (2014). Towards understanding the benefits and challenges of Smart/
renewable DGs. Also, both normal and critical emergency conditions are Micro-Grid for electricity supply system in Nigeria. Renewable and Sustainable Energy
Reviews, 38, 1003–1014.
regarded for the network operation. Unlike the MINLP models and Ebadat-Parast, Mahsa, Nazari, Mohammad Hassan, & Hosseinian, Seyed Hossein (2022).
artificial intelligence-based approaches like genetic algorithm (GA), Distribution system resilience enhancement through resilience-oriented optimal
particle swarm optimization, etc., which do not provide optimal solu­ scheduling of multi-microgrids considering normal and emergency conditions
interlink utilizing multi-objective programming. Sustainable Cities and Society, 76,
tions, the optimality of results in this paper is guaranteed due to con­ Article 103467.
vexification of non-linear relations and making an MIQCP model which Elkadeem, Mohamed R., Elaziz, Mohamed Abd, Ullah, Zia, Wang, Shaorong, &
is solved by “GUROBI” as a powerful solver capable of finding global Sharshir, Swellam W. (2019). Optimal planning of renewable energy-integrated
distribution system considering uncertainties. IEEE Access, 7, 164887–164907.
optimum solutions for convex problems ([Link], Soroudi, 2017). The Galvani, Sadjad, Bagheri, Amir, Farhadi-Kangarlu, Mohammad, & Nikdel, Nazila (2022).
results verified that the presence of DG units cause the need for more A multi-objective probabilistic approach for smart voltage control in wind-energy
HV/MV substations and high-cost feeders. Also, the operation and integrated networks considering correlated parameters. Sustainable Cities and Society,
78, Article 103651.
emission costs are reduced, and as a result, the total cost is decreased Gilani, Mohammad Amin, Dashti, Reza, Ghasemi, Mostafa, Amirioun, Mohammad
with a percent of about 12% by maintaining the resiliency index in an Hassan, & Shafie-khah, Miadreza (2022). A microgrid formation-based restoration
acceptable range. model for resilient distribution systems using distributed energy resources and
demand response programs. Sustainable Cities and Society, 83, Article 103975.
Gritzalis, Dimitris, Theocharidou, Marianthi, & Stergiopoulos, George (2019). Critical
5. Conclusion infrastructure security and resilience. Springer International Publishing, 10, 978.
Haghshenas, Morteza, Hooshmand, Rahmat-Allah, & Gholipour, Mehdi (2022). Power
distribution system resilience enhancement planning against extreme dust storms via
A new convex model based on micro-grid formation is proposed in
pre-and post-event actions considering uncertainties. Sustainable Cities and Society,
this paper for resilience-oriented expansion planning of distribution 78, Article 103626.
systems. Different planning alternatives such as HV/MV substations, Lavorato, Marina, Franco, John F., Rider, Marcos J., & Romero, Rubén (2011). Imposing
main and reserve feeders, and renewable and non-renewable DG units radiality constraints in distribution system optimization problems. IEEE Transactions
on Power Systems, 27(1), 172–180.
have been considered for the conducted planning. The proposed model Lopes, J. A. Pecas, Hatziargyriou, Nikos, Mutale, Joseph, Djapic, Predrag, & Jenkins, N.
is able to reconfigure the network into appropriate micro-grids to (2007). Integrating distributed generation into electric power systems: A review of

16
S. Behzadi and A. Bagheri Sustainable Cities and Society 97 (2023) 104787

drivers, challenges and opportunities. Electric Power Systems Research, 77(9), Shahbazi, Amid, Aghaei, Jamshid, Pirouzi, Sasan, Shafie-khah, Miadreza, &
1189–1203. Catalão, João P. S. (2021). Hybrid stochastic/robust optimization model for resilient
Prathapaneni, Dimple Raja, & Detroja, Ketan P. (2019). An integrated framework for architecture of distribution networks against extreme weather conditions.
optimal planning and operation schedule of microgrid under uncertainty. Sustainable International Journal of Electrical Power & Energy Systems, 126, Article 106576.
Energy, Grids and Networks, 19, Article 100232. Soroudi, Alireza. (2017). Power system optimization modeling in GAMS, 78. Cham,
Rajabzadeh, Mohammad, & Kalantar, Mohsen (2022). Enhance the resilience of Switzerland: Springer.
distribution system against direct and indirect effects of extreme winds using battery Sun, Qirun, Wu, Zhi, Gu, Wei, Zhu, Tao, Zhong, Lei, & Gao, Ting (2021). Flexible
energy storage systems. Sustainable Cities and Society, 76, Article 103486. expansion planning of distribution system integrating multiple renewable energy
Ravadanegh, Sajad Najafi, Jamali, Saeed, & Vaniar, Arash Mohammadi (2022). Multi- sources: An approximate dynamic programming approach. Energy, 226, Article
infrastructure energy systems resiliency assessment in the presence of multi-hazards 120367.
disasters. Sustainable Cities and Society, 79, Article 103687. Wu, Zhi, Liu, Pengxiang, Gu, Wei, Huang, He, & Han, Jun (2018). A bi-level planning
Saberi, Reza, Falaghi, Hamid, Esmaeeli, Mostafa, Ramezani, Maryam, Ashoornezhad, Ali, approach for hybrid AC-DC distribution system considering N-1 security criterion.
& Izadpanah, Reza (2023). Power distribution network expansion planning to Applied energy, 230, 417–428.
improve resilience. IET Generation, Transmission & Distribution. Zakernezhad, Hamid, Nazar, Mehrdad Setayesh, Shafie-khah, Miadreza, & Catalão, João
Sallam, Abdelhay A., & Malik, Om P. (2018). Electric distribution systems. NJ: Wiley. P. S. (2021). Multi-level optimization framework for resilient distribution system
Saravi, Vahid Sabzpoosh, Kalantar, Mohsen, & Anvari-Moghaddam, Amjad (2022). expansion planning with distributed energy resources. Energy, 214, Article 118807.
Resilience-constrained expansion planning of integrated power–gas–heat Zhou, Siyu, Han, Yang, Yang, Ping, Mahmoud, Karar, Lehtonen, Matti,
distribution networks. Applied Energy, 323, Article 119315. Darwish, Mohamed M. F., & Zalhaf, Amr S. (2022). An optimal network constraint-
Seifi, Hossein, & Sadegh Sepasian, Mohammad (2011). Electric power system planning: based joint expansion planning model for modern distribution networks with multi-
Issues, algorithms and solutions. Vol. 49. Berlin: Springer. types intermittent RERs. Renewable Energy.
Seta, Felipe da Silva, Oliveira, Leonardo [Link], & Oliveira, Edimar [Link] (2019).
Comprehensive approach for distribution system planning with uncertainties. IET
Generation, Transmission & Distribution, 13(24), 5467–5477.

17