A Model Predictive Control Strategy for Distribution Grids: Voltage and Frequency Regulation for Islanded Mode Operation

: In the last few years, one of the most important challenges of power technologies has been the integration of traditional energy production systems and distributed energy resources. Large-scale photovoltaic systems and wind farms may decrease the quality of the electrical grid service, mainly due to voltage and frequency peaks and ﬂuctuations. Besides, new functionalities, such as the operation in islanded mode of some portions of the medium-voltage grid, are more and more required. In this respect, a model predictive control for voltage and frequency regulation in interconnected local distribution systems is presented. In the proposed model, each local system represents a collection of intelligent buildings and microgrids with a large capacity in active and reactive power regulation. The related model formalization includes a linear approximation of the power ﬂow equations, based on stochastic variables related to the electrical load and to the production from renewable sources. A model predictive control problem is formalized, and a closed-loop linear control law has been obtained. In the results section, the proposed approach has been tested on the Institute of Electrical and Electronics Engineers(IEEE) 5 bus system, considering multiple loads and renewable sources variations on each local system.


Introduction
The increase of intermittent renewable energy sources (RES) has created instability issues, requiring new controllers for modern smart microgrids [1][2][3]. Actually, RES may cause voltage/frequency peaks and fluctuations that negatively affect the quality of the electrical grid service. On the other hand, RES represents an important alternative to reduce global fossil fuel consumption. The need for new control strategies is also motivated by the increase of distributed small generators that cause millions of new producers potentially involved in the energy market. In such a scenario, in the balancing market, new actors, such as aggregators [4], are present and new functionalities are more and more required. For example, in order to preserve the electrical grid from faults and emergency situations, it may be necessary to reduce the power demand and to operate some portions of the medium-voltage grid in islanded mode. These new set-ups require new fast controllers, which can take into account voltage and frequency models.
In addition, the grid operators, i.e., the transmission system operator (TSO) and distribution system operator (DSO), are forced to deeply change their roles. As an example, the DSO must become active in its activities as a grid manager [5]. The DSO, in particular, is responsible for the normal operation of the distribution grid and its subsystems (i.e., managing power quality and energy losses).
Among others, the new roles for the DSO are: operational security management, dynamic voltage and frequency control, outages response, performing restoration actions, and market operations. In this framework, it becomes hard to define an effective optimal control, which derives from the structure of the power grid, and, specifically, from the presence of several issues: renewable and traditional power production, bidirectional power flows, and stochastic modelling issues (due to uncertainties in forecasting power generated from renewable resources). The main challenge is to define and to solve control problems embedding prediction models for the regulation of frequency, voltage, and active and reactive power flows. Indeed, the complete power flow equations are nonlinear and nonconvex [6]. Therefore, optimization problems can hardly apply them to compute the control in real time. On the other hand, there is the need to represent in detail the physical system to provide reliable control strategies [7].
This work focuses on the optimal control of interconnected local systems in a distribution grid. A local system can represent a group of renewable generators, passive and active loads, or microgrids with a large capacity for active and reactive power regulation. The proposed approach can be used: • To reduce frequency oscillations and to improve inertia response due to a sudden change in the grid (i.e., a variation on a renewable power plant production), and • To manage a portion of the distribution grid that is operated in islanded mode because of demand response programs and/or emergency requirements.
In fact, an increase in the RES penetration can decrease the number of generation units that provide reserve power for primary and secondary control. So, the frequency deviation increases [8]. New equations for inverters' control to regulate frequency and voltage are required when a change in the grid occurs (e.g., a loss of load).
The modelling approach proposed in this paper is based on the concept of virtual inertia [9][10][11]]. An electrical model is defined, including voltage (magnitude and phase), frequency, active and reactive power, and stochastic components. So, the definition of an optimal control problem for frequency and voltage regulation in multi-local systems is introduced. The aim is to define a closed-loop controller for the mitigation of frequency and voltage variations.
The paper is organized as follows. Section 2 reports the state-of-the-art and discusses the innovation of the proposed work. In Sections 3 and 4, the system and the optimal control problem are formalized. Finally, the application to a case study is described in Section 5, while conclusions and future developments are reported in Section 6. The Appendix A reports proofs related to demonstrations not included in the text for sake of readability.

Literature Review
Nowadays, the control of interconnected local distribution systems is a challenging problem. In this respect, several papers focus on the definition and solution of optimization and control problems. A detailed survey [12] on multi-microgrids systems has been recently published, where operational management and control are quoted as challenging problems to be faced. Zheng et al. [13] propose a distributed model predictive control (MPC) approach using discrete-time Laguerre functions for load frequency control in multi-area interconnected power systems, in order to overcome the problems of computational burden and online optimization. The advantage of an MPC approach is that it allows for implementing a multivariable controller that controls the outputs simultaneously by taking into account all the interactions among system variables. Another strength of MPC (not used in our work, where a closed fast feedback is adopted) is that it can handle constraints. Khooban and Niknam [14] propose a heuristic algorithm and a fuzzy logic approach to tune parameters of classical proportional integral (PI) controllers for multi-area load frequency control. However, electrical models are not detailed. In the literature of microgrids, there are several interesting approaches based on MPC but related to the single microgrid, and thus not interconnected with other local systems. In Reference [15], a novel fractional-order model predictive control method is proposed to achieve the optimal frequency control of an islanded microgrid by introducing a fractional order integral cost function into the MPC approach. As in the present work, frequency regulation is attained, even if hereinafter we propose a more detailed system model, where active and reactive powers are present and voltage is regulated.
In Reference [16], a feedback linearization is applied to design a predictive controller for the secondary control of the voltage and frequency in an islanded microgrid. The goal is to improve the performance of a closed-loop nonlinear system. The control scheme is fully distributed, and necessary and sufficient conditions for convergence and stability of the whole system are described. As a similar approach, we model both active and reactive power, but providing a control strategy for the whole optimization problem that can be used under a MPC scheme. In Reference [17], a MPC-based controller uses a simplified voltage prediction model to predict the voltage behavior in an islanded microgrid. Apart from the fact that we consider multiple microgrids, in our model, we use a more detailed electrical model and we derive an optimal control law. The same considerations can be done in comparison with References [18,19].
The approach considered in the present paper gives a strong emphasis on the model of the electrical grid. Specifically, the approach is based on the concept of virtual inertia [20,21]. The problem of virtual inertia control is mainly felt in microgrids and grids where RES are present, as it is necessary to cope with small inertia and uncertainties. In Reference [22], the authors propose an approach to increase the inertia of the photovoltaic (PV) system through inertia emulation, which can be obtained by controlling the charging/discharging of the DC-link capacitor and adjusting the PV generation. Wu et al. [23] consider an islanded microgrid with renewables and storage systems, in which virtual inertia allows power curtailment of RES units and demand response (loads' shifting). In fact, these approaches are also useful for islanded microgrids, and several works are present in the recent literature for frequency and voltage regulation, and for stability [24][25][26][27][28].
In Reference [24], the authors assess that droop control methods are not suitable for renewable energy-based microgrids due to their limited capability for power delivery. They propose a supercapacitor-based DC-link voltage stabilizer and a battery energy storage system-based frequency stabilizer. In particular, as in Reference [25], a fuzzy-based controller is proposed. In Reference [26], the authors highlight the drawbacks of PI controller gains and propose the Grasshopper Optimization Algorithm to optimize the PI controller parameters.
An additional critical issue is represented by the demand of electric vehicles in a smart grid. To face such problem, in Reference [27], distributed controllers are proposed. Other papers [28] focus on improvements of droop controllers and on new approaches for the optimal dispatch in hierarchical microgrids [25]. Specifically, two control modes are considered for the generation units: a voltage control mode, with an active droop control loop, and a power control mode, which allows setting the output power in advance. Some recent works that study virtual inertia are References [29][30][31][32][33]. In Reference [29], the authors study a novel control scheme for virtual inertia converters to operate in unbalanced settings, operating test on real hardware. In Reference [30], an adaptive controller method for a PV-based virtual inertia system is discussed. The approach is studied in the case of the small-signal stability and low-frequency oscillation, which are characteristics of the PV generators. In Reference [31], the authors aim at improving the inertia of a microgrid by reducing the time delay of the communication system to overcome frequency oscillations, finding a compromise with the increase of operational costs.
In Reference [32], the authors propose an approach to find the proper value of inertia constant, together with the frequency droop coefficient of DERs (distributed energy resources) and tuning load frequency controllers' parameters to improve the frequency stability. A multi-objective decision problem is formulated in which the following goals are considered: minimizing the maximum frequency deviation while promptly bringing the frequency back to the desirable region, and minimizing a higher inertia constant that results in a higher cost of energy storages. With a significant difference, here, we propose an optimal control problem to minimize frequency and voltage deviations for multiple microgrids and the optimization problem is analytically solved taking into account both active and reactive power.
In Reference [33], virtual inertia in islanded microgrids is studied. An H-infinite (H∞) robust control approach to the virtual inertial control loop is implemented, taking into account the high penetration of renewables, thus enhancing the robust performance and the stability of the microgrid during contingencies. In the present paper, we propose a control law to minimize frequency and voltage deviations taking into account both active and reactive power under an MPC scheme.
Hereinafter, the electrical grid connecting different local systems is modelled. This paper aims at studying the regulation of voltage and frequency within local systems interconnected through a distribution network. Since these systems have a strong renewable, uncontrolled, and intermittent component, they intrinsically represent a disturbance that, if not controlled, may cause a decrease in the system performance and quality. It is worthwhile to observe that the problem faced is not to be confused with the classic problem of inter-area oscillations between systems connected by a high-voltage grid.
In this work, an MPC approach is proposed, taking into account uncertainties, and deriving an optimal solution that can be used in closed-loop and in real time. In the literature of power systems, similar control approaches can be found [34,35]. Delfino et al. [36] propose a multilevel architecture for interconnected microgrids in which the DC power flow equations are used at a higher level, an analytical solution emulates each microgrid, and the AC power flow equations are used at the lower level. Minciardi and Robba [37] formalize a decision model in which a stochastic optimal control problem is defined and solved. An optimal control strategy can be also found in Reference [38]. The authors consider the problem of load-frequency control in a multi-area power system proposing a hierarchical optimal robust controller. They decompose the overall problem in subsystems, solve the lower level through a control strategy, and finally, iteratively coordinate local systems. However, a detailed representation of the electrical network is not present. Annaswamy and Kiani, in Reference [39], include a discrete-time decentralized linear quadratic regulator (LQR) controller in a transactive architecture for hierarchical smart grids. In the present paper, with respect to References [34,35], a detailed representation of the electrical network is included in order to regulate voltage, and control active and reactive power. With respect to References [36,37], the decision problem has been solved analytically. In addition, References [34][35][36][37][38][39] do not consider frequency variations and equivalent inertia. Other optimal control strategies for the smart grid are present in the literature [40,41] but not in the area of interconnected local systems.
In summary, the main contributions of the proposed paper are: • The statement of the frequency and voltage control problem of distribution grids in islanded mode expressed as an MPC of an affine system affected by additive noise.

•
The inclusion of models for the control of frequency, voltage, and active and reactive power in distribution grids.

•
The derivation of an analytical solution for the stochastic MPC problem based on a two-stage decomposition.

•
An optimal control strategy for the inter-local systems variation problem that can be used either for frequency and voltage regulation or the operation in the islanded mode of portions of the distribution grid for demand response purposes.

The System Model and the Optimization Problem
Consider a distribution power grid ( Figure 1) modeled as an undirected graph G = (A, L), where the set of nodes A = {1, . . . , N} characterizes the set of local systems and L = {1, . . . , b} is the set of links that connect the nodes. Specifically, each local system is interfaced with the distribution grid by a power converter which can perform a virtual inertia control (allowing the inverter to behave as a fictitious synchronous generator) with equivalent inertia time constant H i (i ∈ A) and damping factor KD i (i ∈ A). The set of arcs, L, represents the set of lines that connect the local systems-each line has a proper impedance Z i,j = R i,j + jX i,j (i, j ∈ A), with R i,j (i, j ∈ A) being the line resistance and X i,j ( i, j ∈ A) being the line inductance. All local systems are composed of different devices: distributed controllable resources with high regulation capacity (denoted by G), reactive power-regulating sources (denoted by NG), renewable non-controllable resources (denoted by RES), and passive loads (denoted by L). Each local system can be modeled as an equivalent synchronous generator with equivalent inertia to consider its contribution to the frequency regulation of the whole grid.
Energies 2020, 13, x FOR PEER REVIEW 5 of 30 passive loads (denoted by L). Each local system can be modeled as an equivalent synchronous generator with equivalent inertia to consider its contribution to the frequency regulation of the whole grid.
In the following paragraphs, the state and the control variables are defined, and the optimization problem is formalized. All the quantities (but the phases) and the equations are expressed at discrete time per unit (p.u.) values. As suggested by the literature, this method offers computational simplicity by eliminating units and expressing system quantities as dimensionless ratios. The decision variables are: − : reactive power injected at node i.
: distributed controllable source active power at node i.
: reactive power regulating source reactive power at node i.
Pm : regulation active power injected by local system i.
The parameters are: : load active power at node i.
: load reactive power at node i.
: distributed controllable source reactive power at node i.
: reactive power regulating source active power at node i.
The state variables are: The decision variables detail for a local system is reported in Figure 2. In the following paragraphs, the state and the control variables are defined, and the optimization problem is formalized. All the quantities (but the phases) and the equations are expressed at discrete time per unit (p.u.) values. As suggested by the literature, this method offers computational simplicity by eliminating units and expressing system quantities as dimensionless ratios.
The decision variables are: . . , T − 1: reactive power regulating source reactive power at node i. • Pm i,t : regulation active power injected by local system i.
The parameters are: . . , T − 1: reactive power regulating source active power at node i.
The state variables are: The decision variables detail for a local system is reported in Figure 2. The electrical grid inside each local system i is not considered, and thus only the power balance must be satisfied: where the frequency is related to the active power supply-demand imbalance: where i KD is the local system damping factor. The distribution grid model includes the linearized power flow equations for a medium voltage AC grid, with as the line susceptance. That is: The model described by Equations (5) and (6) is the well-known linearized version of power flow equations, obtained considering the first term of Taylor expansion of the sinusoidal terms of the equations and just considering the voltage difference between nodes [42][43][44]. In the following paragraphs, the linear model described by the previous equations is reformulated in the state space canonical form.
Proposition I: The system model, described by Equations (1)-(6), can be expressed as an affine system affected by additive noise: The electrical grid inside each local system i is not considered, and thus only the power balance must be satisfied: Each local system is associated with a time-varying voltage phasor of magnitude V i,t and phase δ i,t . The discrete-time system dynamics (∆t is the time interval length) are governed by the equations of the virtual inertia controller, linking the node angle, δ i,t , to the frequency difference to the nominal value (ω i,t − ω 0 ): where the frequency is related to the active power supply-demand imbalance: where KD i is the local system damping factor. The distribution grid model includes the linearized power flow equations for a medium voltage AC grid, with G i,j = 1/R i,j (i, j ∈ A) as the line conductance, and B i,j = −1/X i,j (i, j ∈ A) as the line susceptance. That is: The model described by Equations (5) and (6) is the well-known linearized version of power flow equations, obtained considering the first term of Taylor expansion of the sinusoidal terms of the equations and just considering the voltage difference between nodes [42][43][44].
In the following paragraphs, the linear model described by the previous equations is reformulated in the state space canonical form.
Proposition I: The system model, described by Equations (1)-(6), can be expressed as an affine system affected by additive noise: with where A, B, C, andĈ are known matrices, h t is a vector of parameters, andĥ t represents stochastic noise with E ĥ t = 0. Proof of Proposition I: The Proof is shown in the Appendix A.1.
In the day ahead, the DSO plans the production and energy flow based on different objectives, such as the minimization of losses and the capacity of the lines. It is reasonable to consider variations concerning an already predetermined working point. However, given the strong variability of part of the energy sources involved in our problem, control policies are required to mitigate deviations from an initial working point.
The proposed objective function aims at minimizing four different contributions: the square of the difference between the phases of the grid nodes, the quadratic deviation of the frequency from the steady-state value, the voltage quadratic deviation from the desired value, V 0 , and the quadratic deviation of the active regulation power to the set point Pm * : Proposition II: The objective function (8) can be expressed as a stochastic quadratic function with mixed terms: where x t and u t are the state and control vectors respectively,x,û are the tracking parameters, and Q t , R t , F t , Q T are known matrices. Proof of Proposition II: The Proof is shown in the Appendix A.2.

The MPC Scheme
The MPC scheme is defined in Figure 3. The optimization horizon (T) is the length over which the optimization problem runs and (S) is the simulation horizon. Both horizons (T) and (S) can be chosen according to the available information about uncertain parameters/forecasts (current states, forecasted demands, and renewables). In the MPC approach, in the first run, the solution related to the first time interval is considered. Then, uncertain parameters (new information from sensors in the field and new forecasts) are inserted in the procedure and a new optimization problem is run (second run), and so on for the whole duration of the simulation horizon (S). The MPC approach allows for updating new information when available and reducing uncertainties related to renewables' availability and load forecasting.

Explicit Solution of the MPC Problem
The optimization problem is in the form: subject to (s.t.) Equation (7).

Explicit Solution of the MPC Problem
The optimization problem is in the form: subject to (s.t.) Equation (7).

t t t t t t t t T T T
where: The Equations (11)-(13) are a stochastic MPC problem, which is "non-standard," due to the presence of the known input sequence, t μ , in the state equation. In the following claims, it is shown how to solve the MPC problem in a closed-form through a two-stage decomposition. Claim I: Due to its linearity, the decision problem in Equations (11) s.t.  By a change of variables z t = x t −x and e t = u t −û t , the optimization problem is given by: where: The Equations (11)-(13) are a stochastic MPC problem, which is "non-standard," due to the presence of the known input sequence, µ t , in the state equation. In the following claims, it is shown how to solve the MPC problem in a closed-form through a two-stage decomposition.
The optimal control strategy is given by: z * Claim II: The optimal control problem to determine e d, * t , due to its linearity, can be written as: where: Proof of Claim II: The Proof is shown in the Appendix A.4. Proposition III: The optimal control law for the problem in Equations (11)-(13), and taking into account Claim I and Claim II, is given by: with being: Proof of Proposition III: The Proof is shown in the Appendix A.5. As reported in Claim I and Claim II, the original problem can be so decomposed into a stochastic problem and a deterministic discrete tracking problem. The closed-loop optimal control strategy is obtained as follows [46]: Sub-problem 1: Stochastic problem. The optimal control problem to determine e s, * t is given by: s.t. z s t+1 = Az s t + Be s t +Ĉĥ t t = 0, . . . , T − 1 Proposition IV: The solution of the optimal control problem, Equations (36)- (38), is given by the following recursive equations: e s * t = −K t z s * t t = 0, . . . , T − 1 (39) with Proof of Proposition IV: The Proof is shown in the Appendix A.6. Sub-problem 2: Deterministic problem. As reported in Claim II, the optimal control problem to determine e d, * t can be written as: The optimal control strategy for the problem, Equations (43)- (45), is given by with

Case Study Description
To show the effectiveness of the theoretical results presented in the previous section, the proposed control strategy has been applied to a modified IEEE 5-bus system (N = 5) [47], shown in Figure 4. The test network has been adapted to our problem through the following assumptions. Each bus corresponds to a local system, and, keeping the network topology unchanged, the values of cables' resistance, R i,j , and reactance, X i,j , have conveniently rescaled to a distribution grid system (assuming base quantities of 10 MVA for apparent power, V 0 = 15 KV for voltage and ω 0 = 50 Hz for frequency). Table 1 shows the detailed values of the line parameters R i,j, and X i,j for the case study. Energies 2020, 13, x FOR PEER REVIEW 12 of 30  The parameters that characterize each local system, i.e., the inertia constant, Hi, and the equivalent damping factor, KDi, are reported in Table 2. The local system 5 is the slack node, in which the voltage, t V , is fixed to 1 (p.u.) and the phase, t δ , is equal to 0 (rad). The case study can be representative of small local systems connected by a distribution grid (such as interconnected microgrids), in which there is a great presence of generation from renewable resources. The simulation horizon is 60 s, with a discretization step t Δ of 0.5 s. This horizon length has been chosen as a trade-off between the following two aspects: (a) to avoid delays due to the centralized communication from the controller to the local systems, and (b) to enhance the reliability of the short-term forecasting of RES and loads.
To make the approach more effective, the MPC strategy is applied, in which only a part of the optimal control law is applied (e.g., 10 seconds over a T = 60 second horizon), and there is a continuous recalculation of the control law based on new measurements and forecasting. The length of the MPC control horizon is justified by recent developments on the nowcasting of renewable power sources that predict 1 min future PV values [48,49].
At this point, by the means of the local systems' parameters and the distribution network, it is possible to quantify the matrices A and B, defined in the Appendix, Section A1, described by Equations (A19) and (A20), and needed by Equation (7):  The parameters that characterize each local system, i.e., the inertia constant, H i , and the equivalent damping factor, KD i , are reported in Table 2. The local system 5 is the slack node, in which the voltage, V 5,t , is fixed to 1 (p.u.) and the phase, δ 5,t , is equal to 0 (rad). The case study can be representative of small local systems connected by a distribution grid (such as interconnected microgrids), in which there is a great presence of generation from renewable resources. The simulation horizon is 60 s, with a discretization step ∆t of 0.5 s. This horizon length has been chosen as a trade-off between the following two aspects: (a) to avoid delays due to the centralized communication from the controller to the local systems, and (b) to enhance the reliability of the short-term forecasting of RES and loads.
To make the approach more effective, the MPC strategy is applied, in which only a part of the optimal control law is applied (e.g., 10 seconds over a T = 60 second horizon), and there is a continuous recalculation of the control law based on new measurements and forecasting. The length of the MPC control horizon is justified by recent developments on the nowcasting of renewable power sources that predict 1 min future PV values [48,49]. At this point, by the means of the local systems' parameters and the distribution network, it is possible to quantify the matrices A and B, defined in the Appendix A.1, described by Equations (A19) and (A20), and needed by Equation (7): In each MPC step, the optimal control law is the sum of the solution of two different control problems: a stochastic problem defined by Equations (36)- (38), with the optimal solution, Equations (39)- (42), and a deterministic problem, Equations (43)- (45), with the optimal solution, Equations (46)-(53). Both control laws have been recursively defined, i.e., all the parameters that define the control can be calculated offline, the parameters are known, and the forecasts of renewable generation and load.
The variations on the system are caused by the combined variation of load and RES Psc ,i,t on the local systems defined by Equation (A4), represented into the model in the terms h t andĥ t defined by Equations (A21) and (A22) respectively, deterministic and stochastic parts of the noise. Figure 5 reports the values, Psc d i,t , for the whole simulation horizon. Data have been obtained by adapting real power measurements from photovoltaic, residential, and industrial loads, and rescaled to obtain congruent measurements for the distribution networks. The stochastic part of the noise, Psc s i,t~N (0,1), affects the state at each instant t.

Optimal Control Results
Optimal solutions for the state variables in the case with control or without control are reported The power profiles on each local system could be critical for the whole system. Without any upper-level controller but only with the damping given by each local system, the generators could not be able to efficiently compensate (i.e., in a reasonable time, according to the dynamics of the system) the given power variation, causing frequency and voltage variations.
The objective of the proposed controller is to reduce the voltage and frequency variations (during the transient period) and, after the perturbation, bringing them back to their reference value. In a practical context, numerous frequency variations and voltage peaks on the grid could cause the intervention of protection systems, sensitive to variations in voltage or frequency, disconnecting some portions of the network itself, causing a blackout for many users.

Optimal Control Results
Optimal solutions for the state variables in the case with control or without control are reported in Figures 6 and 7. It can be seen how the presence of the proposed control can limit frequency variations (i.e., the power quality remains within acceptable limits, the voltage deviation from the reference value is less than 5%). The dotted line represents the system with no control (only the damping coefficient is considered). In this case, each local system operates without any central controller (like the one developed here), but the system only reacts based on the physics of the system with constant inputs equal to those of the equilibrium point. The bold line represents the controlled system using the developed approach. It is important to note that without the proposed control strategy, the system does not reach a steady state for the whole simulation.
Energies 2020, 13, x FOR PEER REVIEW 15 of 30 Figure 6. Phase optimal results at each local system. Figure 6. Phase optimal results at each local system. Figure 8 shows the voltage pattern for each local system in the presence or absence of the developed control. The control strategy helps in containing the variations affecting this variable, improving the system's power quality. Finally, Figures 9 and 10 show the optimal active and reactive power at each local system.    When there is no control, the objective function value is equal to 0.1098, while when the proposed control strategy is applied, it is 0.0151. Table 3 shows the objective function values and contributions. The test has been developed in MATLAB-SIMULINK (2018b, Mathworks, Natick, Massachusetts, USA), running on a PC Intel i7, 16 GB RAM. The MATLAB tool has been used to calculate a priori the control laws defined recursively in Section 4.

Closed-Loop Stability Analysis
This subsection aims to perform a sensitivity analysis concerning the asymptotic stability properties of the proposed system, with respect to the damping parameter, KD i . The choice of this parameter is fundamental with regards to the design of the control system based on virtual inertia. It is well known from power systems theory that the asymptotic stability of a synchronous machine is directly linked to its damping factor, similarly for the virtual inertia that mimics this behavior.   The test has been developed in MATLAB-SIMULINK (2018b, Mathworks, Natick, Massachusetts, USA), running on a PC Intel i7, 16 GB RAM. The MATLAB tool has been used to calculate a priori the control laws defined recursively in Section 4.

Closed-Loop Stability Analysis
This subsection aims to perform a sensitivity analysis concerning the asymptotic stability properties of the proposed system, with respect to the damping parameter, KDi. The choice of this parameter is fundamental with regards to the design of the control system based on virtual inertia. It is well known from power systems theory that the asymptotic stability of a synchronous machine is directly linked to its damping factor, similarly for the virtual inertia that mimics this behavior.
The first part of Figure 11 shows how different choices of damping factors at each node affect It is noteworthy that for lower values of KD i , some eigenvalues have a norm greater than 1, so the system is unstable. This is consistent with the power system theory. To this end, the asymptotical closed-loop stability of the proposed controller for an infinite time control law defined is analyzed, as in Reference . As can be seen from the second part of Figure 11, the proposed optimal controller also guarantees the asymptotical stability in the cases in which the system was designed as unstable. This aspect is crucial for the optimal control of this kind of problem, as it also fixes design instabilities.

Conclusion and Future Developments
The proposed work introduced an approach based on optimal control for the regulation of voltage and frequency in multi-local systems (such as interconnected microgrids) in a distribution grid. This represents one of the most challenging problems regarding the dynamic aspects of a smart grid.
If not considered and properly mitigated, these aspects can drastically damage the entire system. To reduce these effects, the paper considers accurate grid models and the presence of innovative devices that allow for better regulating the power flows on the grid and, consequently, voltage and frequency at each node. It is worthwhile to observe that the main contribution of this paper is not to assess the controller performance (that is optimal by definition to the linearized model) but to propose a model that can be linearized and controlled by a closed-form control law.
To cope with the variability represented by renewable sources and load, this approach was used in an MPC framework. In this case, just the control law corresponding to the first time interval would be applied, while an update of the law would be performed according to the new forecasts made in the short term. This would allow a further increase in the quality of the system. The IEEE 5 bus system was used to test the control strategy, considering multiple loads and renewable sources variations on each local system. The results showed a strong decrease in the number of variations as well as strong damping of the frequency and voltage peaks in each local system. This greatly increases the power quality of the electrical system. A sensitivity analysis about the asymptotical stability was also given, showing that the proposed control can also prevent instabilities caused by an incorrect design of the damping parameters.

Conclusions and Future Developments
The proposed work introduced an approach based on optimal control for the regulation of voltage and frequency in multi-local systems (such as interconnected microgrids) in a distribution grid. This represents one of the most challenging problems regarding the dynamic aspects of a smart grid.
If not considered and properly mitigated, these aspects can drastically damage the entire system. To reduce these effects, the paper considers accurate grid models and the presence of innovative devices that allow for better regulating the power flows on the grid and, consequently, voltage and frequency at each node. It is worthwhile to observe that the main contribution of this paper is not to assess the controller performance (that is optimal by definition to the linearized model) but to propose a model that can be linearized and controlled by a closed-form control law.
To cope with the variability represented by renewable sources and load, this approach was used in an MPC framework. In this case, just the control law corresponding to the first time interval would be applied, while an update of the law would be performed according to the new forecasts made in the short term. This would allow a further increase in the quality of the system. The IEEE 5 bus system was used to test the control strategy, considering multiple loads and renewable sources variations on each local system. The results showed a strong decrease in the number of variations as well as strong damping of the frequency and voltage peaks in each local system. This greatly increases the power quality of the electrical system. A sensitivity analysis about the asymptotical stability was also given, showing that the proposed control can also prevent instabilities caused by an incorrect design of the damping parameters. Future developments will concern several aspects. First, delays in the communication between the central controller and each local system could be inserted in the optimal control problem.
The particular model of the local systems will be the focus of future investigations. Furthermore, distributed and team robust approaches to overcome the poor robustness of totally centralized control approaches in case of emergencies will be considered.  (5) and (6) in the absence of shunt admittances can be written in the compact form: where K p and K q are known parameters coming from considering the voltage slack bus equal to 1. Equation (1) can be inserted in the frequency state Equation (4): being where Psc d i,t and Psc s i,t are the deterministic and stochastic parts of the net load, respectively. The frequency state equation is given by: being: Substituting Equation (A1) in the frequency state equation, we obtain: then, the voltage can be obtained from the Equation of reactive power (A1). That is: and substituting Equation (A7) in the frequency state Equation (A6) we have: substituting Equation (2) (where the control variable is QNG i,t ) into Equation (A8), we get: Let us pose: then, we obtain: It is important to note that Q TN,t = [QTN i,t , i = 1, . . . , N − 1] is an auxiliary variable that is crucial to define the objective function accurately (as it will be explained in the following). Thus, in the state equation, we have to satisfy the following condition: B ω GB −1 p . Then, we obtain: as regards the other state equation, one can write: , it is possible to write: where: Appendix A.2. Proof of Proposition II As regards the objective function, we consider the four terms in Equation (7): • The first term can be written as δ T t Q δ δ t .

•
The fourth term is given by (P m,t − P * m ) T R p (P m,t − P * m ).
With P * m = [Pm * i , i = 1, . . . , N − 1] where Q δ Q ω and R p are the identity matrix of appropriate dimensions.
The third objective is the most complicated and it is given by (substituting voltage with Equation (A7)): (−B −1 p (Q G,t + Q NG,t + Q RES,t − Q L,t − K q ) − B −1 p Gδ t − I) T R v (−B −1 p (Q G,t + Q NG,t + Q RES,t − Q L,t − K q ) − B −1 p Gδ t − I) Equation (A25) can be written as posing Y t = B −1 p (Q G,t + Q RES,t − Q L,t + K q ) − I: to reformulate Equation (A26), we can write, posing: and posing Q TN,t = −B −1 and thus, posing M = (B −1 p G) 2 , the overall objective function can be expressed as: Equation (A29) can be expressed in a vector form by using the following notation: Appendix A.3. Proof of Claim I Due to its linearity, the system state equation can be written separating the stochastic (z s t+1 ) and deterministic (z d t+1 ) components [45]: z t+1 = z s t+1 + z d t+1 = A z s t + z d t + B e s t + e d t + µ t +Ĉĥ t t = 0, . . . , T − 1 (A33) thus, the system can be so decomposed into a stochastic and deterministic subsystem. That is: z s t+1 = Az s t + Be s t +Ĉĥ t t = 0, . . . , T − 1 z s 0 = 0 (A34) z d t+1 = Az d t + Be d t + µ t t = 0, . . . , T − 1 Also, the cost function (Equation (A36)) can be decomposed into its stochastic and deterministic components represented by the functions J(z s , e s ) and J z d , e d , defined straightforwardly. This means that for each admissible solution (z, e), and the corresponding stochastic (z s , e s ) and deterministic z d , e d components, the following identity holds: J(z, e) = J(z s , e s ) + J z d , e d . (A37) We may assume that λ * t = P t z d1, * t − g t , and we obtain z d1, * t+1 = Az d1, * t − E t P t+1 z d1, * t+1 + E t g t+1 − M t (z d1, * t − r t ) t = 0, . . . , T − 1 (A76) z d1, * t+1 = (I + E t P t+1 ) −1 Az d1, * t + E t g t+1 − M t (z d1, * t − r t ) t = 0, . . . , T − 1 (A77) This equation must hold for all values of the state z * t , which in turn leads to the fact that the coefficient of z t and the rest of the terms in Equation (A77) must individually vanish. That is Finally, the closed-loop optimal control is e d, * k = −K t z d1 t + K g t (B g t+1 + F t r t ) t = 0, . . . , T − 1 (A81) with K t = (R t + B P t+1 B) −1 (B P t+1 A + F t ) t = 0, . . . , T − 1 (A82)