$\Delta t=0.1ms$, $I_1=1A$, $E_{23}=1V$
$L_{12}=L_{23}=0.5mH$
$C_1=C_2=C_3=0.5mF$
$R_1=R_2=R_3=R_{12}=R_{23}=1\Omega $
$V_1(0)=V_2(0)=V_3(0)=0$
$I_{12}(0)=I_{23}(0)=0$
import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(sign=' ')
# Circuit parameters:
R = 1.0
L = 0.5e-3
C = 0.5e-3
I1 = 1.0
E23 = 1
# Incidence matrix
M = np.array([ [1, 0],
[-1, 1],
[0, -1]])
# Number of branches Nb and number of nodes Nn:
(Nn, Nb) = np.shape(M)
# Branches
R_b = np.array([ [R, R]])
L_b = np.array([ [L, L]])
# Nodes
R_n = np.array([ [R, R, R]])
C_n = np.array([ [C, C, C]])
# Total simulation time
T_total = 0.008
# Simulaiton time step
Ts = 0.1e-3
# Number of simulation time steps
npoint = int(np.round(T_total/Ts))
print('Total simulation time: ' + str(T_total))
print('\nSimulation time step: ' + str(Ts))
# Current source vectors for all nodes
H = np.zeros((npoint, Nn))
H[:,0] = I1*np.ones((npoint))
# Voltage source vectors for all branches
E = np.zeros((npoint, Nb))
E[:,1] = E23*np.ones((npoint))
# Create matrix to store solution of branch currents in every simulation steps
I = np.zeros((npoint, Nb))
# Create matrix to store solution of node voltages in every simulation steps
V = np.zeros((npoint, Nn))
# Initial conditions
# Branch inductances: current = 0
I[0,:] = np.zeros((1, Nb))
# Node capacitance: voltage = 0
# Calculate V(Ts/2) based on a given initial conditions V(0)
V_0 = np.zeros((1, Nn))
for j in np.arange(0, Nn):
V[0,j] = (C_n[0,j]*V_0[0,j]/Ts - np.matmul(M[j,:], np.transpose(I[0,:])) + H[0,j]) / (C_n[0,j]/Ts + 1/R_n[0,j])
# Time loop
for i in np.arange(1,npoint):
# Branch loop - calculation of I(i*Ts) based on I((i-1)*Ts) and V((i-1/2)*Ts)
for j in np.arange(0, Nb):
I[i,j] = I[i-1,j] + Ts/L_b[0,j] * (np.matmul(np.transpose(M[:,j]), np.transpose(V[i-1,:])) - R_b[0,j] * I[i-1,j] + E[i-1,j])
# Node loop - calculation of V((i+1/2)*Ts) based on I(i*Ts) and V((i-1/2)*Ts)
for j in np.arange(0, Nn):
V[i,j] = (C_n[0,j]*V[i-1,j]/Ts - np.matmul(M[j,:], np.transpose(I[i,:])) + H[i,j]) / (C_n[0,j]/Ts + 1/R_n[0,j])
print('Initialization:')
print('V1(0.05ms) = ' + str(V[0,0]))
print('V2(0.05ms) = ' + str(V[0,1]))
print('V3(0.05ms) = ' + str(V[0,2]) + '\n')
print('Branch currents at t=0.1ms:')
print('I12(0.1ms) = ' + str(I[1,0]))
print('I23(0.1ms) = ' + str(I[1,1]) + '\n')
print('Node voltages at t=0.15ms:')
print('V1(0.15ms) = ' + str(V[1,0]))
print('V2(0.15ms) = ' + str(V[1,1]))
print('V3(0.15ms) = ' + str(V[1,2]) + '\n')
print('Branch currents at t=0.2ms:')
print('I12(0.2ms) = ' + str(I[2,0]))
print('I23(0.2ms) = ' + str(I[2,1]) + '\n')
# Time vectors
# for branch currents
t_I = np.arange(0, npoint)*Ts
# for node voltages
t_V = (np.arange(0, npoint) + 0.5)*Ts
# Plots
plt.figure(figsize=(8,6))
plt.xlabel('Time [s]')
plt.ylabel('Current [A]')
plt.title('Branch currents')
plt.plot(t_I,I[:,0], t_I, I[:,1], linewidth=2)
plt.xlim([0, (npoint-1)*Ts])
plt.legend(['Current $i_{12}$(t)', 'Current $i_{23}$(t)'], prop={'size': 14})
plt.show()
# Plots
plt.figure(figsize=(8,6))
plt.xlabel('Time [s]')
plt.ylabel('Voltage [V]')
plt.title('Node voltages')
plt.plot(t_V,V[:,0], t_V, V[:,1], t_V,V[:,2], linewidth=2)
plt.xlim([0, (npoint-1)*Ts])
plt.legend(['Voltage $v_{1}$(t)', 'Voltage $v_{2}$(t)', 'Voltage $v_{3}$(t)'], prop={'size': 14})
plt.show()