A Two-Level Simulation-Assisted Sequential Distribution System Restoration Model With Frequency Dynamics Constraints
11 Frequency Dynamics Constrained SequentialService Restoration in Inverter-DominatedMicrogrids
Qianzhi Zhang,
Student Member, IEEE,
Zixiao Ma,
Student Member, IEEE,
Yongli Zhu,
Member, IEEE, and Zhaoyu Wang,
Senior Member, IEEE.
Abstract —This paper proposes a service restoration modelfor unbalanced distribution systems and inverter-dominated mi-crogrids (MGs), in which frequency dynamics constraints aredeveloped to optimize the load pick-up during restoration processand guarantee frequency stability. After extreme events, thedamaged distribution systems can be sectionalized into severalisolated MGs to restore critical loads and tripped distributedgenerations (DGs) by black-start DGs. However, the high pene-tration of inverter-based DGs reduce the system inertia, whichresults in low-inertia issues and large frequency fluctuationduring the restoration. To address this challenge, we formulate asequential service restoration model as a mixed integer linearprogramming (MILP) problem. The proposed MILP modelexplicitly incorporates the frequency response into constraints,by interfacing with transient simulations of inverter-dominatedMGs. Numerical results on a modified IEEE 123-bus systemhave validated that the frequency dynamic performance of theproposed service restoration model are indeed improved.
Index Terms —Frequency dynamics, service restoration, net-work reconfiguration, inverter-dominated microgrids, simulation-based optimization.
I. I
NTRODUCTION E XTREME events can cause severe damages to powerdistribution systems [1], e.g. substation outage, line out-age, generator tripping, load loss, and consequently large-scale system blackouts [2]. During the network and servicerestoration, in order to isolate faults and restore critical loads,a distribution system can be sectionalized into several isolatedmicrogirds (MGs) [3]. Through the MG formation, buses,lines and loads in outage areas can be locally energized bydistributed generations (DGs), where more outage areas couldbe restored and the number of switching operations could beminimized [4]–[9]. In [4], the self-healing mode of MGs isconsidered to provide reliable power supply for critical loadsand restore the outage areas. In [5], a networked MGs-aidedapproach is developed for service restoration, which considersboth dispatchable and non-dispatchable DGs. In [6] and [7],the service restoration problem is formulated as a mixed inte-ger linear program (MILP) to maximize the critical loads to bepicked up while satisfying constraints for MG formation andremotely controlled devices. In [8], the formation of adaptive
This work was supported in part by the U.S. Department of Energy WindEnergy Technologies Office under Grant DE-EE00008956 (
Correspondingauthor: Zhaoyu Wang ).The authors are with the Department of Electrical and Computer Engineer-ing, Iowa State University, Ames, IA 50011 USA (e-mail: [email protected]). multiple MGs is developed as part of the critical servicerestoration strategy. In [9], a sequential service restorationframework is proposed to generate restoration solutions forMGs in the event of large-scale power outages. However, theprevious methods mainly use the conventional synchronousgenerators as the black-start units, and only consider steady-state constraints in the service restoration models, which havelimitations in the following aspects:(1)
An inverter-dominated MG can have low inertia:
Withthe increasing penetration of inverter-based DGs (IBDGs)in distribution systems, such as distributed wind and photo-voltaics (PVs) generations, the system inertia becomes lower[10], [11]. When sudden changes happen, such as DG outputchanging, load reconnecting, and line switching, the dynamicfrequency performance of such low-inertia distribution sys-tems can deteriorate [12]. This issue becomes even worsewhen restoring low-inertia inverter-dominated MGs. Withoutconsidering frequency dynamic constraints, the load pick-upplans may not be implemented in practice.(2)
Frequency dynamics need to be considered:
Previousstudies [13]–[16] have considered the impact of disturbanceson frequency dynamics in the service restoration problemusing different approaches. In [13], the amount of load restoredby DGs is limited by a fixed frequency response rate andmaximum allowable frequency deviation. However, becausethe frequency response rate is pre-determined in an off-linemanner, the impacts of significant load pick-up, topologychange, and load variations may not be fully captured by theoff-line model. In [14], the stability and security constraintsare incorporated into the restoration model. However, themodel has to be solved by meta-heuristic methods due tothe nonlinearity of the stability constraints, which may leadto large optimality gaps. In [15], the transient simulationresults of voltage and frequency are considered to evaluatethe potential MG restoration paths in an off-line manner. In[16], even though a control strategy of real-time frequencyregulation for network reconfiguration is developed, it is notco-optimized with the switch operations.(3)
Grid-forming DGs need to be considered:
In preciousstudies on optimal service restoration, DGs are usually mod-eled as grid-following sources (i.e., PQ sources). However,during the service restoration after a network blackout andloss of connection to the upstream feeder, a grid-formingIBDG will be needed to setup voltage and frequency referencesfor the blackout network [17]. The grid-following IBDGs a r X i v : . [ ee ss . S Y ] J a n are usually considered as PQ sources, which simply injectactive and reactive power based on the control commands.During outages, the grid-following IBDGs will be switchedoff. After outages, the grid-forming IBDGs have the black-startcapability, which can pick up loads after the faults are isolated.Unlike conventional synchronous generators, the IBDGs areconnected with power electronics converters and have norotating mass. Because there is no such P-f coupling in IBDGs,the concept of virtual synchronous generator (VSG) [18], [19]is usually used to emulate the inertia property in an IBDG.To alleviate the frequency fluctuations caused by load pick-up, we establish a MILP-based optimization model with fre-quency dynamics constraints for sequential service restorationto generate sequential actions for remotely controlled switches,restoration status for buses, lines, loads, operation actionsfor grid-forming and grid-following IBDGs, which interactswith the transient simulations of inverter-dominated MGs. Toincorporate the frequency dynamics constraints explicitly inthe optimization formulation, we associate the frequency nadirof the transient simulation with respect to the maximum loadthat a MG can pick up. The main contribution of this paperis two-folded: • A two-level simulation-assisted sequential servicerestoration model is established within a rolling horizonframework, which combines a MILP-based optimizationlevel of service restoration and a transient simulationlevel of inverter-dominated MGs. • Frequency dynamics constraints are developed and ex-plicitly incorporated in the optimization model, to asso-ciate the simulated frequency responses with the decisionvariables of maximum load step at each stage. Theseconstraints help restrict the system frequency drop duringthe transient periods of restoration. Thus, the generatedrestoration solution can be more secure and practical.The reminder of the paper is organized as follows: SectionII presents the overall framework of the proposed servicerestoration model. Section III introduces frequency dynamicsconstrained MILP-based sequential service restoration. SectionIV describes transient simulation of inverter-dominated MGs.Numerical results and conclusions are given in Section V andSection VI, respectively.II. O
VERVIEW OF THE P ROPOSED SERVICE RESTORATION M ODEL
The general framework of the proposed frequency dynamicsconstrained service restoration is shown in Fig. 1, includ-ing an optimization level of MILP-based sequential servicerestoration model and a transient simulation level of th-orderelectromagnetic inverter-dominated MG dynamic model. Afteroutages, the fault-affected areas of the distribution systemwill be isolated. Consequently, each isolated sub-network canbe considered as a MG [20], which can be formed by thevoltage/frequency support from the grid-forming IBDGs, andactive/reactive power supply from the grid-following IBDGs.In the proposed optimization level, each MG will determineits restoration solutions, including optimal service restorationstatus, optimal operation of remotely controlled switches and Time (s) F re qu e n c y ( H z ) t1 t2 t3 Optimization level: MILP sequential service restoration model Simulation level: dynamic model of inverter-dominated MG Restoration solutions
34 5 6
892 1 Stage III
134 5 6 Stage II
Stage I
Frequency nadir L oa d p i c k - up ( k W ) Time (s)
Dynamic model of inverter-dominated MG System frequency response
Fig. 1. The overall framework of the proposed service restoration model withoptimization level and simulation level optimal active/reactive power dispatches of IBDGs. To preventlarge frequency fluctuation due to a large load pick-up, themaximum restorable load for a given period is limited bythe proposed frequency dynamics constraints. In this way,the whole restoration process is divided into multiple stages.As shown in Fig. 1, the information exchanged between theoptimization level and the simulation level are the restorationsolution (by optimization) and MG system frequency nadirvalue (by transient simulation): at each restoration stage, theoptimization level will send the optimal restoration solutionto the simulation level; then, after receiving the restorationsolution, the simulation level will begin to run transientsimulation by the proposed dynamic model of each inverter-dominated MG, and send the frequency nadir value to theoptimization level for next restoration stage.To accurately reflect the dynamic frequency-supporting ca-pacities of grid-forming IBDGs during the service restorationprocess, a rolling-horizon framework is implemented in theproposed service restoration model. More specifically, we repeatedly run the MILP-based sequential service restorationmodel by incorporating the network configuration from thepreceding stage as the initial condition, and then feedbackthe frequency nadir value from the transient simulation tothe frequency dynamics constraints. For each stage: (1) thehorizon length will be fixed (2) then only the restorationsolution of first horizon of the current stage is retainedand transferred to the simulation level, while the remaininghorizons are discarded. (3) This process will keep going untilthe maximum restored load is reached in each MG .III. F
REQUENCY D YNAMICS C ONSTRAINED SERVICERESTORATION
This section presents the mathematical formulation forcoordinating remotely controlled switches, grid-forming and grid-following IBDGs, and the sequential restoration status ofbuses, lines and loads. Here, we consider a three-phase unbal-anced radial distribution system. The three-phase φ a , φ b , φ c are simplified as φ . Define the set Ω L = Ω SW L ∪ Ω NSW L ,where Ω SW L and Ω NSW L represent the set of switchable loadand the set of non-switchable loads, respectively. Define theset Ω G = Ω BS ∪ Ω NBS , where Ω BS and Ω NBS represent theset of grid-forming IBDGs with black start capability and theset of grid-following IBDGs without black start capability,respectively. Define the set Ω K = Ω SW K ∪ Ω NSW K , where Ω SW and Ω NSW represent the set of switchable lines and theset of non-switchable lines, respectively. Define Ω BK as theset of bus blocks, where bus block [9] is a group of busesinterconnected by non-switchable lines and those bus blocksare interconnected by switchable lines. It is assumed that eachbus block can be energized by grid-forming IBDGs. By forcingthe related binary variables of faulted lines to be zeros, eachfaulted area remains isolated during the restoration process. A. MILP-based sequential service restoration Formulation
The objective function (1) aims to maximize the totalrestored loads with priority factor w Li over a rolling horizon [ t, t + T ] as shown below: max (cid:88) t ∈ [ t,t + T ] (cid:88) i ∈ Ω L (cid:88) φ ∈ Ω φ ( w Li x Li,t P Li,φ,t ) (1)where P Li,φ,t and x Li,t are the restored load and restorationstatus of load at t . If the load demand P Li,φ,t is restored,then x Li,t = 1 . T is horizon length in the rolling horizonoptimization problem. In this work, the amount of restoredload is also bounded by frequency dynamics constraints withrespect to frequency response and maximum load step. Moredetails of frequency dynamics constraints are discussed inSection III-B.Constraints (2)-(11) are defined for the operation in eachformed MG during the service restoration process. Constraints(2) and (3) are the active and reactive nodal power balanceconstraints, where P K k,φ,t and Q K k,φ,t are the active and reactivepower flows along line k , and P G i,φ,t and Q G i,φ,t are the poweroutputs of the generators. Constraints (4) and (5) representthe active and reactive power limits of the lines, where thelimits ( P K , max k and Q K , max k ) are multiplied by the line statusbinary variable x Kk,t . Therefore, if a line is disconnected ordamaged x Kk,t = 0 , then constraints (4) and (5) are relaxed,which means that power cannot flow through the line. Twotypes of IBDGs are considered in the proposed model. Duringan outage, the grid-following IBDGs will be switched off.If grid-following IBDGs is connected with an energized bus,then it can be switched on. The grid-forming IBDGs providea reference voltage in the MG, which can energize the busand restore the part of the network that is not damaged ifthe fault is isolated. In constraints (6) and (7), the active andreactive power outputs of the grid-forming IBDGs are limitedby the maximum active and reactive capacities P G , max i and Q G , max i , respectively. Constraints (8) and (9) limit the activeand reactive outputs of the grid-following IBDGs. The grid-following IBDGs constraints are multiplied by binary variable x Gi,t . Consequently, if grid-following IBDG is switched off ortripped x Gi,t = 0 during the restoration, then constraints (8)and (9) will be relaxed. Constraints (10) and (11) calculatethe voltage difference along line k between bus i and bus j ,where U i,φ,t,s is the square of voltage magnitude of bus i . Weuse the big-M method [9] to relax constraints (10) and (11),if lines are damaged or disconnected x Kk,t = 0 . ˆ R k and ˆ X k are the unbalanced three-phase resistance matrix and reactancematrix of line k , which can be referred to [21]. The vector p k,φ represents the phases of line k . Constraint (12) guarantees thatthe voltage is limited within a specified region [ U min i , U max i ],and is set to 0 if the bus is in an outage area x Bi,t = 0 . (cid:88) k ∈ Ω K ( i,. ) P K k,φ,t − (cid:88) k ∈ Ω K ( .,i ) P K k,φ,t = P G i,φ,t − x Li,t P Li,φ,t , ∀ i, φ, t (2) (cid:88) k ∈ Ω K ( i,. ) Q K k,φ,t − (cid:88) k ∈ Ω K ( .,i ) Q K k,φ,t = Q G i,φ,t − x Li,t Q Li,φ,t , ∀ i, φ, t (3) − x Kk,t P K , max k ≤ P K k,φ,t ≤ x Kk,t P K , max k , ∀ k ∈ Ω K , φ, t (4) − x Kk,t Q K , max k ≤ Q K k,φ,t ≤ x Kk,t Q K , max k , ∀ k ∈ Ω K , φ, t (5) ≤ P G i,φ,t ≤ P G , max i , ∀ i ∈ Ω BS , φ, t (6) ≤ Q G i,φ,t ≤ Q G , max i , ∀ i ∈ Ω BS , φ, t (7) ≤ P G i,φ,t ≤ x Gi,t P G , max i , ∀ i ∈ Ω NBS , φ, t (8) ≤ Q G i,φ,t ≤ x Gi,t Q G , max i , ∀ i ∈ Ω NBS , φ, t (9) U i,φ,t − U j,φ,t ≥
2( ˆ R k P K k,φ,t + ˆ X k Q K k,φ,t )+ ( x Kk,t + p k,φ − M, ∀ k, ij ∈ Ω K , φ, t (10) U i,φ,t − U j,φ,t ≤
2( ˆ R k P K k,φ,t + ˆ X k Q K k,φ,t )+ (2 − x Kk,t − p k,φ ) M, ∀ k, ij ∈ Ω K , φ, t (11) x Bi,t U min i ≤ U i,φ,t ≤ x Bi,t U max i , ∀ i, φ, t (12)Constraints (13)-(17) ensure the physical connectionsamong buses, lines, IBDGs and loads during restorationprocess. In constraint (13), grid-following IBDGs will beswitched off x Gi,t = 0 if the connected bus is not energized x Bi,t = 0 . Constraint (14) implies a switchable line can onlybe energized when both end buses are energized. Constraint(15) presents that a non-switchable line can be energized onceone of two end buses is energized. Constraint (16) ensuresthat a switchable load can be energized if the connected busis energized. Constraint (17) allows that a non-switchableload can be immediately energized once the connected busis energized. x Gi,t ≤ x Bi,t , ∀ i ∈ Ω NBS , t (13) x Kk,t ≤ x Bi,t , x
Kk,t ≤ x Bj,t , ∀ k, ij ∈ Ω SW K , t (14) x Kk,t = x Bi,t , x
Kk,t = x Bj,t , ∀ ( i, k ) ∈ Ω NSW K , t (15) x Li,t ≤ x Bi,t , i ∈ Ω SW L , t (16) x Li,t = x Bi,t , ∀ i ∈ Ω NSW L , t (17)Constraints (18)-(20) ensure that each formed MG remainsisolated from each other and each MG can maintain a treetopology during the restoration process. Constraint (18) im-plies that if one bus i is located in one bus block, i ∈ Ω BK ,then the energization status of bus and the corresponding busblock keep the same. Here x BKB,t represents the energizationstatus of bus block BK . To avoid forming loop topology,constraint (19) guarantees that a switchable line cannot beclosed at time t if its both end bus blocks are already energizedat previous time t − . If one bus block is not energized atprevious time t − , then constraint (20) makes sure that thisbus block can only be energized at time t by at most oneof the connected switchable lines. Constraints (21) and (22)make sure that each formed MG has a reasonable restorationand energization sequence of switchable lines and bus blocks.Constraints (21) implies that energized switchable lines canenergize the connected bus block. Constraints (22) requiresthat a switchable line can only be energized at time t , if atleast one of the connected bus block is energized at previoustime t − . x Bi,t = x BKi,t , ∀ i ∈ Ω BK , t (18) ( x BKi,t − x BKi,t − ) + ( x BKj,t − x BKj,t − ) ≥ x Kk,t − x Kk,t − , ∀ k, ij ∈ Ω SW K , t ≥ (19) (cid:88) ki,k ∈ Ω i ( x Kki,t − x Kki,t − ) + (cid:88) ij,j ∈ Ω i ( x Kij,t − x Kij,t − ) ≤ x BKi,t − M, ∀ k, ij ∈ Ω SW K , t ≥ (20) x BKi,t − ≤ (cid:88) ki,k ∈ Ω i ( x Kki,t ) + (cid:88) ij,j ∈ Ω i ( x Kij,t ) , ∀ k, ij ∈ Ω SW K , t ≥ (21) x Kij,t ≤ x BKi,t − + x BKj,t − , ∀ ij ∈ Ω SW K , t ≥ (22) B. Simulation-based Frequency Dynamics Constraints
By considering the dynamics of each isolated inverter-dominated MG during the transitions of network reconfigu-ration and service restoration, constraints (23) and (27) havebeen added here to avoid the potential large frequency devi-ations caused by MG formation and oversized load pick-up.The variable of maximum load step P G,MLSi,t has been appliedin constraint (23) to ensure that the restored load is limited byan upper bound for each restoration stage, as follows: ≤ P G,MLSi,t ≤ P G,MLSi,t − + α (∆ f max − ∆ f measure ) , i ∈ Ω BS , t ≥ (23)where ∆ f max is the user-defined maximum allowable fre-quency drop limit and α is hyper-parameter to be deter-mined later. ∆ f measure is the measured maximum transientfrequency drop from the results of simulation level. Thecoefficient α is related to the virtual inertia characteristics ofthe VSG. This can be shown by the following manipulation: α (∆ f max − ∆ f measure ) = α ∆ t (∆ f max − ∆ f measure )∆ t = α ∆ t ( f − f min − ( f − f nadir ))∆ t = α ∆ t ( f nadir − f min )∆ t ≈ α ∆ t ∆ f ∆ t (24)where ∆ f is transient variation of frequency and ∆ t is thetransient variation time of frequency during the transient sim-ulation. f is the nominal steady-state frequency, e.g. 60Hz. f nadir is the lowest frequency reached during the transientsimulation. f min is the minimum allowable frequency. In thispaper, ∆ f max is set to 0.3 Hz based on [22]. On the otherhand, the active power output of a VSG unit is described as: P V SG = K I d ∆ wdt + K p ∆ w (25)where the inertia emulating characteristic K I [19] can becalculated as below: K I = 2 HP g w (26)In (26), P g is the nominal apparent power of the genera-tor, H is the inertia constant of the mimicked synchronousmachine (unit: seconds), and w is the nominal frequencyof the grid. After per-unit conversion, it is easy to see that α ∆ t = 2 H . Let the estimated ∆ t be 6 seconds and H be0.3 seconds, then α = 0 . . Finally, constraint (27) ensures therestored load and frequency response of the IBDGs do notexceed the user-defined thresholds. − x Gi,t P G,MLSi,t ≤ P G i,φ,t − P G i,φ,t − ≤ x Gi,t P G,MLSi,t , i ∈ Ω BS , φ, t ≥ (27)Note that now P G,MLSi,t is not a constant as in previ-ous literature, but a decision variable who is dynamicallydetermined by (27) during the optimization combining withtransient simulation information.IV. T
RANSIENT S IMULATION OF I NVERTER -D OMINATED
MG F
ORMATION
In optimization level, our target is to pick up the mostload while satisfying a series of constraints. One of theseconstraints should be frequency dynamics constraint which isderived from simulation level. However, due to the differenttime scales and nonlinearity, the conventional dynamic securityconstraints cannot be directly solved in optimization problem,such as Lyapunov theory, LaSalle’s theorem and so on. There-fore, we need a connection variable between the two levels.For this purpose, we assume that the changes of typologiesbetween each two sequential stages can be represented by thechange of picked-up loads P L . The sudden load change of P L results in a disturbance in MGs in the time-scale of simulationlevel. During the transience to the new equilibrium (operationpoint), the system states such as frequency will deviate fromtheir nominal values. Therefore, it is natural to estimate thedynamic security margin with the allowed maximum range ofdeviations.Since the frequency of each inverter-dominated MG ismainly controlled by the grid-forming IBDGs, we can approx-imate the maximum frequency deviation during the transience Coordinate transformation Power calculation Low pass filter Droop controller − − − − + + + + 𝑽 𝑰 𝑽 𝒅 𝑽 𝒒 𝑰 𝒅 𝑰 𝒒 𝑽 𝒅 𝑰 𝒅 + 𝑽 𝒒 𝑰 𝒒 𝑽 𝒅 𝑰 𝒅 - 𝑽 𝒒 𝑰 𝒒 𝑷 𝒎 𝑸 𝒎 𝑸 𝑳 𝑷 𝑳 𝝎 𝒄 𝒔 + 𝝎 𝒄 𝝎 𝒄 𝒔 + 𝝎 𝒄 𝑸 𝑷 𝑫 𝑷 𝑫 𝑸 𝑽 𝜽 𝝎 𝝎 𝒔𝒆𝒕 𝑽 𝒔𝒆𝒕 𝒂𝒃𝒄 𝒅𝒒 Fig. 2. Diagram of studied MG control system. by observing the dynamic response of the grid-forming IBDGsunder sudden load change. Usually, the standard outer droopcontrol together with inner double-loop control structure isadopted for each IBDGs unit. The inner double-loop controlis designed as two PI control loops governing output voltageand frequency of the inverter, respectively. The voltage andfrequency reference signals are provided by the outer droopcontroller. It is worth noting that, the inner control loop isalways designed to operate on much faster time-scale forrapid response, thus allowing us to tuning the inner and outercontrollers separately. By taking advantage of this two-time-scale property inverter can be effectively represented by onlyusing the terminal states and line states of the inverter, suchas terminal angle, frequency, voltage and line current [23].The overall simplified control structure is shown as Fig. 2.The three-phase output voltage V ,abc and current I ,abc aremeasured from the terminal bus of the inverter and transformedinto dq axis firstly. Then, the active and reactive power P and Q are obtained by filtering the calculated power measurements P m and Q m with cut-off frequency ω c as below: ˙ P = ω c ( P m − P ) , (28) ˙ Q = ω c ( Q m − Q ) . (29)In the next step, to integrate in the sudden load change, thedroop characteristic is designed as: ω = ω set − D P ( P − P L ) , (30) V = V set − D Q ( Q − Q L ) . (31)where ω set and V set are the set points of frequency andvoltage controllers, respectively; D P and D Q are P − ω and Q − V droop gains respectively; P L and Q L are the restoredactive and reactive loads, respectively. Finally, a th-orderelectromagnetic MG model augmented by low-pass filter canbe derived and concluded as follows: ˙ P = ω c ( V cos θI d + V sin θI q − P ) , (32) ˙ Q = ω c ( V sin θI d − V cos θI q − Q ) , (33) ˙ θ = ω − ω , (34) ˙ ω = ω c ( ω set − ω + D P ( P − P L )) , (35) ˙ V = ω c ( V set − V + D Q ( Q − Q L )) , (36) ˙ I d = ( V cos θ − V bus − RI d ) /L + ω o I q , (37) ˙ I q = ( V sin θ − RI q ) /L − ω o I d , (38) where θ is phase angle; ω is angular frequency in rad/s ; ω is a fixed angular frequency; I d and I q are dq -axis currents; R and L are aggregate inductance and resistance of connectionsfrom the inverter terminal’s point view, respectively.After the process of fault detection [24] and sub-gridsisolation are finished, the proposed service restoration modelwill begin to work. Each isolated network will begin toform a MG depending on the location of the nearest grid-forming IBDG. The interaction between the proposed transientsimulation and the established optimization problem of servicerestoration is described as follows: (a) Solving the optimal service restoration problem : Givenhorizon length T in each restoration stage, the MILP-basedsequential service restoration problem (1)–(23) and (27) issolved, and the restoration solution is obtained for each formedMG. (b) Transient simulation of inverter-dominated MGs : Afterreceiving restoration solutions of current stage from opti-mization level, the dynamic frequency response is simulatedby (32)–(38) and the frequency nadir is calculated for eachinverter-dominated MG. (c)
Check the progress of service restoration : Go back to(a) to generate the restoration solution with newly obtainedfrequency responses of all MGs for next restoration stage, untilthe maximum service restoration progress is made, e.g. all thecritical loads are energized.V. N
UMERICAL R ESULTS
A modified IEEE 123-bus test system [25] in Fig. 3 is usedto test the performance of the proposed frequency dynam-ics constrained service restoration model. The modified testsystem has equipped with remotely controlled switches, grid-following and grid-forming IBDGs. There are four line faultson lines 1, 19, 54 and 70 are detected, which are shown as thered dashed lines. They are assumed to be persisting during therestoration process until the faulty areas are cleared to maintainthe radial topology and isolate the faulty areas. Consequently,four MGs can be formed for service restoration with grid-forming IBDGs and switches. We demonstrate the effective-ness of our proposed service restoration model through numer-ical evaluations on the following experiments : (i) Comparisonbetween a base case (i.e. without the proposed frequencydynamics constraints) and the case with those constraints. (ii)Cases with the proposed frequency dynamics constraints underdifferent α values (iii) Cases with the proposed frequencydynamics constraints under different values of the horizonlength T . All the case studies are implemented using a PC withIntel Core i7-4790 3.6 GHz CPU and 16 GB RAM hardware.The simulations are performed in MATLAB R2019b, whichintegrates YALMIP Toolbox with IBM ILOG CPLEX 12.9solver for optimization. A. Sequential Service Restoration Results
As shown in (23), the relationship between the maximumload step and the frequency nadir is influenced by the valueof hyper-parameter α in the frequency-dynamics constraints.Therefore, different α values may lead to different service
67 8
25 24
93 9457
54 55
61 60 59 6263 64 873638 37
39 40
98 99 DG Three-phase Single-phase Grid-forming DG
Tap Changer DG DG Two-phase T1 T2 T3 DG DG DG DG DG Grid-following DG
SW1 SW2
Switch
SW3
SW4
SW6
SW7
SW8
SW9
SW10 DG SW11SW12
SW13SW17
SW14
SW15 SW16
Fault DG DG DG SW5
DGDG
DG DG
DG DG DG DG DGDG DG DG DG DG DG DG DG Fig. 3. Modified IEEE 123 node test feeder.
13 2 12 9 14151617 18 828321 20 19222325 24
31 29
34 35
49 48 50
52 53 1198180
76 74
75 7372 91
86 9071
93 94
57 5854
55 5661
60 59 6263 64 873638
37 39 40
41 44
106 107
97 98 99113 DG Three-phase Single-phase Grid-forming DG
DG DG
Two-phase T2 DGDG DGDG DG
Grid-following DG
SW1(1) SW2 (4)
Switch
SW3 (1) SW4 (1)
SW6 (2) SW7 (3)
SW8 (4) SW9 (1)SW10 (2) DG SW11 (2)
SW1
SW13 (4)SW17 (2) SW14 (4)SW15 (5)
SW16 (2)
MG1
MG2 MG4MG3
DG DGDG
SW5 (2)
DGDGDG DG DG DG DG DG DGDGDG DG DG DGDG DG DG DG
Fig. 4. Restoration solutions for the formed MG1-MG4. restoration results. In this case, the horizon length T and thehyper-parameter α are set to 4 and 0.1, respectively.As shown in Fig. 4, the system is partitioned into fourMGs by energizing the switchable lines sequentially, and theradial structure of each MG is maintained at each stage. Insideeach formed MG, the power balance is achieved betweenthe restored load and power outputs of IBDGs. The valuenearby each switch in Fig. 4 represents the stage-id when itcloses. It can be observed that MG2 only needs 3 stages to befully restored, while MG1 and MG3 can restore in 4 stages.However, due to the heavy loading situation, MG4 is graduallyrestored in 5 stages to ensure a relatively smooth frequencydynamics.For each restoration stage, the restored loads and frequency TABLE IR
ESTORED L OADS , F
REQUENCY N ADIR AND C OMPUTATION T IME FOR
MG1-MG4Cases Restored load(kW) Frequency nadir(Hz)MG1( T = 4 and α = 0 . ) Stage 1 280.5 59.7044Stage 2 280.5 59.9992Stage 3 280.5 59.9992Stage 4 346.5 59.9200Stage 5 346.5 59.9989MG2( T = 4 and α = 0 . ) Stage 1 230.0 59.7079Stage 2 360.0 59.8201Stage 3 420.0 59.9146Stage 4 420.0 59.9984Stage 5 420.0 59.9984MG3( T = 4 and α = 0 . ) Stage 1 212.5 59.7116Stage 2 212.5 59.9990Stage 3 212.5 59.9990Stage 4 382.5 59.7656Stage 5 382.5 59.9985MG4( T = 4 and α = 0 . ) Stage 1 192.0 59.7910Stage 2 324.0 59.8541Stage 3 414.0 59.9003Stage 4 570.0 59.8230Stage 5 624.0 59.9364 nadir in MG1-MG4 are shown in Table I. Total 1773 kWof load are restored at the end of the 5 stages. It can beobserved the service restoration actions happened in certainstages rather than in all stages. For example, MG1 restores280.5 kW of load in Stage 1, but it restores no more load untilStage 4. While MG4 takes action on service restoration in eachstage. It is because the sequential service restoration is limitedby operational constraints, among which the maximum loadstep in each stage is again limited by the proposed frequency-dynamics constraints.The comparison of total restored loads with and withoutconsidering the proposed frequency dynamics constraints isshown in Fig. 5. It can be observed that the proposed servicerestoration model restores no more loads between Stage 5 andStage 6, i.e. it can restore all the loads in the first 5 stages.While the base case needs 6 stages to be fully restored. In theearly stages 1 to 3, the restored load of the proposed modelis a little bit less than the base case. It is because duringthe early restoration stages, the proposed service restorationmodel generated a restoration solution that prevents too lowfrequency nadir during transients. The base case restores moreloads at Stage 1 to Stage 3 without considering the limitationon the frequency nadir. However, Stage 4 is a turning pointwhen the proposed model restores more loads than the basecase. Fig. 6 shows the frequency responses of MG4 in Stage 1with and without the frequency dynamics constraints. By thiscomparison, it can be observed that both the rate of change offrequency and frequency nadir are significantly improved.Fig. 7 shows the frequency dynamics of each inverter-dominated MG based on the proposed restoration model. Theresults show that the MG frequency drops when the load is Stage R es o r e d l o a d ( k W ) w/ frequency dynanmics constraintsw/o frequency dynanmics constraints Fig. 5. Total restored load with and without considering frequency dynamicsconstraints.Fig. 6. Frequency responses of MG4 in Stage 1 with and without frequencydynamics constraints. restored. Because the maximum load step is constrained in theproposed MILP-based sequential service restoration model, thefrequency nadir is also constrained. When load is picked upas the frequency drops, the frequency nadir can be effectivelymaintained above the f min threshold. B. Impact of α in Frequency Dynamics Constraints Compared to other MGs, MG4 is heavily loaded with thelargest number of nodes. Based on the results of Fig. 4, MG4needs more stages to be fully restored compared to other MGs.Therefore, MG4 is chosen to test the effect of different α values. The dynamic frequency responses of MG4 with α =0 . , α = 0 . and α = 1 . are shown in Fig. 8. When α = 0 . ,5 stages are required to fully restore all the loads; while only4 restoration stages are needed when α = 0 . or α = 1 . . It t (0.1s) f ( H z ) MG t (0.1s) f ( H z ) MG t (0.1s) f ( H z ) MG t (0.1s) f ( H z ) MG Fig. 7. Dynamic frequency responses of inverter-dominated MG1-MG4. can be also observed that the frequency nadir with α = 0 . is higher than the frequency nadir with α = 0 . or α = 1 . .Hence, there can be a trade-off between dynamic frequencyperformance and restoration performance regarding the choiceof α : too small α may lead to too slow restoration; in turn,too large α may deteriorate the dynamic performance of thesystem frequency in a realistic restoration. C. Impact of Time Horizon
This section shows that different values of the horizonlength T may cause different service restoration results. TableII summarizes the total restored loads and computation timeusing different horizon lengths in the proposed service restora-tion model. On the one side, the restored loads of case with T = 2 and T = 3 are less than that of the cases with T ≥ ,where the total restored load can reach the maximum level.Therefore, the results with small number of horizon length T = 2 and T = 3 are sub-optimal restoration solutions.On the other side, the longer horizon length also leads toheavy computation burden and increase the computation time.Similar to the impact of α , there can be a trade-off between the Fig. 8. Dynamic frequency responses of MG4 with different α .TABLE IIR ESTORED L OADS , F
REQUENCY N ADIR AND C OMPUTATION T IME WITH D IFFERENT H ORIZON L ENGTHS IN
MG4Total restored load (kW) Computation time (s) T = 2 T = 3 T = 4 T = 5 T = 6 computation time and the quality of solution when determiningthe value of T . VI. C ONCLUSION
To improve the dynamic performance of the system fre-quency during service restoration of a unbalanced distributionsystems in an inverter-dominated environment, we propose asimulation-assisted optimization model considering frequencydynamics constraints with clear physical meanings. Resultsdemonstrate that: (i) The proposed frequency dynamics con-strained service restoration model can significantly reduce thetransient frequency drop during MGs forming and servicerestoration. (ii) Other steady-state performance indicators ofour proposed method can rival that of the conventional meth-ods, in terms of the final restored total load and the requirednumber of restoration stages. Investigating on how to choosethe best hyper-parameter α and horizon length T will be thenext research direction. R EFERENCES[1] E. O. of the President, “Economic benefits of increasing electric gridresilience to weather outages,” White House, Tech. Rep., 2020.[2] A. M. Salman, Y. Li, and M. G. Stewart, “Evaluating system reliabilityand targeted hardening strategies of power distribution systems subjectedto hurricanes,”
Reliab. Eng. Syst. Saf. , vol. 144, pp. 319–333, Dec. 2015. [3] H. Haggi, R. R. nejad, M. Song, and W. Sun, “A review of smart gridrestoration to enhance cyber-physical system resilience,” in , 2019, pp. 4008–4013.[4] Z. Wang and J. Wang, “Self-healing resilient distribution systems basedon sectionalization into microgrids,”
IEEE Trans. Power Syst. , vol. 30,pp. 3139–3149, Nov. 2015.[5] A. Arif and Z. Wang, “Networked microgrids for service restorationin resilient distribution systems,”
IET Gener. Transm. Distrib. , vol. 11,no. 14, pp. 3612–3619, Aug. 2017.[6] C. Chen, J. Wang, F. Qiu, and D. Zhao, “Resilient distribution system bymicrogrids formation after natural disasters,”
IEEE Trans. Smart Grid ,vol. 7, no. 2, pp. 958–966, Mar. 2016.[7] S. Yao, P. Wang, and T. Zhao, “Transportable energy storage for moreresilient distribution systems with multiple microgrids,”
IEEE Trans. onSmart Grid , vol. 10, pp. 3331–3341, May 2019.[8] L. Che and M. Shahidehpour, “Adaptive formation of microgrids withmobile emergency resources for critical service restoration in extremeconditions,”
IEEE Trans. Power Syst. , vol. 34, no. 1, pp. 742–753, Jan.2019.[9] B. Chen, C. Chen, J. Wang, and K. L. Butler-Purry, “Sequential servicerestoration for unbalanced distribution systems and microgrids,”
IEEETrans. Power Syst. , vol. 33, pp. 1507–1520, Mar. 2018.[10] Y. Wen, W. Li, G. Huang, and X. Liu, “Frequency dynamics constrainedunit commitment with battery energy storage,”
IEEE Trans. Power Syst. ,vol. 31, no. 6, pp. 5115–5125, Nov. 2016.[11] H. Gu, R. Yan, T. K. Saha, E. Muljadi, J. Tan, and Y. Zhang, “Zonalinertia constrained generator dispatch considering load frequency relief,”
IEEE Trans. Power Syst. , vol. 35, no. 4, pp. 3065–3077, Jul. 2020.[12] Y. Wen, C. Y. Chung, X. Liu, and L. Che, “Microgrid dispatchwith frequency-aware islanding constraints,”
IEEE Trans. Power Syst. ,vol. 34, no. 3, pp. 2465–2468, May 2019.[13] O. Bassey, K. L. Butler-Purry, and B. Chen, “Dynamic modeling ofsequential service restoration in islanded single master microgrids,”
IEEE Trans. Power Syst. , vol. 35, no. 1, pp. 202–214, Jan. 2020.[14] B. Qin, H. Gao, J. Ma, W. Li, and A. Y. Zomaya, “An input-to-statestability-based load restoration approach for isolated power systems,”
Energies , vol. 11, pp. 597–614, Mar. 2018.[15] Y. Xu, C. Liu, K. P. Schneider, F. K. Tuffner, and D. T. Ton, “Microgridsfor service restoration to critical load in a resilient distribution system,”
IEEE Trans. Smart Grid , vol. 9, pp. 426–437, Jan. 2018.[16] Y. Du, X. Lu, J. Wang, and S. Lukic, “Distributed secondary controlstrategy for microgrid operation with dynamic boundaries,”
IEEE Trans.Smart Grid , vol. 10, no. 5, pp. 5269–5285, Sept. 2019.[17] B. K. Poolla, D. Grob, and F. Dorfler, “Placement and implementationof grid-forming and grid-following virtual inertia and fast frequencyresponse,”
IEEE Trans. Power Syst. , vol. 34, pp. 3035–3046, Jul. 2019.[18] K. Y. Yap, C. R. Sarimuthu, and J. M.-Y. Lim, “Virtual inertia-basedinverters for mitigating frequency instability in grid-connected renewableenergy system: A review,”
Appl. Sci. , vol. 9, no. 24, p. 5300, Dec. 2019.[19] H. Bevrani, T. Ise, and Y. Miura, “Virtual synchronous generators: Asurvey and new perspectives,”
Int. J. Electr. Power Energy Syst. , vol. 54,pp. 244–254, Jan. 2014.[20] Y. Kim, J. Wang, and X. Lu, “A framework for load service restorationusing dynamic change in boundaries of advanced microgrids withsynchronous-machine dgs,”
IEEE Trans. Smart Grid , vol. 9, no. 4, pp.3676–3690, Jul. 2018.[21] Q. Zhang, K. Dehghanpour, and Z. Wang, “Distributed CVR in unbal-anced distribution systems with PV penetration,”
IEEE Trans. SmartGrid , vol. 10, no. 5, pp. 5308–5319, Sept. 2019.[22] U. L. S. R. Group, “Underfrequency load shedding assessment reportwestern interconnection,” Western Electricity Coordinating Council,Tech. Rep., 2013.[23] P. Vorobev, P. Huang, M. A. Hosani, J. L. Kirtley, and K. Turitsyn,“High-fidelity model order reduction for microgrids stability assess-ment,”
IEEE Trans. Power Syst. , vol. 33, pp. 874–887, Jan. 2018.[24] Y. Yuan, K. Dehghanpour, F. Bu, and Z. Wang, “Outage detectionin partially observable distribution systems using smart meters andgenerative adversarial networks,”