Releasing version 0.2.5

Releasing the version 0.2.5:
http://sourceforge.net/projects/pythonpowerelec/

Next update

A loooong time since I posted. Made a huge number of changes to the simulator. It is now working for a circuit like a three-phase diode rectified feeding a dc capacitor that is in turn the dc bus for a three-phase 6 switch inverter that has an LC filter. So I feel most of the major bugs have been ironed out.

Anyway, I will be releasing the next version of the simulator tomorrow. But before I end the day, I would like to ask a question to some of the readers of the blog.

Since, documenting this project has been a major pain - this blog is the only thing I have and I am quite glad I did at least that. But this blog too is a mess and I need something a little more systematic. So I was thinking of writing a book. Does anyone know of any book on circuit simulation that actually describes how a simulator can be built? Or if anyone would be interested in reading such a book, what would your expectations be from such a book and why? Out of curiosity, because you are trying something similar, or because you are not happy with existing software?

Would be great to hear any feedback. Please send an email to pythonpowerelectronics@gmail.com.

Voltage measurement

Have been working furiously on modifying the simulator for larger circuits. But found there is a fatal error in voltage measurement. Take a look at the test circuit for a three-phase diode rectifier.

The voltmeter is modelled as a large resistance which does not disturb the current in the remaining circuit. The problem is that the branch containing the voltmeter is a stiff branch that appears in the loop equations. These equations appear randomly. So in the above circuit, when the simulation begins, suppose diodes D1 and D6 conduct. This means voltage vac appears across the capacitor Cdc and the voltmeter Vo. What would need to happen is that the voltmeter loop should be written with Cdc and therefore as Cdc charges, Vo output voltage should increase. But what is Vo appears in the loop with vac, D1 and D6? This is perfectly possible the way the program works now. This means Vo will jump to the voltage vac even though Cdc is charging gradually.

So Vo does not reflect the voltage that appears across the nodes but rather on the voltage impressed on it through the loop it appears on. And this inevitably happens whenever the circuit becomes large and loop equations become more and more random. Also, if the voltmeter were connected across a largely inductive circuit or for that matter over two fairly distant nodes in the circuit. What then?

The only error free way to solve this seems to be to do a nodal analysis with only the voltmeters for the nodes where they exist. This calls for a fairly large additional piece of code that will run at simulation time step.

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]

view raw nodal_analysis_check.py hosted with ❤ by GitHub

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)

view raw freewheelbranches.py hosted with ❤ by GitHub

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

view raw determine_state.py hosted with ❤ by GitHub

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

view raw determinestate_diode.py hosted with ❤ by GitHub

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

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.

Restarting

Been a long time since I posted. Have been busy with building my hardware setup for my research and now they are almost working the way I want. With this done, I now will have to divert some of my time towards this project.

Spent the past week or so reading most of the code all over again. Another task has been portability to Python 3. Have to look into compatibility issues with respect to that. But before that to complete a tested power electronics library.

One of the problems with combining mesh analysis and nodal analysis was that if both are performed whenever an event occurs, freewheeling can be said to occur even when it is not meant to occur. For example, in a three-phase diode bridge rectifier fed by a source with a finite inductance, the turning off of a diode can be seen as cause for freewheeling of the inductor current and causing the other diode in the leg to freewheel.

So, my guess is there comes the need to distinguish between nonlinear element events -  a hard switched event or a soft switched event. A soft switched event is when a device (diode or switch) turns off when the current through it becomes negative. In that case, freewheeling should not happen. However, when a switch is turned off, and the current through it (irrespective of the magnitude) is broken, freewheeling needs to occur if there was an inductor in a branch connected to the node.

That is the immediate next step.

Hardware setup

Been a while since I posted about the simulator. Have been crazy busy with my day job but am happy to announce, that my microgrid setup is finally ready!

Inverter and control board with sensors and microcontroller:

Intelligent power module with dc bus capacitors:

Hopefully I will be able to work on the simulator again before my trip to Japan.

Initializing new stiff equations

Another problem that arises when there are a large number of nonlinear devices is that nonstiff loops become stiff loops but their current has still not become negligible because the diode has turned off at the previous iteration. For example, current through a diode is -0.1A causing the diode to turn off and it's resistance to change to megaohms. When this stiff loops interacts with other stiff loops and the currents of those stiff loops are calculated, all other stiff loops end up having currents of close magnitudes of 0.1 A. Result: entire simulation is disrupted.

Two things need to be done:
1. Need to identify which loops have just become stiff.
2. Calculate those stiff loops before running the ode solver.

The code to identify the newly formed stiff loops is:

	# Check if the system is stiff and recalculate if it is.
	# The stiff loops are isolated to the minimum number of loops
	# So that the system dynamics are captured as far as possible.
	remove_stiffness(sys_mat_e, sys_mat_a, sys_mat_b, dt, branch_params, stiff_ratio, system_loop_map, source_list, [curr_state_vec, next_state_vec], ode_var, loop_stiff_info)
	# This list if of those loops that were nonstiff
	# and have become stiff in the previous calculation
	loop_event=[]
	for c1 in range(len(system_loop_map)):
	is_loop_stiff="no"
	loop_event.append("no")
	for c2 in range(len(system_loop_map[c1])):
	if system_loop_map[c1][c2]=="stiff_forward" or system_loop_map[c1][c2]=="stiff_reverse":
	is_loop_stiff="yes"
	if loop_stiff_info[c1]=="no" and is_loop_stiff=="yes":
	loop_event[c1]="yes"
	else:
	loop_event[c1]="no"
	# This is just a list of which loops are stiff and which are nonstiff
	if is_loop_stiff=="yes":
	loop_stiff_info[c1]="yes"
	else:
	loop_stiff_info[c1]="no"
	# Recalculate the newly formed stiff loops.
	mat_ode_stiff(sys_mat_e, sys_mat_a, sys_mat_b, [curr_state_vec, next_state_vec], sys_mat_u, t-t_ode_prev, ode_var, system_loop_map, node_list, loop_stiff_info, branch_events, loop_event)

view raw init_stiff_loops.py hosted with ❤ by GitHub

And running the stiff loop ode:

	def 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, ode_row):
	""" Solving stiff loops that are algebraic equations only """
	if not (matrix_a.data[ode_row][ode_row]):
	#The variable has no dynamics and can't even be calculated statically!
	# May be due to a redundant loop.
	ode_vectors[4].data[ode_row][0]=0.0
	state_vector[0].data[ode_row][0]=0.0
	state_vector[1].data[ode_row][0]=0.0
	else:
	ode_vectors[4].data[ode_row][0]=0.0
	state_vector[0].data[ode_row][0]=0.0
	state_vector[1].data[ode_row][0]=0.0
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	state_vector[0].data[ode_row][0]+=matrix_b.data[ode_row][0]*input_vector
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	for c2 in range(matrix_b.columns):
	state_vector[0].data[ode_row][0]+=matrix_b.data[ode_row][c2]*input_vector.data[c2][0]
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	for c2 in range(matrix_a.columns):
	if not (ode_row==c2):
	state_vector[0].data[ode_row][0]-=matrix_a.data[ode_row][c2]*state_vector[1].data[c2][0]
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	for c2 in range(matrix_e.columns):
	state_vector[0].data[ode_row][0]-=matrix_e.data[ode_row][c2]*ode_vectors[4].data[c2][0]/dt
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	state_vector[0].data[ode_row][0]=state_vector[0].data[ode_row][0]/matrix_a.data[ode_row][ode_row]
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	return
	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 """
	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_loop_solver.py hosted with ❤ by GitHub

One other thing still needs to be done. What if the newly formed stiff loops themselves are interacting with each other through a stiff branch? Now even this ode solver would fail because they would disrupt each other. The only thing to do would be do perform an upper triangularization of the newly formed stiff loops with respect to the newly formed stiff branches and solve them backwards. Will do this later.

Nodal analysis

There might have been a mistake in the previous current initialization of branch currents. The branches are current sources if:
1. They have inductors and are not stiff.
2. They have voltages and have zero impedance.

If the branches are stiff or are purely resistive with or without a voltage, their currents will need to be calculated by nodal analysis. The code to set the branch currents is:

	# Initializing the branch currents for the next set of nodal analysis
	for c1 in range(len(branch_params)):
	branch_currents[c1]=0.0
	for c1 in range(len(branch_params)):
	# Check if the branch is stiff
	if stiff_ratio[c1]=="no":
	# Check if the branch has inductance
	if branch_params[c1][-1][0][1]:
	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]
	# This second run of the function is to determine the new
	# branch currents because of any possible change in the state
	# of nonlinear devices.
	current_continuity(node_list, branch_params, stiff_ratio, node_voltage, branch_currents, admittance_matrix, source_vector, sys_mat_u)

view raw nodal_current_init.py hosted with ❤ by GitHub

The actual nodal analysis is:

	# Start iterating through the nodes
	for c1 in range(len(list_of_nodes)):
	# Iterate through all the branches.
	for c2 in range(len(branch_info)):
	# Look for the node as the starting or ending node in each branche
	if ((list_of_nodes[c1]==branch_info[c2][0]) or (list_of_nodes[c1]==branch_info[c2][-2])):
	if (list_of_nodes[c1]==branch_info[c2][0]):
	start_node_pos=list_of_nodes.index(branch_info[c2][0])
	# Look for the start node in the short node list
	# If it is a short node, the first short node
	# in that list to which the current short node
	# belongs will be the start node pos.
	# this is because all short nodes are the same.
	for c3 in range(len(shortnode_list)):
	if start_node_pos in shortnode_list[c3]:
	start_node_pos=shortnode_list[c3][0]
	# Mark the position of the other node
	end_node_pos=list_of_nodes.index(branch_info[c2][-2])
	# Again check if it is a short node.
	for c3 in range(len(shortnode_list)):
	if end_node_pos in shortnode_list[c3]:
	end_node_pos=shortnode_list[c3][0]
	# The default direction of current is taken to be
	# away from the starting node - start to end.
	branch_src_dir=1.0
	else:
	start_node_pos=list_of_nodes.index(branch_info[c2][-2])
	for c3 in range(len(shortnode_list)):
	if start_node_pos in shortnode_list[c3]:
	start_node_pos=shortnode_list[c3][0]
	end_node_pos=list_of_nodes.index(branch_info[c2][0])
	for c3 in range(len(shortnode_list)):
	if end_node_pos in shortnode_list[c3]:
	end_node_pos=shortnode_list[c3][0]
	branch_src_dir=-1.0
	# If resistance of the branch is non-zero
	# if branch_info[c2][-1][0][0]:
	# # If the branch has no inductance or is a stiff branch
	# # in which case inductance is negligible compared
	# # to the resistance and the branch is stiff.
	# if (branch_info[c2][-1][0][1]==0.0):
	# mho_matrix[start_node_pos][start_node_pos]+=1.0/branch_info[c2][-1][0][0]
	# mho_matrix[start_node_pos][end_node_pos]-=1.0/branch_info[c2][-1][0][0]
	#
	# # Similarly, any voltage source that exists
	# # will be treated as positive if it forces a
	# # a current away the starting node.
	# # The way branch_info contains these sources is
	# # compatible to the definition
	# # src_vector is on the RHS of the equation,
	# # so thus the -ve sign
	# for c3 in range(len(branch_info[c2][-1][1])):
	# src_vector[start_node_pos]-=branch_src_dir*branch_info[c2][-1][1][c3]*sys_inputs.data[c3][0]/branch_info[c2][-1][0][0]
	# # Check if inductance is non-zero
	# if branch_info[c2][-1][0][1]:
	# # Check if branch is not stiff
	# if branch_stiff[c2]=="no":
	# # Add the current to the src_vector as a current source
	# src_vector[start_node_pos]-=branch_src_dir*br_currents[c2]
	# # Check if is a zero impedance branch
	# if ((branch_info[c2][-1][0][0]==0.0) and (branch_info[c2][-1][0][1]==0.0)):
	# if ((1.0 in branch_info[c2][-1][1]) or (-1.0 in branch_info[c2][-1][1])):
	# src_vector[start_node_pos]-=branch_src_dir*br_currents[c2]
	# Check if branch is not stiff
	if branch_stiff[c2]=="no":
	# Check if inductance is non-zero
	if branch_info[c2][-1][0][1]:
	src_vector[start_node_pos]-=branch_src_dir*br_currents[c2]
	# Check if is a zero impedance branch with a voltage
	elif ((branch_info[c2][-1][0][0]==0.0) and (1.0 in branch_info[c2][-1][1] or -1.0 in branch_info[c2][-1][1])):
	src_vector[start_node_pos]-=branch_src_dir*br_currents[c2]
	# Check if it is a resistive branch
	if ((branch_info[c2][-1][0][0]) and (branch_info[c2][-1][0][1]==0.0)):
	mho_matrix[start_node_pos][start_node_pos]+=1.0/branch_info[c2][-1][0][0]
	mho_matrix[start_node_pos][end_node_pos]-=1.0/branch_info[c2][-1][0][0]
	# Similarly, any voltage source that exists
	# will be treated as positive if it forces a
	# a current away the starting node.
	# The way branch_info contains these sources is
	# compatible to the definition
	# src_vector is on the RHS of the equation,
	# so thus the -ve sign
	for c3 in range(len(branch_info[c2][-1][1])):
	src_vector[start_node_pos]+=branch_src_dir*branch_info[c2][-1][1][c3]*sys_inputs.data[c3][0]/branch_info[c2][-1][0][0]

view raw nodal_analysis.py hosted with ❤ by GitHub

Other than that, the code is quite the same as the previous case. This code was tested with a three phase diode bridge rectifier and was found to be working OK. The previously released code was breaking.

Rewriting the loop finder algorithm

Previously, the system loops were added if they met the following criteria:

1. Did not pass through the same branch or same node more than once.
2. Last node was the same as the first node.

During the loop finder, no check was imposed on whether the loops were independent and so they were several loops which were independent.The way to deal with this was to let the simulation run and after the first run, the dependent loops would become void and be deleted. A simpler method could be used to find the independent loops:

1. Did not pass through the same branch or same node more than once.
2. Last node was the same as the first node.
3. The system branches must be contained in at least one loop.
4. As the loops are added, a loop will be added if and only if it passes through a new branch that the previously found loops have not passed through yet.

The 4th check simplifies the process and at the first iteration, the loops found are independent loops and equal to the number (Branches-Nodes+1).

Here is the code (click on "view raw" to see it in a new window):

	# Determining the loops
	loop_list=[]
	# Loop map is the list of all the node pairs
	# that have branches between them.
	loop_map_list=[]
	for c1 in range(len(branch_map)):
	for c2 in range(c1+1, len(branch_map)):
	if branch_map[c1][c2]:
	if [c1, c2] not in loop_map_list:
	loop_map_list.append([c1, c2])
	# A special check for parallel branches between nodes
	# Add the pair of nodes directly to the loop_list
	for c1 in range(len(loop_map_list)):
	if len(branch_map[loop_map_list[c1][0]][loop_map_list[c1][1]])>1:
	if not check_loop_repeat([loop_map_list[c1], loop_map_list[c1]], loop_list):
	loop_list.append([loop_map_list[c1], loop_map_list[c1]])
	while loop_map_list:
	# loop_iter is the list that iteratively checks the
	# loops for continuity and closing.
	loop_iter=[]
	loop_iter.append([loop_map_list[0]])
	find_loop(branch_map, node_list, loop_list, loop_iter, \
	loop_map_list[0], loop_map_list)
	# Delete all the node pairs (branches found)
	# from the loop_map_list so that new loops only look for
	# new branches.
	for c1 in range(len(loop_list)):
	for c2 in range(len(loop_map_list)-1, -1, -1):
	if (loop_map_list[c2] in loop_list[c1]):
	del loop_map_list[c2]
	elif ([loop_map_list[c2][1], loop_map_list[c2][0]] in loop_list[c1]):
	del loop_map_list[c2]
	# Remove any repitions in loop_list
	for c1 in range(len(loop_list)-1):
	for c2 in range(len(loop_list)-1, c1, -1):
	if loop_list[c1]==loop_list[c2]:
	del loop_list[c2]
	loop_count=len(loop_list)

view raw loop_finder.py hosted with ❤ by GitHub

 The rest of the loop finder code is pretty much the same. The only difference is that the parallel branches are found in the main block above so they have been eliminated from the other functions. Anyway, the rest of the function code is below (click on "view raw" to see it in a new window):

	def loop_copy(loop_inp):
	""" Will return a copy of a loop list
	Used when a change needs to be made."""
	loop_op=[]
	for c1 in range(len(loop_inp)):
	row_vector=[]
	for c2 in range(len(loop_inp[c1])):
	row_vector.append(loop_inp[c1][c2])
	loop_op.append(row_vector)
	return loop_op
	def check_loop_repeat(lp_iter, lp_list):
	""" Will return 1 if the loop already
	exists if the loop_list found so far."""
	# Make a copy of the entire loop list
	# found so far. Just in case.
	list_cmp=loop_copy(lp_list)
	# As a default, loop is not repeated
	# A single instance of repitition will
	# set it to 1.
	lp_repeat=0
	# Go through every loop found
	# so far
	for c1 in range(len(list_cmp)):
	# Make a copy of the loop being checked
	iter_cmp=loop_copy(lp_iter)
	# Move in the reverse direction of the loop
	# This is because elements are deleted
	# as they are found.
	for c2 in range(len(iter_cmp)-1, -1, -1):
	# Check if the element is found or if the
	# symmetrical element is found in any of
	# the loops found so far.
	# Because they are the same.
	if ([iter_cmp[c2][0], iter_cmp[c2][1]] in list_cmp[c1]) or \
	([iter_cmp[c2][1], iter_cmp[c2][0]] in list_cmp[c1]):
	# If so, remove the element
	iter_cmp.remove(iter_cmp[c2])
	# If all element have been deleted
	# it means the loop exists
	if not iter_cmp:
	lp_repeat=1
	return lp_repeat
	def loop_addition(br_map, nd_list, lp_list, curr_lp_iter, lp_update, curr_elem, main_loops):
	""" Take a new element of br_map in any direction.
	Check for a parallel branch at that element.
	Add that element if not already there in the temporary loop."""
	row=curr_elem[0]
	col=curr_elem[1]
	# Check if there is an element
	if br_map[row][col]:
	# Temp list to make copies
	# of loops found
	row_vec=[]
	for item in curr_lp_iter:
	row_vec.append(item)
	# Check if an element has been
	# encountered before
	# not going back and forth that is
	if not (([row, col] in row_vec) or ([col, row] in row_vec)):
	# Add the element found
	lp_update.append(row_vec)
	lp_update[main_loops].append([row, col])
	# Update the main loops counter
	main_loops+=1
	# If an element has not been found
	# or is a duplicate, lp_update
	# won't contain it.
	#all_loops=all_loops+main_loops
	return main_loops
	def loop_closing(br_map, lp_list, nd_list, lp_update, c1):
	""" The check imposed is whether the loop has the same
	elements as any other loops already found in lp_list. """
	# Check if the loop is new even
	# if all nodes have been passed through.
	lp_exist="found"
	if not check_loop_repeat(lp_update[c1], lp_list):
	lp_exist="not_found"
	if lp_exist=="not_found":
	# Add that loop to loop list
	lp_list.append(lp_update[c1])
	# Remove the loop from the temp
	# variable.
	del lp_update[c1]
	return
	def loop_horiz(br_map, nd_list, lp_list, lp_iter, elem, lp_map_list):
	""" Moves horizontally in a loop find.
	Looks for every element in a particular column of br_map.
	Each element is added onto the exiting loops. """
	# lp_list is the loops found.
	# lp_iter is the list of loops as they are being identified.
	# Those that are loops will be added to lp_list.
	no_nodes=len(br_map)
	start_row=elem[0]
	start_col=elem[1]
	# Temp list to make copies
	# of exisiting lists
	# if elements exist on that row
	lp_update=[]
	# Counter for number
	# of elements found in the row
	c2=0
	for c1 in range(len(lp_iter)):
	# Set loop row and counter
	# to the latest element
	# in lp_iter
	loop_row=lp_iter[c1][-1][0]
	loop_column=lp_iter[c1][-1][1]
	# Start from element in next column
	# to end of matrix
	#for c3 in range(loop_column+1, no_nodes):
	for c3 in range(0, no_nodes):
	curr_elem=[loop_row, c3]
	c2=loop_addition(br_map, nd_list, lp_list, lp_iter[c1], lp_update, curr_elem, c2)
	# for c1 in range(len(lp_update)-1, -1, -1):
	any_loop_closes="no"
	c1=len(lp_update)-1
	while any_loop_closes=="no" and c1>=0:
	# End condition
	# Latest element has ending or starting node
	# Same as last element of first branch
	last_elem_frst=br_map[start_row][start_col][0][-1]
	frst_elem_curr=br_map[lp_update[c1][-1][0]][lp_update[c1][-1][1]][0][0]
	last_elem_curr=br_map[lp_update[c1][-1][0]][lp_update[c1][-1][1]][0][-1]
	if (frst_elem_curr==last_elem_frst or \
	last_elem_curr==last_elem_frst):
	loop_closing(br_map, lp_list, nd_list, lp_update, c1)
	any_loop_closes="yes"
	c1-=1
	if any_loop_closes=="no":
	# lp_iter will be the same as lp_update
	lp_iter=[]
	for c1 in range(len(lp_update)):
	lp_iter.append(lp_update[c1])
	else:
	lp_iter=[]
	# lp_iter contains ongoing loops
	# Closed loops are moved to lp_list
	# Broken loops are dropped
	return lp_iter
	def loop_vert(br_map, nd_list, lp_list, lp_iter, elem, lp_map_list):
	""" Moves vertically in a loop find.
	Looks for every element in a particular column of br_map.
	Each element is added onto the exiting loops. """
	# lp_list is the loops found.
	# lp_iter is the list of loops as they are being identified.
	# Those that are loops will be added to lp_list.
	no_nodes=len(br_map)
	start_row=elem[0]
	start_col=elem[1]
	# Temp list to make copies
	# of exisiting lists
	# if elements exist on that column
	lp_update=[]
	# Counter for number
	# of elements found in the column
	c2=0
	for c1 in range(len(lp_iter)):
	# Set loop row and counter
	# to the latest element
	# in lp_iter
	loop_row=lp_iter[c1][-1][0]
	loop_column=lp_iter[c1][-1][1]
	# Start from element from first row
	# to end of column
	for c3 in range(0, no_nodes):
	curr_elem=[c3, loop_column]
	c2=loop_addition(br_map, nd_list, lp_list, lp_iter[c1], lp_update, curr_elem, c2)
	# for c1 in range(len(lp_update)-1, -1, -1):
	any_loop_closes="no"
	c1=len(lp_update)-1
	while any_loop_closes=="no" and c1>=0:
	# End condition
	# Latest element has ending or starting node
	# Same as last element of first branch
	last_elem_frst=br_map[start_row][start_col][0][-1]
	frst_elem_curr=br_map[lp_update[c1][-1][0]][lp_update[c1][-1][1]][0][0]
	last_elem_curr=br_map[lp_update[c1][-1][0]][lp_update[c1][-1][1]][0][-1]
	if (frst_elem_curr==last_elem_frst or \
	last_elem_curr==last_elem_frst):
	loop_closing(br_map, lp_list, nd_list, lp_update, c1)
	any_loop_closes="yes"
	c1-=1
	if any_loop_closes=="no":
	# lp_iter will be the same as lp_update
	lp_iter=[]
	for c1 in range(len(lp_update)):
	lp_iter.append(lp_update[c1])
	else:
	lp_iter=[]
	# lp_iter contains ongoing loops
	# Closed loops are moved to lp_list
	# Broken loops are dropped
	return lp_iter
	def find_loop(br_map, nd_list, lp_list, lp_iter, elem, lp_map_list):
	"""Find the loops from info on branches and
	nodes. The starting point is the first branch in br_map.
	The loops found need not be independent loops."""
	no_nodes=len(br_map)
	# First branch
	start_row=elem[0]
	start_col=elem[1]
	# Move right from that element
	# This is the first element
	# In a general sense, the direction is horiz
	loop_dir="horiz"
	# The termination condition is
	# that there should not be any element
	# in the nd_list. The nodes are deleted
	# as a completed loop contains them.
	# This is to ensure that all the nodes
	# are included in the loops found.
	# To ensure that parallel loops between
	# a few pair of nodes, do not cause
	# loops to be left out, additionally,
	# it is checked whether
	# Loops < Branches - Nodes + 1
	# while (nd_list or lp_count<lp_limit):
	while (lp_iter):
	# Will be executed if we are moving horizontally
	if (loop_dir == "horiz"):
	lp_iter=loop_horiz(br_map, nd_list, lp_list, lp_iter, elem, lp_map_list)
	# Change direction to vertical
	loop_dir="vert"
	# Will be executed if we are moving vertically
	if (loop_dir == "vert"):
	lp_iter=loop_vert(br_map, nd_list, lp_list, lp_iter, elem, lp_map_list)
	# Change direction to horizontal
	loop_dir="horiz"
	return

view raw loops.py hosted with ❤ by GitHub

After a long break

Been a while since I updated this blog. I have been working on the project off and on but I was busy doing all sorts of other things that made documentation difficult if not impossible. I have had to extend my work permit and that means getting together all my documents. I have had to send a paper to a conference in Hiroshima, Japan which meant I had to rush to get the final paper ready. And finally, had to plan my trip to Japan which will be in the month of May and get my visa ready.

There have been some questions about how the simulator should work for nonlinear circuits and I am going to start thinking out aloud on the next few blog posts.

Releasing version 0.2.4

Releasing the next version with those minor changes.

http://sourceforge.net/projects/pythonpowerelec/

Global variables for control code

So far, the only data that are being plotted in the output file are the values measured by the meters. But when writing control code, the user would need to plot control variables to debug the control code. Also, in the case when there are multiple control blocks and some may be cascaded - that is the output of one is the input of the other, this can't be done in the present form as the inputs are only from meters.

To solve these two shortcomings, I thought of creating another category of variables called "VariableStorage". From the name, the primary use of these are to tag certain variables in the control code as variables that can be plotted. Also, these are global variables for all the control code - this means one such global variable defined for one control code can be used anywhere. So the output of one can be used as the input of another.

But, the creation of such global variables is a major threat because some stupid statement somewhere in some control code could be resetting a variable to an arbitrary value. Not sure what protection I need to incorporate to prevent random usage. For now this is that I will do. (click on "view raw" below the code box to see the code in another window):

	if params_from_file[c2][0].lower()=="variablestorage":
	# If it is a variable storage, the dictionary key
	# will be the variable name.
	varstore_name=params_from_file[c2][1].split("=")[1]
	while varstore_name[0]==" ":
	varstore_name=varstore_name[1:]
	while varstore_name[-1]==" ":
	varstore_name=varstore_name[:-1]
	varstore_val=params_from_file[c2][2].split("=")[1]
	while varstore_val[0]==" ":
	varstore_val=varstore_val[1:]
	while varstore_val[-1]==" ":
	varstore_val=varstore_val[:-1]
	try:
	float(varstore_val)
	except:
	pass
	else:
	varstore_val=float(varstore_val)
	varstore_plot=params_from_file[c2][3].split("=")[1]
	while varstore_plot[0]==" ":
	varstore_plot=varstore_plot[1:]
	while varstore_plot[-1]==" ":
	varstore_plot=varstore_plot[:-1]
	if varstore_name not in control_file_variablestorage.keys():
	control_file_variablestorage[varstore_name]=[varstore_val, varstore_plot]
	else:
	print "Multiple definitions of variable storage %s" %varstore_name

view raw variable_store1.py hosted with ❤ by GitHub

Events and simulation time step - continued

A solution to the problem of running certain control algorithms at a faster rate than the simulation is to create a time event list. The closest time event will be the next simulation time instant. This will result in a variable time-step solver and am not sure how this will affect the future network analysis. The user simulation time step therefore becomes the maximum allowable time step. The code is as follows (click on "view raw" below the code box to see the code in a new window):

	# If there has been a system change or time_event
	# or if time is greater than the simulation time step.
	if ((t>=t_ode) or (time_event=="yes")):
	# Update the input values in u matrix
	for c1 in range(len(source_list)):
	component_objects[source_list[c1]].generate_val(source_list, system_loops_copy, sys_mat_e, sys_mat_a, sys_mat_b, sys_mat_u, t, t-t_ode_prev)
	# The ODE solver.
	# Will return the x(k+1) value called
	# as next_state_vec from x(k) value
	# called curr_state_vec
	mat_ode(sys_mat_e, sys_mat_a, sys_mat_b, [curr_state_vec, next_state_vec], sys_mat_u, t-t_ode_prev, ode_var, system_loops_copy, node_list)
	# Update the values of objects based on x(k+1)
	# So far only ammeter affected by this.
	for comps in component_objects.keys():
	component_objects[comps].update_val(system_loop_map, stiff_ratio, sys_mat_e, sys_mat_a, sys_mat_b, next_state_vec, sys_mat_u, branch_params, branch_events)
	# Recalculate the branch currents from the loop currents
	for c1 in range(len(branch_params)):
	branch_params[c1][-1][2]=0.0
	for c2 in range(len(system_loop_map)):
	if system_loop_map[c2][c1]=="forward" or system_loop_map[c2][c1]=="stiff_forward":
	branch_params[c1][-1][2]+=next_state_vec.data[c2][0]
	elif system_loop_map[c2][c1]=="reverse" or system_loop_map[c2][c1]=="stiff_reverse":
	branch_params[c1][-1][2]-=next_state_vec.data[c2][0]
	# x(k)=x(k+1) for next iteration.
	for c1 in range(system_size):
	curr_state_vec.data[c1][0]=next_state_vec.data[c1][0]
	# Save the previous time instant of ODE solver
	t_ode_prev=t
	# If the above code is running because of a
	# change in the system and not because of time
	# greater than t_ode, do not increment the t_ode.
	if time_event=="yes" and t<t_ode:
	time_event="no"
	else:
	t_ode=t+dt
	# Store the measured values.
	if (t>=t_store):
	f.write("%s " %str(t),)
	for c1 in range(len(meter_list)):
	f.write("%s " %component_objects[meter_list[c1]].op_value,)
	#for c1 in range(len(system_loops_copy)):
	# f.write("%s " %next_state_vec.data[c1][0],)
	f.write("\n")
	t_store=t_store+dt_store
	# This time vector will contain all the future time events
	# generated by the control functions and the main ODE solver
	time_vector=[t_ode]
	# Call the control codes only if controlled elements exist.
	if controlled_elements:
	# Call the control functions in the main control programs.
	# Use the eval function to call the functions as string arguments
	for c1 in range(len(control_files)):
	eval("%s(control_file_inputs, control_file_outputs, control_file_staticvars, control_file_timeevents, component_objects, c1, t, time_vector)" %control_functions[c1])
	# The soonest event will be the next time instant.
	time_vector.sort()
	t=time_vector[0]

view raw ode_update.py hosted with ❤ by GitHub

In the control function, one of the last statements updates the time event vector with the time events in the individual control code. So, __control.py will look like this (click on "view raw" below the code box to see the code in a new window):

	import math
	def control_func(interface_inputs, interface_outputs, interface_static, interface_time, circuit_components, pos, t_clock, sys_t_events):
	Vo=circuit_components['9Q'].op_value
	x_tri=interface_static[pos]['x_tri']
	t1=interface_time[pos]['t1']
	import math
	carr_freq=5000.0
	if (x_tri>1.0):
	x_tri=1.0
	if (x_tri<0.0):
	x_tri=0.0
	dt_sample=1.0e-6
	try:
	str(switch_control)
	except:
	switch_control=0.0
	else:
	pass
	if t_clock>=t1:
	x_tri+=(0.5*carr_freq)*dt_sample
	if (x_tri>1.0):
	x_tri=0.0
	if (x_tri<0.5):
	switch_control=1.0
	else:
	switch_control=0.0
	t1=t1+dt_sample
	circuit_components['5C'].control_values[0]=switch_control
	interface_outputs[pos]['5C'][1][2]=switch_control
	interface_static[pos]['x_tri']=x_tri
	interface_time[pos]['t1']=t1
	sys_t_events.append(t1)
	return

view raw control__.py hosted with ❤ by GitHub

The last statement appends the time event t1 to the sys_t_events list. In case they were multiple time events, all of them would be appended to the list.

Events and simulation time step

There is a serious flaw in the way I handle simulation time steps. So far there are only two time events - the simulation time step that is decided by the user and the control time step that is decided by the user in the control file. However, it is assumed that the control time step will be larger than the simulation time step which is completely wrong.

For example, the controller of the buck converter will update the duty cycle once every switching cycle. But the comparison between the duty cycle and the carriers sawtooth waveform has to take place much faster - actually at least 9 times faster than the carrier frequency to get the crossover points accurately.

For now a simple way to resolve this would be to include two time events in the control code. However, at a later stage, the pulse width modulator would have to made into another block where the time step would default to a value appropriate with the switching frequency.

In general, the control philosophy will change. Instead of having the controllers as functions, it might be better to have them as objects. This would enable inputs and outputs to be connected together and also control signals will need to be measured to debug errors.

Releasing version 0.2.3

And with the changes made, releasing the next version of the circuit simulator.

http://sourceforge.net/projects/pythonpowerelec/

Not sure if this will work in every case. So that will be the next step, to test this for all the different types of power electronic converters. Will keep updating this blog as and how new findings are made.

As before, feel free to email me at pythonpowerelectronics@gmail.com if you have any questions.

Solving ODEs

This was an error I am surprised I was even able to locate to the solver because I thought it was something related to the circuit solver. So far, the ODE is solved as an entire state vector at a time. This means, all states are calculated together and when one state needs the information of another state it used the previous value of the state. So even if a more updated value of a state is available, with this method the previous value is used.

This actually defeats the purpose of solving ODEs in a backward manner. That is we solve from the last state to the first state since the matrices have been made upper triangular.

So, the ODE was changed such that a particular state is calculated completely, its value is updated and this updated value is available for the next state to be calculated.

Here is the code (click on "view raw" below the code box to view the code in another window):

	def mat_ode(matrix_e, matrix_a, matrix_b, state_vector, input_vector, dt, ode_vectors, sys_loops, list_of_nodes):
	""" Solves the ODE using Runge Kutta 4th order method. """
	def runge_function4(x_plus_dxdt, ddt_mat, dbydt_order, ode_row):
	""" Defines the function dx/dt=f(x) for 4th order
	Runge Kutta solver. """
	#Calculate d/dt vector in reverse.
	#for c1 in range(matrix_e.rows-1, -1, -1):
	if matrix_e.data[ode_row][ode_row]:
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	ddt_mat[dbydt_order-1].data[ode_row][0]=matrix_b.data[ode_row][0]*input_vector
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	ddt_mat[dbydt_order-1].data[ode_row][0]=0.0 # Added on Dec. 29, 2012.
	for c2 in range(matrix_b.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]+=matrix_b.data[ode_row][c2]*input_vector.data[c2][0]
	if (dbydt_order==2):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*0.5*ddt_mat[dbydt_order-2].data[c3][0]
	if (dbydt_order==3):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*0.5*ddt_mat[dbydt_order-2].data[c3][0]
	if (dbydt_order==4):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-2].data[c3][0]
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*x_plus_dxdt.data[c3][0]
	ddt_mat[dbydt_order-1].data[ode_row][0]=ddt_mat[dbydt_order-1].data[ode_row][0]*dt
	for c3 in range(ode_row+1, matrix_e.columns):
	# for c3 in range(matrix_e.columns):
	if not ode_row==c3:
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_e.data[ode_row][c3]*ode_vectors[4].data[c3][0]
	ddt_mat[dbydt_order-1].data[ode_row][0]=ddt_mat[dbydt_order-1].data[ode_row][0]
	# ddt_mat[dbydt_order-1].data[ode_row][0]=ddt_mat[dbydt_order-1].data[ode_row][0]/matrix_e.data[ode_row][ode_row]
	else:
	if not (matrix_a.data[ode_row][ode_row]):
	#The variable has no dynamics and can't even be calculated statically!
	# May be due to a redundant loop.
	ddt_mat[dbydt_order-1].data[ode_row][0]=0.0
	state_vector[0].data[ode_row][0]=0.0
	state_vector[1].data[ode_row][0]=0.0
	else:
	ddt_mat[dbydt_order-1].data[ode_row][0]=0.0
	state_vector[0].data[ode_row][0]=0.0
	state_vector[1].data[ode_row][0]=0.0
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	state_vector[0].data[ode_row][0]+=matrix_b.data[ode_row][0]*input_vector
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	for c2 in range(matrix_b.columns):
	state_vector[0].data[ode_row][0]+=matrix_b.data[ode_row][c2]*input_vector.data[c2][0]
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	for c2 in range(matrix_a.columns):
	if not (ode_row==c2):
	state_vector[0].data[ode_row][0]-=matrix_a.data[ode_row][c2]*state_vector[1].data[c2][0]
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	state_vector[0].data[ode_row][0]=state_vector[0].data[ode_row][0]/matrix_a.data[ode_row][ode_row]
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	return
	def runge_function5(x_plus_dxdt, ddt_mat, dbydt_order, ode_row):
	""" Defines the function dx/dt=f(x) for higher order
	ODE solver """
	#Calculate d/dt vector in reverse.
	#for c1 in range(matrix_e.rows-1, -1, -1):
	if matrix_e.data[ode_row][ode_row]:
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	ddt_mat[dbydt_order-1].data[ode_row][0]=matrix_b.data[ode_row][0]*input_vector
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	ddt_mat[dbydt_order-1].data[ode_row][0]=0.0 # Added on Dec. 29, 2012.
	for c2 in range(matrix_b.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]+=matrix_b.data[ode_row][c2]*input_vector.data[c2][0]
	if (dbydt_order==2):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-2].data[c3][0]
	if (dbydt_order==3):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*0.5*ddt_mat[dbydt_order-3].data[c3][0]
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*0.5*ddt_mat[dbydt_order-2].data[c3][0]
	if (dbydt_order==4):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-4].data[c3][0]*(14.0/64.0)
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-3].data[c3][0]*(5.0/64.0)
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-2].data[c3][0]*(-3.0/64.0)
	if (dbydt_order==5):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-5].data[c3][0]*(-12.0/96.0)
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-4].data[c3][0]*(-12.0/96.0)
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-3].data[c3][0]*(8.0/96.0)
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-2].data[c3][0]*(64.0/96.0)
	if (dbydt_order==6):
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-5].data[c3][0]*(-9.0/64.0)
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-4].data[c3][0]*(5.0/64.0)
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-3].data[c3][0]*(16.0/64.0)
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*ddt_mat[dbydt_order-2].data[c3][0]*(36.0/64.0)
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_a.data[ode_row][c3]*x_plus_dxdt.data[c3][0]
	ddt_mat[dbydt_order-1].data[ode_row][0]=ddt_mat[dbydt_order-1].data[ode_row][0]*dt
	for c3 in range(matrix_e.columns):
	if not ode_row==c3:
	ddt_mat[dbydt_order-1].data[ode_row][0]-=matrix_e.data[ode_row][c3]*ode_vectors[4].data[c3][0]
	ddt_mat[dbydt_order-1].data[ode_row][0]=ddt_mat[dbydt_order-1].data[ode_row][0]
	else:
	if not (matrix_a.data[ode_row][ode_row]):
	#The variable has no dynamics and can't even be calculated statically!
	# May be due to a redundant loop.
	state_vector[0].data[ode_row][0]=0.0
	state_vector[1].data[ode_row][0]=0.0
	else:
	ddt_mat[dbydt_order-1].data[ode_row][0]=0.0
	state_vector[0].data[ode_row][0]=0.0
	state_vector[1].data[ode_row][0]=0.0
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	state_vector[0].data[ode_row][0]+=matrix_b.data[ode_row][0]*input_vector
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	for c2 in range(matrix_b.columns):
	state_vector[0].data[ode_row][0]+=matrix_b.data[ode_row][c2]*input_vector.data[c2][0]
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	for c2 in range(matrix_a.columns):
	if not (ode_row==c2):
	state_vector[0].data[ode_row][0]-=matrix_a.data[ode_row][c2]*state_vector[1].data[c2][0]
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	state_vector[0].data[ode_row][0]=state_vector[0].data[ode_row][0]/matrix_a.data[ode_row][ode_row]
	state_vector[1].data[ode_row][0]=state_vector[0].data[ode_row][0]
	return
	def trapezoidal_mthd(x_plus_dxdt, ddt_mat, dbydt_order):
	""" Defines the function dx/dt=f(x) for Trapezoidal
	method. """
	#Calculate d/dt vector in reverse.
	for c1 in range(matrix_e.rows-1, -1, -1):
	if matrix_e.data[c1][c1]:
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	ddt_mat[dbydt_order-1].data[c1][0]=matrix_b.data[c1][0]*input_vector
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	ddt_mat[dbydt_order-1].data[c1][0]=0.0 # Added on Dec. 29, 2012.
	for c2 in range(matrix_b.columns):
	ddt_mat[dbydt_order-1].data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0]
	for c3 in range(matrix_a.columns):
	ddt_mat[dbydt_order-1].data[c1][0]-=matrix_a.data[c1][c3]*x_plus_dxdt.data[c3][0]
	ddt_mat[dbydt_order-1].data[c1][0]=ddt_mat[dbydt_order-1].data[c1][0]*dt
	for c3 in range(matrix_e.columns):
	if not c1==c3:
	ddt_mat[dbydt_order-1].data[c1][0]-=matrix_e.data[c1][c3]*ode_vectors[2].data[c3][0]
	ddt_mat[dbydt_order-1].data[c1][0]=ddt_mat[dbydt_order-1].data[c1][0]
	else:
	if not (matrix_a.data[c1][c1]):
	#The variable has no dynamics and can't even be calculated statically!
	# May be due to a redundant loop.
	ddt_mat[dbydt_order-1].data[c1][0]=0.0
	state_vector[0].data[c1][0]=0.0
	state_vector[1].data[c1][0]=0.0
	else:
	ddt_mat[dbydt_order-1].data[c1][0]=0.0
	state_vector[0].data[c1][0]=0.0
	state_vector[1].data[c1][0]=0.0
	try:
	if input_vector.rows:
	pass
	except:
	if not (matrix_b.columns==1):
	print "Input signal has to be a real number and not a vector."
	else:
	state_vector[0].data[c1][0]+=matrix_b.data[c1][0]*input_vector
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	else:
	if not (matrix_b.columns==input_vector.rows):
	print "Dimension of input vector incorrect."
	else:
	for c2 in range(matrix_b.columns):
	state_vector[0].data[c1][0]+=matrix_b.data[c1][c2]*input_vector.data[c2][0]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	for c2 in range(matrix_a.columns):
	if not (c1==c2):
	state_vector[0].data[c1][0]-=matrix_a.data[c1][c2]*state_vector[0].data[c2][0]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	state_vector[0].data[c1][0]=state_vector[0].data[c1][0]/matrix_a.data[c1][c1]
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]
	return
	# Trapezoidal method
	## trapezoidal_mthd(state_vector[0], ode_vectors, 2) #k1=dt*(dx/dt)
	## ode_vectors[2].data[c1][0]=(ode_vectors[0].data[c1][0]+ode_vectors[1].data[c1][0])/2.0
	## state_vector[1].data[c1][0]=state_vector[0].data[c1][0]+ode_vectors[2].data[c1][0]
	## ode_vectors[0].data[c1][0]=ode_vectors[1].data[c1][0]
	# Runge-Kutta 4th order method
	for c1 in range(matrix_e.rows-1, -1, -1):
	runge_function4(state_vector[0], ode_vectors, 1, c1) #k1=dt*(dx/dt)
	runge_function4(state_vector[0], ode_vectors, 2, c1) #k2=dt*d/dt(x+k1/2)
	runge_function4(state_vector[0], ode_vectors, 3, c1) #k3=dt*d/dt(x+k2/2)
	runge_function4(state_vector[0], ode_vectors, 4, c1) #k4=dt*d/dt(x+k3)
	ode_vectors[4].data[c1][0]=(ode_vectors[0].data[c1][0]+2.0*ode_vectors[1].data[c1][0]+2.0*ode_vectors[2].data[c1][0]+ode_vectors[3].data[c1][0])/6.0
	state_vector[1].data[c1][0]=state_vector[0].data[c1][0]+ode_vectors[4].data[c1][0]
	# Runge-Kutta 5th order method
	# for c1 in range(state_vector[1].rows):
	# runge_function5(state_vector[0], ode_vectors, 1, c1) #k1=dt*(dx/dt)
	# runge_function5(state_vector[0], ode_vectors, 2, c1) #k2=dt*d/dt(x+k1)
	# runge_function5(state_vector[0], ode_vectors, 3, c1) #k3=dt*d/dt(x+(k1+k2)/2)
	# runge_function5(state_vector[0], ode_vectors, 4, c1) #k4=dt*d/dt(x+(14k1+5k2-3k3)/64)
	# runge_function5(state_vector[0], ode_vectors, 5, c1) #k4=dt*d/dt(x+(-12k1-12k2+8k3+64k4)/96)
	# runge_function5(state_vector[0], ode_vectors, 6, c1) #k4=dt*d/dt(x+(-9k2+5k3+16k4+36k5)/64)
	# ode_vectors[4].data[c1][0]=(ode_vectors[0].data[c1][0]+2.0*ode_vectors[1].data[c1][0]+2.0*ode_vectors[2].data[c1][0]+ode_vectors[3].data[c1][0])/6.0
	# state_vector[1].data[c1][0]=state_vector[0].data[c1][0]+ode_vectors[4].data[c1][0]
	## ode_vectors[6].data[c1][0]=7.0*ode_vectors[0].data[c1][0]
	## ode_vectors[6].data[c1][0]+=12.0*ode_vectors[2].data[c1][0]
	## ode_vectors[6].data[c1][0]+=7.0*ode_vectors[3].data[c1][0]
	## ode_vectors[6].data[c1][0]+=32.0*ode_vectors[4].data[c1][0]
	## ode_vectors[6].data[c1][0]+=32.0*ode_vectors[5].data[c1][0]
	## ode_vectors[6].data[c1][0]=ode_vectors[6].data[c1][0]/90.0
	## state_vector[1].data[c1][0]=state_vector[0].data[c1][0]+ode_vectors[6].data[c1][0]
	return

view raw ode.py hosted with ❤ by GitHub

Simplifying loop manipulation

Initially, all the code was written with respect to system loops. Also, when manipulating loops either for removing stiffness or for any other reason, the loops themselves were manipulated. And this means the entire loop information, including the branch elements. All this only adds to overhead. Particularly since two other elements have been defined - branch_params which contains information of the branch - the elements, the resistance, inductance, voltages, and current. The other element is the system_loop_map which contains information of which branch exists in each loop.

So why not get rid of system_loops in further manipulations? All that is needed to be done is to extend the definition of system_loop_map. Before system_loop_map had only "stiff", "yes" and "no" as entries. Now it will have:

forward - branch is in the direction of loop
reverse - branch is against loop direction
stiff_forward - stiff branch but same direction
stiff_reverse - stiff branch but opposite direction

With this the entire information of the loop is contained in a far more compact manner. And all the functions were rewritten. For example, the remove_stiffness function is as follows (click on "view raw" below the code box to see the code in another window):

	def remove_stiffness(matrix_e, matrix_a, matrix_b, dt, sys_loops, branch_info, stiff_info, sys_loop_map, src_list):
	""" Reduces the stiffness of the ODE by limiting the
	appearance of stiff branches in the loops."""
	# First to find the stiff loops by locating the first branch
	# that is stiff. After that make all the loops below that loop
	# have no element for that stiff branch.
	# To first arrange the loops such that the first branch
	# that is stiff occurs in the first loop and so on.
	map_loop_count=0
	map_branch_count=0
	# Start with the first branch of the first loop.
	while map_loop_count<len(sys_loop_map):
	if sys_loop_map[map_loop_count][map_branch_count]=="stiff_forward" or sys_loop_map[map_loop_count][map_branch_count]=="stiff_reverse":
	loop_found="yes"
	else:
	loop_found="no"
	# If a loop is not found and the branches are not exhausted.
	while loop_found=="no" and map_branch_count<len(sys_loop_map[0]):
	# Check if the branch is stiff in the loop.
	if sys_loop_map[map_loop_count][map_branch_count]=="stiff_forward" or sys_loop_map[map_loop_count][map_branch_count]=="stiff_reverse":
	loop_found="yes"
	# If not, look at all the remaining loops corresponding to that branch.
	# The idea is that if a branch is stiff, it will appear in at least
	# one loop and so need to look through all loops to make sure it is
	# not missed.
	c1=map_loop_count+1
	while loop_found=="no" and c1<len(sys_loop_map):
	if sys_loop_map[c1][map_branch_count]=="stiff_forward" or sys_loop_map[c1][map_branch_count]=="stiff_reverse":
	for c2 in range(len(sys_loop_map[0])):
	# If a subsequenct loop is found to have the branch as stiff,
	# interchange the rows.
	sys_loop_map[c1][c2], sys_loop_map[map_loop_count][c2] = sys_loop_map[map_loop_count][c2], sys_loop_map[c1][c2]
	loop_found="yes"
	c1+=1
	# If all loops are exhausted and no element is found as stiff
	# for the branch, it means the branch is not stiff
	# or has been accounted in a previous loop.
	# So move on to the next branch.
	if loop_found=="no":
	map_branch_count+=1
	# If a loop has been found, the loops need to be made as upper triangular
	# So all stiff branches in the loops subsequent to the loop will have to eliminated.
	if loop_found=="yes":
	for c1 in range(map_loop_count+1, len(sys_loop_map)):
	if sys_loop_map[c1][map_branch_count]=="stiff_forward" or sys_loop_map[c1][map_branch_count]=="stiff_reverse":
	# This is simply row operation of loops.
	if sys_loop_map[map_loop_count][map_branch_count]==sys_loop_map[c1][map_branch_count]:
	for c2 in range(len(sys_loop_map[0])):
	if sys_loop_map[map_loop_count][c2]=="forward" and sys_loop_map[c1][c2]=="forward":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="reverse" and sys_loop_map[c1][c2]=="reverse":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="reverse" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="forward"
	elif sys_loop_map[map_loop_count][c2]=="forward" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="reverse"
	elif sys_loop_map[map_loop_count][c2]=="stiff_forward" and sys_loop_map[c1][c2]=="stiff_forward":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="stiff_reverse" and sys_loop_map[c1][c2]=="stiff_reverse":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="stiff_reverse" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="stiff_forward"
	elif sys_loop_map[map_loop_count][c2]=="stiff_forward" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="stiff_reverse"
	else:
	for c2 in range(len(sys_loop_map[0])):
	if sys_loop_map[map_loop_count][c2]=="forward" and sys_loop_map[c1][c2]=="reverse":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="reverse" and sys_loop_map[c1][c2]=="forward":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="reverse" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="reverse"
	elif sys_loop_map[map_loop_count][c2]=="forward" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="forward"
	elif sys_loop_map[map_loop_count][c2]=="stiff_forward" and sys_loop_map[c1][c2]=="stiff_reverse":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="stiff_reverse" and sys_loop_map[c1][c2]=="stiff_forward":
	sys_loop_map[c1][c2]="no"
	elif sys_loop_map[map_loop_count][c2]=="stiff_reverse" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="stiff_reverse"
	elif sys_loop_map[map_loop_count][c2]=="stiff_forward" and sys_loop_map[c1][c2]=="no":
	sys_loop_map[c1][c2]="stiff_forward"
	# Increment the loop count and make the branch count equal that
	# So the default is a diagonal stiff element.
	map_loop_count+=1
	map_branch_count=map_loop_count
	# This is to make sure, that stiff loops
	# are connected to an input.
	for c1 in range(len(sys_loop_map)):
	is_loop_stiff="no"
	has_input="no"
	for c3 in range(len(branch_info)):
	#if branch_info[c3][:-1]==sys_loops[c1][c1][c2][:-1]:
	if sys_loop_map[c1][c3]=="stiff_forward" or sys_loop_map[c1][c3]=="stiff_reverse":
	is_loop_stiff="yes"
	# Check if any branch has a non-zero B matrix entry.
	for c4 in range(len(branch_info[c3][-1][1])):
	if branch_info[c3][-1][1][c4]:
	has_input="yes"
	# If loop is stiff and has no input
	if is_loop_stiff=="yes" and has_input=="no":
	c2=0
	flag_input="no"
	# Check all other loops whether they are stiff
	# and whether they have inputs.
	while flag_input=="no" and c2<len(sys_loop_map):
	if not c1==c2:
	is_loop_stiff="no"
	flag_input="no"
	for c4 in range(len(branch_info)):
	if sys_loop_map[c2][c4]=="stiff_forward" or sys_loop_map[c2][c4]=="stiff_reverse":
	is_loop_stiff="yes"
	if is_loop_stiff=="no":
	for c5 in range(len(branch_info[c4][-1][1])):
	if branch_info[c4][-1][1][c5]:
	flag_input="yes"
	c2=c2+1
	# Perform a row operation with a loop that is non-stiff
	# and has an input.
	if is_loop_stiff=="no" and flag_input=="yes":
	c2=c2-1
	#loop_manipulate(sys_loops, c1, c2, "add")
	# Update the sys_loop_map info
	for c3 in range(len(branch_info)):
	if sys_loop_map[c1][c3]=="forward" and sys_loop_map[c2][c3]=="reverse":
	sys_loop_map[c1][c3]="no"
	elif sys_loop_map[c1][c3]=="reverse" and sys_loop_map[c2][c3]=="forward":
	sys_loop_map[c1][c3]="no"
	elif sys_loop_map[c1][c3]=="forward" and sys_loop_map[c2][c3]=="no":
	sys_loop_map[c1][c3]="forward"
	elif sys_loop_map[c1][c3]=="reverse" and sys_loop_map[c2][c3]=="no":
	sys_loop_map[c1][c3]="reverse"
	elif sys_loop_map[c1][c3]=="stiff_forward" and sys_loop_map[c2][c3]=="stiff_reverse":
	sys_loop_map[c1][c3]="no"
	elif sys_loop_map[c1][c3]=="stiff_reverse" and sys_loop_map[c2][c3]=="stiff_forward":
	sys_loop_map[c1][c3]="no"
	elif sys_loop_map[c1][c3]=="stiff_forward" and sys_loop_map[c2][c3]=="no":
	sys_loop_map[c1][c3]="stiff_forward"
	elif sys_loop_map[c1][c3]=="stiff_reverse" and sys_loop_map[c2][c3]=="no":
	sys_loop_map[c1][c3]="stiff_reverse"
	# Get rid of repeating loops by setting all branches "no"
	for c1 in range(len(sys_loop_map)-1):
	for c2 in range(c1+1, len(sys_loop_map)):
	loop_repeats="yes"
	for c3 in range(len(sys_loop_map[c1])):
	if not sys_loop_map[c1][c3]==sys_loop_map[c2][c3]:
	loop_repeats="no"
	if loop_repeats=="yes":
	for c3 in range(len(sys_loop_map[c2])):
	sys_loop_map[c2][c3]="no"
	# Re-initialize the matrices A, B, and E
	matrix_a.zeros(len(sys_loops),len(sys_loops))
	matrix_e.zeros(len(sys_loops),len(sys_loops))
	matrix_b.zeros(len(sys_loops),matrix_b.columns)
	# Recalculate the matrices for the new system loops.
	recalculate_sys_matrices(sys_loop_map, branch_info, matrix_a, matrix_e, matrix_b)
	for c1 in range(matrix_a.rows):
	for c2 in range(matrix_a.columns):
	if abs(matrix_a.data[c1][c2])<1.0e-10:
	matrix_a.data[c1][c2]=0.0
	for c1 in range(matrix_e.rows):
	for c2 in range(matrix_e.columns):
	if abs(matrix_e.data[c1][c2])<1.0e-10:
	matrix_e.data[c1][c2]=0.0
	for c1 in range(matrix_b.rows):
	for c2 in range(matrix_b.columns):
	if abs(matrix_b.data[c1][c2])<1.0e-10:
	matrix_b.data[c1][c2]=0.0
	return

view raw remove_stiffness.py hosted with ❤ by GitHub

And the resultant code is much cleaner, the loop manipulations much easier. This will also change the code in the circuit elements (click on "view raw" below the code box to see the code in another window):

	def update_val(self, sys_loop_map, lbyr_ratio, mat_e, mat_a, mat_b, state_vec, mat_u, sys_branches, sys_events):
	""" This function calculates the actual current in the
	diode branch. With this, the branch voltage is found
	with respect to the existing diode resistance. The diode
	voltage is then used to decide the turn on condition. """
	for c1 in range(len(sys_branches)):
	if csv_tuple(self.pos) in sys_branches[c1]:
	branch_pos=c1
	self.current=0.0
	for c1 in range(len(sys_loop_map)):
	# If diode cathode is after the voltmeter
	# position, it means current is positive.
	if sys_branches[branch_pos].index(self.polrty)>sys_branches[branch_pos].index(csv_tuple(self.pos)):
	# Then check is the loop is aiding or opposing
	# the main loop.
	if sys_loop_map[c1][branch_pos]=="forward" or sys_loop_map[c1][branch_pos]=="stiff_forward":
	self.current+=state_vec.data[c1][0]
	elif sys_loop_map[c1][branch_pos]=="reverse" or sys_loop_map[c1][branch_pos]=="stiff_reverse":
	self.current-=state_vec.data[c1][0]
	else:
	if sys_loop_map[c1][branch_pos]=="forward" or sys_loop_map[c1][branch_pos]=="stiff_forward":
	self.current-=state_vec.data[c1][0]
	elif sys_loop_map[c1][branch_pos]=="reverse" or sys_loop_map[c1][branch_pos]=="stiff_reverse":
	self.current+=state_vec.data[c1][0]
	self.voltage=self.current*self.resistor
	# Diode will turn on when it is forward biased
	# and it was previously off.
	if self.status=="off" and self.voltage>1.0:
	sys_events[branch_pos]="yes"
	self.status="on"
	# Diode will turn off only when current becomes
	# negative.
	if self.status=="on" and self.current<-1.0e-5:
	sys_events[branch_pos]="yes"
	self.status="off"
	if self.status=="off":
	self.resistor=self.resistor_off
	else:
	self.resistor=self.resistor_on
	return

view raw update_value.py hosted with ❤ by GitHub

Diode freewheeling - continued

There were a couple of fatal flaws in the previous code. When I expand the circuit to a buck converter with two diodes connected in parallel as shown below, it doesn't work.

So essentially there are two freewheeling paths. In the previous algorithms, I had considered as current sources:
1. Inductors
2. Voltage sources without any resistance and inductance
3. Branches without any resistance and inductance and also no voltage - short branches or wires.

The 3rd point was wrong. A wire will carry any current and to fix the current through it will be wrong. For now, I think my assumption of a voltage source with zero impedance as a current source is acceptable because the sources should continue supplying what they were during an event to check which branches can change their status.

So, to include the third point, I came up with the concept of shortnode. A supernode will be a node having a zero impedance branch with a voltage. On the other hand, a short node is a node that has a zero impedance branch with no voltage - or a node with a "wire" connected to it. The significance of a short node will be that a "wire" will have two short nodes at the same voltage and in this manner, nodes connected through wires will end up as a group of short nodes. Similar to supernodes, a circuit could have many groups of short nodes. The only difference between a super node and a short node is that when performing KCL, is a node is encountered that is a short node, only the voltage of the first node in that group of short nodes will be included.

To expand on this. If there is a group of four short nodes. All these short nodes are connected by "wires" and so they are at the same voltage. So when performing KCL only one node voltage needs to be included. The reason is that there can be a large number of short nodes and therefore, to write separate equations V1-V2=0 will increase the number of equations well beyond the number of nodes, messing up the KCL.

Once this is done, the KCL is performed. All the node voltages are calculated and all the branch currents too. Except for the short branches or "wires". There is no way to calculate the current through them by the above method as they do not have a resistance. So, for this another KCL is written. The code for this is completely new and so is below (click on "view raw" below the code box to view the code in another window):

	# In the above nodal analysis, the inductor currents and the currents
	# through branches containing a voltage source are considered as
	# current sources. Finally, the branch currents calculated are the
	# currents with resistances or stiff branches where inductance is negligible.
	# But this leaves out the short branches or "wires". There is no means
	# to calculate the current through them as they do not have a resistance.
	# So the next part does that merely by using KCL, currents at a node
	# must be zero.
	# This is a list of the short branches or "wires"
	short_branch_list=[]
	for c1 in range(len(branch_info)):
	if branch_info[c1][-1][0][0]==0.0 and branch_info[c1][-1][0][1]==0.0:
	if ((1.0 not in branch_info[c1][-1][1]) and (-1.0 not in branch_info[c1][-1][1])):
	short_branch_list.append(c1)
	# This is a matrix showing the presence of short branches
	connectivity_short=[]
	# This is a vector with the remaining currents calculated
	# above by nodal analysis at a node.
	connectivity_src=[]
	for c1 in range(len(list_of_nodes)):
	row_vector=[]
	connectivity_src.append(0.0)
	for c1 in range(len(branch_info)):
	row_vector.append(0.0)
	connectivity_short.append(row_vector)
	for c1 in range(len(list_of_nodes)):
	for c2 in range(len(branch_info)):
	# Check if this branch is a short branch
	is_short_node="no"
	if branch_info[c2][-1][0][0]==0.0 and branch_info[c2][-1][0][1]==0.0:
	if ((1.0 not in branch_info[c2][-1][1]) and (-1.0 not in branch_info[c2][-1][1])):
	is_short_node="yes"
	# If it is, then this will apear in the connectivity matrix
	# If the branch originates at the node, it is taken with
	# a minus sign because it is on the other side of the KCL equation.
	if list_of_nodes[c1]==branch_info[c2][0]:
	connectivity_short[c1][c2]=-1.0
	if list_of_nodes[c1]==branch_info[c2][-2]:
	connectivity_short[c1][c2]=1.0
	if is_short_node=="no":
	# If it is not a short branch, add the currents calculated in KCL
	# to the source vector.
	if list_of_nodes[c1]==branch_info[c2][0]:
	connectivity_src[c1]+=br_currents[c2]
	if list_of_nodes[c1]==branch_info[c2][-2]:
	connectivity_src[c1]-=br_currents[c2]
	# Iterate through the short branch list
	c2=0
	# Iterate through each node - each row of connectivity_short.
	for c1 in range(len(list_of_nodes)):
	# Check if the row and branch number of the short branch is 1.0
	# This should be because if it is a short branch, there should
	# be a 1 or a -1 somewhere in that column for short branch number.
	if not connectivity_short[c1][short_branch_list[c2]]:
	# If not, go through all the remaining nodes to check
	# if any subsequent row has a nonzero.
	for c3 in range(c1+1, len(list_of_nodes)):
	if connectivity_short[c3][short_branch_list[c2]]:
	# If so, interchange the two rows.
	for c4 in range(len(connectivity_short[c1])):
	connectivity_short[c1][c4], connectivity_short[c3][c4] = connectivity_short[c3][c4], connectivity_short[c1][c4]
	connectivity_src[c1], connectivity_src[c3] = connectivity_src[c3], connectivity_src[c1]
	# If all the short branches are not exhausted, increment c2.
	if c2<len(short_branch_list)-1:
	c2+=1
	# The above row interchanges will gurantee the first rows
	# equal to the number of short branches will have nonzero
	# entries corresponding to the short branches.
	# Now, a row operation will be performed to make all other rows zero.
	# This will make connectivity_short upper triangular.
	# Once again iterate through the short branch list
	for c1 in range(len(short_branch_list)):
	# Go through all the remaining rows with the aim
	# of making all entries of the short branch zero.
	for c2 in range(c1, len(list_of_nodes)-1):
	if connectivity_short[c2][short_branch_list[c1]]:
	for c3 in range(c1+1, len(list_of_nodes)):
	if connectivity_short[c3][short_branch_list[c1]]:
	short_factor=connectivity_short[c3][short_branch_list[c1]]
	for c4 in short_branch_list:
	connectivity_short[c3][c4]-=connectivity_short[c2][c4]*short_factor/connectivity_short[c2][short_branch_list[c1]]
	connectivity_src[c3]-=connectivity_src[c2]*short_factor/connectivity_short[c2][short_branch_list[c1]]
	# Initialize a list of short branch currents
	short_branch_currents=[]
	for c1 in short_branch_list:
	short_branch_currents.append(0.0)
	# Calculate the short branch currents in reverse as
	# an the matrices were made upper triangular by previous row
	# operation.
	for c1 in range(len(short_branch_list)-1, -1, -1):
	short_branch_currents[c1]=connectivity_src[c1]
	for c2 in range(len(short_branch_list)-1, c1, -1):
	short_branch_currents[c1]-=connectivity_short[c1][short_branch_list[c2]]*short_branch_currents[c2]
	short_branch_currents[c1]=short_branch_currents[c1]/connectivity_short[c1][short_branch_list[c1]]
	# The remaining elements in branch currents are filled.
	# These were zero before.
	for c1 in range(len(short_branch_list)):
	br_currents[short_branch_list[c1]]=short_branch_currents[c1]

view raw short_branch.py hosted with ❤ by GitHub

After this change, the code now works for the circuit shown above.