MODULE 5ΒΆ

Virtual Synchronous Generator - VSGΒΆ

Feedback of the Reactive Power at the GapΒΆ

Synchronous Generator:

Drawing Drawing

The voltage equation of the Synchronous Generator becomes: (1)vΞ±Ξ²=βˆ’rsiΞ±Ξ²βˆ’LsdiΞ±Ξ²dt+eΞ±Ξ²

The EMF voltage vector eΞ±Ξ² is defined as: (2)eΞ±Ξ²=Ο•dΟ‰r[cos⁑(ΞΈr)sin⁑(ΞΈr)]

The eletromagnetic torque Te is: (3)Te=Ο•diΞ±Ξ²T[cos⁑(ΞΈr)sin⁑(ΞΈr)]

The torque and rotor speed relationship are derived From the Newton's second law, as: (4)dΟ‰rdt=1J(Tiβˆ’Teβˆ’DpΟ‰r) (5)dΞΈrdt=Ο‰r where Ti is the input torque, J is the moment of inertia and Dp is a damping coefficient.

The reactive power associated with the EMF Qgap is given by: (6)Qgap=ψdΟ‰riΞ±Ξ²T[βˆ’sin⁑(ΞΈr)cos⁑(ΞΈr)]

and the reactive power reactive exchanged with grid at the PoC Qpoc: (7)Qpoc=(vΞ±Ξ²[01βˆ’10]T)iΞ±Ξ²

Consider the equivalent circuit illustrated in Fig. 1:

The current I and aparent power S from circuit shown in Fig. 2 can be found as:

(8)I=Ecos⁑(Ξ΄)+iEsin⁑(Ξ΄)βˆ’ViXL

(9)S=(Ecos⁑(Ξ΄)+iEsin⁑(Ξ΄))(Esin⁑(Ξ΄)+i(Ecos⁑(Ξ΄)βˆ’V))XL

The real power P and reactive power Q, delivered from the voltage source E to the grid V through the reactance XL, are:

(10)S=P+iQ

(11)P=VXLEsin⁑(δ)

(12)Q=E2XLβˆ’VXLEcos⁑(Ξ΄)

The eletromagnetic torque Te and input torque Ti are:

(13)Te=Pω

(14)Ti=Poωo

The VSG Control Block diagram is:

The nonlinear dynamic equations that describe the behavior of the angular frequency, power angle and VSG flux are:

(15)ddtΟ‰=1J(βˆ’Dp(Ο‰βˆ’Ο‰o)βˆ’Vψsin⁑(Ξ΄)XL+PoΟ‰o)

(16)ddtΞ΄=Ο‰βˆ’Ο‰g

(17)ddtψ=1K(Dq(Eoβˆ’V)+Qo+Ο‰ΟˆVXLcos⁑(Ξ΄)βˆ’Ο‰2ψ2XL)

where the flux loop constant K and moment of inertia J are given from:

(18)J=Ο„fDp

(19)K=τvωgDq

where Ο„f and Ο„v are the time constants associated with the voltage and frequency Droop variables, Dp and Dq, respectively.

In order to simplify the notation let us define some constants, that are:

a11=βˆ’DpJ b11=1J(PoΟ‰o+DpΟ‰o)

a12=βˆ’VJXL b21=βˆ’Ο‰g

a32=βˆ’1KXL b31=1K(Qo+Dq(Eoβˆ’V))

a33=VKXL

then:

(20)ddtΟ‰=a11Ο‰+a12ψsin⁑δ+b11

(21)ddtΞ΄=Ο‰+b12

(22)ddtψ=a32Ο‰2ψ2+a33Ο‰Οˆcos⁑(Ξ΄)+b31

The next step consists in compute the equilibrium point of the nonlinear dynamic equations:

(23)0=a11Ο‰+a12ψsin⁑δ+b11

(24)0=Ο‰+b12

(25)0=a32Ο‰2ψ2+a33Ο‰Οˆcos⁑(Ξ΄)+b31

or:

(26)Ο‰=Ο‰g

(27)Ο‰gψVXLsin⁑(Ξ΄)=Ο‰g(Dp(Ο‰oβˆ’Ο‰g)+PoΟ‰o)

(28)Ο‰gψVXLcos⁑(Ξ΄)=Dq(Vβˆ’Eo)βˆ’Qo+Ο‰g2ψ2XL

From (27) and (28) we have:

(29)(Ο‰gψVXL)2=[Ο‰g(Dp(Ο‰oβˆ’Ο‰g)+PoΟ‰o)]2+(Dq(Vβˆ’Eo)βˆ’Qo+Ο‰g2ψ2XL)2

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:

(30)J=[a11a12Ξ¨ocos⁑(Ξ΄o)a12sin⁑(Ξ΄o)1002a32Ο‰oΞ¨o2+a33Ξ¨ocos⁑(Ξ΄o)βˆ’a33Ξ¨oΟ‰osin⁑(Ξ΄o)2a32Ο‰o2Ξ¨o+a33Ο‰ocos⁑(Ξ΄o)]

Example:ΒΆ

In [12]:
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") )  
 XL =  0.21  Ξ©
 Dp =  4.64
 Dq =  1500.00
In [13]:
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") )  
 Psi =  0.48 T
 Delt =  0.00 Β°
 Ji =  0.05
 K =  28274.33

The Jacobian matrix is given by:

In [14]:
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)     
 J = 
 [[-1.00e+02 -1.06e+04 -0.00e+00]
 [ 1.00e+00  0.00e+00  0.00e+00]
 [-1.12e-02 -0.00e+00 -8.78e+00]]

And the eigenvalues are:

In [15]:
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 ] =  ') )   
 [ λ₁  Ξ»β‚‚  λ₃ ]'=  [[ -8.78 +0.j  ]
                    [-50.  +89.84j]
                    [-50.  -89.84j]]

In order to characterizes the interaction between modes and state variables, let us compute the participation factors:

In [16]:
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" )
 R = 
 [[ 0.00e+00+0.j   -1.00e+00+0.j   -1.00e+00-0.j  ]
 [ 0.00e+00+0.j    4.73e-03+0.01j  4.73e-03-0.01j]
 [ 1.00e+00+0.j   -4.71e-05-0.j   -4.71e-05+0.j  ]]
 
 L = 
 [[ 1.00e-05+0.00e+00j  1.21e-02+4.34e-19j  1.00e+00+0.00e+00j]
 [-5.00e-01-2.78e-01j  0.00e+00-5.88e+01j  0.00e+00-0.00e+00j]
 [-5.00e-01+2.78e-01j  0.00e+00+5.88e+01j  0.00e+00+0.00e+00j]]

 
 Participation Factors: 
 
               λ₁                         Ξ»β‚‚                          λ₃ 
 
|L(0,0)*R(0,0)|= 0.00      |L(1,0)*R(0,1)|= 0.57       |L(2,0)*R(0,2)|= 0.57        Ο‰
|L(0,1)*R(1,0)|= 0.00      |L(1,1)*R(1,1)|= 0.57       |L(2,1)*R(1,2)|= 0.57        Ξ΄
|L(0,2)*R(2,0)|= 1.00      |L(1,2)*R(2,1)|= 0.00       |L(2,2)*R(2,2)|= 0.00        ψ

Static Caracteristics - Voltage DroopΒΆ

In [17]:
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
In [18]:
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)

Static Caracteristics - Frequency DroopΒΆ

In [19]:
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
In [20]:
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)