电力系统潮流计算的计算机算法（二）——牛顿-拉夫逊法的原理

本篇为本科课程《电力系统稳态分析》的笔记。

本篇为这一章的第二篇笔记。上一篇传送门，下一篇传送门。

潮流计算的牛顿-拉夫逊方法

牛顿-拉夫逊法的原理

牛顿-拉夫逊法是求解非线性方程的最好方法之一。先从一维非线性方程式的解来说明他的意义和推导过程。

牛顿-拉夫逊法的意义和推导

设一维方程为：$f(x)=0$。假设其准确解为$x^{(\ast )}$，而$x^{(k)}$是其近似解，他和准确解之间的差值为$\Delta x^{(k)}$，所以有：
$$x^{(\ast )}=x^{(k)}+\Delta x^{(k)}$$

将上式带入方程中，可得：
$$f(x^{(\ast )})=f(x^{(k)}+\Delta x^{(k)})=0$$

将上式在$x^{(k)}$处进行泰勒展开，得到：
$$f(x^{(k)}+\Delta x^{(k)})=f(x^{(k)})+f^{\prime }(x^{(k)})\Delta x^{(k)}+\cdots =0$$

略去上面展开的高次项，可得到线性部分：
$$f(x^{(k)})+f^{\prime }(x^{(k)})\Delta x^{(k)}=0$$

这一个方程就是修正方程，利用他可以解出：
$$\Delta x^{(k)}=-\frac{f(x^{(k)})}{f^{\prime }(x^{(k)})}$$

而$\Delta x^{(k)}$被称为修正量。这样就可以得到迭代下一回的解，它有望更加接近于准确解：
$$x^{(k+1)}=\Delta x^{(k)}+x^{(k)}$$

然后通过相同的接法求出$\Delta x^{(k+1)}$，得到有望更加接近于准确解的$x^{(k+2)}$。

上面就是牛顿拉夫逊法的步骤，这样继续迭代下去，直到满足$|\Delta x^{(k)}|\leqslant \varepsilon $ （精度），所求出来的解就是最终结果。

上面的迭代过程可以整理成下面的迭代形式：
$$\begin{gathered} \Delta x^{(k)}=-\frac{f(x^{(k)})}{f^{\prime }(x^{(k)})} \\ {x}^{(k+1)}=\Delta x^{(k)}+x^{(k)} \\ (k=0,1,\cdots ) \end{gathered}$$
初值$x^{(0)}$需要用户给定。

牛拉法的特点

1. 它是迭代法，是逐步逼近的方法。
2. 修正方程式线性化的方程，体现在把非线性方程在$x^{(k)}$处进行泰勒展开，并略去了高阶小量。
3. 它对于初值是敏感的，如果初始值远离准确解，会导致不收敛。

多变量非线性方程的解

设n维非线性方程组为：
$$\{\begin{array}{l}{f}_1\left(x_1,x_2,\cdots ,x_n\right)=0\\ f_2\left(x_1,x_2,\cdots ,x_n\right)=0\\ ⋮\\ f_n\left(x_1,x_2,\cdots ,x_n\right)=0\end{array}$$

假设初始值给定为$x_1^{(0)},x_2^{(0)},...,x_n^{(0)}$，并令各变量的修正量为$\Delta x_1^{(0)},\Delta x_2^{(0)},\cdots ,\Delta x_n^{(0)}$，那么对于上面的n个方程，在初始值附近展开成泰勒级数，并忽略修正量二次和更高次数的项，得到：
$$\begin{gathered} f_1\left(x_1^{(0)},x_2^{(0)},\cdots x_n^{(0)}\right)+\left(\frac{\partial f_1}{\partial x_1}{|}_0\Delta x_1^{(0)}+\frac{\partial f_1}{\partial x_2}{|}_0\Delta x_2^{(0)}+\cdots +\frac{\partial f_1}{\partial x_n}{|}_0\Delta x_n^{(0)}\right)=0 \\ {f}_2\left(x_1^{(0)},x_2^{(0)},\cdots x_n^{(0)}\right)+\left(\frac{\partial f_2}{\partial x_1}{|}_0\Delta x_1^{(0)}+\frac{\partial f_2}{\partial x_2}{|}_0\Delta x_2^{(0)}+\cdots +\frac{\partial f_2}{\partial x_n}{|}_0\Delta x_n^{(0)}\right)=0 \\ ⋮ \\ {f}_n\left(x_1^{(0)},x_2^{(0)},\cdots x_n^{(0)}\right)+\left(\frac{\partial f_n}{\partial x_1}{|}_0\Delta x_1^{(0)}+\frac{\partial f_n}{\partial x_2}{|}_0\Delta x_2^{(0)}+\cdots +\frac{\partial f_n}{\partial x_n}{|}_0\Delta x_n^{(0)}\right)=0 \end{gathered}$$

其中，$\frac{\partial f_i}{\partial x_j}{|}_0$是函数$f_i(x_1,x_2,\cdots ,x_n)$对自变量$x_j$的偏导数在$(x_1^{(0)},x_2^{(0)},\cdots ,x_n^{(0)})$处的值。

将上方程组表示为矩阵形式，可得：
$$\left[\begin{array}{c}{f}_1(x_1^{(0)},x_2^{(0)},\cdots ,x_n^{(0)})\\ f_1(x_1^{(0)},x_2^{(0)},\cdots ,x_n^{(0)})\\ ⋮\\ f_n(x_1^{(0)},x_2^{(0)},\cdots ,x_n^{(0)})\end{array}\right]=-\left[\begin{array}{ccc}\frac{\partial f_1}{\partial x_1}{|}_0\frac{\partial f_1}{\partial x_2}{|}_0& ...& \frac{\partial f_1}{\partial x_n}{|}_0\\ \frac{\partial f_2}{\partial x_1}{|}_0\frac{\partial f_2}{\partial x_2}{|}_0& ...& \frac{\partial f_2}{\partial x_n}{|}_0\\ ⋮& & \\ \frac{\partial f_n}{\partial x_1}{|}_0\frac{\partial f_n}{\partial x_2}{|}_0& ...& \frac{\partial f_n}{\partial x_n}{|}_0\end{array}\right]\left[\begin{array}{c}\Delta x_1^{(0)}\\ \Delta x_2^{(0)}\\ ⋮\\ \Delta x_n^{(0)}\end{array}\right]$$

使用向量$\mathbf{x},{\mathbf{x}}^{(0)},\Delta {\mathbf{x}}^{(0)}$分别代表变量、变量的初值和修正量所组成的向量，则将上面的矩阵可以化为：
$${\mathbf{J}}^{(0)}\Delta {\mathbf{x}}^{(0)}=-\mathbf{f}({\mathbf{x}}^{(0)})$$

其中，$\mathbf{J}$是雅可比矩阵，为向量$\mathbf{f}(\mathbf{x})$对变量$\mathbf{x}$的一阶偏导数矩阵。

上面这个方程就是修正量的线性方程。从上式中可以解出$\Delta x_1^{(0)},\Delta x_2^{(0)},\cdots ,\Delta x_n^{(0)}$，然后就可以计算出修正后的值：
$$\{\begin{array}{l}{x}_1^{(1)}=x_1^{(0)}+\Delta x_1^{(0)}\\ x_2^{(1)}=x_2^{(0)}+\Delta x_2^{(0)}\\ ⋮\\ x_n^{(1)}=x_n^{(0)}+\Delta x_n^{(0)}\end{array}$$

然后将求出来的$x_1^{(1)},x_2^{(1)},...,x_n^{(1)}$作为新的初值，带入式子重新形成并求解新的修正方程式，这样就可以得到新的修正量$\Delta x_1^{(1)},\Delta x_2^{(1)},\cdots ,\Delta x_n^{(1)}$，然后按照同样的步骤继续。

那么整个过程可以总结为，给定${\mathbf{x}}^{(0)}$：
$$\begin{gathered} {\mathbf{J}}^{(k)}\Delta {\mathbf{x}}^{(k)}=-\mathbf{f}({\mathbf{x}}^{(k)}) \\ {\mathbf{x}}^{(k+1)}={\mathbf{x}}^{(k)}+\Delta {\mathbf{x}}^{(0)} \\ (k=0,1,\cdots ) \end{gathered}$$

迭代的收敛判断依据为：
$$||\mathbf{f}({\mathbf{x}}^{(k)})|{|}_{\infty }=\underset{i}{\max }|f_i({\mathbf{x}}^{(k)})|⩽\epsilon$$

其中$||\mathbf{f}({\mathbf{x}}^{(k)})|{|}_{\infty }$代表的是向量$\mathbf{f}({\mathbf{x}}^{(k)})$无穷范数，即$\mathbf{f}({\mathbf{x}}^{(k)})$中绝对值最大的那个元素，$\epsilon$为所给定的容许误差。

牛拉法主要计算量在于修正方程的求解，下面说明一般线性方程组的基本求解方法。对于n阶线性方程组$\mathbf{A}\mathbf{x}=\mathbf{b}$，通过高斯消去，并记录消去过程中的因子，可以将系数矩阵$\mathbf{A}$分解成两个三角矩阵的乘积，即：
$$\mathbf{A}=\mathbf{L}\mathbf{U}$$

$$L=\left[\begin{array}{cccccc}1& & & & & \\ l_{21}& 1& & & & \\ ⋮& \ddots & \ddots & & & \\ l_{n-1,1}& \cdots & l_{n-1,n-2}& 1& & \\ l_{n,1}& l_{n,2}& \cdots & l_{n,n-1}& 1& \end{array}\right];U=\left[\begin{array}{ccccc}{u}_{11}& u_{12}& \cdots & u_{1,n-1,}& u_{1,n}\\ & u_{22}& \cdots & u_{2,n-1}& u_{2,n}\\ & & \ddots & \ddots & ⋮\\ & & & u_{n-1,n-1}& u_{n-1,n}\\ & & & & u_{nn}\end{array}\right]$$

其中，$\mathbf{L}$是单位下三角矩阵，$\mathbf{U}$是上三角矩阵。

该分解为Doolittle分解，其分解过程为：
$$\begin{gathered} u_{ij}=a_{ij}-\sum _{k=1}^{i-1}{l}_k{u}_{kj}(j=i,\cdots ,n) \\ {l}_{ji}=\left(a_{ji}-\sum _{k=1}^{i-1}{l}_{jk}{u}_{ki}\right)/u_{ii}(j=i+1,\cdots ,n) \\ (i=1,\cdots ,n) \end{gathered}$$

这个分解过程需要乘除的运算次数为$\frac{n^3-n}{3}$。

利用上面的分解，线性方程组的求解就变成了对两个三角方程的相继求解：
$$\begin{gathered} \mathbf{L}\mathbf{y}=\mathbf{b} \\ \mathbf{U}\mathbf{x}=\mathbf{y} \end{gathered}$$

上面第一个式子的求解称为前代，其求解过程为：
$$\begin{gathered} y_1=b_1 \\ {y}_i=b_i-\sum _{j=1}^{i-1}{l}_{ij}{y}_j \\ (i=2,\cdots ,n) \end{gathered}$$

上面第二个式子的求解称为回代，其求解过程为：
$$\begin{gathered} x_n=y_n/u_m \\ {x}_i=\left(y_i-\sum _{j=i+1}^n{u}_{ij}{x}_j\right)/u_i \\ (i=n-1,\cdots ,1) \end{gathered}$$

前代和回代总共需要乘除次数为$n^2$。

当$\mathbf{A}$为对称矩阵的时候，其三角分解师可以简化为：$\mathbf{A}=\mathbf{L}\mathbf{D}{\mathbf{L}}^T$。其中，$\mathbf{L}$是单位下三角矩阵，$\mathbf{D}$是对角矩阵， D = d i a g { d 1 , d 2 , ⋯   , d n } \boldsymbol{D}=\mathrm{diag}\{d_1,d_2,\cdots,d_n\} D=diag{d1​,d2​,⋯,dn​}。

该分解为Choleskey分解，其分解过程为：
$$\begin{gathered} d_i=a_i-\sum _{k=1}^{i-1}{l}_k^2{d}_k \\ {l}_{ji}=\left(a_{ji}-\sum _{k=1}^{i-1}{l}_k{d}_k{l}_{ji}\right)/d_i\left(j=i+1,\cdots ,n\right) \\ (i=1,\cdots ,n) \end{gathered}$$

利用上面的分解，线性方程组的求解就变成了如下两个三角方程的相继求解：
$$\begin{gathered} \mathbf{L}\mathbf{y}=\mathbf{b} \\ {\mathbf{L}}^T\mathbf{x}={\mathbf{D}}^{-1}\mathbf{y} \end{gathered}$$

上面第一个式子的求解前面已经介绍过，为前代。上面第二个式子的求解过程为：
$$\begin{gathered} x_n=y_n/d_n \\ {x}_i=y_i/d_i-\sum _{j=i+1}^n{l}_{ji}{x}_j\left(i=n-1,\cdots ,1\right) \end{gathered}$$