

small (250x250 max)
medium (500x500 max)
Large
Extra Large
large ( > 500x500)
Full Resolution


NATURAL GAS AND ELECTRICITY OPTIMAL POWER FLOW By SEUNGWON AN Bachelor of Engineering Korea Maritime University Pusan, Korea February, 1991 Master of Engineering Oklahoma State University Stillwater, Oklahoma May, 1999 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial ful¯llment of the requirements for the Degree of DOTOR OF PHILOSOPHY May, 2004 NATURAL GAS AND ELECTRICITY OPTIMAL POWER FLOW Thesis Approved: Thesis Advisor Dean of the Graduate College ii Acknowledgment I would like to express my sincere appreciation to my advisor Dr. Thomas W. Gedra for his intelligent supervision, friendship and con¯dence in me. I wish to thank the other members of my doctoral committee, Dr. Ramakumar, Dr. Yen and Dr. Misawa for their assistance throughout this work. I would like to give very special thanks to Dr. Ramakumar for his continuous guidance during my study at OSU. I wish to take this opportunity to express my gratitude to Mr. Lee Clark for his assistance and encouragement in all aspects of my work. I wish to express my sincere gratitude to the School of Electrical and Computer Engineering for providing me with this research opportunity and generous ¯nancial support. I am deeply indebted to all members of the Korean Catholic Community at Still water for their spiritual support. Especially, I would like to thank Mr. Hyunwoong Hong and Mrs. Eunjoo Kang for their brotherhood. I am deeply grateful to my wife, Misook Kim, for her patience and to our twin boys, Clemens and Martino for their existence. My ¯nal appreciation goes to my parents and brothers for their enduring love, support and encouragement. iii TABLE OF CONTENTS 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 AC LOADFLOW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1 Flows on Transmission Systems . . . . . . . . . . . . . . . . . . . . . 6 2.2 Load°ow Problem Statement . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 NewtonRaphson Method . . . . . . . . . . . . . . . . . . . . . . . . 13 3 ECONOMIC DISPACTH . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.1 Lossless Economic Dispatch . . . . . . . . . . . . . . . . . . . . . . . 18 3.2 Economic Dispatch with Transmission Losses . . . . . . . . . . . . . . 22 4 OPTIMAL POWER FLOW AND SOLUTION METHODS . . . . . . . . . 25 4.1 OPF by Newton's Method . . . . . . . . . . . . . . . . . . . . . . . . 27 4.2 PrimalDual InteriorPoint (PDIP) Method . . . . . . . . . . . . . . . 31 5 NATURAL GAS FLOW MODELING . . . . . . . . . . . . . . . . . . . . 39 5.1 Elements of Natural Gas Transmission Network . . . . . . . . . . . . 40 5.2 Network Topology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.3 Matrix Representations of Network . . . . . . . . . . . . . . . . . . . 43 5.4 Flow Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.5 Compressor Horsepower Equation . . . . . . . . . . . . . . . . . . . . 45 5.6 Conservation of Mass Flow . . . . . . . . . . . . . . . . . . . . . . . . 46 6 NATURAL GAS LOADFLOW . . . . . . . . . . . . . . . . . . . . . . . . 49 6.1 Load°ow Problem Statement . . . . . . . . . . . . . . . . . . . . . . . 50 6.2 Load°ow without Compressors . . . . . . . . . . . . . . . . . . . . . . 52 6.3 Load°ow with Compressors . . . . . . . . . . . . . . . . . . . . . . . 53 iv 7 UPFC MODELING FOR STEADYSTATE ANALYSIS . . . . . . . . . . 57 7.1 Operating Principles . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 7.2 Uncoupled Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 7.3 Ideal Transformer Model . . . . . . . . . . . . . . . . . . . . . . . . . 62 7.4 UPFC in a Transmission Line . . . . . . . . . . . . . . . . . . . . . . 65 8 GAS AND ELECTRICITY OPTIMAL POWER FLOW . . . . . . . . . . 69 8.1 Gas and Electric Combined Network . . . . . . . . . . . . . . . . . . 69 8.2 Gas and Electricity Optimal Power Flow . . . . . . . . . . . . . . . . 70 8.2.1 Cost, Bene¯t, and Social Welfare . . . . . . . . . . . . . . . . 71 8.2.2 Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 9 OPTIMAL LOCATION OF A UPFC IN A POWER SYSTEM . . . . . . . 79 9.1 Optimal Power Flow with UPFC . . . . . . . . . . . . . . . . . . . . 79 9.2 FirstOrder Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . 81 9.3 SecondOrder Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . 84 9.4 Estimation of Incremental Value . . . . . . . . . . . . . . . . . . . . . 87 10 RESULTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 10.1 Result of Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . 88 10.2 Result of GEOPF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 10.2.1 Test System I . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 10.2.2 Test System II . . . . . . . . . . . . . . . . . . . . . . . . . . 95 11 CONCLUSIONS AND FUTURE WORK . . . . . . . . . . . . . . . . . . . 101 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 A DERIVATIVES REQUIRED FOR GEOPF . . . . . . . . . . . . . . . . . 108 A.1 Electric Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 A.2 Gas Network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 A.3 Objective Function of GEOPF . . . . . . . . . . . . . . . . . . . . . . 119 B DERIVATIVES REQUIRED FOR UPFC SENSITIVITIES . . . . . . . . . 120 v C INPUT DATA FILE FORMAT . . . . . . . . . . . . . . . . . . . . . . . . 130 C.1 IEEE Common Data Format . . . . . . . . . . . . . . . . . . . . . . . 130 C.2 Generator Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 C.3 Electricity Load Data . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 C.4 Natural Gas Common Data Format . . . . . . . . . . . . . . . . . . . 135 C.5 Natural Gas Load Data . . . . . . . . . . . . . . . . . . . . . . . . . . 137 C.6 Natural Gas Source Data . . . . . . . . . . . . . . . . . . . . . . . . . 138 D MATLAB CODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 D.1 GEOPF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 D.2 Base Case OPF with Sensitivity Analysis . . . . . . . . . . . . . . . . 157 D.3 OPF with a UPFC . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 vi LIST OF FIGURES 1.1 Monthly natural gas price trends at AECO hub in Canada . . . . . . 2 1.2 Average Virginia mine coal price (19752000, Virginia) . . . . . . . . 3 2.1 Bus i and connection to transmission system . . . . . . . . . . . . . . 7 4.1 E®ect of barrier term . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2 Graphical representation of central path . . . . . . . . . . . . . . . . 37 5.1 Pipeline network representation . . . . . . . . . . . . . . . . . . . . . 40 5.2 Graph of the gas network . . . . . . . . . . . . . . . . . . . . . . . . . 41 7.1 General UPFC scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 58 7.2 Proposed UPFC model in a transmission line . . . . . . . . . . . . . . 59 7.3 Phasor diagram of UPFC inputoutput voltages and currents . . . . . 60 7.4 Uncoupled UPFC model in a transmission line. . . . . . . . . . . . . 61 7.5 UPFC ideal transformer model . . . . . . . . . . . . . . . . . . . . . 62 7.6 Simpli¯ed UPFC circuit . . . . . . . . . . . . . . . . . . . . . . . . . 64 7.7 Cascaded transmission line with a UPFC . . . . . . . . . . . . . . . . 66 8.1 Combined natural gas and electricity network . . . . . . . . . . . . . 70 10.1 Diagram of 5bus subset of IEEE 14bus system . . . . . . . . . . . . 88 10.2 Marginal and incremental values for 5bus system . . . . . . . . . . . 89 10.3 Marginal and incremental values for IEEE 14bus system . . . . . . . 90 10.4 Marginal and incremental values for IEEE 30bus system . . . . . . . 91 vii 10.5 Combined gas and electric network . . . . . . . . . . . . . . . . . . . 92 10.6 SW losses due to nonintegrated operation for test system I . . . . . . 96 10.7 Combined gas and electric network . . . . . . . . . . . . . . . . . . . 97 10.8 SW losses due to nonintegrated operation for test system II . . . . . 98 viii LIST OF TABLES 10.1 Real power generation, line loss and incremental value . . . . . . . . . 89 10.2 GEOPF results for test system I . . . . . . . . . . . . . . . . . . . . . 99 10.3 GEOPF results for test system II . . . . . . . . . . . . . . . . . . . . 100 ix CHAPTER 1 INTRODUCTION Restructuring of natural gas and electric industries, both known as network industries, and deregulation of their products have been undertaken in the U.S. to minimize the social ine±ciency incurred by regulated monopolies when natural monopoly has disappeared from their technologies [1]. Due to the successful deregulation of the natural gas market, the wellhead price of natural gas declined by 44 % between 1983 and 1997. Electric power generation technologies utilizing natural gas are generally less capitalintensive and often more technically e±cient than other alternatives. Thus, natural gas is playing an important role in electric power generation since 1990 par tially due to incentives of the Clean Air Act Amendments (CAAA) of 1990. Since 1997, more than 120 GW of new capacity has been added by merchant energy com panies at a lower cost per installed MW and with shorter construction times than in the past [2]. Most of this new generation is fueled by natural gas. However, due to the lack of ¯nancial instruments to lock in future prices to reduce price risk, the large increase in natural gas prices, and excess reserve margins, more than 110 GW of installed capacity belongs to entities that have belowinvestmentgrade credit ratings (i.e., junk bonds) [2]. Thus, we have been experiencing large price volatility in major natural gas markets while price volatility in electric markets has declined steadily over the past three years [2]. In addition, the °oor price for natural gas in summer 1 Figure 1.1: Monthly natural gas price trends at AECO hub in Canada 2003 reached over U.S. $4.00 per MMBTU, as shown in Figure 1.1 (on the basis of 1 U.S $ = 0.75 nominal Canadian $ in December 2003). In the meantime, due to regulatory uncertainties posed by the restructuring of electric industries and environmental concerns, baseload units such as coal¯red power plants have been rarely built since 1990. The amount of electric energy generated by coal¯red generators has been almost constant since that time, resulting in decreases in coal price as illustrated in Figure 1.2. With increasing gas prices and with a diverging price gap between coal and natural gas, generation using natural gas may be no longer competitive with coal¯red generation. Consequently, a signi¯cant amount of planned natural gas generation capacity addition using natural gas has been cancelled or postponed in 2003 [2]. The exposure to price volatility and price risk should be better understood in competitive energy markets to promote e±cient use and expansion of electric trans mission and gas transportation networks, and to facilitate e®ective competition in power generation and gas supply [3]. The largescale construction of gas¯red elec 2 Figure 1.2: Average Virginia mine coal price (19752000, Virginia) tric generation, combined with electric power restructuring, is expected to create an even greater convergence between electricity and natural gas markets. Lack of fuel diversity (i.e., excessive dependence on natural gas with volatile prices, as opposed to coal, especially newer \clean coal" technology) will likely increase volatility of gas and electric prices, even if the price of coal is lower and less volatile. Therefore, viewing these markets separately without recognizing their increasing convergence is myopic at best, and disastrous at worst. Economic forces in energy markets have been driving not only the optimal oper ation of energy networks but also e±cient pricebased system planning. The uni¯ed power °ow controller (UPFC) is one of the most technically promising devices in the Flexible AC Transmission System (FACTS) family [4, 5, 6, 7]. It has the capability to control both voltage magnitude and phase angle, and can also independently provide (positive or negative) reactive power injections. However, it cannot be installed in all possible transmission lines due to its high capital cost. Thus, a need exists for developing a costbene¯t analysis technique to determine if installation of a UPFC 3 would be bene¯cial and the best location to install the UPFC. In principle, determin ing the optimal location for a UPFC is simple. For each possible location, we place a UPFC in the power system model and calculate the cost savings with respect to a base case (with no new UPFC installed). The operating cost at each time throughout the year and for each potential location is determined using an optimal power °ow (OPF) program. However, the computational burden of evaluating this annual value for every possible line is immense because an OPF problem must be solved for each possible UPFC location and at each of several time periods. Therefore, an e±cient screening technique is desired to identify the most promising locations so that at each point in time throughout the year, the exhaustive calculations described above do not have to be carried out for every location that is a candidate for installing the UPFC. In this research, we introduce fundamental natural gas modeling for steadystate analysis, and propose general formulations to solve a combined natural gas and elec tric optimal power °ow (GEOPF). In addition, a screening technique is developed for greatly reducing the computation involved in determining the optimal UPFC lo cation in a large power system. This technique requires running only one optimal power °ow (OPF) to obtain UPFC sensitivities for all possible transmission lines. To implement the screening technique, we develop a new mathematical model of UPFC under steadystate, consisting of an ideal transformer with a complex turns ratio and a variable shunt admittance. We begin by reviewing basic concepts which describe power °ows in power systems. Economic dispatch is discussed next, followed by the introduction of the optimal power °ow concept and its solution methodology. A GEOPF model is described in detail after the introduction of natural gas network modeling. A screening technique based on sensitivity analysis is discussed next and a new UPFC model is described. The GEOPF model has been successfully applied to two test cases: a 15node gas network with (i) a 5bus electric network and (ii) a 9bus electric network. The 4 screening technique (electricity only) has also been implemented in a 5bus system and IEEE 14 and 30bus systems. Discussions of results are presented in Chapter 10. 5 CHAPTER 2 AC LOADFLOW Electric load°ow problems are solved routinely to study power systems under both normal operating conditions and under various contingencies using predicted data or to analyze \what if" scenarios. It is also of great importance in power system planning for future expansion. The principal information obtained from load°ow analysis is the magnitude and phase angle of the voltage at each bus. Using these values, we can calculate real and reactive power °ows in transmission lines and transformers, and line losses. In this chapter, we describe the formulation of the bus admittance (¹Y bus) matrix, which represents a transmission line network, and explain how complex power injec tions into a network are related to the bus voltage magnitudes and angles. Then, we construct the load°ow problem and discuss some of its characteristics. 2.1 Flows on Transmission Systems Consider a network of transmission lines connecting a set of buses in a power trans mission system. We will ¯nd out the relationship between complex power and current injections into the transmission system and the phasor voltages at each bus. Figure 2.1 [8] shows connections to a transmission system at bus i. The doublesubscript notation ¹ Sik and ~Iik, for k 6= i, indicates complex power or phasor current (respec 6 Figure 2.1: Bus i and connection to transmission system [8] tively) °owing from bus i towards bus k along line ik. For i = k, the double subscript notations ¹ Sii and ~Iii indicate complex power and phasor current °owing to ground through the shunt element at bus i. The single subscript notations ¹ Si and ~Ii indicate total complex power and phasor current injection into the transmission system at bus i. Generally, the complex power injected at bus i will consist of power generated by any generator(s) at bus i, minus any power consumed at bus i, or ¹ Si = ¹ SGi ¡ ¹ SDi : (2.1) If a line exists between bus i and bus k, and has a series impedance of ¹zik, the series admittance of the line will be de¯ned as ¹yik = 1 ¹zik : (2.2) If bus i and bus k are not directly connected, then ¹yik = 0. In a similar way, we de¯ne the admittance of any shunt element connecting bus i to ground as ¹yii, which includes transmission line capacitive susceptance jb=2. The current °owing from bus i towards bus k is ~Iik = 1 ¹zik (~V i ¡ ~V k) = ¹yik(~V i ¡ ~V k); (2.3) while the current °owing through the shunt element yii from bus i to ground is ~Iii = ¹yii~V i; (2.4) 7 where ~V i denotes the phasor voltage at bus i. If the total number of buses is n, the total current injected into the power transmission system at bus i is ~Ii = ¹yii ~Vi + Xn k=1 k6=i ¹yik ³ ~V i ¡ ~V k ´ ; (2.5) or, collecting the ~V i terms, ~Ii = ~V i Xn k=1 ¹yik ¡ Xn k=1 k6=i ¹yik ~V k: (2.6) The right side of the above equation motivates the following de¯nition: ¹ Yii = Xn k=1 ¹yik; (2.7) ¹ Yik = ¡¹yik; i 6= k: (2.8) With this de¯nition, we can write ~Ii = Xn k=1 ¹ Yik ~V k: (2.9) In matrix form, the de¯nitions in equations (2.7) and (2.8) can be written as ¹Y bus = 2 66666664 ¹y11 + ¢ ¢ ¢ ¹y1n ¡¹y12 ¢ ¢ ¢ ¡¹y1n ¡¹y21 ¹y21 + ¢ ¢ ¢ ¹y2n ¢ ¢ ¢ ¡¹y2n ... . . . ... ¡¹yn1 ¡¹yn2 ¢ ¢ ¢ ¹yn1 + ¢ ¢ ¢ ¹ynn 3 77777775 : In words, the diagonal element ii of the ¹Y bus matrix consists of the sum of all admit tances connected to bus i (whether they are series admittances connecting to another bus or are shunt admittances), while the o®diagonal element ik is the negative of the series admittance connecting bus i to bus k. In the case that there is no transmission line connecting bus i with bus k, ¹yik = 0 as mentioned above. Likewise, ¹yii = 0 if no shunt element is present between bus i and ground. Note that, since the series admittance ¹yik connecting bus i to bus k is the same as that connecting bus k to bus i, 8 we have ¹yik = ¹yki, so that the matrix ¹Y bus is symmetrical, assuming no phaseshifting elements, which we have implicitly done. If we de¯ne the bus voltage and bus current injection vectors as ~V = 2 66666664 ~V 1 ~V 2 ... ~V n 3 77777775 and ~I = 2 66666664 ~I1 ~I2 ... ~In 3 77777775 ; then we can write equation (2.9) in matrixvector form as ~I = ¹Y bus~V : (2.10) The matrix ¹Y bus is called the bus admittance matrix, since it relates the bus voltages and the bus current injections. The complex power injection at bus i can be written in terms of ¹Y bus matrix elements as ¹ Si = ~V i~I¤ i = ~V i Xn k=1 ¹ Y ¤ ik ~V¤ k ; (2.11) = Vi Xn k=1 VkYike(µi¡µk¡±ik); where ~V i = Vi\µi; ~V k = Vk\µk; and ¹ Yik = Gik + jBik = Yik\±ik: We can split equation (2.12) into real and reactive parts, in terms of the bus voltage magnitudes fVi; i = 1; : : : ; ng and bus voltage angles fµi; i = 1; : : : ; ng. Then the real and imaginary parts of equation (2.12) become Pi(~V ) = Xn k=1 ViVkYik cos(µi ¡ µk ¡ ±ik); (2.12) Qi(~V ) = Xn k=1 ViVkYik sin(µi ¡ µk ¡ ±ik): (2.13) 9 The real and reactive power injections at bus i need to satisfy the following conditions: PGi ¡ PLi = Pi(~V ); i = 1; : : : ; n; (2.14) QGi ¡ QLi = Qi(~V ); i = 1; : : : ; n; (2.15) where PGi = real power generation at bus i; PLi = real power load at bus i; QGi = reactive power generation at bus i; QLi = reactive power load at bus i: We can rewrite equations (2.12), (2.12) and (2.13) in matrix form as ¹S = 2 66666664 ¹ S1 ¹ S2 ... ¹ Sn 3 77777775 = 2 66666664 ~V 1~I¤ 1 ~V 2~I¤ 2 ... ~V n~I¤ n 3 77777775 = diagf~V g~I¤ = diagf~I¤g~V ; (2.16) P = real(¹S ); (2.17) Q = imag(¹S ); (2.18) where diagf~V g denotes the diagonal matrix with the elements of the vector ~V on its diagonal. 2.2 Load°ow Problem Statement Since the load°ow solution is a prerequisite for many other analytical studies, it is necessary to mathematically state the problem to prepare the ground for the estab lishment of its links to economic dispatch and optimal power °ow. The load°ow problem is stated below [8, 9, 10]: 10 ² Given a power system described by a ¹Y bus matrix, and given a subset of bus voltage magnitudes, bus voltage angles, and real and reactive power bus injec tions, ² Determine the other voltage magnitudes and angles and real and reactive power injections. More precisely, two of the four quantities Vi; µi; Pi; and Qi at each bus are speci¯ed, and the other two are to be determined. Each bus is classi¯ed based on the two known quantities, as follows: ² PV bus, or generator bus, or voltagecontrolled bus. For a bus i of this type, we assume that we know the real power injection Pi and the voltage magnitude Vi. This is because the real power generation can be controlled by adjusting the prime mover input, and the voltage magnitude can be controlled by adjusting the generator ¯eld current. However, certain buses without generators may have some means of voltage support at that bus (such as a synchronous condenser, capacitor banks, or a static VAr compensator). ² PQ bus, or load bus. For a nongenerator bus i, we assume that we know the real and reactive power injections Pi and Qi, and the bus voltage magnitude Vi and angle µi are to be determined. In practice, only real power load is known and the reactive power load is calculated based on an assumed power factor by Qi = QGi ¡ PLi p:f:i q 1 ¡ p:f:2 i : We assume that a load is characterized by its constant complex power demand (which does not depend on the bus voltage Vi, as opposed to assuming that the load is a constant current or constant impedance load in which case the complex power demand would depend on the bus voltage Vi.) 11 In fact, solving the load°ow problem with buses of only these two types is not in general possible. The ¯rst reason is that, in the load°ow equations, the bus voltage angles never appear by themselves, but instead appear only as angle di®erences of the form µi ¡ µk. Therefore, adding an angle to every bus voltage angle, will not change the values of real or reactive power injections in equations (2.12) and (2.13). Since phasor voltages are always expressed with respect to some reference voltage, we must decide on some reference phasor and refer all other phase angles to that reference. In fact, therefore, there are only n¡1 angles which in°uence the load°ows. We therefore pick one bus, say bus 1, to serve as our phasor reference, and set µ1 = 0o. Another reason is that solving a load°ow for a system containing only PV and PQ buses would imply that we know the real power injections at every single bus. In fact, we cannot specify all n real power injections, since we do not know all n real power injections until we solve the load°ow. Mathematically, the load°ow is overdetermined if we suppose we know all n real power injections. Specifying the injections at all buses is the same as specifying the real losses of the power system, which we cannot know until the load°ow is solved. Instead, we must pick one bus, and allow the real power injection at that bus to be whatever value is required to solve the load°ow equations. Thus, in addition to the PV and PQ bus types, we have a third bus type: ² Slack Bus, or V µ, or reference, or swing bus. Typically it is a generator bus, and the voltage angle of the slack bus serves as a reference for the angles of all other bus voltages. We assume that we know V1 and µ1, but we do not know P1 and Q1. The usual practice is to set µ1 = 0o and V1 = 1 although other values of V1 are ¯ne too. 12 2.3 NewtonRaphson Method Equations (2.14) and (2.15) are in the form of the vector nonlinear equation y = f (x), which can be solved by NewtonRaphson method [8, 9, 10]. The NewtonRaphson method is based on Taylor series expansion of f (x) about an operating point xo. y = f (xo) + @f @x ¯¯¯¯ x=xo (x ¡ xo) + H:O:T: (2.19) Neglecting the higher order terms in equation (2.19) and solving for x, we have x = xo + · @f @x ¯¯¯¯ x=xo ¸¡1 ¢ (y ¡ f (xo)) : (2.20) The NewtonRaphson method replaces xo by the old value xk and x by the new value xk+1 for the iterative solution as shown below xk+1 = xk + J¡1 ¢ ¢f ; (2.21) where J = @f @x ¯¯¯¯ x=xk ; ¢f = ³ y ¡ f (xk) ´ : Equation (2.21) is repeated until the mismatches ¢f are less than a speci¯ed tolerance, or the algorithm diverges. Now, let us construct a mathematical formulation to solve the AC load°ow prob lem using the NewtonRaphson method. Assume that there are g generators. We assume that bus 1 is the slack bus, and buses 2 through g are PV buses. The n ¡ g buses will be assumed to be PQ buses. Then the problem becomes: Given : µ1 P2 ¢ ¢ ¢ Pg Pg+1 ¢ ¢ ¢ Pn V1 V2 ¢ ¢ ¢ Vg Qg+1 ¢ ¢ ¢ Qn Determine : P1 µ2 ¢ ¢ ¢ µg µg+1 ¢ ¢ ¢ µn Q1 Q2 ¢ ¢ ¢ Qg Vg+1 ¢ ¢ ¢ Vn 13 Note that ¯nding P1 and Q1 through Qg is trivial, once all the voltage magnitudes and angles are known. The di±cult part is to ¯nd n¡1 unknown angles and the n¡g unknown voltage magnitudes. Let us de¯ne the x, y, and f vectors for the load°ow problem as x = 2 64 µ V 3 75 = 2 666666666666664 µ2 ... µn Vg+1 ... Vn 3 777777777777775 ; y = 2 666666666666664 P2 ... Pn Qg+1 ... Qn 3 777777777777775 ; f (x) = 2 666666666666664 P2(x) ... Pn(x) Qg+1(x) ... Qn(x) 3 777777777777775 : Note that there are 2n¡g¡1 nonlinear equations, and there are 2n¡1¡g unknown voltages and angles. So we have exactly the right number of variables to force all the mismatches to zero. We will guess the unknown voltage magnitudes and angles of xk, and compare the calculated values of P(xk) and Q(xk) to the known values of P and Q. Then, we can evaluate the mismatches by ¢f = 2 666666666666664 ¢P2 ... ¢Pn ¢Qg+1 ... ¢Qn 3 777777777777775 = 2 666666666666664 P2 ¡ P2(xk) ... Pn ¡ Pn(xk) Qg+1 ¡ Qg+1(xk) ... Qn ¡ Qn(xk) 3 777777777777775 : (2.22) In this set of equations, the functions Pi(xk) and Qi(xk) are those obtained from equations (2.12) and (2.13), or (2.17) and (2.18), while the numbers Pi and Qi are the known values of real and reactive power injections at the buses where the injections are known. The resulting quantities ¢Pi and ¢Qi are known as the real and reactive mismatch terms, because they represent the di®erence between the known injections and the values calculated based on our guesses. 14 To solve the AC load°ow problem, we need to update the Jacobian at each iter ation. Many authors [9, 10] derived the Jacobian directly from equations (2.12) and (2.13). Even though this method is straightforward, it requires at least two for loop routines in Matlab. To avoid the for loop command, and to improve the speed of each iteration, we present an e±cient technique to construct the Jacobian. Let us de¯ne JSV by JSV = 2 6666666664 @ ¹ S1 @V1 @ ¹ S1 @V2 ¢ ¢ ¢ @ ¹ S1 @Vn @ ¹ S2 @V1 @ ¹ S2 @V2 ¢ ¢ ¢ @ ¹ S2 @Vn ... ... . . . ... @ ¹ Sn @V1 @ ¹ Sn @V2 ¢ ¢ ¢ @ ¹ Sn @Vn 3 7777777775 : The diagonal and o®diagonal elements of JSV are given by @ ¹ Si @Vi = ejµi " Xn k=1 ¹ Yik ~V k #¤ + ~V i ¹ Y ¤ ii e¡jµi ; (2.23) @ ¹ Si @Vk = ~V i ¹ Y ¤ ike¡jµk ; i 6= k (2.24) By noticing similarities in the above two equations, we can rewrite JSV in compact form to be used in Matlab by JSV = diag(~V ) ¤ conj(¹Y bus) ¤ conj(diag(expjµ)) + diag(expjµ :  ¤{zconj(~I}) element to element multiplication ): (2.25) Now, let us de¯ne JSµ by Jsµ = 2 6666666664 @ ¹ S1 @µ1 @ ¹ S1 @µ2 ¢ ¢ ¢ @ ¹ S1 @µn @ ¹ S2 @µ1 @ ¹ S2 @µ2 ¢ ¢ ¢ @ ¹ S2 @µn ... ... . . . ... @ ¹ Sn @µ1 @ ¹ Sn @µ2 ¢ ¢ ¢ @ ¹ Sn @µn 3 7777777775 : The diagonal and o®diagonal elements of JSµ are given by @ ¹ Si @µi = j ¹ Si ¡ j~V i ¹ Y ¤ ii ~V ¤ i ; (2.26) @ ¹ Si @µk = ¡j~V i ¹ Y ¤ ik ~V ¤ k ; i 6= k: (2.27) 15 In a similar way, JSµ can be expressed in compact form by JSµ = ¡jdiag(~V ) ¤ conj(¹Y bus) ¤ conj(diag(~V )) + diag(j¹S ): (2.28) We can split the Jacobian into real and reactive parts as follows: PSµ = real(JSµ ); PSV = real(JSV ); QSµ = imag(JSµ ); QSV = imag(JSV ): Then, the full Jacobian J is constructed as J = 2 64 PSµ PSV QSµ QSV 3 75 : Since the voltage magnitudes at PV buses and the angle at the slack bus are known, their corresponding columns will be truncated from the Jacobian J. In addition, since the real power injection at the slack bus, and the reactive power injections at PV buses and at the slack bus are unknown, their corresponding rows are removed from the Jacobian J. Then, the modi¯ed Jacobian Jm at the kth iteration becomes Jm = 2 666666666666666664 @P2(xk) @µ2 ¢ ¢ ¢ @P2(xk) @µn @P2(xk) @Vg+1 ¢ ¢ ¢ @P2(xk) @Vn ... . . . ... ... . . . ... @Pn(xk) @µ2 ¢ ¢ ¢ @Pn(xk) @µn @Pn(xk) @Vg+1 ¢ ¢ ¢ @Pn(xk) @Vn @Qg+1(xk) @µ2 ¢ ¢ ¢ @Qg+1(xk) @µn @Qg+1(xk) @Vg+1 ¢ ¢ ¢ @Qg+1(xk) @Vn .... . . ... ... . . . ... @Qn(xk) @µ2 ¢ ¢ ¢ @Qn(xk) @µn @Qn(xk) @Vg+1 ¢ ¢ ¢ @Qn(xk) @Vn 3 777777777777777775 : By using the mismatches in equation (2.22) and the modi¯ed Jacobian Jm, we will update the unknown variables x by equation (2.21). Finally, once we have forced all the real and reactive mismatches to reasonably small values, by ¯nding the unknown values of voltage magnitude and angle, we can now go back and ¯nd the unknown real and reactive injections. 16 CHAPTER 3 ECONOMIC DISPACTH The economic dispatch problem is de¯ned as the process of providing the required real power load demand and line losses by allocating generation among a set of online generating units such that total generation cost is minimized [8, 11, 12]. Let C be the generating cost in $/hr, and it is often modelled analytically as a quadratic function of the power generated [8]. The cost function is derived from the generator heat rate curve. Analytically, the heat rate curve with a unit of BTU/kWh is represented in the following form: H(PG) = a PG + b + cPG; (3.1) where PG is the unit's real power generation level measured in kW. Given the heat rate curve, the next important function is the fuel rate, which is simply the rate, in BTU/hr, of consumption of fuel energy. So F(PG) = PGH(PG) = a + bPG + cP2G : (3.2) If the cost of fuel is known as k $/BTU, then the cost function with a unit of $/hr becomes C(PG) = kF(PG) = ka + kbPG + kcP2G = ® + ¯PG + °P2G : (3.3) In classic economic dispatch problems, the essential constraint on the operation is that the sum of the output powers must equal the load demand plus losses. In 17 addition, there are two inequality constraints that must be satis¯ed for each of the generator units. That is, the power output of each unit must be greater than or equal to the minimum power permitted and must also be less than or equal to the maximum power permitted on that particular unit. 3.1 Lossless Economic Dispatch We assume that the utility is responsible for supplying its customers' load. The utility's objective is to minimize the total cost of generation, assuming that trans mission losses are neglected. Then, the ideal economic dispatch problem is stated in terms of minimizing total generation cost subject to satisfying the total load demand. Mathematical formulation of lossless economic dispatch problem can be expressed as min fPG1 ;¢¢¢ ;PGng Xn i=1 Ci(PGi) (3.4) subject to: Xn i=1 PGi = PD (3.5) and Pmin Gi · PGi · Pmax Gi ; i = 1; : : : ; n (3.6) where PD = Xn i=1 PDi is the total system load. We will ignore the constraints in equation (3.6) on generator limits, assuming that these limits are not binding. Considering only constraint (3.5), the Lagrangian is L = Xn i=1 Ci(PGi) ¡ ¸ " Xn i=1 PGi ¡ PD # : (3.7) Di®erentiating with respect to PGi , we obtain the ¯rstorder conditions 0 = @L @PGi = MCi ¡ ¸; i = 1; : : : ; n (3.8) 18 or MCi = ¸; i = 1; : : : ; n where MCi = @Ci @PGi : The quantity ¸, the Lagrange multiplier, associated with the energy balance con straint in equation (3.5), is universally called the \system lambda" and is the price associated with generating slightly more energy. Thus, the criterion for optimal eco nomic distribution of load among n generators is that all generators should generate output at the same marginal operating cost, sometimes called incremental cost. That is @C1 @PG1 = @C2 @PG2 = ¢ ¢ ¢ = @Cn @PGn : (3.9) In the case of quadratic cost functions, the solution to the economic dispatch problem can be calculated analytically, as follows [8]. We have, for each unit, ¸ = MCi = ¯i + 2°iPGi ; (3.10) which can be solved for PGi to give PGi = ¸ ¡ ¯i 2°i : (3.11) The total generation at this ¸ is obtained by summing over all the units. Setting this total equal to the load PD, we obtain PD = Xn i=1 PGi = Xn i=1 ¸ ¡ ¯i 2°i = ¸ Xn i=1 1 2°i ¡ Xn i=1 ¯i 2°i ; (3.12) which can be solved to obtain the ¸ required for a given PD: ¸ = · PD + Pn i=1 ¯i 2°i ¸ Pn i=1 1 2°i (3.13) 19 Finally, this value of ¸ can be substituted back into the expression for PGk to obtain PGk = 1 2°k Pn i=1 1 2°i " PD + Xn i=1 ¯i 2°i # ¡ ¯k 2°k ; (3.14) which gives the dispatch for unit k directly in terms of the total load PD. We de¯ne the \participation factor" for unit k as Kk = 1 2°k Pn i=1 1 2°i : We can see from the de¯nition that Xn i=1 Ki = 1: Generally, the meaning of the participation factor is Kk = dPGk dPD : Thus, if the load were to increase by a small increment (say 1MW), unit k would supply the fraction Kk of the increase. The fact that the participation factors sum to 1 simply means that any increase in load is met exactly by an increase in generation. We can include the generator limit constraints (3.6) in one of two ways. Tradi tionally, we ignore the limits, perform economic dispatch, and check to see if any limits are violated. If the limit for one generator is violated, we set the generator to its limit, remove it from the set of generators included in economic dispatch, and subtract the limit from the total load PD. In other words, we take the generator \o® of economic dispatch" and then proceed to treat it as a negative load1. More formally, we can just include the generation limits in the Lagrangian func tion. Let ¹min i be the Lagrange multiplier for the lower limit of unit i, and ¹max i be 1This may result in \cycling" of the set of active constraints, especially if there are many diverse generators. Other optimization techniques, such as a primaldual interiorpoint (PDIP) method [13] can be used to avoid this problem. 20 the multiplier for the generator's upper limit. Then, the Lagrangian becomes L = Xn i=1 Ci(PGi) ¡ ¸ " Xn i=1 PGi ¡ PD # + Xn i=1 ¹min i £ Pmin Gi ¡ PGi ¤ + Xn i=1 ¹max i £ PGi ¡ Pmax Gi ¤ : (3.15) Now, taking the derivative of L with respect to each PGi , and setting to zero, we obtain 0 = @L @PGi = MCi ¡ ¸ ¡ ¹min i + ¹max i ; i = 1; : : : ; n (3.16) or MCi = ¸ + ¹min i ¡ ¹max i ; i = 1; : : : ; n (3.17) We also have the complementary slackness conditions on the inequality constraints: ¹min i £ Pmin Gi ¡ PGi ¤ = 0; and ¹max i £ PGi ¡ Pmax Gi ¤ = 0: But PGi cannot be at its upper and lower limits simultaneously, so there are only three possibilities for each generator: PGi = Pmax Gi ) ¹max i ¸ 0; ¹min i = 0; Pmin Gi · PGi · Pmax Gi ) ¹max i = 0; ¹min i = 0; PGi = Pmin Gi ) ¹max i = 0; ¹min i ¸ 0: If the limits on some generators are binding, the generators may not be able to generate output at the system incremental cost ¸. Suppose that none of generators reach their limits, and all generators participate in economic dispatch. As the load increases by a small increment, each generator will supply the fraction of the load increase, which is determined by the participation factor. If a generator hits its maximum limit, the generator cannot participate in economic dispatch even though it has a lower incremental cost than the system incremental cost ¸. 21 On the other hand, suppose that the load decreases and one of generators hits its minimum limit. Since the generator must generate the required minimum capacity in order to stay online, the generator cannot participate in economic dispatch for the load decrements. Thus, if the load decreases further, the incremental cost of the generator hit its minimum limit is greater than or equal to the system incremental cost ¸. The optimal criterion for lossless economic dispatch with generator limits can be summarized as [11] PGi = Pmax Gi ) dCi dPGi · ¸; Pmin Gi · PGi · Pmax Gi ) dCi dPGi = ¸; PGi = Pmin Gi ) dCi dPGi ¸ ¸: 3.2 Economic Dispatch with Transmission Losses We now wish to include the e®ect of transmission losses on economic dispatch [8, 12]. In lossless dispatch, the location of individual loads did not matter, but when transmission losses are considered, the solution depends both on the load locations and on the outputs of individual generators around the system. If we know the load PDi , and generation PGi at each bus in the system, we can calculate the real power injections P2; : : : ; Pn using Pi = PGi ¡ PDi . From load°ow, we can then calculate the power injection at the slack bus, P1 = PG1 ¡ PD1 . From this information, we can calculate losses as PL(PG2 ; ¢ ¢ ¢ ; PGn; PD2 ; ¢ ¢ ¢ ; PDn) = Xn i=2 (PGi ¡ PDi) + P1(PG2 ; ¢ ¢ ¢ ; PGn; PD2 ; ¢ ¢ ¢ ; PDn): (3.18) Note that in our formulation, we do not consider PL to be a function of PG1 , since the slack bus generation PG1 is completely determined by the injections at other buses and the slack bus demand PD1 . 22 Now our problem becomes min fPG1 ;¢¢¢ ;PGng Xn i=1 Ci(PGi) (3.19) subject to: Xn i=1 PGi = PD + PL (PG2 ; ¢ ¢ ¢ ; PGn) (3.20) and Pmin Gi · PGi · Pmax Gi ; i = 1; : : : ; n (3.21) where we have suppressed the dependence of PL on the loads fPDig, which are as sumed to be ¯xed. The Lagrangian for the problem now becomes L = Xn i=1 Ci(PGi) ¡ ¸ " Xn i=1 PGi ¡ PD ¡ PL(PG2 ; ¢ ¢ ¢ ; PGn) # : (3.22) Di®erentiating and setting to zero yields 0 = @L @PG1 = MC1 ¡ ¸; and 0 = @L @PGi = MCi ¡ ¸ µ 1 ¡ @PL @PGi ¶ ; i = 2; : : : ; n or ¸ = MC1; and ¸ = MCi 1 1 ¡ @PL @PGi ; i = 2; : : : ; n Finally, let us de¯ne the losspenalty factors L1 = 1; and Li = 1 1 ¡ @PL @PGi ; i = 2; : : : ; n (3.23) Then we can write our condition for economic dispatch as ¸ = MCi ¢ Li; i = 1; : : : ; n (3.24) The idea behind the loss penalty factors can be explained as follows. If increasing a generator's output increases system losses, then that generator'entire increment of 23 generation is not available to the system. So, if a generator's output increases by 1MW but losses increase by 0.1MW, the generator has only provided a net increase in system generation of 0.9MW. Thus, the cost of this increment is not the generator's marginal cost but MC=0:9 instead, or MCi ¢Li. Similarly, an increase in a generator's output could result in a decrease in system losses. Such a generator would have a loss penalty factor Li < 1, since the e®ective cost of incremental generation from such a generator would be less than its actual marginal cost. The loss penalty factors can either be obtained as the result of running a load°ow, or they can be obtained from a quadratic approximation to the loss function, called Bmatrix formula [11]. 24 CHAPTER 4 OPTIMAL POWER FLOW AND SOLUTION METHODS Note that the classic economic dispatch (ED) problem does not strictly take into account power °ows in the transmission system. Typically, one would solve the eco nomic dispatch problem, then the load°ow problem, and repeat the process. The most recent load°ow provides the losspenalty factors for the current operating state, then the ED changes the state slightly to minimize cost, and so on [8]. Provided no transmission lines are overloaded, and no voltages are outside of the allowable range, this works well. However, if the result of the economic dispatch is fed into the load °ow, and the result indicates that transmission line loading is excessive, or that bus voltage magnitudes are outside of the allowable range (typically 95% to 105% of the nominal value), no information is given by the load°ow to indicate how to redispatch generation to alleviate the overloads or restore acceptable voltage levels. Often, the system dispatcher, being extremely familiar with the system, will take some units o® of economic dispatch and manually redispatch them to remove the line overload. There is no guarantee that such a procedure will result in the minimum cost subject to operating limits although sometimes they do fairly well. Similarly, economic dispatch does not consider reactive °ows or bus voltages, whereas load°ow considers these quantities as given or to be found from other voltages and reactive injections [8]. Neither problem considers the adjustment of bus voltages to help to minimize cost, even though reactive power °ows can contribute signi¯cantly 25 to real power losses and thus to overall cost. However, we can reformulate our optimization problem by including economic dispatch, voltage and reactive injections as decision variables, and various operational limits. The optimal power °ow (OPF) combines economic dispatch and load°ow into a single problem, which can be written in very general terms as min Y C(Y ) (4.1) subject to: hi(Y ) = 0; i = 1; ¢ ¢ ¢ ; n; gi(Y ) · 0; i = 1; ¢ ¢ ¢ ; m; where ² Y : a vector of control and state variables, ² C(Y ): an objective function, ² h(Y ): a set of equality constraints, which are power °ow equations, ² g(Y ): a set of inequality constraints, such as voltage limits, generator capacity limits, and line °ow limits. Tapchanging transformers and/or Flexible AC Transmission System (FACTS) de vices also can be incorporated in the OPF problem [9, 14]. The OPF has been used as a tool to improve power system planning, and oper ation by adjusting the objective function and/or the constraints. From the power system planning point of view, the OPF can be used to determine the optimal types, sizes, settings, capital costs, and optimal locations of resources, such as generators, transmission lines, and FACTS devices [14, 15]. Another application of the OPF is to determine various system marginal costs. It can be used to determine shortrun electric pricing (i.e. spot pricing), transmission 26 line pricing, and pricing ancillary services such as voltage support through MVAR support. Since the work in this research is based on solving the proposed optimization problem by a primaldual interiorpoint (PDIP) method using logarithmic barrier function, the PDIP will be discussed in detail after a brief introduction of Newton's method. 4.1 OPF by Newton's Method NewtonRaphson method has been the standard solution algorithm for the economic dispatch and load°ow problems for several decades [11, 12]. Newton's method is a very powerful algorithm because of its rapid convergence near the solution. This property is especially bene¯cial for power system applications because an initial guess close to the solution is easily obtained [16]. For example, voltage magnitude at each bus is presumably near the rated system value, generator outputs can be estimated from historical data, and transformer tap ratios are near 1.0 p.u. during steadystate operation. The solution of the constrained optimization problem stated in (4.1) requires the mathematical formation of the Lagrangian by L(Y; ¸; ¹) = C(Y ) + Xn i=1 ¸ihi(Y ) + X i2A ¹igi(Y ); (4.2) where ¸i is the Lagrange multiplier for the ith equality constraint. Assuming that we know which inequality constraints are binding, and have put them in the set A, then the inequality constraints can now be enforced as equality constraints. Thus the ¹0 is in equation (4.2) have the same property as ¸0 is and they are the Lagrange multipliers for binding inequality constraints. However, we need ¹i ¸ 0 for every i [8]. We can ignore the inequality constraints that are not binding since their ¹0s are 27 known to be zero by complementary slackness condition [8]. That is gi(Y ) · 0 ) ¹i = 0; gi(Y ) = 0 ) ¹i ¸ 0: Therefore, only binding inequality constraints are included in the Lagrangian function (4.2) with corresponding nonzero ¹0s. Solution of a constrained optimization problem can be solved by adjusting control and state variables, and Lagrange multipliers to satisfy the following firstorder necessary optimality conditions : 1) @L @Yi = 0; (4.3) 2) @L @¸i = 0; (4.4) 3) @L @¹i = gi = 0; (4.5) 4) ¹i ¸ 0 and ¹igi = 0: (4.6) The above equations are also called the KarushKuhnTucker (KKT) conditions. Let us de¯ne !(z) = rzL(z) = · @L @Y @L @¸ @L @¹A ¸T = 0; (4.7) where z is a vector of [ Y T ¸T ¹T A ]T , and A represents the binding inequality constraints. To solve the KKT conditions, Newton's method is applied by using the Taylor's series expansion around a current point zp as: !(z) = !(zP ) + @!(z) @z ¯¯¯¯ z=zp ¢ (z ¡ zp) + 1 2 (z ¡ zp)T ¢ @2!(z) @z2 ¯¯¯¯ z=zp ¢ (z ¡ zp) + ¢ ¢ ¢  {z } H.O.T = 0; The current point zp can either be an initial guess in the ¯rst iteration of the computation, or the estimate solution from the prior iteration. Recall that we want !(z) = 0. By ignoring the high order terms (H.O.T) and de¯ning ¢z = z ¡ zp, the 28 above equation can be rewritten as: @!(z) @z ¯¯¯¯ z=zp ¢ ¢z = ¡!(zp): (4.8) The quantity ¢z is the update vector, or the Newton step, and it tells how far and in which direction the variables and multipliers should move from this current point zp to get closer to the solution. Since !(z) is the gradient of the Lagrangian function L(z), equation (4.8) can be written in terms of the Lagrangian function L(z) as: @2L(z) @z2 ¯¯¯¯ z=zp ¢ ¢z = ¡ @L(z) @z ¯¯¯¯ z=zp : (4.9) Or, simply W ¢ ¢z = ¡!(z); (4.10) where W denotes the second order derivatives (or the Hessian matrix) and ! is the gradient, both of the Lagrangian function with respect to z evaluated at the current point. Equation (4.10) can be written in matrix form as 2 6666664 HY JT AT J 0 0 A 0 0 3 7777775 2 6666664 ¢Y ¢¸ ¢¹A 3 7777775 = ¡ 2 666664 rY L r¸L r¹AL 3 777775 ; (4.11) where the Hessian matrix HY and Jacobian matrices J and A are given as follows: HY = @2L(z) @Y 2 ; J = @2L(z) @¸@Y = @h(Y ) @Y ; A = @2L(z) @¹A@Y = @gA(Y ) @Y ; where gA is a set of binding inequality constraints. The Newton step can be obtained by solving (4.10). Then a vector of estimated solution for the next iteration is updated as: zp+1 = zp + ®¢z; (4.12) 29 where ® is usually 1, but can be adjusted to values above or below 1 to speed up convergence or cause convergence in a divergent case. It is important that special attention be paid to the inequality constraints. Equation (4.2) only includes binding inequality constraints being enforced as equality constraints. Thus, after obtaining an updated set of variables and multipliers, a new set of binding inequality constraints (or what we think is the active set) should be determined as follows [12]: ² If the updated ¹0s of the constraint functions in the current active set are zero or have become negative, then the corresponding constraints must be released from the current active set because ¹i < 0 implies that gi = 0 keeps the trial solution at the edge of the feasible region instead of allowing the trial solution to move into the interior of the feasible region. ² If other constraint functions evaluated at the updated variables violate their limits, then those constraints must be included in the new active set. The variable ® may be chosen to prevent constraint violations, but ® < 1 to avoid infeasibility implies that the constraint would otherwise be violated. As a result, if ¹i is positive, continued enforcement will result in an improvement of the objective function, and enforcement is maintained. If ¹i is negative, then enforcement will result in an decrease of the objective function, and enforcement is stopped. Once the active set has been updated, !(zp+1) is checked for convergence. There are several criteria for checking convergence of Newton's method. The convergent tolerance may be set on the maximum absolute value of elements in !(z), or on its norm. If the updated zp+1 does not satisfy the desired convergence criterion, the Newton step calculation is repeated. 30 4.2 PrimalDual InteriorPoint (PDIP) Method One of disadvantages of Newton's method is to identify a set of binding inequality constraints, or active constraints. Among several methods to avoid the di±culty asso ciated with guessing the correct active set, the PDIP method has been acknowledged as one of the most successful [12, 17, 18]. Interior point methods for optimization have been widely known since the publication of Karmarkar's seminal paper in 1984 [19]. Barrier function methods were proposed much earlier in Russia but little attention was paid because the algorithm was so slow in implementation. Later, this method was shown to be equivalent to the interior point methods. Karmarkar's method re sults in numerical illconditioning although this problem is not so bad with the PDIP method. The method uses a barrier function that is continuous in the interior of the feasible set, and becomes unbounded as the boundary of the set is approached from its interior. Two examples of such a function [13, 20] are the logarithmic function, as shown in Figure 4.1 Á(Y ) = ¡ Xm i=1 ln(¡gi(Y )); (4.13) and the inverse of the inequality function Á0(Y ) = Xm i=1 1 ¡gi(Y ) ; (4.14) where gi(Y ) · 0: The barrier method generates a sequence of strictly feasible iterates that converge to a solution of the problem from the interior of the feasible region [13, 20]. To apply the primaldual interior point algorithm to the OPF problem that has equality and inequality constraints, we construct the nonlinear equality and 31 inequalityconstrained optimization problem as min Y f(Y ) (4.15) subject to hi(Y ) = 0; i = 1; ¢ ¢ ¢ ; n; gi(Y ) · 0; i = 1; ¢ ¢ ¢ ; m: To solve the problem, we ¯rst form the logarithmic barrier function as B = f(Y ) ¡ º Xm i=1 ln (¡gi(Y )) ; (4.16) where the parameter º is referred to as the barrier parameter, a positive number that is reduced to approach to zero as the algorithm converges to the optimum. Then we solve a sequence of constrained minimization problems of the form min (Y;¸;º) B(Y; ¸; º) (4.17) subject to hi(Y ) = 0; i = 1; ¢ ¢ ¢ ; n; for a sequence fºkg of positive barrier parameters that decrease monotonically to zero. The solution of this problem by Newton's method requires the formulation of the Lagrangian function L = f(Y ) ¡ Xn i=1 ¸ihi(Y ) ¡ º Xm i=1 ln (¡gi(Y )) : (4.18) Because the barrier term is in¯nite on the boundary of the feasible region, as shown in Fig. 4.1, it acts as a repelling force that drives the current trial solution away from the boundary into the interior of the feasible region. As the barrier parameter º is decreased, the e®ect of the barrier term is diminished, so that the iterates can grad ually approach the constraint boundaries of the feasible region for those constraints which eventually turn out to be binding. 32 Figure 4.1: E®ect of barrier term To solve equation (4.17), we need to ¯nd the gradient of L with respect to Y : rY L = rY f(Y ) ¡ Xn i=1 ¸irY hi(Y ) + º Xm i=1 1 ¡gi(Y ) rY gi(Y ) = rY f(Y ) ¡ Xn i=1 ¸irY hi(Y ) + º Xm i=1 1 si rY gi(Y ) = rY f(Y ) ¡ Xn i=1 ¸irY hi(Y ) + Xm i=1 ¹irY gi(Y ) = rY f(Y ) ¡ hTY ¸ + gT Y ¹; (4.19) where ¡gi(Y ) = si; i = 1; ¢ ¢ ¢ ; m; ¹isi = º; i = 1; ¢ ¢ ¢ ; m; hTY is a matrix consisting of the gradient rY hi(Y ) as columns, and gT Y is a matrix consisting of the gradient rY gi(Y ) as columns, or hTY = · rY h1(Y ) rY h2(Y ) ¢ ¢ ¢ rY hn(Y ) ¸ ; gT Y = · rY g1(Y ) rY g2(Y ) ¢ ¢ ¢ rY gm(Y ) ¸ : Also, ¸ is a vector containing the elements ¸i, and ¹ is a vector containing the elements 33 ¹i. The set of equations we must solve is rY f ¡ Xn i=1 ¸irhi + Xm i=1 ¹irgi = 0; (4.20) hi(Y ) = 0; i = 1; ¢ ¢ ¢ ; n (4.21) gi(Y ) + si = 0; i = 1; ¢ ¢ ¢ ;m (4.22) ¹isi ¡ º = 0; i = 1; ¢ ¢ ¢ ;m (4.23) si > 0; i = 1; ¢ ¢ ¢ ;m (4.24) ¹i > 0; i = 1; ¢ ¢ ¢ ;m (4.25) These equations can be solved using the Newton iterative method. Our four sets of variables for which we must solve are Y; ¸; ¹ and s. The equation for a ¯rst order approximation of a Taylor series of a function, F, whose independent variables are Y; ¸; ¹ and s is F(Y; ¸; ¹; s) »= F(Yo; ¸o; ¹o; so) + @F @Y 4Y + @F @¸ 4¸ + @F @¹ 4¹ + @F @s 4s: (4.26) We want to ¯nd the values of Y; ¸; ¹ and s where the expressions on the left side of the equations we want to solve evaluate to zero. We use NewtonRaphson to do so. Taking the ¯rst order approximation to the Taylor series for each of the four expressions given in equations (4.20)(4.23) and setting them equal to zero (the desired value for each) give us ¡ rY f ¡ hTY¸ + gT Y ¹ ¢ + (r2 Y f ¡ Xn i=1 ¸ir2 Y hi + Xm i=1 ¹ir2 Y gi)¢Y ¡hTY ¢¸ + gT Y¢¹ = 0; (4.27) h + hY¢Y = 0; (4.28) (g + s) + gY¢Y + I¢s = 0; (4.29) (MSe ¡ ºe) + S¢¹ +M¢s = 0; (4.30) 34 where I : an identity matrix; S : a diagonal matrix constructed from (s1; s2; ¢ ¢ ¢ ; sm); M : a diagonal matrix constructed from (¹1; ¹2; ¢ ¢ ¢ ; ¹m); e : a column vector with all elements 1. The above equations can be organized in matrix form as 2 6666666664 HY hTY gT Y 0 hY 0 0 0 gY 0 0 I 0 0 S M 3 7777777775 2 666666664 ¢Y ¢¸ ¢¹ ¢s 3 777777775 = 2 6666666664 ¡rY f + hTY ¸ ¡ gT Y ¹ ¡h ¡g ¡ s ºe ¡MSe 3 7777777775 ; (4.31) or W¢z = ¢F; (4.32) where HY = r2 Y f ¡ Xn i=1 ¸ir2 Y hi + Xm i=1 ¹ir2 Y gi: We initially set º to some relatively large number, such as 10. Starting with an initial guess zo = 2 66666664 Yo ¸o ¹o so 3 77777775 ; we calculate the W matrix and the ¢F vector. If all the elements of ¢F are su± ciently close to zero, we have found the solution zo. Otherwise we must solve for ¢z by ¢z = W¡1¢F: (4.33) 35 Then the original z is updated by z = zo + ®¢z: (4.34) However, following conditions need to be satis¯ed when updating z: ¹i > 0; si > 0: Therefore, when calculating the updated ¹ and s, we must make sure that each ¹i and each si is still strictly greater than zero when ¢¹i and ¢si are added to it, respectively. If adding ¢¹i or ¢si violates this condition, all ¢¹i and all ¢si must be scaled by some factor ® less than one before adding them as in (4.34). Now that we have a new guess for z, the W matrix and ¢F vector are calculated again. If the elements of the ¢F vector are su±ciently close to zero, we have found the z which solves our problem. Otherwise we solve for ¢z and update z again. This process is repeated until we ¯nd the z which makes ¢F very close to zero. After we solve for z, we reduce the barrier parameter º by some factor ·. For example, if · = 0:4, the barrier parameter º is updated by letting the new º be 0.4 times the old º. Using the value of z obtained using the old º as an initial guess zo, we use the iterative procedure described above to solve for the value of z which makes ¢F approximately zero for the new barrier parameter. The process of solving for z, decreasing º, and then solving for z again is repeated until º becomes a very small number, such as 10¡10. This process is illustrated in Figure 4.2. When º gets this small, we have found the z that solves our problem. When solving for z given a particular º, it is not necessary to force ¢F to be zero. What we really want to know is how close we are to the central path. The central path is de¯ned by a sequence of solutions fY (º); ¸(º); ¹(º); s(º)g which make ¢F evaluate to zero for every possible value of º. One way of measuring how close we are to the central path is by checking to see how close each product ¹isi is to the barrier parameter º. One 36 Figure 4.2: Graphical representation of central path proposed way of doing this is by ¯rst calculating the average value of ¹isi, also known as the dualityfactor [13], by bº = Xm i=1 ¹isi m ; (4.35) which evaluates to zero when each product ¹isi is equal to º. Instead of checking to see if ¢F is su±ciently close to zero, we check to see when the following is true: °°°°°°°°°°°°° ¹1s1 ¡ bº ¹2s2 ¡ bº ... ¹msm ¡ bº °°°°°°°°°°°°° 2 · ¿º (4.36) where 0 < ¿ < 1: When this logical statement becomes true, we are su±ciently close to the central path, and we can reduce º. Despite several attractive features of the PDIP method, it has an inherent disad vantage that the size of problem is much bigger than the Newton active set method since the PDIP method includes all inequality constraints. This problem becomes 37 even more involved when we consider a large scale natural gas and power system network. 38 CHAPTER 5 NATURAL GAS FLOW MODELING Natural gas is transported from gas producers to customers at various locations. A typical natural gas transmission system today consists of a large number of gas pro ducers, various customers, storage, many compressor stations, thousands of pipelines, and many other devices such as valves and regulators, including midstream gas pro cessing (between well and pipeline) A typical pipeline network for transporting natural gas consumes a signi¯cant amount of fuel per day to operate compressors pumping natural gas. This is because there is a gas pressure loss due to friction between gas and pipe inner walls. Moreover, energy is lost by heat transfer between gas and its environment. To compensate for these losses of energy and to keep the gas °owing, compressor stations are installed in the network, which consume a part of the transported gas resulting in economic losses1. The purpose of this chapter is to provide the underlying mathematical modeling of natural gas transmission networks. It will include gas °ow equations, compressor horsepower equation, and matrix representations of natural gas networks. 1Losses often refer to leakage of natural gas from the pipeline system, which is negligible if there is no theft problem. We will use the term \losses" only for the gas required to power the compressors. 39 Figure 5.1: Pipeline network representation 5.1 Elements of Natural Gas Transmission Net work Three basic types of entities are considered for the modeling of a natural gas trans mission network: pipelines, compressor stations, both of which are represented by branches, and interconnection points, represented by nodes. For simulation of a network, we assume that nodes represent pipeline connections while branches represent pipelines and compressor stations, which have °ow directions assigned. We also have di®erent types of nodes. A source node represents a gas production or storage facility. A load node represents a place where gas is to be taken out of the system, either for consumption or storage. Each compressor branch de¯nes two more node types technically referred to as suction and discharge nodes. Likewise, each pipeline branch has two types of nodes, a sending end node and a receiving end node. We de¯ne fkij (or fk for simple notation) as the °ow rate through a branch, say number k, which begins at node i and ends at node j. Figure 5.1 shows an example of a natural gas network. It is a simple treestructure pipeline network, which consists of one source at node 1, two loads at nodes T + 1 & K + 1, and several pipeline and compressor branches. wS1 is the amount of gas supplied at the source node, and wLT+1 & wLk+1 are the amounts of gas consumed at the load nodes, respectively. 40 For a given gas network, there are three types of relevant decision variables: the °ow rate through a pipeline, the °ow rate through a compressor, and the pressure at each node. For a compressor, these variables are further restricted by a set of constraints that depend on the operating attributes of the compressor. 5.2 Network Topology Analysis of natural gas pipeline networks is relatively complex, particularly if the network consists of a large number of pipelines, several compressors, and various sup pliers and consumers [21]. Matrix notation is a simple and useful way of representing a network. In gas network analysis, matrices turn out to be the natural way of ex pressing the problem [22]. A gas pipeline system can be described by a set of matrices based on the topology of the network. Consider the gas network represented by the graph in Figure 5.2 [22]. This network consists of one source at node 1, three loads at nodes 2, 3 & 4 and ¯ve pipeline branches ¬, , ®, ¯ & °. For network analysis it is necessary to select at least one reference node. Mathe matically, the reference node is also referred to as an independent node, and all nodal and branch quantities are dependent on it. The pressure at the reference node is Figure 5.2: Graph of the gas network [22] 41 known. A network may contain several pressurede¯ned nodes and these form a set of reference nodes for the network. A gas injection node is a point where gas is injected into the network, which may be positive, negative, or zero. Negative injection represents a load demand for gas from the network. This node may be supplying domestic or commercial consumers, charging gas storages, or even accounting for leakage in the network. A positive injection represents a supply of gas to the network. It may take gas from storage, source or another network. A zero injection is assigned to nodes that do not have a load or source but are used to represent a point of change in the network topology, such as a junction of several branches. Then, the vector of gas injections w at each node for Figure 5.2 is de¯ned as w = · w1 w2 w3 w4 ¸T ; = · wS1 ¡wL2 ¡wL3 ¡wL4 ¸T ; where wi = wSi ¡ wLi ; wSi = gas injection at node i; wLi = gas removed at node i: At steadystate conditions, the total load on the network is balanced by the supply into the network at the source node. To de¯ne the network topology completely, it is necessary to assign a direction to each branch. Each branch direction is assigned arbitrarily and is assumed to be the positive direction of °ow in the branch. If the °ow has a negative value, then the direction of °ow is opposite to the branch direction. References such as Osiadacz's book [22] may be consulted for further detail. 42 5.3 Matrix Representations of Network The interconnection of a network can be described by the branchnodal incidence matrix A. This matrix is rectangular, with the number of rows NN equal to the number of nodes (including reference nodes), and the number of columns NP equal to the number of pipeline and compressor branches in the network. The element Aij of the matrix A corresponds to node i and branch j, and is de¯ned as Aij = 8>>>>< >>>>: +1; if pipeline branch j enters node i; ¡1; if pipeline branch j leaves node i; 0; if pipeline branch j is not connected to node i: For the network in Figure 5.2, the branchnodal incidence matrix is A = 2 66666664 ¡1 ¡1 ¡1 0 0 1 0 0 1 0 0 1 0 ¡1 ¡1 0 0 1 0 1 3 77777775 ; where NN = 4 nodes and NP = 5 branches. One important thing to note is that the sum of all rows in the matrix A becomes zero due to the de¯nition of each element of the matrix A. Thus, the matrix A is always rank1 de¯cient, and therefore the matrix AAT is singular. This problem can be avoided by removing the rows of the matrix A corresponding to knownpressure nodes, which will be discussed in next chapter. 5.4 Flow Equation For isothermal gas °ow in a long horizontal pipeline, say number k, which begins at node i and ends at node j, the general steadystate °ow rate (in standard ft3/hr, or SCF/hr at T0 = 520oR and P0 = 14:65 psia) is often expressed by the following 43 formula [22, 23] derived from energy balance: fkij = Sij £ 3:22 T0 ¼0 s Sij ¡ ¼2 i ¡ ¼2 j ¢ D5 k FkGLkTkaZa ; (5.1) where fkij = pipeline °owrate (SCF/hr); Sij = 8>< >: +1 if ¼i ¡ ¼j > 0; ¡1 if ¼i ¡ ¼j < 0; Fk = pipeline friction factor; Dk = internal diameter of pipe between nodes (inch); G = gas speci¯c gravity (air=1.0, gas=0.6); Lk = pipeline length between nodes (miles); ¼i = pressure at node i (psia); ¼j = pressure at node j (psia); ¼0 = standard pressure (psia); T0 = standard temperature (±R); Tka = average gas temperature (±R); Za = average gas compressibility factor. There are several di®erent °ow equations in use in the natural gas transmission indus try [22, 24]. The di®erences are mainly due to the empirical expression assumed for the friction factor, Fk, and in some cases the rigorous consideration of the deviation of the behavior of natural gas from that of an ideal gas. The change of °ow changing from partial turbulence to full turbulence is referred as the transition region. The Reynolds number, a measure of the ratio of the inertia force on an element of °uid to the viscous force on an element, at which this transition occurs is dependent on the diameter of the pipeline and its roughness, and is typically about 107 [22]. The extent of the transition region is dependent on the system considered, and within this region 44 the frictional resistance depends on both the Reynolds number and the pipe charac teristics. However, in the fullyturbulent °ow region for highpressure networks, the friction factor Fk is strictly dependent on pipeline diameter [22], [25]. That is Fk = 0:032 D1=3 k : (5.2) Then, equation (5.1) becomes fk = fkij = Sij Mk q Sij ¡ ¼2 i ¡ ¼2 j ¢ ; SCF/hr (5.3) where Mk = ² 18:062 T0 D8=3 k ¼0 p GLkTkaZa ; ² = pipeline e±ciency: As indicated in equation (5.3), gas °ow can be determined once ¼i and ¼j are known for given conditions. Equation (5.3), known as Weymouth °ow equation, is most satisfactory for large diameter (¸ 10 inches) lines with high pressures [23]. 5.5 Compressor Horsepower Equation During transportation of gas in pipelines, gas °ow loses a part of its initial energy due to frictional resistance which results in a loss of pressure. To compensate the loss of energy and to move the gas, compressor stations are installed in the network. In general, the nature of the compressor work function is very complex and depends also on consideration such as the number of compressors running within the compressor station, how the compressor units are con¯gured (i.e. in series, in parallel, combina tion of both, etc.), physical properties of the compressor units, and type of compressor unit [22]. One of the most common con¯gurations implemented in practice is that of compressor stations consisting of identical centrifugal compressor units operating in parallel. Centrifugal compressors are versatile, compact, and generally used in the 45 range of 1,000 to 100,000 inlet ft3=min for process and pipeline compression applica tions [24]. In a centrifugal compressor, work is done on the gas by an impeller. Gas is discharged at a high velocity into a di®user. The velocity of gas is reduced and its kinetic energy is converted to static pressure. A key characteristic of the centrifugal compressor is the horsepower consumption, which is a function of the amount of gas that °ows through the compressor and the relative boost ratio between the suction and the discharge pressures. After empirical modi¯cation to account for deviation from ideal gas behavior, the actual adiabatic (zero heat transfer) compressor horsepower equation [23] at To = 60oF (= 520oR) and ¼o = 14:65 psia becomes Hk = Hkij = Bk fk "µ ¼j ¼i ¶Zki(®¡1 ® ) ¡ 1 # ; HP (5.4) where Bk = 3554:58 Tki ´k µ ® ® ¡ 1 ¶ ; fk = °ow rate through compressor (SCF/hr); ¼i = compressor suction pressure (psia); ¼j = compressor discharge pressure (psia); Zki = gas compressibility factor at compressor inlet, Tki = compressor suction temperature (±R); ® = speci¯c heat ratio (cp=cV ); ´k = compressor e±ciency: 5.6 Conservation of Mass Flow The mass°ow balance equation at each node can be written in matrix form as (A + U)f + w ¡ T¿ = 0; (5.5) 46 where A = a branchnodal incidence matrix; Uik = 8>>>>< >>>>: +1; if the kth unit has its outlet at node i; ¡1; if the kth unit has its inlet at node i; 0; otherwise. Tik = 8>< >: +1; if the kth turbine gets gas from node i; 0; otherwise. f = a vector of mass °ow rates through branches; w = a vector of gas injections at each node: In addition to the matrix A, which represents the interconnection of pipelines and nodes, we de¯ne the matrix U, which describes the connection of units (compressors) and nodes. The vector of gas injections w is obtained by w = wS ¡ wL; (5.6) where wS = a vector of gas supplies at each node; wL = a vector of gas demands at each node: Thus, a negative gas injection means that gas is taken out of the network. The matrix T and the vector ¿ represent where gas is withdrawn to power a gas turbine to operate the compressor. So if a gas compressor, say k, between nodes i and j, is driven by a gas¯red turbine, and the gas is tapped from the suction pipeline i, we have the following representation: Tik = +1; Tjk = 0 ; and ¿k = amount tapped: Conversely, if the gas were tapped at the compressor outlet, we would have Tik = 0 ; Tjk = +1; and ¿k = amount tapped: 47 Analytically, we will assume that ¿k can be approximated as ¿k = ®Tk + ¯TkHkij + °TkH2 kij ; (5.7) where Hk = Hkij is the horsepower required for the gas compressor k in equation (5.4). 48 CHAPTER 6 NATURAL GAS LOADFLOW The problem of simulation of a gas network with NN nodes in steady state, known as load°ow, is usually that of computing the values of node pressures and °ow rates in the individual branches for known values of NS source pressures (NS ¸ 1) and of gas injections in all other nodes. Gas load°ow analyses are required operationally whenever signi¯cant changes in demands or supplies are expected to occur. It is also used for system planning pur poses. For example, when a gas¯red generator is located in the gas network, we need to see whether the network has enough capability to carry the required amount of gas to the generator while satisfying various network constraints, such as pressure limits at each node and compressor operation limits. In this chapter, we state the load°ow problem in a general way, and construct a mathematical formation. We present a load°ow problem for a network without compressors. Then we introduce a general load°ow analysis with compressors. We formulate the gas load°ow problem in a way similar to the electric load°ow problem, and include the gas consumption rates at compressor stations. However, we omit gas storage to avoid inter  temporal linkageswe only attempt to solve the gas load°ow at a single timea \snapshot". 49 6.1 Load°ow Problem Statement The gas load°ow problem is stated below: ² Given a natural gas system described by a branchnodal incidence matrix A, a unitnodal incidence matrix U, a gas turbinenodal incidence matrix T, and given a set of gas injections except at the NS knownpressure sources (injections at these nodes initially unknown), and each unit's operating condition (such as the compression ratio, the °ow rate through the compressor, or the suction or discharge pressure), ² determine all other pressures, and calculate the °ow rates in all branches and the gas consumptions at compressor stations. Simply speaking, one of two quantities, nodal pressure ¼i and gas injection wi at each node, and one compressor operating condition are speci¯ed, and other values are to be determined. Speci¯ed quantities are chosen based on the following conditions: Nodes: ² KnownInjection Node: For a node i of this type, we assume that we know a gas injection wi, and the pressure ¼i is to be determined. Generally, source and load nodes, and junctions with no gas injections belong to this node. Elec trically, this is analogous to a \load bus." In fact, solving the load°ow problem with only this type of node is not in general possible. The ¯rst reason is that, in the °ow equation (5.1), the pressures never appear by themselves, but instead appear only as a squaredpressure di®erence of the form ¼2 i ¡ ¼2 j . Therefore, there are only NN ¡ 1 pressures which a®ect the load°ow. We therefore pick NS ¸ 1 nodes to provide reference pressures. These nodes are generally the external gas sources supplying our system. 50 Another reason is that with is that solving a load°ow for a network containing only knowninjection nodes would imply that we know the gas injections at every single node. In fact, we cannot mathematically specify all NN gas injections, as it may not be possible to ¯nd a solution to the load°ow equations. Specifying the injections at all nodes is the same as specifying the gas supplies to gas turbines driving gas compressors, which we cannot know until the load°ow is solved. Instead, we must pick at least one node, allowing the set of gas injection(s) to be whatever is required to solve the load°ow equations. Thus, we have to specify another node type: ² KnownPressure Node: Each is typically one of the source nodes, and the pressures of such nodes serve as references for all other pressures. We assume that we know f¼i; i = 1; : : : ;NSg, but we do not know the corresponding gas injections. Electrically this is analogous to a (possibly distributed) \slack bus." In addition to nodes, the other main components are branches, which connect the nodes. Branches: ² Pipelines: Pipeline °ow modelling has already been discussed in the previous chapter. Other than the physical characteristics of the pipeline, the only vari ables that the °ow fk = fkij on pipeline k depends on are the pressures ¼i and ¼j at the two ends of the pipeline. ² Compressors: The other key component we will model in a gas network is a compressor (also called a unit). The connection between the unit's inlet and outlet nodes is not de¯ned by the branchnodal incidence matrix A, but by U. The compression ratio between the compressor inlet and outlet, and the °ow rate through the compressor are governed by the horsepower equation (5.4), not by the °ow equation (5.1). Compressor data (other than the physical charac teristics of the compressor) can be speci¯ed in several ways [22] for compressor 51 k: relative boost Rk = Rkij = ¼j=¼i, or absolute boost ¼j ¡ ¼i, or mass°ow rate fkij . The inlet pressure ¼i or the outlet pressure ¼j could also be speci¯ed. 6.2 Load°ow without Compressors We will ¯rst consider a gas network with only pipelines to explain the Newtonnodal method. Assume that node 1 is the knownpressure node, and all other nodes are the knowninjection nodes. Since we do not have compressors, the unit's inletoutlet nodes are not de¯ned. Then the load°ow problem is described as follows: Given : ¼1 w2 : : : wNN Determine : w1 ¼2 : : : ¼NN It appears that we have NN ¡ 1 quantities known, and NN ¡ 1 quantities unknown. The set of nodal °ow equations that describes a gas network with only pipelines is given by ¹ w = A1 ¢ f(¹¼; ~¼); (6.1) where A1 = a branchnodal incidence matrix except the knownpressure nodes; ¹ w = ¹ ws ¡ ¹ wL; f = a vector of °ow rates through pipelines: We used ¹¼ to indicate the part of ¼ we know, and ~¼ for unknown part. Likewise for ¹ w. Flow equation (5.3) is used to get each pipeline branch °ow as fkij = Sij Mk q Sij ¡ ¼2 i ¡ ¼2 j ¢ ; SCF/hr (6.2) In the Newtonnodal method [22], an initial approximation is made to the nodal pressures. This approximation is then iteratively corrected until the ¯nal solution is 52 reached. At each iteration the righthand side of equation (6.1) is not equal to zero. The pressures are only approximations of their true values and the °ows calculated from these pressures are not balanced at each node. The imbalance at a node is the nodal error which is a function of all the nodal pressures (except the knownpressure nodes). The Newtonnodal method solves the set of equations (6.1) iteratively until the nodal °ow errors are small enough to be insigni¯cant. The iterative scheme for correcting the approximations to the nodal pressures is ~¼k+1 = ~¼k + ¢~¼k; (6.3) where k = the number of iterations. Term ¢~¼ is computed from the following equa tion: Jk ¢ ¢~¼k = ¡(A1 ¢ f(¹¼; ~¼) ¡ ¹ w)k; (6.4) where the matrix J is the nodal Jacobian matrix and is given by [22] J = ¡A1DAT 1¦1; (6.5) and D = diag µ fkij ¼2 i ¡ ¼2 j ¶ ; k = 1; : : : ;NP ¦1 = diag(~¼): Since the matrix A1DAT 1 is symmetric and positive de¯nite, and ¦1 is a diagonal matrix, NewtonRaphson algorithm can be implemented with a great computational e±ciency [26], [21]. 6.3 Load°ow with Compressors We assume that the pressures at the NS knownpressure nodes are known, and that the injections at the knowninjection nodes are speci¯ed. Also, some operating pa rameter for each compressor, say (for purposes of illustration) the relative boost 53 Rk = Rkij , is speci¯ed, and let's say there are NP branches in the system, of which NC are compressors. We can state the load°ow problem this way: Given : ¼1 ¢ ¢ ¢¼NS wNS+1: : :wNN R1¢ ¢ ¢RNC Find : w1 ¢ ¢ ¢wNS ¼NS+1 : : :¼NN f1 ¢ ¢ ¢ fNC fNC+1: : :fNP It appears that there are NN +NC quantities given, while NN +NP must be found. It is clear that NC < NP : This inequality is strict, unless we have no pipelines at all, only NP compressors connected to each other without any intervening pipelines { a silly situation. So since NN + NC < NN + NP , the system appears to be undetermined. Note, however, that from (6.2), the °ow fk depends only on the pressures ¼i and ¼j of the nodes it connects. Likewise, the horsepower Hk required by the compressor depends only on the °ow fk and the ratio Rk = ¼j=¼i, and therefore only depends on the pressures ¼i and ¼j . The tapo® loss ¿k depends only on Hk and thus on the nodal pressures. So, if we knew f¼i; i = 1; : : : ;NNg, we would know all other quantities we have discussed. But we only know them for i = 1; : : : ;NS. Let's use ¼ to indicate the part of ¼ we know, and ~¼ for the unknown part. Likewise for w and ~ w. The objective, then, is to calculate ~¼, giving us the entire vector ¼, from which all other quantities can be calculated. We can use the massbalance equation (5.5) to write w = ~T ¿ (¼; ~¼) ¡ (~A + ~U )f(¼; ~¼); where ~A , ~ U, and ~T are obtained from A, U, and T as in equation (5.5) after we have removed the rows corresponding to the NS nodes (the knownpressure nodes with unknown injections), that is, we have dropped the ¯rst NS equations of (5.5) 54 corresponding to ~ w. This gives us NN ¡ NS equations (for the elements of w) in NN ¡NS unknowns (the elements of ~¼). At this point, standard NewtonRaphson or other iterative methods can be employed to drive the \mismatch" ¢w (the di®erence between the values of w computed from above and the true values of w given as part of the data) to zero by correcting our current guess for the correct value of ~¼. We have shown that the dimension of the load °ow problem with compressors is NN ¡ NS. However, it is convenient to add new equations and new variables if they make the problem easier to implement. Since the compressor horsepower in equation (5.4) and the gas consumption rate in equation (5.7) are functions of horsepowers, compressor °ow rates, and gas consumption rates as well as the speci¯ed relative boost rates, we will add three extra decision variables for each compressor station. Then the new decision variables become x = · ~¼T HT ¿ T fT C ¸T ; where H = a vector of horsepowers at each compressor, ¿ = a vector of gas consumptions at each compressor, fC = a vector of compressor branch °ow rates. Since there are NN ¡ NS + 3NC decision variables, we need to de¯ne extra 3NC equations. The mass°ow balance equations F1 in matrix form are given by F1 = (~A + ~U ) ¢ 2 64 fC fP (¹¼; ~¼) 3 75 + ¹ w ¡ ~T ¢ ¿ = 0; (6.6) where fP (¹¼; ~¼) is a vector of mass °ow rates through pipelines, which are obtained by °ow equation (6.2). Note that the size of the vector F1 is NN ¡ NS. Assuming that the compression ratio Rkij is given, and gas compressor k located between nodes i and j is driven by a gas¯red turbine, the extra 3NC equations related with compressor 55 operation are F2k = ¿k ¡ ®Tk ¡ ¯TkHkij ¡ °TkH2 kij = 0; k = 1; : : : ;NC (6.7) F3k = Hkij ¡ Bk fCk "µ ¼j ¼i ¶Zki(®¡1 ® ) ¡ 1 # = 0; k = 1; : : : ;NC (6.8) F4k = µ ¼j ¼i ¶ ¡ Rkij = 0; k = 1; : : : ;NC (6.9) Let us de¯ne the vector of error functions F(x) by F(x) = 2 66666664 F1(x) F2(x) F3(x) F4(x) 3 77777775 : Note that there are NN ¡NS +3NC equations and NN ¡NS +3NC unknown decision variables. So we have the right number of variables to force all the mismatches to zero by using NewtonRaphson or other iterative methods. The iteration will be repeated until the mismatches become small enough to be insigni¯cant. 56 CHAPTER 7 UPFC MODELING FOR STEADYSTATE ANALYSIS Modeling UPFC for steadystate analysis has been considered by several researchers and an injection model [6] and an uncoupled model [14] have been proposed. These models can be easily incorporated into steadystate load°ow or optimal power °ow studies. However, these models employ four UPFC control variables that depend on the UPFC input and output currents and voltages and both models require adding two additional buses to the load°ow or OPF problem formulation. The voltage and current relationships between UPFC input and output need to be included explicitly. In addition, a constraint on real power conservation needs to be added, thereby reducing the degrees of freedom four to three. Since electricity prices at the UPFC input and output buses, which are the dual variables (or \shadow" prices) associated with real and reactive power injections, are meaningless when no power is bought or sold at these ¯ctitious buses, their addition to the problem only serves to increase the size of the problem. Therefore, these models are undesirable for our UPFC sensitivity analysis. To overcome these problems, a new steadystate mathematical model for a UPFC, the UPFC ideal transformer model, is proposed. In this model, UPFC control vari ables do not depend on UPFC input and output voltages and currents, and therefore addition of ¯ctitious input and output buses are not necessary. This model is eas ily combined with transmission line models using ABCD twoport representations, 57 which can then be converted to Yparameter representations. Thus, UPFC model is embedded in the ¹Y bus matrix, and so the size of the ¹Y bus matrix is not changed. 7.1 Operating Principles The Uni¯ed Power Flow Controller (UPFC) is one of the most technologically promis ing devices in the FACTS family [4, 5, 6, 7]. It has the capability to control voltage magnitude and phase angle, and can also independently provide (positive or negative) reactive power injections. A UPFC consists of a shunt transformer, a series trans former, power electronic switching devices and a DC link, as shown in Figure 7.1 [7]. Inverter 1 is functionally a static VAR compensator assuming that inverter 2 is not connected. It injects reactive power in the form of current at the shunt transformer, and the current phasor ~IT is in quadrature to the input voltage ~VI . Inverter 2 by itself represents the socalled advanced controllable series compensator (ACSC) assuming that inverter 1 is not connected. It injects reactive power by adding voltage through the series transformer. The injected voltage ~V T is in quadrature to the receiving end Figure 7.1: General UPFC scheme [7] 58 Figure 7.2: Proposed UPFC model in a transmission line current ~Io. Now if we connect inverter 1 to inverter 2 through a DC link, inverter 1 can provide real power to inverter 2. Therefore the UPFC can independently control real and reactive power injections through the series transformer, but the real power injected at the series transformer is provided by the shunt transformer through the DC link. Inverter 1 must provide the real power used by inverter 2 via the DC link, but can also independently inject reactive power (positive or negative) through the shunt transformer. In summary, note that the UPFC conserves real power but can still generate (or sink) reactive power at either transformer or both. 7.2 Uncoupled Model Figure 7.2 shows a basic UPFC model, where the UPFC is located between buses i and k. Each part of the transmission line is represented as an equivalent ¦ circuit. Figure 7.3 shows a phasor diagram illustrating UPFC inputoutput voltage and current re lationships. The injected series voltage ~V T can be resolved into inphase component Vp and quadrature component Vq with respect to the UPFC output current ~Io, and can be expressed as ~V T = (Vp + jVq) ej±o : (7.1) 59 Figure 7.3: Phasor diagram of UPFC inputoutput voltages and currents Since ~V T is dependent on the UPFC output current phase angle ±o, it requires adding an extra bus for the UPFC output terminal. The current IT injected by the shunt transformer contains a real component Ip, which is in phase or in opposite phase with the input voltage. It also has a reactive component Iq, which is in quadrature with the input voltage. Then the injected current ~IT can be written by ~IT = (Ip + jIq) ejµI ; (7.2) where µI is the UPFC input voltage phase angle. Thus, a second extra bus is required for the UPFC input terminal. The UPFC inputoutput voltage and current can be represented by ~V o = ~V I + ~V T = VIejµI + Vpej±o + jVqej±o ; (7.3) ~Io = ~II ¡ ~IT = IIej±I ¡ IpejµI ¡ jIqejµI ; (7.4) where ±I is the UPFC input current phase angle. Then, the complex power injected into the transmission line by the series transformer can be resolved into the real and reactive power in simple form as ST = ~V T ¢ ~I¤ o = V p{¢zI}o PT + jVq{¢zI}o QT : (7.5) 60 Figure 7.4: Uncoupled UPFC model in a transmission line. The inphase voltage Vp is associated with a real power supply and the quadrature voltage Vq with an inductive or capacitive reactance in series with the transmission line. Since the real power PT (which may be negative) is provided by the current Ip in the shunt transformer, we can derive the following relationship: Vp ¢ Io ¡ VI ¢ Ip = 0: (7.6) or the real power input equals the real power output. Due to (7.6), the number of degrees of freedom for the UPFC is reduced to three. Now that two extra buses for the UPFC input and output terminals are added, and the UPFC voltage and current relationships (7.3, 7.4) and real power °ow equation (7.6) are established, we can represent the UPFC uncoupled model as shown in Figure 7.4. The currents injected into the UPFC input and output buses are ~II = µ Yi 2 + 1 Zi ¶ ~V I ¡ 1 Zi ~Vi ; (7.7) ~Io = µ Yk 2 + 1 Zk ¶ ~V o ¡ 1 Zk ~ Vk : (7.8) The magnitudes of the injected voltage VT and current IT are limited by the maximum voltage and current ratings of the inverters and their associated transformers, which need to be included as inequality constraints in OPF. 61 7.3 Ideal Transformer Model Since the UPFC conserves real power, and generates or consumes reactive power, it can be modelled using an ideal transformer and a shunt branch, as shown in Figure 7.5 [27]. The advantage of this model is that the ideal transformer turns ratio and the variable shunt susceptance are independent variables, which are not directly as sociated with the UPFC inputoutput voltages and currents. We de¯ne the UPFC variables as follows: T = transformer voltage magnitude turns ratio (real); Á = phase shifting angle; ½ = shunt susceptance; and the ideal transformer turns ratio can be written by ¹ T = TejÁ: It is important to note that the ideal transformer does not generate real and reactive power, and the reactive power is generated (or consumed) by the shunt admittance only. Since the UPFC inputoutput voltage and current relationship can be expressed Figure 7.5: UPFC ideal transformer model 62 as ~V I = ~V oT\Á; (7.9) ~II = j ¹ T½~V o + 1 ¹ T¤ ~Io; (7.10) the UPFC can be represented by an ABCD matrix as 2 64 ~V I ~II 3 75 = ABCDU ¢ 2 64 ~V o ~Io 3 75 ; (7.11) where ABCDU = 2 64 ¹ T 0 j ¹ T½ 1 ¹ T¤ 3 75 : (7.12) Note that equation (7.11) is not bilateral unless ¹ T = 1\0. Now, we will show that this ideal transformer model represents the UPFC by comparing the complex power injections at the UPFC input and output. Using (7.11), the complex power injection at the UPFC input can be obtained by ¹ SI = ~V I ~I¤ I ; (7.13) = ¹ T~V o µ j ¹ T½~V o + 1 ¹ T¤ ~Io ¶¤ ; = ~V o~I¤ o ¡ jj ¹ Tj2 ¢ j~V oj2½; = ¹ So ¡ jj ¹ Tj2 ¢ j~V oj2½; and the real and reactive power injections can be obtained by PI = Real ¡ ¹ SI ¢ ; QI = Imag ¡ ¹ SI ¢ ; Po = Real ¡ ¹ So ¢ ; Qo = Imag ¡ ¹ So ¢ : Thus, we can derive the following relationships between the UPFC input and output: PI = Po; (7.14) Qo = QI + j ¹ Tj2 ¢ j~V oj2½: (7.15) 63 Figure 7.6: Simpli¯ed UPFC circuit Equations (7.14) and (7.15) mean that the ideal transformer model conserves real power and generates or consumes (for ½ < 0) reactive power. To determine how much real and reactive power is injected in the series and shunt transformers, we will map the complex turns ratio ¹ T in the ideal transformer and the shunt susceptance ½ to the injected voltage ~V T and current ~IT in the UPFC uncoupled model. Since the UPFC input voltage and current are expressed as ~V I = ~V o ¡ ~V T = ~V o Ã 1 ¡ ~V T ~V o ! = ~V oT\Á; (7.16) ~II = ~Io + ~IT = ~Io + µ 1 T \Á ¡ 1 ¶ ~Io + j½~V I ; (7.17) the injected voltage ~V T and current ~IT can be obtained by ~VT = (1 ¡ T\Á)~V o; (7.18) ~IT = µ 1 T \Á ¡ 1 ¶ ~Io + j½~V I : (7.19) Then, the power °ows through each inverter, as shown in Figure 7.6, can be obtained 64 by ¹ S1 = ~V I ~I¤ T ; = ~V oT\Á ·µ 1 T \Á ¡ 1 ¶ ~Io + j½~V oT\Á ¸¤ ; = (1 ¡ T\Á) ¹ So ¡ j½j ¹ Tj2j~V oj2; (7.20) ¹ S2 = ¡~V T ~I¤ o ; = (T\Á ¡ 1)~V o~I¤ o ; = (T\Á ¡ 1) ¹ So: (7.21) Thus, ¹ S1 + ¹ S2 = ¡j½j ¹ Tj2j~V oj2; (7.22) which veri¯es that the UPFC conserves real power and can generate (or consume) reactive power. Since the UPFC is modelled using passive circuit elements only, nonideal UPFC characteristics, such as shunt and series transformer reactances can be easily incor porated into this framework. 7.4 UPFC in a Transmission Line A twoport ABCD matrix is the most convenient method to represent cascaded net works [9]. Let us divide a transmission line between buses i and k with a UPFC into three cascaded networks, a UPFC input transmission line, a UPFC, and a UPFC output transmission line, as shown in Figure 7.7. The UPFC input transmission line , and the UPFC output transmission line are easily represented by twoport ABCD matrices since the transmission lines are modelled using ¦ equivalent circuits. We call ABCDi and ABCDk as the ABCD matrices for each transmission line, and de¯ned 65 Figure 7.7: Cascaded transmission line with a UPFC by ABCDi = 2 64 Ai Bi Ci Di 3 75 and ABCDk = 2 64Ak Bk Ck Dk 3 75 ; where each element is de¯ned by Ai = Di = 1 + YiZi 2 ; Bi = Zi; Ci = Yi(1 + YiZi 4 ); Ak = Dk = 1 + YkZk 2 ; Bk = Zk; Ck = Yk(1 + YkZk 4 ): The ABCD parameters of each transmission line can be obtained after we identify the propagation constant ° and the characteristic impedance Zc. Since we are using IEEE test cases with no knowledge of ° or Zc, we compute these values using the expressions ° = 1 l cosh¡1(1 + Y Z 2 ); (7.23) Zc = Z sinh(°l) ; (7.24) where l is the distance between buses i and k measured in kilometers. Then, assuming the UPFC is installed in x (0 < x < 1), ¦ equivalent circuit values for each section 66 of the transmission line can be found by Zi = Zc sinh(°l ¢ x); Yi = 2 Zi (cosh(°l ¢ x) ¡ 1) ; Zk = Zc sinh(°l ¢ (1 ¡ x)); Yk = 2 Zk (cosh(°l ¢ (1 ¡ x) ¡ 1) : Now, the three cascaded networks are combined to obtain 2 64 ~V i ~Ii 3 75 = ABCDi ¢ ABCDU ¢ ABCDk 2 64 ~V k ¡~Ik 3 75 ; (7.25) = 2 64 Aik Bik Cik Dik 3 75 2 64 ~V k ¡~Ik 3 75 ; where ABCDU is given in (7.12). So Aik = ¹ TAiAk + j ¹ TBiAk½ + 1 ¹ T¤BiCk; Bik = ¹ TAiBk + j ¹ TBiBk½ + 1 ¹ T¤BiDk; Cik = ¹ TCiAk + j ¹ TDiAk½ + 1 ¹ T¤DiCk; Dik = ¹ TCiBk + j ¹ TDiBk½ + 1 ¹ T¤DiDk: By rearranging (7.25) and solving for Ii and Ik, we have 2 64 ~Ii ~Ik 3 75 = ¹Y busik 2 64 ~V i ~V k 3 75 ; (7.26) where ¹Y busik = 2 64 Dik Bik Cik ¡ AikDik Bik ¡ 1 Bik Aik Bik 3 75 : In general, det 0 B@ 2 64 Aik Bik Cik Dik 3 75 1 CA = 1\2Á; (7.27) 67 If Á = 0, that is, complex ¹ T is real, this determinant is one at angle zero, and complex ¹Y busik becomes symmetrical. Note that equation (7.25) represents a bilateral two port network only if ¹ T = 1\0. As seen in (7.26), since the UPFC is embedded in the ¹Y bus matrix, the size of the ¹Y bus matrix is not changed, so UPFC sensitivity analysis can be performed using this ideal transformer model. 68 CHAPTER 8 GAS AND ELECTRICITY OPTIMAL POWER FLOW The purpose of this chapter is to construct a mathematical formulation to solve natural gas and electricity optimal power °ow problems. To do this we present a simple combined gas and electric network, and synthesize an optimal gas °ow (OGF) and an optimal electric power °ow (OPF) into a single natural gas and electricity optimal power °ow (GEOPF) problem. Though this integration of gas and electric system is based upon deterministic prices of gas in source nodes, we will certainly need this GEOPF algorithm for stochastic cases for future studies where the price will be a stochastic variable. For our purposes, gas and electricity storage are neglected. 8.1 Gas and Electric Combined Network For an integrated gas and electric network, even though the gas network and electric network are physically overlapped, we represent the two systems separately. One example is shown in Figure 8.1. Electric generator buses which are coincident with any gas nodes can be used to integrate the gas and electric networks. The generators in the combined nodes are assumed to be driven by gaspowered turbines. 69 Figure 8.1: Combined natural gas and electricity network 8.2 Gas and Electricity Optimal Power Flow The mathematical formulation of the GEOPF can be expressed as min Y C(Y ) ¡ B(Y ) (8.1) subject to: hi(Y ) = 0; i = 1; : : : ; n (8.2) gj(Y ) · 0; j = 1; : : : ;m (8.3) where ² C(Y ) is total cost of system operation, and B(Y ) is the total bene¯t to society and is treated as a negative cost. SW = B(Y )¡C(Y ) is the social welfare, which is the total bene¯t B(Y ) to society, less the cost of combined system operation, and thus is the negative of the net cost of system operation to society. So minimizing net cost C(Y ) ¡ B(Y ) is the same as maximizing social welfare. 70 ² Y is a vector of decision variables (i.e. the voltage magnitude and angle at each bus, real and reactive power generations, real and reactive power consumptions, nodal pressures, °ow rates through pipelines, °ow rates through compressors, gas consumptions by gas turbines, etc.). ² h is a set of equality constraints, such as the real and reactive power balance at each bus, pipeline branch °ow rates, the mass °ow balance at each node, etc.. ² g is a set of inequality constraints, such as line °ow limits, voltage limits, gen eration limits, relative boost limits, nodal pressure limits, etc.. Then the Lagrangian for the GEOPF problem can be written L(Y ) = ¡SW(Y ) ¡ Xn i=1 ¸ihi(Y ) + Xm j=1 ¹jgj(Y ): (8.4) Here, the objective of the GEOPF is to maximize the social welfare, while satisfying n equality constraints and m inequality constraints. 8.2.1 Cost, Bene¯t, and Social Welfare In the \Poolco" environment, central dispatch sets the price of gas or electricity in the system to optimize a networkwide nondiscriminatory service, while consumers and producers contract competitively [8, 1]. The type of problem that the system operator wishes to solve is an optimal power °ow (OPF) problem in electricity and an optimal gas °ow (OGF) problem in natural gas, but doing so simultaneously in an integrated manner. The objective function to be maximized that we will use in the GEOPF problem is social welfare [8]. To understand the concept of the social welfare, we introduce the fact that there are two types of market participants, producers and consumers. A producer is a supplier of some type of goods. Of course, there is a cost to produce the good. Let's say we know the cost of producing the good as a function of the amount 71 of the good that is produced. In our case, the electricity producer is supplying electric power, so its cost is a function of the amount of real power PG it produces. The total cost of real power generation is then (assuming quadratic cost curves) CE = XNG i=1 i=2G CEi(PGi); (8.5) = XNG i=1 i=2G £ ®Gi + ¯GiPGi + °GiP2G i ¤ ; where NG = total number of real power generators; ®Gi ; ¯Gi ; °Gi = cost coe±cients of real power generator i; G = electric nodes with gas¯red generators: For the combined node, it is important to note that the electric generation cost is not included in the equation (8.5). This is because the combined nodes simply transform gas energy into electric energy. For the gas supplier, its production cost is a function of the amount of gas pro duced. Then, the total production cost becomes CG = XNS i=1 ciwSi106 £ GHV; (8.6) where NS = total number of source nodes; ci = gas production cost at source node i ($/MMBTU), GHV = gas gross heating value (BTU/SCF): A consumer purchases a good because he receives some bene¯t from using the good. Let us say that this bene¯t that he/she receives can be measured in dollars and is a quadratic function of the amount of the goods purchased. For the electric 72 power consumer, the bene¯t BE in dollars is a function of the amount of real power PL consumed. The total bene¯t received from the consumption of real power (assuming quadratic bene¯ts) by loads is obtained by BE = NXEL i=1 BEi(PLi); (8.7) = NXEL i=1 £ ¯ELiPLi + °ELiP2 Li ¤ ; where NEL = total number of electric consumers; ¯ELi ; °ELi = bene¯t coe±cients of real power consumer i: For a gas consumer, his bene¯t BG is a function of the amount of gas consumed. Then, the total bene¯t received from the consumption of gas is expressed as BG = NXGL i=1 i=2G BGi(wLi); (8.8) = NXGL i=1 i=2G £ ¯GLiwLi + °GLiw2L i ¤ ; where NGL = total number of gas consumers; ¯GLi ; °GLi = bene¯t coe±cients of gas consumer i: Note again that the gas consumer's bene¯t at the combined node is not included since this node is simply converting gas energy to electric energy. Now we de¯ne the meaning of social welfare. It is the total bene¯t received by society due to the production and consumption of the good. In the case of a combined gas and electric network, the social welfare is de¯ned as SW = BE + BG ¡ CE ¡ CG (8.9) = NXEL i=1 BEi(PLi) + NXGL i=1 i=2G BGi(wLi) ¡ XNG i=1 i=2G CEi(PGi) ¡ XNS i=1 CGi(wsi); 73 and it is the objective function of the GEOPF which we want to maximize, or the negative of that which we wish to minimize. 8.2.2 Constraints The GEOPF is a constrained optimization problem, which must satisfy various equal ity and inequality constraints at the optimum. ² Equality Constraints { Electric Network Kircho®'s laws must be satis¯ed at each electric bus, resulting in the stan dard power °ow (power balance) equations at each bus PGi ¡ PLi = Pi(V; µ); i = 1; : : : ;NB (8.10) QGi ¡ QLi = Qi(V; µ); i = 1; : : : ;NB (8.11) The above two equations are conventional power °ow equations, and tell that the real and reactive power injection at each bus is a function of the system bus voltage magnitudes and phase angles, or Pi(V; µ) and Qi(V; µ). In our notation, V is the vector of system bus voltage magnitudes, and µ is the vector of bus voltage phase angles. The amount of reactive power consumed by the load, QL, is not considered a variable since it will be assumed that the load has a ¯xed power factor where the power factor is de¯ned as p:f: = p PL P2 L + Q2 L : (8.12) This assumption is not critical to our formulation, and is made for simplic ity and concreteness only. Then QL is a linear function of the real power 74 consumed, or QL(PL) = PL p:f: q 1 ¡ p:f:2 : (8.13) The real and reactive injection at bus i can be obtained by Si(V; µ) = ~V i~I¤ i ; i = 1; : : : ;NB (8.14) Pi(V; µ) = Re( ¹ Si); Qi(V; µ) = Im( ¹ Si); where the vector of current injections at each bus is obtained by ~I = ¹Y bus ¢ ~V : { Gas Network There are two sets of steadystate network °ow equations, which must be satis¯ed as equality constraints. The ¯rst equality constraint says that the amount of gas at node i transported from and to other nodes must be the same as that of gas injected at that node. This is known as a mass °ow balance equation. The set of mass °ow equations that describe a gas network is given by (A + U)f + w ¡ T¿ = 0; (8.15) where units of f, w and ¿ are SCF/hr, hence equivalent to conversation of mass. Note that the °ow rates to the compressor inlet and from the compressor outlet are di®erent if the compressor is driven by a gas¯red turbine. ¿k is the amount of gas supplied to gas¯red turbine k measured in SCF/hr. For simplicity, we assume ¿k = ®Tk + ¯TkHk + °TkH2 k ; k = 1; : : : ;NC (8.16) where Hk is the horsepower required for compressor k, and de¯ned by Hk = Bk fk "µ ¼j ¼i ¶Zki(®¡1 ® ) ¡ 1 # : k = 1; : : : ;NC (8.17) 75 The other equality constraint which is described by equation (5.3) is the pipeline gas °ow rate through each pipeline, and it is stated by ¼2 i ¡ ¼2 j = 1 M2 k Sij f2 kij ; i; j = 1; : : : ;NN (8.18) Note that while the nodal °ow equation (8.15) is linear, the pipeline branch °ow equation (8.18) is quadratic, and we need to take into account the fact that the °ow direction must be explicitly accounted for (by Sij . { Gas and Electricity Combined Node This is the location of an actual or planned gas¯red generation facility, which would be located where a gas node and an electric power bus have a common location. The gas node may have other loads or input injections, as may the electrical node the generator is connected to. Let us suppose that the electrical bus in question is i, which is fed gas from gas node ji. Then wji = wji,other (8.19) ¡ ¡ aGi + bGiPGi + cGiP2G i ¢ 1 GHV ; i 2 G where GHV = gas gross heating value (BTU/SCF); aGi ; bGi ; cGi = gas fuel rate coe±cients at node i, G = electric nodes with gas¯red generators. ² Inequality Constraints { Electric Network The system operating constraints include maximum and minimum limits 76 on the voltage magnitude at each bus and the magnitude of the line current °owing on each line. These can be written as Vmini · Vi · Vmaxi ; i = 1; : : : ;NB (8.20) for the voltage magnitude at bus i, and I2 ikmax ¸ ¯¯¯ ~Iik ¯¯¯ 2 = ~Iik ¢ ~I¤ ik; i; k = 1; : : : ;NB (8.21) where ~Iik = ³ ~V i ¡ ~V k ´ ¹yik + ~V i¹yshuntik for the magnitude of the line current from bus i to bus k. Equation (8.21) is called a transmission thermal limit. Since the voltage magnitude at each bus is changed, we use the line current instead of the apparent power °owing on each line for the thermal limit, although an apparent power constraint could be used instead. Also, if voltage drop or the steady state stability limits real power °ow, these constraints could be used instead. There are also limits on the amount of real and reactive power produced by each generator. There may also be a limit on the amount of real power consumed by each load. This gives us an additional set of inequality con straints. PGmini · PGi · PGmaxi ; i 2 NPG (8.22) QGmini · QGi · QGmaxi ; i 2 NQG (8.23) PLmini · PLi · PLmaxi ; i 2 NEL (8.24) The tap magnitude and the tap angle of a regulating transformer, the volt age magnitude and the angle drop along a transmission line, MW inter change transactions, shunt reactors or capacitors, etc. can be also included as inequality constraints if required [11]. 77 { Gas Network The system operating constraints include maximum and minimum limits on the pressure at each node and maximum relative boost limits for each compressor. These can be written as ¼mink · ¼k · ¼maxk ; k = 1; : : : ;NN (8.25) for the pressure ¼i at node i, and 1 · µ ¼j ¼i ¶ k · Rmaxk ; k = 1; : : : ;NC (8.26) for the relative boost rate between the inlet pressure ¼i and the outlet pressure ¼j at compressor station k. There are also limits on the amount of gas supplied at each source node. There may also be limits on the amount of gas consumed by gas consumers. This gives us an additional set of inequality constraints. wminSi · wSi · wmaxSi ; i = 1; : : : ;NS (8.27) wminLi · wLi · wmaxLi ; i = 1; : : : ;NGL (8.28) 78 CHAPTER 9 OPTIMAL LOCATION OF A UPFC IN A POWER SYSTEM This chapter presents a screening technique for greatly reducing the computation involved in determining the optimal location and estimation of the generation cost saving when a UPFC is installed and operated optimally. An technique to obtain the ¯rst and secondorder sensitivities of the generation cost with respect to UPFC con trol parameters is described. The sensitivity analysis attempts to estimate the UPFC value without having to run the OPF with the UPFC several times. This sensitivity technique can be used to optimally locate a UPFC in a large power system by ignoring the transmission lines with low marginal value (MV) and low estimated incremental value (IV), and running a full OPF only for the lines with higher estimated IV to obtain actual IV. Thus, this technique can greatly reduce the computational burden of determining the optimal location of the UPFC in a large power system. 9.1 Optimal Power Flow with UPFC Suppose that a UPFC is installed in transmission line ik. The mathematical formu lation of the OPF with the UPFC can be expressed as min y;xik C(y; xik) (9.1) 79 subject to: hi(y; xik) = 0; i = 1; : : : ; n (9.2) gj(y; xik) · 0; j = 1; : : : ;m (9.3) where ² C(y; xik) is the total generation cost. ² y is a vector of decision variables. ² xik = [ Tik Áik ½ik ]T is a vector of the UPFC control variables in line ik. ² fhi : i = 1; : : : ; ng is the set of equality constraint functions. ² fgj : j = 1; : : : ;mg is the set of inequality constraint functions. We use the UPFC ideal transformer model to construct the equations for OPF with a UPFC. It is important to note that the number of equality constraints is the same as that of the base case OPF with no UPFC. This is because the UPFC control variables do not depend on UPFC input and output voltages and currents, and the UPFC model is embedded in the ¹Y bus matrix, and because we ignore UPFC operation limits. Now, let us construct the Lagrange for the OPF problem as Lo(y; ¸; ¹; xik) = C(y; xik) + Xn i=1 ¸ihi(y; xik) + Xm j=1 ¹jgj(y; xik); (9.4) where ¸i and ¹j are the Lagrange multipliers for the equality and inequality con straints, respectively. To solve the proposed OPF problem with inequality con straints, we use the primaldual interiorpoint method. At the optimum, the last term of (9.4) must satisfy the complementary slackness condition such that ¹jgj = 0 for each j = 1; : : : ;m. Therefore, if an inequality constraint is binding in (9.4), we could treat it as an equality constraint, and we could ignore it if it is not binding. 80 Since we are using interiorpoint methods { not activeset methods { we do not have to distinguish between active and inactive constraints until the OPF problem is solved [20]. This avoids \cycling" behavior in the active set associated with \activeset" methods such as Newton's. Then, to derive the ¯rstorder sensitivities, we rewrite (9.4) as Lo(y; ¸; xik) = C(y; xik) + X j2A ¸jhj(y; xik); (9.5) where A is the set of active constraints, which are known once the basecase OPF is solved. 9.2 FirstOrder Sensitivity Analysis We consider the case where the UPFC is inserted in line ik, and the UPFC ideal transformer model is used for the analysis. The marginal values(MVs) of the UPFC, to be installed in line ik, are simply the amounts by which the total cost of system operation could be changed by allowing a small change of the UPFC control variables in line ik. We can obtain the MVs by assuming that there is a UPFC in line ik, but that the UPFC is not operating. So we add three extra constraints Tik = T; Áik = Á; ½ik = ½ to the original OPF problem, and for simplicity, we denote the constraints as xik = x. Then, the new Lagrangian can be written Lik(y; ¸; xik; ¸x) = C(y; xik) + X j2A ¸jhj(y; xik) + ¸Tx (x ¡ xik); (9.6) where ¸x = [ ¸T ¸Á ¸½ ]T : We de¯ne the function ¸¤ x(x) to be the optimal value of the Lagrange multiplier on the constraint xik = x. Here, we are most interested in ¸¤ x(x = x0), which is associated 81 with the constraints: Tik = 1; Áik = 0; ½ik = 0: This is because the OPF problem when solved with the UPFC control parameters x = x0 yields the same result for y and ¸ as the base case where there is no UPFC in line ik. Using the ¯rstorder conditions for the solution of the OPF problem for xik = x0; @Lik @xik = 0; (9.7) we can solve for ¸¤ x(x0) to obtain ¸¤ x(x0) = " @C(y¤; xik) @xik + X j2A ¸¤ j @hj(y¤; xik) @xik #¯¯¯¯¯ xik=x0 : (9.8) Note that @C(y¤; xik) @xik and ¸¤ j @hj(y¤; xik) @xik ¯¯¯¯ xik=x0 are easy to compute. Equation (9.8) indicates that the marginal value ¸¤ x(x0) can be determined once we know y¤ and ¸¤, which are obtained from the base case OPF with no UPFC. Thus, if we know y¤ and ¸¤, we can obtain the ¯rstorder sensitivities of cost with respect to UPFC control variables x for each possible transmission line by solving only the basecase OPF. Now, we will give the desired interpretation for the Lagrange multiplier ¸x by essentially rederiving the envelope theorem for our case. The Lagrangian in (9.6) can be rewritten as Lik(z; ¸z) = C(z) + ¸Tz hz(z); (9.9) where z = [ y xik ]T , ¸z = [ ¸ ¸x ]T and hz = [ h x ¡ xik ]T . Let the solution of the OPF with the UPFC be denoted by C¤(x), and let the optimal values of z and ¸z be denoted by z¤(x) and ¸¤z (x), respectively. Then C¤(x) = Lik(z¤(x); ¸¤ z(x); x) . (9.10) 82 We di®erentiate C¤(x) with respect to x to obtain d dx C¤(x) = d dx Lik(z¤(x); ¸¤ z(x); x) (9.11) = d dx £ C(z¤(x); x) + ¸¤T z (x)hz(z¤(x); x) ¤ = · @ @z C + ¸¤T z (x) @ @z hz ¸  {z } 0 dz¤(x) dx + £ hTz (z¤(x); x) ¤  {z } 0 d¸¤z (x) dx + @C @x (z¤(x); x) + ¸¤T z (x) @hz @x (z¤(x); x) = @C @x (z¤(x); x) + ¸¤T z (x) @hz @x (z¤(x); x); where the two terms in square brackets are equal to zero by virtue of the KKT ¯rst order optimality conditions. Equation (9.11) indicates that the change in the system generation cost results from the change in C due to the change in x and the change in hz due to the change in x, as the envelope theorem says. The interpretation of ¸x can be obtained by considering the constraint (x ¡ xik) and writing C(z¤(x)) = C(z¤(x); x) (9.12) and h(z¤(x); x) = h (z¤(x)) + [ 0 ¢ ¢ ¢ 0 x ¡ xik ]T : (9.13) Then the envelope theorem [8, 28] says that d dx C¤(x) = ¸¤T x (x): (9.14) For any value of the parameter x, the interpretation of the multiplier ¸¤ x(x) is the marginal change in the total generation cost as the constraint xik = x is changed slightly. If such a change is allowed, it would be made to improve the objective function. The sign of ¸¤ x determines the desired direction of the change in components of xik. Therefore, only the absolute value of the multiplier matters. As we vary x from x0 to the set of optimal values x¤ ik, the multiplier will vary from the marginal 83 value ¸¤ x(x0) to the value ¸¤ x(x¤ ik) = 0. This is because, if we constrain xik to be equal to its unconstrained optimal value, the multiplier for this constraint must be zero. Since the UPFC model is embedded in the ¹Y bus matrix, as explained in (7.26), the ¯rstorder UPFC sensitivity analysis is only associated with complex power in jections at buses i and k, and the thermal limit of transmission line ik if it is binding. If transmission °ow is limited by steadystate stability, such constraints (jPikj · Pmax; or jµi ¡ µkj · µmax) can be included as well. 9.3 SecondOrder Sensitivity Analysis The secondorder sensitivity we wish to obtain is Hx = d2C¤(xik) dx2 ¯¯¯¯ xik=xo = d¸¤ x(xik) dx ¯¯¯¯ xik=xo : (9.15) As in the ¯rst order sensitivity analysis, we arti¯cially insert the UPFC between bus i and bus k, but then constrain xik = x, where we are primarily interested in the case where x = x0, since it represents the base case OPF with no UPFC installed. Rather than including the additional control variable xik in the vector of decision variables y, and the additional multipliers ¸x in the vector of multipliers ¸, we keep these two variables separate for clarity. By the ¯rstorder conditions for optimality, the Lagrange function (9.6) must satisfy the following condition: !(u; x) = ruLik (9.16) = 2 6666666664 ryLik r¸Lik rxikLik r¸xLik 3 7777777775 = 2 6666666664 ryLik h(y; xik) rxikLik (x ¡ xik) 3 7777777775 = 0; 84 where u = 2 66666664 y ¸ xik ¸x 3 77777775 : Note that x occurs only in the term ¸Tx (x ¡ xik) of the Lagrangian in (9.6), so that mixed partials involving y and x, ¸ and x, and xik and x are zero. Thus we have @! @x = 2 66666664 0 0 0 +I 3 77777775 ; (9.17) where I is an identity matrix with the size of xik. As the vector x of the UPFC parameters changes, so will the optimal value of u. Let us denote the optimal value by u¤(x). Then for each x, it must satisfy the condition (9.16): !(u¤(x); x) = 0; 8x (9.18) Since optimality is maintained as we change x, it must satisfy d dx !(u¤(x); x) = 0: (9.19) If we expand (9.19) by the chain rule, then we have @! @u ¢ du¤(x) dx + @! @x = 0. (9.20) The ¯rst term of (9.20) can be expressed as @! @u = Wik(u¤(x); x) (9.21) = 2 6666666664 Hy JT y ryrTx ik Lik 0 Jy 0 Jxik 0 rxikrTy Lik JT xik Hxik ¡I 0 0 ¡I 0 3 7777777775 ; 85 where Hy and Hxik are the Hessians of the Lagrangian, and Jy and Jxik are the Jacobians of constraints in the active set A with respect to y and xik, respectively. Using Wik(u¤(x); x), we can rewrite (9.20) in simple form as Wik du¤(x) dx = ¡ @! @x : (9.22) Then, (9.22) becomes Wik ¢ 2 66666666664 dy¤(x) dx d¸¤(x) dx dx¤ ik(x) dx d¸¤ x(x) dx 3 77777777775 = 2 666666666664 0 0 0 ¡I 3 777777777775 : (9.23) In equation (9.23), we can see that dx¤ ik=dx = I, an obvious result since we have constrained xik = x. The ¯rst two equations can be combined and rearranged to obtain 2 64 dy¤(x) dx d¸¤(x) dx 3 75 = ¡W¡1 ¢ 2 664 ryrTx ik Lik Jxik 3 775 ; (9.24) where W = 2 64 Hy JT y Jy 0 3 75 : Once this is found, the third equation can be rearranged to yield d¸¤ x(x) dx ¯¯¯¯ x=x0 = · rxikrTy Lik JT xik ¸ ¢ 2 64 dy¤(x) dx d¸¤(x) dx 3 75 ¯¯¯¯¯¯¯ x=x0 + Hxik : (9.25) The matrix W is not a®ected by the UPFC since we set x = x0. Thus, to calculate d¸¤ x=dx for x = x0, we already have available W matrix in factored form, so only single forward substitution and back substitution are required in (9.24). 86 9.4 Estimation of Incremental Value using First and SecondOrder Sensitivities Now, we will estimate the incremental value of the UPFC using the ¯rst and second order sensitivities. We approximate ¸x(x) as a linear function of x. Let us expand C(z¤(x); x) at the optimum with respect to x. Then, we have C(z¤(x); x) = C(x) (9.26) = C(x0) + ¸¤T x (x0)¢x + 1 2 ¢xTHx(x0)¢x + H:O:T; where ¢x = x ¡ x0: Neglecting the high order terms in (9.26) and di®erentiating with respect to x yield rxC(x) »= ¸¤ x(x0) + Hx(x0)¢x: (9.27) Since rxC(x) = 0 at the optimum, the estimated optimal value of x is obtained as ^x¤ = ¡H¡1 x (x0)¸¤ x(x0) + x0: (9.28) The estimated incremental value (cIV ) at the estimated optimum ^x¤ is then calculated as cIV = C(x0) ¡ C(^x¤); (9.29) = ¡¸¤T x (x0) (^x¤ ¡ x0) ¡ 1 2 (^x¤ ¡ x0)T Hx(x0) (^x¤ ¡ x0) ; = ¸¤T x (x0)H¡1 x (x0)¸¤ x(x0) ¡ 1 2 ¸¤T x (x0)H¡T x (x0)Hx(x0)H¡1 x (x0)¸¤ x(x0); = 1 2 ¸¤T x (x0)H¡1 x (x0)¸¤ x(x0) ¸ 0: Although IV is by de¯nition nonnegative, nonnegativity of cIV is guaranteed only Hx is positive semide¯nite (which is not guaranteed). The accuracy of cIV depends also on the fact that A does not change in (9.5) since this would change (9.16), the fundamental equation for the sensitivity analysis. 87 CHAPTER 10 RESULTS 10.1 Result of Sensitivity Analysis The proposed sensitivity methods are tested on a 5bus system derived from the IEEE 14bus system, and IEEE 14 and 30bus systems to establish their e®ectiveness. These systems have 7, 20, 41 transmission lines, respectively. Figure 10.1 shows the 5bus system [15]. The system consists of two generators at buses 1 & 2 and one synchronous condenser at bus 3. We assume that generator 2 has higher generation cost than generator 1, and UPFCs are installed in the middle of each transmission line. Loads are assigned such that the current °ow constraint in line 1 is binding (we assume that thermal constraints limit line °ow for this example). The marginal values (j¸T j; j¸Áj and j¸½j), estimated incremental values (cIV ) and Figure 10.1: Diagram of 5bus subset of IEEE 14bus system [15] 88 1 2 3 4 5 6 7 0 2 4 6 8 10 12 Line number MV / Actual IV lT/15 lf/15 lr/5 Estimated IV Actual IV ($/hr) Figure 10.2: Normalized marginal and incremental values for 5bus system. actual incremental values (IV ) for the 5bus system are shown in Figure 10.2. It shows that the lines with high MVs usually produce high estimated IVs. The UPFC location in line 1 produces high MVs for j¸T j and j¸Áj, and the highest estimated IV. In general, it is more economical to locate the UPFC in heavilyloaded high voltage PG1 PG2 Ploss Est. IV Actual IV UPFC Location MW MW MW $/hr $/hr Without UPFC 198.18 71.60 10.51 0 0 Line 1 219.31 51.18 11.64 11.85 6.90 Line 2 224.56 46.19 11.90 7.69 6.38 Line 3 196.85 72.77 10.77 1.13 0.73 Line 4 202.45 67.40 10.99 1.97 1.83 Line 5 204.52 65.44 11.11 2.17 2.07 Line 6 196.85 72.73 10.73 0.75 0.67 Line 7 202.59 67.22 10.95 0.78 0.80 Table 10.1: Real power generation, line loss and total generation cost for 5bus system 89 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 Line number MV lT/15 lf/15 lr/5 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 Line number IV Estimated IV ($/hr) Actual IV ($/hr) Figure 10.3: Marginal and incremental values for IEEE 14bus system. line 1 since it allows more power to be transmitted in underutilized line 2 while preventing line 1 from overloading. Eventually, no further savings due to UPFC operation can be achieved because of the voltage constraint at bus2. Also, reasonable agreement can be observed between the estimated IVs and the actual IVs except line 1 and line 2. The di®erence between the estimated IV and the actual IV is attributed to the fact that some inequality constraints become binding so that the UPFC cannot be fully utilized. Thus, it is important that the set of binding constraints does not change as x changes in order to obtain a reliable estimated IV. (However, we cannot actually verify this without many additional OPF solutions). For the 5bus case, Table 10.1 shows real power generations, transmission line losses and incremental values. Since the generation marginal cost, adjusted for marginal losses at bus 1 is lower than that at bus 2, it is pro¯table to obtain more 90 real power from generator 1 as long as its lossadjusted marginal cost stays lower and no operational limits are reached. The test results of the IEEE 14bus system are shown in Figure 10.3. Similar load conditions as in the 5bus system are used such that the constraint on current in line 1 is binding. As before, the lines with higher MVs produce higher estimated and actual IVs. An important fact to note is that higher voltage lines 1 through 7 are the most suitable locations to install the UPFC. This is because higher voltage lines have lower p.u. impedances and therefore, most of the power will be transferred through those lines and some of the lines may reach their maximum transfer capabilities. The simulation result for the IEEE 30bus system is shown in Figure 10.4. Again, lines with low marginal values tend to yield low estimated and actual incremental values, and the most valuable locations tend to be in the high voltage portions of the 0 5 10 15 20 25 30 35 40 0 5 10 15 20 Line number MV lT/15 lf/15 lr/5 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 Line number IV Estimated IV ($/hr) Actual IV ($/hr) Figure 10.4: Marginal ana incremental values for IEEE 30bus system. 91 system. Figure 10.4 also shows that estimated IVs are quite close to actual IVs. In addition, we see that locations with large incremental values are also those with large marginal values, thus supporting the idea of using only ¯rstorder sensitivity analysis to screen for promising locations for installation of a UPFC. 10.2 Result of GEOPF This section contains case studies that evaluate the capability of the GEOPF in maximizing social welfare associated with electricity generation cost, gas supply cost, and gas and electricity consumer's bene¯t. The coe±cients of quadratic generation cost and bene¯t curves are summarized in Appendix C. Since combined nodes are considered as nodes which simply convert gas energy to electric energy, these nodes are not associated with generation cost and gas consumer's bene¯t. Instead, we will Figure 10.5: Combined gas and electric network 92 use the coe±cients of the fuel rate curve for the gas¯red generator to calculate the amount of gas supplied to the combined node. 10.2.1 Test System I We will ¯rst consider a simple combined network with two real power generators, one of which is a gas¯red generator and the other is a coal¯red generator. Tieline connections in both natural gas and electricity networks are not considered in this test system. The system analyzed is a combined natural gas and electricity network, which consists of a 5bus system and a 15node gas system, as shown in Figure 10.5. The 15node system is a modi¯cation of an example system in [22], plus some made up bene¯t and cost functions. The 5bus electric system obtained by modifying the IEEE 14bus system has two real power generators at buses 1 and 2, one synchronous condenser at bus 3, and four loads at buses 2 through 5, and bus 1 serves as the reference/slack bus in the electric network. The generator at bus 2 is coincident with gas network node 15. We assume that node 1 is the knownpressure node in the gas system. The gas network has ¯fteen nodes consisting of ¯ve loads at nodes 3, 4, 13, 14, and 15, two sources at nodes 1 and 2, and eight compressor inletoutlet nodes 5 through 12. We assume that the compressors in the gas network are driven by gas turbines, and the gas is tapped from the outlet node of the compressor station. Appendix C shows the input data for gas and electricity network represented in Figure 10.5. We assume that generator 1 uses coal with a ¯xed price, and generator 2 uses natural gas with high price variations. To study the impact of the wellhead gas prices at source nodes 1 and 2 to the real power generation at combined node 15, three test conditions with di®erent wellhead gas prices are considered as follows: 93 Case 1: Wellhead gas price at node 1 = $2.07 /(106£ BTU) Wellhead gas price at node 2 = $2.12 /(106£ BTU) Case 2: Wellhead gas price at node 1 = $2.46 /(106£ BTU) Wellhead gas price at node 2 = $2.71 /(106£ BTU) Case 3: (pipeline branch between nodes 12 and 14 removed) Wellhead gas price at node 1 = $2.46 /(106£ BTU) Wellhead gas price at node 2 = $2.71 /(106£ BTU) Wellhead gas prices in cases 2 and 3 are higher than those in case 1. Table 10.2 shows simulation results for the three case studies. In case 1, the gas¯red generator reaches its maximum generation limit due to the low gas price at combined node 15. As noted in Table 10.2, the optimal real power generation PG2 at the combined node is quite sensitive to the wellhead gas prices, and reduced greatly as the wellhead gas prices increase. In case 3, we assume that the pipeline branch between nodes 12 and 14 is unavailable due to a forcedoutage. The gas¯red generator reaches its minimum generation limit, which indicates that contingencies occurring in the gas network may result in the signi¯cant change of electricity generation patterns. Now, let us compare the social welfare of nonintegrated gas and electricity operation with that obtained from the GEOPF for case 2. The social welfare from the GEOPF is $22,405.83/hr, and the gas price at node 15 is $2.82/(106£BTU). For nonintegrated network operation, let's say that there is a broker who purchases gas from node 15, and sells it to the gas¯red generator at bus 2. He
Click tabs to swap between content that is broken into logical sections.
Rating  
Title  Natural Gas and Electricity Optimal Power Flow 
Date  20040501 
Author  An, Seungwon 
Keywords  Nucleotide cofactorbinding, Cytochrome P450 reductase, NADH 
Department  Electrical Engineering 
Document Type  
Full Text Type  Open Access 
Abstract  The objective of this study is to propose an integrated analysis and solution methodology to the natural gas and electricity optimal power flow (GEOPF) problem for the operation of combined natural gas and electric power transmission networks from a central planning point of view. This paper also presents a screening technique for greatly reducing the computation involved in determining the optimal location of a unified power flow controller (UPFC) in a power system. The first and secondorder sensitivities of the generation cost with respect to UPFC control parameters are derived. A new model for an ideal UPFC is also developed. This work reviews fundamental modeling of the natural gas network to be used for the GEOPF, and describes the equality constraints which describe the energy transformation between gas and electric networks at combined nodes (i.e., generators). We also present the formulation of the natural gas loadflow problem, which includes the amount of gas consumed in compressor stations. We also show that the GEOPF returns the social welfare maximizing solution, which is difficult to obtain by nonintegrated operation of the two networks. Since electric energy and gas prices are so volatile, the GEOPF will significantly contribute to achieving optimal system operation for the combined gas and electric network. We have proposed first and secondorder sensitivity techniques to screen for optimal UPFC locations in a large power system by ignoring the transmission lines with low marginal values (MVs), and running a full optimal power flow (OPF) only for the lines with higher MVs to obtain actual incremental value (IV). This technique requires running only one OPF to obtain UPFC sensitivities for all possible transmission lines. To implement a sensitivitybased screening technique for optimally locating a UPFC in a power system, we propose a new UPFC model, which consists of an ideal transformer with a complex turn ratio and a variable shunt admittance. In this model, the UPFC control variables do not depend explicitly on UPFC input and output currents and voltages. Therefore, the problem size is not increased by the addition of the UPFC. 
Note  Dissertation 
Rights  © Oklahoma Agricultural and Mechanical Board of Regents 
Transcript  NATURAL GAS AND ELECTRICITY OPTIMAL POWER FLOW By SEUNGWON AN Bachelor of Engineering Korea Maritime University Pusan, Korea February, 1991 Master of Engineering Oklahoma State University Stillwater, Oklahoma May, 1999 Submitted to the Faculty of the Graduate College of the Oklahoma State University in partial ful¯llment of the requirements for the Degree of DOTOR OF PHILOSOPHY May, 2004 NATURAL GAS AND ELECTRICITY OPTIMAL POWER FLOW Thesis Approved: Thesis Advisor Dean of the Graduate College ii Acknowledgment I would like to express my sincere appreciation to my advisor Dr. Thomas W. Gedra for his intelligent supervision, friendship and con¯dence in me. I wish to thank the other members of my doctoral committee, Dr. Ramakumar, Dr. Yen and Dr. Misawa for their assistance throughout this work. I would like to give very special thanks to Dr. Ramakumar for his continuous guidance during my study at OSU. I wish to take this opportunity to express my gratitude to Mr. Lee Clark for his assistance and encouragement in all aspects of my work. I wish to express my sincere gratitude to the School of Electrical and Computer Engineering for providing me with this research opportunity and generous ¯nancial support. I am deeply indebted to all members of the Korean Catholic Community at Still water for their spiritual support. Especially, I would like to thank Mr. Hyunwoong Hong and Mrs. Eunjoo Kang for their brotherhood. I am deeply grateful to my wife, Misook Kim, for her patience and to our twin boys, Clemens and Martino for their existence. My ¯nal appreciation goes to my parents and brothers for their enduring love, support and encouragement. iii TABLE OF CONTENTS 1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 AC LOADFLOW . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 2.1 Flows on Transmission Systems . . . . . . . . . . . . . . . . . . . . . 6 2.2 Load°ow Problem Statement . . . . . . . . . . . . . . . . . . . . . . . 10 2.3 NewtonRaphson Method . . . . . . . . . . . . . . . . . . . . . . . . 13 3 ECONOMIC DISPACTH . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.1 Lossless Economic Dispatch . . . . . . . . . . . . . . . . . . . . . . . 18 3.2 Economic Dispatch with Transmission Losses . . . . . . . . . . . . . . 22 4 OPTIMAL POWER FLOW AND SOLUTION METHODS . . . . . . . . . 25 4.1 OPF by Newton's Method . . . . . . . . . . . . . . . . . . . . . . . . 27 4.2 PrimalDual InteriorPoint (PDIP) Method . . . . . . . . . . . . . . . 31 5 NATURAL GAS FLOW MODELING . . . . . . . . . . . . . . . . . . . . 39 5.1 Elements of Natural Gas Transmission Network . . . . . . . . . . . . 40 5.2 Network Topology . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 5.3 Matrix Representations of Network . . . . . . . . . . . . . . . . . . . 43 5.4 Flow Equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 5.5 Compressor Horsepower Equation . . . . . . . . . . . . . . . . . . . . 45 5.6 Conservation of Mass Flow . . . . . . . . . . . . . . . . . . . . . . . . 46 6 NATURAL GAS LOADFLOW . . . . . . . . . . . . . . . . . . . . . . . . 49 6.1 Load°ow Problem Statement . . . . . . . . . . . . . . . . . . . . . . . 50 6.2 Load°ow without Compressors . . . . . . . . . . . . . . . . . . . . . . 52 6.3 Load°ow with Compressors . . . . . . . . . . . . . . . . . . . . . . . 53 iv 7 UPFC MODELING FOR STEADYSTATE ANALYSIS . . . . . . . . . . 57 7.1 Operating Principles . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 7.2 Uncoupled Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 7.3 Ideal Transformer Model . . . . . . . . . . . . . . . . . . . . . . . . . 62 7.4 UPFC in a Transmission Line . . . . . . . . . . . . . . . . . . . . . . 65 8 GAS AND ELECTRICITY OPTIMAL POWER FLOW . . . . . . . . . . 69 8.1 Gas and Electric Combined Network . . . . . . . . . . . . . . . . . . 69 8.2 Gas and Electricity Optimal Power Flow . . . . . . . . . . . . . . . . 70 8.2.1 Cost, Bene¯t, and Social Welfare . . . . . . . . . . . . . . . . 71 8.2.2 Constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 9 OPTIMAL LOCATION OF A UPFC IN A POWER SYSTEM . . . . . . . 79 9.1 Optimal Power Flow with UPFC . . . . . . . . . . . . . . . . . . . . 79 9.2 FirstOrder Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . 81 9.3 SecondOrder Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . 84 9.4 Estimation of Incremental Value . . . . . . . . . . . . . . . . . . . . . 87 10 RESULTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 10.1 Result of Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . 88 10.2 Result of GEOPF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 10.2.1 Test System I . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 10.2.2 Test System II . . . . . . . . . . . . . . . . . . . . . . . . . . 95 11 CONCLUSIONS AND FUTURE WORK . . . . . . . . . . . . . . . . . . . 101 BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 A DERIVATIVES REQUIRED FOR GEOPF . . . . . . . . . . . . . . . . . 108 A.1 Electric Networks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108 A.2 Gas Network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115 A.3 Objective Function of GEOPF . . . . . . . . . . . . . . . . . . . . . . 119 B DERIVATIVES REQUIRED FOR UPFC SENSITIVITIES . . . . . . . . . 120 v C INPUT DATA FILE FORMAT . . . . . . . . . . . . . . . . . . . . . . . . 130 C.1 IEEE Common Data Format . . . . . . . . . . . . . . . . . . . . . . . 130 C.2 Generator Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 C.3 Electricity Load Data . . . . . . . . . . . . . . . . . . . . . . . . . . . 134 C.4 Natural Gas Common Data Format . . . . . . . . . . . . . . . . . . . 135 C.5 Natural Gas Load Data . . . . . . . . . . . . . . . . . . . . . . . . . . 137 C.6 Natural Gas Source Data . . . . . . . . . . . . . . . . . . . . . . . . . 138 D MATLAB CODE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 D.1 GEOPF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 D.2 Base Case OPF with Sensitivity Analysis . . . . . . . . . . . . . . . . 157 D.3 OPF with a UPFC . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167 vi LIST OF FIGURES 1.1 Monthly natural gas price trends at AECO hub in Canada . . . . . . 2 1.2 Average Virginia mine coal price (19752000, Virginia) . . . . . . . . 3 2.1 Bus i and connection to transmission system . . . . . . . . . . . . . . 7 4.1 E®ect of barrier term . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 4.2 Graphical representation of central path . . . . . . . . . . . . . . . . 37 5.1 Pipeline network representation . . . . . . . . . . . . . . . . . . . . . 40 5.2 Graph of the gas network . . . . . . . . . . . . . . . . . . . . . . . . . 41 7.1 General UPFC scheme . . . . . . . . . . . . . . . . . . . . . . . . . . 58 7.2 Proposed UPFC model in a transmission line . . . . . . . . . . . . . . 59 7.3 Phasor diagram of UPFC inputoutput voltages and currents . . . . . 60 7.4 Uncoupled UPFC model in a transmission line. . . . . . . . . . . . . 61 7.5 UPFC ideal transformer model . . . . . . . . . . . . . . . . . . . . . 62 7.6 Simpli¯ed UPFC circuit . . . . . . . . . . . . . . . . . . . . . . . . . 64 7.7 Cascaded transmission line with a UPFC . . . . . . . . . . . . . . . . 66 8.1 Combined natural gas and electricity network . . . . . . . . . . . . . 70 10.1 Diagram of 5bus subset of IEEE 14bus system . . . . . . . . . . . . 88 10.2 Marginal and incremental values for 5bus system . . . . . . . . . . . 89 10.3 Marginal and incremental values for IEEE 14bus system . . . . . . . 90 10.4 Marginal and incremental values for IEEE 30bus system . . . . . . . 91 vii 10.5 Combined gas and electric network . . . . . . . . . . . . . . . . . . . 92 10.6 SW losses due to nonintegrated operation for test system I . . . . . . 96 10.7 Combined gas and electric network . . . . . . . . . . . . . . . . . . . 97 10.8 SW losses due to nonintegrated operation for test system II . . . . . 98 viii LIST OF TABLES 10.1 Real power generation, line loss and incremental value . . . . . . . . . 89 10.2 GEOPF results for test system I . . . . . . . . . . . . . . . . . . . . . 99 10.3 GEOPF results for test system II . . . . . . . . . . . . . . . . . . . . 100 ix CHAPTER 1 INTRODUCTION Restructuring of natural gas and electric industries, both known as network industries, and deregulation of their products have been undertaken in the U.S. to minimize the social ine±ciency incurred by regulated monopolies when natural monopoly has disappeared from their technologies [1]. Due to the successful deregulation of the natural gas market, the wellhead price of natural gas declined by 44 % between 1983 and 1997. Electric power generation technologies utilizing natural gas are generally less capitalintensive and often more technically e±cient than other alternatives. Thus, natural gas is playing an important role in electric power generation since 1990 par tially due to incentives of the Clean Air Act Amendments (CAAA) of 1990. Since 1997, more than 120 GW of new capacity has been added by merchant energy com panies at a lower cost per installed MW and with shorter construction times than in the past [2]. Most of this new generation is fueled by natural gas. However, due to the lack of ¯nancial instruments to lock in future prices to reduce price risk, the large increase in natural gas prices, and excess reserve margins, more than 110 GW of installed capacity belongs to entities that have belowinvestmentgrade credit ratings (i.e., junk bonds) [2]. Thus, we have been experiencing large price volatility in major natural gas markets while price volatility in electric markets has declined steadily over the past three years [2]. In addition, the °oor price for natural gas in summer 1 Figure 1.1: Monthly natural gas price trends at AECO hub in Canada 2003 reached over U.S. $4.00 per MMBTU, as shown in Figure 1.1 (on the basis of 1 U.S $ = 0.75 nominal Canadian $ in December 2003). In the meantime, due to regulatory uncertainties posed by the restructuring of electric industries and environmental concerns, baseload units such as coal¯red power plants have been rarely built since 1990. The amount of electric energy generated by coal¯red generators has been almost constant since that time, resulting in decreases in coal price as illustrated in Figure 1.2. With increasing gas prices and with a diverging price gap between coal and natural gas, generation using natural gas may be no longer competitive with coal¯red generation. Consequently, a signi¯cant amount of planned natural gas generation capacity addition using natural gas has been cancelled or postponed in 2003 [2]. The exposure to price volatility and price risk should be better understood in competitive energy markets to promote e±cient use and expansion of electric trans mission and gas transportation networks, and to facilitate e®ective competition in power generation and gas supply [3]. The largescale construction of gas¯red elec 2 Figure 1.2: Average Virginia mine coal price (19752000, Virginia) tric generation, combined with electric power restructuring, is expected to create an even greater convergence between electricity and natural gas markets. Lack of fuel diversity (i.e., excessive dependence on natural gas with volatile prices, as opposed to coal, especially newer \clean coal" technology) will likely increase volatility of gas and electric prices, even if the price of coal is lower and less volatile. Therefore, viewing these markets separately without recognizing their increasing convergence is myopic at best, and disastrous at worst. Economic forces in energy markets have been driving not only the optimal oper ation of energy networks but also e±cient pricebased system planning. The uni¯ed power °ow controller (UPFC) is one of the most technically promising devices in the Flexible AC Transmission System (FACTS) family [4, 5, 6, 7]. It has the capability to control both voltage magnitude and phase angle, and can also independently provide (positive or negative) reactive power injections. However, it cannot be installed in all possible transmission lines due to its high capital cost. Thus, a need exists for developing a costbene¯t analysis technique to determine if installation of a UPFC 3 would be bene¯cial and the best location to install the UPFC. In principle, determin ing the optimal location for a UPFC is simple. For each possible location, we place a UPFC in the power system model and calculate the cost savings with respect to a base case (with no new UPFC installed). The operating cost at each time throughout the year and for each potential location is determined using an optimal power °ow (OPF) program. However, the computational burden of evaluating this annual value for every possible line is immense because an OPF problem must be solved for each possible UPFC location and at each of several time periods. Therefore, an e±cient screening technique is desired to identify the most promising locations so that at each point in time throughout the year, the exhaustive calculations described above do not have to be carried out for every location that is a candidate for installing the UPFC. In this research, we introduce fundamental natural gas modeling for steadystate analysis, and propose general formulations to solve a combined natural gas and elec tric optimal power °ow (GEOPF). In addition, a screening technique is developed for greatly reducing the computation involved in determining the optimal UPFC lo cation in a large power system. This technique requires running only one optimal power °ow (OPF) to obtain UPFC sensitivities for all possible transmission lines. To implement the screening technique, we develop a new mathematical model of UPFC under steadystate, consisting of an ideal transformer with a complex turns ratio and a variable shunt admittance. We begin by reviewing basic concepts which describe power °ows in power systems. Economic dispatch is discussed next, followed by the introduction of the optimal power °ow concept and its solution methodology. A GEOPF model is described in detail after the introduction of natural gas network modeling. A screening technique based on sensitivity analysis is discussed next and a new UPFC model is described. The GEOPF model has been successfully applied to two test cases: a 15node gas network with (i) a 5bus electric network and (ii) a 9bus electric network. The 4 screening technique (electricity only) has also been implemented in a 5bus system and IEEE 14 and 30bus systems. Discussions of results are presented in Chapter 10. 5 CHAPTER 2 AC LOADFLOW Electric load°ow problems are solved routinely to study power systems under both normal operating conditions and under various contingencies using predicted data or to analyze \what if" scenarios. It is also of great importance in power system planning for future expansion. The principal information obtained from load°ow analysis is the magnitude and phase angle of the voltage at each bus. Using these values, we can calculate real and reactive power °ows in transmission lines and transformers, and line losses. In this chapter, we describe the formulation of the bus admittance (¹Y bus) matrix, which represents a transmission line network, and explain how complex power injec tions into a network are related to the bus voltage magnitudes and angles. Then, we construct the load°ow problem and discuss some of its characteristics. 2.1 Flows on Transmission Systems Consider a network of transmission lines connecting a set of buses in a power trans mission system. We will ¯nd out the relationship between complex power and current injections into the transmission system and the phasor voltages at each bus. Figure 2.1 [8] shows connections to a transmission system at bus i. The doublesubscript notation ¹ Sik and ~Iik, for k 6= i, indicates complex power or phasor current (respec 6 Figure 2.1: Bus i and connection to transmission system [8] tively) °owing from bus i towards bus k along line ik. For i = k, the double subscript notations ¹ Sii and ~Iii indicate complex power and phasor current °owing to ground through the shunt element at bus i. The single subscript notations ¹ Si and ~Ii indicate total complex power and phasor current injection into the transmission system at bus i. Generally, the complex power injected at bus i will consist of power generated by any generator(s) at bus i, minus any power consumed at bus i, or ¹ Si = ¹ SGi ¡ ¹ SDi : (2.1) If a line exists between bus i and bus k, and has a series impedance of ¹zik, the series admittance of the line will be de¯ned as ¹yik = 1 ¹zik : (2.2) If bus i and bus k are not directly connected, then ¹yik = 0. In a similar way, we de¯ne the admittance of any shunt element connecting bus i to ground as ¹yii, which includes transmission line capacitive susceptance jb=2. The current °owing from bus i towards bus k is ~Iik = 1 ¹zik (~V i ¡ ~V k) = ¹yik(~V i ¡ ~V k); (2.3) while the current °owing through the shunt element yii from bus i to ground is ~Iii = ¹yii~V i; (2.4) 7 where ~V i denotes the phasor voltage at bus i. If the total number of buses is n, the total current injected into the power transmission system at bus i is ~Ii = ¹yii ~Vi + Xn k=1 k6=i ¹yik ³ ~V i ¡ ~V k ´ ; (2.5) or, collecting the ~V i terms, ~Ii = ~V i Xn k=1 ¹yik ¡ Xn k=1 k6=i ¹yik ~V k: (2.6) The right side of the above equation motivates the following de¯nition: ¹ Yii = Xn k=1 ¹yik; (2.7) ¹ Yik = ¡¹yik; i 6= k: (2.8) With this de¯nition, we can write ~Ii = Xn k=1 ¹ Yik ~V k: (2.9) In matrix form, the de¯nitions in equations (2.7) and (2.8) can be written as ¹Y bus = 2 66666664 ¹y11 + ¢ ¢ ¢ ¹y1n ¡¹y12 ¢ ¢ ¢ ¡¹y1n ¡¹y21 ¹y21 + ¢ ¢ ¢ ¹y2n ¢ ¢ ¢ ¡¹y2n ... . . . ... ¡¹yn1 ¡¹yn2 ¢ ¢ ¢ ¹yn1 + ¢ ¢ ¢ ¹ynn 3 77777775 : In words, the diagonal element ii of the ¹Y bus matrix consists of the sum of all admit tances connected to bus i (whether they are series admittances connecting to another bus or are shunt admittances), while the o®diagonal element ik is the negative of the series admittance connecting bus i to bus k. In the case that there is no transmission line connecting bus i with bus k, ¹yik = 0 as mentioned above. Likewise, ¹yii = 0 if no shunt element is present between bus i and ground. Note that, since the series admittance ¹yik connecting bus i to bus k is the same as that connecting bus k to bus i, 8 we have ¹yik = ¹yki, so that the matrix ¹Y bus is symmetrical, assuming no phaseshifting elements, which we have implicitly done. If we de¯ne the bus voltage and bus current injection vectors as ~V = 2 66666664 ~V 1 ~V 2 ... ~V n 3 77777775 and ~I = 2 66666664 ~I1 ~I2 ... ~In 3 77777775 ; then we can write equation (2.9) in matrixvector form as ~I = ¹Y bus~V : (2.10) The matrix ¹Y bus is called the bus admittance matrix, since it relates the bus voltages and the bus current injections. The complex power injection at bus i can be written in terms of ¹Y bus matrix elements as ¹ Si = ~V i~I¤ i = ~V i Xn k=1 ¹ Y ¤ ik ~V¤ k ; (2.11) = Vi Xn k=1 VkYike(µi¡µk¡±ik); where ~V i = Vi\µi; ~V k = Vk\µk; and ¹ Yik = Gik + jBik = Yik\±ik: We can split equation (2.12) into real and reactive parts, in terms of the bus voltage magnitudes fVi; i = 1; : : : ; ng and bus voltage angles fµi; i = 1; : : : ; ng. Then the real and imaginary parts of equation (2.12) become Pi(~V ) = Xn k=1 ViVkYik cos(µi ¡ µk ¡ ±ik); (2.12) Qi(~V ) = Xn k=1 ViVkYik sin(µi ¡ µk ¡ ±ik): (2.13) 9 The real and reactive power injections at bus i need to satisfy the following conditions: PGi ¡ PLi = Pi(~V ); i = 1; : : : ; n; (2.14) QGi ¡ QLi = Qi(~V ); i = 1; : : : ; n; (2.15) where PGi = real power generation at bus i; PLi = real power load at bus i; QGi = reactive power generation at bus i; QLi = reactive power load at bus i: We can rewrite equations (2.12), (2.12) and (2.13) in matrix form as ¹S = 2 66666664 ¹ S1 ¹ S2 ... ¹ Sn 3 77777775 = 2 66666664 ~V 1~I¤ 1 ~V 2~I¤ 2 ... ~V n~I¤ n 3 77777775 = diagf~V g~I¤ = diagf~I¤g~V ; (2.16) P = real(¹S ); (2.17) Q = imag(¹S ); (2.18) where diagf~V g denotes the diagonal matrix with the elements of the vector ~V on its diagonal. 2.2 Load°ow Problem Statement Since the load°ow solution is a prerequisite for many other analytical studies, it is necessary to mathematically state the problem to prepare the ground for the estab lishment of its links to economic dispatch and optimal power °ow. The load°ow problem is stated below [8, 9, 10]: 10 ² Given a power system described by a ¹Y bus matrix, and given a subset of bus voltage magnitudes, bus voltage angles, and real and reactive power bus injec tions, ² Determine the other voltage magnitudes and angles and real and reactive power injections. More precisely, two of the four quantities Vi; µi; Pi; and Qi at each bus are speci¯ed, and the other two are to be determined. Each bus is classi¯ed based on the two known quantities, as follows: ² PV bus, or generator bus, or voltagecontrolled bus. For a bus i of this type, we assume that we know the real power injection Pi and the voltage magnitude Vi. This is because the real power generation can be controlled by adjusting the prime mover input, and the voltage magnitude can be controlled by adjusting the generator ¯eld current. However, certain buses without generators may have some means of voltage support at that bus (such as a synchronous condenser, capacitor banks, or a static VAr compensator). ² PQ bus, or load bus. For a nongenerator bus i, we assume that we know the real and reactive power injections Pi and Qi, and the bus voltage magnitude Vi and angle µi are to be determined. In practice, only real power load is known and the reactive power load is calculated based on an assumed power factor by Qi = QGi ¡ PLi p:f:i q 1 ¡ p:f:2 i : We assume that a load is characterized by its constant complex power demand (which does not depend on the bus voltage Vi, as opposed to assuming that the load is a constant current or constant impedance load in which case the complex power demand would depend on the bus voltage Vi.) 11 In fact, solving the load°ow problem with buses of only these two types is not in general possible. The ¯rst reason is that, in the load°ow equations, the bus voltage angles never appear by themselves, but instead appear only as angle di®erences of the form µi ¡ µk. Therefore, adding an angle to every bus voltage angle, will not change the values of real or reactive power injections in equations (2.12) and (2.13). Since phasor voltages are always expressed with respect to some reference voltage, we must decide on some reference phasor and refer all other phase angles to that reference. In fact, therefore, there are only n¡1 angles which in°uence the load°ows. We therefore pick one bus, say bus 1, to serve as our phasor reference, and set µ1 = 0o. Another reason is that solving a load°ow for a system containing only PV and PQ buses would imply that we know the real power injections at every single bus. In fact, we cannot specify all n real power injections, since we do not know all n real power injections until we solve the load°ow. Mathematically, the load°ow is overdetermined if we suppose we know all n real power injections. Specifying the injections at all buses is the same as specifying the real losses of the power system, which we cannot know until the load°ow is solved. Instead, we must pick one bus, and allow the real power injection at that bus to be whatever value is required to solve the load°ow equations. Thus, in addition to the PV and PQ bus types, we have a third bus type: ² Slack Bus, or V µ, or reference, or swing bus. Typically it is a generator bus, and the voltage angle of the slack bus serves as a reference for the angles of all other bus voltages. We assume that we know V1 and µ1, but we do not know P1 and Q1. The usual practice is to set µ1 = 0o and V1 = 1 although other values of V1 are ¯ne too. 12 2.3 NewtonRaphson Method Equations (2.14) and (2.15) are in the form of the vector nonlinear equation y = f (x), which can be solved by NewtonRaphson method [8, 9, 10]. The NewtonRaphson method is based on Taylor series expansion of f (x) about an operating point xo. y = f (xo) + @f @x ¯¯¯¯ x=xo (x ¡ xo) + H:O:T: (2.19) Neglecting the higher order terms in equation (2.19) and solving for x, we have x = xo + · @f @x ¯¯¯¯ x=xo ¸¡1 ¢ (y ¡ f (xo)) : (2.20) The NewtonRaphson method replaces xo by the old value xk and x by the new value xk+1 for the iterative solution as shown below xk+1 = xk + J¡1 ¢ ¢f ; (2.21) where J = @f @x ¯¯¯¯ x=xk ; ¢f = ³ y ¡ f (xk) ´ : Equation (2.21) is repeated until the mismatches ¢f are less than a speci¯ed tolerance, or the algorithm diverges. Now, let us construct a mathematical formulation to solve the AC load°ow prob lem using the NewtonRaphson method. Assume that there are g generators. We assume that bus 1 is the slack bus, and buses 2 through g are PV buses. The n ¡ g buses will be assumed to be PQ buses. Then the problem becomes: Given : µ1 P2 ¢ ¢ ¢ Pg Pg+1 ¢ ¢ ¢ Pn V1 V2 ¢ ¢ ¢ Vg Qg+1 ¢ ¢ ¢ Qn Determine : P1 µ2 ¢ ¢ ¢ µg µg+1 ¢ ¢ ¢ µn Q1 Q2 ¢ ¢ ¢ Qg Vg+1 ¢ ¢ ¢ Vn 13 Note that ¯nding P1 and Q1 through Qg is trivial, once all the voltage magnitudes and angles are known. The di±cult part is to ¯nd n¡1 unknown angles and the n¡g unknown voltage magnitudes. Let us de¯ne the x, y, and f vectors for the load°ow problem as x = 2 64 µ V 3 75 = 2 666666666666664 µ2 ... µn Vg+1 ... Vn 3 777777777777775 ; y = 2 666666666666664 P2 ... Pn Qg+1 ... Qn 3 777777777777775 ; f (x) = 2 666666666666664 P2(x) ... Pn(x) Qg+1(x) ... Qn(x) 3 777777777777775 : Note that there are 2n¡g¡1 nonlinear equations, and there are 2n¡1¡g unknown voltages and angles. So we have exactly the right number of variables to force all the mismatches to zero. We will guess the unknown voltage magnitudes and angles of xk, and compare the calculated values of P(xk) and Q(xk) to the known values of P and Q. Then, we can evaluate the mismatches by ¢f = 2 666666666666664 ¢P2 ... ¢Pn ¢Qg+1 ... ¢Qn 3 777777777777775 = 2 666666666666664 P2 ¡ P2(xk) ... Pn ¡ Pn(xk) Qg+1 ¡ Qg+1(xk) ... Qn ¡ Qn(xk) 3 777777777777775 : (2.22) In this set of equations, the functions Pi(xk) and Qi(xk) are those obtained from equations (2.12) and (2.13), or (2.17) and (2.18), while the numbers Pi and Qi are the known values of real and reactive power injections at the buses where the injections are known. The resulting quantities ¢Pi and ¢Qi are known as the real and reactive mismatch terms, because they represent the di®erence between the known injections and the values calculated based on our guesses. 14 To solve the AC load°ow problem, we need to update the Jacobian at each iter ation. Many authors [9, 10] derived the Jacobian directly from equations (2.12) and (2.13). Even though this method is straightforward, it requires at least two for loop routines in Matlab. To avoid the for loop command, and to improve the speed of each iteration, we present an e±cient technique to construct the Jacobian. Let us de¯ne JSV by JSV = 2 6666666664 @ ¹ S1 @V1 @ ¹ S1 @V2 ¢ ¢ ¢ @ ¹ S1 @Vn @ ¹ S2 @V1 @ ¹ S2 @V2 ¢ ¢ ¢ @ ¹ S2 @Vn ... ... . . . ... @ ¹ Sn @V1 @ ¹ Sn @V2 ¢ ¢ ¢ @ ¹ Sn @Vn 3 7777777775 : The diagonal and o®diagonal elements of JSV are given by @ ¹ Si @Vi = ejµi " Xn k=1 ¹ Yik ~V k #¤ + ~V i ¹ Y ¤ ii e¡jµi ; (2.23) @ ¹ Si @Vk = ~V i ¹ Y ¤ ike¡jµk ; i 6= k (2.24) By noticing similarities in the above two equations, we can rewrite JSV in compact form to be used in Matlab by JSV = diag(~V ) ¤ conj(¹Y bus) ¤ conj(diag(expjµ)) + diag(expjµ :  ¤{zconj(~I}) element to element multiplication ): (2.25) Now, let us de¯ne JSµ by Jsµ = 2 6666666664 @ ¹ S1 @µ1 @ ¹ S1 @µ2 ¢ ¢ ¢ @ ¹ S1 @µn @ ¹ S2 @µ1 @ ¹ S2 @µ2 ¢ ¢ ¢ @ ¹ S2 @µn ... ... . . . ... @ ¹ Sn @µ1 @ ¹ Sn @µ2 ¢ ¢ ¢ @ ¹ Sn @µn 3 7777777775 : The diagonal and o®diagonal elements of JSµ are given by @ ¹ Si @µi = j ¹ Si ¡ j~V i ¹ Y ¤ ii ~V ¤ i ; (2.26) @ ¹ Si @µk = ¡j~V i ¹ Y ¤ ik ~V ¤ k ; i 6= k: (2.27) 15 In a similar way, JSµ can be expressed in compact form by JSµ = ¡jdiag(~V ) ¤ conj(¹Y bus) ¤ conj(diag(~V )) + diag(j¹S ): (2.28) We can split the Jacobian into real and reactive parts as follows: PSµ = real(JSµ ); PSV = real(JSV ); QSµ = imag(JSµ ); QSV = imag(JSV ): Then, the full Jacobian J is constructed as J = 2 64 PSµ PSV QSµ QSV 3 75 : Since the voltage magnitudes at PV buses and the angle at the slack bus are known, their corresponding columns will be truncated from the Jacobian J. In addition, since the real power injection at the slack bus, and the reactive power injections at PV buses and at the slack bus are unknown, their corresponding rows are removed from the Jacobian J. Then, the modi¯ed Jacobian Jm at the kth iteration becomes Jm = 2 666666666666666664 @P2(xk) @µ2 ¢ ¢ ¢ @P2(xk) @µn @P2(xk) @Vg+1 ¢ ¢ ¢ @P2(xk) @Vn ... . . . ... ... . . . ... @Pn(xk) @µ2 ¢ ¢ ¢ @Pn(xk) @µn @Pn(xk) @Vg+1 ¢ ¢ ¢ @Pn(xk) @Vn @Qg+1(xk) @µ2 ¢ ¢ ¢ @Qg+1(xk) @µn @Qg+1(xk) @Vg+1 ¢ ¢ ¢ @Qg+1(xk) @Vn .... . . ... ... . . . ... @Qn(xk) @µ2 ¢ ¢ ¢ @Qn(xk) @µn @Qn(xk) @Vg+1 ¢ ¢ ¢ @Qn(xk) @Vn 3 777777777777777775 : By using the mismatches in equation (2.22) and the modi¯ed Jacobian Jm, we will update the unknown variables x by equation (2.21). Finally, once we have forced all the real and reactive mismatches to reasonably small values, by ¯nding the unknown values of voltage magnitude and angle, we can now go back and ¯nd the unknown real and reactive injections. 16 CHAPTER 3 ECONOMIC DISPACTH The economic dispatch problem is de¯ned as the process of providing the required real power load demand and line losses by allocating generation among a set of online generating units such that total generation cost is minimized [8, 11, 12]. Let C be the generating cost in $/hr, and it is often modelled analytically as a quadratic function of the power generated [8]. The cost function is derived from the generator heat rate curve. Analytically, the heat rate curve with a unit of BTU/kWh is represented in the following form: H(PG) = a PG + b + cPG; (3.1) where PG is the unit's real power generation level measured in kW. Given the heat rate curve, the next important function is the fuel rate, which is simply the rate, in BTU/hr, of consumption of fuel energy. So F(PG) = PGH(PG) = a + bPG + cP2G : (3.2) If the cost of fuel is known as k $/BTU, then the cost function with a unit of $/hr becomes C(PG) = kF(PG) = ka + kbPG + kcP2G = ® + ¯PG + °P2G : (3.3) In classic economic dispatch problems, the essential constraint on the operation is that the sum of the output powers must equal the load demand plus losses. In 17 addition, there are two inequality constraints that must be satis¯ed for each of the generator units. That is, the power output of each unit must be greater than or equal to the minimum power permitted and must also be less than or equal to the maximum power permitted on that particular unit. 3.1 Lossless Economic Dispatch We assume that the utility is responsible for supplying its customers' load. The utility's objective is to minimize the total cost of generation, assuming that trans mission losses are neglected. Then, the ideal economic dispatch problem is stated in terms of minimizing total generation cost subject to satisfying the total load demand. Mathematical formulation of lossless economic dispatch problem can be expressed as min fPG1 ;¢¢¢ ;PGng Xn i=1 Ci(PGi) (3.4) subject to: Xn i=1 PGi = PD (3.5) and Pmin Gi · PGi · Pmax Gi ; i = 1; : : : ; n (3.6) where PD = Xn i=1 PDi is the total system load. We will ignore the constraints in equation (3.6) on generator limits, assuming that these limits are not binding. Considering only constraint (3.5), the Lagrangian is L = Xn i=1 Ci(PGi) ¡ ¸ " Xn i=1 PGi ¡ PD # : (3.7) Di®erentiating with respect to PGi , we obtain the ¯rstorder conditions 0 = @L @PGi = MCi ¡ ¸; i = 1; : : : ; n (3.8) 18 or MCi = ¸; i = 1; : : : ; n where MCi = @Ci @PGi : The quantity ¸, the Lagrange multiplier, associated with the energy balance con straint in equation (3.5), is universally called the \system lambda" and is the price associated with generating slightly more energy. Thus, the criterion for optimal eco nomic distribution of load among n generators is that all generators should generate output at the same marginal operating cost, sometimes called incremental cost. That is @C1 @PG1 = @C2 @PG2 = ¢ ¢ ¢ = @Cn @PGn : (3.9) In the case of quadratic cost functions, the solution to the economic dispatch problem can be calculated analytically, as follows [8]. We have, for each unit, ¸ = MCi = ¯i + 2°iPGi ; (3.10) which can be solved for PGi to give PGi = ¸ ¡ ¯i 2°i : (3.11) The total generation at this ¸ is obtained by summing over all the units. Setting this total equal to the load PD, we obtain PD = Xn i=1 PGi = Xn i=1 ¸ ¡ ¯i 2°i = ¸ Xn i=1 1 2°i ¡ Xn i=1 ¯i 2°i ; (3.12) which can be solved to obtain the ¸ required for a given PD: ¸ = · PD + Pn i=1 ¯i 2°i ¸ Pn i=1 1 2°i (3.13) 19 Finally, this value of ¸ can be substituted back into the expression for PGk to obtain PGk = 1 2°k Pn i=1 1 2°i " PD + Xn i=1 ¯i 2°i # ¡ ¯k 2°k ; (3.14) which gives the dispatch for unit k directly in terms of the total load PD. We de¯ne the \participation factor" for unit k as Kk = 1 2°k Pn i=1 1 2°i : We can see from the de¯nition that Xn i=1 Ki = 1: Generally, the meaning of the participation factor is Kk = dPGk dPD : Thus, if the load were to increase by a small increment (say 1MW), unit k would supply the fraction Kk of the increase. The fact that the participation factors sum to 1 simply means that any increase in load is met exactly by an increase in generation. We can include the generator limit constraints (3.6) in one of two ways. Tradi tionally, we ignore the limits, perform economic dispatch, and check to see if any limits are violated. If the limit for one generator is violated, we set the generator to its limit, remove it from the set of generators included in economic dispatch, and subtract the limit from the total load PD. In other words, we take the generator \o® of economic dispatch" and then proceed to treat it as a negative load1. More formally, we can just include the generation limits in the Lagrangian func tion. Let ¹min i be the Lagrange multiplier for the lower limit of unit i, and ¹max i be 1This may result in \cycling" of the set of active constraints, especially if there are many diverse generators. Other optimization techniques, such as a primaldual interiorpoint (PDIP) method [13] can be used to avoid this problem. 20 the multiplier for the generator's upper limit. Then, the Lagrangian becomes L = Xn i=1 Ci(PGi) ¡ ¸ " Xn i=1 PGi ¡ PD # + Xn i=1 ¹min i £ Pmin Gi ¡ PGi ¤ + Xn i=1 ¹max i £ PGi ¡ Pmax Gi ¤ : (3.15) Now, taking the derivative of L with respect to each PGi , and setting to zero, we obtain 0 = @L @PGi = MCi ¡ ¸ ¡ ¹min i + ¹max i ; i = 1; : : : ; n (3.16) or MCi = ¸ + ¹min i ¡ ¹max i ; i = 1; : : : ; n (3.17) We also have the complementary slackness conditions on the inequality constraints: ¹min i £ Pmin Gi ¡ PGi ¤ = 0; and ¹max i £ PGi ¡ Pmax Gi ¤ = 0: But PGi cannot be at its upper and lower limits simultaneously, so there are only three possibilities for each generator: PGi = Pmax Gi ) ¹max i ¸ 0; ¹min i = 0; Pmin Gi · PGi · Pmax Gi ) ¹max i = 0; ¹min i = 0; PGi = Pmin Gi ) ¹max i = 0; ¹min i ¸ 0: If the limits on some generators are binding, the generators may not be able to generate output at the system incremental cost ¸. Suppose that none of generators reach their limits, and all generators participate in economic dispatch. As the load increases by a small increment, each generator will supply the fraction of the load increase, which is determined by the participation factor. If a generator hits its maximum limit, the generator cannot participate in economic dispatch even though it has a lower incremental cost than the system incremental cost ¸. 21 On the other hand, suppose that the load decreases and one of generators hits its minimum limit. Since the generator must generate the required minimum capacity in order to stay online, the generator cannot participate in economic dispatch for the load decrements. Thus, if the load decreases further, the incremental cost of the generator hit its minimum limit is greater than or equal to the system incremental cost ¸. The optimal criterion for lossless economic dispatch with generator limits can be summarized as [11] PGi = Pmax Gi ) dCi dPGi · ¸; Pmin Gi · PGi · Pmax Gi ) dCi dPGi = ¸; PGi = Pmin Gi ) dCi dPGi ¸ ¸: 3.2 Economic Dispatch with Transmission Losses We now wish to include the e®ect of transmission losses on economic dispatch [8, 12]. In lossless dispatch, the location of individual loads did not matter, but when transmission losses are considered, the solution depends both on the load locations and on the outputs of individual generators around the system. If we know the load PDi , and generation PGi at each bus in the system, we can calculate the real power injections P2; : : : ; Pn using Pi = PGi ¡ PDi . From load°ow, we can then calculate the power injection at the slack bus, P1 = PG1 ¡ PD1 . From this information, we can calculate losses as PL(PG2 ; ¢ ¢ ¢ ; PGn; PD2 ; ¢ ¢ ¢ ; PDn) = Xn i=2 (PGi ¡ PDi) + P1(PG2 ; ¢ ¢ ¢ ; PGn; PD2 ; ¢ ¢ ¢ ; PDn): (3.18) Note that in our formulation, we do not consider PL to be a function of PG1 , since the slack bus generation PG1 is completely determined by the injections at other buses and the slack bus demand PD1 . 22 Now our problem becomes min fPG1 ;¢¢¢ ;PGng Xn i=1 Ci(PGi) (3.19) subject to: Xn i=1 PGi = PD + PL (PG2 ; ¢ ¢ ¢ ; PGn) (3.20) and Pmin Gi · PGi · Pmax Gi ; i = 1; : : : ; n (3.21) where we have suppressed the dependence of PL on the loads fPDig, which are as sumed to be ¯xed. The Lagrangian for the problem now becomes L = Xn i=1 Ci(PGi) ¡ ¸ " Xn i=1 PGi ¡ PD ¡ PL(PG2 ; ¢ ¢ ¢ ; PGn) # : (3.22) Di®erentiating and setting to zero yields 0 = @L @PG1 = MC1 ¡ ¸; and 0 = @L @PGi = MCi ¡ ¸ µ 1 ¡ @PL @PGi ¶ ; i = 2; : : : ; n or ¸ = MC1; and ¸ = MCi 1 1 ¡ @PL @PGi ; i = 2; : : : ; n Finally, let us de¯ne the losspenalty factors L1 = 1; and Li = 1 1 ¡ @PL @PGi ; i = 2; : : : ; n (3.23) Then we can write our condition for economic dispatch as ¸ = MCi ¢ Li; i = 1; : : : ; n (3.24) The idea behind the loss penalty factors can be explained as follows. If increasing a generator's output increases system losses, then that generator'entire increment of 23 generation is not available to the system. So, if a generator's output increases by 1MW but losses increase by 0.1MW, the generator has only provided a net increase in system generation of 0.9MW. Thus, the cost of this increment is not the generator's marginal cost but MC=0:9 instead, or MCi ¢Li. Similarly, an increase in a generator's output could result in a decrease in system losses. Such a generator would have a loss penalty factor Li < 1, since the e®ective cost of incremental generation from such a generator would be less than its actual marginal cost. The loss penalty factors can either be obtained as the result of running a load°ow, or they can be obtained from a quadratic approximation to the loss function, called Bmatrix formula [11]. 24 CHAPTER 4 OPTIMAL POWER FLOW AND SOLUTION METHODS Note that the classic economic dispatch (ED) problem does not strictly take into account power °ows in the transmission system. Typically, one would solve the eco nomic dispatch problem, then the load°ow problem, and repeat the process. The most recent load°ow provides the losspenalty factors for the current operating state, then the ED changes the state slightly to minimize cost, and so on [8]. Provided no transmission lines are overloaded, and no voltages are outside of the allowable range, this works well. However, if the result of the economic dispatch is fed into the load °ow, and the result indicates that transmission line loading is excessive, or that bus voltage magnitudes are outside of the allowable range (typically 95% to 105% of the nominal value), no information is given by the load°ow to indicate how to redispatch generation to alleviate the overloads or restore acceptable voltage levels. Often, the system dispatcher, being extremely familiar with the system, will take some units o® of economic dispatch and manually redispatch them to remove the line overload. There is no guarantee that such a procedure will result in the minimum cost subject to operating limits although sometimes they do fairly well. Similarly, economic dispatch does not consider reactive °ows or bus voltages, whereas load°ow considers these quantities as given or to be found from other voltages and reactive injections [8]. Neither problem considers the adjustment of bus voltages to help to minimize cost, even though reactive power °ows can contribute signi¯cantly 25 to real power losses and thus to overall cost. However, we can reformulate our optimization problem by including economic dispatch, voltage and reactive injections as decision variables, and various operational limits. The optimal power °ow (OPF) combines economic dispatch and load°ow into a single problem, which can be written in very general terms as min Y C(Y ) (4.1) subject to: hi(Y ) = 0; i = 1; ¢ ¢ ¢ ; n; gi(Y ) · 0; i = 1; ¢ ¢ ¢ ; m; where ² Y : a vector of control and state variables, ² C(Y ): an objective function, ² h(Y ): a set of equality constraints, which are power °ow equations, ² g(Y ): a set of inequality constraints, such as voltage limits, generator capacity limits, and line °ow limits. Tapchanging transformers and/or Flexible AC Transmission System (FACTS) de vices also can be incorporated in the OPF problem [9, 14]. The OPF has been used as a tool to improve power system planning, and oper ation by adjusting the objective function and/or the constraints. From the power system planning point of view, the OPF can be used to determine the optimal types, sizes, settings, capital costs, and optimal locations of resources, such as generators, transmission lines, and FACTS devices [14, 15]. Another application of the OPF is to determine various system marginal costs. It can be used to determine shortrun electric pricing (i.e. spot pricing), transmission 26 line pricing, and pricing ancillary services such as voltage support through MVAR support. Since the work in this research is based on solving the proposed optimization problem by a primaldual interiorpoint (PDIP) method using logarithmic barrier function, the PDIP will be discussed in detail after a brief introduction of Newton's method. 4.1 OPF by Newton's Method NewtonRaphson method has been the standard solution algorithm for the economic dispatch and load°ow problems for several decades [11, 12]. Newton's method is a very powerful algorithm because of its rapid convergence near the solution. This property is especially bene¯cial for power system applications because an initial guess close to the solution is easily obtained [16]. For example, voltage magnitude at each bus is presumably near the rated system value, generator outputs can be estimated from historical data, and transformer tap ratios are near 1.0 p.u. during steadystate operation. The solution of the constrained optimization problem stated in (4.1) requires the mathematical formation of the Lagrangian by L(Y; ¸; ¹) = C(Y ) + Xn i=1 ¸ihi(Y ) + X i2A ¹igi(Y ); (4.2) where ¸i is the Lagrange multiplier for the ith equality constraint. Assuming that we know which inequality constraints are binding, and have put them in the set A, then the inequality constraints can now be enforced as equality constraints. Thus the ¹0 is in equation (4.2) have the same property as ¸0 is and they are the Lagrange multipliers for binding inequality constraints. However, we need ¹i ¸ 0 for every i [8]. We can ignore the inequality constraints that are not binding since their ¹0s are 27 known to be zero by complementary slackness condition [8]. That is gi(Y ) · 0 ) ¹i = 0; gi(Y ) = 0 ) ¹i ¸ 0: Therefore, only binding inequality constraints are included in the Lagrangian function (4.2) with corresponding nonzero ¹0s. Solution of a constrained optimization problem can be solved by adjusting control and state variables, and Lagrange multipliers to satisfy the following firstorder necessary optimality conditions : 1) @L @Yi = 0; (4.3) 2) @L @¸i = 0; (4.4) 3) @L @¹i = gi = 0; (4.5) 4) ¹i ¸ 0 and ¹igi = 0: (4.6) The above equations are also called the KarushKuhnTucker (KKT) conditions. Let us de¯ne !(z) = rzL(z) = · @L @Y @L @¸ @L @¹A ¸T = 0; (4.7) where z is a vector of [ Y T ¸T ¹T A ]T , and A represents the binding inequality constraints. To solve the KKT conditions, Newton's method is applied by using the Taylor's series expansion around a current point zp as: !(z) = !(zP ) + @!(z) @z ¯¯¯¯ z=zp ¢ (z ¡ zp) + 1 2 (z ¡ zp)T ¢ @2!(z) @z2 ¯¯¯¯ z=zp ¢ (z ¡ zp) + ¢ ¢ ¢  {z } H.O.T = 0; The current point zp can either be an initial guess in the ¯rst iteration of the computation, or the estimate solution from the prior iteration. Recall that we want !(z) = 0. By ignoring the high order terms (H.O.T) and de¯ning ¢z = z ¡ zp, the 28 above equation can be rewritten as: @!(z) @z ¯¯¯¯ z=zp ¢ ¢z = ¡!(zp): (4.8) The quantity ¢z is the update vector, or the Newton step, and it tells how far and in which direction the variables and multipliers should move from this current point zp to get closer to the solution. Since !(z) is the gradient of the Lagrangian function L(z), equation (4.8) can be written in terms of the Lagrangian function L(z) as: @2L(z) @z2 ¯¯¯¯ z=zp ¢ ¢z = ¡ @L(z) @z ¯¯¯¯ z=zp : (4.9) Or, simply W ¢ ¢z = ¡!(z); (4.10) where W denotes the second order derivatives (or the Hessian matrix) and ! is the gradient, both of the Lagrangian function with respect to z evaluated at the current point. Equation (4.10) can be written in matrix form as 2 6666664 HY JT AT J 0 0 A 0 0 3 7777775 2 6666664 ¢Y ¢¸ ¢¹A 3 7777775 = ¡ 2 666664 rY L r¸L r¹AL 3 777775 ; (4.11) where the Hessian matrix HY and Jacobian matrices J and A are given as follows: HY = @2L(z) @Y 2 ; J = @2L(z) @¸@Y = @h(Y ) @Y ; A = @2L(z) @¹A@Y = @gA(Y ) @Y ; where gA is a set of binding inequality constraints. The Newton step can be obtained by solving (4.10). Then a vector of estimated solution for the next iteration is updated as: zp+1 = zp + ®¢z; (4.12) 29 where ® is usually 1, but can be adjusted to values above or below 1 to speed up convergence or cause convergence in a divergent case. It is important that special attention be paid to the inequality constraints. Equation (4.2) only includes binding inequality constraints being enforced as equality constraints. Thus, after obtaining an updated set of variables and multipliers, a new set of binding inequality constraints (or what we think is the active set) should be determined as follows [12]: ² If the updated ¹0s of the constraint functions in the current active set are zero or have become negative, then the corresponding constraints must be released from the current active set because ¹i < 0 implies that gi = 0 keeps the trial solution at the edge of the feasible region instead of allowing the trial solution to move into the interior of the feasible region. ² If other constraint functions evaluated at the updated variables violate their limits, then those constraints must be included in the new active set. The variable ® may be chosen to prevent constraint violations, but ® < 1 to avoid infeasibility implies that the constraint would otherwise be violated. As a result, if ¹i is positive, continued enforcement will result in an improvement of the objective function, and enforcement is maintained. If ¹i is negative, then enforcement will result in an decrease of the objective function, and enforcement is stopped. Once the active set has been updated, !(zp+1) is checked for convergence. There are several criteria for checking convergence of Newton's method. The convergent tolerance may be set on the maximum absolute value of elements in !(z), or on its norm. If the updated zp+1 does not satisfy the desired convergence criterion, the Newton step calculation is repeated. 30 4.2 PrimalDual InteriorPoint (PDIP) Method One of disadvantages of Newton's method is to identify a set of binding inequality constraints, or active constraints. Among several methods to avoid the di±culty asso ciated with guessing the correct active set, the PDIP method has been acknowledged as one of the most successful [12, 17, 18]. Interior point methods for optimization have been widely known since the publication of Karmarkar's seminal paper in 1984 [19]. Barrier function methods were proposed much earlier in Russia but little attention was paid because the algorithm was so slow in implementation. Later, this method was shown to be equivalent to the interior point methods. Karmarkar's method re sults in numerical illconditioning although this problem is not so bad with the PDIP method. The method uses a barrier function that is continuous in the interior of the feasible set, and becomes unbounded as the boundary of the set is approached from its interior. Two examples of such a function [13, 20] are the logarithmic function, as shown in Figure 4.1 Á(Y ) = ¡ Xm i=1 ln(¡gi(Y )); (4.13) and the inverse of the inequality function Á0(Y ) = Xm i=1 1 ¡gi(Y ) ; (4.14) where gi(Y ) · 0: The barrier method generates a sequence of strictly feasible iterates that converge to a solution of the problem from the interior of the feasible region [13, 20]. To apply the primaldual interior point algorithm to the OPF problem that has equality and inequality constraints, we construct the nonlinear equality and 31 inequalityconstrained optimization problem as min Y f(Y ) (4.15) subject to hi(Y ) = 0; i = 1; ¢ ¢ ¢ ; n; gi(Y ) · 0; i = 1; ¢ ¢ ¢ ; m: To solve the problem, we ¯rst form the logarithmic barrier function as B = f(Y ) ¡ º Xm i=1 ln (¡gi(Y )) ; (4.16) where the parameter º is referred to as the barrier parameter, a positive number that is reduced to approach to zero as the algorithm converges to the optimum. Then we solve a sequence of constrained minimization problems of the form min (Y;¸;º) B(Y; ¸; º) (4.17) subject to hi(Y ) = 0; i = 1; ¢ ¢ ¢ ; n; for a sequence fºkg of positive barrier parameters that decrease monotonically to zero. The solution of this problem by Newton's method requires the formulation of the Lagrangian function L = f(Y ) ¡ Xn i=1 ¸ihi(Y ) ¡ º Xm i=1 ln (¡gi(Y )) : (4.18) Because the barrier term is in¯nite on the boundary of the feasible region, as shown in Fig. 4.1, it acts as a repelling force that drives the current trial solution away from the boundary into the interior of the feasible region. As the barrier parameter º is decreased, the e®ect of the barrier term is diminished, so that the iterates can grad ually approach the constraint boundaries of the feasible region for those constraints which eventually turn out to be binding. 32 Figure 4.1: E®ect of barrier term To solve equation (4.17), we need to ¯nd the gradient of L with respect to Y : rY L = rY f(Y ) ¡ Xn i=1 ¸irY hi(Y ) + º Xm i=1 1 ¡gi(Y ) rY gi(Y ) = rY f(Y ) ¡ Xn i=1 ¸irY hi(Y ) + º Xm i=1 1 si rY gi(Y ) = rY f(Y ) ¡ Xn i=1 ¸irY hi(Y ) + Xm i=1 ¹irY gi(Y ) = rY f(Y ) ¡ hTY ¸ + gT Y ¹; (4.19) where ¡gi(Y ) = si; i = 1; ¢ ¢ ¢ ; m; ¹isi = º; i = 1; ¢ ¢ ¢ ; m; hTY is a matrix consisting of the gradient rY hi(Y ) as columns, and gT Y is a matrix consisting of the gradient rY gi(Y ) as columns, or hTY = · rY h1(Y ) rY h2(Y ) ¢ ¢ ¢ rY hn(Y ) ¸ ; gT Y = · rY g1(Y ) rY g2(Y ) ¢ ¢ ¢ rY gm(Y ) ¸ : Also, ¸ is a vector containing the elements ¸i, and ¹ is a vector containing the elements 33 ¹i. The set of equations we must solve is rY f ¡ Xn i=1 ¸irhi + Xm i=1 ¹irgi = 0; (4.20) hi(Y ) = 0; i = 1; ¢ ¢ ¢ ; n (4.21) gi(Y ) + si = 0; i = 1; ¢ ¢ ¢ ;m (4.22) ¹isi ¡ º = 0; i = 1; ¢ ¢ ¢ ;m (4.23) si > 0; i = 1; ¢ ¢ ¢ ;m (4.24) ¹i > 0; i = 1; ¢ ¢ ¢ ;m (4.25) These equations can be solved using the Newton iterative method. Our four sets of variables for which we must solve are Y; ¸; ¹ and s. The equation for a ¯rst order approximation of a Taylor series of a function, F, whose independent variables are Y; ¸; ¹ and s is F(Y; ¸; ¹; s) »= F(Yo; ¸o; ¹o; so) + @F @Y 4Y + @F @¸ 4¸ + @F @¹ 4¹ + @F @s 4s: (4.26) We want to ¯nd the values of Y; ¸; ¹ and s where the expressions on the left side of the equations we want to solve evaluate to zero. We use NewtonRaphson to do so. Taking the ¯rst order approximation to the Taylor series for each of the four expressions given in equations (4.20)(4.23) and setting them equal to zero (the desired value for each) give us ¡ rY f ¡ hTY¸ + gT Y ¹ ¢ + (r2 Y f ¡ Xn i=1 ¸ir2 Y hi + Xm i=1 ¹ir2 Y gi)¢Y ¡hTY ¢¸ + gT Y¢¹ = 0; (4.27) h + hY¢Y = 0; (4.28) (g + s) + gY¢Y + I¢s = 0; (4.29) (MSe ¡ ºe) + S¢¹ +M¢s = 0; (4.30) 34 where I : an identity matrix; S : a diagonal matrix constructed from (s1; s2; ¢ ¢ ¢ ; sm); M : a diagonal matrix constructed from (¹1; ¹2; ¢ ¢ ¢ ; ¹m); e : a column vector with all elements 1. The above equations can be organized in matrix form as 2 6666666664 HY hTY gT Y 0 hY 0 0 0 gY 0 0 I 0 0 S M 3 7777777775 2 666666664 ¢Y ¢¸ ¢¹ ¢s 3 777777775 = 2 6666666664 ¡rY f + hTY ¸ ¡ gT Y ¹ ¡h ¡g ¡ s ºe ¡MSe 3 7777777775 ; (4.31) or W¢z = ¢F; (4.32) where HY = r2 Y f ¡ Xn i=1 ¸ir2 Y hi + Xm i=1 ¹ir2 Y gi: We initially set º to some relatively large number, such as 10. Starting with an initial guess zo = 2 66666664 Yo ¸o ¹o so 3 77777775 ; we calculate the W matrix and the ¢F vector. If all the elements of ¢F are su± ciently close to zero, we have found the solution zo. Otherwise we must solve for ¢z by ¢z = W¡1¢F: (4.33) 35 Then the original z is updated by z = zo + ®¢z: (4.34) However, following conditions need to be satis¯ed when updating z: ¹i > 0; si > 0: Therefore, when calculating the updated ¹ and s, we must make sure that each ¹i and each si is still strictly greater than zero when ¢¹i and ¢si are added to it, respectively. If adding ¢¹i or ¢si violates this condition, all ¢¹i and all ¢si must be scaled by some factor ® less than one before adding them as in (4.34). Now that we have a new guess for z, the W matrix and ¢F vector are calculated again. If the elements of the ¢F vector are su±ciently close to zero, we have found the z which solves our problem. Otherwise we solve for ¢z and update z again. This process is repeated until we ¯nd the z which makes ¢F very close to zero. After we solve for z, we reduce the barrier parameter º by some factor ·. For example, if · = 0:4, the barrier parameter º is updated by letting the new º be 0.4 times the old º. Using the value of z obtained using the old º as an initial guess zo, we use the iterative procedure described above to solve for the value of z which makes ¢F approximately zero for the new barrier parameter. The process of solving for z, decreasing º, and then solving for z again is repeated until º becomes a very small number, such as 10¡10. This process is illustrated in Figure 4.2. When º gets this small, we have found the z that solves our problem. When solving for z given a particular º, it is not necessary to force ¢F to be zero. What we really want to know is how close we are to the central path. The central path is de¯ned by a sequence of solutions fY (º); ¸(º); ¹(º); s(º)g which make ¢F evaluate to zero for every possible value of º. One way of measuring how close we are to the central path is by checking to see how close each product ¹isi is to the barrier parameter º. One 36 Figure 4.2: Graphical representation of central path proposed way of doing this is by ¯rst calculating the average value of ¹isi, also known as the dualityfactor [13], by bº = Xm i=1 ¹isi m ; (4.35) which evaluates to zero when each product ¹isi is equal to º. Instead of checking to see if ¢F is su±ciently close to zero, we check to see when the following is true: °°°°°°°°°°°°° ¹1s1 ¡ bº ¹2s2 ¡ bº ... ¹msm ¡ bº °°°°°°°°°°°°° 2 · ¿º (4.36) where 0 < ¿ < 1: When this logical statement becomes true, we are su±ciently close to the central path, and we can reduce º. Despite several attractive features of the PDIP method, it has an inherent disad vantage that the size of problem is much bigger than the Newton active set method since the PDIP method includes all inequality constraints. This problem becomes 37 even more involved when we consider a large scale natural gas and power system network. 38 CHAPTER 5 NATURAL GAS FLOW MODELING Natural gas is transported from gas producers to customers at various locations. A typical natural gas transmission system today consists of a large number of gas pro ducers, various customers, storage, many compressor stations, thousands of pipelines, and many other devices such as valves and regulators, including midstream gas pro cessing (between well and pipeline) A typical pipeline network for transporting natural gas consumes a signi¯cant amount of fuel per day to operate compressors pumping natural gas. This is because there is a gas pressure loss due to friction between gas and pipe inner walls. Moreover, energy is lost by heat transfer between gas and its environment. To compensate for these losses of energy and to keep the gas °owing, compressor stations are installed in the network, which consume a part of the transported gas resulting in economic losses1. The purpose of this chapter is to provide the underlying mathematical modeling of natural gas transmission networks. It will include gas °ow equations, compressor horsepower equation, and matrix representations of natural gas networks. 1Losses often refer to leakage of natural gas from the pipeline system, which is negligible if there is no theft problem. We will use the term \losses" only for the gas required to power the compressors. 39 Figure 5.1: Pipeline network representation 5.1 Elements of Natural Gas Transmission Net work Three basic types of entities are considered for the modeling of a natural gas trans mission network: pipelines, compressor stations, both of which are represented by branches, and interconnection points, represented by nodes. For simulation of a network, we assume that nodes represent pipeline connections while branches represent pipelines and compressor stations, which have °ow directions assigned. We also have di®erent types of nodes. A source node represents a gas production or storage facility. A load node represents a place where gas is to be taken out of the system, either for consumption or storage. Each compressor branch de¯nes two more node types technically referred to as suction and discharge nodes. Likewise, each pipeline branch has two types of nodes, a sending end node and a receiving end node. We de¯ne fkij (or fk for simple notation) as the °ow rate through a branch, say number k, which begins at node i and ends at node j. Figure 5.1 shows an example of a natural gas network. It is a simple treestructure pipeline network, which consists of one source at node 1, two loads at nodes T + 1 & K + 1, and several pipeline and compressor branches. wS1 is the amount of gas supplied at the source node, and wLT+1 & wLk+1 are the amounts of gas consumed at the load nodes, respectively. 40 For a given gas network, there are three types of relevant decision variables: the °ow rate through a pipeline, the °ow rate through a compressor, and the pressure at each node. For a compressor, these variables are further restricted by a set of constraints that depend on the operating attributes of the compressor. 5.2 Network Topology Analysis of natural gas pipeline networks is relatively complex, particularly if the network consists of a large number of pipelines, several compressors, and various sup pliers and consumers [21]. Matrix notation is a simple and useful way of representing a network. In gas network analysis, matrices turn out to be the natural way of ex pressing the problem [22]. A gas pipeline system can be described by a set of matrices based on the topology of the network. Consider the gas network represented by the graph in Figure 5.2 [22]. This network consists of one source at node 1, three loads at nodes 2, 3 & 4 and ¯ve pipeline branches ¬, , ®, ¯ & °. For network analysis it is necessary to select at least one reference node. Mathe matically, the reference node is also referred to as an independent node, and all nodal and branch quantities are dependent on it. The pressure at the reference node is Figure 5.2: Graph of the gas network [22] 41 known. A network may contain several pressurede¯ned nodes and these form a set of reference nodes for the network. A gas injection node is a point where gas is injected into the network, which may be positive, negative, or zero. Negative injection represents a load demand for gas from the network. This node may be supplying domestic or commercial consumers, charging gas storages, or even accounting for leakage in the network. A positive injection represents a supply of gas to the network. It may take gas from storage, source or another network. A zero injection is assigned to nodes that do not have a load or source but are used to represent a point of change in the network topology, such as a junction of several branches. Then, the vector of gas injections w at each node for Figure 5.2 is de¯ned as w = · w1 w2 w3 w4 ¸T ; = · wS1 ¡wL2 ¡wL3 ¡wL4 ¸T ; where wi = wSi ¡ wLi ; wSi = gas injection at node i; wLi = gas removed at node i: At steadystate conditions, the total load on the network is balanced by the supply into the network at the source node. To de¯ne the network topology completely, it is necessary to assign a direction to each branch. Each branch direction is assigned arbitrarily and is assumed to be the positive direction of °ow in the branch. If the °ow has a negative value, then the direction of °ow is opposite to the branch direction. References such as Osiadacz's book [22] may be consulted for further detail. 42 5.3 Matrix Representations of Network The interconnection of a network can be described by the branchnodal incidence matrix A. This matrix is rectangular, with the number of rows NN equal to the number of nodes (including reference nodes), and the number of columns NP equal to the number of pipeline and compressor branches in the network. The element Aij of the matrix A corresponds to node i and branch j, and is de¯ned as Aij = 8>>>>< >>>>: +1; if pipeline branch j enters node i; ¡1; if pipeline branch j leaves node i; 0; if pipeline branch j is not connected to node i: For the network in Figure 5.2, the branchnodal incidence matrix is A = 2 66666664 ¡1 ¡1 ¡1 0 0 1 0 0 1 0 0 1 0 ¡1 ¡1 0 0 1 0 1 3 77777775 ; where NN = 4 nodes and NP = 5 branches. One important thing to note is that the sum of all rows in the matrix A becomes zero due to the de¯nition of each element of the matrix A. Thus, the matrix A is always rank1 de¯cient, and therefore the matrix AAT is singular. This problem can be avoided by removing the rows of the matrix A corresponding to knownpressure nodes, which will be discussed in next chapter. 5.4 Flow Equation For isothermal gas °ow in a long horizontal pipeline, say number k, which begins at node i and ends at node j, the general steadystate °ow rate (in standard ft3/hr, or SCF/hr at T0 = 520oR and P0 = 14:65 psia) is often expressed by the following 43 formula [22, 23] derived from energy balance: fkij = Sij £ 3:22 T0 ¼0 s Sij ¡ ¼2 i ¡ ¼2 j ¢ D5 k FkGLkTkaZa ; (5.1) where fkij = pipeline °owrate (SCF/hr); Sij = 8>< >: +1 if ¼i ¡ ¼j > 0; ¡1 if ¼i ¡ ¼j < 0; Fk = pipeline friction factor; Dk = internal diameter of pipe between nodes (inch); G = gas speci¯c gravity (air=1.0, gas=0.6); Lk = pipeline length between nodes (miles); ¼i = pressure at node i (psia); ¼j = pressure at node j (psia); ¼0 = standard pressure (psia); T0 = standard temperature (±R); Tka = average gas temperature (±R); Za = average gas compressibility factor. There are several di®erent °ow equations in use in the natural gas transmission indus try [22, 24]. The di®erences are mainly due to the empirical expression assumed for the friction factor, Fk, and in some cases the rigorous consideration of the deviation of the behavior of natural gas from that of an ideal gas. The change of °ow changing from partial turbulence to full turbulence is referred as the transition region. The Reynolds number, a measure of the ratio of the inertia force on an element of °uid to the viscous force on an element, at which this transition occurs is dependent on the diameter of the pipeline and its roughness, and is typically about 107 [22]. The extent of the transition region is dependent on the system considered, and within this region 44 the frictional resistance depends on both the Reynolds number and the pipe charac teristics. However, in the fullyturbulent °ow region for highpressure networks, the friction factor Fk is strictly dependent on pipeline diameter [22], [25]. That is Fk = 0:032 D1=3 k : (5.2) Then, equation (5.1) becomes fk = fkij = Sij Mk q Sij ¡ ¼2 i ¡ ¼2 j ¢ ; SCF/hr (5.3) where Mk = ² 18:062 T0 D8=3 k ¼0 p GLkTkaZa ; ² = pipeline e±ciency: As indicated in equation (5.3), gas °ow can be determined once ¼i and ¼j are known for given conditions. Equation (5.3), known as Weymouth °ow equation, is most satisfactory for large diameter (¸ 10 inches) lines with high pressures [23]. 5.5 Compressor Horsepower Equation During transportation of gas in pipelines, gas °ow loses a part of its initial energy due to frictional resistance which results in a loss of pressure. To compensate the loss of energy and to move the gas, compressor stations are installed in the network. In general, the nature of the compressor work function is very complex and depends also on consideration such as the number of compressors running within the compressor station, how the compressor units are con¯gured (i.e. in series, in parallel, combina tion of both, etc.), physical properties of the compressor units, and type of compressor unit [22]. One of the most common con¯gurations implemented in practice is that of compressor stations consisting of identical centrifugal compressor units operating in parallel. Centrifugal compressors are versatile, compact, and generally used in the 45 range of 1,000 to 100,000 inlet ft3=min for process and pipeline compression applica tions [24]. In a centrifugal compressor, work is done on the gas by an impeller. Gas is discharged at a high velocity into a di®user. The velocity of gas is reduced and its kinetic energy is converted to static pressure. A key characteristic of the centrifugal compressor is the horsepower consumption, which is a function of the amount of gas that °ows through the compressor and the relative boost ratio between the suction and the discharge pressures. After empirical modi¯cation to account for deviation from ideal gas behavior, the actual adiabatic (zero heat transfer) compressor horsepower equation [23] at To = 60oF (= 520oR) and ¼o = 14:65 psia becomes Hk = Hkij = Bk fk "µ ¼j ¼i ¶Zki(®¡1 ® ) ¡ 1 # ; HP (5.4) where Bk = 3554:58 Tki ´k µ ® ® ¡ 1 ¶ ; fk = °ow rate through compressor (SCF/hr); ¼i = compressor suction pressure (psia); ¼j = compressor discharge pressure (psia); Zki = gas compressibility factor at compressor inlet, Tki = compressor suction temperature (±R); ® = speci¯c heat ratio (cp=cV ); ´k = compressor e±ciency: 5.6 Conservation of Mass Flow The mass°ow balance equation at each node can be written in matrix form as (A + U)f + w ¡ T¿ = 0; (5.5) 46 where A = a branchnodal incidence matrix; Uik = 8>>>>< >>>>: +1; if the kth unit has its outlet at node i; ¡1; if the kth unit has its inlet at node i; 0; otherwise. Tik = 8>< >: +1; if the kth turbine gets gas from node i; 0; otherwise. f = a vector of mass °ow rates through branches; w = a vector of gas injections at each node: In addition to the matrix A, which represents the interconnection of pipelines and nodes, we de¯ne the matrix U, which describes the connection of units (compressors) and nodes. The vector of gas injections w is obtained by w = wS ¡ wL; (5.6) where wS = a vector of gas supplies at each node; wL = a vector of gas demands at each node: Thus, a negative gas injection means that gas is taken out of the network. The matrix T and the vector ¿ represent where gas is withdrawn to power a gas turbine to operate the compressor. So if a gas compressor, say k, between nodes i and j, is driven by a gas¯red turbine, and the gas is tapped from the suction pipeline i, we have the following representation: Tik = +1; Tjk = 0 ; and ¿k = amount tapped: Conversely, if the gas were tapped at the compressor outlet, we would have Tik = 0 ; Tjk = +1; and ¿k = amount tapped: 47 Analytically, we will assume that ¿k can be approximated as ¿k = ®Tk + ¯TkHkij + °TkH2 kij ; (5.7) where Hk = Hkij is the horsepower required for the gas compressor k in equation (5.4). 48 CHAPTER 6 NATURAL GAS LOADFLOW The problem of simulation of a gas network with NN nodes in steady state, known as load°ow, is usually that of computing the values of node pressures and °ow rates in the individual branches for known values of NS source pressures (NS ¸ 1) and of gas injections in all other nodes. Gas load°ow analyses are required operationally whenever signi¯cant changes in demands or supplies are expected to occur. It is also used for system planning pur poses. For example, when a gas¯red generator is located in the gas network, we need to see whether the network has enough capability to carry the required amount of gas to the generator while satisfying various network constraints, such as pressure limits at each node and compressor operation limits. In this chapter, we state the load°ow problem in a general way, and construct a mathematical formation. We present a load°ow problem for a network without compressors. Then we introduce a general load°ow analysis with compressors. We formulate the gas load°ow problem in a way similar to the electric load°ow problem, and include the gas consumption rates at compressor stations. However, we omit gas storage to avoid inter  temporal linkageswe only attempt to solve the gas load°ow at a single timea \snapshot". 49 6.1 Load°ow Problem Statement The gas load°ow problem is stated below: ² Given a natural gas system described by a branchnodal incidence matrix A, a unitnodal incidence matrix U, a gas turbinenodal incidence matrix T, and given a set of gas injections except at the NS knownpressure sources (injections at these nodes initially unknown), and each unit's operating condition (such as the compression ratio, the °ow rate through the compressor, or the suction or discharge pressure), ² determine all other pressures, and calculate the °ow rates in all branches and the gas consumptions at compressor stations. Simply speaking, one of two quantities, nodal pressure ¼i and gas injection wi at each node, and one compressor operating condition are speci¯ed, and other values are to be determined. Speci¯ed quantities are chosen based on the following conditions: Nodes: ² KnownInjection Node: For a node i of this type, we assume that we know a gas injection wi, and the pressure ¼i is to be determined. Generally, source and load nodes, and junctions with no gas injections belong to this node. Elec trically, this is analogous to a \load bus." In fact, solving the load°ow problem with only this type of node is not in general possible. The ¯rst reason is that, in the °ow equation (5.1), the pressures never appear by themselves, but instead appear only as a squaredpressure di®erence of the form ¼2 i ¡ ¼2 j . Therefore, there are only NN ¡ 1 pressures which a®ect the load°ow. We therefore pick NS ¸ 1 nodes to provide reference pressures. These nodes are generally the external gas sources supplying our system. 50 Another reason is that with is that solving a load°ow for a network containing only knowninjection nodes would imply that we know the gas injections at every single node. In fact, we cannot mathematically specify all NN gas injections, as it may not be possible to ¯nd a solution to the load°ow equations. Specifying the injections at all nodes is the same as specifying the gas supplies to gas turbines driving gas compressors, which we cannot know until the load°ow is solved. Instead, we must pick at least one node, allowing the set of gas injection(s) to be whatever is required to solve the load°ow equations. Thus, we have to specify another node type: ² KnownPressure Node: Each is typically one of the source nodes, and the pressures of such nodes serve as references for all other pressures. We assume that we know f¼i; i = 1; : : : ;NSg, but we do not know the corresponding gas injections. Electrically this is analogous to a (possibly distributed) \slack bus." In addition to nodes, the other main components are branches, which connect the nodes. Branches: ² Pipelines: Pipeline °ow modelling has already been discussed in the previous chapter. Other than the physical characteristics of the pipeline, the only vari ables that the °ow fk = fkij on pipeline k depends on are the pressures ¼i and ¼j at the two ends of the pipeline. ² Compressors: The other key component we will model in a gas network is a compressor (also called a unit). The connection between the unit's inlet and outlet nodes is not de¯ned by the branchnodal incidence matrix A, but by U. The compression ratio between the compressor inlet and outlet, and the °ow rate through the compressor are governed by the horsepower equation (5.4), not by the °ow equation (5.1). Compressor data (other than the physical charac teristics of the compressor) can be speci¯ed in several ways [22] for compressor 51 k: relative boost Rk = Rkij = ¼j=¼i, or absolute boost ¼j ¡ ¼i, or mass°ow rate fkij . The inlet pressure ¼i or the outlet pressure ¼j could also be speci¯ed. 6.2 Load°ow without Compressors We will ¯rst consider a gas network with only pipelines to explain the Newtonnodal method. Assume that node 1 is the knownpressure node, and all other nodes are the knowninjection nodes. Since we do not have compressors, the unit's inletoutlet nodes are not de¯ned. Then the load°ow problem is described as follows: Given : ¼1 w2 : : : wNN Determine : w1 ¼2 : : : ¼NN It appears that we have NN ¡ 1 quantities known, and NN ¡ 1 quantities unknown. The set of nodal °ow equations that describes a gas network with only pipelines is given by ¹ w = A1 ¢ f(¹¼; ~¼); (6.1) where A1 = a branchnodal incidence matrix except the knownpressure nodes; ¹ w = ¹ ws ¡ ¹ wL; f = a vector of °ow rates through pipelines: We used ¹¼ to indicate the part of ¼ we know, and ~¼ for unknown part. Likewise for ¹ w. Flow equation (5.3) is used to get each pipeline branch °ow as fkij = Sij Mk q Sij ¡ ¼2 i ¡ ¼2 j ¢ ; SCF/hr (6.2) In the Newtonnodal method [22], an initial approximation is made to the nodal pressures. This approximation is then iteratively corrected until the ¯nal solution is 52 reached. At each iteration the righthand side of equation (6.1) is not equal to zero. The pressures are only approximations of their true values and the °ows calculated from these pressures are not balanced at each node. The imbalance at a node is the nodal error which is a function of all the nodal pressures (except the knownpressure nodes). The Newtonnodal method solves the set of equations (6.1) iteratively until the nodal °ow errors are small enough to be insigni¯cant. The iterative scheme for correcting the approximations to the nodal pressures is ~¼k+1 = ~¼k + ¢~¼k; (6.3) where k = the number of iterations. Term ¢~¼ is computed from the following equa tion: Jk ¢ ¢~¼k = ¡(A1 ¢ f(¹¼; ~¼) ¡ ¹ w)k; (6.4) where the matrix J is the nodal Jacobian matrix and is given by [22] J = ¡A1DAT 1¦1; (6.5) and D = diag µ fkij ¼2 i ¡ ¼2 j ¶ ; k = 1; : : : ;NP ¦1 = diag(~¼): Since the matrix A1DAT 1 is symmetric and positive de¯nite, and ¦1 is a diagonal matrix, NewtonRaphson algorithm can be implemented with a great computational e±ciency [26], [21]. 6.3 Load°ow with Compressors We assume that the pressures at the NS knownpressure nodes are known, and that the injections at the knowninjection nodes are speci¯ed. Also, some operating pa rameter for each compressor, say (for purposes of illustration) the relative boost 53 Rk = Rkij , is speci¯ed, and let's say there are NP branches in the system, of which NC are compressors. We can state the load°ow problem this way: Given : ¼1 ¢ ¢ ¢¼NS wNS+1: : :wNN R1¢ ¢ ¢RNC Find : w1 ¢ ¢ ¢wNS ¼NS+1 : : :¼NN f1 ¢ ¢ ¢ fNC fNC+1: : :fNP It appears that there are NN +NC quantities given, while NN +NP must be found. It is clear that NC < NP : This inequality is strict, unless we have no pipelines at all, only NP compressors connected to each other without any intervening pipelines { a silly situation. So since NN + NC < NN + NP , the system appears to be undetermined. Note, however, that from (6.2), the °ow fk depends only on the pressures ¼i and ¼j of the nodes it connects. Likewise, the horsepower Hk required by the compressor depends only on the °ow fk and the ratio Rk = ¼j=¼i, and therefore only depends on the pressures ¼i and ¼j . The tapo® loss ¿k depends only on Hk and thus on the nodal pressures. So, if we knew f¼i; i = 1; : : : ;NNg, we would know all other quantities we have discussed. But we only know them for i = 1; : : : ;NS. Let's use ¼ to indicate the part of ¼ we know, and ~¼ for the unknown part. Likewise for w and ~ w. The objective, then, is to calculate ~¼, giving us the entire vector ¼, from which all other quantities can be calculated. We can use the massbalance equation (5.5) to write w = ~T ¿ (¼; ~¼) ¡ (~A + ~U )f(¼; ~¼); where ~A , ~ U, and ~T are obtained from A, U, and T as in equation (5.5) after we have removed the rows corresponding to the NS nodes (the knownpressure nodes with unknown injections), that is, we have dropped the ¯rst NS equations of (5.5) 54 corresponding to ~ w. This gives us NN ¡ NS equations (for the elements of w) in NN ¡NS unknowns (the elements of ~¼). At this point, standard NewtonRaphson or other iterative methods can be employed to drive the \mismatch" ¢w (the di®erence between the values of w computed from above and the true values of w given as part of the data) to zero by correcting our current guess for the correct value of ~¼. We have shown that the dimension of the load °ow problem with compressors is NN ¡ NS. However, it is convenient to add new equations and new variables if they make the problem easier to implement. Since the compressor horsepower in equation (5.4) and the gas consumption rate in equation (5.7) are functions of horsepowers, compressor °ow rates, and gas consumption rates as well as the speci¯ed relative boost rates, we will add three extra decision variables for each compressor station. Then the new decision variables become x = · ~¼T HT ¿ T fT C ¸T ; where H = a vector of horsepowers at each compressor, ¿ = a vector of gas consumptions at each compressor, fC = a vector of compressor branch °ow rates. Since there are NN ¡ NS + 3NC decision variables, we need to de¯ne extra 3NC equations. The mass°ow balance equations F1 in matrix form are given by F1 = (~A + ~U ) ¢ 2 64 fC fP (¹¼; ~¼) 3 75 + ¹ w ¡ ~T ¢ ¿ = 0; (6.6) where fP (¹¼; ~¼) is a vector of mass °ow rates through pipelines, which are obtained by °ow equation (6.2). Note that the size of the vector F1 is NN ¡ NS. Assuming that the compression ratio Rkij is given, and gas compressor k located between nodes i and j is driven by a gas¯red turbine, the extra 3NC equations related with compressor 55 operation are F2k = ¿k ¡ ®Tk ¡ ¯TkHkij ¡ °TkH2 kij = 0; k = 1; : : : ;NC (6.7) F3k = Hkij ¡ Bk fCk "µ ¼j ¼i ¶Zki(®¡1 ® ) ¡ 1 # = 0; k = 1; : : : ;NC (6.8) F4k = µ ¼j ¼i ¶ ¡ Rkij = 0; k = 1; : : : ;NC (6.9) Let us de¯ne the vector of error functions F(x) by F(x) = 2 66666664 F1(x) F2(x) F3(x) F4(x) 3 77777775 : Note that there are NN ¡NS +3NC equations and NN ¡NS +3NC unknown decision variables. So we have the right number of variables to force all the mismatches to zero by using NewtonRaphson or other iterative methods. The iteration will be repeated until the mismatches become small enough to be insigni¯cant. 56 CHAPTER 7 UPFC MODELING FOR STEADYSTATE ANALYSIS Modeling UPFC for steadystate analysis has been considered by several researchers and an injection model [6] and an uncoupled model [14] have been proposed. These models can be easily incorporated into steadystate load°ow or optimal power °ow studies. However, these models employ four UPFC control variables that depend on the UPFC input and output currents and voltages and both models require adding two additional buses to the load°ow or OPF problem formulation. The voltage and current relationships between UPFC input and output need to be included explicitly. In addition, a constraint on real power conservation needs to be added, thereby reducing the degrees of freedom four to three. Since electricity prices at the UPFC input and output buses, which are the dual variables (or \shadow" prices) associated with real and reactive power injections, are meaningless when no power is bought or sold at these ¯ctitious buses, their addition to the problem only serves to increase the size of the problem. Therefore, these models are undesirable for our UPFC sensitivity analysis. To overcome these problems, a new steadystate mathematical model for a UPFC, the UPFC ideal transformer model, is proposed. In this model, UPFC control vari ables do not depend on UPFC input and output voltages and currents, and therefore addition of ¯ctitious input and output buses are not necessary. This model is eas ily combined with transmission line models using ABCD twoport representations, 57 which can then be converted to Yparameter representations. Thus, UPFC model is embedded in the ¹Y bus matrix, and so the size of the ¹Y bus matrix is not changed. 7.1 Operating Principles The Uni¯ed Power Flow Controller (UPFC) is one of the most technologically promis ing devices in the FACTS family [4, 5, 6, 7]. It has the capability to control voltage magnitude and phase angle, and can also independently provide (positive or negative) reactive power injections. A UPFC consists of a shunt transformer, a series trans former, power electronic switching devices and a DC link, as shown in Figure 7.1 [7]. Inverter 1 is functionally a static VAR compensator assuming that inverter 2 is not connected. It injects reactive power in the form of current at the shunt transformer, and the current phasor ~IT is in quadrature to the input voltage ~VI . Inverter 2 by itself represents the socalled advanced controllable series compensator (ACSC) assuming that inverter 1 is not connected. It injects reactive power by adding voltage through the series transformer. The injected voltage ~V T is in quadrature to the receiving end Figure 7.1: General UPFC scheme [7] 58 Figure 7.2: Proposed UPFC model in a transmission line current ~Io. Now if we connect inverter 1 to inverter 2 through a DC link, inverter 1 can provide real power to inverter 2. Therefore the UPFC can independently control real and reactive power injections through the series transformer, but the real power injected at the series transformer is provided by the shunt transformer through the DC link. Inverter 1 must provide the real power used by inverter 2 via the DC link, but can also independently inject reactive power (positive or negative) through the shunt transformer. In summary, note that the UPFC conserves real power but can still generate (or sink) reactive power at either transformer or both. 7.2 Uncoupled Model Figure 7.2 shows a basic UPFC model, where the UPFC is located between buses i and k. Each part of the transmission line is represented as an equivalent ¦ circuit. Figure 7.3 shows a phasor diagram illustrating UPFC inputoutput voltage and current re lationships. The injected series voltage ~V T can be resolved into inphase component Vp and quadrature component Vq with respect to the UPFC output current ~Io, and can be expressed as ~V T = (Vp + jVq) ej±o : (7.1) 59 Figure 7.3: Phasor diagram of UPFC inputoutput voltages and currents Since ~V T is dependent on the UPFC output current phase angle ±o, it requires adding an extra bus for the UPFC output terminal. The current IT injected by the shunt transformer contains a real component Ip, which is in phase or in opposite phase with the input voltage. It also has a reactive component Iq, which is in quadrature with the input voltage. Then the injected current ~IT can be written by ~IT = (Ip + jIq) ejµI ; (7.2) where µI is the UPFC input voltage phase angle. Thus, a second extra bus is required for the UPFC input terminal. The UPFC inputoutput voltage and current can be represented by ~V o = ~V I + ~V T = VIejµI + Vpej±o + jVqej±o ; (7.3) ~Io = ~II ¡ ~IT = IIej±I ¡ IpejµI ¡ jIqejµI ; (7.4) where ±I is the UPFC input current phase angle. Then, the complex power injected into the transmission line by the series transformer can be resolved into the real and reactive power in simple form as ST = ~V T ¢ ~I¤ o = V p{¢zI}o PT + jVq{¢zI}o QT : (7.5) 60 Figure 7.4: Uncoupled UPFC model in a transmission line. The inphase voltage Vp is associated with a real power supply and the quadrature voltage Vq with an inductive or capacitive reactance in series with the transmission line. Since the real power PT (which may be negative) is provided by the current Ip in the shunt transformer, we can derive the following relationship: Vp ¢ Io ¡ VI ¢ Ip = 0: (7.6) or the real power input equals the real power output. Due to (7.6), the number of degrees of freedom for the UPFC is reduced to three. Now that two extra buses for the UPFC input and output terminals are added, and the UPFC voltage and current relationships (7.3, 7.4) and real power °ow equation (7.6) are established, we can represent the UPFC uncoupled model as shown in Figure 7.4. The currents injected into the UPFC input and output buses are ~II = µ Yi 2 + 1 Zi ¶ ~V I ¡ 1 Zi ~Vi ; (7.7) ~Io = µ Yk 2 + 1 Zk ¶ ~V o ¡ 1 Zk ~ Vk : (7.8) The magnitudes of the injected voltage VT and current IT are limited by the maximum voltage and current ratings of the inverters and their associated transformers, which need to be included as inequality constraints in OPF. 61 7.3 Ideal Transformer Model Since the UPFC conserves real power, and generates or consumes reactive power, it can be modelled using an ideal transformer and a shunt branch, as shown in Figure 7.5 [27]. The advantage of this model is that the ideal transformer turns ratio and the variable shunt susceptance are independent variables, which are not directly as sociated with the UPFC inputoutput voltages and currents. We de¯ne the UPFC variables as follows: T = transformer voltage magnitude turns ratio (real); Á = phase shifting angle; ½ = shunt susceptance; and the ideal transformer turns ratio can be written by ¹ T = TejÁ: It is important to note that the ideal transformer does not generate real and reactive power, and the reactive power is generated (or consumed) by the shunt admittance only. Since the UPFC inputoutput voltage and current relationship can be expressed Figure 7.5: UPFC ideal transformer model 62 as ~V I = ~V oT\Á; (7.9) ~II = j ¹ T½~V o + 1 ¹ T¤ ~Io; (7.10) the UPFC can be represented by an ABCD matrix as 2 64 ~V I ~II 3 75 = ABCDU ¢ 2 64 ~V o ~Io 3 75 ; (7.11) where ABCDU = 2 64 ¹ T 0 j ¹ T½ 1 ¹ T¤ 3 75 : (7.12) Note that equation (7.11) is not bilateral unless ¹ T = 1\0. Now, we will show that this ideal transformer model represents the UPFC by comparing the complex power injections at the UPFC input and output. Using (7.11), the complex power injection at the UPFC input can be obtained by ¹ SI = ~V I ~I¤ I ; (7.13) = ¹ T~V o µ j ¹ T½~V o + 1 ¹ T¤ ~Io ¶¤ ; = ~V o~I¤ o ¡ jj ¹ Tj2 ¢ j~V oj2½; = ¹ So ¡ jj ¹ Tj2 ¢ j~V oj2½; and the real and reactive power injections can be obtained by PI = Real ¡ ¹ SI ¢ ; QI = Imag ¡ ¹ SI ¢ ; Po = Real ¡ ¹ So ¢ ; Qo = Imag ¡ ¹ So ¢ : Thus, we can derive the following relationships between the UPFC input and output: PI = Po; (7.14) Qo = QI + j ¹ Tj2 ¢ j~V oj2½: (7.15) 63 Figure 7.6: Simpli¯ed UPFC circuit Equations (7.14) and (7.15) mean that the ideal transformer model conserves real power and generates or consumes (for ½ < 0) reactive power. To determine how much real and reactive power is injected in the series and shunt transformers, we will map the complex turns ratio ¹ T in the ideal transformer and the shunt susceptance ½ to the injected voltage ~V T and current ~IT in the UPFC uncoupled model. Since the UPFC input voltage and current are expressed as ~V I = ~V o ¡ ~V T = ~V o Ã 1 ¡ ~V T ~V o ! = ~V oT\Á; (7.16) ~II = ~Io + ~IT = ~Io + µ 1 T \Á ¡ 1 ¶ ~Io + j½~V I ; (7.17) the injected voltage ~V T and current ~IT can be obtained by ~VT = (1 ¡ T\Á)~V o; (7.18) ~IT = µ 1 T \Á ¡ 1 ¶ ~Io + j½~V I : (7.19) Then, the power °ows through each inverter, as shown in Figure 7.6, can be obtained 64 by ¹ S1 = ~V I ~I¤ T ; = ~V oT\Á ·µ 1 T \Á ¡ 1 ¶ ~Io + j½~V oT\Á ¸¤ ; = (1 ¡ T\Á) ¹ So ¡ j½j ¹ Tj2j~V oj2; (7.20) ¹ S2 = ¡~V T ~I¤ o ; = (T\Á ¡ 1)~V o~I¤ o ; = (T\Á ¡ 1) ¹ So: (7.21) Thus, ¹ S1 + ¹ S2 = ¡j½j ¹ Tj2j~V oj2; (7.22) which veri¯es that the UPFC conserves real power and can generate (or consume) reactive power. Since the UPFC is modelled using passive circuit elements only, nonideal UPFC characteristics, such as shunt and series transformer reactances can be easily incor porated into this framework. 7.4 UPFC in a Transmission Line A twoport ABCD matrix is the most convenient method to represent cascaded net works [9]. Let us divide a transmission line between buses i and k with a UPFC into three cascaded networks, a UPFC input transmission line, a UPFC, and a UPFC output transmission line, as shown in Figure 7.7. The UPFC input transmission line , and the UPFC output transmission line are easily represented by twoport ABCD matrices since the transmission lines are modelled using ¦ equivalent circuits. We call ABCDi and ABCDk as the ABCD matrices for each transmission line, and de¯ned 65 Figure 7.7: Cascaded transmission line with a UPFC by ABCDi = 2 64 Ai Bi Ci Di 3 75 and ABCDk = 2 64Ak Bk Ck Dk 3 75 ; where each element is de¯ned by Ai = Di = 1 + YiZi 2 ; Bi = Zi; Ci = Yi(1 + YiZi 4 ); Ak = Dk = 1 + YkZk 2 ; Bk = Zk; Ck = Yk(1 + YkZk 4 ): The ABCD parameters of each transmission line can be obtained after we identify the propagation constant ° and the characteristic impedance Zc. Since we are using IEEE test cases with no knowledge of ° or Zc, we compute these values using the expressions ° = 1 l cosh¡1(1 + Y Z 2 ); (7.23) Zc = Z sinh(°l) ; (7.24) where l is the distance between buses i and k measured in kilometers. Then, assuming the UPFC is installed in x (0 < x < 1), ¦ equivalent circuit values for each section 66 of the transmission line can be found by Zi = Zc sinh(°l ¢ x); Yi = 2 Zi (cosh(°l ¢ x) ¡ 1) ; Zk = Zc sinh(°l ¢ (1 ¡ x)); Yk = 2 Zk (cosh(°l ¢ (1 ¡ x) ¡ 1) : Now, the three cascaded networks are combined to obtain 2 64 ~V i ~Ii 3 75 = ABCDi ¢ ABCDU ¢ ABCDk 2 64 ~V k ¡~Ik 3 75 ; (7.25) = 2 64 Aik Bik Cik Dik 3 75 2 64 ~V k ¡~Ik 3 75 ; where ABCDU is given in (7.12). So Aik = ¹ TAiAk + j ¹ TBiAk½ + 1 ¹ T¤BiCk; Bik = ¹ TAiBk + j ¹ TBiBk½ + 1 ¹ T¤BiDk; Cik = ¹ TCiAk + j ¹ TDiAk½ + 1 ¹ T¤DiCk; Dik = ¹ TCiBk + j ¹ TDiBk½ + 1 ¹ T¤DiDk: By rearranging (7.25) and solving for Ii and Ik, we have 2 64 ~Ii ~Ik 3 75 = ¹Y busik 2 64 ~V i ~V k 3 75 ; (7.26) where ¹Y busik = 2 64 Dik Bik Cik ¡ AikDik Bik ¡ 1 Bik Aik Bik 3 75 : In general, det 0 B@ 2 64 Aik Bik Cik Dik 3 75 1 CA = 1\2Á; (7.27) 67 If Á = 0, that is, complex ¹ T is real, this determinant is one at angle zero, and complex ¹Y busik becomes symmetrical. Note that equation (7.25) represents a bilateral two port network only if ¹ T = 1\0. As seen in (7.26), since the UPFC is embedded in the ¹Y bus matrix, the size of the ¹Y bus matrix is not changed, so UPFC sensitivity analysis can be performed using this ideal transformer model. 68 CHAPTER 8 GAS AND ELECTRICITY OPTIMAL POWER FLOW The purpose of this chapter is to construct a mathematical formulation to solve natural gas and electricity optimal power °ow problems. To do this we present a simple combined gas and electric network, and synthesize an optimal gas °ow (OGF) and an optimal electric power °ow (OPF) into a single natural gas and electricity optimal power °ow (GEOPF) problem. Though this integration of gas and electric system is based upon deterministic prices of gas in source nodes, we will certainly need this GEOPF algorithm for stochastic cases for future studies where the price will be a stochastic variable. For our purposes, gas and electricity storage are neglected. 8.1 Gas and Electric Combined Network For an integrated gas and electric network, even though the gas network and electric network are physically overlapped, we represent the two systems separately. One example is shown in Figure 8.1. Electric generator buses which are coincident with any gas nodes can be used to integrate the gas and electric networks. The generators in the combined nodes are assumed to be driven by gaspowered turbines. 69 Figure 8.1: Combined natural gas and electricity network 8.2 Gas and Electricity Optimal Power Flow The mathematical formulation of the GEOPF can be expressed as min Y C(Y ) ¡ B(Y ) (8.1) subject to: hi(Y ) = 0; i = 1; : : : ; n (8.2) gj(Y ) · 0; j = 1; : : : ;m (8.3) where ² C(Y ) is total cost of system operation, and B(Y ) is the total bene¯t to society and is treated as a negative cost. SW = B(Y )¡C(Y ) is the social welfare, which is the total bene¯t B(Y ) to society, less the cost of combined system operation, and thus is the negative of the net cost of system operation to society. So minimizing net cost C(Y ) ¡ B(Y ) is the same as maximizing social welfare. 70 ² Y is a vector of decision variables (i.e. the voltage magnitude and angle at each bus, real and reactive power generations, real and reactive power consumptions, nodal pressures, °ow rates through pipelines, °ow rates through compressors, gas consumptions by gas turbines, etc.). ² h is a set of equality constraints, such as the real and reactive power balance at each bus, pipeline branch °ow rates, the mass °ow balance at each node, etc.. ² g is a set of inequality constraints, such as line °ow limits, voltage limits, gen eration limits, relative boost limits, nodal pressure limits, etc.. Then the Lagrangian for the GEOPF problem can be written L(Y ) = ¡SW(Y ) ¡ Xn i=1 ¸ihi(Y ) + Xm j=1 ¹jgj(Y ): (8.4) Here, the objective of the GEOPF is to maximize the social welfare, while satisfying n equality constraints and m inequality constraints. 8.2.1 Cost, Bene¯t, and Social Welfare In the \Poolco" environment, central dispatch sets the price of gas or electricity in the system to optimize a networkwide nondiscriminatory service, while consumers and producers contract competitively [8, 1]. The type of problem that the system operator wishes to solve is an optimal power °ow (OPF) problem in electricity and an optimal gas °ow (OGF) problem in natural gas, but doing so simultaneously in an integrated manner. The objective function to be maximized that we will use in the GEOPF problem is social welfare [8]. To understand the concept of the social welfare, we introduce the fact that there are two types of market participants, producers and consumers. A producer is a supplier of some type of goods. Of course, there is a cost to produce the good. Let's say we know the cost of producing the good as a function of the amount 71 of the good that is produced. In our case, the electricity producer is supplying electric power, so its cost is a function of the amount of real power PG it produces. The total cost of real power generation is then (assuming quadratic cost curves) CE = XNG i=1 i=2G CEi(PGi); (8.5) = XNG i=1 i=2G £ ®Gi + ¯GiPGi + °GiP2G i ¤ ; where NG = total number of real power generators; ®Gi ; ¯Gi ; °Gi = cost coe±cients of real power generator i; G = electric nodes with gas¯red generators: For the combined node, it is important to note that the electric generation cost is not included in the equation (8.5). This is because the combined nodes simply transform gas energy into electric energy. For the gas supplier, its production cost is a function of the amount of gas pro duced. Then, the total production cost becomes CG = XNS i=1 ciwSi106 £ GHV; (8.6) where NS = total number of source nodes; ci = gas production cost at source node i ($/MMBTU), GHV = gas gross heating value (BTU/SCF): A consumer purchases a good because he receives some bene¯t from using the good. Let us say that this bene¯t that he/she receives can be measured in dollars and is a quadratic function of the amount of the goods purchased. For the electric 72 power consumer, the bene¯t BE in dollars is a function of the amount of real power PL consumed. The total bene¯t received from the consumption of real power (assuming quadratic bene¯ts) by loads is obtained by BE = NXEL i=1 BEi(PLi); (8.7) = NXEL i=1 £ ¯ELiPLi + °ELiP2 Li ¤ ; where NEL = total number of electric consumers; ¯ELi ; °ELi = bene¯t coe±cients of real power consumer i: For a gas consumer, his bene¯t BG is a function of the amount of gas consumed. Then, the total bene¯t received from the consumption of gas is expressed as BG = NXGL i=1 i=2G BGi(wLi); (8.8) = NXGL i=1 i=2G £ ¯GLiwLi + °GLiw2L i ¤ ; where NGL = total number of gas consumers; ¯GLi ; °GLi = bene¯t coe±cients of gas consumer i: Note again that the gas consumer's bene¯t at the combined node is not included since this node is simply converting gas energy to electric energy. Now we de¯ne the meaning of social welfare. It is the total bene¯t received by society due to the production and consumption of the good. In the case of a combined gas and electric network, the social welfare is de¯ned as SW = BE + BG ¡ CE ¡ CG (8.9) = NXEL i=1 BEi(PLi) + NXGL i=1 i=2G BGi(wLi) ¡ XNG i=1 i=2G CEi(PGi) ¡ XNS i=1 CGi(wsi); 73 and it is the objective function of the GEOPF which we want to maximize, or the negative of that which we wish to minimize. 8.2.2 Constraints The GEOPF is a constrained optimization problem, which must satisfy various equal ity and inequality constraints at the optimum. ² Equality Constraints { Electric Network Kircho®'s laws must be satis¯ed at each electric bus, resulting in the stan dard power °ow (power balance) equations at each bus PGi ¡ PLi = Pi(V; µ); i = 1; : : : ;NB (8.10) QGi ¡ QLi = Qi(V; µ); i = 1; : : : ;NB (8.11) The above two equations are conventional power °ow equations, and tell that the real and reactive power injection at each bus is a function of the system bus voltage magnitudes and phase angles, or Pi(V; µ) and Qi(V; µ). In our notation, V is the vector of system bus voltage magnitudes, and µ is the vector of bus voltage phase angles. The amount of reactive power consumed by the load, QL, is not considered a variable since it will be assumed that the load has a ¯xed power factor where the power factor is de¯ned as p:f: = p PL P2 L + Q2 L : (8.12) This assumption is not critical to our formulation, and is made for simplic ity and concreteness only. Then QL is a linear function of the real power 74 consumed, or QL(PL) = PL p:f: q 1 ¡ p:f:2 : (8.13) The real and reactive injection at bus i can be obtained by Si(V; µ) = ~V i~I¤ i ; i = 1; : : : ;NB (8.14) Pi(V; µ) = Re( ¹ Si); Qi(V; µ) = Im( ¹ Si); where the vector of current injections at each bus is obtained by ~I = ¹Y bus ¢ ~V : { Gas Network There are two sets of steadystate network °ow equations, which must be satis¯ed as equality constraints. The ¯rst equality constraint says that the amount of gas at node i transported from and to other nodes must be the same as that of gas injected at that node. This is known as a mass °ow balance equation. The set of mass °ow equations that describe a gas network is given by (A + U)f + w ¡ T¿ = 0; (8.15) where units of f, w and ¿ are SCF/hr, hence equivalent to conversation of mass. Note that the °ow rates to the compressor inlet and from the compressor outlet are di®erent if the compressor is driven by a gas¯red turbine. ¿k is the amount of gas supplied to gas¯red turbine k measured in SCF/hr. For simplicity, we assume ¿k = ®Tk + ¯TkHk + °TkH2 k ; k = 1; : : : ;NC (8.16) where Hk is the horsepower required for compressor k, and de¯ned by Hk = Bk fk "µ ¼j ¼i ¶Zki(®¡1 ® ) ¡ 1 # : k = 1; : : : ;NC (8.17) 75 The other equality constraint which is described by equation (5.3) is the pipeline gas °ow rate through each pipeline, and it is stated by ¼2 i ¡ ¼2 j = 1 M2 k Sij f2 kij ; i; j = 1; : : : ;NN (8.18) Note that while the nodal °ow equation (8.15) is linear, the pipeline branch °ow equation (8.18) is quadratic, and we need to take into account the fact that the °ow direction must be explicitly accounted for (by Sij . { Gas and Electricity Combined Node This is the location of an actual or planned gas¯red generation facility, which would be located where a gas node and an electric power bus have a common location. The gas node may have other loads or input injections, as may the electrical node the generator is connected to. Let us suppose that the electrical bus in question is i, which is fed gas from gas node ji. Then wji = wji,other (8.19) ¡ ¡ aGi + bGiPGi + cGiP2G i ¢ 1 GHV ; i 2 G where GHV = gas gross heating value (BTU/SCF); aGi ; bGi ; cGi = gas fuel rate coe±cients at node i, G = electric nodes with gas¯red generators. ² Inequality Constraints { Electric Network The system operating constraints include maximum and minimum limits 76 on the voltage magnitude at each bus and the magnitude of the line current °owing on each line. These can be written as Vmini · Vi · Vmaxi ; i = 1; : : : ;NB (8.20) for the voltage magnitude at bus i, and I2 ikmax ¸ ¯¯¯ ~Iik ¯¯¯ 2 = ~Iik ¢ ~I¤ ik; i; k = 1; : : : ;NB (8.21) where ~Iik = ³ ~V i ¡ ~V k ´ ¹yik + ~V i¹yshuntik for the magnitude of the line current from bus i to bus k. Equation (8.21) is called a transmission thermal limit. Since the voltage magnitude at each bus is changed, we use the line current instead of the apparent power °owing on each line for the thermal limit, although an apparent power constraint could be used instead. Also, if voltage drop or the steady state stability limits real power °ow, these constraints could be used instead. There are also limits on the amount of real and reactive power produced by each generator. There may also be a limit on the amount of real power consumed by each load. This gives us an additional set of inequality con straints. PGmini · PGi · PGmaxi ; i 2 NPG (8.22) QGmini · QGi · QGmaxi ; i 2 NQG (8.23) PLmini · PLi · PLmaxi ; i 2 NEL (8.24) The tap magnitude and the tap angle of a regulating transformer, the volt age magnitude and the angle drop along a transmission line, MW inter change transactions, shunt reactors or capacitors, etc. can be also included as inequality constraints if required [11]. 77 { Gas Network The system operating constraints include maximum and minimum limits on the pressure at each node and maximum relative boost limits for each compressor. These can be written as ¼mink · ¼k · ¼maxk ; k = 1; : : : ;NN (8.25) for the pressure ¼i at node i, and 1 · µ ¼j ¼i ¶ k · Rmaxk ; k = 1; : : : ;NC (8.26) for the relative boost rate between the inlet pressure ¼i and the outlet pressure ¼j at compressor station k. There are also limits on the amount of gas supplied at each source node. There may also be limits on the amount of gas consumed by gas consumers. This gives us an additional set of inequality constraints. wminSi · wSi · wmaxSi ; i = 1; : : : ;NS (8.27) wminLi · wLi · wmaxLi ; i = 1; : : : ;NGL (8.28) 78 CHAPTER 9 OPTIMAL LOCATION OF A UPFC IN A POWER SYSTEM This chapter presents a screening technique for greatly reducing the computation involved in determining the optimal location and estimation of the generation cost saving when a UPFC is installed and operated optimally. An technique to obtain the ¯rst and secondorder sensitivities of the generation cost with respect to UPFC con trol parameters is described. The sensitivity analysis attempts to estimate the UPFC value without having to run the OPF with the UPFC several times. This sensitivity technique can be used to optimally locate a UPFC in a large power system by ignoring the transmission lines with low marginal value (MV) and low estimated incremental value (IV), and running a full OPF only for the lines with higher estimated IV to obtain actual IV. Thus, this technique can greatly reduce the computational burden of determining the optimal location of the UPFC in a large power system. 9.1 Optimal Power Flow with UPFC Suppose that a UPFC is installed in transmission line ik. The mathematical formu lation of the OPF with the UPFC can be expressed as min y;xik C(y; xik) (9.1) 79 subject to: hi(y; xik) = 0; i = 1; : : : ; n (9.2) gj(y; xik) · 0; j = 1; : : : ;m (9.3) where ² C(y; xik) is the total generation cost. ² y is a vector of decision variables. ² xik = [ Tik Áik ½ik ]T is a vector of the UPFC control variables in line ik. ² fhi : i = 1; : : : ; ng is the set of equality constraint functions. ² fgj : j = 1; : : : ;mg is the set of inequality constraint functions. We use the UPFC ideal transformer model to construct the equations for OPF with a UPFC. It is important to note that the number of equality constraints is the same as that of the base case OPF with no UPFC. This is because the UPFC control variables do not depend on UPFC input and output voltages and currents, and the UPFC model is embedded in the ¹Y bus matrix, and because we ignore UPFC operation limits. Now, let us construct the Lagrange for the OPF problem as Lo(y; ¸; ¹; xik) = C(y; xik) + Xn i=1 ¸ihi(y; xik) + Xm j=1 ¹jgj(y; xik); (9.4) where ¸i and ¹j are the Lagrange multipliers for the equality and inequality con straints, respectively. To solve the proposed OPF problem with inequality con straints, we use the primaldual interiorpoint method. At the optimum, the last term of (9.4) must satisfy the complementary slackness condition such that ¹jgj = 0 for each j = 1; : : : ;m. Therefore, if an inequality constraint is binding in (9.4), we could treat it as an equality constraint, and we could ignore it if it is not binding. 80 Since we are using interiorpoint methods { not activeset methods { we do not have to distinguish between active and inactive constraints until the OPF problem is solved [20]. This avoids \cycling" behavior in the active set associated with \activeset" methods such as Newton's. Then, to derive the ¯rstorder sensitivities, we rewrite (9.4) as Lo(y; ¸; xik) = C(y; xik) + X j2A ¸jhj(y; xik); (9.5) where A is the set of active constraints, which are known once the basecase OPF is solved. 9.2 FirstOrder Sensitivity Analysis We consider the case where the UPFC is inserted in line ik, and the UPFC ideal transformer model is used for the analysis. The marginal values(MVs) of the UPFC, to be installed in line ik, are simply the amounts by which the total cost of system operation could be changed by allowing a small change of the UPFC control variables in line ik. We can obtain the MVs by assuming that there is a UPFC in line ik, but that the UPFC is not operating. So we add three extra constraints Tik = T; Áik = Á; ½ik = ½ to the original OPF problem, and for simplicity, we denote the constraints as xik = x. Then, the new Lagrangian can be written Lik(y; ¸; xik; ¸x) = C(y; xik) + X j2A ¸jhj(y; xik) + ¸Tx (x ¡ xik); (9.6) where ¸x = [ ¸T ¸Á ¸½ ]T : We de¯ne the function ¸¤ x(x) to be the optimal value of the Lagrange multiplier on the constraint xik = x. Here, we are most interested in ¸¤ x(x = x0), which is associated 81 with the constraints: Tik = 1; Áik = 0; ½ik = 0: This is because the OPF problem when solved with the UPFC control parameters x = x0 yields the same result for y and ¸ as the base case where there is no UPFC in line ik. Using the ¯rstorder conditions for the solution of the OPF problem for xik = x0; @Lik @xik = 0; (9.7) we can solve for ¸¤ x(x0) to obtain ¸¤ x(x0) = " @C(y¤; xik) @xik + X j2A ¸¤ j @hj(y¤; xik) @xik #¯¯¯¯¯ xik=x0 : (9.8) Note that @C(y¤; xik) @xik and ¸¤ j @hj(y¤; xik) @xik ¯¯¯¯ xik=x0 are easy to compute. Equation (9.8) indicates that the marginal value ¸¤ x(x0) can be determined once we know y¤ and ¸¤, which are obtained from the base case OPF with no UPFC. Thus, if we know y¤ and ¸¤, we can obtain the ¯rstorder sensitivities of cost with respect to UPFC control variables x for each possible transmission line by solving only the basecase OPF. Now, we will give the desired interpretation for the Lagrange multiplier ¸x by essentially rederiving the envelope theorem for our case. The Lagrangian in (9.6) can be rewritten as Lik(z; ¸z) = C(z) + ¸Tz hz(z); (9.9) where z = [ y xik ]T , ¸z = [ ¸ ¸x ]T and hz = [ h x ¡ xik ]T . Let the solution of the OPF with the UPFC be denoted by C¤(x), and let the optimal values of z and ¸z be denoted by z¤(x) and ¸¤z (x), respectively. Then C¤(x) = Lik(z¤(x); ¸¤ z(x); x) . (9.10) 82 We di®erentiate C¤(x) with respect to x to obtain d dx C¤(x) = d dx Lik(z¤(x); ¸¤ z(x); x) (9.11) = d dx £ C(z¤(x); x) + ¸¤T z (x)hz(z¤(x); x) ¤ = · @ @z C + ¸¤T z (x) @ @z hz ¸  {z } 0 dz¤(x) dx + £ hTz (z¤(x); x) ¤  {z } 0 d¸¤z (x) dx + @C @x (z¤(x); x) + ¸¤T z (x) @hz @x (z¤(x); x) = @C @x (z¤(x); x) + ¸¤T z (x) @hz @x (z¤(x); x); where the two terms in square brackets are equal to zero by virtue of the KKT ¯rst order optimality conditions. Equation (9.11) indicates that the change in the system generation cost results from the change in C due to the change in x and the change in hz due to the change in x, as the envelope theorem says. The interpretation of ¸x can be obtained by considering the constraint (x ¡ xik) and writing C(z¤(x)) = C(z¤(x); x) (9.12) and h(z¤(x); x) = h (z¤(x)) + [ 0 ¢ ¢ ¢ 0 x ¡ xik ]T : (9.13) Then the envelope theorem [8, 28] says that d dx C¤(x) = ¸¤T x (x): (9.14) For any value of the parameter x, the interpretation of the multiplier ¸¤ x(x) is the marginal change in the total generation cost as the constraint xik = x is changed slightly. If such a change is allowed, it would be made to improve the objective function. The sign of ¸¤ x determines the desired direction of the change in components of xik. Therefore, only the absolute value of the multiplier matters. As we vary x from x0 to the set of optimal values x¤ ik, the multiplier will vary from the marginal 83 value ¸¤ x(x0) to the value ¸¤ x(x¤ ik) = 0. This is because, if we constrain xik to be equal to its unconstrained optimal value, the multiplier for this constraint must be zero. Since the UPFC model is embedded in the ¹Y bus matrix, as explained in (7.26), the ¯rstorder UPFC sensitivity analysis is only associated with complex power in jections at buses i and k, and the thermal limit of transmission line ik if it is binding. If transmission °ow is limited by steadystate stability, such constraints (jPikj · Pmax; or jµi ¡ µkj · µmax) can be included as well. 9.3 SecondOrder Sensitivity Analysis The secondorder sensitivity we wish to obtain is Hx = d2C¤(xik) dx2 ¯¯¯¯ xik=xo = d¸¤ x(xik) dx ¯¯¯¯ xik=xo : (9.15) As in the ¯rst order sensitivity analysis, we arti¯cially insert the UPFC between bus i and bus k, but then constrain xik = x, where we are primarily interested in the case where x = x0, since it represents the base case OPF with no UPFC installed. Rather than including the additional control variable xik in the vector of decision variables y, and the additional multipliers ¸x in the vector of multipliers ¸, we keep these two variables separate for clarity. By the ¯rstorder conditions for optimality, the Lagrange function (9.6) must satisfy the following condition: !(u; x) = ruLik (9.16) = 2 6666666664 ryLik r¸Lik rxikLik r¸xLik 3 7777777775 = 2 6666666664 ryLik h(y; xik) rxikLik (x ¡ xik) 3 7777777775 = 0; 84 where u = 2 66666664 y ¸ xik ¸x 3 77777775 : Note that x occurs only in the term ¸Tx (x ¡ xik) of the Lagrangian in (9.6), so that mixed partials involving y and x, ¸ and x, and xik and x are zero. Thus we have @! @x = 2 66666664 0 0 0 +I 3 77777775 ; (9.17) where I is an identity matrix with the size of xik. As the vector x of the UPFC parameters changes, so will the optimal value of u. Let us denote the optimal value by u¤(x). Then for each x, it must satisfy the condition (9.16): !(u¤(x); x) = 0; 8x (9.18) Since optimality is maintained as we change x, it must satisfy d dx !(u¤(x); x) = 0: (9.19) If we expand (9.19) by the chain rule, then we have @! @u ¢ du¤(x) dx + @! @x = 0. (9.20) The ¯rst term of (9.20) can be expressed as @! @u = Wik(u¤(x); x) (9.21) = 2 6666666664 Hy JT y ryrTx ik Lik 0 Jy 0 Jxik 0 rxikrTy Lik JT xik Hxik ¡I 0 0 ¡I 0 3 7777777775 ; 85 where Hy and Hxik are the Hessians of the Lagrangian, and Jy and Jxik are the Jacobians of constraints in the active set A with respect to y and xik, respectively. Using Wik(u¤(x); x), we can rewrite (9.20) in simple form as Wik du¤(x) dx = ¡ @! @x : (9.22) Then, (9.22) becomes Wik ¢ 2 66666666664 dy¤(x) dx d¸¤(x) dx dx¤ ik(x) dx d¸¤ x(x) dx 3 77777777775 = 2 666666666664 0 0 0 ¡I 3 777777777775 : (9.23) In equation (9.23), we can see that dx¤ ik=dx = I, an obvious result since we have constrained xik = x. The ¯rst two equations can be combined and rearranged to obtain 2 64 dy¤(x) dx d¸¤(x) dx 3 75 = ¡W¡1 ¢ 2 664 ryrTx ik Lik Jxik 3 775 ; (9.24) where W = 2 64 Hy JT y Jy 0 3 75 : Once this is found, the third equation can be rearranged to yield d¸¤ x(x) dx ¯¯¯¯ x=x0 = · rxikrTy Lik JT xik ¸ ¢ 2 64 dy¤(x) dx d¸¤(x) dx 3 75 ¯¯¯¯¯¯¯ x=x0 + Hxik : (9.25) The matrix W is not a®ected by the UPFC since we set x = x0. Thus, to calculate d¸¤ x=dx for x = x0, we already have available W matrix in factored form, so only single forward substitution and back substitution are required in (9.24). 86 9.4 Estimation of Incremental Value using First and SecondOrder Sensitivities Now, we will estimate the incremental value of the UPFC using the ¯rst and second order sensitivities. We approximate ¸x(x) as a linear function of x. Let us expand C(z¤(x); x) at the optimum with respect to x. Then, we have C(z¤(x); x) = C(x) (9.26) = C(x0) + ¸¤T x (x0)¢x + 1 2 ¢xTHx(x0)¢x + H:O:T; where ¢x = x ¡ x0: Neglecting the high order terms in (9.26) and di®erentiating with respect to x yield rxC(x) »= ¸¤ x(x0) + Hx(x0)¢x: (9.27) Since rxC(x) = 0 at the optimum, the estimated optimal value of x is obtained as ^x¤ = ¡H¡1 x (x0)¸¤ x(x0) + x0: (9.28) The estimated incremental value (cIV ) at the estimated optimum ^x¤ is then calculated as cIV = C(x0) ¡ C(^x¤); (9.29) = ¡¸¤T x (x0) (^x¤ ¡ x0) ¡ 1 2 (^x¤ ¡ x0)T Hx(x0) (^x¤ ¡ x0) ; = ¸¤T x (x0)H¡1 x (x0)¸¤ x(x0) ¡ 1 2 ¸¤T x (x0)H¡T x (x0)Hx(x0)H¡1 x (x0)¸¤ x(x0); = 1 2 ¸¤T x (x0)H¡1 x (x0)¸¤ x(x0) ¸ 0: Although IV is by de¯nition nonnegative, nonnegativity of cIV is guaranteed only Hx is positive semide¯nite (which is not guaranteed). The accuracy of cIV depends also on the fact that A does not change in (9.5) since this would change (9.16), the fundamental equation for the sensitivity analysis. 87 CHAPTER 10 RESULTS 10.1 Result of Sensitivity Analysis The proposed sensitivity methods are tested on a 5bus system derived from the IEEE 14bus system, and IEEE 14 and 30bus systems to establish their e®ectiveness. These systems have 7, 20, 41 transmission lines, respectively. Figure 10.1 shows the 5bus system [15]. The system consists of two generators at buses 1 & 2 and one synchronous condenser at bus 3. We assume that generator 2 has higher generation cost than generator 1, and UPFCs are installed in the middle of each transmission line. Loads are assigned such that the current °ow constraint in line 1 is binding (we assume that thermal constraints limit line °ow for this example). The marginal values (j¸T j; j¸Áj and j¸½j), estimated incremental values (cIV ) and Figure 10.1: Diagram of 5bus subset of IEEE 14bus system [15] 88 1 2 3 4 5 6 7 0 2 4 6 8 10 12 Line number MV / Actual IV lT/15 lf/15 lr/5 Estimated IV Actual IV ($/hr) Figure 10.2: Normalized marginal and incremental values for 5bus system. actual incremental values (IV ) for the 5bus system are shown in Figure 10.2. It shows that the lines with high MVs usually produce high estimated IVs. The UPFC location in line 1 produces high MVs for j¸T j and j¸Áj, and the highest estimated IV. In general, it is more economical to locate the UPFC in heavilyloaded high voltage PG1 PG2 Ploss Est. IV Actual IV UPFC Location MW MW MW $/hr $/hr Without UPFC 198.18 71.60 10.51 0 0 Line 1 219.31 51.18 11.64 11.85 6.90 Line 2 224.56 46.19 11.90 7.69 6.38 Line 3 196.85 72.77 10.77 1.13 0.73 Line 4 202.45 67.40 10.99 1.97 1.83 Line 5 204.52 65.44 11.11 2.17 2.07 Line 6 196.85 72.73 10.73 0.75 0.67 Line 7 202.59 67.22 10.95 0.78 0.80 Table 10.1: Real power generation, line loss and total generation cost for 5bus system 89 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 10 12 Line number MV lT/15 lf/15 lr/5 0 2 4 6 8 10 12 14 16 18 20 0 2 4 6 8 Line number IV Estimated IV ($/hr) Actual IV ($/hr) Figure 10.3: Marginal and incremental values for IEEE 14bus system. line 1 since it allows more power to be transmitted in underutilized line 2 while preventing line 1 from overloading. Eventually, no further savings due to UPFC operation can be achieved because of the voltage constraint at bus2. Also, reasonable agreement can be observed between the estimated IVs and the actual IVs except line 1 and line 2. The di®erence between the estimated IV and the actual IV is attributed to the fact that some inequality constraints become binding so that the UPFC cannot be fully utilized. Thus, it is important that the set of binding constraints does not change as x changes in order to obtain a reliable estimated IV. (However, we cannot actually verify this without many additional OPF solutions). For the 5bus case, Table 10.1 shows real power generations, transmission line losses and incremental values. Since the generation marginal cost, adjusted for marginal losses at bus 1 is lower than that at bus 2, it is pro¯table to obtain more 90 real power from generator 1 as long as its lossadjusted marginal cost stays lower and no operational limits are reached. The test results of the IEEE 14bus system are shown in Figure 10.3. Similar load conditions as in the 5bus system are used such that the constraint on current in line 1 is binding. As before, the lines with higher MVs produce higher estimated and actual IVs. An important fact to note is that higher voltage lines 1 through 7 are the most suitable locations to install the UPFC. This is because higher voltage lines have lower p.u. impedances and therefore, most of the power will be transferred through those lines and some of the lines may reach their maximum transfer capabilities. The simulation result for the IEEE 30bus system is shown in Figure 10.4. Again, lines with low marginal values tend to yield low estimated and actual incremental values, and the most valuable locations tend to be in the high voltage portions of the 0 5 10 15 20 25 30 35 40 0 5 10 15 20 Line number MV lT/15 lf/15 lr/5 0 5 10 15 20 25 30 35 40 0 5 10 15 20 25 Line number IV Estimated IV ($/hr) Actual IV ($/hr) Figure 10.4: Marginal ana incremental values for IEEE 30bus system. 91 system. Figure 10.4 also shows that estimated IVs are quite close to actual IVs. In addition, we see that locations with large incremental values are also those with large marginal values, thus supporting the idea of using only ¯rstorder sensitivity analysis to screen for promising locations for installation of a UPFC. 10.2 Result of GEOPF This section contains case studies that evaluate the capability of the GEOPF in maximizing social welfare associated with electricity generation cost, gas supply cost, and gas and electricity consumer's bene¯t. The coe±cients of quadratic generation cost and bene¯t curves are summarized in Appendix C. Since combined nodes are considered as nodes which simply convert gas energy to electric energy, these nodes are not associated with generation cost and gas consumer's bene¯t. Instead, we will Figure 10.5: Combined gas and electric network 92 use the coe±cients of the fuel rate curve for the gas¯red generator to calculate the amount of gas supplied to the combined node. 10.2.1 Test System I We will ¯rst consider a simple combined network with two real power generators, one of which is a gas¯red generator and the other is a coal¯red generator. Tieline connections in both natural gas and electricity networks are not considered in this test system. The system analyzed is a combined natural gas and electricity network, which consists of a 5bus system and a 15node gas system, as shown in Figure 10.5. The 15node system is a modi¯cation of an example system in [22], plus some made up bene¯t and cost functions. The 5bus electric system obtained by modifying the IEEE 14bus system has two real power generators at buses 1 and 2, one synchronous condenser at bus 3, and four loads at buses 2 through 5, and bus 1 serves as the reference/slack bus in the electric network. The generator at bus 2 is coincident with gas network node 15. We assume that node 1 is the knownpressure node in the gas system. The gas network has ¯fteen nodes consisting of ¯ve loads at nodes 3, 4, 13, 14, and 15, two sources at nodes 1 and 2, and eight compressor inletoutlet nodes 5 through 12. We assume that the compressors in the gas network are driven by gas turbines, and the gas is tapped from the outlet node of the compressor station. Appendix C shows the input data for gas and electricity network represented in Figure 10.5. We assume that generator 1 uses coal with a ¯xed price, and generator 2 uses natural gas with high price variations. To study the impact of the wellhead gas prices at source nodes 1 and 2 to the real power generation at combined node 15, three test conditions with di®erent wellhead gas prices are considered as follows: 93 Case 1: Wellhead gas price at node 1 = $2.07 /(106£ BTU) Wellhead gas price at node 2 = $2.12 /(106£ BTU) Case 2: Wellhead gas price at node 1 = $2.46 /(106£ BTU) Wellhead gas price at node 2 = $2.71 /(106£ BTU) Case 3: (pipeline branch between nodes 12 and 14 removed) Wellhead gas price at node 1 = $2.46 /(106£ BTU) Wellhead gas price at node 2 = $2.71 /(106£ BTU) Wellhead gas prices in cases 2 and 3 are higher than those in case 1. Table 10.2 shows simulation results for the three case studies. In case 1, the gas¯red generator reaches its maximum generation limit due to the low gas price at combined node 15. As noted in Table 10.2, the optimal real power generation PG2 at the combined node is quite sensitive to the wellhead gas prices, and reduced greatly as the wellhead gas prices increase. In case 3, we assume that the pipeline branch between nodes 12 and 14 is unavailable due to a forcedoutage. The gas¯red generator reaches its minimum generation limit, which indicates that contingencies occurring in the gas network may result in the signi¯cant change of electricity generation patterns. Now, let us compare the social welfare of nonintegrated gas and electricity operation with that obtained from the GEOPF for case 2. The social welfare from the GEOPF is $22,405.83/hr, and the gas price at node 15 is $2.82/(106£BTU). For nonintegrated network operation, let's say that there is a broker who purchases gas from node 15, and sells it to the gas¯red generator at bus 2. He 



A 

B 

C 

D 

E 

F 

I 

J 

K 

L 

O 

P 

R 

S 

T 

U 

V 

W 


