Beyond Binary Controllability: Reachability Margins for Power Systems via Scalable Screening

Introduction and Paper Contributions

Power systems operate around economic dispatch and optimal power flow (OPF) setpoints. Their dynamics are studied locally through linearized versions of the nonlinear differential-algebraic equation (DAE) models. For a multi-machine power network, DAEs couple dynamic states such as frequencies with algebraic, network power flow constraints and variables (such as voltages and angles) [1, 2, 3]. Much of the literature removes the exact dynamic-algebraic coupling to obtain ODE models, simplifying analysis and enabling known control methods primarily developed for ODEs.

Paper’s objective and scope. The paper’s objective is to study this screening problem: given an OPF setpoint, a schedule of load and renewable injections, a linearized DAE model, a control horizon of $T$ seconds, bounds on the control input deviations, and an arbitrary terminal target $x^{*}_{\mathrm{T}}$, a prescribed value of the deviation $\Delta x(T)$, compute a scalable margin $\delta\le0$ that quantifies the difficulty of this engineering problem. A negative $\delta$ certifies that the target is unreachable and equals the negative distance to the terminal reachable set, while $\delta= 0$ means the target is not flagged. The idea is presented in Figure 1.

For multi-machine power network models, the Kalman controllability test is numerically unreliable even when the PBH test certifies controllability [4], and models controllable in the rank sense remain practically uncontrollable along weak directions where the required control effort is prohibitive [5]. In this context, absolute controllability is not very meaningful, as some state targets are not interesting from the power system side and hence their reachability is not relevant. However, this investigation becomes pertinent in renewables-rich operation, because the operating point, the local Jacobians, and

the effective control authority all vary with dispatch and network stress. With more renewables, synchronous machines and their controls are displaced, which reduces the available control authority and typically shrinks the reachable space under bounded inputs [6]. For the problem we study in this paper, some terminal state requests may remain easy to realize, while others become difficult because the horizon $T$ can be short, the controller bounds are tight, or the algebraic network closure prevents the full terminal state from matching the request. A useful screening tool therefore needs to depend explicitly on the inputs described in Figure 1.

The central modeling point is that the terminal target lives on the full DAE-consistent manifold. Reduced ODE models remain useful, but they hide failure mechanisms carried by algebraic variables such as voltages and network-balance quantities. In contrast, the target in this paper may simultaneously request acceptable recovery of frequencies, internal EMF variables, and other dynamic states. The resulting tool herein screens frequency control and voltage control reachability in one scalable computation. The relevant literature is reviewed next.

Brief literature review. The relevant literature spans DAE system theory, scalable reachability computation, and power system applications. The studies [7, 8] characterize reachable and controllable subspaces through matrix-pencil structure and DAE indices. The authors in this extensive survey [9] systematize the different controllability notions that arise once consistency constraints are present, and the study [10] develops reachable set estimation for singular systems with explicit treatment of admissible initial conditions. The study [11] treats high-index linear DAE reachability by separating slow and algebraic dynamics through consistency enforcing constructions. These contributions explain how DAEs differ from ODEs and why algebraic constraints must be handled explicitly, but they do not formulate a target-specific screening metric. In the scalable reachability analysis more broadly, the study [12] propagates uncertain linear reachable sets with zonotopes (convex polytopes). The authors in [13] develop support function propagation for linear systems by repeatedly maximizing linear functionals rather than enumerating full sets. The work [14] uses conservative linearization to bound nonlinear reachable sets, and the study [15] combines zonotopes to improve computational scalability for large systems. These works make clear why support functions are important: they encode the directional action of a high-dimensional reachable set, support separating-hyperplane arguments, and avoid explicit construction of a complicated set in every computation. Our work here, instead of propagating an entire reachable tube, uses the support function of the terminal reachable set to compute an exact target-to-set shortfall for a prescribed terminal request.

In power systems, reachable set ideas mainly support transient stability assessment and formal verification. The study [16] investigates transient stability assessment through reachable set computations. The study [17] investigates reachability analysis for nonlinear differential-algebraic systems and shows how set-propagation ideas adapt to DAE structure. The dissertation [18] develops formal and compositional reachability analysis for power systems under uncertainty, including decomposition strategies for scalability, conservative set propagation around operating points, and reachable set based verification. Another dissertation [19] extends this work and develops transient stability assessment, region of attraction estimation, LPV-based approximations, and reachability-assisted control design for large power networks; see also studies [20, 21] that delineate the contributions of the dissertation [19]. These studies compute (or over-approximate) reachable sets to certify safety, stability, or control performance. They do not derive an exact Euclidean distance from a prescribed full terminal target to the terminal reachable set of a linearized power system DAE model. Our recent work [22] paper develops a scalable, target-aware reinterpretation of the PBH controllability test by using dual infeasibility vectors to certify unreachable terminal states. This study extends our work with a focus on DAEs, power systems, and continuous-time (CT) dynamic models compared to the discrete-time (DT) ODEs in [22].

Paper Contributions. The paper contributions are:

• We formulate terminal reachability directly on the DAE-consistent manifold, so the target may include both differential and algebraic power system variables.

• For index-1 linearized power system DAEs with bounded controls, we derive an explicit support function representation of the full terminal reachable set and a screening margin $\delta$ that equals the negative Euclidean distance to that set when the target is infeasible. This quantifies the difficulty in steering the prescribed system to the target.

• We present three computational routes tailored to screening: a quadrature-based conic approximation, a projected subgradient method, and a bundle-type method. The methods are endowed with different properties and scale nicely. Thorough case studies are included to corroborate the theory.

Section II introduces the model and the screening problem. Section III derives the margin and its geometric meaning. Section IV discusses how to compute the metric efficiently. Section V presents case studies and Section VI concludes the paper.

Symbol	Typical contents	Role in screening
$x_d$	rotor angles, frequencies, internal EMFs, governor states, inverter-controller states	dynamic quantities that evolve over $[0,T]$
$x_a$	bus voltages, phase angles, nodal balance variables, algebraic network quantities	algebraic quantities that must close consistently at the terminal time
$u$	AGC/secondary-control actions, inverter setpoints, FACTS/HVDC modulation, controllable load actions	bounded actuation used to move the terminal state
$x_T^\star$	selected terminal values from both $x_d$ and $x_a$	screening target to be tested for reachability
$\xi_T$	target after removing the free response	forced displacement that the controls must generate

System Model and Problem Formulation

We start from the nonlinear DAE model [1, 23]

generator physics: $\dot{x}_d(t)=f_d(x_d(t),x_a(t),u(t)),$

network power flow: $0=f_a(x_d(t),x_a(t),u(t)),$

where $x_d(t)\in\mathbb{R}^{n_d}$ collects the differential states, $x_a(t)\in\mathbb{R}^{n_a}$ collects the algebraic states, and $u(t)\in\mathbb{R}^{m}$ is the control input. Around an equilibrium $(\bar{x}_d,\bar{x}_a,\bar{u})$, the linearized model is

or

with

Tab. I summarizes the main variables in the power system model and the screening tool. We note the following: (i) in many classical models $B_a=0$ because generators’ control inputs do not appear in the algebraic power flow constraints, but we retain $B_a$ because several grid control strategies (i.e., demand response or renewables curtailment) act directly on the algebraic constraints. (ii) Likewise, loads and renewable injections do not appear as separate variables in (2) as their influence is absorbed into the operating point and the Jacobian matrices. There is no loss of generality here, as the model can still accommodate their deviations. The math is simpler without them, and hence they are omitted.

Assumption 1 (Index-1 regularity). Matrix $A_{aa}$ is nonsingular, and the descriptor pair $(E,A)$ is regular and impulse-free.

Under Assumption 1, which is common in the power systems literature and is satisfied for a host of standard benchmark models near regular operating points [3, 24], the algebraic state is given by

Substituting into (2a) yields reduced dynamics

with $A_r:=A_{dd}-A_{da}A_{aa}^{-1}A_{ad}$ and $B_r:=B_u-A_{da}A_{aa}^{-1}B_a$. Matrix $C_r$ is the linearized algebraic-recovery map: it tells us how a perturbation in the differential state instantaneously propagates to the algebraic variables through the power flow equations. Given symmetric controller limits $\bar{u}\in\mathbb{R}^{m}_{>0}$ on the deviations of the generators’ control actions $u(t)$ (recall that $u(t)$ is the deviation in the setpoints), we use the admissible input class $U_T:=\{u\in L_\infty([0,T],\mathbb{R}^m)\mid\lvert u(t)\rvert\le\bar{u}\ \text{a.e. on }[0,T]\}$. By (2), the algebraic closure holds for all $t$; at the terminal time it is determined by the admissible endpoint input value. We therefore introduce an explicit terminal control variable $u_T\in\mathbb{R}^m$ satisfying $\lvert u_T\rvert\le\bar{u}$. This avoids committing the theory to a specific endpoint convention. If the implementation is continuous at $T$, one identifies $u_T=u(T)$; if a sampled implementation uses a pre-jump endpoint convention, one identifies $u_T=u(T^-)$. The analysis below only needs the admissible endpoint value $u_T$. For any admissible endpoint control value $u_T$, define

We also use the admissible terminal manifold $\mathcal{M}_T:=\bigcup_{\lvert u_T\rvert\le\bar{u}}\mathcal{M}_T(u_T)$. Let $x_{d,0}$ be a consistent initial differential state and let $x_T^\star=[x_{d,T}^\star,x_{a,T}^\star]^\top$ be the desired terminal state, or a terminal state value whose (un)reachability we seek to quantify. The free terminal response is defined as $x_T^{\mathrm{fr}}$ and the forced terminal operator induced by $(u(\cdot),u_T)$ is defined as $\Phi_T(u,u_T)$. We also define $P:=[I\ \ C_r]$ and $G:=[0\ \ D_r]$. We can now write

Accordingly,

By construction, $\mathcal{R}_T^{\mathrm{full}}\subseteq\mathcal{M}_T$. Therefore, any target outside $\mathcal{M}_T$ is immediately impossible because it violates terminal algebraic consistency for every admissible endpoint control. We assume throughout that $\mathcal{R}_T^{\mathrm{forc}}$ is nonempty, compact, and convex. For a desired terminal target $x_T^\star$, let $\xi_T:=x_T^\star-x_T^{\mathrm{fr}}$ define the target after removing the free response. Then $x_T^\star\in\mathcal{R}_T^{\mathrm{full}}$ if and only if $\xi_T\in\mathcal{R}_T^{\mathrm{forc}}$. The problem formulation is given next.

Margin Formulation on DAE Manifolds

In this section, we present the simple theoretical result that allows us to compute the margin $\delta$. To do so, we need to define support functions. Readers familiar with support functions can skip the next paragraph.

For a nonempty set $S\subset\mathbb{R}^n$, the support function $h_S(y):=\sup_{x\in S}y^\top x$ gives the maximum projection of the set along the direction $y$. The support function is a convex function that maps a vector to the maximum inner product over that set, essentially measuring the distance of the set’s supporting hyperplane from the origin. For example, if $y$ places weights on frequencies, internal EMF states, voltages, phase angles, and other terminal coordinates, then $h_{\mathcal{R}_T^{\mathrm{forc}}}(y)$ is the largest terminal displacement that the admissible controls can produce along that weighted combination of quantities.

This interpretation is standard in reachability analysis, where support functions summarize high-dimensional reachable sets through linear functionals, support separating hyperplane arguments, and enable scalable outer approximations without explicit set gridding [13, 15]. This same idea makes the screening problem meaningful and scalable. Rather than constructing the full terminal reachable set, the paper works with its directional action and recovers an exact target-to-set shortfall from that information. Having defined support functions and introduced the model and the terminal manifold, the next step is to characterize this directional action explicitly for the index-1 power system DAE. The following lemma provides this result.

The proof is given in Appendix A. Targets outside $\mathcal{M}_T$ can be screened out immediately. The next result handles the nontrivial case $x_T^\star\in\mathcal{M}_T$, where algebraic consistency is possible but reachability is not guaranteed. This result addresses Problem 1.

The proof is given in Appendix B. Theorem 1 shows that the screening margin is a one-sided control shortfall metric. It returns zero when the target is not certified as unreachable and becomes negative precisely when the target lies outside the full terminal reachable set. The direction $y^\star$ indicates where the terminal request most strongly exceeds the attainable directional capacity. In power system terms, the dominant entries of $y^\star$ localize the shortfall at specific machines (frequency entries) or buses (voltage and angle entries), and $x_T^{\mathrm{fr}}+\xi_T-\lvert\delta\rvert y^\star$ is the nearest attainable terminal state, hence a minimally relaxed target. Typical requests include full or fractional recovery of a disturbance and redispatch toward a neighboring equilibrium, with $T$ set by the secondary-control window and $\bar{u}$ by governor and inverter headroom.

The proof is given in Appendix C. Corollary 1 formalizes the sense in which the metric gets worse for a more unreachable instance, i.e., if a design change shrinks the terminal reachable set, then the margin cannot improve for the same terminal request.

Computing the Margin

The simple results in the previous section reduce screening to repeated evaluation of the support function in Lemma 1 and the optimization problem in Theorem 1. In this section, we develop three approaches to solve the optimization. Approach I approximates the exact margin by quadrature rules and then solves a finite-dimensional convex conic problem. Approach II keeps the CT model intact and applies projected subgradient ascent to the dual objective. Approach III uses the same CT dual formulation but stabilizes the search with a bundle-type model.

Approach I: Quadrature and Conic Optimization

Approach I targets a high-fidelity approximation of (6). The key observation is that the support function in Lemma 1 separates into a CT dynamic integral and a finite-dimensional terminal algebraic term. We therefore approximate only the integral by Gauss–Legendre quadrature [25]. Let $\{\hat{\zeta}_j,\hat{w}_j\}_{j=1}^{N_q}$ be the Gauss–Legendre nodes and weights on $[-1,1]$, the roots of the degree-$N_q$ Legendre polynomial and their quadrature weights, and define, for $j=1,\ldots,N_q$, the following terms

Then

and the margin is approximated by

Introducing slack variables $t_j\in\mathbb{R}^m$ and $s\in\mathbb{R}^m$ yields:

This is a finite-dimensional second order cone program (SOCP) and may therefore be solved globally by standard conic optimization solvers. It is always feasible since $y=0$, $t_j=0$, and $s=0$ satisfy all constraints and its objective is bounded below by $-\|\xi_T\|_2$, so a minimizer exists. Its main cost comes from repeatedly evaluating the matrix exponentials embedded in the quadrature matrices $Q_j$ and from the number $N_q$ of quadrature nodes. This formulation can become expensive if many operating points or many target requests must be screened.

Approach II: Projected Subgradient Method

The first approach, while benchmark worthy, can be computationally challenging for large-scale power networks. To that end, the second approach works directly with the CT dual and never discretizes the DAE model at the formulation level. Define

Then Theorem 1 $\delta(x_T^\star)=-\max_{\|y\|_2\leq1}J(y)$. The objective $J(y)$ is concave but generally nonsmooth because of the $\ell_1$-norm terms. Consequently, one natural computational route is a projected subgradient method, in which a subgradient ascent step is taken on $J(y)$ and the iterate is then projected back onto the feasible set $\{y\in\mathbb{R}^{n_d+n_a}\mid\|y\|_2\leq1\}$; see, e.g., [26]. A second option is to use a bundle-type method [27], which retains several past subgradients and function values to construct a stabilized piecewise-linear model of the nonsmooth objective and often yields more reliable practical progress than plain subgradient steps; see [27]. This section focuses on a projected subgradient method. This route keeps the DAE model continuous in time and moves all optimization to the finite-dimensional variable $y$. For a given $y$, let $v(\tau)$ solve the adjoint reduced dynamics

Then $v(\tau)=e^{A_r^{\mathsf T}(T-\tau)}P^{\mathsf T}y$, and a valid extreme control is

with any subgradient choice at zero. The extreme terminal point is

so a valid subgradient of $J(y)$ is

In the numerical implementation, $J(y)$ and $g(y)$ are evaluated with $N_q$ quadrature nodes, so the projected subgradient routine is indexed by the triplet $(N_q,K,\alpha_k)$, where $K$ is the iteration budget.

A projected subgradient ascent step is therefore

where $\Pi_{\|\cdot\|_2\leq1}$ denotes Euclidean projection onto the unit ball and $\alpha_k>0$ is a stepsize [26]. Standard diminishing stepsize rules such as $\alpha_k=\alpha_0/\sqrt{k+1}$ keep the iteration stable in the nonsmooth setting. If $\{y^k\}_{k=0}^K$ denotes the generated iterates, a practical estimate of the screening margin is $\delta_{\mathrm{sg}}^{(K)}(x_T^\star):=-\max_{0\leq k\leq K}J(y^k)$. Whenever $J(y^k)>0$ for some iterate, the corresponding estimate is strictly negative and terminal infeasibility is certified.

Approach III: Bundle-Type Method

While projected subgradient ascent uses only the most recent subgradient information, a bundle-type method [27] retains several past function values and subgradients in order to build a local piecewise-linear model of the concave nonsmooth objective $J(y)$. This is often more stable in practice, especially when the dual landscape contains many changes in active signs induced by the $\ell_1$-norm terms. The key idea is to replace the true objective locally by the best available collection of affine upper models and then maximize this model with an explicit proximal regularization over the Euclidean unit ball.

Suppose that at bundle points $y^{(i)}$, $i=1,\ldots,p$, we have already computed objective values $J(y^{(i)})$ and valid subgradients $g^{(i)}\in\partial J(y^{(i)})$. Since $J(y)$ is concave, every subgradient gives the global affine upper model

The resulting piecewise-linear model is

Rather than maximizing $\widehat{J}_p$ directly, the method stabilizes the search around a current center $y^c$ by solving

where $\rho>0$ is a proximal parameter. Problem (15) is always feasible, for example $y=y^c$ with $\eta=\widehat{J}_p(y^c)$, and eliminating $\eta$ leaves a strongly concave maximization over a compact set, so its solution exists and is unique. The optimizer of (15) produces a new trial point. At that point, one evaluates $J(y)$ and the subgradient (12) using exactly the same adjoint computation as in Approach II, appends the new cut to (15), and repeats. The resulting estimate is

In practice, the bundle is truncated to at most $p_{\max}$ cuts, while $J$ and $g$ are evaluated at each trial point with a composite trapezoidal rule on a uniform grid with $N_q$ nodes. The ensuing case studies tests all of these methods.

Numerical Experiments

We apply the results of Lemma 1 and Theorem 1 to linearized DAE models obtained from the IEEE 9, 39, 57, and 200 bus benchmarks under several renewable penetration and load levels. For each renewable penetration and load scaling pair, renewables enter as active power injections that offset part of the demand at the larger load buses; we solve the resulting nonlinear equilibrium, form the corresponding descriptor linearization $(E,A,B)$, set $\bar{u}$ as a fraction of the equilibrium input magnitudes that shrinks with penetration and loading, generate the post-step initial condition from a small load increase and a simultaneous renewable drop, compute the corresponding free terminal state $x_T^{\mathrm{F}}$ over $T=10$ s, and then construct two DAE-consistent full state targets along a weak direction $y_{\mathrm{cert}}$ obtained from a Gramian or PBH analysis of $(A_r,B_r)$, placed at 0.5 and 1.25 times the directional support value. The latter is unreachable by construction since the dual objective in (8) is positive at $y_{\mathrm{cert}}$, and all screening studies below use these certified unreachable targets. An implementation of the three approaches, together with the linearized benchmark models, is available at [28]. Table III lists the algorithmic parameters.

Appr.	Settings
I	$N_q=40$, SOCP solver tolerance $10^{-8}$
II	$N_q=80$, $K=300$, $\alpha_k=\alpha_0/\sqrt{k+1}$ with $\alpha_0=1.0$, backtracking factor $\beta=0.5$, backtracking cap 8, improvement tolerance $10^{-10}$, normalized subgradient, and stall restart.
III	$N_q=80$, $K=120$, proximal parameter $\rho=2.0$, bundle cap $p_{\max}=15$, and optimization tolerance $10^{-8}$.

Method Comparison

We begin by comparing the three methods developed in Section IV across all benchmark systems. For each operating point listed in Table II, the same full state terminal request $x_T^*$ is screened using

Case	$(n,m)$	LoadScale	Ops.	$\lVert\bar{u}\rVert_\infty$
IEEE 9	$(36,6)$	1.00, 1.10, 1.20	9	0.040 to 0.202
IEEE 39	$(138,20)$	1.00, 1.05, 1.10	7	0.865 to 4.274
IEEE 57	$(156,14)$	1.00, 1.05, 1.10	9	0.090 to 0.564
IEEE 200	$(694,98)$	1.00, 1.05, 1.10	9	0.077 to 0.510

Approaches I–III and the resulting margins are recorded together with the corresponding computational times. Figure 2 compares the margins obtained from the three computational approaches across all benchmark systems and operating points. Approach I serves as the conic benchmark because it directly solves (9) through the resulting SOCP (10). As expected, margins obtained with Approaches II and III remain close to those of Approach I across all systems, indicating that the nonsmooth dual optimization in (11) captures the same separating geometry described in Theorem 1. This confirms that the directional shortfall for power systems identified by the dual variable $y^*$ is robust with respect to the numerical approach used to compute it. The underlying limitation is therefore structural (induced by actuator placement, network topology, and algebraic coupling) rather than some artifact of numerical discretization.

Convergence and Scalability

Table IV summarizes the benchmark-relative convergence and computational scalability of the three screening methods. Approach I is treated as the benchmark reference because for a sufficiently accurate quadrature rule it solves the conic approximation of the screening problem to high accuracy. Approaches II and III are interpreted as first-order approximations to this benchmark rather than competing exact solvers. For each sampled infeasible instance, the certificate produced by Approach I is denoted by $\delta^{\mathrm{I}}$. Approaches II and III are declared to have reached the benchmark once their best-so-far certificate enters the tolerance band

Case	I [s]	II to I [s] (speedup)	II med. iters (hit rate)	III to I [s] (speedup)	III med. iters (hit rate)
IEEE 9	0.757	0.005 (156.10×)	6.0 (100%)	0.136 (5.58×)	3.0 (100%)
IEEE 39	0.760	0.221 (3.44×)	100.0 (20%)	0.315 (2.41×)	15.0 (100%)
IEEE 57	0.622	0.359 (1.73×)	273.0 (20%)	0.103 (6.05×)	5.0 (100%)
IEEE 200	18.172	3.456 (5.26×)	123.0 (80%)	0.329 (55.19×)	5.0 (40%)

with $\varepsilon_{\mathrm{abs}}=10^{-5}$ and $\varepsilon_{\mathrm{rel}}=10^{-2}$. This is appropriate here because the objective is not exact iterate matching but practical agreement with the benchmark certificate at a scale commensurate with the certificate magnitude. The reported hit rate is the fraction of sampled instances that enter this band before the prescribed iteration cap. The results show a clear separation between the two first-order schemes. Approach III reaches the benchmark tolerance band in very few iterations on the small and medium systems: the median count is 3 iterations for IEEE 9 and 5 iterations for IEEE 57, with 100% hit rate in both cases. Even on IEEE 39, Approach III reaches the benchmark within 15 iterations at 100% hit rate. The corresponding wall clock times remain below the benchmark runtime in all three systems. Approach II is also faster than Approach I in wall clock time, but its convergence is markedly less reliable under the same caps, hit rate drops to 20% on IEEE 39 and IEEE 57, and the median benchmark-reaching iteration count is much larger when convergence occurs.

For the large system (IEEE 200), Approach III attains the benchmark band in median 5 iterations and achieves a speedup exceeding one order of magnitude in wall clock time although its hit rate drops to 40% under the fixed cap $K_{\mathrm{III}}=120$. Approach II attains a higher hit rate of 80% on this case but requires substantially more iterations and remains slower than Approach III. This indicates that the bundle-type updates in Approach III are substantially more effective at extracting a high-quality screening certificate early whereas the projected subgradient method can require many more iterations to achieve the same benchmark-relative accuracy. Overall, these results support the simple intended message of the paper. Approach I provides a benchmark-quality certificate, while Approaches II and III provide scalable approximations whose practical value is determined by how rapidly they approach that benchmark. Among the two first-order methods, Approach III gives the strongest tradeoff between runtime and benchmark-relative accuracy, especially when only a small iteration budget is available.

Effect of Operating-Point Stress on Reachability

We now study how operating-point stress affects the reachability margin. Stress is induced through simultaneous changes in renewable penetration and loading level, which modify both the Jacobian matrices and the effective control directions. Figure 3 shows the change of the screening margin as the operating point becomes more stressed. Across all tested systems, the margin becomes increasingly negative as renewable penetration and loading increase. This is consistent with a shrinking reachable set along the tested directions.

From a control-theoretic standpoint, stressed operating points modify the reduced dynamics matrix $A_r$ and reduce the effectiveness of the control matrix $B_r$ along critical dynamic directions. In power system terms, increased renewable penetration replaces synchronous inertia with inverter-based resources which weakens electromechanical coupling and reduces controllability of frequency-related modes. These combined effects reduce the directional support function in (4) which leads directly to more negative values of $\delta$. The resulting trend hence provides a quantitative interpretation of reduced controllability in stressed grids, not all terminal frequency or voltage recovery targets remain reachable within the prescribed time horizon and control envelope.

Influence of Control Limits and Terminal Demand

We now investigate how tightening the control envelope or increasing the terminal demand affects the reachability margin. The terminal request is parameterized along a reference direction $x_T(\alpha)=x_T^{\mathrm{fr}}+\alpha s_T^{\mathrm{ref}}$ (where $s_T^{\mathrm{ref}}$ is the displacement direction) so that increasing $\alpha$ corresponds to requesting larger forced deviations from the free terminal response. Figure 4 shows the margin obtained as the terminal request magnitude increases. For all systems, the margin decreases as $\alpha$ grows. This matches the geometric interpretation of Theorem 1 as the requested terminal displacement moves farther from the reachable set, the Euclidean shortfall grows monotonically.

A similar trend is observed when the control limits $\bar{u}$ are reduced. Tightening the control envelope shrinks the reachable set through the scaling matrix $\bar{R}_u$ in (4). This entails smaller control limits correspond to reduced governor authority, limited inverter headroom, or curtailed demand-response participation. Each of these conditions reduces the ability of the system to counteract disturbances within the finite horizon. These observations show that the screening margin captures both actuator limitations and target severity in a unified quantity. This enables operators to assess whether proposed frequency or voltage recovery objectives remain achievable under realistic actuator constraints without performing repeated time-domain simulations.

Case	I [s]	III [s]	CORA [s]	mean $\delta^{\mathrm{I}}$	mean $\delta^{\mathrm{III}}$	mean $\delta_{\mathrm{CORA}}^{\mathrm{out}}$
IEEE 9	0.841	1.427	13.476	-0.0279	-0.02789	-0.02716
IEEE 39	0.658	2.024	–	-1.433	-1.384	–
IEEE 200	1.806	4.961	–	-0.0631	-0.0584	–

Comparison with CORA

We compare the proposed screening method with CORA [29]. Since CORA natively propagates reachable sets for continuous dynamic models, we perform the comparison on the same reduced index-1 model used in (4)–(5), with the same initial state $x_{d,0}$, horizon $T$, and dynamic terminal target $x^{\star}_{d,T}$. In CORA, the input box is represented as a zonotope and the terminal reachable set is computed as an outer approximation $\widehat{\mathcal{R}}_{d,T}$ of the reduced dynamics. When the returned final set is represented as a zonotope, we denote its center by $c_T$ and its generator matrix by $\Gamma_T$, so that

Its support function is $h_{\widehat{\mathcal{R}}_{d,T}}(y)=y^\top c_T+\|\Gamma_T^\top y\|_1$. We then wrap the CORA output into the same target-specific margin geometry used in this paper by defining the outer set margin

If $\delta_{\mathrm{CORA}}^{\mathrm{out}}<0$, then $x^\star_{d,T}\notin\widehat{\mathcal{R}}_{d,T}$ and the reduced dynamic target is certified unreachable; if $\delta_{\mathrm{CORA}}^{\mathrm{out}}\ge0$, the result is inconclusive because $\widehat{\mathcal{R}}_{d,T}$ is an outer approximation. The corresponding algebraic variables may then be recovered afterward through the linearized constraint map $x_a=C_rx_d+D_ru_T$. Results of this comparison are reported in Tab. V.

Since $\widehat{\mathcal{R}}_{d,T}$ is an outer approximation of the exact terminal set of the same reduced model, $\delta_{\mathrm{CORA}}^{\mathrm{out}}\ge\delta^{\mathrm{I}}$ up to solver tolerances, and in Tab. V this ordering holds with gaps of 1.5% to 8.8% of $|\delta^{\mathrm{I}}|$ on IEEE 9. The cost of the explicit set representation grows quickly with dimension. With CORA’s adaptive algorithm at error tolerance $10^{-2}$, the returned IEEE 9 terminal sets carry 6,865 to 34,057 zonotope generators at $n_d=12$, and the same computation failed on IEEE 39 and IEEE 200. This is expected since tightness and tractability depend strongly on the time step, the zonotope order, and the reduction settings. In that sense, while CORA is a powerful general purpose reachability toolbox, it is less suitable for repeated screening of realistic power system models than the present margin computation which avoids explicit reachable set construction altogether.

Paper Summary and Future Work

The paper develops a DAE reachability screening framework for power systems in which the central object is a margin that uses the operating point, the horizon $T$, the control limits $\bar{u}$, and the linearized descriptor model $(E,A,B)$ to assess how far a demanding terminal target is from feasibility. The resulting margin supports rapid comparison of stressed and renewable-rich operating conditions before heavier reachable set construction or nonlinear simulation is carried out. Future extensions include operational path constraints on $[0,T]$ and using the formulation to place PMUs, actuation nodes, or HVDC lines to enlarge the terminal reachable set under heavy renewables or stressed conditions.

Proof of Lemma 1

Every element of $\mathcal{R}^{\mathrm{full}}_T$ is generated by the terminal relation $x(T)=x_T^{\mathrm{fr}}+\Phi_T(u,u_T)$, so it necessarily satisfies the algebraic closure induced by some admissible endpoint value $u_T$; hence $\mathcal{R}^{\mathrm{full}}_T\subseteq\mathcal{M}_T$. By definition of the support function and of $\Phi_T(u,u_T)$,

Let $q(\tau):=B_r^\top e^{A_r^\top(T-\tau)}P^\top y$ and $r:=D_r^\top y_a$. Because the admissible set is a product of independent componentwise box constraints in time and at the endpoint, the maximization separates into a pointwise maximization for $u(\tau)$ and a finite-dimensional maximization for $u_T$. For every $\tau\in[0,T]$,

with a maximizer given componentwise by $u^\star(\tau)=\bar{R}_u\operatorname{sign}(q(\tau))$. Since $q(\tau)$ is continuous on $[0,T]$, $u^\star(\tau)$ is measurable and satisfies $|u^\star(\tau)|\le\bar{u}$ almost everywhere, so $u^\star\in\mathcal{U}_T$ and the supremum over $\mathcal{U}_T$ equals the integral of these pointwise suprema. Likewise,

with maximizer $u_T^\star=R_u\operatorname{sign}(r)$. Substituting these expressions yields

which is exactly (4). Finally, since $\mathcal{R}_T^{\mathrm{full}}=x_T^{\mathrm{fr}}\oplus\mathcal{R}_T^{\mathrm{forc}}$,

which proves (5).

Appendix B: Proof of Theorem 1

Lemma 1 gives the explicit expression for $h_{R_T^{\mathrm{forc}}}(y)$, so (6) is a computable optimization problem. The remaining argument is geometric. Since $\mathcal{R}_T^{\mathrm{forc}}$ is compact and convex, the projection $p^\star:=\Pi_{\mathcal{R}_T^{\mathrm{forc}}}(\xi_T)$ exists. If $x_T^\star\in\mathcal{R}_T^{\mathrm{full}}$, then $\xi_T\in\mathcal{R}_T^{\mathrm{forc}}$. Hence $h_{R_T^{\mathrm{forc}}}(y)\ge y^{\mathsf{T}}\xi_T$ for every $y$, so $\delta(x_T^\star)\ge0$. Choosing $y=0$ in (6) gives $\delta(x_T^\star)\le0$, hence $\delta(x_T^\star)=0$.

If $x_T^\star\notin\mathcal{R}_T^{\mathrm{full}}$, then $\xi_T\notin\mathcal{R}_T^{\mathrm{forc}}$. Define

By the projection theorem,

so $h_{R_T^{\mathrm{forc}}}(y^\star)={y^\star}^{\mathsf{T}}p^\star$. Therefore,

For any $\lVert y\rVert_2\le1$,

Hence

Rewriting this equality yields the second line of (7). Finally, since $\mathcal{R}_T^{\mathrm{full}}=x_T^{\mathrm{fr}}\oplus\mathcal{R}_T^{\mathrm{forc}}$ and $\xi_T=x_T^\star-x_T^{\mathrm{fr}}$, translation invariance gives

which completes the proof.

Appendix C: Proof of Corollary 1

If $\bar{x}_T\in\mathcal{R}_{T,1}^{\mathrm{full}}$, then $\bar{x}_T\in\mathcal{R}_{T,2}^{\mathrm{full}}$ by inclusion, and Theorem 1 gives $\delta_1(\bar{x}_T)=\delta_2(\bar{x}_T)=0$. If $\bar{x}_T\in\mathcal{R}_{T,2}^{\mathrm{full}}\setminus\mathcal{R}_{T,1}^{\mathrm{full}}$, then $\delta_2(\bar{x}_T)=0$ while $\delta_1(\bar{x}_T)<0$, so the claim again holds. If $\bar{x}_T\notin\mathcal{R}_{T,2}^{\mathrm{full}}$, then $\bar{x}_T$ lies outside both sets. Since distance to a set decreases when the set enlarges,

Applying Theorem 1 to both instances yields

The magnitude statement follows immediately.