MODULE 5

Small Signal Analyzes of the Inductive Droop

Integrator in the Reactive Power - Voltage Loop

Grid-Forming Inverter with Inductive Droop Control and Integrator in the Reactive Power - Voltage Loop:

Drawing Drawing

The real power p and reactive power q, delivered from the voltage source e to the grid V through the reactance XL, are:

(1)p=VXLesin(δ)

(2)q=VXLecos(δ)V2XL

where δ is the power angle. The active and reactive power are usually filtered, here this filter is implemented as:

(3)ddtpf=apf+ap

(4)ddtqf=apf+aq

where a defines the low-pass filter cutoff frequency.

The Droop Law is:

(5)ed=(EoV)kvkp(qfQo)

(6)ω=ωom(pfPo)

where ed is the input of the integrator of q channel and m, kp and kv corresponds to the Droop coefficients. Therefore, these equations can be rewritten as:

(7)qf=Qokqed+kv(EoV)kq

(8)pf=ωo+mPoωm

By assuming that Po, Qo, Eo and ωo are constant, it can be considered that:

(9)ddted=kpddtqf

(10)ddtω=mddtpf

and by rearanging the above equations, it is possible to obtain:

(11)ddted=kq(aqf+aq)=kq[a(Qokded+kv(EoV)kq)+a(VXLecos(δ)V2XL)]

(14)ddtω=m(apf+ap)=m[a(ωo+mPoωm)+a(VXLesin(δ))]

(13)ddted=aedkqVaXLecos(δ)+akqQo+kqV2aXL+akv(EoV)

(14)ddtω=aωamVXLesin(δ)+aωo+amPo

Note that ω is the Droop control variable that by integration will define the angular position of the voltage sinthesized by the inverter, that is:

(15)ddtθ=ω

Since we are consider the operation in sinusoidal steady state, the θ is given by:

(16)θ=ωrt+δ

therefore:

(17)ddtθ=ωr+ddtδ

where:

(18)ddtδ=ωωr

The nonlinear dynamic equations that govern the grid connected converter with droop are described by:

(19)ddtδ=ωωr

(20)ddtω=aωamVXLesin(δ)+a(ωo+mPo)

(21)ddted=aedkqVaXLecos(δ)+akqQo+kqV2aXL+akv(EoV)

(22)ddte=ed

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

a21=a b11=ωr

a22=mVaXL b21=a(ωo+mPo)

a32=kqVaXL b31=akqQo+kqV2aXL+akv(EoV)

then:

(23)ddtδ=ω+b11

(24)ddtω=a21ω+a22esin(δ)+b21

(25)ddted=a21ed+a32ecos(δ)+b31

(26)ddte=ed

Let us start by computing the equilibrium point of the nonlinear dynamic equations:

(27)0=ω+b11ω=b11

(28)0=a21ω+a22esin(δ)+b21esin(δ)=b21+a21b11a22

(29)0=a21ed+a32ecos(δ)+b31ecos(δ)=b31a32

(30)ed=0

From (28) and (29) we have:

(31)e2=(b21+a21b11a22)2+(b31a32)2

or:

(32)e=(b21+a21b11a22)2+(b31a32)2

The load angle can be found as:

(33)δ=atan2(b31a32,b21+a21b11a22)

and the inverter reference frequency:

(34)ω=b11

Example:

In [66]:
pi = 3.14159265359
V  = 220.         # Grid Voltage
wr = 2.*pi*60.    # Grid Frequency
L  = 1.e-3        # Inductance 
XL = wr*L         # Inductive Reactance 
a  = 2.*pi*30.         # Low-pass filters
Zb = V**2/33000.       # Based Impedance
m  = 2.*pi*3./33.33e3  # Droop coefficient m
kv = 20.               # Droop coefficient kv
kq = kv*n              # Droop coefficient kq
print (" XL = ", format(XL, ".2f") , " \u03A9")  
print (" XL = ", format(XL/Zb, ".2f"), " pu ")  
 XL =  0.38  Ω
 XL =  0.26  pu 

The reference voltages and frequency will be expressed as a function of the ative and reative power to be dispached to the grid:

In [67]:
Qo = 33e3 * 0.
Eo = V

Po = 33e3 * 1.
wo = wr 

print (" Eo = ", format(Eo, ".2f"),  " V")
print (" wo = ", format(wo, ".2f"),  " rad/s" )
print (" Po = ", format(Po, ".2f"),  "   W" )
print (" Qo = ", format(Qo, ".2f"),  "   VAr" )
 Eo =  220.00  V
 wo =  376.99  rad/s
 Po =  33000.00    W
 Qo =  0.00    VAr
In [68]:
a21 = - a
a22 = - (a*m*V)/XL
a32 = - (a*kq*V)/XL
b11 = - wr
b21 = a*(wo + m*Po)
b31 = a*kq*Qo + (kq*a*V**2)/XL + a*kv*(Eo - V)

print (" a21 = ", format(a21, ".2f"))
print (" a22 = ", format(a22, ".2f"))
print (" a32 = ", format(a32, ".2f"))
print (" b11 = ", format(b11, ".2f"))
print (" b21 = ", format(b21, ".2f"))
print (" b31 = ", format(b31, ".2f"))
 a21 =  -188.50
 a22 =  -62.21
 a32 =  -1452.15
 b11 =  -376.99
 b21 =  74579.03
 b31 =  319471.95
In [69]:
e = (((-b21 + a21*b11)/a22)**2 + (-b31/a32)**2)**(1/2)

print (" e = ", format(e, ".2f"), "  V")
 e =  227.15   V
In [70]:
import math

delt = math.atan2( (-b21 + a21*b11)/a22 , -b31/a32)

Delt = delt*180./pi

print (" Delt = ", format(Delt, ".2f"), "\u00b0")                  
 Delt =  14.42 °
In [71]:
p = (V/XL)*e*math.sin(delt)

q = (V/XL)*e*math.cos(delt) - (V**2)/XL

print (" p = ", format(p/1000, ".2f") , "kW")        
print (" q = ", format(q/1000, ".2f") , " kVAr")                  
 p =  33.00 kW
 q =  0.00  kVAr

The Jacobian matrix is given by:

In [72]:
import numpy as np 
from scipy import sqrt, exp, cos, matrix, vstack, hstack, zeros

np.set_printoptions(precision=2)

J = matrix([[0, 1, 0, 0], [(a22*e*math.cos(delt)), a21, 0, (a22*math.sin(delt))], [(-a32*e*math.sin(delt)), 0, a21, a32*math.cos(delt)], [0, 0, 1, 0]])

print (" J = \n", J)                              
 J = 
 [[ 0.00e+00  1.00e+00  0.00e+00  0.00e+00]
 [-1.37e+04 -1.88e+02  0.00e+00 -1.55e+01]
 [ 8.21e+04  0.00e+00 -1.88e+02 -1.41e+03]
 [ 0.00e+00  0.00e+00  1.00e+00  0.00e+00]]
In [73]:
import scipy.linalg as la

evals, evecs = la.eig(J)    
    
Eig_V_J = evals  

print (" [ \u03BB\u2081  \u03BB\u2082  \u03BB\u2083  \u03BB\u2084]'= ",  np.array2string(Eig_V_J.reshape((-1, 1)), prefix=' [ \u03BB\u2081  \u03BB\u2082  \u03BB\u2083 \u03BB\u2084] =  ') )   
 [ λ₁  λ₂  λ₃  λ₄]'=  [[  -8.39 +0.j  ]
                      [ -94.25+68.55j]
                      [ -94.25-68.55j]
                      [-180.11 +0.j  ]]

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

In [74]:
evals, evecs = la.eig(J)   

T = evals.real 

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                          \u03BB\u2084 \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"), "      |L(3,0)*R(0,3)|=", format(np.abs(L[3][0]*R[0][3]), ".2f"), "       \u03B4")
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"), "      |L(3,1)*R(1,3)|=", format(np.abs(L[3][1]*R[1][3]), ".2f"), "       \u03C9")
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"), "      |L(3,2)*R(2,3)|=", format(np.abs(L[3][2]*R[2][3]), ".2f"), "       ed" )
print("|L(0,3)*R(3,0)|=", format(np.abs(L[0][3]*R[3][0]), ".2f"), "     |L(1,3)*R(3,1)|=", format(np.abs(L[1][3]*R[3][1]), ".2f"), "      |L(2,3)*R(3,2)|=", format(np.abs(L[2][3]*R[3][2]), ".2f"), "      |L(3,3)*R(3,3)|=", format(np.abs(L[3][3]*R[3][3]), ".2f"), "       e" )
 R = 
 [[-1.51e-04+0.00e+00j -1.02e-03-7.40e-04j -1.02e-03+7.40e-04j
   7.06e-06+0.00e+00j]
 [ 1.26e-03+0.00e+00j  1.47e-01-7.63e-17j  1.47e-01+7.63e-17j
  -1.27e-03+0.00e+00j]
 [-9.93e-01+0.00e+00j -9.89e-01+0.00e+00j -9.89e-01-0.00e+00j
   1.00e+00+0.00e+00j]
 [ 1.18e-01+0.00e+00j  6.86e-03+4.99e-03j  6.86e-03-4.99e-03j
  -5.55e-03+0.00e+00j]]
 
 L = 
 [[ 6.03e+01-1.14e-13j  3.35e-01+0.00e+00j  4.96e-02-8.74e-19j
   8.94e+00-1.12e-16j]
 [-4.91e-13+6.81e+02j  3.44e+00+4.73e+00j  4.37e-03+6.01e-03j
  -1.25e-14+8.67e-01j]
 [-3.86e-13-6.81e+02j  3.44e+00-4.73e+00j  4.37e-03-6.01e-03j
  -1.23e-14-8.67e-01j]
 [ 5.99e+01-5.80e-14j  7.14e+00-4.24e-16j  1.06e+00-5.06e-19j
   8.87e+00-7.38e-17j]]


 Participation Factors: 
 
               λ₁                         λ₂                         λ₃                          λ₄ 
 
|L(0,0)*R(0,0)|= 0.01      |L(1,0)*R(0,1)|= 0.86       |L(2,0)*R(0,2)|= 0.86       |L(3,0)*R(0,3)|= 0.00        δ
|L(0,1)*R(1,0)|= 0.00      |L(1,1)*R(1,1)|= 0.86       |L(2,1)*R(1,2)|= 0.86       |L(3,1)*R(1,3)|= 0.01        ω
|L(0,2)*R(2,0)|= 0.05      |L(1,2)*R(2,1)|= 0.01       |L(2,2)*R(2,2)|= 0.01       |L(3,2)*R(2,3)|= 1.06        ed
|L(0,3)*R(3,0)|= 1.06      |L(1,3)*R(3,1)|= 0.01       |L(2,3)*R(3,2)|= 0.01       |L(3,3)*R(3,3)|= 0.05        e

Static Caracteristics - Voltage Droop

In [75]:
PN = 100
j = np.arange(start=1, stop=101, step=1)
DV= 22.

Vj= 220 - DV + j*(2*DV)/PN

wr = 2.*pi*60.
L  = .5e-3
XL = wr*L 
n  = 22./(33.33e3)
m  = (2.*pi*3.)/(33.33e3)
a  = 2.*pi*12.
wo = wr
Eo = 220.
Po = 0.
Qo = 0.
kv = 1.
kq = n*kv

b11 = -wr
a21 = -a
b21 = a*(wo + m*Po)
a22j = - a*m*Vj/XL
b31j = a*kq*Qo + kq*(Vj**2)*a/XL + a*kv*(Eo - Vj)
a32j = - kq*Vj*a/XL

ej = sqrt(((-b21 + a21*b11)/a22j)**2 + (-b31j/a32j)**2)

deltj = np.arctan2( (-b21 + a21*b11)/a22j , -b31j/a32j)

qj = (Vj/XL)*ej*np.cos(deltj) - (Vj**2)/XL

pj = (Vj/XL)*ej*np.sin(deltj)
In [76]:
import matplotlib.pyplot as plt

plt.figure(figsize=(12, 4
                   ))
plt.subplot(121)
plt.plot(Vj, qj/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{qj}{33000}$', rotation=0) 
plt.xlabel('$\dfrac{Vj}{220}$') 

plt.subplot(122)
plt.plot(Vj, pj/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{pj}{33000}$', rotation=0) 
plt.xlabel('$\dfrac{Vj}{220}$') 

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 [77]:
PN = 100
j  = np.arange(start=1, stop=101, step=1)
DW = 2.*pi*3.

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

L  = 1.e-3
XLj = wrj*L 
n  = 22./(33.33e3)
m  = (2.*pi*3.)/(33.33e3)
a  = 2.*pi*12.
wo = 2.*pi*60.
Eo = 220.
Po = 0.
Qo = 0.
kv = 1.
kq = n*kv

b11j = -wrj
a21 = -a
b21 = a*(wo + m*Po)
a22j = - a*m*V/XLj
b31j = a*kq*Qo + kq*(V**2)*a/XLj + a*kv*(Eo - V)
a32j = - kq*V*a/XLj

ej = sqrt(((-b21 + a21*b11j)/a22j)**2 + (-b31j/a32j)**2)

deltj = np.arctan2( (-b21 + a21*b11j)/a22j , -b31j/a32j)

qj = (V/XLj)*ej*np.cos(deltj) - (V**2)/XLj

pj = (V/XLj)*ej*np.sin(deltj)
In [78]:
import matplotlib.pyplot as plt

plt.figure(figsize=(12, 4
                   ))
plt.subplot(121)
plt.plot(wrj/(2*pi*60), qj/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{qj}{33000}$', rotation=0) 
plt.xlabel('$\dfrac{\u03C9 rj}{2 \u03C0 60}$') 

plt.subplot(122)
plt.plot(wrj/(2*pi*60), pj/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{pj}{33000}$', rotation=0) 
plt.xlabel('$\dfrac{\u03C9 rj}{2 \u03C0 60}$') 

plt.subplots_adjust(left=0.125,
                    bottom=0.1, 
                    right=0.9, 
                    top=0.9, 
                    wspace=0.4, 
                    hspace=0.35)