from sympy import *
import numpy as np
from tabulate import tabulate
import pandas as pd
import matplotlib.pyplot as plt
import SymMNA
from IPython.display import display, Markdown, Math, Latex
init_printing()10 Initial Conditions
Initial conditions for capacitors and inductors are essential in transient circuit analysis as they represent the energy stored in these components immediately prior to a change in the circuit, such as a switch opening or closing at time \(t=0\).
Because capacitors store energy in an electric field and inductors store energy in a magnetic field, they possess a fundamental property: the voltage across a capacitor cannot change instantaneously, and the current through an inductor cannot change instantaneously.
Mathematically, this is expressed as \(v_C(0^-) = v_C(0^+)\) and \(i_L(0^-) = i_L(0^+)\), where \(0^-\) is the instant just before the change and \(0^+\) is the instant just after.
To find these initial values, one typically analyzes the circuit for \(t<0\) (often assuming it has been in a steady-state for a long time) by replacing the capacitor with an open circuit (since DC current is zero) and the inductor with a short circuit (since DC voltage is zero) and then solving for the capacitor voltage and inductor current at \(t=0^-\).
These continuous values then carry over to \(t=0^+\) and serve as the necessary constants for solving the differential equations that describe the circuit’s transient response for \(t>0\).
10.1 Capacitor Initial Conditions
C_IC_test_jig_1_net_list = '''
* C_IC_test_jig_1.asc
C1 2 1 1e-6
V1 1 0 1
R1 2 0 1e3
'''The MNA equations are generated from the function SymMNA.smna.
report, network_df, i_unk_df, A, X, Z = SymMNA.smna(C_IC_test_jig_1_net_list)The code below assembles the network equations from the MNA matrices and displays the equations.
# Put matrices into SymPy
X = Matrix(X)
Z = Matrix(Z)
# put the matricies into equation form
NE_sym = Eq(A*X,Z)
# free symbols are entered as SymPy variables
var(str(NE_sym.free_symbols).replace('{','').replace('}',''))
# the element values are put into a dictionary
element_values = SymMNA.get_part_values(network_df) # get dummy element vales from netlist
# display the equations
temp = ''
for i in range(shape(NE_sym.lhs)[0]):
temp += '${:s} = {:s}$<br>'.format(latex(NE_sym.rhs[i]),latex(NE_sym.lhs[i]))
Markdown(temp)\(0 = C_{1} s v_{1} - C_{1} s v_{2} + I_{V1}\)
\(0 = - C_{1} s v_{1} + v_{2} \left(C_{1} s + \frac{1}{R_{1}}\right)\)
\(V_{1} = v_{1}\)
SOlve the equations and display the results.
U_sym = solve(NE_sym,X)
temp = ''
for i in U_sym.keys():
temp += '${:s} = {:s}$<br>'.format(latex(i),latex(U_sym[i]))
Markdown(temp)\(v_{1} = V_{1}\)
\(v_{2} = \frac{C_{1} R_{1} V_{1} s}{C_{1} R_{1} s + 1}\)
\(I_{V1} = - \frac{C_{1} V_{1} s}{C_{1} R_{1} s + 1}\)
Initial conditions for the capacitor \(C_1\) is set by making \(V_1\) a step voltage equal to the value of the initial voltage condition. \(1/s\) is substituted for the value of \(V_1\).
t = symbols('t',positive=True) # t > 0Substitute \(1/s\) for \(V_1\)
element_values[V1] = 1/s #laplace_transform(Heaviside(t), t, s)[0]
NE_ic = NE_sym.subs(element_values)
# display the equations
temp = ''
for i in range(shape(NE_ic.lhs)[0]):
temp += '${:s} = {:s}$<br>'.format(latex(NE_ic.rhs[i]),latex(NE_ic.lhs[i]))
Markdown(temp)\(0 = I_{V1} + 1.0 \cdot 10^{-6} s v_{1} - 1.0 \cdot 10^{-6} s v_{2}\)
\(0 = - 1.0 \cdot 10^{-6} s v_{1} + v_{2} \cdot \left(1.0 \cdot 10^{-6} s + 0.001\right)\)
\(\frac{1}{s} = v_{1}\)
Solve the network equations and display the results.
U_ic = solve(NE_ic,X)
temp = ''
for i in U_ic.keys():
temp += '${:s} = {:s}$<br>'.format(latex(i),latex(U_ic[i]))
Markdown(temp)\(v_{1} = \frac{1}{s}\)
\(v_{2} = \frac{1}{s + 1000.0}\)
\(I_{V1} = - \frac{1}{1000.0 s + 1000000.0}\)
The voltage on node 2 is symplified with the chain of operators applied to the expression; nsimplify(), simplify(), expand() and together(). This helps SymPy solve the inverse Laplace transform.
node_2_voltage_s = U_ic[v2].nsimplify().simplify().expand().together()
node_2_voltage_s\(\displaystyle \frac{1}{s + 1000}\)
node_2_voltage_t = inverse_laplace_transform(node_2_voltage_s, s, t)
node_2_voltage_t\(\displaystyle e^{- 1000 t}\)
func_node_2_voltage_t = lambdify(t, node_2_voltage_t) Plot the voltage at node 2.
x_axis = np.linspace(0, 10e-3, 2000, endpoint=True)
plt.title('Node 2 voltage')
plt.plot(x_axis*1e3, func_node_2_voltage_t(x_axis),'-b')
plt.ylabel('v(t), volts')
plt.xlabel('time, msec')
plt.grid()
plt.show()
Plot the current in the capacitor, this is the same as \(I_{V1}\).
V1_current_s = U_ic[I_V1].nsimplify().simplify().expand().together()
V1_current_s\(\displaystyle - \frac{1}{1000 \left(s + 1000\right)}\)
V1_current_t = inverse_laplace_transform(V1_current_s, s, t)
V1_current_t\(\displaystyle - \frac{e^{- 1000 t}}{1000}\)
func_V1_current_t = lambdify(t, V1_current_t) The plot below shows the node voltages versus time.
x_axis = np.linspace(0, 10e-3, 2000, endpoint=True)
plt.title('V1 current')
plt.plot(x_axis*1e3, func_V1_current_t(x_axis),'-b')
plt.ylabel('i(t), amps')
plt.xlabel('time, msec')
plt.grid()
plt.show()
10.2 Inductor Initial Conditions
- need to update values
- swap polarity of I1, change \(L_1\) value to 1
In the LTSpice simulation \(I_1=0\) which since there is no current folowing is an open circuit. The initial conditions for \(L_1\) is set by .ic I(L1)=-1e-3, notice the negative sign.
L_IC_test_jig_1_net_list = '''
* L_IC_test_jig_1.asc
R1 1 0 1e3
L1 1 0 1
I1 0 1 1
'''The MNA equations are generated from the function SymMNA.smna.
report, network_df, i_unk_df, A, X, Z = SymMNA.smna(L_IC_test_jig_1_net_list)The code below assembles the network equations from the MNA matrices and displays the equations.
# Put matrices into SymPy
X = Matrix(X)
Z = Matrix(Z)
# put the matricies into equation form
NE_sym = Eq(A*X,Z)
# free symbols are entered as SymPy variables
var(str(NE_sym.free_symbols).replace('{','').replace('}',''))
# the element values are put into a dictionary
element_values = SymMNA.get_part_values(network_df) # get dummy element vales from netlist
# display the equations
temp = ''
for i in range(shape(NE_sym.lhs)[0]):
temp += '${:s} = {:s}$<br>'.format(latex(NE_sym.rhs[i]),latex(NE_sym.lhs[i]))
Markdown(temp)\(I_{1} = I_{L1} + \frac{v_{1}}{R_{1}}\)
\(0 = - I_{L1} L_{1} s + v_{1}\)
Solve the equations and display the results.
U_sym = solve(NE_sym,X)
temp = ''
for i in U_sym.keys():
temp += '${:s} = {:s}$<br>'.format(latex(i),latex(U_sym[i]))
Markdown(temp)\(v_{1} = \frac{I_{1} L_{1} R_{1} s}{L_{1} s + R_{1}}\)
\(I_{L1} = \frac{I_{1} R_{1}}{L_{1} s + R_{1}}\)
The initial conditions for the inductor is set by the current source \(I_1\). A step function with a value of the step is used. The steady state value, a DC level, of the step is shorted by the inductor and the voltage across the inductor from this contribution is 0.
t = symbols('t',positive=True) # t > 0element_values[I1] = 1/s #laplace_transform(1e-3*Heaviside(t), t, s)[0]
NE_ic = NE_sym.subs(element_values)# display the equations
temp = ''
for i in range(shape(NE_ic.lhs)[0]):
temp += '${:s} = {:s}$<br>'.format(latex(NE_ic.rhs[i]),latex(NE_ic.lhs[i]))
Markdown(temp)\(\frac{1}{s} = I_{L1} + 0.001 v_{1}\)
\(0 = - 1.0 I_{L1} s + v_{1}\)
Solve the network equations and display the results.
U_ic = solve(NE_ic,X)
temp = ''
for i in U_ic.keys():
temp += '${:s} = {:s}$<br>'.format(latex(i),latex(U_ic[i]))
Markdown(temp)\(v_{1} = \frac{1000.0}{s + 1000.0}\)
\(I_{L1} = \frac{1000.0}{s^{2} + 1000.0 s}\)
The voltage on node 2 is symplified with the chain of operators applied to the expression; nsimplify(), simplify(), expand() and together(). This helps SymPy solve the inverse Laplace transform.
temp = U_ic[v1].nsimplify().simplify().expand().together()
temp\(\displaystyle \frac{1000}{s + 1000}\)
node_1_voltage = inverse_laplace_transform(temp, s, t)
node_1_voltage\(\displaystyle 1000 e^{- 1000 t}\)
func_node_1_voltage = lambdify(t, node_1_voltage) The plot below shows the node voltages versus time.
x_axis = np.linspace(0, 10e-3, 2000, endpoint=True)
plt.title('Node 1 voltage')
plt.plot(x_axis*1e3, func_node_1_voltage(x_axis),'-b')
plt.ylabel('v1(t), volts')
plt.xlabel('time, msec')
plt.grid()
plt.show()
Find the currrent in the inductor.
L1_current_s = U_ic[I_L1].nsimplify().simplify().expand().together()
L1_current_s\(\displaystyle \frac{1000}{s \left(s + 1000\right)}\)
L1_current_t = inverse_laplace_transform(L1_current_s, s, t)
L1_current_t\(\displaystyle 1 - e^{- 1000 t}\)
func_L1_current_t = lambdify(t, L1_current_t) The plot below shows inductor current versus time.
x_axis = np.linspace(0, 10e-3, 2000, endpoint=True)
plt.title('Inductor current')
plt.plot(x_axis*1e3, func_L1_current_t(x_axis),'-b')
plt.ylabel('L1 current, amps')
plt.xlabel('time, msec')
plt.grid()
plt.show()
Inductor current should approach zero. Can’t determine the inductor current since it’s being driven by \(I_1\) and \(I_1\) is a step function which has a final value of 1.
10.3 Circuit 2
A circuit with initial conditions consisting of a capacitor with an inital voltage and an inductor with initial current is analysed. The circuit in Figure 10.1 has seven branches and four nodes. Capacitor, C1, has an initial voltage at t(0) and the inductor, L1, has an initial current at t(0). The voltage source V2 represents the initial voltage on the capacitor and current source, I1, represents the initial current flowing in the inductor. The circuit has a DC voltage source, V1. There a current controlled current source, F1, that is controlled by the current in V1.
The circuit was drawn using LTSpice and the netlist was pasted into the code. The Heaviside function is used to define the initial current and voltage on the inductor and capacitor. The Heaviside function is named after Oliver Heaviside, who made significant contributions to electrical engineering.
The following Python modules are used.
t = symbols('t',positive=True) # t > 010.3.1 Load the netlist
The netlist generated by LTSpice is pasted into the cell below and some edits were made to remove the inductor series resistance and the independent sources are set to their DC values.
net_list = '''
V1 1 0 1
R2 3 2 5
C1 3 4 1
V2 4 0 1
L1 1 2 3
I1 1 2 1
F1 2 0 V1 2
'''The MNA equations are generated from the function SymMNA.smna.
report, network_df, i_unk_df, A, X, Z = SymMNA.smna(net_list)The code below assembles the network equations from the MNA matrices and displays the equations.
# Put matrices into SymPy
X = Matrix(X)
Z = Matrix(Z)
# put the matricies into equation form
NE_sym = Eq(A*X,Z)
# free symbols are entered as SymPy variables
var(str(NE_sym.free_symbols).replace('{','').replace('}',''))
# the element values are put into a dictionary
element_values = SymMNA.get_part_values(network_df) # get dummy element vales from netlist
# display the equations
temp = ''
for i in range(shape(NE_sym.lhs)[0]):
temp += '${:s} = {:s}$<br>'.format(latex(NE_sym.rhs[i]),latex(NE_sym.lhs[i]))
Markdown(temp)\(- I_{1} = I_{L1} + I_{V1}\)
\(I_{1} = I_{F1} - I_{L1} + \frac{v_{2}}{R_{2}} - \frac{v_{3}}{R_{2}}\)
\(0 = - C_{1} s v_{4} + v_{3} \left(C_{1} s + \frac{1}{R_{2}}\right) - \frac{v_{2}}{R_{2}}\)
\(0 = - C_{1} s v_{3} + C_{1} s v_{4} + I_{V2}\)
\(V_{1} = v_{1}\)
\(V_{2} = v_{4}\)
\(0 = - I_{L1} L_{1} s + v_{1} - v_{2}\)
\(0 = I_{F1} - I_{V1} f_{1}\)
As shown above MNA generated many equations and these would be difficult to solve by hand and a symbolic soultion would take a lot of computing time. The equations are displace in matrix notation.
The symbols generated by the Python code are extraced by the SymPy function free_symbols and then declared as SymPy variables.
The SymPy Heaviside function is used to define the initial current and voltage on the inductor and capacitor.
element_values[V1] = laplace_transform(Heaviside(t), t, s)[0]
element_values[V2] = laplace_transform(-0.2*Heaviside(t), t, s)[0]
element_values[I1] = laplace_transform(0.1*Heaviside(t), t, s)[0]
NE_ic = NE_sym.subs(element_values)
# display the equations
temp = ''
for i in range(shape(NE_ic.lhs)[0]):
temp += '${:s} = {:s}$<br>'.format(latex(NE_ic.rhs[i]),latex(NE_ic.lhs[i]))
Markdown(temp)\(- \frac{0.1}{s} = I_{L1} + I_{V1}\)
\(\frac{0.1}{s} = I_{F1} - I_{L1} + 0.2 v_{2} - 0.2 v_{3}\)
\(0 = - 1.0 s v_{4} - 0.2 v_{2} + v_{3} \cdot \left(1.0 s + 0.2\right)\)
\(0 = I_{V2} - 1.0 s v_{3} + 1.0 s v_{4}\)
\(\frac{1}{s} = v_{1}\)
\(- \frac{0.2}{s} = v_{4}\)
\(0 = - 3.0 I_{L1} s + v_{1} - v_{2}\)
\(0 = I_{F1} - 2.0 I_{V1}\)
Solve the network equations and display the results.
U_ic = solve(NE_ic,X)
temp = ''
for i in U_ic.keys():
temp += '${:s} = {:s}$<br>'.format(latex(i),latex(U_ic[i]))
Markdown(temp)\(v_{1} = \frac{1}{s}\)
\(v_{2} = \frac{13.0 s^{2} + 53.0 s + 10.0}{10.0 s^{3} + 50.0 s^{2} + 10.0 s}\)
\(v_{3} = \frac{- 2.0 s^{2} - 7.0 s + 10.0}{10.0 s^{3} + 50.0 s^{2} + 10.0 s}\)
\(v_{4} = - \frac{0.2}{s}\)
\(I_{V1} = \frac{- s - 4.0}{10.0 s^{2} + 50.0 s + 10.0}\)
\(I_{V2} = \frac{3.0 s + 12.0}{10.0 s^{2} + 50.0 s + 10.0}\)
\(I_{L1} = \frac{- s - 1.0}{10.0 s^{3} + 50.0 s^{2} + 10.0 s}\)
\(I_{F1} = \frac{- s - 4.0}{5.0 s^{2} + 25.0 s + 5.0}\)
The voltage on node 2 is symplified with the chain of operators applied to the expression; nsimplify(), simplify(), expand() and together(). This helps SymPy solve the inverse Laplace transform.
temp = U_ic[v2].nsimplify().simplify().expand().together()
temp\(\displaystyle \frac{13 s^{2} + 53 s + 10}{10 s \left(s^{2} + 5 s + 1\right)}\)
node_2_voltage = inverse_laplace_transform(temp, s, t)
node_2_voltage\(\displaystyle \left(\frac{\sqrt{21} \cdot \left(3 e^{\frac{t \left(5 - \sqrt{21}\right)}{2}} + \sqrt{21} e^{\frac{t \left(5 - \sqrt{21}\right)}{2}} - 3 e^{\frac{t \left(\sqrt{21} + 5\right)}{2}} + \sqrt{21} e^{\frac{t \left(\sqrt{21} + 5\right)}{2}}\right)}{140} + e^{5 t}\right) e^{- 5 t}\)
func_node_2_voltage = lambdify(t, node_2_voltage) The voltage on node 3 is obtained in a sumular way.
temp = U_ic[v3].nsimplify().simplify().expand().together()
temp\(\displaystyle \frac{- 2 s^{2} - 7 s + 10}{10 s \left(s^{2} + 5 s + 1\right)}\)
node_3_voltage = inverse_laplace_transform(temp, s, t)
node_3_voltage\(\displaystyle \left(\frac{\sqrt{21} \left(- 2 \sqrt{21} e^{\frac{t \left(5 - \sqrt{21}\right)}{2}} + 9 e^{\frac{t \left(5 - \sqrt{21}\right)}{2}} - 2 \sqrt{21} e^{\frac{t \left(\sqrt{21} + 5\right)}{2}} - 9 e^{\frac{t \left(\sqrt{21} + 5\right)}{2}}\right)}{70} + e^{5 t}\right) e^{- 5 t}\)
func_node_3_voltage = lambdify(t, node_3_voltage) The plot below shows the node voltages versus time.
x = np.linspace(0, 10, 2000, endpoint=True)
plt.title('Node voltages vs time')
plt.plot(x, func_node_2_voltage(x),label='v2(t)')
plt.plot(x, func_node_3_voltage(x),label='v3(t)')
plt.ylabel('v(t), volts')
plt.xlabel('time, sec')
plt.legend()
plt.grid()
plt.show()
10.4 Summary
In this chapter, newtork equations were solved which had initial conditions.
Note
Still investigating the use of a step function for the inductor’s initial conditions since the results don’t look correct.