# Nonlinear equation solvers in SPICE simulators

We have an assignment at the Parallel processing class, the target is to implement a Non-linear Equations solver on cuda based on Newton Raphson method and to interface this solver with an application that deals with nonlinear set of equations. We wanted to interface our solver with circuit simulators. We have picked up an open source simulator and every time the simulator is performing a DC operating point simulation it will invoke our cuda code. At this point we wanted to compare the performance of our solver against solvers implemented in other circuit simulators such as

• LTspice
• Ngspice
• Qucs

And also other software solvers e.g. Matlab optimization toolbox

We have tested these solvers against a circuit [which should map to a huge set of nonlinear eqns.]

The circuit netlist is generated by a script where the number of the nodes is given. We have formulated the set of nonlinear equations governing the above circuit as the following

The unknowns x[i] in this set of equations are the voltage nodes and the current at each resistor

We have managed to write this as a matlab function to test this circuit against matlab nonlinear solver algorithms featured in the optimization toolbox.

function F = non_linear_diode(X)
% Len(x) is always even 2*d
d = length(X)/2;
F = zeros(1, d*2);
i = 1e-3;     % Current source magnitude
r = 50;       % Resistors value
c1 = 1e-15;   % Diodes I_s
c2 = 0.0258;  % Diodes N*V_th

F(1) = X(1) - X(2) - i*r;
F(d) = X(d) - X(d+1) - X(2*d)*r;
F(d+1) = i - c1*(exp(X(2)/c2) -1) - X(d+2);

for ii = 2:(d -1)
F(ii) = X(ii) - X(ii+1) - X(ii+d)*r;
F(ii+d) = X(ii+d) - c1*(exp(X(ii+1)/c2) -1) - X(ii + d + 1);
end
F(d*2) = X(2*d) - c1*(exp(X(d+1)/c2) -1);
end


We have also written a script to generate Spice netlist for this problem

def gen_ckt(num):
ret = ""
for i in range(1, num):
ret += 'R'+str(i)+" "+str(i)+" "+str(i+1)+" 50\n"
ret += 'D'+str(i)+" "+str(i+1)+" 0 DI1N4004\n"


Where DI1N4004 is our diode model defined in the netlist. When testing the above solvers against the problem with 70,000 nodes i.e. 140,000 equations and unknowns

• Matlab runs out of memory
• Qucs takes forever
• All the spice based solvers somehow solved this problem in less than 2 seconds

We actually have no idea how spice solvers managed to avoid this memory problem, and even when testing this problem against fewer number of unknowns e.g. 3,000 the spice solvers always outperform matlab and qucs. Though as mentioned in [1] [2] [3] spice uses the damped Newton-Raphson approach to solve circuits with nonlinear components which is the same as all the solvers mentioned above

• Matlab Fsolve : dogleg method [Newton + Trust-region + steepest decent] [4]
• Qucs : damped Newton-Raphson [5]

My questions are

• How spice managed to solve such a huge system of nonlinear equations quickly and without running out of memory ? is it exploiting the circuit structure making use of the repeated elements ?
• Is this example fair enough ? i mean do we need to consider more practical circuits ? and if so can any one give an example for a circuit(s) where DC operating point simulation might be the bottleneck of the simulation time ?
• In [6] the authors mentioned ISCAS85 benchmark circuits should we consider these circuits in our tests ?
• Is the DC operating point the simulation type we should be targeting ? i mean should we focus to interface our solver to other types of simulations e.g. the transient analysis ?

References

1 : http://www.ni.com/white-paper/5808/en/

2 : http://www.ecircuitcenter.com/SpiceTopics/Overview/Overview.htm

3 : http://www.electronicdesign.com/boards/taking-peek-under-hood-your-spice-circuit-simulation-engine

4 : https://www.mathworks.com/help/optim/ug/fsolve.html

5 : http://qucs.sourceforge.net/tech/node16.html

6 : http://www.mos-ak.org/bucharest/presetnations/Lannutti_MOS-AK_Bucharest.pdf

• I don't have time at the moment to answer, but this book can give you a good idea of how circuit simulators do their magic, I highly recommend it: amazon.com/Circuit-Simulation-Farid-N-Najm/dp/0470538716/… – John D Aug 17 '17 at 21:36
• You might consider asking (a subset?) of this question on computational science – The Photon Aug 18 '17 at 0:42
• This isn't really an electronics question... SPICE isn't even just limited to electronics AFAIK. Shouldn't it be on Math.SE or SciComp.SE or CS.SE or something else? – user541686 Aug 18 '17 at 5:16
• @Mehrdad i totally agree with ThePhoton only part of the question should be asked on cse.se but considering the rest of the questions i think this is the best place to ask them. – Elbehery Aug 18 '17 at 12:35

How spice managed to solve such a huge system of nonlinear equations quickly and without running out of memory ?

A good SPICE very likely defaults (or knows when to switch to) a sparse matrix solver, since large circuits typically produce sparse matrices (each node connects to only a small fraction of the branches) it would be an obvious optimization to use (or have available) a sparse matrix solver in a SPICE. Even Nagel's original 1975 report (Thesis?) on SPICE discusses using sparse matrix methods.

Matlab certainly has a sparse matrix solver available, but you probably have to invoke it explicitly.

Qucs may not have this capability, or it might not be implemented particularly well, because it's a relatively raw open source project and its developers might not have got to the point of testing it on anything bigger than a toy problem.

(Hat tip to @jonk for the link to the Nagel report)

Is this example fair enough ? i mean do we need to consider more practical circuits ?

I would think you want to demonstrate your solver on a wide variety of different types of circuits. You probably want to consider circuits that are known to produce ill-conditioned matrices. Positive feedback circuits also often produce difficulties for non-linear solvers.

and if so can any one give an example for a circuit(s) where DC operating point simulation might be the bottleneck of the simulation time?

I would expec this to be common in any ill-conditioned circuit when setting up an AC simulation.

Is the DC operating point the simulation type we should be targeting ? i mean should we focus to interface our solver to other types of simulations e.g. the transient analysis ?

The other main simulation types (ac and transient) only require linear solvers. The AC simulation explicitly is about small variations about the operating point, so that the circuit can can be considered linear by perturbation theory. The transient solver linearizes the circuit at each time step, but re-calculates the local linear equivalent circuit for each time step. So if you're trying to demonstrate a non-linear solver, the DC solver is the one to demonstrate.

• Hat tip back to you for finding an online version of Nagel's paper. (I think it was a Thesis, but then morphed into a department memorandum of some kind.) – jonk Aug 18 '17 at 0:19
• @jonk, I think you gave me the link to it a couple weeks (months?) ago. – The Photon Aug 18 '17 at 0:33
• I think I gave the reference itself, because I have a complete copy of it in paper form that I got from Berkeley, directly, well more than a decade ago. Perhaps it wasn't you that found the link. But I definitely remember I didn't know about the web version until I saw it posted here on June 20th, this year. (Grabbed it, of course, as I like being able to search for words.) – jonk Aug 18 '17 at 1:01
• These days, sparse solvers are used for almost everything, including dense matrices. But looking at the OP's circuit diagram the system matrices is probably banded with a very small bandwidth - it might even by tridiagonal. If you ignored that fact in your Matlab code, you could expect speedup factors of 1000 or bigger with very little effort, even using 1960s-era numerical algorithms, not modern ones that are harder to understand! – alephzero Aug 18 '17 at 3:43
• @Elbehery, to get better performance from Matlab you may have to be very careful about never making it store the matrix as a dense matrix, or making it convert between dense and sparse representations. Unfortunately I don't know enough about Matlab to give more detailed advice. – The Photon Aug 28 '17 at 21:49