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.