基于Newton-Raphson法的含二极管非线性电路仿真实现问题排查

我正在开发一款基于**Modified Nodal Analysis（MNA）**的电路仿真器，目标是支持二极管这类非线性元件。为求解非线性方程，我实现了Newton-Raphson方法，但当前输出结果与其他仿真器差异极大，需要排查问题。

实现代码 #

// 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;
}

仿真结果对比 #

• 其他两款仿真器结果：
  • 第一款：Vd = 739 mV，Id = 70.74 mA
  • 第二款：Vd = 765 mV，Id = 70 mA
• 我的实现输出：Vd = 175 mV，Id = 2.14 µA

补充推导与验证 #

基于基尔霍夫电压定律（KVL）和基尔霍夫电流定律（KCL）推导的节点方程如下：

KCL at node 1: I + (v1 - v2)/R1 = 0 // I是电压源支路电流
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

我的代码运行得到的解向量：

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

该解仅满足第一和最后一个方程，无法匹配所有方程。使用专业仿真工具得到的正确解向量：

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

此解代入所有方程均成立，与推导的方程组完全匹配。

核心排查方向 #

1. MNA方程组与残差向量构建错误 #

代码中残差向量r1/r2/r3的逻辑完全不符合KCL的物理意义：

• 残差的本质是节点电流不平衡量，即所有流入节点的电流之和减去流出节点的电流之和，当前计算逻辑完全偏离这一规则。
• 未正确处理电压源的约束：节点1电压应固定为vs=5V，不应作为迭代变量；同时MNA需引入电压源电流作为额外变量，当前3变量模型维度不足，无法完整描述电路。

2. 雅可比矩阵元素错误 #

雅可比矩阵的每个元素是残差对对应节点电压的偏导数，当前矩阵存在多处错误：

• 节点2的残差包含v2/R2项，其偏导数应为1/R2，但当前j22中仅存在-1/r2 -g，符号和项数均错误。
• 二极管导纳g的位置与符号不符合MNA的互导纳规则，自导纳和互导纳的增减逻辑错误。

3. 变量更新与线性求解错误 #

Newton-Raphson的更新公式为x_new = x_old - J^{-1} * R(x_old)，需确认：

• 克莱姆法则求解线性方程组时的符号是否正确，当前dv1/dv2/dv3的计算可能存在符号偏差。
• 变量更新时v1 -= dv1的符号是否匹配求解结果的逻辑。

4. 初始值设置错误 #

节点1电压由电压源直接约束，初始值不应设为0.0，这会导致迭代方向完全错误，甚至无法收敛到正确解空间。

内容的提问来源于stack exchange，提问作者Abdo21