Resilient Transmission Grid Design: AC Relaxation vs. DC approximation
Harsha Nagarajan, Russell Bent, Pascal Van Hentenryck, Scott Backhaus, Emre Yamangil
11 Resilient Transmission Grid Design: ACRelaxation vs. DC approximation
Harsha Nagarajan † , Russell Bent † , Pascal Van Hentenryck ‡ , Scott Backhaus † , Emre Yamangil †† Center for Nonlinear Studies, Los Alamos National Laboratory, NM, United States. ‡ Department of Industrial Operations Engineering, University of Michigan, Ann Arbor, United States.Contact: [email protected]
Abstract —As illustrated in recent years (SuperstormSandy, the Northeast Ice Storm of 1998, etc.), extremeweather events pose an enormous threat to the electric powertransmission systems and the associated socio-economic sys-tems that depend on reliable delivery of electric power.Besides inevitable malfunction of power grid components, de-liberate malicious attacks can cause high risks to the service.These threats motivate the need for approaches and methodsthat improve the resilience of power systems. In this paper,we develop a model and tractable methods for optimizingthe upgrade of transmission systems through a combinationof hardening existing components, adding redundant lines,switches, generators, and FACTS and phase-shifting devices.While many of these controllable components are includedin traditional design (expansion planning) problems, weuniquely assess their benefits from a resiliency point ofview. More importantly, perhaps, we evaluate the suitabilityof using state-of-the-art AC power flow relaxations versusthe common DC approximation in resilience improvementstudies. The resiliency model and algorithms are tested on amodified version of the RTS-96 (single area) system.
I. I
NTRODUCTION T HE modern electrical system is designed for trans-portation of large amounts of power from sources ofsupply to distant points of demand. Within these systems,the underlying high-voltage transmission networks playa vital role in achieving this mission. However, whentransmission networks are exposed to extreme event con-ditions, the ability to deliver power is degraded becauseof physical damage to overhead transmission lines andtowers. One example of such events are ice storms. Duringan ice storm, transmission towers can fail due to legbuckling and lines can fail due to the combined stressof ice accumulation and wind [1], [2], [3].When such events occur on large scales, outages andimpacts can be extreme. For example, in the winter of1998, an ice storm in northeastern North America toppledover 1000 transmission towers and 30,000 wooden utilitypoles. Over 5 million people were without power and theeconomic impacts were estimated at $2.6 billion [4]. Thus,given the potential social and economic impacts of theseevents, it is important to consider how to upgrade the de-sign of transmission systems to improve their performanceunder these conditions. In our preliminary work [5], we formulated the OptimalResilient Grid Design problem for Transmission systems(ORGDT) as a two-stage mixed-integer stochastic opti-mization problem and developed an algorithm to solvethis problem. The first stage selects from a set of potentialupgrades to the network. The second stage evaluates thenetwork performance benefit of the upgrades with a set ofdamage scenarios sampled from a stochastic distributionof events of concern. One of the important questions thatwas unanswered in this early work centered on the use ofrecently developed convex relaxations of AC power flowphysics [6] in resiliency planning. Here, we address thisconcern and convincingly show that using DC approxima-tion of power flow physics leads to severely suboptimaland infeasible solutions under certain operating conditionsof the grid. These results establish the necessity of usingmodels of AC physics in resilient design problems andanalysis.Like [5], we adopt the methods discussed in [3] inorder to sample realistic damage scenarios for transmis-sion systems. The ORGDT upgrade options include: a)Build new lines, b) Build switches, FACTS devices andPhase-shifting transformers (PST) to provide operationalflexibility, c) Harden existing lines to lower the proba-bility of damage, and d) Build new generation facilities.Minimal network (resiliency) performance is measured bysatisfying a minimum fraction of critical and non-criticalloads. We use the exact algorithms developed in [5] tocompare solutions found using the DC approximation andthe convex quadratic relaxations of [6].
Literature Review
Our recent paper [5] is the mostclosely related work. This paper develops the ORGDTmodel and algorithm used here and is based on theresilient distribution system design work developed in [7].Another important area of related work is interdictionmodeling and optimization. Here, the goal is to operateor design a system to make it as resilient as possibleto an adversary who can damage up to k elements [8],[9], [10], [11]. These models are a generalization of ourmodel when k is chosen to bound a worst-case event.However, given their min-max structure, these modelsare computationally challenging and are solvable only for a r X i v : . [ m a t h . O C ] M a r small k . Instead, we exploit the probabilistic nature of theadversary and we are able to address larger problems. Inaddition, existing interdiction models also do not generallyinclude AC physics. Closely related to interdiction andour model is reference [12]. They develop a model forhardening and upgrading transmission systems that isresilient to all possible physical damages bounded by acertain budget (i.e. weighted k ). Instead of developingan attacker-defender model as in interdiction modeling,they include each possible damage as a scenario (likeour model). Given the combinatorial number of possiblescenarios, they develop heuristics for reducing the numberof scenarios that are included. Unlike our approach, theyfocus only on the DC approximation and leave the de-velopment of decomposition algorithms for future work.A third area of related work is stochastic transmissionand generation expansion planning, where a recent surveydescribes some of the state-of-the-art [13].Overall, most papers in network design and expansionplanning use the linearized DC model and few studiesconsider FACTS devices and PSTs, although they mayhave significant benefits. Some notable exceptions includethe use of PSTs in network expansion [14], which uses agenetic algorithm over the DC model. See also the recentwork in [15] for the optimal placement of these devicesto avoid congestion.The key contributions of this paper include: • A detailed case study and extensions of the ORGDTmodel (from [5]) that compares DC approximationbased solutions and solutions based on state-of-the-art AC power flow relaxations. • An approach for recovering AC feasible solutionsfrom solutions found using the DC approximationand AC relaxations. • A detailed expansion model with appropriate gener-ator capacity modeling, FACTS devices and PST de-vices, and cutting-plane-based algorithmic enhance-ments that was omitted from [5].The rest of this paper is organized as follows. In sectionII we describe the ORGDT model. Section III describesour decomposition algorithms and Section IV our casestudies. We conclude with Section V.II. ORGDT O
PTIMIZATION M ODEL N OMENCLATURE
Parameters and notation N set of nodes (buses) E set of edges (lines and transformers) S set of disaster scenarios D s set of edges that are inoperable during s ∈ S i imaginary number constant | . | absolute value of the input argument c xij cost to build a line ( i , j ) if the line does not already exist c τij cost to build a switch on line ( i , j ) c tij cost to harden a line ( i , j ) c δij cost to build a FACTS device on line ( i , j ) c γij cost to build a PST device on line ( i , j ) c zpi cost of real generation capacity at bus ic ui cost to build a generation facility at bus iT ij apparent power thermal limit on line ( i , j ) L set of buses whose loads are critical G ij + i B ij admittance of line ( i , j ) ¯ G ij + i ¯ B ij modified admittance of line ( i , j ) due to top transformer R ij + i X ij impedance of line ( i , j ) R ij + i ¯ X ij modified impedance of line ( i , j ) due to top transformer θ u phase angle difference limit θ M big-M value on phase angle difference limit φ u phase shift limit v li , v ui lower and upper bound on voltage at bus i , respectively dp i + i dq i AC power demand at bus igp ui + i gq ui max. existing AC generation capacity at bus izp ui + i zq ui max. AC generation capacity that can be built at bus ilp cr + i lq cr fraction of critical AC power loads that must be served lp ncr + i lq ncr fraction of non-critical AC power loads that must beserved Binary variables x ij determines if line ( i , j ) is built τ ij determines if line ( i , j ) has a switch built t ij determines if line ( i , j ) is hardened u i determines the generation capacity built at bus iδ ij determines if FACTS device is built on line ( i , j ) γ ij determines if PST device is built on line ( i , j ) x sij determines if line ( i , j ) is used during disaster s ∈ S τ sij determines if switch ( i , j ) is used during disaster s ∈ S t sij determines if line ( i , j ) is hardened during disaster s ∈ S δ sij determines if FACTS device on line ( i , j ) is used duringdisaster s ∈ S γ sij determines if PST device on line ( i , j ) is used duringdisaster s ∈ S Continuous variables θ si phase angle at bus i during disaster s ∈ S φ sij phase angle shift introduced by PST on line ij duringdisaster s ∈ S v si voltage at bus i during disaster s ∈ S l sij current magnitude squared ( | I sij | ) on line ( i , j ) duringdisaster s ∈ S p sij + i q sij AC power flow on line ( i , j ) during disaster s ∈ S zp i + i zq i determines the capacity for AC generation at bus izp si + i zq si determines the capacity for AC generation at bus i duringdisaster s ∈ S gp si + i gq si AC power generated at bus i during disaster s ∈ S lp si + i lq si AC power load delivered at bus i during disaster s ∈ S yp si + i yq si determines the fraction of AC power load served at bus i during disaster s ∈ S We formulate the ORGDT as a two-stage mixed-integernonlinear program, P AC . The first-stage variables (with-out superscript s ) specify the new infrastructure enhance-ments and the second-stage variables (with superscript s )describe the operation of invested infrastructure to serveloads for each damage scenario s ∈ S . P AC := min (cid:88) ij ∈E (cid:16) c xij x ij + c τij τ ij + c tij t ij + c δij δ ij + c γij γ ij (cid:17) + (cid:88) i ∈N ( c ui u i + c zpi zp i ) (1a)s.t. x sij ≤ x ij ∀ ij ∈ D s , s ∈ S (1b) x sij = x ij ∀ ij / ∈ D s , s ∈ S (1c) τ sij ≤ τ ij , t sij ≤ t ij ∀ ij ∈ E , s ∈ S (1d) δ sij ≤ δ ij , γ sij ≤ γ ij ∀ ij ∈ E , s ∈ S (1e) zp si ≤ zp i , zq si ≤ zq i ∀ i ∈ N , s ∈ S (1f) ≤ zp i ≤ zp ui u i , | zq i | ≤ zq ui u i ∀ i ∈ N (1g) zp i ≥ | zq i | ∀ i ∈ N (1h) ( x s , τ s , t s , δ s , γ s , zp s , zq s , u ) ∈ Q AC ( s ) ∀ s ∈ S (1i) x ij , τ ij , t ij , δ ij , γ ij , u i ∈ { , } ∀ ij ∈ E , i ∈ N (1j) In P AC , Eq. (1a) minimizes the total upgrade cost, whichincludes new lines, switches, hardening, FACTS, phaseshifters, generators and real power capacity, respectively.Constraints (1b)-(1f) link the first-stage (construction)decisions with second-stage variables in Q AC ( s ) . Con-straints (1g) and (1h) model generation capacity upgrades.Without loss of generality, reactive power capacity islimited to half of the real power capacity. Constraint (1h)ensures that the generators are not built purely for reactivepower support. Constraint (1i) states that the mixed-integervector ( x s , τ s , t s , δ s , γ s , zp s , zq s , u ) ∈ Q AC ( s ) is anAC feasible transmission network for scenario s subjectto the constraints of Q AC ( s ) .The constraints of Q AC ( s ) describe the AC-power flowequations and budget constraints on resiliency options.Ohm’s law is modeled in constraints (2)-(5). For a giventopology of the network. ˜ x s , flow on line ij is forcedto zero when ˜ x sij = 0 . Kirchhoff’s current law (flowconservation) is given in constraints (12) and (13). Con-straints (9)-(11) model the connection between currentmagnitude ( l sij ) and the power loss equations as discussedin [16]. Operational constraints (14)-(17) represent thethermal limits, phase angle difference limits and voltagebounds, respectively. Constraint (19) models the damagedlines of the scenario s ∈ S , i.e., a line is inoperablewhen damaged and unhardened. Constraint (20) definesthe topology for scenario s . A switch is available foroperation only if the line is active. Constraints (22) boundthe real and reactive power capacities that can be built(if u i = 1 ) at bus i ∈ N . Constraints in (23) and (24)represent total real and reactive power capacities includingthe existing generators of the network. Finally, constraint(25) captures the continuous power shedding subject to theresiliency constraints in (26)-(29). User-defined resiliencyparameters ( lp cr , lq cr , lp ncr , lq ncr ) ensure that a minimumfraction of critical and non-critical loads is served duringevery damage scenario s ∈ S . Q AC ( s ) = { x s , τ s , t s , δ s , γ s , zp s , zq s , u : On/Off AC power flow equations: p sij = ˜ x sij ( ˜ G ij v s i − ˜ G ij v si v sj cos( θ sij ) − ˜ B ij v si v sj sin( θ sij )) ∀ ij ∈ E , (2) q sij = ˜ x sij ( − ˜ B ij v s i + ˜ B ij v si v sj cos( θ sij ) − ˜ G ij v si v sj sin( θ sij )) ∀ ij ∈ E , (3) p sji = ˜ x sij ( ˜ G ij v s j − ˜ G ij v si v sj cos( θ sij )+ ˜ B ij v si v sj sin( θ sij )) ∀ ij ∈ E , (4) q sji = ˜ x sij ( − ˜ B ij v s j + ˜ B ij v si v sj cos( θ sij ) + ˜ G ij v si v sj sin( θ sij )) ∀ ij ∈ E , (5) θ sij = θ si − θ sj + φ sij ∀ ij ∈ E , (6) ˜ G ij = ( ¯ G ij − G ij ) δ sij + G ij ∀ ij ∈ E , (7) ˜ B ij = ( ¯ B ij − B ij ) δ sij + B ij ∀ ij ∈ E , (8) p sij + p sji = R ij l sij ∀ ij ∈ E , (9) q sij + q sji = X ij l sij ∀ ij ∈ E , (10) p s ij + q s ij = l sij v s i ∀ ij ∈ E , (11) gp si − lp si = (cid:88) ij ∈E p sij + (cid:88) ji ∈E p sij ∀ i ∈ N , (12) gq si − lq si = (cid:88) ij ∈E q sij + (cid:88) ji ∈E q sij ∀ i ∈ N , (13) Operational limits and topology constraints: p sij + q sij ≤ ˜ x sij T ij ∀ ij ∈ E , (14) p sji + q sji ≤ ˜ x sij T ij ∀ ij ∈ E , (15) | θ sij | ≤ ˜ x sij θ u + (1 − ˜ x sij ) θ M ∀ ij ∈ E , (16) v li ≤ v si ≤ v u ∀ i ∈ N , (17) − φ u γ sij ≤ φ sij ≤ φ u γ sij ∀ ij ∈ E , (18) x sij = t sij ∀ ij ∈ D s (19) ˜ x sij = x sij − τ sij ≥ ∀ ij ∈ E (20) δ sij ≤ ˜ x sij , γ sij ≤ ˜ x sij ∀ ij ∈ E , (21) Generation constraints ∀ i ∈ N :0 ≤ zp si ≤ zp ui u i , | zq si | ≤ zq ui u i , zp si ≥ | zq si | (22) ≤ gp si ≤ gp ui + zp si , (23) ( gq li − | zq si | ) ≤ gq si ≤ ( gq ui + | zq si | ) , (24) Resilience constraints: lp si = yp si dp i , lq si = yq si dq i ∀ i ∈ N , (25) (cid:88) i ∈L lp si ≥ lp cr (cid:88) i ∈L dp i , (26) (cid:88) i ∈N\L lp si ≥ lp ncr (cid:88) i ∈N\L dp i (27) (cid:88) i ∈L lq si ≥ lq cr (cid:88) i ∈L dq i (28) (cid:88) i ∈N\L lq si ≥ lq ncr (cid:88) i ∈N\L dq i (29) x s , τ s , t s , δ s , γ s ∈ { , } ; 0 ≤ yp s , yq s ≤ } A. Additional Investment Options
In this paper, we also include FACTS and PST devicesand as investment options. These devices are often usefulfor addressing congestion in overloaded transmission sys-tems, a situation that can occur during major disruptions.In addition, these devices are cost-effective [17] andmay reduce resiliency costs significantly by replacing theneed for new transmission lines or hardening the existingdamaged lines.
Series Compensators (FACTS)
We model series com-pensation with reactance reduction. As described in thenomenclature, δ sij indicates if a compensation device isused on line ij during scenario s ∈ S . To model thisdisjunction we introduce ( ˜ G ij , ˜ B ij ) , a tuple representingthe following conditional expression: ( ˜ G ij , ˜ B ij ) = (cid:40) ( ¯ G ij , ¯ B ij ) , if used ( δ sij = 1)( G ij , B ij ) , otherwise ( δ sij = 0) This expression is modeled as constraints (7) and (8).Further, this disjunction is captured in the modified powerflow equations (2)-(5). The non-convexity of this disjunc-tion is discussed later in this section.
Phase-Shifting Transformers (PSTs)
PSTs are devicesthat adjust voltage phase angles. As described in thenomenclature, γ sij indicates if a phase shifter is used online ij during scenario s ∈ S . This disjunction is describedin the following conditional expression: θ sij = θ si − θ sj + φ sij , − φ u ≤ φ sij ≤ φ u . if installed ( γ sij = 1) ,θ si − θ sj , otherwise ( γ sij = 0) . This expression is modeled through constraints (6) and(18).
B. Convex Relaxations for On/Off AC Power Flow
Numerous convex relaxations of AC power flow equa-tions have been proposed in the literature to obtain tightlower bounds on the original nonconvex formulation.These relaxations have various trade-offs in complexity,quality, and performance. However, very few relaxationmethods can be applied in the context of transmissionswitching problems as the addition of binary variablesincreases the complexity of the problems significantly.Hence, given the mixed-integer nature of ORGDT, whichleads to an even more complex formulation with variousresiliency options, we extend the recently developed con-vex quadratic relaxations of the power flow equations toORGDT [6]. For completeness, we discuss the necessaryconvex relaxations for multilinear and trigonometric ex-pressions which appear in nonconvex constraints (2)-(5)and (11).
Multilinear expressions
Given any two variables x i , x j ∈ [ x i x i ] × (cid:2) x j x j (cid:3) , the McCormick relaxation is used tolinearize the bilinear product x i x j by introducing a newvariable (cid:99) x ij ∈ (cid:104) x i , x j (cid:105) MC . The feasible region of thisvariable is defined by equations (30). This relaxation isexact if one variable is binary [18]. (cid:99) x ij ≥ x i x j + x j x i − x i x j (30a) (cid:99) x ij ≥ x i x j + x j x i − x i x j (30b) (cid:99) x ij ≤ x i x j + x j x i − x i x j (30c) (cid:99) x ij ≤ x i x j + x j x i − x i x j (30d) Multilinear expressions of the form x i x j . . . x p are re-laxed using a sequential bilinear approach as follows: (cid:104)(cid:104) x i , x j (cid:105) MC , . . . , x p (cid:105) MC . Quadratic terms
Given a variable x i ∈ [ x i x i ] , a second-order conic relaxation is applied to x i by introducing anew variable (cid:98) x i ∈ (cid:104) x i (cid:105) MC − q , such that, (cid:98) x i ≥ x i (31a) (cid:98) x i ≤ ( x i + x i ) x i − x i x i (31b) By applying the above approaches and disjunctive con-vex relaxations for the trigonometric functions with on/offvariables [6], we outer-approximate the non-convex power flow constraints. For brevity, we focus on constraint (2)but similar relaxations hold for constraints (3)-(5). p sij =( G ij − G ij ) (cid:16)(cid:100) δxv sij − (cid:100) δwc sij (cid:17) + ( B ij − B ij ) (cid:100) δws sij + G ij (cid:0)(cid:99) xv sij − (cid:99) wc sij (cid:1) − B ij (cid:99) ws sij ∀ ij ∈ E where (cid:99) xv sij , (cid:99) wc sij , (cid:99) ws sij , (cid:100) δxv sij , (cid:100) δwc sij and (cid:100) δws sij satisfythe following constraints: (cid:98) v si ∈ (cid:104) v si (cid:105) MC − q , (cid:99) xv sij ∈ (cid:104) ˜ x sij , (cid:98) v si (cid:105) MC (32a) (cid:98) w sij ∈ (cid:104) v si , v sj (cid:105) MC (32b) (cid:99) wc sij ∈ (cid:104) (cid:98) w sij , (cid:98) cs sij (cid:105) MC , (cid:99) ws sij ∈ (cid:104) (cid:98) w sij , (cid:99) sn sij (cid:105) MC (32c) (cid:100) δwc sij ∈ (cid:104) δ sij , (cid:99) wc sij (cid:105) MC , (cid:100) δws sij ∈ (cid:104) δ sij , (cid:99) ws sij (cid:105) MC (32d) (cid:98) cs sij ≤ ˜ x sij − (cid:18) − cos( θ u )( θ u ) (cid:19) ( θ sij + (˜ x sij − θ M ) ) (32e) ˜ x sij cos( θ u ) ≤ (cid:98) cs sij ≤ ˜ x sij (32f) (cid:99) sn sij ≤ cos( θ u / θ sij + ˜ x sij (sin( θ u / − θ u / θ u / (32g) + (1 − ˜ x sij )(cos( θ u / θ M + 1) (cid:99) sn sij ≥ cos( θ u / θ sij − ˜ x sij (sin( θ u / − θ u / θ u / (32h) − (1 − ˜ x sij )(cos( θ u / θ M + 1)˜ x sij sin( − θ u ) ≤ (cid:99) sn sij ≤ ˜ x sij sin( θ u ) (32i) Finally, we relax constraint (11) into a second-order conicconstraint of the form p s ij + q s ij ≤ l sij ˆ v si . C. DC approximation
Within the expansion planning literature, the DC ap-proximation is often used for designing power systems.For comparison purposes, we develop a version of theORGDT problem based on the DC model in (33). Q DC ( s ) = { x s , τ s , t s , δ s , γ s , zp s , u : On/Off DC power flow equations: p sij ≤ − ˜ B ij (cid:16) θ sij + θ M (1 − ˜ x sij ) (cid:17) ∀ ij ∈ E , (33a) p sij ≥ − ˜ B ij (cid:16) θ sij − θ M (1 − ˜ x sij ) (cid:17) ∀ ij ∈ E , (33b) ˜ G ij = ( ¯ G ij − G ij ) δ sij + G ij ∀ ij ∈ E , (33c) ˜ B ij = ( ¯ B ij − B ij ) δ sij + B ij ∀ ij ∈ E , (33d) p sij + p sji = 0 , θ sij = θ si − θ sj + φ sij ∀ ij ∈ E , (33e) gp si − lp si = (cid:88) ij ∈E p sij ∀ i ∈ N , (33f) Operational limits and topology constraints: | p sij | ≤ ˜ x sij T ij ∀ ij ∈ E , (33g) | θ sij | ≤ ˜ x sij θ u + (1 − ˜ x sij ) θ M ∀ ij ∈ E , (33h) − φ u γ sij ≤ φ sij ≤ φ u γ sij ∀ ij ∈ E , (33i) x sij = t sij ∀ ij ∈ D s (33j) ˜ x sij = x sij − τ sij ≥ ∀ ij ∈ E (33k) Generation and demand constraints ∀ i ∈ N :0 ≤ zp si ≤ zp ui u i , ≤ gp si ≤ gp ui + zp si , (33l) (cid:88) i ∈L lp si ≥ lp cr (cid:88) i ∈L dp i , (33m) (cid:88) i ∈N\L lp si ≥ lp ncr (cid:88) i ∈N\L dp i (33n) lp si = yp si dp i ; 0 ≤ yp si ≤ , x s , τ s , t s ∈ { , }} D. AC-feasibility formulation
It is important to note that the convex relaxationsintroduced in section II-B and the DC approximationof (33) might violate the nonconvex AC-equations de-scribed in Eqs. (2)-(5), (11). In this section, we developa model for recovering primal feasible solutions to assesthe quality of these approximations. This model fixes allthe build decisions from the second stage solutions forevery disaster s ∈ S and solves for the feasibility of thenonconvex equations based on a primal-dual interior pointmethod [19]. Formally, let ( ˜x s , t s , δ s , γ s , zp s , zq s , u ) bethe build decisions obtained in the first stage. Under theDC model, at bus i , we assume that zq si is equal to halfof the real power capacity, that is | zq si | = zp si / . Wethen minimize a weighted sum of the violations to theresilience constraints which serves as a proxy to measurethe distance to AC feasibility. p := min M ( λ pcr + λ qcr ) + λ pncr + λ qncr s.t. Given ( ˜x s , t s , δ s , γ s , zp s , zq s , u )Constraints (2) − (6) , (9) − (17) , (22) − (25) ,λ pcr ≥ (cid:32) lp cr (cid:88) i ∈L dp i (cid:33) − (cid:88) i ∈L lp si , (34a) λ pncr ≥ lp ncr (cid:88) i ∈N\L dp i − (cid:88) i ∈N\L lp si (34b) λ qcr ≥ (cid:32) lq cr (cid:88) i ∈L dq i (cid:33) − (cid:88) i ∈L lq si (34c) λ qncr ≥ lq ncr (cid:88) i ∈N\L dq i − (cid:88) i ∈N\L lq si (34d) λ pcr , λ qcr , λ pncr , λ qncr ≥ (34e) where M is used to weight satisfaction of critical load andand λ pcr , λ qcr , λ pncr , λ qncr represent the load shed variableson critical and non-critical load constraints (up to theresilience criteria) as shown in constraints (34a)-(34d).III. A LGORITHMS
In this section we discuss algorithms we use to solvethe ORGDT. The ORGDT is a two-stage Mixed-IntegerQuadratically Constrained Program (MIQCP) with a blockdiagonal structure. In order to exploit this structure, wegeneralize the scenario-based decomposition (SBD) tech-niques of [7] to solve the ORGDT. In the remainder of thispaper, let P AC ( S (cid:48) ) denote the ORGDT with the scenarioset S (cid:48) ⊆ S and σ denote the vector of constructionvariables x ij , τ ij , t ij , γ ij , δ ij for all ij ∈ E and u i , zp i for all i ∈ N . SBD can be applied to ORGDT after thefollowing key observation: Observation III.1.
The second stage variables do notappear in the objective function. Therefore any optimalfirst stage solution based on a subset of the second stagesubproblems that is feasible for the remaining scenarios,is an optimal solution for the original problem.
The SBD approach is outlined in Algorithm 1. At a highlevel, Algorithm 1 solves problems with iteratively largersets of scenarios until a solution is obtained that is feasiblefor all scenarios. The algorithm takes as input the set ofscenarios and an initial scenario to consider, S (cid:48) . Line 2solves the ORGDT on S (cid:48) , where P AC ( S (cid:48) ) and σ ∗ are usedto denote the problem and solution respectively. Line 3then evaluates σ ∗ on the remaining scenarios in S \S (cid:48) . Thefunction l : P (cid:48) ( s, σ ∗ ) → R + , is an infeasibility measurethat is 0 if the problem is feasible for scenario s , positiveotherwise. This function is implemented by solving thefeasibility subproblem reliability constraints, i.e. total andcritical demands are satisfied. We use a version of (34)where M = 1 and the non convex constraints are replacedwith their convex counterparts. This function prices thecurrent solution over s ∈ S \ S (cid:48) . If all prices are 0,then the algorithm terminates with solution σ ∗ (lines 4-5).Otherwise, the algorithm adds the scenario with the worstinfeasibility measure to S (cid:48) (line 7). Algorithm 1:
Scenario-Based Decomposition input:
A set of disasters S and let S (cid:48) = S ; while S \ S (cid:48) (cid:54) = ∅ do σ ∗ ← Solve P AC ( S (cid:48) ) ; I ← (cid:10) s , s . . . s |S\S (cid:48) | (cid:11) s ∈ S \ S (cid:48) : l ( P AC (cid:48) ( s i , σ ∗ )) ≥ l ( P AC (cid:48) ( s i +1 , σ ∗ )) ; if l ( P AC (cid:48) ( I (0) , σ ∗ )) ≤ then return σ ∗ ; else S (cid:48) ← S (cid:48) ∪ I (0) ; return σ ∗ Remark III.2.
We observed that the LP-relaxation forthe ORGDT is very loose. To overcome this issue, weaugmented every iteration of Algorithm 1 with the previousoptimum objective value as a lower bound for the currentiteration.
Cutting-plane algorithm for quadratic constraints:
Even though optimization theory guarantees that the set ofconvex inequalities in ORGDT can be solved efficiently,several numerical experiments demonstrated that it waschallenging to solve even moderately sized problems usingstate-of-the-art MIQCP solvers (CPLEX/Gurobi). Eitherthe solver convergence was very slow or it terminated with“numerical trouble”. In order to circumvent this issue, weadopted a cutting-plane approach. Let P miq be a generalMIQCP which contains the following quadratic constraint(rotated second-order cone): x T x ≤ yz, y ≥ , z ≥ , ( x , y, z ) ∈ R | x | +2 (35)Algorithm 2 outlines the cutting plane procedure tosolve P miq . The key idea behind this algorithm is to solvea finite sequence of MILPs to obtain an optimal solutionfor the original MIQCP. Let P mil , a MILP, represent the original problem without the quadratic constraints. Thesolution in step 3 of Algorithm 2 is then guaranteed tobe a lower bound to P miq . This lower bound is tightenedfurther for every infeasible solution (step 4) by adding avalid cutting plane (step 5). This cutting plane is valid asit is an outer-approximation of the original feasible convexset. This procedure is repeated until the solution obtainedis feasible (satisfies step 4), and hence optimal, for P miq . Algorithm 2:
Cutting plane algorithm for MIQCPs Notation : f ( x , y ) = x T x y Input : P mil , tol > , L = ∅ ˆ x , ˆ y, ˆ z ← Solve P mil ( x , y, z ) while ˆ x T ˆ x > ˆ y ˆ z + tol do Augment L with the following cut: f (ˆ x , ˆ y )+ (cid:80) | x | i =1 ∂f ( x ,y ) ∂x i ( x i − ˆ x i )+ ∂f ( x ,y ) ∂y ( y − ˆ y ) ≤ z ˆ x , ˆ y, ˆ z ← Solve (cid:0) P mil ( x , y, z ) ∪ L (cid:1) IV. C
ASE S TUDY
In this section, we discuss three numerical studies forthe ORGDT based on the modified single area IEEE RTS-96 system [20]. We first discuss the computational benefitsassociated with using the SBD algorithm. We secondcompare solutions based on the DC approximation withsolutions based on the QC relaxation. Finally, we examinethe benefits of FACTS and PST devices. All modeling wasdone using JuMP [21] on a computer with 32 threads,a 2.6GHz Intel 64 bit processor, 25.6MB L3 cache and64GB of memory. Gurobi 6.5.2 was used for solving theMIQCPs and MILPs (optimality gap 0.1%) and Knitro9.1.0 was used for solving the nonlinear programs.
Test system
We use a modified IEEE single-area RTS-96 system that has 24 buses including 17 load buses, 38transmission lines and 32 conventional generators [20].The total installed capacity of the existing generators is3405 MW. The total load in the system is 2850 MW.We labeled 1740 MW as critical loads. The network wasspatially placed in an area of 4250 miles . Since this testsystem does not have the co-ordinate data for every bus,we approximately evaluated them based on the lengths oflines [20]. The admittance, impedance and apparent powerthermal limit values on lines are from the standard testcase. We assume that FACTS, PSTs and new generationcapacities, if chosen, can be built anywhere in the network.The remaining parameters are outlined in Table I. Scenario generation
The damage scenarios are basedon line failure probabilities that follow a multivariateGaussian distribution with a mean placed at the center ofthe network. A Bernoulli trial was applied to every linefor chosen percentiles (% damage) to generate the randomscenarios. Empirically, we observed that 20 scenarios weresufficient to represent the features of the distribution.
TABLE IP
ARAMETERS FOR THE TEST SYSTEM c xij , c tij , c τij $1.35m/mile, $5000/mile, $1000 c δij , c γij $50000, $50000 c ui , c zpi $0.1m, $0.817m/MW v li , v ui , φ u , ¯ X ij ◦ , . X ij lp cr , lq cr , lp ncr , lq ncr A. Computational Performance of SBD
Table II compares the computational time of SBD withthe computational time of solving the full model using aMIQCP (Gurobi) solver. Here, we focus on the problemwhere only new lines, hardened lines, and distributedgenerations are available as design choices. On average,SBD is 17x faster when θ u = 15 ◦ and 26x faster when θ u = 45 ◦ . We highlight in the 80% damage case thatafter 20 hours the MIQCP solver could only provide a 7%optimality gap, whereas SBD terminated with an optimalsolution after 12 minutes. It is important to note that whenFACTS and PST devices are included as design options,the MIQCP solver fails to find any feasible solution after10 hours of computation. TABLE IIC
OMPARISON OF WALL TIME ( SEC .) OF SBD
AND
MIQCP
SOLVERSON PROBLEMS WITHOUT
FACTS
AND
PST
DEVICES . *
INDICATESTHAT THE BEST GAP WAS
HE FIRST COLUMN INDICATES THELEVEL OF DAMAGE . T
HE SECOND / FIFTH AND FOURTH / SEVENTHCOLUMNS INDICATE THE WALL TIMES OF
SBD
AND
MIQCP,
RESPECTIVELY . T
HE THIRD AND SIXTH COLUMNS SHOW THENUMBER OF ITERATIONS FOR
SBD θ u = 15 ◦ θ u = 45 ◦ Damage SBD Iters Full SBD Iters Full90% 634.8 3 5008.2 2419.6 4 13729.280% 303.7 5 1508.7 730.1 5 72886.3 ∗
70% 417.7 4 6566.9 485.1 4 2339.060% 188.3 2 8853.8 157.6 3 3537.550% 90.2 2 6156.7 53.1 1 9372.8
B. AC vs. DC power flow models
In this section we scale the loads of the model between0.8 and 1.25 at .05 intervals. We define the gap betweenthe solutions based on the QC relaxation and solutionsbased on the DC approximation as ζ = (cid:32) opt( P AC ) − opt( P DC )opt( P AC ) (cid:33) ∗ where opt( P ACo ) and opt( P DCo ) are optimal upgradecosts for QC relaxation and DC approximation models,respectively. As shown in table III, DC-approximation-based solutions are always lower bounds to solutionsbased on QC relaxations. In the case of small phase angledifferences ( θ u = 15 ), the maximum and average ζ is atthe maximum damage case (90%). When the phase angledifferences are increased, the value of ζ increases dramat-ically. This behavior in the solutions is not unexpected. DC approximations tend to perform poorly when there arelarge phase angle differences as the assumptions behindthe approximation no longer hold. In Figure 1, we expandthe results of the 90% damage row of Table III and seehow the DC approximation consistently under estimatesthe design costs.
TABLE IIIE
ACH ROW DESCRIBES THE AVERAGE , MIN , AND MAX ζ ACROSSLOAD SCALING VALUES BETWEEN
AND
X FOR A SPECIFIEDDAMAGE LEVEL . θ u = 15 ◦ θ u = 45 ◦ Dam. avg min max avg min maxgap (%) gap (%) gap (%) gap (%) gap (%) gap (%)90% 9.7 0 35.4 22.4 5.5 47.380% 9.2 0 25.3 22.1 11.7 37.670% 7.0 0 29.3 24.0 15.4 55.760% 5.6 0 29.4 24.9 7.4 43.150% 2.1 0 11.3 35.9 0 53.1 . . . . . . . . . · Load scale T o t a l up g r a d ec o s t( d o ll a r s ) QC relaxationDC approximation (a) 90% damage, θ u = 15 o . . . . . . . . · Load scale T o t a l up g r a d ec o s t( d o ll a r s ) QC relaxationDC approximation (b) 90% damage, θ u = 45 o Fig. 1. Comparison of total upgrade cost for the QC relaxation with theDC approximation
AC feasibility analysis ζ is not the only metric forcomparing the solutions based on the DC and QC models.It is also important to understand distance to AC feasibility(section II-D) in each scenario. For a given scenario s , wecalculate the percent apparent-load shed with: µ = (cid:114) ( λpcr )2 + ( λqcr )2 (cid:115)(cid:18)(cid:80) i ∈L dpi (cid:19) (cid:18)(cid:80) i ∈L dqi (cid:19) , (cid:114) ( λpncr )2 + ( λqncr )2 (cid:118)(cid:117)(cid:117)(cid:116)(cid:32)(cid:80) i ∈N\L dpi (cid:33) (cid:32)(cid:80) i ∈N\L dqi (cid:33) × for critical and noncritical loads, respectively.Table IV summarizes the average and maximum µ for each scenario for the solution found with the DCapproximation for a load scale of 1.0. The QC relaxation istight and no (additional) load shedding is required. Thoughthe DC approximation is considered a good approximationwhen phase angle differences are small, we still observeload sheds of up to 3.7% of critical loads and 6.7% ofnon-critical loads in the maximum damage case. Again,the DC approximation’s performance degrades further forlarger phase angle differences with a maximum µ of 4.5%(critical) and 9.4% (non-critical). Figure 2 further expandson µ for each scenario in the 90% damage case.To summarize, OGRDT solutions based on the DCapproximation perform reasonable well only on systemswith small phase angle differences and limited damage.When the grid is damaged heavily, DC approximationperforms poorly, even at smaller θ u values. At larger θ u values, the DC approximation drastically under approxi-mates the required upgrades, and results in solutions thatseverely violate resilience criteria based on load shedding.However, these issues can be addressed by applying ACrelaxations in lue of DC approximations. TABLE IVA
VERAGE AND MAXIMUM µ ( CRITICAL , NON - CRITICAL ) FOR DC SOLUTIONS AND A LOAD SCALE OF θ u = 15 ◦ θ u = 45 ◦ Damage avg. µ max. µ avg. µ max µ
90% 0.3, 0.7 3.7, 6.7 0.5, 1.6 4.5, 8.980% 0.1, 1.2 0.3, 7.3 0.2, 3.0 1.7, 9.270% 0.1, 0.1 0.3, 2.7 0.15, 2.5 0.9, 9.460% 0.1, 0.1 0.3, 0.2 0.13, 2.4 0.3, 9.050% 0, 0 0, 0.03 0, 1.8 0, 7.5 P e r ce n t o f a pp a r e n t - l oa d s h e d ( µ ) (a) 90% damage, θ u = 15 o P e r ce n t o f a pp a r e n t - l oa d s h e d ( µ ) non-critical (QC)critical (QC)non-critical (DC)critical (DC) (b) 90% damage, θ u = 45 o Fig. 2. Values of µ for the QC relaxation and DC approximation for90% damage scenarios when the load scale is 1.0. C. Benefits of FACTS and PST devices
In this section, we summarize the benefits of includ-ing FACTS and PST device options in the ORGDT. Toquantify the benefits, we define the following: ψ = (cid:32) opt( P AC ) − opt( P AC − dev )opt( P AC ) (cid:33) × FACTS devices
Table V summarizes the benefits of in-cluding FACTS devices to achieve advantages of operatingtransmission grids at larger θ u values without the stabilityissues. Overall, we observe that savings is a much as17.2% when θ u = 15 ◦ . The benefit of FACTS devices issmaller when θ u is larger, since the operating points areless congested and the need for such devices is reduced.We use Figure 3 to expand on ψ for the 90% damagecase. TABLE VA
VERAGE AND MAXIMUM ( PARENTHESIS INDICATE LOAD SCALEVALUE FOR THE MAX ) ψ ACROSS LOAD SCALES BETWEEN
AND
WHEN
FACTS
DEVICES ARE INCLUDED IN THE
OGRDT. θ u = 15 ◦ θ u = 45 ◦ Damage avg. ψ max. ψ avg. ψ max. ψ
90% 9.3 31.6 (1.25) 5.6 20.7 (1.25)80% 9.6 31.3 (1.2) 5.5 20.2 (1.25)70% 10.5 28.0 (1.2) 4.5 14.6 (1.25)60% 10.3 29.7 (1.2) 4.6 28.0 (1.25)50% 17.2 33.2 (1.2) 11.8 29.4 (1.25)
In Figure 3, it is interesting to note the (sometimes) non-increasing trend in ψ as load scale increases. Intuitively, Load scale T o t a l up g r a d ec o s t( $ ) QC objective (without FACTS) QC objective (with FACTS)
Load scale .10 ψ (%) ψ (a) 90% damage, θ u = 15 ◦ (b) 90% damage, θ u = 45 ◦ Fig. 3. Comparison of total upgrade costs and ψ for optimal upgradesolutions with and without FACTS devices. this indicates that upgrades required for certain load scalesare also sufficient for further increases in load. More im-portantly, perhaps, the OGRDT solution with and withoutFACTS devices is sometimes identical. This phenomenaindicates that at certain load scales, FACTS devices cannoteliminate the need for other upgrades (line hardening, newlines, etc.) and once those upgrades are installed, the needfor FACTS devices is eliminated. Of course, as the loadscale increases, FACTS devices are once again required. PST devices
Generally speaking, we only found PSTdevices to be of benefit when the network is congested andhas larger load scales. Under the 90% damage cases, θ u =15 ◦ and load scales greater than 1.20, we found solutionswith savings of around 10%. Though the improvement in ψ when PSTs are included is not as impressive as whenFACTS devices are included, we believe that the savingsare non trival when considering large-scale grids.V. C ONCLUSIONS
In this paper, we formulated, modeled, and developed amodel of ORGDT. Our contributions include a comparisonof solutions to the ORGDT with the DC approximationand the QC relaxation that demonstrates the necessity ofAC power flow modeling in resilience studies. We alsodeveloped an approach for recovering fully AC feasiblesolutions that indicate that solutions based on the QCrelaxation tend to be tight in practice. Finally, we havealso shown that a small set of strategically located PSTand FACTS devices as design options can reduce theneed for other design considerations like line hardeningor distributed generation.There remain a number of interesting future directionsfor work in this area. First, it will be important to scalethe approaches of the paper to larger, more realistictransmission grids. An important idea in this area is tolimit the upgrades to those that are deemed practical bysubject matter experts. Second, it will be important tointroduce tighter and improved models of FACTS and PSTdevices (such as continuous set points). Finally, it willbe interesting to consider= approaches for using FACTSand PST devices to improve the restoration process oftransmission grids, as seen in [22]. R
EFERENCES[1] E. Brostrom, “Ice Storm Modelling in Transmission System Relia-bility Calculations,” Ph.D. dissertation, Royal Institute of Technol-ogy, 2007.[2] J. Eidinger and L. Kempner, “Reliability of Transmission Towersunder Extreme Wind and Ice Loading,”
G&E Engineering SystemsInc. and Bonneville Power Administration , 2013.[3] S. Wang, “Influence of ice accretions on mechanical and electricalperformance of overhead transmission lines,”
Journal of Electricaland Mechanical Engineering , vol. 1, no. 1, 2010.[4] E. Mills, “Electric Grid Disruptions and Extreme Weather,” in
US Disaster Reanalysis Workshop National Climatic Data Center ,Asheville, NC, 2012.[5] H. Nagarajan, E. Yamangil, R. Bent, P. V. Hentenryck, and S. Back-haus, “Optimal resilient transmission grid design,” in , June 2016, pp. 1–7.[6] H. Hijazi, C. Coffrin, and P. Van Hentenryck, “Convex quadratic re-laxations for mixed-integer nonlinear programs in power systems,”
Mathematical Programming Computation , pp. 1–47, 2014.[7] E. Yamangil, R. Bent, and S. Backhaus, “Resilient upgrade ofelectrical distribution grids,” in
Twenty-Ninth AAAI Conference onArtificial Intelligence , 2015.[8] R. L.-Y. Chen and C. A. Phillips, “k-edge failure resilient networkdesign,”
Electronic Notes in Discrete Mathematics , vol. 41, no. 0,pp. 375 – 382, 2013.[9] R. L.-Y. Chen, A. Cohn, N. Fan, and A. Pinar, “Contingency-RiskInformed Power System Design,”
IEEE Transactions on PowerSystems , vol. 29, no. 5, pp. 2087–2096, Sep. 2014.[10] J. Salmeron, K. Wood, and R. Baldick, “Worst-case interdictionanalysis of large-scale electric power grids,”
IEEE Transactions onPower Systems , vol. 24, no. 1, pp. 96–104, 2009.[11] A. Delgadillo, J. Arroyo, and N. Alguacil, “Analysis of ElectricGrid Interdiction with Line Switching,”
IEEE Transactions onPower Systems , vol. 25, no. 2, pp. 633–641, 2010.[12] N. Nezamoddini, S. Mousavian, and M. Erol-Kantarci, “A riskoptimization model for enhanced power grid resilience againstphysical attacks,”
Electric Power Systems Research , vol. 143, pp.329–338, 2017.[13] V. Krishnan, J. Ho, B. Hobbs, A. Liu, J. McCalley, M. Shahideh-pour, and Q. Zheng, “Co-Optimization of Electricity Transmissionand Generation Resources for Planning and Policy Analysis: Re-view of Concepts and Modeling Approaches,”
Energy Systems ,2015.[14] C. T. Miasaki, E. Franco, and R. A. Romero, “Transmission net-work expansion planning considering phase-shifter transformers,”
Journal of Electrical and Computer Engineering , vol. 2012, 2012.[15] V. Frolov, S. Backhaus, and M. Chertkov, “Efficient algorithm forlocating and sizing series compensation devices in large powertransmission grids: Ii. solutions and applications,”
New Journal ofPhysics , vol. 16, no. 10, p. 105016, 2014.[16] M. Baran and F. F. Wu, “Optimal sizing of capacitors placed on aradial distribution system,”
IEEE Trans. on power Delivery , vol. 4,no. 1, pp. 735–743, 1989.[17] R. Baldick and R. O’Neill, “Estimates of comparative costs foruprating transmission capacity,”
Power Delivery, IEEE Trans. on ,vol. 24, no. 2, pp. 961–969, April 2009.[18] H. Nagarajan, M. Lu, E. Yamangil, and R. Bent, “Tighteningmccormick relaxations for nonlinear programs via dynamic mul-tivariate partitioning,” in
International Conference on Principlesand Practice of Constraint Programming . Springer, 2016, pp.369–387.[19] S. J. Wright,
Primal-dual interior-point methods . SIAM, 1997.[20] P. Wong, P. Albrecht, R. Allan, R. Billinton, Q. Chen, C. Fong,S. Haddad, W. Li, R. Mukerji, D. Patton et al. , “The ieee reliabilitytest system-1996. a report prepared by the reliability test systemtask force of the application of probability methods subcommittee,”
Power Systems, IEEE Trans. on , vol. 14, pp. 1010–1020, 1999.[21] I. Dunning, J. Huchette, and M. Lubin, “Jump: A modelinglanguage for mathematical optimization,” arXiv:1508.01982 , 2015.[22] T. Mak, C. Coffrin, P. Van Hentenryck, I. Hiskens, and D. Hill,“Power System Restoration Planning with Standing Phase Angleand Voltage Difference Constraints,” in