Thursday, March 21, 2013

Version 0.1.2 released.

Released the next minor patch:
http://sourceforge.net/projects/pythonpowerelec/files/
Posted by at No comments:

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.


Posted by at No comments:

Sunday, March 10, 2013

Stiff systems - II

Stupid mistake. The algorithm wasn't totally wrong. Still is wrong but not completely wrong. The problem was in the loop interaction.

Take this sample circuit:


The large resistors are Resistor_V2 (at 4D and of 500000.0 ohm) and Resistor_V1 (at 12 J and of 10000000.0 ohm).

The matrices turned out to be:
-----------------------------------------------------------------------------------------------------------
500000.0  500000.0  0.0  500000.0  500000.0
500000.0  500010.1  0.0  -500000.1  500000.1
0.0  0.0  10000010.0  0.0  -10000000.0
500000.0  -500000.1  0.0  500020.1  -499990.1
500000.0  500000.1  -10000000.0  -499990.1  10500010.1

0.0  0.0  0.0  0.0  0.0
0.0  0.0001  0.0  -0.0001  0.0001
0.0  0.0  0.1  0.0  0.0
0.0  -0.0001  0.0  0.2001  0.0999
0.0  0.0001  0.0  0.0999  0.1001

1.0
0.0
0.0
0.0
0.0
-----------------------------------------------------------------------------------------------------------

For loops defined as:
 -----------------------------------------------------------------------------------------------------------
1D 2D 3D 4D 5D 6D 6C 6B 6A 5A 4A 3A 2A 1A 1B 1C 1D

1H 1G 1F 1E 1D 2D 3D 4D 5D 6D 6E 6F 6G 6H 5H 4H 3H 2H 1H

9L 10L 11L 12L 12K 12J 12I 12H 11H 10H 9H 9I 9J 9K 9L

1H 1G 1F 1E 1D 2D 3D 4D 5D 6D 6E 6F 6G 6H 6I 6J 6K 6L 5L 4L 3L 2L 1L 1K 1J 1I 1H

1H 1G 1F 1E 1D 2D 3D 4D 5D 6D 6E 6F 6G 6H 7H 8H 9H 10H 11H 12H 12I 12J 12K 12L 11L 10L 9L 8L 7L 6L 5L 4L 3L 2L 1L 1K 1J 1I 1H
-----------------------------------------------------------------------------------------------------------

Which is wrong. Because, the interaction between loop 2 and loop 4 is -500000.1 in the matrix A but should be positive from the sense of the loops above. The problem is in the code that defines the "forward" and "reverse" directions.

This is the new code (click on "view raw" to see the code in a new window):


# A temporary list that stores the nodes
# common to two loops
nodes_in_loop=[]
# A array of all the loops of the system
# including common branches between loops.
system_loops=[]
for c1 in range(loop_count):
row_vector=[]
for c2 in range(loop_count):
row_vector.append([])
system_loops.append(row_vector)
for c1 in range(len(loop_branches)):
# The diagonal elements of system_loops
# will be the loops themselves.
for c2 in range(len(loop_branches[c1])):
system_loops[c1][c1].append(loop_branches[c1][c2])
# The system_loops array will be symmetric
for c2 in range(c1+1, len(loop_branches)):
# Find the nodes in loop1
for c3 in range(len(node_list)):
if node_list[c3] in loop_branches[c1]:
nodes_in_loop.append(node_list[c3])
# Find out the nodes common to loop1
# and loop2.
for c3 in range(len(nodes_in_loop)-1, -1, -1):
if nodes_in_loop[c3] not in loop_branches[c2]:
del nodes_in_loop[c3]
# If there are two or more nodes common
# between loop1 and loop2, there are
# common elements.
if len(nodes_in_loop)>1:
comm_seg_in_loop=comm_elem_in_loop(loop_branches[c1], loop_branches[c2])
# The list of common branches between
# loop c1 and loop c2
sys_loop_off_diag=comm_branch_in_loop(comm_seg_in_loop, loop_branches[c1])
for c3 in range(len(sys_loop_off_diag)):
#for c4 in range(len(sys_loop_off_diag[c3])):
system_loops[c1][c2].append(sys_loop_off_diag[c3])
system_loops[c2][c1].append(sys_loop_off_diag[c3])
# Determining the direction of common
# branches between loops c1 and c2
for c3 in range(len(sys_loop_off_diag)):
st_node=sys_loop_off_diag[c3][0]
end_node=sys_loop_off_diag[c3][-1]
# Map the start and end nodes of common
# branch c3 to loops c1 and c2
loop1_st_node=loop_branches[c1].index(st_node)
loop1_end_node=loop_branches[c1].index(end_node)
loop2_st_node=loop_branches[c2].index(st_node)
loop2_end_node=loop_branches[c2].index(end_node)
# This check is because the first node of a
# loop is the same as the last node.
# So this is to make sure, the correct start
# and end nodes are mapped, the difference between
# the indices is compared to the length of the common
# branch. If not equal, this means, the last node of
# a loop needs to be taken instead of the first node.
if (abs(loop1_end_node-loop1_st_node)!= (len(sys_loop_off_diag[c3])-1)):
if (len(loop_branches[c1])-1-loop1_st_node==(len(sys_loop_off_diag[c3])-1)):
loop1_end_node=len(loop_branches[c1])-1
if (len(loop_branches[c1])-1-loop1_end_node==(len(sys_loop_off_diag[c3])-1)):
loop1_st_node=len(loop_branches[c1])-1
if (abs(loop2_end_node-loop2_st_node)!= (len(sys_loop_off_diag[c3])-1)):
if (len(loop_branches[c2])-1-loop2_st_node==(len(sys_loop_off_diag[c3])-1)):
loop2_end_node=len(loop_branches[c2])-1
if (len(loop_branches[c2])-1-loop2_end_node==(len(sys_loop_off_diag[c3])-1)):
loop2_st_node=len(loop_branches[c2])-1
if (loop1_st_node<loop1_end_node):
loop1_dir="ahead"
else:
loop1_dir="back"
if (loop2_st_node<loop2_end_node):
loop2_dir="ahead"
else:
loop2_dir="back"
if loop1_dir=="ahead" and loop2_dir=="ahead":
sys_loop_off_diag[c3].append("forward")
if loop1_dir=="back" and loop2_dir=="back":
sys_loop_off_diag[c3].append("forward")
if loop1_dir=="ahead" and loop2_dir=="back":
sys_loop_off_diag[c3].append("reverse")
if loop1_dir=="back" and loop2_dir=="ahead":
sys_loop_off_diag[c3].append("reverse")
nodes_in_loop=[]
view raw loop_interac_dir.py hosted with ❤ by GitHub

The problem is with the judgement of length of branch with respect to the start and end nodes. The length of the branch will the number of elements in the branch. Which will be (end_node_st_node+1). This is exact block where there was the error:
-----------------------------------------------------------------------------------------------------------
if (abs(loop1_end_node-loop1_st_node)!= (len(sys_loop_off_diag[c3])-1)):
    if (len(loop_branches[c1])-1-loop1_st_node==(len(sys_loop_off_diag[c3])-1)):
        loop1_end_node=len(loop_branches[c1])-1
    if (len(loop_branches[c1])-1-loop1_end_node==(len(sys_loop_off_diag[c3])-1)):
        loop1_st_node=len(loop_branches[c1])-1
if (abs(loop2_end_node-loop2_st_node)!= (len(sys_loop_off_diag[c3])-1)):
    if (len(loop_branches[c2])-1-loop2_st_node==(len(sys_loop_off_diag[c3])-1)):
        loop2_end_node=len(loop_branches[c2])-1
    if (len(loop_branches[c2])-1-loop2_end_node==(len(sys_loop_off_diag[c3])-1)):
        loop2_st_node=len(loop_branches[c2])-1
-----------------------------------------------------------------------------------------------------------


So, the row manipulations to remove the stiff elements from the matrix produce the following matrices:
-----------------------------------------------------------------------------------------------------------
a=
500000.0  0.0  0.0  0.0  0.0
0.0  10.1  0.0  0.0999999999767  0.0999999999767
0.0  0.0  10500010.0  0.0  0.0
0.0  0.0999999999767  0.0  20.1  10.1
0.0  0.0999999999767  10.0  10.1  10.0999999996

e=
0.0  0.0  0.0  0.0  0.0
0.0  0.0001  0.0  0.0001  0.0001
0.0  0.0  0.0  0.0  0.0
0.0  0.0001  0.0  0.2001  0.1001
0.0  0.0001  0.1  0.1001  0.1001

b=
1.0
-1.0
1.0
-1.0
-1.0
-----------------------------------------------------------------------------------------------------------


The isolation has taken place. Because the two large resistors have been isolated to two separate rows. But the loop resistance between the 10000000.0 ohm resistor and the source is wrong. It should be 10000010.1 instead of 10500010.0.

Well, it doesn't blow up like before. So an improvement.

Releasing this as a new patch anyway (v0.1.2). Check on the SourceForge site soon.
Posted by at No comments:

Friday, March 8, 2013

Stiff systems

Changing the title at least to reflect the focus of the problem.

Change in strategy. Make a collection of branches and for each branch have the entries for A, E and B matrices. So this will give an indication of the L, R and V values in every branch.

What is a stiff branch?

  1. It is not a branch with negligible or zero inductance. If the resistance is small enough, it could be a load resistance or even a branch resistance.
  2. It is a branch with a large resistance and negligible or zero inductance.
So, how do know this? It is relative. If all the resistances in the circuit are in a narrow range - say 0.1 to 1000, the system is not stiff. How large a resistance is again depends on the voltage level - for a 120 V system, 100 kiloohm would be large, for a 120 kV system, it would be nothing. But to avoid a per unit system and still classify a system as stiff will mean I need to arrange them in an order according to the resistances. Which is not a difficult task.

After this, I need to calculate the L/R ratio just in case though a very large inductance accompanying a large resistance is not very practical. A large resistance simply means an open circuit or a non-conducting medium but a large inductance means a huge core and a large number of turns in the winding. The worst I can think of in a huge transformer that has a large magnetising inductance but this will have a low winding resistance and a large core loss resistor across it. So it won't be a stiff system. On the contrary it will be a system with a large time constant which is actually the case because when you energize one of these monsters, the magnetizing current can take a long time to decay.

The smallest resistance and inductance would really be the measure of a circuit because they would be the parasitics. So suppose I take these as a base and try to calculate ratios of other elements with respect to it. So a parasitic resistance of 0.005 ohm for a 120 V circuit would be acceptable wire resistance. And a 10 ohm load resistor would result in  a ratio of 10/0.005=2000 which is quite OK. But a 1 megaohm resistor would result in a ratio 2.0e+8 which should raise an alarm.

Taking such a ratio is itself a form of per unit, just that the per unit is arbitrary and if a user puts ridiculous values and an equally impractical simulation time step, a circuit will classified as stiff.
Posted by at No comments:

Thursday, March 7, 2013

Circuit Solver - IX

Too quick to celebrate. One of the reason why I thought things were working was because the system was small and it wasn't too easy to miss out branches. The row operations to get rid of the large resistances and isolate the stiff equations was correct. But the subsequent matrix manipulation to make the matrix symmetric was wrong.

The consider the diagonal element of a manipuled row to be the sum of the offdiagonal elements is no way right. Simply because the diaognal element may already contain a branch present in any of the off diagonal elements and so the effective resistance/inductance may be larger.

I was trying to do this because I didn't want to do another parameter read from the objects. But looks like that may be the only way. Or at least can't be completely avoided. Need to figure out how to minimize the object read.
Posted by at No comments:

Wednesday, March 6, 2013

Circuit Solver - VIII

Documenting what I have been doing for the past two weeks won't be easy. Anyway, the main focus of my problem has been to be able to solve any type of differential equations. When I say any type, this means a combination of differential equations that are stiff and non-stiff. Basically a stiff equation is one that has a small time constant and therefore needs a small simulation time step to capture the dynamics effectively. Something I don't want to do because reducing the time step too much will cause the simulation to be very slow.

So to describe the process from the outer layers of the code to the inner layers.

First major change. The way the current loops of the system and their interaction is decribed. A loop will have the beginning and the ending node as the same and will be a continuous list. The problem is if loop manipulations have to be done - that is differences between loops have to be found in terms of branch segments, it is a little difficult.
So all loops and common branch segments are defined in terms of branches that are elements between two nodes. This is the smallest denomination. Also, every one of them has a direction. If a branch is in a loop and is the diagonal element in system_loops, it has the default direction as "forward". If it shows the interaction between loops, it was have the necessary direction.

So before the ODE solver, a pre-processing has been done to break up all elements in system_loops into lists of branches. This is the code (click on "View Raw" below the code box to see it in a new window):


for c1 in range(len(system_loops)):
for c2 in range(len(system_loops)):
if c1==c2:
row1_branch_list=[]
row1_node_list=[]
for c3 in range(len(node_list)):
if node_list[c3] in system_loops[c1][c2]:
row1_node_list.append(system_loops[c1][c2].index(node_list[c3]))
row1_node_list.sort()
for c3 in range(len(row1_node_list)-1):
seg_vector=system_loops[c1][c2][row1_node_list[c3]:row1_node_list[c3+1]+1]
seg_vector.append("forward")
row1_branch_list.append(seg_vector)
seg_vector=system_loops[c1][c2][row1_node_list[c3+1]:len(system_loops[c1][c2])]
seg_vector.append("forward")
row1_branch_list.append(seg_vector)
system_loops[c1][c2]=row1_branch_list
else:
row2_branch_list=[]
row2_node_list=[]
for c3 in range(len(system_loops[c1][c2])):
for c4 in range(len(node_list)):
if node_list[c4] in system_loops[c1][c2][c3]:
row2_node_list.append(system_loops[c1][c2][c3].index(node_list[c4]))
row2_node_list.sort()
for c4 in range(len(row2_node_list)-1):
seg_vector=system_loops[c1][c2][c3][row2_node_list[c4]:row2_node_list[c4+1]+1]
seg_vector.append(system_loops[c1][c2][c3][-1])
row2_branch_list.append(seg_vector)
row2_node_list=[]
system_loops[c1][c2]=row2_branch_list
view raw branches_break.py hosted with ❤ by GitHub

The next major change is in the ODE solver. Three steps here:

  1. Check for stiff equations - L/R ratio and isolate every stiff branch to a single equation by row operations.
  2. With each row operation on the matrices A, B and E, perform the same row operations on the system_loops. So remove elements in common and add new elements.
  3. Since, the objective has been to isolate a stiff loop and approximate it to a static equation without a d/dt, the interactions with other loops has to be eliminated and these interactions will have to be readjusted.

Now in detail. Here is the code (click on "View Raw" below the code box to see it in a new window):


# Getting rid of elements with small time constants
# that appear multiple times in a circuit.
for c1 in range(matrix_a.rows):
if matrix_a.data[c1][c1]:
if (matrix_e.data[c1][c1]/matrix_a.data[c1][c1])<0.1*dt:
for c2 in range(c1+1, matrix_a.columns):
if matrix_a.data[c1][c2]:
# Check is an off-diagonal element has a
# small L/R ratio.
if (matrix_e.data[c1][c2]/matrix_a.data[c1][c2])<0.1*dt:
# Check if the diagonal element is row c1 has
# smaller L/R than diagonal element in row c2
# The idea is to arrange the equations as per
# time constants. So manipulate the larger time
# constants with respect to the smaller time constants
# and ignore the ones that are below a threshold.
if (matrix_e.data[c1][c1]/matrix_a.data[c1][c1])<(matrix_e.data[c2][c2]/matrix_a.data[c2][c2]):
# If row c2 has larger time constant,
# so remove the common element from row c2.
loop_manipulate(sys_loops, list_of_nodes, c2, c1)
for c3 in range(matrix_a.columns):
matrix_a.data[c2][c3]-=matrix_a.data[c1][c3]
matrix_e.data[c2][c3]-=matrix_e.data[c1][c3]
for c3 in range(matrix_b.columns):
matrix_b.data[c2][c3]-=matrix_b.data[c1][c3]
else:
loop_manipulate(sys_loops, list_of_nodes, c1, c2)
for c3 in range(matrix_a.columns):
matrix_a.data[c1][c3]-=matrix_a.data[c2][c3]
matrix_e.data[c1][c3]-=matrix_e.data[c2][c3]
for c3 in range(matrix_b.columns):
matrix_b.data[c1][c3]-=matrix_b.data[c2][c3]
view raw loop_row_op.py hosted with ❤ by GitHub

Check if a diagonal element has a small L/R ratio (time constant). Now this judgement is purely a guess and I have chosen it as 0.1*SimulationStep. But may have to be changed. If that is the case, this equation has to be isolated. So, the interaction between this loop and other loops, must be eliminated and first we start with the branch that is causing the problems. So look for off diagonal elements in the same row with small time constants.

Suppose an off diagonal element [c1, c2] has a small time constant. This means the branch between loops c1 and c2 has a small time constant. And this therefore means that loop c2 will also have this branch and will also have a small time constant. So, the way to go is to check which loop has a smaller time constant - c1 or c2 by looking at the diagonal elements c1,c1 and c2,c2. The smaller one remains as it is and the larger one is manipulated to remove this stiff branch. This is by changing row c2 to c2-c1.

Look through the entire matrix this way. The objective is then that only the stiff branches will be in their own isolated loops and the stiff elements removed from all other loops.


Next - loop manipulations by changing the elements.


def loop_manipulate(system_loops, list_nodes, row1, row2):
"""Does a row1=row1-row2 for the system loops."""
# A branch segment is that between 2 nodes.
# A loop will contain several branch segments
# in continuity and the first and last nodes will
# be the same. These are diagonal elements in sys_loops.
# The off diagonal elements are interactions between
# loops and may be branch segments without continuity.
# First check in row2 for branches
# not found in row1. Add those changing the direction
# - forward becomes reverse and vice versa.
for c1 in range(len(system_loops)):
# Iterate through branch segments in row2, c1 element.
for c2 in range(len(system_loops[row2][c1])):
# Default flag is no. A single occurance changes the flag.
branch_found="no"
# Go through every branch segment in row1, c1 element.
# Comparing branches in row1, row2 and same column c1
for c3 in range(len(system_loops[row1][c1])):
if system_loops[row2][c1][c2]==system_loops[row1][c1][c3]:
branch_found="yes"
# If not found, add the branch segment
if branch_found=="no":
# Empty segment
br_vector=[]
# The last element is "forward" or "reverse"
for c3 in range(len(system_loops[row2][c1][c2])-1):
br_vector.append(system_loops[row2][c1][c2][c3])
# Check for direction, and reverse it.
if system_loops[row2][c1][c2][-1]=="forward":
br_vector.append("reverse")
else:
br_vector.append("forward")
system_loops[row1][c1].append(br_vector)
# Deletion in row 1 is done after addition.
# The addition of elements has been done with
# direction reversed and so they can't be found.
for c2 in range(len(system_loops[row2][c1])):
# Always go in opposite direction while
# deleting elements.
# Iterate through every branch segment
# of row2, c1
for c3 in range(len(system_loops[row1][c1])-1, -1, -1):
# Iterate through every branch segment of row1, c1
# and delete as soon as duplicate is found.
if system_loops[row1][c1][c3]==system_loops[row2][c1][c2]:
del system_loops[row1][c1][c3]
return
view raw loop_elem_row.py hosted with ❤ by GitHub

The reason for doing this is that when the branch currents are measured, they are measured with respect to the branches present in the loops. So with the above row operations, the loops have changed, so the actual loop descriptions must also change. So when we do c2=c2-c1, we need to find out which elements in c1 are present in c2 and delete them from c2 and which elements are present in c1 but not in c2 and add them reversing their direction.


Last step - readjust the matrices. Here is the code (click on "View Raw" below the code box to see it in a new window):


# The aim is to isolate the loops that
# have low time constants and are difficult
# to integrate. So these equations will be
# converted to static equations and will
# contain only the branch that is stiff
# and a voltage source with the resitance between
# them.
for c1 in range(matrix_e.rows):
# Check if there is a resistance.
# Should be but check anyway
if matrix_a.data[c1][c1]:
# Check if time constant is smaller than
# simulation time step.
if (matrix_e.data[c1][c1]/matrix_a.data[c1][c1])<0.1*dt:
# Check if this equation is related
# to a source. There should be a
# non zero element in B.
flag_input="no"
for c2 in range(matrix_b.columns):
if matrix_b.data[c1][c2]:
flag_input="yes"
# If no connection to input,
# perform row operation with another
# row that has an input.
while flag_input=="no":
c2=0
for c3 in range(matrix_b.columns):
if matrix_b.data[c2][c3]:
flag_input="yes"
if flag_input=="yes":
for c3 in range(matrix_a.columns):
# Make all off diagonal
# elements of A zero in the stiff
# equation. So there should be only
# one large total resistance in the
# stiff branch.
if not c1==c3:
matrix_a.data[c1][c3]=0.0
matrix_a.data[c1][c3]+=matrix_a.data[c2][c3]
matrix_e.data[c1][c3]+=matrix_e.data[c2][c3]
for c3 in range(matrix_b.columns):
matrix_b.data[c1][c3]+=matrix_b.data[c2][c3]
c2=c2+1
for c1 in range(matrix_a.rows):
if matrix_a.data[c1][c1]:
# Look for the stiff loops.
if (matrix_e.data[c1][c1]/matrix_a.data[c1][c1])<0.1*dt:
for c2 in range(matrix_a.columns):
# Remove all elements in E.
# No d/dt terms.
matrix_e.data[c1][c2]=0.0
# Due to the row operations, to connect
# the stiff branch to an input source,
# Some off diagonal terms may contain
# resitances. These will be part of the
# total resistance between the branch and
# the source.
# So add all the offdiagonal elements
# in the stiff loop row to the diagonal
# element. And then set them to zero to
# remove interactions with other loops.
if not c1==c2:
matrix_a.data[c1][c1]+=matrix_a.data[c1][c2]
matrix_a.data[c1][c2]=0.0
# The stiff loop would originally have been
# interacting with other loops. After the row
# operations, the stiff elements would have
# disappeared from the other loops. But there
# may still be some interaction elements
# remaining in the off diagonal terms.
# My assumption is to isolate the stiff loop.
# But the interaction with other loops which
# will be non zero elements in the same column
# as the stiff diagonal element will need to be
# translated into interation between non-stiff
# loops. In simple terms, if X and Y are related
# to Z, then they are also related to each other.
# Make a list of all the rows that have non
# zero elements in the same column as the stiff
# loop.
interac_stiff=[]
for c2 in range(matrix_a.rows):
if not c2==c1:
if matrix_a.data[c2][c1]:
interac_stiff.append(c2)
# The interaction between the non stiff loops is
# as follows.
# If X is aiding Z and Y is also aiding Z,
# X and Y are aiding each other.
# If X is opposing Z and Y is aiding Z,
# X and Y are opposing each other.
if len(interac_stiff)>1:
for c2 in range(len(interac_stiff)-1):
elem1=matrix_a.data[interac_stiff[c2]][c1]
elem2=matrix_a.data[interac_stiff[c2+1]][c1]
if (elem1>0 and elem2>0):
comm_elem=elem1
if (elem1<0 and elem2<0):
comm_elem=-elem1
if (elem1>0 and elem2<0):
comm_elem=-elem1
if (elem1<0 and elem2>0):
comm_elem=elem1
matrix_a.data[interac_stiff[c2]][interac_stiff[c2+1]]=comm_elem
matrix_a.data[interac_stiff[c2+1]][interac_stiff[c2]]=comm_elem
matrix_a.data[interac_stiff[c2]][interac_stiff[c2]]+=abs(comm_elem)
matrix_a.data[interac_stiff[c2+1]][interac_stiff[c2+1]]+=abs(comm_elem)
view raw matrix_re.py hosted with ❤ by GitHub

This was actually a bit tricky. The first part is the stiff equation must be connected to an input source because interactions with all other loops has been eliminated. To eliminate the interaction with other loops, the off diagonal elements in matrix A are set to zero in the stiff loop.

Then check if any of the elements in the row of matrix B is nonzero. If yes, leave it as it is. If not, add the stiff loop with another nonstiff loop that is connected to the input. This will cause one of the elements of B matrix in the stiff row as non zero. At the same time, some of the off diagonal elements in the stiff row will now be nonzero as these are the resistances in the nonstiff loop with which the addition occurred. These resistances are essentially the resistances between the source and the stiff branch. So add all those to the diagonal element in A and then get rid of them. So now the diagonal element of A contains the entire resistance as a loop between the source and the stiff branch and is also linked to the source. So the stiff dynamical equation has been approximated to a static isolated equation.

Next, remove all the elements in matrix E - eliminating all d/dt terms and making it a static equation. The last part is to account for the interaction elements that still linger between the stiff loops and non stiff loops. So if Z is the stiff loop and X and Y are nonstiff loops that interact with Z in that they have nonzero elements in the same column as Z, this means that X and Y will interact with each other. This interaction needs to be captured as we want to isolated loop Z. So look for the elements and move them to the offdiagonal entries between X and Y with the appropriate sign.
Posted by at No comments:

Tuesday, March 5, 2013

Circuit Solver - VII

I think the problem has been solved. Updated the latest code on my Source Forge site:
https://sourceforge.net/projects/pythonpowerelec/files/

Need to close my eyes now. So details tomorrow.
Posted by at No comments:

Friday, March 1, 2013

Circuit Solver - VI

Last post had a mistake in the logic.

By performing row operations, the solution set does not change. Also, the way the ammeter currents are calculated is dependent on the loops. And the loops are described in system_loops and have not changed.

So essentially, the current through the large resistor will be a sum of loops i1, i3 and i4. And the accuracy of this will depend on the simulation time step. Which is why it is wrong. Currently, the current through the large resistor has a peak of 0.04Amps whereas it should be in microAmps.

So now what? This means any change I make to the resistors has to be reflected back to the system_loops. And finally, when the component currents are calculated they have to know that the system has been approximated. This is going to take some time.
Posted by at No comments:
Subscribe to: Posts (Atom)