Halite - an analog circuit simulator in ~1k lines of code

[mystran]

(or a whole lot less if you strip out all the comments)

So, I was writing some new circuit simulation code to see how complicated it would be to build a little test-bed for this kind of stuff and somewhere around the time I got the basic solver working, it became clear that it actually didn't take much code at all. Since I've kept telling people that automating MNA is easy to do, I figured I could publish this as an example of just how easy it can actually be.

So, I took the code, cleaned it a bit, stripped out everything that wasn't essential to the simulation (eg. solver optimisations, graphical plotting) and concatenated it all into a single file C++ file. The result can be found here:

https://gist.github.com/signaldust/74ce ... 4a19ce8365

It's tested with OSX clang, but should work with any decent C++ compiler on any platform really.

The circuit is defined in the bottom of the file in main(), so it takes a recompile to change the circuit or the simulation parameters; there's no front-end parser.

The circuit included is a simple 2-transistor amplifier (one for voltage amp, another for buffering output) that appears to work, even though some analog wizard will undoubtedly tell me my circuit is very badly designed. I really just wrote it as a test for the BJT model and tweaked it until I got (close to) the maximum voltage gain (I'd guess it's about 10 or so) that wouldn't clip to supply or ground with the .5Vpp PWM input, so take it for what it's worth.

The actual simulator tries to be somewhat user-friendly with it's output (and wastes quite a few lines of code doing so), with pretty-printing of the MNA matrix (eg. should be easy to see what's going on) and some output formatting (eg. should be easy to import into spreadsheets for plotting).

Calling buildSystem() will initialise the system for DC analysis (ie. "infinite time-step") and calling simulateTick() directly will find the DC operating point. Calling setTimeStep() will setup transient analysis (eg. simulateTick() will then advance by one time-step). If you did the DC analysis, then transient analysis will start from the DC solution, otherwise it'll start with zero-state everywhere.

The supported components are resistors, caps, voltage sources (fixed and callback function), voltage probes, pn-diodes and npn/pnp-transistors (with Ebers-Moll model; see comments for caveats). The code for the implemented components should make it fairly obvious how to add more (or change the existing models).

The capacitors and pn-junctions are solved on their own rows. This is explained in the comments and not really useful with this "dumb" solver, but with models built like this you can separate the dynamic and non-linear parts of the system and only iterate/step a reduced system (something I stripped out from this code to simplify the presentation), which would be a huge win if you wanted to do real-time simulation.

There is very little (read: none) error checking in this thing (beyond printing "nan nan nan" when stuff starts failing), but hopefully the code is easy enough to follow that it serves as a demonstration of how circuit sim with MNA doesn't need to be very complicated at all.

Update: While there's no plotting built-in (because dragging in enough dependencies to make things works makes everything too complicated), here's a picture of what a plot of selected traces from the example circuit could look like (observe how the DC operating point allows the transient analysis to start from a sane initial condition):

halite-example-plot.png

[mystran]

If someone's too lazy to scroll all the way to the bottom of the file, here's the example "main()" code so you can get a hang of exactly what I'm talking about here:

Code: Select all

int main()
{
    /*
        BJT test circuit
        this is kinda good at catching problems :)
    */
    NetList * net = new NetList(8);

    // node 1 is +12V
    net->addComponent(new Voltage(+12, 1, 0));

    // bias for base
    net->addComponent(new Resistor(100e3, 2, 1));
    net->addComponent(new Resistor(11e3, 2, 0));

    // input cap and function
    net->addComponent(new Capacitor(.1e-6, 2, 5));
    net->addComponent(new Resistor(1e3, 5, 6));
    net->addComponent(new Function(fnGen, 6, 0));

    // collector resistance
    net->addComponent(new Resistor(10e3, 1, 3));

    // emitter resistance
    net->addComponent(new Resistor(680, 4, 0));

    // bjt on b:2,c:3,e:4
    net->addComponent(new BJT(2,3,4));

    // output emitter follower (direct rail collector?)
    net->addComponent(new BJT(3, 1, 7));
    net->addComponent(new Resistor(100e3, 7, 0));

    net->addComponent(new Probe(7, 0));

    net->buildSystem();

    // get DC solution (optional)
    net->simulateTick();
    net->setTimeStep(1e-4);
    net->printHeaders();

    for(int i = 0; i < 25; ++i)
    {
        net->simulateTick();
    }
    net->printHeaders();

    net->dumpMatrix();
};

And before you ask, yes .. you could put the rest of the simulator into a header, but I wanted a single-file thingie.

[stratum]

After staring at the code for a few minutes I have found that I need to read http://qucs.sourceforge.net/tech/node14.html

For the time being, I'm trying to understand the differences between nodal analysis and its modified version.

Thanks for sharing.

[mystran]

Oh and some additional observations: if you call dumpMatrix() before any simulateTick() calls, you'll get the matrix in the order it was built. This will probably make it easier to understand what is going on if you're not entirely fluent in linear algebra and find the row-swaps confusing.

If on the other hand you call dumpMatrix() after simulateTick(), it will be printed after the row-swaps from partial-pivoting by LU (ie. in the order it was actually solved) and you will see '#' entries where the original matrix had a zero, but LU filled it with something (where as the dot-entries are zero even after LU). The row-swapped system is obviously fully equivalent, but the MNA stamps can be harder to see as the rows and node-vector entries are in different orders (since the node-vector stays in the original order).

In the past some people were concerned that LU-pivoting would cause problems for real-time implementation, but it turns out that if you modify luFactor to print a message when it swaps rows, you'll notice that this rarely happens after the matrix has been solved once. In fact the order is surprisingly stable considering we will always swap rows here for a "better" pivot no matter how small the improvement. With an additional condition that requires some minimum improvement (eg. new pivot is order of magnitude larger), the order should practically never change once chosen.

[mystran]

For the time being, I'm trying to understand the differences between nodal analysis and its modified version.

In regular nodal analysis, all the "unknown" variables in the circuit are voltages and all the right-hand side values are currents. The problem is, this only works when you can solve the voltages from the currents.

For some components like voltage sources, we can't possibly solve for voltage (it's fixed!) and the unknown is the current instead. For other components like diodes, we could probably do it either way in theory, but solving the current from the voltage works much better in practice (except when it doesn't, see footnote).

So in "modified nodal analysis" we use regular nodal analysis solving for the nodal voltages, but add additional unknowns to the system (eg. for currents) whenever this is necessary or convenient (and my code actually goes a bit further down the "convenient" path than what you would find in text-book descriptions of MNA).

Now strictly speaking we could also use regular nodal analysis and just use gyrators to transform eg. current sources into voltage sources, but since we still would need to add extra nodes and the resulting equivalent circuits end up having stamps that match the MNA stamps, it sort of just makes sense to call the unknowns "current" (or whatever else you want to solve for) and the method "modified" and be done with it.

----
Footnote: So above I noted that with diodes solving i=f(v) doesn't always work. What can happen is that during Newton-iteration the linearised system solves into a voltage that is unrealistically large (and the current would exceed floating-point range). In this case it's necessary to figure out what current resulted in this (non-sense) voltage and use that current to get a more realistic estimate of the voltage. Unfortunately trying to use current as the unknown runs into the opposite problem with low voltages, so really we have to pick one or the other on per-iteration basis either way. So at the end of the day it comes down to the fact using voltage as the unknown just gives us better numerical behaviour.

[mystran]

Add a plot for those wondering what the example circuit is doing when you run the simulation.

[mystran]

If anyone wonders where this is heading, my current optimiser (work-in-progress, we can do better) claims it can solve the system for this particular circuit (not including the actual Newton evaluation functions, just the LU) with 35 multiply-adds and 4 divisions per iteration (including the final one that halite doesn't count), another 9 multiply-adds per time-step to advance the simulation (eg. solve cap-state and v:probe for output) and finally another 41 multiply-adds if we want the full solution (not needed by the simulation itself).

Obviously a lot more work needs to be done when the time-step is changed and these numbers actually grow pretty fast with circuit complexity, but .. just something to think about anyway.

[stratum]

For the time being, I'm trying to understand the differences between nodal analysis and its modified version.

In regular nodal analysis, all the "unknown" variables in the circuit are voltages and all the right-hand side values are currents. The problem is, this only works when you can solve the voltages from the currents.

For some components like voltage sources, we can't possibly solve for voltage (it's fixed!) and the unknown is the current instead. For other components like diodes, we could probably do it either way in theory, but solving the current from the voltage works much better in practice (except when it doesn't, see footnote).

So in "modified nodal analysis" we use regular nodal analysis solving for the nodal voltages, but add additional unknowns to the system (eg. for currents) whenever this is necessary or convenient (and my code actually goes a bit further down the "convenient" path than what you would find in text-book descriptions of MNA).

Now strictly speaking we could also use regular nodal analysis and just use gyrators to transform eg. current sources into voltage sources, but since we still would need to add extra nodes and the resulting equivalent circuits end up having stamps that match the MNA stamps, it sort of just makes sense to call the unknowns "current" (or whatever else you want to solve for) and the method "modified" and be done with it.

----
Footnote: So above I noted that with diodes solving i=f(v) doesn't always work. What can happen is that during Newton-iteration the linearised system solves into a voltage that is unrealistically large (and the current would exceed floating-point range). In this case it's necessary to figure out what current resulted in this (non-sense) voltage and use that current to get a more realistic estimate of the voltage. Unfortunately trying to use current as the unknown runs into the opposite problem with low voltages, so really we have to pick one or the other on per-iteration basis either way. So at the end of the day it comes down to the fact using voltage as the unknown just gives us better numerical behaviour.

I remember that problem from this paper http://www.simulanalog.org/statevariable.pdf

Oddly, it doesn't reference or mention any MNA related paper. Maybe generic realtime circuit simulation wasn't feasible with the CPU's of the day and they didn't even search for it.

[mystran]

I remember that problem from this paper http://www.simulanalog.org/statevariable.pdf

Oddly, it doesn't reference or mention any MNA related paper. Maybe generic realtime circuit simulation wasn't feasible with the CPU's of the day and they didn't even search for it.

No idea. I don't think I ever deciphered a single simulanalog paper successfully.

[Music Engineer]

thank you very much! circuit emulation is one of the many things that i don't know much about but i think i should. i will take a closer look as soon as time permits! ...and hopefully will learn something

[stratum]

This tutorial appears to be good as an introduction.

https://www.swarthmore.edu/NatSci/echee ... /MNA2.html

It explains what MNA is, additionally how we can just 'stamp' the components to a matrix (the 'observations' section). (At first this wasn't clear to me at all)

[mystran]

It explains what MNA is, additionally how we can just 'stamp' the components to a matrix (the 'observations' section). (At first this wasn't clear to me at all)

Yup.

One final observation is that while you need a "special ground" to reference your solution against (otherwise there would be an infinite number of solutions), you can actually do what I did and save some conditionals by treating it as "node zero" as far as the stamping process goes and simply skip over it in the solver (which is why the LU factorisation starts from index 1).

edit: actually another useful observation is that you can relax many of the "text-book rules" quite a bit and just treat the matrix more or less like a generic set of linear equations.. like my capacitor formulation probably looks weird on first sight, but it's really just the standard TDF2 integrator embedded directly into the MNA matrix (because why not)

[Mayae]

Thanks for sharing, looks really cool. Will also take a closer look when time permits.

[mystran]

In the spirit of what I posted in the SEM thread for OTA, here's a basic opamp for Halite using the same diode-clamping strategy to keep things supply limited, so it's safe to use even without feedback (although for comparator applications you might want to add a pole for some time-delay). It probably clips a bit "too nice" when compared to your average real-life opamp (and goes unreasonably close to supply if you push it very hard), but should behave realistically enough to actually work most of the time (except as far as the gain is unreasonable at high frequencies which might cause problems and can be fixed by adding a pole and the output current comes from "thin air" as it doesn't actually draw any current from the supply).

It's a bit spammy in terms of internal nodes (with useless values) so skipping the printing of the internal nodes might be a good idea (which is easy to do by adding a flag into MNANodeInfo but for some reason I didn't include that in the original code).

While this thing has 6 internal nodes (to keep it easy to understand and one of which is purely to calculate the actual output current for printing purposes and could be removed), only two of these need iteration, so it might be useful for real-time situations where you still want supply limits (although I suppose for opamps just using atan() or something would also be an option).

Code: Select all

    struct OPA : Component<5, 6>
    {
        // diode clamps
        JunctionPN  pnPP, pnNN;

        double g;

        // pins: out, in+, in-, supply+, supply-
        OPA(int vOut, int vInP, int vInN, int vPP, int vNN)
        {
            pinLoc[0] = vOut;
            pinLoc[1] = vInP;
            pinLoc[2] = vInN;
            pinLoc[3] = vPP;
            pinLoc[4] = vNN;

            // the DC voltage gain
            g = 10e3;

            // any sort of reasonable diode will do
            double is = 6.734e-15;
            double n = 1.24;
            initJunctionPN(pnPP, is, n);
            initJunctionPN(pnNN, is, n);

            linearizeJunctionPN(pnPP, 0);
            linearizeJunctionPN(pnNN, 0);
        }

        bool newton(MNASystem & m)
        {
            return newtonJunctionPN(pnPP, m.b[nets[5]].lu)
                && newtonJunctionPN(pnNN, m.b[nets[6]].lu);
        }

        void stamp(MNASystem & m)
        {
            // What we want here is a high-gain VCVS where
            // we then bypass to rails if we get too close.
            //
            // Here it's important not to have series resistance
            // for thee diodes, otherwise we can still exceed rails.
            //
            // NOTE: the following ignores supply currents
            //
            //      0   1   2   3  4   5   6   7   8    9
            //    vout in+ in- V+ V- vd+ vd- id+ id- iout
            // 0 |  .   .   .   .  .   .   .  +1  -1   +1 | vout
            // 1 |  .   .   .   .  .   .   .   .   .    . | in+
            // 2 |  .   .   .   .  .   .   .   .   .    . | in-
            // 3 |  .   .   .   .  .   .   .   .   .    . | v+
            // 4 |  .   .   .   .  .   .   .   .   .    . | v-
            // 5 |  .   .   .   .  . gpp   .  -1   .    . | vd+  = i0pp
            // 6 |  .   .   .   .  .   . gnn   .  -1    . | vd-  = i0nn
            // 7 | -1   .   .  +1  .  +1   .  ~0   .    . | id+  = +1.5
            // 8 | +1   .   .   . -1   .  +1   .  ~0    . | id-  = +1.5
            // 9 | -1  +g  -g   .  .   .   .   .   .   ro | iout
            //
            // We then add one useless extra row just to add
            // the currents together, so one can plot the real
            // current that actually get pushed into vOut

            // output currents
            m.stampStatic(+1, nets[0], nets[7], "+1");
            m.stampStatic(-1, nets[0], nets[8], "-1");
            m.stampStatic(+1, nets[0], nets[9], "+1");

            // output feedback
            m.stampStatic(-1, nets[7], nets[0], "-1");
            m.stampStatic(+1, nets[8], nets[0], "+1");
            m.stampStatic(-1, nets[9], nets[0], "-1");

            // voltage input
            m.stampStatic(+g, nets[9], nets[1], "+G");
            m.stampStatic(-g, nets[9], nets[2], "-G");

            // supply voltages
            m.stampStatic(+1, nets[7], nets[3], "+1");
            m.stampStatic(-1, nets[8], nets[4], "-1");

            // voltage drops from the supply, should be slightly
            // more than the drop voltage drop across the diodes
            m.b[nets[7]].g += 1.5;
            m.b[nets[8]].g += 1.5;

            // diode voltages to currents
            m.stampStatic(+1, nets[7], nets[5], "+1");
            m.stampStatic(-1, nets[5], nets[7], "-1");

            m.stampStatic(+1, nets[8], nets[6], "+1");
            m.stampStatic(-1, nets[6], nets[8], "-1");

            double ro = 10;

            // the series resistance for the diode clamps
            // needs to be small to handle the high gain,
            // but still use something just slightly non-zero
            double rs = gMin;

            // series resistances for diodes
            // allow "gmin" just for pivoting?
            m.stampStatic(rs, nets[7], nets[7], "~0");
            m.stampStatic(rs, nets[8], nets[8], "~0");

            // series output resistance
            m.stampStatic(ro, nets[9], nets[9], "ro");

            // Finally (not show above) put a resistor between
            // the input pins.. this could be more as well.
            double ri = 50e6;
            m.stampConductor(1. / ri, nets[1], nets[2], "ri");

            // junctions
            m.A[nets[5]][nets[5]].gdyn.push_back(&pnPP.geq);
            m.A[nets[5]][nets[5]].txt += "gm:V+";
            m.A[nets[6]][nets[6]].gdyn.push_back(&pnNN.geq);
            m.A[nets[6]][nets[6]].txt += "gm:V-";

            m.b[nets[5]].gdyn.push_back(&pnPP.ieq);
            m.b[nets[5]].txt += "i0:V+";
            m.b[nets[6]].gdyn.push_back(&pnNN.ieq);
            m.b[nets[6]].txt += "i0:V-";

            m.nodes[nets[5]].name = "v:opa+";
            m.nodes[nets[6]].name = "v:opa-";
            
            m.nodes[nets[7]].name = "i:opa+";
            m.nodes[nets[7]].type = MNANodeInfo::tCurrent;

            m.nodes[nets[8]].name = "i:opa-";
            m.nodes[nets[8]].type = MNANodeInfo::tCurrent;

            m.nodes[nets[9]].name = "i:opaG";
            m.nodes[nets[9]].type = MNANodeInfo::tCurrent;

            // this is useless as far as simulation goes
            // it's just for getting a nice current value
            m.stampStatic(+1, nets[10], nets[7], "+1");
            m.stampStatic(-1, nets[10], nets[8], "-1");
            m.stampStatic(+1, nets[10], nets[9], "+1");
            m.stampStatic(+1, nets[10], nets[10], "+1");

            m.nodes[nets[10]].name = "i:opa";
            m.nodes[nets[10]].type = MNANodeInfo::tCurrent;
        }
    };

Turning this into the proposed OTA model is left as an exercise for the reader.

[Richard_Synapse]

Spend some time fiddling with your code, gotta say this is dope!

Tried a Fuzz Face circuit just for fun, and it seemed to render fine without any problems

Richard