MSP Simulation Example - Resistive Companion

Sample Circuit

$R_1$=$1 \Omega$, $R_2$=$2 \Omega$
$L_1$=$20 mH$, $L_2$=$100 mH$, $L_3$=$50 mH$
$V_{in}(t)$=$10 V sin(\omega t)$

Circuit and Simulation Setup

In [1]:
import numpy as np
import ipywidgets as widget
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings('ignore')
%matplotlib inline
np.set_printoptions(sign=' ')

# Circuit parameters
R1 = 1
R2 = 2
L1 = 20e-3
L2 = 100e-3
L3 = 50e-3
Emax = 10

G1 = 1/R1
G2 = 1/R2

#default start values for Ts and T_total
Ts=1e-7
T_total=0.05

Euler Backward Integration Method

In [2]:
#function to perform the euler backward integration 
def int_EB(Ts, T_total, npoint):
    
    GL1_back = Ts/(L1)
    GL2_back = Ts/(L2)
    GL3_back = Ts/(L3)

    # Conductance matrix
    Gn = np.array([ [G1+GL1_back,       -GL1_back,                      0],
                    [-GL1_back,     GL1_back+GL2_back+GL3_back,       -GL3_back],
                    [0,                -GL3_back,                  GL3_back+G2]])

    # Node voltage vector
    vn_back = np.zeros((3,npoint))

    # Current vector in the 3 inductances
    in_back = np.zeros((3,npoint))

    # Intial conditions
    # Voltage source at E(t=0)=0;
    # i_L(t=0)=0  -> e1 = E(t=0) =0; e2 = 0; e3 = 0;
    in_back[:,0] = np.zeros(3)
    vn_back[:,0] = np.zeros(3) 

    # Enter Loop
    # tic %Start a timer
    for i in np.arange(1,npoint):
        # Update source vector
        AL1_back = in_back[0,i-1]
        AL2_back = in_back[1,i-1]
        AL3_back = in_back[2,i-1]
        E = Emax*np.sin(2*np.pi*60*i*Ts)
        Jn = np.array([E/R1-AL1_back, AL1_back-AL2_back-AL3_back, AL3_back])

        # Matrix inversion and solution of the equation G*e=J;
        vn_back[:,i] = np.linalg.solve(Gn, Jn)

        # post step
        in_back[0,i] = AL1_back+GL1_back*(vn_back[0,i]-vn_back[1,i])
        in_back[1,i] = AL2_back+GL2_back*(vn_back[1,i])
        in_back[2,i] = AL3_back+GL3_back*(vn_back[1,i]-vn_back[2,i])

    # toc %stop the timer 

    
    return in_back

Euler Forward Integration Method

In [3]:
#function to perform the euler forward integration 
def int_EF(Ts, T_total, npoint):
    
    # Large resistor added in parallel with inductors (3 current sources coincide at a node)
    Gadd = 1e-4
    # Conductance matrix
    Gn = np.array(  [[G1+Gadd,       -Gadd,           0],
                    [-Gadd,         3*Gadd,        -Gadd],
                    [   0,           -Gadd,       Gadd+G2]])

    # Node voltage vector
    vn_forw = np.zeros((3,npoint))

    # Current vector in the 3 inductances
    in_forw = np.zeros((3,npoint))

    # Intial conditions
    # Voltage source at E(t=0)=0;
    # i_L(t=0)=0  -> e1 = E(t=0) =0; e2 = 0; e3 = 0;
    in_forw[:,0] = np.zeros(3)
    vn_forw[:,0] = np.zeros(3) 

    #Enter Loop
    # tic % Start a timer
    for i in np.arange(1,npoint):
        # Update source vector
        AL1_forw = in_forw[0,i-1]+(vn_forw[0,i-1]-vn_forw[1,i-1])*Ts/L1
        AL2_forw = in_forw[1,i-1]+(vn_forw[1,i-1])*Ts/L2
        AL3_forw = in_forw[2,i-1]+(vn_forw[1,i-1]-vn_forw[2,i-1])*Ts/L3
        E = Emax*np.sin(2*np.pi*60*i*Ts)
        Jn = np.array([E/R1-AL1_forw, AL1_forw-AL2_forw-AL3_forw, AL3_forw])

        # Matrix inversion and solution of the equation G*e=J 
        vn_forw[:,i] = np.linalg.solve(Gn, Jn)

        # post step
        in_forw[0,i]= AL1_forw
        in_forw[1,i]= AL2_forw
        in_forw[2,i]= AL3_forw
        
    
    return in_forw

Trapezoidal Integration Method

In [4]:
#function to perform the trapezoidal integration 
def int_TR(Ts, T_total, npoint):
    
    GL1_trap = Ts/(2*L1)
    GL2_trap = Ts/(2*L2)
    GL3_trap = Ts/(2*L3)

    # Conductance matrix
    Gn = np.array(  [[G1+GL1_trap,           -GL1_trap,                         0],
                     [-GL1_trap,     GL1_trap+GL2_trap+GL3_trap,           -GL3_trap],
                     [0,                    -GL3_trap,                   GL3_trap+G2]])

    # Node voltage vector
    vn_trap = np.zeros((3,npoint))

    # Current vector in the 3 inductances
    in_trap = np.zeros((3,npoint))

    # Intial conditions
    # Voltage source at E(t=0)=0
    # i_L(t=0)=0  -> e1 = E(t=0) =0; e2 = 0; e3 = 0;
    in_trap[:,0] = np.zeros(3)
    vn_trap[:,0] = np.zeros(3)

    #Enter Loop
    #tic
    for i in np.arange(1,npoint):
        # Update source vector
        AL1 = in_trap[0,i-1]+(vn_trap[0,i-1]-vn_trap[1,i-1])*Ts/(2*L1)
        AL2 = in_trap[1,i-1]+(vn_trap[1,i-1])*Ts/(2*L2)
        AL3 = in_trap[2,i-1]+(vn_trap[1,i-1]-vn_trap[2,i-1])*Ts/(2*L3)
        E = Emax*np.sin(2*np.pi*60*i*Ts)
        Jn = np.array([E/R1-AL1, AL1-AL2-AL3, AL3])

        # Matrix inversion and solution of the equation G*e=J
        vn_trap[:,i] = np.linalg.solve(Gn, Jn)

        # post step
        in_trap[0,i]= AL1+GL1_trap*(vn_trap[0,i]-vn_trap[1,i])
        in_trap[1,i]= AL2+GL2_trap*(vn_trap[1,i])
        in_trap[2,i]= AL3+GL3_trap*(vn_trap[1,i]-vn_trap[2,i])
    #toc
    
    return in_trap

Result Comparison for different Integration Methods

In [5]:
#Plot function
def plot_fun(time_step, npoint):
    
    #Number of subplots (1*3) and size of the figure where the subplots will appear
    fig, ax = plt.subplots(1,3, figsize=(16, 5))
    
    #Calculate time values
    t = np.arange(0, npoint)*Ts
    
    #Transform time values to milliseconds for better representation
    t= np.multiply(t, 10e+3)
              
    #Plot results E. Backward
    ax[0].set_title('E. Backward')
    ax[0].set_xlim([0, (npoint-1)*Ts*10e+3])
    ax[0].set_xlabel('Time [ms]')
    ax[0].set_ylabel('Current [A]')
    ax[0].plot(t,in_back[0,:], 'r', t, in_back[2,:], 'r:')
    leg = ax[0].legend(['$I_{L1}$(t)','$I_{L2}$(t)'])
    
    #Plot results E. Forward
    ax[1].set_title('E. Forward')
    #Set axis limite and transform to milliseconds
    ax[1].set_xlim([0, (npoint-1)*Ts*10e+3])
    ax[1].set_xlabel('Time [ms]')
    ax[1].set_ylabel('Current [A]')
    ax[1].plot(t, in_forw[0,:] ,'g', t, in_forw[2,:], 'g:', linewidth=2)
    leg = ax[1].legend(['$I_{L1}$(t)','$I_{L2}$(t)'])
        
    #Plot results Trapezoidal
    ax[2].set_title('Trapezoidal')
    ax[2].set_xlim([0, (npoint-1)*Ts*10e+3])
    ax[2].set_xlabel('Time [ms]')
    ax[2].set_ylabel('Current [A]')
    ax[2].plot(t,in_trap[0,:], 'b', t, in_trap[2,:], 'b:', linewidth=2)
    leg = ax[2].legend(['$I_{L1}$(t)','$I_{L2}$(t)'])
    
    
#Slider function to forward the selected values to the simulation functions
def slider_fun(time_step):
    
    global in_back , in_forw , in_trap, npoint 
    
    if time_step==0 or T_total==0:
        print('\n \n \n' + 'Simulation time step is set to: ' + str(time_step) + ' seconds')
        print('Total simulation time: ' + str(T_total) + ' seconds' + '\n \n')
    
    else:
        print('\n \n \n' + 'Simulation time step is set to: ' + str(time_step) + ' seconds')
        print('Total simulation time: ' + str(T_total) + ' seconds' + '\n \n')
    
        npoint = int(np.round(T_total/time_step))

        in_back=int_EB(time_step, T_total, npoint)

        in_forw=int_EF(time_step, T_total, npoint)

        in_trap=int_TR(time_step, T_total, npoint)
        
        plot_fun(time_step, npoint)

    return time_step, T_total

#Values range for the simulation steps
values=[round(i*10**-7, 7) for i in range(101)]

#It is important to set the argument continous_update to 'false' so that the slider_fun is called after the user finishes dragging the slider. 
#Otherwise the function will be continuously called for several values in the dragging interval and the expected results will take alot of time to appear   
slider= widget.interactive(slider_fun, time_step=widget.SelectionSlider(description="Time step", continuous_update=False, options=[("%g"%i,i) for i in values], layout=widget.Layout( width='50%', height='50px',)));

display(slider)

Critical time step:

$7,6 \cdot 10^{-6} s < T_{critical} $

In [ ]: