Unified Diagnostics for Quantifying AC Operating-Point Robustness Under Injection and Topological Uncertainties with Regime Changes

DRAFT FOR JOURNAL SUBMISSION, FEB 2026 1 Uniﬁed Diagnostics for Quantifying A C Operating-Point Rob ustness Under Injection and T opological Uncertainties with Re gime Changes Laurent ¸ iu Lucian Anton* 1 , Marija Ili ´ c 2 Massachusetts Institute of T echnolo gy , 77 Massachusetts A ve, Cambridge , MA 02139, United States Abstract —In the presence of uncertainties in load, generation, and network topology , power system planning must reﬂect oper - ational conditions, while operations require situational awareness over credible uncertainty sets. Diverse methods exist to screen, analyze, embed, and propagate operational uncertainty in power ﬂow and optimal power ﬂow settings. While complementary , these appr oaches offer only partial insight into how physical constraints, controls, and economic interactions shape the ro- bustness of steady-state operating-points under uncertainty . By formulating operating-point rob ustness as a post-solution phys- ical response problem around a solved A C optimal power ﬂow (A C-OPF) equilibrium, this paper provides a uniﬁed framework for assessing rob ustness under injection and topological uncer - tainty without re-optimization. W e construct a primal physical response mapping that accounts for connectivity changes, active power redistribution, generator saturation including P V → P Q transitions, and A C network propagation, and introduce quasi- duals that offer a geometric interpr etation of shadow prices for off-optimal equilibria. Using these r esponse mappings, we develop a suite of deterministic screening procedures that generalize N − k contingency analysis to include cost vulnerability C − k , and local analogs N + δ ( k ) and C + δ ( k ) , deﬁned through sensitivity- normalized margins and risk tolerances. W e extend the frame- work to include probabilistic scr eening for distrib ution-based and moment-based uncertainties, and introduce sequentially-pruned mixture modeling and α -stressed regime constructions to manage combinatorial branching from probabilistic regime transitions. W e present a case study on the Puerto Rican bulk power system that demonstrates how violation and vulnerability diagnostics can be integrated with geospatial data to produce visualizations that enhance situational awareness for operational and planning applications. Index T erms —Operating-Point Robustness, P ost-Solution Screening, Physical Uncertainty Propagation, A C Optimal Power Flow , Power System Reliability I . I N T RO D U C T I O N Modern po wer systems operate under gro wing uncertainty , jointly driv en by imperativ es to decarbonize generation, elec- trify end-use, and maximize network utilization. As variable renew able capacity expands, dispatchable thermal capacity continues to retire. These ﬂeet changes, shaped by both economics and policy incentiv es, vary by technology and region and can reduce ﬂexibility in balancing supply and demand [1, 2, 3]. At the same time, electriﬁcation of heating and transport, combined with rapid data-center and industrial * Corresponding author . 1 ORCID: 0000-0001-7129-2256 2 ORCID: 0000-0002-6835-197X growth, is contributing to sustained demand increases across sectors [4, 5]. T ransmission and distribution networks, already under economic and regulatory pressure, are now operated with reduced margins while increasingly exposed to climate- related risks [6]. Under these stressors, ev en small differences between mod- eled and realized operating conditions, arising from forecast errors, modeling inaccuracies, deviations in generator or DER output, or topology changes, can push a system to ward ther- mal, v oltage, or stability limits. This motiv ates not only re- newed interest in constraint margins and uncertainty propaga- tion, but also a systematic vie w of operating-point robustness to realistic uncertainties in injections and network parameters. A diverse set of fast screening tools exists, each addressing speciﬁc uncertainties but providing only partial visibility into how likely an operating point remains within limits over different uncertainty sets. Traditional security criteria, which perform systematic screening of robustness, are formulated in terms of discrete disturbance sets such as N − 1 and selected N − k contingencies, and most probabilistic or risk- based extensions remain tied to these contingency lists. Gi ven the heterogeneous uncertainties modern power systems must contend with, operators and planners beneﬁt from lightweight, interpretable diagnostics that generalize these notions across varied uncertainty sets, whether discrete, local, or probabilistic. A. Bac kgr ound W e substantiate this narrativ e in detail. W e ﬁrst high- light recent operator reports that document operational trends affecting constraint margins, motiv ating renewed attention to operating-point robustness. W e then provide a high-lev el ov erview of analytical tools developed to capture the effects of injection and topological uncertainty on steady-state equilibria, both to summarize the state of the art and to motiv ate a uniﬁed framew ork that builds on these mature foundations to deliv er new system-wide diagnostics of operating-point robustness ov er meaningful uncertainty sets. 1) Operational T r ends on Constraint Margins: Across the United States, reports from system operators speak to the combined pressures of increasing uncertainties in generation capabilities, load proﬁles, and system states, while operating with reduced margins. NY -ISO warns that as more fossil generators retire, new carbon-free resources are not being added fast enough to DRAFT FOR JOURNAL SUBMISSION, FEB 2026 2 keep pace with expected demand growth, raising concerns about whether the aging ﬂeet will be able to provide essential grid reliability services under shrinking margins [5]. ISO-NE reports “dramatic changes” in its energy mix, shifting from coal and oil to natural gas, with an interconnection queue dominated by wind and batteries. With steady gro wth expected from increased heating and transport electriﬁcation, a shift from summer to winter peaking patterns brings concerns over energy adequacy and deli verability on both the grid and gas pipeline systems [7, 8]. Both PJM and MISO caution that dispatchable resources are increasingly being replaced by weather-dependent gener- ation with materially different capabilities, raising concerns ov er operational challenges associated with limited-duration storage, extreme-weather exposure, and large single-site load additions [9, 10]. CAISO reports gro wing mismatches between supply and demand resulting in wind and solar curtailment, increased ramping requirements, mentions heightened climate- related threats to grid infrastructure, and that battery storage now plays an important role in po wer balancing [11]. ERCO T highlights signiﬁcant projected load growth, in- creased thermal generation retirement, rapid growth in both transmission-connected and distrib uted wind, solar, and en- ergy storage de velopment [12]. SPP similarly comments on shrinking reserve margins that are e xpected to fall belo w its 15% Planning Reserve Margin by 2027, resulting from thermal retirements and expiring contracts, while expressing an increased reliance on demand response [13]. Puerto Rico represents perhaps the most acute example of these operational challenges. Operated by LUMA Energy , the grid remains fragile, contending with compounding climate- related damages, while simultaneously working tow ards rapid grid decarbonization and modernization. Recurrent blackouts, limited ﬂexible reserves, narrow operating margins, and fund- ing constraints continue to impact progress tow ards reliable, resilient operation [14, 15]. The reports abov e collectiv ely point to increasing reliability concerns and gro wing vulnerability to uncertain realizations in generation, demand, or topology under tighter operating margins. These trends motiv ate renewed attention to how close an operating point lies to its physical, regulatory , and stability limits, and how quickly such mar gins are consumed under uncertainty . W e ne xt survey the body of analytical methods that characterize how injection and topological perturbations propagate through steady-state network equations. 2) Appr oaches to Incorporating Uncertainty in Steady-State P ower Flow and Dispatch: Operational uncertainty has long been addressed through methods that either embed uncertainty directly into dispatch formulations or provide post-solution assessments of how perturbations propagate through steady- state netw ork equations. These two classes form the foundation of existing notions of operating-point robustness. A substantial body of work incorporates uncertainty di- rectly into optimal power ﬂow (OPF) formulations. Stochastic, chance-constrained, robust, and distributionally robust OPF treat injections or model parameters as random v ariables or uncertainty sets, enforcing operational limits with a speciﬁed conﬁdence level or across worst-case realizations. Chance- constrained OPF (CC-OPF) requires that constraints hold with probability at least 1 − ϵ [16], while stochastic OPF (S- OPF) often constrains or optimizes a conditional value-at- risk (CV aR). For topological uncertainty , robust OPF (R-OPF) enforces feasibility across speciﬁed contingency sets, and distributionally robust OPF (DR-OPF) optimizes performance under the worst-case distribution within a deﬁned uncertainty set [17, 18]. In practice, market processes such as security- constrained unit commitment and economic dispatch (SC-UC and SC-ED) model contingencies and reserve requirements to handle uncertainty prior to real-time operation. These approaches provide ex ante guarantees by embedding un- certainty , rather than assessing the robustness of a speciﬁc realized AC operating point in post. Despite the v ariety of tools for embedding uncertainty in dispatch, realized system states often occur that were not considered in predictiv e formulations. A complementary suite of methods focus on post-solution screening of operating points, to speciﬁc contingencies. Screening for likely post- contingency constraint violations is often performed using lin- earized distrib ution factors that consider incremental changes in po wer transfers (PTDF), line outages (LODF), line additions (LADF), and outage transfers (O TDF) [19]. These methods provide local sensitivity information, treat uncertainty sources one at a time, and often rely on linear approximations that are accurate within a limited neighborhood of the base case. Probabilistic load ﬂo w (PLF) methods propagate uncertainty in injections or parameters to state variables such as voltages and branch ﬂows. Monte Carlo sampling estimates violation frequencies by repeatedly solving A C power ﬂows under sampled conditions [20]. Analytical PLF approaches propagate input moments or cumulants through linearized mappings to approximate state probability distributions [21, 22, 23]. Point estimation methods (PEMs) approximate these distributions using a small number of deterministic points deriv ed from input moments [24]. PLF has also been formulated to con- sider constrained load ﬂo w and correcti ve control actions to minimize violations [25]. Robustness to steady-state voltage stability limits is com- monly assessed through homotopy-based continuation power ﬂow (CPF), which traces loadability curves and identiﬁes the bifurcation point that deﬁnes maximum loadability [26], beyond which voltage collapse can occur . Steady-state stability margins are thus deﬁned as a distance to this point. In contrast, PLF extensions have used linearized sensiti vities to estimate the probability of collapse and identify correctiv e control actions [27]. These methods each provide important insights into voltage phenomena b ut focus on a speciﬁc degrees of freedom rather than broad robustness assessments. A related class of methods considers risk-based security assessment, where the probability of a contingency and the sev erity of its consequences are combined to compute e xpected risk indices [28]. These methods can be interpreted as proba- bilistic generalizations of N − 1 screening, since the underlying uncertainty sets are enumerated contingency scenarios rather than continuous injection or topology uncertainty models. As such, risk-based approaches extend deterministic formulations, DRAFT FOR JOURNAL SUBMISSION, FEB 2026 3 rather than providing a uniform PFL-style treatment of uncer- tainty around deﬁned operating points. T aken together , these approaches offer valuable yet frag- mented notions of operating-point robustness, each deﬁned with respect to particular physical mappings and uncertainty sets that differ substantially in structure, dimensionality , and operational meaning. None provide a uniﬁed way to assess how a solved A C-OPF equilibrium responds to heterogeneous injection and topology uncertainty . As a result, e xisting notions of rob ustness remain dispersed across tools de veloped for different purposes, making it difﬁcult to interpret what it means for an operating point to be “robust, ” and with respect to which uncertainties. This motiv ates the development of a uniﬁed post-solution framework that ev aluates robustness with respect to clearly deﬁned and meaningful uncertainty sets around an operating point. 3) T owar d a Uniﬁed F rame work for Operating-P oint Ro- bustness: A uniﬁed view of operating-point robustness can be deﬁned relativ e to a known network state and a speciﬁed uncertainty set, and should consider both physical feasibility and optimality . An A C optimal power ﬂow (A C-OPF) solution provides such a reference, as it speciﬁes a feasible A C oper- ating point that satisﬁes nonlinear power balance equations and operational limits while optimizing a chosen objective, often cost. The associated primal and dual information jointly characterize physical mar gins and economic tightness, mak- ing AC-OPF a natural generalized reference formulation for propagating injection and topological perturbations in a uniﬁed robustness frame work. While an AC-OPF solution fully speciﬁes base operating points, it does not by itself deﬁne ho w the system responds to physical perturbations. Changes in steady-state equilibria should reﬂect layered activ e and reactive po wer control mech- anisms, on the timescale of interest. AC power ﬂow (A C- PF) and many PLF formulations assume that acti ve po wer imbalances are absorbed by a single slack bus [29], with distributed slack treatments considered in speciﬁc settings [30]. These are typically interpreted as feedback control to fre- quency deviations, as in automatic generation control (AGC), in contrast to feed-forward prev entiv e or proportional reserve allocations that balance predicted imbalances. Reactiv e power control introduces further complexities through control regime transitions for voltage-controlling generators when reactive limits are reached. These saturation effects are critical for realistic steady-state voltage behavior but poorly represented in PLF studies [31]. A uniﬁed robustness frame work must therefore encode these acti ve and reacti ve control layers, rather than relying solely on a single slack-b us or static bus types. Perturbations applied to an A C-OPF equilibrium with spec- iﬁed control logic generally driv e the system toward a ne w steady-state that is no longer optimal with respect to the orig- inal objective and may violate operational limits. In this off- optimal re gime, the Karush–K uhn–T ucker (KKT) conditions no longer hold and an exact dual formulation no longer exists. Nev ertheless, local projections of constraint sensitivities can be interpreted as quasi-dual signals that approximate how tight constraints are economically in off-optimal conditions, enabling robustness assessment from both primal and quasi- dual perspectives. Once a physical response mapping is speciﬁed, operating- point robustness can be assessed with respect to arbitrary perturbations in injections and topology . The various notions of robustness re viewed in Section I-A2 consider different un- certainty sets and a base operating point. A uniﬁed frame work should therefore admit these div erse uncertainty sets within a common formulation, so that con ventional criteria like N − 1 security become special cases. In many practical settings, detailed probability distributions for uncertain injections or parameters are unav ailable, and only limited statistical information can be estimated reliably . In such cases, distribution-free inequalities such as Cantelli bounds provide conserv ativ e guarantees on constraint violation probabilities based solely on means and variances. These bounds have been used to derive tractable reformulations in CC-OPF problems [32], b ut, to the best of our knowledge, have not been adopted as post-solution diagnostics for A C-OPF operating points, using PLF techniques. A uniﬁed rob ustness framew ork should accommodate both distributional models and moment-based bounds, enabling conservati ve robustness assessments when uncertainty information is incomplete. A uniﬁed formulation should thus generalize traditional se- curity notions by quantifying exposure to constraint violations, erosion of margins, and potential constraint release under deterministic or probabilistic perturbations. In practice, these diagnostic layers can be aggregated into system-lev el visu- alizations that summarize robustness across entire uncertainty sets, extending traditional single-contingenc y views into multi- scenario assessments. Such visualizations support situational awareness by revealing spatial patterns in state margins and violations, relev ant to both operational and planning contexts. B. Contrib utions This paper dev elops a uniﬁed post-solution frame work for assessing the robustness of A C-OPF operating points under injection and topology uncertainties. In doing so, we: • Formalize operating-point robustness of a solved A C- OPF equilibrium by unifying deterministic contingenc y analysis, local perturbation screening, and probabilistic uncertainty modeling within a common framework. • De velop a control- and constraint-aware steady-state re- sponse mapping that incorporates regime changes result- ing from islanding, generator control layers, and satu- ration effects when propagating injection and topology perturbations through the network. • Propose a quasi-dual construction that extends the inter- pretation of shadow prices to off-optimal steady-states reached under small perturbations. • Deﬁne deterministic and probabilistic robustness metrics for post-solution screening, including margin-based indi- cators, violation and release probabilities, and moment- based bounds for when uncertainty information is limited. • Introduce probabilistic methods for handling combina- torial, uncertainty-induced regime changes, including a sequentially pruned mixture model for distributional un- DRAFT FOR JOURNAL SUBMISSION, FEB 2026 4 certainty and an α -stressed construction for moment- based uncertainty . • De velop system-wide diagnostic layers that unify primal, quasi-dual, and probabilistic information, enabling inter- pretable visual assessments of robustness over speciﬁed uncertainty sets. • Demonstrate the frame work on the Puerto Rican bulk power system, illustrating its ability to visualize system- wide vulnerabilities. The rest of the paper is organized as follows. Section II dev elops the framework for steady-state A C physical response. Sections III and IV introduce the proposed deterministic and probabilistic robustness screening tools. Section V demon- strates these tools on the Puerto Rican bulk power system, with visualizations. Section VI summarizes the ﬁndings and outlines directions for future work. I I . S T E A D Y - S TA T E P H Y S I C A L R E S P O N S E T O U N C E RT A I N T Y A R O U N D AC - O P F E Q U I L I B R I A Assessing robustness requires understanding how an AC system in equilibrium responds to small changes in injections or network parameters. An AC-OPF solution provides an operating point that is optimal over the feasible region deﬁned by the constraint set. Primal variables encompass bus voltages, branch ﬂows, and generator dispatch, equality constraints determine the feasible manifold and inequality constraints determine the feasible region(s) on that manifold. T o quantify how constraint margins are af fected, we model the steady-state response of the system to variations in activ e and reactive injections and in network admittances. From an A C-OPF equilibrium, we construct a primal perturbation map- ping that considers active power redistrib ution via participation factors and generator limits, reactiv e power saturation via P V → P Q regime changes, and propagation through A C network equations around the base point. W e then introduce a quasi-dual construction that interprets how the cost gradient aligns with the local constraint manifold when the system settles at a physically feasible but off-optimal equilibrium. A. A C-OPF Primal F ormulation and Notation W ithout loss of generality , we consider an AC-OPF for- mulation for activ e po wer economic dispatch. Let N , G , D , and B be the sets of buses, generators, loads, and branches, respectiv ely . Let Y bus be the admittance matrix, deﬁned by considering all network components including series and shunt compensation. Y bus has complex elements Y ij = | Y ij | ∠ δ ij = G ij + j B ij between buses i and j , that can be represented either using admittance magnitude | Y ij | and angle δ ij , or conductance G ij and susceptance B ij . Deﬁne apparent po wer as S = P + j Q , with acti ve po wer P , and reacti ve po wer Q , and bus voltage V = | V | ∠ θ , with v oltage magnitude | V | and angle θ . At each bus i , deﬁne P i = P g ∈G i P g − P d ∈D i P d and Q i = P g ∈G i Q g − P d ∈D i Q d as the net activ e and reactiv e power injection giv en generators G i ⊆ G and loads D i ⊆ D connected at bus i . Let S ij = V i ( V i − V j ) ∗ Y ∗ ij and θ ij = θ i − θ j be the apparent power ﬂow and angle difference across a branch connecting buses i and j . Branches are taken to be either lines or transformers. Lastly , assume at most quadratic functions f ( P g ) = a g P 2 g + b g P g + c g in the objectiv e function. W e then have: min { P g , Q g , | V i | , θ i } X g ∈G f ( P g ) (1) s . t . P i = | V i | X j ∈N | Y ij || V j | cos( θ ij − δ ij ) ∀ i ∈ N (2) Q i = | V i | X j ∈N | Y ij || V j | sin( θ ij − δ ij ) ∀ i ∈ N (3) P min g ≤ P g ≤ P max g , ∀ g ∈ G (4) Q min g ≤ Q g ≤ Q max g , ∀ g ∈ G (5) V min i ≤ | V i | ≤ V max i , ∀ i ∈ N (6) | S ij | ≤ S max ij , ∀ ( i, j ) ∈ B (7) θ min ij ≤ θ ij ≤ θ max ij , ∀ ( i, j ) ∈ B (8) θ i 0 = 0 , for the reference bus. (9) W e represent (2) and (3) with | Y ij | and ∠ δ ij instead of G ij and B ij , for compactness. Fixed voltage setpoints can be modeled with tight box limits in (6) or replaced with equality constraints V i = V set i . Branch limits may constrain P ij rather than S ij . W e lastly note that many other objective functions are possible for AC-OPF , including those that assign costs for reactiv e power compensation or emissions [33]. B. Modeling Steady-State AC Physical Response How a system responds to changes in inputs depends on the timescale of interest. Here we assume a continuous and stable transition between steady-state A C equilibria, without triggering protection or discrete switching ev ents. W e model the physical response of a system in two stages. First we model distributed activ e power balancing, to generate a combined disturbance vector with a minimized imbalance. This combined response is then used to parameterize a stan- dard A C-PF problem. T ogether , these stages characterize the physical response of the network using both generator controls and power balance equations (2)–(3). 1) A C P ower Flow Bus Conventions: In traditional AC-PF , two of four bus variables ( P , Q , | V | , and θ ) are assigned by a giv en dispatch, thereby labeling the buses by what is ﬁx ed. The standard bus deﬁnitions are P V -b uses, where generator active power output is ﬁxed and the bus voltage is locally controlled via generator reactiv e power output; P Q -buses, where the activ e and reactiv e power injections are assigned ﬁxed, often at load buses; and a slack bus, which both serves as a voltage reference bus, and also the reserve generator bus that absorbs all mismatches after A C-PF is solved for all other buses [34]. 2) F ixed and Non-Contr ollable Elements: W e treat all loads as ﬁxed P Q elements. Non-dispatchable generation (e.g., renew able resources without acti ve control) or generation not participating in acti ve po wer balancing are represented as negati ve loads. Thus, buses containing only ﬁxed injections are modeled as P Q -buses. DRAFT FOR JOURNAL SUBMISSION, FEB 2026 5 3) Generator Active P ower Response: When deﬁning a PF problem from an OPF state, active power redistribution must be deﬁned. W e distinguish four power balancing con ventions: 1) Single-slac k: a simpliﬁed representation, where one large generator absorbs all power mismatches; 2) Primary contr ol: distributed compensation ov er partici- pating generators via frequency-droop response; 3) Secondary control: distributed compensation via auto- matic generator control (AGC) response, acting on slower timescales and restoring nominal frequency; 4) T ertiary contr ol: redispatch via a new OPF solve, repre- senting an economic re-optimization rather than a physi- cal response, and not considered here. The single-slack formulation is the standard AC-PF rep- resentation, with a designated slack generator at the refer- ence bus absorbing activ e power imbalances. Both primary and secondary frequency control admit a distributed slack interpretation using participation factors. Primary control acts ov er seconds via gov ernor po wer–frequency droop to restore power balance, while secondary control acts ov er minutes via A GC to restore nominal frequency and eliminate area control error (A CE) [35, 36]. In both cases, the resulting steady-state redistribution corresponds to automated feedback in response to frequency de viations, rather than preplanned, feed-forward correctiv e redispatch conditioned on speciﬁc disturbance sce- narios. Apart from how participation factors are determined, these control layers are algebraically equiv alent at equilibrium and deﬁne the activ e power response model adopted here. W e model all controllable generators as responsive ele- ments with hard capacity limits on acti ve and reactiv e power , represented in (4) and (5). Let G r ⊆ G be the set of dis- patchable generators participating in automated activ e power balancing. W e consider the generalized case, where more than one generator can participate in active power balance via participation factors, α g , where P g ∈G r α g = 1 . Generators whose injections are directly perturbed are excluded from the balancing set G r . Assume a net perturbation in active power injections, dP inj = P i ∈N p dP i where N p is the set of buses with perturbed injections. The response from each participating generator is then dP g = − α g dP inj . (10) This deﬁnes a local feasibility range for perturbations, or necessitates at least one generator without hard constraint limits. Further , the response from each participating generator should, in principle, respond to the total imbalance including incremental losses, dP losses . Ho wever , since the incremental losses are not kno wn a priori, these must be implicitly handled. Here, we assume that the incremental losses are picked up by a designated generator g 0 after dP inj is distributed among participating generators. If g 0 is at the reference bus, this coincides with the numerical slack bus in A C-PF formulations. a) Generator Limit Considerations: W e consider active power operational limits using an iterati ve residual-imbalance method. Let dP ( ℓ ) res denote the residual net-injection imbalance at iteration ℓ , initialized as dP (0) res = dP inj . For each partici- pating generator g ∈ G ( ℓ ) r , the predicted response is P pred , ( ℓ ) g = P g + dP ( ℓ ) g = P g − α ( ℓ ) g dP ( ℓ ) res . (11) If one or more generators are predicted to violate limits, they are removed from the balancing set and clamped at the corresponding limit. The realized output is P ( ℓ +1) g =      P pred , ( ℓ ) g , if P min g < P pred , ( ℓ ) g < P max g , P min g , if P pred , ( ℓ ) g ≤ P min g , P max g , if P pred , ( ℓ ) g ≥ P max g . (12) Let G ( ℓ ) s denote the set of generators clamped at their limits up to iteration ℓ , and deﬁne the newly clamped set ∆ G ( ℓ +1) s = G ( ℓ +1) s \ G ( ℓ ) s . For each s ∈ ∆ G ( ℓ +1) s , assign α ( ℓ +1) s = 0 and renormalize the remaining participation factors such that P g ∈G ( ℓ +1) r α ( ℓ +1) g = 1 . The residual imbalance is then dP ( ℓ +1) res = dP ( ℓ ) res + X s ∈ ∆ G ( ℓ +1) s  P ( ℓ +1) s − P s  . (13) This is repeated until the residual power is balanced ( dP ( ℓ +1) res = 0 ) or no more participating generators remain ( G ( ℓ +1) r = ∅ ). 4) Generator Reactive P ower Response: Activ e power bal- ancing through participation factors uses system-wide fre- quency deviations as a shared feedback signal. Reactiv e power response lacks an analogous uni versal imbalance signal. In- stead, voltage regulation is implemented through local excita- tion control, while wider-area coordination such as pilot-bus Automatic V oltage Control (A VC) is imposed by supervisory control based on zonal measurements and setpoint adjustments [37]. W e therefore adopt the standard AC-PF conv ention of representing all voltage-controlled buses as P V b uses, for which the bus v oltage magnitude is held ﬁxed while the generator reactive output satisﬁes Q min g ≤ Q g ≤ Q max g [34]. If a perturbation causes a generator’ s reactiv e output to reach one of its limits, we assume saturation and ﬁx Q g = Q min g or Q max g , thereby losing v oltage regulation at the corresponding bus and conv erting it from P V to P Q operation [38, 39]. Once saturated, the generator is assumed to supply a ﬁxed reactiv e injection, and the bus voltage is no longer directly controlled. W e refer to this loss of voltage regulation as a Q -saturation regime change, discussed further in Section II-B7. 5) A C Network Response: W e no w describe how the system state responds to perturbations in power injections and net- work topology . Local perturbations are modeled as continuous variations, while larger disturbances are treated as discrete ev ents. T opological changes are represented through variations in network admittance elements Y ij . Although physically implemented via discrete actions such as tap adjustments or switching of transformers and ﬂexible AC transmission system (F ACTS) devices, they are commonly approximated as continuous controls for sensitivity analysis. Generator and branch contingencies are permitted as discrete e vents. Let N p ⊆ N denote the set of b uses where acti ve power injections are perturbed, and N r ⊆ N represent the buses DRAFT FOR JOURNAL SUBMISSION, FEB 2026 6 where generators in G r are connected. After redistribution in response to dP inj = P i ∈N p dP (inj) i , the resulting acti ve power perturbation at bus i is dP i =      dP (inj) i , ∀ i ∈ N p , dP ( l ) i , ∀ i ∈ N r \ ( N p ∩ N r ) , 0 , ∀ i / ∈ ( N p ∪ N r ) . (14) The updated acti ve po wer perturbation vector d P = [ dP i ] ⊤ is used to deﬁne the new dispatch condition for the perturbed state, and is propagated through the A C network equations to obtain perturbed bus state for P , Q , | V | , and θ . Let ρ = ( u , Y bus ) collect the kno wn parameters. The vector u = [ P Q ] ⊤ contains the speciﬁed active and reactiv e po wer injections, while y = [ θ | V | ] ⊤ collects the corresponding voltage variables. In particular, P is speciﬁed and θ is solved for at all non-slack buses N \ { 0 } , while Q is speciﬁed and | V | is solved for at P Q buses N P Q . Deﬁne bus power mismatches as F ( y ; ρ ) =  P ( y ; Y bus ) − P Q ( y ; Y bus ) − Q  = G ( y ; Y bus ) − u . (15) For a perturbation d ρ = ( d u , d Y bus ) , the total dif ferential of each component F k is dF k = X m ∂ F k ∂ y m dy m + X m ∂ F k ∂ u m du m + N X i =1 N X j =1 ∂ F k ∂ Y ij d Y ij . (16) Stacking all k giv es d F = [ dF k ] ⊤ . Compactly , d F = F y d y + F u d u + F Y : d Y bus , (17) with F y = ∂ F ∂ y , F u = ∂ F ∂ u , F Y = ∂ F ∂ Y bus . (18) W e take the tensor contraction F Y : d Y bus ov er a minimal basis spanning the admissible perturbation subspace of Y bus , using only indices corresponding to perturbed branch admit- tances. Since F ( y ; ρ ) = G ( y ; Y bus ) − u , we then have F y = ∂ G ∂ y = J , F u = − I , F Y = ∂ G ∂ Y bus , (19) with J being the Jacobian. Thus, in general form, d F = J ( y ; ρ ) d y − d u + F Y ( y ; ρ ) : d Y bus . (20) At any equilibrium, there are no po wer imbalances. Then F ( y ; ρ ) = 0 , deﬁnes a smooth manifold of feasible operating points in the joint space of voltage v ariables and parameters. A perturbation ( d u , d Y bus ) induces displacements ( d y , d ρ ) that remain on this manifold if and only if the total differential d F ( y ; ρ ) v anishes. a) F irst-Or der Appr oximation: Consider a perturbation ∆ ρ that induces a change ∆ y such that the system settles to a ne w equilibrium. Since both the base and perturbed operating points satisfy the power -ﬂow equations, we hav e F ( y + ∆ y ; ρ + ∆ ρ ) − F ( y ; ρ ) = 0 . (21) T aylor expanding about ( y 0 , ρ 0 ) , we get J 0 ∆ y − ∆ u + F Y 0 : ∆ Y bus + O ( ∥ ∆ y , ∆ ρ ∥ 2 ) = 0 , (22) where J 0 = J ( y 0 ; ρ 0 ) , F Y , 0 = F Y ( y 0 ; ρ 0 ) , and O ( ∥ ∆ y , ∆ ρ ∥ 2 ) represents all higher-order terms. T o the ﬁrst- order , ∆ y ≈ J − 1 0 (∆ u − F Y , 0 : ∆ Y bus ) . (23) This ﬁrst-order mapping describes the local sensitivity of bus voltages and angles to small injection and topology pertur- bations, assuming the Jacobian remains nonsingular and b us deﬁnitions remain unchanged. 6) Monitor ed Quantities and F irst-Or der Response: Con- sider a monitored quantity , x ( y ; Y bus ) , such as a bus state variable, generator output or branch ﬂow . Under small pertur- bations in ∆ u and ∆ Y bus , the ﬁrst-order change in x is ∆ x = ∂ x ∂ y     0 ∆ y + ∂ x ∂ Y bus     0 : ∆ Y bus . (24) Substituting (23) giv es a practical decomposition into an injection term and a topology term: ∆ x =  ∂ x ∂ y     0 J − 1 0  ∆ u +  ∂ x ∂ Y     0 − ∂ x ∂ y     0 J − 1 0 F Y , 0  : ∆ Y bus = L x ∆ u + T x : ∆ Y bus . (25) All deri vati ves and sensiti vities are e valuated at the ﬁxed operating point ( y 0 , ρ 0 ) and deﬁne numerical linear operators once the base point is speciﬁed. While useful to characterize directional sensitivities, perturbed monitored quantities x ′ can also be ev aluated from the perturbed steady-state solution, y ′ . 7) Deterministic Regime Changes: The network response abov e assumes generators control modes are ﬁxed and the net- work topology remains well-deﬁned. In practice, perturbations in injections or topology may trigger discrete regime changes associated with generator saturation or network islanding. These ev ents introduce piecewise-smoothness to the feasible manifold, and can be treated sequentially . W e consider regime changes in the following order: 1) network admittance changes related to islanding, 2) generator activ e po wer saturation, then ﬁnally 3) generator reacti ve power saturation resulting in P V → P Q changes in bus deﬁnitions. P -saturation changes are detailed in Section II-B3. W e de- scribe the remaining two classes of regime changes below . a) Y bus Re gime Chang es: Discrete perturbations in branch admittances may induce changes in network connec- tivity when admittances are driv en to zero. This effecti vely remov es the connection and can electrically island buses. T o av oid ill-conditioned or singular sensitivity matrices, the activ e system connected to the original slack bus must be reduced by remo ving these islanded components, which can then be handled separately based on modeling assumptions (e.g., ride-through control protocols or assumed blackout). Since generator control actions and P V → P Q transitions depend on connectivity , Y bus regime changes must be resolved prior to any generator saturation sequencing. DRAFT FOR JOURNAL SUBMISSION, FEB 2026 7 b) Q -Saturation Regime Changes: Following resolution of active po wer and network-admittance regime changes, reac- tiv e power limits are enforced through sequential P V → P Q transitions. Beginning from the solved base state, we compute the incremental change in voltages ∆ y in the current regime. From this, the ﬁrst-order change in reactive power at each PV bus is obtained as ∆ Q i = ∂ Q i ∂ y ∆ y + ∂ Q i ∂ Y bus : ∆ Y bus . (26) The predicted reactiv e power Q pred i = Q i + ∆ Q i is then distributed to each generator g at b us i and compared to its limits [ Q min g i , Q max g i ] . Generators for which Q pred g i lies outside this range are ﬂagged as in violation. If no violations are detected, the perturbation does not induce a regime change and the iteration ends. Otherwise, each violating generator is assigned a normalized ev ent ratio t lim g i = Q lim g i − Q g i ∆ Q g i , (27) where t lim g i measures how far along the disturbance direction each output would change before its Q limit is reached. The generator with the smallest positiv e t lim g i is considered the ﬁrst to reach a limit under the applied perturbation, and is clamped to the associated limit. If no more voltage-controlling generators exist at i , the bus is reclassiﬁed from P V to P Q and system equations are updated accordingly . The Jacobian and associated partial deriv ativ es are reconstructed for the ne w regime, and the perturbation is remapped to obtain an updated state response. This process is repeated until no further Q -limit violations are detected. 8) Discussion on Proximity to Stability Limits: V oltage stability phenomena are dynamical with nonlinear evolution, impacted by control interactions. In steady-state, voltage insta- bility is associated with loss of regularity of the power-ﬂo w equations, which manifests as a singular Jacobian. Stability margins can be computed using CPF along prescribed load- ability paths [26]. Implicit proximity measures can be deriv ed from Jacobian-based indicators such as the smallest singular values, which provide a direction-independent measure of conditioning [40]. Moreover , Jacobian singular value decom- position further insight by identifying combinations of power - balance equation imbalances for which small coordinated violations require large changes in bus voltage magnitudes and angles to restore equilibrium, often referred to as weak modes [40]. These modes are intrinsic properties of the operating point and network structure and describe patterns of imbalance to which the system is vulnerable, independent of any speciﬁc disturbance direction being applied. Assessing proximity to stability limits along a speciﬁed ∆ ρ is conceptually closest to CPF methods. In principle, an analogous continuation procedure could be constructed by iterativ ely scaling ∆ ρ , recomputing the power ﬂo w solution, and tracing the ev olution of monitored quantities such as voltage magnitudes and reactiv e injections. As an initial fast screening, one could use the local linearized change ∆ x in response to ∆ ρ as a crude indicator for loss of regularity . If ∆ x > ε reg or ∥ ∆ x ∥ > ε reg for some tolerance ε reg , the equilibrium response may be ill-conditioned along ∆ ρ , and can trigger more detailed assessments. C. Of f-Optimal Quasi-Dual F ormulation At an A C-OPF optimum, the KKT conditions imply that the cost gradient lies in the span of the gradients of the activ e constraints. The associated Lagrange multipliers ν = ( λ , µ , σ ± , τ ± , η ± , π ± , ζ ) represent shadow prices associated with the power balance equations (2)–(3) and activ e opera- tional limits (4)–(9). These multipliers are well deﬁned only at optimality , where the stationarity condition holds exactly . Follo wing a disturbance, the system settles to an equilibrium that is physically-feasible with respect to equality constraints (2) and (3), b ut possibly not operationally-feasible with respect to inequality constraints (4)–(9). Since no re-optimization occurred, the KKT conditions no longer hold and exact mul- tipliers do not exist. Still, it is useful to retain multiplier-lik e quantities that describe how the cost gradient interacts with the local constraint manifold. This motiv ates an off-optimal quasi-dual construction, which extends the notion of shadow prices by projecting the cost gradient onto the span of active constraint gradients [41]. Let the primal decision vector be z = ( P G , Q G , V N , θ N ) , (28) with associated multipliers ν ≥ 0 . The pre-disturbance opti- mal equilibrium ( z ⋆ , ν ⋆ ) satisﬁes all KKT conditions. A small perturbation ( d u , d Y bus ) modiﬁes injections and network pa- rameters to ρ ′ = ρ ⋆ + d ρ = ( u ⋆ + d u , Y ⋆ bus + d Y bus ) , (29) resulting in a new physically-feasible off-optimal state z ′ = z ⋆ + d z . (30) 1) Constraint Manifold and Cost Gr adient: Around any state z ′ , the activ e constraints deﬁne a smooth manifold in the space of decision v ariables. The directions normal to that surface are gi ven by the constraint gradients, and the directions along the surface are those that keep constraints satisﬁed to ﬁrst order . Let the corresponding constraint-gradient matrix be C ( z ′ ) =  ∇ z h ( z ′ ) ∇ z g A ( z ′ )  , (31) where h ( z ′ ) = 0 collects the equality constraints, and g A ( z ′ ) ≤ 0 collects the activ e and possibly violated inequality constraints at z ′ . The cost gradient at z ′ is ∇ z f ( z ′ ) =     ∇ P G f ( P ′ G ) ∇ Q G f ( Q ′ G ) ∇ | V | N f ( | V | ′ N ) ∇ θ N f ( θ ′ N )     =     ∇ P G f ( P ′ G ) 0 Q G 0 | V | N 0 θ N     , (32) since only P G enters the objectiv e (1). At z ⋆ the KKT stationarity condition can be written as ∇ z f ( z ⋆ ) + C ( z ⋆ ) ⊤ ν ⋆ = 0 , (33) DRAFT FOR JOURNAL SUBMISSION, FEB 2026 8 which means ∇ z f ( z ⋆ ) is a linear combination of constraint normals. At an A C feasible but suboptimal z ′ ∇ z f ( z ′ ) + C ( z ′ ) ⊤ ν  = 0 for any ν , (34) and the nonzero residual is the component of the cost gradient that still points along the tangent directions, indicating possible cost reduction without violating constraints. 2) Quasi-Duals fr om Least-Squar es Pr ojection: W e seek multiplier-lik e coefﬁcients that make (33) hold as closely as possible at z ′ . T o do this, we perform a linear least-squares ﬁt of the cost gradient by a combination of constraint normals. Deﬁne the quasi-dual vector at z ′ as ˜ ν ( z ′ ) = arg min ν A   ∇ z f ( z ′ ) + C ( z ′ ) ⊤ ν   2 2 . (35) The quantity inside the norm is the stationarity residual , which vanishes only at a true optimum. This formulation does not impose sign restrictions on the components of ˜ ν . Instead, the quasi-duals are deﬁned purely by the geometry of the cost gradient relativ e to the local constraint manifold. Minimizing (35) leads to the normal equations C ( z ′ ) C ( z ′ ) ⊤ ˜ ν ( z ′ ) = − C ( z ′ ) ∇ z f ( z ′ ) , (36) and if C ( z ′ ) C ( z ′ ) ⊤ is inv ertible, the closed-form solution is ˜ ν ( z ′ ) = −  C ( z ′ ) C ( z ′ ) ⊤  − 1 C ( z ′ ) ∇ z f ( z ′ ) . (37) Each component ˜ ν k of ˜ ν ( z ′ ) measures the signed projection of the unconstrained cost gradient onto the corresponding constraint normal. Its magnitude indicates ho w strongly that constraint shapes nearby cost-reducing directions. If ˜ ν k > 0 , the cost-reducing direction points into the constraint, so the constraint resists further cost reduction and acts as an econom- ically acti ve limit. If ˜ ν k < 0 , relaxing the constraint would reduce cost, indicating that it is not economically limiting at z ′ , e ven if violated. Because quasi-duals depend on the local activ e set, they may vary discontinuously and are best interpreted as local geometric indicators. I I I . D E T E R M I N I S T I C C O N S T R A I N T R O B U S T N E S S W e now characterize robustness with respect to opera- tional constraints for deterministic perturbations. Let moni- tored quantity x hav e bounds L min ≤ x ≤ L max . W e deﬁne the associated constraint margin as M lim x = s ( L lim − x ) , with s = ( +1 , if L lim = L max , − 1 , if L lim = L min , (38) so that M lim x > 0 indicates a respected limit for either bound. Except where needed, we drop “lim” superscripts for clarity , with L and M x understood to either bound unless speciﬁed. A. Discr ete Physical Robustness ( N − k ) For a k th -order discrete disturbance ∆ ρ , the system is ev aluated at the perturb state z ′ and the feasibility of all constraints is checked. If M x < 0 for any x , then a constraint is violated, and the perturbation is labeled as critical. A critical subset of violating quantities { x viol } can be created for each critical perturbation for further consideration. This corresponds directly to the classical N − k screening procedure used in operational planning, where a system is “ N − k secure” if, for all contingencies in the screening set, a feasible solution exists that respects all limits [42]. Whether as a post-solution screen or as part of an R-OPF formulation, it is also possible to allow for corrective actions between pre and post disturbance states. Here, we consider automated generator actions deﬁned in Sections II-B3 and II-B4. B. Local Physical Robustness ( N + δ k ) Let d ρ be any k th -order perturbation direction, meaning any simultaneous variation of k components of d ρ . The directional sensitivity of x with respect to d ρ is dx d ρ = ∂ x ∂ y d y d ρ + ∂ x ∂ ρ . (39) The directional change in the margin M x is dM x d ρ = d  s ( L − x )  d ρ = − s dx d ρ . (40) W e deﬁne the sensitivity-normalized margin as R x = M x dM x /d ρ = − sM x dx/d ρ , M x > 0 . (41) The index R x ev aluated along d ρ is ﬁnite when dM x /d ρ  = 0 and is inﬁnite when the constraint is insensitive in ﬁrst order . Its magnitude | R x | = M x / | dM x /d ρ | represents the normalized distance to the limit in units of perturbation mag- nitude, while its sign indicates whether the mar gin shrinks or expands along d ρ . Negati ve values identify shrinking margins, while positiv e values indicate expanding margins that are directionally robust ev en when M x is small. W e can use R x to create a local analog of N − k security . If − ε R < R lim x < 0 for a risk tolerance ε R > 0 , we classify the pair ( x, lim) as locally vulnerable. W e can thus say that an operating point is locally N + δ k robust if, for e very considered k th -order d ρ , all R x satisfy R x / ∈ [ − ε R , 0) . a) Choice of Disturbance Basis: For contingencies, the direction of d ρ is ﬁxed by the loss of individual components. In contrast, the local formulation admits inﬁnitely many per- turbation directions within a neighborhood of ρ , requiring the selection of a representati ve direction set for rob ustness screening. Lacking additional information, one may consider a cartesian-aligned basis with equal signed ± variations to individual P , Q , G , or B components, or a parameter-aligned basis that scales perturbations along or orthogonal to P + j Q or G + j B . In principle, any basis may be chosen. For example, local directions may be constructed to represent dominant modes of measurement or prediction error . DRAFT FOR JOURNAL SUBMISSION, FEB 2026 9 R x also depends on the normalization of d ρ . Directional deriv ati ves are deﬁned with respect to unit-length directions, and are dimensionless. For screening purposes, howe ver , it may be useful to retain a physical scaling of ρ , such as a unit step as a 10 MW injection change or a 1 % variation in line admittance. In this case, quantities such as dM x /d ρ represent changes in margin per physically meaningful disturbance. Unit-normalized and physically-scaled ρ can thus be used to address different operational questions. C. Discr ete Economic Robustness ( C − k ) For a discrete ∆ ρ , the system is ev aluated at the perturbed operating point z ′ and the associated quasi-duals ˜ ν ( z ′ ) are computed. Economic rob ustness is then assessed by comparing ˜ ν ( z ′ ) against a prescribed cost-tolerance threshold ε ν > 0 , either element-wise or with a chosen norm. If any element or the chosen-norm exceeds ε ν , the perturbation is labeled as economically vulnerable. This construction mirrors classical N − k screening on the dual side. For a set of discrete k -order disturbances, the system is said to be C − k economically robust if for all admissible k th - order perturbations, ˜ ν ( z ′ ) or its elements are less than ε ν . C − k robustness thus provides an economic analog to physical N − k security , characterizing whether post-disturbance operating points remain economically well-conditioned. D. Local Economic Robustness ( C + δ k ) Section III-B introduced the sensiti vity-normalized margin R x which measures, to ﬁrst order , how far a constraint remains feasible under an inﬁnitesimal disturbance direction d ρ . A natural question is whether an analogous quantity can be deﬁned on the dual side. In principle, one may consider a dual margin M ν x = ν max x − ν x , and the associated sensiti vity- normalized margin R ν x = M ν x dν x /d ρ , M ν x ≥ 0 , (42) which would represent the ﬁrst-order change in ν x associated with x tow ard a shadow-price limit ν max x . Howe ver , two issues arise. First, many markets do not consider explicit price caps, making ν max x undeﬁned. Second, computing dν x /d ρ requires differentiating the quasi-dual system (36), with respect to the disturbance direction d ρ . This introduces deri vati ves of both the active-constraint matrix C ( z ′ ) and the cost gradient ∇ z f ( z ′ ) , producing a second-order system that is compara- tiv ely more expensi ve to e valuate. A simple alternati ve is to characterize ho w strongly each activ e constraint contributes to the change in marginal cost along d ρ . Using the standard KKT decomposition ∇ z f ( z ) ≈ C ⊤ ν . (43) The directional cost sensitivity then satisﬁes d f d ρ ≈ ν ⊤ C d z d ρ = X a ∈A ν a dg a d ρ , (44) where dg a /d ρ is the directional change of the constraint function induced by the state response. Equiv alently , this directional change satisﬁes dg x /d ρ = − dM x /d ρ . W e deﬁne the shadow-price volatility index (SPVI) for a constrained monitored quantity as ς x = ν x dg x d ρ = − ν x dM x d ρ , (45) which measures how strongly the binding constraint on x contributes to the local cost sensiti vity along d ρ . T o turn this into a screening tool, we introduce a volatility tolerance ε SPVI > 0 and classify each activ e constraint k as locally cost-robust along ρ if | ς x ( ρ ) | < ε SPVI , and cost-vulnerable otherwise. Similar to N + δ k physical robustness screening, a system is said to be C + δ k economically robust with respect to a family of k th -order disturbance directions if every active constraint is cost-robust. I V . P RO B A B I L I S T I C C O N S T R A I N T R O B U S T N E S S The screening tools in Section III assess how operating margins change in response to a speciﬁed perturbation or direction. In practice, forecast errors and operational variability are often random. It is thus valuable to characterize not only the distance to constraints but also the probabilities of constraint violation. Here we build on PLF foundations to extend the deterministic framew ork to random perturbations, using the ﬁrst-order mappings deriv ed in Section II. W e ﬁrst consider perturbations within a ﬁxed network and control regime, such that the linearized response is well deﬁned. W e then address regime changes in Section IV -F. A. Linear Primal Uncertainty Mapping In the deterministic setting, perturbations in injections ∆ u and network parameters ∆ Y bus appear as separate objects. In the probabilistic setting, howe ver , uncertain injections and topology parameters in general act as a joint disturbance with possible cov ariance. W e therefore collect all perturbations into a single real-v alued stacked disturbance vector , ∆ ρ =   ∆ u v ec(∆ G bus ) v ec(∆ B bus )   , (46) where ∆ Y bus = ∆ G bus + j ∆ B bus . As in Section II-B5, uncertainty in Y bus is restricted to a minimal basis that spans the support of admissible network uncertainty such that ∆ Y bus = K X k =1 δ y k E k , (47) where { E k } K k =1 identify the nonzero perturbation direc- tions permitted by the uncertainty model, and δ Y = [ δ y 1 , . . . , δ y K ] ⊤ is the associated vector of random coefﬁ- cients. This basis is not assumed orthogonal or unique, and serves only to identify the support of admissible network per- turbations. Separating real and imaginary components yields v ec(∆ G bus ) = E G δ G , v ec(∆ B bus ) = E B δ B , (48) where E G and E B are ﬁxed embedding matrices, and δ G , δ B collect the real-v alued network uncertainty parameters. DRAFT FOR JOURNAL SUBMISSION, FEB 2026 10 For a scalar monitored quantity x , the linearized response can then be written as ∆ x = L x ∆ u + v ec( T x,G ) ⊤ v ec(∆ G bus ) + vec( T x,B ) ⊤ v ec(∆ B bus ) , = L x ∆ u + S x,G δ G + S x,B δ B , (49) where T x,G = ∂ x/∂ G bus and T x,B = ∂ x/∂ B bus are e valu- ated at the operating point associated with the current regime, and S x,G = v ec( T x,G ) ⊤ E G and S x,B = v ec( T x,B ) ⊤ E B . Collecting terms yields the compact representation ∆ x = S x ∆ ρ , (50) with S x =  L x v ec( T x,G ) ⊤ v ec( T x,B ) ⊤  . If ∆ ρ has ﬁnite ﬁrst and second moments, E [∆ ρ ] = µ ρ , Co v (∆ ρ ) = Σ ρ , (51) then by linearity in (50), the induced moments of ∆ x are µ x = E [∆ x ] = S x µ ρ , (52) σ 2 x = V ar(∆ x ) = S x Σ ρ S ⊤ x . (53) B. Moment-Based Constraint Mar gins Consider x with limits L min ≤ x ≤ L max , base value x ⋆ , and margin M ⋆ x ≥ 0 . Under a random perturbation ∆ ρ , the linearized response ∆ x = S x ∆ ρ results in a constraint violation whenever s ∆ x ≥ M ⋆ x . (54) For a binding constraint with M ⋆ x = 0 , x is released when s ∆ x < 0 . (55) T o support robustness screening, we deﬁne two moment-based margins. First, we deﬁne the σ -normalized margin as R σ,x = M ⋆ x − s µ x σ x , (56) which measures the distance between the limit and the mean drift of x , in standard deviations of ∆ x . W e next deﬁne the one-sided upper-tail α -quantile [43] of s ∆ x by q α,x such that Pr( s ∆ x ≤ q α,x ) = α ⇐ ⇒ Pr( s ∆ x ≥ q α,x ) = 1 − α. (57) W e deﬁne the α -conﬁdence margin as M α,x = M ⋆ x − q α,x , (58) so that M α,x > 0 implies Pr( s ∆ x ≤ M ⋆ x ) ≥ α . Using the ﬁrst two moments, we have q α,x = s µ x + ξ α σ x . Then M α,x = M ⋆ x − s µ x − ξ α σ x . (59) a) Choice of One-sided Quantile: The formulation above applies for any uncertain disturbance with ﬁnite ﬁrst and sec- ond moments. When only these moments are known, Cantelli’ s inequality [44] provides a distribution-free upper bound on the one-sided tail of s ∆ x . For any α ∈ (0 , 1) , Pr  s ∆ x ≥ s µ x + r α 1 − α σ x  ≤ 1 − α. (60) Equiv alently , with probability at least α , s ∆ x ≤ s µ x + ξ ( C ) α σ x , ξ ( C ) α = r α 1 − α . (61) When the distribution of ∆ x is known, the one-sided score ξ α may be ev aluated directly from the corresponding cumulativ e distribution function (CDF). If ∆ ρ is modeled as jointly normal, then s ∆ x is normally distrib uted and ξ ( N ) α = Φ − 1 ( α ) , (62) where Φ( · ) is the standard normal CDF . For example, ξ ( N ) 0 . 95 = 1 . 645 , recovers the normal 95% one-sided quantile. C. Distrib ution-F r ee Bounds on Constraint V iolation and Re- lease Pr obabilities Often, not enough data exist to quantify distrib utional un- certainties reliably . When only the ﬁrst and second moments are av ailable, we can nev ertheless obtain conservati ve bounds on violation and release probabilities. 1) Upper Bound on V iolation Pr obabilities: When sµ x < M ⋆ x , we can obtain an upper bound on the probability of constraint violation. Applying Cantelli’ s inequality giv es Pr  s ∆ x ≥ M ⋆ x  ≤ σ 2 x σ 2 x + ( M ⋆ x − s µ x ) 2 = 1 1 + R 2 σ,x . (63) The resulting bound remains valid for any underlying distribu- tion with the same moments, whether non-Gaussian, asymmet- ric, or heavy-tailed. The condition s µ x < M ⋆ x indicates that the mean drift remains inside the feasible region. If s µ x ≥ M ⋆ x , (63) reduces to the trivial bound Pr( s ∆ x ≥ M ⋆ x ) ≤ 1 , since it must also hold when σ x = 0 . 2) Lower Bound on Release Pr obabilities: For binding constraints when s µ x < 0 , we can calculate a lower bound on the release probability . Applying Cantelli’ s inequality giv es Pr( s ∆ x < 0) = 1 − Pr( s ∆ x ≥ 0) (64) ≥ 1 − 1 1 + R 2 σ,x = R 2 σ,x 1 + R 2 σ,x , (65) with M ⋆ x = 0 and thus R σ,x = − s µ x /σ x . If s µ x ≥ 0 , we get the trivial bound Pr( s ∆ x < 0) ≥ 0 . D. Constraint V iolation and Release Pr obabilities for Distur- bance Distributions When the joint distribution of the uncertain disturbance vector ∆ ρ is deﬁned, the afﬁne map (50) induces well- deﬁned distrib utions for ∆ x . This enables exact ev aluations of constraint violation and release probabilities. DRAFT FOR JOURNAL SUBMISSION, FEB 2026 11 1) V iolation Pr obabilities for Gaussian Uncertainty: A common assumption in PLF and OPF studies is that ∆ ρ follows a multi variate Gaussian distrib ution [45]. Under this assumption, the afﬁne response ∆ x is itself Gaussian, with mean µ x and variance σ 2 x as in (52) and (53). Under Gaussian uncertainty , the violation probability is Pr( s ∆ x ≥ M ⋆ x ) = Pr  s ∆ x − s µ x σ x ≥ M ⋆ x − s µ x σ x  = 1 − Φ( R σ,x ) . (66) 2) Release Pr obabilities for Gaussian Uncertainty: For binding constraints, the release probability is Pr( s ∆ x < 0) = Pr  s ∆ x − s µ x σ x < − s µ x σ x  = Φ( R σ,x ) . (67) E. Linear Quasi-Dual Uncertainty Mapping In Section II-C, deterministic quasi-duals are constructed using an already calculated z ′ by projecting the cost gradient onto the span of constraint normals activ e at that point. In the probabilistic setting, howe ver , z ′ is random, and both the activ e constraint set and the local constraint geometry vary across realizations. As a result, ˜ ν ( z ′ ) is not directly amenable to linearization with respect to z ′ ov er the full uncertainty support. T o obtain a tractable approximation, we instead prop- agate uncertainty from ∆ ρ into the quasi-dual space using a linearization ev aluated at a representativ e operating point, and restrict the resulting quasi-dual variations to domains where the corresponding constraints are applicable. For compactness, deﬁne the projected cost gradient q ( z ) = C ( z ) ∇ z f ( z ) , (68) and the normal matrix H ( z ) = C ( z ) C ( z ) ⊤ , (69) so that, within a region with a ﬁxed set of activ e constraints, the quasi-duals satisfy ˜ ν ( z ) = − H ( z ) − 1 q ( z ) . (70) Let ¯ z = E [ z ′ ] ≈ z ⋆ + E [∆ z ] (71) denote the mean perturbed state obtained from the linear primal uncertainty mapping. While the active set is ﬁxed, the quasi-dual mapping (70) is smooth. The dif ferential is d ˜ ν = − d  H − 1  q − H − 1 d q . (72) Using the identity d  H − 1  = − H − 1 ( d H ) H − 1 and substi- tuting q = − H ˜ ν , we obtain d ˜ ν = − H − 1 d q − H − 1 ( d H ) ˜ ν . (73) Next, differentiating (68)–(69) giv es d H = d ( CC ⊤ ) = ( d C ) C ⊤ + C ( d C ) ⊤ , (74) d q = d ( C ∇ z f ) = ( d C ) ∇ z f + C d ( ∇ z f ) . (75) Both d C and d ( ∇ z f ) depend linearly on d z through second deriv ati ves (constraint and cost Hessians). Writing the differ - ential columnwise, for each state component z k , ∂ ν ∂ z k = − H − 1 " ∂ q ∂ z k + ∂ H ∂ z k ν # , (76) with ∂ H ∂ z k = ∂ C ∂ z k C ⊤ + C  ∂ C ∂ z k  ⊤ , (77) ∂ q ∂ z k = ∂ C ∂ z k ∇ z f + C ∂ 2 f ∂ z ∂ z k . (78) Collecting (76) ov er all k gi ves the compact form d ˜ ν = K ( ¯ z ) ˜ ν d z , (79) where K ( ¯ z ) ν is ev aluated at z = ¯ z and uses the activ e set A ( ¯ z ) . Using the linearized state response ∆ z = S z ∆ ρ , we get ∆ ˜ ν ≈ K ˜ ν ∆ ρ , K ˜ ν = K ( ¯ z ) ˜ ν S z , (80) which describes the linear propagation of uncertainty into the quasi-dual space under the assumption of ﬁxed constraint geometry . Under this approximation, the ﬁrst and second moments of ∆ ν follow directly from those of ∆ ρ : E [∆ ˜ ν ] = K ˜ ν µ ρ , (81) Co v (∆ ˜ ν ) = K ˜ ν Σ ρ K ⊤ ˜ ν . (82) The linear mapping (80) is valid only within regions where the activ e set remains unchanged, the regime remains un- changed, and the ﬁrst-order approximation holds. W e do not attempt to quantify the size of this region or the probability of active-set changes here. 1) Inequality Quasi-Dual Domain Restriction: Quasi-duals associated with inequality constraints are physically mean- ingful only when the corresponding constraints are binding or violated. W e then restrict each component of ∆ ˜ ν to an admissible domain deﬁned by the corresponding constraint margin. Deﬁne the admissible domain of ˜ ν x as D x := { ∆ ρ : s ∆ x ≥ M ⋆ x } . (83) Since ∆ x = S x ∆ ρ , the domain D x corresponds to a half- space in ∆ ρ . The truncated quasi-dual variation associated with constraint x is then deﬁned as ∆ ¯ ν x = ( ∆ ˜ ν x , if ∆ ρ ∈ D x , 0 , otherwise . (84) Equiv alently , truncation occurs at the hyperplane s ∆ x = M ⋆ x , which serves as the activ ation threshold for the corresponding quasi-dual component. D x is ev aluated using the same linear response ∆ x . When the distribution of ∆ x is known, the probability of activ ation Pr( D x ) can be computed directly . When only ﬁrst and second moments are av ailable, conserv ative bounds on Pr( D x ) can be obtained using distribution-free inequalities, as discussed in Section IV -C. DRAFT FOR JOURNAL SUBMISSION, FEB 2026 12 F . Pr obabilistic Re gime Changes Sections II-B3–II-B7 resolve generator and network dis- continuities through regime changes associated with activ e power limits, network connectivity , and reacti ve po wer limits. These regime classes dif fer in how a unique ordering of ev ents can be deﬁned. In the deterministic setting, all detected activ e power saturation ev ents are treated as simultaneous, and higher-order network connectivity changes are likewise modeled as a single event by deﬁnition of N − k contingencies. Reactiv e po wer saturation admits a natural ordering using directional voltage sensitivities and predicted ov ershoots to deﬁne a unique sequence of Q -limit acti vations as generators transition from P V to P Q operation. For uncertain disturbances, these deterministic sequences are only a subset of the possible regime ev olutions. Dif ferent orderings and combinations of regime changes may occur with non-negligible probability , based on how probability masses interact with operational limits and connectivity thresholds. Ignoring e vent order , each generator with potential acti ve- or reactiv e power saturation admits three possible states (lower limit, interior , or upper limit), while each candidate discon- nection admits two states (connected or disconnected). The number of regime patterns is thus bounded by N pattern ≤ 3 N P + N Q 2 N Y , (85) where N P and N Q denote the numbers of generators with possible P - and Q -saturation e vents, and N Y the number of branches with possible disconnection e vents. If e vent order is also considered, the combinatorial space grows substantially . Allo wing any subset of the N P + N Q + N Y candidate events to occur in arbitrary sequence giv es at most N sequence ≤ N P + N Q + N Y X k =0 ( N P + N Q + N Y )! ( N P + N Q + N Y − k )! (86) sequences, which grows superexponentially with system size [46]. Exhaustive enumeration quickly becomes intractable. T o address this, we introduce two order-reduction methods. For distrib utional uncertainties, we construct a sequentially- pruned regime tree that enumerates physically consistent and probabilistically rele vant re gime ev olutions while preserving exact probability mass. When only moments are kno wn or used, we construct an α -stressed regime for e valuating con- servati ve bounds on violation and release probabilities. 1) Sequential-Pruning Mixtur e Representation Under Known Disturbance Distributions: Assume a known distribution for ∆ ρ . W e construct a regime tree where each tree node n represents a regime R n deﬁned by a conditioning region in disturbance space, together with a regime-local conditional disturbance ∆ ρ ( n ) ∼ ∆ ρ | R n and an associated linear response model S ( n ) . An edge from parent node m to child node n corresponds to a single regime-change event and carries conditional probability p n | m = Pr( R n | R m ) , so that each node has unconditional probability mass γ n = γ m p n | m , with γ 0 = 1 for the root node. T ree construction proceeds sequentially by regime class. At each non-terminal node m , a single candidate regime-change Fig. 1: Hierarchical composition of the sequentially-pruned regime tree, illustrating the ordered attachment of sub-trees. ev ent is selected and tested, inducing a binary split of R m into two complementary outcomes. If both resulting child probabilities e xceed a prescribed pruning threshold ε p , both children are added as nodes. If only one outcome exceeds ε p , the parent node is updated in place to represent the surviving conditional regime, and the complementary low- probability outcome is collapsed into a probability atom to the corresponding limit ev ent. This path-compression strategy preserves probability mass while minimizing tree gro wth. A node is terminal for a gi ven regime class when no further candidate events or controllable elements remain. The full regime tree is assembled by attaching sub-trees in the same regime-class order as in Section II-B7. A ∆ Y bus sub-tree is constructed ﬁrst to resolv e connectivity changes. T o each of its terminal nodes, a ∆ P sub-tree is attached to represent acti ve power saturation. Finally , ∆ Q sub-trees are attached to capture reactiv e power saturation and P V → P Q transitions. Each leaf of the resulting composite tree corresponds to a realized sequence of regime changes and accumulates the conditional probabilities along its path from the root. This structure is illustrated in Figure 1. Each node implicitly represents a conditioning region R n deﬁned by the intersection of branching ev ents along its path. The resulting tree does not enumerate all possible regime orderings, nor does it attempt to represent exact regime-local probability densities. It provides an accounting of physically consistent regime evolutions whose conditional probability exceeds the prescribed tolerance. a) ∆ Y bus -T r ee Construction: Let ∆ Y ( m ) bus denote the regime-local admittance perturbation at node m . Using a minimal basis that spans the support of admissible network uncertainty , ∆ Y ( m ) bus = P K m k =1 δ y ( m ) k E ( m ) k . Uncertainty in ∆ Y ( m ) bus is then characterized by the distribution of δ ( m ) Y , with conditional mean and cov ariance E [ δ ( m ) Y | R m ] = µ ( m ) Y , Co v ( δ ( m ) Y | R m ) = Σ ( m ) Y . (87) For a candidate branch ( i, j ) , let Y pred , ( m ) ij = Y ij + ∆ Y ( m ) ij , with ∆ Y ( m ) ij = ∆ G ( m ) ij + j∆ B ( m ) ij . A near-zero admittance (disconnection) event is deﬁned as DRAFT FOR JOURNAL SUBMISSION, FEB 2026 13 E ( m ) ij =     Y pred , ( m ) ij    2 ≤ Y 2 min  , (88) with   Y pred , ( m ) ij   2 =  G ij + ∆ G ( m ) ij  2 +  B ij + ∆ B ( m ) ij  2 . Each connectivity event corresponds to a quadratic inequality in the jointly random variables (∆ G ( m ) ij , ∆ B ( m ) ij ) . Because these variables may be correlated with other perturbation elements, the resulting probability does not admit a closed- form expression and is ev aluated via numerical integration. The disconnection probability at node m is p ( m ) ij = Pr  E ( m ) ij | R m  , and the candidate set of connecti vity regime ev ents is C ( m ) Y = n ( i, j )    p ( m ) ij > ε Y o , (89) for some signiﬁcance threshold ε Y . If C ( m ) Y = ∅ , no con- nectivity branching occurs at node m . Otherwise, the branch ( c i, j ) ∈ C ( m ) Y with the largest disconnection probability p ( m ) ˆ ij is selected for branching. T w o candidate outcomes are deﬁned: E ( m ) disc = n   Y pred , ( m ) ij   2 ≤ Y 2 min o , and (90) E ( m ) conn = n   Y pred , ( m ) ij   2 > Y 2 min o , (91) with regime-conditional edge probabilities p disc | m = p ( m ) ˆ ij , p conn | m = 1 − p ( m ) ˆ ij . (92) The corresponding node probabilities are γ disc = γ m p disc | m and γ conn = γ m p conn | m . If both resulting edge probabilities exceed the pruning threshold ε p , two child nodes are created. In the disconnection child, the branch admittance is clamped to zero, Y (disc) ˆ ij = 0 . The corresponding basis directions for ∆ G ˆ ij and ∆ B ˆ ij are remov ed from ∆ Y ( m ) bus , and the coef ﬁcient vector δ ( m ) Y is reduced. The mean and covariance of the remaining coef- ﬁcients are updated by restriction. If islanding occurs, the activ e system is reduced to the reachable sub-network. In the connected child, the network topology is preserved, and δ ( m ) Y is conditioned on E ( m ) conn . This yields a truncated distrib ution for (∆ G ( m ) ˆ ij , ∆ B ( m ) ˆ ij ) while retaining all basis directions. If only one outcome yields a node probability exceeding ε γ , the parent regime is updated in place according to the survi ving outcome. The complementary probability mass is collapsed to the corresponding connectivity event. If disconnection sur- viv es, the basis reduction and network update are applied as needed. If connection surviv es, the admittance uncertainty is truncated accordingly . Let B m denote the inde x set of basis directions acti ve at node m , and let C m ⊆ B m denote those correspond- ing to branch admittances that have already been tested for disconnection. A connectivity basis direction is considered resolved once its associated branch has been selected and ev aluated for branching. A node m is considered terminal when ∀ ( i, j ) ∈ B m \ C m , p ( m ) ij ≤ ε Y , or when no further uncertain admittance directions remain. b) ∆ P -T r ee Construction: At node m , let ∆ P ( m ) inj ∈ R |N | denote the activ e power injection components of the regime-local disturbance ∆ ρ ( m ) . Deﬁne the aggregation vector h = 1 |N | , so that the residual activ e power imbalance is ∆ P ( m ) res = h ⊤ ∆ P ( m ) inj = X i ∈N ∆ P ( m ) i , (93) with ∆ P (0) res = ∆ P inj as in Section II-B3. If ∆ P ( m ) inj has conditional mean µ ( m ) P and covariance Σ ( m ) P , then µ ( m ) res = h ⊤ µ ( m ) P , (94) ( σ ( m ) res ) 2 = h ⊤ Σ ( m ) P h . (95) All participating generators respond to the same scalar imbalance. For generator g , ∆ P ( m ) g = − α ( m ) g ∆ P ( m ) res with α ( m ) g > 0 , and the signed margin to an active po wer limit is M ( m ) g , lim = s  P lim − ( P g + ∆ P ( m ) g )  . Substituting the participation response yields the scalar threshold condition M ( m ) g , lim ≤ 0 ⇐ ⇒ s ∆ P ( m ) res ≤ s τ ( m ) g , lim , (96) with τ ( m ) g , lim = ( P g − P lim ) /α ( m ) g . Thus, all acti ve power limit ev ents correspond to parallel half-spaces in the scalar ∆ P ( m ) res . For each participating generator and limit, the conditional violation probability is p ( m ) g , lim = Pr  M ( m ) g , lim ≤ 0    R m  , (97) and the candidate set of activ e power regime events is C ( m ) P =  ( g , lim) : p ( m ) g , lim > ε P  , (98) for threshold ε P . If C ( m ) P = ∅ , no active po wer branching occurs at node m . Otherwise, the e vent (ˆ g, c lim) with the largest p ( m ) g , lim is selected for branching. Branching is induced by conditioning on the complementary ev ents E ( m ) lim =  s ∆ P ( m ) res ≤ s τ ( m ) ˆ g , c lim  , E ( m ) int =  E ( m ) lim  c , (99) with conditional probabilities p lim | m = p ( m ) ˆ g , c lim , p int | m = 1 − p ( m ) ˆ g , c lim . (100) If both resulting edge probabilities exceed the pruning threshold ε p , two child nodes are created. In the interior child, ∆ P ( m ) res is truncated to E ( m ) int . In the limit child, generator ˆ g is clamped at P c lim ˆ g , removed from the participation set, and participation factors are renormalized. The post-clamping imbalance is ∆ P (lim) res ∼  ∆ P ( m ) res | E ( m ) lim  + α ( m ) ˆ g τ ( m ) ˆ g , c lim , (101) with subsequent branching using this updated distrib ution. If only one outcome exceeds ε p , the parent node is updated in place to represent the surviving regime, with the com- plementary probability mass collapsed into the corresponding saturation e vent. Here, ∆ P ( m ) res remains unchanged; it is only conditioned upon branching. DRAFT FOR JOURNAL SUBMISSION, FEB 2026 14 c) ∆ Q -T r ee Construction: The construction of the ∆ Q - tree follows the same sequential-pruning logic as the ∆ P -tree, differing in the structure of the reactive-po wer response and the physical consequences of saturation. At node m , the regime-local disturbance ∆ ρ ( m ) induces a linearized reactiv e-power response ∆ Q ( m ) = S ( m ) Q ∆ ρ ( m ) , where S ( m ) Q depends on the current network topology , bus-type assignments, and set of voltage-controlling generators. For each generator q operating in P V mode, the predicted reacti ve- power output is Q pred , ( m ) q = Q q + ∆ Q ( m ) q , with signed margin M ( m ) q , lim = s  Q lim − Q pred , ( m ) q  . V iolation probabilities p ( m ) q , lim , candidate sets C ( m ) Q , ev ent selection, and branching e vents E ( m ) lim = { M ( m ) q , lim ≤ 0 } and E ( m ) int = ( E ( m ) lim ) c are deﬁned the same way as in the ∆ P -tree. In the limit child, generator ˆ q is clamped at Q c lim ˆ q , and the bus is reclassiﬁed from P V to P Q if no other P V generators remain at that bus. The linearized network equations are reconstructed, and a new sensitivity mapping S (lim) Q is formed before further branching. In the interior child, the generator remains in P V mode and the disturbance distribution is conditioned on E ( m ) int . If only one outcome exceeds ε p , the parent node is updated in place to represent the survi ving regime. If the limit outcome surviv es, the P V → P Q transition is enforced and the comple- mentary probability mass is collapsed into the violated limit. If the interior outcome survi ves, the disturbance distribution is truncated accordingly . Branching continues until no remaining P V generators exhibit Q -limit violation probabilities exceeding threshold ε Q . d) Mixtur e Repr esentation and Evaluation of Leaf Mo- ments and Pr obabilities: The sequentially-pruned tree pro- duces a ﬁnite set of leaf regimes {R ℓ } L ℓ =1 with associated unconditional weights γ ℓ satisfying P L ℓ =1 γ ℓ = 1 . Each leaf ℓ corresponds to a fully speciﬁed operating regime and an explicit conditioning region R ℓ obtained as the intersection of all edge e vents along the unique root-to-leaf path. In leaf ℓ , the regime-local disturbance is deﬁned by the conditional law ∆ ρ | R ℓ , and deterministic regime edits which act on additi ve perturbation coordinates through projection, redistribution, or clamping operations. Since these operations are linear or constant-shift transformations on ﬁrst-order de- viation variables, their effect on the surviving disturbance can be represented by an af ﬁne transformation ∆ ρ ( ℓ ) = A ℓ ∆ ρ + b ℓ , ∆ ρ ∼ (∆ ρ | R ℓ ) , (102) where A ℓ encodes projection, redistribution, or reparameteri- zation, and b ℓ accounts for ﬁxed offsets introduced by clamp- ing. Conditioning itself and the collapse of low-probability complementary outcomes into probability atoms act at the lev el of probability measures and are handled separately such that af ﬁne maps are applied only to the propagated random component. The leaf-local moments are then µ ( ℓ ) ρ = E [∆ ρ ( ℓ ) ] = A ℓ µ ρ | ℓ + b ℓ , (103) Σ ( ℓ ) ρ = Co v (∆ ρ ( ℓ ) ) = A ℓ Σ ρ | ℓ A ⊤ ℓ , (104) where µ ρ | ℓ = E [∆ ρ | R ℓ ] and Σ ρ | ℓ = Co v (∆ ρ | R ℓ ) . Proba- bility mass accounting is handled entirely by the conditioning region R ℓ , while the deterministic edits ( A ℓ , b ℓ ) only modify the effecti ve disturbance representation within the leaf. For any monitored scalar quantity x , leaf ℓ contains a deterministic linear response map ∆ x ( ℓ ) = S ( ℓ ) x ∆ ρ ( ℓ ) , with S ( ℓ ) x computed from the physical model in that regime. Con- sequently , the induced leaf-local moments are µ ( ℓ ) x = E [∆ x ( ℓ ) ] = S ( ℓ ) x µ ( ℓ ) ρ , (105) σ 2 , ( ℓ ) x = V ar(∆ x ( ℓ ) ) = S ( ℓ ) x Σ ( ℓ ) ρ S ( ℓ ) ⊤ x . (106) Deﬁne the leaf-wise violation probability as p ( ℓ ) viol ,x := Pr  s ∆ x ( ℓ ) ≥ M ⋆ x  , (107) with p ( ℓ ) rel ,x = 1 − p ( ℓ ) viol ,x . The overall probabilities are obtained by mixture aggregation over leaves, p viol ,x = L X ℓ =1 γ ℓ p ( ℓ ) viol ,x , p rel ,x = L X ℓ =1 γ ℓ p ( ℓ ) rel ,x . (108) Lastly , in each node regime, we can compute ¯ z ( l ) for ev aluating quasi-duals and deﬁne mixture-model quasi-duals ov er the tree over the union or intersection of activ e constraint sets, as desired. W e omit these details for bre vity . e) Combinatorial Reduction: Sequential pruning con- structs a reduced binary regime tree whose growth is tied to the number of regime ev ents realized along each path. In the worst case, the number of leaf regimes is bounded by R ≤ 2 N P + N ∆ Y + N Q . (109) In practice, the realized tree is far smaller . Only events with conditional violation probabilities e xceeding prescribed tolerances are admitted as branching candidates, suppressing growth driven by distribution tails. Conditioning at each interior node further reduces the likelihood of do wnstream ev ents, while clamped outcomes rapidly diminish the remain- ing probability mass. Moreov er , each P V → P Q transition permanently removes a generator from subsequent reactive- power consideration. If additional reduction is required, further approximations may be applied. For example, pruning can be based on cumu- lativ e node probability γ n rather than conditional probabilities p n | m , or enforced through structural limits on tree breadth, depth, or total node count. Such choices trade enforcement of hard constraints on generators for computation time. W e recommend always checking for islanding. 2) Appr oximation of Distribution-F ree Cantelli Bounds in an α -Str essed Regime: The sequentially-pruned regime-tree abov e is feasible when the distrib ution of ∆ ρ is speciﬁed. When only ﬁrst and second moments of ∆ ρ are available, neither the node probabilities γ ℓ nor the conditional proba- bilities p n | m can be identiﬁed uniquely , since many distinct distributions can share the same mean and cov ariance but induce different probabilities of regime transitions. T o obtain a tractable, distribution-free alternati ve, we con- struct a single α -str essed re gime . Using only ﬁrst and second DRAFT FOR JOURNAL SUBMISSION, FEB 2026 15 moments we identify which regime changes are unavoidable at a prescribed conﬁdence lev el α , and deterministically enforce those changes prior to ev aluating constraint margins. All remaining constraints are then assessed within this single α -stressed regime using Cantelli bounds. This construction yields conservati ve guarantees while avoiding combinatorial regime enumeration. This method can be applied even if the distribution of ∆ ρ is known. Starting from the base operating point and disturbance moments ( µ ρ , Σ ρ ) , the α -stressed regime approximation pro- ceeds in three stages, follo wing the same ordering of regime classes as the deterministic and mixture-based constructions. Each stage applies an α -conﬁdence test to all candidate ev ents of that class simultaneously . All tests that cannot be certiﬁed at conﬁdence α are conserv atively enforced. The system model, disturbance representation, and linearized response maps are then updated deterministically before proceeding to the next stage. a) ∆ Y bus α -T est: For each branch ( i, j ) with uncertain admittance, let Y pred ij =  G ij + ∆ G ij B ij + ∆ B ij  (110) be the predicted admittance vector . Its mean and covariance, ( µ ij , Σ ij ) , are obtained as afﬁne projections of ( µ ρ , Σ ρ ) . Disconnection is deﬁned by the ev ent ∥ Y pred ij ∥ 2 ≤ Y min . W ith only ﬁrst and second moments, this quadratic ev ent cannot be ev aluated exactly . Instead, we create a distrib ution- free α -conﬁdence certiﬁcate. Using the expectation of a quadratic form [47], we hav e E     Y pred ij − µ ij    2 2  = tr( Σ ij ) . (111) Applying Markov’ s inequality [48] to ∥ Y pred ij − µ ij ∥ 2 2 giv es Pr     Y pred ij − µ ij    2 2 ≥ tr( Σ ij ) 1 − α  ≤ 1 − α, (112) or equiv alently , Pr     Y pred ij − µ ij    2 ≤ r α,ij  ≥ α, r α,ij = r tr( Σ ij ) 1 − α . (113) Thus, with conﬁdence at least α , the random vector Y pred ij lies inside the Euclidean ball of radius r α,ij around its mean. By the triangle inequality [49], ∥ µ ij ∥ 2 =    Y pred ij − ( Y pred ij − µ ij )    2 ≤ ∥ Y pred ij ∥ 2 + ∥ Y pred ij − µ ij ∥ 2 , = ⇒ ∥ Y pred ij ∥ 2 ≥ ∥ µ ij ∥ 2 − ∥ Y pred ij − µ ij ∥ 2 . (114) Combining (114) with the α -conﬁdence ev ent in (113), with probability at least α , ∥ Y pred ij ∥ 2 ≥ ∥ µ ij ∥ 2 − r α,ij . (115) W e thus deﬁne the α -conﬁdence connectivity margin M α,Y ij = ∥ µ ij ∥ 2 − q tr( Σ ij ) 1 − α − Y min . (116) If M α,Y ij > 0 , branch ( i, j ) is certiﬁed to remain connected with conﬁdence at least α . Otherwise, disconnection cannot be excluded at le vel α and the branch is removed. All branches with M α,Y ij ≤ 0 are simultaneously discon- nected. The network is then reduced to the subgraph reachable from the slack bus. Corresponding admittance uncertainty directions are remov ed from ∆ ρ , and the disturbance mo- ments ( µ ρ , Σ ρ ) are updated by restriction to the remaining degrees of freedom. The network model and linearization are reconstructed before proceeding. b) ∆ P α -T est: Let ∆ P res = h ⊤ ∆ P inj denote the net activ e power imbalance, with mean µ res and v ariance σ 2 res . For each participating generator g with participation factor α g , the predicted response is ∆ P g = − α g ∆ P res . For each activ e power limit, deﬁne the signed margin M ⋆ g , lim as in the deterministic case. Using the moment-based α -conﬁdence margin, M α,g , lim = M ⋆ g , lim − s µ g − ξ α σ g , all generators for which M α,g , lim ≤ 0 are clamped at the corresponding limit. Let G ( P ) sat denote the set of clamped generators. These are remov ed from the participation set, and participation factors are renormalized ov er the remaining generators. The residual imbalance is updated afﬁnely by subtracting the deterministic contributions of the clamped units, which shifts µ res but lea ves σ res unchanged. The activ e power control model and linear response mappings are then updated. c) ∆ Q α -T est: For each generator operating in P V mode, the linearized reactive power response is ∆ Q q = S q ∆ ρ , with moments ( µ q , σ q ) computed from the current disturbance moments. For each reactive power limit, deﬁne the α -conﬁdence mar gin M α,q , lim = M ⋆ q , lim − s µ q − ξ α σ q . All generators with M α,q , lim ≤ 0 are clamped at their reactive limits and reclassiﬁed from P V to P Q mode. The bus-type assignment is updated, the linearized network equations are reconstructed, and the sensitivities are recomputed. d) Evaluation in the α -Str essed Re gime.: After applying all three classes of α -tests, the system is left in a single α -stressed operating regime with a fully speciﬁed network topology , generator participation set, and bus-type classiﬁca- tion. The disturbance representation ∆ ρ has updated moments ( µ ( α ) ρ , Σ ( α ) ρ ) , and all linear response mappings S ( α ) x are re- computed consistently with this regime. For any x , ∆ x = S ( α ) x ∆ ρ , with moments µ ( α ) x = S ( α ) x µ ( α ) ρ , ( σ ( α ) x ) 2 = S ( α ) x Σ ( α ) ρ S ( α ) ⊤ x . (117) V iolation and release probabilities are then bounded using Cantelli’ s inequality with these moments, yielding distribution- free guarantees that are valid for all disturbance distributions consistent with the speciﬁed ﬁrst and second moments. G. Rob ustness Scr eening and Meaningful Uncertainty Models Sections IV -A–IV -F de velop a framew ork that is distribution-agnostic. In practice, system-wide screening requires disturbance models that reﬂect credible operational uncertainty . Here we brieﬂy revie w commonly used disturbance models. DRAFT FOR JOURNAL SUBMISSION, FEB 2026 16 (a) N − 1 screening with preset voltage setpoints. (b) N − 1 screening with optimized voltage setpoints. Fig. 2: N − 1 contingenc y screening on Puerto Rican bulk power system model, around a base operating point representing peak load conditions with dispatch optimized using A C OPF economic dispatch with (a) preset and (b) optimized voltage setpoints. a) Independent Uncertainties: Independent disturbances across injections and branches take the block diagonal form ∆ ρ ∼ P  µ , diag (Σ 2 1 , . . . , Σ 2 N G + N Y )  (118) and hav e been explored extensi vely in both PLF and CC-OPF contexts [50, 51, 52]. Such models deﬁne Cartesian uncertainty sets and naturally support k th -order screening analyses, anal- ogous to N − k screening. For kno wn distrib utions, a system is N + P ( µ , Σ ) k robust if all admissible k th -order ∆ ρ satisfy prescribed violation and release probability tolerances. When only ﬁrst and second moments are av ailable, an analogous N + ( µ , Σ ) k notion applies. b) Load and Rene wable Uncertainty: Load, wind, and solar photovoltaic uncertainties exhibit strong temporal and spatial dependence dri ven by shared factors such as weather conditions, time of day , and geographic proximity . The litera- ture captures these correlations using joint models, including copula-based methods, mixed PDF/CDF formulations, Gaus- sian mixture models, and sampling approaches [53, 54, 55]. These joint models induce correlated disturbance directions in ∆ ρ , which can be sampled or approximated within the proposed screening frame work without modiﬁcation. c) Line P ar ameter and T opology Uncertainty: T ransmis- sion uncertainty arises from errors in line and transformer parameters, operating conditions, and network switching or outages. In the PLF literature, line parameters are treated as uncertain inputs to the po wer-ﬂow equations using prob- abilistic, bounded, or scenario-based models, while topology uncertainty is handled through discrete contingencies or out- age probabilities [45]. Here, continuous parameter uncertainty enters through perturbations of Y bus under ﬁxed topology , whereas discrete topology changes are handled explicitly through the regime and mixture constructions of Section IV -F. V . C A S E S T U DY : D E M O N S T R A T I O N O N P U E RT O R I C A N B U L K P O W E R S Y S T E M T o demonstrate ho w the proposed screening metrics can be communicated at scale, we apply them to a publicly av ailable A C power ﬂow model of Puerto Rico’ s bulk po wer system (BPS) ﬁrst created by [56] and visualize results geospatially . The model has 12 generators, 960 loads, 385 buses and 861 lines representing the 38 kV , 115 kV and 230 kV net- work, deriv ed from public data sources. W e implement minor modiﬁcations to limits, setpoints, and branch impedances as described in [57]. The 38 kV sub-transmission lev el is approx- imated using Kron reduction. As such, all reported violations, contingencies, and screening visualizations are restricted to the 115 kV and 230 kV BPS layer . Puerto Rico’ s BPS spans the island with major load centers including San Juan, Ponce, Mayag ¨ uez, Arecibo, and Caguas, while legac y thermal generation is largely concentrated near coastal ports. Bulk transfers commonly trav erse the mountainous interior [57]. W e consider a peak loading scenario, which is solved in [57] with an AC OPF economic dispatch formulation, using the SmartGridz solver [58] with both preset and adjustable voltage setpoints at P V buses. All generators are participating in P -balancing, in proportion to their MW capacities. First, we perform standard N − 1 screening on both cases with the settings mentioned, considering generator outages and outages on 115 kV branches and abov e. Figure 2 shows violation frequencies on buses, branches, and generators in both cases, illustrated geospatially (generators are aggregated to corresponding buses) and as a barplot for the top 10 components with the most limit violations. The side-by-side comparison shows that by optimizing voltage setpoints, the number of violations system-wide decreases and stress is more ev enly distributed across components. W e next extend the analysis to compare the array of screen- ing tools presented in this paper , on the optimized voltage DRAFT FOR JOURNAL SUBMISSION, FEB 2026 17 (a) Limit violation ( M x < 0 ) frequencies from N − 1 screening. (b) Quasi-dual vulnerabilities ( ˜ ν x > 0 ) from C − 1 screening. (c) Marginal vulnerabilities ( R x ∈ [ − 10 , 0) ) from N + δ 1 screening. (d) SPVI vulnerabilities ( ς x > 0 ) from C + δ 1 screening. (e) Max violation probabilities from N + N ( µ , Σ ) 1 screening. (f) Max violation probabilities from N + ( µ , Σ ) 1 screening. Fig. 3: Geospatial visualizations of system-wide screening metrics on the Puerto Rican bulk power system under a peak- load operating point solved by A C OPF with optimized voltage setpoints. Panels (a)–(d) sho w deterministic N − 1 , C − 1 , N + δ (1) , and C + δ (1) screenings, with the ten components exhibiting the highest frequencies highlighted. Panels (e)–(f) present probabilistic single-contingency screening under Gaussian disturbances using the sequentially pruned mixture model and the α -stressed regime method, respectiv ely . For probabilistic panels, the 100 largest violation probabilities are displayed, with generator probabilities aggregated to their corresponding b uses. DRAFT FOR JOURNAL SUBMISSION, FEB 2026 18 setpoint case. In particular , we perform: 1) N − 1 screening as above. 2) C − 1 screening via the N − 1 perturbed states, marking all positive inequality quasi-duals as vulnerable. 3) N + δ 1 screening, using aligned disturbance directions on the N − 1 component set, and ε R = 10 . 4) C + δ 1 screening, using the N + δ 1 directions, marking all positive inequality SPVI indices as vulnerable. 5) N + N ( µ , Σ ) 1 screening, using multiv ariate Gaussian dis- tributions with zero means and 0.02% of base injections and admittances as variances on the N − 1 component set. W e set Y min = 1 e − 5 and an equal signiﬁcance threshold of 5% for all condition tests. W e dra w 10,000 Monte Carlo samples for each instantiation of ∆ ρ ( m ) . 6) N + ( µ , Σ ) 1 screening, using the same moments as in the N + N ( µ , Σ ) 1 screening. W e consider α = 95% conﬁdence tests and Cantelli bounds. The results are shown in Figure 3 as a dashboard. All maps presented in this paper are screenshots of interactiv e maps that identify which uncertainties correspond to which exceedance ﬂags or probabilities per component. Each screening result highlights dif ferent aspects of system robustness. N − 1 and C − 1 screens (Figures 3a–3b) correspond to classical contingency analyses with regime transitions. They identify components prone to limit violations under outages and rev eal how frequently active constraints exert positive economic pressure via quasi-duals. In contrast, N + δ 1 and C + δ 1 screens (Figures 3c–3d) pro vide local analogs that highlight components which are directionally vulnerable to violations or shadow-price ampliﬁcation. The deterministic results are giv en using frequency of threshold exceedance, but could be sho wn instead using continuous measure such as mean or median magnitudes of exceedance. W e implement this with the probabilistic results in Figure 3e and Figure 3f by sho wing violation probabilities. Each plot sho ws, giv en uncertainty sets with the same ﬁrst and second moments, the maximum probability of violation per component over the uncertainty sets. As expected, the results align, with the α -stressed re gime method systematically providing more conservati ve estimates. Interestingly , across all screening methods, the system was found to be most vulnerable to b us voltage violations. V I . C O N C L U S I O N This paper presented a uniﬁed diagnostic frame work for assessing the robustness of AC operating points, by perturbing primal and dual states deriv ed from A C OPF solutions. W e dev elop a control- and constraint-aware steady-state response model to assess how operating points change in response to injection and topology perturbations, while considering regime changes dri ven by connectivity loss, generator saturation, and P V → P Q transitions. Propagating disturbances through a common response mapping enables consistent deterministic and probabilistic robustness metrics that generalize classical N − 1 screening to include economic indicators, local neigh- borhoods, and uncertain inputs. T o interpret economic stress away from optimality , we introduced quasi-dual variables as local geometric analogs of shadow prices, which are nonzero for activ e constraints, sim- ilar to dual variables deriv ed from optimization. Unlike dual variables, quasi-duals admit both positi ve and negati ve values at off-optimal operating points, distinguishing constraints that activ ely oppose cost reduction from those that are merely activ e or violated. This provides a complementary economic diagnostic that separates geometric proximity to limits from economic relev ance. The framework was extended to probabilistic robustness assessment under both distribution-based and moment-based uncertainty models. T o manage combinatorial regime changes under uncertainty , we dev eloped a sequentially pruned tree for mixture models when distrib utions are kno wn, and an α -stressed regime for when only ﬁrst and second moments are av ailable. These methods enable tractable estimation of violation and release probabilities while remaining consistent with physical regime transitions. A case study on a model of the Puerto Rican b ulk po wer sys- tem demonstrated ho w the proposed diagnostics can be used to increase situational awareness of system robustness. The sample visualizations selected from the broader set of results illustrated ho w different screening metrics emphasize com- plementary aspects of system stress, including contingency- driv en violations, directional vulnerability , economic tightness, and probabilistic risk, providing insight beyond feasibility and optimality alone. Sev eral directions for future work remain, including ex- tensions to more detailed steady-state voltage control logic, including A VC and P Q → P V transitions, higher-order response mappings, and improved probabilistic quasi-dual estimation. More broadly , the proposed framework provides a systematic means of auditing the robustness of A C-OPF solutions across scenario-based operating conditions, serving as a post-optimal diagnostic layer alongside primal and dual solutions for planners and operators. D E C L A R A T I O N O F C O M P E T I N G I N T E R E S T The authors declare that they hav e no competing ﬁnancial interests or personal relationships that could have appeared to inﬂuence the work reported in this paper . D A TA A V A I L A B I L I T Y The base po wer ﬂow model and A C OPF solutions used in the case study are available in the following GitHub repository: https://github .com/lauanton/puerto rico robustness case study. Accessed on: February 10, 2026. A C K N O W L E D G M E N T S The research presented in this manuscript was funded by the MIT Martin Family Society of Fello ws for Sustainability . The case study uses a publicly av ailable model of the Puerto Rican system described in [57]. The authors thank Rupamathi Jaddiv ada, Director of Innov ation at SmartGridz , for her support and initial guidelines on using the SmartGridz A C optimal power ﬂow solver . DRAFT FOR JOURNAL SUBMISSION, FEB 2026 19 R E F E R E N C E S [1] International Energy Agency (IEA), “W orld energy out- look 2024, ” tech. rep., International Energy Agency (IEA), Oct. 2024. [2] D. A yres and L. Zamora, “Renewable Power Generation Costs in 2023, ” tech. rep., International Renewable En- ergy Agency (IRENA), Abu Dhabi, U AE, 2024. ISBN: 978-92-9260-621-3. [3] North American Electric Reliability Corporation (NERC), “In verter-based resource performance issues public report, ” tech. rep., North American Electric Reliability Corporation (NERC), 2023. Level 2 alert, March 2023. [4] M. I. Fernandez, Y . I. Go, D. M. W ong, and W .-G. Fr ¨ uh, “Revie w of challenges and ke y enablers in energy systems tow ards net zero target: Renew ables, storage, buildings, & grid technologies, ” Heliyon , vol. 10, no. 23, p. e40691, 2024. [5] Ne w Y ork Independent System Operator (NYISO), “2025 power trends: The new york iso annual grid and markets report, ” tech. rep., Ne w Y ork Independent System Operator (NYISO), 2025. Annual state-of-the- grid and markets report. [6] International Energy Agency (IEA), “Climate resilience – power systems in transition, ” tech. rep., International Energy Agency (IEA), 2021. Analysis report, published 12 April 2021. [7] ISO New England (ISO-NE), “2024 regional electricity outlook, ” tech. rep., ISO New England, 2024. [8] ISO New England (ISO-NE), “Regional electricity out- look, ” tech. rep., ISO New England (ISO-NE), 2025. Accessed: October 29, 2025. [9] PJM Interconnection (PJM), “Energy transition in pjm: Resource retirements, replacements, and risks, ” tech. rep., PJM Interconnection (PJM), 2023. Special Report. [10] Midcontinent Independent System Operator (MISO), “2024 reliability imperativ e report, ” tech. rep., Midconti- nent Independent System Operator (MISO), 2024. Living Document: Updated February 2024, Accessed: October 30, 2025. [11] California Independent System Operator (CAISO), “Managing the ev olving grid, ” tech. rep., California In- dependent System Operator (CAISO), 2025. Accessed: October 29, 2025. [12] Electric Reliability Council of T e xas (ERCO T), “Report on existing and potential electric system constraints and needs, ” tech. rep., Electric Reliability Council of T exas (ERCO T), 2024. [13] Southwest Power Pool (SPP), “2025 spp summer re- source adequacy report, ” tech. rep., Southwest Power Pool (SPP), 2025. [14] LUMA Energy, “Luma’ s report: Q4 metrics sho w lack of funding leading to longer and more frequent outages, ” tech. rep., LUMA Ener gy, 2025. Accessed: October 30, 2025. [15] LUMA Energy, “New resource adequacy report: Lack of reliable po wer supply from generators will pose ongoing threat to service reliability in 2025, ” tech. rep., LUMA Energy, 2024. Accessed: October 30, 2025. [16] D. Bienstock, M. Chertkov , and S. Harnett, “Chance- Constrained Optimal Power Flow: Risk-A ware Network Control under Uncertainty , ” SIAM Review , vol. 56, no. 3, pp. 461–495, 2014. Publisher: Society for Industrial and Applied Mathematics. [17] A. Esteban-P ´ erez and J. M. Morales, “Distributionally robust optimal power ﬂow with contextual information, ” Eur opean Journal of Operational Researc h , vol. 306, no. 3, pp. 1047–1058, 2023. [18] Y . Guo, K. Baker , E. Dall’Anese, Z. Hu, and T . Summers, “Stochastic optimal po wer ﬂow based on data-dri ven distributionally robust optimization, ” in 2018 Annual American Contr ol Confer ence (ACC) , pp. 3840–3846, 2018. [19] Po werW orld Corporation, “Linear analysis techniques in powerw orld simulator, ” tech. rep., PowerW orld Corpora- tion, Champaign, IL, 2015. Accessed: October 31, 2025. [20] A. L. da Silv a, R. Allan, S. Soares, and V . Arienti, “Probabilistic load ﬂo w considering network outages, ” IEE Pr oceedings C (Generation, T ransmission and Dis- tribution) , vol. 132, pp. 139–145, 1985. [21] B. Borko wska, “Probabilistic load ﬂo w , ” IEEE T rans- actions on P ower Apparatus and Systems , vol. P AS-93, no. 3, pp. 752–759, 1974. [22] J. Dopazo, O. Klitin, and A. Sasson, “Stochastic load ﬂows, ” IEEE T ransactions on P ower Apparatus and Systems , vol. 94, no. 2, pp. 299–309, 1975. [23] J. Usaola, “Probabilistic load ﬂow in systems with wind generation, ” IET Generation, T ransmission & Distrib u- tion , vol. 3, pp. 1031–1041, 2009. [24] G. Plattner , H. F . Semlali, and N. Kong, “ Analysis of probabilistic load ﬂow using point estimation method to ev aluate the quantiles of electrical networks state variables, ” CIRED , vol. 2017, pp. 2087–2091, 2017. [25] T . Karakatsanis and N. Hatziargyriou, “Probabilistic con- strained load ﬂow based on sensiti vity analysis, ” IEEE T ransactions on P ower Systems , vol. 9, no. 4, pp. 1853– 1860, 1994. [26] V . Ajjarapu and C. Christy , “The continuation power ﬂow: a tool for steady state voltage stability analysis, ” IEEE T r ansactions on P ower Systems , v ol. 7, no. 1, pp. 416–423, 1992. [27] N. Hatziargyriou and T . Karakatsanis, “Probabilistic load ﬂow for assessment of voltage instability , ” IEE Pr o- ceedings - Generation, T ransmission and Distribution , vol. 145, pp. 196–202, 1998. [28] R. Billinton and R. N. Allan, Reliability Evaluation of P ower Systems . New Y ork, NY : Springer , 1 ed., 1996. [29] P . Chen, Z. Chen, and B. Bak-Jensen, “Probabilistic load ﬂow: A revie w , ” in 2008 Thir d International Confer ence on Electric Utility Der e gulation and Restructuring and P ower T echnologies , pp. 1586–1591, 2008. [30] D. D. Le, A. Berizzi, C. Bov o, E. Ciapessoni, D. Cirio, A. Pitto, and G. Gross, “ A probabilistic approach to power system security assessment under uncertainty , ” in 2013 IREP Symposium Bulk P ower System Dynamics and DRAFT FOR JOURNAL SUBMISSION, FEB 2026 20 Contr ol - IX Optimization, Security and Contr ol of the Emer ging P ower Grid , pp. 1–7, 2013. [31] A. Dimitrovski and K. T omsovic, “Slack bus treatment in load ﬂow solutions with uncertain nodal powers, ” in 2004 International Conference on Probabilistic Methods Applied to P ower Systems , pp. 532–537, 2004. [32] A. M. Hou and L. A. Roald, “Chance constraint tuning for optimal power ﬂo w , ” in 2020 International Confer- ence on Pr obabilistic Methods Applied to P ower Systems (PMAPS) , pp. 1–6, 2020. [33] Z. W ang, A. Y ounesi, M. V . Liu, G. C. Guo, and C. L. Anderson, “A C Optimal Power Flow in Power Systems W ith Renew able Ener gy Integration: A Revie w of Formulations and Case Studies, ” IEEE Access , v ol. 11, pp. 102681–102712, 2023. [34] J. D. Glov er, M. S. Sarma, and T . J. Overbye, P ower System Analysis and Design, F ifth Edition, SI . Global Engineering, 5th ed., 2012. [35] P . Kundur , N. Balu, and M. Lauby , P ower System Stabil- ity and Contr ol . EPRI po wer system engineering series, McGraw-Hill Education, 1994. [36] M. Ili ´ c and J. Zaborszky , Dynamics and Contr ol of Large Electric P ower Systems . Wile y , 2000. [37] W . Y an, W . Cui, W .-J. Lee, J. Y u, and X. Zhao, “Pilot- bus-centered automatic v oltage control with high pen- etration lev el of wind generation, ” IEEE T ransactions on Industry Applications , vol. 52, no. 3, pp. 1962–1969, 2016. [38] B. Stott, “Revie w of load-ﬂow calculation methods, ” Pr oceedings of the IEEE , vol. 62, no. 7, pp. 916–929, 1974. [39] Y . Lei, R. W ang, T . Li, Q. T ang, Y . W ang, and J. Li, “Modeling pv/pq switching in security constrained opti- mal power ﬂo w , ” in 2019 IEEE Inno vative Smart Grid T echnologies - Asia (ISGT Asia) , pp. 84–88, 2019. [40] A. V enkataramana, Computational T echniques for V olt- age Stability Assessment and Contr ol . Po wer Electronics and Power Systems, Boston, MA: Springer US, 2007. [41] J. Nocedal and S. J. Wright, Numerical Optimization . Springer , second edition ed., 2006. [42] C. M. Da vis and T . J. Ov erbye, “Multiple element contin- gency screening, ” IEEE T ransactions on P ower Systems , vol. 26, no. 3, pp. 1294–1301, 2011. [43] M. J. Evans and J. S. Rosenthal, Pr obability and Statis- tics: The Science of Uncertainty . W . H Freeman, 2003. [44] I. R. Sav age, “Probability inequalities of the tchebychef f type, ” 1952. Report Number: NBS Report 1744. [45] B. R. Prusty and D. Jena, “ A critical re view on prob- abilistic load ﬂo w studies in uncertainty constrained power systems with photov oltaic generation and a new approach, ” Renewable and Sustainable Ener gy Reviews , vol. 69, pp. 1286–1302, 2017. [46] G. Casella and R. L. Berger , Statistical Inference . Cen- gage Learning, second edition ed., 2002. [47] D. A. Kendrick, Stochastic Contr ol for Economic Models . McGraw-Hill, second edition ed., 1981. [48] M. Huber , “Halving the bounds for the markov , cheby- shev , and chernoff inequalities using smoothing, ” The American Mathematical Monthly , vol. 126, no. 10, pp. 915–927, 2019. [49] M. A. Khamsi and W . A. Kirk, An Intr oduction to Metric Spaces and F ixed P oint Theory . W iley-Interscience, ﬁrst edition ed., 2001. [50] Z. Hu and X. W ang, “ A probabilistic load ﬂow method considering branch outages, ” IEEE T r ansactions on P ower Systems , vol. 21, no. 2, pp. 507–514, 2006. [51] F . A. Hossain, M. Rokonuzzaman, N. Amin, J. Zhang, M. K. Mishu, W .-S. T an, M. R. Islam, and R. B. Roy , “Probabilistic load ﬂow–based optimal placement and sizing of distributed generators, ” Energies , vol. 14, no. 23, 2021. [52] M. V elay , M. V inyals, Y . Besanger, and N. Retiere, “Dis- tributed chance-constrained optimal power ﬂow based on primary frequenc y control, ” in e-Ener gy ’18: Proceedings of the Ninth International Conference on Futur e Energy Systems , (New Y ork, NY , USA), p. 366–374, Association for Computing Machinery , 2018. [53] Y . Zhang, J. W ang, and X. W ang, “Revie w on probabilis- tic forecasting of wind power generation, ” Renewable and Sustainable Energy Reviews , vol. 32, pp. 255–270, 2014. [54] M.-L. Kloubert, “Probabilistic load ﬂo w approach con- sidering dependencies of wind speed, solar irradiance, electrical load and energy exchange with a joint proba- bility distribution model, ” Ener gies , vol. 13, no. 7, 2020. [55] Z. W ang, S. Bu, J. W en, and C. Huang, “ A comprehen- siv e revie w on uncertainty modeling methods in mod- ern po wer systems, ” International J ournal of Electrical P ower & Energy Systems , vol. 166, p. 110534, 2025. [56] M. Ilic, R. O. Salcedo Ulerio, E. G. Corbett, E. Austin, M. Shatz, and E. R. Limpaecher , “ A Framework for Evaluating Electric Power Grid Improvements in Puerto Rico, ” tech. rep., MIT Lincoln Laboratory , 2019. [57] L. L. Anton, “A C Optimal Po wer Flow for Physically and Economically Informed Grid Decarbonization, ” thesis, Massachusetts Institute of T echnology , 2024. [58] SmartGridz, “SmartGridz User Manual, ” tech. rep., SmartGridz, 2025.

Unified Diagnostics for Quantifying AC Operating-Point Robustness Under Injection and Topological Uncertainties with Regime Changes

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment