Implementing Newton-Raphson Method for Nonlinear Circuit Simulation with Diode Components

Question

I am working on a circuit simulator based on Modified Nodal Analysis (MNA), aiming to support nonlinear components such as diodes. To solve the nonlinear equations, I have implemented the Newton-Raphson method. However, I am facing issues with my implementation. Here is a simplified version of the code I am using to solve the below circuit:

// Function to calculate the diode current using the Shockley diode equation
double diodeEquation(double vd, double is, double vt) {
    return is * (exp(vd / vt) - 1.0);
}

// Function to calculate the diode current derivative w.r.t. voltage using the Shockley diode equation
double diodeDerivative(double vd, double is, double vt) {
    return is * exp(vd / vt) / vt;
}

int main() {
    // Step 1: Set up the MNA equations
    // Circuit components
    double vs = 5.0;  // Voltage source value
    double r1 = 10.0; // Resistor R1 value
    double r2 = 20.0; // Resistor R2 value
    double r3 = 30.0; // Resistor R3 value

    // Diode parameters for the Shockley diode equation
    double is = 1e-12;  // Saturation current
    double vt = 25.85e-3; // Thermal voltage

    // Step 2: Initialize variables
    double v1 = 0.0; // Initial guess for node voltage V1
    double v2 = 0.0; // Initial guess for node voltage V2
    double v3 = 0.0; // Initial guess for node voltage V3
    double v1_old, v2_old, v3_old; // Previous values for convergence check
    double tol = 1e-6; // Tolerance for convergence
    int maxIterations = 100; // Maximum number of iterations
    int iterations = 0; // Iteration counter
    bool converged = false; // Convergence flag

    // Step 4: Perform the Modified Newton-Raphson iteration
    while (!converged && iterations < maxIterations) {
        // Save the previous node voltages for convergence check
        v1_old = v1;
        v2_old = v2;
        v3_old = v3;

        // Calculate the diode current and its derivative
        double vd = v2 - v3;
        double id = diodeEquation(vd, is, vt);
        double g = diodeDerivative(vd, is, vt);

        // Formulate the MNA equations
        double i1 = (v1 - vs) / r1;
        double i2 = (v2 - v3) / r2;
        double i3 = v3 / r3;
        double i4 = id;

        // Construct the Jacobian matrix
        double j11 = 1.0 / r1;
        double j12 = -1.0 / r1;
        double j13 = 0.0;
        double j21 = 1.0 / r2;
        double j22 = -1.0 / r2 - g;
        double j23 = g;
        double j31 = 0.0;
        double j32 = g;
        double j33 = -g;
        
        // Construct the residual vector
        double r1 = i1 - i2;
        double r2 = -i1 + i2 - i3;
        double r3 = -i2 + i3 - i4;

        // Solve the linear system of equations to obtain the update
        double detJ = j11 * (j22 * j33 - j32 * j23) - j12 * (j21 * j33 - j31 * j23) + j13 * (j21 * j32 - j31 * j22);
        double dv1 = ((j22 * j33 - j32 * j23) * r1 - (j12 * j33 - j32 * j13) * r2 + (j12 * j23 - j22 * j13) * r3) / detJ;
        double dv2 = (-(j21 * j33 - j31 * j23) * r1 + (j11 * j33 - j31 * j13) * r2 - (j11 * j23 - j21 * j13) * r3) / detJ;
        double dv3 = ((j21 * j32 - j31 * j22) * r1 - (j11 * j32 - j31 * j12) * r2 + (j11 * j22 - j21 * j12) * r3) / detJ;

        // Update the solution vector
        v1 -= dv1;
        v2 -= dv2;
        v3 -= dv3;

        // Check for convergence
        double norm = sqrt(pow(v1 - v1_old, 2) + pow(v2 - v2_old, 2) + pow(v3 - v3_old, 2));
        if (norm < tol) {
            converged = true;
        }

        iterations++;
    }

    // Step 5: Output the final solution
    if (converged) {
        std::cout << "Converged to solution:\n";
        std::cout << "V1 = " << v1 << " V\n";
        std::cout << "V2 = " << v2 << " V\n";
        std::cout << "V3 = " << v3 << " V\n";
    }

    return 0;
}

I have tested the same circuit in two other circuit simulators, The first simulator reports Vd = 739 mV and Id = 70.74 mA, and the second simulator reports Vd = 765 mV and Id = 70 mA. However, in my implementation, the output is Vd = 175 mV and Id = 2.14 µA.

I would appreciate any insights into what might be wrong with my implementation of the Newton-Raphson method for solving the circuit equations, especially with regards to the diode component.

Edit:

After considering the suggestions from EvgsupportsModerator, I have formulated the equations based on Kirchhoff's Voltage Law (KVL) and Kirchhoff's Current Law (KCL). The equations for each node are as follows:

KCL at node 1: I + (v1 - v2)/R1 = 0 // I is the current through the voltage source
KCL at node 2: (v2 - v1)/R1 + v2/R2 + Is*(exp((v2 - v3)/vt) - 1) = 0
KCL at node 3: -Is*(exp((v2 - v3)/vt) - 1) + v3/R3 = 0
KVL at node 1: v1 - E = 0

I obtained the following vector solution when I execute the previous code that I provide:

v1 = 5
v2 = 3.33142
v3 = 0.00861309
I = -0.166858

However, when I substituted these values back into the system, I found that they do not satisfy all of the equations. They only satisfy the first and last equations.

To validate my calculations, I simulated the same circuit using the ngspice software, and it provided a different vector solution:

v1 = 5
v2 = 2.808430
v3 = 2.362064
I = -0.219157

Interestingly, when I substituted these values into the system, it appeared to be a correct solution. Therefore, it seems that the values obtained from the ngspice simulation align better with the given equations.

Answer 1

First of all, you have only TWO free parameters, not four. For network problems (water flow, electric circuits, ... ) these are most conveniently taken as the driving potentials at nodes, since currents/flows depend only on the difference between these potentials and the resistances in the connectors. In this problem, the best two parameters are V2 and V3 (voltages at nodes 2 and 3).

Secondly, a diode is a highly NON-LINEAR device. Thus, Newton-Raphson is not necessarily a good way to solve this problem. (By hand, I used single-point iteration, since the inverted diode formula gives a logarithm, which changes relatively slowly.)

In the end I did manage to solve your problem by Newton-Raphson, but I had to do two key things. (1) I needed a decent initial guess (to find which I basically took the diode as offering no resistance). (2) I had to UNDER-RELAX the update step. (This is a classic technique in non-linear problems such as computational fluid mechanics.)

The solver that I have used here is a very cheap one (2x2 matrices only). You can use your library one for bigger problems.

#include <iostream>
#include <vector>
#include <cmath>
using namespace std;

void solve( vector<vector<double>> A, vector<double> &x, const vector<double> &rhs );

int main()
{
   // Convergence
   const int MAXITER = 1000;
   const double TOLERANCE = 1.0e-8;
   const double UNDER_RELAX = 0.1;

   // Resistors
   const double R1 = 10.0, R2 = 20.0, R3 = 30.0; // resistances in ohms

   // Diode
   const double Is = 2.52e-9;                    // saturation current in A (weird value, though)
   const double Vt = 0.026;                      // thermal potential in V

   // Fixed voltages
   const double EMF = 5;                         // fixed by the battery
   const double EARTH = 0.0;                     //

   // Initial guess for free parameters
   int N = 2;
   double V2 = 25.0 / 11.0;                      // potential at node 2 in V
   double V3 = V2;                               // potential at node 3 in V

   double I1, I2, I3;                            // currents through resistors
   double Id;                                    // current for diode
   vector<double> f( N );                        // net current out of node
   vector<vector<double>> dfdx( N, vector<double>(N) );   // derivatives of f

   // Solve f(x) = 0   by  f -> f(x*) + dfdx.dx = 0   (matrix Newton-Raphson)
   bool converged = false;
   for ( int iter = 1; iter <= MAXITER; iter++ )
   {
      // Currents through resistors
      I1 = ( EMF - V2   ) / R1;
      I2 = ( V2 - EARTH ) / R2;
      I3 = ( V3 - EARTH ) / R3;

      // Current through diode
      double expfunc = exp( ( V2 - V3 ) / Vt );
      Id = Is * ( expfunc - 1.0 );     // current from diode

      // Net current from nodes 2, 3
      f[0] = -I1 + I2 + Id;
      f[1] =       I3 - Id;

      // Check convergence
      double normsq = 0.0;
      for ( int i = 0; i < N; i++ ) normsq += f[i] * f[i];
      double error = sqrt( normsq );
      converged = error < TOLERANCE;
      if ( converged ) break;

      // Derivative of f with respect to parameters V2 and V3
      dfdx[0][0] = 1.0 / R1 + 1.0 / R2 + ( Is / Vt ) * expfunc;
      dfdx[0][1] =                     - ( Is / Vt ) * expfunc;
      dfdx[1][0] =                     - ( Is / Vt ) * expfunc;
      dfdx[1][1] = 1.0 / R3            + ( Is / Vt ) * expfunc;

      vector<double> rhs(N);                     
      for ( int i = 0; i < N; i++ ) rhs[i] = -f[i];

      // Solve for increments in parameters
      vector<double> x(N);              
      solve( dfdx, x, rhs );

      // Update parameters
      V2 += UNDER_RELAX * x[0];
      V3 += UNDER_RELAX * x[1];

      cout << "iter: " << iter << "    V2: " << V2 << "   V3: " << V3 << "    error: " << error << '\n';
   }

   if ( !converged ) cout << "Not converged\n";
   cout << "V2: " << V2      << '\n'
        << "V3: " << V3      << '\n'
        << "Vd: " << V2 - V3 << '\n'
        << "I1: " << I1      << '\n'
        << "I2: " << I2      << '\n'
        << "I3: " << I3      << '\n'
        << "Id: " << Id      << '\n';
}


void solve( vector<vector<double>> A, vector<double> &x, const vector<double> &b )
// !!!!!!! This solves a 2x2 system only - replace with something better !!!!!!!!!!!!!!!!!
{
   double det = A[0][0] * A[1][1] - A[0][1] * A[1][0];
   x[0] = (  A[1][1] * b[0] - A[0][1] * b[1] ) / det;
   x[1] = ( -A[1][0] * b[0] + A[0][0] * b[1] ) / det;
}

Final iterations and output:

.......
.......
iter: 149    V2: 2.80885   V3: 2.36018    error: 1.34202e-08
iter: 150    V2: 2.80885   V3: 2.36018    error: 1.20782e-08
iter: 151    V2: 2.80885   V3: 2.36018    error: 1.08704e-08
V2: 2.80885
V3: 2.36018
Vd: 0.44867
I1: 0.219115
I2: 0.140442
I3: 0.0786726
Id: 0.0786726