Optimal Power Flow

Division of Economics and Business Working Paper Series
A Primer on Optimal Power Flow: Theory, Formulation, and Practical Examples

Stephen Frank Steen Rebennack
Working Paper 2012-14 http://econbus.mines.edu/working-papers/wp201214.pdf Colorado School of Mines Division of Economics and Business 1500 Illinois Street Golden, CO 80401
October 2012
c 2012 by the listed authors. All rights reserved.
Colorado School of Mines Division of Economics and Business Working Paper No. 2012-14 October 2012 Title: A Primer on Optimal Power Flow: Theory, Formulation, and Practical Examples
Author(s): Stephen Frank Department of Electrical Engineering & Computer Science Colorado School of Mines Golden, CO 80401-1887 [email protected] Steen Rebennack Division of Economics and Business Colorado School of Mines Golden, CO 80401-1887 [email protected]
ABSTRACT The set of optimization problems in electric power systems engineering known collectively as Optimal Power Flow (OPF) is one of the most practically important and well-researched subelds of constrained nonlinear optimization. OPF has enjoyed a rich history of research, innovation, and publication since its debut ve decades ago. Nevertheless, entry into OPF research is a daunting task for the uninitiatedboth due to the sheer volume of literature and because OPFs familiarity within the electric power systems community has led authors to assume a great deal of prior knowledge that readers unfamiliar with electric power systems may not possess. This primer provides a practical introduction to OPF from an Operations Research perspective; it describes a complete and concise basis of knowledge for beginning OPF research. The primer is tailored for the Operations Researcher who has experience with optimization but little knowledge of Electrical Engineering. Topics include power systems modeling, the power ow equations, typical OPF formulations, and data exchange for OPF. AMS subject classications: 90-01, 90C26, 90C30, 90C90 Keywords: power ow, optimal power ow, electric power systems analysis, electrical engineering, nonlinear programming, optimization, operations research.
This work is supported by the National Science Foundation Graduate Research Fellowship under Grant No. DGE-1057607.
2
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49
STEPHEN FRANK AND STEFFEN REBENNACK
1. Introduction. The set of optimization problems in electric power systems engineering known collectively as Optimal Power Flow (OPF) is one of the most practically important and well-researched subelds of constrained nonlinear optimization. In 1962, Carpentier [8] introduced OPF as an extension to the problem of optimal economic dispatch (ED) of generation in electric power systems. Carpentiers key contribution was the inclusion of the electric power ow equations in the ED formulation. Today, the dening feature of OPF remains the presence of the power ow equations in the set of equality constraints. In general, OPF includes any optimization problem which seeks to optimize the operation of an electric power system (specically, the generation and transmission of electricity) subject to the physical constraints imposed by electrical laws and engineering limits on the decision variables. This general framework encompasses dozens of optimization problems for power systems planning and operation [11, 36, 37]. As illustrated in Figure 1.1, OPF may be applied to decision making at nearly any planning horizon in power systems operation and controlfrom long-term transmission network capacity planning to minute-by-minute adjustment of real and reactive power dispatch [14, 32, 37]. To date, thousands of articles and hundreds of textbook entries have been written about OPF. In its maturation over the past ve decades, OPF has served as a practical proving ground for many popular nonlinear optimization algorithms, including gradient methods [4, 10, 22], Newton-type methods [29], sequential linear programming [3, 27], sequential quadratic programming [7], and both linear and nonlinear interior point methods [15, 30, 31]. These OPF algorithms, among others, are reviewed in several surveys [1719, 36], including one that we recently published [11, 12]. Although OPF spans both Operations Research and Electrical Engineering, the accessibility of the OPF literature is skewed heavily toward the Electrical Engineering community. Both conventional power ow (PF) and OPF have become so familiar within the electric power systems community that the recent literature assumes a great deal of prior knowledge on the part of the reader. For example, while conducting our recent survey we found that few papers even include a full OPF formulation, much less explain the particulars of the objective function or constraints. Even introductory textbooks [26, 32, 37] require a strong background in power systems analysis, specically regarding the form, construction, and solution of the electrical power ow equations. Although many Electrical Engineers have this prior knowledge, an Operations Researcher likely will not. We believe this accessibility gap has been detrimental to the involvement of the Operations Research community in OPF research; our impression is that most OPF papers continue to be published in engineering journals by Electrical Engineers. What is missing from the literatureand what we provide in this primeris a practical introduction to OPF from an Operations Research perspective. The goal of this primer is to describe the tool set required to formulate, solve, and analyze a practical OPF problem. Other introductory texts for OPF focus heavily on optimization theory and tailored solution algorithms. In contrast, this primer places an emphasis on the theory and mechanics of the OPF formulationthe least documented aspect of OPFusing practical, illustrative examples. Because we have written this primer for the Operations Researcher, we assume that the reader has signicant experience with nonlinear optimization but little or no background in Electrical Engineering. Specically, this paper requires a solid understanding of
A PRIMER ON OPTIMAL POWER FLOW
Planning Horizon Years Months Weeks Days Hours Minutes Real-Time
asin Incre Decr
del g Mo
n lutio Reso
easin g
Unce rtain
ty
Long-term Generation Scheduling Capacity Expansion Planning Real Power Dispatch Reactive Power Dispatch Reactive Power Planning
Unit Commitment SecurityConstrained Unit Commitment
Economic Dispatch Classic Optimal Power Flow
Automatic Generation Control
SecurityConstrained Unit Commitment
Optimal Reactive Power Flow
Voltage Control
Fig. 1.1. Optimization in power system operation via incremental planning. Long-term planning procedures make high level decisions based on coarse system models. Short-term procedures ne tune earlier decisions, using detailed models but a more limited decision space. Bold text indicates planning procedures which incorporate variants of optimal power ow.
50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
linear algebra [16], complex number theory [13, 16], analysis of dierential equations in the frequency domain [16], and linear and nonlinear optimization theory and application [20, 24]. Readers who also have a working knowledge of electrical circuits [21] and electric power systems analysis [14] will nd the development of the power ow equations familiar. Other readers may wish to expand their understanding by consulting a good power systems text such as [14] or [32]. However, prior familiarity with power ow is not strictly required in order to follow the development presented in this primer. The primer begins with a guide to OPF notation in 2, including notational dierences between electric power systems engineers and Operations Researchers. Sections 3 and 4 introduce the modeling of electric power systems and the power ow equations; these fundamental topics are omitted in most other introductory OPF texts. Building upon the previous sections, 5 discusses OPF formulations; this section includes full formulations for several of the decision processes shown in Figure 1.1. 6 provides a descriptive guide to two common le formats for exchanging PF and OPF data. Finally, 7 concludes the primer.
4
67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115
2. Notation. In order to remain consistent with the existing body of OPF literature, this primer uses notation that follows the conventions of electric power systems engineering rather than the Operations Research community. Where there are significant dierences, we have added clarifying remarks. 2.1. General. Throughout this primer, italic roman font (A) indicates a variable or parameter, bold roman font (A) indicates a set, and a tilde over a symbol (a) indicates a phasor quantity (complex number). Letter case does not dierentiate variables from parameters; a given quantity may be a variable in some cases and a parameter in others. Subscripts indicate indices. Symbolic superscripts are used as qualiers to dierentiate similar variables, while numeric superscripts indicate mathematical operations. For example, the superscript L dierentiates load power P L from net power P , but P 2 indicates (net) power squared. We use the following general notation for optimization formulations: u vector of control variables (independent decision variables) x vector or state variables (dependent decision variables) vector of uncertain parameters f (u, x) objective function (scalar) g (u, x) vector function of equality constraints h(u, x) vector function of inequality constraints The symbols e and j represent mathematical constants: e Eulers number (the base of the natural logarithm), e 2.71828 j the imaginary unit or 90 operator, j = 1 Remark 2.1. This diers from the common use in Operations Research of e as the unit vector and j as an index. In Electrical Engineering, j , rather than i or , designates the imaginary unit. For this reason, we avoid the use of j as an index throughout this primer. In the examples, electrical units are specied where applicable using regular roman font. The unit for a quantity follows the numeric quantity and is separated by a space; for example 120 V indicates 120 Volts. The following electrical units are used in this primer: V Volt (unit of electrical voltage) A Ampere (unit of electrical current) W Watt (unit of real electrical power) VA Volt-Ampere (unit of apparent electrical power) VAR Volt-Ampere Reactive (unit of reactive electrical power) 2.2. Dimensions, Indices, and Sets. The following dimensions and indices are used in the OPF formulations within this primer: N total number of system buses (nodes) L total number of system branches (arcs) M number of system PQ buses i, k indices corresponding to system buses and branches c contingency case index t time period index Remark 2.2. We use L to indicate the number of system branches because B is reserved for the bus susceptance matrix. Remark 2.3. System branches are indexed as arcs between buses. For example, the branch between buses i and k is denoted by (i, k ) or ik . Remark 2.4. In the optimization community, c typically refers to a vector of objective function coecients. In this primer, however, we use upper case C for

116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
objective function coecients and reserve lower case c for the contingency case index of security-constrained economic dispatch as described in 5.2.1. There is no standard set notation within the OPF literature. (Many authors do not use sets in their formulations.) For convenience, however, we adopt the following sets in this primer: N set of system buses (nodes) L set of system branches (arcs) M set of load (PQ) buses G set of controllable generation buses H set of branches with controllable phase-shifting transformers K set of branches with controllable tap-changing transformers Q set of planned locations (buses) for new reactive power sources C set of power system contingencies for contingency analysis T set of time-periods for multi-period OPF Remark 2.5. For clarity, we use H and K to represent sets of controllable phase-shifting and tap-changing transformers rather than S (often used to designate sources or scenarios) and T (often used to designate time periods). The letters H and K otherwise have no special association with phase-shifting and tap-changing transformers. 2.3. Electrical Quantities. In power systems analysis, electrical quantities are represented in the frequency domain as phasor quantities (complex numbers). Complex numbers may be represented as a single complex variable, as two real-valued variables in rectangular form a + jb, or as two real-valued variables in polar form c ; all of these notations are found in the OPF literature. (Complex number notation is explained in more detail in 3.2.) Here, we document the usual symbols and relationships used for the electrical quantities; some notational exceptions exist in the literature. Zik Rik Xik yik gik bik
Sh yik Sh gik bSh ik S yi S gi bS i
135 136 137 138 139 140 141 142
143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163
Yik Yik ik Gik
2.3.1. Admittance. complex impedance of branch ik resistance of branch ik (real component of Zik ) reactance of branch ik (imaginary component of Zik ) Zik = Rik + jXik complex series admittance of branch ik series conductance of branch ik (real component of yik ) series susceptance of branch ik (imaginary component of yik ) yik = 1/Zik = gik + jbik complex shunt admittance of branch ik Sh shunt conductance of branch ik (real component of yik ) Sh shunt susceptance of branch ik (imaginary component of yik ) Sh Sh Sh yik = gik + jbik complex shunt admittance at bus i S shunt conductance at bus i (real component of yi ) S shunt susceptance at bus i (imaginary component of yi ) S S S yi = gi + jbi complex ik th element of the bus admittance matrix magnitude of ik th element of the bus admittance matrix angle of ik th element of the bus admittance matrix conductance of ik th element of the bus admittance matrix (real component of
6
164 165 166 167 168 169 170 171
Yik ) susceptance of ik th element of the bus admittance matrix (imaginary component of Yik ) Yik = Yik ik = Gik + jBik Remark 2.6. Note the distinction between lowercase y , g , and b and uppercase Y , G, and B : the former represents the values corresponding to individual system branch elements, while the latter refers to the admittance matrix which models the interaction of all system branches. Bik Vi Vi i Ei Fi 2.3.2. Voltage. complex (phasor) voltage at bus i voltage magnitude at bus i voltage angle at bus i real component of complex voltage at bus i imaginary component of complex voltage at bus i Vi = Vi i = Ei + jFi 2.3.3. Current. complex (phasor) current injected at bus i magnitude of current injected at bus i complex (phasor) current in branch ik , directed from bus i to bus k magnitude of current in branch ik
172 173 174 175 176 177 178
179 180 181 182 183
Ii Ii Iik Iik PiL QL i L Si
184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
2.3.4. Power. load (demand) real power at bus i load (demand) reactive power at bus i load (demand) complex power at bus i L Si = PiL + jQL i G Pi generator (supply) real power at bus i QG generator (supply) reactive power at bus i i G Si generator (supply) complex power at bus i G Si = PiG + jQG i Pi net real power injection at bus i (Pi = PiG PiL ) L Qi net reactive power injection at bus i (Qi = QG i Qi ) Si net complex power injection at bus i G L Si = Si Si = Pi + jQi Remark 2.7. The phasor indicator is omitted for complex power S , as S is always understood to be a complex quantity. In the literature, the indicator is also often omitted for V , I , y , and Y , but we include it here to disambiguate the complex quantities from their associated (real-valued) magnitudes. ik Tik 2.3.5. Other. phase shift of phase-shifting transformer in branch ik tap ratio of tap-changing transformer in branch ik
201 202 203
204 205 206 207 208
3. Fundamental Concepts. Scholarly literature discussing OPF assumes a working knowledge of power systems models and electrical concepts, many of which may be unfamiliar to the Operations Researcher. In this section, we summarize several fundamental concepts required for power systems analysis and the development of the power ow equations as presented in 4.
Bus 2
2) (1, ch n a Br
Branch (2,4)
Bus 4
Br an
ch an r B
ch (1, 3)
Branch (3,5) Bus 3 Bus 5
Fig. 3.1. Bus and branch indices in an example 5-bus electrical network.
209 210 211 212 213
3.1. System Representation. Electric power systems are modeled as a network of electrical nodes (buses) interconnected via admittances (branches) that represent transmission lines, cables, transformers, and similar power systems equipment. Buses are referenced by node with index i N, while branches are referenced as arcs between nodes (i, k ) L, where i, k N. Example 3.1. The network in Figure 3.1 has N = 5 buses and L = 6 branches, with corresponding sets N = {1, 2, 3, 4, 5} and L = {(1, 2), (1, 3), (2, 4), (3, 4), (3, 5), (4, 5)} .
214 215 216 217 218 219 220 221
For power ow analysis, the electric power system is analyzed in the frequency domain under the assumption of sinusoidal steady-state operation. At sinusoidal steady-state, all voltage and current waveforms are sinusoids with xed magnitude, frequency, and phase shift, and all system impedances are xed. Under these conditions, the dierential equations governing power system operation reduce to a set of complex algebraic equations involving the phasor representation of the system electrical quantities. This algebraic representation is much easier to solve. 3.2. Phasor Quantities. Steady-state sinusoidal voltages and currents can be fully characterized by their magnitude and phase shift using phasors. A phasor transforms a sinusoidal time-domain signal into a complex exponential in the frequency domain, using the relationship c sin (2f t + ) Time Domain cej Frequency Domain
The frequency f of the signal is xed and therefore omitted from the phasor notation. The phasor cej may be written as c in polar coordinates or as a + jb in rectangular
Branch (4,5)
Bus 1
(3
) ,4
coordinates, where, according to Eulers formula, a = c cos , b = c sin , a2 + b2 , b = arctan . a c=

222 223
Freitag and Busam [13] provide an overview of complex number theory, while OMalley [21, Ch. 11] discusses the use of phasor quantities in Electrical Engineering. Remark 3.2. In Electrical Engineering, voltage and current phasors are expressed as root-mean-square (RMS) quantities rather than peak quantities. This is done so that frequency domain power calculations yield the correct values without the need for an additional scaling factor. For a sinusoid, the RMS magnitude is 1/ 2 times the peak magnitude. Thus, in Electrical Engineering, a time-domain voltage waveform VPk sin (2f t + ) has the frequency domain phasor VPk j e . 2
224
Example 3.3. The time domain voltage waveform 120 2 sin(377t 30) V represents the standard outlet voltage in the United States with the angle referenced, for instance, to the high side of a utility distribution transformer. This voltage has the frequency domain phasor V = 120ej 30 V. In polar coordinates, this phasor is V = 12030 V, while in rectangular coordinates it is V = 120 cos 30 + j 120 sin 30 V, 103.9 j 60.0 V.
225 226 227 228 229 230 231 232
Note that the RMS voltage magnitude of 120 V is the familiar quantity. 3.3. Complex Power. The product of voltage and current is electrical power. In AC power systems, however, instantaneous electrical power uctuates as the voltage and current magnitudes change over time. For frequency domain analysis, power systems engineers use the concept of complex power to characterize these time domain power uctuations. Complex power is a phasor quantity consisting of real power and reactive power. Real power represents real work, that is, a net transfer of energy from source to load.
9
v (t )
Total Current t Quadrature
i (t ) = i d(t ) + i q(t ) i d(t ) i q(t )
Direct
Instantaneous Power
Real t Power
v (t )i (t ) v (t )i d(t ) v (t )i q(t )
Reactive Power
Fig. 3.2. Conceptual illustration of real and reactive power using time domain waveforms. In the gure, current i(t) lags voltage v (t) by 30.
233 234 235 236 237 238 239 240
Reactive power, on the other hand, represents circulating energyan cyclic exchange of energy that averages zero net energy transfer over time. Real power transfer occurs when voltage and current are in phase, while reactive power transfer occurs when voltage and current are 90 out of phase (that is, orthogonal). An arbitrary AC current i(t) can be represented by the sum of direct current id (t) (in phase with the voltage) and quadrature current iq (t) (orthogonal to the voltage). Direct current produces real power and quadrature current reactive power, as illustrated in Figure 3.2. By convention, reactive power is considered positive when current lags voltage. Therefore, complex power S can be computed from1 S = V I = P + jQ
241 242 243
and consists of orthogonal components P (real power) and Q (reactive power). The magnitude of complex power, |S |, is called the apparent power and is often used to specify power systems equipment and transmission line ratings. Complex and
1 Here and elsewhere in this primer, the symbol denotes complex conjugation rather than an optimal value. This use is typical in Electrical Engineering and consistent with most OPF literature.
10
244 245 246 247
apparent power have units of Volt-Amperes (VA), real power has units of Watts (W), and reactive power has units of Volt-Amperes Reactive (VAR). Remark 3.4. Reactive power is sometimes called imaginary power, both because it does not perform real work and because it is the imaginary part of S . Example 3.5. A small electrical appliance draws 2 A30 from a 120 V0 source. The complex power draw of the appliance is S = (120 V0) (2 A30) , = 240 VA 30, 207.8 W + j 120 VAR.
= (120 V 2 A) (0 + 30) ,
248 249
The apparent power draw of the appliance is 240 VA, the real power is 207.8 W, and the reactive power is 120 VAR. We can verify that the real power is correct by examining the average power in the time domain. The voltage and current waveforms are v (t) = 120 2 sin(377t) V, i(t) = 2 2 sin(377t 30) A, respectively. The instantaneous power is p(t) = v (t)i(t) = 480 sin(377t) sin(377t 30) W. Or, by using trigonometric identities, p(t) = 240 cos(30) 240 cos (2 (377t) 30) W. Over time, the average power is pAvg =
0
1 60
240 cos(30) 240 cos (2 (377t) 30) dt W,
= 240 cos(30) 207.8 W,

250 251 252 253 254 255 256 257 258 259 260 261
which is identical to the real power P computed in the frequency domain. Both real and reactive power aect power systems operation and are therefore modeled in PF and OPF. OMalley [21, Ch. 15] and other introductory circuits texts provide a more complete overview of complex power. 3.4. The Per-Unit System. Electric power systems quantities are usually expressed as a ratio of the actual quantity to a reference, or base, quantity; this practice is called the per-unit system. Per-unit quantities are unitless and are labeled using a designation, if any, of p.u.. Nearly all OPF literature assumes a working knowledge of per-unit on the part of the reader, but this assumption is rarely stated explicitly. Indeed, power systems texts frequently mix per-unit and SI (metric) units, for instance, reporting voltage in per-unit and power in MW. As a result, per-unit can be a signicant source of confusion when working with practical OPF formulations. In power systems analysis, base quantities are given in the SI system (Volts, Amperes, Watts, Ohms, etc.), while per-unit quantities are dimensionless. The perunit value of an SI quantity x on a given base xBase is x . xpu = xBase

262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280
11
Correct interpretation of the SI value of a per-unit quantity requires knowledge of the base quantity. For example, a power of 0.15 p.u. on a 10 MVA base is equal to 1.5 MW, but 0.15 p.u. on a 1000 MVA base is equal to 150 MW.2 All calculations that can be performed in the SI system can also be performed in per-unit. However, (i) per-unit and SI quantities cannot be mixed in calculations, and (ii) all per-unit calculations must be performed on a consistent set of bases. With a proper selection of system bases, the per-unit system has several advantages over the SI system of measurement, most notably 1. The use of per-unit eliminates the need to distinguish between single-phase and three-phase electrical quantities; 2. The use of per-unit eliminates the need to apply voltage and current scaling factors at the majority of system transformers; 3. The use of per-unit automatically adjusts for the phase shift of three-phase transformers (Wye-Delta or Delta-Wye); 4. Per-unit quantities have consistent magnitudes on the order of 1.0, which improves the numerical stability of power ow calculations; and 5. The per-unit system is easier to interpret at a glance. (For example, per-unit voltage should always lie within the approximate range 0.951.05 p.u., regardless of the SI voltage.) Once two system bases are specied, the others are xed exactly. In power ow analysis, the voltage and power bases are specied, VBase = Line-to-line Voltage, SBase = Three-phase Power, and the remaining three-phase system bases are calculated according to SBase , IBase = 3 VBase VBase ZBase = = 3 IBase 3 IBase = YBase = VBase
2 VBase , SBase
SBase 1 = . 2 VBase ZBase
281 282 283 284 285 286
SBase is constant throughout the power system, but VBase is distinct for each system bus. For mathematical convenience, power systems engineers typically set SBase to one of 10, 100, or 1000 MVA and select VBase as the nominal line-to-line voltage at each bus. When VBase is selected in this way, the voltage ratio of most system transformers becomes 1:1 in per-unit, simplifying the development of the system admittance matrix; see 4.1. Example 3.6. 12.47 kV is a common distribution voltage level in the United States. If a nominally 12.47 kV feeder is operated at 0.95 p.u., then its voltage level in SI units is VSI = 0.95 12.47 kV, = 11.85 kV.
2 In per-unit, real power (W), reactive power (VAR), and apparent power (VA) share a common base with units of VA.
12
Similarly, a nominally 12.47 kV feeder operating at 13.2 kV (another common distribution voltage level) has a per-unit voltage of Vpu = 13.2 kV , 12.47 kV = 1.059 p.u.
287 288 289 290 291 292 293
In OPF, a common convention is to specify source and load power in SI units, indicate the system power base, and specify all other quantities directly in per-unit; see 6. The power values must be converted to per-unit prior to evaluating the power ow equations, but usually no other conversions are necessary. Glover, Sarma, and Overbye [14] provide additional discussion of the per-unit system, including instructions for base conversions, relationships for single-phase bases, and practical examples. 4. The Power Flow Equations. In this section, we develop the power ow equations and present the mechanics of their construction. The steady-state operation of an alternating current (AC) electrical network is governed by the matrix equation I =YV where I = I1 , . . . , IN is an N -dimensional vector of phasor currents injected at each system bus, V = V1 , . . . , VN is an N -dimension vector of phasor voltages at each system bus, and Y11 . Y = . . ... .. . ... Y1N . . . (4.1)
YN 1
YN N
294 295 296 297 298 299
is the N N complex bus admittance matrix. In practical power systems, I represents the current supplied by generators and demanded by loads, while Y models transmission lines, cables, and transformers. Traditionally, the elements of Y were considered constant, but in newer OPF formulations Y may contain both constants and control (decision) variables. The voltages V are state variables which fully characterize the system operation for a given matrix Y . In power systems analysis, it is more convenient to work with power ows than currents because (i) injected powers are independent of system voltage angle while injected currents are not, and (ii) working directly with power allows straightforward computation of required electrical energy. Therefore, power systems engineers transform (4.1) into the complex power ow equation S =V YV
(4.2)
13
Bus i ei V
e I i yiS e
yik e ykS e
e I k
Bus k
ek V
Fig. 4.1. Example two-bus network illustrating the denitions of bus voltage and injected current.
where S = P + jQ is a vector of complex power injections at each bus and denotes element-wise vector multiplication. At each bus, the total injected power is the dierence between the generation and the load,
G L Si = Si Si ,
L Qi = QG i Qi .
300 301
Pi = PiG PiL ,
Typically, load real and reactive power are xed while generation real and reactive power are control variables with minimum and maximum limits. 4.1. The Admittance Matrix. The bus admittance matrix Y forms the core of the power ow equations. OPF data generally does not give Y directly, and therefore we summarize the mechanics of its construction here. The elements of Y are derived from the application of Ohms law, Kirchos current law (KCL), and Kirchos voltage law (KVL) to a steady-state AC electrical network; OMalley [21] provides a concise summary of these electrical laws. At each bus i, Ii is the net current owing out of the bus through all connected branches, that is, the current injected from outside sources (such as connected generators or loads) required to satisfy KCL. From Ohms law and KVL,
S,Total Ii = Vi yi +
k:(i,k)L
302 303 304
Vi Vk yik +
k:(k,i)L
Vi Vk yki
(4.3)
S,Total where yik is the admittance of branch (i, k ) and yi is the total shunt admittance from bus i to neutral. In matrix form, (4.1) is equivalent to (4.3) when the elements of Y are dened as
Yii = Yik =
305 306 307 308 309 310
Admittances directly connected to bus i Admittances directly connected between bus i and bus k (4.4)
Typically, only a single branch (i, k ) connects bus i to bus k , in which case the odiagonal elements become Yik = Yki = yik . If there is no connection between buses i and k , Yik = 0. Thus, Y is sparse, having dimension N N but only N + 2L nonzero entries. In this section, we rst document the types of branch elements used in power ow analysis and then develop a general expression for the entries of Y that satises (4.4). Example 4.1. Figure 4.1 shows an example network consisting of two buses i and k and a single branch (i, k ) between them. Branch (i, k ) has series admittance
14
Rik Bus i bik 2

Sh
Xik Bus k
Sh
bik 2
Fig. 4.2. branch model for cables and transmission lines.
yik , and each bus also has a shunt admittance. For this network, writing (4.3) for each bus yields the matrix equation Ii Ik
311
S yik + yi yik
yik S yik + yk
Vi Vk
4.1.1. Cables and Transmission Lines. Power cables and transmission lines are modeled as branch circuits (Figure 4.2). The line characteristics are specied as a series impedance Rik + jXik and a branch shunt admittance jbSh ik , which is sometimes given as line charging reactive power. The branch series admittance yik for inclusion in Y is yik = gik bik
312 313 314 315
Xik Rik 1 = 2 2 j R2 + X 2 , Rik + jXik Rik + Xik ik ik Rik = 2 2 , Rik + Xik Rik = 2 2 . Rik + Xik
(4.5)
Branch shunt susceptance bSh ik is related to but distinct from the net shunt susceptance at buses i and k . Specically, the branch shunt susceptance bSh ik is divided into two S at each end of the line, as equal parts and added to bus shunt susceptances bS and b i k illustrated in Figure 4.2. For short lines, branch shunt susceptance is usually omitted. 4.1.2. Transformers. Most power systems transformers have nominal turns ratios, that is, the voltage ratio across the transformer exactly equals change in system voltage base across the transformer (a 1:1 voltage ratio in per-unit, with no phase shift). Because the per-unit system automatically accounts for the turns ratio, the branch model for a transformer with nominal turns ratios is identical to the branch circuit for a transmission line (Figure 4.2). However, in power ow analysis, transformer branches are almost always modeled with zero shunt susceptance and often with zero series resistance as well. 4.1.3. O-Nominal Transformers. Any transformer that does not have exactly a 1:1 voltage ratio in per-unit is an o-nominal transformer. This category includes xed-tap transformers with o-nominal turns ratios, tap-changing transformers, and phase-shifting transformers. O-nominal transformers require modied entries in Y to account for the additional voltage magnitude or phase angle change relative to the nominal case. Proper modeling of o-nominal transformers is a key
316 317 318 319 320 321 322 323
324 325 326 327 328 329
15
Bus i ei V
e I ik
eik I' ei V' a:1
yik e
e I ki
Bus k
ek V
Fig. 4.3. O-nominal transformer branch model.
330 331
skill in developing practical OPF algorithms. Unfortunately, this topic is neglected in most introductory OPF texts. Figure 4.3 displays the standard model for o-nominal transformers found in practical power ow and OPF software. In the model, bus i is the tap bus and bus k is the impedance bus or Z bus. The transformer turns ratio in per-unit is a:1, where a is a complex exponential consisting of magnitude T and phase shift , a = T ej ,
332 333 334 335
such that Vi = aVi and Iik = Iik /a . Selecting T = 1 and = 0 yields the nominal turns ratio. In OPF, either T or (or both) may be a control variable: controllable T models an on-load tap changer, while controllable models a phase shifting transformer. In order to include the eects of o-nominal transformer (i, k ) in Y , partial admittance matrix entries for branch (i, k ) are required such that Iik Iki = Yii Yki Yik Ykk Vi Vk .
Using the turns ratio denitions given above and Ohms law, the expression for current Iik is developed as follows: Iik = 1 I , a ik 1 = yik Vi Vk , a 1 1 = yik Vi Vk , a a 1 1 = yik Vi yik Vk . aa a
(4.6)
Similarly, Iki = yik Vk Vi , 1 = yik Vk Vi , a 1 = yik Vi + yik Vk . a
(4.7)
16
In matrix form, expressions (4.6) and (4.7) become 1 1 y y ik Iik Vi ik a = aa1 Iki Vk yik yik a
336 337 338 339 340 341 342
(4.8)
The diering expressions for Iik and Iki in (4.8) indicate the importance of the dierence between the tap bus i and the Z bus k ; reversing the two leads to considerable error. When constructing the full admittance matrix Y , the relationships of (4.8) must be preserved. If Y has previously been constructed according to (4.4) with o-nominal voltage ratios neglected, then the following correction procedure is required for each o-nominal branch (i, k ): 1. The partial diagonal term corresponding to branch (i, k ) in Yii is divided by 2 aa = |a| , as given by the replacement procedure
Old Yii Yii yik +
1 yik , aa
343 344 345
where yik is the uncorrected partial diagonal admittance term. 2. The partial diagonal admittance term corresponding to branch (i, k ) in Ykk remains unchanged. 3. The o-diagonal admittance matrix entry Yik is divided by a , as given by the replacement procedure Yik 1 yik . a
4. The o-diagonal admittance matrix entry Yki is divided by a, as given by the replacement procedure 1 Yki yik . a
346 347 348 349 350 351 352
When the procedure is complete, the eects of the o-nominal turns ratios are included directly in Y . For notational convenience, this correction procedure is written with respect to an already constructed admittance matrix Y . In practice, the corrections are made in the initial construction of Y , as given in 4.1.4, rather than performed as a replacement procedure. Remark 4.2. A transformer with o-nominal magnitude only (real valued a) leaves Y a symmetric matrix, but a phase-shifting transformer (complex a) does not. 4.1.4. Construction Equations for Admittance Matrix. In general, any of the branch elements described above can be represented by a series admittance Sh yik , a shunt admittance yik , and a complex turns ratio (nominal or o-nominal) jik aik = Tik e . Using (4.4) and the correction procedure given in 4.1.3, the entries of Y become
S Yii = yi +
1 |aik | k:(i,k)L
2
1 Sh yik + yik 2
+
k:(k,i)L
1 Sh yki + yki , 2
(4.9) (4.10)
Yik =
k:(i,k)L
1 yik a ik
k:(k,i)L
1 yki , aik
i=k
17
Table 4.1 Branch impedance data for Example 4.3. All quantities except phase angles are given in perunit. Dots indicate nominal voltage ratios and phase angles. Series Resistance Rik 0.000 0.023 0.006 0.020 0.000 0.000 Series Reactance Xik 0.300 0.145 0.032 0.260 0.320 0.500 Shunt Susceptance bSh ik 0.000 0.040 0.010 0.000 0.000 0.000 Voltage Ratio Tik 0.98 Phase Angle ik 3.0
From Bus i 1 1 2 3 3 4
To Bus k 2 3 4 4 5 5
where a = 1 for any branch with a nominal turns ratio. Equations (4.9)(4.10) can also be separated until real and imaginary parts using the denition Y = G + jB and the identity aik = Tik (cos ik + j sin ik ),
S Gii = gi + k:(i,k)L
1 2 Tik
1 Sh gik + gik 2
+
k:(k,i)L
1 Sh gki + gki , 2
(4.11)
Gik =
k:(i,k)L
1 (gik cos ik bik sin ik ) Tik
k:(k,i)L
1 (gki cos ki + bki sin ki ) , Tki 1 2 Tik 1 bik + bSh 2 ik +

k:(k,i)L
i=k 1 bki + bSh , 2 ki
(4.12) (4.13)
Bii = bS i +
k:(i,k)L
Bik =
k:(i,k)L
1 (gik sin ik + bik cos ik ) Tik
353 354 355 356
k:(k,i)L
1 (gki sin ki + bki cos ki ) Tki
i = k.
(4.14)
Example 4.3. Table 4.1 provides a set of branch data for the 5-bus example system of Figure 3.1. Note that branch (3, 4) is a phase-shifting transformer and branch (4, 5) has an o-nominal voltage ratio. In addition to the branch data, bus 2 has a shunt susceptance of j 0.30 pu and bus 3 has a shunt conductance of 0.05 pu. To compute the admittance matrix for this system, we rst compute the series admittance yik of each branch using (4.5). For example, the series admittance of branch (1, 3) is y13 = 0.023 0.145 j 1.067 j 6.727. 0.0232 + 0.1452 0.0232 + 0.1452 y12 0.000 j 3.333,
The remaining branches have series admittances y24 5.660 j 30.189,
y35 0.000 j 3.125,
y34 0.294 j 3.824,
18 and
y45 0.000 j 2.000.

357
Verication of these values is left as an exercise for the reader. Next, we construct Y using (4.9)(4.10). For example, diagonal element Y33 consists of summing the admittances of branches (1, 3), (3, 4), and (3, 5), plus the contributions of the shunt conductance at bus 3 and the shunt susceptance of branch (1, 3). y34 and y35 have o-nominal turns ratios a34 = 1.0ej 3.0 0.999 j 0.052 and a35 = 0.98ej 0.0 = 0.980. (Note that a34 a 34 = 1.0 if rounding errors are neglected.) Therefore, the full expression for Y33 is Y33 (1.067 j 6.727) + j 1.41 j 13.78. 0.294 j 3.824 j 3.125 0.04 + + 0.05, 2 (0.999 j 0.052) (0.999 + j 0.052) 0.9802
An example o-diagonal element is Y34 , which from (4.10) is Y34
0.294 j 3.824 0.09 + j 3.83. (0.999 + j 0.052)

0 0 0.00 + j 3.19 . 0.00 + j 2.00 0.00 j 5.13
The full admittance matrix is

1.07 j 10.04 0.00 + j 3.33 Y 1.07 + j 6.73 0 0 0.00 + j 3.33 5.66 j 33.22 0 5.66 + j 30.19 0 1.07 + j 6.73 0 1.41 j 13.78 0.49 + j 3.80 0.00 + j 3.19 0 5.66 + j 30.19 0.09 + j 3.83 5.95 j 36.01 0.00 + j 2.00
358 359 360 361 362 363 364 365
Verication of the remaining matrix elements is left as an exercise for the reader.
4.2. AC Power Flow. For a given Y , equation (4.2) can be decomposed into a set of equations for the real and reactive power injections by evaluating the real and imaginary parts of S , respectively. This process yields a pair of equations called the AC power ow equations, which can be written in several equivalent forms depending on whether the voltages and admittance matrix elements are expressed in polar or rectangular coordinates. In the literature, the most common forms of the AC power ow equations are (in order) 1. Selection of polar coordinates for voltage, Vi = Vi i , and rectangular coordinates for admittance, Yik = Gik + jBik :
N
Pi (V, ) = Vi
k=1 N
Vk Gik cos (i k ) + Bik sin (i k ) Vk Gik sin (i k ) Bik cos (i k )
(4.15)
Qi (V, ) = Vi
k=1
(4.16)
19
2. Selection of polar coordinates for voltage, Vi = Vi i , and polar coordinates for admittance, Yik = Yik ik :
N
Pi (V, ) = Vi
k=1 N
Vk Yik cos (i k ik ) Vk Yik sin (i k ik )
(4.17)
Qi (V, ) = Vi
k=1
(4.18)
3. Selection of rectangular coordinates for voltage, Vi = Ei + jFi , and rectangular coordinates for admittance, Yik = Gik + jBik :
N
Pi (E, F ) =
k=1 N
Gik (Ei Ek + Fi Fk ) + Bik (Fi Ek Ei Fk ) Gik (Fi Ek Ei Fk ) Bik (Ei Ek + Fi Fk )
(4.19)
Qi (E, F ) =
k=1
366 367 368 369 370 371 372 373 374 375 376 377
(4.20)
Power systems texts [14, 32, 37] provide exact derivations of these three forms of the AC power ow equations.3 Each form of the equations involves real-valued quantities only. However, all forms are equivalent and give the exact solution to the power ow under the assumptions outlined in 3.1. From an OPF perspective, there is little dierence between the selection of polar or rectangular coordinates for the admittance matrix. Rectangular coordinates are more common in practice because they facilitate the use of certain approximations in fast-decoupled solution methods for conventional PF [37]. These approximations are also useful in the development of the DC power ow equations; see 4.3. Rectangular coordinates also facilitate the inclusion of transformer voltage ratios and phase angles as decision variables. However, neither of these advantages strongly aects the AC power ow equations as used in most OPF formulations. The more important distinction is the choice of polar or rectangular coordinates for voltage. The advantage of voltage polar coordinates is that constraints on the voltage magnitude can be enforced directly, Vi Vimin , Vi Vimax .
In voltage rectangular coordinates, on the other hand, voltage magnitude limits require the functional inequality constraints
2 + F 2 V min , Ei i i 2 + F 2 V max . Ei i i
Similarly, if the voltage magnitude is xed (for instance at a PV bus; see 4.4), then in polar coordinates Vi can be replaced with a constant value. In rectangular coordinates, however, a xed voltage magnitude requires the equality constraint
2 + F2 = V . Ei i i
3 A fourth formselection of rectangular coordinates for voltage and polar coordinates for admittanceis theoretically possible but has no advantages for practical use.
20
Table 4.2 Comparison of the selection of voltage polar coordinates versus voltage rectangular coordinates for the power ow equations. Bold entries indicate the superior characteristic. Polar Coords. Variable limit Variable elimination by substitution N +M 1 Trigonometric Trigonometric Trigonometric Rectangular Coords. Nonlinear functional inequality constraint Nonlinear functional equality constraint 2N 2 Quadratic Linear Constant
Voltage magnitude limit Fixed voltage magnitude # of variables in conventional PF Nature of PF equations 1st derivative of PF equations 2nd derivative of PF equations
378 379 380 381
Thus, for a xed voltage magnitude, use of voltage polar coordinates leads to a reduction of variables while use of voltage rectangular coordinates leads to an increase in (non-linear, non-convex) equality constraints. For this reason, polar coordinates are preferred both for conventional PF and most OPF formulations. There is, however, one compelling reason to use voltage rectangular coordinates: expressing voltage in rectangular coordinates eliminates trigonometric functions from the power ow equations. The resulting power ow equations (4.19)(4.20) are quadratic, which presents several advantages [30]: 1. The elimination of trigonometric functions speeds evaluation of the equations. 2. The 2nd order Taylor series expansion of a quadratic function is exact; this yields an eciency advantage in higher-order interior-point algorithms for OPF. 3. The Hessian matrix for a quadratic function is constant and need be evaluated only once. This simplies the application of Newtons method to the KKT conditions of the OPF formulation. In some cases, these computational advantages outweigh the disadvantages associated with enforcing voltage magnitude constraints. Table 4.2 summarizes the dierences between the two voltage coordinate choices. Example 4.4. Using the admittance matrix developed in Example 4.3, we can write the real and reactive power ow equations for any bus in the 5-bus example system. From (4.15), the real power injection at bus 1 is
382 383 384 385
386 387 388 389 390 391
392 393 394
P1 (V, ) = V1
k=1 1.07V12
Vk G1k cos (1 k ) + B1k sin (1 k ) , cos (1 1 ) 1.07V1 V3 cos (1 3 ) 10.04V12 sin (1 1 )
1.07V12 1.07V1 V3 cos (1 3 ) + 3.33V1 V2 sin (1 2 ) + 6.73V1 V3 sin (1 3 ).
+ 3.33V1 V2 sin (1 2 ) + 6.73V1 V3 sin (1 3 ),
21
Similarly, from (4.16), the reactive power injection at bus 1 is

5
Q1 (V, ) = V1
k=1 1.07V12
Vk G1k sin (1 k ) B1k cos (1 k ) , sin (1 1 ) 1.07V1 V3 sin (1 3 ) + 10.04V12 cos (1 1 )
1.07V1 V3 sin (1 3 ) + 10.04V12 3.33V1 V2 cos (1 2 ) 6.73V1 V3 cos (1 3 ).

395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411
3.33V1 V2 cos (1 2 ) 6.73V1 V3 cos (1 3 ),
Evaluation of the remaining buses is left as an excercise for the reader. 4.3. DC Power Flow. The AC power ow equations are nonlinear. For conventional PF, this nonlinearity requires the use of an iterative numerical method; for OPF it implies both a nonlinear formulation and non-convexity in the feasible region. In order to simplify the system representation, power systems engineers have developed a linear approximation to the power ow equations. This approximation is called DC power ow.4 The conventional development of the DC power ow equations requires several assumptions regarding the power system [26, 37]: 1. All system branch resistances are approximately zero, that is, the transmission system is assumed to be lossless. As a result, all ik = 90 and all Gik = 0. 2. The dierences between adjacent bus voltage angles are small, such that sin(i k ) i k and cos(i k ) 1. 3. The system bus voltages are approximately equal to 1.0. This assumption requires that there is sucient reactive power generation in the system to maintain a level voltage prole. 4. Reactive power ow is neglected. Applying these assumptions to (4.15) produces the DC power ow equation
N
Pi ( )
412 413 414 415 416
k=1
Bik (i k )
(4.21)
Under normal system operating conditions, DC power ow models real power transfer quite accurately. It has been successfully used in many OPF applications that require rapid and robust solutions. However, the assumptions required for DC power ow can lead to signicant errors for stressed systems. The exact equation for branch power transfer is Pik = gik Vi2 gik Vi Vk cos (i k ) bik Vi Vk sin (i k ), (4.22)
417
cf. (4.15), while the DC power ow approximation is Pik bik sin (i k ). (4.23)
418 419 420 421
The bik term dominates the exact expression because Vi2 Vi Vk cos (i k ) and therefore the rst two terms in (4.22) largely cancel. We observe that (4.23) overestimates the magnitude of the branch power transfer (4.22) if
4 The DC power ow is so named because the equations resemble the power ow in a direct current (DC) network. However, the DC power ow equations still model an AC power system.
22
422 423 424 425 426 427 428 429 430 431 432 433 434 435
(i) The bus voltages at either end of the branch are depressed relative to the assumed value of 1.0 p.u., or (ii) The angle dierence between the buses is too large. Observation (ii) follows from the relationship |sin (i k )| |i k |. Depressed voltages and larger than normal angle dierences are common in stressed power systems. In particular, large dierences in voltage in dierent areas of the system can lead to signicant error [28]. Therefore, the DC power ow equations should not be used for OPF under stressed system conditions unless they have been carefully evaluated for accuracy in the system under test. Example 4.5. Consider a transmission line from bus i to bus k with admittance 0.05 j 2.0. Let Vi = 0.95 0 and Vk = 0.90 20. (These numbers do not represent normal operation, but are plausible for a stressed power system. Operating voltages as low as 0.9 p.u. are allowable in emergency conditions, and angle dierences of up to 30 can occur on long, heavily loaded transmission lines.) The exact power transfer for this line is Pik = 0.05 0.952 0.05 0.95 0.90 cos (0 + 20) 2.0 0.95 0.90 sin (0 + 20), = 0.590 p.u. The approximate power transfer is Pik 2.0 sin (0 + 20), 0.684 p.u.
436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461
The error in the approximate power transfer is 16%; most of this error is attributable to the voltage dierence. Even under normal operation, the approximation of a lossless transmission network can also lead to signicant errors in generator scheduling, branch power ow estimates, and marginal fuel cost estimates. Power transfer errors for certain critically loaded branches can be much higher than the average branch error. Therefore, in practical DC power ow models an estimate of the losses must be reintroduced using approximate methods, especially if the network is large [28]. For further discussion regarding the advantages and disadvantages of DC power ow, including loss approximation methods, we refer the interested reader to [25, 28]. Throughout the rest of this primer, we use the AC power ow equations. 4.4. Solution Methods for Conventional PF. Many practical OPF algorithms incorporate aspects of conventional PF solution methods. Therefore, a basic understanding of these methods is helpful when reviewing OPF literature. Here, we discuss the solution of the AC power ow equations with voltage polar coordinates. The solution method for voltage rectangular coordinates is similar; Zhu [37] provides a good summary. Each system bus has four variables (real power injection Pi , reactive power injection Qi , voltage magnitude Vi , and voltage angle i ) and is governed by two equations: either (4.15)(4.16) or (4.17)(4.18). Thus, a unique solution to the conventional PF requires xing the values of two out of four variables at each bus. Remark 4.6. Even though the power ow equations are nonlinear, there exists only one physically meaningful solution for most power systems models given an equal number of equations and unknowns. Other solutions may exist mathematically, but have no realistic physical interpretation. (An example would be any solution which returns a negative voltage magnitude, as magnitudes are by denition nonnegative.)
A PRIMER ON OPTIMAL POWER FLOW Table 4.3 Power system bus types and characteristics for conventional power ow. Bus Type # of buses in system Known quantities Unknown quantities # of equations in conventional PF Slack 1 , V P, Q 0 PQ M P, Q , V 2 PV N M 1 P, V , Q 1
23
462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480
In conventional PF, all system buses are assigned to one of three bus types: Slack Bus At the slack bus, the voltage magnitude and angle are xed and the power injections are free. There is only one slack bus in a power system. Load Bus At a load bus, or PQ bus, the power injections are xed while the voltage magnitude and angle are free. There are M PQ buses in the system. Voltage-Controlled Bus At a voltage controlled bus, or PV bus, the real power injection and voltage magnitude are xed while the reactive power injection and the voltage angle are free. (This corresponds to allowing a local source of reactive power to regulate the voltage to a desired setpoint.) There are N M 1 PV buses in the system. Assigning buses in this way establishes an equal number of equations and unknowns. Table 4.3 summarizes the three bus types. If all voltage magnitudes and angles in the system are known, then the power injections are fully determined. Solving the power ow therefore requires determining N 1 voltage angles (corresponding to the PQ and PV buses) and M voltage magnitudes (corresponding to the PQ buses only). This is done by solving N + M 1 simultaneous nonlinear equations with known right hand side values. This equation set consists of the real power injection equation (4.15) at each PV and PQ bus and the reactive power injection equation (4.16) at each PQ bus. Newtons method is commonly used to solve this system. The 1st order Taylor series approximation about the current estimate of V and yields P P V P Q Q V , Q V J , V (4.24)
where J is the Jacobian matrix of the system. At each iteration, the mismatches in the power ow equations are Pi = PiG PiL Pi (V, ) , QG i QL i (4.25) (4.26)
Qi =
481 482 483 484 485 486
Qi (V, ) .
Newtons method consists of iteratively solving (4.24) for the and V required to correct the mismatch in the power ow equations computed from (4.25)(4.26). Newtons method is locally quadratically convergent. Therefore, given a suciently good starting point, the method reliably nds the correct solution to the PF equations. Newtons method for conventional PF is described more fully in power systems texts [14, 32, 37].
24
487 488 489 490 491 492 493 494 495 496 497 498
In OPF, the decision variables are often partitioned into a set of control (independent) variables u and a set of state (dependent) variables x [7, 10]. At each search step, the OPF algorithm xes u and derives x by solving a conventional PF. When this method is used, the Jacobian matrix J plays several important roles: 1. It provides the linearization of the power ow equations required for successive linear programming (SLP) OPF algorithms, 2. It provides sensitivities in the power ow injections with respect to the state variables, 3. It provides a direct calculation of portions of the Hessian matrix of the Lagrangian function in OPF (see [29]), and 4. It is therefore often used to improve computational eciency in computing the KKT conditions for the Lagrangian function. In DC power ow, there is no distinction between PV and PQ buses because all voltage magnitudes are considered to be 1.0. As with the AC power ow, the slack bus angle is xed. Because the DC power ow equations are linear, they may be solved directly for the voltage angles using = B 1 P.
499
(This is simply a solved matrix representation of (4.21).) 4.5. Practical Considerations. In our experience, there are several practical and computational aspects of OPF stemming from the power ow equations that can cause confusion. One of these is the use of the per-unit system, which is discussed in 3.4. We discuss a few others here. 4.5.1. Degrees versus Radians. Power systems engineers usually report angles in degrees, including in data les for OPF (see 6). For computation, these angles must be converted to radians, for two reasons: 1. Nearly all optimization software and descriptive languagesincluding AMPL and GAMSimplement trigonometric functions in radians, not degrees. 2. Even when using the DC power ow equations (which require no trigonometric function evaluations), radians must be used. If degrees are used, the powers computed from DC power ow will have a scaling error of 180/ . Power ow software typically handles these conversions transparently, accepting input and giving output in degrees. Thus, it can be dicult to remember that generalpurpose optimization software requires an explicit conversion. 4.5.2. System Initialization. In both conventional power ow and OPF, the convergence of the power ow equations depends strongly on the selection of a starting point. Given a starting point far from the correct solution, the power ow equations may converge to a meaningless solution, or may not converge at all. In the absence of a starting point, standard practice is to initialize all voltage magnitudes to 1.0 p.u. and all voltage angles to zero; this is called a cold start or a at start. The alternative is a hot start, in which the voltages and angles are initialized to the solution of a pre-solved power ow. Hot starts are often used in online OPF to minimize computation time and ensure that the search begins from the current system operating condition. 4.5.3. Decoupled Power Flow versus Decoupled OPF. In practical power systems, real power injections are strongly coupled to voltage angles and reactive power injections are strongly coupled to voltage magnitudes. Conversely, real power
500 501 502 503
504 505 506 507 508 509 510 511 512 513 514
515 516 517 518 519 520 521 522 523 524
25
injections are weakly coupled to voltage magnitudes and reactive power injections are weakly coupled to voltage angles. This feature has led to the development of decoupled solution methods for the power ow equations [14, 37]. The most basic decoupling method is to use a set of approximate Taylor series expansions of the form P , Q Q V. V P
525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546
This allows the use of separate Newton updates for and V with correspondingly smaller matrices; this is a signicant computational advantage. Although decoupled power ow uses an approximate update method, it still uses exact real and reactive power mismatches P and Q from (4.25)(4.26) and updates both V and at each iteration. Decoupled power ow therefore is locally convergent to the exact solution to the power ow. However, because of the approximated Jacobian matrix, more iterations are required for convergence [14]. Zhu [37] discusses several decoupled power ow variants in detail. Decoupled OPF also takes advantage of the strong P - and Q-V relationships by formulating a real subproblem and a reactive subproblem. The optima of the subproblems are assumed to be independent. Unlike decoupled power ow, however, decoupled OPF solves the subproblems sequentially rather than simultaneously: the real subproblem solves for the optimal values of P and while holding Q and V constant, and the reactive subproblem solves for the optimal values of Q and V while holding P and constant [9,29]. Decoupled OPF is therefore distinctly dierent from decoupled power ow in that the decoupled OPF solution is inexact. The error is a function of the accuracy of the decoupling assumptions; these assumptions should therefore be evaluated for accuracy if a decoupled OPF approach is considered. Remark 4.7. In the OPF literature, it is not always clear whether decoupled OPF is in use or whether a decoupled power ow procedure is used within the solution algorithm for a coupled OPF. Because of the implications for the OPF solution quality, the careful reader should try to discern which is the case. 5. Optimal Power Flow. Broadly speaking, any power systems optimization problem which includes the power ow equations in the set of equality constraints is an OPF problem. Thus, the term OPF now encompasses an extremely wide variety of formulations, many with tailored solution methods [11]. Most of these variants, however, build upon the classic formulation of Carpentier [8] and Dommel and Tinney [10]. (This is so common that most OPF papers omit the core of the formulation entirely, focusing only on novel enhancements or algorithmic development.) Here, we rst present the classical formulation and then briey discuss several common extensions. 5.1. Classical Formulation. The classical OPF formulation of Dommel and Tinney is an extension of economic dispatch (ED): its objective is to minimize the total cost of electricity generation while maintaining the electric power system within safe operating limits. The power system is modeled as a set of buses N connected by a set of branches L. Controllable generators are located at a subset G of the system buses. The operating cost of each generator is a (typically quadratic) function of its real output power: Ci PiG . The objective is to minimize the total cost of generation.
547 548 549 550 551 552 553 554 555
556 557 558 559 560 561 562
26
The classical form of the formulation is min

i G
Ci PiG , Pi (V, ) = PiG PiL QG i i N, i G, i N,
(5.1) (5.2) (5.3) (5.4) (5.5) (5.6) (5.7)
s.t.
Qi (V, ) = QL i G,min G,max G Pi P i Pi G,min G,max QG Qi i Qi Vimin Vi Vimax min max i i i
i N,
i G, i N.
In (5.2)(5.3), Pi (V, ) and Qi (V, ) represent the power ow equations in polar formeither (4.15)(4.16) or (4.17)(4.18). The vector of control variables (independent decision variables) is
G u = PiG :iG , Qi:iG
and the vector of state variables (dependent decision variables) is x = (2 , . . . , N , V2 , . . . , VN ) .

563 564
The voltage magnitude and angle at the system slack bus (by convention, bus 1) are xed, usually to V1 = 1.00. If the system contains controllable phase-shifting or tap-changing transformers, then the corresponding phase angles and tap ratios are introduced into the set of control variables. The control variable vector u becomes
G u = PiG :iG , Qi:iG , ik:ikH , Tik:ikK ,
where H and K are the sets of branches with controllable-phase shifting transformers and tap-changing transformers, respectively. Since and T alter the elements of admittance matrix Y , the left hand sides of (5.2) and (5.3) become functions of and T : Pi (V, , , T ) and Qi (V, , , T ), respectively. The formulation is also augmented with bound constraints on the phase angles
max min ik ik ik
ik H
(5.8)
and the tap ratios

min max Tik Tik Tik
ik K.
(5.9)
Although not considered in the earliest papers, more recent OPF formulations also consider branch current limits. Unlike the previous bounds, the branch current limits require functional inequality constraints. By Ohms law, the current magnitude in branch ik is Iik = Vi Vk yik , (5.10)
A PRIMER ON OPTIMAL POWER FLOW Table 5.1 Bus data for Example 5.2. All quantities are given in per-unit. Dots indicates zero values. Load Real Power PiL 0.900 0.239 Load Reactive Power QL i 0.400 0.129 Min. Bus Voltage Vimin 1.00 0.95 0.95 0.95 0.95 Max. Bus Voltage Vimax 1.00 1.05 1.05 1.05 1.05
27
Bus i 1 2 3 4 5
Table 5.2 Generator data for Example 5.2. All quantities are given in per-unit. Min. Generator Real Power PiG,min 0.10 0.05 Max. Generator Real Power PiG,max 0.40 0.40 Min. Generator Reactive Power ,min QG i 0.20 0.20 Max. Generator Reactive Power ,max QG i 0.30 0.20
Bus i 1 3 4
where yik is the magnitude of the branch admittance. Thus, we can constrain the branch current using
max Vi Vk yik Iik ,

565 566 567 568 569 570 571 572 573 574 575 576
(Vi cos i Vk cos k ) + (Vi sin i Vk sin k ) (Vi cos i Vk cos k ) + (Vi sin i Vk sin k )
2 2
max Iik , yik max (Iik ) 2 yik 2
ik L. (5.11)
Rather than bounding the square of the current as is given in (5.11), many formulations bound the total real and reactive power ow entering the line. However, (5.11) gives a more exact representation of the true constraint, which is technically a maximum current, not a maximum power. Remark 5.1. For branches with o-nominal turns ratios, the tap bus voltage Vi and angle i in (5.11) must be corrected for the o-nominal turns ratio. In this case, Vi is replaced by Vi = Vi /Tik and i is replaced by i = i ik . Even without variable phase angles, variable tap ratios, or branch current limits, the classical OPF formulation is dicult to solve. The power ow constraints (5.2) (5.3) are both nonlinear and non-convex, and the presence of trigonometric functions complicates the construction of approximations. For this reason, OPF problems have historically been solved using tailored algorithms rather than general purpose solvers. Example 5.2. We now develop the classical OPF formulation for the 5-bus example system rst presented in Figure 3.1. For this example, the branch impedance data are as given in Table 4.1, except that we assign 34 and T35 to be decision variables representing a phase shifting transformer and an on-load tap changer, respectively. 34 and T35 have limits 30.0 34 30.0
28 and
0.95 T35 1.05. Consider the bus data (voltage limits, load, and generation) given in Tables 5.1 and 5.2. The system power base is 100 MW. Given this data, the sets dening the formulation are: N = {1, 2, 3, 4, 5} , L = {(1, 2), (1, 3), (2, 4), (3, 4), (3, 5), (4, 5)} ,
G = {1, 3, 4} , H = {(3, 4)} , and K = {(3, 5)} .
The three generator cost functions, in thousands of dollars, are

G G C1 P1 = 0.35P1 , G G G C3 P3 = 0.20P3 + 0.40 P3 G G C4 P4 = 0.30P4 + 0.50
577
, ,
G 2 P4
where the PiG are expressed in per-unit. To develop the full formulation, it is rst necessary to re-write Y from Example 4.3 to explicitly include 34 and T35 . Let a34 = cos 34 + j sin 34 and a35 = T35 .
(Note that a34 a 34 = 1.0, 1/a34 = cos 34 j sin 34 = a34 , and 1/a34 = cos 34 + j sin 34 = a34 .) Then, using (4.9)(4.10) and simplifying,
Y33 = 1.41 j 10.53 j
1 2 3.13, T35
Y43 = 0.29 cos 34 + 3.82 sin 34 + j (3.82 cos 34 + 0.29 sin 34 ) , 1 Y35 = j 3.13, T35 and Y53 = j
578
Y34 = 0.29 cos 34 3.82 sin 34 + j (3.82 cos 34 0.29 sin 34 ) ,
1 3.13. T35
Y44 , Y55 , and the remaining matrix entries are unchanged.
29
Bus 1 is the system slack bus, and therefore V1 is xed to 1.00.0. To construct the formulation, we round all numerical values to two decimal places. (This rounding does not aect model feasibility because sucient degrees of freedom exist in the state variables.) Following (5.1)(5.7), (5.8), and (5.9), the full formulation is
min s.t.
G G G 0.35P1 + 0.20P3 + 0.40 P3 G P1
G G + 0.30P4 + 0.50 P4
0 = 5.66V22 5.66V2 V4 cos (2 4 )

G P3 = 1.41V32 1.07V3 cos (3 )
= 1.07 1.07V3 cos (3 ) + 3.33V2 sin (2 ) + 6.73V3 sin (3 ) , + 3.33V2 sin (2 ) + 30.19V2 V4 sin (2 4 ) ,
+ 6.73V3 sin (3 ) + (3.82 cos 34 0.29 sin 34 ) V3 V4 sin (3 4 ) 3.13 V3 V5 sin (3 5 ) , + T35 G P4 0.900 = 5.95V42 5.66V4 V2 cos (4 2 ) + 30.19V4 V2 sin (4 2 ) + (0.29 cos 34 + 3.82 sin 34 ) V4 V3 cos (4 3 )
+ (0.29 cos 34 3.82 sin 34 ) V3 V4 cos (3 4 )
+ (3.82 cos 34 + 0.29 sin 34 ) V4 V3 sin (4 3 ) + 2.00V4 V5 sin (4 5 ) , 3.13 V5 V3 sin (5 3 ) + 2.00V5 V4 sin (5 4 ) , 0.239 = T35 QG 1 = 10.04 1.07V3 sin (3 ) 3.33V2 cos (2 ) 6.73V3 cos (3 ) , 0 = 5.66V2 V4 sin (2 4 ) + 33.22V22 3.33V2 cos (2 ) 30.19V2 V4 cos (2 4 ) , 3.13 2 T35 V32 6.73V3 cos (3 ) 3.13 V3 V5 cos (3 5 ) , T35
QG 3
= 1.07V3 sin (3 ) + (0.29 cos 34 3.82 sin 34 ) V3 V4 sin (3 4 ) + 10.53 +
(3.82 cos 34 0.29 sin 34 ) V3 V4 cos (3 4 ) QG 4 0.940 = 5.66V4 V2 sin (4 2 ) + (0.29 cos 34 + 3.82 sin 34 ) V4 V3 sin (4 3 )
(3.82 cos 34 + 0.29 sin 34 ) V4 V3 cos (4 3 ) 2.00V4 V5 cos (4 5 ) , 3.13 0.129 = 5.13V52 V5 V3 cos (5 3 ) 2.00V5 V4 cos (5 4 ) , T35
+ 36.01V42 30.19V4 V2 cos (4 2 )
30

G 0.10 P3 0.40,
180.0 i 180.0,
579 580 581
0.95 Vi 1.05,
0.95 T35 1.05,
30.0 34 30.0, i {2, 3, 4, 5} ,
0.20 QG 4 0.20,
0.20 QG 3 0.30,
G 0.05 P4 0.40,
i {2, 3, 4, 5} .
Voltage angles 1 , 2 , 3 , and 4 are restricted to one full sweep of the unit circle. G Slack bus generator powers P1 and QG 1 are unrestricted, and branch current limits are neglected. For this formulation, the vector of control variables is
G G G G G u = P1 , P3 , P4 , QG 1 , Q3 , Q4 , 34 , T35
and the vector of state variables is x = (2 , 3 , 4 , 5 , V2 , V3 , V4 , V5 ) . The optimal solution for this formulation is V2 0.981,
G P1 QG 1
2 12.59, 0.947, 0.387,
V3 0.957,
G P3 QG 3
3 1.67, 0.192, 0.127,
V4 0.968,
G P4 QG 4
4 13.86, 0.053, 0.200,
V5 0.959,
5 9.13,
34 12.38,
T35 0.95,
with objective function value 0.4016596. If the controllable phase-shifting and tapchanging transformers are instead xed to 34 = 3.0 and T35 = 0.98, the optimal solution becomes V2 0.983,
G P1 QG 1
582 583 584 585 586 587 588 589
2 7.50, 0.946, 0.249,
V3 0.964,
G P3 QG 3
3 4.22, 0.195, 0.072,
V4 0.970,
G P4 QG 4
4 8.20, 0.058, 0.200,
V5 0.950,
5 8.64,
with objective function value 0.4041438, a cost increase of approximately 0.6%. To obtain the optimal solution for Example 5.2, we implemented three versions of the classical formulation (5.1)(5.9) in the GAMS modeling language. The three versions each use a dierent form of the power ow equations: (i) polar voltage coordinates with rectangular admittance coordinates (4.15)(4.16), (ii) polar voltage coordinates with polar admittance coordinates (4.17)(4.18), and (iii) rectangular voltage coordinates with rectangular admittance coordinates (4.19)(4.20). The model is publicly available in the GAMS model library [1].5 For the example, the model yielded
5 Note to the reviewer: the GAMS model is attached at the end of the article. The model is not yet available for download in the GAMS model library as stated here, but will be nalized and made available when this paper is published.

590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606
31
identical optimal solutions using three local nonlinear solvers, SNOPT, MINOS, and CONOPT, and veried as globally optimal using the global solver LINDOGlobal. 5.2. Special Applications. Besides the classical ED formulation, several other OPF variants are common in both industry and research. These include securityconstrained economic dispatch (SCED), security-constrained unit commitment (SCUC), optimal reactive power ow (ORPF), and reactive power planning (RPP). 5.2.1. Security-Constrained Economic Dispatch. Security-constrained economic dispatch (SCED), sometimes referred to as security-contrained optimal power ow (SCOPF), is an OPF formulation which includes power system contingency constraints [4]. A contingency is dened as an event which removes one or more generators or transmission lines from the power system, increasing the stress on the remaining network. SCED seeks an optimal solution that remains feasible under any of a pre-specied set of likely contingency events. SCED formulations typically have the same objective function and decision variables u as the classic formulation.6 However, they introduce NC additional sets of state variables x and accompanying sets of power ow constraints, where NC is the number of contingencies. SCED can be expressed in a general way as min s.t. f (u, x0 ), g0 (u, x0 ) = 0, h0 (u, x0 ) 0, gc (u, xc ) = 0 hc (u, xc ) 0 (5.12) c C, c C,
607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628
where C = {1, . . . , NC } is the set of contingencies to consider. Each contingency has a distinct admittance matrix Yc with less connectivity than the original system. Apart from the contingency index, f , g , and h are dened as the objective function, equality contstraints, and inequality constraints in the classical OPF formulation of 5.1, respectively. In other words, for each contingency c C, the post-contigency power ow must remain feasible for the original decision variables u: (i) The power ow equations must have a solution, (ii) The contingency state variables xc must remain within limits, and (iii) Any inequality constraints, such as branch ow limits, must be satised. Remark 5.3. Typically, the limits on the contingency-dependent state variables xc and other functional inequality constraints are relaxed for the contingency cases compared to the base case. For example, system voltages are allowed to dip further during an emergency than under normal operating conditions. The relaxation of system limits is justied because operation under a contingency is temporary: when a contingency occurs, operators immediately begin re-conguring the system to return all branches and buses to normal operating limits. SCED is a restriction of the classic OPF formulation: for the same objective function, the optimal solution to SCED will be no better than the optimal solution without considering contingencies. The justication for the restriction is that SCED mitigates the risk of a system failure (blackout) should one of the contingencies occur. SCED has interesting connections to other areas of optimization. The motivation for SCED is theoretically similar to that of Robust Optimization (RO) [6], although
6 In SCED, the slack bus real and reactive power are treated as state variables because they must be allowed to change for each contingency in order for the system to remain feasible.
32
629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648
RO typically addresses continuous uncertain parameters rather than discrete scenarios. Additionally, because the constraints are separable for a xed u, SCED lends itself well to parallelization and decomposition algorithms [23]. 5.2.2. Security-Constrained Unit Commitment. In electric power systems operation, unit commitment (UC) refers to the scheduling of generating units such that total operating cost is minimized. UC diers from ED in that it operates across multiple time periods and schedules the on-o status of each generator in addition to its power output. UC must address generator startup and shutdown time and costs, limits on generator cycling, ramp rate limits, reserve margin requirements, and other scheduling constraints. UC is a large-scale, multi-period, mixed-integer nonlinear programming (MINLP) problem. Many UC formulations relax certain aspects of the problem in order to obtain a mixed-integer linear program (MILP) insteadfor instance by using linearized cost functions. If the power ow equations are added to the UC problem, the formulation becomes security-constrained unit commitment (SCUC). In SCUC, a power ow is applied at each time period to ensure that the scheduled generation satises not only the scheduling constraints but also system voltage and branch ow limits. In other words, SCUC ensures that the UC algorithm produces a generation schedule that can be physically realized in the power system. Because of its complexity, research on SCUC has accelerated only with the advent of faster computing capabilities. In SCUC, we introduce a time index t T and a set of binary control variables wit to the OPF formulation. Each wit indicates whether or not generator i is committed for time period t. The modied formulation becomes min
t T i G G SU wit Ci Pit + Ci wit (1 wi,t1 ) SD + Ci (1 wit ) wi,t1 ,
(5.13) i N, t T, i G, t T, i G, t T, i N, t T, i N, t T, (5.14) (5.15) (5.16) (5.17) (5.18) (5.19) (5.20) (5.21) (5.22) (5.23) (5.24)
s.t.
G L Pit (Vt , t ) = Pit Pit
Qit (Vt , t ) =
QG it
min max i it i max min ik ikt ik min max Tik Tikt Tik max Iikt (Vt , t ) Iik
G wit PiG,min Pit wit PiG,max ,min G,max wit QG QG it wit Qi i Vimin Vit Vimax
QL it
i N, t T,
ik K, t T, PiUp PReserve ik L, t T, i G, t T, t T.
ik H, t T,
PiDown
i G
649 650 651 652 653
wit PiG,max
G Pit
G Pi,t 1 i G
G Pit
The objective function (5.13) includes terms for unit startup costs C SU and shutdown costs C SD in addition to the generation costs in each time period. The generation limits (5.16)(5.17) are modied such that uncommitted units must have zero real and reactive power generation. Current limit constraint (5.22) is a compact expression of (5.11) with an added time index. Constraint (5.23) species positive and negative

654 655 656 657 658 659 660 661 662 663
33
generator ramp limits P Up and P Down , respectively; these are physical limitations of the generators. Constraint (5.24) requires a spinning reserve margin of at least PReserve ; sometimes this constraint is written such that PReserve is a fraction of the total load in each time period. The SCUC formulation (5.13)(5.24) is one of many possible formulations. Some formulations include more precise ramp limits and startup and shutdown characteristics; others include constraints governing generator minimum uptime and downtime. Because of the scale and presence of binary decision variables, SCUC is one of the most dicult power systems optimization problems. Zhu [37, ch. 7] and Bai and Wei [5] provide more discussion of SCUC, including detailed formulations. 5.2.3. Optimal Reactive Power Flow. Optimal reactive power ow (ORPF), also known as reactive power dispatch or VAR control, seeks to optimize the system reactive power generation in order to minimize the total system losses. In ORPF, the system real power generation is determined a priori, from the outcome of (for example) a DC OPF algorithm, UC, or other form of ED. A basic ORPF formulation is min s.t. P1 , Pi (V, ) = PiG PiL L Qi (V, ) = QG i Qi ,min G,max QG QG i Qi i Vimin Vi Vimax min max i i i max min ik ik ik min max Tik Tik Tik max Iik (V, ) Iik i N, i G, i N, i N, i N, (5.25) (5.26) (5.27) (5.28) (5.29) (5.30) (5.31) (5.32) (5.33)
ik K, ik L.
ik H,
The vector of control variables is u = P1 , QG i:iG , ik:ikH , Tik:ikK ,

664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679
while the vector of state variables x = (, V ) is identical to the classical formulation. In ORPF, all real power load and generation is xed except for the real power at the slack bus, P1 . Minimizing P1 is therefore equivalent to minimizing total system loss. One motivation for using ORPF is the reduction of the variable space compared to fully coupled OPF [9]; another is the ability to reschedule reactive power to optimally respond to changes in the system load without changing the system real power setpoints. Many interior point algorithms for OPF have focused specically on ORPF [11]. Zhu [37, ch. 10] discusses several approximate ORPF formulations and their solution methods. 5.2.4. Reactive Power Planning. Reactive power planning (RPP) extends the ORPF problem to the optimal allocation of new reactive power sourcessuch as capacitor bankswithin a power system in order to minimize either system losses or total costs. RPP modies ORPF to include a set of possible new reactive power sources; the presence or absence of each new source is modeled with a binary variable. The combinatorial nature of installing new reactive power sources has inspired many papers which apply heuristic methods to RPP [12].
34
A basic RPP formulation which minimizes total costs is min s.t. C 1 ( P1 ) +

i Q Install wi Ci ,
(5.34) i N, (5.35) (5.36) (5.37) (5.38) (5.39) (5.40) (5.41) (5.42) (5.43)
Pi (V, ) = PiG PiL QG i
Qi (V, ) = + QNew QL i i G,min G,max Qi QG Q i i New,min ,max New Qi wi Qi wi QNew i Vimin Vi Vimax min max i i i max min ik ik ik min max Tik Tik Tik max Iik (V, ) Iik
i G, i Q,
i N,
i N, ik H, ik L. i N,
ik K,
Install where Ci represents the capital cost of each new reactive power source i Q; New Qi is the amount of reactive power provided by each new reactive power source, New,min New,max subject to limits Qi and Qi ; and wi is a binary variable governing the decision to install each new reactive power source. The modied vector of control variables is New u = P1 , QG i:iG , wi:iQ , Qi:iQ , ik:ikH , Tik:ikK .
680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702
(5.44)
Some variants of RPP also include real power dispatch in the decision variables or include multiple load scenarios. By necessity, RPP optimizes with respect to uncertain future conditions typically reactive power requirements for worst-case scenarios. This uncertainty, together with the problem complexity, make RPP a very challenging optimization problem [34]. Zhang et al. [34, 35] review both formulations and solution techniques for RPP. 6. Data Exchange. Two common formats for the exchange of power ow and OPF case data are the IEEE Common Data Format [33] and the MATPOWER Case Format [38]. A number of publicly available test cases for OPF are distributed in one or both of these two formats [2, 38]. This section summarizes the structure of these formats and their relationship to the classical OPF formulation; the goal is to assist the reader in interpreting and applying available published data. 6.1. The IEEE Common Data Format. The IEEE Common Data Format (CDF) was rst developed in order to standardize the exchange of PF case data among utilities [33]. It has since been used to archive and exchange power systems test case data for the purpose of testing conventional PF and OPF algorithms. The format includes sections, or cards,7 for title data, bus data, branch data, loss zone data, and interchange data. Only the title, bus, and branch data are relevant for classical OPF as described in this primer. The full specication for the IEEE CDF can be found in [33] and an abbreviated description is available at [2]. Each IEEE CDF data card consists of plain text with elds delimited by character column. The title data card is a single line which includes summary information for
7 Originally,
the CDF data was exchanged among utilities by mail on paper card media.
35
Table 6.1 Field specication for IEEE Commmon Data Format bus data. The sixth column maps the eld to an index, parameter, or variable used in the classical OPF formulation given in 5.1. (Some elds are used indirectly via inclusion in Y .) Field 1 2 3 4 5 6 7 8 9 10 11 12 13 14 Columns 14 617 1920 2123 26 2833 3440 4149 5058 5967 6875 7783 8590 9198 Field Name Bus Number Bus Name Bus Area Loss Zone Numbera Bus Type Voltage Magnitude Voltage Angle Load Real Power Load Reactive Power Gen. Real Power Gen. Reactive Power Base Voltagea Desired Voltage Max. Reactive Power or Max. Voltage Magnituded Min. Reactive Power or Min. Voltage Magnituded Bus Shunt Conductance Bus Shunt Susceptance Remote Bus Number Data Type Integer Text Integer Integer Integer Numeric Numeric Numeric Numeric Numeric Numeric Numeric Numeric Numeric Numeric Numeric Numeric Numeric Numeric Integer Units Quantity in OPF i (bus index)
p.u. deg. MW MVAR MW MVAR kV p.u. MVAR p.u. MVAR p.u. p.u. p.u.
Specialb Vi i PiL QL i PiG QG i Vi (Specialc ) ,max QG i Vimax ,min QG i Vimin S gi bS i
15
99106
16 107114 17 115122 18 124127 a Optional eld b 0=PQ, 1=PQ (within voltage limits), 2=PV (within VAR limits), 3=Swing c Indicates target voltage magnitude for voltage-controlled (PV) buses d Gives reactive power limits if bus type is 2, voltage limits if bus type is 1
703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720
the case, including the power base SBase in MVA. The bus and branch data cards follow, beginning with the characters BUS DATA FOLLOWS and BRANCH DATA FOLLOWS, respectively, and ending with the ag characters -999. Each line within the card gives the data for a single bus or branch. Tables 6.1 and 6.2 list the IEEE CDF eld specications for bus and branch data, respectively. The elds include a mixture of SI and per-unit quantities. Conversion of all quantities to per-unit is required prior to use in an OPF formulation. Nominalvalued and unused elds in the data have zero entries. This quirk of the specication requires some caution in processing the data; for example, a value of 0.0 in the branch voltage ratio eld should be interpreted as a nominal tap ratio (T = 1.0). The IEEE CDF format is adapted for the compact exchange of system control data rather than OPF data. The eld structure therefore has several limitations: 1. Some IEEE CDF elds specify nal variable values (for instance, voltages Vi ) for conventional power ow. For OPF, these elds should be understood as a feasible or near feasible starting point rather than an optimal solution. (Due to rounding, the reported solution may not be strictly feasible.) 2. The elds bus type (bus eld 5) and branch type (branch eld 6) specify system control methods, and are therefore of limited use in OPF. However, the bus
36
Table 6.2 Field specication for IEEE Commmon Data Format branch data. The sixth column maps the eld to an index, parameter, or variable used in the classical OPF formulation given in 5.1. (Some elds are used indirectly via inclusion in Y .) Field 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Columns 14 69 1112 1315 17 19 2029 3039 4149 5155 5761 6367 6972 74 7782 8490 9197 Field Name Tap Bus Number Z Bus Number Line Areaa Loss Zone Numbera Circuit Number Branch Type Branch Resistance Branch Reactance Branch Shunt Susceptance Line Rating 1a Line Rating 2a Line Rating 3a Control Bus Number Side Voltage Ratio Phase Angle Min. Voltage Tap or Min. Phase Angled Max. Voltage Tap or Max. Phase Angled Tap Step Size or Phase Angle Step Sized Min. Voltage or Min. MVar Transfer or Min. MW Transfere Max. Voltage or Max. MVar Transfer or Max. MW Transfere Data Type Integer Integer Integer Integer Integer Integer Numeric Numeric Numeric Numeric Numeric Numeric Integer Integer Numeric Numeric Numeric Numeric Numeric Numeric Numeric Numeric Numeric Numeric Numeric Numeric Numeric Units Quantity in OPF i (from bus index) k (to bus index)
p.u. p.u. p.u. MVA MVA MVA
Specialb Rik Xik bSh ik max c Iik
p.u. deg. p.u. deg. p.u. deg. p.u. deg. p.u. MVar MW p.u. MVar
Tik ik min Tik min ik max Tik max ik
18
98104
19
105111
20
113119
21
120126
Numeric MW eld b 0=Transmission line, 1=Fixed T and , 2=Controllable T and xed (voltage control), 3=Controllable T and xed (MVAR control), 4=Fixed T and controllable c Conversion to per-unit current (using rated branch voltage) is required d Gives voltage tap limits or step if branch type is 2 or 3, phase angle limits or step if branch type is 4 e Gives voltage limits if branch type is 2, MVAR limits if branch type is 3, MW limits if branch type is 4
a Optional

721 722 723 724 725 726 727 728 729 730 731 732 733
37
and branch types govern the interpretation of certain other elds in the IEEE CDF, as described in the table footnotes. For example, for PQ buses, bus elds 14 and 15 give voltage limits Vimax and Vimin , respectively. For PV buses, these same elds ,max ,min instead give reactive power generation limits QG and QG , respectively. i i 3. For IEEE CDF elds which depend on the bus and branch types, the data are sucient for conventional PF but incomplete for OPF. For example, the IEEE CDF lacks voltage limits for PV buses and reactive power generation limits at PQ buses; the eld structure prevents these data from being available. The user must supply (or assume) values for the incomplete data. 4. The IEEE CDF lacks other data required for OPF, including generator real power limits and cost data. Given these limitations, publicly archived IEEE CDF case data is most useful for obtaining the network structure and associated bus and branch admittance data. 6.2. MATPOWER Case Format. MATPOWER [39] is an open-source software package for MATLAB8 including functions for both conventional PF and OPF. The MATPOWER case format is a set of standard matrix structures used to store power systems case data and closely resembles the IEEE CDF. The format is described in detail in the MATPOWER manual [38]. MATPOWER case data consists of a MATLAB structure with elds baseMVA, bus, branch, gen, and gencost. baseMVA is a scalar giving the system power base SBase in MVA. The remaining elds are matrices. Like the IEEE CDF, the MATPOWER case structure uses a mixture of SI and per-unit quantities and species nominal-valued branch tap ratios as 0 instead of 1.0. Tables 6.3, 6.4, and 6.5 describe the bus, branch, and gen matrices. The gencost matrix has the same number of rows as the gen matrix, but the column structure provides a exible description of the generator cost function. Column 1 species the type of cost model: 1 for piecewise linear or 2 for polynomial. Columns 2 and 3 give the generator startup and shutdown costs. The interpretation of column numbers 4 and greater depends on the type of cost model: For a piecewise linear cost model, column 4 species the number of coordinate pairs n of the form (P, C ) that generate the piecewise linear cost function. The next 2n columns, beginning with column 5, give the coordinate pairs (P0 , C0 ), . . . , (Pn1 , Cn1 ), in ascending order. The units of C are $/hr and the units of P are MW. For a polynomial cost model, column 4 species the number n of polynomial cost coecients. The next n columns, beginning with column 5, give the cost coecients Cn1 , . . . , C0 in descending order. The corresponding polynomial cost model is Cn1 P n1 + . . . + C1 P + C0 . The units are such that the cost evaluates to dollars $/hr for power given in MW. If gencost is included, then MATPOWER case data contains nearly all the information necessary to formulate the classical OPF problem as described in 5.1. However, MATPOWER makes no provision for including transformer tap ratios or phase shifting transformer angles in the set of decision variables; therefore, limits on these variables are not present in the data structure. The user must supply limits for these controls if they exist in the formulation. 7. Conclusion. In this primer, we have addressed the basic, practical aspects of Optimal Power Flow formulations. For the reader interested in learning more, particu8 MATLAB
734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765
766 767
is a popular technical computing environment produced by The MathWorks, Inc.
38
Table 6.3 Field specication for bus data matrix in MATPOWER case data (input elds only). The fth column maps the eld to an index, parameter, or variable used in the classical OPF formulation given in 5.1. (Some elds are used indirectly via inclusion in Y .) Column Field Description Data Type Units 1 Bus Number Integer 2 Bus Type Integer 3 Load Real Power Numeric MW 4 Load Reactive Power Numeric MVAR 5 Bus Area Integer 6 Bus Shunt Conductance Numeric MWb 7 Bus Shunt Susceptance Numeric MVARb 8 Voltage Magnitude Numeric p.u. 9 Voltage Angle Numeric deg. 10 Base Voltage Numeric kV 11 Loss Zone Integer 12 Max. Voltage Magnitude Numeric p.u. 13 Min. Voltage Magnitude Numeric p.u. a 1=PV, 2=PQ, 3=Swing, 4=Isolated b Specied as a MW or MVAR demand for V = 1.0 p.u. Quantity in OPF i (bus index) Speciala PiL QL i
S gi bS i Vi i
Vimax Vimin
Table 6.4 Field specication for branch data matrix in MATPOWER case data (input elds 111 only). The fth column maps the eld to an index, parameter, or variable used in the classical OPF formulation given in 5.1. (Some elds are used indirectly via inclusion in Y .) Column Field Description Data Type Units Quantity in OPF 1 Tap Bus Number Integer i (from bus index) 2 Z Bus Number Integer k (to bus index) 3 Branch Resistance Numeric p.u. Rik 4 Branch Reactance Numeric p.u. Xik 5 Branch Shunt Susceptance Numeric p.u. bSh ik max c 6 Line Rating (Long-term) Numeric MVA Iik 7 Line Rating (Short-term) Numeric MVA 8 Line Rating (Emergency) Numeric MVA 9 Voltage Ratio Numeric p.u. Tik 10 Phase Angle Numeric deg. ik 11 Branch Status Binary a Conversion to per-unit current (using rated branch voltage) is required
768 769 770 771 772 773 774 775 776 777 778
larly regarding optimization algorithms than have been used for OPF, we recommend any of the following: 1. Read the classical papers on OPF, for instance [4, 10, 27, 29]. These papers provide a detailed discussion of the foundations of OPF and provide context for more recent work. 2. Review textbooks which describe the OPF problem [32,37]. These textbooks provide clear, detailed formulations and also provide lists of relevant references. 3. Review the survey papers on OPF from the past several decades, for instance [11,12,1719]. Reading the older surveys prior to the more recent ones provides insight into how OPF has developed over time. 4. Experiment with the GAMS OPF formulations provided to accompany Ex-
39
Table 6.5 Field specication for generator data matrix in MATPOWER case data (input elds 110 only). The fth column maps the eld to an index, parameter, or variable used in the classical OPF formulation given in 5.1. Column Field Description Data Type Units Quantity in OPF 1 Bus Number Integer i (generator index) 2 Gen. Real Power Numeric MW PiG 3 Gen. Reactive Power Numeric MVAR QG i ,max 4 Max. Reactive Power Numeric MVAR QG i G,min 5 Min. Reactive Power Numeric MVAR Qi 6 Voltage Setpoint Numeric p.u. 7 Gen. MVA Basea Numeric MVA 8 Generator Statusb Binary 9 Max. Real Power Numeric MW PiG,max 10 Min. Real Power Numeric MW PiG,min a Defaults to system power base S Base b 0 indicates generator out of service (remove from OPF formulation)
779 780 781 782 783 784 785 786 787
ample 5.2, which are available in the GAMS model library [1]. Alternatively, install and experiment with the OPF capabilities available in MATPOWER [39]. Either software will provide insight into the practical challenges of OPF. The material presented in this primer should provide a sucient foundation for understanding the content of the references cited in this list. In recent years, OPF has become one of the most widely researched topics in electric power systems engineering. We hope that this primer encourages a similar level of engagement within the Operations Research community, particularly in the development of new, ecient OPF algorithms. Acknowledgment. We thank Kathryn Schumacher of the University of Michigan for her valuable comments and suggestions during the drafting of this primer.
REFERENCES [1] GAMS Model Library. Internet site. Available: http://gams.com/modlib/modlib.htm. [2] Power Systems Test Case Archive. Internet site. Available: http://www.ee.washington.edu/ research/pstca/. [3] O. Alsac, J. Bright, M. Praise, and B. Stott, Further developments in LP-based optimal power ow, IEEE Transactions on Power Systems, 5 (1990), pp. 697711. [4] O. Alsac and B. Stott, Optimal load ow with steady-state security, IEEE Transactions on Power Apparatus and Systems, PAS-93 (1974), pp. 745751. [5] X. Bai and H. Wei, Semi-denite programming-based method for security-constrained unit commitment with operational and optimal power ow constraints, IET Generation, Transmission & Distribution, 3 (2009), pp. 182197. [6] D. Bertsimas, D.B. Brown, and C. Caramanis, Theory and applications of robust optimization, SIAM Review, 53 (2011), pp. 464501. [7] R.C. Burchett, H.H. Happ, and K.A. Wirgau, Large scale optimal power ow, IEEE Transactions on Power Apparatus and Systems, 101 (1982), pp. 37223732. [8] J. Carpentier, Contribution to the economic dispatch problem, Bulletin de la Soci et e Fran caise des Electriciens, 8 (1962), pp. 431447. [9] G.C. Contaxis, C. Delkis, and G. Korres, Decoupled Optimal Load Flow Using Linear or Quadratic Programming, IEEE Transactions on Power Systems, PWRS-I (1986), pp. 17. [10] H.W. Dommel and W.F. Tinney, Optimal power ow solutions, IEEE Transactions on Power Apparatus and Systems, 87 (1968), pp. 18661876.
788 789
790
791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810
40
811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872
[11] S. Frank, I. Steponavice, and S. Rebennack, Optimal Power Flow: A Bibliographic Survey I, Formulations and Deterministic Methods, Energy Systems, tba (2012), p. tba. , Optimal Power Flow: A Bibliographic Survey II, Non-Deterministic and Hybrid Meth[12] ods, Energy Systems, tba (2012), p. tba. [13] Eberhard Freitag and Rolf Busam, Complex Analysis, Springer, 2 ed., 2009. [14] J.D. Glover, M.S. Sarma, and T.J. Overbye, Power System Analysis and Design, Cengage Learning, Stamford, CT, 4 ed., 2008. [15] S. Granville, Optimal reactive dispatch through interior point method, IEEE Transactions on Power Systems, 9 (1994), pp. 136146. [16] M. Greenberg, Advanced Engineering Mathematics, Prentice Hall, Upper Saddle River, NJ, 2 ed., 1998. [17] M. Huneault and F.D. Galiana, A survey of the optimal power ow literature, IEEE Transactions on Power Systems, 6 (1991), pp. 762770. [18] J.A. Momoh, M.E. El-Hawary, and R. Adapa, A Review of Selected Optimal Power Flow Literature to 1993 Part II: Newton, Linear Programming and Interior Point Methods, IEEE Transactions on Power Systems, 14 (1999), pp. 105111. [19] J. A. Momoh, M.E. El-Hawary, and R. Adapa, A Review of Selected Optimal Power Flow Literature to 1993 Part I: NonLinear and Quadratic Programming Approaches, IEEE Transactions on Power Systems, 14 (1999), pp. 105111. [20] J. Nocedal and S. Wright, Numerical Optimization, Springer, New York, NY, 2nd ed., 2006. [21] J. OMalley, Schaums Outline of Basic Circuit Analysis, Schaums Outline Series, McGrawHill, New York, NY, 2 ed., 2011. [22] J. Peschon, D.W. Bree, and L.P. Hajdu, Optimal Power-Flow Solutions for Power System Planning, Proceedings of the IEEE, 6 (1972), pp. 6470. [23] W. Qiu, A.J. Flueck, and F. Tu, A new parallel algorithm for security constrained optimal power ow with a nonlinear interior point method, in IEEE Power Engineering Society General Meeting, 2005, pp. 24222428. [24] R.L. Rardin, Optimization in Operations Research, Prentice Hall, Upper Saddle River, NJ, 1 ed., 1997. [25] N.S. Rau, Issues in the path toward an RTO and standard markets, IEEE Transactions on Power Systems, 18 (2003), pp. 435443. , Optimization Principles: Practical Applications to the Operation and Markets of the [26] Electric Power Industry, Wiley-IEEE Press, Hoboken, NJ, 2003. [27] B. Stott and E. Hobson, Power System Security Control Calculation Using Linear Programming Parts I and II, IEEE Transactions on Power Apparatus and Systems, PAS-97 (1978), pp. 17131731. [28] B. Stott, J. Jardim, and O. Alsac, DC power ow revisited, IEEE Transactions on Power Systems, 24 (2009), pp. 12901300. [29] D.I. Sun, B. Ashley, B. Brewer, A. Hughes, and W.F. Tinney, Optimal power ow by Newton approach, IEEE Transactions on Power Apparatus and Systems, 103 (1984), pp. 2864 2880. [30] G.L. Torres and V.H. Quintana, An Interior-Point Method For Nonlinear Optimal Power Flow Using Voltage Rectangular Coordinates, IEEE Transactions on Power Systems, 13 (1998), pp. 12111218. [31] L.S. Vargas, V.H. Quintana, and A. Vannelli, A tutorial description of an interior point method and its applications to security-constrained economic dispatch, IEEE Transactions on Power Systems, 8 (1993), pp. 13151323. [32] A.J. Wood and B.F. Wollenberg, Power Generation, Operation, and Control, John Wiley & Sons, Inc., New York, NY, 2 ed., 1996. [33] Working Group on a Common Format for Exchange of Solved Load Flow Data, Common Format For Exchange of Solved Load Flow Data, IEEE Transactions on Power Apparatus and Systems, PAS-92 (1973), pp. 19161925. [34] W. Zhang, F. Li, and L.M. Tolbert, Review of Reactive Power Planning: Objectives, Constraints, and Algorithms, IEEE Transactions on Power Systems, 22 (2007), pp. 21772186. [35] W. Zhang and L.M. Tolbert, Survey of Reactive Power Planning Methods, in IEEE Power Engineering Society General Meeting, vol. 2, June 2-16 2005, pp. 14301440. [36] X.-P. Zhang, Fundamentals of Electric Power Systems, John Wiley & Sons, Inc., 2010, ch. 1, pp. 152. [37] J. Zhu, Optimization Of Power System Operation, Wiley-IEEE, Piscataway, NJ, 2009. nchez, MATPOWER 4.1 Users Manual, Power Sys[38] R.D. Zimmerman and C.E. Murillo-Sa tems Engineering Research Center (PSERC), 2011. nchez, and R.J. Thomas, MATPOWER: Steady-State [39] R.D. Zimmerman, C.E. Murillo-Sa

873 874
41
Operations, Planning, and Analysis Tools for Power Systems Research and Education, IEEE Transactions on Power Systems, 26 (2011), pp. 1219.

Optimal Power Flow

Uploaded by

Copyright:

Available Formats

Optimal Power Flow

Uploaded by

Document Information

Original Description:

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Optimal Power Flow

Uploaded by

Copyright:

Available Formats

Division of Economics and Business Working Paper Series

A Primer on Optimal Power Flow: Theory, Formulation, and Practical Examples

c 2012 by the listed authors. All rights reserved.

STEPHEN FRANK AND STEFFEN REBENNACK

A PRIMER ON OPTIMAL POWER FLOW

Planning Horizon Years Months Weeks Days Hours Minutes Real-Time

asin Incre Decr

Unit Commitment SecurityConstrained Unit Commitment

Economic Dispatch Classic Optimal Power Flow

Automatic Generation Control

SecurityConstrained Unit Commitment

Optimal Reactive Power Flow

STEPHEN FRANK AND STEFFEN REBENNACK

A PRIMER ON OPTIMAL POWER FLOW

135 136 137 138 139 140 141 142

Yik Yik ik Gik

STEPHEN FRANK AND STEFFEN REBENNACK

172 173 174 175 176 177 178

179 180 181 182 183

Ii Ii Iik Iik PiL QL i L Si

201 202 203

204 205 206 207 208

A PRIMER ON OPTIMAL POWER FLOW

Branch (3,5) Bus 3 Bus 5

209 210 211 212 213

214 215 216 217 218 219 220 221

STEPHEN FRANK AND STEFFEN REBENNACK

coordinates, where, according to Eulers formula, a = c cos , b = c sin , a2 + b2 , b = arctan . a c=

A PRIMER ON OPTIMAL POWER FLOW

Total Current t Quadrature

i (t ) = i d(t ) + i q(t ) i d(t ) i q(t )

233 234 235 236 237 238 239 240

241 242 243

STEPHEN FRANK AND STEFFEN REBENNACK

240 cos(30) 240 cos (2 (377t) 30) dt W,

= 240 cos(30) 207.8 W,

A PRIMER ON OPTIMAL POWER FLOW

SBase 1 = . 2 VBase ZBase

281 282 283 284 285 286

STEPHEN FRANK AND STEFFEN REBENNACK

287 288 289 290 291 292 293

294 295 296 297 298 299

A PRIMER ON OPTIMAL POWER FLOW

STEPHEN FRANK AND STEFFEN REBENNACK

Rik Bus i bik 2

Fig. 4.2. branch model for cables and transmission lines.

316 317 318 319 320 321 322 323

324 325 326 327 328 329

A PRIMER ON OPTIMAL POWER FLOW

eik I' ei V' a:1

Fig. 4.3. O-nominal transformer branch model.

332 333 334 335

Similarly, Iki = yik Vk Vi , 1 = yik Vk Vi , a 1 = yik Vi + yik Vk . a

STEPHEN FRANK AND STEFFEN REBENNACK

343 344 345

A PRIMER ON OPTIMAL POWER FLOW

1 (gik cos ik bik sin ik ) Tik

1 (gki cos ki + bki sin ki ) , Tki 1 2 Tik 1 bik + bSh 2 ik +

i=k 1 bki + bSh , 2 ki