wrestling math and electronics while following "from circuit to code"

[antto]

a bunch of months ago i watched the "from circuit to code" video (from Andy (from Cytomic)) and this is the first thing that began to make some sense to me (i hope you haven't forgotten that i'm not good at math)
so i tried to follow it with something very simple which i understand DSP-wise, a RC lowpass.
i did it by hand and it seemed to work... so i rewatched the video a pile of times (i kinda couldn't follow once the diodes got into the picture) ... then i went to read elsewhere about NA and MNA, but basically if there are too many components (and by that i mean like more than 3) i better NOT do it by hand because i'm slow and i'll mess it up and then i will only waste time and not get anywhere

so instead, i tried to make a tool that automates the hand-work.
i wrote a CLI program that reads a netlist and spews out the node equations... then i formatted the equations so they can be directly (pretty much) pasted into wxMaxima, which outputs voltage equations
and after a lot of butchering, it worked - at least resistors and capacitors..
unfortunately, in the beginning, i made the prediction that this probably won't work at all and i kept it simple and assumed that components have 2 pins.

i had to leave it and a pile of time passed, i kinda forgot about it and in what state it was.
some months later, i "discovered" it again and wanted to ask here but i had forgotten about it so hard that i couldn't even formulate a question.. i had to re-learn the stuff i forgot and see where i had left it
in the process i fixed some things, broke some other things, i really wanted to add a proper potentiometer to it but that required 3 pins.
so i butchered it heavy and now it can do multi-pin components

i went back to the "circuit to code" video and tried hard to figure out the scary diode stuff
after a long time of head<->wall "physical interaction", i think i got the "hpf-lpf-diode-clipper" circuit from the video working!

and now that i could have more than 2 pins, i really really wanted to have transistors, because there is this one pesky circuit involving a PNP transistor that has kept me up for long days and nights...

i went to the qucs pages which have technical explanations of the transistor... yeah... no way i can do this. just the mere list of parameters is 1 page long o_O

i tried other places, there was nothing simple anywhere.
it's 2026, so i gave up and asked the chatbot - it gladly and very happily synthesized a transistor for me!
and of course, it didn't work (or i didn't know how to make it work).
i spent a pile of times trying different things in that direction (asking again, but more clearly, etc..) every time the result was different code and often it looked obviously broken.
i gave up with that.. but at least in the process, i learned some useful keywords and terminology about what i'm trying to do (and what i'm NOT doing).
well, i hated that i couldn't get anything working out of this, and looked in a few other places and saw some familiar chunks and stuff, so i went back to the qucs page, looked at it and i said to myself "i spent so much time elsewhere, this page now looks small, i could try my best to write down the math from it on the side into a text editor, figure out what goes where, and see if i can make something out of it, and if i can omit some stuff easily (since that page talks about the gummel-poon model, which now i know is not bad, but for audio it might be even overkill, it depends)
well, i did that, and the only unclear thing was the actual contributions to the node equations, since that page gives you some table stuff, not equations
well i asked the chatbot about that, it kinda helped with that, the table is a matrix, i wasn't sure about the polarities of the "ieq" variables but i just threw it in and.... it didn't explode immediately...

so maybe i was getting close to something now.
after some more code butchering, i got a very dumb circuit with an NPN to at least not explode and give some values that are not too far from what i get in falstad.
so i tried a PNP circuit and x_x it exploded of course, turns out the "ieq" had to be inverted, but there were other problems too.
couldn't figure out what's wrong.. i went to see how falstad's transistor works, and hah - saw some familiar stuff (but with slightly different variable names), and i think it uses the "classic SPICE gummel-poon" model, with the "go" and "gm" conductances instead of the "gmf" "gmr".
but i saw some things that i clearly didn't have - one was a mechanism that limits the speed with which vbe and vbc can "move" versus the previous time (so this requires memory) ... i added something like that and it was still exploding
i was looking elsewhere in my code for bugs too and noticed that in my test circuit in the code i feed the oscillator (a ramp) backwards, it's supposed to be going from 11 to 5V but i had it going from 5 to 11V (upward), flipped that and:



looks very familiar
but the transistor is still not quite working properly because it explodes with an upward going ramp.
also, i'm not sure if i've done things properly...

okay, well, i just want to say a huge thank you to KVR people, Andy, mystran, aciddose, kunn, karrikuh, Robin, Urs, and the folks from #musicdsp Hideki, trip-, mike, jazzdude, ashafq, i'm probably forgetting a bunch, forgive me

so this CLI program right now works kinda in 2 steps with some manual hand-work inbetween (but tbh i think it's not bad):
- it reads a netlist, writes node equations in wxMaxima-friendly format into an output file.
- you're supposed to copy/paste them into wxMaxima (and if you know what you're doing you might actually modify some stuff, because i'm not good with Maxima either).
- Maxima solve()s the equations and gives voltage equations and optimize()s them, then this is printed as one line...
- you feed back that one line into the input netlist and re-run the program
- it picks up the voltage equations and fixes the variable names..
so it adds a C++ class with component declarations, init, and a tick() function

initial netlist:

Code: Select all

# circuit-to-code HP+LP+diode-clipper
#float_type: float
#tick_args: const float VS1
VS1 v1 v0
R1 v1 v2, 2.2k
C1 v3 v2, 470n
C2 v3 v0, 10n
R2 v3 v0, 6.8k
D1 v3 v0, isat=1e-15
D2 v0 v3, isat=1e-15

my program vomits:

Code: Select all

kill(all);
R1: 2.2*1000;
gr1: 1.0/R1;
C1: 470/1000000000;
gc1: (2*C1)/dt;
C2: 10/1000000000;
gc2: (2*C2)/dt;
R2: 6.8*1000;
gr2: 1.0/R2;
/* ----------- */
v0: 0;
/* ----------- */
eq0: 0 = +c1_g*(v3-v2)-c1_ieq+c2_g*(v3-v0)-c2_ieq+r2_g*(v3-v0)-d1_gd*(v0-v3)-d1_ieq+d2_gd*(v3-v0)+d2_ieq;
eq1: 0 = -r1_g*(v1-v2)-c1_g*(v3-v2)+c1_ieq;
res: solve([eq0, eq1], [v3, v2]), numer;
res2: ratsimp(res);
res2: factor(res2);
res2: expand(res2);
res2: radcan(res2);
res2: horner(res2);
res2: float(res2);
res3: optimize(res2);
res3: subst(["^"=pow],res3);
display2d:false; negsumdispflag:true;
/* set_display('none); */
printf(false,"#maxima_opt: ~s",res3);
/* ----------- */

i intentionally added the known passive elements to the top so that (if you want to) you can completely solve the result with Maxima, like, if you wanna know what voltage comes out of a voltage divider or something like that, but in general i don't give those values to Maxima, just the stuff from v0 down...

you get something like this back, copy it and add it to the netlist file:

Code: Select all

#maxima_opt: block([%1,%2,%3,%4,%5],%1:(d2_gd+d1_gd+c2_g+c1_g)*r1_g,%2:1/((r1_g+c1_g)*r2_g+%1+c1_g*d2_gd+c1_g*d1_gd+c1_g*c2_g),%3:c1_g*c2_ieq,%4:c1_g*d1_ieq,%5:-1.0*c1_g*d2_ieq,[[v3 = %2*(c1_g*r1_g*v1+((-1.0*d2_ieq)+d1_ieq+c2_ieq+c1_ieq)*r1_g+%5+%4+%3),v2 = %2*((r1_g*r2_g+%1)*v1-1.0*c1_ieq*r2_g+%5-1.0*c1_ieq*d2_gd+%4-1.0*c1_ieq*d1_gd+%3-1.0*c1_ieq*c2_g)]])"

and then my program makes this:

Code: Select all

struct circuit
{
    float dt;
    static constexpr size_t _max_iterations = 128;
    size_t _iter;
    as_bzzzt::resistor<float> r1;
    as_bzzzt::capacitor<float> c1;
    as_bzzzt::capacitor<float> c2;
    as_bzzzt::resistor<float> r2;
    as_bzzzt::diode<float, 1e-15f, 0.026f> d1;
    as_bzzzt::diode<float, 1e-15f, 0.026f> d2;
    // ----
    float v3;
    float v2;
    static constexpr float v0 = 0.0;
    float v1;
    void init(const float Fs)
    {
        dt = 1.f/Fs;
        r1.set_value(2200.f);
        c1.init(Fs, 4.7e-07f, 0);
        c2.init(Fs, 1e-08f, 0);
        r2.set_value(6800.f);
        v3 = 0; // TODO initial value!
        v2 = 0; // TODO initial value!
    }
    void tick(const float VS1)
    {
        v1 = (v0+VS1);
        _iter = 0;
        while (_iter < _max_iterations)
        {
            const float _v3 = v3;
            const float _v2 = v2;
            d1.tick(v0-v3);
            d2.tick(v3-v0);
            // voltage equations from external math solver:
            const float _tmp1 = (d2.gd+d1.gd+c2.g+c1.g)*r1.g;
            const float _tmp2 = 1.f/((r1.g+c1.g)*r2.g+_tmp1+c1.g*d2.gd+c1.g*d1.gd+c1.g*c2.g);
            const float _tmp3 = c1.g*c2.ieq;
            const float _tmp4 = c1.g*d1.ieq;
            const float _tmp5 = -1.f*c1.g*d2.ieq;
            v3 = _tmp2*(c1.g*r1.g*v1+((-1.f*d2.ieq)+d1.ieq+c2.ieq+c1.ieq)*r1.g+_tmp5+_tmp4+_tmp3);
            v2 = _tmp2*((r1.g*r2.g+_tmp1)*v1-1.f*c1.ieq*r2.g+_tmp5-1.f*c1.ieq*d2.gd+_tmp4-1.f*c1.ieq*d1.gd+_tmp3-1.f*c1.ieq*c2.g);
            // end of equations
            ++_iter;
            static constexpr float _quiet = 1e-06;
            if      (std::abs(_v3-v3) > _quiet) { continue; }
            else if (std::abs(_v2-v2) > _quiet) { continue; }
            else { break; }
        }
        c1.save(v3-v2);
        c2.save(v3-v0);
    }
    /*
        VS1 v1 v0
        R1 v1 v2, 2.2k
        C1 v3 v2, 470n
        C2 v3 v0, 10n
        R2 v3 v0, 6.8k
        D1 v3 v0, isat=1e-15
        D2 v0 v3, isat=1e-15
    */
};

after the potentiometer, and especially after the first tries with the transistor, i decided to move the actual runtime implementations of the components into a .hpp instead of emitting them from the CLI tool, and yes, i didn't have a better name for this so i called it "bzzzt"

Code: Select all

template<typename T>
struct resistor
{
    T g;    //!< conductance
    inline void set_value(const T ohms) { g = 1 / ohms; }
};
template<typename T>
struct capacitor
{
    T g;    //!< conductance
    T ieq;  //!< current memory variable
    inline void init(const T Fs, const T farads, const T _ieq = 0) { ieq = _ieq; set_value(farads, Fs); }
    inline void set_value(const T farads, const T Fs) { g = Fs*2*farads; }
    inline void save(const T vd)
    {
        ieq += 2*(g*vd-ieq);
    }
};

template<typename T, T R_total, T R_min = 0.5>
struct rvar // variable resistor
{
    T g;
    static constexpr T R_add = R_total-R_min;
    /// wiper goes 0 to 1 (and resistance grows)
    inline void set_wiper(const T wiper) { g = _calc_vr_g<T, R_add, R_min>(wiper); }
};
template<typename T, T R_total, T R_min = 0.5>
struct rpot // potentiometer
{
    T g_a, g_b;
    static constexpr T R_add = R_total-(R_min+R_min);
    /// wiper goes 0 to 1 (from A to B)
    inline void set_wiper(const T wiper) { g_a = _calc_vr_g<T, R_add, R_min>(wiper); g_b = _calc_vr_g<T, R_add, R_min>(1-wiper); }
};
template <typename T, T R_off=1e+6, T R_on = 0.5>
struct sw_spst : public rvar<T, R_off, R_on> // SPST switch
{
    using base_t = rvar<T, R_off, R_on>;
    static constexpr T lut_g[2] = { _calc_vr_g<T, R_off-R_on, R_on>(0), _calc_vr_g<T, R_off-R_on, R_on>(1) };
    inline void set(bool closed) { base_t::g = lut_g[closed]; }
};
template <typename T, T R_off=1e+6, T R_on = 0.5>
struct sw_spdt : public rpot<T, R_off, R_on> // SPDT switch, terminals A COM B
{
    using base_t = rpot<T, R_off, R_on>;
    static constexpr T lut_g[2] = { _calc_vr_g<T, R_off-R_on, R_on>(0), _calc_vr_g<T, R_off-R_on, R_on>(1) };
    /// 0 == open; 1 == A<-COM; 2 == COM->B; 3 == A<-COM->B (fully closed)
    inline void set(uint8_t mask)
    {
        base_t::g_a = lut_g[(mask   )&0x01];
        base_t::g_b = lut_g[(mask>>1)&0x01];
    }
};

template
<
    typename T,
    T isat = 1e-15, // saturation/reverse current
    T vt = 26e-3 // diode thermal voltage
>
struct diode
{
    static constexpr T vt_inv = static_cast<T>(1)/vt;
    T g; //!< diode conductance
    T ieq;
    inline void tick(const T vd)
    {
        // ed = exp(vd/vt);
        // id = isat*ed-isat;
        // gd = isat*ed/vt;
        // ieq = id-gd*vd;
        const T ed = std::exp(vd*vt_inv);
        const T ed_isat = ed*isat;
        g  = ed_isat*vt_inv;
        ieq = ed_isat-isat-g*vd;
    }
};

[BertKoor]

[mystran]

but basically if there are too many components (and by that i mean like more than 3) i better NOT do it by hand because i'm slow and i'll mess it up and then i will only waste time and not get anywhere

Congratulations, you are not the only one that has reached this conclusion.

[torgover]

https://github.com/hal0zer0/melange

[reubenwbetts]

AUTO-ADMIN: Non-MP3, WAV, OGG, SoundCloud, YouTube, Vimeo, Twitter and Facebook links in this post have been protected automatically. Once the member reaches 5 posts the links will function as normal.

I'm in the same boat, I've written a JUCE solver that parses a netlist and solves it automatically by LU decomp, which matches SPICE for linear circuits, and mostly seems to match SPICE for diodes (although getting the SPICE diode model to actually match mine is somewhat of a challenge), but porting the whole thing into matlab to verify properly is going to be somewhat arduous so haven't got round to it. My main problem is implementing the Eber-molls model, no matter what I do it seems to become unstable/not match SPICE at all. Theres a few possible explanations, but I'm at a loss as to which one it is:

A) my transistor parameters are wrong - I've been unable to work out what an appropriate forward and backward alpha for my current sources are - it could also be the wrong ideality factor for IS term in the diodes?

B) the situation im putting it under is inherently unstable, where I'm driving the base with a pulse(0-15v) that's a higher voltage than the collector (10v). So the transistor is going from closed -> active -> saturation in an incredibly short time.

C) My solver is not correct - it performs the full LU decomposition of the A*V = B matrices and linearisation-process for each newton iteration, and recalculates the diode currents based of the previous iterations calculation of the voltage across the diode. The other thing to note is it linearises the Shockley equation to get its conductance and current, and stamps them into the a matrix and b vector:

Code: Select all (#)

double v_diff = diode->vd;
        
double id = is * exp(v_diff/(n*vt));
double gd = (is/(n*vt)) * exp(v_diff/(n*vt)); // conductance stamps into a_matrix        
double ieq = id - gd * v_diff; // current stamps in b_vector

after stamping these, it then does the LU decomp of the entire circuit, and the voltage difference calculated from that is used to re-evaluate. At the moment I'm just running this 50 times (I don't care about efficiency really I just want it to work lol).

Again my maths isn't particularly great so trying to work out all of this has been a massive headache. I don't think(?) I'm doing newtons method properly, as I don't think I'm calculating the Jacobian etc... but again I'm not even sure if thats the case as my understanding of this is pretty limited.

[mystran]

My main problem is implementing the Eber-molls model, no matter what I do it seems to become unstable/not match SPICE at all. Theres a few possible explanations, but I'm at a loss as to which one it is:

Have you tried the "Qucs trick" (ie. a trick explained in Qucs technical manuals, here: https://qucs.sourceforge.net/tech/node1 ... 0000000000).

The TL;DR is that if you try to solve the operating point for NR linearization from the voltage alone, there's a problem when the linear solver gives you a voltage that's too large, because the exponential gets much larger than it should (and potentially overflows to infinity). On the other hand, if you solve from current, then we run into a similar problem with reverse bias where the slope is close to zero.

So... the solution is to do both: if the predicted voltage from the linear solver is below a threshold (the optimal threshold is pretty much at the knee of the diode curve) then use the predicted voltage directly, but if it's above the knee, solve for the current instead as predicted by the linear solver (which we can get from the predicted voltage and the previous guess linearization) and then use this predicted current with the Ebers-Moll equation backwards to solve for a more sensible voltage that you can use for NR without things blowing up.

ps. If it sounds like first taking a logarithm only to feed it back into an exponent seems expensive... trust me, you'll more than make up for the performance penalty from much faster convergence even where it wouldn't blow up.

[Music Engineer]

https://github.com/hal0zer0/melange

"Melange is a circuit compiler. Give it a SPICE netlist — or a KiCad schematic — and it emits optimized, real-time-safe DSP code. It handles the math so you can focus on your product."

This sounds amazing! It actually even sounds too good to be true! Can it also emit the ODE system? I suppose, it must at some point generate it as an intermediate step anyway, right? I'd personally be most interested in an automated "circuit -> math" translation and would perhaps rather do the "math -> code" translation myself. This is because reading circuit diagrams is outside my skillset - but implementing ODE solvers is kinda within (although I'm not claiming to be good at it - I just know how it generally works). That, together with a library of famous audio circuits, would be a great resource for learning about circuit modeling.

[torgover]

Melange does have an in-memory intermediate representation, and it set up so that the IR can be compiled down into "pretty much whatever", in terms of language. At the moment it's just Rust as the proof-of-concept language, once that's solid it'll do other forms of output (c++, matlab, etc)

[torgover]

Note - still early alpha, no release builds, no announcements, incomplete docs, etc

[antto]

melange - interesting, but i've already written a small bunch of code here

i wanted to add the ability for "series resistance", i tried to cut some corners, but it seems the best way is to literally put hidden resistors into the netlist, that involves creating "internal nodes" and re-routing the original netlist

well, i almost did that now (need to re-wire the resistance value to come from the RS parameter of the (for example) diode instead of being an independant literal)

also, since the parameters especially on the transistor are so many (and i don't even have half of those of actual spice), i moved them into a param struct, and i added "presets" to my netlist, it declares a named preset with parameters

Code: Select all

!preset N4148 pinout=AK, IS=4.352E-9, N=1.906, RS=0.6458

D1 net3 net4, &N4148, TNOM=25

wxMaxima - you might have noticed the "res3: subst(["^"=pow],res3);" chunk earlier - this tries to substitute the "^" power operator of Maxima with "pow()" and when i tested it on some small expression - it seemed to work and even preserved the operator precedence - great! but then some time later it broke something else - it replaced some divisions with pow(xxx,-1).

in any case, i guess my scheme with the intermediate "manual steps" of tossing around math expressions into an external program (wxMaxima in my case) is kinda annoying during the developement of the actual tool, but on the other hand it gives the opportunity for an actual user to interfere and do a better job with the math (if he knows math), which is probably a good thing

Robin: is this maybe what you're looking for? i'm not sure what you mean exactly with "ODE".. or is that more than merely the KCL node equations but also the other equations for the components which allow full computations inside a program like wxMaxima (or other), because in my case i do not want to give the full diode equations to wxMaxima because i don't want to solve the voltage there but at runtime.. yet i can give some of the values to it like the Resistor value and the Voltage value (for fixed voltage rails)

Have you tried the "Qucs trick" (ie. a trick explained in Qucs technical manuals, here: https://qucs.sourceforge.net/tech/node1 ... 0000000000).

mmm, so that's where it comes from, thank you, i need to check it closer

[Marvinh]

very impressive antto and torgover never occurred to me we have to treat v0 v1 separately the + - I always think of it as a single input.

This part here is just for cleaning up "Solve"?
res2: ratsimp(res);
res2: factor(res2);
res2: expand(res2);
res2: radcan(res2);
res2: horner(res2);

[antto]

well, the output from solve() tends to have tons of repeating stuff

these bunch of function calls there are a very simple attempt (and a hint) that some math simplification and/or optimization steps can be plugged at that point. i actually don't know if they help or make things worse, and it might actually be hard to say.
i think a given set of Maxima "post-processing" steps might actually give decent results for one particular circuit, and not so great results for another, but i'm not sure.

i've used ratsimp() expand() and optimize() before, the others i got by describing what i'm doing to a chatbot and he suggested some of those extra ones

the optimize() output for the otherwise "simple" circuit which generated the waveform from my first picture looks like this:



if you think that's a whole screen, well, the output of solve() before that is several screens

btw: you can also probably see that i've been also wrestling with Maxima's tendancy to format the equations in fancy "math notation" which is not cool for copy/pasting into "code-grade math" ...
i saw there's some special "code-gen" plugin in maxima that is supposedly specifically made to help with this kind of thing, but i couldn't figure out how to use it, so i came up with the "print()" function and at least that just writes the whole thing into one line that you must transfer back to the netlist to complete the code-gen. then i fix all the math constants to the target float/double format. Maxima seems to sometimes write "1" other times "1.0", and if my output float type is "float" i need all of these rewritten as "1.f" because on Cortex-M4F for example, the FPU is only single-precision and any accidental "1.0" literal can turn the whole expression into double-precision and that goes above the FPU's head and the performance falls to the floor.

[mystran]

well, the output from solve() tends to have tons of repeating stuff

Yes, this is generally the wrong way to go about it unless your circuit is very simple.

The 'lu_factor' and 'lu_backsub' from the 'linearalgebra' package might give you slightly more useful output, but even that doesn't really make the computational structure entirely obvious.

To avoid repetition, the best approach (IMHO) is to just give up on the math packages entirely and simply write a little code-gen that takes the MNA matrix in symbolic form, this can be as simple as opaque (as far as the code-gen cares) strings vs. canonically zero entries (usually that's enough to find workable pivots, though if the system isn't too large, you can even experiment with pivoting manually to find a solution that leads to less fill-in; note that optimal minimization of fill-in is NP-hard) and perform LU factorization by dumping C++ code for the arithmetic operations directly.

The resulting code might still have a bit of redundancies and such, but the C++ compiler's optimizer can usually clean up the lower level optimizations pretty well.

[Archit3ch]

i really wanted to add a proper potentiometer to it but that required 3 pins.

You can make a potentiometer as 2 variable-resistors.

...and if i can omit some stuff easily (since that page talks about the gummel-poon model, which now i know is not bad, but for audio it might be even overkill, it depends)

Start with Ebers-Moll for audio.

okay, well, i just want to say a huge thank you to KVR people, Andy, mystran, aciddose, kunn, karrikuh, Robin, Urs, and the folks from #musicdsp Hideki, trip-, mike, jazzdude, ashafq, i'm probably forgetting a bunch, forgive me

We have an #analog-modeling channel in The Audio Programmer Discord. There are people posting about their solver attempts.

Code: Select all

struct circuit
{
    void init(const float Fs)
    {
        dt = 1.f/Fs; // don't do fixed rate integration!!
    }
    void tick(const float VS1)
    {
        while (_iter < _max_iterations) // don't quit early! this accepts wrong answers
        {
            static constexpr float _quiet = 1e-06; // abstol only? needs reltol as well
        }
    }
};

See my comments in the code.

Code: Select all

template<typename T>
struct capacitor
{
    inline void save(const T vd)
    {
        ieq += 2*(g*vd-ieq); // hardcoded trapezoidal integrator. not bad, but this should be abstracted out as an integrator choice
    }
};

See my comment in the code.

https://github.com/hal0zer0/melange

Looks vibecoded. There are plenty of mistakes, scaling problems, etc.

To be fair, it looks better than "everyone's first circuit solver".

Again my maths isn't particularly great so trying to work out all of this has been a massive headache. I don't think(?) I'm doing newtons method properly, as I don't think I'm calculating the Jacobian etc... but again I'm not even sure if thats the case as my understanding of this is pretty limited.

You should be able to match SPICE (accounting for different models, added gmin, etc.). Calculating the Jacobian is important.

Have you tried the "Qucs trick" (ie. a trick explained in Qucs technical manuals, here: https://qucs.sourceforge.net/tech/node1 ... 0000000000).

The TL;DR is that if you try to solve the operating point for NR linearization from the voltage alone, there's a problem when the linear solver gives you a voltage that's too large, because the exponential gets much larger than it should (and potentially overflows to infinity). On the other hand, if you solve from current, then we run into a similar problem with reverse bias where the slope is close to zero.

Yeah, everyone with hand-rolled solvers runs into this eventually.

I use industrial-grade solvers. You can pretty much throw a netlist at them to brute force. No need to resort to hacks for convergence. But it can still be important for efficiency.

"Melange is a circuit compiler. Give it a SPICE netlist — or a KiCad schematic — and it emits optimized, real-time-safe DSP code. It handles the math so you can focus on your product."

This sounds amazing! It actually even sounds too good to be true!

It's hardly novel! Livespice came out in 2013. Very interesting quasi-symbolic approach, too...

Can it also emit the ODE system? I suppose, it must at some point generate it as an intermediate step anyway, right? I'd personally be most interested in an automated "circuit -> math" translation and would perhaps rather do the "math -> code" translation myself. This is because reading circuit diagrams is outside my skillset - but implementing ODE solvers is kinda within (although I'm not claiming to be good at it - I just know how it generally works).

Yup! I've settled on having an intermediate stage for symbolic math as well.

My steps are roughly:
1) Pantelides index reduction
2) Symbolic simplification (sympy/Mathematica)
3) Manual re-writes (very system-dependent, not automated)
4) Schur linear/nonlinear split

That, together with a library of famous audio circuits, would be a great resource for learning about circuit modeling.

I should write a book.

Melange does have an in-memory intermediate representation, and it set up so that the IR can be compiled down into "pretty much whatever", in terms of language. At the moment it's just Rust as the proof-of-concept language, once that's solid it'll do other forms of output (c++, matlab, etc)

IR is great for the last step (=code generation).

But for symbolic expressions, you want something optimized for them. Even then, symbolic-solve speed will explode for large systems.

in any case, i guess my scheme with the intermediate "manual steps" of tossing around math expressions into an external program (wxMaxima in my case) is kinda annoying during the developement of the actual tool, but on the other hand it gives the opportunity for an actual user to interfere and do a better job with the math (if he knows math), which is probably a good thing

Yes, as I said. You want to have that manual step.

...and if my output float type is "float" i need all of these rewritten as "1.f" because on Cortex-M4F for example, the FPU is only single-precision and any accidental "1.0" literal can turn the whole expression into double-precision and that goes above the FPU's head and the performance falls to the floor.

Uh, none of that is gonna run on a Cortex-M4F.

With the possible exception of the simplest of simple circuits.

Yes, this is generally the wrong way to go about it unless your circuit is very simple.

For very simple circuits, compute the sparse symbolic inverse.

To avoid repetition, the best approach (IMHO) is to just give up on the math packages entirely and simply write a little code-gen that takes the MNA matrix in symbolic form, this can be as simple as opaque (as far as the code-gen cares) strings vs. canonically zero entries (usually that's enough to find workable pivots, though if the system isn't too large, you can even experiment with pivoting manually to find a solution that leads to less fill-in; note that optimal minimization of fill-in is NP-hard) and perform LU factorization by dumping C++ code for the arithmetic operations directly.

I believe Andy does something similar.

IMO, start from the circuit equations in DAE form instead of MNA. Then you run the 2 steps I wrote above (Mathematica + manual symbolic re-writes). Then you do the Schur/W codegen. In my case it's Julia, not C++.

For pivoting, I just run the simulation. Telemetry has given me optimal pivots so far. Kinda cheating that one can use runtime information to apply compile-time optimization for the next pass. PGO is a valid optimization?

The resulting code might still have a bit of redundancies and such, but the C++ compiler's optimizer can usually clean up the lower level optimizations pretty well.

Yup. I'm using non-C++ LLVM. It does a good job of picking up repeat computation at the lowest level after codegen.

[mystran]

Have you tried the "Qucs trick" (ie. a trick explained in Qucs technical manuals, here: https://qucs.sourceforge.net/tech/node1 ... 0000000000).

[...]

Yeah, everyone with hand-rolled solvers runs into this eventually.

I use industrial-grade solvers. You can pretty much throw a netlist at them to brute force. No need to resort to hacks for convergence. But it can still be important for efficiency.

You will run into this problem whether you "hand roll" your solver or not (at least if you use a SPICE-like solver with NR) and the suggested fix is not a "hack" as much as a very elegant solution to make sure that NR always converges optimally without having to resort to actual hacks such as clamping the forward voltage, falling back to a small time-step, etc.

Yes, as I said. You want to have that manual step.

I think you should really not have any manual steps at all if you can help it.

Manual steps only serve to introduce errors.

Yes, this is generally the wrong way to go about it unless your circuit is very simple.

For very simple circuits, compute the sparse symbolic inverse.

This only works well [edit: in the sense of giving you reasonably efficient code] for very simple circuits and even then LU usually gives better code.

IMO, start from the circuit equations in DAE form instead of MNA. Then you run the 2 steps I wrote above (Mathematica + manual symbolic re-writes). Then you do the Schur/W codegen. In my case it's Julia, not C++.

Well, the MNA form effectively is nothing but the original DAE with NR linearizations substituted for non-linearities. I would sometimes go further and also substitute the integrators directly into the matrix, because this can actually lead to slightly better code (mainly because you need to solve less dimensions when the solves produces new state variables directly).

You might think this means I can no longer easily switch between integrators? Well, worry not... 'cos if you full automate the whole netlist->MNA->code path (or even schematic->netlist->MNA->code; turning a graphical schematic into a netlist is easy, the only difficult part is writing a workable schematic editor) then ... building a different MNA matrix to substitute a different integrator is no longer a big deal.

For pivoting, I just run the simulation. Telemetry has given me optimal pivots so far. Kinda cheating that one can use runtime information to apply compile-time optimization for the next pass. PGO is a valid optimization?

Any optimization is valid in my book... although if you just bruteforce every combinatino that's a lot of combinations to try even for a fairly small system. It should perhaps also be pointed out that partial pivots alone won't usually give you the optimal solution.

That said... in practice I've come to the conclusion that the difference between optimal and "whatever works" pivots is usually not that huge, because usually it's the NR iteration that takes most of the time and the submatrix you need to iterate tends to end up dense anyway (at least if you make sure to minimize the dimensions, but that's a bigger when than trying to keep things sparse).

ps. Also my two cents is this: if you're using NR, then your first priority with regards to optimization should always be trying to minimize the number of NR iterations you need to take and minimizing the number of dimensions you actually need to solve during your NR iteration. If you can improve on either of these by making the rest of the circuit a bit more complicated, it's pretty much always worth it.