MODULE 5

Virtual Synchronous Generator - VSG

Feedback of the Reactive Power at the PoC

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(TiTeDpω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αβ[0110]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=E2XLVXLEcos(δ)

The eletromagnetic torque Te and input torque Ti are:

(13)Te=Pω

(14)Ti=Poωo

Block diagram of VSG implementation for feedback of the reactive power at the PoC:

The 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(EoV)+QoVωψcos(δ)XL+V2XL)

From the dynamic equations, the equilibrium point can be obtained:

(18)0=1J(Dp(ωωo)Vψsin(δ)XL+Poωo)

(19)0=ωωg

(20)0=1K(Dq(EoV)+QoVωψcos(δ)XL+V2XL)

or:

(21)ω=ωg

(22)ωgψVXLsin(δ)=ωg(Dp(ωoωg)+Poωo)

(23)ωgψVXLcos(δ)=Dq(EoV)+Qo+V2XL

From (22) and (23) we have:

(24)(ωgψVXL)2=[ωg(Dp(ωoωg)+Poωo)]2+(Dq(EoV)+Qo+V2XL)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:

(11)J=[DpJVoψocos(δo)JXLVosin(δo)JXL100Voψocos(δo)KXLVoψoωosin(δo)KXLVoωocos(δo)KXL]

Example:

In [50]:
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 [51]:
import numpy as np 

Vo = V
a  = ((wg*V)/XL)**2
b  = wg*(Dp*(wo-wg) + Po/wo)
c  = Dq*(Eo - V) + Qo + (V**2)/XL

Psio  = np.sqrt((b**2 + c**2)/a)
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.50 T
 Delt =  0.00 °
 Ji =  0.05
 K =  28274.33

The Jacobian matrix is given by:

In [52]:
import math
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], [-Vo*Psio*math.cos(Delto)/(K*XL), Vo*wo*Psio*math.sin(Delto)/(K*XL), -Vo*wo*math.cos(Delto)/(K*XL)]])

print (" J = \n", J)     
 J = 
 [[-1.00e+02 -1.10e+04 -0.00e+00]
 [ 1.00e+00  0.00e+00  0.00e+00]
 [-1.80e-02  0.00e+00 -1.37e+01]]
In [53]:
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 ] =  ') )   
 [ λ₁  λ₂  λ₃ ]'=  [[-13.65 +0.j  ]
                    [-50.  +92.09j]
                    [-50.  -92.09j]]

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

In [54]:
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.55e-03+0.01j  4.55e-03-0.01j]
 [ 1.00e+00+0.j   -6.69e-05-0.j   -6.69e-05+0.j  ]]
 
 L = 
 [[ 2.51e-05-1.36e-20j  2.02e-02+1.73e-18j  1.00e+00+0.00e+00j]
 [-5.00e-01-2.71e-01j  3.08e-15-5.96e+01j  0.00e+00-0.00e+00j]
 [-5.00e-01+2.71e-01j -3.08e-15+5.96e+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 [55]:
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)
cj = Dq*(Eo - Vj) + Qo + (Vj**2)/XL


Psi1j = np.sqrt((b**2 + cj**2)/aj)

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.sin(Delt1j))/XL
Q2j = (E1j*Vj*np.cos(Delt1j) - Vj**2)/XL
In [57]:
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.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{Q2j}{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 [39]:
PN = 10
j  = np.arange(start=0, stop=21, step=1)
DW = 2.*pi*3.

V  = 220
wrj = 2.*pi*(60. - 6.) + (j/PN)*2.*DW
wg = wrj  

L  = 1.e-3
XLj = wrj*L 

aj = ((wg*V)/XLj)**2
b  = wg*(Dp*(wo-wg) + Po/wo)
cj = Dq*(Eo - V) + Qo + (V**2)/XLj


Psi1j = sqrt((b**2 + cj**2)/aj)

Delt1j = np.arctan2( wg*(Dp*(wo-wg)+Po/wo) , Dq*(V-Eo)-Qo+(wg**2*Psi1j**2)/XLj)


E1j = Psi1j*wg 
P1j = V*E1j*np.sin(Delt1j)/XLj
Q1j = (E1j**2 - E1j*V*np.sin(Delt1j))/XLj
Q2j = (E1j*V*np.cos(Delt1j) - V**2)/XLj
In [41]:
import matplotlib.pyplot as plt

plt.figure(figsize=(12, 4
                   ))
plt.subplot(121)
plt.plot(wrj/(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(wrj/(2*pi*60), Q2j/33000, color="blue", label="${Q_{poc}}$")
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{Q2j}{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)