Table Of ContentSolution 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
