Document Sample
voltage Powered By Docstoc
					                              Summary of Papers

1. P. Sauer and M. Pai, “Power System Steady-State Stability and the Load Flow Jacobian,” IEEE
   Transactions on Power Systems, Vol. 5, No. 4, Nov. 1990
2. V. Ajjarapu and C. Christy, “The Continuation Power Flow: A Tool for Steady-State Voltage
   Stability Analysis,” IEEE Transactions on Power Systems, Vol. 7, No. 1, Feb., 1992.
3. S. Greene, I. Dobson, and F. Alvarado, “Sensitivity of the Loading Margin to Voltage Collapse
   with Respect to Arbitrary Parameters,” IEEE Transactions on Power Systems, Vol. 12, No. 1,
   Feb. 1997, pp. 232-240.
4. S. Greene, I. Dobson, and F. Alvarado, “Contingency Ranking for Voltage Collapse via
   Sensitivities from a Single Nose Curve,” IEEE Transactions on Power Systems, Vol. 14, No. 1,
   Feb. 1999, pp. 262-272.

                     Voltage Security

Voltage security is the ability of the system to maintain
adequate and controllable voltage levels at all system load buses.
The main concern is that voltage levels outside of a specified
range can affect the operation of the customer’s loads.

Voltage security may be divided into two main problems:
1. Low voltage: voltage level is outside of pre-defined range.
2. Voltage instability: an uncontrolled voltage decline.

You should know that
• low voltage does not necessarily imply voltage instability
• no low voltage does not necessarily imply voltage stability
• voltage instability does necessarily imply low voltage
 There have been several individuals that have significantly
 progressed the field of voltage security. These include:

 • Ajjarapu from ISU

 • Van Cutsem: See the book by Van Cutsem and Vournas.

 • Alvarado, Dobson, Canizares, & Greene:

There are a couple other texts that provide good treatments of
the subject:
• Carson Taylor: “Power System Voltage Stability”
• Prabha Kundur: “Power System Stability & Control”

Our treatment of voltage security will proceed as follows:

 • Voltage instability in a simple system
 • Voltage instability in a large system
 • Brief treatment of bifurcation analysis
 • Continuation power flow (path following) methods
 • Sensitivity methods

     Voltage stability in a simple system

 Consider the per-phase equivalent of a very simple three
 phase power system given below:
V1           V2

                         Node 1                           Node 2
                             +      I                     +

                          V1                               V2

                               _                          _

                                   S12          SD=-S12            5
              Z  R  jX  Y  G  jB                             Note B>0

                    S12  P  jQ12

P | V1 |2 G  | V1 || V2 | G cos(1   2 ) | V1 || V2 | B sin(1   2 )

Q12 | V1 |2 B  | V1 || V2 | B cos(1   2 ) | V1 || V2 | G sin(1   2 )

Let G=0. Then….
                P | V1 || V2 | B sin(1   2 )

               Q12 | V1 |2 B  | V1 || V2 | B cos( 1   2 )
Now we can get SD=PD+jQD=-(P21+jQ21) by

• exchanging the 1 and 2 subscripts in the previous equations.
• negating
  PD   P21   | V1 || V2 | B sin( 2  1 )
              | V1 || V2 | B sin(1   2 )

  QD  Q21   | V2 |2 B  | V1 || V2 | B cos( 2  1 )
                 | V2 |2 B  | V1 || V2 | B cos( 1   2 )

Define 12 =1- 2
                        PD | V1 || V2 | B sin 12
                        QD   | V2 |2 B  | V1 || V2 | B cos12   7
Define:  is the power factor angle of the load, i.e.,

                    V2  I
Then we can also express SD as:
 S D  V2 I * | V2 || I | e j
               | V2 || I | (cos  j sin  )
                                        sin 
               | V2 || I | cos (1  j       )
                PD (1  j tan  )                  Note that phi, and
                                                    therefore beta, is
Define =tan. Then                                 positive for lagging,
                                                    negative for leading.

           S D  PD  jQ D  PD (1  j )                             8
   So we have developed the following equations….

                  PD | V1 || V2 | B sin 12
                  QD   | V2 |2 B  | V1 || V2 | B cos12

               S D  PD  jQ D  PD (1  j )

 Equating the expressions for PD and for QD, we have:
                              QD  PD    | V2 |2 B  | V1 || V2 | B cos12
PD | V1 || V2 | B sin 12          PD   | V2 |2 B | V1 || V2 | B cos12

             Square both equations and add them to get…..

 PD 2  ( PD   | V2 |2 B) 2 | V1 |2 | V2 |2 B 2 (sin 2 12  cos 2 12 )
  PD  ( PD   | V2 |2 B) 2 | V1 |2 | V2 |2 B 2

 Manipulation yields:

| V |          2 PD 
                                                         
       2 2                       2         PD
   2                    | V1 |  | V2 |  2 1   2  0

                 B                        B

Note that this is a quadratic in |V2|2. As such, it has the solution:
                                                                 1/ 2
         | V1 |   PD  | V1 | PD  PD
                  2                4
                                                 2 
| V2 | 
                                    | V1 | 
            2      B   4       B  B              
Let’s assume that the sending end voltage is |V1|=1.0 pu
and B=2 pu. Then our previous equation becomes:

             1  PD  1  PD ( PD  2  )1 / 2
   | V2 |2 
                          % pf = 0.97 lagging
                          pdn=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.78];
                          v2n=sqrt((1-beta.*pdn - sqrt(1-pdn.*(pdn+2*beta)))/2);
                          pdp=[0.78 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0];
                          v2p=sqrt((1-beta.*pdp + sqrt(1-pdp.*(pdp+2*beta)))/2);
                          pd1=[pdn pdp];
    You can make          v21=[v2n v2p];
                          % pf = 1.0
    the P-V plot using    pdn=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.99];
                          v2n=sqrt((1-beta.*pdn - sqrt(1-pdn.*(pdn+2*beta)))/2);
    the following         pdp=[0.99 0.9 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0];
                          v2p=sqrt((1-beta.*pdp + sqrt(1-pdp.*(pdp+2*beta)))/2);

    matlab code.          pd2=[pdn pdp];
                          v22=[v2n v2p];
                          % pf = .97 leading
                          pdn=[0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3];
                          v2n=sqrt((1-beta.*pdn - sqrt(1-pdn.*(pdn+2*beta)))/2);
                          pdp=[1.3 1.2 1.1 1.0 0.9 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0];
                          v2p=sqrt((1-beta.*pdp + sqrt(1-pdp.*(pdp+2*beta)))/2);
                          pd3=[pdn pdp];
                          v23=[v2n v2p];
       Plots of the previous equation for different power factors


                  Real power loading, PD                            12
Some comments regarding the PV curves:
• Each curve has a maximum load. This value is typically
   called the maximum system load or the system loadability.
• If the load is increased beyond the loadability, the voltages will
  decline uncontrollably.
• For a value of load below the loadability, there are two
  voltage solutions. The upper one corresponds to one that can be
  reached in practice. The lower one is correct mathematically, but I
  do not know of a way to reach these points in practice.
• In the lagging or unity power factor condition, it is clear that the
  voltage decreases as the load power increases until the loadability.
  In this case, the voltage instability phenomena is detectable, i.e.,
  operator will be aware that voltages are declining before the
  loadability is exceeded.
• In the leading case, one observes that the voltage is flat, or perhaps
   even increasing a little, until just before the loadability. Thus, in
   the leading condition, voltage instability is not very detectable.
   The leading condition occurs during high transfer conditions when the
   load is light or when the load is highly compensated.                   13
                     QV Curves
   We consider our simple (lossless) system again, with the equations
                   PD | V1 || V2 | B sin 12
                   QD   | V2 |2 B  | V1 || V2 | B cos12
   Now, again assume that V1=1.0, and for a given value of PD
   and V2, compute 12 from the first equation, and then Q from the
   second equation. Repeat for various values of V2 to obtain a QV
   curve for the specified real load PD.

                You can make the P-V plot using
b=1.0;          the following matlab code.


                              The curve on the next page illustrates….                   14
       Q-V Curve


                   QD   15
1. Draw the PV-curve for the following cases, and for each, determine the loadability.
     a. B=2, |V1|=1.0, pf=0.97 lagging
     b. B=2, |V1|=1.0, pf=0.95 lagging
     c. B=2, |V1|=1.06, pf=0.97 lagging
     d. B=10, |V1|=1.0, pf=0.97 lagging
   Identify the effect on loadability of power factor, sending-end voltage, and line reactance.
2. Draw the QV-curves for the following cases, and for each, determine the maximum QD.
     a. B=1, |V1|=1.0, PD=0.1
     b. B=1, |V1|=1.0, PD=0.2
     c. B=1, |V1|=1.06, PD=0.1
     d. B=2, |V1|=1.0, PD=0.1
  Identify the effect on maximum QD of real power demand, sending-end voltage, and line

 Some comments regarding the QV Curves
• In practice, these curves may be drawn with a power flow program
      1. modeling at the target bus a synchronous condenser (a
         generator with P=0) having very wide reactive limits
      2. Setting |V| to a desired value
      3. Solving the power flow.
      4. Reading the Q of the generator.
      5. Repeat 2-4 for a range of voltages.
• QV curves have one advantage over PV curves:
  They are easier to obtain if you only have a power flow (standard
     power flows will not solve near or below the “nose” of PV curves
     but they will solve completely around the “nose” of QV curves.)

         Voltage instability in a large system:
         Influential factors:
             • Load modeling
             • Reactive power limits on generators
             • Loss of a circuit
             • Availability of switchable shunt devices
Two important ideas on which understanding of the above
influences rest:
  1. Voltage instability occurs when the reactive power supply
      cannot meet the reactive power demand of the network.
     •   Transmission line loading is too high
     •   Reactive sources (generators) are too far from load centers
     •   Generator terminal voltages are too low.
     •   Insufficient load reactive compensation
  2. Reactive power cannot be moved very far in a network
     (vars do not travel!).
 Implication: The SYSTEM can have a var surplus but experience
 voltage instability if a local area has a var deficiency.             18
  Load modeling
In analyzing voltage instability, it is necessary to consider the network
under various voltage profiles.

Voltage stability depends on the level of current drawn by the loads.

The level of current drawn by the loads can depend on the voltage seen
by the loads.

Therefore, voltage instability analysis requires a model of how the
load responds to load variations.

Thus, load modeling is very influential in voltage instability analysis.

Exponential load model
A typical load model for a load at a bus is the exponential model:
                                          
          V                      V 
   P  P0  
          V               Q  Q0  
                                   V 
           0                      0
where the subscript 0 indicates the initial operating conditions.
The exponents  and  are specific to the type of load, e.g.,
                                                      
Incandescent lamps            1.54                     -
Room air conditioner          0.50                     2.5
Furnace fan                   0.08                     1.6
Battery charger               2.59                     4.06
Electronic compact florescent1.0                       0.40
Conventional florescent       2.07                     3.21
   Polynomial load model
The ZIP model is a special case of the more general exponential model,
given by a sum of 3 exponential models with specified subscripts:
         V 2     V                          V 2    V     
P  P0  p1    p2  p3 
            V                        Q  Q0 q1    q2  q3 
                                                  V 
         0
                    V0    
                                               0
                                                          V0    
    p1  p2  p3  1.0                        q1  q2  q3  1.0
 where again the subscript 0 indicates the initial operating conditions.

 Usually, values p2 and q2 are the largest.

            So this model is composed of three components:

            • constant impedance component (p1, q1)
            • constant current component (p2, q2)
            • constant power component (p3,, q3)
Effect of Load modeling
Understanding the effect of each component on voltage instability
depends on understanding two ideas:

1. Voltage instability is alleviated when the demand reduces. This
   is because I reduces and I2X reactive losses in the circuits reduce.

2. Since voltage instability causes voltage decline, alleviation of
   voltage instability results if demand reduces with voltage decline.
   This gives the key to understanding the effect of load modeling.

       • constant impedance load (p1) is GOOD since demand
         reduces with square of voltage.
       • constant current load (p2) is OK since demand reduces
         with voltage.
       • Constant power load (p3) is BAD since demand does
          not change as voltage declines.                            22
Some considerations in load modeling
The effects of voltage variation on loads, and thus of loads on
voltage instability, cannot be fully captured using exponential or
polynomial load models because of the following three aspects.

• Induction motor stalling/tripping
• Thermostatic load recovery
• Load tap changers

Thermostatic load recovery
Heating load is the most common type of thermostatic load, and it is
one for which we are all quite familiar. Although much heating is
done with natural gas as the primary fuel, heating systems always
contain some electric components as well, e.g., the fans.

When voltage drops, thermostatic loads initially decrease in power
consumption. But after voltages remain low for a few minutes, the
load regulation devices (thermostats) will start the loads or will
maintain them for longer periods so that more of them are on at the
same time. This is referred to as thermostatic load recovery, and it
tends to exacerbate voltage problems at the high voltage level.

Induction motor stalling/tripping
Three phase induction motors comprise a significant portion of
the total load and so its response to voltage variation is important,
especially since it has a rather unique response.

Consider the steady-state induction motor per-phase equivalent
              Za=R1+jX1                 X’2

    V1                Zb=                              R’2+R’2(1-s)/s
                     Rc//jXm                              =R’2 / s

 Induction motor stalling/tripping
 The current is given by:                           Vth
                                  I '2 
                                         Z th  ( R'2 / s)  jX '2
                         Zb                                         Z a Zb
 where       Vth  V1                  and     Z th  Z a // Z b 
                      Z a  Zb                                     Z a  Zb
Under normal conditions, the slip is typically very small. In this case,
R’2/s >> R’2, and I’2 is small.

But as voltage V1 decreases, the electromagnetic torque developed decreases
as well, the motor slows down. Ultimately, the motor may stall. In this case,
s=1, causing R’2/s = R’2. Thus, one sees that the current I’2 is much larger
for stalled conditions than for normal conditions. Because of X1 and X’2 of
the induction motor, the large “stall” current represents a large reactive load.

Large motors have undervoltage tripping to guard against this, but smaller
motors may not.                                                            26
     Tap changers:

Load tap changers (LTC, OLTC, ULTC, TCUL) are transformers that
connect the transmission or subtransmission systems to the distribution
systems. They are typically equipped with regulation capability that
allow them to control the voltage on the low side so that voltage
deviation on the high side is not seen on the low side.
 V1 and t are                  t:1
 given in pu.   HV side   V1            V1/t   LV side

In per unit, we say that the tap is t:1, where
• t may range from 0.85-1.15 pu
• a single step may be about 0.005 pu (5/8%=0.00625 is very common)
• a change of one step typically requires about 5 seconds.
• there is a deadband of 2-3 times the tap step to prevent excessive tap change.
Under low voltage conditions at the high side, the LTC will decrease t
in order to try and increase V1/t.                                    27
 Tap changers:

Thus, as long as the LTC is regulating (not at a limit), a voltage
decline on the high side does not result in voltage decline at the
load, in the steady-state, so that even if the load is constant Z,
it appears to the high side as if it is constant power. So a simple
load model for voltage instability analysis is constant power!

There are two qualifications to using such a simple model:
1. “Fast” voltage dips are seen at the low side (since LTC
action typically requires minutes), and if the dip is low enough,
induction motors may trip, resulting in an immediate decrease in
load power.
2. Once the LTC hits its limit (minimum t), then the low side
voltage begins to decline, and it becomes necessary to model
the load voltage sensitivity.
   Generator capability curve:

                                 Field current limit due to field heating,
                                 enforced by overexcitation limiter on If.
                                             Armature current limit due to
Typical                                      armature heating, enforced by
approximation                                operator control of P and If.
used in power
flow programs.                                  P


                             Limit due to steady-state instability (small
                             internal voltage E gives small |E||V|Bsin), and
                             due to stator end-region heating from induced
                             eddy currents.                                29
      Effect of generator reactive power limits:

1. Voltage instability is typically preceded by generators hitting their upper reactive
    limit, so modeling Qmax is very important to analysis of voltage instability.
2. Most power flow programs represent generator Qmax as fixed. However, this
   is an approximation, and one that should be recognized. In reality, Qmax is not fixed.
   The reactive capability diagram shows quite clearly that Qmax is a function of P and
   becomes more restrictive as P increases. A first-order improvement to fixed Qmax
   is to model Qmax as a function of P.
3. Qmax is set according to the Over-eXcitation Limiter (OXL). The field circuit has a
   rated steady-state field current2 If-max, set by field circuit heating limitations. Since
   heating is proportional to  I f dt , we see that smaller overloads can be tolerated
   for longer times.         overload
                               tim e
   Therefore, most modern OXLs are set with a time-inverse characteristic:
4. As soon as the OXL acts to limit If, then no             2.0
   further increase in reactive power is possible.                      OXL characteristic
   When drawing PV or QV curves, the                     If
   action of a generator hitting Qmax, will              Irated
   manifest itself as a sharp discontinuity                 1.0
                                                                10 Overload time (sec)  30
   in the curve.
Effect of OXL action on PV curve:

                                        One generator hits reactive limit
                                                        No reactive
                                         o              limits modeled


                 Note: Georgia Power Co. models its loadability
                       limit at point x, not point o.

 Loss of a circuit

     Compare reactive losses with and without second circuit
     Assume both circuits have reactance of X.

        I/2                                I

        I/2               P                                    P

Qloss=(I/2)2X+ (I/2)2X=I2X/2                Qloss=I2X

Implication: Loss of a circuit will always increase reactive losses
            in the network. This effect is compounded by the fact
            that losing a circuit also means losing its line charging
            capacitance.                                          32
Kundur, on pp. 979-990, has an excellent example which illustrates
many of the aforementioned effects. The illustration was done using
a long-term time domain simulation program (Eurostag).

Influence of switched shunt capacitors

      I                                  I

                       P                                     P


                           Without                   With capacitor

                                             (demand)            34
 But, shunt compensation has some drawbacks:

 • It produces reactive power in proportion to the square of the
   the voltage, therefore when voltages drop, so does the reactive
   power supplied by the capacitor.

 • It has a maximum compensation level beyond which stable
   operation is not possible (See pg. 972 of Kundur, and next slide).

(A synchronous condenser and an SVC do not have these 2 drawbacks)
  • It results in a flatter PV curve and therefore makes voltage
    instability less detectable. Therefore, as the load grows in areas
    lacking generation, more and more shunt compensation is used to
    keep voltages in normal operating ranges. By so doing, normal
    operating points progressively approach loadability.
    V1=1.0                                         Each QV curve/Capacitor characteristic
                         V2                        intersection shows the operating point. Note
                                                   that for the first three operating points, a
                                PL                 small increase in Q-comp (indicated by
                                QL=0               arrows) results in voltage increase, but for
                                                   the last operating point (950), more Q-comp
S=|V2|2B*Sbase                                     results in a voltage decrease.
with |V2|=1.0            675 Mvar 450 Mvar    300 Mvar
                   950 Mvar


  QV-curves drawn
  using synchronous
  condensor approach.

      1600      1400   1200   1000   800     600     400      200

                          Capacitive Mvars                                                  36
Bifurcation analysis (ref: A. Gaponov-Grekhov, “Nonlinearities in action” and
also Van Cutsem & Vournas, “Voltage stability of electric power systems.”)
A bifurcation, for a dynamic system, is an acquisition of a new
quality by the motion the dynamic system, caused by small changes
in its parameters. A power system that has experienced a bifurcation
will generally have corresponding motion that is undesirable.
Consider representing the dynamics of the power system as:
                             x  F ( x, y , p )
                                                         Eqts. 1
                             0  G ( x, y , p )
A differential-algebraic system (DAS):
Here x represents state variables of the system (e.g., rotor angles, rotor
speed, etc), y represents the algebraic variables (bus voltage magnitudes
& voltage angles), and p represents the real and reactive power injections
at each bus. The function F represents the differential equations for the
generators, and the function G represents the power flow equations. 37
Types of bifurcations

      There are at least two types of bifurcation:
      • Hopf: two eigenvalues become purely imaginary:
      a birth of oscillatory or periodic motion.
      • Saddle node: a disappearance of an equilibrium state.
      The stable operating equilibrium coalesces with an unstable
      equilibrium and disappears. The dynamic consequence of a
      generic saddle node bifurcation is:
              a monotonic decline in system variables.
       So we think it is the saddle node bifurcation that causes
       voltage instability.

  The unreduced Jacobian:
      The Jacobian matrix of eqts. 1 is

                     F X             FY 
                  J                    
                     G X             GY 
and it is referred to as the unreduced Jacobian of the DAS, where

                       x 
                                x 
                        0   J y 
                                                   Eqt. 2

                                

    The reduced Jacobian:
       We may reduce eq. 2 by eliminating the variable y
               x   F X
                                    F Y  x 
                0   G                 y 
                  X               GY   
This means we need to force the top right hand submatrix to 0, which we can
do by multiplying the bottom row by -FYGY-1 and then adding to the top row.

       x    F X  F Y G Y 1 G X 0  x 
       0                                
                                       G Y  y 
         G X
  This results in:
                                    1
                   x  F X  F Y G Y G X  x           
  So that the reduced Jacobian matrix is a Schur’s complement:
                         A  F X  F Y GY G X
 Fact:The conditions for a saddle node bifurcation are
 1. Equilibrium:         x  F ( x, y , p )
                         0  G ( x, y , p )       F X           FY 
 2. Singularity of the unreduced Jacobian J  
                                                  G X           GY 
     det(J)=0.
  Implication: The stability of an equilibrium point of the DAS depends on
  the eigenvalues of the unreduced Jacobian J. The system will experience a
  SNB as parameter p increases when J has a zero eigenvalue.
 Fact: The determinant of a Schur’s complement times the determinant of
 GY gives the determinant of the original matrix: det(J)=det(A)*det(GY)
 if GY is nonsingular.
1. If GY is nonsingular, then singularity of A implies singularity of J
   so that we may analyze eigenvalues of A to ascertain stability.
2. The fact that GY may be nonsingular, yet A singular, means that
   load flow convergence is not a sufficient condition for voltage 41
     Singularity of load flow Jacobian:

So voltage instability analysis using only a load flow Jacobian may
yield optimistic results when compared to results from analysis of A,
that is, stable points may not be really stable.

However, I believe that it is true that unstable points identified using
the load flow Jacobian will be really unstable (Schur’s complement
does not support this, however, because it is only valid if GY is

Note: Sauer and Pai, 1990, provide an in-depth analysis of the relation
between singularity of GY and singularity of J, and shows some special
cases for which singularity of GY implies singularity of J.

    Singularity of load flow Jacobian:

 So we will assume that load flow Jacobian analysis provides an upper
 bound on stability.

 Fact: The bifurcation point of the load flow Jacobian corresponds to
 the “turn-around point” (i.e., the “nose” point) of a P-V or Q-V curve
 drawn using a power flow program.
     This can be proven using an optimization approach.
     See pp. 218-220 of the text by Van Cutsem and Vournas.

We have previously denoted the power flow equations as G(x,y,p)=0,
but now we denote them as G(y,p)=0, without the dependence on the
state variables x (which relate to the machine modeling and include,
minimally,  and  of each machine).                                 43
So we turn our effort to identifying the saddle node bifurcation
(SNB) for the power flow Jacobian matrix.

The Jacobian can reach an SNB in many ways. For example,
• increase the impedance in a key tie line
• increase the generation level at a generator with weak transmission, while
  decreasing generation at all other generators.
• increase the load at a single bus
• increase the load at all buses.                              |V|
In all cases, we are looking for the “nose” point of the                            
V- curve ( is the parameter that is being increased.)
Most applications focus on the last method (increase load at all buses).
Key questions here are:
• “direction” of increase: are bus loads increased proportionally, or in some other way?
• dispatch policy: how do the generators pick up the load increase ?

We will assume proportional load increase with “governor” load flow
(generators pick up in proportion to their rating)
   Define: critical point - the operating conditions, characterized
           by a certain value of , beyond which operation is not

Question 1:
What can cause the critical point to differ    |V|
from the SNB point ?                                            

Question 2:                                              
How can knowledge of the critical point provide a security
Question 3:
Does the P-V curve provide a forecast of the system trajectory ?

Solution approaches to finding *, the value of  corresponding to SNB.

   Approach 1: Search for * using some iterative search procedure.
            1. i=1
            2. Using (i), solve power flow using Newton-Raphson.
                Here, we iteratively solve G(y,p)=0. At each step,
                we must solve for y in the eqt: GY y = p
            3. If solved,
                       (i+1)= (i)+ .
                       go to 2
               else if not solved,
                       *= (i+1)
            4. End

   But big problem: as  gets close to *, GY becomes ill-conditioned.
   This means that step 2 is no longer feasible.                    46
Approach 2: Use the continuation power flow (CPF).

                   Predictor step

                   Corrector step

                                           No.     Select
                      Pass   * ?                continuation

   The predictor step:
   The power flow equations are functions of the bus voltages and
   bus angles and the bus injections:
                         0  G ( y, p )
 Augment the power flow equations so that they are functions of 
 (dependence on p is carried through the dependence on ).
         pp0          0  G ( y,  )
   Now recognize that y                   0  G( ,V ,  )
                                  so that
                            V 
If we want to compute the change in the power flow equations dG
due to small changes in the variables , V, and ,
• that move us closer to the loadability point
• as we move from one solution i to another “close” solution i+1, then
dG= G((i),V(i),(i))- G((i+1),V(i+1),(i+1)) = 0 – 0 = 0          48
                           dG      dG      dG
                      dG     d     dV     d
                           d      dV      d
Here, each set of partial derivatives are evaluated at a particular set of operating
conditions. If the power flow equations are linear with the 3 sets of variables in the
region between the old solution and the (close) new one, the following is satisfied:

                  dG      dG      dG
             dG     d     dV     d  0
                  d      dV      d                                      d 
                                  Eq. 3        G         GV           d V   0
                                                                    G   
                                                                          d 
                                                                          
BUT, we have added one unknown, , to the power flow problem without adding a
corresponding equation, i.e., in G(,V,)=0, there are are N equations but N+1
variables, so that in eq. 3, the matrix [G GV, G], has N rows (the number of eqts
being differentiated) and N+1 columns (the number of variables for which each eqt is
differentiated). So we need another equation in order to solve this. What to do ? 49
The answer to this can be found by identifying how we will be using using the
solution to eqt. 3. Note the solution corresponding to the “new” point is:

               (i 1, p )   (i )   d '         Here the “p” indicates
               (i 1, p )   (i )   
                              V   dV '
                                                        that this is the
              V                                        “predicted” point.
               (i 1, p )   (i )   d ' 
                               
     If we define  to be the “step size,” then we can rewrite this as
    (i 1, p )   (i )    d 
    (i 1, p )   (i )     dV  where  d  '     d 
   V              V     
                                          d V '   d V 
    ( i 1 p )    (i ) 
                               d                  
                           
                                           d ' 
                                                     d 
                                                       
We call the update vector (with the differentials) the
“tangent” vector, denoted by t.
                        d 
                   t   dV 
                        
                        d 
                        

This vector provides the direction to move in order
to find a new solution (i+1,p) from the old one (i).
We can think of this in terms of the following picture…..

      Tangent vector



Note: In specifying a direction using an n-dimensional vector, only n-1 of the
elements are constrained - one element can be chosen to be any value we like.

      For example, consider a 2-dimensional vector….

                                             x2=x1tan(30) so:
                                             - the direction is specified by
                                               selecting x1=1, x2=0.5774,
                  Direction                  - the direction is specified by
                   = 30o
                              x1               selecting x1=0.5, x2=0.2246.

            So we can set one of the tangent vector elements to
            any value we like, then compute the other elements.
            This provides us with our other equation….                     53
    Suppose that we set the k-th parameter in the tangent vector
    to be 1.0. Then our equation given as eq. 3 can be augmented
    to become:
                                               d 
                   G         GV       G     0 
                                            dV   1 
                              ek             d   
                                               

 e k  [0      0     ...        0       1     0    ...    0]
To select , we would have:

  e k  [0     0         ...        0     0    0    ...      1]   54
The parameter for which we select k is called the continuation
parameter, and it can be any load level (or group of load levels),
or it can be a voltage magnitude. Initially, when the solution is
far from the nose, the continuation parameter is typically .
             (i 1, p )   (i )     d 
             (i 1, p )   (i ) 
            V              V    dV 
                                         
             (i 1, p )   (i )     d 
                                     
                (i 1, p )
             y            y         t
                                   (i )

 The parameter  is called the step size, and it can be selected
 using various techniques. The simplest of these is to just
 set it to a constant.
Corrector step
 Note, however, that the predicted point will satisfy the
 power flow equations only if the power flow equations are
 linear, which they are not.

 So our point needs correction. This leads to the corrector step.

 There are two different approaches for performing the
 corrector step.

     Approach a: Perpendicular intersection method.

     Approach b: Parameterization method

        Approach a: perpendicular intersection

        Here, we find the intersection between the power flow
        equations (the PV curve) and a plane that is perpendicular to
        the tangent vector.

                                        |V|                   y(i)       t
 Solve simultaneously,
 for y(i+1)                                                                     y(i+1,p)
                    (i 1)
 0  G( y                    , )                                     y(i+1)

 y   (i 1)
                y (i 1, p )  t  0
The last equation says the inner
(dot) product of 2  vectors is zero.

 Use Newton-Raphson to solve the above (requires only 1-3 iterations since we have
 good starting point). If no convergence, cut step size () by half and repeat.          57
   Approach b: Parameterization

   Here, the corrector step is performed by fixing the continuation
   parameter and then solving the power flow equations.

                       |V|                                       y(i)      t
Solve simultaneously,
for y(i+1)
                                      Vertical corrections
 G ( y  )
         (i 1)                       correspond to a                   y(i+1)
                0                  load-continuation
                                      parameter, horizontal
  yk (i 1)   
                                    to a voltage correction.

Here, yk(i+1) is the continuation parameter; it corresponds to the k-th element in the tangent
vector. The parameter  is the value to which yk is set, which would be the value found in
the predictor step. Again, you can solve this using Newton-Raphson. If no convergence,
cut step size () by half and repeat.                                                            58
Detection of critical point:

  We will know that we have surpassed the critical
  point when the sign of d in the tangent vector
  becomes negative, because it is at this point where
  the loading reaches a maximum point and begins
  to decrease.
           |V|                       increasing

                                     decreasing

                                                           59
    Selection of continuation parameter:

      The continuation parameter is selected from among 
      and the state variables in y according to the one that is
      changing the most with . This will be the parameter that
      has the largest element in the tangent vector.
      • relatively unstressed conditions (far from nose): generally 
      • relatively stressed conditions (close to nose): generally the
        voltage magnitude of the weakest bus, as it changes a great
        deal as  is changed, when we are close to *.

                  (i 1, p )   (i )    d 
                  (i 1, p )   (i )     dV 
Typically, yk
is going to be   V              V     
                  (i 1, p )   (i )    d 
one of these.
                                         
Selection of continuation parameter:

  The continuation parameter is selected from among 
  and the state variables in y according to the one that is
  changing the most with . This will be the parameter that
  has the largest element in the tangent vector.
  • relatively unstressed conditions (far from nose): generally ,
    for parameterization, it looks like:
|V|                            y(i)

                                 y(i+1)              Here,  is fixed.

                                                                        61
   • relatively stressed conditions (close to nose): generally the
     voltage magnitude of the weakest bus. Here, the voltage being
     plotted is chosen as the continuation parameter.

     |V|                           y(i)

                                          y(i+1)              Here, |V| is fixed.

“Essentially, a variable is fixed as a parameter (the voltage), and
the parameter () is treated as a variable. This process of selecting
a variable to fix is sometimes called the parameterization step.”
                -Scott Greene, Ph.D. dissertation, 1998.              62
      A central question:

      How does the continuation technique alleviate the ill-
      conditioning problem experienced by a regular power flow ?

Refer to the solutions procedures for the two corrector approaches.
     Perpendicular interesection             Parameterization
     Solve simultaneously,                   Solve simultaneously,
     for y(i+1)                              for y(i+1)

      0  G( y (i 1) ,  )                  G ( y (i 1)  )
                                                             0
      y   (i 1)
                     y (i 1, p )  t  0    yk (i 1)   
                                                             
 In both cases, we use Newton-Raphson to solve, so we need to obtain the
 Jacobian. But the Jacobian is slightly different than in normal power flow.
The Jacobian of the power flow equations is just Gy, but the
Jacobian of the equations in the two corrector approaches will
have an extra row and column.
                       G y     G xk 
                                    
                       C y
                               C xk 
  Here, C is the additional equation, and xk is the selected
  continuation parameter.

   This addition of a row and column to the Jacobian has the
   effect of improving the conditioning so that the previously
   singular points can in fact be obtained. In other words, the
   additional row and column provides that this Jacobian is
   nonsingular at * where the standard Jacobian is singular.
Known codes for continuation methods:

1. Claudio Canizarres at University of Waterloo: C-code
   See http://www.power.uwaterloo.ca/~claudio/claudio.html
   UWPFLOW is a research tool that has been designed to calculate local bifurcations related to system
   limits or singularities in the system Jacobian. The program also generates a series of output files that
   allow further analyses, such as tangent vectors, left and right eigenvectors at a singular bifurcation
   point, Jacobians, power flow solutions at different loading levels, voltage stability indices, etc
2. I have Matlab code that does it – from Scott Greene.
3. Venkataramana Ajjarapu (ISU): Fortran code

 Calculation of sensitivities for voltage instability analysis

What is a sensitivity ?

       It is the derivative of an equation with respect to a variable.
       It shows how parameter 1 changes with parameter 2.

It is: exact when parameter 2 depends linearly on parameter 1.
It is approximate when parameter 2 depends nonlinearly on parameter 1,

               but it is quite accurate if it is only
               used close to where it is calculated.

   Consider the system characterized by G(y). Then        y   y*
   is the sensitivity of the equation G with respect to y,
   evaluated at y*.

                                       Slope is G/y
                                       evaluated at y*.


                         y y*
It’s usefulness is that once it is calculated, it can be used to
QUICKLY evaluate f(y) from G(y)G(y*)+ (G/y|y*)y,
Consider parameter p: we desire to obtain the sensitivity of
G(y,p) to p. Typical parameters p would be a bus load, a bus
power factor, or a generation level.

Very important to distinguish between
   • voltage sensitivities

   • voltage instability sensitivities

What is the difference between them in terms of
• what they mean ?
• how to compute them ?

               Sensitivities for bus voltage

            These we compute at the current operating condition.

            For a given continuation parameter, they can be obtained
            from the first predictor step in the continuation power flow.
                                 Recall that this provides us with  d 
                                 the tangent vector, given by:  t   dV 
                                                                     
                                                                     d 
                                                                     
                                 The tangent vector is the vector of
                                sensitivities with respect to a small
                                 change in , so the portion of the vector
                                 designated as dV is exactly the voltage
 operating point                 sensitivities.                         69
 Sensitivities for voltage instability
 Here, it is important to realize that the measure of voltage instability,
        the loading margin, depends on an operating condition
              different from the present operating condition.

                                  The implication is that we must look at sensitivities
                                       of the loading margin, not of the voltage.

             Loading margin                     So we want the sensitivities
                                                evaluated at this point, i.e.,
                                                the SNB point.

 operating point                                                                    70
  Derivation of loading margin sensitivities at SNB point.

Let S be the vector of real and reactive load powers,
and k be the direction of load increase.

               S  S0   k
 Also, define L as the loading margin (a scalar), so that
 the load powers resulting in the SNB point are given by:

               S  S 0  Lk
We desire to find the sensitivity of the loading margin L to a
change in the parameter p. We denote this sensitivity by Lp.
Consider the system characterized by
                                                  We want the sensitivity of
                                                  the loading margin to p.
                          G(y,S, p) = 0

Assumption: the system has a SNB at (y*,S*, p*), i.e., :

1. G(y*,S*, p*) = 0 (an equilibrium point)

2. Gy(y*,S*, p*) is singular (zero eigenvalue), and
   w is a left eigenvector of Gy(y*,S*, p*), corresponding
   to the zero eigenvalue so that (by definition of the left
                        wT Gy(y*,S*, p*) =0 wT=0
   Note that Gy(y*,S*, p*), being singular, cannot be inverted, but
   but we can compute it (Gy y*,S*, p*) and its eigenvectors.

3. wT GS(y*,S*, p*)  0
The points (y,S, p) satisfying numbers 1 and 2 correspond to SNB

    and we can obtain a curve of such points by
    varying p about its nominal value p*.

Linearization of this curve about the SNB point results in

                  G y y  G S * S  G p p  0
                         *                            *
  where the notation |* indicates the derivatives are evaluated at the SNB point.
  Pre-multiplication by the left eigenvector w results in:
           w G y y  w G S * S  w G p p  0
              T                  T                    T
                     *                                       *

 By #2 on the previous slide, the first term in the above is zero. So...
                  w G S * S  w G p p  0
                     T                T
                                                                     Eqt. *

 Now recall the relation of the load powers to the loading margin….

                         S  S 0  Lk
                          S  L k
     Substituting this expression for the load powers into eqt. *,

w G  * Lk  w G p p  0  L w G  * k  w G p p
 T                  T                             T                 T
                          *                                                  *

      And the loading margin sensitivity to parameter p is:
                                                          So p may be, for example,
                          L  w G p *
                                                          real power load at a bus (to
                    Lp     T                           detect the most effective load
                       *  p w G S k                      shedding) or reactive power
                                                          at a bus (to determine where
                                                      *                              74
                                                          To site a shunt cap).
                                                                w Gp
 Some comments about computing Lp                     Lp         T
• The left eigenvector w must be computed for the               w GS k
  Jacobian Gy evaluated at the SNB point.

• You only need to compute w and GS once, independent of
  how many sensitivities you need. Methods to compute the left eigenvector
  w include QR or inverse iteration.

• The vector of derivatives with respect to the parameter p, which is Gp, is
  typically sparse. For example, if you want to compute the sensitivity to a
  bus power, then there would be only 1 non-zero entry in Gp.

• The matrix of derivatives with respect to the load powers, GS, using constant
  power load models, is a diagonal matrix with ones in the rows corresponding
  to load buses. This is because a particular load variable would ONLY occur
  in the equation corresponding to the bus where it is located, and for these
  equations, these variables appear linearly with 1 as coefficient.

     Some comments about extensions

• Multiple sensitivities may be computed using Gp (a matrix) instead of Gp (a vector).
  In this case, the result is a vector.
                                                     w Gp

                                           Lp         T
                                                     w GS k
• Getting multiple sensitivities can be especially attractive when we want to find
  the sensitivity to several simultaneous changes. One good example is to find the
  sensitivity to changes in multiple loads.

• A special case of this is to find the sensitivity to changes at ALL loads, which is
  very typical, given a particular loading direction k . Then
                                                        Lall loads *   ki L pi
 • A sensitivity to a line outage may be obtained by letting p contain elements
   corresponding to the outaged line parameters.
   Some comments about extensions
• A sensitivity to a line outage may be obtained by letting p contain elements
  corresponding to the outaged line parameters: R (series conductance), X (series
  reactance), and B (line charging). Then use the multiple parameter approach.
                                                        Zpq=R + jX
                  w Gp

         Lp         T
                   w GS k                   p                                    q

             L  L p * p                                                jB
  • Here, p = [R X B]T.

  • Note that p is NOT SMALL ! Therefore L may have considerable error.
    For that reason, this one needs to be careful about using this approach to
    compute the actual loading margins following contingencies.

  • However, it certainly can be used for RANKING contingencies. One might
    consider having a “quick approximation” and a “long exact” risk calculation.77
 Some comments about alternatives

• Greene, et al., also propose a quadratic sensitivity which requires
  calculation of a second order term Lpp . This is used together
  with the linear sensitivity according to
                    L  L p p  L pp (p ) 2
                            *    2    *
It requires significantly more computation but can provide greater
accuracy over a larger range of p.

• Invariant Subspace Parametic Sensitivity (ISPS) by Ajjarapu.
    – based on differential-algebraic model
    – provides sensitivities at ANY point on the P-V curve

Shared By: