Circuit Simulation Part IV: Nonlinear circuits and implementation

In the last part we presented the major part of the UI, but we still hadn’t talked about the implementation of solving. In this post we will talk about how to actually solve non-linear equations both in theory and in implementation.

Non-linear circuits

Since last time, I took a pass at cleaning up the code and I decided to implement a few more components (inductors and diodes). Inductors were relatively easy since it is still linear, but the diode presented more problems. Let’s go through the theory of non-linear solves.

We started in Part I laying out how to devise an equation per KCL node. Each equation becomes a row in the matrix and then we just solve a linear system to get the final values. That’s all good and fine, but it only works for linear components like resistors, inductors, capacitors and simple source terms like voltages and currents. For things like diodes and transistors which have crap loads of annoying non-linear exponential crap in them, we are kinda screwed. So, we must turn to nonlinear solving techniques, in particular Newton-Raphson iteration which root solves non linear equations by solving linear approximations and iterating. Let’s do a brief overview to see what we’ll need. Let’s suppose the ith equation
$\mathbf{F}_i(\mathbf{x}) = 0$
But since $\mathbf{F}$ this is nonlinear. So using the taylor approximation we can write an approximation to the function as
$\mathbf{F}_i(\mathbf{x}) \approx \mathbf{F}_i(\mathbf{x}^{n}) + \left. \frac{\partial \mathbf{F}}{\partial \mathbf{x}} \right|_{\mathbf{x}^n}\delta\mathbf{x}$
where $\delta\mathbf{x} = \mathbf{x}^{n+1} - \mathbf{x}^n$ and $\partial\mathbf{F}/\partial\mathbf{x}|_{\mathbf{x}^n}$ is the Jacobian at the current iterate.
We still want the root to the approximation so we need to now solve
$\mathbf{F}_i(\mathbf{x}^{n}) + \left. \frac{\partial \mathbf{F}}{\partial \mathbf{x}} \right|_{\mathbf{x}^n}\delta\mathbf{x} = 0$
or specifically
$\left. \frac{\partial \mathbf{F}}{\partial \mathbf{x}} \right|_{\mathbf{x}^n}\delta\mathbf{x} = -\mathbf{F}_i(\mathbf{x}^{n})$
Colloquially this is just a matrix equation $\mathbf{A x}=\mathbf{b}$ . However, every Newton iteration we need to form a new A matrix. Thus given a current solution state $\mathbf{x}_i$ we can go to a new one by following the steps:

Solve $\left( \left. \frac{\partial \mathbf{F}}{\partial \mathbf{x}}\right|\mathbf{x}^n \right) \mathbf{\delta x}= -\mathbf{F}( \mathbf{x}^n )$
Update the solution to be $\mathbf{x}^{n+1} = \mathbf{x}^n + \mathbf{\delta x}$

On each time step we just repeat 1 and 2 as many times as needed for convergence given some norm. I just use check if $\mathbf{\delta x}$ is small for now. So to implement this, we need a way to form the matrix, a way to compute derivatives and a way to solve matrices.

Implementation

Building a matrix

Conceptually, we need to form a matrix every time step and every Newton step. This requires two steps, identifying all the unknowns (which defines matrix dimensions), followed by determining the coefficients of the matrices. For various reasons (that probably are invalid and I will have to fix later) I thought it would be good to use separate objects for the solve and the UI. Thus, in the last segment we had a function that took “Symbol” objects and produced a final description like:

["R",["n1","n2"],1000]
["C",["n2","GND"],0.0001]
["V",["n1","GND"],"square,0.1,-1,1"]

We will go through the list and make “solver” objects now. It looks like this:

function createCircuitFromNetList(names,components){
    circuit=new Circuit();
    // Make a list of all the voltage nodes to graph
    namesToGraph=[]
    for(var n in names) namesToGraph.push(n);
    // Take visual components and turn them into simulation components
    for(var n in components){
        var node=components[n];
        switch(node[0]){
            case "R": circuit.addR(node[1][0],node[1][1],node[2]);break;
            case "C": circuit.addC(node[1][0],node[1][1],node[2]);break;
            case "V": circuit.addVoltage(node[1][0],node[1][1],node[2]);break;
            case "D": circuit.addD(node[1][0],node[1][1],node[2]);break;
            case "L": circuit.addI(node[1][0],node[1][1],node[2]);break;
        }
    }
    // solve and graph the circuit
    var foo=circuit.transient(3.0,.005,namesToGraph);
    graph(0,2,namesToGraph,foo);
}

Now essentially, we have connected up our UI part to our solver part. From here on out we are going to talk about the parts of the solver. It is possible to use those parts without a UI in a sense, so we basically have some compartmentalization!

Let’s look at what addR and addVoltage do because they are interesting:

Circuit.prototype.addVoltage=function(name1,name2,value){
    n1=this.allocNode(name1)
    n12=this.allocNode("sc"+this.num)
    n2=this.allocNode(name2)
    this.components.push(new Voltage(n1,n12,n2,value))
}
Circuit.prototype.addR=function(name1,name2,R){
    var n1=this.allocNode(name1)
    var n2=this.allocNode(name2);
    this.components.push(new Resistor(n1,n2,R))
}

The allocNode assigns that named node a number in the unknown vector. Once that is done, we make the component object using only the numbers. Each object knows how to contribute to the matrix and which values it contributes to. So in the example above we’d probably assign the matrix unknowns as [n1,n2,GND,sc3]. And remember these unknowns are going to actually be the delta in those values rather than the values themselves. Let’s look at how Resistor is implemented:

function Resistor(node1,node2,R){
    this.node1=node1
    this.node2=node2
    this.oneOverR=1.0/R
}
Resistor.prototype.matrix=function(dt, time, system, vPrev, xOld){
    system.addToMatrix(this.node1,this.node1,this.oneOverR);
    system.addToMatrix(this.node1,this.node2,-this.oneOverR);
    system.addToMatrix(this.node2,this.node1,-this.oneOverR);
    system.addToMatrix(this.node2,this.node2,this.oneOverR);
    system.addToB(this.node1,-vPrev[this.node1]*this.oneOverR+vPrev[this.node2]*this.oneOverR);
    system.addToB(this.node2,+vPrev[this.node1]*this.oneOverR-vPrev[this.node2]*this.oneOverR);
}

The constructor caches the nodes and the conductance (inverse of resistance). The matrix function is where all the magic happens. The function is called from the global build matrix function for all components. Several arguments are provided which are important to be able to add to the matrix. dt is the time step, time is the current time we are simulating at (good for voltage sources), system is an interface for adding values to the matrix, vPrev is the last newton iteration value ( $x^{n-1}$ in notation above). Why am I using v instead of x? Because I was stupid when I started writing this and thought everything would be a voltage, but clearly that can’t work. system provides two functions to use addToMatrix() and addToB that add stuff to the A matrix and the b right-hand-side, respectively. Everything we do to the matrix is opposite and equal and thus we always end up with a symmetric matrix, which is to be expected from a physical system. Also note that it is the linearity of the derivative that lets us differentiate the contribution to the full nodal equations for just the one component we are considering here.

Let’s dig into the math of the derivatives for a second. The contribution to the current of a resistor to a node will be I=V/R based on ohm’s law. So suppose we have a resistor for $n_1$ to $n_2$ . Then we have a contribution of $n_2/R-n_1/R$ on n1’s net current and $n1/R-n2/R$ net current. Differentiating with respect to each variable yields on each equation gives the Jacobian. These nodal equation derivatives only contribute to the rows corresponding to the KCL current and the only columns that are non-zero are the ones where the variable is n1 or n2. The derivative w.r.t. to either n1 or n2 gives 1/R (this is why we cached the conductance). Thus if this was the only equation our partial Jacobian matrix would look like
$\left(\begin{array}{ccccc} 0& \vdots & 0& \vdots & 0 \\ \cdots & 1/R & \cdots & -1/R & \cdots \\ 0& \vdots & 0& \vdots & 0\\ \cdots & -1/R & \cdots & 1/R & \vdots \\ 0 & \vdots & 0& \vdots & 0\\ \end{array} \right)$
You can think of each component producing its own Jacobian matrix and just adding them all. Implementing the code like that would be very slow, so instead, we have this addToMatrix function do the same thing for us. On the right hand side we just evaluate the ohm’s law contributions with respect to the vPrev values. Thus our final system contribution from just the resistance in this particular case is
$\left(\begin{array}{cccc} 1/R & -1/R & 0 & 0 \\ -1/R & 1/R & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ \end{array}\right)\left(\begin{array}{c} n_1 \\ n_2 \\ \textrm{GND} \\ \textrm{sc3} \\ \end{array}\right) = \left( \begin{array}{c} -\textrm{vPre}_0/R+\textrm{vPre}_1/R \\ \textrm{vPre}_0/R-\textrm{vPre}_1/R \\ 0 \\ 0 \\ \end{array}\right)$ .
This is consistent with the code presented above.

Now that we’ve built the matrix for the resistor, let’s build it for the voltage source. The voltage has three matrix equations it contributes to. The two current equations for the nodes it is connected to and its own current which is unknown. From our example we have that $n_1-n_2=v(t)$ where v is the voltage source that only depends on time. But we also have the unknown sc3 which contributes to the KCL current equation for n1 and GND. Thus this function

Voltage.prototype.matrix=function(dt,time,system,vprev,vold){
    // the constraint that the voltage from node1 to node2 is 5
    var v=5; // default to DC of constant 5 v
    if(this.type=="sin"){ 
    	v=(Math.sin(time*2*Math.PI/this.period)+1)*(this.max-this.min)*.5+this.min
    }
    // constrain the voltage
    system.addToMatrix(this.nodeCurrent,this.node1,1);
    system.addToMatrix(this.nodeCurrent,this.node2,-1);
    system.addToB(this.nodeCurrent,v-(vprev[this.node1]-vprev[this.node2]));
    // nodeCurrent is the current through the voltage source. it needs to be added to the KVL of node 1 and node2
    system.addToMatrix(this.node1,this.nodeCurrent,-1);
    system.addToMatrix(this.node2,this.nodeCurrent,1);
    system.addToB(this.node1,vprev[this.nodeCurrent]);
    system.addToB(this.node2,-vprev[this.nodeCurrent]);
}

yields on our example the matrix
$\left(\begin{array}{cccc} 0 & 0 & 0 & -1 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 1 & -1 & 0 & 0 \\ \end{array}\right)\left(\begin{array}{c} n_1 \\ n_2 \\ \textrm{GND} \\ \textrm{sc3} \\ \end{array}\right) = \left( \begin{array}{c} 0 \\ 0\\ 0 \\ v(t) \\ \end{array}\right)$ .
It is somewhat weird that GND has all zeroes after my lecture about equal and opposite and symmetry. Well, the upshot is that GND is set to be zero, so we don’t let the matrix or right hand side get any contributions there (it is in the implementation of addToMatrix and addToB). At the end of the matrix construction we set the on-diagonal element for the gnd node row to be 1 which means that we enforce the equation gnd=0.

Each component has a similar implementation. Some are more complicated and need vOld as well (capacitors need to computer voltage derivatives). Since everybody has a definition, we have a simple buildMatrix function that makes the whole matrix:

Circuit.prototype.buildMatrix=function(vPrev,vOld,dt,time){
    var system=new System(this.num,this.ground);
    this.sys=system;
    for(var comp=0;comp<this.components.length;comp++){
        this.components[comp].matrix(dt,time,system,vPrev,vOld);
    }
}

Solving a matrix

Now we have constructed the matrix in our System object (which you need to look at the source code for). It basically consists of a array of arrays member called A and a array member consisting of b. We can solve the matrix equation using Gaussian elimination which we have just implemented in a straightforward way using partial pivoting (column swapping). Nothing exciting here, you can go look at wikipedia for more information on Gaussian elimination. That being said, using a dense matrix and gaussian elimination may not be the most efficient solver here. At least using a sparse LU or other optimized direct solver would be better. You could use LAPACK or SparseLU or whatever if you were in C/C++/Fortran. I am just in it for a quick and fast demo, but many opportunities for improvement are available which would be essential to make large circuits (1000’s of components) really efficient.

The Newton Solve

The full loop is implemented in transient() and it basically follows these steps:

function transient(maxTime,dt){
    vOld=vector(N,0); // N entries that are all zero
    vPrev=vector(N,0); // N entries that are all zero
    while(time<maxTime){
        vPrev.copyVector(vOld);
        for(var newton=0;newton<maxNewtonSteps;newton++){
            var matrix=buildMatrix(vPrev,vOld,dt,time);
            var difference=solve(matrix); /
            var maxDiff=difference.maxNorm();
            vPrev.addVector(difference)
            if(maxDiff < diffTolerance) break; // converged!
        }
        vOld=vPrev; // copy values
        time+=dt;
    }
    return vOld;
}

This is basically everything that is needed. There are several things that can go wrong. If you don’t specify a gnd, the solve will not work. If you make a loop of capacitors it probably won’t work. If you don’t have everything connected in a consistent way, it won’t work. There is very little error checking right now. It would be good to check if ground is connected. It would be good to check if the matrix is non-singular on every newton step.

A Diode

Before we finish let’s go through how the diode is implemented. This is a very important part because it is usually not talked about in simple circuit solvers that only do linear elements. A diode can be modeled by the Shockley equation that relates current to voltage
$I=I_S \left( e^{\frac{V_D}{n V_T}} - 1\right)$
where $I_S$ is the saturation current of the diode, $V_T$ is the thermal voltage (26 mV at room temperature according to wikipedia) and $n$ varies between 1 and 2 to represent how ideal the diode is. Suppose we have nodal voltages n1 and n2, then we have $I_S \left( e^{\frac{n_1-n_2}{n V_T}} - 1\right)$ . The derivatives with respect to n1 and n2 are
$\pm \frac{I_S}{n V_T} e^{\frac{n_1-n_2}{n V_T}}$ , respectively.

All this fun can be encoded into a component function that looks like this:

Diode.prototype.matrix=function(dt, time, system, vPrev, vOld){
	var n1=this.node1;
	var n2=this.node2;
	var denomInv=1./this.n*this.VT;
	var value=this.IS*(Math.exp((vPrev[n1]-vPrev[n2])*denomInv)-1)
	var deriv=this.IS*(Math.exp((vPrev[n1]-vPrev[n2])*denomInv))*denomInv;
	system.addToMatrix(this.node1,this.node1,deriv);
	system.addToMatrix(this.node1,this.node2,-deriv);
	system.addToMatrix(this.node2,this.node1,-deriv);
	system.addToMatrix(this.node2,this.node2,deriv);
	system.addToB(this.node1,-value);
	system.addToB(this.node2,value);
}

Using this we can now make a bridge rectifier circuit plus filter to turn AC into DC currents with minimal ripple.

Graphing

The transient() function populates arrays for each value it is requested to collect. These are then send to the graph() routine which is basically a html5 canvas draw routine that puts those on the graphs on the right. See the code for more details, nothing super exciting there. I would like to make the graphing selective and have the ability to graph some quantities on the same graph window, but this suffices for now. I would also like to extend the simulator to produce currents for every element as well. That will likely require some rework to how the matrices are created which will make the union find from the last post less important.

Try it yourself

You can now try the program online at http://www.andyselle.com/circuit/.
There is still a lot more to do to make this an awesome circuit simulator, but it’s basically working. Yay! Here’s another example I did with some weird non-linear inductor behavior