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

view raw compute_loops.py hosted with ❤ by GitHub

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.