Here is a code for a simple model of three phase two level inverter with constant DC voltage source and three phase RL load.

import numpy as np
from datetime import datetime
from scipy.integrate import solve_ivp
from scipy import signal

#time axis
starting_time = 0 #[s]
step_time = .2e-6 #[s]
ending_time = 50e-3 #[s]
#PWM modulator
m = 0.85 #[-] modulation index
f_e = 1e3 #[Hz] electrical frequency
f_sw = 40e3 #[Hz] switching frequency
Vdc = 800 #[V] DC-link voltage
R = 0.5 #[ohm]
L = 400e-6 #[H] self-inductance
K = 0 #inductances coupling coefficient, <1
M = K*L #[H] mutual inductance


def ref(time, amplitude, frequency, phase):
    ref = amplitude*np.sin(2*np.pi*frequency*time + phase)
    return ref

def triangle(time, frequency):
    triangle = signal.sawtooth(2*np.pi*frequency*time, 0.5)
    return triangle

def SPWM(mod, carrier):
    if mod >= carrier:
        d = 1
    if mod < carrier:
        d = 0
    return d

def system(t, y):
    mod1 = ref(t, m, f_e, 0)
    mod2 = ref(t, m, f_e, 2/3*np.pi)
    mod3 = ref(t, m, f_e, 4/3*np.pi)
    carrier = triangle(t, f_sw)
    d1 = SPWM(mod1, carrier)
    d2 = SPWM(mod2, carrier)
    d3 = SPWM(mod3, carrier)
    v1_gnd = d1*Vdc
    v2_gnd = d2*Vdc
    v3_gnd = d3*Vdc
    vgnd_n = -(v1_gnd + v2_gnd + v3_gnd)/3
    v1 = v1_gnd + vgnd_n 
    v2 = v2_gnd + vgnd_n 
    v3 = v3_gnd + vgnd_n 
    v12 = v1 - v2 
    v23 = v2 - v3 
    v13 = v1 - v3
    i1 = 1/(L*L+L*M-2*M*M)*((L+M)*y[0]-M*y[1]-M*y[2])
    i2 = 1/(L*L+L*M-2*M*M)*((L+M)*y[1]-M*y[0]-M*y[2])
    i3 = 1/(L*L+L*M-2*M*M)*((L+M)*y[2]-M*y[0]-M*y[1])
    dl1dt = -R*i1+v1
    dl2dt = -R*i2+v2
    dl3dt = -R*i3+v3

    return dl1dt, dl2dt, dl3dt

dt1 = datetime.now()

time = np.arange(starting_time,ending_time,step_time)
time = np.round(time,9)

sol = solve_ivp(system, [0,ending_time], [0,0,0], t_eval = time, method = 'LSODA', max_step = 1/50/f_sw, rtol = 1e-6, atol = 1e-6)
t = sol.t
l1 = sol.y[0]
l2 = sol.y[1]
l3 = sol.y[2]
i1 = 1/(L*L+L*M-2*M*M)*((L+M)*l1-M*l2-M*l3)
i2 = 1/(L*L+L*M-2*M*M)*((L+M)*l2-M*l1-M*l3)
i3 = 1/(L*L+L*M-2*M*M)*((L+M)*l3-M*l1-M*l2)

dt2 = datetime.now()

I need to consider an high switching frequency (max 40kHz), so sampling should be around 1/40e3/50 to get good results.

At the moment, with my PC I'm only able to simulate 50ms in around 1 min. I was wondering if there's a way to optimize this problem and speed up the simulation. Ideally I would like to be able to simulate minutes of inverter operation in reasonable time.

I tested that using a for loop with fixed time step and explicit integration method speeds up the simulation (not that much). Is there something better?



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.