Version 0.2.1 released

Released the next version of the simulator with diode model:
http://sourceforge.net/projects/pythonpowerelec/

As always, for questions or comments, please send me an email at pythonpowerelectronics@gmail.com

Diode model - loops and manipulations

For any circuit, the number of independent loops is Branches-Nodes+1. The loop finder function returns excess loops for all circuits except maybe the simplest circuits. Initially, I had removed the excess loops and assumed they were linear combinations of the previous loops. But when the circuit has a ladder structure, it may so happen that some of the essential loops appear right at the end and get thrown out.

This is what I found when I was trying to figure out why the three-phase diode bridge rectifier was not working. The only way to ensure that essential loops are not deleted is to let all of them be to begin with. With row operations, some loops will be eliminated. As and how loops are eliminated, they should be completely deleted and the system should be reduced.

So, first the loop finder function (click on "View Raw" below the code box to see the code in a new window):

	def find_loop(br_map, nd_list, lp_list, lp_iter, elem, lp_count, lp_limit):
	"""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_count, lp_iter=loop_horiz(br_map, nd_list, lp_list, lp_iter, elem, lp_count)
	# Change direction to vertical
	loop_dir="vert"
	# Will be executed if we are moving vertically
	if (loop_dir == "vert"):
	lp_count, lp_iter=loop_vert(br_map, nd_list, lp_list, lp_iter, elem, lp_count)
	# Change direction to horizontal
	loop_dir="horiz"
	return lp_count

view raw loop_finder.py hosted with ❤ by GitHub

The difference is in the terminating condition:
if loop_iter:
As long as loop_iter is searching, let it search. So let it add as many valid loops as possible.
The number of excess loops can be pretty huge (x6).

Next comes the main circuit_solver.py. In this, another matrix has been conceived called the system_loop_map. This is to indicate which branches are stiff so as to eliminate stiff branches from as many loops as possible. The code for this has been put together in one block (click on "View Raw" below the code box to see the code in a new window):

	# Stiff ratio indicates whether a branch
	# is stiff.
	stiff_ratio=[]
	for c1 in range(len(branch_params)):
	stiff_ratio.append("no")
	max_res=abs(branch_params[0][-1][0][0])
	for c1 in range(1, len(branch_params)):
	if branch_params[c1][-1][0][0]>max_res:
	max_res=branch_params[c1][-1][0][0]
	# Calculates the time constants of the loops
	# from the diagonal elements as L/R ratios.
	for c1 in range(len(branch_params)):
	if branch_params[c1][-1][0][0]:
	if abs(branch_params[c1][-1][0][1]/branch_params[c1][-1][0][0])<0.1*dt:
	if branch_params[c1][-1][0][0]/max_res > dt:
	stiff_ratio[c1]="yes"
	system_loop_map=[]
	for c1 in range(len(system_loops_copy)):
	br_vector=[]
	for c2 in range(len(branch_params)):
	br_vector.append("no")
	system_loop_map.append(br_vector)
	for c1 in range(len(system_loops_copy)):
	for c3 in range(len(branch_params)):
	for c2 in range(len(system_loops_copy[c1][c1])):
	if branch_params[c3][:-1]==system_loops_copy[c1][c1][c2][:-1]:
	if stiff_ratio[c3]=="yes":
	system_loop_map[c1][c3]="stiff"
	else:
	system_loop_map[c1][c3]="yes"

view raw system_loop_map.py hosted with ❤ by GitHub

So essentially, for every loop there is minimal information about every branch in the circuit - if it exists and if it does, is the branch stiff.

	# The concept here is to minimize the number of times
	# a stiff branch appears in loops. To do so a sys_loop_map
	# has been constructed. This has three options "stiff", "yes"
	# and "no". Using row operations, this sys_loop_map is made
	# into a upper triangular matrix only with respect to the
	# stiff branches.
	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][2]=="stiff":
	is_loop_stiff="yes"
	# Check if there is another loop after c1 (>c1)
	# which has a stiff branch "sooner" than c1 has.
	# All this is mere nomenclature.
	# If so, exchange the loops for the triangularization.
	if is_loop_stiff=="yes":
	c2=0
	while c2<len(branch_info) and sys_loop_map[c1][c2]!="stiff":
	for c3 in range(c1+1, len(sys_loop_map)):
	if sys_loop_map[c3][c2]=="stiff":
	c2=len(branch_info)-1
	for c4 in range(len(branch_info)):
	sys_loop_map[c1][c4], sys_loop_map[c3][c4] = sys_loop_map[c3][c4], sys_loop_map[c1][c4]
	sys_loops[c1][c1], sys_loops[c3][c3] = sys_loops[c3][c3], sys_loops[c1][c1]
	c2=c2+1
	# Look for the first stiff branch
	c2=0
	while c2<len(branch_info) and sys_loop_map[c1][c2]!="stiff":
	c2=c2+1
	if c2<len(branch_info):
	# Eliminate that stiff branch from subsequent loops.
	for c3 in range(c1+1, len(sys_loop_map)):
	if sys_loop_map[c3][c2]=="stiff":
	# Loop manipulations may cause loops
	# to change lenghts rapidly. So the try/except
	# statmements are used to avoid exceeding indices.
	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 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
	# 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"
	elif sys_loop_map[c1][c4]=="stiff" and sys_loop_map[c3][c4]=="stiff":
	sys_loop_map[c3][c4]="no"
	elif sys_loop_map[c1][c4]=="stiff" and sys_loop_map[c3][c4]=="no":
	sys_loop_map[c3][c4]="stiff"
	# This is to make sure, that stiff loops
	# are connected to an input.
	for c1 in range(len(sys_loops)):
	is_loop_stiff="no"
	has_input="no"
	for c2 in range(len(sys_loops[c1][c1])):
	# Check if any of the branches are stiff
	for c3 in range(len(branch_info)):
	if branch_info[c3][:-1]==sys_loops[c1][c1][c2][:-1]:
	if stiff_info[c3]=="yes":
	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_loops):
	if not c1==c2:
	is_loop_stiff="no"
	flag_input="no"
	for c3 in range(len(sys_loops[c2][c2])):
	for c4 in range(len(branch_info)):
	if branch_info[c4][:-1]==sys_loops[c2][c2][c3][:-1]:
	if stiff_info[c4]=="yes":
	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")

view raw stiffness_removal.py hosted with ❤ by GitHub

So essentially, I use system_loops_map to make the system upper triangular as far as stiff branches are concerned. The next block I am not sure if it is needed. Whether it is necessary to make sure a stiff loop is connected to the input. I'll test it and try to get rid of it. For some reason, it looks like an ugly code block.

The last part of to reduce the size of the system by getting rid of redundant loops (click on "View Raw" below the code box to see the code in a new window):

	# 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]
	# This part if to eliminate the empty (redundant) loops
	# based on zero rows in A, E and B matrices.
	# At the same time, A, B and E matrices are resized.
	for c1 in range(len(sys_loops)-1, -1, -1):
	empty_loop="yes"
	for c2 in range(matrix_a.columns):
	if matrix_a.data[c1][c2]:
	empty_loop="no"
	for c2 in range(matrix_e.columns):
	if matrix_e.data[c1][c2]:
	empty_loop="no"
	for c2 in range(matrix_b.columns):
	if matrix_b.data[c1][c2]:
	empty_loop="no"
	if empty_loop=="yes":
	# If a loop is empty, the column c1 is first deleted
	# and then the row c1 is deleted.
	for c2 in range(len(sys_loops)-1, -1, -1):
	for c3 in range(len(sys_loops[c2][c1])-1, -1, -1):
	del sys_loops[c2][c1][c3]
	del sys_loops[c2][c1]
	for c2 in range(len(sys_loops[c1])-1, -1, -1):
	del sys_loops[c1][c2]
	del sys_loops[c1]
	# Same for the matrices A and E
	for c2 in range(matrix_a.rows-1, -1, -1):
	del matrix_a.data[c2][c1]
	del matrix_e.data[c2][c1]
	matrix_a.columns-=1
	matrix_e.columns-=1
	for c2 in range(matrix_a.columns-1, -1, -1):
	del matrix_a.data[c1][c2]
	del matrix_e.data[c1][c2]
	del matrix_a.data[c1]
	del matrix_e.data[c1]
	matrix_a.rows-=1
	matrix_e.rows-=1
	for c2 in range(matrix_b.columns-1, -1, -1):
	del matrix_b.data[c1][c2]
	del matrix_b.data[c1]
	matrix_b.rows-=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

view raw sys_red.py hosted with ❤ by GitHub

Diode model - part I

Where do I begin? Massive number of changes. So first, the diode class (clink on "View Raw" below the code box to see the code in a new window):

	class Diode:
	""" Diode class. Contains functions to initiliaze
	the diode according to name tag, unique cell position,
	update system matrix on each iteration. """
	def __init__(self, diode_index, diode_pos, diode_tag):
	""" Constructor to initialize value.
	Also, takes in the identifiers -
	index (serial number), cell position and tag. """
	self.type="Diode"
	self.number=diode_index
	self.pos=diode_pos
	self.tag=diode_tag
	self.has_voltage="yes"
	self.diode_level=120.0
	self.current=0.0
	self.voltage=0.0
	self.polrty=[-1, -1]
	self.resistor_on=0.01
	self.status="off"
	def display(self):
	print "Diode is ",
	print self.tag,
	print " located at ",
	print self.pos,
	print " with cathode 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."""
	diode_params=["Diode"]
	diode_params.append(self.tag)
	diode_params.append(self.pos)
	diode_params.append("Voltage level (V) = %f" %self.diode_level)
	if self.polrty==[-1, -1]:
	# Looking for a default value of polarity
	# in the neighbouring cells
	self.diode_elem=csv_tuple(self.pos)
	if self.diode_elem[0]>0:
	if ckt_mat[self.diode_elem[0]-1][self.diode_elem[1]]:
	self.polrty=[self.diode_elem[0]-1, self.diode_elem[1]]
	if self.diode_elem[1]>0:
	if ckt_mat[self.diode_elem[0]][self.diode_elem[1]-1]:
	self.polrty=[self.diode_elem[0], self.diode_elem[1]-1]
	if self.diode_elem[0]<len(ckt_mat)-1:
	if ckt_mat[self.diode_elem[0]+1][self.diode_elem[1]]:
	self.polrty=[self.diode_elem[0]+1, self.diode_elem[1]]
	if self.diode_elem[1]<len(ckt_mat)-1:
	if ckt_mat[self.diode_elem[0]][self.diode_elem[1]+1]:
	self.polrty=[self.diode_elem[0], self.diode_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!!! Diode polarity should be in the same branch as the diode. Check diode at %s" %self.pos
	print "!"*50
	print
	diode_params.append("Cathode polarity towards (cell) = %s" %csv_element(self.polrty))
	x_list.append(diode_params)
	return
	def get_values(self, x_list, ckt_mat):
	""" Takes the parameter from the spreadsheet."""
	self.diode_level=float(x_list[0].split("=")[1])
	# Choosing 1 micro Amp as the leakage current that
	# is drawn by the diode in off state.
	self.resistor_off=self.diode_level/1.0e-6
	self.resistor=self.resistor_off
	diode_polrty=x_list[1].split("=")[1]
	# Convert the human readable form of cell
	# to [row, column] form
	while diode_polrty[0]==" ":
	diode_polrty=diode_polrty[1:]
	self.polrty=csv_tuple(diode_polrty)
	if not ckt_mat[self.polrty[0]][self.polrty[1]]:
	print "Polarity incorrect. Branch does not exist at %s" %csv_element(self.polrty)
	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 diode."""
	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 diode
	to the branch list """
	if csv_tuple(self.pos) in sys_branch:
	sys_branch[-1][0][0]+=self.resistor
	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 diode 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
	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. """
	# 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 diode cathode polarity is after the diode
	# 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
	for c1 in range(len(sys_branches)):
	if csv_tuple(self.pos) in sys_branches[c1]:
	branch_pos=c1
	if self.status=="off" and self.voltage>0.9:
	sys_events[branch_pos]="yes"
	self.status="on"
	if self.status=="on" and self.current<-1.0e-5:
	sys_events[branch_pos]="yes"
	self.status="off"
	if self.current<-1.0e-5:
	self.status="off"
	if self.status=="off":
	self.resistor=self.resistor_off
	else:
	self.resistor=self.resistor_on
	return

view raw diode_class.py hosted with ❤ by GitHub

Similar to many of the others. In the parameter specification, the user needs to enter the voltage level of the diode and the polarity in terms of where the cathode is.

The diode resistance changes with the status whether "ON" or "OFF". The diode is ON when it is forward biased beyond a certain threshold voltage. It turns "OFF" when current becomes negative. The ON resistance, forward bias threshold voltage and the ON drop voltage can be made user defined parameters.