Efficient methods for dealing with nonlinear elements in MNA solvers

[Hemmar]

I'm following along the tutorial here (DAFx17 Tutorial 3: Dave Berners - Modeling Circuits with Nonlinearities in Discrete Time): https://www.youtube.com/watch?v=7Npx0eaSxow&t=2455s and have found it very instructive. (specific timestamp for this method is 40:55)

Dave describes a technique for avoiding needing to invert the full admittance matrix with an initial example given for a resistor whose value might change in time, though the idea is developed further to apply to the linearised contribution of a diode (which I may have further questions on later). Essentially the (unknown) current is placed on the RHS (source vector) then determined from voltage differences of the solution (I think).

Screenshot 2022-12-15 at 14.27.40.png

I have a couple of questions:

1. to solve for voltages, surely I need to put something for I_R into my source vector in order to solve my linear system? But I can't work out I_R until I solve for node voltages (can I?) so I_R remains unknown until after the linear system solve? I can see how I could do this analytically (see below) but I can't always expect analytic expressions (e.g. for a very complex system)
2. are the a_n, a_m, b_n, b_m trivially found from the remaining MNA A matrix numerically/by inspection and if so how? I can't quite connect the dots here!

My attempt at using this method, for a voltage divider with R1/R2 (with R1 variable) is shown below. I've done this symbolically, but am struggling to see how to do it for arbitrary (more complex) examples. Basically I'd like to know how to do this "programatically".

Screenshot 2022-12-15 at 14.22.25.png

[mystran]

You don't want to invert anything.

You want to solve Ax=b and it's faster to directly solve Ax=b than to invert A. I'd argue it's easier too.

The weapon of choice is usually LU, which also has the nice property that you can do a partial decomposition (ie. essentially block LU step) and only afterwards add additional terms to the not-yet-decomposed bottom-right block, which allows you to iterate a smaller subsystem simply by virtue of ordering all the iterated nodes bottom-right. Entirely mechanical stuff, no symbolic computation necessary at all.

The initial values of the iterated values are called "initial guess" and really any sensible value will usually do, although a better initial guess will converge faster with less iterations. In theory Newton is not guaranteed to converge unless the initial guess is close enough, but with well designed component models "close enough" should ideally cover the whole operating range, because otherwise you risk random runtime blowups. For transient analysis (ie. simulation over time) choosing the final values from the previous timestep as the "initial guess" for the current timestep is a reasonable choice that only gets better when you probably end up oversampling at least a bit to avoid aliasing. For the first time-step, you can usually compute values for zero voltage or something.

My tiny example circuit solver doesn't perform the sub-block optimization thing (it just brute-forces the whole system; node sorting and partial decomposition might add some 100-200 lines and when I wrote it I wanted to keep it to the basics only) and realistically LU itself can be made much faster with abstract interpretation and code-generation (so compiler can optimize the code for a particular circuit directly), but ... this simplified solver illustrates how a "full" solver can be built entirely mechanically simply by stamping terms into a matrix and solving with a standard numerical algorithm (ie. LU).

There are few other other sensible variations possible... but inverting matrices is not sensible except as purely mathematical notation for illustrating that we solve Ax=b for x.

[Hemmar]

Thanks for the response! I've actually been pouring over halite for the past few weeks and found it incredibly instructive, just the right level of abstraction to learn from but still be functional/performant so many thanks for that!

I realise my wording was a bit clumsy - I should clarify that I'm not explicitly inverting but using LU to solve. And from the halite thread (and your comment above) I've a good understanding of how you can structure the problem to minimise the impact of "re-factorising" if some entries change (nonlinear pieces, user parameters). So in a sense I have at least one good solution to this problem "banked". At this stage I'm just curious about this alternative approach.

As I understood it, Dave's method offers an alternative approach whereby arbitrary conductances can be moved to the RHS and solved for per timestep with relatively little cost (inverting 2x2 subproblem where as and bs are known at compile time). You still solve the problem with remaining conductances sensibly (LU decomposition of what he calls Y once at the start). What's not clear to me is if this subproblem can be easily constructed by looking at the admittance matrix (good) or if a symbolic solution to the whole MNA problem is required (bad, impractical beyond a certain size). For simplicity of discussion I'm happy to assume we're only solving circuits with resistors and sources (as once that is cracked, the rest, caps diodes follows naturally from Norton equivalents etc) - so we don't need to think about Newton etc yet.

EDIT: For a bit more context, I'm interested in transient analysis. I'm working on a python script (with heavy reference to Halite) that generates C code to solve small circuits symbolically. This works well but as expected, doesn't generalise to larger circuits. I also have a brute force solver, that solves numerically by LU, which inefficiently refactorises the whole thing if the matrix changes, and am learning about various strategies to improve upon the brute force.

[mystran]

What's not clear to me is if this subproblem can be easily constructed by looking at the admittance matrix (good) or if a symbolic solution to the whole MNA problem is required (bad, impractical beyond a certain size).

There are probably several approaches, but probably the easiest way is to construct it as a by-product of how LU works.

Basically if you look at https://en.wikipedia.org/wiki/Block_LU_decomposition you'll see that the bottom-right block has the form D-CA^-1B after the block LU step. What this means is that the decomposition process (block or not) doesn't care about the contents of D at all until you get far enough to start decomposing that sub-block. The practical consequence then is that you can partially decompose before you even add the D terms into the partially decomposed matrix.

So.. all you need to do is sort your nodes in such a way that all those terms you want to iterate end up in the D block. Then you decompose the rest of it and D-CA^-1B is the sub-system you iterate, where every iteration can substitute different D. Partial forward-substitution of the RHS works too.

It actually makes sense to do this twice: order all the resistors for nodes you don't care about at runtime first, then nodes like capacitors or variable resistors that need per-sample updates in either the matrix or the RHS, then the iterated stuff last. Then do a partial decomposition of the "don't care at runtime" nodes and you can throw them away completely before you even do code-generation.

As for generating fast code... basically give each non-zero element (you want to track zero elements anyway to avoid pivoting on them) of the matrix a separate variable name (so that there's no explicit array in generated code; compiler will optimize better), then run LU as "abstract interpretation" where rather than actually computing stuff, you print out C code that would perform the computation. In a sense, you're basically unrolling the whole thing so that there are no branches whatsoever (well, except for the inner iteration loop) and compiler can perform constant folding and other optimizations like register alloc as if it was straight line code... but your code generator doesn't need to understand what any of the elements mean, it just needs to know what goes into which phase (ie. precompute vs. per-sample vs. iteration).

[Hemmar]

OK thanks - I'm probably flogging a dead horse trying to work out the specifics of his method (if I do work it out I'll update this thread). You've given me plenty of good advice here, so I'll focus on implementing this instead. I'll try and post some example code back here if/when I get it working. Thanks again for the advice.

[mystran]

The one thing I want to add is that you absolutely can solve these things symbolically if you want to gain insight into what is going on, but I like to advocate a purely numerical approach simply because (1) it's easier to keep track of numerical values than long and unwieldy symbolic expressions and (2) this is just as efficient anyway.

[mystran]

After watching the video a couple of times, I honestly can't see what the exact idea is... because I just don't see how you can keep the quadratic convergence of Newton if you only include the diode current in the circuit equation and if you lose the quadratic convergence then you almost certainly lose performance 'cos even an "expensive" direct solution (ie. LU) to a low-dimensional system (ie. usually one dimension per pn-junction) is pretty fast compared to evaluating the non-linear functions. The quadratic convergence of Newton is almost never worth giving up, unless you're aiming for a strategy that doesn't iterate at all.

The strategy where you iterate without changing the admittance matrix (except possibly once per timestep) makes sense if solving Ax=b is very expensive compared to evaluating the non-linearities multiple times.. but that's just not typically true for relatively small circuits that we would hope to simulate in realtime.

I don't know... I don't think the video presentation is all that clear about exactly what is suggested and what the performance characteristics of an actual implementation actually are.

[Hemmar]

Thanks again for your advice - I tried this on a couple of problems with some success and had a follow up question. First I implemented this on a small Schmitt trigger circuit (falstad example) and had no issues:

schmitt_maths.png

However I'm now trying to think about how to solve this more generically. The second example is a diode ladder filter. I can automatically apply row/column swaps (remembering to apply to x and b as appropriate). The code snippet below shows the starting matrix M (with 'N' marking the nonlinear elements, i.e. those that will change each iteration) and a sample reorganisation and partitioning into a block matrix form M -> [[A, B], [C, D]]. The issue is that while M is full rank, the submatrices A and D aren't necessarily, and if they aren't then solution will fail. So my question(s) are:

* is there an intelligent algorithm to swap rows and cols and choose the submatrix split points so that the two submatrices A and D remain invertable / are full rank, whilst retaining this desirable property of partitioning the change and non-changing values. I appear to have chanced upon it for the simpler problem but cannot for the harder
* am I thinking about this the right way? I'm primarily working with symbolic solutions at the moment FYI

Note in this example the first col of D is zero so no surprise it fails, but even with tinkering to have filled rows/cols I still end up with D.det() == 0 / no solution etc.

Code: Select all

[XX                X     ]    [NN     |     X   X       ]
[XXX                     ]    [NNN    |                 ]
[ XX                     ]    [ NNN   |                 ]
[   XX      X       X    ]    [  NNN  |                 ]
[   XXX                  ]    [   NNN |                 ]
[    XNX      N          ]    [    NNN|                 ]
[     XX             X   ]    [     NN|X   X            ]
[       XX               ]    [-------------------------]
[       XXXXX            ]    [       | X              X]
[        XX           X  ]    [       | XXXX           X]
[        X X           X ]    [       | XX           X  ]
[   X    X  XX           ]    [       | X X           X ]
[           XN    N     X]    [      X| X  X   X        ]
[     N       NN         ]    [X      |     X       X   ]
[             NNN        ]    [       |      XX  X      ]
[              NNN       ]    [       |      XX         ]
[               NNN      ]    [       |    X   XX  X    ]
[            N   NN      ]    [X      |        XX       ]
[X                       ]    [       |      X   XX     ]
[  XX                    ]    [       |          X      ]
[    X  X                ]    [       |       XX        ]
[         X              ]    [       |         X      X]
[          X             ]    [       |  X              ]
[           X            ]    [       |   X             ]
                              [       |    X            ]

[mystran]

* is there an intelligent algorithm to swap rows and cols and choose the submatrix split points so that the two submatrices A and D remain invertable / are full rank, whilst retaining this desirable property of partitioning the change and non-changing values. I appear to have chanced upon it for the simpler problem but cannot for the harder

I'm not sure if there is a good algorithm (other than reduce LU row-by-row with partial pivoting and treat the remaining block as inner loop when you can't find a pivot from a linear row; whether or not naive partitioning works without fancy pivoting though might depend on how you design the component stamps), but as far as I can see in your simple example you need A-B*D.inv*C and D be invertible, but A doesn't necessarily need to be.

That said, be careful with idealized components as sometimes these can lead to circuits (depending on how they are connected) where multiple variables are linearly dependent no matter what. When this happens adding some small parasitic series resistances is usually the easiest thing to do.

[Z1202]

That said, be careful with idealized components as sometimes these can lead to circuits (depending on how they are connected) where multiple variables are linearly dependent no matter what. When this happens adding some small parasitic series resistances is usually the easiest thing to do.

Would you mind giving a simple illustrating example, if possible. I'm not sure what exactly you are referring to, however it reminds me of a discussion with a colleague who claimed that you almost always have redundant equations in linear circuit analysis. Turned out he prefers to use voltage differences . OTOH, I at least don't remember such problems when dealing with voltage potentials directly. But maybe I don't remember

[mystran]

That said, be careful with idealized components as sometimes these can lead to circuits (depending on how they are connected) where multiple variables are linearly dependent no matter what. When this happens adding some small parasitic series resistances is usually the easiest thing to do.

Would you mind giving a simple illustrating example, if possible. I'm not sure what exactly you are referring to, however it reminds me of a discussion with a colleague who claimed that you almost always have redundant equations in linear circuit analysis.

I don't think there is a general problem with redundant equations as such, though whether you actually need to solve all variables is another thing. I think any issues with inability to solve a circuit are almost always situations where for one reason or another, one accidentally violates model assumptions. In a sense, all cases I can think of are "degenerate" but the concern I'm trying to raise is that sometimes one has to be careful about what exactly is "degenerate" as far as idealized components are concerned, because this can happen by accident (well, at least by empirical experience).

Probably the canonical example of an ill-defined idealized circuit would be floating capacitors. Let's take a trivial circuit with two idealized capacitors and try so solve V1 (pardon my ascii art):

GND ---]C1[--- V1 ---]C2[--- GND.

There is no DC solution for V1 here, because it's floating. The voltages over C1 and C2 must be same magnitude and opposite polarity, but unless we fix one of them, any choice of V1 can satisfy the circuit laws and transient analysis from initial conditions where voltages over C1 and C2 are different would have to solve for infinite current (well, sort of; it's a Dirac spike, so we might or might not be able to numerically integrate, but trapezoidal method in particular is not L-stable). If we do something like extend our capacitor model with a tiny bit of parallel leakage conductance, then these issues mostly just go away: now the DC solution is V1=0 and transient analysis is stable as long as we leak more than we accumulate rounding errors.

The point I'm trying to make here is that by assuming idealized capacitors with no leakage and no wire resistance, we've built a "circuit" that is unphysical. We can fix this by modifying the circuit (merge the two capacitors into one) such that the assumption for the idealized components hold, or we can fix this by using components that are less ideal... but if something like this happens in a more complex circuit by accident, it's not always obvious exactly where the problem is, rather you just have a matrix with no solution.

[Z1202]

It seems to me you're talking of a slightly more elaborate case of connecting capacitors to voltage sources or connecting inductors to current sources, which are clear violations of the idealized model, of course, but I wonder, whether there are any real circuits which would want to do that? So does this occur in analysis of actual real circuits?

[Hemmar]

Would you mind giving a simple illustrating example, if possible.
 

Definitely, it's hard to talk about in the abstract. Few comments first:

• I work with sympy to symbolically solve the Mx = b MNA problem
• I then autogenerate C++ code (and in my case use in VCV Rack plugins)
• The diode ladder code works as I expect, but I solve for everything every Newton step which isn't necessary
• I'm aware of other optimisations (simd for the exponentials), here I'm thinking about more high level stuff
• not shown is the actual linearisation of the diode/newton step, but that is straightforward

Link to a stripped down jupyter notebook with the details / motivation / what I've tried so far etc:
https://gist.github.com/hemmer/e24957a1 ... d6130ac0ba

The reason I'm surprised it doesn't neatly partition, is that the diode ladder and the preceding amplification appear to be minimally coupled subcircuits. Is there a straightforward way to think about subcircuits in MNA? Can I treat abstract away the diode ladder and treat the whole thing as a two input ports, one output port black box that I solve seperately?

[Hemmar]

That said, be careful with idealized components as sometimes these can lead to circuits (depending on how they are connected) where multiple variables are linearly dependent no matter what. When this happens adding some small parasitic series resistances is usually the easiest thing to do.

Interesting I'll watch out for that, though in this case the full matrix problem can be solved without issues so I don't think it's that. I also tried using unique symbols for resistors R1, R2 rather than R100k R33k in case these somehow exactly cancelled / lead to linear dependencies, but same result.

[mystran]

It seems to me you're talking of a slightly more elaborate case of connecting capacitors to voltage sources or connecting inductors to current sources, which are clear violations of the idealized model, of course, but I wonder, whether there are any real circuits which would want to do that? So does this occur in analysis of actual real circuits?

Why is this a clear violation? Take any electronics textbook and they ask you to solve capacitor mazes. It seems reasonable to expect you could also simulate them.

That said, I really wish I could give a simple, practical example that was slightly less obvious (and ideally without capacitors), but sadly it's been quite a while so I would have to search for one. If I recall correctly though, there is nothing special about integrators in particular, but rather a general problem when you introduce dependent sources (eg. ideal opamps, simple Ebers-Moll transistors, etc), where if you connect several of them together with no additional constraints then you might get linear dependencies unless you follow additional rules of what can be connected where. If these violations are "clear" then perhaps this is a non-issue.