Robust Distributed Secondary Voltage Restoration Control of AC Microgrids under Multiple Communication Delays

: This paper focuses on the robust distributed secondary voltage restoration control of AC microgrids (MGs) under multiple communication delays and nonlinear model uncertainties. The problem is addressed in a multi-agent fashion where the generators’ local controllers play the role of cooperative agents communicating over a network and where electrical couplings among generators are interpreted as disturbances to be rejected. Communications are considered to be affected by heterogeneous network-induced time-varying delays with given upper-bounds and the MG is subjected to nonlinear model uncertainties and abrupt changes in the operating working condition. Robustness against uncertainties is achieved by means of an integral sliding mode control term embedded in the control protocol. Then, the global voltage restoration stability, despite the communication delays, is demonstrated through a Lyapunov-Krasovskii analysis. Given the delays’ bounds, and because the resulting stability conditions result in being non-convex with respect to the controller gain, then a relaxed linear matrix inequalities-based tuning criteria is developed to maximize the controller tuning, thus minimizing the restoration settling-time. By means of that, a criteria to estimate the maximal delay margin tolerated by the system is also provided. Finally, simulations on a faithful nonlinear MG model, showing the effectiveness of the proposed control strategy, are further discussed.


Introduction
Microgrids (MGs) are small-scale power systems consisting of localized grouping of Distributed Generators (DGs), storage systems and loads. In general, MGs operate either in islanded, autonomous mode or connected to the main power system. Recently the MG control system has been standardized into three layers [1]. The inner layer is called "Primary Control". It generates the actual command for the DG power converters' DC-AC interfaces. Commonly, it is implemented in a decentralized way and consists of the droop power control terms which aim to establish a desired power sharing among DGs and the inner voltage and current control loops. The "Secondary Control" layer enables the compensation of the frequency and voltage deviations introduced by the droop power control terms. This is because inverter-based DGs have no inertia [2]. Lastly, to optimize the DGs injected power from/to the main grid, a "Tertiary Control" layer is further designed. It allows us to adapt online the droop coefficients accordingly with the energy costs and the demand constraints imposed by the energy provider and vendors. It is worth mentioning that, among the three control layers, the secondary one is that of more practical interest to make compliant the integration of renewable resources within the AC power transmission paradigm. Moreover, due to the possibility of temporarily modifying the DGs' frequencies and output voltages to certain set-points while preserving the power sharing it would also be useful to perform the seamless and safe transition of the MG from islanded to the (i, j), are recounted in [19]. Therein, under the assumption of Markovian distributed delays and globally known set-points, a new stochastic small-signal MG model and a secondary controller is proposed. Small-signal analysis of an MG equipped with a time-delayed secondary voltage control is also studied in [20]. Therein, an analytical formula is developed to determine the individual delay margin for each link with respect to different sets of controller gains and leader pinning conditions. Lastly, we further mention [21], which proposes a time-delay tolerant distributed observer to track the average voltage of all DGs. Then, thanks to this estimation the implementation of lower level decentralized restoration strategies is enabled. Notice now that, although all the previously mentioned works account for the presence of different types of delays, all consist of a linear control protocol, thus they suffer the same limit of [7]. Indeed, they can only provide stability features in a local sense since their results are obtained by means of small signal analysis. Moreover, since they are linear, they cannot also face-off the presence of model non-linearities and/or disturbances.
Summarizing from the previous discussion it can be concluded that the control protocols that account robustness features against model uncertainties and disturbances are not able to face-off communication delays, whereas those designed to deal with the presence of communication delays lack robust control features because they are usually control protocols. Moreover, among the delay-oriented strategies, only [20] proposes a criterion for the secondary control tuning, whereas the other limits their analysis to provide necessary conditions for the closed loop stability within a local domain.

Statement of Contributions
Thus motivated, this paper proposes a robust voltage secondary control strategy that is able to perform the exact distributed tracking of the voltage set-point despite the absence of feedback-linearisation terms, and the presence of heterogeneous network-induced timevarying delays. Let us further note that the proposed control protocol is integral with discontinuous derivative thus well-suited to fed the inner control layers. It consists of two terms: (1) a distributed Integral SMC [22] component that aims to enforce each agent (DG) to behave as a reference unperturbed dynamic; (2) an ad-hoc designed distributed averaging control term aiming to globally, exponentially, achieve voltage restoration despite the non-uniform time-varying communication delays. A Lyapunov-Krasovskii analysis that used the simpler averaging problems, without [23], and with leaders [24], is also provided, and an LMI tuning criterion is established. Differently to the State-of-the-Art: (a) at the same time, both robust control concepts and the presence of multiple non-uniform time-varying communication delays are considered; (b) under the assumption that the power flowing within the MG is bounded, conditions for the robust global exponential stability of the system are established; (c) a new stability analysis for the MG's secondary control problem is proposed; (d) an off-line LMI-based algorithm aimed to find the maximal controller gain that, under the given delays bounds, guarantees the control objective is given along with criteria to estimate the corresponding maximal tolerated delay margin; (e) the robustification method in [22] is extended to the presence of delays among the agents' communications. Preliminary results of this research have been presented in [25]. However, although the controller structure here proposed is the same as in [25], this paper differs from that in the following aspects: • Here a much more faithful MG modelization expressed in the d-q-reference frame is considered, whereas in our past work DGs were simply modelled as second-order coupled oscillators. • The tuning criteria in our past work involves non-convex optimization which can thus hardly be solved in practice. On the other hand, our novel Algorithm 1 is based on the solution of a standard LMI problem, solvable by means of standard interior-point solvers.

•
Here Lemma 1 proofs that matrix Φ in Algorithm 1 is always negative definite, whereas in our past work this was an operating assumption; • The proposed control protocol is now tested on a realistic MG model which accounts, along with time-varying delays, 3-phase PWM controlled power bridges, noisy measurements with a realistic Signal-to-Noise Ratio, and unplanned faults.
• Define the gain decrement 0 < ∆k k , the edges weight (b) Solve for η > 0 and the positive definite matrices P, Q l , Q g , R g , W l , W g ∈ R 2N×2N the LMI problem (c) If the given solution satisfies (3) then: Break Else: k = k − ∆k EndIf.

Paper Organization
This manuscript is organized as follows-in Section 2 the nonlinear modelization of a MG for secondary control purposes is given. The voltage secondary control problem, along with the proposed robust control strategy and the main paper results are illustrated in Section 3. In Section 4, the robustness of the proposed controller is confirmed by means of simulations on a faithful 3-phase nonlinear MG model accounting IGBT powerconverters, time-varying communication delays, noisy measurements and uncertainties, and unplanned faults. Lastly, in Section 5, concluding remarks and possible hints for future research are provided.

Mathematical Notation
The set of real numbers is denoted by R. Let A = [a ij ] ∈ R n×n , with n > 0, its transpose is A , whereas a ij ∈ R denotes its entry in position (i, j). Assume A is symmetric, A 0 (A 0) denotes A positive (semi-)definite. Let A 2 and A F = trace{A A} be the Euclidean and Frobenius norms, then it holds that A 2 ≤ A F . Let B ∈ R n×n , A ⊗ B = [a ij B] ∈ R 2n×2n denote the Kronecker product between A and B, with i and j = 1, 2, . . . , n. The sign operator sign(a) is understood in the Filippov sense [22], such that sign(a) = 1 if a > 0, −1 if a < 0, otherwise it is set-valued as follows sign(a) ∈ [−1, 1]. I n denotes the n-dimensional identity matrix. 1 n and 0 n are the all 1, and all 0, n-dimensional column vectors. Let x ∈ R n , and A 0, and h > 0, the Jensen inequality for the integral terms Lemma 2.1 in [26] states that

Nomenclature
For the sake of clarity, the list of variables of the considered nonlinear model of MG, expressed in the d-q (direct-quadrature) reference frame, is provided below. The reader is referred to Figure 1 for a graphical explanation of the meaning of some variables. The considered nomenclature for the DG variables is as follows: ω com : Speed of the common rotating reference frame ω i : Speed of the rotating reference frame of the i-th DG δ i : Angle between the local and common reference frames υ ni , ω ni : Voltage and frequency secondary control commands; υ ki , i ki : 3-ph voltages and currents at node k of the i-th DG υ kdi , υ kqi : d-q voltages of the i-th DG at node k i kdi , i kqi : d-q currents of the i-th DG at node k k = l, o, b: Denote the input "k = l", output "k = o" and bus node "k = b"of a DG ω 0 , υ 0 : Desired voltage and frequency secondary restoration values P i , Q i : Active and reactive powers dc-components at the output node of the i-th DG υ * odi , υ * oqi : d-q voltage set-points of the i-th DG ψ di , ψ qi : d-q voltage error's integral of the i-th DG i * ldi , i * lqi : d-q current set-points of the i-th DG φ di , φ qi : d-q current error's integral of the i-th DG

Microgrid (MG) Modeling for Secondary Control Design
An MG is a cyber-physical system in which physical and software components are deeply intertwined. At the physical level, an MG is a geographically distributed power system consisting of DGs and loads connected by transmission lines. On the other hand, the supervisory and control algorithms, and in particular the secondary control layer, run over a communication infrastructure as shown in Figure 2. An inverter-based DG includes a 3-ph power converter in which the DC side is connected to a dc power source (e.g., photovoltaic panels [27], fuel cell system [28] or a wind turbine [8]), while its AC side is connected to the 3-phase power grid by the series of the coupling and the output filters, see Figure 1.
The Primary Control layer consists of three nested feedback control loops that aim to adjust the power, voltages, and currents injected within the main grid. The behaviour of a primary controlled DG is generally nonlinear. This is because the current and voltage control loops are expressed in the d-q local DG reference frame, rotating at the frequency ω i [29]. The Secondary Control layer aims at generating the commands/references for the inner Primary Control layer. In particular, through the signals ω ni and υ ni , the local DG frequencies ω i and the direct voltage components υ odi should be restored to their respective set-points ω re f and υ re f . In general, the phase angle δ i of each DG i is expressed with respect to a specifically selected DG considered as the common reference frame. Let ω com be the rotating speed of the common reference frame, then one has that The local frequency and voltage droop primary control characteristics are as follows where m i and n i represent the droop power coefficients, which aim to mimic the generator's inertia. Then, P i (t) and Q i (t) denote the DC components of active and reactive power flows at the DG outputs. A low-pass process with cut-off frequency ω ci ω i ω ni is designed to obtain P i (t) and Q i (t) as next Finally, the inner voltage and current primary control loops consists of cascaded PI controllers combined with feed-forward compensation terms. The voltage control is as follows whereas the inner current controller satisfies the next relations whereω is specified as the nominal angular frequency of the MG. Lastly, L f i and C f i denote the inductance and capacitance of the coupling filter drawn in Figure 1, and k pvi , k ivi , k f vi and k pci , k ici , k f ci denotes, respectively, the proportional, integral and feed-forward gains of those controllers. Finally, the voltage and current dynamics associated with the local LC filter and the output connector are where, in accordance with Figure 1, υ bdi (t), υ bqi (t) denote the d-q voltages at the connection bus. Let υ * ldi (t) and υ * lqi (t) in (20) and (21) express the smooth signals presented by inner current primary control loop, whereas υ ldi (t) and υ lqi (t) in (22) and (25) represent the effective non-smooth dispatched voltages of the inverter at the AC port after a high frequency PWM process. It is worth mentioning that, since the switching frequencies of power bridges are very high, the inverter dynamics can safely be ignored by comparing the MG dynamics. Hence, according to the averaging principle, (20) and (21) can be replaced for υ ldi and υ lqi into (22) and (23), more details can be found in [30,31].
Henceforth, the following reasonable assumption is assumed to hold.

Assumption 1.
Consider the primary controlled MG (8)- (21). Let a priori known bounds Π P and Π Q exist, such that, for each DG the active and reactive powers at the DG's output node "o" are assumed to be bounded by Remark 1. Assumption 1 is justified because the power flowing in the lines and/or absorbed by the load is bounded everywhere due to: (a) the passive behaviour of loads and lines; (b) the bounded operating range of the inverter voltages and currents due to their physical limits; (c) the presence of protection apparatus and the inner voltage and current primary control loops. This keeps the power flows within pre-specified ranges as discussed for instance in [11,30,32].
As discussed in [8,11,15,31], the cascaded voltage and current control loops in, resp., (14)- (17) and (18)-(21), are theoretically justified based on time-scale separation techniques. In particular because of the filter inductance L f i in (22), and because the inverse of the current integral gain 1/k ici in (18) can reasonably be considered small enough, then from (18), (20), (22) the following the Singular Perturbation analysis takes place where = max{L f i , 1/k ici } is the so-called perturbation parameter. Now according to the singular perturbation method and setting = 0, the quasi steady-state behaviour of the current fast dynamics can be obtained as Hence, at the quasi-steady mode the current control induces the inverter to send its set-point (16) almost immediately despite the inherent connection between the d and q component. Likewise one also obtains that i lqi (t) ≡ i * lqi (t). From (16), (24) and (29) the reduced-order relation between υ odi and the control input υ ni iṡ For later use, we express (30) in the next augmented form with whereυ n i is the actual voltage secondary control command anḋ plays as a matched unknown term accounting perturbation, parameter uncertainties, unmodelled dynamics, as well as chances in the MG working-point. Let us further note that due to Assumption 1, and the inherent boundedness of the MG variables enforced by the Primary Control, which makes the DG's frequency and voltages smooth, with bounded derivatives, then the existence of a sufficiently large constant W i such that (33) holds if fairly reasonable.

MG's Network Models
To achieve the Secondary Control objectives in a distributed way, DGs are assumed to be provided with computing and communication facilities embedded within their local controllers. Thus, each DG plays the role of autonomous agent over a networked cyberphysical system. Agents communicate over a network in which topology is encoded by a real-weighted graph G N = (V, E N ), where V = 1, 2, . . . , N is the set of local secondary controllers (or without loss of generality simply DGs), whereas E N is the set of the available communication links between DGs, and A = [α ij ] ∈ R N×N is the so-called adjacency matrix, where α ij = 1 if (i, j) ∈ E N , namely if DG i can receive data delivered by DG j, otherwise α ij = 0. Accordingly with the leader-follower paradigm, we further assume the existence of a root node, labelled by 0 in the augmented graph G N+1 (0 ∪ V, E N+1 ), which aims to provide the secondary control set-points to a non-empty set of DGs. Accordingly with Appendix A, where preliminary results on Graph Theory are provided, in the remainder, L and L N+1 denote the Laplacian matrices associated with G N and G N+1 . In the remainder, communications will be assumed to be affected by the network-induced time-varying delays. Further details on that will be given in the next section.

Voltage Secondary Control Design
Let the desired voltage set-point provided by the virtual node 0 be denoted by In the absence of a secondary control loop, and by assuming the voltage set-point globally available to each DG, thus by letting υ ni ≡ υ 0 , it can be noted that, due to the droop-power terms in (10) and (11) and the voltage drops in (30), at the resulting equilibrium υ di (t) will be smaller than υ 0 , thus restoration is needed. Our control objective can thus be summarized as follows where υ 0 is available only to a subset of DGs. To globally solve the secondary control problem (34) in a distributed way and by accounting network-induced delays in the communication links the next control protocol is proposeḋ where s i (t) is the desired integral sliding mode control manifold designed as follows and where m i > 0 and k > 0 are some constant gains to be designed, and α ij denotes whether DG j is enabled to communicate with DG i,

Remark 2.
It is worth remarking that the DG's output voltage derivatives are not available for measurements. However, they can easily be estimated and then used for output feedback purposes by means of differentiators implemented within each local controller. Following [33], and thanks to the finite-time convergence properties provided by the well-known Levant's exact differentiator, where if C i large enough, it results that, after a finite interval of time the estimationυ odi (t) coincides with the true valueυ odi (t). This is the reason why in (35)-(38) the quantitiesυ odi (t) are treated as if they were known.
Remark 3. According to Figure 2, the DGs dynamics are coupled with each other due to the interconnecting power lines and are also affected by the local utilities. Although (9)- (27) do not explicitly show this coupling, these effects can be encoded by the bus voltages υ bi [31]. Let us further note that, differently to [9,10], where these voltages were considered to be known and then compensated using feedback-linearization/feedforward approaches, in the present work they are supposed to be unknown. It is also worth mentioning that the considered model is strongly nonlinear because the primary control is designed in the d-q reference frame. Thus, due to its inherent nonlinear behaviour and because of, to the best of our knowledge, all the available time-delay tolerant secondary controllers are linear, then they can provide stable results only in the local sense. On the contrary, here, by means of the proposed robust control protocol, we can account for both the presence of communication delays and model nonlinearities/uncertainties without any model approximation, thus providing stability in the large sense.

Remark 4. The design (35)-(38)
is inspired by the robustification method proposed in [22] for arbitrary distributed consensus and optimization problems. Differently to that, here also the presence of communication delays is accounted for. The discontinuous term consists of an ISMC control term that aims to enforce each DG to behaves as a nominal multi-agent system, namely such thaẗ υ i = k · u i . Then, u i (t) is properly designed to fulfil (34) under the given communication topology and delays.
For later use, the collection of delays affecting the available communication links between the leader (node 0) and the generic DG j will be denoted by the set T (t) = {τ 1 , τ 2 , . . . , τ q }, where each element τ l (t) ∈ T (t) is such that Then, we indicate with S(t) = {σ 1 , σ 2 , . . . , σ m } the collection of delays affecting the oriented communications among DGs, and where σ g (t) ∈ S(t) is such that Note that the indexes m ≤ N(N − 1) and q ≤ N equal their maximum values only if G N+1 is a complete graph, and all the delays, for a given time t, are different. The common assumption of slowly-varying, bounded delays is made in [15,23,24,34,35]. Assumption 2. Let a priori known bounds τ l , σ g , δ l , δ g > 0 exist, such that, for each delay τ l ∈ T and σ g ∈ S, To correlate the data received from a DG for feedback purposes we further assume that: Assumption 3. Each communication delay τ ij (t) is detectable for all t ≥ 0, and (i, j) ∈ G N+1 .

Remark 5.
Communication protocols used to include the data-packet timestamp, thus the requirement of detectable delays is costless [36]. Once the delay is detected by means of local buffers each controller is enabled to retrieve its own state at that time and performs (38). These averaging laws are common in networked applications with multiple delays [15,23,24,37], and others.
The main results of the present paper are now presented: Theorem 1. Consider an MG of N DGs in which voltage dynamics satisfy (31), under the voltage control (35)- (38). Let the communication topology G N be connected and undirected. Let node 0 be a root node over the augmented graph G N+1 (0 ∪ V, E N+1 ). Let the communications subjected to heterogeneous time-varying delays τ ij (t) with given upper-bounds τ l , σ g , and δ l < 1 and δ g < 1 as in (43). Let Assumptions 1, 2, 3, be in force. Let m i > W i in (33). If the matrix inequality problem (1)-(6) admits a feasible solution with respect to k > 0, η > 0, and the symmetric positive definite matrices P, Q l , Q g , R g , W l , W g ∈ R 2N×2N , then the voltage restoration objective (34) can be globally uniformly asymptotically achieved.

Proof of Theorem 1. See the Appendix D.
Remark 6. The matrix inequality problem (1)- (6) is strongly non-linear with respect to k, thus standard convex optimization solvers will fail in finding their solutions because of the quadratic dependence from k of (5), and (6). However, for a given k, (1)-(6) degenerates into a set of LMIs, thus is solvable by means of standard convex optimization solvers, as for instance the well-known feasp MATLAB's LMI solver which is based on the Nesterov and Nemirovski's Polynomial Projective Method. The following Algorithm 1 is thus proposed to provide a method to find systematically, if it exists, the maximal feasible gain k for (35)- (38) such that, for the given delays' bounds, the fulfilment of the control objective is guaranteed.
Before presenting Algorithm 1 the next instrumental lemma is provided.

Proof. See Appendix C.
Corollary 1. Assume the condition of Theorem 1 to be satisfied. Hence a k = k exists such that, for the given delays' bounds the control objectives can be fulfilled. Then, Algorithm 1 will find such a feasible gain value k ; (b) LetĤ = ∑ m g=1 W g + ∑ q l=1 W l and k = k , an estimation of the maximum delay tolerated by the control system is given by Proof. Consider the LMI (5). Independently from the given k > 0 and η > 0, the delays' bounds, and the positive definite matrices W l and W g , there will always exist a sufficiently large positive definite matrix Q l 0 such that −(1 − δ l )Q l + η(A l (k) ∑ m g=1 σ g W g + ∑ q l=1 τ l W l A l (k)) ≺ 0, ∀ l. It follows that (5) is always feasible. The same can be concluded for (6) and ∀ g. Consider now the LMIs (3) and (4), which are coupled to each other by means of the matrices R g 0, and which play as dummy variables. Thus by substituting the right-hand side of (4) into (3), and by lettingQ = (∑ q l=1 Q l + ∑ m g=1 Q g + η · A 0 ∑ m g=1 σ g W g + ∑ q l=1 τ l W l A 0 ), the following LMI takes place: Φ P + PΦ +Q ≺ 0. Thus, because of Lemma 1, provided in Appendix B, Φ is Hurwitz, andQ is, by construction, positive definite, then a matrix P 0 can be found, satisfying the resulting Lyapunov equation Φ P + PΦ +Q = − I 2N , ∀ > 0. Thus, (3) is also always feasible. It follows that, if a feasible gain k exists, by means of Algorithm 1, the feasibility problem associated with (1)-(6) can be numerically studied iteratively for different k, until such a value k = k is found. If the actual k does not correspond with a feasible solution then k is decreased by ∆k > 0, and the routine is repeated. On the other hand, if k becomes smaller than or equal to zero, it means that for the given delays' bounds, ∆k andk, a feasible k does not exist. This confirms sentence (a).
Once a feasible solution for (1)-(6) is found, a lower estimation of the maximum delay tolerated τ by the system can be easily derived as follows: (1) Let us first replace such a solution on (1)-(6); (2) Then, by substituting τ (∑ m g=1 W g + ∑ q l=1 W l ) to ∑ m g=1 σ g W g + ∑ q l=1 τ l W l , a set of inequality in which the only unknown variable is τ takes place; (3) Finally, by solving for τ , by applying norms to the numerator and denominator and because · 2 ≤ · F , after manipulations, (44) is obtained.

Remark 7.
The achievement of the voltage restoration condition (34) through the proposed robust control protocol (35)- (38) does not affect the stability of the DGs' frequencies. In fact, because of the presence of the droop power control (9)-(10), which guarantees the MG synchronism and stability to load changes/variations, and because the achievement of (34) simply implies some smooth changes on P i (t) in accordance with (12), now υ odi → υ 0 . Then, the achievement of the voltage restoration will simply modify the actual frequency equilibria, which will correspond, following the analysis in [2], to the next quantity Let us further note that, because at each time t the total power in the MG is time-invariant, then the current's steady equilibria will change accordingly, as well. The results shown in Figure 3 corroborate these statements.

Test Rig Design
The proposed voltage control is tested on a realistic MG model developed by means of the Simscape Power Systems TM MATLAB/Simulink toolbox. Figure 2 depicts the schematics of the considered inverter-based MG. Each DG model comprises a 3-ph An IGBT bridge with 10 kVA of rated power provided by a 800 V dc-source. PWM Generators with 2 kHz carrier are employed to control the switching devices. In accordance with Equations (9)-(21) and (35)- (38), the primary and secondary controls are described in the local rotating reference frames by means of d-q Parks' Transformations. The dynamics of the LC filters (22)- (25), and of the output connectors (26) and (27), and the 3-ph power lines are instead formulated in the abc frame by using 3-ph Series RLC Branches. Loads consists of 3-ph Parallel RLC Loads. The current and voltage outputs of the inner voltage and current PIs, as well as the PWM modules, are saturated in accordance with the DGs' rated powers, resp., 380V ph−ph , 32A. The MG and DGs parameters are listed in Table 1. Lastly, to show the robustness of the proposed control rule to sudden unplanned events, a 3-ph to ground fault on Line 3 of Figure 2 will be scheduled throughout the simulation. Then, due to the surge of transient currents generated by the faults and according to the delays of protection devices, 10 ms later, two circuit breakers located at the branch buses of DG 3 and DG 4 will isolate the MG in two portions accordingly with the fault location.
The Runge-Kutta fixed-step solver is employed in the MATLAB environments to carry out the simulation with sampling time T s = 2 µs. The primary and secondary controllers have been discretized using a sampling stepT s = 500 µs. It is also worth mentioning that all the measurements used for feedback control are affected by realistic measurement noise. In particular, the MG currents and voltages measurements were first converted into the 4 ÷ 20 mA wired protocol with power transmission equal to 0.2 W, and were then corrupted by an Additive White Gaussian Noise with a realistic Signal-to-Noise Ratio of 90 dB [11].
The communication topology is chosen to be connected and undirected as in Figure 2. Only DG 1 can access the voltage set-point υ 0 from the virtual node 0. Communication delays are time-varying such that, for each oriented link (i, j) ∈ E N+1 ,τ ij (t) is randomly uniformly distributed within the open-interval (−1, 1). The delays' bounds in (43) are τ l = σ g = 0.03 s, and d l = d g = 0.999, ∀ l and g.

Voltage Secondary Control Implementation
Algorithm 1 has been implemented in the MATLAB environment, and the LMI problem (1)-(6) built by means of the lmiedit symbolic interface. The implementation of step (b) takes instead the advantage of the feasp LMI solver. Algorithm 1 is initialized withk = 5 and ∆k = 0.01. For the given delays' bounds and the communication topology G N+1 , it results in an optimal maximal gain k = k = 3.42, which corresponds to a delay stability margin of τ = 0.0781s. On the other hand, the discontinuous gain m i has been set equal to 140. It is also worth mentioning that to reduce the so-called chattering effect the sign operator in (35) is approximated as a sigmoidal function, namely sign(a) ≈ a/(|a| + ε) with ε = 0.01, and where a ∈ R denotes the operator argument. Finally the differentiators gains in (40) are chosen to be large enough. In particular, we set C i = 1000.

Case Study
The system is tested for 75s, and the delays are randomly initialized as next τ ij (0) ∈ [0.5τ , τ ]. The list of events scheduled throughout the test is displayed: • Step 1 (t = 0 s): Only the primary control is active with ω ni = 2π50 Hz, υ ni = 220 V RMS (per phase rms); • Step 2 (t = 5 s): The voltage control with υ ni as in (35)- (38) and υ 0 = 220 V RMS is activated; • Step 3 (t = 15/25 s): An additional load (P L3 , Q L3 ) is connected/disconnected at the bus 3; • Step 4 (t = 35 s): A 3-ph to ground fault occurs on the Line 3; • Step 5 (t = 35.1 s): Over-current protection devices isolate Line 3, thus DG4 results electrically isolated; • Step 6 (t = 45 s): The set-point for the voltage Secondary Control is changed to υ 0 = 225 V RMS ; Consider Figure 4a. In the first five seconds, when only the primary control is active, all the DG voltages are less than the reference value of 220 V RMS . It is evident that restoration is needed. Once at t = 5 s the proposed voltage control is activated, promptly the DGs' output voltages are globally, asymptotically, restored to the desired setpoint, despite the communication delays. Then, the MG operating working points change after the connection/disconnection of the additional load (P L3 , Q L3 ) by means of a 3-ph breaker. This verifies that the proposed controller is robust against unexpected changes of demands. Moreover, at t = 35 s also a 3-ph to ground fault is triggered on the line between DG 3 and DG 4. Then, 10 ms later that line is isolated and the MG is sectioned into two sub-MGs, one consisting of DG 1, DG 2, and DG 3 and the respective loads, and another composed only of DG 4 and the local load (P L4 , Q L4 ). From Figure 4a it can be observed that although the occurrence of such a critical event at t = 35 s, and the fact that now DG 4 is isolated, all the DG's voltages remain at the desired nominal setpoint v ref = 220 V RMS . Then, at t = 45 s such a setpoint is changed to υ 0 = 225 V RMS . Consequently, all the DG's voltages converge to it. Finally, the time evolution of the voltage secondary control signals is shown in Figure 4b. The results justify that the proposed scheme provides a satisfactory performance and smooth control signals. Following what was stated in Remark 7, and with the aim of showing that the achievement of the voltage restoration by means of the proposed controller does not affect the MG stability, in Figure 3a the DGs' frequencies measured by means of 3-phase Phase-Locked Loop sensors deployed at the DGs' output nodes are drawn. Then, for the sake of completeness, the temporal response of the DGs' RMS currents is also further provided in Figure 3b.

Conclusions
Most of the existing time-delay tolerant MG control protocols are linear. Thus, because an MG is a strongly nonlinear system, all those approaches can only provide stability features in the local sense under additional model approximation or linearization procedures. Our approach instead improves the current State-of-the-Art by providing stability features in the large sense by means of a co-design control method, which inherits both the robustness features of the sliding mode control paradigm and the resilience of the LMI-based controller against network induced time-varying communication delays. Future research activities will be devoted to removing the assumption of slowly varying delay and will consider more general communication topologies. Moreover, extensions of the proposed method to the frequency restoration problem are currently under investigation.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Graph Theory
A directed graph G N (V, E N ) is a topological tool used to describe pairwise interactions between objects, abstractly referred as nodes, or agents in the multi-agent research area. The node set is V = {1, 2, . . . , N}, whereas the edges set E N ⊆ {V × V } describe the node's pairwise interactions. Let (i, j) ∈ E N , we refer to i; and j as the tail and head of the edge. The adjacency matrix of G N is A = [α ij ] ∈ R N×N , and α ij denotes the edge weight. If (i, j) / ∈ E then α ij = 0, otherwise α ij = 1. A directed path is an alternating sequence of nodes and edges with both endpoints of an edge appearing adjacent to it in the sequence. Node i is a root node if it can reached from any other vertex by traversing a directed path. If there exists a root node then G N is said to have a spanning tree. Let B = diag{b 1 , . . . , b n }, with b i = ∑ N j=1 α ij , ∀ i ∈ V, the Laplacian matrix of G N is L = B − A. Since by construction L1 N = 0 N , then 1 N is an eigenvector of L associated with a 0 eigenvalue. If (i, j) ∈ E =⇒ (j, i) ∈ E , G N is undirected, and it results that L = L . An undirected graph is connected if it has not disconnected nodes.

Appendix C. Proof of Lemma 1
Since node 0 is a root node over the graph G N+1 , from Proposition 1, its corresponding Laplacian matrix has a simple 0 eigenvalue, whereas the rest N − 1 have all positive real part, and L is the Laplacian matrix associated with the reduced graph G N (V, E N ). It is easy note that F is diagonally dominant, namely | f ii | ≥ ∑ j =i f ij ∀ i, and symmetric. Thus, the eigenvalues of F are real and non-negative. Let us now define x = 1 N x 0 , with x 0 ∈ R and note that the following continued equality holds Thus, because of F is invertible, we can conclude that −F is Hurwitz. Let us now consider the matrices A l (k), andÂ g (k) in (1), and matrix Φ(k) in (2). By inspection it results that Let D be the matrix which contains the generalized eigenvectors of F as columns, it results that where λ i > 0. Let us now note that the Jordan Canonical representation associated with matrix (A6) satisfies Consider, now the eigenpolynomial of Φ, we have From (A9), we can immediately conclude that the eigenvalues of Φ(k) will be strictly negative for all λ i , and k > 0. This conclude the proof of Lemma 1.

Appendix D. Proof of Theorem 1
Let us consider the sliding manifold (36), provided below for the clarities' sake By substituting the proposed voltage secondary control (35)- (38) into the DG's voltage dynamics (31), one has For the discontinuous differential Equations (A11) and (A12), their solutions can be understood in the Filippov sense. Namely, as the solution of an appropriate differential inclusion, the existence of which is guaranteed and for which its absolute continuity is in force. The reader is referred to [22] for a comprehensive account of the necessary notions of non-smooth analysis. Particularly, following Theorem 3 in [22], namely by differentiating the locally Lipschitz Lyapunov functional V = ∑ N i=1 |s i (t)|, it results that, except for Lebesgue measure zero points, which can be disregarded, it yieldṡ Thus, thanks to (33), and m i > Γ P i , it results that condition s i =ṡ i = 0 is time-invariant for all t ≥ 0. By substituting u i (t) in (38), and s i =ṡ i = 0, into (A13), the following secondo-order closed-loop local dynamic takes place The point now is to find which conditions (A14) have to meet to asymptotically achieve the control objective (34). Let us first define the following voltage restoration error vector e(t) = e 1 (t), e 2 (t), . . . , e N (t) with e i (t) = υ odi (t) − υ 0 υ odi (t) −υ 0 (A15) Consider now (1), by differentiating e i (t), after some algebraic manipulations, one deriveṡ Now, in order to provide a compact state-space representation of the networked error dynamics associated with the vectors (A15), and accordingly with the notation introduced for the delays in (41) and (42), where q = card{T (t)} and m = card{S (t)}, it resultṡ For stability analysis purposes, and on the basis of the so-called Leibniz-Newton formula we introduce the following transformations Hence, from (A18), the networked closed-loop dynamic in (A17) can be recast as nexṫ where P, Q l , Q g , W l , W g ∈ R 2N×2N are constant, symmetric, and positive definite matrices to be determined and η is a positive scalar. Let τ = max l,g {τ l , σ g } > 0, and let α(e(t)) = e(t) Pe(t), β(e(t)) = e(t) Pe(t) + e(t − σ g (t)) Q g e(t − σ g (t))(1 −σ g (t)) Because of Assumption 2, and by adding to the right-hand side of (A24) the next identically zero quantity m ∑ g=1 σ g e(t) R g e(t) − m ∑ g=1 σ g e(t) R g e(t) = 0, where R g 0 is a matrix to be determined, (A24) can be upper-estimated as nexṫ and after substituting (A26), and (A28), the inequality (A25) can further be upper-estimated as followsV (t) ≤ρ(t) Σρ(t) + ηė(t) Hė(t) − q ∑ l=1 e(t − τ l (t)) Q l e(t − τ l (t))(1 − d l ) Define now the following augmented state vector ξ(t) = e(t) , e(t − τ 1 (t)) , . . . , e(t − τ q (t)) , e(t − σ 1 (t)) , . . . , e(t − σ m (t)) , then, by substituting (A17) into the second term of (A29), and after lengthy manipulations, (A29) can finally be recast as nexṫ where Θ is an upper triangular block matrix in the from and which diagonal blocks take the following form η Q l + A l H A l , l = 1, 2, . . . , q, η Q g +Â g HÂ g , g = 1, 2, . . . , m.
Hence, if the LMIs (3)-(6) are satisfied, then the matrices Σ in (A30) and Θ in (A33) are negative definite. Therefore, it results thatV(t) < 0 and thus also condition (A3) of Theorem 2 is satisfied. Moreover, by choosing α(s) as in (A21), it follows that lim s→∞ α(s) = +∞, and hence the error vector e(t) globally uniformly converges to zero, which it further implies the voltage restoration achievement. This concludes the proof.