Solution Manual for Control and Dynamics in Power Systems and Microgrids Lingling Fan January 30, 2017 1 1 Chapter 2 Dynamic Simulation 1. A voltage source is serving an RLC series connected circuit. Let R = 0.01Ω, L = 0.01 H, C = 0.001 F. The compensation degree of the system is X /X , approximately 70.36%. Find its c L current response for a step response of the voltage source, and a sinusoidal 60 Hz input (amplitude 1 V)of the voltage source. Use Laplace transformation to find the current in Laplace domain and current in time domain. Solution: We will first find the transfer function from the source voltage to the current by imple- menting the impedance model. The impedance model of the circuit is R+Ls+ 1 . Therefore, the Cs transfer function is I(s) 1 = . (1) V (s) R+Ls+ 1 s Cs Next, weconsidertwocases. Case1, thesourcevoltageisastepresponseanditsLaplacetransform is 1. Case 2, the source voltage is sinusoidal. s Case 1: V (s) = 1 for a step response. Therefore, we may find the current’s Laplace transform s s as 1 Cs 0.001 100 I(s) = = = (2) sLCs2+RCs+1 10−5s2+10−5s+1 s2+s+105 A few algebraic manipulations will lead to 100 I(s) = . (3) (s+0.5)2+3162 Since the inverse Laplace transform for b is eatsin(bt), we then find i(t) in the following (s−a)2+b2 expression. i(t) = 0.3162e−0.5tsin(316t) A (4) Case 2: v (t) = sin(377t). Its Laplace transform is V (s) = 377 . We may find the current’s s s s2+3772 Laplace transform as 377 0.001s 37700s I(s) = = (5) s2+377210−5s2+10−5s+1 (s2+3772)(s2+s+105) The denominator has four roots: ±j377 and −0.5±j316. In order to find i(t), we need to carry out inverse Laplace transform. The technique used is Partial Fraction Expansion to split up a complicated fraction into forms that are in the Laplace transform table. For the above I(s), its fraction expansion expression is k s+k k s+k 1 2 3 4 I(s) = + (6) s2+3772 s2+s+105 where k s are real numbers. i Tofindk andk , let(5)and(6)bothbemultipliedby(s2+3772). Thenweevaluatedtheright 1 2 side terms of the two equations at ±j377. Note that k3s+k4 (s2+377) will be 0 if it is evaluated s2+s+105 2 at s = ±j377. Therefore the following equations will be found. (cid:12) (cid:12) (cid:12) 37700s (cid:12) k1s+k2(cid:12)(cid:12) = s2+s+105(cid:12)(cid:12) (7a) −j377 −j377 (cid:12) (cid:12) (cid:12) 37700s (cid:12) k1s+k2(cid:12)(cid:12) = s2+s+105(cid:12)(cid:12) (7b) j377 j377 =⇒ k = −0.8948,k = 3.0187 1 2 Similarly, we can find k = 0.8948 and k = −2.1239. 3 4 The inverse Laplace transform is k k i(t) = k cos(377t)+ 2 sin(377t)+k e−0.5tsin(316t)+ 4 e−0.5tcos(316t) 1 3 377 316 (8) = −0.8948cos(377t)+0.008sin(377t)+0.8948e−0.5tsin(316t)−0.0067e−0.5tcos(316t) The above expression can be further arranged as follows. i(t) = 0.8948sin(377t−1.5618)+0.8948e−0.5tsin(316t−0.0075) (9) Remarks: We may examine the above expression and know that the steady-state response of the current is 0.8948sin(377t − 1.5618). This time-domain function can also be obtained by phasor-based calculation. The impedance at 60 Hz or 377 rad/s is Z = R+j377L−j 1 = 1.1175ej1.5618. 377C Therefore the current phasor magnitude should be the voltage source magnitude divided by |Z| = 1.1175 and the phase shift should be −1.5618 radian. The steady-state current should be 1 i(t) = sin(377t−1.5618) = 0.8948sin(377t−1.5618). 1.1175 2. Use MATLAB linear system analysis tools to define a linear system for the above RLC circuit. Treat the voltage source as the input while the current as the output. Give a set of Bode plots of the system by varying R. Notate the plot properly. Use MATLAB function step to examine the dynamic response of the current with a step response of the voltage source. Use Matlab function lsim to examine the dynamic response of the current with a sinusoidal input. Solution: The MATLAB codes are as follows. R = [0.01, 0.05, 0.1] ; L = 0.01; C = 0.001; s = tf(’s’); for i=1: length(R) plant(i) = 1/(R(i) + L*s + 1/(C*s)); bode(plant(i)); hold on; end figure step(plant(1)); grid on; t=0:0.001:1; u = sin(377*t); y = lsim(plant(1), u, t, 0); figure plot(t,y); xlabel(’Time (s)’); grid on; 3 Bode Diagram R = 0.01 30 R = 0.05 Magnitude (dB)1200 R =0.1 0 90 Phase (deg)−44550 −90 102.4 102.5 102.6 Frequency (rad/s) (a) Bode plots for three different Rs. Step Response 0.4 0.3 0.2 mplitude 0.01 A −0.1 −0.2 −0.3 −0.4 0 2 4 6 8 10 12 Time (seconds) (b) Problem 2: step output for a step input. 2 1.5 1 0.5 0 −0.5 −1 −1.5 −2 0 0.2 0.4 0.6 0.8 1 Time (s) (c) Problem 2: lsim output for a sinusoidal input. Figure 1: Problem 2. 4 The generated figures are shown in Figs. 1a, 1b, and 1c. 3. For the above RLC circuit, build a two-order state-space model. The state variables are the current and the voltage across the capacitor. Use MATLAB function ode to simulate the dynamic response of the current for a step input and a sinusoidal input. Solution: The state-space model has been built in (2.19). The MATLAB codes to conduct the two simulation case studies are as follows. R = 0.1 ; L = 0.01; C = 0.001; dxdt_step = @(t, x, R, L, C)[(1-R*x(1)-x(2))/L; x(1)/C] ; dxdt_sin = @(t, x, R, L, C)[(sin(377*t)-R*x(1)-x(2))/L; x(1)/C] ; [T1,y1]= ode23(@(t,x)dxdt_step(t, x, R, L, C), [0 1],[0; 0]); [T2,y2]= ode23(@(t,x)dxdt_sin(t, x, R, L, C), [0 1],[0; 0]); figure plot(T1,y1); grid on; xlabel(’Time (s)’); figure plot(T2,y2); grid on; xlabel(’Time (s)’); 2 5 capacitor voltage capacitor voltage 1.5 1 0 0.5 current 0 −0.5 −5 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 Time (s) Time (s) (a) (b) Figure 2: Problem 3: (a) Dynamic response for a step input; (b) Dynamic response for a sinusoidal input. 2 Chapter 3 Frequency Control 1. Use parameters in Example 1 Fig. 3.29. For a single generator load serving system, derive the linear system model and build the model in MATLAB/Simulink. Find the droop to make ∆ω = −0.2 for ∆P = 0.1. L • Find the bandwidth of the system with only primary frequency control. • Providethedynamicsimulationofthesystemfrequencyduetoastepresponseofloadincrease 0.1. • Specify ∆P to bring ∆ω back to zero. c 5 • Design the secondary frequency control to bring the system frequency back to nominal when load varies. Choose a gain of the integration block to make the frequency return to nominal in less than 100 seconds. Find the value of ∆P at steady-state and see if this value matches c your previous calculation. Provide the dynamic responses of ∆ω and ∆P . c • Increasing the gain of the integration block will lead to instability. Find the marginal value of the gain. Designate an analysis procedure to analytically find the marginal gain. Hint: you can use root loci method. Solution: a) The droop parameter R is not given. The electromechanical dynamics is given as 1 = 10 . This indicates that H = 0.5 and D = 0.1. At steady-state, we know that 2Hs+D1 10s+1 1 ∆ω = −1 ∆P . Apply the condition: ∆P = 0.1, ∆ω = −0.2. We may find D +1/R = 0.5. D1+1/R L L 1 Therefore, R = 1/(0.5−0.1) = 2.5. The bandwidth of the system can be found by examining the closed-loop transfer function. The closed-loop transfer function is ∆ω − 1 G = = s+0.1 ∆P 1+ 1 1 1 L s+0.1Rs+1 s+1 (10) = (s+0.1)(s+1)+1/R s+1 = s2+1.1s+0.5 Since the steady-state gain of G is 2. The bandwidth can be found by examining the Bode plot of √ G/2 and the corresponding frequency when the gain is 1/ 2. The Bode plot of G/2 is shown in Fig. 3. Bode Diagram 10 System: untitled1 Frequency (rad/s): 0.812 B) 0 Magnitude (dB): −3 d e ( ud −10 nit g Ma −20 −30 0 g) e d e ( −45 s a h P −90 10−2 10−1 100 101 Frequency (rad/s) Figure 3: Bode plot of G/2. From the Bode plot, we can find the bandwidth is 0.812 rad/s. b) and c) are examined by simulation. ∆P should match ∆P to bring ∆ω back to zero. c L The Simulink model is shown in Fig. 4. The dynamic simulation results are shown in Fig. 5. The simulation results show that given a 0.1 increase in the load, there is 0.2 pu decrease in the 6 Figure 4: Simulink model for Problem 1. 0.1 ∆ P L 0.05 ∆ P c 0 0 5 10 15 20 0.1 0 ∆ω −0.1 −0.2 −0.3 0 5 10 15 20 Time (s) Figure 5: Dynamic simulation results. At t = 1 s, ∆P changes from 0 to 0.1. At t = 10 s, ∆P L c changes from 0 to 0.1. frequency. Further,thefrequencycanbebroughtbacktonominalbyincreasingthereferencepower to match the load. d) and e) examine secondary frequency control design and identify the limit of the gain of the integral controller. The Simulink model is shown in Fig. 6 and the simulation results are shown in Fig. 7. It can be seen that when k is 0.2, the secondary frequency control can bring the frequency back i to nominal. ∆P will eventually become 0.1 to match the load change. This observation confirms c our computation in a). When k is 0.6, the system loses stability. The magnitude of the oscillation i keeps increasing. We may use root locus method to show the limit of k . The open-loop transfer i 7 Figure 6: Simulink model for a system with secondary frequency control. 0.2 0.15 ∆ Pc ki=0.2 00..23 ∆ P ki=0.6 0.1 ∆ P c L 0.1 ∆ P 0.05 L 0 0 0 5 10 15 20 −0.1 0 5 10 15 20 0.2 0.2 0.1 0.1 ∆ω 0 ∆ω 0 −0.1 −0.1 −0.2 −0.2 0 5 10 15 20 0 5 10 15 20 Time (s) Time (s) Figure 7: Dynamic simulation results for a system with secondary frequency control. function is first obtained by decoupling the secondary frequency control loop and it is: 1 1 1 ∆ω 1 G(cid:48) = = s+1s+0.1 s∆P s1+ 1 1 1 c s+1s+0.1R (11) 1 = s(s2+1.1s+0.5) The root loci for G(cid:48) are shown in Fig. 8. It can be seen that the gain cannot exceed 0.559. 2. UsetheparametersinExample3. InMatlab/Simulink, buildthelinearizedmodelforatwo-area interconnection system. Each area consists of a generator with a load. A load change of 100 MW occurs in Area 1. • DemonstratetheACEsteady-statevalues,steady-statefrequencyinHz,changeintie-lineflow in MW for each area through time-domain simulation for the system with primary frequency control only. Observe if the values match the calculation. • Design secondary frequency control for each area and demonstrate that the frequency and tie-line power flow will go back to nominal. Please present simulation results. Solution: This problem is for students to repeat the example presented in Section 3.5.3. 8 Root Locus 2 1.5 −1s) 1 d n o ec 0.5 0.559 s s ( Axi 0 y ar n −0.5 gi a m I −1 −1.5 −2 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 Real Axis (seconds−1) Figure 8: Root loci for G(cid:48). 3 Chapter 4 Synchronous Generator Models 1. A round-rotor generator (X = 1.0, r = 0.1) is synchronized to a bus whose voltage is 1∠00. At s synchronization i = 1000A (actual). The generator is then adjusted until S = 0.8+j0.6. (S is F G G the power supplied to the generator bus.) a. Find i and the generator efficiency (assuming no generator loss except I2r). F b. With the same i , what is the maximum active power the generator can deliver. F Solution: E is proportional to the rotor excitation current i . At the synchronization condition, a F the generator is sending out 0 power and the stator current is also 0. Therefore, E = V , the a a internal voltage is equal to the terminal voltage. So we find that when E is 1, i is 1000 A. We a F can then look at the second operation condition and use phasor diagram or circuit analysis to find E . Then from E , we find i . a a F For the second operating condition, the generator is sending out complex power S = 0.8+j0.6 G while the terminal bus voltage is 1∠00. Thus the stator current I can be found as a (cid:18)S (cid:19)∗ I = G = 0.8−j0.6 = 1∠−36.870. a V a Based on the Thevenin equivalent circuit, we can find E as a E = V +(r+jX )I . a a s a Thus E = 1.0+(0.1+j1.0)(0.8−j0.6) = 1.836∠23.770. a a) When E = 1, i is 1000 A. Therefore, when E = 1.836, i will be 1836 A. a F a F The generator efficiency can be found by computing the real power loss on the stator. The stator current magnitude is 1. Then the loss is I2r = 0.1. The efficiency is 0.8 = 88.9%. 0.8+0.1 9 b) If i does not change, then E will be kept at 1.86. The angle δ may change to have a F a different output power. The expression of the complex power to the grid is as follows. (cid:18) (cid:19)∗ E −V ∗ a a S = V I = V G a a a r+jX s (12) V E V2 = a ae−j(δ−θz)− a ejθz |Z| |Z| (cid:112) where |Z| = r2+X2 = 1.005 and θ = cos−1 r = 84.290. s z |Z| Take the real part of S , we have G V E V2 P = a a cos(δ−θ )− a cosθ . g z z |Z| |Z| When δ = θ , P will be maximized and z G V E V2 1.0×1.836 1.02 Pmax = a a − a cosθ = − cos84.290 = 1.7279. G |Z| |Z| z 1.005 1.005 2. Considerasalient-polegeneratordeliveringpowerthroughashorttransmissionlinetoaninfinite bus. V = 1∠00, E = 1.4. The active power delivered to the infinite bus is 0.6. We are given the ∞ a generator reactances X = 1.6 and X = 1.0 and the line reactance X = 0.4. Neglect resistances d q L and find E and I . a a Solution: A most straightforward approach for this problem is to have the active power ex- pression in terms of δ and E ready. We know that a (cid:32) (cid:33) E V 1 1 V2 P = a ∞ sinδ+ − ∞ sin(2δ). X(cid:101)d X(cid:101)q X(cid:101)d 2 where X(cid:101)d = Xd+XL, X(cid:101)q = Xq +XL. Therefore, we can find δ from that equation. 1.4×1 (cid:18) 1 1 (cid:19) V2 0.6 = sinδ+ − ∞ sin(2δ) 2.0 1.4 2.0 2 (13) 0.6 = 0.7sinδ+0.1071sin(2δ) We may use MATLAB’s function fzero to find roots for a nonlinear equation 0.7sin(x) + 0.1071sin(2x)−0.6. An initial guess should be given as the second argument of this function. Here we may select 0 as the initial guess. fzero(@(x) 0.7*sin(x)+0.1071*sin(2*x)-0.6, 0) The angle is δ = 0.7812rad = 44.760. Based on the phasor diagram, we may then find the dq-axis currents. (cid:40) Ea = V∞cosδ−X(cid:101)dIad (14) 0 = −V∞sinδ+X(cid:101)qIaq (cid:40) I = cos0.7812−1.4 = −0.3450 =⇒ ad 2 (15) I = sin0.7812 = 0.5030 aq 1.4 10