Synchronous Generator:
![]() |
![]() |
The voltage equation of the Synchronous Generator becomes:
The EMF voltage vector is defined as:
The eletromagnetic torque is:
The torque and rotor speed relationship are derived From the Newton's second law, as: where is the input torque, is the moment of inertia and is a damping coefficient.
The reactive power associated with the EMF is given by:
and the reactive power reactive exchanged with grid at the PoC :
Consider the equivalent circuit illustrated in Fig. 1:
The current and aparent power from circuit shown in Fig. 2 can be found as:
The real power and reactive power , delivered from the voltage source to the grid through the reactance , are:
The eletromagnetic torque and input torque are:
The VSG Control Block diagram is:
The nonlinear dynamic equations that describe the behavior of the angular frequency, power angle and VSG flux are:
where the flux loop constant and moment of inertia are given from:
where and are the time constants associated with the voltage and frequency Droop variables, and , respectively.
In order to simplify the notation let us define some constants, that are:
then:
The next step consists in compute the equilibrium point of the nonlinear dynamic equations:
or:
From (27) and (28) we have:
where can be found from the real positive solution of (29).
The local qualitative behavior of the VSG can be analyzed from the eigenvalues of the Jacobian matrix, that is:
pi = 3.14159265359
V = 220. # Grid Voltage
wo = 2.*pi*60. # Fundamental Frequency
L = 5.7e-4 # Inductance
XL = wo*L # Inductive Reactance
wg = 2.*pi*60. # Grid Frequency
Eo = 220. # Reference Voltage
Po = 33e3 * 0. # Active Power
Qo = 33e3 * -1. # Reactive Power
Dp = (33000/wo)*(1/(0.05*wo)) # Droop Coefficient for Frequency
Dq = (33000/(0.1*220)) # Droop Coefficient for Voltage
print (" XL = ", format(XL, ".2f") , " \u03A9")
print (" Dp = ", format(Dp, ".2f") )
print (" Dq = ", format(Dq, ".2f") )
import math
import numpy as np
Vo = V
a = ((wg*V)/XL)**2
b = (wg*(Dp*(wo-wg) + Po/wo))**2
c = Dq*(V - Eo) - Qo
d = wg**2/XL
Psio = np.sqrt(((a/2)+np.sqrt(a**2-4*a*c*d-4*b*d**2)/2-c*d)/(d**2))
Delto = np.arctan2( (wg*(Dp*(wo-wg)+Po/wo)) , (Dq*(V-Eo)-Qo+((wg**2)*(Psio**2))/XL))
Tf = 0.01 # Frequency Time Constant
Tv = 0.05 # Voltage Time Constant
Ji = Dp*Tf # Moment of Inertia
K = Tv*wg*Dq # Flux Loop Constant
Delt = Delto*180/pi
print (" Psi = ", format(Psio, ".2f") , "T")
print (" Delt = ", format(Delt, ".2f") , "\u00b0")
print (" Ji = ", format(Ji, ".2f") )
print (" K = ", format(K, ".2f") )
The Jacobian matrix is given by:
import numpy as np
from scipy import sqrt, exp, cos, matrix, vstack, hstack, zeros
np.set_printoptions(precision=2)
J = matrix([[-Dp/Ji, -Vo*Psio*math.cos(Delto)/(Ji*XL), -Vo*math.sin(Delto)/(Ji*XL)], [1, 0, 0], [-2*(Psio**2)*wo/(K*XL)+Vo*Psio*math.cos(Delto)/(K*XL), -Vo*wo*Psio*math.sin(Delto)/(K*XL), -2*Psio*(wo**2)/(K*XL)+Vo*wo*math.cos(Delto)/(K*XL)]])
print (" J = \n", J)
And the eigenvalues are:
import scipy.linalg as la
evals, evecs = la.eig(J)
Eig_V_J = evals
print (" [ \u03BB\u2081 \u03BB\u2082 \u03BB\u2083 ]'= ", np.array2string(Eig_V_J.reshape((-1, 1)), prefix=' [ \u03BB\u2081 \u03BB\u2082 \u03BB\u2083 ] = ') )
In order to characterizes the interaction between modes and state variables, let us compute the participation factors:
evals, evecs = la.eig(J)
T = evals # Vector with the eigenvalues of J
R = evecs # Matrix whose columns are the Right eigenvectors of J
L = np.linalg.inv(R) # Matrix whose rows are the left eigenvecotors of J
print (" R = \n", R)
print (" \n L = \n", L)
print ("\n \n Participation Factors: ")
print (" \n \u03BB\u2081 \u03BB\u2082 \u03BB\u2083 \n ")
print("|L(0,0)*R(0,0)|=", format(np.abs(L[0][0]*R[0][0]), ".2f"), " |L(1,0)*R(0,1)|=", format(np.abs(L[1][0]*R[0][1]), ".2f"), " |L(2,0)*R(0,2)|=", format(np.abs(L[2][0]*R[0][2]), ".2f"), " \u03C9")
print("|L(0,1)*R(1,0)|=", format(np.abs(L[0][1]*R[1][0]), ".2f"), " |L(1,1)*R(1,1)|=", format(np.abs(L[1][1]*R[1][1]), ".2f"), " |L(2,1)*R(1,2)|=", format(np.abs(L[2][1]*R[1][2]), ".2f"), " \u03B4")
print("|L(0,2)*R(2,0)|=", format(np.abs(L[0][2]*R[2][0]), ".2f"), " |L(1,2)*R(2,1)|=", format(np.abs(L[1][2]*R[2][1]), ".2f"), " |L(2,2)*R(2,2)|=", format(np.abs(L[2][2]*R[2][2]), ".2f"), " \u03C8" )
PN = 10
j = np.arange(start=0, stop=11, step=1)
Vj= 220 - 22 + j*(44)/PN
aj = ((wg*Vj)/XL)**2
b = (wg*(Dp*(wo-wg) + Po/wo))**2
cj = Dq*(Vj- Eo) - Qo
d = (wg**2)/XL
Psi1j = np.sqrt(((aj/2)+np.sqrt(aj**2-4*aj*cj*d-4*b*d**2)/2-cj*d)/(d**2))
Delt1j = np.arctan2( (wg*(Dp*(wo-wg)+Po/wo)) , (Dq*(Vj-Eo)-Qo+((wg**2)*(Psi1j**2))/XL))
E1j = Psi1j*wg
P1j = Vj*E1j*np.sin(Delt1j)/XL
Q1j = (E1j**2 - E1j*Vj*np.cos(Delt1j))/XL # Ggap
Q2j = (E1j*Vj*np.cos(Delt1j) - Vj**2)/XL # Qpoc
import matplotlib.pyplot as plt
plt.figure(figsize=(12, 4
))
plt.subplot(121)
plt.plot(Vj, P1j/33000)
plt.xlim([198, 242])
plt.xticks([198, 209, 220, 231, 242],
["0.9", "0.95", "1", "1.05", "1.1"])
plt.ylim([-1.1, 1.1])
plt.yticks([-1,-0.5, 0, 0.5, 1],
["-1", "-0.5", "0.", "0.5", "1"])
plt.grid(True)
plt.ylabel('$\dfrac{P1j}{33000}$', rotation=0)
plt.xlabel('$\dfrac{Vj}{220}$')
plt.subplot(122)
plt.plot(Vj, Q2j/33000, color="blue", label="${Q_{poc}}$")
plt.plot(Vj, Q1j/33000, color="red", label="${Q_{gap}}$")
plt.xlim([198, 242])
plt.xticks([198, 209, 220, 231, 242],
["0.9", "0.95", "1", "1.05", "1.1"])
plt.ylim([-4.1, 1.1])
plt.yticks([-4, -3, -2, -1, 0, 1],
["-4", "-3", "-2", "-1", "0.", "1"])
plt.grid(True)
plt.ylabel('$\dfrac{Q}{33000}$', rotation=0)
plt.xlabel('$\dfrac{Vj}{220}$')
plt.legend(loc="upper right", frameon=False)
plt.subplots_adjust(left=0.125,
bottom=0.1,
right=0.9,
top=0.9,
wspace=0.4,
hspace=0.35)
Qo = 0
PN = 20
j = np.arange(start=0, stop=PN+1, step=1)
DW = 2.*pi*3.
V = 220
wgj = 2.*pi*(60. - 3.) + (j/PN)*2.*DW
L = 5.7e-4
XLj = wgj*L
aj = ((wgj*V)/XLj)**2
bj = (wgj*(Dp*(wo-wgj) + Po/wo))**2
c = Dq*(V- Eo) - Qo
dj = (wgj**2)/XLj
Psi1j = np.real(sqrt(((aj/2)+np.sqrt(aj**2-4*aj*c*dj-4*bj*dj**2)/2-c*dj)/(dj**2)))
Delt1j = np.arctan2( (wgj*(Dp*(wo-wgj)+Po/wo)) , (Dq*(V-Eo)-Qo+((wgj**2)*(Psi1j**2))/XLj))
E1j = Psi1j*wgj
P1j = V*E1j*np.sin(Delt1j)/XLj
Q1j = (E1j**2 - E1j*V*np.cos(Delt1j))/XLj # Ggap
Q2j = (E1j*V*np.cos(Delt1j) - V**2)/XLj # Qpoc
import matplotlib.pyplot as plt
plt.figure(figsize=(12, 4
))
plt.subplot(121)
plt.plot(wgj/(2*pi*60), P1j/33000)
plt.xlim([0.95, 1.05])
plt.xticks([0.95, 0.975, 1, 1.025, 1.05],
["0.95", "0.975", "1", "1.025", "1.05"])
plt.ylim([-1.1, 1.1])
plt.yticks([-1,-0.5, 0, 0.5, 1],
["-1", "-0.5", "0.", "0.5", "1"])
plt.grid(True)
plt.ylabel('$\dfrac{P1j}{33000}$', rotation=0)
plt.xlabel('$\dfrac{\u03C9 rj}{2 \u03C0 60}$')
plt.subplot(122)
plt.plot(wgj/(2*pi*60), Q2j/33000, color="blue", label="${Q_{poc}}$")
plt.plot(wgj/(2*pi*60), Q1j/33000, color="red", label="${Q_{gap}}$")
plt.xlim([0.95, 1.05])
plt.xticks([0.95, 0.975, 1, 1.025, 1.05],
["0.95", "0.975", "1", "1.025", "1.05"])
plt.ylim([-3.1, 1.1])
plt.yticks([-4, -3, -2, -1, 0, 1],
["-4", "-3", "-2", "-1", "0.", "1"])
plt.grid(True)
plt.ylabel('$\dfrac{Q}{33000}$', rotation=0)
plt.xlabel('$\dfrac{\u03C9 rj}{2 \u03C0 60}$')
plt.legend(loc="upper right", frameon=False)
plt.subplots_adjust(left=0.125,
bottom=0.1,
right=0.9,
top=0.9,
wspace=0.4,
hspace=0.35)