1 s2.0 S0142061521009042 Main

Electrical Power and Energy Systems 137 (2022) 107674
Contents lists available at ScienceDirect
International Journal of Electrical Power and Energy Systems

journal homepage: www.elsevier.com/locate/ijepes
An MILP formulation for the optimum operation of AC microgrids with

hierarchical control
Youthanalack Vilaisarn a, *, Majid Moradzadeh b, Morad Abdelaziz b, Jérôme Cros a
a
Department of Electrical and Computer Engineering, Université Laval, Québec, Canada
b
School of Engineering, University of British Columbia, Kelowna, Canada
A R T I C L E I N F O A B S T R A C T
Keywords: An accurate linear model for the optimal operation of the microgrid is a key enabling solution for many real time
Distributed generation and stochastic applications. Nonetheless, many of the existing microgrid’s optimal power flow models fall short
Hierarchical control in considering the details of the microgrid hierarchical control, which can lead to infeasible and suboptimal
AC microgrid
solutions. In this work, a new mixed integer linear programing (MILP) model for the optimal operation of three
Linear optimum power flow
phase microgrids with hierarchical control is proposed. First, the nonlinear formulation is derived to account for
the microgrid philosophy of operation in both grid connected and islanding modes of operation. Subsequently, a
linearization approach is adopted to the derived nonlinear model to transform it into the required MILP
formulation. The effectiveness of the proposed model is validated through several case studies considering
balanced and unbalanced test systems. Additionally, a detailed error assessment is conducted to justify the
proposed model.
upstream.
In order to optimize and control the operation of the microgrid, the
1. Introduction hierarchical control structure has been widely adopted in the literature
[4]. This structure employs several control layers, wherein each of these
Power distribution systems are currently undergoing a major tran layers has its own functionality and response time. A higher-level layer
sition towards higher penetration of renewable energy recourses along is typically employed for economic purposes and manages the DG units
with an increased application of energy storage technologies. As such, and ESSs commitment and dispatching decisions as well as the load
these systems are transforming from their traditionally passive structure curtailment decisions. Lower-level and intermediary-level layers are on
into being active distribution systems (ADSs) [1,2]. This transition to the other hand deployed for controlling the DG units power production
wards the ADS structure enhances the distribution system’s operators and for maintaining the microgrid’s voltage and frequency within
(DSOs) ability to operate and manage the distribution of electricity by bounds [5]. Still, the operation of the microgrid and its hierarchical
offering new possibilities and flexibilities. In this context, microgrids control structure are conceptually different in both grid-connected and
(MGs) have been put forward as attractive new concept. Microgrids are islanding modes of operation. In the grid-connected mode, the voltage
electric regions within the ADS comprising distributed generation (DG) and the frequency of the microgrid are imposed by the upstream grid at
units and possibly energy storage systems (ESSs). With their superior the point of common coupling (PCC) and any mismatch between the
ability of operating independently (i.e., islanding mode) and as a part of microgrid’s local generation and demand is supplied by the upstream
their upstream grids (i.e., grid connected mode), microgrids have grid. As such, in the grid-connected mode of operation, the DG units are
become a key enabling solution for future smart ADSs [3,4]. Specifically, controlled to inject constant amounts of power prespecified by the high-
a microgrid can be disconnected from the main grid and operated in the level layer of the hierarchical control structure [6]. On the other hand, in
islanding mode due to intentional (scheduled) or unintentional events the islanding mode of operation, the DG units cannot be controlled to
(unscheduled). In this islanding mode of operation, the microgrid’s local supply some prespecified constant amount of power and need to follow
generation and energy storage resources are controlled to supply its the load variations in order to ensure a stable system operation [78]. To
demand. On the other hand, a microgrid can also be operated in the grid this end, in the islanding mode, primary droop control is deployed as the
connected mode to provide ancillary services to the main grid in its
* Corresponding author.
E-mail address: youthanalack.vilaisarn.1@ulaval.ca (Y. Vilaisarn).
https://doi.org/10.1016/j.ijepes.2021.107674
Received 6 February 2021; Received in revised form 18 September 2021; Accepted 26 September 2021
Available online 30 November 2021
0142-0615/© 2021 Elsevier Ltd. All rights reserved.
Y. Vilaisarn et al. International Journal of Electrical Power and Energy Systems 137 (2022) 107674
Nomenclature time i
PLi,Ψ,t Active power demand phase Ψ located at bus i, time t
Abbreviation
PLi,Ψpeak Annual peak active demand of bus i
ADS Active Distribution System
QWT 3Φ
Three phase reactive power generated by ith wind units at
DDG Dispatchable Distributed Generation i,t
DER Distributed Energy Resources time i

DG Distributed Generation QLi,Ψ,t Reactive power demand phase Ψ located at bus i, time t
ESS Energy Storage System QLi,Ψpeak Annual peak reactive demand of bus i
GC Grid connected mode RC Certain radiation point, typically set to 0.15 kW/m2
IS Islanding mode
RSTD Solar radiation in standard test conditions, typically
ISW Isolated Switch
specified as 1000 W/m2
MG Microgrid 0
SOCi Initial state of charge of the ESS unit i
MILP Mixed Integer Linear Problem
Ticurt Maximum hours of load curtailment of bus i
MINLP Mixed Integer Nonlinear Problem
V Minimum voltage violation of the bus i
OPF Optimum Power flow
PCC Point of Common Coupling V Maximum voltage violation of the bus i
PV Photovoltaic Energy Vi* Voltage reference for primary control of the ith DDG or the
SOC State of Charge of the ESS ith ESS with primary droop control mode with a supervisory
SRF-PLL Synchronous Reference Frame Phase Locked-Loop secondary control layer
TR Transformer Vi** Voltage reference for secondary control of the ith DDG or
WT Wind Turbine Energy the ith ESS with primary droop control mode with a
supervisory secondary control layer
Sets and indices GC
Vspec Nominal specified voltage during grid-connected mode
ΩT Set of segments or hours in the planning horizon e.g. ΩT = m
bt Binary input parameter indicates the operation mode of the
{1, 2, 3, ⋯, 24}
microgrid at time t
Ωx Set of buses type × , where × can be {BUS, LIN, PCC, DDG,
bESm Binary input parameter indicates the operation mode of the
WDG, PV or ESS} i
R Continuous variable or real number ith ESS

Z Integer variable or binary number ki ω
Integral gain for frequency restoration for DDG or for ESS
with primary droop control mode with a supervisory
Parameter secondary control layer
E ES ES
Minimum/maximum allowable stored energy of the ESS mpi Active droop gain of the ith DDG unit or the ith ESS with
i /Ei
unit i primary droop control mode with a supervisory secondary
br control layer
Ik Maximum ampacity of the branch k
DDG DDG nqi Reactive droop gain of the ith DDG unit or the ith ESS with
P /Pi Minimum/maximum allowable for active power
i primary droop control mode with a supervisory secondary
generated by ith DDG unit
control layer
DDG DDG
Q i /Qi Minimum/maximum allowable for reactive power rtm Solar irradiance in kW/m2
generated by ith DDG unit vci , vco , vra Cut-in, cut-out, and the rated speed of the wind DG,
ES
Q ES
i /Qi Minimum/maximum allowable for reactive power of ith respectively
ESS unit vmt Wind speed in m/s
DDG
Si Capacity rating of the DDG unit i ω* Nominal angular frequency of the system
PV
Si Capacity rate of the PV DG located at bus i ηES
i Efficiency of the energy storage
TR σ curt Cost of load curtailment [$/kWh]
Si Capacity rate of the transformer located at the substation
WT
σ fuel Natural gas price [$/Thousands of cubic feet]
Si Capacity rate of the wind DG located at bus i σ PCC Market energy price at time t [$/kWh]
t
APLt Annual peak load factor according to seasonal pattern Δt Step time (i.e., 1 segment = 1 h)
BΨβ
ij Element of the susceptance matrix denote the connection Nicurt Limit on the number of load shedding switchings for each
of bus i phase Ψ to bus j phase β load at bus i
CRDDG
i Consumption rate of the natural gas of DDGi [scf/kWh]
Variables for OPF
GΨβ Element of the conductance matrix denote the connection
ij EES
i,t Energy SOC of the ESS unit i at time t
of bus i phase Ψ to bus j phase β
PDDG
i,t
3Φ
Three phase active power generated by the ith DDG unit at
PCVI Penalty coefficient of voltage deviation [$/kV]
time t
PFDDG 1st,4th
Minimum operating power factor allowable of the 1st/ PCC 3Φ
i Pi,t Three phase active power imported from the main grid
4th quadrant of the DDG unit i
across the ith PCC of the MG at time t
PFTR 1st,4th
Minimum operating power factor allowable of the 1st/
i Pch
i,t
3Φ
Active power charge by ith ESS unit at time t
4th quadrant of the transformer i
Pdch 3Φ
Active power discharge by ith ESS unit at time t
PFPV i Power factor set of the PV unit i i,t
Pdch Active power phase Ψ discharged by ith ESS unit at time t
PFWT i Power factor set of the wind unit i i,Ψ,t
PLi,t 3Φ Three phase active power demand located at bus i, time t Pi,Ψ,t Active power injected phase Ψ to bus i at time t
th PDDG Active power phase Ψ generated by the ith DDG unit at time
PPV
i,t
3Φ
Three phase active power generated by i PV units at time i i,Ψ,t
t
WT 3Φ
Pi,t Three phase active power generated by ith wind units at
PPCC
i,Ψ,t Active power phase Ψ imported from the main grid across
2
the ith PCC of the MG at time t Binary variables for OPF

Pbr
k,Ψ,t Active power phase Ψ flowing across the transmission bES
i,t Binary variable indicates the status of ESS i.e. 1: discharge
branch k at time t and 0: charge
QES
i,t
3Φ
Three reactive power generate/receive of the ith ESS unit, bcurt
i,t Binary variable representing load curtailment (1: load
at time t shaded and 0: load fed)
QES Reactive power phase Ψ generated by ith ESS unit at time t zcurt
i,t Binary variable representing a linearization of absolute
i,Ψ,t ⃒ ⃒
Qi,Ψ,t Re-active power injected phase Ψ to bus i at time t ⃒
term ⃒bcurt curt ⃒
i,t − bi,t− 1 ⃒
QDDG
i,Ψ,t Reactive power phase Ψ generated by the ith DDG unit at
Function
time t
CCURT Cost of energy curtailment during islanding mode at time t
QPCC
i,Ψ,t Three phase reactive power imported from the main grid t
CDDG Cost of natural gas consumption by the DDG units at time t
across the ith PCC of the MG at time t t
CPCC Cost of energy imported from the substation across the PCC
Qbr
k,Ψ,t Reactive power phase Ψ flowing across the transmission t
of the microgrid at time t
branch k at time t
CVD
t Cost of voltage deviation at time t
Vi,Ψ,t Voltage magnitude phase Ψ, of the bus i at time t
CVI
t
Lin
Linearized version of CVD
t
uv,int
i,Ψ,t Integrator output of secondary voltage restoration 3Φ 3Φ
XMINLP ,XMILP Vector containing the decision variables for the MINLP
considering phase Ψ of ith DDG unit or the ith ESS operating
and MILP problem formulation
with secondary control at time t
δi,Ψ,t Voltage magnitude angle of phase Ψ, bus i at time t
δMG
Ψ,t Phase angle of the bus where the frequency control act
Fig. 1. Basic scheme of the microgrid.
lower-level control layer in order to control the DG units to follow the work, the primary droop control has been considered and the microgrid
load variations. Secondary control is deployed as the intermediary-level operation was optimized for minimizing the distribution losses. Addi
control layer responsible for compensating the steady-state errors in the tionally, an MINLP model for the optimal operation of unbalanced
system voltage and frequency derived by the primary control. The sec islanded droop based microgrid has also been proposed in [6]. This
ondary control hence ensures that the microgrid’s voltage and frequency formulation was also linearized and transferred into a MILP model that
are maintained within bounds and that the microgrid is ready for re- can be solved with linear commercial solvers. While the work in [11,6]
synchronization with the main grid when moving back from the has overcome earlier drawbacks pertaining to the consideration of the
islanding mode to the grid-connected mode. Finally, tertiary control is primary droop control layer, still these models were restricted to the
adopted as the higher-level control layer used for ensuring the economic islanded mode of operation and did not consider the secondary control
viability of the microgrid operation by optimizing the microgrid’s layer. Recently, a generalized model has been introduced in [12] to
commitment and dispatch decisions. consider both the grid-connected and islanded modes of operation of an
The optimization of the operation of the microgrid resources and unbalanced microgrid. However, the secondary control layer has not
loads, performed by the higher-level control layer of the microgrid’s been considered in this work. On the other hand, as the secondary
hierarchal control structure, is based on an AC optimum power flow control layer change the DG units operating points and their active and
(ACOPF) problem. In this context, the considered ACOPF problem seeks reactive power injections, overlooking the secondary control in the
to optimize the microgrid’s operation subject to the system’s operational microgrid ACOPF problem formulation can lead to infeasible and sub
constraints. Several research works have recently attempted to accu optimal operating points [13].
rately formulate and solve the ACOPF problem for microgrid systems. In this sense, seeking to overcome the mentioned limitations of the
For instance, a bi-level model has been proposed in [9] to optimally existing literature, in this work we seek to develop an ACOPF formula
schedule the microgrid operation. However, this model was restricted to tion for AC microgrids with hierarchical control to simultaneously
the grid connected mode and the ability of the system to operate in consider: 1) both grid-connected and islanded modes of operation; 2)
islanded mode was not considered. The mathematical formulation pro balanced and unbalanced operation of microgrids; 3) the primary and
posed in [10] has considered the islanded mode of the microgrid. Yet, secondary control layers of the DG units in the islanded mode; 4) general
the primary and secondary control layers were not considered, which models of dispatchable resources, non-dispatchable resources and the
can lead to infeasible solutions. In [11], an ACOPF for inverter based energy storage; 5) the time coupling constraints of the resources and
microgrids operating in the islanded mode has been proposed. In this storage located in the microgrid. To this end, first we develop a new
3
Fig. 2. (a) Basic scheme of the microgrid (b) Steady-state model of DDG unit’s operation with primary and secondary control under IS mode.
MINLP formulation for the optimal operation of AC microgrids with ( )

∑
hierarchical control. Subsequently, the proposed MINLP formulation is CtPCC = Sbase σPCC
t PPCC
i,t
3Φ
Δt, ∀t ∈ ΩT (2)
linearized and a MILP formulation is extracted, which can be solved by i∈ΩPCC
commercial optimization solvers. The remainder of this paper is orga ( )

nized as follows: First, the proposed MINLP formulation to model the ∑
CtDDG = Sbase CRDDG PDDD 3Φ fuel
σ Δt, ∀t ∈ ΩT (3)
optimal operation of the unbalanced AC microgrid is presented in Sec i∈ΩDDG
i i,t
tion 2. Subsequently, proposed MINLP model is linearized and the

proposed MILP model is introduced in Section 3. The numerical simu (
∑
)
lation results are given in Section 4. Finally, Section 5 concludes the CtCURT = Sbase σcurt bcurt load L 3Φ
i,t wi Pi,t Δt, ∀t ∈ ΩT (4)
paper. i∈ΩBUS
( [ ])
∑ ∑ ⃒⃒( ) GC ⃒
2. MINLP problem formulation CtVD =Vbase PCVD ⃒ 1− bmt Vspec
⃒
+bmt V ** − Vi,Ψ,t ⃒ Δt, ∀t ∈ΩT
i∈ΩBUS Ψ∈F
In this section, the problem for the optimal operation of a three phase (5)
microgrid is formulated as an MINLP problem. The proposed formula
tion accounts for both islanded and grid connected modes of operation,
wherein the operation mode is defined by the status of the island 2.2. Constraints
isolating switch (ISW) located at the microgrid’s PCC as shown in Fig. 1.
btm denotes the status of ISW along the planning horizon Ω T; if btm = 0, The inequality (6) is adopted to have the voltage magnitude of
the microgrid operates in grid-connected mode, and if btm = 1, the different buses restricted to voltage magnitude constraints. Eq. (7) en
microgrid operates in islanding mode. This allows the problem to sures the quality of power supply by limiting the load curtailment
consider the objective function and constraints pertaining to the duration. Then, inequality (8) is used to limit the number of load
microgrid’s operation mode. In this work, the load and the generation shedding switching cycles for each load at bus i. Finally, the inequality
are modelled using the method described in Appendix A [14]. The (9) ensures bcurt
i,t activate only for islanded mode.
objective function and the related constraints for the MINLP model are V ≤ Vi,Ψ,t ≤ V, ∀i ∈ ΩBUS , ∀Ψ ∈ F , ∀t ∈ ΩT (6)
mathematically formulated in the following subsections.
∑
bcurt curt
i,t ≤ Ti , ∀i ∈ ΩBUS (7)
2.1. Objective function t∈ΩT
∑⃒⃒ ⃒
curt ⃒
The objective function given by (1) is comprised of four terms aiming ⃒bcurt curt
i,t − bi,t− 1 ⃒ ≤ 2Ni , ∀i ∈ ΩBUS (8)
to minimize the total operation cost over the planning horizon ΩT . The t∈ΩT
first term in (1) represents the cost of power exchange between the
0 ≤ bcurt m
i,t ≤ bt , ∀i, k ∈ ΩBUS , ∀t ∈ ΩT (9)
microgrid and the main grid as determined by (2). The second term,
given by (3), refers to the fuel cost consumed by the DDG units [15]. In the following subsections, the constraints corresponding to the
These two terms are restricted for the grid connected mode of the modeling of the microgrid, and the components located inside and
microgrid. The third term represents the cost of load curtailment during outside the microgrid (i.e., DDG, ESS, and substation) are described.
islanding mode of the microgrid, given by (4) [6]. The last term char
acterizes the voltage violation index used for improving the voltage 2.2.1. Power flow constraints
quality of the microgrid during both operation modes by (5). The three-phase power flow model is used to model the microgrid
∑[( )( ) ] given by (10)–(13). The model takes into account the mutual inductance
min CtPCC + CtDDG 1 − bmt + CtCURT bmt + CtVI (1)
X 3Φ
MINLP t∈Ω T
and interphase capacitance among the different phase [16].
4
Fig. 3. Feasible operation region of the DDG units in IS and GC mode.
⎛ ⎞
⎜ + PPCC PDDG ⎟
⎜ i,Ψ,t i,Ψ,t ⎟
⎛ ) ⎞ ⎜
( ⎜
⎛ ⎞ ⎟
⎟
∑ ∑ GΨβ
ij cos
δi,Ψ,t − δi,β,t ⎜ WT
Pi,t 3Φ PV 3Φ
+ Pi,t ⎟
Pi,Ψ,t = Vi,Ψ,t Vi,Ψ,t ⎝ ( ⎠ ⎜
) =⎜+⎝ ⎠/3 ⎟,
⎟ ∀i ∈ ΩBUS , ∀Ψ ∈ F , ∀t ∈ ΩT (10)
j∈ΩBUS β∈F +BΨβ
ij sin δi,Ψ,t − δi,β,t ⎜ +Pdch
i,t
3Φ
− ηESi Pchi,t
3Φ
⎟
⎜ ( ) ⎟
⎜ ⎟
⎝ − PLi,Ψ,t 1 − bcurt
t
⎠
⎛ ⎞
⎛ ⎞ ⎜ +QPCC DDG
⎟
( )
⎜ i,Ψ,t + Qi,Ψ,t ⎟
∑ ∑ − BΨβ
ij cos δi,Ψ,t − δi,β,t ⎜ ⎟
Qi,Ψ,t = Vi,Ψ,t Vi,Ψ,t ⎝ ( ) ⎠ WT 3Φ
= ⎜ +(Qi,t ES 3Φ
+ Qi,t )/3 ⎟, ∀i ∈ ΩBUS , ∀Ψ ∈ F , ∀t ∈ ΩT (11)
j∈ΩBUS β∈F +GΨβ
ij sin δi,Ψ,t − δi,β,t
⎜
⎝ ( ) ⎟
⎠
− QLi,Ψ,t 1 − bcurt
t
⎛ ( ) ⎞
Ψβ
∑( ) gi(k)j(k) cos δi(k),Ψ,t − δj(k),Ψ,t
Pbr
k,Ψ,t = Vi(k),Ψ,t Vi(k),β,t − Vj(k),β,t ⎝ Ψβ ( ) ⎠, ∀k ∈ ΩLIN , ∀Ψ ∈ F , ∀t ∈ ΩT (12)
β∈F +bi(k)j(k) sin δi(k),Ψ,t − δj(k),Ψ,t
⎛ ( ) ⎞
∑( ) gΨβ
i(k)j(k) sin δi(k),Ψ,t − δj(k),Ψ,t
Qbr
k,Ψ,t = Vi(k),Ψ,t Vi(k),β,t − Vj(k),β,t ⎝ Ψβ ( ) ⎠, ∀k ∈ ΩLIN , ∀Ψ ∈ F , ∀t ∈ ΩT (13)
β∈F − bi(k)j(k) cos δi(k),Ψ,t − δj(k),Ψ,t
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
( ) 2 ( ) 2 and reactive power flowing along the branch k in the microgrid [17]. It
(14)
br
Pbr
k,Ψ,t + Qbr
k,Ψ,t ≤ I k Vi(k),Ψ,t , ∀k ∈ ΩLIN , ∀Ψ ∈ F , ∀t ∈ ΩT is worth noting that, i(k) and j(k) denote bus i located at the upstream of
branch k and bus j for the downstream of the branch k, respectively. Eq.
The constraints (10) and (11) defines the power injected to the ith bus (14) is the quadratic constraints representing the limitation over the
of the microgrid which are the well-known node-based three phase transmission line ampacity.
power flow equations. Then, (12) and (13) define respectively the active
5
variable in the proposed problem. This allows us to consider uv,int to be a

variable to be determined in the proposed optimization problem. As
such, the primary droop and secondary control are considered with the
DDG model in the proposed model. Let fbus denotes the ith bus where the
secondary frequency control acts. The DDG unit considering three phase
microgrid can be modeled by (19)–(23):
⎧[ ( ) ]
⎪
⎪ ω MG p m m DDG 3Φ
⎪ − ki δΨ,t − π/2 /mi bt + (1 − bt )Pi,t
⎪ /3,Ψ = a
⎪
⎨[ ( ) ]
p
PDDG
i,Ψ,t = − kiω δMG m m DDG 3Φ
Ψ,t + π /6 /mi bt + (1 − bt )Pi,t /3,Ψ = b ,
⎪
⎪ [ ( ) ]
⎪
⎪
⎪ ω MG p m m
⎩ − ki δΨ,t − 7π/6 /mi bt + (1 − bt )Pi,t DDG 3Φ
/3,Ψ = c
∀i ∈ ΩDDG ,∀t ∈ ΩT (19)

[( ) q] m
QDDG
i,Ψ,t = Vi* − Vi,Ψ,t + uv,int m DDG
i,Ψ,t /ni bt + (1 − bt )Qi,t
3Φ
/3,
∀i ∈ ΩDDG , ∀Ψ ∈ F , ∀t ∈ ΩT (20)
Fig. 4. Principle of the ESS units. ∑
PDDG
i,t
3Φ
= PDDG
i,Ψ,t , ∀i ∈ ΩDDG , ∀t ∈ ΩT (21)
2.2.2. DDG operation constraints
Ψ∈F
The DDG units are typically interfaced with power electric con ∑
QDDG 3Φ
= QDDG
i,Ψ,t , ∀i ∈ ΩDDG , ∀t ∈ ΩT (22)
verters and passive output filters [17]. In this study, the hierarchical i,t
Ψ∈F
control structure is employed to control the operation of DDG units.
Fig. 2(a) illustrates the basic scheme of the DDG unit with hierarchical Vi,Ψ,t bmt = Vi** bmt , ∀i ∈ ΩDDG , ∀Ψ ∈ F , ∀t ∈ ΩT , (23)
control. The control structure involves multi-layers i.e., internal
voltage/current control, primary droop control, and secondary control The constraint (19) and (20) defines the active and reactive power
layers. The synchronous reference frame phase-locked loop (SRF-PLL) is generated by the DDG units. The first terms of the constraints refer to the
used to measure the frequency of the microgrid [13]. More detail about secondary control mode of DDGs in the islanding operation mode, while
the DDG’s control structure and their functionality can be found in their second terms represent the PQ mode for the grid connected oper
[5,18]. ation mode. The binary input parameter btm is used to activate the terms
In this work, the steady state model is adopted for dealing with the of (19) and (20) according to the microgrid operation mode. In this
operation constraints of the DDG units [17,13]. Fig. 2(b) shows the work, the active and reactive droop gains of the DDG units are deter
steady state model for the DDG units operated with the primary droop mined using the capacity-based method and expressed by (24) and (25).
and secondary controls. This model is sufficient for determining the By using this definition, the load in the system can be shared in pro
electrical variables of the PCC bus where the DDG unit is connected, portion to the capacity of the DDG units [6,12].
while respecting to the DDG unit’s internal circuit and the control
(24)
DDG
mpi = Δω/Pi , ∀i ∈ ΩDDG
structure.
As can be seen in the Fig. 2(b), the operation of the DDG unit, their
(25)
DDG
nqi = ΔV/2Qi , ∀i ∈ ΩDDG
internal circuits and its control structure can be modeled as a voltage
source where the frequency and the voltage are controlled via the droop Then, the constraints (21) and (22) defines the three-phase active
control and restored through the secondary control. Following this and reactive power flowing to the bus where the DDG is connected. The
definition, the DGG considering the primary droop and secondary con constraints (23) define the voltage magnitude restored by the secondary
trol in the steady state can be modeled as expressed by (15) and (16): control of the DDGs.
With the interfaced power electronic converter, the DDG units is
ω = ω* − mp PDDG + uω (15)
allowed to control the reactive power by producing to or receiving from
the microgrid [18,19]. Fig. 3 shows the feasible operating region of the
V = V * − nq QDDG + uv (16)
DDG units and the following constraints defines this region accurately.
where the output signals uω and uv are delivered through a proportional- P DDG
≤ PDDG 3Φ DDG
≤ Pi , ∀i ∈ ΩDDG , ∀t ∈ ΩT (26)
integral (PI) controller of the secondary layer and dependent to the
i i,t
measured frequency ωm and voltage V m as expressed by:

(27)
DDG DDG
∫ Q i ≤ QDDG
i,t
3Φ
≤ Qi , ∀i ∈ ΩDDG , ∀t ∈ ΩT
uω = kpω (ω* − ωm ) + kiω (ω* − ωm )dt (17) ( ) ( )
− PDDG
i,t
3Φ
tg cos− 1 PF DDG
i
4th
≤ QDDG
i,t
3Φ
≤ PDDG
i,t
3Φ
tg cos− 1 PF DDG
i
1st
,
∫ ∀i ∈ ΩDDG , ∀t ∈ ΩT (28)
uv = kpv (V ** − V m ) + kiv (V ** − V m )dt (18)
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
( )2 ( )2
(29)
DDG
As the proposed problem is formulated to model the steady state PDDG
i,t
3Φ
+ QDDG i,t
3Φ
≤ Si , ∀i ∈ ΩDDG , ∀t ∈ ΩT
condition, the proportional term for (17) and (18) can be eliminated (in
steady state: ω* − ωm = 0 and V ** − V m = 0) [13]. Nonetheless, the same where (26) and (27) limit the active and reactive power while the
approach could not be applied to the integrator part of the (17) and (18). constraint (28) guarantees the operation of the DDG unit respecting to
Let uω,int denotes the integrator part of (17) and uv,int for (18), as such the minimum power factor set. Finally, (29) outlines the boundary of the
uω = uω,int and uv = uv,int . As explained in [13], uω,int can be determined feasible region for limiting the active and reactive power with respect to
using uω,int = − kiω (δMG − π/2) under steady state condition. On the other the rating of the DDG. More details regarding the derivation for the
hand, in the steady state condition the voltage of the voltage-controlled operation feasible region can be found in the Appendix C.
bus V m is fixed by V ** = V m and there is no need to consider V m to be a
6
2.2.3. Energy storage operational constraints ( ) ( ) ∑ PCC

In this study, a general model is employed for the ESSs that is not PPCC
i,t
3Φ
1 − bmt = 1 − bmt Pi,Ψ,t , ∀i ∈ ΩPCC , ∀t ∈ ΩT (37)
limited to any particular energy storage technology. The ESSs are
Ψ∈F
interfaced with power conversion systems (PCSs) to control the ESSs ( ) ( ) ∑ PCC
QPCC 3Φ
1 − bmt = 1 − bmt Qi,Ψ,t , ∀i ∈ ΩPCC , ∀t ∈ ΩT (38)
active and reactive powers in accordance with the connected network i,t
Ψ∈F
requirements [20]. Accordingly, the ESS can be utilized to provide
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
)2̅
ancillary services such as voltage support, frequency regulation, grid ( )2 (
(39)
PCC 3Φ PCC 3Φ TR
stabilization, etc. In this work, the ESSs are supposed to operate at four- Pi,t + Qi,t ≤ Si , ∀i ∈ ΩPCC , ∀t ∈ ΩT
quadrant areas. To this end, the ESS active and reactive powers are
( ) ( )
limited via the quadratic function as illustrated in Fig. 4. PCC
− Pi,t 3Φ
tg cos− 1 PFTR
i
4th
≤ QPCC
i,t
3Φ
≤ PPCC
i,t
3Φ
tg cos− 1 PF TR
i
1st
,
Considering three phase microgrid, the ESS, therefore, can be model ∀i ∈ ΩPCC , ∀t ∈ NT (40)
using [20,21].
( ) ( )
Vi,Ψ,t 1 − bt m GC
= Vspec 1 − bmt , ∀i ∈ ΩPCC , ∀Ψ ∈ F , ∀t ∈ ΩT (41)
(30)
ES
0 ≤ Pdch
i,t
3Φ
≤ Pi bES
i,t , ∀i ∈ ΩEES , ∀t ∈ ΩT
( ) ⎧ ( )
⎪
⎪ 0 1 − bmt , Ψ = a
(31)
ES
0 ≤ Pch
i,t
3Φ
≤ Pi 1 − bES
i,t , ∀i ∈ ΩEES , ∀t ∈ ΩT ( ) ⎨ ( )
PCC
δi,Ψ,t 1 − bt = − 2π/3 1 − bmt , Ψ = b ,
m
∀i ∈ ΩPCC , ∀t ∈ ΩT (42)
⎪
⎪ ( )
⎩ +2π/3 1 − b , Ψ = c
m
(32)
ES ES
Q i ≤ QES
i,t
3Φ
≤ Qi , ∀i ∈ ΩESS , ∀t ∈ ΩT t
√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
)̅
( ) 2 ( 2 3. Linearize formulation
(33)
ES
Pdch
i,t
3Φ
− Pch
i,t
3Φ
+ QES
i,t
3Φ
≤ Si , ∀i ∈ ΩESS , ∀t ∈ ΩT
Solving the proposed MINLP model is computationally expensive
⎧ ( )
⎪ ES
⎨ SOC0i Ei + ηES ch 3Φ due to the nonlinear terms found in (1)–(37). In this section, the
i Pi,t − Pdch 3Φ
Δt, t = 1
nonlinear terms are linearized, and the proposed MILP model is derived.
i,t
ES
Ei,t = ( ) , ∀i ∈ ΩEES , ∀t ∈ ΩT
⎪ ES
⎩ Ei,t− ES ch 3Φ
− Pdch 3Φ
1 + ηi Pi,t i,t Δt, t ≥ 2
3.1. Linearization approach
(34)
(a) Linearization of Absolute Value Operator
(35)
ES ES ES
E i ≤ Ei,t ≤ Ei , ∀i ∈ ΩEES , ∀t ∈ ΩT
Eqs. (5) and (8) include absolute value operators which make the
(36)
ES
≥ SOC0i Ei , model non-linear. In order to linearize (5), the two following assump
ES
Ei,t ∀i ∈ ΩEES , t = NT
tions are made:
The ESS charging and discharging powers are restricted by their
rated power using (30)–(31). The output apparent power of the energy
(1) During the grid connected mode of the microgrid, the DDG units
storage, connected through an inverter, is restricted by the maximum
ES are operated in PQ mode and no DDG operates in PV mode in the
apparent power supplied/consumed of the inverter, denoted as (33). Si system similar to the work in [10,12]. Thus, there is no other type
The ESS stored energy in each time segment, EES i,t , is updated based on its GC
of resources can specify the voltage rather than Vspec imposed by
charging/discharging energy flows and the efficiency of the ESS (34). the upstream network. In this case, the voltage magnitude for all
The efficiency ηES the round trip efficiency of the energy storage systems the buses located at downstream of the bus which interlinked
[22,23]. Although the stored energy level of the ESS is not affected by with the upstream power grid (i.e., bus #1) are lower than VspecGC
.
the consumed/supplied reactive power of the connected inverter, the
(2) During the islanded mode of the microgrid, the DDG units are
amount of active power is dependent on it due to the line current limi
operated with secondary control layer. Then, a voltage restora
tations (33)–(34). The energy stored in the ESS is limited to its rated
tion reference V** is identically set for all DDG units allocated in
capacity and minimum stored energy (35). The ESSs are assumed to be
ES
the test system. In this case, the voltage of the system attempt to
initially charged up to SOC0i Ei and it must remain equal to or greater follow V**.
than the same amount of energy at the end of the planning horizon (36).
In this work, the ESS units are assumed to be operated in power- As there is no bus in the microgrid that its voltage magnitude could
controlled mode (PQ mode) formulated by (30)–(36) similar to the overcome Vspec GC
or V**, the absolute of voltage deviation can be omitted.
works in [6,12,20]. Nevertheless, as the ESS units are interfaced with
Hence, the linearized voltage deviation index (5) can be modified as:
power electronics inverters, the ESSs units are capable of operating in ( [ ])
droop mode and accordingly with the secondary control. A case study is ∑ ∑( ) spec m **
provided in subsection 4.5 to investigate the ESS operation with primary CtVI Lin =Vbase PCVI 1− bmt VGC +bt V − Vi,Ψ,t Δt, ∀t ∈ΩT
droop and secondary control layer.
i∈ΩBUS Ψ∈F
(43)
2.2.4. Power exchange constraints Nevertheless, if the aforementioned assumptions made to linearize
The constraints represent the power exchange between the micro (5) are not applicable, the method introduced in the Appendix D can be
grid, and the main-grid can be given by (37)–(42). The constraints (37) applied.
and (38) defines the active and reactive three phase power at the PCC of The same linearization method described in Appendix D can be
microgrid. Then, the power imported from the main grid were enforced adopted to linearize the absolute value operator employed in (8). In this
their limit by the rating of the transformer located in the substation and case, a new binary variable zcurti is introduced for representing the term
⃒ ⃒
given by the quadratic Eq. (39). The constraint (40) limits the operation ⃒ curt ⃒
⃒bi,t − bcurt
i,t− 1 ⃒. Then, the inequality constraints employed for linearizing
of the transformer above the minimum power factor allowable. More
over, (41) and (42) specified the voltage magnitude and their phase (8) can be expressed by:
∑
angle at the PCC which imposed by the main grid during grid connected zcurt curt
(44)
i,t ≤ 2Ni , ∀i ∈ ΩBUS
mode. t∈ΩT
7
nbus
∑ ( )
HijΨβ δi,Ψ,t − δj,β,t ≈ HijΨβ δΨ,t
j , ∀i, j ∈ ΩBUS , ∀Ψ, β ∈ F , ∀t ∈ ΩT (51)
j=1
where H can be replaced by G or B depending on (10) and (11); Cosine

and sine functions within (49) and (50) have already been eliminated by
Approximation #1. In this work, the frequency of the system is restored
to the nominal value due to the DDGs’ secondary control and its effect on
the lines impedances can be neglected. By applying Approximations #1
and #2, the three-phase power flow model can be linearized as:
⎛ ⎞
⎛ ⎞
Ψ,a:c Ψ,a:c a:ct,s,se
⎜ G ◦R1 .V j ⎟
⎜ ⎝ ij ⎠ ⎟
⎜ Ψ,a:c Ψ,a:c Ψt,s,se ⎟
n
∑ bus ⎜ +Gij .R2 δj ⎟
⎜ ⎟
Pi,Ψ,t = ⎜ ⎛ ⎞ ⎟
⎜ Ψ,a:c a:ct,s,se ⎟
j=1 ⎜ − BΨ,a:c
ij ◦R 1 .δ j ⎟
⎜ +⎝ ⎠⎟
⎝ Ψ,a:c Ψ,a:c ⎠
+Bij .R2
⎛ ⎞
⎜ PPCC + PDDG ⎟
⎜ i,Ψ,t i,Ψ,t ⎟
⎜ ⎛ ⎞ ⎟
⎜ WT 3Φ PV 3Φ ⎟
⎜ Pi,t + Pi,t ⎟
=⎜
⎜ + ⎝
dch 3Φ ch 3Φ
⎠/3 ⎟,
⎟ ∀i ∈ ΩBUS ,∀Ψ ∈ F ,∀t ∈ ΩT
Fig. 5. Polygon piecewise linearization plan for the quadratic function ⎜ +Pi,t − η P
ESi i,t ⎟
⎜ ⎟
approximation. ⎜ ( ) ⎟
⎝ − PLi,Ψ,t 1 − bcurt
t
⎠
bcurt curt curt

i,t − bi,t− 1 ≤ zi,t , ∀i ∈ ΩBUS , ∀t ∈ ΩT (45) (52)
⎛ ⎞
− bcurt curt curt
i,t + bi,t− 1 ≤ zi,t , ∀i ∈ ΩBUS , ∀t ∈ ΩT (46) ⎛ t,s,se
⎞
⎜ BΨ,a:c ◦RΨ,a:c .V a:c ⎟
⎜ − ⎝ ij 1 j
⎠ ⎟
i,t to be “1 when
The constraints (45) and (46) are adopted to force zcurt ′′ ⎜ t,s,se ⎟
nbus ⎜
∑ +BΨ,a:c .RΨ,a:c δΨj ⎟
there is a transition between load shedding and power supply for the ⎜ ij 2 ⎟
Qi,Ψ,t = ⎜ ⎛ ⎞ ⎟
⎜ ⎟
load in the bus i. Otherwise, zcurt is forced to be zero using (44)–(46).
t,s,se
i
j=1 ⎜ − GΨ,a:c
ij ◦RΨ,a:c
1 .δa:c
j ⎟
⎜ +⎝ ⎠ ⎟
⎝ ⎠
+GΨ,a:c .RΨ,a:c
(b) Linearization of power flow and line flow
ij 2
⎛ ⎞
In (10) and (11), with the production of voltage magnitude and ⎜ +QPCC DDG
i,Ψ,t + Qi,Ψ,t ⎟
⎜ ⎟
trigonometric i.e., sines and cosines the non-linearity has been intro ⎜ WT 3Φ ES 3Φ
+ Qi,t )/3
⎟
=⎜ +(Qi,t ⎟, ∀i ∈ ΩBUS , ∀Ψ ∈ F , ∀t ∈ ΩT
duced in these equations. In this case, two linear approximation is ⎜
⎝ ( ) ⎟
⎠
employed for these equations: − Qi,Ψ,t 1 − bcurt
L
t
Approximation #1: This approximation technique deals with the

(53)
non-linearity introduced by the trigonometric sine and cosine functions
and the deviation between the phase angle θΨβ as follows [16]: where V a:c a b c a:c
= [δaj δbj δcj ]. In this work, “
j = [Vj Vj Vj ] and δj ” and “.”
◦
ij
⎧
⎨ √̅̅̅ ( 1, if Ψ = β) respectively represent the Hadamard and scalar multiplications. RΨ,a:c
1
( )
cos δi,Ψ,t − δj,β,t ≈ − 1/2 −
√
3 /2 δ − δj,Ψ,t , if β lag to Ψ ,
̅̅̅ ( i,Ψ,t ) and RΨ,a:c
2 are the row vector of rotation based on cosine and sines of the
⎩
− 1/2 + 3 /2 δi,Ψ,t − δj,Ψ,t , if β lead to Ψ phase deviation, respectively. These vectors are as follows:
∀i, j ∈ ΩBUS , ∀Ψ, β ∈ F , ∀t ∈ ΩT (47) ⎡ ⎤
1 − 1/2 − 1/2
⎧
a:c,a:c
R1 ⎣
= − 1/2 1 − 1/2 ⎦ (54)
δ − δ , if Ψ=β
( ) ⎨ i,Ψ,t ( j,Ψ,t ) √̅̅̅ − 1/2 − 1/2 1
sin δi,Ψ,t − δj,β,t ≈ − 1/2 δi,Ψ,t − δj,Ψ,t + √3
̅̅̅
/2, if β lag to Ψ ,
⎩ ( ) ⎡ √̅̅̅ √̅̅̅ ⎤
− 1/2 δi,Ψ,t − δj,Ψ,t − 3 /2, if β lead to Ψ
√0̅̅̅ 3 /2 − √̅̅̅3 /2
∀i,j∈ΩBUS ,∀Ψ,β∈F ,∀t ∈ΩT (48) Ra:c,a:c
2 = ⎣− √̅̅̅3 /2 √0̅̅̅ 3 /2 ⎦ (55)
3 /2 − 3 /2 0
Approximation 2: This approximation technique deals with the sec
ond order of the voltage as well as the relation between the phase angle The same approximation techniques are also applied to the nonlinear
deviation and the susceptance term as follows [16]: constraints (12)–(13) and they are linearized as follows:
( ) ( ) ⎛ ⎞ ⎛ ⎞
Vi,Ψ,t Vj,β,t HijΨβ cos δi,Ψ,t − δj,β,t ≈ Vj,β,t HijΨβ cos δi,Ψ,t − δj,β,t , t,s,se t,s,se
⎜ GΨ,a:c ◦RΨ,a:c .V a:c ⎟ ⎜ GΨ,a:c ◦RΨ,a:c .V a:c ⎟
(49)
1 1
∀i, j ∈ ΩBUS , ∀Ψ, β ∈ F , ∀t ∈ ΩT ⎜ ij i
⎟ ⎜ ij j
⎟
⎜ t,s,se ⎟ ⎜ Ψ,a:c Ψt,s,se ⎟
Pbr
k,Ψ,t = − ⎜ +GΨ,a:c
ij .RΨ,a:c
2 δΨi Ψ,a:c
⎟ + ⎜ +Gij .R2 δj ⎟,
( ) ( ) ⎜ ⎟ ⎜ ⎟
⎝ t,s,se ⎠ ⎝ Ψ,a:c a:ct,s,se ⎠
Vi,Ψ,t Vj,β,t HijΨβ sin δi,Ψ,t − δj,β,t ≈ HijΨβ sin δi,Ψ,t − δj,β,t , ∀i, j ∈ ΩBUS , ∀Ψ, β − BΨ,a:c
ij ◦RΨ,a:c
1 .δa:c
i − BΨ,a:c
ij ◦R 1 .δ j
∈ F , ∀t ∈ ΩT
∀k ∈ ΩLIN , ∀Ψ ∈ F , ∀t ∈ ΩT (56)
(50)
8
Table 1
Summarize of piecewise linear constraint formation for each equipment.
Equipment operating Set of piecewise X Y Z set belong formulation
(quadrant) (name) (elements selected) #
DDG 1st/2nd LDDG

i { L1st 4th
n , Ln }
DDG
Pi,t 3Φ DDG
Qi,t 3Φ DDG
Si , ∀i ∈ ΩDDG (66)
ESS full LESS

i { L1st 2nd 3rd 4th
n , Ln , Ln Ln } Pdch
i,t
3Φ
− Pch
i,t
3Φ Qi,t
ES 3Φ ES
Si , ∀i ∈ ΩESS (67)
TR (PCC) full LPCC

n , Ln , Ln Ln } PPCC
i,t
3Φ
PPCC
i,t
3Φ TR
Si , ∀i ∈ ΩPCC (68)
LINE full Lbr
n , Ln , Ln Ln } Pbr
k,Ψ,t Qbr
k,Ψ,t
br
Ik Vi(k),Ψ,t , ∀k ∈ ΩLIN , ∀i ∈ ΩBUS (69)
where n = 1, 2, ⋯, npw ; Ψ ∈ F , ∀t ∈ NT .
Table 2
The characteristic of the DDG, WTs, PVs, ESS, transformer, and transmission line.
33-bus (balanced) 25-bus (unbalanced)
No. Bus Smax CRDDG

i
Mode PF Emax/Emin Imax No. Bus Smax CRDDG
i
Mode PF Emax/Emin Imax
(kVA) (cf/kWh) (GC/IS) (kWh) (A/Ψ) (kVA) (cf/kW) (GC/IS) (kWh) (A/Ψ)
DDGs 1 2 1500 11.105 PQ/2nd – – – 1 13 2000 11.148 PQ/2nd – – –

2 9 900 7.806 PQ/2nd – – – 2 19 1400 11.105 PQ/2nd – – –
3 8 1700 11.165 PQ/2nd – – – 3 25 2500 11.165 PQ/2nd – – –
4 22 500 7.316 PQ/2nd – – – – – – – – – – –
5 25 2000 11.418 PQ/2nd – – – – – – – – – – –
WTs 1 18 100 – PQunity 1.00 – – 1 12 100 – PQunity 1.00 – –
2 24 150 – PQunity 1.00 – – 2 22 200 – PQlag 0.95 – –
3 27 100 – PQlag 0.95 – – – – – – – – – –
4 32 100 – PQlag 0.95 – – – – – – – – – –
PVs 1 24 100 – PQunity 1.00 – – 1 5 150 – PQunity 1.00 – –
2 29 200 – PQunity 1.00 – – 2 15 150 – PQunity 1.00 – –
3 31 200 – PQunity 1.00 – – – – – – – – – –
ESSs 1 21 350 – PQ – 550/160 – 1 5 250 – PQ – 500/100
2 25 200 – PQ – 400/80 – 2 17 350 – PQ – 700/140
3 33 200 – PQ – 400/80 – – – – – – – – –
TR – 1 5300 – – – – – – – 5000 – – – – –
Line 1–8 – – – – – – 500 1–6 – – – – – – 700
9–32 – – – – – – 200 7–24 – – – – – – 350
Table 3
General parameter information.
Parameter Value Parameter Value Parameter Value
P DDG DDG
/Pi 0/Smax fbus /kiω 2/4 Q ES ES − 0.6Smax/+0.6Smax
i i /Qi
Q DDG DDG
/Qi − 0.2Smax/+0.6Smax PFDDG
i
1st
/PFDDG
i
4th 0.4/0.4 ηES
i
95 %
i
V /V* ** 1.01/1.05 σ fuel 0.007 $/Mcf [32] PFTR

i
1st
/PFTR
i
4th 0.75/0.75
⎛ ⎞ ⎛ ⎞ {
(59)
st st st
Ψ,a:c Ψ,a:c a:c t,s,se Ψ,a:c Ψ,a:c a:c t,s,se
Ln2nd = − α1n X + β1n Y − γ 1 Z ≤ 0 , n = 1, 2, ⋯, npw
⎜ − Bij ◦R1 .V i ⎟ ⎜ − Bij ◦R1 .V j ⎟
⎜ t,s,se ⎟ ⎜ t,s,se ⎟
Qbr = − ⎜ Ψ,a:c
⎜ − Bij .R2 δi
Ψ,a:c Ψ ⎟ + ⎜ − Bij .R
⎟ ⎜
Ψ,a:c Ψ,a:c Ψ
δj ⎟,
⎟
{
(60)
k,Ψ,t 2 th th th
⎝ t,s,se ⎠ ⎝ t,s,se ⎠ Ln3rd = − α4n X − β4n Y − γ4 Z ≤ 0 , n = 1, 2, ⋯, npw
− GΨ,a:c ◦RΨ,a:c .δa:c − GΨ,a:c ◦RΨ,a:c .δa:c
ij 1 i ij 1 j (57)
{
∀k ∈ ΩLIN , ∀Ψ ∈ F , ∀t ∈ ΩT (61)
th th th
Ln4th = α4n X − β4n Y − γ4 Z ≤ 0 , n = 1, 2, ⋯, npw
(c) Linearization of quadratic function where npw is the number of piecewise constraints for each quadrant. The
npw can be set high to improve the accuracy. Nonetheless, higher npw
The piecewise linearization technique is applied for the nonlinear increases the number of constraints and the computation burden is
quadratic constraints (14), (29), (33) and (39). Specifically in this work, accordingly increased. The piecewise constraint coefficients for each
multi-segment polygon region piecewise linearization technique is quadrant can be calculated as following:
employed [24–26]. This technique focuses on forming multiple con ( st th ( st th ) st th st th )
(62)
st th
straints representing the perimeters of circle for restricting the quadratic α1n ,4 = cos φ1n ,4 + n1 ,4 − 1 φ1sel,4 /n1pw,4 , n = 1, 2, ⋯, npw
summation. Fig. 5 illustrates the polygon piecewise technique. ( st th ( st th ) st th st th )
Let’s Ln and denote the nth piecewise constraints of the (63)
st ,4th
{1st,2nd,3rd,4th}
β1n = sin φ1n ,4 + n1 ,4 − 1 φ1sel,4 /n1pw,4 , n = 1, 2, ⋯, npw
quadrants {1st, 2nd, 3rd, 4th}. According to Fig. 5, each segment con
straints Ln can be formulated by: (64)
{1st,2nd,3rd,4th} st ,4th st th st th
γ1 = cos(φ1sel,4 /2n1pw,4 ), n = 1, 2, ⋯, npw
{ 1st 1st
(58)
1st st
Ln = αn X + βn Y − γ 1 Z ≤ 0 , n = 1, 2, ⋯, npw
9
Fig. 6. Total active power profile for the 25-bus test system: (a) SCEN 1 (b) SCEN 2 (c) SCEN 3.
Table 4
Optimization results obtained by MILP for 33-bus and 25-bus in difference scenarios.
Optimization goals 33-bus (balanced) 25-bus (unbalanced)
SCEN 1 SCEN 2 SCEN 3 SCEN 1 SCEN 2 SCEN 3
Total CPCC ($/24 h) 1123.22 0.00 805.45 993.63 0.00 706.89

Total CDDG ($/24 h) 0.00 5109.92 1676.71 0.00 4802.60 1603.03
Total CCURT ($/24 h) 0.00 149.94 145.35 0.00 829.60 414.80
Total CVD ($/24 h) 401.80 94.95 301.81 158.24 56.99 124.32
Total cost ($/24 h) 1525.03 5354.81 2929.32 1151.87 5689.20 2849.04
min{Vi,Ψ,t} (p.u.) 0.9716 1.0124 0.9719 0.9753 1.0131 0.9753
Avg. Runtime (s) 4.85 58.78 34.34 23.89 63.57 68.04
φ10
st ,4th st th st th
= φ1sel,4 /2n1pw,4 , n = 1, 2, ⋯, npw (65) 4. Numerical results
By applying the piecewise constraints to the quadratic functions of In this section, the well-known balanced 33-bus and unbalanced 25-
Section II and respecting to the operation characteristic of the equip bus test system have been chosen for testing the proposed MINLP and
ment, the formation of linear piecewise constraint associated with each MILP model. The characteristic of these test system can be found in [27].
equipment is summarized in Table 1. The planning horizon is set to 24 h representing a typical day. The daily
market electricity price and the percentage of annual peak load are
3.2. MILP model taken from NYISO [28] and [29]. The wind speed and solar irradiance
data are given in Appendix B. The load curtailment cost and the penalty
The proposed MILP model can be formulated by the set of linearized coefficient of the voltage violation are obtained from [30] and [31]. The
constraints and objective function as following: characteristic of the DDG, WTs, PVs, ESS, transformer, and transmission
∑ [( )( ) ] line are summarized in the Tables 2 and 3, respectively. The weighing
CtPCC + CtDDG 1 − bmt + CtCURT bmt + CtVI Lin (70)
coefficient of the load at ith bus and its maximum curtailment duration
min
X 3Φ
MILP t∈Ω
are assumed to be proportional to their annual peak and are summarized
DDG
subject to: in Table 8 of Appendix B. The parameters of the adopted piecewise

st th st th
linearization technique is set to φ1sel = φ4sel = 90◦ and n1pw = n4pw = 15
• Voltage and load curtailment constraints: (6), (7), (44), (45), (46), pieces which is high enough to obtain an acceptable accuracy [26]. The
(9). maximum and minimum voltage of the buses is set to 1.05 and 0.95 per
• Power flow constraints: (52), (53), (56), (57), (69). units, respectively.
• DDG operational constraints: (19), (20), (21), (22), (23), (26), (27), To study the effectiveness of the proposed ACOPF model considering
(28), (66). the microgrid operation in both grids connected (GC) and islanding (IS),
• ESS operational constraints: (30), (31), (32), (67), (34), (35), (36). three scenarios (SCEN 1, SCEN 2 and SCEN 3) are considered. First, in
• Power exchange constraints: (37), (38), (68), (40), (41), (42). SCEN 1 the microgrid is assumed to operates in grid connected mode all
along the considered time horizon. Then, SCEN 2 assume that the
The details of the decision variables of the proposed model are also microgrid operates only in islanding mode. Finally, SCEN 3 assume that
summarized in Appendix E. the microgrid is disconnected from the substation by opening the switch
ISW (i.e., btm = 1) between 8:00 and 14:00. In this work, the proposed
MILP model and the reference MINLP model have been built in MATLAB
10
Table 5
Comparison of results obtained by MILP to those obtained by NPL.
Error indices 33-bus (balanced) 25-bus (unbalanced)
NLP MILP Error NLP MILP Error
Total energy DDG (MWh) 23.094 22.858 1.02% 20.721 20.548 0.84%
Total energy curtailment (kWh) 44.1 44.1 0.00% 122 122 0%
Total energy losses (MWh) 1.907 1.930 1.20% 1.632 1.676 2.68%
min{Vi,Ψ,t} (p.u.) 0.9703 0.972 0.17% 0.9790 0.9753 0.38%
min{δi,Ψ,t} (deg.) − 0.50 − 0.496 0.14% − 120.66 − 120.87 0.17%
max{δi,Ψ,t} (deg.) 90.418 90.446 0.03% 211.31 211.36 0.02%
max{|Ibrk,Ψ,t|} (kA) 0.322 0.321 0.27% 0.640 0.643 0.35%
Avg. Runtime (min) 137.22 0.572 612.01 1.134
Fig. 9. Daily voltage magnitude profile of bus #25 obtained by the MILP2nd
Fig. 7. Voltage profile at 7:00 (IS mode) and at 8:00 (GC mode) for the un and MILPdroop for the 33-bus test system SCEN 2.
balanced 33-bus test case obtained by NPL and MILP.
environment and solved using “Intlingprog” and “Fmincon” solvers,

respectively. All the simulation are implemented on a PC computer
equipped with 3.40 GHz Quad-Core Intel(R) Core™ i7-6700 CPU with
32 GB of RAM.
4.1. Optimum operation of the microgrid
The MG optimal operation obtained by solving the proposed MILP Fig. 10. Voltage magnitude profile with and without secondary control for the
model for the three aforementioned scenarios (i.e., SCEN 1, SCEN 2 and 33-bus test system for SCEN 2 at 12:00 am.
SCEN 3) is presented. Fig. 6 shows the total three phase active power for
the generation, storage and load for the planning horizon. Table 4 shows
a comparison between SCEN 1, SCEN 2 and SCEN 3 in terms of costs and
voltage violation.
Due to the availability of the main grid during SCEN 1, the demand
can be fully supplied by the grid and the resources located in the test
system. Thus, there is no load curtailment in this scenario compared to
other two scenarios. On the other hand, 0.38% of energy demand and
0.19% of energy demand has to be respectively curtailed for the sce
narios SCEN 2 and SCEN 3, not only for balancing the power in the Fig. 11. Daily power flowing on branch #24 obtained by MILP2nd and MILP
droop
for the 33-bus test system for SCEN 2.
system but also for satisfying the secondary control of the DDG units
such as voltage restoration layer. The cost associated with the power
exchange with the main grid is zero in SCEN 2 as the microgrid is iso
lated from the main grid during this scenario. The PV and WT renewable
resources account for 3.91% of the total generated energy for the three
scenarios. In the grid connected mode, the DDG units are turned off due
to the high cost of fuel consumption compared to the cost of energy
imported from the main grid. Thus, no active power is generated by the
DDG units and its reactive power is accordingly bounded to zero by (28).
In other words, as it is not economical to generate active power by DDGs
Fig. 12. Daily system frequency obtained by the MILP2nd and MILPdroop for the
during grid connected operation mode, their reactive power generation
33-bus test system SCEN 2.
will accordingly be bounded to zero by (28). Hence, no reactive power is
Fig. 8. Voltage profile at 7:00 (IS mode) and at 8:00 (GC mode) for the unbalanced 25-bus test case obtained by NPL and MILP: (a) Phase A, (b) Phase B, (c) Phase C.
11

Fig. 13. Number of load shedding switching state cycle.
Fig. 14. Profile of load shedding at bus #19.
Fig. 15. Daily profile for (a) binary indicate charging/discharging state of ESS localized at bus #13 (b) voltage for the bus #13.
generated. As the DDG units are generating no reactive power in this

case, the cost of voltage deviation during the SCEN 1 is high. For the
same reason, the cost of voltage deviation in the SCEN 2 and SCEN 3 is
lower than that in SCEN 1.
Nevertheless, the total operation cost during SCEN 1 is always lower
than that in SCEN 2 and SCEN 3. That is why the IS mode is usually
operated during emergencies and/or when the main grid is not avail
able. In term of the minimum voltage, they are limited the minimum
Fig. 16. Active discharged power of ESS localized on bus #13 for both case ESS voltage magnitude constraints.
with and without considering droop + 2nd control.
12
Table 6
Optimization results for CASE2nd and CASEP. Vi,ta bmt = Vi,tb bmt = Vi,tc bmt , ∀i ∈ ΩDDG , ∀t ∈ ΩT , (75)
Indicators CASE2nd CASEP Finally, both MILP OPF models, with and without the secondary
Total C DDG
($/24 h) 5113.15 5113.15 control, are used to solve the 33-bus test system considering SCEN 2. The
Total CCURT ($/24 h) 0.00 0.00 results are derived and the effect of primary droop and secondary con
Total cost ($/24 h) 5188.24 5188.24 trols are evaluated accordingly. For the ease of notation, we denote the
min{Vi,Ψ=A,t} (p.u.) 1.014 1.013
proposed MILP model in which both primary and secondary control
min{Vi,Ψ=A,t} (p.u.) 1.05 1.05
min{δi,Ψ=A,t} (deg.) 88.000 88.113 layer have been considered as MILP2nd. Similarly, we denote the
max{δi,Ψ=A,t} (deg.) 90.779 90.609 modified model with only droop primary control layer as MILPdroop.
max{|Ibrk,Ψ=A,t|} (kA) 0.102 0.099 First, with the application of secondary control layer the voltage
Total energy loss (MWh) 0.912 0.896 magnitudes at the buses where the DDG units are connected are
restored. Fig. 9 shows the voltage magnitude profile for bus #25 where a
4.2. Accuracy evaluation DDG is connected.
In SCEN 2, the microgrid is operated in IS mode along the planning
For evaluating the accuracy, the results obtained by the proposed horizon 24 h. As can be seen in Fig. 9, with MILP2nd the voltage resto
MILP model are compared to those obtained by the original MINLP ration layer of secondary control is activated and the voltage magnitude
model. For the sake of the comparison with the solutions obtained by the at the bus #25 is restored to the secondary voltage control reference
proposed MILP, the MINLP model is changed to a nonlinear problem (V ** = 1.05p.u.). Nonetheless, the result derived by the MILPdroop, show
(NLP) by setting all binary variables according to those obtained by that the voltage magnitude at bus #25 follows the primary control
solving the proposed MILP model [6]. The third scenario (SCEN 3) is voltage droop with a reference V * = 1.01p.u.. Additionally, Fig. 10
used and the solution obtained via NLP and MILP as well as the relative shows the voltage profile for all the system buses during the at 12:00 am
error obtained by the solving the models are provided in Table 5. in SCEN 2.
Table 5 shows that the relative errors in the calculation of the total As can be seen in Fig. 10, the voltage magnitude profile derived by
energy produced by the DDG units are around 1%. In terms of power the MILP2nd follows V ** at the DDG units connection points, while the
loss, the maximum absolute error in the lines active power losses by lines voltage profile obtained using the MILPdroop model is always maintained
is lower than 3 × 10− 3 per unit. Nonetheless, the error is accumulated below V * by following the primary droop control. Fig. 11 shows the
when calculating the total energy losses along the considered time ho power flowing in the branch #24 i.e., |P24 |. This branch is connected
rizon and lead to 1.20% and 2.68% for 33-bus and 25-bus test systems, between bus #24 and #25. As can see, the power flows obtained by
respectively. MILPdroop are different compared to those obtained by MILP2nd.
On the other hand, the relative error of the voltage magnitude, the The secondary control has also been used for maintaining the system
voltage phase angle, and the current is lower than 0.5% which validates frequency to its nominal value and make the microgrid ready for syn
the accuracy of the proposed MILP model. Figs. 7 and 8 show the voltage chronization and reconnection to the main grid. The daily frequency of
profile at time t = 7:00 (the MG operate in GC mode) and t = 8:00 (the the microgrid system is shown in Fig. 12. The system frequency derived
MG operate in IS mode) obtained by proposed MILP and NPL for the 33- by the proposed MILP2nd model can be directly specified by ωt = 1.00p.
bus and 25-bus test systems, respectively. As Table 5 shows, the pro u. or e.g., 60 Hz in this work. Nonetheless, in the MILPdroop model with
posed MILP is computationally superior compared to NLP which is a only the primary droop control adoption, the system frequency is a
crucial for real-time operation and stochastic planning applications. variable to be found [12]. In this case, this frequency is dependent on the
active droop characteristic P − ω and the normal angular frequency, set
4.3. Primary droop and secondary control impact evaluation to ω* = 1.00p.u. in this work.
Given that the results of the MILP2nd closely match those of the
As discussed in the literature, the secondary control has been adop detailed MINLP model as depicted in Table 5, the results of the com
ted for compensating the steady-state error derived by the primary parison between MILP2nd OPF model and MILPdroop OPF model
droop control layer and to restore the microgrid’s voltage and frequency demonstrate that the MILPdroop OPF formulation is not sufficient for
[13]. In this subsection, the impact of considering the secondary control representing the actual behavior of the microgrid with hierarchal con
layer in the OPF formulation on the microgrid’s operation is evaluated. trol i.e., when a secondary control layer is employed.
To this end, the constraints (19)–(23) of the proposed MILP model
representing the operation of DDGs with secondary control are modified 4.4. Evaluation of load shedding switching state cycle
and transformed into the well-known droop-based power generation
model as expressed in (71)–(75) [12], i.e., only primary control. Then, In this section, in order to evaluate the impacts of constraints (44)–
the set of decision variables, included in Table 10 of Appendix E, are (46) employed to limit the number of load shedding cycles on the
modified by eliminating a set of uvi,Ψ,t and introducing the system fre microgrid operation, we consider a case in which it is assumed that the
quency ωt as variables to be solved [12]. DDGs #3 to #5 are disactivated. This will force a large portion of load to
be shed in order to balance the generation and demand. Fig. 13 shows
PDDG 3Φ
= [(ω* − ωt )/mpi ]bmt + (1 − bmt )PDDG 3Φ
/3, ∀i ∈ ΩDDG , ∀t ∈ ΩT
i,t i,t
the total number of switching cycles for each load i within 24 h for 33-
(71) bus test system in SCEN 2 with DDGs #3 to #5 disactivated.
[( ) ] As can be seen in Fig. 13, the maximum number of switching cycles
QDDG
i,t
3Φ
= Vi* − Vi,ta /nqi bmt +(1 − bmt )QDDG
i,t
3Φ
/3, ∀i ∈ ΩDDG ,∀t ∈ ΩT for the case without considering (44)–(46) are 5.5 cycles. It this neces
(72) sary to note that, 0.5 cycle represents a single switching state, e.g., on to
off or vice versa. Two more cases where the number of cycles is limited
∑
PDDG 3Φ
= PDDG ∀i ∈ ΩDDG , ∀t ∈ ΩT (73) by Ncurt = 1 and Ncurt = 3 are also provided. In these cases, the number
i,t i,Ψ,t ,
Ψ∈F of cycles has been limited to Ncurt = 1 and Ncurt = 3 cycles. The loads
∑ energy demand has been shed from the network by 26.95%, 27.27%,
QDDG
i,t
3Φ
= QDDG
i,Ψ,t , ∀i ∈ ΩDDG , ∀t ∈ ΩT (74) and 26.47% for the cases without (44)–(46), with Ncurt = 1, and with
Ψ∈F
Ncurt = 3, respectively. The bcurt
19,t for these cases is shown in Fig. 14. As
can be seen, the load in bus #19 is more often switched between the
13
shedding and supply states in the “no limit” case compared to other ∑
Pdch 3Φ ES m ESm
bi,t bt bi = bES m ESm
i,t bt bi Pdch
i,Ψ,t , ∀i ∈ ΩESS , ∀t ∈ ΩT (78)
cases. This situation might not be satisfying for the consumers/operators i,t
Ψ∈F
as there is no control over the load shedding cycles. Nonetheless, with
∑
(44)–(46) applied, the number of cycles is limited. ES
Qi,t 3Φ ES m ESm
bi,t bt bi = bES m ESm
QES (79)
i,t bi,t bi i,Ψ,t , ∀i ∈ ΩESS , ∀t ∈ ΩT
Ψ∈F
4.5. Evaluation of ESS with primary droop and secondary control

Vi,Ψ,t bES m ESm
i,t bt bi = Vi** bES m ESm
i,t bt bi , ∀i ∈ ΩESS , ∀Ψ ∈ F , ∀t ∈ ΩT , (80)
(a) Formulation for ESS with primary droop and secondary control
These set of constraints are derived from those implemented for the
DDG units operating with secondary control. Nonetheless, bES i,t intro
With the power electronics inverter interface, the ESS has the ability
duced in (76)–(80) to ensure that the constraints are activated only
to be operated in droop control mode with a supervisory secondary
control layer. In this mode, the ESS units participate with the dis when the ESS is in discharging state i.e., bES i,t = 1. Moreover, the binary
patchable resources microgrid to maintain the system frequency and the input parameter bm t introduced in this set of constraints guaranties that
voltage to the nominal value. The ESS operates in primary droop control the ESS operates in PQ mode during grid-connected operation mode.
mode with a supervisory secondary control layer when it discharges, The binary input parameter bESm i represents the operation mode for the
while during charging mode they act like a load and are as such operated ESS and is set to “0′′ for PQ mode and “1” for droop control mode with a
in the PQ mode. For considering the ESS operation in primary droop supervisory secondary control layer. The active and reactive droop gains
control mode with a supervisory secondary control layer, a new set of of the ESS units are determined using the capacity-based method similar
constraints expressed by (76)–(80) are added to the constraints (30)– to those used for the DDG unit. As can be seen in (76)–(80), there are
(36), representing the operation of the ESS in PQ mode: nonlinear terms in this set of constraints due to the multiplication be
⎧[ ( ) ] tween continuous variables and the binary variable bES i,t . In order to
⎪
⎪ ω MG p ES m ESm
⎪ − ki δΨ,t − π/2 /mi bi,t bt bi , Ψ = a resolve this nonlinearity, the linearization approach for the multiplica
⎪
⎪
⎨[ ( ) ] tion between continuous and binary variables introduced in Appendix D
p
Pdch ES m ESm
− kiω δMG ES m ESm
Ψ,t + π /6 /mi bi,t bt bi is adopted. Finally, as Pdch
i,Ψ,t and Qi,Ψ,t are introduced for each phase, the
i,Ψ,t bi,t bt bi = ,Ψ = b , ES
⎪
⎪ [ ( ) ]
⎪
⎪
⎪ ω MG p ES m ESm
⎩ − ki δΨ,t − 7π/6 /mi bi,t bt bi , Ψ = c power flow constraints are needed to be modified from (10) and (11) to
(81) and (82), while respecting to bES i,t , bt and bi
m ESm
as follows:
∀i ∈ ΩESS , ∀t ∈ ΩT (76)
[( ) q ] ES m ESm
QES ES m ESm
i,Ψ,t bi,t bt bi = Vi* − Vi,Ψ,t + uv,int
i,Ψ,t /ni bi,t bt bi , ∀i ∈ ΩESS , ∀Ψ
∈ F , ∀t ∈ ΩT
(77)
⎛ ⎞
⎜ PPCC DDG dch ES m ESm ⎟
⎜ i,Ψ,t + Pi,Ψ,t + Pi,Ψ,t bi,t bt bi ⎟
⎛ ⎞ ⎜ ⎛ ⎞ ⎟
( ) ⎜ ⎟
∑ ∑ GΨβ ⎜ PWT 3Φ
+ PPV 3Φ
⎟
ij cos δi,Ψ,t − δi,β,t i,t i,t
Pi,Ψ,t = Vi,Ψ,t Vi,Ψ,t ⎝ ( ) ⎠=⎜
⎜ +
⎜
⎝ ( ) ⎟ ⎟
⎠/3 ⎟, ∀i ∈ ΩBUS , ∀Ψ ∈ F , ∀t ∈ ΩT (81)
j∈ΩBUS β∈F +BΨβ
ij sin δi,Ψ,t − δi,β,t
⎜
⎜ +Pdch
i,t
3Φ
1 − bES m ESm
i,t bt bi − ηESi Pch
i,t
3Φ ⎟
⎟
⎜ ( ) ⎟
⎜ ⎟
⎝ − PLi,Ψ,t 1 − bcurt t ⎠
⎛ ⎞
⎜ ⎟
⎜ +QPCC + QDDG ES ES m ESm
i,Ψ,t + Qi,Ψ,t bi,t bt bi
⎟
⎜ ⎛ i,Ψ,t ⎞ ⎟
⎛ ( ) ⎞ ⎜ ⎟
⎜ QWT 3Φ ⎟
∑ ∑ − BΨβ
ij cos δi,Ψ,t − δi,β,t ⎜ ⎜ i,t ⎟
Qi,Ψ,t = Vi,Ψ,t Vi,Ψ,t ⎝ ) ⎠=⎜ ( )⎟ ⎟, ∀i ∈ ΩBUS , ∀Ψ ∈ F , ∀t ∈ ΩT (82)
Ψβ
( ⎜ +⎝ +QEs 3Φ 1 − bES bm bESm ⎠/3 ⎟
j∈ΩBUS β∈F +G ij sin δ i,Ψ,t − δi,β,t ⎜ i,t i,t t i ⎟
⎜ ⎟
⎜ ( ) ⎟
⎜ − Qi,Ψ,t 1 − bcurt
L ⎟
⎝ t ⎠
14
(b) Evaluation of ESS with primary droop and secondary control microgrid with hierarchical control structure was proposed. The model
considers the microgrid philosophy of operation during GC and IS in a
In this subsection, the operation of the ESS in primary droop control single comprehensive model. The energy storage system and the inter
mode with a supervisory secondary control layer is evaluated. A new mittent resources such as WTs and PVs were considered in the proposed
ESS is added into microgrid with the capability to operate in both the model. Additionally, the proposed MINLP model was linearized and the
power control mode or in primary droop control mode with a supervi MILP model of the problem was obtained which can be solved efficiently
sory secondary control layer. As such, two new case studies, CASEP and with commercial solvers. The obtained solutions can guarantee the
CASE2nd, are introduced for the proposed MILP model. In CASEP, the ESS optimal operation of the microgrid respecting to their philosophy of
is operated power control PQ mode while in CASE2nd it is operated in operation during both operation mode. Numerical simulation studies
droop control mode with a supervisory secondary control layer. Both were conducted to prove the effectiveness of the proposed models. The
case studies are used to solve the 33-bus test system considering SCEN2. error assessment was also conducted for validating the proposed model.
The ESS is allocated to bus #13 of the test system with the following The results showed that the proposed MILP model achieves the optimum
ES ES
characteristic: Si = 1.2MVA, Ei = 4.5MWh, E ES results with significantly lower computational time compared to that in
i = 0.9MWh, and
ηES = 95%. the MINLP. Such performance is crucially needed in the real-time and
i
stochastic planning application for the microgrid.
Fig. 15 shows the charging/discharge state of the new ESS allocated
at bus #13 as well as the voltage for this bus. As can see at a particular
period of time i.e., 12 h00-13 h00, the ESS operated in discharging mode
CRediT authorship contribution statement
for both cases i.e., the ESS with and without the primary droop control
mode with a supervisory secondary control layer. As expected, the ESS
Youthanalack Vilaisarn: Conceptualization, Data curation, Formal
operating in droop control mode with a supervisory secondary control
analysis, Funding acquisition, Investigation, Methodology, Resources,
layer can maintain the voltage refers to their nominal voltage setting
Software, Validation, Writing – original draft, Writing – review & edit
V ** . Nonetheless, with ESS in PQ mode the voltage for the bus #13 in
ing. Majid Moradzadeh: Conceptualization, Data curation, Formal
this case is below V ** of the DDG units. For the rest of period, the ESS for analysis, Funding acquisition, Investigation, Methodology, Resources,
both cases were operated similarly. Software, Validation, Writing – original draft, Writing – review & edit
Fig. 16 shows the active power profile for a particular period of ESS ing. Morad Abdelaziz: Conceptualization, Data curation, Formal anal
operating in discharging mode i.e., 12 h00-13 h00. When the ESS oper ysis, Funding acquisition, Investigation, Methodology, Resources,
ates in PQ mode, the active power has been scheduled to inject its power Software, Validation, Writing – original draft, Writing – review & edit
into the network according to the optimal operation of the microgrid. ing. Jérôme Cros: Conceptualization, Data curation, Formal analysis,
Nonetheless, when the ESS operates in droop control mode with a su Funding acquisition, Investigation, Methodology, Resources, Software,
pervisory secondary control layer, the ESS injects the specific amount of Validation, Writing – original draft, Writing – review & editing.
active power by Pdch dch
13,A,t=12 = 566.81kW and P13,A,t=13 = 559.96kW cor
responding to the droop settings under δMG MG
A,t=12 = 89.932 and δA,t=13 =
◦
89.933 , respectively.
◦
Declaration of Competing Interest
The overall optimization results have been summarized in Table 6.
The results show that the total cost is not dependent on the operating The authors declare that they have no known competing financial
mode of the ESS. For the same reason, no load is curtailed in both cases. interests or personal relationships that could have appeared to influence
Nonetheless, the minimum voltage for CASE2nd is slightly higher than the work reported in this paper.
that in CASEP due to primary droop and the secondary control adopted.
Finally, the maximum branch current and total energy loss were higher Acknowledgments
in CASE2nd compared to those obtained in CASEP.
This work was partially supported by Programme Canadien de
5. Conclusion Bourses de la Francophonie (PCBF) and NSERC (Canada). The author
Youthanalack Vilaisarn especially thanks the PCBF for their support of
An MINLP model for the optimal operation problem of unbalanced his Doctoral Program studies.
Appendix A
In this work, the wind units (WTs) and photovoltaic units (PVs) are interfaced with the power electronic converter and able to track a MPPT power
for the generation. As such, the WTs and PVs can be modeled as the sources of static injection power into the microgrid where their amount of power
generated are respected to the capacity windspeed and irradiance, respectively. In this way the power generated via these resources can be calculated
using the static model [14] and describe as following:
⎧
⎪
⎪ 0, vmt < vci , vmt ≥ vco
⎨ ( m ) ra
(83)
WT 3Φ WT
Pi,t = Si PF i vt − v /(v − v ), vci ≤ vmt ≤ vra
WT ci ci
, ∀i ∈ ΩWT , ∀t ∈ NT
⎪
⎪
⎩ WT WT ra m co
Si × PFi , v < vt < v
15
then, the reactive power generated by the WTs can be calculated respecting to their specified setting power factor. The formulation used for this
calculation are illustrated as following:
( ( ))
QWT
i,t
3Φ
= PWT
i,t
3Φ
tg cos− 1 PF WT
i , ∀i ∈ ΩWT , ∀t ∈ NT (84)
next, the PVs assumed controlled in the power factor unity mode in this work. As such, there is no reactive power participated via the PV units. The
static model used for the PV units can be:
⎧
⎪ PV ( m )2 STD C m C
⎪
⎪ Si rt /(R − R ), rt ≤ R
⎨
PV ( m )2
PV 3Φ
Qi,t = Si rt /R , STD
R < rt < RSTD ,
C m ∀i ∈ ΩPV , ∀t ∈ NT (85)
⎪
⎪
⎪
⎩S ,
PV m STD
i rt ≥ R
In term of the demand, the load are modeled by considering the pattern percentage of the annual peak demands [29]. In this way, load at the time t
can be calculated though the multiplication of their peak with the annual peak load factor along the planning horizon.
PLi,Ψ,t = (APLt /100) × PLi,Ψpeak , ∀i ∈ ΩBUS , Ψ ∈ F , ∀t ∈ NT (86)
QLi,Ψ,t = (APLt /100) × QLi,Ψpeak , ∀i ∈ ΩBUS , Ψ ∈ F , ∀t ∈ NT (87)
Appendix B
See Tables 7 and 8.
Table 7
Daily information data input for the 33-bus and 25-bus test system.
Time σPCC
t
APLt wt rt Time σPCC
t
APLt wm
t rm
t Time σPCC
t
APLt wm
t rm
t
(h) ($/MWh) (%) (m/s) (kW/m2) (h) ($/MWh) (%) (m/s) (kW/m2) (h) ($/MWh) (%) (m/s) (kW/m2)
1:00 15.63 64.0 8.1 0 9:00 12.10 88.0 6.2 0.592 17:00 18.48 95.0 2.2 0.629
2:00 15.50 60.0 7.6 0 10:00 12.87 95.0 7.2 0.752 18:00 20.03 95.0 0.5 0.368
3:00 14.63 58.0 7.6 0 11:00 14.01 98.0 7.1 0.868 19:00 21.85 93.0 2.2 0
4:00 13.28 55.0 6.9 0 12:00 14.88 100.0 7.5 0.893 20:00 20.92 92.0 2.6 0
5:00 12.88 55.0 6.9 0 13:00 15.33 98.0 7.1 0.916 21:00 20.06 92.0 4.8 0
6:00 12.31 57.0 7.0 0.391 14:00 16.21 100.0 4.8 0.206 22:00 19.17 93.0 6.2 0
7:00 12.07 64.0 7.0 0.212 15:00 17.21 100.0 3.8 0.390 23:00 17.41 88.0 4.8 0
8:00 11.85 76.0 7.3 0.371 16:00 17.26 96.0 3.3 0.117 24:00 16.15 72.0 5.1 0
Table 8
load weight and maximum duration of load curtailment 33-bus and 25-bus test system.
33-bus (balanced) 25-bus (unbalanced)
Bus wload
i
Tcurt
i Bus wload
i
Ticurt Bus wload
i
Ticurt Bus wload
i
Tcurt
i Bus wload
i
Ticurt Bus wload
i
Ticurt
(h) (h) (h) (h) (h) (h)
1 0 24 12 0.016 21 23 0.024 19 1 0.000 24 12 0.047 15 23 0.056 15

2 0.027 19 13 0.016 21 24 0.113 4 2 0.000 24 13 0.033 17 24 0.033 17
3 0.024 19 14 0.032 18 25 0.113 4 3 0.033 17 14 0.047 15 25 0.056 15
4 0.032 18 15 0.016 21 26 0.016 21 4 0.047 15 15 0.124 4
5 0.016 21 16 0.016 21 27 0.016 21 5 0.037 17 16 0.037 17
6 0.016 21 17 0.016 21 28 0.016 21 6 0.037 17 17 0.037 17
7 0.054 14 18 0.024 19 29 0.032 18 7 0.000 24 18 0.037 17
8 0.054 14 19 0.024 19 30 0.054 14 8 0.037 17 19 0.056 15
9 0.016 21 20 0.024 19 31 0.04 16 9 0.056 15 20 0.033 17
10 0.016 21 21 0.024 19 32 0.057 14 10 0.033 17 21 0.037 17
11 0.012 21 22 0.024 19 33 0.016 21 11 0.042 17 22 0.047 15
16
Appendix C
The operation range of the DDGs shown in Fig. 3 has been derived with respect to the feasible operating regions of DDGs [19]. In this work, it is
assumed the DDGs are generating units and allowed to generate active power and to exchange reactive power using the interfaced power electronic
devices. Fig. 17 shows the procedure executed for extracting the DDGs’ feasible operation region. Table 9 summarizes the relation between the
Fig. 17. Derivation of DDG feasible operation region.
17
Table 9
Constraints and parameters related to the determination of DDG’s feasible operation region.
Constraint # Fig. 17 Parameters
DDG DDG DDG DDG
P i Pi Q i Qi PFDDG
i
1st
PFDDG
i
4th
Si
DDG
(26) (a) 0 [0, R] – – – – –

(27) (b) – – [ − R, 0] [0, R] – – –
(28) (c) – – – – [0, R] [0, R] –
(29) (d) – – – – – – [0, R]
constraints, Fig. 17, and their relative parameters for determining the DDG’s feasible operation region.
The DDG i active power P DDGi in (26) is set to as low as zero as the DDGs are supposed to be operating as a source of active power. As such, the left-
DDG
hand side of the operating region is excluded as shown in Fig. 17(a)–(d). The Pi in (26), on the other hand, can be set according to the DDGs
DDG
maximum allowable power generation to adjust the right-hand side operating region as illustrated by Fig. 17(a). Nonetheless, Pi is usually set to
DDG
Si for allowing the DDGs to generate their active power proportional to their size. Fig. 17(b) illustrates the impact of the reactive power limits,
DDG
Q DDG
i and Qi given by (27), on the DDGs feasible operation region. They restrict the DDGs’ reactive power generated to and received from the main
grid, respectively. Hence, the DDGs feasible operating region is restricted. The power factor PFDDG
i
1st
and PFDDG
i
4th
in (28) limit the DDG’s minimum
generated active power proportional to its generated and received reactive power, respectively. As shown in Fig. 17(c), PFDDG i
1st
and PFDDG
i
4th
adjust
the angle of the radius located in the first and fourth quadrant of the DDG’s operating circle. Finally, Fig. 17(d) demonstrates the impact of the
DDG
quadratic power constraint (29). As can be seen, increasing Si expands the perimeter of the DDGs feasible operation circle. Finally, the feasible
operation region for the DDGs can be determined by the constraints (26)–(29) as shown in Fig. 17(e).
Appendix D
(a) Linearization of Absolute Value Operator
In order to linearize the nonlinear expression z = |y| in the problem with objective function Min|y|, where y is a continuous variable and the |.| is the
absolute value operator, the auxiliary constraints (89)–(90) are considered and the problem is reformulated as follows [33]:
Minz (88)
y≤z (89)
− y≤z (90)
where (88)–(89) perform the function of the absolute value operator for y ≥ 0, while (88) and (90) do the same for y ≤ 0.
(b) Linearization of the production between continuous and binary variables
Let z = A × b be the expression of the production between the continuous variable (A) and the binary variable (b), if A is bounded between zero and
the big M i.e., A, then, a representing variable z can be introduced into the problem formulation subjected to (91)–(94) [21]:
z≤A×b (91)
z≤A (92)
z ≥ A − (1 − b)A (93)
0≤z≤A (94)
Appendix E
The decision variables of the model are formed by the combination of the groups of vectors of control and dependent variables which are illustrated
in the Table 10. Here it is worth noting that the decision vector is the same for both IS and GC operation modes. However, the relevant decision
variables, within the decision vector, are activated/disactivated by using the binary input parameter bm t in the proposed problem (see Table 10).
18
Table 10
Control and dependent variables of the MINLP and MILP model for minimizing their objective functions.
Group Variables Type Set X3Φ
MILP X3Φ
MINLP
Operation mode
# ∀i, k ∈ IS GC
1 Vi,Ψ,t ,Ψ∈F , ∀t ∈ ΩT R ΩBUS ✓ ✓ ✓ ✓

2 δi,Ψ,t ,Ψ∈F , ∀t ∈ ΩT R ΩBUS ✓ ✓ ✓ ✓
3 bcurt
i,t , ∀t ∈ ΩT Z ΩBUS ✓ ✓ ✓
4 zcurt
i,t , ∀t ∈ ΩT Z ΩBUS ✓ ✓
5 PDDG
i,t
3Φ , ∀t ∈ ΩT R ΩDDG ✓ ✓ ✓ ✓
6 PDDG
i,Ψ,t
,Ψ∈F , ∀t ∈ ΩT R ΩDDG ✓ ✓ ✓ ✓
7 QDDG
i,t
3Φ , ∀t ∈ ΩT R ΩDDG ✓ ✓ ✓ ✓
8 QDDG
i,Ψ,t
,Ψ∈F , ∀t ∈ ΩT R ΩDDG ✓ ✓ ✓ ✓
9 uv,int ,Ψ∈F , ∀t ∈ ΩT R ΩDDG ✓ ✓ ✓
i,Ψ,t
10 Pdch
i,t
3Φ , ∀t ∈ ΩT R ΩESS ✓ ✓ ✓ ✓
11 Pch
i,t
3Φ , ∀t ∈ ΩT R ΩESS ✓ ✓ ✓ ✓
12 QES
i,t
3Φ , ∀t ∈ ΩT R ΩESS ✓ ✓ ✓ ✓
13 EES
i,t
, ∀t ∈ ΩT R ΩESS ✓ ✓ ✓ ✓
14 bES
i,t
, ∀t ∈ ΩT Z ΩESS ✓ ✓ ✓ ✓
15 PPCC
i,t
3Φ , ∀t ∈ ΩT R ΩPCC ✓ ✓ ✓
16 PPCC
i,Ψ,t
,Ψ∈F , ∀t ∈ ΩT R ΩPCC ✓ ✓ ✓
17 QPCC
i,t
3Φ , ∀t ∈ ΩT R ΩPCC ✓ ✓ ✓
18 QPCC
i,Ψ,t
,Ψ∈F , ∀t ∈ ΩT R ΩPCC ✓ ✓ ✓
19 Pbr
k,Ψ,t
,Ψ∈F , ∀t ∈ ΩT R ΩLIN ✓ ✓ ✓ ✓
20 Qbr
k,Ψ,t
,Ψ∈F , ∀t ∈ ΩT R ΩLIN ✓ ✓ ✓ ✓
References [15] Abdelaziz MMA, Farag HE, El-Saadany EF. Optimum Reconfiguration of Droop-
Controlled Islanded Microgrids. IEEE Trans Power Syst 2016;31(3):2144–53.
https://doi.org/10.1109/TPWRS.2015.2456154.
[1] Abdelaziz MMA. OpenCL-Accelerated Probabilistic Power Flow for Active
[16] Wang Yi, Zhang N, Li H, Yang J, Kang C. Linear three-phase power flow for
Distribution Networks. IEEE Trans Sustain Energy 2018;9(3):1255–64. https://doi.
unbalanced active distribution networks with PV nodes. CSEE J Power Energy Syst
org/10.1109/TSTE.2017.2781148.
2017;3(3):321–4.
[2] Abdelaziz MMA. Effect of Detailed Reactive Power Limit Modeling on Islanded
[17] Abdelaziz MMA, Farag HE, El-Saadany EF, Mohamed Y-R. A Novel and Generalized
Microgrid Power Flow Analysis. IEEE Trans Power Syst 2016;31(2):1665–6.
Three-Phase Power Flow Algorithm for Islanded Microgrids Using a Newton Trust
https://doi.org/10.1109/TPWRS.2015.2412690.
Region Method. IEEE Trans Power Syst 2013;28(1):190–201. https://doi.org/
[3] Rodrigues YR, Abdelaziz M, Wang L. D-PMU based secondary frequency control for
10.1109/TPWRS.2012.2195785.
islanded microgrids. IEEE Trans Smart Grid 2020;11(1):857–72. https://doi.org/
[18] Guerrero JM, Chandorkar M, Lee T-L, Loh PC. Advanced control architectures for
10.1109/TSG.516541110.1109/TSG.2019.2919123.
intelligent microgridspart i: Decentralized and hierarchical control. IEEE Trans Ind
[4] Abdelaziz MMA, Farag HE. An Enhanced Supervisory Control for Islanded
Electron 2013;60(4):1254–62. https://doi.org/10.1109/TIE.2012.2194969.
Microgrid Systems. IEEE Trans Smart Grid 2016;7(4):1941–3. https://doi.org/
[19] Macedo LH, Franco JF, Rider MJ, Romero R. Optimal Operation of Distribution
10.1109/TSG.2016.2553580.
Networks Considering Energy Storage Devices. IEEE Trans Smart Grid 2015;6(6):
[5] Rocabert J, Luna A, Blaabjerg F, Rodríguez P. Control of power converters in AC
2825–36.
microgrids. IEEE Trans Power Electron 2012;27(11):4734–49. https://doi.org/
[20] Mohseni-Bonab SM, Kamwa I, Moeini A, Rabiee A. Voltage Security Constrained
10.1109/TPEL.2012.2199334.
Stochastic Programming Model for Day-Ahead BESS Schedule in Co-Optimization
[6] Vergara PP, Lopez JC, Rider MJ, da Silva LCP. Optimal Operation of Unbalanced
of T&D Systems. IEEE Trans Sustain Energy 2020;11:391–404. https://doi.org/
Three-Phase Islanded Droop-Based Microgrids. IEEE Trans Smart Grid 2019;10(1):
10.1109/TSTE.2019.2892024.
928–40. https://doi.org/10.1109/TSG.2017.2756021.
[21] Moradzadeh M, Abdelaziz MMA. A New MILP Formulation for Renewables and
[7] Olivares DE, Mehrizi-Sani A, Etemadi AH, Canizares CA, Iravani R, Kazerani M,
Energy Storage Integration in Fast Charging Stations. IEEE Trans Transp Electrif
et al. Trends in microgrid control. IEEE Trans Smart Grid 2014;5(4):1905–19.
2020;6(1):181–98. https://doi.org/10.1109/TTE.668731610.1109/
https://doi.org/10.1109/TSG.2013.2295514.
TTE.2020.2974179.
[8] Wu D, Tang F, Dragicevic T, Vasquez JC, Guerrero JM. A Control Architecture to
[22] Alsaidan I, Khodaei A, Gao W. A Comprehensive Battery Energy Storage Optimal
Coordinate Renewable Energy Sources and Energy Storage Systems in Islanded
Sizing Model for Microgrid Applications. IEEE Trans Power Syst 2018;33(4):
Microgrids. IEEE Trans Smart Grid 2015;6(3):1156–66. https://doi.org/10.1109/
3968–80. https://doi.org/10.1109/TPWRS.5910.1109/TPWRS.2017.2769639.
TSG.2014.2377018.
[23] Moradzadeh M, Abdelaziz MMA. A Stochastic Optimal Planning Model for Fully
[9] Wu X, Wang X, Qu C. A hierarchical framework for generation scheduling of
Green Stand-Alone PEV Charging Stations. IEEE Trans Transp Electrif 2021;7(4):
microgrids. IEEE Trans Power Deliv 2014;29(6):2448–57. https://doi.org/
2356–75. https://doi.org/10.1109/TTE.2021.3069438.
10.1109/TPWRD.2014.2360064.
[24] Yang Z, Zhong H, Xia Q, Bose A, Kang C. Optimal power flow based on successive
[10] Vergara PP, López JC, da Silva LCP, Rider MJ. Security-constrained optimal energy
linear approximation of power flow equations. IET Gener Transm Distrib 2016;10
management system for three-phase residential microgrids. Electr Power Syst Res
(14):3654–62. https://doi.org/10.1049/gtd2.v10.1410.1049/iet-gtd.2016.0547.
2017;146:371–82. https://doi.org/10.1016/j.epsr.2017.02.012.
[25] Yang Z, Bose A, Zhong H, Zhang N, Xia Q, Kang C. Optimal Reactive Power
[11] Riva Sanseverino E, Nguyen Quang N, Di Silvestre ML, Guerrero JM, Li C. Optimal
Dispatch With Accurately Modeled Discrete Control Devices: A Successive Linear
power flow in three-phase islanded microgrids with inverter interfaced units. Electr
Approximation Approach. IEEE Trans Power Syst 2017;32(3):2435–44. https://
Power Syst Res 2015;123:48–56. https://doi.org/10.1016/j.epsr.2015.01.020.
doi.org/10.1109/TPWRS.2016.2608178.
[12] Vergara PP, Rey JM, Lopez JC, Rider MJ, da Silva LCP, Shaker HR, et al.
[26] Yang Z, Zhong H, Bose A, Zheng T, Xia Q, Kang C. A Linearized OPF Model with
A Generalized Model for the Optimal Operation of Microgrids in Grid-Connected
Reactive Power and Voltage Magnitude: A Pathway to Improve the MW-Only DC
and Islanded Droop-Based Mode. IEEE Trans Smart Grid 2019;10(5):5032–45.
OPF. IEEE Trans Power Syst 2018;33(2):1734–45. https://doi.org/10.1109/
https://doi.org/10.1109/TSG.516541110.1109/TSG.2018.2873411.
TPWRS.5910.1109/TPWRS.2017.2718551.
[13] Agundis-Tinajero G, Segundo-Ramírez J, Visairo-Cruz N, Savaghebi M,
[27] Vilaisarn Y, Abdelaziz M. An inversion-free sparse Z power flow algorithm for
Guerrero JM, Barocio E. Power flow modeling of islanded AC microgrids with
large-scale droop controlled islanded microgrids. Int J Electr Power Energy Syst
hierarchical control. Int J Electr Power Energy Syst 2019;105:28–36. https://doi.
2020;121:106048. https://doi.org/10.1016/j.ijepes.2020.106048.
org/10.1016/j.ijepes.2018.08.002.
[28] New York Independent System Opearator; n.d. [Online]. Available: http://www.
[14] Abdelaziz M, Moradzadeh M. Monte-Carlo simulation based multi-objective
nyiso.com.
optimum allocation of renewable distributed generation using OpenCL. Electr
Power Syst Res 2019;170:81–91. https://doi.org/10.1016/j.epsr.2019.01.012.
19
[29] Barani M, Aghaei J, Akbari MA, Niknam T, Farahmand H, Korpas M. Optimal [31] Samadi Gazijahani F, Salehi J. Optimal Bilevel Model for Stochastic Risk-Based
Partitioning of Smart Distribution Systems Into Supply-Sufficient Microgrids. IEEE Planning of Microgrids Under Uncertainty. IEEE Trans Ind Informatics 2018;14(7):
Trans Smart Grid 2019;10(3):2523–33. https://doi.org/10.1109/ 3054–64. https://doi.org/10.1109/TII.2017.2769656.
TSG.516541110.1109/TSG.2018.2803215. [32] U.S. Energy Information Administration; 2020. [Online]. Available: https://www.
[30] Gazijahani FS, Salehi J. Robust Design of Microgrids with Reconfigurable Topology eia.gov/dnav/ng/hist/n3035us3m.htm.
under Severe Uncertainty. IEEE Trans Sustain Energy 2018;9(2):559–69. https:// [33] Dimitris Bertsimas. Introduction to Linear Optimization; 1997.
doi.org/10.1109/TSTE.2017.2748882.
20

1 s2.0 S0142061521009042 Main

Uploaded by

Copyright:

Available Formats

1 s2.0 S0142061521009042 Main

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

1 s2.0 S0142061521009042 Main

Uploaded by

Copyright:

Available Formats

Electrical Power and Energy Systems 137 (2022) 107674

Contents lists available at ScienceDirect

International Journal of Electrical Power and Energy Systems

An MILP formulation for the optimum operation of AC microgrids with

DER Distributed Energy Resources time i

R Continuous variable or real number ith ESS

the ith PCC of the MG at time t Binary variables for OPF

Fig. 1. Basic scheme of the microgrid.

MINLP formulation for the optimal operation of AC microgrids with ( )

commercial optimization solvers. The remainder of this paper is orga­ ( )

tion 2. Subsequently, proposed MINLP model is linearized and the

Fig. 3. Feasible operation region of the DDG units in IS and GC mode.

variable in the proposed problem. This allows us to consider uv,int to be a

∀i ∈ ΩDDG ,∀t ∈ ΩT (19)

measured frequency ωm and voltage V m as expressed by:

2.2.3. Energy storage operational constraints ( ) ( ) ∑ PCC

where H can be replaced by G or B depending on (10) and (11); Cosine

bcurt curt curt

Approximation #1: This approximation technique deals with the

(quadrant) (name) (elements selected) #

DDG 1st/2nd LDDG

ESS full LESS

TR (PCC) full LPCC

No. Bus Smax CRDDG

DDGs 1 2 1500 11.105 PQ/2nd – – – 1 13 2000 11.148 PQ/2nd – – –

V /V* ** 1.01/1.05 σ fuel 0.007 $/Mcf [32] PFTR

SCEN 1 SCEN 2 SCEN 3 SCEN 1 SCEN 2 SCEN 3

Total CPCC ($/24 h) 1123.22 0.00 805.45 993.63 0.00 706.89

subject to: in Table 8 of Appendix B. The parameters of the adopted piecewise

NLP MILP Error NLP MILP Error

environment and solved using “Intlingprog” and “Fmincon” solvers,

4.1. Optimum operation of the microgrid

Fig. 13. Number of load shedding switching state cycle.

Fig. 14. Profile of load shedding at bus #19.

generated. As the DDG units are generating no reactive power in this

4.5. Evaluation of ESS with primary droop and secondary control

PLi,Ψ,t = (APLt /100) × PLi,Ψpeak , ∀i ∈ ΩBUS , Ψ ∈ F , ∀t ∈ NT (86)

QLi,Ψ,t = (APLt /100) × QLi,Ψpeak , ∀i ∈ ΩBUS , Ψ ∈ F , ∀t ∈ NT (87)

See Tables 7 and 8.

1 0 24 12 0.016 21 23 0.024 19 1 0.000 24 12 0.047 15 23 0.056 15

Fig. 17. Derivation of DDG feasible operation region.

(26) (a) 0 [0, R] – – – – –

(a) Linearization of Absolute Value Operator

(b) Linearization of the production between continuous and binary variables

1 Vi,Ψ,t ,Ψ∈F , ∀t ∈ ΩT R ΩBUS ✓ ✓ ✓ ✓

You might also like

commercial optimization solvers. The remainder of this paper is orga ( )