Releasing version 0.2.2

Finally, releasing the next version 0.2.2. The example of the buck converter is there in testckt16.csv with the control code in control.py. And the description of the circuit parameters in testckt16_params.csv and the controls in control_desc.csv.

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

Well, well, well. It is the 31st of December 2013. I started this project a little more than a year ago and frankly I never thought I would get this far. So far everything has been coded to be automatic and without hardcoding any particular circuits which I feared I might have to. If all goes well, I will add other equipment like transformers, motors. But first I need to simulate power electronic converters of every type which just the library that I have. I suppose a lot of errors will creep up but I think the main concept is there and should work out.


And finally a big thank you to my blog readers. I see a lot of hits to the blog and also received some very nice comments. Thank you so much for reading my blog and I hope that soon I'll have a real simulator that every one of you can use on any OS. This has been a great year 2013 and I hope 2014 is just as good if not better.

Happy New Year to all of you and I wish you all the success for an awesome 2014!

I'll take a break for a few days and post again when I recover from my hangover :)

Diode freewheeling - nodal analysis

Where do I begin? Been coding continuously for the past three days and now have to find where I am.

As posted before, the only way to solve the problem of getting the diode to freewheel is by a nodal analysis. Essentially, the current through an inductor should not change instantaneously. If it does, it is bad design and I'll have to figure that out later. But in standard converter circuits, the inductor is expected to freewheel through a diode. As in the example of the buck converter.

So what I have coded is a function to perform nodal analysis take inductors as current sources. So I perform a dc analysis on a snapshot of the circuit at an instant and find out where the currents would flow if the inductor current were to remain constant.

Take the snapshot when the switch turns off with the inductor La having a non-zero current. So this current will flow through the switch and the diode almost equally since these are two high resistance paths.

With these currents calculated through nodal analysis, I try to figure out if any of the devices in the circuit has the capability to change its state. Since only nonlinear devices can do so, the remaining devices will have empty functions. The diode and switch will have the function called determine_state. Here is the code for the switch (click on "View Raw" at the bottom of the code box to see code in another 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

And here is the code for the nodal analysis (click on "View Raw" at the bottom of the code box to see code in another window)

	def current_continuity(list_of_nodes, branch_info, branch_stiff, nd_voltage, br_currents, mho_matrix, src_vector, sys_inputs):
	"""Following an event, to find the branch currents through nodal analysis.
	The objective is to determine if certain devices must change their status
	to maintain continuity of inductor currents."""
	# Node voltages. There will be one reference node
	for c1 in range(len(nd_voltage)):
	nd_voltage[c1]=0.0
	# The admittance matrix. As it is more or less
	# a dc analysis of a snapshot of the circuit
	# only resistances are considered.
	for c1 in range(len(mho_matrix)):
	for c2 in range(len(mho_matrix)):
	mho_matrix[c1][c2]=0.0
	# The source vector. This will contain
	# currents sources - inductor currents,
	# currents through zero impedance branches
	# and finally also voltages across zero
	# impedance branches
	for c1 in range(len(src_vector)):
	src_vector[c1]=0.0
	# The concept is of nodal analysis.
	# If a branch has only a resistance the currents
	# are expressed as a difference of node voltages divided
	# by the resistance.
	# If a branch has an inductance and it is non stiff,
	# it is treated as a current source. This is the main
	# feature as the inductor current must not change
	# instantaneously.
	# If it has zero impedance, it is treated as a current
	# source. This was done later as the currents through
	# these branches can't be accounted for by nodal analysis
	# and these currents aren't suppose to change anyway.
	# It is only the nonlinear devices that are essentially
	# variable resistances that are supposed to have changing
	# currents.
	# Start iterating through the nodes
	for c1 in range(len(list_of_nodes)):
	# Look for the node as the starting or ending node in branches
	for c2 in range(len(branch_info)):
	if ((list_of_nodes[c1]==branch_info[c2][0]) or (list_of_nodes[c1]==branch_info[c2][-2])):
	# Mark the position of the other node
	if (list_of_nodes[c1]==branch_info[c2][0]):
	end_node_pos=list_of_nodes.index(branch_info[c2][-2])
	# The default direction of current is taken to be
	# away from the starting node - start to end.
	branch_src_dir=1.0
	else:
	end_node_pos=list_of_nodes.index(branch_info[c2][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[c1][c1]+=1.0/branch_info[c2][-1][0][0]
	mho_matrix[c1][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[c1]-=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[c1]-=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)):
	src_vector[c1]-=branch_src_dir*br_currents[c2]
	# There are some rows that are exact negatives of the other for the reason, that
	# they are essentially reference nodes. So these need not be included
	for c1 in range(len(mho_matrix)-1):
	for c2 in range(c1+1, len(mho_matrix)):
	rows_the_same="yes"
	for c3 in range(len(mho_matrix[c1])):
	if (mho_matrix[c1][c3]+mho_matrix[c2][c3]):
	rows_the_same="no"
	if rows_the_same=="yes":
	src_vector[c2]=0.0
	for c3 in range(len(mho_matrix[c2])):
	mho_matrix[c2][c3]=0.0
	# Finding the supernodes in the system
	# Supernodes are essentially those nodes that
	# have a zero impedance branch connected to them
	# These branches may or may not have a voltage source.
	# When there is a supernode, KCL continues to the other node
	# of the zero impedance branch. So essentially the two end
	# nodes of a zero impedance branch have their KCL equations
	# added and their node voltages are expressed by an equation
	# depending on whether they have a voltage source.
	# The complete list of supernodes
	supernode_list=[]
	for c1 in range(len(list_of_nodes)):
	# Each small segment of supernodes
	current_supernode=[]
	# Check if the node has been found as a supernode before
	node_found="no"
	for c2 in range(len(supernode_list)):
	if c1 in supernode_list[c2]:
	node_found="yes"
	# If not add it as a supernode
	# This is the first node in the current
	# running list. If no other nodes are found,
	# it will be deleted at the next iteration
	# as the node is not a supernode in that case.
	if node_found=="no":
	current_supernode.append(c1)
	# Iterate through the remaining nodes
	for c2 in range(len(list_of_nodes)):
	# Look for the nodes in both the
	# current list and the complete list
	node_found="no"
	if c2 in current_supernode:
	node_found=="yes"
	if node_found=="no":
	for c3 in range(len(supernode_list)):
	if c2 in supernode_list[c3]:
	node_found="yes"
	# If the node has not been found, check if it is connected
	# to any of the existing nodes in the current list
	# by branches with zero impedance.
	if node_found=="no":
	for c3 in range(len(branch_info)):
	# Check if a branch has zero impedance.
	if ((branch_info[c3][-1][0][0]==0.0) and (branch_info[c3][-1][0][1]==0.0)):
	# If the current node is the first node in the current
	# branch being examined, check if any of the nodes
	# in the current supernode list are the ending nodes.
	# If so append the current node to the current super node list
	if (list_of_nodes[c2]==branch_info[c3][0]):
	for c4 in current_supernode:
	if list_of_nodes[c4]==branch_info[c3][-2]:
	current_supernode.append(c2)
	# If the current node is the ending node in the current
	# branch being examined, check if any of the nodes
	# in the current supernode list are the starting nodes.
	# If so append the current node to the current super node list
	if (list_of_nodes[c2]==branch_info[c3][-2]):
	for c4 in current_supernode:
	if list_of_nodes[c4]==branch_info[c3][0]:
	current_supernode.append(c2)
	# There is another possibility
	# A node could be a supernode, but is connected only
	# to another node which appears later in the node list
	# So, the above calculation will not work as it looks
	# for a node connected to existing supernodes by zero
	# impedance branches.
	# So essentially it is a rerun to pickup any nodes
	# that were missed the first time.
	# Iterate through existing list of current supernodes.
	for c2 in current_supernode:
	# Then look through the entire list of nodes.
	for c3 in range(len(list_of_nodes)):
	# Check it is has been found as a supernode
	# in the current list
	node_found="no"
	if c3 in current_supernode:
	node_found="yes"
	# Check it has been found elsewhere.
	if node_found=="no":
	for c4 in range(len(supernode_list)):
	if c3 in supernode_list[c4]:
	node_found="yes"
	# Check if the current node is connected to any of the nodes
	# in the current supernode list by a zero impedance branch.
	if node_found=="no":
	for c5 in range(len(branch_info)):
	if ((branch_info[c5][-1][0][0]==0.0) and (branch_info[c5][-1][0][1]==0.0)):
	if (list_of_nodes[c3]==branch_info[c5][0]):
	for c6 in current_supernode:
	if (list_of_nodes[c6]==branch_info[c5][-2]):
	current_supernode.append(c3)
	if (list_of_nodes[c3]==branch_info[c5][-2]):
	for c6 in current_supernode:
	if (list_of_nodes[c6]==branch_info[c5][0]):
	current_supernode.append(c3)
	# Sort the list of current supernodes
	# The reason is so that the KCL equations can be added
	# and the first node will contain the sum.
	current_supernode.sort()
	# # If the current super node list has more than
	# # one node, there is a zero impedance branch
	# # and so it is super node list.
	if len(current_supernode)>1:
	supernode_list.append(current_supernode)
	# For every supernode list,
	# add the KCL equations at subsequnt supernodes to the
	# first node in the supernode list. Then set the other
	# equations to zero. Also, with the source vector.
	for c1 in range(len(supernode_list)):
	row1=supernode_list[c1][0]
	for c2 in range(1, len(supernode_list[c1])):
	row2=supernode_list[c1][c2]
	src_vector[row1]+=src_vector[row2]
	src_vector[row2]=0.0
	for c3 in range(len(mho_matrix[0])):
	mho_matrix[row1][c3]+=mho_matrix[row2][c3]
	mho_matrix[row2][c3]=0.0
	# Make a list of the zero rows.
	zero_rows=[]
	for c1 in range(len(mho_matrix)):
	is_row_zero="yes"
	for c2 in range(len(mho_matrix[c1])):
	if mho_matrix[c1][c2]:
	is_row_zero="no"
	if is_row_zero=="yes":
	zero_rows.append(c1)
	# Move the zero rows to the end of the matrix by row interchanges.
	for c1 in zero_rows:
	for c2 in range(len(mho_matrix)-1, c1, -1):
	if c2 not in zero_rows:
	src_vector[c1], src_vector[c2] = src_vector[c2], src_vector[c1]
	for c3 in range(len(mho_matrix[c1])):
	mho_matrix[c1][c3], mho_matrix[c2][c3] = mho_matrix[c2][c3], mho_matrix[c1][c3]
	# Now to add equations that describe the voltage difference
	# between two nodes connected by a zero impedance branch
	# The starting row if the size of the addmitance matrix
	# minus the zero rows.
	supernode_row=len(mho_matrix)-len(zero_rows)
	# Iterate through the supernode lists
	for c1 in range(len(supernode_list)):
	# In every supernode list, take a supernode
	# and express its voltage with respect to every
	# other supernode in that list
	for c2 in range(len(supernode_list[c1])-1):
	for c3 in range(c2+1, len(supernode_list[c1])):
	# Find out which supernode is the start and which is
	# the end node of the branch.
	for c4 in range(len(branch_info)):
	if ((branch_info[c4][0]==list_of_nodes[supernode_list[c1][c2]]) and (branch_info[c4][-2]==list_of_nodes[supernode_list[c1][c3]])):
	# Confirm whether the branch is a zero impedance branch
	if (branch_info[c4][-1][0][0]==0.0 and branch_info[c4][-1][0][1]==0.0):
	# This is mere V1-V2=Vbranch
	# The plus and minus signs are for different cases of
	# starting node and ending node.
	mho_matrix[supernode_row][supernode_list[c1][c2]]=1.0
	mho_matrix[supernode_row][supernode_list[c1][c3]]=-1.0
	for c5 in range(len(branch_info[c4][-1][1])):
	src_vector[supernode_row]+=branch_info[c4][-1][1][c5]*sys_inputs.data[c5][0]
	# Increment the supernode row pointer
	supernode_row+=1
	if ((branch_info[c4][-2]==list_of_nodes[supernode_list[c1][c2]]) and (branch_info[c4][0]==list_of_nodes[supernode_list[c1][c3]])):
	if (branch_info[c4][-1][0][0]==0.0 and branch_info[c4][-1][0][1]==0.0):
	mho_matrix[supernode_row][supernode_list[c1][c2]]=-1.0
	mho_matrix[supernode_row][supernode_list[c1][c3]]=1.0
	for c5 in range(len(branch_info[c4][-1][1])):
	src_vector[supernode_row]-=branch_info[c4][-1][1][c5]*sys_inputs.data[c5][0]
	supernode_row+=1
	# Now to solve the equation AX=B
	# Look for diagonal elements - check if they are zero
	# If so, look in that same column in subsequent rows
	# if there is a non zero element and exchange them
	# Later, make the matrix upper triangular.
	# Using row manipulations, make all the elements
	# below a diagonal row zero.
	for c1 in range(len(mho_matrix)):
	if not mho_matrix[c1][c1]:
	for c2 in range(c1+1, len(mho_matrix)):
	if mho_matrix[c2][c1]:
	src_vector[c1], src_vector[c2] = src_vector[c2], src_vector[c1]
	for c3 in range(len(mho_matrix[c1])):
	mho_matrix[c1][c3], mho_matrix[c2][c3] = mho_matrix[c2][c3], mho_matrix[c1][c3]
	if mho_matrix[c1][c1]:
	for c2 in range(c1+1, len(mho_matrix)):
	if mho_matrix[c2][c1]:
	src_vector[c2]-=src_vector[c1]*mho_matrix[c2][c1]/mho_matrix[c1][c1]
	for c3 in range(len(mho_matrix[c1])):
	mho_matrix[c2][c3]-=mho_matrix[c2][c1]*mho_matrix[c1][c3]/mho_matrix[c1][c1]
	# The last row of the manipulated admittance matrix will
	# be zero as it is the reference node. So taking this
	# voltage to be zero, calculate all the other node voltages
	# from second last to first.
	for c1 in range(len(mho_matrix)-2, -1, -1):
	if mho_matrix[c1][c1]:
	nd_voltage[c1]=src_vector[c1]
	for c2 in range(c1+1, len(mho_matrix[c1])):
	nd_voltage[c1]-=mho_matrix[c1][c2]*nd_voltage[c2]
	nd_voltage[c1]=nd_voltage[c1]/mho_matrix[c1][c1]
	# Calculate the branch currents as I=(Vnode1-Vnode2-Vsource)/Rbranch
	for c1 in range(len(branch_info)):
	if branch_stiff=="yes" or branch_info[c1][-1][0][1]==0.0:
	start_node=list_of_nodes.index(branch_info[c1][0])
	end_node=list_of_nodes.index(branch_info[c1][-2])
	if branch_info[c1][-1][0][0]:
	br_currents[c1]=(nd_voltage[start_node]-nd_voltage[end_node])/branch_info[c1][-1][0][0]
	for c2 in range(len(branch_info[c1][-1][1])):
	br_currents[c1]+=branch_info[c1][-1][1][c2]*sys_inputs.data[c2][0]/branch_info[c1][-1][0][0]
	return

view raw nodal_analysis.py hosted with ❤ by GitHub

Yes, it is a monster. Anyway, the basic concept is:

1. If a branch is not stiff and has an inductor, it becomes a current source.
2. If it has no impedance, the nodes on either end of the branch become super nodes.
3. If it only a resistance or is a stiff branch, it is expressed in KCL
(Vnode1-Vnode2)/Rbranch=Ibranch

4. Also if a branch has zero impedance, the current through it is treated as a current source.
5. Because of this, the only thing to be calculated are the resistive branches which is what needs to be done because that is where the non-linear devices will be.

6. The super node was a little tough to handle. When two nodes are connected by a branch with zero impedance, the current through that branch is tough to define. Therefore, the KCLs at the two nodes can be added up to form one KCL equation. And the second equation will be the voltage of the branch connecting these nodes.

7. In general, several nodes can be connected together to form a super node. Even worse, several groups of super nodes can exist.

8. Once the KCL equations are established, they need to reduced to upper triangular form and solved backwards. One node will become the reference node and will have zero voltage.

9. The branch currents will then be calculated from these node voltages. and there comes the next part of the story.


Because the nodal analysis could have changed the nature of the circuit by changing the state of nonlinear devices, it may so happen that the loop currents calculated before are not correct. So we need to recalculate the loop currents from the branch currents that the nodal analysis has given. The code for that is in compute_loop_currents (click on "View Raw" at the bottom of the code box to see code in another 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":
	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]=="yes":
	# 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]=="yes":
	does_branch_occur="yes"
	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 loop that exists in a loop
	if sys_loop_map[c1][c2]=="yes":
	# Look at the remaining nonstiff loops.
	for c3 in nonstiff_loops:
	if not c1==c3:
	# 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]=="yes":
	# Basically take a running counter with loop c1
	# starting from the last branch and keep checking
	# whether a branch exists. This is because as
	# manipulations take place, branches disappear
	# so an exception mght be thrown.
	c4=len(sys_loops[c1][c1])-1
	while c4>=0:
	try:
	sys_loops[c1][c1][c4]
	except:
	pass
	else:
	c5=len(sys_loops[c3][c3])-1
	while c5>=0:
	try:
	sys_loops[c3][c3][c5]
	except:
	pass
	else:
	if sys_loops[c1][c1][c4][:-1]==sys_loops[c3][c3][c5][:-1]:
	if branch_info[c2][:-1]==sys_loops[c1][c1][c4][:-1]:
	# If the branches are in the same sense, subtract them
	# or else add them.
	if sys_loops[c1][c1][c4][-1]==sys_loops[c3][c3][c5][-1]:
	loop_manipulate(sys_loops, c3, c1, "diff")
	else:
	loop_manipulate(sys_loops, c3, c1, "add")
	c5=c5-1
	c4=c4-1
	# 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)
	# Update the sys_loop_map info
	for c4 in range(len(branch_info)):
	if sys_loop_map[c1][c4]=="yes" and sys_loop_map[c3][c4]=="yes":
	sys_loop_map[c3][c4]="no"
	elif sys_loop_map[c1][c4]=="yes" and sys_loop_map[c3][c4]=="no":
	sys_loop_map[c3][c4]="yes"
	# 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]
	for c2 in range(len(sys_loops[loop_row][loop_row])):
	if (sys_loops[loop_row][loop_row][c2][:-1]==branch_info[single_branches[c1]][:-1]):
	if sys_loops[loop_row][loop_row][c2][-1]=="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 is to recalculate the sys_loops matrix.
	# First the diagonal loops are recalculated.
	# Then the off-diagonal (interactions)
	readjust_sys_loops(sys_loops, branch_info, stiff_info)
	# 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 A, B and E
	for c1 in range(len(sys_loops)):
	for c2 in range(len(sys_loops)):
	for c3 in range(len(sys_loops[c1][c2])):
	for c4 in range(len(branch_info)):
	if sys_loops[c1][c2][c3][:-1]==branch_info[c4][:-1]:
	if c1==c2:
	matrix_a.data[c1][c2]+=branch_info[c4][-1][0][0]
	matrix_e.data[c1][c2]+=branch_info[c4][-1][0][1]
	if sys_loops[c1][c2][c3][-1]=="forward":
	for c5 in range(matrix_b.columns):
	matrix_b.data[c1][c5]+=branch_info[c4][-1][1][c5]
	else:
	for c5 in range(matrix_b.columns):
	matrix_b.data[c1][c5]-=branch_info[c4][-1][1][c5]
	else:
	if sys_loops[c1][c2][c3][-1]=="forward":
	matrix_a.data[c1][c2]+=branch_info[c4][-1][0][0]
	matrix_e.data[c1][c2]+=branch_info[c4][-1][0][1]
	else:
	matrix_a.data[c1][c2]-=branch_info[c4][-1][0][0]
	matrix_e.data[c1][c2]-=branch_info[c4][-1][0][1]
	# Make sure that the stiff loops have no inductances
	# so that they are treated as static equations.
	for c1 in range(len(sys_loops)):
	for c2 in range(len(sys_loops[c1][c1])):
	for c3 in range(len(branch_info)):
	if sys_loops[c1][c1][c2][:-1]==branch_info[c3][:-1]:
	if stiff_info[c3]=="yes":
	for c4 in range(matrix_e.columns):
	matrix_e.data[c1][c4]=0.0
	return

view raw compute_loops.py hosted with ❤ by GitHub

This function looks to isolate branches to single loops and once that happens, the current of that branch becomes the loop current. If you are lucky, and in the case of the buck converter I was, the loops will be such that every non stiff loop will have a non stiff branch in that loop only and no other. Otherwise, a row operation is performed to isolate branches.

Anyway, long post but it is long overdue. The buck converter is now working. I am not sure if this will work in its present form for all power electronic converters. Need to test it thoroughly. But I'll release the next version soon.

Getting the diode to freewheel - II

To begin with, today marks the completion of one year of blogging on this project. The project has come pretty far though not as far as I had hoped.

The problem with the approach that I have been using is that it is based exclusively on loop analysis. This breaks almost every time with every modification to the logic I try to make for every variation of the circuit. Nothing universal seems to be coming out. I have been spending the entire week drawing circuits, writing loops, scratching them out and repeating over and over again.

So need to change my approach. How about I mix in some nodal analysis as well? The loop analysis gives me all the loop currents for any one iteration. The control changes the circuit. Now, to determine how the currents will be rearranged, nodal analysis seems to be the only solution.

I am thinking about it this way:

1. At the end of an iteration, mark which loops were non-stiff and therefore may have had inductors carrying a non-negligible current.

2. When an event occurs, we need to rearrange the circuit, the loops, the loop currents and also decide any freewheeling effects.

3. At this point, the one circuit law that should be obeyed is - the current through an inductor cannot change instantaneously. Using this as a starting point, apply nodal analysis to calculate the current through all the branches.

4. Based on the preliminary calculations, check if the branch currents result in any changes in circuit topology - diodes freewheeling or IGBTs conducting. This seems a bit complicated - will need to expand on this.

5. Now that circuit topology is updated, another nodal analysis will update the branch currents. Using these updated branch currents, the loop currents will need to be calculated for the next round of loop analysis.


Only worry is that these tasks seem computationally heavy. Even if they are executed only when events occur, every effort will need to be made to make them simpler.

Getting the diode to freewheel

The diode must freewheel when:

1. An inductor (for now, but in general an energy storage element) contains non-negligible current because a negligible current could simply mean it is in a loop that is stiff but finds itself only in stiff loops.

2. For this current to break would be either bad design which means undue voltage stress on devices as stored energy has nowhere to go or the current freewheels through a diode (or any other device such as an IGBT with an anti-parallel diode).

And this is essentially how the algorithm develops:

1. After one iteration is over, suppose the ideal switch in the buck converter turns off.

2. The next iteration will throw up an event. So new loop manipulations will be done, an effort will be made to remove stiffness, and this effort will fail because there are two parallel stiff branches.

3. The check to be done - is there any inductor in all the loops which has a non-negligible current but now finds itself only in stiff loops? Right now I will restrict the search to inductors but I will expand the criteria later if needed.

4. If so, check if it can freewheel through a device capable of freewheeling. Now this is the tough part.

5. Start with the loop that contains the inductor but has a non-negligible current. This is because the inductor could have been in multiple loops to begin with and some of them may have been stiff in the first place.

6. Check if any of the devices in the loop are capable of freewheeling. This will be done by running a function within all devices but will be empty for all devices that can't freewheel and will contain code for those than can freewheel. The function will check if the loop current and direction are such that the device will start conducting if the current were to flow through it. If so, the device state changes.

7. Now check again if loop containing the inductor is stiff. Why? Because there may be two diodes connected in series. So repeat the above procedure.

8. There could be another possibility. Suppose there are two diodes connected in series in the buck converter and they have high resistance snubbers connected in parallel with them. So now we won't know how the loops are formed. It could be that the inductor appears in a loop with the snubbers and not the diodes. The diodes appear in a loop with the snubbers. So the above algorithm will give up because the inductor loop never encounters a diode.

9. So if upon encountering a stiff branch, it so turns out that the stiff branch does not contain an element that can freewheel, check if that branch exists in other loops. If so, perform a loop manipulation and bring those other loops into the inductor loop. Now repeat the above process.

10. When will the algorithm give up? Not sure how this will work. I will have to run different cases and check. There may be times when it goes into an infinite loop.

Now I will code this. Not sure if I can finish it tonight. Probably will post again tomorrow.
 

Ideal switch

It is winter again and I am buried in snow in Toronto. I am not much into winter activities and this means only coding can stop me for going crazy. So I am back again to this project. After weeks of dithering, pretending to review code and thinking of future strategies, I finally gotten around to taking this project further.

Like I said before, the diode model is not complete because the diode doesn't know when and where to freewheel. To highlight this, I made the ideal switch model and tried to simulate a buck converter. Here's the code for the ideal switch (click on view raw below the code box to see the code in another window):

	class Switch:
	""" Ideal switch class. Contains functions to initiliaze
	the switch according to name tag, unique cell position,
	update system matrix on each iteration. """
	def __init__(self, switch_index, switch_pos, switch_tag):
	""" Constructor to initialize value.
	Also, takes in the identifiers -
	index (serial number), cell position and tag. """
	self.type="Switch"
	self.number=switch_index
	self.pos=switch_pos
	self.tag=switch_tag
	self.has_voltage="yes"
	self.switch_level=120.0
	self.current=0.0
	self.voltage=0.0
	self.polrty=[-1, -1]
	self.resistor_on=0.01
	self.status="off"
	self.control_tag=["Control"]
	self.control_values=[0.0]
	def display(self):
	print "Switch is ",
	print self.tag,
	print " located at ",
	print self.pos,
	print " with negative polarity towards %s" %(csv_element(self.polrty))
	return
	def ask_values(self, x_list, ckt_mat, sys_branch):
	""" Writes the values needed to the spreadsheet."""
	switch_params=["Switch"]
	switch_params.append(self.tag)
	switch_params.append(self.pos)
	switch_params.append("Voltage level (V) = %f" %self.switch_level)
	if self.polrty==[-1, -1]:
	# Looking for a default value of polarity
	# in the neighbouring cells
	self.switch_elem=csv_tuple(self.pos)
	if self.switch_elem[0]>0:
	if ckt_mat[self.switch_elem[0]-1][self.switch_elem[1]]:
	self.polrty=[self.switch_elem[0]-1, self.switch_elem[1]]
	if self.switch_elem[1]>0:
	if ckt_mat[self.switch_elem[0]][self.switch_elem[1]-1]:
	self.polrty=[self.switch_elem[0], self.switch_elem[1]-1]
	if self.switch_elem[0]<len(ckt_mat)-1:
	if ckt_mat[self.switch_elem[0]+1][self.switch_elem[1]]:
	self.polrty=[self.switch_elem[0]+1, self.switch_elem[1]]
	if self.switch_elem[1]<len(ckt_mat)-1:
	if ckt_mat[self.switch_elem[0]][self.switch_elem[1]+1]:
	self.polrty=[self.switch_elem[0], self.switch_elem[1]+1]
	else:
	for c1 in range(len(sys_branch)):
	if csv_tuple(self.pos) in sys_branch[c1]:
	if not self.polrty in sys_branch[c1]:
	print
	print "!"*50
	print "ERROR!!! Switch polarity should be in the same branch as the switch. Check switch at %s" %self.pos
	print "!"*50
	print
	switch_params.append("Negative polarity towards (cell) = %s" %csv_element(self.polrty))
	switch_params.append("Name of control signal = %s" %self.control_tag[0])
	print switch_params
	x_list.append(switch_params)
	return
	def get_values(self, x_list, ckt_mat):
	""" Takes the parameter from the spreadsheet."""
	self.switch_level=float(x_list[0].split("=")[1])
	# Choosing 1 micro Amp as the leakage current that
	# is drawn by the switch in off state.
	self.resistor_off=self.switch_level/1.0e-6
	self.resistor=self.resistor_off
	switch_polrty=x_list[1].split("=")[1]
	# Convert the human readable form of cell
	# to [row, column] form
	while switch_polrty[0]==" ":
	switch_polrty=switch_polrty[1:]
	self.polrty=csv_tuple(switch_polrty)
	if not ckt_mat[self.polrty[0]][self.polrty[1]]:
	print "Polarity incorrect. Branch does not exist at %s" %csv_element(self.polrty)
	self.control_tag[0]=x_list[2].split("=")[1]
	while self.control_tag[0][0]==" ":
	self.control_tag[0]=self.control_tag[0][1:]
	return
	def transfer_to_sys(self, sys_loops, mat_e, mat_a, mat_b, mat_u, source_list):
	""" The matrix A in E.dx/dt=Ax+Bu will be updated by the
	resistor value of the switch."""
	for c1 in range(len(sys_loops)):
	for c2 in range(c1, len(sys_loops)):
	# Updating the elements depending
	# on the sense of the loops (aiding or opposing)
	for c3 in range(len(sys_loops[c1][c2])):
	# Check if current source position is there in the loop.
	if csv_tuple(self.pos) in sys_loops[c1][c2][c3]:
	# Add current source series resistor
	# if branch is in forward direction
	if sys_loops[c1][c2][c3][-1]=="forward":
	mat_a.data[c1][c2]+=self.resistor
	else:
	# Else subtract if branch is in reverse direction
	mat_a.data[c1][c2]-=self.resistor
	# Because the matrices are symmetric
	mat_a.data[c2][c1]=mat_a.data[c1][c2]
	# If the positive polarity appears before the voltage position
	# it means as per KVL, we are moving from +ve to -ve
	# and so the voltage will be taken negative
	if sys_loops[c1][c1][c2].index(self.polrty)>sys_loops[c1][c1][c2].index(csv_tuple(self.pos)):
	if sys_loops[c1][c1][c2][-1]=="forward":
	mat_b.data[c1][source_list.index(self.pos)]=-1.0
	else:
	mat_b.data[c1][source_list.index(self.pos)]=1.0
	else:
	if sys_loops[c1][c1][c2][-1]=="forward":
	mat_b.data[c1][source_list.index(self.pos)]=1.0
	else:
	mat_b.data[c1][source_list.index(self.pos)]=-1.0
	return
	def transfer_to_branch(self, sys_branch, source_list):
	""" Update the resistor info of the switch
	to the branch list """
	if csv_tuple(self.pos) in sys_branch:
	sys_branch[-1][0][0]+=self.resistor
	# For the switch forward voltage drop.
	if csv_tuple(self.pos) in sys_branch:
	if sys_branch.index(self.polrty)>sys_branch.index(csv_tuple(self.pos)):
	sys_branch[-1][1][source_list.index(self.pos)]=-1.0
	else:
	sys_branch[-1][1][source_list.index(self.pos)]=1.0
	return
	def generate_val(self, source_lst, sys_loops, mat_e, mat_a, mat_b, mat_u, t, dt):
	""" The switch forward drop voltage is updated
	in the matrix u in E.dx/dt=Ax+Bu ."""
	if self.status=="on":
	mat_u.data[source_lst.index(self.pos)][0]=0.7
	else:
	mat_u.data[source_lst.index(self.pos)][0]=0.0
	return
	def update_val(self, sys_loops, lbyr_ratio, mat_e, mat_a, mat_b, state_vec, mat_u, sys_branches, sys_events):
	""" This function calculates the actual current in the
	switch branch. With this, the branch voltage is found
	with respect to the existing switch resistance. The switch
	voltage is then used to decide the turn on condition. """
	# Local variable to calculate the branch
	# current from all loops that contain
	# the current source branch.
	act_current=0.0
	for c1 in range(len(sys_loops)):
	for c2 in range(len(sys_loops[c1][c1])):
	if csv_tuple(self.pos) in sys_loops[c1][c1][c2]:
	# If switch negative polarity is after the switch
	# position, the current is positive.
	if sys_loops[c1][c1][c2].index(self.polrty)>sys_loops[c1][c1][c2].index(csv_tuple(self.pos)):
	# Then check is the loop is aiding or opposing
	# the main loop.
	if sys_loops[c1][c1][c2][-1]=="forward":
	act_current+=state_vec.data[c1][0]
	else:
	act_current-=state_vec.data[c1][0]
	else:
	if sys_loops[c1][c1][c2][-1]=="forward":
	act_current-=state_vec.data[c1][0]
	else:
	act_current+=state_vec.data[c1][0]
	self.current=act_current
	self.voltage=self.current*self.resistor
	# Identifying the position of the switch branch
	# to generate events.
	for c1 in range(len(sys_branches)):
	if csv_tuple(self.pos) in sys_branches[c1]:
	branch_pos=c1
	# Switch will turn on when it is forward biased
	# and it is gated on.
	if self.control_values[0]>=1.0 and self.voltage>1.0:
	if self.status=="off":
	sys_events[branch_pos]="yes"
	self.status="on"
	# Switch will turn off when gated off or
	# when current becomes negative.
	if self.control_values[0]==0.0 or self.current<0.0:
	if self.status=="on":
	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 switch.py hosted with ❤ by GitHub

This is the circuit for the buck converter.

And this is the control code - actually just an open loop with a constant 50% duty ratio  (click on view raw below the code box to see the code in another window):

	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=10.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+10.0e-6

view raw buck_control.py hosted with ❤ by GitHub

And on execution, the current (measured by Ammeter_La) through the inductor La is:

And as can be seen, the diode does not free wheel at all. The current rises when the switch turns on and drops to zero when it is turned off.

The strategy to overcome this will be my next blog entry as it is a fairly detailed algorithm.