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.