System Dynamics
System Dynamics
Home

Impedance Modeling in Python

Systematically solving for transfer functions using impedance methods is straightforward with Python. The following example illustrates the process.

Example

Consider the linear graph of an electronic system, below. Use impedance methods to derive the transfer functions from inputs VS and IS to outputs vC2 and iR1.

Linear graph.
Figure : Linear graph.

The normal tree is shown, below, along with contours to be used for continuity equations.

Normal tree.
Figure : Normal tree.

Let’s continue in Python. Load the SymPy package as follows:

import sympy as sp

First, define the symbolic variables as follows:

s = sp.symbols("s", complex=True)
R1, R2, R3, L1, L2, C1, C2 = sp.symbols(
    "R1, R2, R3, L1, L2, C1, C2", positive=True
)
VR1, VR2, VR3, VL1, VL2, VC1, VC2, \
IR1, IR2, IR3, IL1, IL2, IC1, IC2 = sp.symbols(
    "VR1, VR2, VR3, VL1, VL2, VC1, VC2," \
    "IR1, IR2, IR3, IL1, IL2, IC1, IC2", complex=True
)
VS, IS = sp.symbols("VS, IS", complex=True)
ZR1, ZR2, ZR3, ZL1, ZL2, ZC1, ZC2 = sp.symbols(
    "ZR1, ZR2, ZR3, ZL1, ZL2, ZC1, ZC2", complex=True
)

Define the unknown variables and the input and output variables for the requested transfer functions as follows:

vars_input = [VS, IS]
vars_unknown = [
    VR1, VR2, VR3, VL1, VL2, VC1, VC2,
	IR1, IR2, IR3, IL1, IL2, IC1, IC2
]
input_indices = [0, 1]  # VS, IS
output_indices = [6, 7]  # VC2, IR1

Write the elemental equations as follows:

eqs_ele = [
    VR1 - IR1 * ZR1,  # == 0
    VR2 - IR2 * ZR2,  # == 0
    VR3 - IR3 * ZR3,  # == 0
    VL1 - IL1 * ZL1,  # == 0
    VL2 - IL2 * ZL2,  # == 0
    VC1 - IC1 * ZC1,  # == 0
    VC2 - IC2 * ZC2,  # == 0
]

Write the continuity equations as follows:

eqs_con = [
    IC1 - IL1 + IS + IR2,  # == 0
    IR1 - IL1,  # == 0
    IC2 - IS - IR2,  # == 0
    IR3 - IS - IR2,  # == 0
]

Write the compatibility equations as follows:

eqs_com = [
    VL1 + VR1 - VS + VC1,  # == 0
    VL2 - VC2,  # == 0
    VR2 - VC1 + VC2 - VR3,  # == 0
]

Write the impedances as a substitution dictionary as follows:

impedance_rules = {
	ZR1: R1,
	ZR2: R2,
	ZR3: R3,
	ZL1: L1 * s,
	ZL2: L2 * s,
	ZC1: 1 / (C1 * s),
	ZC2: 1 / (C2 * s),
}

Define the output and input variables as follows:

Solve the linear algebraic system of equations as follows:

eqs = eqs_ele + eqs_con + eqs_com
sol = sp.solve(eqs, vars_unknown, dict=True)

We now have the solution for all unknown variables in terms of the input variables and the impedances. Let’s print out the solution for the two output variables as follows:

for outie in output_indices:
	print(sol[0][vars_unknown[outie]])
(IS*ZC1*ZC2*ZR2 + IS*ZC2*ZL1*ZR2 + IS*ZC2*ZR1*ZR2 + VS*ZC1*ZC2)/(ZC1*ZC2 + ZC1*ZL1 + ZC1*ZR1 + ZC1*ZR2 - ZC1*ZR3 + ZC2*ZL1 + ZC2*ZR1 + ZL1*ZR2 - ZL1*ZR3 + ZR1*ZR2 - ZR1*ZR3)
(IS*ZC1*ZR2 + VS*ZC1 + VS*ZC2 + VS*ZR2 - VS*ZR3)/(ZC1*ZC2 + ZC1*ZL1 + ZC1*ZR1 + ZC1*ZR2 - ZC1*ZR3 + ZC2*ZL1 + ZC2*ZR1 + ZL1*ZR2 - ZL1*ZR3 + ZR1*ZR2 - ZR1*ZR3)

Note that the output variables are solved in terms of both of the input variables and the impedances. Let’s generate a matrix of transfer functions in the form $$H_{ij}(s) = \frac{Y_i(s)}{U_j(s)},$$ where Yi(s) is the ith output variable and Uj(s) is the jth input variable.

H = sp.zeros(len(vars_unknown), len(vars_input))
for i, un_ in enumerate(vars_unknown):
	for j, in_ in enumerate(vars_input):
		# Create dict to sub all other input variables to zero
		sub_dict = {an_in_: 0 for an_in_ in vars_input if an_in_ != in_}
		H[i, j] = sol[0][un_].subs(sub_dict).subs(impedance_rules)
		H[i, j] = H[i, j].simplify()
		H[i, j] = H[i, j].expand(numer=True).expand(denom=True)
		H[i, j] = sp.rcollect(H[i, j], s)  # Collect terms in s

Now print the transfer functions of interest as follows:

for i, j in zip(output_indices, input_indices):
	print(f"H_{i}{j}(s) = Y_{i}(s) / U_{j}(s)")
	print(H[i, j])
H_60(s) = Y_6(s) / U_0(s) = 
VS/(s**3*(C1*C2*L1*R2 - C1*C2*L1*R3) + s**2*(C1*C2*R1*R2 - C1*C2*R1*R3 + C1*L1 + C2*L1) + s*(C1*R1 + C2*R1 + C2*R2 - C2*R3) + 1)
H_71(s) = Y_7(s) / U_1(s) = 
C2*IS*R2*s/(s**3*(C1*C2*L1*R2 - C1*C2*L1*R3) + s**2*(C1*C2*R1*R2 - C1*C2*R1*R3 + C1*L1 + C2*L1) + s*(C1*R1 + C2*R1 + C2*R2 - C2*R3) + 1)

Online Resources for Section 12.4

No online resources.