MODULE 5¶

Small Signal Analyzes of the Inductive Droop¶

Grid-Forming Inverter with Inductive Droop Control:

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

$p = \frac{V}{X_L} e \sin(\delta) \tag{1}$

$q = \frac{V}{X_L} e \cos(\delta) - \frac{V^2}{X_L} \tag{2}$

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

$\frac{d}{dt}p_f = - a p_f + a p \tag{3}$

$\frac{d}{dt}q_f = - a p_f + a q \tag{4}$

where $a$ defines the low-pass filter cutoff frequency.

The Droop Law is:

$\omega = \omega_o - m(p_f - P_o) \tag{5}$

$e = E_o - n(q_f - Q_o) \tag{6}$

where $n$ and $m$ corresponds to the Droop coefficients. Therefore, these equations can be rewritten as:

$p_f = \frac{\omega_o + m P_o - \omega}{m} \tag{7}$

$q_f = \frac{E_o + n Q_o - e}{n} \tag{8}$

By assuming that $P_o$ , $Q_o$ , $E_o$ and $\omega_o$ are constant, it can be considered that:

$\frac{d}{dt} e = -n \frac{d}{dt} q_f \tag{9}$

$\frac{d}{dt} \omega = -m \frac{d}{dt} p_f \tag{10}$

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

$\frac{d}{dt} e = -n (-a q_f + a q) = -n \biggl[ -a \biggl(\frac{E_o + n Q_o - e}{n}\biggr) + a \biggl(\frac{V}{X_L} e \cos(\delta) - \frac{V^2}{X_L}\biggr)\biggr] = -a e - \frac{n V a}{X_L} e \cos(\delta) + a E_o + a n Q_o + \frac{n V^2 a}{X_L} \tag{11}$

$\frac{d}{dt} \omega = -m (-a p_f + a p) = -m \biggl[ -a \biggl(\frac{\omega_o + m P_o - \omega}{m}\biggr) + a \biggl(\frac{V}{X_L} e \sin(\delta)\biggr)\biggr] = -a \omega - \frac{a m V}{X_L} e \sin(\delta) + a \omega_o + a m P_o \tag{12}$

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

$\frac{d}{dt} \theta = \omega \tag{13}$

Since we are consider the operation in sinusoidal steady state, the $\theta$ is given by:

$\theta = \omega_r t + \delta \tag{14}$

therefore:

$\frac{d}{dt} \theta = \omega_r + \frac{d}{dt} \delta \tag{15}$

where:

$\frac{d}{dt} \delta = \omega - \omega_r \tag{16}$

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

$\frac{d}{dt} \delta = \omega - \omega_r \tag{16}$

$\frac{d}{dt} \omega = -a \omega - \frac{a m V}{X_L} e \sin(\delta) + a (\omega_o + m P_o) \tag{17}$

$\frac{d}{dt} e = -a e - \frac{n V a}{X_L} e \cos(\delta) + a (E_o + n Q_o) + \frac{n V^2 a}{X_L} \tag{18}$

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

$a_{21} = -a$ $\quad$ $\quad$ $\quad$ $b_{11} = - \omega_r$

$a_{22} = - \frac{m V a}{X_L}$ $\quad$ $\quad$ $b_{21} = a ( \omega_o + m P_o)$

$a_{32} = - \frac{n V a}{X_L}$ $\quad$ $\quad$ $b_{31} = a(E_o + n Q_o) + \frac{n V^2 a}{X_L}$

then:

$\frac{d}{dt} \delta = \omega + b_{11} \tag{19}$

$\frac{d}{dt} \omega = a_{21} \omega + a_{22} e \sin(\delta) + b_{21} \tag{20}$

$\frac{d}{dt} e = a_{21} e + a_{32} e \cos(\delta) + b_{31} \tag{21}$

Let us start by computing the equilibrium point of $(20)$ :

$0 = \omega + b_{11} \quad \quad \qquad \, \, \rightarrow \quad \quad \omega = - b_{11} \tag{22}$

$0 = a_{21} \omega + a_{22} e \sin(\delta) + b_{21} \quad \rightarrow \quad e \sin(\delta) = \frac{- b_{21} + a_{21} b_{11}}{a_{22}} \tag{23}$

$0 = a_{21} e + a_{32} e \cos(\delta) + b_{31} \quad \rightarrow \quad e \cos(\delta) = \frac{- b_{31} + a_{21} E}{a_{32}} \tag{24}$

From $(23)$ and $(24)$ we have:

$e^2 = \biggl(\frac{-b_{21} + a_{21} b_{11}}{a_{22}}\biggr)^2 + \biggl(\frac{-b_{31} - a_{21} e}{a_{32}}\biggr)^2 \tag{25}$

The quadratic equation described in $(25)$ has two solutions. Let us call them as $e_1$ and $e_2$ , that are:

$e_1 = \frac{a^2_{22} a^2_{32} \biggl(\frac{\sqrt{2 a^3_{21} b_{11} b_{21} - a^4_{21} b^2_{11} + a^2_{21} a^2_{32} b^2_{11} - 2 a^2_{21} b^2_{21} - 2 a_{21} a^2_{32} b_{11} b_{21} + a^2_{22} b^2_{31} + a^2_{32} b^2_{21}}}{a_{22} a_{32}}-\frac{a_{21} b_{31}}{a^2_{32}}\biggr)}{a^2_{21} a^2_{22} - a^2_{22} a^2_{32}} \tag{26}$

$e_2 = - \frac{a^2_{22} a^2_{32} \biggl(\frac{\sqrt{2 a^3_{21} b_{11} b_{21} - a^4_{21} b^2_{11} + a^2_{21} a^2_{32} b^2_{11} - 2 a^2_{21} b^2_{21} - 2 a_{21} a^2_{32} b_{11} b_{21} + a^2_{22} b^2_{31} + a^2_{32} b^2_{21}}}{a_{22} a_{32}}+\frac{a_{21} b_{31}}{a^2_{32}}\biggr)}{a^2_{21} a^2_{22} - a^2_{22} a^2_{32}} \tag{27}$

The load angle can be found as:

$\delta_{1,2} = \text{atan2} \biggl( \frac{-b_{31} - a_{21} e_{1,2} }{a_{32}} \text{,} \frac{-b_{21} + a_{21} b_{11}}{a_{22}} \biggr) \tag{28}$

and the inverter reference frequency:

$\omega = - b_{11} \tag{29}$

Example:¶

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 ")

 XL =  0.75  Ω
 XL =  0.51  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" )

 Eo =  220.00  V
 wo =  376.99  rad/s
 Po =  0.00    W
 Qo =  0.00    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"))

 a21 =  -75.40
 a22 =  -12.44
 a32 =  -7.26
 b11 =  -376.99
 b21 =  28424.46
 b31 =  18184.97

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")

 e1 =  266.89   V
 e2 =  220.00   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"))

 Delt_1 =  180.00
 Delt_2 =  0.00

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")

 p1 =  0.00    kW
 q1 =  -142.07 kVAr
 p2 =  0.00    kW
 q2 =  -0.00   kVAr

The Jacobian matrix for the first equilibrium point, described as $e_1$ and $\delta_1$ , 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)

 J = 
 [[ 0.00e+00  1.00e+00  0.00e+00]
 [ 3.32e+03 -7.54e+01 -1.52e-15]
 [ 2.37e-13  0.00e+00 -6.81e+01]]

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 ] =  ') )

 [ λ₁  λ₂  λ₃ ]'=  [[  31.16+0.j]
                    [-106.56+0.j]
                    [ -68.14+0.j]]

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" )

 R = 
 [[ 3.21e-02  9.38e-03  3.99e-19]
 [ 9.99e-01 -1.00e+00 -2.72e-17]
 [ 7.14e-17 -7.15e-17  1.00e+00]]
 
 L = 
 [[ 2.41e+01  2.26e-01 -3.47e-18]
 [ 2.41e+01 -7.74e-01 -3.07e-17]
 [ 0.00e+00 -7.15e-17  1.00e+00]]

 
 Participation Factors: 
 
               λ₁                         λ₂                          λ₃ 
 
|L(0,0)*R(0,0)|= 0.77      |L(1,0)*R(0,1)|= 0.23       |L(2,0)*R(0,2)|= 0.00        δ
|L(0,1)*R(1,0)|= 0.23      |L(1,1)*R(1,1)|= 0.77       |L(2,1)*R(1,2)|= 0.00        ω
|L(0,2)*R(2,0)|= 0.00      |L(1,2)*R(2,1)|= 0.00       |L(2,2)*R(2,2)|= 1.00        e

Let us now consider the second equilibrium point, $e_2$ and $\delta_2$ , 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)

 J = 
 [[ 0.00e+00  1.00e+00  0.00e+00]
 [-2.74e+03 -7.54e+01 -0.00e+00]
 [ 0.00e+00  0.00e+00 -8.27e+01]]

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 ] =  ') )

 [ λ₁  λ₂  λ₃ ]'=  [[-37.7 +36.28j]
                    [-37.7 -36.28j]
                    [-82.66 +0.j  ]]

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" )

 R = 
 [[ 0.01+0.01j  0.01-0.01j  0.  +0.j  ]
 [-1.  +0.j   -1.  -0.j    0.  +0.j  ]
 [ 0.  +0.j    0.  -0.j    1.  +0.j  ]]
 
 L = 
 [[ 0. -37.73j -0.5 -0.52j  0.  -0.j  ]
 [ 0. +37.73j -0.5 +0.52j  0.  +0.j  ]
 [ 0.  +0.j    0.  +0.j    1.  +0.j  ]]

 
 Participation Factors: 
 
               λ₁                         λ₂                          λ₃ 
 
|L(0,0)*R(0,0)|= 0.72      |L(1,0)*R(0,1)|= 0.72       |L(2,0)*R(0,2)|= 0.00        δ
|L(0,1)*R(1,0)|= 0.72      |L(1,1)*R(1,1)|= 0.72       |L(2,1)*R(1,2)|= 0.00        ω
|L(0,2)*R(2,0)|= 0.00      |L(1,2)*R(2,1)|= 0.00       |L(2,2)*R(2,2)|= 1.00        e

Static Caracteristics - Voltage Droop¶

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)

Static Caracteristics - Frequency Droop¶

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)