1 Introduction↩︎

Distributed energy resources (DERs) operate in different functional modes depending on their interface and control objectives. Grid-following (GFL) DERs synchronize with the grid and inject prescribed power, whereas grid-forming (GFM) DERs establish and regulate voltage and frequency, particularly in islanded or weak-grid conditions [1]. Primary droop control enables decentralized operation of DERs [2]; however, by itself, it generally does not restore frequency to its nominal value, does not guarantee accurate reactive power sharing under heterogeneous impedances, and does not explicitly enforce device-level operational constraints [3]–[5]. These limitations motivate secondary control mechanisms [5]–[7] that compute voltage adjustments for GFM-DERs and active-power references for GFL-DERs through constrained optimization problems.

This paper studies the closed-loop dynamics induced by a sampled-data optimization-based secondary controller for DER networks. The controller is designed to coordinate voltage regulation, reactive power sharing, and frequency support while respecting local operational constraints. Its broader architecture, distributed implementation, design rationale, plug-and-play operation, privacy-preserving information exchange, and performance validation are presented in a separate companion article. In contrast, the present paper is devoted to the mathematical analysis of the resulting closed-loop dynamics induced by the distributed secondary controller. The controller model and equations needed for the analysis are specified explicitly in Section 2; consequently, the results developed here are self-contained and do not rely on external implementation details. The closed-loop system is mathematically challenging because the physical inverter dynamics evolve continuously, while the secondary controller is updated only at sampling instants. At each sampling time, constrained optimization problems are instantiated using measured voltage, reactive-power, and frequency signals. The computed optimizer outputs are then applied through interpolation over the next sampling interval. Hence, the resulting dynamics are not those of a standard continuous-time droop-controlled system, but of a nonlinear sampled-data system with optimization in the loop. Existing analyses of secondary control in power systems primarily focus on small-signal, local behavior and continuous-time models [6], [7]. As a result, they do not directly provide large-signal guarantees for the closed-loop voltage and power dynamics and are not applicable to the setting considered here.

We first formulate the closed-loop sampled-data model induced by the optimization-based secondary controller, capturing nonlinear power-flow coupling, filtered power measurements, constrained optimization updates, and interpolation-based actuation. We then provide a self-contained stability analysis of the resulting closed-loop model. The main results establish that i) the DER voltage and filtered reactive-power dynamics remain within an explicit forward-invariant set, ii) a positive steady state exists in this set, iii) the steady state achieves voltage regulation and equal per-unitized reactive power sharing, and iv) the frequency dynamics are input-to-state practically stable with respect to secondary control signals and DER voltages. Together, these results provide rigorous closed-loop guarantees for sampled-data optimization-based secondary control beyond local small-signal analysis.

2 Control Setup↩︎

We consider a network of DERs operating in either GFM or GFL mode. Both GFM- and GFL-DERs operate within a hierarchical control architecture [8] comprising three layers: zero-level, primary-level, and secondary-level control. The GFM-DERs provide voltage and frequency regulation through primary-level droop control, while the GFL-DERs inject active power through primary-level \(\mathrm{PQ}\) dispatch [8].

Although primary droop enables decentralized operation, it generally does not achieve the desired steady-state objectives [5]. In particular, \(P\)-\(\omega\) droop induces a steady-state frequency deviation from the nominal value, while unequal line impedances prevent accurate reactive power sharing under \(Q\)-\(V\) droop and cause GFM-DER terminal voltages to deviate from nominal values. Thus, a secondary-level controller is needed to restore frequency, regulate voltages, and enforce reactive power sharing while respecting device-level constraints [5].

The secondary controller considered in this paper modifies the GFM-DER voltage references and the GFL-DER active-power references using sampled measurements and constrained distributed optimization. Next, we discuss the model of this secondary controller. The subsequent sections analyze the resulting nonlinear sampled-data closed-loop dynamics. References for \(P\)-\(\omega\)/\(Q\)-\(V\) droop control law [9] for GFM-DER \(i\) are given by \[\begin{align} \omega_i &= \textstyle \overline{\omega} - r_{\omega_i} \overline{P}_i, \;\;\overline{P}_i = \frac{1}{(\tau_{P_i} s+1)}P_i, \tag{1}\\ V_i &= \textstyle \overline{V}- r_{V_i}\overline{Q}_i + U_i^\star,\;\;\overline{Q}_i = \frac{1}{(\tau_{Q_i} s+1)}Q_i. \tag{2} \end{align}\] Here, \(\overline{\omega}\), \(\omega_i, \overline{V}, V_i\) are the nominal and measured frequency and voltages, respectively, and \(r_{\omega_i}, r_{V_i} \in \mathbb{R}_{>0}\) are the droop coefficients. The instantaneous active and reactive power \(P_i\) and \(Q_i\), respectively, are given by \[\begin{align}P_i &= \textstyle G_{ii} V_i^2 - \sum_{k\in N_i} B_{ik} V_i V_{k} (\delta_i - \delta_{k}) - P_i^\mathrm{gfl} \tag{3},\\Q_i & = \textstyle -B_{ii} V_i^2 + \sum_{k\in N_i} B_{ik} V_i V_{k},\tag{4} \end{align}\] where, the complex admittance, \(Y_{ik} = Y_{ki} \in \mathbb{C}\), between buses \(i \in \{1,2,\dots,n\}\) and \(k\in \{1,2,\dots,n\}\) are represented as \(Y_{ik} := -j B_{i k} = -j B_{ki} =: Y_{ki}\) with \(B_{i k} = B_{ki} < 0\). The set of neighbors of a bus \(i \in \{1,2,\dots,n\}\) is defined as \(N_i:= \{k| k\in \{1,2,\dots,n\}, k\neq i, Y_{ik} \neq 0 \}\) and \(Y_{ii} = G_{ii} + j (B^{\mathrm{sh}}_{i} + \sum_{k\in N_i} B_{ik}) := G_{i i} + j B_{ii}\) where, \(G_{i i} > 0\) is the shunt conductance, \(B^{\mathrm{sh}}_{i} < 0\) is the shunt susceptance and \(B_{i i} < 0\). The quantities \(\overline{P}_i\) and \(\overline{Q}_i\) in 1 2 are the averaged active and reactive power injections of GFM-DER \(i\), obtained by passing the instantaneous active and reactive power \(P_i\) and \(Q_i\) through low-pass filters with time constants, \(\tau_{P_i}, \tau_{Q_i} \in \mathbb{R}_{>0}\) respectively.

Signals, \(U_i^\star(t), i \in \{1,2,\dots, n\}\) in 2 are determined using secondary-level controller via the following design law, \[\begin{align} \label{eq:U95i95design95law}U_i^\star(t) := u_i(t) - \beta_{Q_i} r_{V_i} \overline{Q}_i(t) + \beta_{V_i} (\overline{V}- V_i(t)), \end{align}\tag{5}\] where \(\beta_{V_i}, \beta_{Q_i} \in \mathbb{R}\) are design hyper-parameters. With the compensation signals \(U_i^\star\) in 5 the reference in 2 becomes \[\begin{align} \label{eq:voltage95reference} V_i &= \textstyle \overline{V}- r_{V_i}\frac{(1 + \beta_{Q_i})}{(1 + \beta_{V_i})} \overline{Q}_i + \frac{1}{(1 + \beta_{V_i})}u_i, \end{align}\tag{6}\] where, \(u_i(t), t\in(t_s,t_{s+1}]\), is updated from the last value \(u_i(t_s)\) via the first-order hold like interpolation as follows² \[u_i(t) =\Bigg\{\begin{align} & u_i(t_s), t \in (t_s, t_s+ \Delta_s),\\ & \textstyle u_i(t_s) + \frac{(t - t_s- \Delta_s)(x_{s,i} - u_{i}(t_s))}{\Delta}, t \in (t_s+ \Delta_s, t_{s+1}] \end{align}\] for all \(i\in\{1,\dots,n\}\). Here \(x_{s,i}\) denotes the solution of the following problem instantiated with measurements \([V_i(t_s), \;\overline{Q}_i(t_s)] \in \mathbb{R}^2, i \in \{1,\dots, n\}\): \[\label{eq:gfm95optimization} \begin{align} &x_s:= \textstyle \mathop{\mathrm{argmin}}_{x_{s,1}, \dots, x_{s,n}} \sum_{i=1}^n \;\frac{1}{2}(x_{s,i} - \alpha_i (t_s))^2 \\ subject to \;x_{s,i} &= x_{s,j}, \;\forall i,j\in \{1,\dots,n\}, \\ x_{s,i} &\in D_i^{V_\Delta}(t_s),\;\forall i \in \{1,\dots,n\} \\ \alpha_i (t_s) &:= (1+\beta_{Q_i})(\overline{V}- V_i(t_s)) \\ &- \beta_{V_i} r_{V_i} \overline{Q}_i (t_s) \;\forall i \in \{1,\dots,n\}. \end{align}\tag{7}\] Here, \(D_i^{V_\Delta}(t_s) := \big\{ \zeta \big | | \zeta - (1+\beta_{Q_i}) r_{V_i} \overline{Q}_i (t_s) + \beta_{V_i} (\overline{V}- V_i(t_s)) | \leq V_\Delta \big\}\) with a prescribed parameter \(V_\Delta \in \mathbb{R}_{>0}\).

In addition, the secondary-level controller designs active power references \(P_i^\star\), \(i \in \{1,\dots,n\}\), for the primary-level \(\mathrm{PQ}\)-dispatch controllers of the GFL-DERs as \[\begin{align} \label{eq:Pi95star95design95law} P_i^\star(t) := P_i^{\min} + p_i(t)(P_i^{\max} - P_i^{\min}). \end{align}\tag{8}\] Here, \(P_i^{\min}\) and \(P_i^{\max}\) denote the minimum and maximum active power generation capacities of GFL-DER \(i\), respectively, and \(p_i(t),t\in(t_s,t_{s+1}]\), are updated from \(p_i(t_s)\) as \[p_i(t) =\Bigg\{\begin{align} & p_i(t_s), t \in (t_s, t_s+ \Delta_s),\\ & \textstyle p_i(t_s) + \frac{(t - t_s- \Delta_s)(y_{s,i} - p_i(t_s))}{\Delta}, t \in (t_s+ \Delta_s, t_{s+1}] \end{align}\] for all \(i\in\{1,\dots,n\}\), where \(y_{s,i}\) denotes the solution of \[\label{eq:gfl95optimization} \begin{align} y_s:= & \textstyle \mathop{\mathrm{argmin}}_{y_{s,1},\dots, y_{s,n}} \sum_{i=1}^{n} \;\frac{1}{2}(P_{y_{i}}^\star)^2 + b_i P_{y_{i}}^\star + c_i \\ &subject to \; 0 \leq y_{s,i} \leq 1 \;\forall i \in \{1,\dots, n\}, \\ &\textstyle \sum_{i=1}^{n} P_{y_{i}}^\star = \rho_\mathrm{d}(t_s) \\ & P_{y_{i}}^\star := P_i^{\min} + y_{s,i} (P_i^{\max} - P_i^{\min}) \;\forall i \in \{1,\dots, n\}, \end{align}\tag{9}\] where \(y_s=[y_{s,1},\dots,y_{s,n}]\in\mathbb{R}^{n}\) denotes the solution of the optimization problem, and \(P_i^{\min}\) and \(P_i^{\max}\) are the active power limits of GFL-DER \(i\), as in 8 and \(\rho_\mathrm{d}(t_s)\) is the total active power requirement. The GFL-DER power set-point estimate in 9 is parametrized as \(P_{y_i}^\star:=P_i^{\min}+y_{s,i}(P_i^{\max}-P_i^{\min})\), with \(0\leq y_{s,i}\leq 1\), ensuring that the resulting active power set-points satisfy capacity constraints. The solution to 9 determines the GFL-DER power reference \(P_i^\star\) in 8 . The GFL-DER output \(P_i^\mathrm{gfl}\) tracks \(P_i^\star\) with negligible, assumed zero, error [10] and therefore, \[\begin{align} \label{eq:gfl95power} \textstyle P_i^\mathrm{gfl} (t) = P_i^\star(t) = P_i^{\min} + p_i(t)(P_i^{\max} - P_i^{\min}) \;\forall t. \end{align}\tag{10}\] Collectively, problems 7 and 9 can be formulated as \[\label{eq:opt95prob95equivalent} \begin{align} & \textstyle \mathop{\mathrm{minimize}}_{\zeta_1, \zeta_2, \dots, \zeta_{n_\zeta}} \quad \sum_{i=1}^{n_\zeta} f_i(\zeta_i)\\ & subject to\zeta_i = \zeta_j\;\forall i, j\in \{1,2,\dots,n_\zeta\}, \\ &\zeta_i \in \mathcal{X}_i, \;\overline{H}_i \zeta_i = \overline{h}_i, \;\underline{H}_i \zeta_i \leq \underline{h}_i \; \forall i \in \{1,\dots,n_\zeta\}. \end{align}\tag{11}\] Here, each DER maintains a local decision variable \(\zeta_i\). The objective \(f_i:\mathbb{R}^{n_\zeta}\rightarrow\mathbb{R}\), constraint set \(\mathcal{X}_i\subseteq\mathbb{R}^{n_\zeta}\), equality constraints \(\overline{H}_i\zeta_i=\overline{h}_i\), and inequality constraints \(\underline{H}_i\zeta_i\leq\underline{h}_i\) are local to DER \(i\) and can be enforced in a decentralized manner. The consensus constraints \(\zeta_i=\zeta_j\) couple the local decisions and require network-level coordination. The developed secondary-level controller, therefore, solves 11 using the distributed discrete-time \(\boldsymbol{\texttt{DC-DistADMM}}\) algorithm developed in [11]. The \(\boldsymbol{\texttt{DC-DistADMM}}\) algorithm has a geometric rate of convergence; after \(\theta\) iterations, \[\begin{align} \label{eq:admm95iter95comp} \|\zeta_i^{(\theta)} - \zeta^\star \|^2 \leq \Upsilon (0.75)^{\theta}, \quad for all \;i, \end{align}\tag{12}\] where \(\Upsilon\) is a known constant determined from problem data. Hence, an accuracy of \(\epsilon\) requires only \(\theta_\epsilon=O(\log(1/\epsilon))\) iterations, and therefore the sampling interval \(\Delta_s\) in the sampled-data secondary controller can be determined according to the chosen accuracy \(\epsilon\).

3 Stability Analysis↩︎

Using 1 2 , GFM-DER \(i\) closed-loop dynamics are, \[\begin{align} \dot{\delta}_i &= \omega_i, \quad \tau_{P_i}\dot{\omega_i} = -(\omega_i - \overline{\omega}) - r_{\omega_i}P_i, \tag{13}\\ \tau_{Q_i}\dot{V_i} &= -(V_i - \overline{V}) - r_{V_i}Q_i + U_i^\star + \tau_{Q_i}\dot{U}_i^\star. \tag{14} \end{align}\]

3.1 Analysis of GFM-DER Voltage Dynamics Loop↩︎

Substituting 5 in 14 we get, \[\begin{align} \label{eq:V95dot95expanded} \tau_{Q_i}\dot{V_i} &= -(V_i - \overline{V}) - r_{V_i}Q_i + U_i^\star + \tau_{Q_i}\dot{U}_i^\star \nonumber \\ & = -(V_i - \overline{V}) - r_{V_i}Q_i + u_i - \beta_{Q_i} r_{V_i} \overline{Q}_i \\ &+ \beta_{V_i} (\overline{V}- V_i) + \tau_{Q_i}\dot{u}_i - \tau_{Q_i}\beta_{Q_i} r_{V_i} \dot{Q}^\mathrm{avg}_i - \tau_{Q_i}\beta_{V_i} \dot{V}_i \nonumber. \end{align}\tag{15}\] Let \(\tilde{\beta}_{Q_i} := (1 + \beta_{Q_i}), \;\tilde{\beta}_{V_i} := (1 + \beta_{V_i}), \;\tilde{\beta}_i := \tilde{\beta}_{Q_i}/\tilde{\beta}_{V_i} \;\forall i.\) Substituting \(\tau_{Q_i} \dot{\overline{Q}}_i = Q_i - \overline{Q}_i\) in 15 we get for all \(i\), \[\begin{align} \label{eq:V95dot95expanded95intermediate} \tau_{Q_i}(1 + \beta_{V_i})\dot{V_i} &= -(1+\beta_{V_i})(V_i - \overline{V}) - r_{V_i}Q_i \nonumber \\ &- \beta_{Q_i} r_{V_i} \overline{Q}_i - \beta_{Q_i} r_{V_i} (Q_i - \overline{Q}_i) + u_i + \tau_{Q_i}\dot{u}_i \nonumber \\ &= -\tilde{\beta}_{V_i}(V_i - \overline{V}) - \tilde{\beta}_{Q_i}r_{V_i}Q_i + u_i + \tau_{Q_i}\dot{u}_i, \end{align}\tag{16}\] Therefore, using 16 we have for all \(i \in \{1,2,\dots,n\}\), \[\begin{align} \label{eq:V95dot95expanded95final} \dot{V_i} = \textstyle -\frac{1}{\tau_{Q_i}}(V_i - \overline{V}) - \frac{\tilde{\beta}_i r_{V_i}}{\tau_{Q_i}}Q_i + \frac{1}{\tau_{Q_i} \tilde{\beta}_{V_i}}u_i + \frac{1}{\tilde{\beta}_{V_i}}\dot{u}_i. \end{align}\tag{17}\] Let \(V:=[V_1,\dots,V_n]^\top\in\mathbb{R}^n, \overline{Q}(V):=\overline{Q} =[\overline{Q}_1,\dots,\) \(\overline{Q}_n]^\top\in\mathbb{R}^{n}\). We establish an invariance result for the closed-loop voltage dynamics around the nominal operating point \(\mathbf{\overline{V}}:=\overline{V}\mathbf{1}_{n}\) and \(\widehat{Q}:=Q(\mathbf{\overline{V}})\). The result provides an explicit forward-invariant set whose size is determined by the network parameters and controller design, characterizing closed-loop boundedness and performance.

Define \(u(t) := [u_1(t), \dots, u_{n}(t)]^\top \in \mathbb{R}^{n}\), and \(\dot{u}(t) := [\dot{u}_1(t), \dots, \dot{u}_{n}(t)]^\top \in \mathbb{R}^{n}\), \(\overline{u}_i :=\tilde{\beta}_{V_i}\tilde{\beta}_i r_{V_i} \widehat{Q}_i\), for all \(i\), \(\overline{u}:= [\overline{u}_1, \dots, \overline{u}_{n}]^\top \in \mathbb{R}^{n}\), \(d_i := |B_{i i}| + \textstyle \sum_{k\in N_i}|B_{i k}|\) for all \(i\), and \(d_B := \sqrt{\sum_{i=1}^{n} d_i^2}\). Fix \((R_V,R_q)>0\), define \(\Omega(R_V,R_q) := \{(\widetilde{V},\widetilde{q}):\|\widetilde{V}\|_2\le R_V,\;\|\widetilde{q}\|_2\le R_q\}, L_Q:= d_B \bigl(2\|\overline{V}\|_2+R_V\bigr)\), and \[\label{eq:constants} \begin{align} B_1 &= \max_i |1+\beta_{Q_i}|+\sqrt{n}\max_i |\beta_{V_i}|, \\ B_2 &= \max_i |\beta_{V_i}r_{V_i}|+\sqrt{n}\max_i |(1+\beta_{Q_i})r_{V_i}|, \\ B_0 &= \sqrt{n}V_\Delta + B_2\|\widehat{Q}\|_2 + \|\overline{u}\|_2, \\ X_\star &:= \textstyle B_0+B_1R_V+B_2R_q, \; D_\star:=\frac{2X_\star}{\Delta}. \end{align}\tag{18}\] Let \(\tau_{\min}:=\min_i \tau_{Q_i}, \tau_{\max}:=\max_i \tau_{Q_i}, \tilde{\beta}_{V,\min}:=\min_i \tilde{\beta}_{V_i},\boldsymbol{\tau}_{\beta r} := \text{diag}(\frac{\tilde{\beta}_i r_{V_i}}{\tau_{Q_i}}), \boldsymbol{\tau}_{\beta_V} := \text{diag}(1/(\tilde{\beta}_{V_i} \tau_{Q_i})),\) \(\boldsymbol{\beta}_V := \text{diag}(1/\tilde{\beta}_{V_i}), \boldsymbol{\tau} := \text{diag}(1/\tau_{Q_i}),\) and \(K_u \le \frac{1}{3\tilde{\beta}_{V,\min}\tau_{\min}^2}, K_d\le \frac{1}{3\tilde{\beta}_{V,\min}\tau_{\min}}.\) Choose \(\varepsilon_{r,1},\varepsilon_{r,2},\varepsilon_u,\varepsilon_d>0\), and \(\varepsilon_q\in(0,2/\tau_{\max})\). Define \(\varepsilon_r:=\varepsilon_{r,1}+\varepsilon_{r,2}\), and \(d_r(R_V) := \textstyle \frac{1}{18\varepsilon_{r,1}} \|\boldsymbol{\tau}_{\beta r}\boldsymbol{\tau} \widehat{Q}\|_2^2 + \frac{L_Q^2}{18\varepsilon_{r,2}} \|\boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}\mathbf{\overline{V}}\|_2^2, + \frac{1}{3} | \mathbf{\overline{V}}^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}\widehat{Q}|, c_V := \textstyle \frac{1}{3\tau_{\max}^2}-\frac{\varepsilon_r+\varepsilon_u+\varepsilon_d}{2}\). \(\varepsilon_{r,1},\varepsilon_{r,2},\varepsilon_u,\varepsilon_d>0\) be chosen such that \(c_V > 0\). Further, let \(\nu>0, \varepsilon_q > 0\) be such that \(0<\nu< \frac{2\varepsilon_q\tau_{\min}^2\,c_V}{L_Q^2}.\) Let \(a_V := \textstyle c_V-\nu\frac{L_Q^2}{2\varepsilon_q\tau_{\min}^2}, \;a_q := \nu\Bigl(\frac{1}{\tau_{\max}}-\frac{\varepsilon_q}{2}\Bigr), M_\Psi := \textstyle \max\left\{\frac{1}{6\tau_{\min}},\frac{\nu}{2}\right\}, \; m_\Psi:=\min\left\{\frac{1}{6\tau_{\max}},\frac{\nu}{2}\right\}, \alpha := \textstyle \frac{\min\{a_V,a_q\}}{M_\Psi}, \; c_\star := \textstyle d_r(R_V) +\frac{K_u^2}{2\varepsilon_u}X_\star^2 +\frac{K_d^2}{2\varepsilon_d}D_\star^2, \rho_\Omega := \textstyle \min\left\{\frac{R_V^2}{6\tau_{\max}}, \frac{\nu R_q^2}{2} \right\}.\) We divide the proof into several steps.

Step 1: Centered dynamics. Since \(V = \widetilde{V} + \mathbf{\overline{V}}, \overline{Q}=\widetilde{q} + \widehat{Q},\) the voltage dynamics 15 can be written as \[\begin{align}\dot{\widetilde{V}} = -\boldsymbol{\tau} \widetilde{V} -\boldsymbol{\tau}_{\beta r}\bigl(Q(V) - \widehat{Q}\bigr) + \boldsymbol{\tau}_{\beta_V}(u - \overline{u}) +\boldsymbol{\beta}_V \dot{u}. \label{eq:centered95voltage95dyn95proof} \end{align}\tag{19}\] Also, using \(\tau_{Q_i}\dot{\overline{Q}}_i = -\overline{Q}_i + Q_i,\) \[\label{eq:centered95filter95dyn95proof} \dot{\widetilde{q}} = -\boldsymbol{\tau} \widetilde{q} + \boldsymbol{\tau}\bigl(Q(V) - \widehat{Q}\bigr).\tag{20}\]

Step 2: Bound the optimizer \(x_s\). Because of consensus constraints \(x_{s,i} = x_{s,j}\), every feasible solution has the form \(x_s= c_{x_s} \mathbf{1}_{n}\), for some scalar \(c_{x_s} \in \mathbb{R}\). Therefore, \(\|x\|_2 = \sqrt{n} |c_{x_s}|.\) Next, from 7 , \(\alpha_i (t_s) = (1+\beta_{Q_i})\overline{V}- (1+\beta_{Q_i})V_i - \beta_{V_i}r_{V_i}Q_i^{\mathrm{avg}}(t_s).\) Thus, \[\begin{align} \|\alpha \|_2 \le & \max_i |1+\beta_{Q_i}| \|\mathbf{\overline{V}} - V\|_2 + \max_i |\beta_{V_i}r_{V_i}| \|\overline{Q}\|_2. \end{align}\] Let \(L_\alpha^V := \max_i |1+\beta_{Q_i}|, L_\alpha^Q := \textstyle \max_i |\beta_{V_i}r_{V_i}|.\) Then, \(\|\alpha\|_2 \le L_\alpha^V \|\mathbf{\overline{V}} - V\|_2 + L_\alpha^Q \|\overline{Q}\|_2.\) Define, \(\overline{\alpha} := \frac{1}{n}\sum_{i=1}^{n} \alpha_i\), and by Cauchy-Schwarz inequality, \[\begin{align} |\bar \alpha| = \textstyle \frac{1}{n} |\mathbf{1}_{n}^\top \alpha| & \textstyle \le \frac{1}{\sqrt{n}}\|\alpha\|_2 \leq \textstyle \frac{L_\alpha^V}{\sqrt n}\|\mathbf{\overline{V}} - V\|_2 + \frac{L_\alpha^Q}{\sqrt n}\|\overline{Q}\|_2. \end{align}\] Now consider sets \(D_i^{V_\Delta}\) in 7 , with the centers given by \(m_i = (1+\beta_{Q_i})r_{V_i}Q_i^{\mathrm{avg}}(t_s) - \beta_{V_i}(\overline{V}- V_i (t_s))\) Therefore, \(|m_i| \le |\beta_{V_i}||\overline{V}- V_i (t_s)| + |(1+\beta_{Q_i})r_{V_i}||Q_i^{\mathrm{avg}}(t_s)|.\) Taking the maximum over \(i\) yields \(\|m\|_\infty \le m_V\|V(t_s) - \overline{V}\|_2 + m_Q\|\overline{Q}(t_s)\|_2\), where \(m_V := \max_i |\beta_{V_i}|, m_Q := \textstyle \max_i |(1+\beta_{Q_i})r_{V_i}|.\) Therefore, the solution \(x_s= c_{x_s} \mathbf{1}_{n}\) is bounded, \(|c_{x_s}| \le |\bar\alpha| + \|m\|_\infty + V_\Delta.\) Substituting the previously derived bounds gives \(|c_{x_s}| \leq \left(\frac{L_\alpha^V}{\sqrt{n}} + m_V \right) \|V(t_s) - \mathbf{\overline{V}}\|_2 + V_\Delta + \textstyle \left(\frac{L_\alpha^Q}{\sqrt{n}} + m_Q \right) \|\overline{Q}(t_s)\|_2\) Multiplying by \(\sqrt{n}\) and using \(\|x_s\|_2=\sqrt{n} |c_{x_s}|\), we conclude \[\begin{align} \label{eq:bound95on95xs} \|x_s\|_2 \le B_1\|\tilde{V}(t_s)\|_2 + B_2\|\overline{Q}(t_s)\|_2 + \sqrt{n}V_\Delta, \end{align}\tag{21}\] where \(B_1,B_2\) are exactly those given in 18 . Hence, by the triangle inequality, \[\begin{align} \|x_s- \overline{u}\|_2 &\le \|x_s\|_2+\|\overline{u}\|_2 \le \sqrt{n}V_\Delta + B_1\|\widetilde{V}(t_s)\|_2 \nonumber \\ &+ B_2\|\widetilde{q}(t_s)\|_2 + B_2\|\widehat{Q}\|_2 + \|\overline{u}\|_2 \nonumber\\ &= B_0 + B_1\|\widetilde{V}(t_s)\|_2+B_2\|\widetilde{q}(t_s)\|_2. \label{eq:x95minus95unom95bound95proof} \end{align}\tag{22}\] Now assume \((\widetilde{V}(t_s),\widetilde{q}(t_s))\in \mathcal{S}_{\rho_\Omega}\), then \(\|\widetilde{V}(t_s)\|_2\le R_V, \;\|\widetilde{q}(t_s)\|_2\le R_q,\) and so 22 yields \[\label{eq:x95minus95unom95bound95region95proof} \|x_s- \overline{u}\|_2 \le B_0+B_1R_V+B_2R_q = X_\star.\tag{23}\]

Step 3: Explicit bounds on \(u - \overline{u}\) and \(\dot{u}\). For \(t\in[t_s,t_s+ \Delta_s)\), the interpolation law gives \(u(t) = u(t_s)\). And for \(t\in[t_s+ \Delta_s,t_{s+1})\), the interpolation law gives \[\begin{align} u(t) = \textstyle \Bigl(1-\frac{t - t_s- \Delta_s}{\Delta}\Bigr)u(t_s) + \frac{t - t_s- \Delta_s}{\Delta}x_s. \end{align}\] Subtracting \(\overline{u}\) for \(t \in [t_s, t_{s+1})\), we have \(u(t) - \overline{u}\) \(= \begin{cases} u(t_s) - \overline{u}\\ \big(1-\frac{t-t_s-\Delta_s}{\Delta}\big)(u(t_s) - \overline{u}) + \frac{t-t_s-\Delta_s}{\Delta}(x_s - \overline{u}). \end{cases}\)

Since both coefficients are nonnegative and sum to one, \[\begin{align} \|u(t) - \overline{u}\|_2 &\le \textstyle \Bigl(1-\frac{t-t_s-\Delta_s}{\Delta}\Bigr)\|u(t_s)- \overline{u}\|_2 \nonumber \\ & \textstyle + \frac{t-t_s-\Delta_s}{\Delta}\|x_s- \overline{u}\|_2 \nonumber\\ & \le \max \bigl\{\|u(t_s) - \overline{u}\|_2, \|x_s- \overline{u}\|_2 \bigr\}.\label{eq:u95convex95bound95proof} \end{align}\tag{24}\] Assume \(\|u(0) - \overline{u}\|_2\le X_\star\). Then, using 23 , we prove by induction over the sampling intervals \(t \in [t_s, t_{s+1})\) that \[\label{eq:u95minus95unom95global95bound95proof} \|u(t) - \overline{u}\|_2\le X_\star.\tag{25}\] as long as \((\widetilde{V}(t),\widetilde{q}(t))\in \mathcal{S}_{\rho_\Omega}\). Indeed, it is true at \(t=0\) by assumption. If it holds at \(t=t_s\), then 24 and 23 imply it holds for all \(t\in[t_s,t_{s+1})\). Also, \(\dot{u}(t)= \begin{cases} 0 & t\in[t_s+ \Delta_s)\\ \frac{x_s- u(t_s)}{\Delta} & t\in[t_s+ \Delta_s,t_{s+1}) \end{cases}\).

\[\begin{align} Hence, \;\|\dot{u}(t)\|_2 & \textstyle \le \frac{1}{\Delta}\bigl(\|x_s- \overline{u}\|_2+\|u(t_s) -\overline{u}\|_2\bigr) \nonumber\\ & \textstyle \le \frac{1}{\Delta}(X_\star+X_\star) = \frac{2X_\star}{\Delta} = D_\star. \label{eq:udot95bound95proof} \end{align}\tag{26}\] Therefore, on \(\Omega(R_V,R_q)\), \[\label{eq:u95udot95explicit95bounds95proof} \|u(t) - \overline{u}\|_2\le X_\star, \quad \|\dot{u}(t)\|_2 \le D_\star.\tag{27}\]

Step 4: Local Lipschitz bound on \(Q(V) - \widehat{Q}\). From the quadratic structure of the reactive-power map, \(\|Q(V) - Q(W)\|_2 \le d_B\bigl(\|V\|_2+\|W\|_2\bigr)\|V-W\|_2.\) Apply this with \(W = \mathbf{\overline{V}}\). Then, \(\|Q(V) - \widehat{Q}\|_2 \le d_B\bigl(\|V\|_2+\|\mathbf{\overline{V}}\|_2\bigr)\|V - \overline{V}\|_2.\) Since \(V=\widetilde{V}+\mathbf{\overline{V}}, \; \|V\|_2\le \|\widetilde{V}\|_2+\|\mathbf{\overline{V}}\|_2,\) and on \(\Omega(R_V,R_q)\) we have \(\|\widetilde{V}\|_2\le R_V\), it follows that \(\|V\|_2\le R_V+\|\mathbf{\overline{V}}\|_2.\) Hence \[\begin{align} \|Q(V) - \widehat{Q}\|_2 &\le d_B\bigl(R_V+\|\mathbf{\overline{V}}\|_2+\|\mathbf{\overline{V}}\|_2\bigr)\|\widetilde{V}\|_2 \nonumber\\ & = d_B\bigl(2\|\mathbf{\overline{V}}\|_2+R_V\bigr)\|\widetilde{V}\|_2 = L_Q\|\widetilde{V}\|_2. \label{eq:LQ95bound95proof} \end{align}\tag{28}\]

Step 5: Derivative of the voltage Lyapunov function. Define \(\Phi(\widetilde{V}):=\frac{1}{6} \widetilde{V}^\top \boldsymbol{\tau} \widetilde{V}.\) Then \(\dot{\Phi} = \dot{\widetilde{V}}^\top \frac{1}{6} \boldsymbol{\tau} \widetilde{V} + \widetilde{V}^\top \frac{1}{6} \boldsymbol{\tau} \dot{\widetilde{V}} = \frac{1}{3} \widetilde{V}^\top \boldsymbol{\tau} \dot{\widetilde{V}}\), since \(\boldsymbol{\tau}\) is diagonal and symmetric. Using 19 , \[\begin{align} \dot{\Phi} &= \textstyle -\frac{1}{3} \widetilde{V}^\top \boldsymbol{\tau}^2 \widetilde{V} -\frac{1}{3} \widetilde{V}^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}(Q - \widehat{Q}) \nonumber\\ & \textstyle +\frac{1}{3} \widetilde{V}^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta_V}(u - \overline{u}) +\frac{1}{3} \widetilde{V}^\top \boldsymbol{\tau}\boldsymbol{\beta}_V\dot{u}. \label{eq:Phi95dot95centered95intermediate95proof} \end{align}\tag{29}\] To make the reactive term symmetric, write \[\begin{align} &\dot{\Phi} = \textstyle -\frac{1}{3} \widetilde{V}^\top \boldsymbol{\tau}^2 \widetilde{V} +\mathcal{T}_Q +\frac{1}{3} \widetilde{V}^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta_V}(u - \overline{u}) +\frac{1}{3} \widetilde{V}^\top \boldsymbol{\tau}\boldsymbol{\beta}_V \dot{u}, \tag{30} \\where \;\mathcal{T}_Q & \textstyle := -\frac{1}{3} \widetilde{V}^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}(Q - \widehat{Q}) \nonumber \\ &\textstyle = -\frac{1}{6} \widetilde{V}^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}(Q-\widehat{Q}) -\frac{1}{6} (Q - \widehat{Q})^\top \boldsymbol{\tau}_{\beta r}\boldsymbol{\tau} \widetilde{V}. \tag{31} \end{align}\]

We now estimate \(\mathcal{T}_Q\) using the improved centered decomposition. Since \(\widetilde{V} = V - \mathbf{\overline{V}}\), we expand 31 as \[\begin{align} &\textstyle -\frac{1}{6} (V - \mathbf{\overline{V}})^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}(Q - \widehat{Q}) -\frac{1}{6} (Q - \widehat{Q})^\top \boldsymbol{\tau}_{\beta r}\boldsymbol{\tau} (V - \mathbf{\overline{V}}) \nonumber\\ & = \textstyle -\frac{1}{6} V^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}Q -\frac{1}{6} Q^\top \boldsymbol{\tau}_{\beta r}\boldsymbol{\tau} V +\frac{1}{6} V^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r} \widehat{Q} \nonumber\\ &\textstyle+\frac{1}{6} (\widehat{Q})^\top \boldsymbol{\tau}_{\beta r}\boldsymbol{\tau} V +\frac{1}{6} (\mathbf{\overline{V}})^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}Q +\frac{1}{6} Q^\top \boldsymbol{\tau}_{\beta r}\boldsymbol{\tau} \mathbf{\overline{V}}\nonumber\\ & \textstyle -\frac{1}{6} (\mathbf{\overline{V}})^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}\widehat{Q} -\frac{1}{6} (\widehat{Q})^\top \boldsymbol{\tau}_{\beta r}\boldsymbol{\tau} \mathbf{\overline{V}}. \label{eq:TQ95expanded95full95proof} \end{align}\tag{32}\] Rearranging terms, we have \[\begin{align} \mathcal{T}_Q &= \textstyle -\frac{1}{6} V^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}Q -\frac{1}{6} Q^\top \boldsymbol{\tau}_{\beta r}\boldsymbol{\tau} V +\frac{1}{3} (\widehat{Q})^\top \boldsymbol{\tau}_{\beta r}\boldsymbol{\tau} \widetilde{V} \nonumber\\ & \textstyle+\frac{1}{3} (\mathbf{\overline{V}})^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}Q. \label{eq:TQ95rearranged95proof} \end{align}\tag{33}\] Using \(Q(V)=\mathrm{diag}(V)\mathcal{B} V\), and defining \(\Sigma(V):=\mathrm{diag}\left(\tilde{\beta}_i r_{V_i}V_i/2\tau_{Q_i}^2\right)\), the first two terms in 33 become \(-\frac{1}{3} V^\top \bigl(\Sigma(V)\mathcal{B}+\mathcal{B}^\top \Sigma(V)\bigr)V.\) Note that the matrix \(\Sigma(V) \mathcal{B} + \mathcal{B}^\top \Sigma(V)\) is positive semi-definite (see Lemma 1), and therefore \(\textstyle -\frac{1}{3} V^\top \bigl(\Sigma(V)\mathcal{B}+\mathcal{B}^\top \Sigma(V)\bigr)V \le 0.\) Hence, \[\begin{align} \mathcal{T}_Q & \textstyle \le \frac{1}{3} (\widehat{Q})^\top \boldsymbol{\tau}_{\beta r}\boldsymbol{\tau} \widetilde{V} +\frac{1}{3} (\mathbf{\overline{V}})^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}Q. \label{eq:TQ95after95PSD95drop95proof} \end{align}\tag{34}\] We now bound the two remaining terms. For the first term, Young’s inequality gives \[\begin{align} \left| \textstyle \frac{1}{3} (\widehat{Q})^\top \boldsymbol{\tau}_{\beta r}\boldsymbol{\tau} \widetilde{V} \right| & \textstyle \le \frac{\varepsilon_{r,1}}{2}\|\widetilde{V}\|_2^2 + \frac{1}{18\varepsilon_{r,1}} \|\boldsymbol{\tau}_{\beta r}\boldsymbol{\tau} \widehat{Q}\|_2^2. \label{eq:TQ95term195bound95proof} \end{align}\tag{35}\] For the second term, using 28 , \[\begin{align} & \left| \textstyle \frac{1}{3} (\mathbf{\overline{V}})^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}Q \right| \le \textstyle \frac{1}{3} \|\boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}\mathbf{\overline{V}}\|_2 \|Q - \widehat{Q}\|_2 + \left| \textstyle \frac{1}{3} \mathbf{\overline{V}}^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}\widehat{Q} \right| \nonumber \\ & \textstyle \le \frac{1}{3} \|\boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}\mathbf{\overline{V}}\|_2 L_Q\|\widetilde{V}\|_2 + \left| \textstyle \frac{1}{3} \mathbf{\overline{V}}^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}\widehat{Q} \right| \nonumber\\ & \textstyle \le \frac{\varepsilon_{r,2}}{2}\|\widetilde{V}\|_2^2 + \frac{L_Q^2}{18\varepsilon_{r,2}} \|\boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}\mathbf{\overline{V}}\|_2^2 + \left| \textstyle \frac{1}{3} \mathbf{\overline{V}}^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta r}\widehat{Q} \right|. \label{eq:TQ95term295bound95proof} \end{align}\tag{36}\] Combining 34 , 35 , and 36 , we get \[\begin{align} \mathcal{T}_Q & \le \textstyle \frac{\varepsilon_r}{2}\|\widetilde{V}\|_2^2+d_r(R_V), \label{eq:TQ95final95bound95proof} \end{align}\tag{37}\] where \(\varepsilon_r=\varepsilon_{r,1}+\varepsilon_{r,2}\). Next, the quadratic term satisfies \[\begin{align} \textstyle -\frac{1}{3} \widetilde{V}^\top \boldsymbol{\tau}^2 \widetilde{V} & \textstyle \le -\frac{1}{3\tau_{\max}^2}\|\widetilde{V}\|_2^2, \label{eq:tau295dissipation95proof} \end{align}\tag{38}\] because the smallest eigenvalue of \(\boldsymbol{\tau}^2\) is \(1/\tau_{\max}^2\). For the control term involving \(u - \overline{u}\), using 27 and Young’s inequality, \[\begin{align} \textstyle \frac{1}{3} \widetilde{V}^\top \boldsymbol{\tau}\boldsymbol{\tau}_{\beta_V}(u - \overline{u}) &\le K_u\|\widetilde{V}\|_2\|u - \overline{u}\|_2 \le \textstyle \frac{\varepsilon_u}{2}\|\widetilde{V}\|_2^2 \nonumber \\ &\textstyle + \frac{K_u^2}{2\varepsilon_u}\|u - \overline{u}\|_2^2 \le \frac{\varepsilon_u}{2}\|\widetilde{V}\|_2^2 + \frac{K_u^2}{2\varepsilon_u}X_\star^2. \label{eq:u95term95bound95proof} \end{align}\tag{39}\] Similarly, for the \(\dot{u}\) term, \[\begin{align} \textstyle \frac{1}{3} \widetilde{V}^\top \boldsymbol{\tau}\boldsymbol{\beta}_V\dot{u} & \textstyle \le K_d\|\widetilde{V}\|_2\|\dot{u}\|_2 \le \frac{\varepsilon_d}{2}\|\widetilde{V}\|_2^2 + \frac{K_d^2}{2\varepsilon_d}\|\dot{u}\|_2^2 \nonumber\\ &\textstyle \le \frac{\varepsilon_d}{2}\|\widetilde{V}\|_2^2 + \frac{K_d^2}{2\varepsilon_d}D_\star^2. \label{eq:udot95term95bound95proof} \end{align}\tag{40}\] Substituting 37 , 38 , 39 , and 40 into 30 yields \[\begin{align} \dot{\Phi} &\textstyle \le -\left( \frac{1}{3\tau_{\max}^2} -\frac{\varepsilon_r+\varepsilon_u+\varepsilon_d}{2} \right)\|\widetilde{V}\|_2^2 +d_r(R_V) \nonumber\\ & \textstyle + \frac{K_u^2}{2\varepsilon_u}X_\star^2 + \frac{K_d^2}{2\varepsilon_d}D_\star^2 \nonumber\\ &= \textstyle -c_V\|\widetilde{V}\|_2^2 +d_r(R_V) + \frac{K_u^2}{2\varepsilon_u}X_\star^2 + \frac{K_d^2}{2\varepsilon_d}D_\star^2. \label{eq:Phi95dot95final95proof} \end{align}\tag{41}\]

Step 6: Derivative of the filter-energy term. From 20 , we have \(\frac{1}{2}\frac{d}{dt}\|\widetilde{q}\|_2^2 = \widetilde{q}^\top \dot{\widetilde{q}} = -\widetilde{q}^\top \boldsymbol{\tau}\widetilde{q} + \widetilde{q}^\top \boldsymbol{\tau}(Q-\widehat{Q}).\) Since \(\textstyle \boldsymbol{\tau} \succeq \frac{1}{\tau_{\max}}I, \; \|\boldsymbol{\tau}\|\le \frac{1}{\tau_{\min}},\) we obtain \(\frac{1}{2}\frac{d}{dt}\|\widetilde{q}\|_2^2 \le -\frac{1}{\tau_{\max}}\|\widetilde{q}\|_2^2 + \frac{1}{\tau_{\min}}\|\widetilde{q}\|_2\|Q - \widehat{Q}\|_2.\) Using 28 , \(\frac{1}{2}\frac{d}{dt}\|\widetilde{q}\|_2^2 \le -\frac{1}{\tau_{\max}}\|\widetilde{q}\|_2^2 + \frac{L_Q}{\tau_{\min}}\|\widetilde{q}\|_2\|\widetilde{V}\|_2\). Applying Young’s inequality with parameter \(\varepsilon_q\), \(\frac{L_Q}{\tau_{\min}}\|\widetilde{q}\|_2\|\widetilde{V}\|_2 \le \frac{\varepsilon_q}{2}\|\widetilde{q}\|_2^2 + \frac{L_Q^2}{2\varepsilon_q\tau_{\min}^2}\|\widetilde{V}\|_2^2.\) \[\begin{align}Thus, \;\textstyle \frac{1}{2}\frac{d}{dt}\|\widetilde{q}\|_2^2 & \textstyle \le -\Bigl(\frac{1}{\tau_{\max}}-\frac{\varepsilon_q}{2}\Bigr)\|\widetilde{q}\|_2^2 + \frac{L_Q^2}{2\varepsilon_q\tau_{\min}^2}\|\widetilde{V}\|_2^2. \label{eq:q95energy95final95proof} \end{align}\tag{42}\]

Step 7: Derivative of the composite Lyapunov function and invariance. Define \(\Psi(\widetilde{V},\widetilde{q}) = \Phi(\widetilde{V})+\frac{\nu}{2}\|\widetilde{q}\|_2^2.\) Multiplying 42 by \(\nu\) and adding to 41 , we get \(\dot{\Psi} \le -\left( c_V-\nu\frac{L_Q^2}{2\varepsilon_q\tau_{\min}^2} \right)\|\widetilde{V}\|_2^2 -\nu\Bigl(\frac{1}{\tau_{\max}}-\frac{\varepsilon_q}{2}\Bigr)\|\widetilde{q}\|_2^2 + c_\star,\) where \(c_\star= d_r(R_V) +\frac{K_u^2}{2\varepsilon_u}X_\star^2 +\frac{K_d^2}{2\varepsilon_d}D_\star^2.\) Due to the choice of \(\nu\), \(a_V= c_V-\nu\frac{L_Q^2}{2\varepsilon_q\tau_{\min}^2}>0.\) Also, since \(\varepsilon_q<2/\tau_{\max}\), \(a_q= \nu\Bigl(\frac{1}{\tau_{\max}}-\frac{\varepsilon_q}{2}\Bigr)>0.\) Therefore, \[\begin{align} \dot{\Psi} &\le -a_V\|\widetilde{V}\|_2^2-a_q\|\widetilde{q}\|_2^2+c_\star. \label{eq:Psi95dot95with95aVaQ95proof} \end{align}\tag{43}\] Next, by definition of \(M_\Psi\), \(\Psi(\widetilde{V},\widetilde{q}) \le M_\Psi\bigl(\|\widetilde{V}\|_2^2+\|\widetilde{q}\|_2^2\bigr).\) Hence, \(a_V\|\widetilde{V}\|_2^2+a_q\|\widetilde{q}\|_2^2 \ge \frac{\min\{a_V,a_q\}}{M_\Psi}\Psi = \alpha \Psi.\) Substituting this into 43 , we obtain \[\begin{align} \dot{\Psi} &\le -\alpha\Psi+c_\star. \label{eq:Psi95dot95scalar95proof} \end{align}\tag{44}\] Since, \(c_\star < \alpha \rho_\Omega\), whenever \(\Psi=\rho_\Omega\), \(\dot{\Psi}\le -\alpha\rho_\Omega+c_\star<0.\) The estimates leading to 44 are valid whenever \((\widetilde{V}(t),\widetilde{q}(t))\in \Omega(R_V,R_q)\), and for almost all \(t\), since \(u(t)\) is piecewise affine. Define the first exit time \(T^\star := \inf\left\{ t\ge 0: (\widetilde{V}(t),\widetilde{q}(t))\notin \Omega(R_V,R_q) \right\}.\) If no such time exists, then \(T^\star=\infty\). On the interval \([0,T^\star)\), the trajectory belongs to \(\Omega(R_V,R_q)\). Therefore, all bounds derived above apply, and 44 holds for almost all \(t\in[0,T^\star)\): \(\dot{\Psi}(t)\le -\alpha\Psi(t)+c_\star.\) By the comparison lemma, for all \(t\in[0,T^\star), \Psi(t) \le e^{-\alpha t}\Psi(0) + \frac{c_\star}{\alpha}\bigl(1-e^{-\alpha t}\bigr).\) If \(\Psi(0)\le \rho_\Omega\), then using \(c_\star<\alpha\rho_\Omega.\) Then, for every \(t\in[0,T^\star)\), \[\Psi(t) \le e^{-\alpha t}\rho_\Omega + \frac{c_\star}{\alpha}\bigl(1-e^{-\alpha t}\bigr) \le \rho_\Omega .\] Thus, \(\Psi(t)\le \rho_\Omega, \forall t\in[0,T^\star)\). Next, by definition of \(\Psi\), \(\Psi(\widetilde{V},\widetilde{q}) = \textstyle \frac{1}{6} \widetilde{V}^\top \boldsymbol{\tau} \widetilde{V} + \frac{\nu}{2}\|\widetilde{q}\|_2^2.\) Finally, by definition of \(m_\Psi\), \(\textstyle \Psi(\widetilde{V},\widetilde{q}) \ge \frac{1}{6\tau_{\max}}\|\widetilde{V}\|_2^2, \Psi(\widetilde{V},\widetilde{q}) \ge \frac{\nu}{2}\|\widetilde{q}\|_2^2.\) Therefore, for all \(t\in[0,T^\star)\), \(\|\widetilde{V}(t)\|_2^2 \le 6\tau_{\max}\Psi(t) \le 6\tau_{\max}\rho_\Omega \le R_V^2,\) and \(\|\widetilde{q}(t)\|_2^2 \le \frac{2}{\nu}\Psi(t) \le \frac{2}{\nu}\rho_\Omega \le R_q^2.\) Hence, \((\widetilde{V}(t),\widetilde{q}(t))\in \Omega(R_V,R_q), \; \forall t\in[0,T^\star).\) Suppose, for contradiction, that \(T^\star<\infty\). Since the state trajectory is continuous, taking the limit \(t \to T^\star\) gives \(\|\widetilde{V}(T^\star)\|_2\le R_V, \; \|\widetilde{q}(T^\star)\|_2\le R_q.\) Thus \[(\widetilde{V}(T^\star),\widetilde{q}(T^\star)) \in \Omega(R_V,R_q),\] which contradicts the definition of \(T^\star\) as the first exit time from \(\Omega(R_V,R_q)\). Therefore, \(T^\star=\infty\). Consequently, \(\Psi(t)\le \rho_\Omega, \;\forall t\ge 0\), and hence \(\|\widetilde{V}(t)\|_2\le R_V, \; \|\widetilde{q}(t)\|_2\le R_q, \; \forall t\ge 0.\) Therefore, \(\Omega(R_V,R_q)\) is forward invariant for all trajectories starting in the sublevel set \(\mathcal{S}_{\rho_\Omega} := \{(\widetilde{V},\widetilde{q}):\Psi(\widetilde{V},\widetilde{q})\le \rho_\Omega\}.\) Moreover, \(\mathcal{S}_{\rho_\Omega}\subseteq \Omega(R_V,R_q)\), and \(\mathcal{S}_{\rho_\Omega}\) is forward invariant. This completes the proof.

We next establish the existence of a steady state within the invariant region. Since the invariance confines trajectories to a compact set, fixed-point arguments can be used to characterize equilibrium behavior. We show that the closed-loop voltage dynamics admit a steady-state operating point consistent with the network power flow and controller structure, and that this equilibrium lies within the invariant region.

Define \(d_{\max} := \max_i \textstyle \left(|B_{ii}|+\sum_{k\in N_i}|B_{ik}|\right), M_Q := d_{\max}(\overline{V}\) \(+ R_V)^2, M_u := B_1\sqrt{n} R_{V} + B_2\sqrt{n}M_Q + \sqrt{n}V_\Delta\), where, \(B_1,B_2\) are as in 18 \(\tilde{\beta}_{Q,\max}:=\max_i |\tilde{\beta}_{Q_i}|, r_{V,\max}:=\max_i |r_{V_i}|\). Let \[\label{eq:selfmap95condition} \textstyle \frac{M_u+\tilde{\beta}_{Q,\max}r_{V,\max}M_Q}{\tilde{\beta}_{V,\min}} \le R_V.\tag{45}\]

We proceed in several steps.

Step 1: \(\mathcal{K}_V\) is a nonempty compact convex positive set. Since \(R_{V}>0\), the set \(\mathcal{K}_V\) is nonempty. Because it is a product of closed bounded intervals, it is compact; it is convex. Since \(\overline{V}> R_{V} >0\). Hence, every \(V\in\mathcal{K}_V\) is componentwise positive: \(V_i\in [\overline{V}-R_{V},\overline{V}+R_{V}] \subset (0,\infty), \;i=1,\dots,n.\)

Step 2: Map \(T\) is continuous. By definition 4 , \(Q(V)\) is continuous for any \(V\) and the solution \(x_s\) of 7 is continuous in \((V,Q(V))\) (see Lemma 2). Therefore, the composition \(V \mapsto x_s\bigl(V,Q(V)\bigr)\) is continuous. Since each component \(T_i\) is obtained from continuous operations, \(T\) is continuous.

Step 3: Uniform bounds for \(Q(V)\) on \(\mathcal{K}_V\). Let \(V\in \mathcal{K}_V\). Since each component satisfies \(|V_i|\le \overline{V}+R_{V}\), we have \(\|V\|_\infty\le \overline{V}+R_{V}\). For each \(i\), using the reactive-power relation 4 , we obtain \(|Q_i(V)| \le \left(|B_{ii}|+\sum_{k\in N_i}|B_{ik}|\right)\|V\|_\infty^2 \le d_{\max}(\overline{V}+R_{V})^2 = M_Q.\) \[\label{eq:Qi95bound95proof}Hence,|Q_i(V)|\le M_Q, \quad i=1,\dots,n.\tag{46}\] which yields \(\|Q(V)\|_2 \le \sqrt{n}M_Q.\) Using 21 , we have \(\|x_s\|_2 \le B_1\|\tilde{V}\|_2 + B_2\|\overline{Q}\|_2 + \sqrt{n}V_\Delta \le \left(B_1\sqrt{n} R_{V} + B_2\sqrt{n}M_Q + \sqrt{n}V_\Delta \right):= M_u.\)

Step 4: \(T(\mathcal{K}_V)\subseteq \mathcal{K}_V\). Fix \(V\in\mathcal{K}_V\). For each \(i\), using 45 , and 46 , we get \[\begin{align} |T_i(V)-\overline{V}| &= \textstyle \left| \frac{x_{s,i}(V,Q(V))-\tilde{\beta}_{Q_i}r_{V_i}Q_i(V)}{\tilde{\beta}_{V_i}} \right| \le \\ &\textstyle \frac{ |x_{s,i}(V,Q(V))| + |\tilde{\beta}_{Q_i}| |r_{V_i}| |Q_i(V)| }{\tilde{\beta}_{V_i}} \le \frac{M_u+\tilde{\beta}_{Q,\max}r_{V,\max}M_Q}{\tilde{\beta}_{V,\min}} \le R_{V} \end{align}\] where the last inequality is exactly 45 . Therefore \(\overline{V}-R_{V}\le T_i(V)\le \overline{V}+R_{V}, \;i=1,\dots,n,\) which shows that \(T(V)\in \mathcal{K}_V\). Since \(V\in\mathcal{K}_V\) is arbitrary, we conclude \(T(\mathcal{K}_V)\subseteq \mathcal{K}_V.\)

Step 5: Existence of a fixed point. The set \(\mathcal{K}_V\) is nonempty, compact, and convex, and since \(T:\mathcal{K}_V\to\mathcal{K}_V\) is continuous, Brouwer’s fixed-point theorem [12] implies that there exists \(V^\star\in\mathcal{K}_V\) such that \(T(V^\star)=V^\star\).

Step 6: Construction of the steady state. Define \(Q^\star:=Q(V^\star), \; (\overline{Q})^\star:=Q^\star, \; u^\star:=x_s(V^\star,Q^\star)\). We claim that the constant signals \(V(t) \equiv V^\star, \; Q(t)\equiv Q^\star, \;\overline{Q}(t) \equiv Q^\star, \;u(t)\equiv u^\star,\) solve the voltage/filter/controller subsystem. First, since \(u(t) = u^\star = x_s(V^\star, Q^\star) = c^\star \mathbf{1}_{\mathrm{gfm}}\), for some constant \(c^\star\), for all \(t\), we have \(\dot{u}_i(t)=0, i=1,\dots,n.\) Second, since \((\overline{Q})^\star=Q^\star\), the averaging filter \(\tau_{Q_i}\dot{(\overline{Q})}^\star = Q^\star_i - (\overline{Q})^\star\) yields \(\dot{(\overline{Q})}_i^\star(t)=0, i=1,\dots,n.\) It remains to be verified that the voltage equation holds. The fixed-point identity \(T(V^\star)=V^\star\) means that, for each \(i\), \(V_i^\star = \overline{V}+ \frac{u_i^\star - \tilde{\beta}_{Q_i}r_{V_i}Q_i^\star}{\tilde{\beta}_{V_i}}.\) Equivalently, \(u_i^\star = \tilde{\beta}_{V_i}(V_i^\star - \overline{V})+\tilde{\beta}_{Q_i}r_{V_i}Q_i^\star.\) Using the relation \(\tilde{\beta}_i=\tilde{\beta}_{Q_i}/\tilde{\beta}_{V_i}\), we may rewrite this as \(0 = -\frac{1}{\tau_{Q_i}}(V_i^\star-\overline{V}) -\frac{\tilde{\beta}_i r_{V_i}}{\tau_{Q_i}}Q_i^\star +\frac{1}{\tau_{Q_i}\tilde{\beta}_{V_i}}u_i^\star.\) Since \(\dot{u}_i=0\), it implies \(\dot{V}_i = 0, i=1,\dots,n\) for all \(i\). Therefore, the voltage equation 17 achieves a steady-state.

Step 7: Positivity. Since \(V^\star\in \mathcal{K}_V\), \(V_i^\star\ge \overline{V}-R_{V}>0, i=1,\dots,n.\) Thus, the steady state voltage is positive.

Hence, the constructed constant signals form a positive steady state of the voltage/filter/controller subsystem.

Having established steady-state existence, we characterize its relation to the secondary-level objectives. We show that the equilibrium induced by the proposed distributed controller coordinates GFM-DERs to achieve voltage regulation and equal per-unitized reactive power sharing.

By Theorem 2, the closed-loop voltage dynamics admit a steady-state operating point. At this equilibrium, the measurements entering 7 are constant, so \(\alpha_i(t_s)\) and \(D_i^{V_\Delta}(t_s)\) are fixed at all steady-state sampling instants. Hence, problem 7 has the same optimizer, denoted \(x^\mathrm{ss}\), at every such instant. By the sampled-data update law, after one update beyond steady state, \(u_i(t_s)=x_i^\mathrm{ss}\) for all sufficiently large \(t_s\) and all \(i\in{1,\dots,n}\). Consequently, \(x_i^\mathrm{ss}-u_i(t_s)=0\), and the interpolation law gives \(\dot{u}_i(t)=0\) for all sufficiently large \(t\). Since \(U_i^\star\) differs from \(u_i\) only by steady-state constant terms, \(\dot{U}_i^\star=0\) for all \(i\in{1,\dots,n}\) at steady state. Setting \(\dot{V}_i = 0\) in 14 and using \(\dot{U}_i^\star = 0, \overline{Q}_i^{\mathrm{ss}} = Q_i^{\mathrm{ss}}\) gives \[\begin{align} \label{eq:Vi95steady95state}(V_i^{\mathrm{ss}} - \overline{V}) = u_i^{\mathrm{ss}} + \beta_{V_i} (\overline{V}- V_i^{\mathrm{ss}}) - \tilde{\beta}_{Q_i} r_{V_i}Q_i^{\mathrm{ss}}. \end{align}\tag{47}\] Let \(x^\mathrm{ss}\) be \(\boldsymbol{\texttt{DC-DistADMM}}\) solution with accuracy \(\varepsilon/4\). Then, \(|V_i^{\mathrm{ss}} - \overline{V}| = \big|u_i^{\mathrm{ss}} + \beta_{V_i} (\overline{V}- V_i^{\mathrm{ss}})- \tilde{\beta}_{Q_i} r_{V_i}Q_i^{\mathrm{ss}}\big| \leq |x^\mathrm{ss}+ \beta_{V_i} (\overline{V}- V_i^{\mathrm{ss}}) - \tilde{\beta}_{Q_i} r_{V_i}Q_i^{\mathrm{ss}}| + |u_i^{\mathrm{ss}} - x^\mathrm{ss}| \leq V_\Delta + \textstyle \sqrt{\Upsilon (0.75)^{\theta_{\varepsilon}}} \leq V_\Delta + \varepsilon/2\), where, we used 12 . From 47 , \[\begin{align} \label{eq:steady95state} (1 + \beta_{V_i})(V_i^{\mathrm{ss}} - \overline{V}) + (1 + \beta_{Q_i}) r_{V_i}Q_i^{\mathrm{ss}} = u_i^{\mathrm{ss}}, \;\forall i. \end{align}\tag{48}\] Since the steady-state inputs \(u_i^\mathrm{ss}\) are obtained from the \(\boldsymbol{\texttt{DC-DistADMM}}\) solution \(x_i^\mathrm{ss}\) [11], the agreement error satisfies \(|u_i^{\mathrm{ss}} - u_j^{\mathrm{ss}}| \leq \textstyle \sqrt{\Upsilon (0.75)^{\theta_{\varepsilon}}} \leq \varepsilon\) for all \(i,j\). Thus, with suitable design parameters, the secondary-level inputs 5 achieve the objectives of tasks \(\mathcal{T}_1\) and \(\mathcal{T}_2\) at steady state:

• Let \(\beta_{V_i} = \beta_{V} = - 1 + \mu\), \(0< \mu \ll 0.01\), \(\beta_{Q_i} = \beta_{Q} = 0\), for all GFM-DERs. Then from 48 , \(|r_{V_i} Q_i^{\mathrm{ss}} - r_{V_j} Q_j^{\mathrm{ss}}| \approx |u_i^{\mathrm{ss}} - u_j^{\mathrm{ss}}| \leq \varepsilon \;\forall i,j\). Hence, equal per-unitized reactive power sharing is achieved among all GFM-DERs.

• Let \(\beta_{V_i} = \beta_{V} = 0, \beta_{Q_i} = \beta_{Q} = - 1 + \mu\), \(0< \mu \ll 0.01\), for all GFM-DERs \(i\). Then from 48 , \(|V_i^{\mathrm{ss}} - V_j^{\mathrm{ss}}| = | (\overline{V}+ u_i^{\mathrm{ss}}) - (\overline{V}+ u_j^{\mathrm{ss}}) | \approx |u_i^{\mathrm{ss}} - u_j^{\mathrm{ss}}| \leq \varepsilon\) for all GFM-DERs \(i,j\). Thus, the GFM-DERs buses have similar voltage magnitudes in steady state. Moreover, if \(\beta_{V_i}=\beta_V=0\) and \(\beta_{Q_i}=\beta_Q=-1+\mu\), with \(0<\mu\ll0.01\), for all GFM-DERs \(i\), then the solution \(x_i^\mathrm{ss}\) of 7 , and hence \(u_i^\mathrm{ss}\), satisfies \(u_i^\mathrm{ss}\approx0\). Therefore, the steady-state voltage magnitude satisfies \(V_i=V_{\mathrm{nom}}+u_i^\mathrm{ss}\approx\overline{V}\).

Therefore, with appropriate design choices, the secondary-level controller achieves equal reactive power sharing and voltage regulation among GFM-DERs at steady state.

3.2 Analysis of GFM-DER Frequency Control Loop↩︎

Next, we establish the stability properties of frequency dynamics 13 . Let \(\delta=[\delta_1,\dots,\delta_{n}]^\top\in\mathbb{R}^{n}, \omega=[\omega_1,\dots,\omega_{n}]^\top\in\mathbb{R}^{n}\) be the system states. Dynamics 13 , using 3 and 10 , can be written compactly as \[\begin{align} \label{eq:phase95omega95closedloop95vector} \dot{\begin{bmatrix} \delta \\ \omega \end{bmatrix}} &= \mathbf{A} \begin{bmatrix} \delta \\ \omega \end{bmatrix} + \mathbf{B} \begin{bmatrix} 0 \\ \Omega \end{bmatrix} + \mathbf{C} \begin{bmatrix} 0 \\ \boldsymbol{p} \end{bmatrix} + \mathbf{D} \begin{bmatrix} 0 \\ \widehat{V} \end{bmatrix}, \end{align}\tag{49}\] where, \(\Omega = [\Omega_1, \dots, \Omega_{n}]^\top \in \mathbb{R}^{n}\) with \(\Omega_i := (\overline{\omega} + r_{\omega_i} P_i^{\min})/\tau_{P_i}\), \(\boldsymbol{p} = [p_1, \dots, p_{n}]^\top \in \mathbb{R}^{n}\), \(\widehat{V} = [V^2_1, \dots, V^2_{n}]^\top \in \mathbb{R}^{n}\) and matrices \(\mathbf{A},\mathbf{B}, \mathbf{C}, \mathbf{D}\) are \[\begin{align} \mathbf{A} &:= \begin{bmatrix}0_{n} \; \mathbb{I}_{n} \\ \mathbf{G}_{\omega} \; \boldsymbol{\tau}_P \end{bmatrix} \mathbf{B}:= \begin{bmatrix} 0_{n} \;\;0_{n} \\ 0_{n} \;\;\mathbb{I}_{n} \end{bmatrix} \mathbf{C} := \begin{bmatrix} 0_{n} \;\;0_{n} \\ 0_{n} \; \;\boldsymbol{P} \end{bmatrix} \nonumber \\ \mathbf{D} &:= \begin{bmatrix} 0_{n} & 0_{n} \\ 0_{n} & \text{diag}((-G_{i i} r_{\omega_i})/\tau_{P_i}) \end{bmatrix}, \label{eq:ABC95matrix95omega95delta} \end{align}\tag{50}\] where, \(\boldsymbol{\tau}_{P} := \text{diag}\left(-1/\tau_{P_i}\right), [\mathbf{G}_\omega]_{i j} := -\frac{r_{\omega_i}B_{i j} V_{i} V_{j}}{\tau_{P_i}},\) \([\mathbf{G}_\omega]_{ii} := -\sum_{j} [\mathbf{G}_\omega]_{i j}, \boldsymbol{P} := \text{diag}(r_{\omega_i}(P_i^{\max} - P_i^{\min})/\tau_{P_i})\). Note that matrix \(\mathbf{G}_\omega\) has a marginal mode associated with the vector \(\mathbf{1}_{n}\). Since the phase angles are invariant under uniform shifts, we study the stability of 49 in relative coordinates. Let \(\xi\in\mathbb{R}^{n}_{>0}\) denote a normalized left null vector of \(\mathbf{G}_\omega\), i.e., \(\xi^\top \mathbf{G}_\omega=0\) and \(\xi^\top\mathbf{1}_{n}=1\), and define the projection matrix \(\Pi:=\mathbb{I}_{n} -\mathbf{1}_{n}\xi^\top\). The projected phase and frequency variables are given by \(\tilde{\delta}:=\Pi\delta, \; \tilde{\omega}:=\Pi\omega.\) Because \(\mathbf{G}_\omega\mathbf{1}_{n}=0\), we have \(\mathbf{G}_\omega\delta=\mathbf{G}_\omega\tilde{\delta}\). Applying the projection to 49 yields \[\begin{align} \dot{\tilde{\delta}} &= \tilde{\omega},\\ \dot{\tilde{\omega}} &= \mathbf{G}_\omega\tilde{\delta} + \Pi\boldsymbol{\tau}_P \tilde{\omega} + \Pi\boldsymbol{\tau}_P\mathbf{1}_{n}\omega_\xi + \Pi\Omega + \Pi\boldsymbol{P}\boldsymbol{p}\\ &+ \Pi\operatorname{diag}\left(-G_{ii}r_{\omega_i}/\tau_{P_i}\right)\widehat{V}, \end{align}\] where \(\omega_\xi:=\xi^\top\omega\) denotes the weighted average frequency component. In particular, if the active-power filter time constants are identical, i.e., \(\tau_{P_i}=\tau_P\) for all \(i\), then \(\boldsymbol{\tau}_P=-(1/\tau_P)\mathbb{I}_{n}\) and the term \(\Pi\boldsymbol{\tau}_P\mathbf{1}_{n}\omega_\xi\) vanishes. In this case, the projected dynamics close in \((\tilde{\delta},\tilde{\omega})\) as \[\label{eq:centered95freq95delta} \begin{align} \dot{\tilde{\delta}} &= \tilde{\omega}, \\ \dot{\tilde{\omega}} &= \textstyle \mathbf{G}_\omega\tilde{\delta} -\frac{1}{\tau_P}\tilde{\omega} + \Pi\left[\Omega + \boldsymbol{P}\boldsymbol{p}+\operatorname{diag}\left(\frac{-G_{ii}r_{\omega_i}}{\tau_{P_i}}\right)\widehat{V}\right]. \end{align}\tag{51}\]

Thus, the marginal absolute-angle direction is removed, and the stability analysis can be carried out on the disagreement subspace. The projected dynamics 51 , with \(\tilde{x}:= [\tilde{\delta}\;\tilde{\omega}]^\top\), can be expressed in a compact form as \[\label{eq:projected95compact} \begin{align} \dot{\tilde{x}} &= \mathbf{A}_\Pi \tilde{x} + \mathbf{B}_\Pi d(t), \\ \mathbf{A}_\Pi &:= \begin{bmatrix} 0_{n} & \mathbb{I}_{n} \\ \mathbf{G}_\omega & -\frac{1}{\tau_P}\mathbb{I}_{n} \end{bmatrix}, \; \mathbf{B}_\Pi := \begin{bmatrix} 0_{n}\\ \Pi \end{bmatrix}, \\ d(t) &:= \textstyle \Omega + \boldsymbol{P}\boldsymbol{p}(t) + \operatorname{diag}\left(\frac{-G_{ii}r_{\omega_i}}{\tau_P}\right)\widehat{V}(t). \end{align}\tag{52}\] We have the following result.

First note that \(\mathcal{S}\) is invariant under 52 . Indeed, if \(\tilde{x}\in\mathcal{S}\), then \(\xi^\top\dot{\tilde{\delta}} = \xi^\top\tilde{\omega}=0,\) and using \(\xi^\top\mathbf{G}_\omega=0\) and \(\xi^\top\Pi=0\), \(\xi^\top\dot{\tilde{\omega}} = \xi^\top\mathbf{G}_\omega\tilde{\delta} - \frac{1}{\tau_P}\xi^\top\tilde{\omega} + \xi^\top\Pi d(t) =0.\) Hence trajectories initialized in \(\mathcal{S}\) remain in \(\mathcal{S}\). Note that under the assumptions in Lemma 3, \(\mathbf{A}_\Pi\) is Hurwitz, and thus for any symmetric positive definite matrix \(\mathbf{Q}\) on \(\mathcal{S}\), there exists a symmetric positive definite matrix \(\mathbf{M}\) on \(\mathcal{S}\) such that \(\mathbf{A}_\Pi^\top\mathbf{M}+\mathbf{M}\mathbf{A}_\Pi=-\mathbf{Q}\) on \(\mathcal{S}\). Consider the Lyapunov function \(W(\tilde{x})=\tilde{x}^\top \mathbf{M}\tilde{x}.\) Along trajectories of the projected system, \(\dot{W} = \tilde{x}^\top \left( \mathbf{A}_\Pi^\top\mathbf{M} + \mathbf{M}\mathbf{A}_\Pi \right) \tilde{x} + 2\tilde{x}^\top \mathbf{M}\mathbf{B}_\Pi d(t) = -\tilde{x}^\top\mathbf{Q}\tilde{x} + 2\tilde{x}^\top \mathbf{M}\mathbf{B}_\Pi d(t).\) Applying Young’s inequality, with \(0 < \eta < 1,\) we get \(2\tilde{x}^\top \mathbf{M}\mathbf{B}_\Pi d(t) \leq \eta\lambda_{\min}(\mathbf{Q})\|\tilde{x}\|^2 + \frac{\|\mathbf{M}\|^2}{\eta\lambda_{\min}(\mathbf{Q})} \|\Pi d(t)\|^2.\) Thus, \(\dot{W} \leq -(1-\eta)\lambda_{\min}(\mathbf{Q})\|\tilde{x}\|^2 + \frac{\|\mathbf{M}\|^2}{\eta\lambda_{\min}(\mathbf{Q})} \|\Pi d(t)\|^2.\) Since, for all \(\tilde{x} \in \mathcal{S}\) we have, \(\lambda_{\min}(\mathbf{M})\|\tilde{x}\|^2 \leq W(\tilde{x}) \leq \lambda_{\max}(\mathbf{M})\|\tilde{x}\|^2,\) there exist constants \(c_1>0\) and \(c_2>0\) such that \(\dot{W} \leq -c_1 W + c_2\|\Pi d(t)\|^2.\) Applying the comparison lemma gives \(W(t) \leq e^{-c_1 t}W(0) + \frac{c_2}{c_1} \sup_{0\leq s\leq t}\|\Pi d(s)\|^2.\) Using the quadratic bounds on \(W\) yields \(\|\tilde{x}(t)\| \leq \textstyle \sqrt{\frac{\lambda_{\max}(\mathbf{M})}{\lambda_{\min}(\mathbf{M})}} e^{\frac{-c_1}{2}t}\|\tilde{x}(0)\| + \sqrt{\frac{c_2}{c_1\lambda_{\min}(\mathbf{M})}} \sup_{0\leq s\leq t}\|\Pi d(s)\|\). Thus, \(\kappa := \sqrt{\frac{\lambda_{\max}(\mathbf{M})}{\lambda_{\min}(\mathbf{M})}}, \lambda = \frac{c_1}{c_2}, \gamma := \sqrt{\frac{c_2}{c_1\lambda_{\min}(\mathbf{M})}}\). Hence, projected dynamics are input-to-state stable with respect to the projected input \(\Pi d\).

We prove the result in three steps. First, we show that the projected state remains in the disagreement subspace. Second, we apply the ISS estimate for the projected dynamics on this subspace from Theorem 4. Third, we combine this estimate with the scalar weighted-average frequency dynamics.

For any initial condition \((\delta(0),\omega(0))\), projected variables \((\tilde{\delta}(0), \tilde{\omega}(0))\) satisfy \(\xi^\top\tilde{\delta}(0) = \xi^\top\Pi\delta(0) =0, \xi^\top\tilde{\omega}(0)= \xi^\top\Pi\omega(0) =0.\) Thus, \(\tilde{x}(0)= \begin{bmatrix} \tilde{\delta}(0) \;\; \tilde{\omega}(0) \end{bmatrix} \in\mathcal{S}.\)

As shown in Theorem 4, \(\mathcal{S}\) is invariant under the projected dynamics 52 . Therefore, if \(\tilde{x}(0)\in\mathcal{S}\), then \(\tilde{x}(t)\in\mathcal{S}\) for all \(t\geq 0\). Hence, the projected trajectory evolves entirely on the disagreement subspace, where the restriction of \(\mathbf{A}_\Pi\) in 52 is Hurwitz by Lemma 3.

Consequently, the ISS estimate on \(\mathcal{S}\) applies to the projected system. Thus, from Theorem 4, there exist constants \(\kappa\geq 1\), \(\lambda>0\), and \(\gamma>0\) such that \[\|\tilde{x}(t)\| \leq \kappa e^{-\lambda t}\|\tilde{x}(0)\| + \gamma\sup_{0\leq s\leq t}\|\Pi d(s)\|, \quad t\geq 0,\] where, \(\Pi d\) is the projected input. Since \(\tilde{\omega}\) is a component of \(\tilde{x}\), it follows that \[\|\tilde{\omega}(t)\| \leq \|\tilde{x}(t)\| \leq \kappa e^{-\lambda t}\|\tilde{x}(0)\| + \gamma\sup_{0\leq s\leq t}\|\Pi d(s)\|.\]

It remains to analyze the weighted-average frequency component. Define \(\omega_\xi(t):=\xi^\top\omega(t), \; e_\xi(t):=\omega_\xi(t)-\omega_{\rm nom}\). Using the original frequency dynamics \(\dot{\omega} = \mathbf{G}_\omega\delta -\frac{1}{\tau_P}\omega + d(t)\), where, \(d(t)\) is as in 52 . Multiplying by \(\xi^\top\), we obtain \[\begin{align} \dot{\omega}_\xi & \textstyle = \xi^\top\dot{\omega} = \xi^\top\mathbf{G}_\omega\delta -\frac{1}{\tau_P}\xi^\top\omega + \xi^\top d(t). \end{align}\] Since \(\xi^\top\mathbf{G}_\omega=0\), this reduces to \(\dot{\omega}_\xi = -\frac{1}{\tau_P}\omega_\xi + \xi^\top d(t).\) Subtracting \(\omega_{\rm nom}\) from both sides gives \[\dot{e}_\xi = \textstyle -\frac{1}{\tau_P}e_\xi + \left( \xi^\top d(t)-\frac{1}{\tau_P}\omega_{\rm nom} \right).\] Define, \(d_\xi(t):= \xi^\top d(t)-\frac{1}{\tau_P}\omega_{\rm nom}\). Then \(\dot{e}_\xi = -\frac{1}{\tau_P}e_\xi+d_\xi(t)\). By the variation-of-constants formula, \[\textstyle e_\xi(t) = e^{-t/\tau_P}e_\xi(0) + \int_0^t e^{-(t-s)/\tau_P}d_\xi(s)ds.\] Taking absolute values yields \[\begin{align} |e_\xi(t)| &\leq \textstyle e^{-t/\tau_P}|e_\xi(0)| + \int_0^t e^{-(t-s)/\tau_P}|d_\xi(s)|ds \\ & \textstyle \leq e^{-t/\tau_P}|e_\xi(0)| + \sup_{0\leq s\leq t}|d_\xi(s)| \int_0^t e^{-(t-s)/\tau_P}ds \\ & \textstyle \leq e^{-t/\tau_P}|e_\xi(0)| + \tau_P\sup_{0\leq s\leq t}|d_\xi(s)|. \end{align}\] Finally, decompose the original frequency vector as \(\omega = \Pi\omega+\mathbf{1}_n\xi^\top\omega = \tilde{\omega}+\mathbf{1}_n\omega_\xi.\) Therefore, \[\omega-\omega_{\rm nom}\mathbf{1}_n = \tilde{\omega} + \mathbf{1}_n(\omega_\xi-\omega_{\rm nom}) = \tilde{\omega} + \mathbf{1}_n e_\xi.\] Using the triangle inequality, \[\begin{align} \|\omega(t)-\omega_{\rm nom}\mathbf{1}_n\| &\leq \|\tilde{\omega}(t)\| + \|\mathbf{1}_n e_\xi(t)\| \\ & = \|\tilde{\omega}(t)\| + \sqrt{n}\|e_\xi(t)\|. \end{align}\] Substituting the bounds on \(\tilde{\omega}(t)\) and \(e_\xi(t)\) gives \[\begin{align} \|\omega(t)-\omega_{\rm nom}\mathbf{1}_n\| &\leq \kappa e^{-\lambda t}\|\tilde{x}(0)\| + \gamma \sup_{0\leq s\leq t}\|\Pi d(s)\|\\ &+ \sqrt{n}e^{-t/\tau_P}|e_\xi(0)| + \sqrt{n}\tau_P \sup_{0\leq s\leq t}|d_\xi(s)|. \end{align}\] This proves the stated bound on the frequency. If, in addition, \(\sup_{t\geq 0}\|\Pi d(t)\|\leq \Delta_d, \sup_{t\geq 0}|d_\xi(t)|\leq \Delta_\xi\), then the preceding estimate immediately implies \[\begin{align} \|\omega(t)-\omega_{\rm nom}\mathbf{1}_n\| &\leq \kappa e^{-\lambda t}\|\tilde{x}(0)\| + \sqrt{n}e^{-t/\tau_P}|e_\xi(0)| \\ &+ \gamma\Delta_d + \sqrt{n}\tau_P\Delta_\xi. \end{align}\] Thus, the frequency deviation from nominal synchronization is input-to-state practically stable with respect to the disagreement input \(\Pi d\) and the weighted average-frequency mismatch \(d_\xi\).

4 Conclusion↩︎

This paper presented a mathematical stability analysis of a sampled-data optimization-based secondary controller for networks of inverter-interfaced DERs. The analysis treated the controller model as given and studied the nonlinear closed-loop dynamics induced by sampled measurements, constrained optimization updates, and interpolation-based actuation. For the GFM-DER voltage loop, we established large-signal boundedness of the voltage and filtered reactive-power dynamics within a certified operating region. We also characterized positive steady-state operating points and showed how the optimization-induced consensus condition connects voltage regulation with equal per-unitized reactive power sharing. For the phase-frequency dynamics, we used a projected representation to remove the marginal absolute-angle mode and established input-to-state stability with respect to active-power mismatch. These results provide closed-loop guarantees for optimization-based secondary control beyond small-signal or purely continuous-time analyses. Future work will focus on relaxing some of the technical assumptions used in the analysis, including identical active-power filter time constants and exact tracking by the inner GFL-DER control loop.

4.1 Supporting Lemmas↩︎

For any \(x \in \mathbb{R}^{n}\) we have \(x^\top (\Sigma \mathcal{B} + \mathcal{B}^\top \Sigma)x = 2x^\top \Sigma \mathcal{B} x\), since \(x^\top \Sigma \mathcal{B} x\) is a scalar and equal to its transpose. Since, \(V_i > 0\), for all \(i\), thus, \([\Sigma]_{i i} > 0\) for all \(i\). Let \(y = \Sigma^{1/2} x\) (well-defined since \(\Sigma\) is a positive diagonal matrix). Then \(x = \Sigma^{-1/2} y\) and \(2 x^\top \Sigma \mathcal{B} x = 2 y^\top ( \Sigma^{-1/2} \mathcal{B} \Sigma^{-1/2} ) y.\) Define \(\nu := \Sigma^{-1/2} \mathcal{B} \Sigma^{-1/2}\). Note that \(\nu = \sigma^\top \mathcal{B} \sigma\) with \(\sigma = \Sigma^{-1/2}\), so \(\nu\) is congruent to \(\mathcal{B}\), and \(\nu \succeq 0\) if and only if \(\mathcal{B} \succeq 0\). Therefore, \(\Sigma \mathcal{B} + \mathcal{B}^\top \Sigma \succeq 0 \iff \mathcal{B} \succeq 0\), which is indeed the case due to the Gershgorin Circle Theorem.

Due to the consensus constraints, every feasible point has the form \(x_s=c\mathbf{1}_{n}\). Hence, the optimization reduces to the scalar problem \(\min_{c\in\mathcal{I}(V,Q)} \|c\mathbf{1}_{n}-\alpha(V,Q)\|_2^2\), where \(\mathcal{I}(V,Q) := \textstyle \bigcap_{i=1}^{n} [m_i(V,Q)-V_\Delta, m_i(V,Q)+V_\Delta]\) and \(m_i(V,Q) := (1+\beta_{Q_i})r_{V_i}Q_i^{\mathrm{avg}}(t_s) - \beta_{V_i}(\overline{V}- V_i (t_s)).\) Writing \(\bar\alpha(V,Q):= \frac{1}{n}\mathbf{1}^\top\alpha(V,Q),\) we have \(\|c\mathbf{1}-\alpha(V,Q)\|_2^2 = n(c-\bar\alpha(V,Q))^2 + \|\alpha(V,Q)-\bar\alpha(V,Q)\mathbf{1}\|_2^2.\) Therefore, the optimizer is \(c_s(V,Q)=\Pi_{\mathcal{I}(V,Q)}(\bar\alpha(V,Q))\). Let \(\ell(V,Q):=\max_i\{m_i(V,Q)-V_\Delta\}, r(V,Q):=\min_i\{m_i(V,Q)+V_\Delta\}.\) Then \(\mathcal{I}(V,Q)=[\ell(V,Q),r(V,Q)]\). Since each \(m_i(V,Q)\) is continuous, \(\ell\) and \(r\) are continuous. Moreover, \(\bar\alpha\) is continuous because \(\alpha\) is continuous. Thus \(c_s(V,Q) = \min\{\max\{\bar\alpha(V,Q),\ell(V,Q)\},r(V,Q)\}\) is continuous. Hence, the solution \(x_s(V,Q)=c_s(V,Q)\mathbf{1}_{n}\) is continuous on the domain where \(\mathcal{I}(V,Q)\) is nonempty.

Define the positive diagonal matrix \(R_\omega := \operatorname{diag}(r_{\omega_1},\dots,r_{\omega_{n}})\). Next, define the symmetric weighted Laplacian \(L_V\) associated with the GFM-DER interaction graph by \([L_V]_{ij} = B_{ij}V_i V_j, \;i\neq j, \;and \; [L_V]_{ii} = -\sum_{j\neq i}B_{ij}V_i V_j.\) Since \(B_{ij}=B_{ji}<0\) on every edge and \(V_i>0\), \(L_V\) is a symmetric positive semidefinite weighted Laplacian. Since the GFM-DER interaction graph is connected, \(\operatorname{ker}(L_V)=\operatorname{span}\{\mathbf{1}_{n}\}.\) Using the definition of \(\mathbf{G}_\omega\), we can write \(\mathbf{G}_\omega = -\frac{1}{\tau_P}R_\omega L_V\). Because \(R_\omega\succ 0\), the matrix \(R_\omega L_V\) is similar to a symmetric positive semidefinite matrix. Indeed, \(R_\omega L_V = R_\omega^{1/2} \left( R_\omega^{1/2}L_VR_\omega^{1/2} \right) R_\omega^{-1/2}\). Thus, \(R_\omega L_V\) has real nonnegative eigenvalues. Consequently, \(\mathbf{G}_\omega\) has real nonpositive eigenvalues. Since the graph is connected, \(\mathbf{G}_\omega\) has one zero eigenvalue corresponding to the uniform-angle direction and all remaining eigenvalues are strictly negative. That is, \(\lambda_1(\mathbf{G}_\omega)=0, \;\lambda_k(\mathbf{G}_\omega)<0, \;k=2,\dots,n.\) The zero eigenvalue corresponds to the direction \(\mathbf{1}_{n}\), which is removed by the projection. Hence, on the disagreement subspace \(\mathcal{S}\), only the modes associated with \(\lambda_k(\mathbf{G}_\omega)<0\) remain.

To show that this implies the Hurwitz property of \(\mathbf{A}_\Pi\) on the disagreement subspace, we now lift the modal properties of \(\mathbf{G}_\omega\) to the second-order phase-frequency dynamics. Since \(\mathbf{G}_\omega\) is similar to a symmetric matrix, it is diagonalizable and has real eigenvalues. Let \(v_k\) be an eigenvector of \(\mathbf{G}_\omega\) associated with a disagreement eigenvalue \(\lambda_k<0\), so that \(\mathbf{G}_\omega v_k=\lambda_k v_k.\) Consider an eigenpair \(\left(s, \begin{bmatrix} \phi\\ \psi \end{bmatrix} \right)\) of \(\mathbf{A}_\Pi\) on the disagreement subspace. Then \(\begin{bmatrix} 0_{n} & \mathbb{I}_{n}\\ \mathbf{G}_\omega & -\frac{1}{\tau_P}\mathbb{I}_{n} \end{bmatrix} \begin{bmatrix} \phi\\ \psi \end{bmatrix} = s \begin{bmatrix} \phi\\ \psi \end{bmatrix}.\) The first block row gives \(\psi=s\phi\). Substituting this relation into the second block row gives \(\mathbf{G}_\omega\phi-\frac{1}{\tau_P}\psi=s\psi.\) Using \(\psi=s\phi\), we obtain \(\mathbf{G}_\omega\phi-\frac{s}{\tau_P}\phi=s^2\phi\), or equivalently, \(\mathbf{G}_\omega\phi = \left(s^2+\frac{s}{\tau_P}\right)\phi\). Thus, for each eigenvalue \(\lambda_k\) of \(\mathbf{G}_\omega\), the corresponding eigenvalues \(s\) of \(\mathbf{A}_\Pi\) satisfy \(s^2+\frac{1}{\tau_P}s-\lambda_k=0.\) On the disagreement subspace, \(\lambda_k<0\). Hence, \(\frac{1}{\tau_P}>0, \;-\lambda_k>0.\) Therefore, by the second-order Routh-Hurwitz criterion, the polynomial \(s^2+\frac{1}{\tau_P}s-\lambda_k\). has both roots in the open left-half complex plane. Consequently, every eigenvalue of \(\mathbf{A}_\Pi\) associated with a disagreement mode has a strictly negative real part. The only eigenvalue of \(\mathbf{G}_\omega\) that is not strictly negative is \(\lambda_1=0\), which corresponds to the uniform-angle direction \(\mathbf{1}_{n}\). For this mode, the characteristic equation becomes \(s^2+\frac{1}{\tau_P}s=0,\) whose roots are \(s=0, s=-\frac{1}{\tau_P}.\) The zero root corresponds to the absolute-angle mode, which is removed by projection onto \(\mathcal{S}\). Therefore, no marginal mode remains on \(\mathcal{S}\), and all eigenvalues of \(\mathbf{A}_\Pi\) restricted to \(\mathcal{S}\) have strictly negative real parts. Hence, \(\mathbf{A}_\Pi\) restricted to \(\mathcal{S}\) is Hurwitz.

References↩︎

---

1. \(^{\dagger}\)Department of Mechanical Science and Engineering, University of Illinois at Urbana-Champaign, IL, USA {vkhatana}@{illinois.edu}, \(^{\ddagger}\)Department of Electrical Engineering, Indian Institute of Science, Karnataka, India {schakraborty@iisc.ac.in}, \(^{\wr}\)Department of Electrical and Computer Engineering, University of Minnesota, MN, USA {murtis@umn.edu}. The research conducted with the support of the United States Department of Energy via grant DE-CR\(0000040\).↩︎

2. Sampling instants satisfy \(t_{s+1}:=t_s+\Delta_s+\Delta\) for \(s\in\{0,1,2,\dots\}\)↩︎