Stiff Systems - III

Several changes made to the entire codes. In order to solve every type of stiff system effectively, the only was seems to be to completely recalculate not only the A, B and E matrices but the KVL loops and their interactions.

At the highest level, I defined another layer called the branch. Currently, it is another list, but it may be an object depending on the later use. I have defined the branch as follows (click on "view raw" to see the code in a new window):

	### This list contains every branch and
	### at the end contains the elements that
	### would go into A, B and E matrices.
	branch_params=[]
	for c1 in range(len(system_loops)):
	# Check if a branch has already been added to
	# branch_params.
	for c2 in range(len(system_loops[c1][c1])):
	branch_found="no"
	for c4 in range(len(branch_params)):
	if system_loops[c1][c1][c2][:-1]==branch_params[c4][:-1]:
	branch_found="yes"
	# Matrix B will have the number
	# of entries equal to the number
	# of sources.
	if branch_found=="no":
	source_add=[]
	for c3 in range(len(source_list)):
	source_add.append(0.0)
	branch_add=system_loops[c1][c1][c2][:-1]
	params_add=[]
	params_add.append([0.0,0.0])
	params_add.append(source_add)
	branch_add.append(params_add)
	branch_params.append(branch_add)
	### This is an additional layer of check to make sure
	### that none of the branches in branch_params is the
	### reverse of another.
	for c1 in range(len(branch_params)-1, -1, -1):
	# Start from the last branch and compute
	# the reverse.
	# Make a copy
	check_br_reverse=[]
	for c2 in range(len(branch_params[c1])-1):
	check_br_reverse.append(branch_params[c1][c2])
	# Reverse the copy
	for c2 in range(len(check_br_reverse)/2):
	check_br_reverse[c2], check_br_reverse[len(check_br_reverse)-1-c2] = \
	check_br_reverse[len(check_br_reverse)-1-c2], check_br_reverse[c2]
	# Iterate through the branches.
	c2=0
	while c2<len(branch_params):
	# Check if the reverse of the branch c1
	# is equal to the current branch
	if check_br_reverse==branch_params[c2][:-1]:
	# In this case, two things need to be done
	# Find out where the branch occurs in
	# system_loops and reverse those branches
	# Second, delete that branch from branch_params
	# There will be two instances of the branch
	# One with a forward appended at the end
	check_br_fwd=[]
	for c3 in range(len(branch_params[c1][:-1])):
	check_br_fwd.append(branch_params[c1][c3])
	check_br_fwd.append("forward")
	# The other with a reverse appended at the end.
	check_br_rev=[]
	for c3 in range(len(branch_params[c1][:-1])):
	check_br_rev.append(branch_params[c1][c3])
	check_br_rev.append("reverse")
	# Look through the entire system_loops matrix
	for c3 in range(len(system_loops)):
	for c4 in range(len(system_loops)):
	# Check if the branch with forward is found.
	if (check_br_fwd in system_loops[c3][c4]):
	# Mark the position
	ex_br_pos=system_loops[c3][c4].index(check_br_fwd)
	# Mark the direction
	ex_br_dir=system_loops[c3][c4][ex_br_pos][-1]
	# Delete the branch
	del system_loops[c3][c4][ex_br_pos]
	# Add the earlier branch found branch_params[c2]
	# But reverse the direction of the branch just
	# deleted.
	if ex_br_dir=="forward":
	new_br=[]
	for c5 in range(len(branch_params[c2])-1):
	new_br.append(branch_params[c2][c5])
	new_br.append("forward")
	system_loops[c3][c4].append(new_br)
	else:
	new_br=[]
	for c5 in range(len(branch_params[c2])-1):
	new_br.append(branch_params[c2][c5])
	new_br.append("reverse")
	system_loops[c3][c4].append(new_br)
	# Repeat the process with the branch with reverse
	if (check_br_rev in system_loops[c3][c4]):
	ex_br_pos=system_loops[c3][c4].index(check_br_rev)
	ex_br_dir=system_loops[c3][c4][ex_br_pos][-1]
	del system_loops[c3][c4][ex_br_pos]
	if ex_br_dir=="forward":
	new_br=[]
	for c5 in range(len(branch_params[c2])-1):
	new_br.append(branch_params[c2][c5])
	new_br.append("forward")
	system_loops[c3][c4].append(new_br)
	else:
	new_br=[]
	for c5 in range(len(branch_params[c2])-1):
	new_br.append(branch_params[c2][c5])
	new_br.append("reverse")
	system_loops[c3][c4].append(new_br)
	# Delete the latest branch because it is a reverse
	del branch_params[c1]
	c2+=1

view raw branch_fed.py hosted with ❤ by GitHub

There was a problem where some branches were the reverse of the other and these were difficult to handle when loops were being manipulated. So I added the code to delete the redundant branches and also remove them from the system_loops matrix.

Another concept planned for the next stage was that each branch can have an "event" flag. That is if there is a change in the parameters like a diode starting or stopping conduction, there will be an event generated. When an event is generated on a branch, an event is generated to the entire system and only then will the system matrices A, B and E be recalculated. Otherwise, leave them as it is and just solve the equations for new inputs. Even in a power electronics circuit, this may result in significant savings as a change may happen once in 20-50 microseconds so that with a time step of 1 microsecond, the computation can be 20-50x faster if all the recalculation are left out.

The branches are updated as follows (click on "view raw" to see the code in a new window):

	# Now transfer values from the branch to the matrices A,B and E.
	# Look where each branch appears and transfer the info.
	for c1 in range(len(system_loops_copy)):
	for c2 in range(len(system_loops_copy)):
	for c3 in range(len(system_loops_copy[c1][c2])):
	for c4 in range(len(branch_params)):
	if system_loops_copy[c1][c2][c3][:-1]==branch_params[c4][:-1]:
	if c1==c2:
	sys_mat_a.data[c1][c2]+=branch_params[c4][-1][0][0]
	sys_mat_e.data[c1][c2]+=branch_params[c4][-1][0][1]
	if system_loops_copy[c1][c2][c3][-1]=="forward":
	for c5 in range(len(source_list)):
	sys_mat_b.data[c1][c5]+=branch_params[c4][-1][1][c5]
	else:
	for c5 in range(len(source_list)):
	sys_mat_b.data[c1][c5]-=branch_params[c4][-1][1][c5]
	else:
	if system_loops_copy[c1][c2][c3][-1]=="forward":
	sys_mat_a.data[c1][c2]+=branch_params[c4][-1][0][0]
	sys_mat_e.data[c1][c2]+=branch_params[c4][-1][0][1]
	else:
	sys_mat_a.data[c1][c2]-=branch_params[c4][-1][0][0]
	sys_mat_e.data[c1][c2]-=branch_params[c4][-1][0][1]

view raw branch_update.py hosted with ❤ by GitHub

Also, each class has another method called "transfer_to_branch".

The really huge change has been in the solver.py file. The basic concept is described here. The file has comments as far as possible.

1. The way to find out if a loop is stiff is to check the L/R ratio. If the L/R ratio of the diagonal element is low, the offdiagonal elements are checked. If they are also stiff, a row manipulation is done with the lower loops so that the stiff elements are removed from the lower loops.

2. Once the stiff elements have been isolated to individual loops by row operations, the loops are now recalculated.

3. Once again the system matrices A, B and E are then calculated from the new loops and now the ODE is solved.