An Adversarial Model for Attack Vector Vulnerability Analysis on Power and Gas Delivery Operations

An Adv ersarial Mo del for A ttac k V ector V ulnerabilit y Analysis on P o w er and Gas Deliv ery Op erations Ignacio Losada Carre ˜ no † , Anna Scaglione † , Anatoly Zlotnik ‡ , Deep jyoti Dek a ‡ and Kaarthik Sundar ‡ † Arizona State Univ ersity , T emp e, AZ, USA, ‡ Los Alamos National Lab oratory , Los Alamos, NM, USA Octob er 10, 2019 Abstract P ow er systems often rely on natural gas pip eline net works to supply fuel for gas-fired gen- eration. Mark et inefficiencies and a lack of formal co ordination b et ween the wholesale p ow er and gas deliv ery infrastructures may magnify the broader impact of a cyb er-attac k on a natural gas pip eline. In this study w e present a mo del that can b e used to quantify the impact of cyb er-attac ks on electricit y and gas delivery op erations. W e model activ ation of cyb er-attac k v ectors that attempt to gain access to pip eline gas compressor con trols using a con tinuous-time Mark ov chain ov er a state space based on the gas op erator Industrial Control System firew all zone partition. Our approach ev aluates the op erating states and decision-making in the net- w orks using physically realistic and op erationally represen tative mo dels. W e summarize these mo dels, the sequence of analyses used to quan tify the impacts of a cyber-incident, and prop ose a Mon te Carlo sim ulation approac h to quan tify the resulting effect on the reliability of the bulk p o w er system by the increase in operational cost. The metho dology is applied to a case study of interacting p o wer, gas, and cyber test netw orks. Keyw ords: Con tin uous-time Marko v c hain, cyb er-ph ysical Systems, optimal p o wer flow, Purdue mo del, transient gas flow 1 In tro duction Reliance on natural gas for electricity generation has increased significantly in many regions of the w orld [1]. Low prev ailing gas prices [2], increasing p enetration of renew able sources [3], and the retiremen t of coal p o w er plan ts [4] hav e led natural gas to be the main source for energy in U.S. and Europ ean p o wer grids. This creates an in terdep endence b et ween p o wer and gas delivery systems that comp ounds their reliabilit y and securit y risks [5]. It is b ecoming increasingly problematic for generation asset managers to pro cure natural gas during cold w eather ev ents [6] for base load and ancillary services, which ma y lead to load curtailments [7, 8]. Critically , the use of natural gas as a generation fuel source has altered the traditional status quo in which gas consumption is relatively steady and predictable throughout the da y [9, 10]. This transition comp els the need for analysis to ols that incorp orate transient pip eline flo w mo deling 1 [11 – 13]. A deficiency in formal co ordination b et ween pow er and gas delivery op erations motiv ates the fo cus of our study on mo deling the comp ounded impacts that may b e caused by a malicious cyb er attack er that seeks to maximize op erational disruption across netw orks. In some areas of the United States, a compressor station ma y b e critical to supply natural gas to fuel up to 7,000 MW [8] of com bined electricit y generation capacit y . An unplanned outage affecting these resources w ould disrupt grid op erations, and th us it is crucial to understand the risk and consequence of a meditated attack on suc h an asset. In general, cyb er-securit y threats to energy deliv ery systems app ear in the form of unauthorized commands (e.g. op ening a circuit break er in the p o wer system, a control v alve in the natural gas net work, or by changing relay settings), or of obfuscation of situational a wareness (e.g. false data injection attac ks that lead to sub-optimal op erator decisions). In the former case, an attack er who is willing to disrupt op erations may use a kno wn or zero-day vulnerabilit y to gain root access to a h uman machine in terface (HMI). Op en commands could b e sen t to a break er and cause other rela ys to trip more lines if remedial actions are not taken in time. In the latter case, malicious agents ma y target communication links that are used to deliv er data pac kets through the sup ervisory con trol and data acquisition (SCADA) netw ork. Compromising a single communication link ma y affect control ov er multiple measuremen ts, because measuremen t data is multiplexed at the remote terminal unit (R TU). Net work security analysis in volv es the sp ecification of attack graphs and trees. Attac k graphs ma y o ccur in v arious forms, and attack v ectors are challenging to mo del. There exists a s ubstan tial literature on analysis of attac k graphs using Bay esian net works [14 – 16]. These studies use securit y metrics based on graph theory [17] or vulnerabilit y metrics [18, 19] suc h as the common vulnerability scoring system (CVSS) to study vectors of attac k. Ho wev er, these previous studies do not consider the importance of the evolving physical state and con trol protocols of the system, and do not examine the physical damage that could result in an energy deliv ery net work. W e argue that in a cyb er-ph ysical systems (CPS), cyb er security is b est understo o d in the con text of mo dels that simulate the underlying physical evolution and real-time op erational control. Ph ysical and economic metrics can then quantify the impacts of a giv en attack v ector on an industrial control system (ICS). Our main con tribution is a modeling approac h that incorp orates these considerations when ev aluating cyb er-attac k impacts on infrastructure systems. This study presents a sto c hastic adv ersarial mo del for quantifying vulnerability of in terdep en- den t gas-electricit y deliv ery infrastructures to op erational disruptions caused by cyb er-attac ks. Our approac h is inspired b y a previous study [20], which examined attack v ectors for p o w er systems as a Marko v decision pro cess (MDP) that represents attac kers attempting to p enetrate a corp orate net work through the internet, and whic h is used to em ulate the b eha vior of a strategic agen t. In that mo del, attack ers ma y gain access to a p ow er system control ro om from where actions (i.e. op ening a circuit breaker) may b e p erformed. Our study fo cuses on actions by non-strategic at- tac kers to gain access and con trol of a gas pip eline netw ork that pro vides fuel for generation that constitutes a significan t p ortion of the generation mix of a dep enden t p o wer grid. W e examine ho w cyb er-ph ysical attacks on gas pip eline comp onen ts can affect gas flow dynamics, sub ject to pip eline engineering and op erating constrain ts, in order to quantify impacts on the ability of the pip eline to service gas-fired p o wer plan ts and, in turn, the consequences on electric p o w er systems securit y . Our approach is intended to enable metho ds to uncov er vulnerabilities that arise from the in terdep endence b et w een these energy systems. In particular, w e seek to represent attac k vectors that p enetrate a gas cyb er-net work by acquiring credentials that allo w the attac ker to p enetrate 2 the cyb er-physical infrastructure to attain direct con trol of op erating p oin ts of pip eline gas com- pressors. The ability of a malicious cyb er-attack er to acces s and control system comp onen ts is represen ted using a finite-horizon continous-time Mark ov c hain (CTMC) [21, 22] that includes b oth the physical and cyb er states. W e prop ose that an energy managemen t system (EMS) contingency analysis mo dule could explore the likelihoo d of cyb er-compromise even ts using a state-space ap- proac h, which could could b e used offline in a N-k probabilistic outage mo del. Our formalism could b e extended to include count ermeasures, such as in trusion detection systems (IDS), honeyp ots, tamp er-proof authentication, restrictive firewall rules, etc. W e build a benchmark test case based on the Purdue netw ork architecture mo del [23]. Finally , we attempt to capture a more realistic attac ker behavior b y incorporating the p o w er and gas dynamics into the cyb er mo deling. This manuscript is organized as follo ws. Section I I presen ts the standard for hierarc hical se- cure net work architectures and describ es in teractions among cyb er no des and the ph ysical system. W e describ e a cyb er-gas-pow er w orkflow and the exchange of information among mo dules in an adv ersarial mo del for attack vectors on p o w er and gas deliv ery op erations. Sections I I I and IV, resp ectiv ely , detail the transient gas and optimal p o wer flo w mo dels used to ev aluate op erational impacts of attac k vectors. Section V presen t the results of a computational study , and w e conclude in Section VI. 2 Cyb er mo deling using Mark o v Chains The cyber-top ology of an ICS is represen ted b y a set of firewall rules that limit the ability of third parties to control op erations. The firewalls and the rules that control access b et ween zones define a graph in whic h no des are the zones that are protected by firewalls and where eac h edge represents the p ossibilit y of devices from a one zone to p enetrate the firew all that defends the next zone. W e refer to this graph as the cyb er-top olo gy . A cyb er top ology contains zones that can include the enterprise network , the manufacturing zone (e.g. the EMS), and the ar e a zone (e.g. the SCADA net work). The Purdue mo del [23] sets the current standard for firew alls that secure architectures in an ICS. All IT systems are lo cated in the enterprise zone, whic h in the case of energy delivery systems may include applications suc h as capacit y planning, maintenance scheduling, e-mail, phone, and printing services, as well as VPN and corp orate internet connections [24]. Because of securit y concerns, the Purdue mo del restricts communications b et ween the enterprise and manufacturing zones to the so-called demilitarized zone (DMZ). In p o wer and natural gas transmission systems, the man ufacturing zone is responsible for controlling op erations, and applications in this zone w ould include the plan t historian, EMS, and file serv ers, whic h are core to the reliabilit y and securit y of the p o w er systems. The con trols resulting from solving mathematical problems in the man ufacturing zone (e.g. a unit commitment and/or economic dispatch in the p o wer system) are then passed to the area zone where R TUs or programmable logic con trollers (PLCs) send commands in real-time to actuators on the system. Information from metering devices, the status of switching devices, and alerts are displa yed in workstations by h uman mac hine in terfaces (HMIs). W e prop ose that in teractions b et ween cyb er-no des and physical states of an energy system can b e mo delled using a con tinuous-time finite-state Marko v pro cess (i.e. chain), with states corresp onding to penetration of the attac ker into a sp ecific zone of the ICS architecture. In this wa y , the CTMC mo del states and state transition matrix corresp ond to the cyb er-top ology of the ICS architecture. 3 2.1 Con tin uous-time Marko v chain mo del of attack vectors W e express the Purdue mo del arc hitecture as an undirected graph G c = ( V c , E c ) where V c denotes the set of cyb er assets and E c are the set of edges or communication links. That is, an edge connecting t wo no des means that bi-directional traffic is allo w ed b et ween the no des. Cyb er attack ers are mo deled as agen ts that follow a sto c hastic pro cess { x ( t ) : t ≥ 0 } within a CTMC [21, 22]. A CTMC is represented as a tuple ( S , P ( t ) , Q, π ( t ) , λ ) defined on a contin uous time interv al T = [0 , T ] with an exp onen tial holding time distribution T i ≈ E xp ( λ i ) λ i ≥ 0 ∀ i ∈ S , where S = { 1 , 2 , . . . , |S |} denotes the discrete state space. There exists a one-to-one mapping b et w een ev ery state s ∈ S and its corresp onding cyb er asset i ∈ V c . In addition, full control of the device is granted to the agen t when a state from the Mark ov c hain is realized. That is, the access to the area zone implies control o ver the substation actuators. Additionally , π ( t ) =  p ( x ( t ) = 1) , p ( x ( t ) = 2) , . . . , p ( x ( t ) = |S | )  is the v ector of probabilities of b eing at a giv en state at time t . The matrix P ( t ) ∈ [0 , 1] |S |×|S | denotes the transition probabilit y matrix and is defined for all times t ∈ T . The set of communication links E c pro vides the structure to P ( t ), therefore the probabilit y of transitioning b et ween t w o states where no communications are allo w ed is 0. Q ∈ R |S |×|S | is the transition rate matrix whic h accoun ts for holding times. The transition rate matrix Q can b e obtained from the jump chain probability as Q ij =  λ i p ij if i 6 = j , − λ i if i = j, (1) where p ij is the em b edded Marko v jump probability b et w een the i -th and j -th states. Note that p ij 6 = [ P ( t )] ij , b ecause the system ev olves with time. Here p ij denotes the probabilit y of an attac k er to b ypass a firew all rule. The state space con tains a recurrent absorbing state that corresp onds to the detection of the attack er b y the op erator, and the remaining states are transient. Therefore, there is a unique stationary p oin t distribution. W e determine p ij according to a security score that rates the lik eliho od of a given asset vulnerability to b e exploited. Security scores can b e found in ICS-CER T [25] under advisories [20]. The holding time parameters λ i , ∀ i ∈ S , dep end on the lev el of activity in the netw ork. Also, calculating P ( t ) for the contin uous-time case is more c hallenging than in the discrete-time problem where P ( t ) = P t for t ≥ 1. The transition probabilit y dynamics ev olve according to d dt P ( t ) = P ( t ) Q, (2) with initial condition P (0) = I (i.e., no transitions are allow ed at the initial state). The equation (2) denotes Kolmogorov’s forw ard differential equation [26]. The solution to (2) is P ( t ) = e tQ ≡ ∞ X n =0 ( tQ ) n n ! . (3) F or conv enience, it is assumed that initially the attac ker resides on the in ternet, i.e., π (0) = (1 , 0 , . . . , 0). Attac ks may o ccur as a series of escalating vulnerabilit y exploitations in the commu- nication infrastructure, which may lead the attack er to the con trol ro om. A ttack ers ma y control the area zone (e.g. a substation R TU or PLC) and may alter con trols sen t to the field devices and c hange the setp oin ts. Softw are of v arious cyb er assets suc h as R TUs or PLC devices are pro vided b y v endors and may con tain kno wn vulnerabilities to an attack er [25]. The attack er may also c hange the factory settings of a sensing device (i.e. a relay) pro vided by the v endor. This relay ma y send an op en signal command to the circuit break er due to bad mo dified settings. The consequences 4 ma y include a redistribution of flo ws due to a change in the system top ology , or an imbalance in the system if the generator setp oint is c hanged. Additionally , the HMI may b e displa ying false measuremen ts, and may induce operators to b elieve that the system is in a different op erating state for whic h sub optimal actions might b e taken. The impacts of a cyb er attac k may compound the consequences of flaws in system design (e.g. resulting from inadequate transmission or capacity planning) exp erienced during extreme even ts. The range of p ossible actions tak en b y a malicious agen t with con trol access to an energy system (i.e. the discrete range of generator/compressor setp oin ts, the status of a circuit break er/transformer, or the discrete in terv al of relay settings) are represented as states in the CTMC. In other words, the CTMC state space is c haracterized by the cyb er and ph ysical top ology (i.e. the differen t hi- erarc hical levels corresp onding to the Purdue mo del, and substations) as well as set of actions an attac ker may take to disrupt ph ysical op erations. An example of the CTMC graph is shown in Fig. 1 When the net work gets attac k ed, i.e. when the Marko v c hain lands in a state that has the ability to control op erations in real-time, a sequence of problems are solved to ev aluate the op erational impact of suc h attac ks as presented in the following section. In this pap er, w e consider tw o ICS, namely the gas and electric grid interactions in the presence of a cyb er-attac k to the gas netw ork. 2.2 Mo delling the interactions b et w een p o wer, cyb er and gas Electricit y and natural gas deliv ery op erations rely on differen t practices to ensure reliabilit y , re- siliency , and securit y . Op erator decisions can maintain normal functions of the system ev en given congestion and unexp ected conditions. Therefore, the mo delling of ph ysics and real-world prac- tices is paramount in assessing the impact of a cyber even t. In the electric p ow er sector, op erations t ypically rely on the solution of t wo core optimization problems: 1) securit y-constrained unit com- mitmen t (SCUC), whic h is used for day-ahead scheduling of generators and reserve pro curement; and 2) a DC optimal p o wer flo w (DCOPF), which is used to re-dispatch fast reacting units (e.g. gas turbines) in near real-time to comp ensate for im balances arising from imp erfect day-ahead forecasts. While the unit commitmen t problem in p o wer systems is a well established framework (see e.g. [27 – 30] and “security-constrained” v ariations of the problem [31, 32]), gas pip eline op er- ations usually dep end on lab or-in tensive heuristic approaches. Physical gas flo w control is often reactiv e and ma y be sub-optimal. The use of optimization for pipeline op erational sc heduling under transien t conditions is an active area of research [9, 11, 12, 33]. The w orkflow of the mathematical problems we solve to study the impacts of cyb er attac ks is presen ted in Fig. 1. The SCUC pro vides the sc hedules and natural gas nominations for p o w er plan ts in the day-ahead markets, and the DCOPF in real-time re-dispatc hes generation as needed. F or simplicit y , w e mo del synchronous 24-hour mark et perio ds, although this does not hold in practice in man y regions. Because pip eline op erations rely en tirely on the judgement of exp erienced op erators rather than on optimization-based decision supp ort, we use transien t optimal gas flow (TOGF) as a b est c ase proxy for the actions of gas con trol op erators to determine compressor op erating setp oin ts. T o ev aluate the op erational state of the pip eline system after a con tingency , we use a transien t gas flow (TGF) initial b oundary v alue problem (IBVP) sim ulation. In the cyb er domain, w e mo del attack v ectors through one CTMC, as is illustrated in Fig. (1), that corresp onds to the cyb er-topology of the natural gas net w ork. In this study , w e only consider the case of attacks that attempt to access and con trol compressor stations, but the mo del is easily extensible to attacks on other system components as well as on the p o w er system. When an attack er causes a compressor 5 Figure 1: Summary of p o w er, cyb er, and gas netw ork mo delling. On the left side, w e show the agen ts mo delled in the problem that include op erators, SCADA comm unications, the attac ker, and the gov erning physics for the p ow er grid and gas pip eline netw ork. In the cen ter of the figure, eac h blo c k denotes our mo del for the decisions or action taken ov er time in sequence aligned v ertically with the v arious agen ts that are op erating in the coupled cyb er-ph ysical system. The times when the attac ker penetrates the system, enters the con trol area zone, causes a contingency , and is detected are denoted b y t 0 , t 1 , t 2 , and t 3 , resp ectively . On the right, we show an example of a Mark ov c hain where the state space includes each of the areas in the Purdue mo del. Eac h area zone is coupled with a ph ysical gas no de of a 4-no de pip eline net work test case (i.e. a substation) where HMIs, PLCs and R TUs are lo cated. Dashed edges connect cyb er and ph ysical no des. Dotted edges denote the transitions to the detection state. Below the CTMC, a sample case of in terdep enden t pow er-gas net works is sho wn. The blo c ks that are gray ed-out are not explicitly co ded in our numerical model. 6 con tingency , as observ ed in a transition of the CTMC to an appropriate state, we supp ose that the b oost ratio of the compressor drops to unity , so that downstream pip eline pressure may drop b elo w minim um levels. If the pressure is sufficiently reduced, the pip eline op erator will detect the even t and issue a w arning to the p o wer system op erator. When pressure is not restored to adequate v alues, gas supplies to p o wer plan ts will b e curtailed. This will require shutdo wn of gas-fired generators, whic h ma y b e providing op erational reserv es to the p o wer grid. In resp onse, the p o w er system will re-solve a DCOPF to dispatch spinning reserves and clear the contingency . Ho wev er, p o w er system reliabilit y practices that require the N-1 security criteria could b e insufficient if multiple gas-fired generators rely on natural gas deliv ered b y a pip eline that depends on a single compressor, and th us, ma y not main tain sufficient line pressure [34]. Also, op erational generation reserves may not b e av ailable when the con tingency o ccurs if they are pro cured from gas-fired p o wer plants. Ultimately , such a c hain of ev ents will result in p o wer load curtailmen t. In the follo wing sections, w e describ e the conceptual and mathematical approac hes, as well as numerical metho ds, that w e use to simulate the decisions that pow er grid and gas pip eline op erators w ould tak e for active managemen t of in teracting systems under attac k. 3 Mo deling of gas system op erations Historically , consumption of natural gas has b een slowly v arying and highly correlated to the w eather [35]. Natural gas pip eline managers hav e used simple capacity-based mo dels to ev aluate the abilit y of natural gas pip elines to supp ort transp ort con tracts in the da y-ahead mark et [12]. The increase of natural gas usage for electricity generation has created significant intra-da y fluctuations in gas consumption. Gas pip eline op erators ma y rely on op erating pro cesses that assume steady ratable use of natural gas [9, 36 – 38], which do not accurately capture large fuel offtak e ramps caused b y the use of gas-fired generators to balance the v ariabilit y and uncertain t y in net-load [39]. This comp els dev elopment and use of mo dels that capture transien t gas flo w phenomena. 3.1 Mo delling transient gas dynamics in a net work Recen t studies ha ve prop osed new transien t sim ulation and optimization tools that mo del in tra-day transien t gas dynamics in large-scale pip elines, their effect on pip eline capacit y , and the conse- quences for dep enden t electric p ow er systems [13]. W e consider the case of a gas netw ork where gas can b e injected, compressed or withdrawn. An example of a single pip eline with a compressor station is sho wn in Fig. 2. A natural gas pip eline netw ork can be represen ted as a directed graph G g = ( V g , E g ) where V g represen ts the set of nodes i ∈ V g where gas can b e withdra wn or injected. E g represen ts the set of gas pip elines e = ( i, j ) ∈ E g that connect tw o no des i, j where gas can b e compressed with a ratio of outlet to inlet pressures α e > 1, or transp orted (i.e., α e = 1). Also, w e denote F i = { ( i, j ) | ( i, j ) ∈ E g and i ∈ V g ∀ j } , T i = { ( j, i ) | ( j, i ) ∈ E g and i ∈ V g ∀ i } as the set of edges connected to the i -th “from” or “to” node, resp ectiv ely . W e describ e the evolution of compressible gas along a pip eline using a simplification of the 1-D Euler equations. W e apply a non-dimensional transformation to the Euler equations in a single pipe following prior w ork [9], to yield ∂ t ρ e ( t, x e ) | {z } Accumulation of mass + ∂ x φ e ( t, x e ) | {z } Net mass flux = 0 , (4) 7 Figure 2: A pip eline with a compressor station ∂ t φ e ( t, x e ) | {z } Inertial forces + ∂ x ρ e ( t, x e ) | {z } Pressure forces = − λ e ` 2 D e φ e ( t, x e ) | φ e ( t, x e ) | ρ e ( t, x e ) | {z } F riction forces , (5) where (4) represen ts the con tinuit y equation (conserv ation of mass) and (5) is the momen tum conserv ation equation. The righ t hand side of (5) represents the Darcy-W eisbach equation where λ is the Darcy friction factor. ρ e and φ e denote the gas density and flux in the direction of motion ( x ), resp ectiv ely . These v ariables are defined for the domain [0 , L e ] × [0 , T ] b eing L e and T the length of the pip e and time horizon. D is the diameter of the pip e and c denotes the sp eed of the sound defined b y the ideal gas law as p = c 2 ρ where p is the pressure. It should b e noted that (4) and (5) assume slo w dynamics ( v ≤ 15 m/s ), neglect gravitational forces (horizontal pip es), and assume no temp erature gradien ts (long pip es and high pressures). Density compatibilit y conditions related to pip e pressures are giv en by ρ e ( t, 0) = α e ( t ) ρ i ( t ) ∀ i : e ∈ F i , (6) ρ e ( t, L e ) = ρ i ( t ) ∀ i : e ∈ T i . (7) Constrain ts on no dal and compressor outlet densities are ρ i ≤ ρ i ( t ) ≤ ρ i ∀ i ∈ V g , (8) α e ( t ) ρ i ( t ) ≤ ρ i ∀ i : e ∈ F i , (9) where ρ i , ρ i denote the lo w er and upp er density limits, resp ectively . Similarly , we imp ose limits on compressor ratios, 1 ≤ α e ( t ) ≤ α e ∀ e ∈ E . (10) where α e is the compressor b o ost ratio upp er limit. The gas flo w balance equation at a no de is defined as d g i ( t ) = X e ∈T i A e φ e ( t, L e ) − X e ∈F i A e φ e ( t, 0) , ∀ i ∈ V g , (11) where d g i denotes the net gas injection at the i-th no de and A e is the cross-sectional area of the pip e e . Similarly we impose limits on gas injections or withdra wals with the constraints d g i ≤ d g i ( t ) ≤ d g i ∀ i ∈ V g (12) 8 where d g i and d g i are the minim um and maximum withdraw als at the no de. In addition, the energy consumed by the compressor stations is giv en as J G = X ∀ e ∈E A e Z T 0 | φ e ( t, 0) | η e ( α ( γ − 1) /γ e ( t ) − 1) dt, (13) where η e is a scaling factor that includes compressor efficiency and γ is the gas specific heat capacity ratio [40]. W e use a reduced con trol system model where the spatial dimension is discretized using a lump ed-elemen t sc heme [9]. Letting N = |V g | and E = |E g | , w e define the vector of no dal densities ρ =  ρ L 1 , . . . , ρ L N  T , the vectors of mass flux at the b eginning and end of the pip e φ 0 =  φ 0 1 , . . . , φ 0 E  T , φ L =  φ L 1 , . . . , φ L E  T , and injections or withdraw als d g =  d g 1 , . . . , d g N  T . In addition, let Φ = 1 2 ( φ 0 + φ L ). W e define the weigh ted incidence matrix B ∈ R N × E b y B ij =    1 edge j en ters no de i, i.e., j = i − α j edge j leav es no de i, i.e., j = i + 1 0 else (14) where A = sign( B ) is the incidence matrix. Similarly , w e define A s and B s as the incidence matrix corresp onding to the no des with known pressure, and B d and A d denote the remaining ro ws. W e define the diagonal matrices Λ, K and X with diagonal entries Λ ee = L e , K ee = `λ e /D e and X ee = π D 2 e / 4 defined ∀ e ∈ E g . In matrix-v ector notation, the system dynamics (4), (5), and (11) no w b ecome | A d | X Λ   B T d   ˙ ρ = 4 ( A d X Φ − d g ) − | A d | X Λ   B T s   ˙ s (15) ˙ Φ = − Λ − 1  B T s s + B T d ρ  − K g  φ,   B T s   s +   B T d   ρ  (16) where the function g : R M × R M + → R M is defined by g j ( x, y ) = x j | x j | /y j . It should b e noted that (15) is a result of combining (4) and (11). Refer to [9, 41] for the full deriv ation of (15) and (16). W e use the ab ov e state space mo del to formulate and solv e the following fundamental gas transp ort problems. 3.1.1 T ransient optimal gas flo w (TOGF) W e solve and optimal transien t control problem where compressor ratios are minimized sub ject to the transient gas flo w dynamic constrain ts and operating limits as follows min α j , ∀ j compressor energy J G (13) s.t. system dynamics (15) − (16) pressure limits (8) − (9) compressor limits (10) supply and demand limits (12) (17) 3.1.2 T ransient gas flo w (TGF) W e solve the differential equations (15)-(16) as an implicit Differential Algebraic Equation (D AE). The DAEs are structured as follows 0 = F   ˙ ρ ˙ Φ  ,  ρ Φ  ,  α d g   = F ( ˙ x |{z} gradient , x |{z} state , u |{z} controls ) (18) where ρ and Φ are state v ariables, x ∈ R N + E , and α, d g are the control v ariables u ∈ R |V d g | + |E c g | . 9 4 Mo deling of p o w er system op erations A p o w er system netw ork is represented as a graph G p = ( V p , E p ) where V p is the set of electric buses and E p is the set of branc hes. Also, let N p = |V p | and E p = |E p | . W e formulate a SCUC problem for day-ahead sc heduling and a DC-OPF for generation dispatch in real-time. 4.1 Securit y-Constrained Unit Commitmen t (SCUC) North American Energy Reliabilit y Corp oration (NERC) standards require p ow er systems to func- tion reliably given an outage of any p o wer system asset [42]. The SCUC form ulation presented b elo w guarantees that the pow er system is secure against a generator outage contin gency . W e solv e an optimization problem defined on a discrete time in terv al T = { 1 , 2 , . . . , N } that optimizes op erational cost: J uc P = P T gt c g + u T gt c nl g + v T gt c su g + w T gt c sd g + r T gt c r g , ∀ t ∈ T , (19) where P gt ∈ R N p is the vector of p o wer injections, u gt ∈ { 0 , 1 } N p and v gt , w gt ∈ [0 , 1] N p are the v ector of commitmen t, start-up and shut-do wn binary v ariables, resp ectiv ely . W e define r gt ∈ R N p as the v ector of reserv e capacities. Similarly , c g , c nl g , c su g , c sd g , c r g ∈ R N p are the v ector of generator, no-load, start-up, shut-do wn, and capacity reserv e costs, resp ectiv ely . W e use the subscript t to index the elements of P g as provided in T . The p o wer flo w thermal limits are P ` ≤ P T D F R  P gt − d p t  ≤ P ` , ∀ t ∈ T , (20) where d p t ∈ R N p is the v ector of p ow er withdraw als and P ` and P ` are the minim um and maximum p o w er line thermal limits, respectively . P T DF R ∈ R E p × N p denotes the p o wer transfer distribution factor with resp ect to the reference bus, and can b e calculated as P T DF = H B 0− 1 where H ∈ R E p × ( N p − 1) is the p o wer flow Jacobian and B 0 ∈ R ( N p − 1) × ( N p − 1) denotes the susceptance matrix excluding the row and column that corresp ond to the reference bus. Constraints on generator limits are given b y diag( P g ) u gt + r gt ≤ P gt ≤ diag( P g ) u gt − r gt , ∀ t ∈ T , (21) where P g and P g are the minim um and maximum generator limits, resp ectiv ely . The op erator diag( . ) returns a diagonal matrix with the elemen ts of the v ector argumen t. Generation m ust balance system load at all times, as enforced b y 1 T P gt = 1 T d p t , ∀ t ∈ T , (22) and we impose ramp rate limits on the generation: P gt − P g , t − 1 ≤ diag( R hr g ) u g , t − 1 + diag( R su g ) v gt , ∀ t ∈ T , (23) P g , t − 1 − P gt ≤ diag( R hr g ) u gt + diag( R sd g ) w gt , ∀ t ∈ T , (24) where R hr g , R su g , R sd g ∈ R N p are the hourly , start-up, and sh ut-down ramping rates. The commit- men t constraints are expressed using binary v ariables, and are given as t X s = t − t u g +1 v gs ≤ u gt , ∀ t ∈ { t u g , . . . , N } , (25) 10 t X s = t − t d g +1 w gs ≤ 1 − u gt , ∀ t ∈ { t d g , . . . , N } , (26) v gt − w gt = u gt − u g , t − 1 , ∀ t ∈ T , (27) where t u g and t d g are the minim um up and down times of generator g . Allo cation of reserve capacity is guaranteed b y imp osing the constraints 1 T r gt ≥ 1 T (0 . 07 d p t ) , ∀ t ∈ T , (28) 1 T r gt ⊗ 1 ≥ P gt + r gt , ∀ t ∈ T , (29) r gt ≤ diag( R g ) u gt , ∀ t ∈ T . (30) In other w ords, the system must b e able to allo cate enough capacit y to cov er the outage of any single generator or 7% of the load at any giv en time, whichev er is larger. The SCUC optimization problem is min P gt op erational cost J uc P (19) s.t. p o w er flow limits (20) generator limits (21) p o w er balance (22) ramping limits (23) − (24) commitmen t constraints (25) − (27) reserv e capacity constrain ts (28) − (30) (31) 4.2 DC Optimal P o w er Flow (DC-OPF) T o em ulate economic re-dispatch in real-time electricity mark ets, we use a DC-OPF form ulation ∀ t ∈ T = { 1 , 2 , . . . , N } where the commitmen t v ariables and reserves are passed as parameters. The optimization problem is min P gt P T gt c g s.t. p o w er flow limits (20) generator limits (21) p o w er balance (22) ramping limits (23) − (24) (32) 4.3 Ph ysical coupling b et w een p o w er and gas infrastructures W e model the interconnections b et w een the wholesale gas and electricity sectors created b y gas-fired generators using quadratic heat curv es of the form d g i [ t ] = β ( a · P 2 it + b · P it + c ) ∀ i ∈ V p ∩ V g , ∀ t ∈ T , (33) where a , b , and c are heat rate co efficien ts and β : MMBTU/MWh → kg/s is a conv ersion factor. 11 T able 1: Coupling of pip eline no des to gas-fired generators. Gas No de Generator ID T otal Capacity (MW) 6 4,20,21 255 8 5,22,23 255 12 11,24,25 255 13 12,26,27 255 18 13,28,29 255 19 15,31 231 24 36,37,38 228 25 39,40 152 5 Sim ulation Case Study and Results W e presen t a study to sho wcase the numerical results that can b e obtained by concatenating the mo dules describ ed in Fig. 1. One realization of the w orkflow is a sto chastic pro cess where the nonlinear evolution in the state of the energy netw orks dep ends on the time in terv al during which the Mark ov chain is in a state that represen ts compromised gas compressor op erations. In order to quan tify the distribution of the sto c hastic pro cess and the resulting effect on energy delivery reliabilit y , w e p erform Mon te Carlo simulation of the stochastic process and examine its statistics. T est System: The cyb er-gas-p o w er numerical analysis is done on the IEEE-118 test case [43] for the p o w er system coupled with a 25-no de gas test case from [11] with fiv e compressor stations. W e model the coupling of gas-fired generators in the gas net work as shown in T able 1. The generator heat rate coefficients are sho wn in T able 2. The natural gas cyb er-net w ork follo ws the Purdue mo del discussed in Sec. 2, resulting in the CTMC shown in Fig. 1. W e use one absorbing state to denote the attack detection state. Attac ks occur as a series of firew all vulnerability exploitations that lead to compressor station contingencies and are mo delled by the transition probabilities. W e assume that compressor stations with higher impact in the a radial gas net w ork to ha v e higher vulnerabilit y scores, thus, reflecting the fact that some areas will b e more difficult to access. T able 2: Gas-fired generator heat rates [44] a b c Com bustion turbine 4.46 -9.95 15.11 Steam turbine 1.7 -3.87 12 Com bined cycle 5.8 -11.2 12.87 Implemen tation: W e use Julia/JuMP v0.19.2 to solve the optimization problems. JuMP is used to call the non-linear solver Ip opt [45] to solv e the non-con v ex, non-linear TOGF, and Gurobi is called for the SCUC and DCOPF whic h are mixed-integer (MIP) and linear (LP) programs, resp ectiv ely . The ph ysical flows in the natural gas net work are mo delled using the TGF IBVP presen ted in Section 3. The DAEs are solved by using the solv er ID A in Sundials. IDA uses a v ariable-order, v ariable-co efficien t Bac kward Differen tiation F ormula (BDF) with a mo dified New- ton iteration metho d and is called through the Julia pack age “DifferentialEquations.jl”. The initial conditions are obtained by solving a steady-state optimization problem, i.e. we set ˙ ρ = 0 and 12 ˙ Φ = 0 in (17), for whic h we use an absolute tolerance of 10 − 5 and a relativ e tolerance of 10 − 6 . W e obtain 10,000 samples of the Mark ov c hain o ver a 24-hour op erational da y through the Monte Carlo sim ulation. The total pro cessing time is less than 30 minutes using an In tel i7 4-core pro cessor 2018 Macb o ok Pro. Sim ulation w orkflow: First, w e compute a baseline SCUC and DCOPF solution for the pow er system, which results in a gas load at the no de. Figure 3 illustrates one instance of the physical system states at gas no de 6 that services gas-fired generators 4, 20 and 21 with a com bined capacit y of 255 MW. Note that the nominal pressure at that no de is just ov er 500 psi, which is the minim um pressure bound in the net work. A TOGF is then used to compute the gas compressor profile and no dal gas pressure under normal op erations (sho wn in blue in the figure). The CTMC is used to compute the timing of a compressor outage in the form of contingency and reco very times, b et ween whic h an upstream gas compressor pressure bo ost ratio drops to α e ≡ 1. W e then p erform a TGF sim ulation to show the resulting effect on the no dal pressure (shown in orange in the figure). When the pressure drops b elo w minim um lev els, w e suppose that the gas-fired generator is curtailed, whic h results in an acceptable no dal pressure (shown in green in the Figure). Ho wev er, the pro duction of the curtailed p ow er plant m ust b e replaced with more exp ensiv e resources, which increases the economic cost of pro duction. W e th us quantify the impacts of cyb er-attac ks b y the resulting cost of reserves and b y the p ercen tage of electric load shed. The range of outcomes that could result from a contingency in a natural gas compressor station is given in Fig. 4. The 25-node tree net work heavily relies on compressor 1 to sustain a minim um pressure. A con tingency in compressor 1 will result in natural gas fuel curtailments to all gas- fired generators. As a consequence, the p o w er system op erator will re-dispatc h its most exp ensiv e generator to balance the system, deplo y reserv es and sometimes curtail p o wer load with results in an increase in cost up to 150%. In fact, such an outcome w ould be exp ected for actual pip elines, whic h ha ve top ological tree structures. Another scenario may inv olv e an attac ker targeting compressor 2, whic h provides natural gas to all generators connected to gas no des 6 and 8. The p o w er system is not complian t with even ts that cause multiple contingencies and thus suc h an even t ma y compromise p o w er system reliabilit y . In Fig. 5, we examine how curtailed gas loads correlate to the increase in op erational cost for p o w er systems, and also examine the p o wer load curtailment with resp ect to the op erational cost increase and curtailed gas loads. W e see a strong correlation b et ween the op erational cost and curtailed gas load ( R = 0 . 9813), and also observe that there exist cases where the increase in cost also results in the highest p o wer load shedding. 6 Conclusion This study provides a detailed simulation framew ork to rigorously analyze the impacts of malicious cyb er-attac ks on the physical and operating states of interdependent natural gas and electric pow er transmission infrastructures. W e incorp orate the p ow er and natural gas dynamics in to the cyb er mo deling to b etter capture the consequences of attac ks on an energy system. The prop osed mo del- ing approach can b e extended to accommo date any ICS with any firewall p olicy in place. Through a n umerical study inv olving cyb er attac ks restricted to gas compressor contingencies in a pip eline, w e sho w that suc h an attack can cause significan t cost increase and disruption in the op erations of the p ow er grid. Our results suggest that N-1 security standards ma y not suffice to guarantee reli- able p o wer system op eration. Additionally , co ordinated pow er and gas delivery operations w ould increase the reliability of b oth sectors. Going forward, w e will extend our attack construction to 13 Figure 3: Sample of the cyb er-gas-p o w er simulation workflo w sho wn for a no de in the gas test net work at no de 6. T op: compressor ratio at station 3; center: no dal gas pressure; b ottom: gas load at the no de. TOGF solutions without c on tingency (blue), given a contingency (orange), and correction by generator curtailmen t (green). 14 Figure 4: Compressor con tingency b o xplot. The y-axis sho ws the p o w er op erational cost increase as a p ercen tage of the base case without con tingencies more p oten t adversaries that can target not just gas compressor op erations, but other system com- p onen ts such as control v alves, storage fields, or p o wer system equipment. Our future work will also extend the mo deling of cyb er-attack er b eha vior by emplo ying a Marko v decision pro cess that in volv es observed system states, rather than a single Marko v c hain only . 7 Ac kno wledgemen ts This study was carried out as part of the “Dynamical Mo deling, Estimation, and Optimal Con trol of Electrical Grid-Natural Gas T ransmission Systems” Pro ject for the D.O.E. Office of Electricit y Adv anced Grid Research and Developmen t program. W ork conducted at Los Alamos National Lab oratory w as done under the auspices of the National Nuclear Securit y Administration of the U.S. Department of Energy under Contract No. 89233218CNA000001. It was also supp orted in part b y the Director, Office of Electricit y Deliv ery and Energy Reliabilit y , Cyb ersecurit y for Energy Deliv ery Systems program, of the U.S. Department of Energy , under contract DOE0000780. An y opinions, and findings expressed in this material are those of the authors and do not necessarily reflect those of the sp onsors. References [1] US Energy Information Administration, “US Energy F acts Explained,” 2018. [Online]. Av ailable: https://www.eia.go v/energyexplained/?page=us { } energy { } home [2] ——, “Natural Gas Prices,” 2019. [Online]. Av ailable: h ttps://www.eia.go v/dnav/ng/ng { } pri { } sum { } dcu { } nus { } m.htm [3] A. Zerv os and R. Adib, “Renew ables 2018 Global Status Rep ort,” Renew able Energy Policy Netw ork for the 21st Cen tury , T ech. Rep., 2018. [Online]. Av ailable: h ttp://www.ren21.net/status- of- renewables/global- status- rep ort/ 15 Figure 5: T op: Operational p ow er cost increase vs. curtailed gas loads. Center: P ow er load shed vs. curtailed gas loads. Bottom: Po wer load shed vs. op erational p o wer cost increase 16 [4] S. F easter, “Record Drop in U . S . Coal-Fired Capacity Lik ely in 2018,” Institute for Energy Economics and Financial Analysis, T ech. Rep. Octob er, 2018. [5] North American Electric Reliability Co orp oration, “Reliability Guideline. Gas and Electrical Op erational Co ordination Considerations,” T ech. Rep. September, 2018. [6] PJM In terconnection, “Analysis of operational even ts and mark et impacts during the January 2014 cold weather ev ents,” T ech. Rep. January , 2014. [On- line]. Av ailable: h ttps://www.p jm.com/ { ∼ } /media/library/rep orts- notices/w eather- related/ 20140509- analysis- of- op erational- even ts- and- mark et- impacts- during- the- jan- 2014- cold- weather- ev en ts. ashx [7] ISO New England, “2018 Regional Electricity Outlo ok,” T ec h. Rep., 2018. [Online]. Av ailable: h ttps://www.iso- ne.com/static- assets/documents/2018/02/2018 { } reo.pdf. [8] ——, “Op erational F uel-Security Analysis,” T ech. Rep., 2018. [9] A. Zlotnik, M. Chertko v, and S. Backhaus, “Optimal con trol of transient flo w in natural gas net works,” in Pr o c. IEEE Conf. on De cision and Contr ol , vol. 54rd IEEE, no. Cdc. IEEE, 2015, pp. 4563–4570. [10] A. V. Zlotnik, P . A. Ruiz, A. M. Rudk evich, S. N. Backhaus, R. G. Carter, and J. D. T aft, “Grid architecture at the gas-electric interface,” GRID Mo dernization Lab oratory Consortium, U.S. Department of Energy , T ech. Rep., 2017. [11] A. Zlotnik, L. Roald, S. Bac khaus, M. Chertko v, and G. Andersson, “Con trol policies for op erational co ordination of electric p o wer and natural gas transmission systems,” pp. 7478– 7483, 2016. [12] A. V. Zlotnik, A. M. Rudkevic h, E. Goldis, P . A. Ruiz, M. Caramanis, R. G. Carter, S. N. Bac khaus, R. T ab ors, and D. Baldwin, “Economic Optimization of Intra-Da y Gas Pip eline Flo w Schedules using T ransien t Flow Mo dels,” in Pip eline Simulation Inter est Gr oup 2017 A nnual Me eting , Atlan ta, Georgia, 2017. [13] A. Zlotnik, L. Roald, S. Backhaus, M. Chertko v, and G. Andersson, “Co ordinated Scheduling for Interdependent Electric Po w er and Natural Gas Infrastructures,” IEEE T r ansactions on Power Systems , vol. 32, no. 1, pp. 600–610, 2017. [14] N. Poolsappasit, R. Dewri, and I. Ray , “Dynamic securit y risk management using ba yesian attac k graphs,” IEEE T r ansactions on Dep endable and Se cur e Computing , vol. 9, no. 1, pp. 61–74, 2011. [15] Y. Zhang, L. W ang, Y. Xiang, and C.-W. T en, “Po w er system reliability ev aluation with scada cyb ersecurit y considerations,” IEEE T r ansactions on Smart Grid , vol. 6, no. 4, pp. 1707–1721, 2015. [16] E. Miehling, M. Rasouli, and D. T eneketzis, “Optimal defense policies for partially observ able spreading processes on bay esian attac k graphs,” in Pr o c e e dings of the Se c ond ACM Workshop on Moving T ar get Defense . ACM, 2015, pp. 67–76. 17 [17] O. V uko vic et al. , “Net work-a ware mitigation of data in tegrity attac ks on p o w er system state estimation,” IEEE Journal on Sele cte d A r e as in Communic ations , v ol. 30, no. 6, pp. 1108–1118, 2012. [18] W. Bo yer and M. McQueen, “Ideal based cyb er securit y technical metrics for con trol systems,” in International Workshop on Critic al Information Infr astructur es Se curity . Springer, 2007, pp. 246–260. [19] P . E. Blac k, K. Scarfone, and M. Souppay a, “Cyb er security metrics and measures,” Wiley Handb o ok of Scienc e and T e chnolo gy for Homeland Se curity , pp. 1–15, 2008. [20] K. R. Da vis, C. M. Davis, S. A. Zonouz, R. B. Bobba, R. Berthier, L. Garcia, and P . W. Sauer, “A cyb er-ph ysical mo deling and assessment framew ork for p o w er grid infrastructures,” IEEE T r ansactions on smart grid , v ol. 6, no. 5, pp. 2464–2475, 2015. [21] W. J. Anderson, Continuous-time Markov chains: A n applic ations-oriente d appr o ach . Springer Science & Business Media, 2012. [22] G. G. Yin and Q. Zhang, Continuous-time Markov chains and applic ations: a singular p ertur- b ation appr o ach . Springer, 2012, vol. 37. [23] L. Obregon, “Secure Arc hitecture for Industrial Con trol Systems,” SANS Institute, T ech. Rep., 2019. [24] P . Didier, F. Macias, J. Harstad, R. Antholine, S. A. Johnston, S. Piy evsky , M. Schillace, G. Wilcox, D. Zaniewski, and S. Zuponcic, “Con verged plan t wide ethernet (cp we) design and implemen tation guide,” Cisc o Systems and R o ckwel l Automation , 2011. [25] “Industrial control systems, advisories and rep orts,” https://www.us- cert.go v/ics, accessed: 2019-09-02. [26] A. Kolmogoroff, “ ¨ Ub er die analytischen metho den in der wahrsc heinlichk eitsrec hnung,” Mathematische Annalen , vol. 104, no. 1, pp. 415–458, Dec 1931. [Online]. Av ailable: h ttps://doi.org/10.1007/BF01457949 [27] W. F. Tinney and C. E. Hart, “P ow er flow solution by newton’s method,” IEEE T r ans. Power App ar atus and Sys. , no. 11, pp. 1449–1460, 1967. [28] R. Kerr, J. Scheidt, A. F ontanna, and J. Wiley , “Unit commitment,” IEEE T r ansactions on Power App ar atus and Systems , no. 5, pp. 417–421, 1966. [29] A. I. Cohen and V. R. Sherk at, “Optimization-based metho ds for op erations scheduling,” Pr o c e e dings of the IEEE , vol. 75, no. 12, pp. 1574–1591, 1987. [30] A. J. W o o d, B. F. W ollenberg, and G. B. Shebl´ e, Power gener ation, op er ation, and c ontr ol . John Wiley & Sons, 2013. [31] J. J. Sha w, “A direct method for security-constrained unit commitmen t,” IEEE T r ans. Power Sys. , vol. 10, no. 3, pp. 1329–1342, 1995. 18 [32] A. Monticelli, M. Pereira, and S. Gran ville, “Securit y-constrained optimal p o wer flo w with p ost-con tingency correctiv e resc heduling,” IEEE T r ansactions on Power Systems , vol. 2, no. 1, pp. 175–180, 1987. [33] R. Z. R ´ ıos-Mercado and C. Borraz-S´ anc hez, “Optimization problems in natural gas transp orta- tion systems: A state-of-the-art review,” Applie d Ener gy , v ol. 147, pp. 536–555, 2015. [34] K. Sundar, H. Nagara jan, M. Lubin, L. Roald, S. Misra, R. Bent, and D. Biensto c k, “Unit commitmen t with n-1 security and wind uncertaint y ,” in 2016 Power Systems Computation Confer enc e (PSCC) . IEEE, 2016, pp. 1–7. [35] S. Edwardson et al. , “The next dimension of flow profiles,” in PSIG Annual Me eting , p. 1617. [36] P . W ong and R. Larson, “Optimization of natural-gas pip eline systems via dynamic program- ming,” IEEE T r ansactions on A utomatic Contr ol , vol. 13, no. 5, pp. 475–481, 1968. [37] B. Rothfarb, H. F rank, D. Rosenbaum, K. Steiglitz, and D. J. Kleitman, “Optimal design of offshore natural-gas pipe line systems,” Op er ations r ese ar ch , vol. 18, no. 6, pp. 992–1020, 1970. [38] C. Luongo, W. Y eung, and B. Gilmour, “Optimizing the op eration of gas transmission net- w orks,” in Computers in engine ering 1991 , 1991. [39] K. P ambour, R. Sopgwi, B.-M. Ho dge, and C. Brancucci, “The v alue of day-ahead co ordination of p o wer and natural gas netw ork op erations,” Ener gies , v ol. 11, no. 7, p. 1628, 2018. [40] E. S. Menon, Gas pip eline hydr aulics . Crc Press, 2005. [41] A. V. Zlotnik, S. N. Backhaus, M. Chertko v, and S. Dy achenk o, “Mo del Reduction and Opti- mization of Natural Gas Pip eline Dynamics,” in ASME Dynamic Systems and Contr ol Con- fer enc e , vol. 14. Los Alamos National Lab, 2015. [42] North American Electric Reliabilit y Corp oration , “Reliabilit y Guideline, Metho ds for Establishing IR OLs,” 2018. [Online]. Av ailable: h ttps://www.eia.gov/energy explained/ ?page=us { } energy { } home [43] R. Christie, “P ow er systems test case archiv e.” [Online]. Av ailable: https://labs.ece.u w.edu/ pstca/ [44] D. Lew et al. , “Impacts of wind and solar on fossil-fueled generators,” National Renew able Energy Lab.(NREL), Golden, CO (United States), T ech. Rep., 2012. [45] A. W¨ ach ter and L. T. Biegler, “On the implemen tation of an interior-point filter line-search algorithm for large-scale nonlinear programming,” Mathematic al pr o gr amming , v ol. 106, no. 1, pp. 25–57, 2006. 19