Version 0.1.2 released.

Released the next minor patch:
http://sourceforge.net/projects/pythonpowerelec/files/

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):



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):




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.

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):



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.

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.

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.

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):



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):



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.



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):



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.

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.

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.