Monday, July 21, 2014

Revised nodal analysis

When must nodal analysis take place?
  1. When there are inductors with sufficient energy and that energy might be abruptly interrupted.
  2. When does an inductor have sufficient energy? When it has been in at least one non stiff loop in the previous iteration. This way, the circuit is power independent.
The code is below (click on "view raw" below the code box to see the code in a new window):

new_branch_events="yes"
while new_branch_events=="yes":
# Initialize currents in branches
# This is used to determine whether any
# inductor current change event occurs.
for c1 in range(len(branch_currents)):
branch_currents[c1]=0.0
# Set the branch currents equal to the inductor current
# and equal to any node with zero impedance
# The concept is that currents through resistance may need
# to be calculated to determine if their state will change
# This is specifically for freewheeling
for c1 in range(len(branch_params)):
# Check if branch has inductance and is not stiff
# In this case it behaves like a current source for
# the nodal analysis.
# Check if the branch is stiff
if stiff_ratio[c1]=="no":
# Check if the branch has inductance
if branch_params[c1][-1][0][1]:
#if branch_params[c1][-1][2]>0.001:
branch_currents[c1]=branch_params[c1][-1][2]
# Check if it is a zero impedance branch with a voltage
elif ((branch_params[c1][-1][0][0]==0.0) and (1.0 in branch_params[c1][-1][1] or -1.0 in branch_params[c1][-1][1])):
branch_currents[c1]=branch_params[c1][-1][2]
# Making a list of all the branches that have inductors
# but are not stiff.
inductor_list=[]
for c1 in range(len(branch_params)):
if branch_params[c1][-1][0][1]:
if stiff_ratio[c1]=="no":
inductor_list.append(c1)
# Determining is nodal analysis is required.
# Check if the inductor appears in a loop
# Check if any other branch of the loop is stiff.
# If that branch has become stiff in the previous iteration
# it is not counted as a stiff branch as that branch has
# experienced an event that is making the loop stiff.
# If at least one non stiff loop exists with the inductor
# the inductor has been in a non stiff loop and a nodal
# analysis needs to be performed
nodal_analysis_reqd="no"
if inductor_list:
for c1 in range(len(inductor_list)):
inductor_loop_stiff="yes"
for c2 in range(len(system_loop_map)):
if system_loop_map[c2][inductor_list[c1]]=="forward" or system_loop_map[c2][inductor_list[c1]]=="reverse":
loop_stiff="no"
for c3 in range(len(system_loop_map[c2])):
if system_loop_map[c2][c3]=="stiff_forward" or system_loop_map[c2][c3]=="stiff_reverse":
if branch_events[c3]=="no":
loop_stiff="yes"
if loop_stiff=="no":
inductor_loop_stiff="no"
if inductor_loop_stiff=="no":
nodal_analysis_reqd="yes"
# This function is to determine if there is a an event.
# In this case the event is - does the current through any inductor
# change instantaneously.
if nodal_analysis_reqd=="yes":
current_continuity(node_list, branch_params, stiff_ratio, node_voltage, branch_currents, admittance_matrix, source_vector, sys_mat_u, branch_events, "det_state")
# The branch currents from the above function
# contain the currents calculated from nodal analysis
# This is compared from the currents from loop analysis
# If there is a difference, there is an event.
continuity_event="no"
for c1 in range(len(branch_params)):
if branch_params[c1][-1][2]:
if (abs(branch_params[c1][-1][2]-branch_currents[c1])/branch_params[c1][-1][2]>0.01):
continuity_event="yes"
# Determine the state of nonlinear devices
# Using the currents from nodal analysis,
# it will check if any of the devices will start
# or stop conducting.
if continuity_event=="yes":
for comps in component_objects.keys():
component_objects[comps].determine_state(branch_currents, branch_params, branch_events)
# Check if there is a difference in the branch events as opposed
# to the previous nodal analysis. Keep performing nodal analysis
# until there are no new branch events
new_branch_events="no"
for c1 in range(len(branch_params)):
if not branch_events[c1]==branch_events_prev[c1]:
new_branch_events="yes"
# If there has been a new branch event, update the branch parameters
if new_branch_events=="yes":
for c1 in range(len(branch_params)):
branch_params[c1][-1][0][0]=0.0
branch_params[c1][-1][0][1]=0.0
for c2 in range(len(source_list)):
branch_params[c1][-1][1][c2]=0.0
for c1 in range(len(branch_params)):
for c2 in range(len(branch_params[c1][:-1])):
try:
comp_pos=csv_element(branch_params[c1][c2])
component_objects[comp_pos]
except:
pass
else:
component_objects[comp_pos].transfer_to_branch(branch_params[c1], source_list)
# Store the branch events
for c1 in range(len(branch_params)):
branch_events_prev[c1]=branch_events[c1]



The process is as follows:

  1. In the beginning, a new branch event is assumed by default whenever a branch event occurs.
  2. Look if there are any inductors that have previously been in non stiff loops.
  3. If so, nodal analysis is required.
  4. Following nodal analysis, check if there is a continuity event with branch currents from nodal analysis being different from the branch currents from previous iteration of loop analysis.
  5. If so, execute the determine_state function to check if any of the elements will change their state. If any of the elements do change state, the respective branches will have branch events.
  6. Check if there are any branch events that had not occurred in the previous iteration. If so, go back to step 2.
  7. This will repeat until no new branch events occur.

In nodal analysis, there has been one significant change. In order to perform nodal analysis, it is essential to calculate the node voltages. So the entire circuit has to be solved. However, the branch currents will be calculated only if:
  1. At a node, check if there has been an branch which has a "hard" event.
  2. A hard event is when a switch turns on or when it turns off while carrying a non negligible current. When a switch or a diode turns off when the current goes negative, that is a soft event. When a diode turns on when it is forward biased, it is a soft event.
  3. If a single branch incident at a node has a "hard" event, the currents of all the branches incident at the node will have to calculated. If not, the branch currents remain at their initial values - either zero or a current it is an inductor.
The code is below (click on "view raw" below the code box to see the code in a new window):
if func_purpose=="det_state":
# Need to differentiate between branches that are hard switched and
# those that are soft switched. If a node does not have a single branch
# that is hard switched, the branch currents of branches
# incident at that node do not get recalculated
freewheel_branches=[]
for c1 in range(len(list_of_nodes)):
hard_switched_branch="no"
for c2 in range(len(branch_info)):
if branch_info[c2][0]==list_of_nodes[c1] or branch_info[c2][-2]==list_of_nodes[c1]:
if br_events[c2]=="hard":
hard_switched_branch="yes"
# If the node is a short node, perform the
# same check with branches incident on all
# the short nodes of the same group
if hard_switched_branch=="no":
for c2 in range(len(shortnode_list)):
if c1 in shortnode_list[c2]:
for c3 in range(len(branch_info)):
for c4 in range(len(shortnode_list[c2])):
if branch_info[c3][0]==list_of_nodes[shortnode_list[c2][c4]] or branch_info[c3][-2]==list_of_nodes[shortnode_list[c2][c4]]:
if br_events[c3]=="hard":
hard_switched_branch="yes"
if hard_switched_branch=="yes":
for c2 in range(len(branch_info)):
if branch_info[c2][0]==list_of_nodes[c1] or branch_info[c2][-2]==list_of_nodes[c1]:
if c2 not in freewheel_branches:
freewheel_branches.append(c2)
else:
if func_purpose=="calc_currents":
freewheel_branches=[]
for c1 in range(len(branch_info)):
freewheel_branches.append(c1)



There has been another tag created called function_purpose. When determining state, function purpose will only calculate branch currents at nodes where there has been a "hard" event. But while calculating currents before the next loop analysis, all branch currents will be recalculated.

Determining switch and diode state

The code for determining switch state is below (click on "view raw" below the code box to view the code in a new window):


def determine_state(self, br_currents, sys_branches, sys_events):
""" Determines the state of the switch following an event
where the continuity of current through an inductor is
about to be broken."""
# Mark the position of the switch in sys_branches
for c1 in range(len(sys_branches)):
if csv_tuple(self.pos) in sys_branches[c1]:
branch_pos=c1
# Since branch current direction is by default considered
# positive when flowing away from the starting node
# If the branch current is negative, with the switch cathode
# closer towards the starting node, current direction is
# positive
## if br_currents[branch_pos]*self.resistor<-1.0:
## if sys_branches[branch_pos].index(self.polrty)<sys_branches[branch_pos].index(csv_tuple(self.pos)):
## if self.status=="off" and self.control_values[0]>0.0:
## sys_events[branch_pos]="yes"
## self.status="on"
# If the current direction is reverse, switch can never conduct
if br_currents[branch_pos]<0.0:
if sys_branches[branch_pos].index(self.polrty)>sys_branches[branch_pos].index(csv_tuple(self.pos)):
if self.status=="on":
sys_events[branch_pos]="yes"
self.status="off"
## if br_currents[branch_pos]*self.resistor>1.0:
## if sys_branches[branch_pos].index(self.polrty)>sys_branches[branch_pos].index(csv_tuple(self.pos)):
## if self.status=="off" and self.control_values[0]>0.0:
## sys_events[branch_pos]="yes"
## self.status="on"
if br_currents[branch_pos]>0.0:
if sys_branches[branch_pos].index(self.polrty)<sys_branches[branch_pos].index(csv_tuple(self.pos)):
if self.status=="on":
sys_events[branch_pos]="yes"
self.status="off"
# Update the value of resistance
if self.status=="off":
self.resistor=self.resistor_off
else:
self.resistor=self.resistor_on
return

The code for the determining diode state is below (click on "view raw" below the code box to view the code in a new window):


def determine_state(self, br_currents, sys_branches, sys_events):
""" Determines the state of the diode following an event
where the continuity of current through an inductor is
about to be broken."""
# Mark the position of the diode in the branches list
for c1 in range(len(sys_branches)):
if csv_tuple(self.pos) in sys_branches[c1]:
branch_pos=c1
# Since branch current direction is by default considered
# positive when flowing away from the starting node
# If the branch current is negative, with the diode cathode
# closer towards the starting node, current direction is
# positive
if br_currents[branch_pos]*self.resistor<-1.0:
#if br_currents[branch_pos]<-0.03:
if sys_branches[branch_pos].index(self.polrty)<sys_branches[branch_pos].index(csv_tuple(self.pos)):
if self.status=="off":
sys_events[branch_pos]="yes"
self.status="on"
# If current direction is reverse, diode can never conduct
if br_currents[branch_pos]<0.0:
if sys_branches[branch_pos].index(self.polrty)>sys_branches[branch_pos].index(csv_tuple(self.pos)):
if self.status=="on":
sys_events[branch_pos]="yes"
self.status="off"
if br_currents[branch_pos]*self.resistor>1.0:
#if br_currents[branch_pos]>0.03:
if sys_branches[branch_pos].index(self.polrty)>sys_branches[branch_pos].index(csv_tuple(self.pos)):
if self.status=="off":
sys_events[branch_pos]="yes"
self.status="on"
if br_currents[branch_pos]>0.0:
if sys_branches[branch_pos].index(self.polrty)<sys_branches[branch_pos].index(csv_tuple(self.pos)):
if self.status=="on":
sys_events[branch_pos]="yes"
self.status="off"
# Update the value of resistance
if self.status=="off":
self.resistor=self.resistor_off
else:
self.resistor=self.resistor_on
return
This function is called after performing nodal analysis to check if the inductor current has caused an non linear device to change state. The difference between the code for switch and diode is that a switch can only turn off during nodal analysis while a diode can both turn on and turn off. My guess is that the switch should turn on only in the update_value function. If for some reason, the switch were to turn off during the nodal analysis because the current through it was negative due to an inductor, the switch remains off and only a diode can turn on to freewheel the current.

The reason for doing this will be the next detailed blog post.

Stiff loop evaluation

Been a long time since I have posted. One of the reasons has been that I have made change after change in the solver. With every change, the circuit that I am working on gets fixed but a previously tested one breaks. Even now I am not totally sure until I test them all. But I will start posting code anyway.

To begin with, the stiff loop evaluation. The concept until now has been:
  1. Sort out the loops into stiff loops and non-stiff loops.
  2. With stiff loops, there are those loops that were non-stiff until the previous iteration and became stiff. So these will have a non-negligible current because a diode may have turned off causing the current to be slightly negative.
  3. So the currents of these stiff loops will have to be re-evaluated to brought back to values that correspond to their being stiff sloops with large resistances. A failure to do so will cause these loops to disrupt the calculation of the other stiff loops.
  4. Simple way to do this is to have another function. The code is below.


This is the function that calls the stiff equation solver which has not changed from before (click "view raw" link below the code box to see code in a new window):


def mat_ode_stiff(matrix_e, matrix_a, matrix_b, state_vector, input_vector, dt, ode_vectors, sys_loop_map, list_of_nodes, loop_info_stiff, branch_info_event, loop_events):
""" Solves thoses loops that are stiff """
# Initalize the loops that became stiff to zero
# so that they do not disrupt other newly formed
# stiff loops in the next calculation
for c1 in range(len(sys_loop_map)):
if loop_events[c1]=="yes":
state_vector[0].data[c1][0]=0.0
state_vector[1].data[c1][0]=0.0
for c2 in range(3):
for c1 in range(len(loop_info_stiff)-1, -1, -1):
if loop_events[c1]=="yes":
stiff_equation_solver(matrix_e, matrix_a, matrix_b, state_vector, input_vector, dt, ode_vectors, sys_loop_map, list_of_nodes, loop_info_stiff, branch_info_event, c1)
for c1 in range(len(loop_info_stiff)-1, -1, -1):
if loop_info_stiff[c1]=="yes":
stiff_equation_solver(matrix_e, matrix_a, matrix_b, state_vector, input_vector, dt, ode_vectors, sys_loop_map, list_of_nodes, loop_info_stiff, branch_info_event, c1)
return
view raw stiff_eqns.py hosted with ❤ by GitHub
There is a bit of guessing here as well. Out of "n" loops, let us suppose "m" loops (m<n) have turned stiff in the previous iteration. The remove_stiffness function makes the stiff loops into a upper triangular matrix with the first stiff branch of every stiff loop being present only in that loop. However, as the number of stiff branches can be larger than the number of stiff loops, the last few stiff branches may appear in multiple loops. As an example, take a look at these loops:

stiff_forward no no no no no no stiff_reverse reverse reverse no no reverse no reverse reverse no reverse
no stiff_forward no no no no no stiff_forward forward forward no no forward no forward forward no forward
no no stiff_forward no no no no stiff_reverse no reverse no no no reverse reverse reverse reverse no
no no no stiff_reverse no no no stiff_reverse no reverse no no no reverse reverse reverse reverse no
no no no no stiff_reverse no no stiff_reverse no no forward forward no no reverse reverse no reverse
no no no no no stiff_reverse no stiff_forward no no reverse reverse no no forward forward no forward
no no no no no no stiff_reverse stiff_reverse no no no no no no no no no no

All the loops are stiff and the 8th branch appears in all the stiff loops. And for all you know, this 8th branch could be the branch that became recently stiff. So there is a possibility that these newly formed loops could disrupt the calculation of other new formed loops. So thus there is the initialization of newly formed stiff loop currents to zero.

The current of a stiff loop is decided only by the voltage in that loop and the resistance of that loop. The loop has no dynamics di/dt. So repeatedly calculating the loops will not change the loop currents unless the currents of one of the other stiff loops has changed. Which is why if the newly formed stiff loops are recalculated a number of times (thrice for now), the values should stabilize. In any case without a di/dt, they should not blow up out of bounds.

Following this, when it comes to calculating loop currents from branches, there were two parts to that function. The code is below (click "view raw" link below the code box to see code in a new window):


def compute_loop_currents(matrix_e, matrix_a, matrix_b, dt, sys_loops, branch_info, stiff_info, sys_loop_map, src_list, state_vector):
""" Calculates the loop currents after an event from branch currents."""
# Make a list of the nonstiff loops
# Because stiff loops will essentially have a negligible
# current.
nonstiff_loops=[]
for c1 in range(len(sys_loop_map)):
# Check if loop c1 has a stiff branch
is_loop_stiff="no"
for c2 in range(len(branch_info)):
if sys_loop_map[c1][c2]=="stiff_forward" or sys_loop_map[c1][c2]=="stiff_reverse":
is_loop_stiff="yes"
if is_loop_stiff=="no":
nonstiff_loops.append(c1)
# To begin with find out which of the branches
# occur only in one loop. This will make
# computations easier as those branch currents
# will automatically become the initial values
# for the loop currents in the next loop analysis.
# A list of loops and branches with this info
# It is a corresponding mapping.
single_loops=[]
single_branches=[]
# Iterate through the nonstiff loops.
for c1 in nonstiff_loops:
# Iterate through the branches
for c2 in range(len(branch_info)):
# Check if the corresponding system map is a "yes"
# Which means branch exixts.
if sys_loop_map[c1][c2]=="forward" or sys_loop_map[c1][c2]=="reverse":
# Look for the branch in all other nonstiff loops
# A single occurance elsewhere will mean loop exists
# So start with the default "no"
does_branch_occur="no"
for c3 in nonstiff_loops:
if not c1==c3:
if sys_loop_map[c3][c2]=="forward" or sys_loop_map[c3][c2]=="reverse":
does_branch_occur="yes"
#print does_branch_occur, c1, c2
#print "Check"
if does_branch_occur=="no":
# Next check is both the loop and the branch have not been
# found before to prevent clashes.
if ((c1 not in single_loops) and (c2 not in single_branches)):
single_loops.append(c1)
single_branches.append(c2)
# If the number of loops where a branch occurs in only one
# loop is less than the number of nonstiff loops, row
# manipulations will have to be done to the remaining loops
# so that one branch current can be taken as the loop current
if (len(single_loops)<len(nonstiff_loops)):
# Iterate through the nonstiff loops
for c1 in nonstiff_loops:
# Make sure the loop has not been found.
# Don't manipulate those which have been isolated.
if c1 not in single_loops:
# Look through the branches
for c2 in range(len(branch_info)):
# Again check if the branch has not been isolated.
if c2 not in single_branches:
# If a loop exists. So this will find the
# first branch that exists in a loop
if sys_loop_map[c1][c2]=="forward" or sys_loop_map[c1][c2]=="reverse":
# Look at the remaining nonstiff loops.
for c3 in nonstiff_loops:
if ((c1!=c3) and (c3 not in single_loops)):
# If branch exists in any other nonstiff loop
# Need to manipulate the loops and for that
# need to find out the sense in which the branches
# are in the loops.
if sys_loop_map[c3][c2]=="forward" or sys_loop_map[c3][c2]=="reverse":
if not sys_loop_map[c1][c2]==sys_loop_map[c3][c2]:
for c4 in range(len(sys_loop_map[c1])):
if sys_loop_map[c1][c4]=="forward" and sys_loop_map[c3][c4]=="reverse":
sys_loop_map[c3][c4]="no"
elif sys_loop_map[c1][c4]=="reverse" and sys_loop_map[c3][c4]=="forward":
sys_loop_map[c3][c4]="no"
elif sys_loop_map[c1][c4]=="forward" and sys_loop_map[c3][c4]=="no":
sys_loop_map[c3][c4]="forward"
elif sys_loop_map[c1][c4]=="reverse" and sys_loop_map[c3][c4]=="no":
sys_loop_map[c3][c4]="reverse"
else:
for c4 in range(len(sys_loop_map[c1])):
if sys_loop_map[c1][c4]=="forward" and sys_loop_map[c3][c4]=="forward":
sys_loop_map[c3][c4]="no"
elif sys_loop_map[c1][c4]=="reverse" and sys_loop_map[c3][c4]=="reverse":
sys_loop_map[c3][c4]="no"
elif sys_loop_map[c1][c4]=="forward" and sys_loop_map[c3][c4]=="no":
sys_loop_map[c3][c4]="reverse"
elif sys_loop_map[c1][c4]=="reverse" and sys_loop_map[c3][c4]=="no":
sys_loop_map[c3][c4]="forward"
# If both loop and branch have not been found add them
if ((c1 not in single_loops) and (c2 not in single_branches)):
single_loops.append(c1)
single_branches.append(c2)
# Now to set the loop currents equal to branch currents.
# Take every loop and branch in the single_loop and single_branch lists
# Since they are a one-to-one mapping, just equate the state vectors to
# to the branch currents.
for c1 in range(len(single_loops)):
loop_row=single_loops[c1]
if sys_loop_map[loop_row][single_branches[c1]]=="forward":
state_vector[0].data[single_loops[c1]][0]=branch_info[single_branches[c1]][-1][2]
state_vector[1].data[single_loops[c1]][0]=branch_info[single_branches[c1]][-1][2]
else:
state_vector[0].data[single_loops[c1]][0]=-branch_info[single_branches[c1]][-1][2]
state_vector[1].data[single_loops[c1]][0]=-branch_info[single_branches[c1]][-1][2]
# This computation is to set the initial values of loop currents
# for stiff loops. Since the stiff loops have been made upper triangular
# the first occurrance of a stiff branch will mean the loop current is
# equal to the stiff branch.
# for c1 in range(len(sys_loop_map)):
# # Go through all the stiff loops.
# if c1 not in nonstiff_loops:
# stiff_branch_found="no"
# c2=0
# # Check if a branch is stiff.
# while stiff_branch_found=="no" and c2<len(sys_loop_map[c1]):
# if sys_loop_map[c1][c2]=="stiff_forward" or sys_loop_map[c1][c2]=="stiff_reverse":
# stiff_branch_found="yes"
# # At the first occurrance of stiff branch, initalize
# # the loop current equal to the stiff branch current.
# if sys_loop_map[c1][c2]=="stiff_forward":
# state_vector[0].data[c1][0]=branch_info[c2][-1][2]
# state_vector[1].data[c1][0]=branch_info[c2][-1][2]
# elif sys_loop_map[c1][c2]=="stiff_reverse":
# state_vector[0].data[c1][0]=-branch_info[c2][-1][2]
# state_vector[1].data[c1][0]=-branch_info[c2][-1][2]
# Now iterate through all the loops previous to this stiff loop.
#for c3 in range(c1):
#if c3 not in nonstiff_loops:
# If the same stiff branch is found in them, the loop currents of
# those loops will have to deducted. This is because the sum of loop
# currents is the branch current. So if the stiff branch occurs in
# multiple loops, it is the sum of all the loop currents that will
# be the current through the stiff branch.
#if sys_loop_map[c3][c2]=="stiff_forward":
#state_vector[0].data[c1][0]-=state_vector[0].data[c3][0]
#state_vector[1].data[c1][0]-=state_vector[1].data[c3][0]
#elif sys_loop_map[c3][c2]=="stiff_reverse":
#state_vector[0].data[c1][0]+=state_vector[0].data[c3][0]
#state_vector[1].data[c1][0]+=state_vector[1].data[c3][0]
# c2+=1
for c1 in range(len(sys_loop_map)-1, -1, -1):
loop_exists="no"
for c2 in range(len(sys_loop_map[c1])):
if not sys_loop_map[c1][c2]=="no":
loop_exists="yes"
if loop_exists=="no":
del sys_loop_map[c1]
# Re-initialize the matrices A, B, and E
matrix_a.zeros(len(sys_loop_map),len(sys_loop_map))
matrix_e.zeros(len(sys_loop_map),len(sys_loop_map))
matrix_b.zeros(len(sys_loop_map),matrix_b.columns)
# Recalculate the system matrices for the new loops.
recalculate_sys_matrices(sys_loop_map, branch_info, matrix_a, matrix_e, matrix_b)
return

As can be seen, the entire block of code that would compute currents of stiff loops has been commented out. The reason is that if the currents of stiff loops are being calculated from the stiff equation solver, there is no need to repeat the process. On the contrary repeating the process has the capacity to disrupt the loop current values.