Making the simulator faster

I was almost ready releasing the next version of the circuit simulator. The only drawback with this version has been that it is a bit slow. Reasons for this are the repeated loop manipulations that are being performed - restricting the stiff loops, computing the loop currents every time there is a switch turning on or off etc.

The past couple of days I tried a nodal analysis approach in comparison with the previous loop analysis approach. I didn't find much improvement in the speed of simulation. However, the matrices are much more unstable and the simulation time step needs to be much smaller. There probably are improvements I could have made to reduce this instability, but I would rather go with the loop analysis method.

To speed up the simulations with loop analysis, the only method I can think of is to maintain a logbook of all the changes taking place in the circuit and how this impacts the loops. This is under the assumption that there will be a finite number of changes possible and therefore a lookup table can be used to get rid of all the loop manipulations. So the simulation will be slow in the beginning but will speed up later.

Anyway, I am taking a break for the next four days to hang out with friends who are beginning to wonder if I am still alive. Seasons greetings to those who celebrate Christmas and for those who don't (and that includes me) get drunk and be merry anyway. See you before new year when I hope I will release my next version.

Diode state determination

Diodes and switches are high value resistances when off and low value resistances when on. The way of detecting when the turn on is to measure the voltage across them at every simulation cycle and change their status when the forward voltage exceeds 1 Volt. For a device in off state this might mean a current of 0.5 microAmps.

When running the simulation of a three-phase inverter, I found that the diodes were trying to turn on at every simulation instant even though the switch across them was conduction. On the contrary, the diodes that were across switches that were on and conducting was significantly higher and therefore were turning on.

To describe this look at this switch/diode combination when the switch is on.

The switch is a part of a loop with a current i1. The diode is off and forms a loop with the switch. By mere observation, it is clear that because of its association with i1, i2 will be large enough for the diode to turn on. There is absolutely nothing wrong with the loop, but the result is blatantly wrong.

This leads me to think about why do I need to use loop analysis to determine the current though stiff elements? The loop i2 could be written in several different ways. If it is written with respect to the dc bus voltage, it would result in a reverse current that would keep the diode off.

The same goes for any stiff element. By making an element stiff, it ensures that it will not disrupt the currents in the circuit. So the effect of the diode can be removed by making it a high resistance. But if the current through the diode matters and in this case it does, loop analysis is not an accurate way to determine the current through the diode.

Since nodal analysis is available and is being used to determine freewheeling, why not use it to determine the current through all stiff elements after the loop analysis ODE is solved.

Here is the code:

Excessive loops

I started working on the case of a three-phase inverter feeding an inductor-capacitor filter with the objective for designing a controller to regulate the voltage across the capacitor. In a nutshell simulate a UPS.

I ran into a number of problems which will be separate posts. I am still not sure if these problems have been completely solved.

The first problem. On simulating in open loop i.e giving balanced sine wave modulating signals I expected smooth balanced current and voltage waveforms. Instead all were badly distorted. The reason was that the circuit has 48 branches and 32 nodes but the number of loops calculated was 20. This isn't so bad because the loop finder works in random ways and could end up with extra loops. The assumption was that these extra loops would be linear combinations of the other loops and would be eliminated by loop manipulations. However, on looking at the system matrices in the differential equations I found 18 loops. So one loop was extra and solving this loop messed up the simulation.

So why did this loop linger? The loop manipulations are done only for stiff loops. This is to restrict stiff branches to the minimum number of loops. But if there are excess loops that are nonstiff to begin with, they may never be eliminated. So, to make sure that the loops are exact, on forming the system loops for the first time, row operations are performed to convert them to an upper triangular form. This will ensure that loops that are linear combinations are eliminated at the beginning of the simulation.

Here's the code:

Releasing version 1.0.1

I am releasing version 1.0.1 following a bug I discovered when I tried to simulate four-wire systems with a neutral wire. The problem lies with the solver and the sequence in which it solves equations. The solution has been to ensure that differential equations get solved before static equations rather than vice versa.

http://www.pythonpowerelectronics.com/downloads.html

 

Three phase grid connected inverter

I posted the next case study - a three phase inverter connected to a three phase grid in current control mode. Check out:

http://www.pythonpowerelectronics.com/blog.html

I have not included the complete control code and will not do so from now on as even though I am releasing Python Power Electronics open source, the control code may be protected legally by someone else. If you would like to write your own control code, check out the tutorials in the documentation section of the website.

PPE Website

I haven't posted for the past week because I have been launching the Python Power Electronics Website. Check it out:

http://www.pythonpowerelectronics.com/

It probably will keep changing for a while until I finalize it. Suggestions are of course very welcome and can be sent to pythonpowerelectronics@gmail.com.

I'll continue this blog but the case studies that I was posting will be posted on the website. This will be more of an informal discussion from now onwards.

PLL with notch filter

The previous post was of a PLL that works when the voltages are balanced. When the voltages are unbalanced, it doesn't work as well. This is because the negative sequence component in the voltage appears as a double frequency (120 Hz) component in the "d" and "q" components of the transformed voltages is hard to suppress. The PI controller can remove higher frequency harmonics and noise present but the 120 Hz component generally gets through.

To describe this concept through simulations, I'll make the "b" phase component of voltage to be zero.

The green line is the "b" phase voltage. This is a fairly extreme case of unbalance and so this is chosen as a benchmark to design a notch filter. The PLL output without a notch filter is:

The above figure shows the double frequency creeping into all control signals and in the above case the angular frequency. Additionally, the double frequency can also be seen in then "d" and "q" components of the transformed voltage as below.

Now that the problem has been described, the notch filter has been inserted before the PI controller. Instead of feeding the "q" component of the transformed voltage to the PI controller, it is passed through a notch filter that removes the 120 Hz component and the filtered voltage is fed to the controller.

The voltage "vq" is the same as above, but the output of the notch filter is:

The 120 Hz component has been attenuated significantly. The output of the PLL is now:

With this done, the next step is to simulate a grid connected inverter feeding a controlled current into the grid.

Questions or comments about this? Email to pythonpowerelectronics@gmail.com 

Phase locked loop

In order to move on to simulating a grid connected inverter, the next step is to implement a Phase Locked Loop (PLL). The method used is the synchronous transfer frame reference alignment. In more detail, starting with a frequency that could be close to what could be the normal frequency, using a Proportional Integral (PI) controller, we adjust the frequency until the signal being locked on to is aligned with the reference frame. So with the frequency being adjusted, calculate the d and q axis components in the rotating reference frame, and adjust the frequency until the q component is zero.

The circuit is pretty simple:

There is an unbalanced three phase resistor bank connected across a three-phase voltage source. Well it doesn't really matter, because we just want the voltage source. The parameters of the circuit are:

Run-time parameters:

-------------------------------------------------------------------------------------------------------------------------
>>>
CSV file containing the network layout --> grid_pll

Enter parameters in file grid_pll_params.csv. When done, close the file and press y and enter to continue -> y

**************************************************
Voltmeter is  Vb  located at  22O  with positive polarity towards 21O
Voltmeter is  Va  located at  22L  with positive polarity towards 21L
Voltage Source is  Vcn of 120.000000 V(peak), 60.000000 Hz(frequency), -240.000000 (degrees phase shift) and 0.000000 dc offset  located at  16D  with positive polarity towards 16E
Resistor is  Rc = 0.100000  located at  16G
Resistor is  Rbc = 1000000.000000  located at  11I
Resistor is  Rb = 0.100000  located at  8G
Voltage Source is  Vbn of 120.000000 V(peak), 60.000000 Hz(frequency), -120.000000 (degrees phase shift) and 0.000000 dc offset  located at  8D  with positive polarity towards 8E
Voltmeter is  Vc  located at  19S  with positive polarity towards 18S
Voltage Source is  Van of 120.000000 V(peak), 60.000000 Hz(frequency), 0.000000 (degrees phase shift) and 0.000000 dc offset  located at  1D  with positive polarity towards 1E
Resistor is  Ra = 0.100000  located at  1G
Resistor is  Rab = 1000000.000000  located at  4I
**************************************************

Enter the control files. Omit the .py extension and just leave spaces between files --> pll_3ph
Enter control parameters in the following files -->
pll_3ph_desc.csv
When ready press y and enter to continue -> y
Enter control code in the following files -->
pll_3ph.py
When ready press y and enter to continue -> y

**************************************************
Meters are in the following sequence: ['19S', '22L', '22O']
Control variables to be plotted are in the following sequence
['control_signal1', 'control_signal3', 'control_signal2', 'control_signal4']
**************************************************
Simulation time step in seconds (recommended <= 10.0e-6) --> 10.0e-6
Duration of simulation in seconds --> 1.0
Data storage time step in seconds ( >= 0.000010) --> 10.0e-6
Output data file name (.dat extension will be added) --> ckt_out


Data stored in ckt_out.dat. Check output log file output_log.txt for description of output
35.8410000801
-------------------------------------------------------------------------------------------------------------------------


In this case, there is a control file called pll_3ph even though there is no controllable element. VoltageSource, Resistor and Voltmeter are the only components in the circuit. In this case, the control file does not perform any control but is merely a signal processing code that generated control_signal1, control_signal2, control_signal3, control_signal4.

I will skip the details of the control and move on to the results. The control_signals 1,2,3,4 generated by the control code are
control_signal1 --> volt_d i.e d component of the voltage in the rotating reference frame
control_signal2 --> volt_q i.e d component of the voltage in the rotating reference frame
control_signal3 --> omega i.e angular frequency of rotating reference frame and the system
control_signal4 --> phase_angle of the sinusoidal signals


So first and foremost, frequency:

 In steady state tracks 377 rad/s which is 2*pi*60.
Next, volt_q:

As stated above in the principle of operation of the PLL, in steady state the q component of voltage in the rotating reference frame will be zero. This is when the rotating frame is at the same frequency as that of the source.

Next, volt_d:

 The d component of the voltage settles to 147V. This is the Root Mean Square (RMS) value of the Line to Line (L-L) voltage of the system. There are several forms of transformations that use several multiplying factors that could case d component to have all kinds of values - peak of Line to Line, peak of phase etc.

The last is the phase angle. This needs to be plotted with the a component of the voltage to give you an idea of the alignment.

 The way the component VoltageSource generates voltage is as a sine template with respect to time. As seen above, the phase angle starts at 0 from the peak of what is the "a" phase voltage. This is because the PLL uses cosine and sine templates of the phase angle to align the rotating reference frame. The cosine of the above phase angle will therefore turn out to be a cosine function corresponding to the sine waveform of the voltage.

The next step is including a notch filter to make sure the PLL works even when the grid is unbalanced. After which I can simulate a grid connected inverter.

 Questions or comments? Send an email to pythonpowerelectronics@gmail.com.

Three-phase standalone inverter

So lets start the simulation series with a standalone three-phase inverter - basically the ends of the filter shorted together. This is the screenshot of the .csv file.

The parameters of the inverter are in this screenshot

A few things to note:

1. The voltage source is just 10V because with a short circuit, that's all I need to drive a large current. A dc voltage is specified by making the ac peak 0 and adding a dc offset. The reason is that the peak automatically means an ac signal of frequency specified in the next column.
2. The switches S1 to S6 are shown to have control signals as s1gate to s6gate. However, when the parameter file is initialized, these values are arbitrary and may even be the same. The control signals can't be the same! The problem is not from the simulator point of view but for the control code. Each switch needs a unique control tag. So change as needed.

Before continuing, these are the main inputs the simulator wants from the user:

------------------------------------------------------------------------------------------------------------------------------

$ python circuit_solver.py

CSV file containing the network layout --> standalone_inv

Enter parameters in file standalone_inv_params.csv. When done, close the file and press y and enter to continue -> y 

**************************************************

Resistor is Cdc = 0.001100 located at 5C

Resistor is Ra = 0.100000 located at 26AA

Resistor is Rc = 0.100000 located at 12AA

Voltage Source is Cdc of 0.000000 V(peak), 0.000000 Hz(frequency), 0.000000 (degrees phase shift) and 10.000000 dc offset located at 9C with positive polarity towards 8C

Ammeter is Ac located at 12U with positive polarity towards 12V

Ammeter is Aa located at 26U with positive polarity towards 26V

Inductor is Lc = 0.001000 located at 12X

Inductor is La = 0.001000 located at 26X

Diode is D2 located at 18H with cathode polarity towards 17H

Diode is D4 located at 18L with cathode polarity towards 17L

Diode is D6 located at 18P with cathode polarity towards 17P

Switch is S4 located at 17J with negative polarity towards 18J

Switch is S6 located at 17N with negative polarity towards 18N

Voltmeter is Vdc located at 11A with positive polarity towards 10A

Inductor is Lb = 0.001000 located at 19X

Switch is S1 located at 6F with negative polarity towards 7F

Diode is D1 located at 6H with cathode polarity towards 5H

Switch is S2 located at 17F with negative polarity towards 18F

Switch is S3 located at 6J with negative polarity towards 7J

Diode is D3 located at 6L with cathode polarity towards 5L

Switch is S5 located at 6N with negative polarity towards 7N

Diode is D5 located at 6P with cathode polarity towards 5P

Ammeter is Ab located at 19U with positive polarity towards 19V

Resistor is Rb = 0.100000 located at 19AA

************************************************** 

Enter the control files. Omit the .py extension and just leave spaces between files --> currcont_3ph

Enter control parameters in the following files -->

currcont_3ph_desc.csv

Enter control code in the following files -->

currcont_3ph.py

When ready press y and enter to continue -> y

**************************************************

Meters are in the following sequence: ['12U', '19U', '26U', '11A']

Control variables to be plotted are in the following sequence

['control_signal1', 'control_signal3', 'control_signal2', 'control_signal4']

**************************************************

Simulation time step in seconds (recommended <= 10.0e-6) --> 1.0e-6

Duration of simulation in seconds --> 0.1

Data storage time step in seconds ( >= 0.000001) --> 1.0e-6

Output data file name (.dat extension will be added) --> ckt_out 

Data stored in ckt_out.dat. Check output log file output_log.txt for description of output

258.522240877

------------------------------------------------------------------------------------------------------------------------------

A quick description:

1. Execute circuit_solver.py. Make sure all the simulator files and .csv files are in the same folder.
2. It asks you for the circuit layout which is a csv file. In this case it is standalone_inv.csv but the .csv is dropped in the input as that is added automatically.
3. If the network layout has just been created, the simulator will read all the components in the layout and will insert them into the parameter file (shown above) with default parameters. These default parameters are almost definitely not the ones you would like to use, so open the standalone_inv_params.csv file, make changes, save it and close it.
4. It prints out the parameters for one last confirmation.
5. Next it asks for control files. Now, these need not be control files but can also be signal processing or any processing file like heat/loss calculation etc. So you do not need to have controlled components to have a control file. Add these control files.
6. Again for each control file the simulator creates a default "descriptor" file as a .csv file. More on this later. The default will have one input (usually a meter). If you would like to use another input or have more inputs, change or copy that row and make changes. There will be one output which will be a controlled entity if there is one. There is one time event. And there is a static variable. And there is one variable storage. More on these later.
7. It asks you to make sure control code is ready in the control files.
8. Then outputs the manner in which the output is saved in the data file. Here the first column is always the time so the other columns start from 2 onwards.
9. Then the next columns will be the stored variables. More on these later.
10. Next inputs are the time step (1 microsec), duration of simulation run (0.1 sec), storage time step (1 microsec but can be higher) and the name of the output data file.
11. On completion tells you how long the simulation took to run.

Now about the descriptor file currcont_3ph_desc.csv

1. In this case, the inputs are the three ammeter outputs. The next column are the names of the variables by which these ammeter outputs will be available in the control code. The simulator makes these meter outputs available with these variable names so you can directly use in the control code. Will be shown.
2. The outputs are shown next. The first three columns are from the parameter file. The next two are specific to the control file. The desired variable name will be assigned to the control outputs at the end so the user can assign any value to the control outputs and they will be passed on. And finally an initial output.
3. The next set are the static variables, these are the variables whose values need to be stored from one iteration to the next. Any variable not declared here but used in the control function is a local variable to that control function. It will exist only in that particular iteration and will cease to exist when the control is executed. So for any values that need to be stored between iterations like itegrals in PI controllers etc, use static variables. In general unless sure that you don't need to store a variable, declare it as a static variable.
4. Two time events are generated - t1, and tcarr. t1 is the time step of the control algorithm. tcarr is the time step of generation of the carrier wave and comparison with the modulation to generate the control signals.
5. The last are the control_signals 1 to 4. These are control signals that can be plotted for debugging like modulation signals, integrals, error signals etc. A note - these control signals are common to all control functions in a simulation. So declaring the same control signal for two control functions is a violation. Also, this is a way to share control signals between control functions. Similar to routing signals between blocks.

Finally the control code:



The control code shows the closed loop control using a PI controller with d-q transformation. Take a moment to find out which variables have not been declared static and the reason why. For that matter, all the variables used can be declared static. Python does not need variables to be declared and variables come into existence when they are used for the first time. So for example dt_sample=100.0e-6 is when dt_sample is declared.


A few results. First one is the current waveforms. The next one is the tracking performance. The reference for d is 10 and q is 0. From waveforms, the actual currents are begin to track the references.

Anyway, that's all for now folks. A long post after a long time.

Continuing with the simulator

Been ages since I posted. Lots of changes - new job, moved to another city and a host of others. Programming has been happening slow and steady. As for this blog and the project in general, I am thinking of taking a slightly different course.

So far have been posting code. From now on I will post simulation results. The purpose is that I would like to create a running document with applications rather than just code. This will be of more use to power electronics engineers than mere code based on the simulator.

So keep visiting this page :)