重温牛顿-拉夫逊法

牛顿-拉夫逊法是数值计算中最经典的方法之一。对电力系统工程师而言，它既是潮流计算的基础，也是很多参数辨识、最优化和非线性方程求解程序的核心工具。本文分三部分讨论：先做一点历史考据；再用标量与向量两种形式重述公式推导；最后给一个简洁而工程上更稳妥的 C++ 例子，并讨论牛顿法与潮流计算之间的关系。

1 历史考据

牛顿法最早可以追溯到艾萨克·牛顿在《流数法》中的工作。牛顿当年并没有直接写成现代教科书中那种整齐的迭代公式，而是通过具体代数问题不断修正近似值。后来的拉夫逊（Joseph Raphson）在 1690 年出版的著作中，把这种思路写得更代数化、更便于重复应用，因此后人通常把该方法称为牛顿-拉夫逊法。

这个命名其实很合理：思想源头主要在牛顿，算法表达更接近现代形式则有赖于拉夫逊的整理。

2 从切线法到一维牛顿迭代

先看最简单的一维情形。要求解

若当前迭代点为 $x_k$，把 $f(x)$ 在 $x_k$ 附近做一阶泰勒展开：

牛顿法的核心思想是：不用真的解原方程，而是去解当前点的线性近似方程。 令右端近似为 0，可得

因此

这个公式有非常明确的几何意义：以 $(x_k,f(x_k))$ 处的切线代替原函数，切线与 $x$ 轴的交点作为下一次迭代点。

2.1 它为什么收敛得快

若初值足够接近真解，且函数在真解附近足够光滑，且 $f’(x^\star)\neq 0$，则牛顿法具有局部二次收敛性。所谓“二次”，是指误差大致满足

这意味着一旦进入真解附近，正确数字位数会迅速翻倍。这正是牛顿法在工程计算中极具吸引力的根本原因。

2.2 它为什么有时会失败

牛顿法并不总是“必胜”。典型失败原因有：

1. 初值离真解太远；
2. 导数接近 0，导致步长异常放大；
3. 函数存在多个根，算法被吸引到别的根；
4. 线性化在当前点附近代表性太差。

所以，牛顿法是一种局部强方法，而不是全局鲁棒方法。工程实现中，常常需要配合阻尼、步长限制、信赖域或良好的初值。

3 多维情形：真正适合电力系统的写法

对非线性方程组

其中

在当前点 $x_k$ 处一阶线性化：

这里

是雅可比矩阵。令线性化后的右端为 0，得到牛顿校正方程：

这才是工程上最推荐的写法。需要特别强调的是：程序中不应显式求逆 $J^{-1}$，而应直接求解线性方程组。 原因有三点：

1. 求逆的计算量更大；
2. 数值误差通常更差；
3. 对稀疏矩阵尤其不合适。

电力系统潮流程序、状态估计程序、最优潮流程序，几乎都遵循这个思路。

4 一个二维非线性方程组例子

考虑方程组：

写成向量形式：

雅可比矩阵为：

于是每一步都要求解

然后更新

5 一个更稳妥的 C++ 实现

下面给出一个简单的 C++ 例子。与很多教材不同，这里不显式写逆矩阵，而是直接解 $2\times 2$ 线性方程组。对小系统而言，这样已经足够清楚；对大系统，则会进一步替换为稀疏矩阵分解。

这段程序的工程意义主要有三点：

1. 先判断残差，再判断步长；
2. 对雅可比行列式过小的情况显式报错；
3. 不显式求逆矩阵。

#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>

struct State {
    double x;
    double y;
};

double f1(double x, double y) {
    return std::sin(x) + std::exp(y) - 1.0;
}

double f2(double x, double y) {
    return x + std::cos(y) - 1.0;
}

bool newton_step(State& s, double tol = 1e-12) {
    const double J11 = std::cos(s.x);
    const double J12 = std::exp(s.y);
    const double J21 = 1.0;
    const double J22 = -std::sin(s.y);

    const double rhs1 = -f1(s.x, s.y);
    const double rhs2 = -f2(s.x, s.y);

    const double det = J11 * J22 - J12 * J21;
    if (std::abs(det) < tol) {
        throw std::runtime_error("Jacobian is singular or nearly singular.");
    }

    // 解 J * delta = rhs
    const double dx = ( rhs1 * J22 - J12 * rhs2) / det;
    const double dy = (-rhs1 * J21 + J11 * rhs2) / det;

    s.x += dx;
    s.y += dy;

    return std::sqrt(dx * dx + dy * dy) < 1e-10;
}

int main() {
    State s{0.1, -0.1};  // 初值可根据经验或图形分析给定

    std::cout << std::fixed << std::setprecision(12);

    for (int k = 0; k < 20; ++k) {
        const double r1 = f1(s.x, s.y);
        const double r2 = f2(s.x, s.y);

        std::cout << "iter " << std::setw(2) << k
                  << "  x = " << s.x
                  << "  y = " << s.y
                  << "  ||F|| = " << std::sqrt(r1 * r1 + r2 * r2)
                  << '\n';

        if (std::sqrt(r1 * r1 + r2 * r2) < 1e-12) {
            std::cout << "Converged.\n";
            return 0;
        }

        if (newton_step(s)) {
            std::cout << "Converged by step size.\n";
            return 0;
        }
    }

    std::cout << "Iteration limit reached.\n";
    return 0;
}

6 牛顿法与潮流计算的关系

电力系统潮流计算的本质，就是把节点有功、无功不平衡方程写成一个非线性方程组，然后用牛顿法反复消除不平衡量。

若把未知量写成电压相角与电压幅值的组合，便可形成标准潮流方程。每一步迭代都做两件事：

1. 根据当前状态计算功率不平衡向量；
2. 用雅可比矩阵求出状态修正量。

从算法骨架上说，潮流计算与前述二维例子没有本质区别；区别只在于：

• 维数更高；
• 矩阵更稀疏；
• 变量分块更明显；
• 需要处理 PQ、PV、平衡节点以及无功限值等工程约束。

也正因为如此，真正工业级的潮流程序，关键不只在“会写牛顿法”，而在于能否高效地构造雅可比矩阵、调用稀疏因子分解、在 PV/PQ 切换和不良初值下仍保持鲁棒。

7 牛顿法与梯度法的关系

牛顿法还可以从最优化角度理解。对一维函数 $f(x)$，若希望最小化它，在当前点二阶展开：

令这个二次近似的导数为零，可得

这就是一维牛顿法。因此，从某种意义上说，牛顿法是在当前点上利用了二阶信息的“最优局部步长”方法。相比只利用一阶信息的梯度法，它在局部通常更快，但也更依赖 Hessian 或 Jacobian 的质量。

8 结语

牛顿-拉夫逊法之所以经典，不是因为它公式漂亮，而是因为它抓住了一个极其有效的思想：复杂的非线性问题，不是直接硬解，而是不断在当前位置求解一个线性近似问题。

这套思想对于电力系统工程师非常重要。潮流、状态估计、参数辨识、最优潮流，很多看似不同的问题，在数值计算层面都共享同样的骨架。真正理解牛顿法，往往不是会背公式，而是理解三件事：

1. 线性化到底近似了什么；
2. Jacobian 为什么决定收敛性；
3. 工程程序为什么从不真正去算逆矩阵。