3  SMNA Example

Figure 3.1: Example circuit

3.1 Introduction

This chapter walks through the Python code used to generate and solve the circuit network equations. Figure 3.1 is the schematic for the circuit used in this example. The analysis procedure first requires a circuit net list, which can be generated by hand with a text editor. For small circuits, such as those assigned as homework problems, this is not difficult since the number of nodes and components is small. College textbook problems are usually meant by the authors to be solved by hand with pencil and paper. In this example LTSpice was used to draw the schematic and label the components and nodes. Most schematic capture programs have the ability to export a SPICE net list, which then can be pasted into the procedure described here.

3.2 Circuit description

The circuit in Figure 3.1 is a 2nd order band pass filter with magnetic coupling. The netlist generated by LTSpice is shown below.

V1 1 0 AC 1
R1 3 1 1k
R4 2 0 10k
C1 3 0 0.01µ
C2 2 0 0.01µ
L1 4 0 1µ
L2 5 0 1µ
R2 4 3 0.5
R3 2 5 0.5
K1 L1 L2 0.15

The LTSpice netlist requires some editing. The line in the netlist that defines the independent AC source, V1, needs to be formatted as if it were a DC source. For symbolic analysis we are only concerned with the label for the source at this time. Later, in this example, an AC analysis will be performed and \(j \omega\) will be substituted for the Laplace variable, \(s\). Also, the suffixes that SPICE allows in the component values e.g., k and \(\mu\), need to be replaced by their corresponding multiplication factors as shown below.

V1 1 0 1
R1 3 1 1e3
R4 2 0 10e3
C1 3 0 0.01e-6
C2 2 0 0.01e-6
L1 4 0 1e-6
L2 5 0 1e-6
R2 4 3 0.5
R3 2 5 0.5
K1 L1 L2 0.15

The following Python modules are used:

import os
from sympy import *
import numpy as np
from scipy import signal
import matplotlib.pyplot as plt
import pandas as pd
from tabulate import tabulate
from IPython.display import display, Markdown, Math, Latex
init_printing()

3.3 Symbolic MNA code

A count of the component types are initialized to zero.

# initialize variables
num_rlc = 0 # number of passive elements
num_ind = 0 # number of inductors
num_v = 0    # number of independent voltage sources
num_i = 0    # number of independent current sources
i_unk = 0  # number of current unknowns
num_opamps = 0   # number of op amps
num_vcvs = 0     # number of controlled sources of various types
num_vccs = 0
num_cccs = 0
num_ccvs = 0
num_cpld_ind = 0 # number of coupled inductors

3.4 Read the net list and preprocess it

The circuit netlist is pasted into the code cell below. A new line character is required at the end of each line and the triple quotes in the code cell below preserve the line breaks.

net_list = '''
V1 1 0 1
R1 3 1 1e3
R4 2 0 10e3
C1 3 0 0.01e-6
C2 2 0 0.01e-6
L1 4 0 1e-6
L2 5 0 1e-6
R2 4 3 0.5
R3 2 5 0.5
K1 L1 L2 0.15
'''

The code cell below performs the following operations:

  • split the net list into a list of lines at the line breaks
  • remove blank lines and comments
  • convert first letter of element name to uppercase
  • removes extra spaces between entries
content = net_list.splitlines()

content = [x.strip() for x in content]  #remove leading and trailing white space
# remove empty lines
while '' in content:
    content.pop(content.index(''))

# remove comment lines, these start with a asterisk *
content = [n for n in content if not n.startswith('*')]
# remove other comment lines, these start with a semicolon ;
content = [n for n in content if not n.startswith(';')]
# remove SPICE directives, these start with a period, .
content = [n for n in content if not n.startswith('.')]
# converts 1st letter to upper case
#content = [x.upper() for x in content] <- this converts all to upper case
content = [x.capitalize() for x in content]
# removes extra spaces between entries
content = [' '.join(x.split()) for x in content]
# display the cleaned up netlist
for i in content:
    print(i)
V1 1 0 1
R1 3 1 1e3
R4 2 0 10e3
C1 3 0 0.01e-6
C2 2 0 0.01e-6
L1 4 0 1e-6
L2 5 0 1e-6
R2 4 3 0.5
R3 2 5 0.5
K1 l1 l2 0.15

3.4.1 Process each line in the netlist

line_cnt = len(content) # number of lines in the netlist
branch_cnt = 0  # number of branches in the netlist
# check number of entries on each line, count each element type
for i in range(line_cnt):
    x = content[i][0]
    tk_cnt = len(content[i].split()) # split the line into a list of words

    if (x == 'R') or (x == 'L') or (x == 'C'):
        if tk_cnt != 4:
            print("branch {:d} not formatted correctly, {:s}".format(i,content[i]))
            print("had {:d} items and should only be 4".format(tk_cnt))
        num_rlc += 1
        branch_cnt += 1
        if x == 'L':
            num_ind += 1
    elif x == 'V':
        if tk_cnt != 4:
            print("branch {:d} not formatted correctly, {:s}".format(i,content[i]))
            print("had {:d} items and should only be 4".format(tk_cnt))
        num_v += 1
        branch_cnt += 1
    elif x == 'I':
        if tk_cnt != 4:
            print("branch {:d} not formatted correctly, {:s}".format(i,content[i]))
            print("had {:d} items and should only be 4".format(tk_cnt))
        num_i += 1
        branch_cnt += 1
    elif x == 'O':
        if tk_cnt != 4:
            print("branch {:d} not formatted correctly, {:s}".format(i,content[i]))
            print("had {:d} items and should only be 4".format(tk_cnt))
        num_opamps += 1
    elif x == 'E':
        if (tk_cnt != 6):
            print("branch {:d} not formatted correctly, {:s}".format(i,content[i]))
            print("had {:d} items and should only be 6".format(tk_cnt))
        num_vcvs += 1
        branch_cnt += 1
    elif x == 'G':
        if (tk_cnt != 6):
            print("branch {:d} not formatted correctly, {:s}".format(i,content[i]))
            print("had {:d} items and should only be 6".format(tk_cnt))
        num_vccs += 1
        branch_cnt += 1
    elif x == 'F':
        if (tk_cnt != 5):
            print("branch {:d} not formatted correctly, {:s}".format(i,content[i]))
            print("had {:d} items and should only be 5".format(tk_cnt))
        num_cccs += 1
        branch_cnt += 1
    elif x == 'H':
        if (tk_cnt != 5):
            print("branch {:d} not formatted correctly, {:s}".format(i,content[i]))
            print("had {:d} items and should only be 5".format(tk_cnt))
        num_ccvs += 1
        branch_cnt += 1
    elif x == 'K':
        if (tk_cnt != 4):
            print("branch {:d} not formatted correctly, {:s}".format(i,content[i]))
            print("had {:d} items and should only be 4".format(tk_cnt))
        num_cpld_ind += 1
    else:
        print("unknown element type in branch {:d}, {:s}".format(i,content[i]))

3.5 Parser

The parser performs the following operations.

  1. puts branch elements into data frame
  2. counts number of nodes

data frame labels:

  • element: type of element
  • p node: positive node
  • n node: negative node, for a current source, the arrow point terminal, LTSpice puts the inductor phasing dot on this terminal
  • cp node: controlling positive node of branch
  • cn node: controlling negative node of branch
  • Vout: Op Amp output node
  • value: value of element or voltage
  • Vname: voltage source through which the controlling current flows. Need to add a zero volt voltage source to the controlling branch.
  • Lname1: name of coupled inductor 1
  • Lname2: name of coupled inductor 2
# build the pandas data frame
df = pd.DataFrame(columns=['element','p node','n node','cp node','cn node',
    'Vout','value','Vname','Lname1','Lname2'])

# this data frame is for branches with unknown currents
df2 = pd.DataFrame(columns=['element','p node','n node'])

3.5.1 Functions to load branch elements into data frame and check for gaps in node numbering

# loads voltage or current sources into branch structure
def indep_source(line_nu):
    tk = content[line_nu].split()
    df.loc[line_nu,'element'] = tk[0]
    df.loc[line_nu,'p node'] = int(tk[1])
    df.loc[line_nu,'n node'] = int(tk[2])
    df.loc[line_nu,'value'] = float(tk[3])

# loads passive elements into branch structure
def rlc_element(line_nu):
    tk = content[line_nu].split()
    df.loc[line_nu,'element'] = tk[0]
    df.loc[line_nu,'p node'] = int(tk[1])
    df.loc[line_nu,'n node'] = int(tk[2])
    df.loc[line_nu,'value'] = float(tk[3])

# loads multi-terminal elements into branch structure
# O - Op Amps
def opamp_sub_network(line_nu):
    tk = content[line_nu].split()
    df.loc[line_nu,'element'] = tk[0]
    df.loc[line_nu,'p node'] = int(tk[1])
    df.loc[line_nu,'n node'] = int(tk[2])
    df.loc[line_nu,'Vout'] = int(tk[3])

# G - VCCS
def vccs_sub_network(line_nu):
    tk = content[line_nu].split()
    df.loc[line_nu,'element'] = tk[0]
    df.loc[line_nu,'p node'] = int(tk[1])
    df.loc[line_nu,'n node'] = int(tk[2])
    df.loc[line_nu,'cp node'] = int(tk[3])
    df.loc[line_nu,'cn node'] = int(tk[4])
    df.loc[line_nu,'value'] = float(tk[5])

# E - VCVS
# in SymPy E is the number 2.718, replacing E with Ea otherwise, sympify() errors out
def vcvs_sub_network(line_nu):
    tk = content[line_nu].split()
    df.loc[line_nu,'element'] = tk[0].replace('E', 'Ea')
    df.loc[line_nu,'p node'] = int(tk[1])
    df.loc[line_nu,'n node'] = int(tk[2])
    df.loc[line_nu,'cp node'] = int(tk[3])
    df.loc[line_nu,'cn node'] = int(tk[4])
    df.loc[line_nu,'value'] = float(tk[5])

# F - CCCS
def cccs_sub_network(line_nu):
    tk = content[line_nu].split()
    df.loc[line_nu,'element'] = tk[0]
    df.loc[line_nu,'p node'] = int(tk[1])
    df.loc[line_nu,'n node'] = int(tk[2])
    df.loc[line_nu,'Vname'] = tk[3].capitalize()
    df.loc[line_nu,'value'] = float(tk[4])

# H - CCVS
def ccvs_sub_network(line_nu):
    tk = content[line_nu].split()
    df.loc[line_nu,'element'] = tk[0]
    df.loc[line_nu,'p node'] = int(tk[1])
    df.loc[line_nu,'n node'] = int(tk[2])
    df.loc[line_nu,'Vname'] = tk[3].capitalize()
    df.loc[line_nu,'value'] = float(tk[4])

# K - Coupled inductors
def cpld_ind_sub_network(line_nu):
    tk = content[line_nu].split()
    df.loc[line_nu,'element'] = tk[0]
    df.loc[line_nu,'Lname1'] = tk[1].capitalize()
    df.loc[line_nu,'Lname2'] = tk[2].capitalize()
    df.loc[line_nu,'value'] = float(tk[3])

# function to scan df and get largest node number
def count_nodes():
    # need to check that nodes are consecutive
    # fill array with node numbers
    p = np.zeros(line_cnt+1)
    for i in range(line_cnt):
        # need to skip coupled inductor 'K' statements
        if df.loc[i,'element'][0] != 'K': #get 1st letter of element name
            p[df['p node'][i]] = df['p node'][i]
            p[df['n node'][i]] = df['n node'][i]

    # find the largest node number
    if df['n node'].max() > df['p node'].max():
        largest = df['n node'].max()
    else:
        largest =  df['p node'].max()

    largest = int(largest)
    # check for unfilled elements, skip node 0
    for i in range(1,largest):
        if p[i] == 0:
            print('nodes not in continuous order, node {:.0f} is missing'.format(p[i-1]+1))

    return largest

3.5.2 Load circuit netlist into the data frames

# load branch info into data frame
for i in range(line_cnt):
    x = content[i][0]

    if (x == 'R') or (x == 'L') or (x == 'C'):
        rlc_element(i)
    elif (x == 'V') or (x == 'I'):
        indep_source(i)
    elif x == 'O':
        opamp_sub_network(i)
    elif x == 'E':
        vcvs_sub_network(i)
    elif x == 'G':
        vccs_sub_network(i)
    elif x == 'F':
        cccs_sub_network(i)
    elif x == 'H':
        ccvs_sub_network(i)
    elif x == 'K':
        cpld_ind_sub_network(i)
    else:
        print("unknown element type in branch {:d}, {:s}".format(i,content[i]))

29 Nov 2023: When the D matrix is built, independent voltage sources are processed in the data frame order when building the D matrix. If the voltage source followed element L, H, F, K types in the netlist, a row was inserted that put the voltage source in a different row in relation to its position in the Ev matrix. This would cause the node attached to the terminal of the voltage source to be zero volts.
Solution - The following block of code was added to move voltage source types to the beginning of the net list data frame before any calculations are performed.

# Check for the position of voltage sources in the data frame.
source_index = [] # keep track of voltage source row number
other_index = [] # make a list of all other types
for i in range(len(df)):
    # process all the elements creating unknown currents
    x = df.loc[i,'element'][0]   #get 1st letter of element name
    if (x == 'V'):
        source_index.append(i)
    else:
        other_index.append(i)

df = df.reindex(source_index+other_index,copy=True) # reorder the data frame
df.reset_index(drop=True, inplace=True) # renumber the index
# count number of nodes
num_nodes = count_nodes()

# Build df2: consists of branches with current unknowns, used for C & D matrices
# walk through data frame and find these parameters
count = 0
for i in range(len(df)):
    # process all the elements creating unknown currents
    x = df.loc[i,'element'][0]   #get 1st letter of element name
    if (x == 'L') or (x == 'V') or (x == 'O') or (x == 'E') or (x == 'H') or (x == 'F'):
        df2.loc[count,'element'] = df.loc[i,'element']
        df2.loc[count,'p node'] = df.loc[i,'p node']
        df2.loc[count,'n node'] = df.loc[i,'n node']
        count += 1

3.7 Initialize matrices

V = zeros(num_nodes,1)
I = zeros(num_nodes,1)
G = zeros(num_nodes,num_nodes)  # also called Yr, the reduced nodal matrix
s = Symbol('s')  # the Laplace variable

# count the number of element types that affect the size of the B, C, D, E and J arrays
# these are element types that have unknown currents
i_unk = num_v+num_opamps+num_vcvs+num_ccvs+num_ind+num_cccs
# if i_unk == 0, just generate empty arrays
B = zeros(num_nodes,i_unk)
C = zeros(i_unk,num_nodes)
D = zeros(i_unk,i_unk)
Ev = zeros(i_unk,1)
J = zeros(i_unk,1)

Debugging notes: Is it possible to have i_unk == 0, what about a network with only current sources? This would make B = 0. See test_14 and test_15.

3.8 G matrix

The G matrix is n by n, where n is the number of nodes. The matrix is formed by the interconnections between the resistors, capacitors and VCCS type elements. In the original paper G is called Yr, where Yr is a reduced form of the nodal matrix excluding the contributions due to voltage sources, current controlling elements, etc. In python row and columns are: G[row, column]

# G matrix
for i in range(len(df)):  # process each row in the data frame
    n1 = df.loc[i,'p node']
    n2 = df.loc[i,'n node']
    cn1 = df.loc[i,'cp node']
    cn2 = df.loc[i,'cn node']
    # process all the passive elements, save conductance to temp value
    x = df.loc[i,'element'][0]   #get 1st letter of element name
    if x == 'R':
        g = 1/sympify(df.loc[i,'element'])
    if x == 'C':
        g = s*sympify(df.loc[i,'element'])
    if x == 'G':   #vccs type element
        g = sympify(df.loc[i,'element'].lower())  # use a symbol for gain value

    if (x == 'R') or (x == 'C'):
        # If neither side of the element is connected to ground
        # then subtract it from the appropriate location in the matrix.
        if (n1 != 0) and (n2 != 0):
            G[n1-1,n2-1] += -g
            G[n2-1,n1-1] += -g

        # If node 1 is connected to ground, add element to diagonal of matrix
        if n1 != 0:
            G[n1-1,n1-1] += g

        # same for for node 2
        if n2 != 0:
            G[n2-1,n2-1] += g

    if x == 'G':    #vccs type element
        # check to see if any terminal is grounded
        # then stamp the matrix
        if n1 != 0 and cn1 != 0:
            G[n1-1,cn1-1] += g

        if n2 != 0 and cn2 != 0:
            G[n2-1,cn2-1] += g

        if n1 != 0 and cn2 != 0:
            G[n1-1,cn2-1] -= g

        if n2 != 0 and cn1 != 0:
            G[n2-1,cn1-1] -= g

G  # display the G matrix

\(\displaystyle \left[\begin{matrix}\frac{1}{R_{1}} & 0 & - \frac{1}{R_{1}} & 0 & 0\\0 & C_{2} s + \frac{1}{R_{4}} + \frac{1}{R_{3}} & 0 & 0 & - \frac{1}{R_{3}}\\- \frac{1}{R_{1}} & 0 & C_{1} s + \frac{1}{R_{2}} + \frac{1}{R_{1}} & - \frac{1}{R_{2}} & 0\\0 & 0 & - \frac{1}{R_{2}} & \frac{1}{R_{2}} & 0\\0 & - \frac{1}{R_{3}} & 0 & 0 & \frac{1}{R_{3}}\end{matrix}\right]\)

3.9 B Matrix

The B matrix is an n by m matrix with only 0, 1 and -1 elements, where n = number of nodes and m is the number of current unknowns, i_unk. There is one column for each unknown current. The code loop through all the branches and process elements that have stamps for the B matrix:

  • Voltage sources (V)
  • Op Amps (O)
  • CCVS (H)
  • CCCS (F)
  • VCVS (E)
  • Inductors (L)

The order of the columns is as they appear in the netlist. CCCS (F) does not get its own column because the controlling current is through a zero volt voltage source, called Vname and is already in the net list.

# generate the B Matrix
sn = 0   # count source number as code walks through the data frame
for i in range(len(df)):
    n1 = df.loc[i,'p node']
    n2 = df.loc[i,'n node']
    n_vout = df.loc[i,'Vout'] # node connected to Op Amp output

    # process elements with input to B matrix
    x = df.loc[i,'element'][0]   #get 1st letter of element name
    if x == 'V':
        if i_unk > 1:  #is B greater than 1 by n?, V
            if n1 != 0:
                B[n1-1,sn] = 1
            if n2 != 0:
                B[n2-1,sn] = -1
        else:
            if n1 != 0:
                B[n1-1] = 1
            if n2 != 0:
                B[n2-1] = -1
        sn += 1   #increment source count
    if x == 'O':  # op amp type, output connection of the Op Amp goes in the B matrix
        B[n_vout-1,sn] = 1
        sn += 1   # increment source count
    if (x == 'H') or (x == 'F'):  # H: ccvs, F: cccs,
        if i_unk > 1:  #is B greater than 1 by n?, H, F
            # check to see if any terminal is grounded
            # then stamp the matrix
            if n1 != 0:
                B[n1-1,sn] = 1
            if n2 != 0:
                B[n2-1,sn] = -1
        else:
            if n1 != 0:
                B[n1-1] = 1
            if n2 != 0:
                B[n2-1] = -1
        sn += 1   #increment source count
    if x == 'E':   # vcvs type, only ik column is altered at n1 and n2
        if i_unk > 1:  #is B greater than 1 by n?, E
            if n1 != 0:
                B[n1-1,sn] = 1
            if n2 != 0:
                B[n2-1,sn] = -1
        else:
            if n1 != 0:
                B[n1-1] = 1
            if n2 != 0:
                B[n2-1] = -1
        sn += 1   #increment source count
    if x == 'L':
        if i_unk > 1:  #is B greater than 1 by n?, L
            if n1 != 0:
                B[n1-1,sn] = 1
            if n2 != 0:
                B[n2-1,sn] = -1
        else:
            if n1 != 0:
                B[n1-1] = 1
            if n2 != 0:
                B[n2-1] = -1
        sn += 1   #increment source count

# check source count
if sn != i_unk:
    print('source number, sn={:d} not equal to i_unk={:d} in matrix B'.format(sn,i_unk))

B   # display the B matrix

\(\displaystyle \left[\begin{matrix}1 & 0 & 0\\0 & 0 & 0\\0 & 0 & 0\\0 & 1 & 0\\0 & 0 & 1\end{matrix}\right]\)

3.10 C matrix

The C matrix is an m by n matrix with only 0, 1 and -1 elements (except for controlled sources). The code is similar to the B matrix code, except the indices are swapped. The code loops through all the branches and process elements that have stamps for the C matrix:

  • Voltage sources (V)
  • Op Amps (O)
  • CCVS (H)
  • CCCS (F)
  • VCVS (E)
  • Inductors (L)

3.10.1 Op Amp elements

The Op Amp element is assumed to be an ideal Op Amp and use of this component is valid only when used in circuits with a DC path (a short or a resistor) from the output terminal to the negative input terminal of the Op Amp. No error checking is provided and if the condition is violated, the results likely will be erroneous. Chen (2018) and Fakhfakh, Tlelo-Cuautle, and Fernandez (2012) were consulted during the debugging of the Op Amp stamp.

# find the the column position in the C and D matrix for controlled sources
# needs to return the node numbers and branch number of controlling branch
def find_vname(name):
    # need to walk through data frame and find these parameters
    for i in range(len(df2)):
        # process all the elements creating unknown currents
        if name == df2.loc[i,'element']:
            n1 = df2.loc[i,'p node']
            n2 = df2.loc[i,'n node']
            return n1, n2, i  # n1, n2 & col_num are from the branch of the controlling element

    print('failed to find matching branch element in find_vname')
# generate the C Matrix
sn = 0   # count source number as code walks through the data frame
for i in range(len(df)):
    n1 = df.loc[i,'p node']
    n2 = df.loc[i,'n node']
    cn1 = df.loc[i,'cp node'] # nodes for controlled sources
    cn2 = df.loc[i,'cn node']
    n_vout = df.loc[i,'Vout'] # node connected to op amp output

    # process elements with input to B matrix
    x = df.loc[i,'element'][0]   #get 1st letter of element name
    if x == 'V':
        if i_unk > 1:  #is B greater than 1 by n?, V
            if n1 != 0:
                C[sn,n1-1] = 1
            if n2 != 0:
                C[sn,n2-1] = -1
        else:
            if n1 != 0:
                C[n1-1] = 1
            if n2 != 0:
                C[n2-1] = -1
        sn += 1   #increment source count

    if x == 'O':  # Op Amp type, input connections of the Op Amp go into the C matrix
        # C[sn,n_vout-1] = 1
        if i_unk > 1:  #is B greater than 1 by n?, O
            # check to see if any terminal is grounded
            # then stamp the matrix
            if n1 != 0:
                C[sn,n1-1] = 1
            if n2 != 0:
                C[sn,n2-1] = -1
        else:
            if n1 != 0:
                C[n1-1] = 1
            if n2 != 0:
                C[n2-1] = -1
        sn += 1   # increment source count

    if x == 'F':  # need to count F (cccs) types
        sn += 1   #increment source count
    if x == 'H':  # H: ccvs
        if i_unk > 1:  #is B greater than 1 by n?, H
            # check to see if any terminal is grounded
            # then stamp the matrix
            if n1 != 0:
                C[sn,n1-1] = 1
            if n2 != 0:
                C[sn,n2-1] = -1
        else:
            if n1 != 0:
                C[n1-1] = 1
            if n2 != 0:
                C[n2-1] = -1
        sn += 1   #increment source count
    if x == 'E':   # vcvs type, ik column is altered at n1 and n2, cn1 & cn2 get value
        if i_unk > 1:  #is B greater than 1 by n?, E
            if n1 != 0:
                C[sn,n1-1] = 1
            if n2 != 0:
                C[sn,n2-1] = -1
            # add entry for cp and cn of the controlling voltage
            if cn1 != 0:
                C[sn,cn1-1] = -sympify(df.loc[i,'element'].lower())
            if cn2 != 0:
                C[sn,cn2-1] = sympify(df.loc[i,'element'].lower())
        else:
            if n1 != 0:
                C[n1-1] = 1
            if n2 != 0:
                C[n2-1] = -1
            vn1, vn2, df2_index = find_vname(df.loc[i,'Vname'])
            if vn1 != 0:
                C[vn1-1] = -sympify(df.loc[i,'element'].lower())
            if vn2 != 0:
                C[vn2-1] = sympify(df.loc[i,'element'].lower())
        sn += 1   #increment source count

    if x == 'L':
        if i_unk > 1:  #is B greater than 1 by n?, L
            if n1 != 0:
                C[sn,n1-1] = 1
            if n2 != 0:
                C[sn,n2-1] = -1
        else:
            if n1 != 0:
                C[n1-1] = 1
            if n2 != 0:
                C[n2-1] = -1
        sn += 1   #increment source count

# check source count
if sn != i_unk:
    print('source number, sn={:d} not equal to i_unk={:d} in matrix C'.format(sn,i_unk))

C   # display the C matrix

\(\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 1 & 0\\0 & 0 & 0 & 0 & 1\end{matrix}\right]\)

3.11 D matrix

The D matrix is an m by m matrix, where m is the number of unknown currents.

m = i_unk = num_v+num_opamps+num_vcvs+num_ccvs+num_ind+num_cccs

Stamps that affect the D matrix are: inductor, ccvs and cccs
inductors: minus sign added to keep current flow convention consistent

Coupled inductors notes: 12/6/2017 doing some debugging on with coupled inductors; LTSpice seems to put the phasing dot on the neg node when it generates the netlist. This Python code uses M for mutual inductance, LTSpice uses k for the coupling coefficient. Inductors in LTSpice have a \(20m\Omega\) series resistance that needs to be set to zero.

# generate the D Matrix
sn = 0   # count source number as code walks through the data frame
for i in range(len(df)):
    n1 = df.loc[i,'p node']
    n2 = df.loc[i,'n node']
    #cn1 = df.loc[i,'cp node'] # nodes for controlled sources
    #cn2 = df.loc[i,'cn node']
    #n_vout = df.loc[i,'Vout'] # node connected to op amp output

    # process elements with input to D matrix
    x = df.loc[i,'element'][0]   #get 1st letter of element name
    if (x == 'V') or (x == 'O') or (x == 'E'):  # need to count V, E & O types
        sn += 1   #increment source count

    if x == 'L':
        if i_unk > 1:  #is D greater than 1 by 1?
            D[sn,sn] += -s*sympify(df.loc[i,'element'])
        else:
            D[sn] += -s*sympify(df.loc[i,'element'])
        sn += 1   #increment source count

    if x == 'H':  # H: ccvs
        # if there is a H type, D is m by m
        # need to find the vn for Vname
        # then stamp the matrix
        vn1, vn2, df2_index = find_vname(df.loc[i,'Vname'])
        D[sn,df2_index] += -sympify(df.loc[i,'element'].lower())
        sn += 1   #increment source count

    if x == 'F':  # F: cccs
        # if there is a F type, D is m by m
        # need to find the vn for Vname
        # then stamp the matrix
        vn1, vn2, df2_index = find_vname(df.loc[i,'Vname'])
        D[sn,df2_index] += -sympify(df.loc[i,'element'].lower())
        D[sn,sn] = 1
        sn += 1   #increment source count

    if x == 'K':  # K: coupled inductors, KXX LYY LZZ value
        # if there is a K type, D is m by m
        vn1, vn2, ind1_index = find_vname(df.loc[i,'Lname1'])  # get i_unk position for Lx
        vn1, vn2, ind2_index = find_vname(df.loc[i,'Lname2'])  # get i_unk position for Ly
        # enter sM on diagonals = value*sqrt(LXX*LZZ)

        D[ind1_index,ind2_index] += -s*sympify('M{:s}'.format(df.loc[i,'element'].lower()[1:]))  # s*Mxx
        D[ind2_index,ind1_index] += -s*sympify('M{:s}'.format(df.loc[i,'element'].lower()[1:]))  # -s*Mxx

# display the The D matrix
D

\(\displaystyle \left[\begin{matrix}0 & 0 & 0\\0 & - L_{1} s & - M_{1} s\\0 & - M_{1} s & - L_{2} s\end{matrix}\right]\)

3.12 V matrix

The V matrix is an n by 1 matrix formed of the node voltages, where n is the number of nodes. Each element in V corresponds to the voltage at the node.

Maybe make small v’s v_1 so as not to confuse v1 with V1.

# generate the V matrix
for i in range(num_nodes):
    V[i] = sympify('v{:d}'.format(i+1))

V  # display the V matrix

\(\displaystyle \left[\begin{matrix}v_{1}\\v_{2}\\v_{3}\\v_{4}\\v_{5}\end{matrix}\right]\)

3.13 J matrix

The J matrix is an m by 1 matrix, where m is the number of unknown currents. >i_unk = num_v+num_opamps+num_vcvs+num_ccvs+num_ind+num_cccs

# The J matrix is an mx1 matrix, with one entry for each i_unk from a source
#sn = 0   # count i_unk source number
#oan = 0   #count op amp number
for i in range(len(df2)):
    # process all the unknown currents
    J[i] = sympify('I_{:s}'.format(df2.loc[i,'element']))

J  # diplay the J matrix

\(\displaystyle \left[\begin{matrix}I_{V1}\\I_{L1}\\I_{L2}\end{matrix}\right]\)

3.14 I matrix

The I matrix is an n by 1 matrix, where n is the number of nodes. The value of each element of I is determined by the sum of current sources into the corresponding node. If there are no current sources connected to the node, the value is zero.

# generate the I matrix, current sources have n2 = arrow end of the element
for i in range(len(df)):
    n1 = df.loc[i,'p node']
    n2 = df.loc[i,'n node']
    # process all the passive elements, save conductance to temp value
    x = df.loc[i,'element'][0]   #get 1st letter of element name
    if x == 'I':
        g = sympify(df.loc[i,'element'])
        # sum the current into each node
        if n1 != 0:
            I[n1-1] -= g
        if n2 != 0:
            I[n2-1] += g

I  # display the I matrix

\(\displaystyle \left[\begin{matrix}0\\0\\0\\0\\0\end{matrix}\right]\)

3.15 Ev matrix

The Ev matrix is mx1 and holds the values of the independent voltage sources.

# generate the E matrix
sn = 0   # count source number
for i in range(len(df)):
    # process all the passive elements
    x = df.loc[i,'element'][0]   #get 1st letter of element name
    if x == 'V':
        Ev[sn] = sympify(df.loc[i,'element'])
        sn += 1

Ev   # display the E matrix

\(\displaystyle \left[\begin{matrix}V_{1}\\0\\0\end{matrix}\right]\)

3.16 Z matrix

The Z matrix holds the independent voltage and current sources and is the combination of 2 smaller matrices I and Ev. The Z matrix is (m+n) by 1, n is the number of nodes, and m is the number of independent voltage sources. The I matrix is n by 1 and contains the sum of the currents through the passive elements into the corresponding node (either zero, or the sum of independent current sources). The Ev matrix is m by 1 and holds the values of the independent voltage sources.

Z = I[:] + Ev[:]  # the + operator in python concatenates the lists
Z  # display the Z matrix

\(\displaystyle \left[ 0, \ 0, \ 0, \ 0, \ 0, \ V_{1}, \ 0, \ 0\right]\)

3.17 X matrix

The X matrix is an (n+m) by 1 vector that holds the unknown quantities (node voltages and the currents through the independent voltage sources). The top n elements are the n node voltages. The bottom m elements represent the currents through the m independent voltage sources in the circuit. The V matrix is n by 1 and holds the unknown voltages. The J matrix is m by 1 and holds the unknown currents through the voltage sources

X = V[:] + J[:]  # the + operator in python concatenates the lists
X  # display the X matrix

\(\displaystyle \left[ v_{1}, \ v_{2}, \ v_{3}, \ v_{4}, \ v_{5}, \ I_{V1}, \ I_{L1}, \ I_{L2}\right]\)

3.18 A matrix

The A matrix is (m+n) by (m+n) and will be developed as the combination of 4 smaller matrices, G, B, C, and D.

n = num_nodes
m = i_unk
A = zeros(m+n,m+n)
for i in range(n):
    for j in range(n):
        A[i,j] = G[i,j]

if i_unk > 1:
    for i in range(n):
        for j in range(m):
            A[i,n+j] = B[i,j]
            A[n+j,i] = C[j,i]

    for i in range(m):
        for j in range(m):
            A[n+i,n+j] = D[i,j]

if i_unk == 1:
    for i in range(n):
        A[i,n] = B[i]
        A[n,i] = C[i]
    A[n,n] = D[0] # added 1/7/2024 while debugging source free circuit with one inductor

A  # display the A matrix

\(\displaystyle \left[\begin{matrix}\frac{1}{R_{1}} & 0 & - \frac{1}{R_{1}} & 0 & 0 & 1 & 0 & 0\\0 & C_{2} s + \frac{1}{R_{4}} + \frac{1}{R_{3}} & 0 & 0 & - \frac{1}{R_{3}} & 0 & 0 & 0\\- \frac{1}{R_{1}} & 0 & C_{1} s + \frac{1}{R_{2}} + \frac{1}{R_{1}} & - \frac{1}{R_{2}} & 0 & 0 & 0 & 0\\0 & 0 & - \frac{1}{R_{2}} & \frac{1}{R_{2}} & 0 & 0 & 1 & 0\\0 & - \frac{1}{R_{3}} & 0 & 0 & \frac{1}{R_{3}} & 0 & 0 & 1\\1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\0 & 0 & 0 & 1 & 0 & 0 & - L_{1} s & - M_{1} s\\0 & 0 & 0 & 0 & 1 & 0 & - M_{1} s & - L_{2} s\end{matrix}\right]\)

3.19 generate the circuit equations

NE_sym = Eq(A*Matrix(X),Matrix(Z))

Display the network 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 = I_{V1} + \frac{v_{1}}{R_{1}} - \frac{v_{3}}{R_{1}}\)
\(0 = v_{2} \left(C_{2} s + \frac{1}{R_{4}} + \frac{1}{R_{3}}\right) - \frac{v_{5}}{R_{3}}\)
\(0 = v_{3} \left(C_{1} s + \frac{1}{R_{2}} + \frac{1}{R_{1}}\right) - \frac{v_{4}}{R_{2}} - \frac{v_{1}}{R_{1}}\)
\(0 = I_{L1} - \frac{v_{3}}{R_{2}} + \frac{v_{4}}{R_{2}}\)
\(0 = I_{L2} - \frac{v_{2}}{R_{3}} + \frac{v_{5}}{R_{3}}\)
\(V_{1} = v_{1}\)
\(0 = - I_{L1} L_{1} s - I_{L2} M_{1} s + v_{4}\)
\(0 = - I_{L1} M_{1} s - I_{L2} L_{2} s + v_{5}\)

3.20 Symbolic solution

# turn the free symbols into SymPy variables
var(str(NE_sym.free_symbols).replace('{','').replace('}',''))

\(\displaystyle \left( C_{1}, \ v_{2}, \ C_{2}, \ L_{1}, \ I_{L1}, \ I_{V1}, \ R_{3}, \ v_{3}, \ v_{4}, \ v_{5}, \ s, \ R_{4}, \ R_{1}, \ I_{L2}, \ V_{1}, \ R_{2}, \ L_{2}, \ v_{1}, \ M_{1}\right)\)

Solve the equations.

U_sym = solve(NE_sym,X)

Display the solution.

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{M_{1} R_{4} V_{1} s}{C_{1} C_{2} L_{1} L_{2} R_{1} R_{4} s^{4} + C_{1} C_{2} L_{1} R_{1} R_{3} R_{4} s^{3} + C_{1} C_{2} L_{2} R_{1} R_{2} R_{4} s^{3} - C_{1} C_{2} M_{1}^{2} R_{1} R_{4} s^{4} + C_{1} C_{2} R_{1} R_{2} R_{3} R_{4} s^{2} + C_{1} L_{1} L_{2} R_{1} s^{3} + C_{1} L_{1} R_{1} R_{3} s^{2} + C_{1} L_{1} R_{1} R_{4} s^{2} + C_{1} L_{2} R_{1} R_{2} s^{2} - C_{1} M_{1}^{2} R_{1} s^{3} + C_{1} R_{1} R_{2} R_{3} s + C_{1} R_{1} R_{2} R_{4} s + C_{2} L_{1} L_{2} R_{4} s^{3} + C_{2} L_{1} R_{3} R_{4} s^{2} + C_{2} L_{2} R_{1} R_{4} s^{2} + C_{2} L_{2} R_{2} R_{4} s^{2} - C_{2} M_{1}^{2} R_{4} s^{3} + C_{2} R_{1} R_{3} R_{4} s + C_{2} R_{2} R_{3} R_{4} s + L_{1} L_{2} s^{2} + L_{1} R_{3} s + L_{1} R_{4} s + L_{2} R_{1} s + L_{2} R_{2} s - M_{1}^{2} s^{2} + R_{1} R_{3} + R_{1} R_{4} + R_{2} R_{3} + R_{2} R_{4}}\)
\(v_{3} = \frac{C_{2} L_{1} L_{2} R_{4} V_{1} s^{3} + C_{2} L_{1} R_{3} R_{4} V_{1} s^{2} + C_{2} L_{2} R_{2} R_{4} V_{1} s^{2} - C_{2} M_{1}^{2} R_{4} V_{1} s^{3} + C_{2} R_{2} R_{3} R_{4} V_{1} s + L_{1} L_{2} V_{1} s^{2} + L_{1} R_{3} V_{1} s + L_{1} R_{4} V_{1} s + L_{2} R_{2} V_{1} s - M_{1}^{2} V_{1} s^{2} + R_{2} R_{3} V_{1} + R_{2} R_{4} V_{1}}{C_{1} C_{2} L_{1} L_{2} R_{1} R_{4} s^{4} + C_{1} C_{2} L_{1} R_{1} R_{3} R_{4} s^{3} + C_{1} C_{2} L_{2} R_{1} R_{2} R_{4} s^{3} - C_{1} C_{2} M_{1}^{2} R_{1} R_{4} s^{4} + C_{1} C_{2} R_{1} R_{2} R_{3} R_{4} s^{2} + C_{1} L_{1} L_{2} R_{1} s^{3} + C_{1} L_{1} R_{1} R_{3} s^{2} + C_{1} L_{1} R_{1} R_{4} s^{2} + C_{1} L_{2} R_{1} R_{2} s^{2} - C_{1} M_{1}^{2} R_{1} s^{3} + C_{1} R_{1} R_{2} R_{3} s + C_{1} R_{1} R_{2} R_{4} s + C_{2} L_{1} L_{2} R_{4} s^{3} + C_{2} L_{1} R_{3} R_{4} s^{2} + C_{2} L_{2} R_{1} R_{4} s^{2} + C_{2} L_{2} R_{2} R_{4} s^{2} - C_{2} M_{1}^{2} R_{4} s^{3} + C_{2} R_{1} R_{3} R_{4} s + C_{2} R_{2} R_{3} R_{4} s + L_{1} L_{2} s^{2} + L_{1} R_{3} s + L_{1} R_{4} s + L_{2} R_{1} s + L_{2} R_{2} s - M_{1}^{2} s^{2} + R_{1} R_{3} + R_{1} R_{4} + R_{2} R_{3} + R_{2} R_{4}}\)
\(v_{4} = \frac{C_{2} L_{1} L_{2} R_{4} V_{1} s^{3} + C_{2} L_{1} R_{3} R_{4} V_{1} s^{2} - C_{2} M_{1}^{2} R_{4} V_{1} s^{3} + L_{1} L_{2} V_{1} s^{2} + L_{1} R_{3} V_{1} s + L_{1} R_{4} V_{1} s - M_{1}^{2} V_{1} s^{2}}{C_{1} C_{2} L_{1} L_{2} R_{1} R_{4} s^{4} + C_{1} C_{2} L_{1} R_{1} R_{3} R_{4} s^{3} + C_{1} C_{2} L_{2} R_{1} R_{2} R_{4} s^{3} - C_{1} C_{2} M_{1}^{2} R_{1} R_{4} s^{4} + C_{1} C_{2} R_{1} R_{2} R_{3} R_{4} s^{2} + C_{1} L_{1} L_{2} R_{1} s^{3} + C_{1} L_{1} R_{1} R_{3} s^{2} + C_{1} L_{1} R_{1} R_{4} s^{2} + C_{1} L_{2} R_{1} R_{2} s^{2} - C_{1} M_{1}^{2} R_{1} s^{3} + C_{1} R_{1} R_{2} R_{3} s + C_{1} R_{1} R_{2} R_{4} s + C_{2} L_{1} L_{2} R_{4} s^{3} + C_{2} L_{1} R_{3} R_{4} s^{2} + C_{2} L_{2} R_{1} R_{4} s^{2} + C_{2} L_{2} R_{2} R_{4} s^{2} - C_{2} M_{1}^{2} R_{4} s^{3} + C_{2} R_{1} R_{3} R_{4} s + C_{2} R_{2} R_{3} R_{4} s + L_{1} L_{2} s^{2} + L_{1} R_{3} s + L_{1} R_{4} s + L_{2} R_{1} s + L_{2} R_{2} s - M_{1}^{2} s^{2} + R_{1} R_{3} + R_{1} R_{4} + R_{2} R_{3} + R_{2} R_{4}}\)
\(v_{5} = \frac{C_{2} M_{1} R_{3} R_{4} V_{1} s^{2} + M_{1} R_{3} V_{1} s + M_{1} R_{4} V_{1} s}{C_{1} C_{2} L_{1} L_{2} R_{1} R_{4} s^{4} + C_{1} C_{2} L_{1} R_{1} R_{3} R_{4} s^{3} + C_{1} C_{2} L_{2} R_{1} R_{2} R_{4} s^{3} - C_{1} C_{2} M_{1}^{2} R_{1} R_{4} s^{4} + C_{1} C_{2} R_{1} R_{2} R_{3} R_{4} s^{2} + C_{1} L_{1} L_{2} R_{1} s^{3} + C_{1} L_{1} R_{1} R_{3} s^{2} + C_{1} L_{1} R_{1} R_{4} s^{2} + C_{1} L_{2} R_{1} R_{2} s^{2} - C_{1} M_{1}^{2} R_{1} s^{3} + C_{1} R_{1} R_{2} R_{3} s + C_{1} R_{1} R_{2} R_{4} s + C_{2} L_{1} L_{2} R_{4} s^{3} + C_{2} L_{1} R_{3} R_{4} s^{2} + C_{2} L_{2} R_{1} R_{4} s^{2} + C_{2} L_{2} R_{2} R_{4} s^{2} - C_{2} M_{1}^{2} R_{4} s^{3} + C_{2} R_{1} R_{3} R_{4} s + C_{2} R_{2} R_{3} R_{4} s + L_{1} L_{2} s^{2} + L_{1} R_{3} s + L_{1} R_{4} s + L_{2} R_{1} s + L_{2} R_{2} s - M_{1}^{2} s^{2} + R_{1} R_{3} + R_{1} R_{4} + R_{2} R_{3} + R_{2} R_{4}}\)
\(I_{V1} = \frac{- C_{1} C_{2} L_{1} L_{2} R_{4} V_{1} s^{4} - C_{1} C_{2} L_{1} R_{3} R_{4} V_{1} s^{3} - C_{1} C_{2} L_{2} R_{2} R_{4} V_{1} s^{3} + C_{1} C_{2} M_{1}^{2} R_{4} V_{1} s^{4} - C_{1} C_{2} R_{2} R_{3} R_{4} V_{1} s^{2} - C_{1} L_{1} L_{2} V_{1} s^{3} - C_{1} L_{1} R_{3} V_{1} s^{2} - C_{1} L_{1} R_{4} V_{1} s^{2} - C_{1} L_{2} R_{2} V_{1} s^{2} + C_{1} M_{1}^{2} V_{1} s^{3} - C_{1} R_{2} R_{3} V_{1} s - C_{1} R_{2} R_{4} V_{1} s - C_{2} L_{2} R_{4} V_{1} s^{2} - C_{2} R_{3} R_{4} V_{1} s - L_{2} V_{1} s - R_{3} V_{1} - R_{4} V_{1}}{C_{1} C_{2} L_{1} L_{2} R_{1} R_{4} s^{4} + C_{1} C_{2} L_{1} R_{1} R_{3} R_{4} s^{3} + C_{1} C_{2} L_{2} R_{1} R_{2} R_{4} s^{3} - C_{1} C_{2} M_{1}^{2} R_{1} R_{4} s^{4} + C_{1} C_{2} R_{1} R_{2} R_{3} R_{4} s^{2} + C_{1} L_{1} L_{2} R_{1} s^{3} + C_{1} L_{1} R_{1} R_{3} s^{2} + C_{1} L_{1} R_{1} R_{4} s^{2} + C_{1} L_{2} R_{1} R_{2} s^{2} - C_{1} M_{1}^{2} R_{1} s^{3} + C_{1} R_{1} R_{2} R_{3} s + C_{1} R_{1} R_{2} R_{4} s + C_{2} L_{1} L_{2} R_{4} s^{3} + C_{2} L_{1} R_{3} R_{4} s^{2} + C_{2} L_{2} R_{1} R_{4} s^{2} + C_{2} L_{2} R_{2} R_{4} s^{2} - C_{2} M_{1}^{2} R_{4} s^{3} + C_{2} R_{1} R_{3} R_{4} s + C_{2} R_{2} R_{3} R_{4} s + L_{1} L_{2} s^{2} + L_{1} R_{3} s + L_{1} R_{4} s + L_{2} R_{1} s + L_{2} R_{2} s - M_{1}^{2} s^{2} + R_{1} R_{3} + R_{1} R_{4} + R_{2} R_{3} + R_{2} R_{4}}\)
\(I_{L1} = \frac{C_{2} L_{2} R_{4} V_{1} s^{2} + C_{2} R_{3} R_{4} V_{1} s + L_{2} V_{1} s + R_{3} V_{1} + R_{4} V_{1}}{C_{1} C_{2} L_{1} L_{2} R_{1} R_{4} s^{4} + C_{1} C_{2} L_{1} R_{1} R_{3} R_{4} s^{3} + C_{1} C_{2} L_{2} R_{1} R_{2} R_{4} s^{3} - C_{1} C_{2} M_{1}^{2} R_{1} R_{4} s^{4} + C_{1} C_{2} R_{1} R_{2} R_{3} R_{4} s^{2} + C_{1} L_{1} L_{2} R_{1} s^{3} + C_{1} L_{1} R_{1} R_{3} s^{2} + C_{1} L_{1} R_{1} R_{4} s^{2} + C_{1} L_{2} R_{1} R_{2} s^{2} - C_{1} M_{1}^{2} R_{1} s^{3} + C_{1} R_{1} R_{2} R_{3} s + C_{1} R_{1} R_{2} R_{4} s + C_{2} L_{1} L_{2} R_{4} s^{3} + C_{2} L_{1} R_{3} R_{4} s^{2} + C_{2} L_{2} R_{1} R_{4} s^{2} + C_{2} L_{2} R_{2} R_{4} s^{2} - C_{2} M_{1}^{2} R_{4} s^{3} + C_{2} R_{1} R_{3} R_{4} s + C_{2} R_{2} R_{3} R_{4} s + L_{1} L_{2} s^{2} + L_{1} R_{3} s + L_{1} R_{4} s + L_{2} R_{1} s + L_{2} R_{2} s - M_{1}^{2} s^{2} + R_{1} R_{3} + R_{1} R_{4} + R_{2} R_{3} + R_{2} R_{4}}\)
\(I_{L2} = \frac{- C_{2} M_{1} R_{4} V_{1} s^{2} - M_{1} V_{1} s}{C_{1} C_{2} L_{1} L_{2} R_{1} R_{4} s^{4} + C_{1} C_{2} L_{1} R_{1} R_{3} R_{4} s^{3} + C_{1} C_{2} L_{2} R_{1} R_{2} R_{4} s^{3} - C_{1} C_{2} M_{1}^{2} R_{1} R_{4} s^{4} + C_{1} C_{2} R_{1} R_{2} R_{3} R_{4} s^{2} + C_{1} L_{1} L_{2} R_{1} s^{3} + C_{1} L_{1} R_{1} R_{3} s^{2} + C_{1} L_{1} R_{1} R_{4} s^{2} + C_{1} L_{2} R_{1} R_{2} s^{2} - C_{1} M_{1}^{2} R_{1} s^{3} + C_{1} R_{1} R_{2} R_{3} s + C_{1} R_{1} R_{2} R_{4} s + C_{2} L_{1} L_{2} R_{4} s^{3} + C_{2} L_{1} R_{3} R_{4} s^{2} + C_{2} L_{2} R_{1} R_{4} s^{2} + C_{2} L_{2} R_{2} R_{4} s^{2} - C_{2} M_{1}^{2} R_{4} s^{3} + C_{2} R_{1} R_{3} R_{4} s + C_{2} R_{2} R_{3} R_{4} s + L_{1} L_{2} s^{2} + L_{1} R_{3} s + L_{1} R_{4} s + L_{2} R_{1} s + L_{2} R_{2} s - M_{1}^{2} s^{2} + R_{1} R_{3} + R_{1} R_{4} + R_{2} R_{3} + R_{2} R_{4}}\)

The transfer function can be simplified with the help of SymPy and the coefficients of the Laplace term can be collected. The transfer function at node 2 is:

H_sym = (U_sym[v2]/U_sym[v1]).simplify().collect(s)
H_sym

\(\displaystyle \frac{M_{1} R_{4} s}{R_{1} R_{3} + R_{1} R_{4} + R_{2} R_{3} + R_{2} R_{4} + s^{4} \left(C_{1} C_{2} L_{1} L_{2} R_{1} R_{4} - C_{1} C_{2} M_{1}^{2} R_{1} R_{4}\right) + s^{3} \left(C_{1} C_{2} L_{1} R_{1} R_{3} R_{4} + C_{1} C_{2} L_{2} R_{1} R_{2} R_{4} + C_{1} L_{1} L_{2} R_{1} - C_{1} M_{1}^{2} R_{1} + C_{2} L_{1} L_{2} R_{4} - C_{2} M_{1}^{2} R_{4}\right) + s^{2} \left(C_{1} C_{2} R_{1} R_{2} R_{3} R_{4} + C_{1} L_{1} R_{1} R_{3} + C_{1} L_{1} R_{1} R_{4} + C_{1} L_{2} R_{1} R_{2} + C_{2} L_{1} R_{3} R_{4} + C_{2} L_{2} R_{1} R_{4} + C_{2} L_{2} R_{2} R_{4} + L_{1} L_{2} - M_{1}^{2}\right) + s \left(C_{1} R_{1} R_{2} R_{3} + C_{1} R_{1} R_{2} R_{4} + C_{2} R_{1} R_{3} R_{4} + C_{2} R_{2} R_{3} R_{4} + L_{1} R_{3} + L_{1} R_{4} + L_{2} R_{1} + L_{2} R_{2}\right)}\)

The numerator and denominator can be extracted from the expression.

num_sym, denom_sym = fraction(H_sym)

Display the numerator:

num_sym

\(\displaystyle M_{1} R_{4} s\)

Display the denominator:

denom_sym

\(\displaystyle R_{1} R_{3} + R_{1} R_{4} + R_{2} R_{3} + R_{2} R_{4} + s^{4} \left(C_{1} C_{2} L_{1} L_{2} R_{1} R_{4} - C_{1} C_{2} M_{1}^{2} R_{1} R_{4}\right) + s^{3} \left(C_{1} C_{2} L_{1} R_{1} R_{3} R_{4} + C_{1} C_{2} L_{2} R_{1} R_{2} R_{4} + C_{1} L_{1} L_{2} R_{1} - C_{1} M_{1}^{2} R_{1} + C_{2} L_{1} L_{2} R_{4} - C_{2} M_{1}^{2} R_{4}\right) + s^{2} \left(C_{1} C_{2} R_{1} R_{2} R_{3} R_{4} + C_{1} L_{1} R_{1} R_{3} + C_{1} L_{1} R_{1} R_{4} + C_{1} L_{2} R_{1} R_{2} + C_{2} L_{1} R_{3} R_{4} + C_{2} L_{2} R_{1} R_{4} + C_{2} L_{2} R_{2} R_{4} + L_{1} L_{2} - M_{1}^{2}\right) + s \left(C_{1} R_{1} R_{2} R_{3} + C_{1} R_{1} R_{2} R_{4} + C_{2} R_{1} R_{3} R_{4} + C_{2} R_{2} R_{3} R_{4} + L_{1} R_{3} + L_{1} R_{4} + L_{2} R_{1} + L_{2} R_{2}\right)\)

Solve for the roots of the denominator.

#solve(denom_sym,s)

The code above was taking a long time, so it was interrupted and commented out.

3.20.1 Built a python dictionary of element values

element_value_keys = []
element_value_values = []

for i in range(len(df)):
    if df.iloc[i]['element'][0] == 'F' or df.iloc[i]['element'][0] == 'E' or df.iloc[i]['element'][0] == 'G' or df.iloc[i]['element'][0] == 'H':
        element_value_keys.append(var(df.iloc[i]['element'].lower()))
        element_value_values.append(df.iloc[i]['value'])

    else:
        element_value_keys.append(var(df.iloc[i]['element']))
        element_value_values.append(df.iloc[i]['value'])

element_values = dict(zip(element_value_keys, element_value_values))

The mutual inductance between L1 and L2 is calculated from the coupling coeeficient.

\(M = k\sqrt{L_1L_2}\)

# calculate the coupling constant from the mutual inductance
element_values[M1] = element_values[K1]*np.sqrt(element_values[L1] * element_values[L2])
print('mutual inductance, M1 = {:.9f}'.format(element_values[M1]))
mutual inductance, M1 = 0.000000150
element_values

\(\displaystyle \left\{ C_{1} : 1.0 \cdot 10^{-8}, \ C_{2} : 1.0 \cdot 10^{-8}, \ K_{1} : 0.15, \ L_{1} : 1.0 \cdot 10^{-6}, \ L_{2} : 1.0 \cdot 10^{-6}, \ M_{1} : 1.5 \cdot 10^{-7}, \ R_{1} : 1000.0, \ R_{2} : 0.5, \ R_{3} : 0.5, \ R_{4} : 10000.0, \ V_{1} : 1.0\right\}\)

3.21 Numeric solution

Substitute the element values into the equations and solve for unknown node voltages and currents. Need to set the current source, I1, to zero.

NE = NE_sym.subs(element_values)

Display the equations with numeric values.

temp = ''
for i in range(shape(NE.lhs)[0]):
    temp += '${:s} = {:s}$<br>'.format(latex(NE.rhs[i]),latex(NE.lhs[i]))

Markdown(temp)

\(0 = I_{V1} + 0.001 v_{1} - 0.001 v_{3}\)
\(0 = v_{2} \cdot \left(1.0 \cdot 10^{-8} s + 2.0001\right) - 2.0 v_{5}\)
\(0 = - 0.001 v_{1} + v_{3} \cdot \left(1.0 \cdot 10^{-8} s + 2.001\right) - 2.0 v_{4}\)
\(0 = I_{L1} - 2.0 v_{3} + 2.0 v_{4}\)
\(0 = I_{L2} - 2.0 v_{2} + 2.0 v_{5}\)
\(1.0 = v_{1}\)
\(0 = - 1.0 \cdot 10^{-6} I_{L1} s - 1.5 \cdot 10^{-7} I_{L2} s + v_{4}\)
\(0 = - 1.5 \cdot 10^{-7} I_{L1} s - 1.0 \cdot 10^{-6} I_{L2} s + v_{5}\)

Solve for voltages and currents in terms of Laplace variable s and displaying the transfer function.

U = solve(NE,X)
H = U[v2]/U[v1]
H

\(\displaystyle \frac{6.0 \cdot 10^{20} s}{391.0 s^{4} + 443010000.0 s^{3} + 8.0144391 \cdot 10^{16} s^{2} + 4.44114 \cdot 10^{22} s + 4.0022001 \cdot 10^{30}}\)

The numerator and denominator can be extracted from the expression.

num, denom = fraction(H)

Display the numerator:

num

\(\displaystyle 6.0 \cdot 10^{20} s\)

Display the denominator:

denom

\(\displaystyle 391.0 s^{4} + 443010000.0 s^{3} + 8.0144391 \cdot 10^{16} s^{2} + 4.44114 \cdot 10^{22} s + 4.0022001 \cdot 10^{30}\)

The roots of the denominator can be obtained with the solve function.

solve(denom,s)

\(\displaystyle \left[ -321649.557914015 - 10842889.0903087 i, \ -321649.557914015 + 10842889.0903087 i, \ -244859.393492634 - 9323420.66838014 i, \ -244859.393492634 + 9323420.66838014 i\right]\)

3.22 AC analysis

Solve equations a frequency of 1.491MHz or \(\omega\) equal to 9.3682292e6 radians per second, s = 9.3682292e6j.

NE_1rad_per_s = NE.subs({s:9.3682292e6j})

Display the equations with numeric values.

temp = ''
for i in range(shape(NE_1rad_per_s.lhs)[0]):
    temp += '${:s} = {:s}$<br>'.format(latex(NE_1rad_per_s.rhs[i]),latex(NE_1rad_per_s.lhs[i]))

Markdown(temp)

\(0 = I_{V1} + 0.001 v_{1} - 0.001 v_{3}\)
\(0 = v_{2} \cdot \left(2.0001 + 0.093682292 i\right) - 2.0 v_{5}\)
\(0 = - 0.001 v_{1} + v_{3} \cdot \left(2.001 + 0.093682292 i\right) - 2.0 v_{4}\)
\(0 = I_{L1} - 2.0 v_{3} + 2.0 v_{4}\)
\(0 = I_{L2} - 2.0 v_{2} + 2.0 v_{5}\)
\(1.0 = v_{1}\)
\(0 = - 9.3682292 i I_{L1} - 1.40523438 i I_{L2} + v_{4}\)
\(0 = - 1.40523438 i I_{L1} - 9.3682292 i I_{L2} + v_{5}\)

U_1rad_per_s = solve(NE_1rad_per_s,X)

table_header = ['unknown', 'mag','phase, deg']
table_row = []

for name, value in U_1rad_per_s.items():
    table_row.append([str(name),float(abs(value)),float(arg(value)*180/np.pi)])

print(tabulate(table_row, headers=table_header,colalign = ('left','decimal','decimal'),tablefmt="simple",floatfmt=('5s','.6f','.6f')))
unknown         mag    phase, deg
---------  --------  ------------
v1         1.000000      0.000000
v2         0.101262    -21.012442
v3         0.101577     -5.111761
v4         0.101240     -2.443264
v5         0.101378    -18.330736
I_V1       0.000899   -179.423099
I_L1       0.009469    -89.691641
I_L2       0.009486   -111.073602

3.23 AC Sweep

Looking at node 2 voltage, calculate the transfer function.

H = U[v2]/U[v1]
H

\(\displaystyle \frac{6.0 \cdot 10^{20} s}{391.0 s^{4} + 443010000.0 s^{3} + 8.0144391 \cdot 10^{16} s^{2} + 4.44114 \cdot 10^{22} s + 4.0022001 \cdot 10^{30}}\)

Get the coefficients of the numerator and denominator polynomials.

num, denom = fraction(H) #returns numerator and denominator

# convert symbolic to numpy polynomial
a = np.array(Poly(num, s).all_coeffs(), dtype=float)
b = np.array(Poly(denom, s).all_coeffs(), dtype=float)
system = (a, b) # system for circuit

Use the SciPy bode function to calculate the magnitude and phase response of the transfer function.

x = np.logspace(6, 6.5, 1000, endpoint=True)*2*np.pi
w, mag, phase = signal.bode(system, w=x) # returns: rad/s, mag in dB, phase in deg

Plot the results with Matplotlib.

fig, ax1 = plt.subplots()
ax1.set_ylabel('magnitude, dB')
ax1.set_xlabel('frequency, Hz')

plt.semilogx(w/(2*np.pi), mag,'-k')    # Bode magnitude plot

ax1.tick_params(axis='y')
plt.grid()

# instantiate a second y-axes that shares the same x-axis
ax2 = ax1.twinx()
plt.semilogx(w/(2*np.pi), phase,':',color='b')  # Bode phase plot

ax2.set_ylabel('phase, deg',color='b')
ax2.tick_params(axis='y', labelcolor='b')
#ax2.set_ylim((-5,25))

plt.title('Magnitude and phase response')
plt.show()

print('peak: {:.2f} dB at {:.3f} MHz'.format(mag.max(),w[np.argmax(mag)]/(2*np.pi)/1e6,))
peak: -19.89 dB at 1.490 MHz

3.24 Summary

In this chapter a walk through of the Python code that implements the symbolic MNA procedure was presented. The network equations in symbolic form were solved and the filter’s symbolic transfer function was obtained. Then the component values were substituted into the network equations and solved again. The node voltages were obtained at a discrete frequency and then an AC sweep was done for the transfer function.