Releasing version 0.2.5

Releasing the version 0.2.5:
http://sourceforge.net/projects/pythonpowerelec/

Next update

A loooong time since I posted. Made a huge number of changes to the simulator. It is now working for a circuit like a three-phase diode rectified feeding a dc capacitor that is in turn the dc bus for a three-phase 6 switch inverter that has an LC filter. So I feel most of the major bugs have been ironed out.

Anyway, I will be releasing the next version of the simulator tomorrow. But before I end the day, I would like to ask a question to some of the readers of the blog.

Since, documenting this project has been a major pain - this blog is the only thing I have and I am quite glad I did at least that. But this blog too is a mess and I need something a little more systematic. So I was thinking of writing a book. Does anyone know of any book on circuit simulation that actually describes how a simulator can be built? Or if anyone would be interested in reading such a book, what would your expectations be from such a book and why? Out of curiosity, because you are trying something similar, or because you are not happy with existing software?

Would be great to hear any feedback. Please send an email to pythonpowerelectronics@gmail.com.

Voltage measurement

Have been working furiously on modifying the simulator for larger circuits. But found there is a fatal error in voltage measurement. Take a look at the test circuit for a three-phase diode rectifier.

The voltmeter is modelled as a large resistance which does not disturb the current in the remaining circuit. The problem is that the branch containing the voltmeter is a stiff branch that appears in the loop equations. These equations appear randomly. So in the above circuit, when the simulation begins, suppose diodes D1 and D6 conduct. This means voltage vac appears across the capacitor Cdc and the voltmeter Vo. What would need to happen is that the voltmeter loop should be written with Cdc and therefore as Cdc charges, Vo output voltage should increase. But what is Vo appears in the loop with vac, D1 and D6? This is perfectly possible the way the program works now. This means Vo will jump to the voltage vac even though Cdc is charging gradually.

So Vo does not reflect the voltage that appears across the nodes but rather on the voltage impressed on it through the loop it appears on. And this inevitably happens whenever the circuit becomes large and loop equations become more and more random. Also, if the voltmeter were connected across a largely inductive circuit or for that matter over two fairly distant nodes in the circuit. What then?

The only error free way to solve this seems to be to do a nodal analysis with only the voltmeters for the nodes where they exist. This calls for a fairly large additional piece of code that will run at simulation time step.

Revised nodal analysis

When must nodal analysis take place?

1. When there are inductors with sufficient energy and that energy might be abruptly interrupted.
2. When does an inductor have sufficient energy? When it has been in at least one non stiff loop in the previous iteration. This way, the circuit is power independent.

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

The process is as follows:

1. In the beginning, a new branch event is assumed by default whenever a branch event occurs.
2. Look if there are any inductors that have previously been in non stiff loops.
3. If so, nodal analysis is required.
4. Following nodal analysis, check if there is a continuity event with branch currents from nodal analysis being different from the branch currents from previous iteration of loop analysis.
5. If so, execute the determine_state function to check if any of the elements will change their state. If any of the elements do change state, the respective branches will have branch events.
6. Check if there are any branch events that had not occurred in the previous iteration. If so, go back to step 2.
7. This will repeat until no new branch events occur.

In nodal analysis, there has been one significant change. In order to perform nodal analysis, it is essential to calculate the node voltages. So the entire circuit has to be solved. However, the branch currents will be calculated only if:

1. At a node, check if there has been an branch which has a "hard" event.
2. A hard event is when a switch turns on or when it turns off while carrying a non negligible current. When a switch or a diode turns off when the current goes negative, that is a soft event. When a diode turns on when it is forward biased, it is a soft event.
3. If a single branch incident at a node has a "hard" event, the currents of all the branches incident at the node will have to calculated. If not, the branch currents remain at their initial values - either zero or a current it is an inductor.

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

There has been another tag created called function_purpose. When determining state, function purpose will only calculate branch currents at nodes where there has been a "hard" event. But while calculating currents before the next loop analysis, all branch currents will be recalculated.

Determining switch and diode state

The code for determining switch state is below (click on "view raw" below the code box to view the code in a new window):



The code for the determining diode state is below (click on "view raw" below the code box to view the code in a new window):


This function is called after performing nodal analysis to check if the inductor current has caused an non linear device to change state. The difference between the code for switch and diode is that a switch can only turn off during nodal analysis while a diode can both turn on and turn off. My guess is that the switch should turn on only in the update_value function. If for some reason, the switch were to turn off during the nodal analysis because the current through it was negative due to an inductor, the switch remains off and only a diode can turn on to freewheel the current.

The reason for doing this will be the next detailed blog post.

Stiff loop evaluation

Been a long time since I have posted. One of the reasons has been that I have made change after change in the solver. With every change, the circuit that I am working on gets fixed but a previously tested one breaks. Even now I am not totally sure until I test them all. But I will start posting code anyway.

To begin with, the stiff loop evaluation. The concept until now has been:

1. Sort out the loops into stiff loops and non-stiff loops.
2. With stiff loops, there are those loops that were non-stiff until the previous iteration and became stiff. So these will have a non-negligible current because a diode may have turned off causing the current to be slightly negative.
3. So the currents of these stiff loops will have to be re-evaluated to brought back to values that correspond to their being stiff sloops with large resistances. A failure to do so will cause these loops to disrupt the calculation of the other stiff loops.
4. Simple way to do this is to have another function. The code is below.

This is the function that calls the stiff equation solver which has not changed from before (click "view raw" link below the code box to see code in a new window):


There is a bit of guessing here as well. Out of "n" loops, let us suppose "m" loops (m<n) have turned stiff in the previous iteration. The remove_stiffness function makes the stiff loops into a upper triangular matrix with the first stiff branch of every stiff loop being present only in that loop. However, as the number of stiff branches can be larger than the number of stiff loops, the last few stiff branches may appear in multiple loops. As an example, take a look at these loops:

stiff_forward no no no no no no stiff_reverse reverse reverse no no reverse no reverse reverse no reverse
no stiff_forward no no no no no stiff_forward forward forward no no forward no forward forward no forward
no no stiff_forward no no no no stiff_reverse no reverse no no no reverse reverse reverse reverse no
no no no stiff_reverse no no no stiff_reverse no reverse no no no reverse reverse reverse reverse no
no no no no stiff_reverse no no stiff_reverse no no forward forward no no reverse reverse no reverse
no no no no no stiff_reverse no stiff_forward no no reverse reverse no no forward forward no forward
no no no no no no stiff_reverse stiff_reverse no no no no no no no no no no

All the loops are stiff and the 8th branch appears in all the stiff loops. And for all you know, this 8th branch could be the branch that became recently stiff. So there is a possibility that these newly formed loops could disrupt the calculation of other new formed loops. So thus there is the initialization of newly formed stiff loop currents to zero.

The current of a stiff loop is decided only by the voltage in that loop and the resistance of that loop. The loop has no dynamics di/dt. So repeatedly calculating the loops will not change the loop currents unless the currents of one of the other stiff loops has changed. Which is why if the newly formed stiff loops are recalculated a number of times (thrice for now), the values should stabilize. In any case without a di/dt, they should not blow up out of bounds.

Following this, when it comes to calculating loop currents from branches, there were two parts to that function. The code is below (click "view raw" link below the code box to see code in a new window):



As can be seen, the entire block of code that would compute currents of stiff loops has been commented out. The reason is that if the currents of stiff loops are being calculated from the stiff equation solver, there is no need to repeat the process. On the contrary repeating the process has the capacity to disrupt the loop current values.

Restarting

Been a long time since I posted. Have been busy with building my hardware setup for my research and now they are almost working the way I want. With this done, I now will have to divert some of my time towards this project.

Spent the past week or so reading most of the code all over again. Another task has been portability to Python 3. Have to look into compatibility issues with respect to that. But before that to complete a tested power electronics library.

One of the problems with combining mesh analysis and nodal analysis was that if both are performed whenever an event occurs, freewheeling can be said to occur even when it is not meant to occur. For example, in a three-phase diode bridge rectifier fed by a source with a finite inductance, the turning off of a diode can be seen as cause for freewheeling of the inductor current and causing the other diode in the leg to freewheel.

So, my guess is there comes the need to distinguish between nonlinear element events -  a hard switched event or a soft switched event. A soft switched event is when a device (diode or switch) turns off when the current through it becomes negative. In that case, freewheeling should not happen. However, when a switch is turned off, and the current through it (irrespective of the magnitude) is broken, freewheeling needs to occur if there was an inductor in a branch connected to the node.

That is the immediate next step.

Hardware setup

Been a while since I posted about the simulator. Have been crazy busy with my day job but am happy to announce, that my microgrid setup is finally ready!

Inverter and control board with sensors and microcontroller:

Intelligent power module with dc bus capacitors:

Hopefully I will be able to work on the simulator again before my trip to Japan.

Initializing new stiff equations

Another problem that arises when there are a large number of nonlinear devices is that nonstiff loops become stiff loops but their current has still not become negligible because the diode has turned off at the previous iteration. For example, current through a diode is -0.1A causing the diode to turn off and it's resistance to change to megaohms. When this stiff loops interacts with other stiff loops and the currents of those stiff loops are calculated, all other stiff loops end up having currents of close magnitudes of 0.1 A. Result: entire simulation is disrupted.

Two things need to be done:
1. Need to identify which loops have just become stiff.
2. Calculate those stiff loops before running the ode solver.

The code to identify the newly formed stiff loops is:


And running the stiff loop ode:



One other thing still needs to be done. What if the newly formed stiff loops themselves are interacting with each other through a stiff branch? Now even this ode solver would fail because they would disrupt each other. The only thing to do would be do perform an upper triangularization of the newly formed stiff loops with respect to the newly formed stiff branches and solve them backwards. Will do this later.

Nodal analysis

There might have been a mistake in the previous current initialization of branch currents. The branches are current sources if:
1. They have inductors and are not stiff.
2. They have voltages and have zero impedance.

If the branches are stiff or are purely resistive with or without a voltage, their currents will need to be calculated by nodal analysis. The code to set the branch currents is:



The actual nodal analysis is:



Other than that, the code is quite the same as the previous case. This code was tested with a three phase diode bridge rectifier and was found to be working OK. The previously released code was breaking.

Rewriting the loop finder algorithm

Previously, the system loops were added if they met the following criteria:

1. Did not pass through the same branch or same node more than once.
2. Last node was the same as the first node.

During the loop finder, no check was imposed on whether the loops were independent and so they were several loops which were independent.The way to deal with this was to let the simulation run and after the first run, the dependent loops would become void and be deleted. A simpler method could be used to find the independent loops:

1. Did not pass through the same branch or same node more than once.
2. Last node was the same as the first node.
3. The system branches must be contained in at least one loop.
4. As the loops are added, a loop will be added if and only if it passes through a new branch that the previously found loops have not passed through yet.

The 4th check simplifies the process and at the first iteration, the loops found are independent loops and equal to the number (Branches-Nodes+1).

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



 The rest of the loop finder code is pretty much the same. The only difference is that the parallel branches are found in the main block above so they have been eliminated from the other functions. Anyway, the rest of the function code is below (click on "view raw" to see it in a new window):

After a long break

Been a while since I updated this blog. I have been working on the project off and on but I was busy doing all sorts of other things that made documentation difficult if not impossible. I have had to extend my work permit and that means getting together all my documents. I have had to send a paper to a conference in Hiroshima, Japan which meant I had to rush to get the final paper ready. And finally, had to plan my trip to Japan which will be in the month of May and get my visa ready.

There have been some questions about how the simulator should work for nonlinear circuits and I am going to start thinking out aloud on the next few blog posts.

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

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



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


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

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



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

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



After this change, the code now works for the circuit shown above.