Distributed Optimization for Reactive Power Sharing and Stability of Inverter-Based Resources Under Voltage Limits

Reactive power sharing and containment of voltages within limits for inverter-based resources (IBRs) are two important, yet coupled objectives in ac networks. In this article, we propose a distributed control technique to simultaneously achieve these objectives. Our controller consists of two components: a purely local nonlinear integral controller which adjusts the IBR voltage setpoint, and a distributed primal-dual optimizer that coordinates reactive power sharing between the IBRs. The controller prioritizes the voltage containment objective over reactive power sharing at all points in time; excluding the IBRs with saturated voltages, it provides reactive power sharing among all the IBRs. Considering the voltage saturation and the coupling between voltage and angle dynamics, a formal closed-loop stability analysis based on singular perturbation theory is provided, yielding practical tuning guidance for the overall control system. To validate the effectiveness of the proposed controller for different case studies, we apply it to a low-voltage microgrid and a microgrid adapted from the CIGRE medium-voltage network benchmark, both simulated in the MATLAB/Simulink environment.

This dependency causes an inherent trade-off between precise reactive power sharing and individual bus voltage regulation, motivating a significant volume of research work on this topic.Different centralized, decentralized, and distributed voltage and reactive power control techniques have been proposed for IBRs.The distributed techniques have attracted significant attention in power system control, especially for large-scale integration of IBRs [3].Compared to their centralized counterparts, they rely on the exchange of information only between neighboring IBRs.In addition, they show better performance and accuracy than decentralized solutions, such as the droop control technique [5].Therefore, it seems that the real-time distributed techniques can be a viable strategy in many situations [3], [6].

A. Literature Review and Research Gaps
Distributed voltage and reactive power sharing control of IBRs has been studied in many papers.Some works have focused solely on the voltage regulation task and have not considered the reactive power sharing problem.The only objective in these papers is to regulate the voltages of the IBRs to a setpoint.This setpoint may be constant, or may, be updated by an external controller.For example, in [7]- [12] and some references therein, assuming that only a few IBRs can directly access the voltage setpoint, a leader-follower consensus algorithm is used for the IBRs to follow this setpoint which is considered a virtual leader.
Conversely, in other lines of research, reactive power sharing is considered as the main objective, and the voltage control requirements are either neglected or discussed only briefly.For example, in [13]- [19], distributed consensus algorithms are used to ensure a proportional reactive power sharing among the IBRs, regardless of the impacts of the controllers on the voltages.However, voltage regulation and reactive power sharing are both important, yet coupled and conflicting; therefore, they should be considered simultaneously.
Simultaneous reactive power sharing and voltage regulation has also been studied.In [20], a leader-follower consensusbased control is proposed for reactive power sharing and voltage tracking problems, where the voltage setpoint is given by a critical bus voltage regulator.Different versions of this scheme are studied and proposed in [21]- [24].A somewhat similar controller is proposed in [25], where unlike in [20], [22]- [24] it is assumed that all the IBRs can directly access the voltage setpoint.These controllers, however, use a single integrator for achieving both the objectives; therefore, the accuracy of reactive power sharing and voltage regulation highly depends on the choice of control gains.The existence of the trade-off between the two objectives is discussed in [24], [25] as well.In [25], tuning of the control gains is suggested as a possible solution for dealing with this trade-off, while in [24], the issue is left as an open problem.
Another combination of control objectives is precise reactive power sharing and average voltage regulation [26]- [35].In this approach, instead of the individual voltages of the IBRs, their estimated average is regulated at a setpoint.To this end, in [26]- [28], the leader-follower consensus algorithm is used for average voltage regulation.Based on the leader-less consensus algorithm, a controller is proposed in [29] where the voltage setpoint of each IBR is corrected by two terms providing average voltage regulation and accurate reactive power sharing, separately.Similarly, two other approaches are proposed in [30], [31], but power sharing is achieved by adjusting the droop coefficient [30] or by changing the virtual impedance [31].To improve the voltage profile and accuracy under input disturbances, some modified controllers are also proposed in [32]- [35].
While the above-mentioned control schemes can provide average voltage regulation, they may result in large deviations in the individual voltages of the IBRs, violating the limits provided by grid standards, e.g., IEEE 1547 [36].Therefore, in many applications, constraining the individual voltages within limits (voltage containment) seems to be a more practical objective [34], [35], [37], [38].In [34], [37], the problem is formulated as an optimization problem, and some controllers based on the primal-dual gradient method are developed.However, these methods require knowledge of the grid model and exchange a relatively large amount of information among the IBRs.In [38], along with a consensus-based control for reactive power sharing, a leader-follower voltage containment controller is proposed to force the voltages into a safe band imposed by some minimum and maximum "leader" IBRs.However, the accuracy of reactive power sharing and voltage containment in this method relies on the selection of the right leaders; i.e., one must already know which IBRs take voltages closer to the minimum and maximum limits and select them as the leader IBRs.In another attempt to bound the voltages, [35] introduces a voltage variance estimation and control loop to the scheme of [29].However, in this method, one "special" IBR is left out of the reactive power sharing task so that the other units can reach simultaneous accurate reactive power sharing and bounded voltages.
Summarizing, we have observed the following research gaps.The works in [7]- [19] have studied either regulation of the individual IBR voltages or reactive power sharing, but not both, while none of the papers in [7]- [33] have considered the operational IBR voltage limits.The accuracy of reactive power sharing and voltage regulation/containment under the proposed controllers in [20]- [25], [34], [35], [37], [38], depends strongly on the choice of control parameters such that if not properly designed, even when steady-state reactive power sharing under voltage limits is possible, these controllers may not provide it.Finally, a rigorous study of stability and synchronization of the power network considering the coupling between angle and voltage dynamics is absent in the above works.

B. Contributions
To address the observed research gap, we propose a distributed control scheme for IBRs to simultaneously achieve voltage containment and reactive power sharing.The main contributions made in this paper are as follows.C1) In our proposed method, we make use of a distributed primal-dual optimizer to generate a globally-unique setpoint to be tracked by a purely local nonlinear integral controller that regulates the IBR's reactive power and tunes its voltage setpoint.This architecture allows maintaining the user-defined voltage constraints, not only in steady state but at all points in time while ensuring that the reactive power demand is shared among the IBRs with a high accuracy.If the above-described reactive power sharing is not possible due to saturation of the voltages, then our controller excludes only the IBRs with saturated voltages from the reactive power sharing task and allows the other IBRs, which are operating away from the voltage limits, to reach a high-accuracy reactive power sharing; i.e., the controller prioritizes voltage containment over reactive power sharing but does not punish all the IBRs.C2) We analyze the system's steady state using graph theory and state its properties.Considering the coupling between voltage and angle dynamics and the voltage saturation, we rigorously study the stability of the system using the Lyapunov method (as recommended in [2]).To this aim, we consider a timescale separation between the dynamics of the primal-dual optimizer and the voltage-angle dynamics, conduct a singular perturbation analysis, and find the stability conditions.We also provide some practical insights into the selection of the control parameters based on the IEEE 1547 standard [36] and the stability analysis.C3) To validate our findings, we adapt the proposed scheme to two test systems, simulated in the MATLAB/Simulink environment.One of the systems is based on a subnetwork of the CIGRE benchmark medium-voltage distribution network.
Our first attempt to address the observed research gap was presented in [39], where we introduced a preliminary version of our control architecture.In this paper, we extend the work in [39] in the following ways.First, we include a leakage term in the local integrator channel to provide anti-wind-up action.Second, we reformulate the selection of the integrator setpoint as an optimization problem.Third, we add a formal stability proof and a parameter selection guideline.Finally, we add a new simulation case study based on a low-voltage microgrid.
The rest of the paper is structured as follows.Section II contains the system modeling and problem statement.In Section III, we introduce our proposed control scheme.We conduct steady state and stability analyses of the closed-loop system in Section IV, where we also provide parameter selection guidelines.In Section V, we present and discuss the simulation results for different case studies.Finally, Section VI concludes the paper.

II. SYSTEM MODELING, POWER SHARING DEFINITION, AND DROOP CONTROL BEHAVIOR
A. Inverter-Based Electric Power Network Under the hierarchical control policy [2], [3], the innermost control loops of inverters are tasked with controlling the LC filter's inductor current and capacitor voltage by generating proper switching signals (see Fig. 1).While different inner loop designs have been proposed, all are designed to act very fast, such that the subsystem denoted by red dashed lines in Fig. 1 has a high bandwidth [3]; virtual impedance control can also be embedded in this subsystem to provide additional decoupling between active and reactive powers and improve system performance [40].For example, in our simulation case High-Bandwidth Subsystem Fig. 1.An inverter-based resource (IBR), governed by the primary control in (1).The high-bandwidth subsystem is denoted by dashed lines.The detailed low-level control structure used in this paper can be found in, e.g., [40].
studies in Section V, we use the cascaded control structure described in [40] and references therein.We also assume that the IBRs use the well-known droop control [40] or equivalently virtual synchronous machine (VSM) control technique [41] as their primary controller, which operates slowly compared with the internal control loops.
In a multi-vendor power system, however, the detailed structure and dynamics of the fast internal controllers are not easily accessible.Therefore, for high-level control design and stability studies, it is preferable to use a simplified generic model for each primary-controlled IBR [42]- [44].For the ith IBR we will use the model where θ i and ω i are the phase angle and angular frequency of the IBR, V i and V set i are the IBR voltage and its setpoint, and ω nom and V nom are the nominal frequency and voltage.The state variables Ω i and v i are the frequency and voltage deviations induced by droop (VSM) controllers in (1b) and (2b), respectively.The constants τ Ω and τ v are the frequency and voltage time constants, respectively.The constants m ω i and m V i are the IBR's frequency and voltage droop coefficients, respectively.The apparent power S rated i is the rated capacity of the IBR while P i and Q i are respectively its active and reactive power injections, which are related to θ = (θ 1 , . . ., θ n ) and V = (V 1 , . . ., V n ) through the following power flow equations where θ ij = θ i − θ j is the phase difference between IBRs i and j; G ij and B ij are the elements of the network's reduced conductance and susceptance matrices [45,Ch. 6.4].

B. Power Sharing Definition and Review of Droop Control
In this subsection, we define power sharing among IBRs and review the steady-state behavior of droop control.As notation, for any variable x, let x denote its steady-state value.
Definition 1.The microgrid system (1)-(3) achieves reactive power sharing if Qi /S rated i = Qj /S rated j = α Q for some α Q .We define active power sharing similarly, using P instead of Q.
According to the droop control in (1b) and (2b) we have Since steady-state frequency is global, for every i and j we have Ωi = Ωj .Thus, following the conventional droop control design criteria [40], by selecting equal frequency droop coefficients for the IBRs, i.e., m ω i = m ω j , we have Pi /S rated i = Pj /S rated j for every i and j, i.e., the frequency droop controller (1b) enforces proportional active power sharing.However, since vi = vj for every i and j does not necessarily hold, selecting m V i = m V j does not guarantee Qi /S rated i = Qj /S rated j ; i.e., the voltage droop controller (2b) cannot enforce reactive power sharing in the same way.In what follows, we propose a distributed control scheme to provide reactive power sharing considering the IBRs voltage limits.

III. PROPOSED CONTROLLER
In this section, we introduce our proposed controller.The controller consists of two subsystems that will be introduced separately: a) a nonlinear leaky integral controller for regulating reactive power ratios of the IBRs and maintaining the voltage limits, and b) a distributed optimizer for obtaining the optimal setpoint for this integrator.

A. Integral Reactive Power Regulation Under Voltage Limits
Let V min i and V max i denote minimum and maximum the desired operational voltage limits for IBR i, with average value ) and maximum allowable deviation ) from that average.In place of the conventional voltage controller (2), we propose the nonlinear integral controller where v i is the state variable of the integrator (4b) with time constant τ v , and where β > 0 is sufficiently small.The variable λ i is a setpoint for the utilization ratio Q i /S rated i , obtained by the optimizer, which will be subsequently described in (8), in the next subsection.The non-negative function ρ i (v i ) is a nonlinear leakage coefficient, defined as The main ideas behind the controller (4) are as follows.
• Since tanh is bounded between −1 and 1, (4a) ensures that at all points in time.In other words, voltage containment is achieved by construction. 1  • The first term in (4b) provides integral action for the utilization ratio Q i /S rated i to track the provided setpoint λ i .The (small) term β∆ i tanh(v i /∆ i ) provides damping, which will assist in our subsequent stability analysis.

B. Distributed Optimization of the Integrator Setpoint λ i
In (4b), λ i acts as a setpoint for the utilization ratio Q i /S rated i .By Definition 1, reactive power sharing will be achieved if the equilibrium values λi are equal, i.e., if λi = λj for all IBRs i and j.We now discuss the optimal selection λi for this setpoint and introduce a distributed algorithm for its online computation.
Following the above discussion, the optimal setpoint selection λi can be formulated via the following optimization problem (5b) We will be seeking a distributed online solution to this optimization problem.To this end, we assume that the IBRs can exchange information over a communication network modeled with an undirected (bidirectional) and connected communication graph; see Appendix A for more info on graph theory.With a ij denoting the elements of the adjacency matrix and N i the set of neighbours of IBR i, the problem ( 5) is equivalent to where k > 0. The constraint (6b) implies that λi = λj for all i and j.We define the Lagrangian associated with the problem (6) as L( λ1 , ζ1 , . . ., λn , ζn ) = C( λ1 , . . ., λn ) + n i=1 ζi zi , where C( λ1 , . . ., λn ) is the total cost function used in (6a) and ζi is the Lagrange multiplier associated with the constraint zi = 0.The problem ( 6) is a quadratic minimization program with linear constraints; hence, Slater's condition holds, and KKT conditions provide necessary and sufficient conditions for optimality [46].In other words, λi , ζi , and zi are optimal if and only if they satisfy the KKT conditions [46, Ch. 5.5] The solution ( λi , ζi ) of ( 7) can be computed in a distributed manner via the so-called primal-dual dynamics [47] where λ i and ζ i are now dynamic state variables which are exchanged between neighboring IBRs in real-time.The parameters τ p and τ d are the primal and dual dynamics time constants, which for our purposes are tunable gains.
To summarize the overall control architecture: the subsystem (8) generates the setpoint λ i to be tracked by the regulator (4b), while the regulator (4b) generates the voltage setpoint (4a) that is saturated within limits; (4b) also provides an anti-wind-up function through the leakage term ρ i (v i )v i when necessary.The general scheme of the proposed controller is shown in Fig. 2.

IV. STEADY-STATE, STABILITY ANALYSIS, AND CONTROLLER GAIN SELECTION
The closed-loop system consists of the angle dynamics (1), the voltage controller ( 4) and ( 8), and the power grid model (3).
In this section, we analyze the steady state of the closed-loop system, study its stability, and state its properties.We also give some insights on the selection of the control parameters.

A. System Steady State and its Properties
We begin by writing the system dynamics in a compact form.Let x = col(x 1 , . . ., x n ) denote the column vector composed of elements x 1 , . . ., x n .We define m = diag(m ω 1 , . . ., m ω n ), S = diag(S rated 1 , . . ., S rated n ), and ρ(v) = diag(ρ 1 (v 1 ), . . ., ρ n (v n )) as well.With this, we can write the differential equations of (1), (4), and (8) in the compact form where L is the Laplacian matrix of the communication graph defined in Appendix A, We can also compactly write the power flow equations (3) and the voltage (4a) as Out first result describes equilibrium points of ( 9)- (10).
Lemma 1 (Steady State).Consider the system (9)-( 10) and suppose that the ac power network has a synchronization frequency of ω syn .Then any steady state of the system satisfies where where Proof.If the ac network has a synchronization frequency of ω syn , then we have θ = ω syn 1 n .Setting this equality together with ẋ = 0 n for any other variable x in (9), we can simply derive the steady-state equations (11).Next, we prove (12).
According to (11a), we have Ω = (ω syn − ω nom )1 n .Multiplying this equation by 1 n , we get ω syn − ω nom = 1 n 1 n Ω and hence Ω = 1 n 1 n 1 n Ω.On the other hand, from (11b) we have Ω = −m S −1 P , where we used m ω i = m for all i.Using the last two equations, we can derive (12a).By connectivity of the communication graph, every solution of equation (11e) has the form λ = α Q 1 n for some α Q ∈ R [48,Ch. 6].Since the graph is undirected, we also have 1 n L = 0 n [48,Ch. 6]; multiplying (11d) by 1 n and using this property we get 1 n λ = 1 n S −1 Q. Setting λ = α Q 1 n in this equation, we can finally arrive at (12b).
Based on Lemma 1, we can now state several practically important properties of the steady state enforced by our controller.
Proposition 1 (Steady State Properties).Consider a steady state as given in (11), let N denote the set of all the IBRs, and define the set of voltage-saturated IBRs as If m ω i = m for all i ∈ N , then the steady state described by Lemma 1 has the following properties: 1) Active Power Sharing and Frequency Regulation: Active power sharing is achieved among all the IBRs and the microgrid's synchronization frequency is ω syn = ω nom − m 1 n 1 n S −1 P .2) Voltage Containment: The steady-state voltages are all in the safe range, i.e., Vi i.e., the IBRs achieve reactive power sharing with a small error proportional to β. 4) Partial Reactive Power Sharing: i.e., only the IBRs that are not in N sat achieve the described almost accurate reactive power sharing, and for the IBRs that belong to N sat , the sharing accuracy decreases (deteriorates) as ρ i (v i ) increases.
Proof.According to (12a) and (11b), we have S −1 P = α P 1 n and hence Pi /S rated i = Pj /S rated j = α P , for every i and j, which according to Definition 1 underlines that active power sharing is achieved among all the IBRs.Inserting (12a) into (11a), we can also write ω syn = ω nom − m 1 n 1 n S −1 P , which proves property 1.The second property is obvious, as we have −1 < tanh(•) < 1. Next, we prove properties 3 and 4.
Inserting (12b) into (11c), and considering ∆tanh ∆ −1 v = V − V , we have Now if N sat = ∅, then for all i we have ρ i (v i ) = 0. From (13b), we can therefore write which proves property 3. Using the definition of the set N set , we can similarly write which, according to Definition 1, proves property 4.

B. Stability Analysis
We want to analyze the stability for the system (9).To this end, we first take some steps to simplify the system dynamics and then analyze stability for the simplified version of (9).
As our controller will maintain voltages within limits around their nominal values, it seems reasonable to use a linearized power flow model to describe the network behavior around the operating point.
Assumption 1. Around a nominal operating point, the power flow equations in (10) can be approximated by where J P θ and J P V (resp.J Q θ and J Q V ) are the n × n Jacobian matrices of f P (θ, V ) (resp.f Q (θ, V )) with respect to θ and V at the linearization point, respectively; w P and w Q are the corresponding intercepts of the linear functions.The matrices J P θ and J Q θ each have an eigenvalue at 0 with corresponding right eigenvector 1 n .
We next reduce the order of the system and transform it into relative coordinates, which allows us to leverage singular perturbation analysis [49,Ch. 11] and find the stability conditions.
1) Model Reduction and Coordinate Transformation: As they are tunable control parameters, we can make the following assumption about the time constants τ Ω , τ p , τ d , τ v in (9).Assumption 2. We have τ Ω , τ p < < τ v and τ p < < τ d .
According to the low-pass filters (9b) and (9d), we have Under Assumption 2, the terms τ Ω Ω and τ p λ can be viewed as some negligible parasitic effects; therefore, the system dynamics are mainly governed by (9a), (9c), and (9e).Indeed, one may apply singular perturbation theory to rigorously reduce the order of the system dynamics to the dynamics of θ, v, and ζ (for an example, see [50]); instead, we omit the details and simply eliminate the left-hand sides of equations (9b) and (9d) and consider Ω = −mS −1 P and (I n + kL)λ = −Lζ + S −1 Q.Therefore, considering the linearized power flow equations in Assumption 1, the system (9) reduces to where We now want to study the stability of the steady state for the system (15) via singular perturbation analysis.In particular, the analysis in [49,Theorem 11.3] requires a system evolving on Euclidean space and an exponentially stable fixed point for the fast dynamics.In order to satisfy these requirements, we instead analyze a version of the system (15) By Assumption 1 and connectivity of the communication graph, the matrices J P θ , J Q θ , and L satisfy [48,Ch. 6].Using T 1 n in (16a) and these properties, we compute that for and write the system dynamics (15) in the new coordinates as According to (16), the first columns of T P θ , T Q θ , and T ζ are all zeros, which means the first elements of x θ and x ζθ av and ζ av -do not influence the dynamics in ( 17) at all.Therefore, (17) is the interconnection of the two cascaded subsystems, given by where their components are given in Appendix B. It should be noted that to obtain ( 18)-( 19) from the dynamics ( 17), we have used the properties Clearly, the dynamics of r θ , r ζ , and v do not depend on θ av and ζ av .Therefore, the steady states of ( 17) and hence (15) are stable, if and only if the steady state of ( 18) is stable.In what follows, we discover this.
2) Timescale Separation and Singular Perturbation Analysis: We are now interested in studying the stability of the steady state of the system (18) using the idea of timescale separation by considering (18a)-(18b) as the slow dynamics and (18c) as the fast dynamics.The following theorem states the stability conditions under these considerations.
Theorem 1 (Exponential Stability for ( 18)).Suppose that the linear matrix inequality in the variables P θ and D v has a solution, where P θ is symmetric, D v is diagonal, and Q is where Then, there exists ε > 0 such that for all τ d < ε τ v the steady state of the system (18) is exponentially stable.
Proof.We consider ε = τ d /τ v small and (18c) as the fast dynamics; therefore, the velocity ṙζ ∝ (1/ε) can be large when ε is small and r ζ in (18c) may rapidly converge to a root of In other words, the subsystem (18c) may quickly achieve a quasi-steady state, where r . We now define the error between the actual r ζ and this quasi-steady state as We can therefore write (18) as the singular perturbation problem below [49,Ch. 11].
where R new vθ and R new vV are given in (20c) and We want to examine the stability of the steady state of ( 21) by examining the reduced system (22a)-(22b) and boundary-layer system (22c), given as where t = t/ε is a stretched timescale with t the time.From Appendix B, we have By the connectivity of the communication graph and the definitions of K, I r , and T , one can establish that the matrix R ζ is negativedefinite; hence, there exists a symmetric matrix P y 0 such that P y R ζ +R ζ P y ≺ 0. With the matrix P y and the matrices P θ 0 and D v 0 in (20a), we now take the following Lyapunov candidates for the slow (22a)-(22b) and fast (22c) dynamics.
where rθ = r θ − rθ , ṽ = v − v, and ỹ = y − ȳ with rθ , v, and ȳ the steady states in (22), and h(τ ) the following function which is element-wise strictly increasing in τ .Time derivatives of S s and S f are Inserting the dynamics ( 22) into (25), we have where δ(ṽ) = ρ(v)(v) − ρ(v)v, η = col(r θ , h(ṽ)) and Q is as given in (20b).The functions δ(ṽ) and h(ṽ) are both elementwise increasing with respect to ṽ; therefore, we have If the matrix inequality (20a) holds, we can write where α s > 0 is the smallest eigenvalue of −(Q + Q ).We can also write where α f > 0 is the smallest eigenvalue of −(P y ζ + R ζ P y ).Using (27a)-(27c), we can now bound the derivatives in (26) as which, together with the fact that S s and S f are both positivedefinite and radially unbounded, show exponential stability of the steady states of the reduced and boundary-layer dynamics.Now, we can say that the singularly perturbed system (21) satisfies all the assumptions of [49, Theorem 11.3]; therefore, there exists ε > 0 such that for all ε < ε or equivalently τ d < ε τ v , the steady state of ( 21) is exponentially stable.

C. Intuition on Parameter Selection
The tunable control parameters in (1), ( 2), (4), and ( 8) are m ω i , m V i , τ Ω , τ p , τ v , τ d , k, and β.Following the standard droop control design [40], we select the droop coefficients m V i = ∆ i and m ω i = m = 2π∆f max for all i, where ∆f max is the maximum steady-state frequency deviation.In what follows, we introduce the impacts and limitations of the remaining parameters and propose a selection procedure.1) (τ p , τ Ω ): According to Assumption 2, our stability analysis is based on τ Ω , τ p < < τ v and τ p < < τ d ; therefore, the smaller τ Ω and τ p , the more reliable our stability analysis.We suggest starting the selection procedure by selecting a small time constant for the low-pass filter (8a), e.g., τ p = 0.01s.One can, however, select a larger τ p for better filtering, if required.On the other hand, according to (1a), decreasing the frequency time constant Here, it should be noted that, in practice, a small τ d requires fast (low-latency) inter-IBR data transmissions.But, selecting a large τ d leads to a large τ v , which in turn causes slower regulation of the IBR voltage and its reactive power.Therefore, while selecting the control parameters, we should consider the practical standards, e.g., IEEE 1547 [36], on this matter.For example, according to [36,Ch. 5.3] the response time for voltage-reactive power control, depending on the mode and application, varies between 1 to 10 seconds.Selecting τ Ω , τ d ∈ [0.1s, 1s] and following the above suggestions, we get τ v ∈ [1s, 10s] which lies in this acceptable range.3) (k, β): Clearly, β helps solvability of the linear matrix inequality (20a); it increases the eigenvalues of −(Q + Q ) and hence the convergence rate α s in (28).But, according to Proposition 1, it degrades the steady-state reactive power sharing.Therefore, we suggest selecting a desired k > 0 first2 and then selecting a small β such that: i) the linear matrix inequality (20a) has a solution, and ii) for every IBR i, the value β∆ i /V i is an acceptable upper bound of the error | Qi /S rated i − α Q |.

V. CASE STUDIES AND SIMULATION RESULTS
To verify the effectiveness of the proposed controller, we applied it to a low-voltage 5-bus meshed microgrid system, simulated in MATLAB/Simscape Electrical software environment.I.The nominal voltage and frequency of the grid are 220-V (RMS) and 50-Hz, respectively.As shown in Fig. 3, the microgrid consists of five local loads energized by five IBRs.Each IBR feeds its corresponding main bus/load via an output connector.Table I shows the electrical and control specifications of the system.It is to note that, in our simulations, we adapted the detailed model of the inverters and internal control loops from [40].

A. Case Study 1: Activation and Load Change
We assume that the droop controllers in (1b) and (2b) control the system before activating the proposed controller.According to Fig. 4(a)-(b), we can see that the voltages deviate from the nominal value, and the IBRs do not share the reactive power proportionally.After activation of the controller at t = 10s, the IBRs start changing their voltages so that their reactive power ratios become equal and, at the same time, their voltages maintain within limits (0.95, 1.05) [p.u.].These results are in line with Properties 2 and 3 in Proposition 1.
At t = 25s, the load at bus number 5 decreases by 80%, which means to keep the proportional sharing, the 5th IBR must feed the other loads instead of the lost 80% local load.Therefore, in a collaborative effort to reach an agreement on a new equal power ratio, this IBR increases its voltage, and the other ones reduce their voltages until the IBRs 1, 2, and 4 reach an agreement on Q i /S rated i .However, the 3rd and 5th IBRs fail to join this agreement because their voltages are already saturated at the minimum and maximum limits, respectively.However, they have come close to the point of agreement and stayed there.We can observe their effort in reaching a consensus with the other IBRs in Fig. 4(d) and Fig. 4(f).After t = 25s, the 5th IBR keeps integrating and increasing v 5 /∆ 5 to increase its voltage to the maximum.However, as the voltage is saturated using the tanh function, the voltage does not change much.Therefore, at t ≈ 26s, when v 5 > 3∆ 5 , the leakage coefficient ρ 5 takes a positive value to prevent the integrator wind-up.Meanwhile, the 3rd IBR also keeps integrating but decreasing h 3 /∆ 3 , until its voltage gets saturated at t ≈ 34s and ρ 3 also takes a positive value.These results are in line with properties 2 and 4 in Proposition 1.After restoring the lost load at t = 40s, the IBRs re-achieve proportional reactive power sharing, and their leakage coefficients are all restored to zero.Fig. 4(c) and Fig. 4(e) show the primal-dual variables; we can see that thanks to the dual variables, the primal variable always converges to the average of the reactive power ratios (see (12b)), no matter if the voltages are saturated or not.Therefore, they are treated as a globallycommon variable like frequency (see Fig. 6(h)) and used as a reliable reference for the reactive power ratios of the IBRs.We can also see the active and frequency responses in Fig 4(g)-(h), reflecting the impacts of the voltage-reactive power controller on the frequency-active power dynamics.

B. Case Study 2: CIGRE Benchmark and Voltage Level Shift
We also applied our controller to a system based on the Subnetwork 1 of the European medium-voltage (20-kV,50-Hz) distribution network benchmark, provided by CIGRE Task Force C6.04.02 [52].Except for the following modifications, all the specifications of the test system are the same as the original network.Based on the Task Force recommendation, the simulated microgrid is an isolated 9-bus subnetwork of the CIGRE system composed of buses number 3 to 11.All the distributed generators at each bus are lumped into one single dispatchable IBR governed by the proposed controller.To meet the maximum load demand in the islanded microgrid the rated power of the IBRs are all increased by 50%.The control gains are similar to the previous case study, but the voltage limits are (0.98, 1.02) [p.u.].The simulation results are shown in Fig. 5.At t = 10s, the controller is activated and the voltages and reactive powers are controlled properly.At t = 20s, we shift the voltage level by setting the new limits (1.01, 1.05) [p.u.].At t = 30s and t = 40s, we disconnect and connect back the residential loads at buses 6 and 8.The results highlight that under the proposed method, we can shift the voltage level in a controlled way while keeping reactive power sharing at different voltage levels.

VI. CONCLUSION
Voltage regulation and reactive power sharing in power systems are two highly coupled control objectives.This coupling is because reactive power flow between two nodes depends more strongly on their voltage differences than the absolute values of the voltages.We proposed a nonlinear controller based on a hyperbolic tangent function and a distributed primal-dual optimizer.The controller provides the IBRs with acceptable reactive power sharing while keeping their voltages within some user-defined limits.We also found stability conditions for the system considering the voltage-angle couplings, under timescale separation between the voltage and optimizer dynamics.The numerical simulations, followed by a proposed parameter selection guideline, indicated a promising performance from the proposed method in controlling the voltage level of the network and achieving reactive power sharing among the IBRs.

APPENDIX A COMMUNICATION NETWORK MODEL AND GRAPH THEORY
An inter-IBR data network can be modeled by an undirected graph where the IBRs and communication links are considered its nodes and edges, respectively.Let G = (N , E, A) be a graph with N = {1, ..., n}, E ⊆ N × N , and A = [a ij ] ∈ R n×n being its node set, edge set, and adjacency matrix, respectively.If the nodes i and j directly exchange data, they are neighbors, meaning that (i, j) ∈ E and (j, i) ∈ E, and a ij = a ji = 0; otherwise, a ij = a ji = 0. Let N i = {j | (j, i) ∈ E} and d i = j∈Ni a ij be the neighbor set and in-degree associated with node i, respectively.Laplacian matrix of G is defined as L = D − A, where D = diag{d i }.A walk (or path) from node i to node j is an ordered sequence of nodes such that any pair of consecutive nodes in the sequence is an edge of the graph.A graph is connected if there exists a walk between any two nodes [48].
in which the system is transformed into relative coordinates.Let us now define the change of coordinates x θ = T θ = [θ av r θ ] and x ζ = T ζ = [ζ av r ζ ] , where θ av and ζ av are respectively the average values of the elements in θ and ζ, the vectors r θ and r ζ belong to R n−1 , and the transformation matrix T ∈ R n×n is

Fig. 4 .
Fig. 4. Simulation results for Case Study 1; (a) reactive power ratios Q i /S rated i , (b) voltages V i , (c) primal variables λ i , (d) normalized integrator states v i /∆ i , (e) dual variables ζ i , (f) leakage coefficients ρ i (v i ), (g) active power ratios P i /S rated i , and frequencies f i = 1 2 ω i /π.

Fig. 5 .
Fig. 5. Simulation results for Case Study 2; (a) voltages and (b) reactive power ratios Q i /S rated i .
[36,n as the Rate of Change of Frequency (RoCoF), which should be limited in practice (see, for example,[36, Table  21]).Thus, we suggest selecting τ Ω = m ω i /(2πRoCoF ), where RoCoF is the maximum withstandable initial RoCoF after a step change in P i from 0 to S rated i or vice versa.2) (τ v , τ d ): Following the discussion in the previous step, to make our stability analysis more reliable we suggest selecting τ d > > τ p and τ v > > τ Ω , τ p , for example, τ d ≥ 10τ p and τ v ≥ max{10τ Ω , 10τ p }. On the other hand, from Theorem 1, one should select τ d < ε τ v .Since the exact value of ε is not easily available, we select τ v as large as possible and τ d as small as possible, for example, we suggest selecting τ v ≥ 10τ d .Combining these suggestions, we get τ d ≥ 10τ p and τ v ≥ max{10τ Ω , 10τ d }.

TABLE I CONTROL
AND ELECTRIC SPECIFICATIONS FOR THE LV SYSTEM IN FIG.3