Grid-Forming Inverter with Inductive Droop Control:
![]() |
![]() |
![]() |
The real power and reactive power , delivered from the voltage source to the grid through the reactance , are:
where is the power angle. The active and reactive power are usually filtered, here this filter is implemented as:
where defines the low-pass filter cutoff frequency.
The Droop Law is:
where and corresponds to the Droop coefficients. Therefore, these equations can be rewritten as:
By assuming that , , and are constant, it can be considered that:
and by rearanging the above equations, it is possible to obtain:
Note that is the Droop control variable that by integration will define the angular position of the voltage sinthesized by the inverter, that is:
Since we are consider the operation in sinusoidal steady state, the is given by:
therefore:
where:
The nonlinear dynamic equations that govern the grid connected converter with droop are described by:
In order to simplify the notation let us define some constants, that are:
then:
Let us start by computing the equilibrium point of :
From and we have:
The quadratic equation described in has two solutions. Let us call them as and , that are:
The load angle can be found as:
and the inverter reference frequency:
pi = 3.14159265359
V = 220. # Grid Voltage
wr = 2.*pi*60. # Grid Frequency
L = 2.e-3 # Inductance
XL = wr*L # Inductive Reactance
n = 11./(33.33e3) # Droop coefficient n
m = 2.*pi*3./33.33e3 # Droop coefficient m
a = 2.*pi*12. # Low-pass filters
Zb = V**2/33000. # Based Impedance
print (" XL = ", format(XL, ".2f") , " \u03A9")
print (" XL = ", format(XL/Zb, ".2f"), " pu ")
The reference voltages and frequency will be expressed as a function of the ative and reative power to be dispached to the grid:
Qo = 33e3 * 0.
Eo = V
Po = 33e3 * 0.
wo = wr
Qo = ((V**2 + (Po/(V/XL))**2)**(1/2) - Eo)/n
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" )
a21 = - a
a22 = - (a*m*V)/XL
a32 = - (a*n*V)/XL
b11 = - wr
b21 = a*(wo + m*Po)
b31 = a*(Eo + n*Qo) + (n*a*V**2)/XL
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"))
e1 = (a22**2)*(a32**2) * (( (2*(a21**3)*(b11)*(b21) - (a21**4)*(b11**2) + (a21**2)*(a32**2)*(b11**2) - (a21**2)*(b21**2) - 2*a21*(a32**2)*b11*b21 + (a22**2)*(b31**2) + (a32**2)*(b21**2))**(1/2)) /(a22*a32) - ((a21*b31)/(a32**2))) / ((a21**2)*(a22**2) - (a22**2)*(a32**2))
e2 = - (a22**2)*(a32**2) * (( (2*(a21**3)*(b11)*(b21) - (a21**4)*(b11**2) + (a21**2)*(a32**2)*(b11**2) - (a21**2)*(b21**2) - 2*a21*(a32**2)*b11*b21 + (a22**2)*(b31**2) + (a32**2)*(b21**2))**(1/2)) /(a22*a32) + ((a21*b31)/(a32**2))) / ((a21**2)*(a22**2) - (a22**2)*(a32**2))
print (" e1 = ", format(e1, ".2f"), " V")
print (" e2 = ", format(e2, ".2f"), " V")
import math
delt1 = math.atan2( XL*(wo + (m*Po) - wr) / (V*m) , (V + ((XL*(Eo + n*Qo - e1))/(V*n))))
Delt1 = delt1*180./pi
delt2 = math.atan2( XL*(wo + (m*Po) - wr) / (V*m) , (V + ((XL*(Eo + n*Qo - e2))/(V*n))))
Delt2 = delt2*180./pi
print (" Delt_1 = ", format(Delt1, ".2f"))
print (" Delt_2 = ", format(Delt2, ".2f"))
p1 = (V/XL)*e1*math.sin(delt1)
q1 = (V/XL)*e1*math.cos(delt1) - (V**2)/XL
p2 = (V/XL)*e2*math.sin(delt2)
q2 = (V/XL)*e2*math.cos(delt2) - (V**2)/XL
print (" p1 = ", format(p1/1000, ".2f") , " kW")
print (" q1 = ", format(q1/1000, ".2f") , "kVAr")
print (" p2 = ", format(p2/1000, ".2f"), " kW")
print (" q2 = ", format(q2/1000, ".2f"), " kVAr")
The Jacobian matrix for the first equilibrium point, described as and , is given by:
import numpy as np
from scipy import sqrt, exp, cos, matrix, vstack, hstack, zeros
np.set_printoptions(precision=2)
e = e1
delt = delt1
J = matrix([[0, 1, 0], [(a22*e*math.cos(delt)), a21, (a22*math.sin(delt))], [(-a32*e*math.sin(delt)), 0, (a21 + a32*math.cos(delt))]])
print (" J = \n", J)
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 ] = ') )
Note that if an eigenvalue has a positive real part, this indicates that the considered equilibrium point is unstable.
In order to characterizes the interaction between modes and state variables, let us compute the participation factors:
evals, evecs = la.eig(J)
T = evals.real
R = evecs.real # 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"), " \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"), " \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"), " e" )
Let us now consider the second equilibrium point, and , and calculate the Jacobian matrix, as:
e = e2
delt = delt2
J = matrix([[0, 1, 0], [(a22*e*math.cos(delt)), a21, (a22*math.sin(delt))], [(-a32*e*math.sin(delt)), 0, (a21 + a32*math.cos(delt))]])
print (" J = \n", J)
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 ] = ') )
Note that these eigenvalues indicate that the considered equilibrium point is stable.
In order to characterizes the interaction between modes and state variables, let us compute the participarion factors:
np.set_printoptions(linewidth=np.inf)
evals, evecs = la.eig(J)
T = evals
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"), " \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"), " \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"), " e" )
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 = 33000.
Qo = 0.
b11 = -wr
a21 = -a
b21 = a*(wo + m*Po)
a22j = - a*m*Vj/XL
b31j = a*(Eo + n*Qo) + n*(Vj**2)*a/XL
a32j = - n*Vj*a/XL
e2j = - (a22j**2)*(a32j**2) * (( (2*(a21**3)*(b11)*(b21) - (a21**4)*(b11**2) + (a21**2)*(a32j**2)*(b11**2) - (a21**2)*(b21**2) - 2*a21*(a32j**2)*b11*b21 + (a22j**2)*(b31j**2) + (a32j**2)*(b21**2))**(1/2)) /(a22j*a32j) + ((a21*b31j)/(a32j**2))) / ((a21**2)*(a22j**2) - (a22j**2)*(a32j**2))
delt2j = np.arctan2( XL*(wo + (m*Po) - wr) / (Vj*m) , (Vj + ((XL*(Eo + n*Qo - e2j))/(Vj*n))))
qj = (Vj/XL)*e2j*np.cos(delt2j) - (Vj**2)/XL
pj = (Vj/XL)*e2j*np.sin(delt2j)
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)
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 = .1e-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.
b11j = - wrj
a21 = - a
b21 = a*(wo + m*Po)
a22j = - a*m*V/XLj
b31j = a*(Eo + n*Qo) + n*(V**2)*a/XLj
a32j = - n*V*a/XLj
e2j = - (a22j**2)*(a32j**2) * (( (2*(a21**3)*(b11j)*(b21) - (a21**4)*(b11j**2) + (a21**2)*(a32j**2)*(b11j**2) - (a21**2)*(b21**2) - 2*a21*(a32j**2)*b11j*b21 + (a22j**2)*(b31j**2) + (a32j**2)*(b21**2))**(1/2)) /(a22j*a32j) + ((a21*b31j)/(a32j**2))) / ((a21**2)*(a22j**2) - (a22j**2)*(a32j**2))
delt2j = np.arctan2( XLj*(wo + (m*Po) - wrj) / (V*m) , (V + ((XLj*(Eo + n*Qo - e2j))/(V*n))))
qj = (V/XLj)*e2j*np.cos(delt2j) - (V**2)/XLj
pj = (V/XLj)*e2j*np.sin(delt2j)
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)