MPC-Based Black Start and Restoration for Resilient DER-Rich Electric Distribution System

For past several years the resiliency of the power grid is severely challenged by extreme events. Black start and restoration scheme (BS&RS) is critical to enhance distribution grid resiliency. Variability in the power supply capacity of distributed energy resources (DERs) and load demand stand against the efficiency of sequential black start restoration of an unbalanced distribution system. An efficient black start restoration scheme must consider forecasted generation and demand profile to assess resource participation and generate valid switching sequence during sequential restoration of distribution grid. This work proposes a novel model predictive control (MPC) based efficient BS&RS strategy for reconfigurable distribution grid with DERs. The proposed algorithm generates optimal switching sequence by coordinating smart switches, black start DERs (BS-DERs) and energy storage systems (ESSs). Active power-frequency droop, reactive power-voltage droop and ramp rate constraints are assigned to BS-DERs for practical purposes. The non-convexity of sequential service restoration problem due to operational constraints of distribution grid, line switches and dynamics of ESS are linearized and solved as a mixed integer linear programming problem (MILP). The proposed scheme is also extended to incorporate errors associated with forecasting, DER capacity and load demand. Results of IEEE 123-bus and 1069-bus systems indicate effectiveness of the proposed approach.


Indices m, n
Index of buses. mn Index of lines. Parameters α h , p 1 , p 2 Weighting factors. P Der n Forecasted active power capacity of DERs at bus n (p.u.). P n /Q n Forecasted active/ reactive power demand at bus n (p.u.). P max mn /Q max mn /S max mn Maximum active/ reactive/ apparent power capacity of line mn (p.u.).
The associate editor coordinating the review of this manuscript and approving it for publication was Cuo Zhang .

R n
Ramp rate of BS-DER connected to n th bus. r P mn (x P mn ) Resistance (Reactance) of line mn (p.u.) between phases P mn .

I. INTRODUCTION
Distribution grid is going through transitions due to increasing integration of inverter connected distributed energy resources (DERs). Large-scale DERs are replacing fossil fuel based plants to decarbonize electricity generation. However these fossil fuel plants have generators with sufficient inertia to retain dynamic stability and robustness of power system. Increase in the percentage of low inertia renewable generators make the power system prone to larger disturbances which may even result in blackouts. Typical reasons behind blackouts are weather related outages, equipment failures, damages caused by vegetation, accidents and vandalism. Efficient and fast acting black start services can enhance system resiliency to blackouts. Therefore, transmission grids within North American Electricity Reliability Corporation (NERC) require to be equipped with black start generators and services [1]. Black start generators such as hydro, thermal, gas turbine and nuclear generators in transmission system implicitly strengthen distribution system resiliency [2]. A grid forming DER operates as BS-DER only when it has sufficient generation to support the restoration process and it can supply power to a power grid without any external power supply [2]. These generators are also capable of regulating voltage and frequency of the restored power grid within acceptable limits [3], [4]. The increase in the percentage of non conventional sources in the power grid calls for extensive research in black start capabilities of DERs for restoration of distribution islands of an isolated system.
Black start restoration is achieved in a sequential manner to maintain adequate reserve till the non black start DERs or generators are restored. Also, with the increase in DERs in the power grid, restoration algorithms need to consider inter-temporal control and operational constraints to identify a feasible restoration solution [5]. These inter-temporal and operational constraints include ramp rate constraints of dispatchable DERs, dynamics of ESS / battery state of charge (SOC) and need for adjustment of the protection system settings with changes in the power grid topology. Sequential restoration schemes can incorporate all these constraints in identifying optimum switching sequence.
Predominantly, black start restoration of a distribution grid is a ''bottom-up'' restoration scheme where BS-DERs form isolated islands to restore loads [6]. In this sequential restoration scheme, the boundaries of energized islands increase, reducing the size of de-energized parts of the distribution grid. This process is repeated until the distribution grid is restored completely. These islands are interconnected further to form larger networks. Depending on the size of a distribution grid, the restoration process can take hours to days. To accelerate the black start restoration process, multiple BS-DERs located at different parts of the distribution grid coordinate and operate parallelly [7].
Distribution grid restoration problem is solved using genetic algorithm (GA), particle swarm optimization (PSO) [8], fuzzy multi-objective [9], branch-and bound (B&B) [10], etc. Though these heuristic approaches are widely used in solving mixed integer problems, optimality of the solution is not guaranteed. With the introduction of convex relaxation, restoration problem can be formulated as mixed integer linear programming (MILP), mixed integer quadratic programming (MIQP) and mixed integer second order cone programming problem (MISOCP), etc [11]. MILP guarantees optimality in a well defined problem. Therefore, a MILP based black start restoration algorithm is proposed in this work.
For systematic black start restoration, it is also imperative to include irregularities associated with future generation from DERs and load demand of the distribution grid. Due to the ability of model predictive control (MPC) to take into account the above factors [12], MPC based restoration schemes outperform conventional restoration strategies [13]. MPC has been applied to power system for economic dispatch in the presence of DERs, energy management system of smart houses, microgrid energy management [14], etc. In [15], one step MPC based restoration algorithm is proposed. Another one step restoration algorithm reported in [16] incorporates MPC to minimize forecast errors or sudden fluctuations in the input variables. One step MPC based restoration algorithms suffer from scalability issues due to increase in the number of binary decision variables with the increase in the size of the system. Therefore, large scale one step MPC based restoration algorithms require to compute possible restoration paths offline [15]. To the best of our knowledge, no formulation exists in literature for MPC based sequential restoration scheme. Hence, a finite horizon MPC based sequential black start and restoration algorithm is proposed in this work.
A transmission grid black start restoration scheme is proposed in [6] but the literature on black start restoration scheme for distribution grid is limited. A sequential bus block restoration scheme is proposed in [17]. In this approach, distribution grid is divided into multiple bus blocks at the beginning of the sequential restoration scheme. The algorithm generates restoration sequence with single black start DER in the distribution grid. A similar approach is proposed with multiple black start DERs in [5], but the efficiency of restoration scheme to incorporate variability in generation and demand or maintain the frequency of the system within allowable limits are not discussed. The above formulation is extended to capture frequency dynamics of a system in [18]. Another efficient MILP based sequential black start restoration algorithm is proposed in [19]. The algorithm utilizes commodity flow model to embody multiple black start DERs in generating switching sequence. However, it does not consider unbalanced grid conditions, variability in load demand, generation capacity and dynamics of ESS. The number of lines repaired or energized at an interval needs to be specified in the above algorithm. This is difficult to achieve as the number of lines energized by BS-DERs is determined by several different factors such as ramp rate of generators, generation and demand profile of that interval etc. The augmented presence of DERs make it important to consider the effect of irregularities of the black start DER generation in distribution grid restoration and develop necessary control mechanism. As MPC is dependent on forecasted generation and demand profile, therefore it is also crucial to include forecasting errors in the formulation. Most of the above approaches do not include uncertainty or forecasting error in their study. This paper is primarily focused on finding restoration sequence for operator's decision support which requires deterministic solution. Consequently in this work, a worst case deterministic MPC based sequential restoration scheme is proposed. Similar approaches can be found in [20] and [21]. To bridge the existing research gap, key contributions of this work are summarized below: 1) Proposed a parallel restoration framework with gridfarming, grid following and black start Distributed Energy Resources (DERs) in Power Grids using Model Predictive Control (MPC) and Mixed-Integer Linear Programming (MILP). The framework utilizes a centralized leader-follower architecture, where Grid-Forming DERs or blackstart-DERs (BS-DERs) act as leaders, and Grid-Following DERs act as followers [22], [23], [24]. The MPC-based receding horizon restoration scheme autonomously terminates upon full line energization and load restoration. 2) In this study, a novel system model is introduced, encompassing the integration of three-phase power flow equations, necessary conditions for radiality in distribution grids, and a comprehensive set of constraints specifically tailored for Distributed Energy Resources (DERs). The proposed model incorporates practical considerations, including DER ramp rate constraints, active power-frequency constraints, and reactive power-voltage droop constraints, which are critical for accurately representing the operational constraints. Additionally, a Quasi-Static Charging and Discharging Mechanism is incorporated into the model, enabling an accurate representation of the behavior of Energy Storage Systems (ESSs). 3) In order to account for the inherent uncertainty associated with DERs capacities to supply power and load demand forecasts, an extension is proposed for the deterministic worst-case Model Predictive Control (MPC) based sequential restoration scheme for robust restoration.
The rest of the paper is structured as follows: In Section II we introduced a centralized black start restoration problem. Several case studies, discussions and comparisons are provided in Section III. The concluding remarks of the study are in Section IV.

II. MATHEMATICAL FORMULATION
To simplify the BS&RS formulation, following assumptions are considered: 1) This study is applicable for all different types of BS-DERs such as solar, wind, hydro, biomass, etc. These non dispatchable DERs are assumed to be with ESSs to operate as dispatchable DERs. Detailed representation of battery integrated DER will be future scope of this work.
2) There is no actual demand in a grid during complete blackouts. Therefore, forecasted load demand profile is utilized in the restoration scheme to restore the power grid. A fractional load restoration scheme is employed, which is feasible with smart meters and demand response scheme.
3) The proposed BS&RS can be used for both complete and partial restoration of the distribution grid. Accordingly, droop controllers are introduced to keep the frequency and voltage levels of the system within allowable limits while generating switching sequence.
In the proposed MPC based BS&RS formulation, forecasted active load demand (P n ) is used for the prediction of other decision variables based on the mathematical model. VOLUME 11, 2023 Generation forecast is also included in the study to introduce variability in the power supply capacity (P Der n ) of BS-DERs. The manipulated variables derived through the online optimization are output power of BS-DERs, status of lines and output of ESSs.
The mathematical model for the proposed MPC base BS&RS is discussed in this section. A schematic of a distribution grid is shown in Fig. 1. Presence of different types of DERs in a distribution grid are also marked in the figure. The outage of substation creates complete blackout in the distribution grid. BS-DER activated at bus n initiates the restoration process. The distribution grid is restored by sequentially closing line switches β kl , β op and β pq . Lets assume at time t = 1, line switch β kl is closed allowing power supply from BS-DER to loads connected to buses l and q. Restoration of bus q enables grid following DER connected to bus q to take part in the restoration process by supplying power. Similarly all other open switches of the distribution grid are closed in the successive time intervals.

A. OBJECTIVE FUNCTION
The multi-objective optimization problem (1a) intends to minimize mismatch between the forecasted and restored active load demand (1b) and maximize the number of switches energized at an interval (1c). The first term in the objective function (1b) intends to minimize the difference between the forecasted load demand (P n ) and the restored active load power (P n ) over H . This allows to prioritize restoration of larger loads at interval h over horizon H . The second term maximizes the number of switches energized at an interval i.e. N R in a weighted manner. Here, α h represents weights assigned to the number of switches energized at different intervals of the prediction horizon (H ).
In this study p 1 and p 2 are set at 1. α h is 10 for all time intervals.

B. POWER BALANCE AND VOLTAGE EQUATIONS
The active power, reactive power and voltage equations of an unbalanced radial power grid are as follows [25].
Ignoring line losses, constraints (2a), (2b), (2c) and (2d) represent the power balance equations of n th bus. The active and reactive power flow (p.u.) through a line mn is represented by P mn and Q mn , respectively. The active and reactive power injections at N Der are indicated by P Der n and Q Der n , respectively. Since the DERs are all inverter connected source, the reactive power injections of DERs are limited by active power injections (P Der n ) and apparent power capacities of inverters, S max n (2d). Constraint (2e) defines the voltage difference equation for line mn. The bi-linear term (e.g., β mn V n ) of (2e) is linearized using Big-M method. Expressions for r P mn andx P mn are provided in (3), (4) and (5) [5].

C. DROOP CONTROLLED BS-DER
Along with restoration, the grid forming or BS-DERs also regulate the voltage and frequency levels of the energized islands. Constraints (6a) and (6b) represent active powerfrequency (P Der n −ω n ) and reactive power-voltage (Q Der n −V n ) droop controllers of BS-DERs. Based on the output power of DERs (P Der n , Q Der n ) the ω n and V n of DER integrated buses are determined by the droop controllers [26]. Here K P and K Q are the droop gains. P Der ref , Q Der ref , V ref and ω ref are the active power, reactive power, voltage and frequency references of the droop controllers, respectively. As it is a static study, the changes in the N Der frequency / voltage with changes in DER active power / reactive power injections are assumed to be instantaneous.

D. CONNECTIVITY TO BS-DERs AND SEQUENTIAL RESTORATION
Due to the presence of redundant lines, the topology of a reconfigurable distribution grid is not radial (Fig.1). Therefore, maintaining radiality of a partially restored grid during sequential restoration process becomes critical.
Single commodity flow model is utilized to ensure connectivity of all energized buses to BS-DERs (7a). Constraint (7a) ensures that every non-grid forming restored or energized bus of the distribution grid is connected to a bus with BS-DER. In this centralized leader-follower architecture all grid forming or BS-DERs operate as leaders and grid following DERs operate as followers [22]. Here, f is an auxiliary variable and holds no significance related to the actual power flow of the distribution grid. Variable f is bounded by (7b), making f zero only if a line is de-energized. The status of buses and lines are represented by binary decision variables v n and β mn , respectively. The status of β mn is dependent on the status of v m and v n (7c). The proposed BS&RS allows flexibility to the number of lines energized (N R ) at an interval (7d). A line energized at a time interval remains energized during the successive time intervals (7e). The number of islands in a distribution grid after complete restoration is determined by the number of BS-DERs in the distribution grid [27], [28]. Therefore, to maintain radiality of the distribution grid during sequential restoration, the number of energized lines at any time interval equals the number of BS-DERs and deenergized buses (7f) [19]. To restore the least priority buses or zero injection buses at the last time interval of sequential restoration, constraint (7g) is introduced. The frequencies of two adjacent buses (m and n) remain same if line mn is energized. However, if mn is de-energized the frequencies of m and n are determined by the the BS-DERs they are connected to (7h).

E. ENERGY STORAGE SYSTEM (ESS)
The efficiency of the MPC based BS&RS is assessed in the presence of quasi-static equations of ESS.
ESS is only charged or discharged if its terminal bus is energized (8a) and (8b). Constraint (8c) represents the dynamics of the SOC of ESS. η ch b and η dch b are the charging and discharging efficiencies of the ESS, respectively. Two binary variables γ ch ∈ {0, 1} |N b | and γ dch ∈ {0, 1} |N b | are introduced to make the charging and the discharging process of the restored ESS complementary (8d). The charging and discharging power P ch n (k+h|k) and P dch n (k+h|k), respectively of an ESS are limited by its capacity (8e-8f). The updated active power balance equations of N b are as per (8g).

F. CAPACITY CONSTRAINTS
Constraints (9a)-(9i) keep the variables of the above formulation within allowable limits. At any instant, the active power injection by BS-DER is limited by its minimum (P min ) and forecasted (P Der ) active power supply capacity (9a). The VOLUME 11, 2023 BS-DERs are assumed to be with ramp rate constraints, where (9b) are the upward and downward ramp rate constraints of active power injections by BS-DERs. The line flow limits are provided by (9c) and (9d). The active and reactive load demand at any bus is limited by the forecasted load demand of that bus based on (9e) and (9f), respectively. The bus voltages and frequencies of the distribution grid are limited by (9g) and (9h), respectively. All bus voltages and frequencies get non zero values only if the buses are energized. The maximum number of lines energized at an interval are limited by the total number of lines in the distribution grid (9i).

G. PROPOSED MPC BASED BLACK START AND RESTORATION SCHEME
A schematic of the proposed MPC based BS&RS is provided in Fig. 2. The algorithm is initiated during blackouts or in the presence of de-energized islands or widespread outages to sequentially restore a distribution grid. At each prediction horizon k = 1, 2, .., the input measurements are obtained from the physical layer. These measurements include the current state of controlled decision variables such as status of line switches, BS-DER set points, stored energy in ESS etc. The forecasted load profile and generation profile are also obtained. The optimization problem with objective function (1a) and constraints (2a)-(9i) is solved for every time interval of the prediction horizon. The solution gives the optimal switching sequence, power injections of DERs and SOC of ESSs over the time intervals [k, H ].
The first sample of the optimal control actions is applied to the system. The process is repeated for the shifted horizon (i.e., k = k + 1), till the grid is completely restored.

III. RESULTS AND ANALYSIS
Several numerical experiments are conducted in this section to evaluate the efficacy of the proposed MPC based BS&RS algorithm. The algorithm is coded in MATLAB-2020b [29], and solved using Gurobi Optimizer toolbox version 9.1.2 [30]. The study is carried out on a computer with Intel(R) Core(TM) i9-6200U CPU processor running at 2.30 GHz and 16 GB of RAM.

A. IEEE 123-BUS SYSTEM DESCRIPTION
The feeder is with 127 switchable lines and 123 buses. The rated power and rated voltage of the system are 1 MVA and 4.16 kV, respectively [31]. The base frequency is 60 Hz. For restoration of the distribution grid from complete blackout, the feeder is equipped with three BS-DERs and four grid following ESSs. The BS-DERs connected to buses 125, 123 and 60 are indicated as BS-1, BS-2 and BS-3, respectively. These DERs are assumed to be photovoltaic (PV) generators equipped with ramp rate constraints and droop controllers.  The duration of an interval ( t) is assumed to be 15 minutes, with time index h ∈ [1, 2, 3, 4, 5].

B. CASE STUDIES 1) MPC BASED BS&RS
In this section the proposed MPC based restoration algorithm discussed in Section II, is utilized for restoration of a distribution grid from complete blackout. To achieve this, the proposed algorithm coordinates three BS-DERs and four ESSs for sequential restoration of the distribution grid. At every time index h for a given k, all the topology, power balance, voltage and capacity constraints are satisfied. In this study, it is assumed that all DERs operate at their rated capacities and load multiplying factor varies in-between 0.5 to 1.1 p.u. The steps of the receding horizon MPC based sequential restoration algorithm are provided in Table 1. As listed in the table, the algorithm is iterated by shifting the prediction horizon (k). In this case study at k = 1 time index (h) is assumed to be 5 and h is reduced with the increase in k [32]. Fig. 3 illustrates the sequence of energized lines in the proposed MPC based BS&RS. In this approach, the number of lines energized at an interval are decided by the algorithm. At every time interval the optimal switching sequence retains connectivity of energized buses to one of the BS-DERs and maintains radial structure of the restored network. Taking into account the irregularities in the forecasted load demand, the algorithm identifies the lines to be energized to maximize load restoration at each time interval. At the final interval, all buses are energized but few switches remain open by the algorithm to enforce radial structure.
Due to brevity, the voltage and frequency of the distribution grid during the first four time intervals are provided in Fig. 4 and 5, respectively. It is evident from Fig. 4 and 5 that the number of buses operating at zero voltage and frequency level (buses not energized) are decreasing in successive time intervals.
The active power injections of BS-DERs are depicted in   lines, respectively. The increase in generation over the subsequent time intervals are evident from Fig. 6. The charging and discharging power and SOC of all four ESSs are illustrated in Fig. 7a and b, respectively. It is evident from Fig. 3 that only ESS-2 is restored during first time interval, however no charging / discharging takes place during this interval (see Fig. 7). ESS-1, 2 and 3 start charging during interval 2:00-3:00. ESS-1 and ESS-2 discharge during interval 3:00-4:00. However, ESS-3 holds its charge till the end of interval 3:00-4:00 and discharges during interval 4:00-5:00 (see Fig. 7b). All the three ESSs in the system participate in the restoration process and discharge completely at the end of interval 4:00-5:00. ESS-4 is also restored at interval 3:00-4:00.
The time varying forecasted active load demand and restored active load during sequential restoration are illustrated in Fig. 8. Figure shows the decrease in the load not  supplied over successive time intervals. However, due to limitations in the power supply capacity of the DERs the load demand is not supplied completely at the final interval.
The proposed receding horizon MPC based BS&RS algorithm does not require offline computation to determine the length of the restoration horizon. It also does not require to specify the number of lines to be switched in during an interval.

2) EFFICIENCY OF MPC IN ENHANCING DISTRIBUTION GRID RESILIENCY
Black start generators are expected to ramp up faster to reduce restoration time and unsupplied demand [7]. To assess the efficiency of BS-DERs in supporting the rapid change in load demand, the MPC based black start restoration scheme is compared with time series restoration scheme proposed in [19]. This case study assumes sequential restoration of a power grid from complete black out with sudden large change in forecasted load demand. For this case study, the prediction horizon is divided into 5 time indices. The maximum number of lines energized / repaired at an interval in time series restoration scheme is assumed to be 127 5 ≈ 26. Here, 127 is the number of lines in the distribution grid and 5 is the restoration horizon. As depicted in Fig. 9c, at interval 3:00-4:00 the load multiplying factor representing the forecasted demand changes from 0.5 to 1.1. The ramp in the generation of BS-DERs with the proposed MPC based BS&RS and time series BS&RS are provided in Fig. 9a and Fig. 9b, respectively. 69184 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE. Restrictions apply.   The efficiency of these two algorithms in restoration of distribution grid is quantified using an index, load not supplied (LNS) (10).
As stated before (section II-G) to take control action at time index (h) the MPC algorithm assesses the response of the system over the entire prediction horizon (H). Therefore, to support the increase in the future load demand in the proposed MPC based algorithm, the BS-DERs start ramping up ahead of time. As a result, compared to the time series black start scheme, the generators support more loads and LNS also reduces in the proposed MPC based algorithm (see Fig. 9d). The lines energized at different time intervals of the MPC based restoration algorithm and time series restoration scheme are listed in Table 2.

3) EFFECT OF SOLAR GENERATION AND LOAD PROFILE ON BS&RS
The effect of variability in generation capacity of BS-DERs and load demand in the sequential restoration scheme are assessed in this section. All BS-DERs are following a forecasted power generation capacity profile illustrated in Fig. 10. Therefore, the power supply capacity of the BS-DERs change at different hours of a day. The load multiplying factor representing the forecasted load demand is also depicted in Fig. 10. In this study, the proposed black start restoration scheme is evaluated at different hours of a day. The number of energized lines at different intervals of a horizon (N R ) are provided in Table 3. These energized lines and its locations are variable and are determined by several different factors such as ramp rate constraints of the DERs, forecasted PV generation profile, forecasted load demand, SOC of ESS, and α h . The efficiency of the sequential load restoration algorithm is illustrated in Fig. 11. As depicted in Fig. 11, at hours (9)-554 (10), (13)-14), (16)(17) and (17)(18), unmet load demand exist due to reduction in solar generation and increase in demand. With the rapid decline in solar generation at hour 17-18, the LNS increases at the intervals of H. It is evident VOLUME 11, 2023    from the study that the algorithm maintains connectivity and radiality of the distribution grid but fails to restore loads completely at different time intervals of a day. The presence of larger ESSs can improve the resiliency of the distribution grid in the presence of PV generator based BS-DERS.

C. DETERMINISTIC WORST-CASE BS&RS
Due to the intermittency in DER generation and errors associated with forecasting load demand, forecasted DER capacity and load demand at an interval are expressed as followŝ In (11), σ represents peak value of the forecasting error and it is assumed to be constant over the entire H . Here, (11a) and (11b) represent the forecasted power supply capacity of DER and active load demand at n th node, taking forecasting errors into account. As this study is suited for power system operator's decision support, a deterministic solution is desired. Therefore, a deterministic worst case MPC based BS&RS scheme is proposed in this section. For the deterministic worst case MPC, four boundary scenarios are considered [34].
A study is carried out for three different values of example uncertainties, σ = 0%, ±20%, ±30%. The effect of forecasting errors of generation and load demand can adversely affect the system's voltage and frequency profiles due to increased mismatch in DER power supply capacity and load demand. The effect of different percentages of forecasting errors on voltage and frequency of the restored system are depicted in Fig.12. It is evident from Fig.12 that with the increase in the σ value, the voltage of the restored system dropped below 0.92 p.u. due to limited apparent power supply capacities of inverter based BS-DERs. However, the effect of σ on the system frequency is not significant. Because, even during the worst case scenario, the active power generation is limited by the ramp rate constraints and the change in system frequency is limited by the ramp rates and droop gains of BS-DERs.

D. SCALABILITY AND REAL TIME CAPABILITY
Scalability is challenging in MPC based algorithms due to increase in the size of the problem with increase in the prediction horizon. Due to the mixed integer nature of the problem and associated computational burden, the feasibility for scalability and real time implementation is assessed in this section. The proposed algorithm (P1) is studied for black starting modified IEEE 123-bus and 1069-bus systems for different prediction horizons.

1) 1069-BUS SYSTEM DESCRIPTION
1069-bus system has four feeders, 1323 lines and 1317 buses (see Fig. 13) [35]. There are 176 switches in the distribution grid. All the switchable lines are shown in Fig. 13 Table 4. For improved decision support, it is suggested to increase the prediction horizon of MPC based algorithms [36]. It is evident from the study that the computation time increases with increase in the prediction horizon and size of the system which makes the real time implementation challenging. To keep computation time smaller than t = 15 minutes, it is suggested to select the prediction horizon as 3.

IV. DISCUSSIONS AND CONCLUSION
This work presents a MPC aided black start service restoration strategy for enhancing resiliency of an unbalanced distribution system with DERs. Based on the fluctuations in the future generation and load demand, the proposed  approach allows real-time control actions with the aim of identifying optimal switching sequence and DER set point to sequentially restore a distribution grid. Case studies show that the proposed MPC based sequential restoration scheme improves distribution grid resiliency compared to existing sequential restoration schemes. Though the algorithm is studied in this work for restoration of a distribution grid from complete blackout, it can be applied for restoration from partial or large scale outages. The novelty of the algorithm also lies in incorporating frequency regulation in sequential service restoration scheme. The MPC based sequential restoration scheme does not require offline study to determine the duration of the restoration process because of receding horizon technique. Due to the presence of power flow equations, constraints associated to reconfiguration and dynamics of ESSs, the problem is non convex in nature. The switching sequence, DERs set points and status of SOC of the problem are achieved at each time interval by solving a linearized convex model of the original problem. The proposed scheme is extended to incorporate errors associated with forecasting of generation capacity and load demand. The algorithm is able to identify valid switching sequence considering uncertainty of forecasted generation and demand. The proposed method performs efficiently when applied to IEEE 123-bus and 1069bus systems. The study also finds a trade-off between the accuracy of the solution and real time implementation. Future work will comprise distributed learning based restoration approaches.

APPENDIX A BS-DER PARAMETERS
See Table 5.