3-phase VSI (PWM)

A sinusoidal pulse width modulated three-phase inverter given below is used to inject power to the three-phase balanced $400\,$V, $50\,$Hz AC voltage source. The DC side of the inverter is maintained at $800\,$V. The A phase modulating voltage is $m_a(t) = 0.816\,\sin(100\,\pi\,t + 30^{\circ})$. The inverter is connected to the grid through a series $RL$ circuit with $R=0.1\,\Omega$ and $L=0.01\,$H. Determine the following.

1. RMS values of the fundamental components of phase voltage and load current
2. active power injected to AC voltage source
3. phase difference between fundamental current through and voltage across the AC source
4. average current drawn from the DC supply?

from IPython.display import Image
Image(filename =r'VSI_3ph_6_fig_1.png', width=650)

# run this cell to view the circuit file.
%pycat VSI_3ph_6_orig.in

We now replace the strings such as \$Vdc, \$L, with the values of our choice by running the python script given below. It takes an existing circuit file VSI_3ph_6_orig.in and produces a new circuit file VSI_3ph_6.in, after replacing \$Vdc, \$L, etc. with values of our choice.

import gseim_calc as calc
s_Vdc = "800"
s_L = "0.01"
s_R = "0.1"
s_ea = "325"
s_sin = "0.816"
s_f_sin = "50"
s_f_carrier = "2e3"
s_dt_min = "0.1e-6"
s_dt_nrml = "5e-6"

phi_a_sin = 30.0
phi_b_sin = phi_a_sin - 120.0
phi_c_sin = phi_a_sin - 240.0

s_phi_a_sin = ("%11.4E"%(phi_a_sin)).strip()
s_phi_b_sin = ("%11.4E"%(phi_b_sin)).strip()
s_phi_c_sin = ("%11.4E"%(phi_c_sin)).strip()

l = [
  ('$Vdc', s_Vdc),
  ('$L', s_L),
  ('$R', s_R),
  ('$A_ea', s_ea),
  ('$A_sin', s_sin),
  ('$f_sin', s_f_sin),
  ('$f_carrier', s_f_carrier),
  ('$dt_min', s_dt_min),
  ('$dt_nrml', s_dt_nrml),
  ('$phi_a_sin', s_phi_a_sin),
  ('$phi_b_sin', s_phi_b_sin),
  ('$phi_c_sin', s_phi_c_sin),
]
calc.replace_strings_1("VSI_3ph_6_orig.in", "VSI_3ph_6.in", l)
print('VSI_3ph_6.in is ready for execution')

VSI_3ph_6.in is ready for execution

Execute the following cell to run GSEIM on VSI_3ph_6.in.

import os
import dos_unix
# uncomment for windows:
#dos_unix.d2u("VSI_3ph_6.in")
os.system('run_gseim VSI_3ph_6.in')

get_lib_elements: filename gseim_aux/xbe.aux
get_lib_elements: filename gseim_aux/ebe.aux
Circuit: filename = VSI_3ph_6.in
main: i_solve = 0
main: calling solve_trns
mat_ssw_1_ex: n_statevar: 3
Transient simulation starts...
i=0
i=1000
i=2000
i=3000
i=4000
i=5000
i=6000
i=7000
i=8000
solve_ssw_ex: ssw_iter_newton=0, rhs_ssw_norm=1.2501e+01
Transient simulation starts...
i=0
i=1000
i=2000
i=3000
i=4000
i=5000
i=6000
i=7000
i=8000
solve_ssw_ex: ssw_iter_newton=1, rhs_ssw_norm=8.8402e-12
Transient simulation starts...
i=0
i=1000
i=2000
i=3000
i=4000
i=5000
i=6000
i=7000
i=8000
solve_ssw_ex: ssw_iter_newton=2, rhs_ssw_norm=1.4696e-13
Transient simulation starts...
i=0
i=1000
i=2000
i=3000
i=4000
i=5000
i=6000
i=7000
i=8000
solve_ssw_ex: ssw_iter_newton=3, rhs_ssw_norm=4.9917e-14
solve_ssw_ex: calling solve_ssw_1_ex for one more trns step
Transient simulation starts...
i=0
i=1000
i=2000
i=3000
i=4000
i=5000
i=6000
i=7000
i=8000
solve_ssw_1_ex over (after trns step for output)
solve_ssw_ex ends, slv.ssw_iter_newton=3
GSEIM: Program completed.

0

The circuit file (VSI_3ph_6.in) is created in the same directory as that used for launching Jupyter notebook. The last step (i.e., running GSEIM on VSI_3ph_6.in) creates data files called VSI_3ph_6_1.dat, etc. in the same directory. We can now use the python code below to compute/plot the various quantities of interest.

import numpy as np
import matplotlib.pyplot as plt 
import gseim_calc as calc
from setsize import set_size

slv = calc.slv("VSI_3ph_6.in")

i_slv = 0
i_out = 0
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u1 = np.loadtxt(filename)
t1 = u1[:, 0]

col_sa = slv.get_index(i_slv,i_out,"sa")
col_sb = slv.get_index(i_slv,i_out,"sb")
col_sc = slv.get_index(i_slv,i_out,"sc")
col_t  = slv.get_index(i_slv,i_out,"t" )
col_g1 = slv.get_index(i_slv,i_out,"g1")
col_g2 = slv.get_index(i_slv,i_out,"g2")
col_g3 = slv.get_index(i_slv,i_out,"g3")
col_g4 = slv.get_index(i_slv,i_out,"g4")
col_g5 = slv.get_index(i_slv,i_out,"g5")
col_g6 = slv.get_index(i_slv,i_out,"g6")

# since we have stored two cycles, we need to divide the last time point
# by 2 to get the period:

T = t1[-1]/2

i_out = 1
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u2 = np.loadtxt(filename)
t2 = u2[:, 0]

col_v_ab = slv.get_index(i_slv,i_out,"v_ab")
col_v_bc = slv.get_index(i_slv,i_out,"v_bc")
col_v_ca = slv.get_index(i_slv,i_out,"v_ca")
col_v_an = slv.get_index(i_slv,i_out,"v_an")
col_v_bn = slv.get_index(i_slv,i_out,"v_bn")
col_v_cn = slv.get_index(i_slv,i_out,"v_cn")

i_out = 2
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u3 = np.loadtxt(filename)
t3 = u3[:, 0]

col_ILa  = slv.get_index(i_slv,i_out,"ILa")
col_ILb  = slv.get_index(i_slv,i_out,"ILb")
col_ILc  = slv.get_index(i_slv,i_out,"ILc")
col_IS   = slv.get_index(i_slv,i_out,"IS")
col_P_ea = slv.get_index(i_slv,i_out,"P_ea")
col_P_eb = slv.get_index(i_slv,i_out,"P_eb")
col_P_ec = slv.get_index(i_slv,i_out,"P_ec")
col_P_Ra = slv.get_index(i_slv,i_out,"P_Ra")
col_P_Rb = slv.get_index(i_slv,i_out,"P_Rb")
col_P_Rc = slv.get_index(i_slv,i_out,"P_Rc")

l_IS   = calc.avg_rms_2(t3, u3[:,col_IS  ], T, 2.0*T, 1.0e-4*T)
l_P_ea = calc.avg_rms_2(t3, u3[:,col_P_ea], T, 2.0*T, 1.0e-4*T)
l_P_Ra = calc.avg_rms_2(t3, u3[:,col_P_Ra], T, 2.0*T, 1.0e-4*T)

print('average power delivered to 3-phase source:', "%11.4E"%(-3.0*l_P_ea[1][0]))
print('average power delivered to 3-phase load:', "%11.4E"%(3.0*l_P_Ra[1][0]))
print('average DC supply current:', "%11.4E"%l_IS [1][0])

color1='green'
color2='crimson'
color3='cornflowerblue'
color4='goldenrod'
color5='blue'

fig, ax = plt.subplots(3, sharex=False)
plt.subplots_adjust(wspace=0, hspace=0.0)

set_size(5.5, 6, ax[0])

for i in range(3):
    ax[i].set_xlim(left=T*1e3, right=2.0*T*1e3)
    ax[i].grid(color='#CCCCCC', linestyle='solid', linewidth=0.5)

ax[0].plot(t1*1e3, u1[:,col_t ], color=color4, linewidth=1.0, label="$t$")
ax[0].plot(t1*1e3, u1[:,col_sa], color=color1, linewidth=1.0, label="$sa$")
ax[0].plot(t1*1e3, u1[:,col_sb], color=color2, linewidth=1.0, label="$sb$")
ax[0].plot(t1*1e3, u1[:,col_sc], color=color3, linewidth=1.0, label="$sc$")

ax[1].plot(t2*1e3, u2[:,col_v_an], color=color1, linewidth=0.7, label="$V_{an}$")

ax[2].plot(t2*1e3, u2[:,col_v_ab], color=color2, linewidth=0.7, label="$V_{ab}$")

ax[2].set_xlabel('time (msec)', fontsize=12)

for i in range(3):
    ax[i].legend(loc = 'lower right',frameon = True, fontsize = 10, title = None,
      markerfirst = True, markerscale = 1.0, labelspacing = 0.5, columnspacing = 2.0,
      prop = {'size' : 12},)

plt.tight_layout()
plt.show()

filename: VSI_3ph_6_1.dat
filename: VSI_3ph_6_2.dat
filename: VSI_3ph_6_3.dat
average power delivered to 3-phase source:  2.5049E+04
average power delivered to 3-phase load:  4.3052E+02
average DC supply current:  3.1904E+01

import numpy as np
import matplotlib.pyplot as plt 
import gseim_calc as calc
from setsize import set_size

slv = calc.slv("VSI_3ph_6.in")

i_slv = 0
i_out = 0
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u1 = np.loadtxt(filename)
t1 = u1[:, 0]

col_sa = slv.get_index(i_slv,i_out,"sa")
col_sb = slv.get_index(i_slv,i_out,"sb")
col_sc = slv.get_index(i_slv,i_out,"sc")
col_t  = slv.get_index(i_slv,i_out,"t" )
col_g1 = slv.get_index(i_slv,i_out,"g1")
col_g2 = slv.get_index(i_slv,i_out,"g2")
col_g3 = slv.get_index(i_slv,i_out,"g3")
col_g4 = slv.get_index(i_slv,i_out,"g4")
col_g5 = slv.get_index(i_slv,i_out,"g5")
col_g6 = slv.get_index(i_slv,i_out,"g6")

# since we have stored two cycles, we need to divide the last time point
# by 2 to get the period:

T = t1[-1]/2

i_out = 1
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u2 = np.loadtxt(filename)
t2 = u2[:, 0]

col_v_ab = slv.get_index(i_slv,i_out,"v_ab")
col_v_bc = slv.get_index(i_slv,i_out,"v_bc")
col_v_ca = slv.get_index(i_slv,i_out,"v_ca")
col_v_an = slv.get_index(i_slv,i_out,"v_an")
col_v_bn = slv.get_index(i_slv,i_out,"v_bn")
col_v_cn = slv.get_index(i_slv,i_out,"v_cn")

i_out = 2
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u3 = np.loadtxt(filename)
t3 = u3[:, 0]

col_ILa  = slv.get_index(i_slv,i_out,"ILa")
col_ILb  = slv.get_index(i_slv,i_out,"ILb")
col_ILc  = slv.get_index(i_slv,i_out,"ILc")
col_IS   = slv.get_index(i_slv,i_out,"IS")
col_P_ea = slv.get_index(i_slv,i_out,"P_ea")
col_P_eb = slv.get_index(i_slv,i_out,"P_eb")
col_P_ec = slv.get_index(i_slv,i_out,"P_ec")
col_P_Ra = slv.get_index(i_slv,i_out,"P_Ra")
col_P_Rb = slv.get_index(i_slv,i_out,"P_Rb")
col_P_Rc = slv.get_index(i_slv,i_out,"P_Rc")

color1='green'
color2='crimson'
color3='cornflowerblue'
color4='goldenrod'
color5='blue'

fig, ax = plt.subplots(2, sharex=False)
plt.subplots_adjust(wspace=0, hspace=0.0)

set_size(5.5, 5, ax[0])

for i in range(2):
    ax[i].set_xlim(left=T*1e3, right=2.0*T*1e3)
    ax[i].grid(color='#CCCCCC', linestyle='solid', linewidth=0.5)

ax[0].plot(t3*1e3, u3[:,col_ILa], color=color1, linewidth=1.0, label="$i_{La}$")
ax[0].plot(t3*1e3, u3[:,col_ILb], color=color2, linewidth=1.0, label="$i_{Lb}$")
ax[0].plot(t3*1e3, u3[:,col_ILc], color=color3, linewidth=1.0, label="$i_{Lc}$")

ax[1].plot(t3*1e3, u3[:,col_IS ], color=color5, linewidth=0.7, label="$i_{dc}$")

ax[1].set_xlabel('time (msec)', fontsize=12)

for i in range(2):
    ax[i].legend(loc = 'lower right',frameon = True, fontsize = 10, title = None,
      markerfirst = True, markerscale = 1.0, labelspacing = 0.5, columnspacing = 2.0,
      prop = {'size' : 12},)

plt.tight_layout()
plt.show()

filename: VSI_3ph_6_1.dat
filename: VSI_3ph_6_2.dat
filename: VSI_3ph_6_3.dat

import math
import numpy as np
import matplotlib.pyplot as plt 
import gseim_calc as calc
from setsize import set_size

slv = calc.slv("VSI_3ph_6.in")

i_slv = 0
i_out = 0
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u1 = np.loadtxt(filename)
t1 = u1[:, 0]

col_sa = slv.get_index(i_slv,i_out,"sa")
col_sb = slv.get_index(i_slv,i_out,"sb")
col_sc = slv.get_index(i_slv,i_out,"sc")
col_t  = slv.get_index(i_slv,i_out,"t" )
col_g1 = slv.get_index(i_slv,i_out,"g1")
col_g2 = slv.get_index(i_slv,i_out,"g2")
col_g3 = slv.get_index(i_slv,i_out,"g3")
col_g4 = slv.get_index(i_slv,i_out,"g4")
col_g5 = slv.get_index(i_slv,i_out,"g5")
col_g6 = slv.get_index(i_slv,i_out,"g6")

# since we have stored two cycles, we need to divide the last time point
# by 2 to get the period:

T = t1[-1]/2

i_out = 1
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u2 = np.loadtxt(filename)
t2 = u2[:, 0]

col_v_ab = slv.get_index(i_slv,i_out,"v_ab")
col_v_bc = slv.get_index(i_slv,i_out,"v_bc")
col_v_ca = slv.get_index(i_slv,i_out,"v_ca")
col_v_an = slv.get_index(i_slv,i_out,"v_an")
col_v_bn = slv.get_index(i_slv,i_out,"v_bn")
col_v_cn = slv.get_index(i_slv,i_out,"v_cn")

i_out = 2
filename = slv.l_filename_all[i_slv][i_out]
print('filename:', filename)
u3 = np.loadtxt(filename)
t3 = u3[:, 0]

col_ILa  = slv.get_index(i_slv,i_out,"ILa")
col_ILb  = slv.get_index(i_slv,i_out,"ILb")
col_ILc  = slv.get_index(i_slv,i_out,"ILc")
col_IS   = slv.get_index(i_slv,i_out,"IS")
col_P_ea = slv.get_index(i_slv,i_out,"P_ea")
col_P_eb = slv.get_index(i_slv,i_out,"P_eb")
col_P_ec = slv.get_index(i_slv,i_out,"P_ec")
col_P_Ra = slv.get_index(i_slv,i_out,"P_Ra")
col_P_Rb = slv.get_index(i_slv,i_out,"P_Rb")
col_P_Rc = slv.get_index(i_slv,i_out,"P_Rc")

# compute Fourier coeffs:

t_start = T 
t_end = 2.0*T

n_fourier = 50

coeff_v_an, thd_v_an = calc.fourier_coeff_1C(t2, u2[:,col_v_an], 
    t_start, t_end, 1.0e-8, n_fourier)

coeff_IS, thd_IS = calc.fourier_coeff_1C(t3, u3[:,col_IS], 
    t_start, t_end, 1.0e-8, n_fourier)

coeff_ILa, thd_ILa, coeff_a_ILa, coeff_b_ILa = calc.fourier_coeff_2A(
  t3, u3[:,col_ILa], t_start, t_end, 1.0e-8, n_fourier)

theta1 = math.atan2(-coeff_b_ILa[1], coeff_a_ILa[1])*180.0/math.pi

# if ILa was written as k * sin (w*t + theta), what would theta be?
theta = theta1 + 90.0
print('phase of ILa w.r.t. V_ea:', "%11.4E"%theta)

print("I_La fundamental: RMS value: ", "%11.4E"%(coeff_ILa[1]/np.sqrt(2.0)))
print("V_an fundamental: RMS value: ", "%11.4E"%(coeff_v_an[1]/np.sqrt(2.0)))

x = np.linspace(0, n_fourier, n_fourier+1)

y_v_an  = np.array(coeff_v_an)
y_IS    = np.array(coeff_IS  )
y_ILa   = np.array(coeff_ILa )

fig, ax = plt.subplots(3, sharex=False)
plt.subplots_adjust(wspace=0, hspace=0.0)
grid_color='#CCCCCC'

set_size(6, 5, ax[0])

delta = 10.0
x_major_ticks = np.arange(0.0, (float(n_fourier+1)), delta)
x_minor_ticks = np.arange(0.0, (float(n_fourier+1)), 1.0)

for i in range(3):
    ax[i].set_xlim(left=-1.0, right=float(n_fourier))
    ax[i].set_xticks(x_major_ticks)
    ax[i].set_xticks(x_minor_ticks, minor=True)
    ax[i].grid(visible=True, which='major', axis='x', color=grid_color, linestyle='-', zorder=0)

ax[0].set_ylabel('$V_{an}$',fontsize=14)
ax[1].set_ylabel('$i_{dc}$',fontsize=14)
ax[2].set_ylabel('$i_{La}$',fontsize=14)

ax[2].set_xlabel('N', fontsize=14)

bars1 = ax[0].bar(x, y_v_an, width=0.3, color='red',   label="$V_{an}$", zorder=3)
bars2 = ax[1].bar(x, y_IS  , width=0.3, color='green', label="$i_{dc}$", zorder=3)
bars3 = ax[2].bar(x, y_ILa , width=0.3, color='blue',  label="$i_{La}$", zorder=3)

plt.tight_layout()
plt.show()

filename: VSI_3ph_6_1.dat
filename: VSI_3ph_6_2.dat
filename: VSI_3ph_6_3.dat
phase of ILa w.r.t. V_ea:  1.6386E+01
I_La fundamental: RMS value:   3.7871E+01
V_an fundamental: RMS value:   2.3080E+02

This notebook was contributed by Prof. Nakul Narayanan K, Govt. Engineering College, Thrissur. He may be contacted at nakul@gectcr.ac.in.