ECE 530 – Analysis Techniques for Large-Scale Electrical Systems

1. ECE 530 – Analysis Techniques for Large-Scale Electrical Systems Lecture 21: SVD, Numeric Solution of Differential Equations Prof. Tom Overbye Dept. of Electrical and Computer Engineering University of Illinois at Urbana-Champaign overbye@illinois.edu

2. Announcements • HW 7 is due Tuesday, November 19 • Special guest lecture by Dr. James Weber from PowerWorld Corporation on November 14

3. Example Power System SVD Application Goal • Enhance Real Time Market Monitoring tools to detect Market Power Potential. Approach • Combine Economic and Engineering models to determine when market participants may have an advantageous position in the market. Impact • Improve the efficiency and reliability of the electric power grid Full algorithm in B.C. Lesieutre, K.M. Rogers, T.J. Overbye, A.R. Borden, "A Sensitivity Approach to Detection of Local Market Power Potential," IEEE Trans. on Power Systems, Nov. 2011, pp. 1980-1988

4. Locational Advantage and Market Power Can suppliers exploit locational advantages? If so - can we identify when this is possible? A supplier with the ability to raise prices without changing dispatch could increase revenues with no change in cost.

5. Motivation: Dispatch Sensitivity M is called the Dispatch Sensitivity Matrix Price Perturbation Vector Dispatch Vector Identify sets of non-substitutable generators: They may adjust price without affecting dispatch The price perturbation vector is not unique. For (m-1) constraints there are m independent vectors that satisfy this equation. For example, a vector of all ones, meaning the price is uniformly changed a each generator.

6. SVD and Null Space • We need to determine the null space of the M matrix • Sis a matrix whose columns correspond to each binding constraint in an OPF/SCOPF, and whose rows correspond to the generators. The entries are then just the sensitivities – how a change in a particular generator will affect the particular constraint. • This is a portion of the LP OPF (or SCOPF) tableau • S is “augmented” with a column of all ones • By doing an SVD on S the columns of V corresponding form the orthonormal basis of the null space of M, which is what we need

7. Example System • The lines connecting buses (2, 6) and buses (5, 7) are constrained, resulting in higher LMPs at buses 6 and 7. • Can already see suggestive patterns in the elements of S

8. Algorithm • For the set of binding line constraints, determine the sensitivity of the line flow (or other constraint) with respect to each of the generators of interest • Augment with a column of all ones • Computationally identical to determining the tableau for an LP OPF • Assume there are (m-1) constraints, n generators, with n >> m • Use a Singular Value Decomposition to obtain an orthonormal basis matrix • the SVD computational order is mn2 so this step is quite quick

9. Algorithm • Do cluster analysis to group generators • Quality Threshold (QT) - does not require a number of clusters do not need to be specified, is deterministic, but is order n2 and requires specifying a distance threshold • K-means - much faster but requires user to specify number of clusters a priori, clusters depend on starting point • Coupling index clustering algorithm we are developing, number of clusters does not need to be specified, allows items to be in multiple clusters • Special implementation of single-linkage agglomerative clustering uses kd-trees and graph theory concepts, order nlogn, requires no parameters, provides a hierarchy which may be sliced at a desired granularity and also studied in a research mode • Perform eigenvector analysis and visualize results

10. Cluster by Rows Use rows of B to initially group generators. For the seven-bus example- Any sensible clustering algorithm appears to work well.

11. Refine Grouping Consider the problem to maximize the magnitudes of the entries associated with clustered elements, while minimizing the rest. A small number of non-zero elements * Mostly zero-valued elements Optimize… *Subscript i denotes elements of cluster i and –i denotes all others

12. Refine Grouping This is an eigenvalue problem with the solution that x is the eigenvector associated with the largest eigenvalue of (BiTBi-B-iTB-i). Solve this problem for each cluster. For the seven-bus example- Examine the result. If there are large entries outside the cluster, add these to the “make-large” group, and repeat previous step.

13. Switching Gears • The analysis we've done so far in the class has been associated with power system static analysis • Determining characteristics of the power system quasi-steady state equilibrium • Now we're going to do a brief coverage of techniques for analysis of power system dynamics, with fuller coverage detailed next semester in 576 • Appropriate models depend on time period of interest • Faster dynamics can be represented as algebraic constraints • Slower dynamics can be represented as constants

14. Lightning Propagation Switching Surges Stator Transients and Subsynchronous Resonance Transient Stability Governor and Load Frequency Control Boiler/Long-Term Dynamics 10-7 10-5 10-3 0.1 10 103 105 Time (Seconds) Power System Time Frames Voltage Stability Power Flow Image source: P.W. Sauer, M.A. Pai, Power System Dynamics and Stability, 1997, Fig 1.2, modified

15. Power Grid Disturbance Example Figures show the frequency change as a result of the sudden loss of a large amount of generation in the Southern WECC Time in Seconds Frequency Contour

16. Frequency Response for Gen. Loss • In response to rapid loss of generation, in the initial seconds the system frequency will decrease as energy stored in the rotating masses is transformed into electric energy • Solar PV has no inertia, and for most new wind turbines the inertia is not seen by the system • Within seconds governors respond, increasing power output of controllable generation • Solar PV and wind are usually operated at maximum power so they have no reserves to contribute

17. Solution Considerations • In ECE 530 we introduce several solution methods that are more fully considered in ECE 576 • A wide variety of different solution methods are possible, with different classes of problems (such as power system transient stability) having customized solutions • There is a balance between the problem to be solved and the solution method • Can we bound the dynamics considered, with fast dynamics represented as algebraic constraints, and slow as constants

18. Differential Algebraic Equations • Many problems, including many in the power area, can be formulated as a set of differential, algebraic equations (DAE) of the form • A power example is transient stability, in which f represents (primarily) the generator dynamics, and g (primarily) the bus power balance equations • We'll initially consider the simpler problem of just

19. Ordinary Differential Equations (ODEs) • Assume we have a problem of the form • This is known as an initial value problem, since the initial value of x is given at some time t0 • We need to determine x(t) for future time • Initial value, x0, must be either be given or determined by solving for an equilibrium point, f(x) = 0 • Higher-order systems can be put into this first order form • Except for special cases, such as linear systems, an analytic solution is usually not possible – numerical methods must be used

20. Equilibrium Points • An equilibrium point x* satisfies • An equilibrium point is stable if the response to a small disturbance remains small • This is known as Lyapunov stability • Formally, if for every e > 0, there exists a d = d(e) > 0 such that if x(0) – x* < d, then x(t) – x* < e for t  0 • An equilibrium point has asymptotic stability if there exists a d > 0 such that if x(0) – x* < d, then

21. Power System Application • A typical power system application is to assume the power flow solution represents an equilibrium point • Back solve to determine the initial state variables, x(0) • At some point a contingency occurs, perturbing the state away from the equilibrium point • Time domain simulation is used to determine whether the system returns to the equilibrium point

22. Initial value Problem Examples

23. Numerical Solution Methods • Numerical solution methods do not generate exact solutions; they practically always introduce some error • Methods assume time advances in discrete increments, called a stepsize (or time step), Dt • Speed accuracy tradeoff: a smaller Dt usually gives a better solution, but it takes longer to compute • Numeric roundoff error due to finite computer word size • Key issue is the derivative of x, f(x) depends on x, the value we are trying to determine • A solution exists as long as f(x) is continuously differentiable

24. Numerical Solution Methods • There are a wide variety of different solution approaches, we will only touch on several • One-step methods: require information about solution just at one point, x(t) • Forward Euler • Runge-Kutta • Multi-step methods: make use of information at more than one point, x(t), x(t-Dt), x(t-D2t)… • Adams-Bashforth • Predictor-Corrector Methods: implicit • Backward Euler

25. Error Propagation • At each time step the total round-off error is the sum of the local round-off at time and the propagated error from steps 1, 2 , … , k − 1 • An algorithm with the desirable property that local round-off error decays with increasing number of steps is said to be numerically stable • Otherwise, the algorithm is numerically unstable • Numerically unstable algorithms can nevertheless give quite good performance if appropriate time steps are used • This is particularly true when coupled with algebraic equations

26. Forward Euler’s Method • The simplest technique for numerically integrating such equations is known as the Euler's Method (sometimes the Forward Euler's Method) • Key idea is to approximate • In general, the smaller the Dt, the more accurate the solution, but it also takes more time steps

27. Euler’s Method Algorithm

28. Euler’s Method Example 1

29. Euler’s Method Example 1, cont’d

30. Euler’s Method Example 2

31. Euler's Method Example 2, cont'd

32. Euler's Method Example 2, cont'd Since we know fromthe exactsolution thatx1 is boundedbetween -1 and 1, clearly themethod isnumericallyunstable

33. Euler's Method Example 2, cont'd Below is a comparison of the solution values for x1(t) at time t = 10 seconds

34. Second Order Runge-Kutta Method • Runge-Kutta methods improve on Euler's method by evaluating f(x) at selected points over the time step • Simplest method is the second order method in which • That is, k1 is what we get from Euler's; k2 improves on this by reevaluating at the estimated end of the time step

35. Second Order Runge-Kutta Algorithm t = 0, x(0) = x0, Dt = step size While t tfinal Do k1 =Dt f(x(t)) k2 =Dt f(x(t) + k1) x(t+Dt) = x(t) + ( k1 + k2) t = t + Dt End While

36. RK2 Oscillating Cart • Consider the same example from before the position of a cart attached to a lossless spring. Again, with initial conditions of x1(0) =1 and x2(0) = 0, the analytic solution is x1(t) = cos(t) • With Dt=0.25 at t = 0

37. RK2 Oscillating Cart

38. Comparison • The below table compares the numeric and exact solutions for x1(t) using the RK2 algorithm

39. Comparison of x1(10) for varying Dt • The below table compares the x1(10) values for different values of Dt; recall with Euler's with Dt=0.1 was -1.41 and with 0.01 was -0.8823

40. RK2 Versus Euler's • RK2 requires twice the function evaluations per iteration, but gives much better results • With RK2 the error tends to vary with the cube of the step size, compared with the square of the step size for Euler's • The smaller error allows for larger step sizes compared to Eulers

41. Fourth Order Runge-Kutta • Other Runge-Kutta algorithms are possible, including the fourth order

42. RK4 Oscillating Cart Example • RK4 gives much better results, with error varying with the time step to the fifth power

43. Multistep Methods • Euler's and Runge-Kutta methods are single step approaches, in that they only use information at x(t) to determine its value at the next time step • Multistep methods take advantage of the fact that using we have information about previous time steps x(t-Dt), x(t-2Dt), etc • These methods can be explicit or implicit (dependent on x(t+Dt) values; we'll just consider the explicit Adams-Bashforth approach

44. Multistep Motivation • In determining x(t+Dt) we could use a Taylor series expansion about x(t) (note Euler's is just the first two terms on the right-hand side)