Robust Two-Layer Control of DC Microgrids With Fluctuating Constant Power Load Demand

This article proposes a cascaded control structure for the regulation of dc microgrids (MGs). It is well known that the negative impedance characteristics of a constant power load (CPL) adversely affect the stability of the network and can cause problems, such as voltage collapse or damage the electronic components. To mitigate this, we propose a two-layer control structure, where at the inner layer, the proposed controller achieves fast tracking of the supplied reference points and ultimate boundedness of the trajectories in a desired set. The outer layer generates the inner-layer reference points, accounts for system constraints, and introduces robustness of the voltage dynamics to unknown perturbations of the CPL demand. For the first time, an investigation of the nonlinear geometric behavior of the CPL is carried out to derive necessary conditions that ensure the boundedness of the network dynamics and feasible regulation to a desired equilibrium set. Finally, control Lyapunov functions are formulated to prove the stability and estimate the region of attraction of the closed-loop dynamics. A simulated scenario of a meshed MG network is presented to confirm the validity of the results.


Robust Two-Layer Control of DC Microgrids With
Fluctuating Constant Power Load Demand Grigoris Michos , Pablo Rodolfo Baldivieso-Monasterios , George C. Konstantopoulos , Member, IEEE, and Paul A. Trodden , Member, IEEE Abstract-This article proposes a cascaded control structure for the regulation of dc microgrids (MGs).It is well known that the negative impedance characteristics of a constant power load (CPL) adversely affect the stability of the network and can cause problems, such as voltage collapse or damage the electronic components.To mitigate this, we propose a two-layer control structure, where at the inner layer, the proposed controller achieves fast tracking of the supplied reference points and ultimate boundedness of the trajectories in a desired set.The outer layer generates the inner-layer reference points, accounts for system constraints, and introduces robustness of the voltage dynamics to unknown perturbations of the CPL demand.For the first time, an investigation of the nonlinear geometric behavior of the CPL is carried out to derive necessary conditions that ensure the boundedness of the network dynamics and feasible regulation to a desired equilibrium set.Finally, control Lyapunov functions are formulated to prove the stability and estimate the region of attraction of the closed-loop dynamics.A simulated scenario of a meshed MG network is presented to confirm the validity of the results.Index Terms-Constraint control, CPL, microgrids, nonlinear systems, robustness.

I. INTRODUCTION
T HE concept of a microgrid (MG) brought a paradigm change to the architecture of conventional power networks [1].The traditionally centralized structure has shifted to geographically decentralized clusters that are able to operate both in a grid-connected setting and isolated, known as an islanded mode.The MG structure can be found in both ac and dc applications; however, in many cases, the use of a dc structure is often preferred because it provides higher efficiency and Grigoris Michos, Pablo Rodolfo Baldivieso-Monasterios, and Paul A. Trodden are with the Department of Automatic Control and Systems Engineering, University of Sheffield, S1 3JD Sheffield, U.K. (e-mail: gmichos1@sheffield.ac.uk; p.baldivieso@sheffield.ac.uk; p.trodden@ sheffield.ac.uk).
George C. Konstantopoulos is with the Department of Electrical and Computer Engineering, University of Patras, 26500 Rion, Greece (email: g.konstantopoulos@ece.upatras.gr).
At the heart of the dc MG lies the bidirectional dc/dc converter, achieving integration of renewable energy sources in the network.These devices are utilized to achieve MG voltage/current regulation and power flow control and ensure that normal network operation is retained in the presence of external disturbances.One of the main challenges of dc/dc converter control is the instability caused by the nonlinearities of constant power loads (CPLs), [3], [4].CPLs introduce negative impedance characteristics to the dynamics, which result in unstable equilibrium points.A plethora of studies has tried to address this problem, with the majority adopting a linear control approach; see, for example, [5], [6], or [7].The main drawback of linear approaches is the use of small-signal model analysis, which only guarantees local stability in an area around a desired operating point, where the linearization remains valid.In the presence of large load demand fluctuations, the nonlinearity of the CPL can drive the system outside the linearization region and destabilize the system.In light of this issue, more sophisticated control techniques have been developed that include the nonlinearity inside the control formulation.In [8], an adaptive backstepping technique was proposed for MGs feeding CPLs, where the estimation of the load demand is carried through a Kalman filter.A sliding-mode controller was presented in [9], where the sliding surface is designed as a linear combination of the voltage and current tracking errors.The stability of a boost converter feeding a CPL is demonstrated and the control scheme regulates the voltage to a desired reference point.Other methods using passivity theory [10], [11] or output regulation [12] have also been proposed.Enhancing the robustness of dc networks comprising buck converters to unknown CPL demand is also a topic of rising interest.Passivity theory was used in [13] and [14] to establish robust stability of the network equilibria to perturbations of the load demand.This work was extended in [15] to achieve global robust stability in the case of constantimpedance-current loads.It has been shown that the adoption of a more advanced technique can enlarge the region of attraction and provide stronger stability guarantees with respect to changes in the load demand.However, the majority of the proposed techniques do not take system constraints into consideration while the behavior induced by the CPL is not extensively studied.A question that arises is whether the behavior of the CPL can be included in the design of a unified constrained MG controller that guarantees stability and boundedness of the voltage and current states in the presence of unknown load demand perturbations.
The need to adopt a control scheme that is able to satisfy constraint sets arises naturally in MG control, most commonly in the form of actuation limits, current and voltage network capacities, or in order to prevent damage to electronic components.Therefore, the system is required to operate within a predefined operating range, which can be translated to the form of constraints on both the current and the voltage of each dc/dc converter.One of the most popular and effective techniques for achieving control in the presence of constraints is model predictive control (MPC).There is rich literature behind robust MPC approaches, the majority of which follow a Tube MPC formulation.The term "Tube MPC" refers to a collection of control approaches that bound the trajectories of the uncertain system within a sequence of sets and regulate this sequence to desired terminal sets [16].One of the most famous approaches dates back to the work of Mayne et al. [17], which standardized the Tube MPC for linear systems.First, the uncertain system is decomposed into nominal and error dynamics; then, an approximation of the minimum robust positive invariant (mRPI) set is calculated to bound the error while the uncertain system is driven by regulating the nominal state trajectory in conjunction with feedback control on the error dynamics.The calculation of the mRPI set approximation requires an explicit form of the integral flow of the system, i.e., the solution of the ordinary differential equation describing the system dynamics [18].Calculating this would present a challenging task in the nonlinear setting, as an analytical solution to the dynamics may not even exist.A few approaches have focused on the nonlinear case, for example, in [19], a linearization around each point in the horizon was proposed while a feedback linearization was utilized in [20].A tube nonlinear MPC (NMPC) was proposed in [21]; however, the calculation of the "restricted" nominal constraints is carried through simulations and lacks an analytical approach.The contraction theory is used to construct the tubes in [22] for a design of a distributed NMPC for dynamically decoupled subsystems.A method to optimize the tube size online was proposed in [23] by exploiting the structure of the adopted boundary layer sliding controller.An approach that constructs positive invariance sets for globally Lipschitz systems was proposed in [24], where the control action relies on computing a quadratic Lyapunov function for the system.However, imposing a globally Lipschitz condition restricts the scope of possible applications of this control method and may result in a conservative controller.It is evident from the literature that a robust NMPC scheme is a subject open for investigation and is often reliant on the specific form of the dynamics.According to the authors' best knowledge, despite its distinct advantages, the extension of such a scheme to MGs is still in its early stages.In a previous work of the authors, a robust NMPC scheme for islanded dc MGs was developed in [25], where the inherent robustness properties of the nominal MPC are exploited to ensure the recursive feasibility of the optimal control problem.It is shown that the system demonstrates a degree of robustness, by bounding the fluctuations of the power demand among the sampling intervals.However, only linear loads were considered in this study.This work was extended to include CPLs in the network in [26], where tools from economic MPC were used to guarantee the recursive feasibility of the proposed control scheme.
In this study, we propose a unified approach for the MG regulation, i.e., we include both a current and a voltage controller.We extend our previous work by relaxing the bound on the load demand fluctuations and investigating its effect on the rest of the network.More specifically, we propose a robust NMPC scheme for regulating the voltage of the network in the presence of persisting, unknown, power demand fluctuations.We exploit the dynamic structure of the converter system in order to parameterize a control law that bounds the system dynamics in a desired region.As will be shown, this results in a Tube MPC performance but in a nonlinear setting, where the bound on the error is chosen based on the interplay of the control and drift vector fields.A preliminary version of this was recently proposed by the authors in [27], where local approximation techniques were used to show local input to state stability of the dynamics and establish a tube-like behavior of the system trajectories.This article significantly extends this work by employing candidate Lyapunov functions of the original nonlinear dynamics as opposed to using linearization techniques.In addition, we study the geometric effect of the CPL on the network dynamics to characterize a positive invariant set under the solution of the voltage dynamics.Contrary to our previous work, this is done in a global sense, by studying the behavior of the network equilibria.Thus, we avoid possible approximation errors due to linearization techniques and reduce the conservativeness of the controller.More specifically, the main contributions of this work are as follows.
1) C1: We propose a modified version of the state-limiting PI, introduced in [28], in order to bound the converter current and introduce overcurrent converter protection during transient performance.Contrary to the literature, instead of a saturation unit, we employ a nonlinear parameterization of the input and show that by the proposed control law, the system flow is smooth and ultimately bounded in a predefined set, thus, facilitating the stability analysis of the overall system.2) C2: For the first time, we investigate the geometric behavior of the CPL and exploit the effect it has on the network to construct a robust positive invariant set of the system dynamics.Contrary to the work in [27], this set is constructed using the displacement of the voltage equilibria due to the power demand instead of employing local linearizations.Then, we employ energy-like functions to establish ultimate boundedness and specify an explicit bound of the error trajectories, as opposed to relying on local approximations of the system's region of attraction.This way, we can characterize a positive invariant set in a global sense, provide stronger theoretical guarantees, and reduce conservativeness by allowing larger unknown load demand perturbations.3) C3: We propose the adoption of an MPC control scheme to drive the nonlinear nominal dynamics to desired reference voltages.We expand the results from the work in [27], by deducing a necessary condition of the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
augmented nominal control law that guarantees stability of the cascaded dynamics and recursive feasibility of the optimal control problem.Furthermore, we include a result on the estimation of the region of attraction for the nominal dynamics.4) C4: We introduce a control Lyapunov function to prove the stability of the inner current dynamics and estimate the region of attraction of the respective equilibrium points.Then, we prove the stability of the cascaded dynamics and establish stronger stability guarantees compared to the linearization approaches, as these do not require explicit knowledge of the load demand and are not affected by linearization errors.The rest of this article is organized as follows.In Section II, we formulate the problem of controlling the dc MG network.In Section III, we introduce the current controller and prove the ultimate boundedness of the current dynamics.In Section IV, we formulate the voltage control that guarantees constraint satisfaction and robustness to constant power demand perturbations.Then, Section V provides the stability results of the cascaded structure, and in Section VI, we demonstrate the application of the proposed control scheme with a simulation scenario.Finally, Section VII concludes this article.

A. Notation
The notation |A| denotes the cardinality of A. The closure, boundary, and interior of a set A are denoted as cl(A), ∂A, and int(A), respectively.The Pontryagin difference of polytopes A and B is defined as A B = {a ∈ R n : a + b ∈ A ∀b ∈ B}.For a vector x ∈ R n , the notation [x] denotes a diagonal matrix on R n×n , where A topological space M is called a smooth manifold if for any p ∈ M , there exist an open set in a neighborhood of p, O ⊂ M and diffeomorphism ψ : O → R n .The combination (ψ, O) is called a coordinate chart of M .The tangent space on a point p ∈ M is denoted as T p M , and a vector field X is a map assigning a vector X(p) ∈ T p M at each point p ∈ M .A Riemannian metric G represents the inner product of the tangent space of M at x, where G : An MG can be seen as an undirected and connected graph G = (M, E) where the set of nodes M represents a collection of power converters and local loads; the set of edges E ⊆ M × M  defining the MG topology is characterized by the node-edge matrix B ∈ R |E|×|M| , which for edge ε = (i, j) ∈ E involving nodes i and j can be defined as [B] ei = 1 if node i is the source of e ∈ E, and [B] ej = −1 if node j is its sink, and zero otherwise.The weights of the edge ε, representing the admittances of the lines, are collected in the line admittances matrix R −1 .
We also use the following well-known result that we include here for reference.
Theorem 1 (Nagumo's theorem): Let ẋ = f (x), where the map f is at least once continuously differentiable and the solutions exist inside an open set O ⊆ R n .Then, the closed subset S ⊂ O is a positive invariant under the flow of the system if and only if for all x ∈ ∂S, where d(•, •) denotes the Euclidean distance.

II. PROBLEM FORMULATION
In this section, we investigate the system of an islanded meshed dc MG composed of n number of power converters, where each ith converter is connected to local CPL, with i ∈ M = [1, 2, . . ., n].Our aim is to introduce robustness of the system operation to load demand fluctuations and ensure that the system dynamics are restricted in a predefined operational range.It is common to represent an islanded dc MG as a meshed, connected, undirected graph, see Fig. 1.Each node models an interlinking dc/dc buck converter, integrating an energy source to the power network.The power consumption of the network is modeled by CPLs, connected to the output capacitor of each converter.This is represented in the converter circuit diagram as a controlled current source, see Fig. 2. In order to simplify the analysis, we only consider the pure resistive component of the lines.This is a common approach in dc MGs as it has been shown that the same stability results can be obtained by omitting the line inductance, see, for example, [29].Therefore, we will consider the Kirchhoff model of the network dynamics given by The capacitance and inductance of each node are collected into diagonal matrices L and C, respectively.r denotes a diagonal matrix collecting the node parasitic resistances, and E denotes the input voltage of the converter.The states of the system are the capacitor voltage v ∈ R n and the inductor current i ∈ R n while the duty ratio of the switching mechanism ν ∈ [0, 1] is regarded as the input to the system.
The connection between the nodes of the graphs can be represented by the weighted adjacency matrix A(R) ∈ R n×n , where a ij = R −1 ij , with R −1 ij the admittance of the line between nodes i and j, and a ij = 0 if the edge (i, j) is not incident.The full topology of the network is represented by the Laplacian matrix L = [A(R)1 n ] − A(R).Therefore, the output network current can be modeled as where P ∈ R n is the power demand vector.To allow the development of a robust control strategy, we assume that the load demands are bounded in a set P ⊂ R n .In addition, the system is subjected to operational constraints, i.e., desired regions of voltage and current operation, denoted as X and I, respectively.In order to ensure a smooth region of operation, we invoke the following assumption on the structure of the constraint sets.Assumption 1: The constraint sets X, I, and load demand P are polytopic and compact sets, where the origin is within the interior of P and I.
It is common to require the inner current and output voltage of a converter to operate in a predefined range, usually as a bounded deviation from a rated value.The above assumption translates this requirement into a mathematical notion, and ensures that the load demand acquires a finite maximum value.This property is used to ensure the convexity of the optimization problem that is introduced in Section IV.In addition, a consequence of this assumption is that the origin is not necessarily included in the voltage constraint set.Thus, the voltage dynamics are smooth over the constraint set, which allows us to consider the respective state subspace as a Riemannian manifold M embedded in R n .The compactness of the constraint set also implies that any closed subset of X is also compact.As a result, any subspace of the metric space (X, d R ) has the Heine-Borel property, i.e., the metric space is compact and complete.This is a particularly useful property that will assist us in the following sections to establish convergence of the system flow.
Linearity of (2) with respect to the load allows us to separate the load current into a nominal current depending on an a-priori known nominal load P and an uncertain one parameterized by a deviation δP = P − P .Then, we define the set of deviations from the nominal load demand as a disturbance set W := {δP ∈ R n : P + δP ∈ P }.Following Assumption 1, the disturbance set W inherits the same properties of P .The control objective is to introduce robustness to the system with respect to fluctuations of the load demand δP from the known, constant value P .To this aim, in the following sections, we first introduce the boundedness of each converter's inner current state i i .Then, we decompose the true voltage dynamics into a nominal state, where the load demand is P at all times, and an error between the true and the nominal voltage emanating from the fluctuations of the load demand δP from the nominal value P .Finally, we introduce conditions on the choice of the tuning parameters such that the distance between the trajectory of the nominal and the true voltages is bounded, the network satisfies the desired operational constraints at all times, and the stability of the overall system is guaranteed.

III. INNER CURRENT CONTROL
The current controller regulates the input current to a desired setpoint by adjusting the duty ratio ν.It is common in the literature to use some form of a saturated controller in order to achieve overcurrent protection of each power converter.To this aim, we employ a modified version of the state-limiting PI, first introduced in [28].This ensures smoothness of the dynamics, and thus avoiding performance degradation or instabilities caused by traditional saturated controllers, as highlighted in [30].Moreover, it allows us to formulate an analytic procedure to choose the control parameters that achieve the desired transient performance and boundedness properties.The parameterization of the duty ratio ν is given by resulting in closed-loop current dynamics where î is the reference current, M, k p , k I ∈ R n×n are the tuning parameters, and σ is the controller introduced integrator state.As shown in the following proposition, this control parameterization restricts the flow of the closed-loop system in a control invariant subset of the constraint set k p,i is the control invariant for the current dynamics (4).
, the derivative vanishes at point |σ i | = 1, i.e., this point represents an equilibrium point of the integrator dynamics.Leveraging on the continuity property of ( 4) and the form of E σ,i , the trajectory Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
cannot escape the ball defined by |σ i | = 1, which leads to a contradiction in the existence of a point σ i (τ 2 ) outside this ball.Now, we proceed with the boundedness of the current dynamics.We consider the C 1 energy function of the inductor E c,i = 1 2 Li 2 i , with time derivative Following [31,Th. 4.18], the flow of the system φ c (i 0 i , t) for i 0 i ≤ M i k p,i is ultimately bounded with bound I max,i = M i k p,i .Therefore, the desired behavior of the inner current loop can be achieved by adjusting the tuning parameters according to the desired bound I max,i = M i k p,i .The stability of the system will be investigated in Section V, where we will discuss the stability properties of the overall cascaded dynamics.In addition, we note that the choice of the inner current bound needs to be made according to the required voltage constraint set X, in order to satisfy the requirement ν i ∈ (0, 1) for all (i i , σ i ) ∈ C i and v i ∈ X i .

IV. VOLTAGE CONTROL
In this section, we formulate the robust control scheme for the voltage node dynamics.To this aim, we separate the dynamics into a nominal and an uncertain part.In the following, we will show that the uncertain part of the dynamics can be bounded within a positive invariant set, leveraging on the geometric properties of the closed-loop vector field.Then, the nominal dynamics are shown to act as a driving subsystem, regulating the true system to some desired reference behavior.We will use the notion of control sets to guarantee constraint satisfaction of the nominal subsystem, which, by an appropriate parameterization of the constraint sets, implies constraint satisfaction of Y = X × C by the uncertain dynamics.In this study, we assume that the current operates at a faster time-scale than the voltage.This is a common assumption in the literature, see, for example, [32] or [33], as it allows us to study each of the dynamic components separately and, thus, significantly simplifies the analysis.This leads to a cascaded structure where the input current in (1b) represents a piecewise constant reference and is considered as the input to the voltage dynamics.Following [34,Th. 4.4] and [35, Appendix A], the reachability, and consequently convergence, properties are preserved under the replacement of integrable smooth controls with piecewise constant controls.Having established the boundedness of the solution of (4), we formalize the time-scale separation by invoking the following assumption.
Assumption 2 (Time-scale separation): The network parameters satisfy Remark 1: As a consequence of Assumption 1, the load demand P i is bounded and, thus, one can enforce the time-scale separation of the node dynamics by an appropriate selection of the tuning parameters, as shown in Assumption 2 where we provide explicit time constants derived using singular perturbation analysis.We refer the interested reader to [31,Ch. 11] for a detailed analysis.

A. Formulation and Boundedness of the Error Dynamics
Following Assumption 2, we begin by writing the node voltage uncertain system and define the error term e = v − z, which separates the nominal part of the dynamics from the uncertain.Our aim is to use the nominal state z as a driving state of the system's flow φ v (•) while restricting the effect of the disturbance such that the distance between the flow of the uncertain system and the nominal satisfy d R (φ v , φ z ) ≤ ε for some ε > 0. We define the control policy where u is the nominal control policy that will be defined later.Then, isolating the nominal part of ( 5), results in Therefore, the evolution of the error dynamics can be described by the ordinary differential equation where for a constant power demand deviation δP s and converged nominal dynamics to a point ẑ, we have Many studies have been focused on the existence of real solutions for the voltage equilibrium map usually deriving a necessary condition to be satisfied by the system parameters, see, for example, [36].It is shown that, if the necessary conditions are met, the equilibrium map is a diffeomorphism where a highand a low-voltage solution exist, see [37] or [38].Note that γ(ê) = 0 is the result of the state transformation v = e + z, similar results with the voltage v equilibrium map can be obtained for ê = γ −1 (0).Therefore, the discontinuity of the true voltage dynamics is translated to the critical point e = −z of the error state space.The second solution of γ(ê) = 0 is then the steady-state displacement caused by the load demand perturbations δP between the true voltage trajectory and the nominal.
In the following, we intend to use the imposed properties of the vector field γ(•) in an area around the origin to restrict the error flow φ e in a positive invariant set.Considering the solution of ê = γ −1 (0) closest to the origin, the maximum displacement due to the load demand can be found by Next, we define the polytope where P = [I n − I n ] and ξ(e m ) = [1 n e m 1 n e m ] , that is, the set N is an n-dimensional "box" around the origin.The question that arises is under which conditions does the system guarantee the desired theoretic properties, i.e., the existence of a unique solution of the error dynamics in N , as well as positive invariance of N .The former is addressed in the following result.Lemma 1 (Lipschitz continuity): Consider the error dynamics given in (8) and bounded positive nominal voltage z ≤ z ≤ z.Given a constant power load (CPL) demand P s ∈ { P } ⊕ W and a scalar positive definite matrix K, the map γ(•) is Lipschitz continuous on N , with Lipschitz constant Proof: We begin by finding the Lie derivative for a constant power demand Due to the fact that z and e m are positive values, the matrix 2 − norm is given by We define and note that the resulting matrix within the norm is symmetric.Then, for a symmetric matrix, the largest singular value is given by its spectral radius, i.e., By the properties of the Laplacian matrix, we have λ n ≥ • • • ≥ λ 2 ≥ λ 1 = 0, and thus, λ max (−L) = 0.The above takes the form This establishes the desired upper bound on the norm of the Lie derivative, i.e., Lie(γ) 2 ≤ L. Therefore, for two infinitesimally different e 1 , e 2 ∈ N and a ξ ∈ (e 1 , e 2 ), an application of the mean value theorem states Therefore, the vector field is Lipschitz continuous on N with Lipschitz constant K.
Having established the Lipschitz continuity of the system flow inside N , we now are able to investigate the stability properties of the equilibrium ê ∈ N .This is demonstrated in the following result.

Proposition 2 (Stability and positive invariance of the error dynamics):
For bounded local nominal dynamics z i ≤ z i ≤ z i , the error dynamics of every node i ∈ N , admit an asymptotically stable equilibrium point in N , if the local feedback gain satisfies with P m i = max(| Pi + δP i |) and δP i ∈ W i .In addition, the set N is a positive invariant under the flow of the system φ e (e 0 , t), for all t > 0 and e 0 ∈ N .
Proof: For any equilibrium point ê ∈ N and respective constant load demand P , the resulting Jacobian matrix takes the form In order for the equilibrium point ê to be stable, the respective Jacobian matrix J needs to be Hurwitz, i.e., to have negative eigenvalues.Using Lyapunov theory, it suffices to investigate the largest eigenvalue of the matrix where we have used the fact that λ max (−L) = 0.In order to ensure that a similar condition holds for any ê ∈ N and any P ∈ P , we note that holds for all ê ∈ N .Therefore, we require ) 2 where we have used the fact that z i ≤ z i ≤ z i , where the lower bound is positive.Following the definition of N , for some point at ∂N , we have that e i = e m for i ∈ [1, . . ., n] and e j ≤ e m for j = i.We consider the quadratic function V = 1 2 e Ce and our aim is to show that V | e∈∂N is negative outside the set N .Since e i ≤ e m for all i ∈ [1, . . ., n], it suffices to show Vi | e i ∈∂N i is nonpositive in both cases of e i = e m i and e i = −e m i , i.e., we require to be nonpositive, where N i = {j ∈ M : L ij = 0, i = j} denotes the set of neighbors of the ith node, i.e., the nodes of the network where there exists a direct line connecting the two nodes.We bifurcate the analysis in two cases: (a) e i = e m and (b) e i = −e m .In case of (a), we need the term inside the parenthesis to be negative.Since e m ≥ e j , the diagonal Laplacian term dominates the off-diagonal elements.Comparing the lower bound of the feedback K ii , it can be seen that when the gain obtains its lower value, it is not necessarily true that −K ii e i + Pi z i ≤ 0, unless the deviation between the maximum load demand and the nominal is "large enough."In order to ensure that the desired result holds in every case, we need to Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
impose an additional condition on the feedback Then, the desired inequality Vi | e i =e m i ≤ 0 holds as the expression inside the parenthesis is a summation of strictly negative elements.In case of (b), we note that, by construction, at the point e = −1 n e m , the velocity is given by +K ii e m + Pi when K ii , z i obtain their respective lower bound and the load demand P i is maximized.Then, for any e j ≥ −e m , P i ≤ P m i and z i ≥ z i , we have that and, therefore, we obtain the desired inequality Therefore, the convex nature of set N allows us to conclude that the velocity vector assigned to every point e ∈ ∂N by the map γ : N → T N is subtangential to the set N and by direct application of Nagumo's theorem, the set N is a positive invariant under the flow of the error dynamics with φ e : N → N .The invariance of set N combined with the asymptotic stability of every equilibrium point in N implies the following result.
Corollary 1: The set N is within the region of attraction of system (8).
Remark 2: The above results reveal the immediate effect of the feedback gain K on the degree of conservativeness of N .Higher values of K result in higher magnitudes of the velocity introduced by the linear part and, therefore, a decrease in the maximum displacement of the equilibrium point from the origin.Therefore, there exists a subset of N parameterized by K, Ñ (K) ⊂ N , such that Ñ is a positive invariant and contains the limit of the flow for any e 0 ∈ N , i.e., lim t→∞ φ e (e 0 , t) ∈ Ñ (K).However, while a high feedback gain would diminish the effect of the load demand on the voltage dynamics, this has a negative effect on the converter current steady state, as can be seen by the reference current parameterization î = −Ke + u.Higher values of the feedback introduce larger deviations of the reference current î from the desired reference u, revealing an inherent tradeoff in the choice of the system parameters.
The results are illustrated in Fig. 3, where the set S(P ) is a polytope similar to N but parameterized by an equilibrium point corresponding to some load demand with P i < P m i .It is seen that the vector field drives the flow on the boundary of S(P ) while N is a positive invariant.The above results also suggest that the set Ω = N ⊕ ẑ is a positive invariant for the uncertain system.A question that arises is that under which conditions, the uncertain state flow φ v , parametrized by the control input î, is driven in Ω while v ∈ X is satisfied at all times.One way to guarantee this is to prove the existence of a robust control invariant set (RCI) for the uncertain voltage dynamics within the constraint set.To this end, we recall the following definition from Blanchini [39], adapted to our setting.

Definition 1 (Robust Control Invariant Set):
The set R is said to be robust control invariant for the system (5), if there exists a control î ∈ C such that for all initial states v 0 ∈ R and any w ∈ W , the flow of the system satisfies v = φ(v 0 , î, w, T ) ∈ R, for all T > 0.
Using this notion of RCI sets, we require Following the previous analysis, we recall that Ω = N ⊕ ẑ with N being a positive invariant set; therefore, we can formulate the above question in the nominal setting, by defining the nominal constraint sets as The above indicates that the input constraint set is "shrinking" with larger choices of the control gain K while, on the other hand, the nominal state constraint set Z approaches the size of the original constraint set X.As it was also discussed in Remark 2, this is the expected tradeoff that needs to be considered when choosing the tuning parameter K.The aim of this interplay is to proportionally distribute the high-transients effect of the CPL between the current and the voltage trajectories, in order to satisfy the operational constraints.Furthermore, we note that the computation of the positive invariant set can be done offline and only utilizes the information regarding the bound on the load demand; this significantly simplifies the development of a control scheme for the nominal voltage dynamics and substantially reduces conservativeness stemming from the need of instantaneous load measurements during operation.

B. Formulation of the Nominal Voltage Controller
The problem now becomes of choosing references within the nominal constraint set, i.e., ẑ ∈ Z, and proving the existence of a control invariant set for the nominal dynamics.While the former can be trivially satisfied, the latter requires the association of a control policy with the candidate control invariant set.A common way to resolve this is by employing an MPC control scheme.In MPC, the control policy is generated by solving a finite horizon optimal control problem subject to the system constraints and dynamics.The cost function is often adopted as a quadratic function penalizing the deviation of the system current state from a reference point, thus achieving regulation of the flow to a desired equilibrium.At each time instant, a sequence of control actions is generated and the first element is used as an input to the system while the rest of the sequence is discarded.This process is repeated in the next time instant, thus achieving the receding horizon implementation.Efficient techniques have been developed to solve the continuous-time counterpart of the optimal control problem that usually involves an approximation of the solutions using a numerical solver.Some of the most commonly adopted techniques are the interior point method or sequential quadratic programming (SQP), see [40,Ch. 10] for a detailed analysis.
One way to achieve a stabilizing and recursively feasible controller is by adopting a positive definite cost, compact constraint sets, and invariant terminal ingredients [17], [41].Noting Assumption 1, it is straightforward to satisfy the first two requirements.However, the CPL destabilizes every fixed point ẑ ∈ Z.To overcome this, we define where η o is the optimal control policy generated by solving the MPC problem.Thus, the nominal dynamics take the form We provide the following result to establish the asymptotic stability of the terminal dynamics.Proposition 3 (Stability of the terminal dynamics): Considering the nominal system (13), there exists a δ > 0 such that the fixed point ẑ is an asymptotically stable equilibrium with a region of attraction A = {z ∈ M : z − ẑ ≤ δ}, if and only if To investigate the stability of the fixed point ẑ, we linearize (13) and find the Jacobian matrix The point ẑ is asymptotically stable if and only if the Jacobian is a Hurwitz matrix.Exploiting the properties of the Laplacian matrix we can define a worst case scenario at λ max (−L) = 0.Then, in combination with the fact that the rest of the terms are symmetric matrices, we can deduce the scalar condition which results in the required inequality i for all ẑi ∈ Z i .Shifting the axes to the desired point ẑ, the new dynamics can be described by ż = J z + g(z) where z = z − ẑ and for some ε > 0 and δ > 0, g(z) satisfies g(z) ≤ ε z in some neighborhood of the origin, i.e., z ≤ δ.Then, let V (z) = z P z be a Lyapunov candidate function, where P is a solution of PJ + J P = −Q with Q 0. Note that existence and uniqueness of P is guaranteed due to J being a Hurwitz matrix [31,Th. 4.6].The time derivative of V (•) results in V = z P(J z + g(z)) + (z J + g (z))P z (14) where we have used the fact that g(z) ≤ ε z and z Qz > λ min (Q) z 2 .Now, for ε < λ min (Q) 2 P , the function V (z) = z P z is a Lyapunov function for the shifted system as V ≤ 0, and thus, there exists a δ > 0 such that the region of attraction of Corollary 2: Any level set Z 0 satisfying Z 0 = βA for 0 ≤ β ≤ 1 is a positive invariant for the dynamics (13).
Following the above result, we formalize the definition of the optimization problem as where (•) is a positive definite function.The resulting control is applied in a receding horizon fashion, where, at some time t 1 , the first element of the resulting optimal control sequence is used as input to the system and the process is repeated at the next sampling instant t 2 = t 1 + T .Note that a direct consequence of Proposition 3 and Corollary 2 is that choosing Z f = Z 0 and J f = 1 2 (z f − ẑ) P(z f − ẑ) leads to recursive feasibility of (17) provided that the problem is feasible at some initial time t 0 .In addition, any sublevel set of the stabilizing terminal cost function is control invariant for the nominal dynamics.For completeness purposes, we will provide here a sketch of the recursive feasibility result and we refer the interested reader to [40, Ch. 5] for a detailed analysis.

Proposition 4 (Recursive feasibility of OCP):
Let the problem ( 17) be feasible at some initial state z 0 and time t 0 > 0.Then, the problem remains feasible for all t > t 0 .
Proof: Let z(z 0 , t 0 ) and η(z 0 , t 0 ) denote the feasible solutions of the optimal control problem, respectively, with initial state z 0 and time t 0 .In addition, let z(φ z (z 0 , t 0 + T ), t 0 + T ) and η(φ z (z 0 , t 0 + T ), t 0 + T ) be the candidate state and control sequences at the next immediate sampling instant, i.e., at time t 0 + T .Assuming no additive uncertainties to the system, following Proposition 3 and Corollary 2, the candidate solutions at time t 0 + T are formulated as the tails of the ones at time t 0 where the predicted state at the end of the horizon satisfies z f (t 0 + T + t f ) ∈ Z f .Positive invariance of Z f implies that z(φ z (z 0 , t 0 + T ), t 0 + T ) and η(φ z (z 0 , t 0 + T ), t 0 + T ) are feasible and the problem is recursively feasible for all t > t 0 .

V. STABILITY ANALYSIS OF THE CASCADED DYNAMICS
Following the analysis of the previous sections, the problem has been shifted to regulating the following cascaded dynamics: The stability of a cascaded structure has been thoroughly investigated in the literature, see [42].The conventional procedure is to separate the dynamics into a driving subsystem and a driven one, where the state of the former is considered an input to the latter.Then, the asymptotic stability of the overall dynamics follows from the asymptotic stability of the driving subsystem and the asymptotic stability or some form of boundedness of the driven dynamics.Given a desired reference nominal voltage z ss , we define the equilibrium set for (19) with state vector x = (i, σ, z, e) as We begin the stability analysis by showing the asymptotic stability of the node current.

Theorem 2 (Lyapunov stability of the driving dynamics):
For every node i ∈ M, the C 1 function W is a control Lyapunov function for the subsystem (19a), (19b), and the subsystem is asymptotically stable with equilibrium point described by (20)  Proof: Following the definition of W, it can be seen that it is positive definite.Therefore, all we need to show is the negative definiteness of the first derivative.Following the assumption on the time-scale separation of the dynamics, the input current is constant, i.e., îi = ī.Then, the time derivative of W i results in Following La Salle's invariance theorem and noting that the derivative vanishes only at i ∈ E, the driving subsystem is asymptotically stable with respect to the equilibrium set E. In addition, following Proposition.1, the region of attraction for the current equilibria is the set C. The next theorem combines the aforementioned results to establish the stability of the cascaded dynamics under the proposed control scheme.
Proof: Asymptotic stability of the driving subsystem follows from Theorem 2. In addition, following Proposition 4, the dynamics described by (19c) are also asymptotically stable, where the proof follows the common approach of exploiting the recursive feasibility properties of the optimal control problem and using the cost function as a Lyapunov function for the system.Finally, boundedness of (19d) follows from Proposition 2. Thus, the cascaded dynamics are an interconnection of Lyapunov stable driving subsystems (19a), (19b), (19c), and a bounded-driven subsystem (19d).Then, according to Isidori et al. [42], the equilibrium points in (20) are asymptotically stable for the cascaded dynamics (19).

VI. SIMULATIONS
In this section, we demonstrate the proposed control scheme in a simulated scenario of a seven-node meshed network, see Fig. 4. We require the nominal voltage to reach given references while satisfying the "tighter" constraint sets; hence, the uncertain voltage to always remain within the original constraint set X.The node voltage evolution is depicted in Fig. 5, where each node voltage is always contained within the respective set N i .The current trajectories of Node 1, 3, and 7, along with the respective generated references, are shown in Fig. 6.It is seen that both the trajectories and the references are contained within the desired constraint set, validating the analysis developed in the previous sections.In addition, we note that the fluctuations of the Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.) respectively.The constrained region is represented with the black solid lines ( ) and the voltage references by black dashed lines ( ).The voltage trajectories are within the respective S i ( ) at all times.
current references are a result of the varying load demand, which also affects the error dynamics.Here, the choice to allow larger current fluctuations is made; however, as it was demonstrated by the previous sections, one may choose to limit the fluctuations of the input currents by enlarging the set N , hence allowing larger deviations between the nominal and uncertain node voltages.
The rated voltage is set to v * = 100 V.The current loop control parameters are chosen as k p,i = 600, M i = 9000 resulting in a maximum current I max,i = 17 A. The nominal power demand is Pi = 500 W and we bound the maximum deviation at |δP i | ≤ 500, see Fig. 7.We furthermore parameterize the voltage constraint set as X i = {v i ∈ R n : 97.9 ≤ v i ≤ 102.6}.A quadratic cost functional is chosen of the form (z − ẑ, η − η) = z − ẑ 2 + η − η 2 and the problem is solved using the SQP solver and "ICLOCS" MATLAB toolbox provided in [43].Finally, we choose the voltage control parameters as K i = 50 and K z,i = 4.

VII. CONCLUSION
In this study, we propose a robust control scheme that ensures stability of the closed-loop system in the presence of fluctuating load demand in a meshed dc MG.As it was shown, a "tube" behavior naturally arises from the interplay of the control policy with the rest of the dynamics, allowing us to constrain the uncertain system trajectory within a predefined positive invariant set centered on a generated nominal trajectory.Furthermore, an inherent tradeoff was revealed between the choice of the tube size and the availability of the current values that reside within the operational constraints and satisfy the load demand.We employed an MPC scheme to generate the nominal voltage trajectory and include the constraint sets in the control design process.To ensure the recursive feasibility of the receding horizon optimization problem, the nominal control action was parameterized by a feedback gain that ensures positive invariance of the terminal ingredients.Finally, a demonstration of the results was given by a simulation scenario on a seven-node network with local CPLs.
In future versions, the potential of extending the current approach to plug-and-play scenarios will be investigated.Additionally, we will focus on enlarging the feasibility region of the optimization problem by allowing the constraint sets to be time-varying.

Manuscript received 9
March 2023; accepted 24 May 2023.Date of publication 16 June 2023; date of current version 1 March 2024.This work was supported by the Research Committee of the University of Patras via "C.CARATHEODORY" program under Grant 81359.Recommended by Associate Editor Joshua A. Taylor.(Corresponding author: Grigoris Michos.)

Fig. 2 .
Fig. 2. Node circuit modeled as a dc/dc buck converter connected in parallel to a local CPL.

Fig. 3 .
Fig. 3. Vector field and positive invariant set N for a two-node system.(a) Given a constant power demand P m , any solution starting in N will remain in N and converge to ∂N.(b) For a power demand P ∈ P , the region of the dynamic map shrinks within the set S(P ) ⊂ N .