牛顿法潮流计算直角坐标系雅各比矩阵推导

用牛顿拉夫逊法线性化了的潮流功率方程组可表达为：

$\Delta f=J\Delta x$

其中 $J$ 称为函数 $f$ 的雅各比矩阵。在各种电力系统分析的教材中，雅各比矩阵的形成方法都是根据功率方程: $S=UY^{*}U^{*}$

将每个节点的注入功率表达式求出，然后再分成P,Q两个分量的函数表达式，再对节点电压的实部和虚部或幅值和相角求偏导数，并得出一长串复杂的偏导数表达式。有时候还需要考察对这个雅各比矩阵的记忆情况，那么多元素记忆得死去活来。

但是，这部分内容应该是计算机潮流计算的介绍，上面那种推导方法是否适合计算机计算呢？答案肯定是否定的，不仅没有性能优势，还造成了巨大的编程压力。要知道现在大部分编程环境都可以直接进行复数计算，支持复数类型了。即使没有支持，也可以自己建立复数类型来计算。

但这还不是主要的，书上所描述的雅各比矩阵推导过程和计算方法更不是计算机所采用的主流方法。现在的线性代数计算库十分健全，完全可以用矩阵的计算来得出雅各比矩阵。然而书上的作法相当于把矩阵运算的一部分手动完成了，变成了极为繁杂的公式，而且这种繁琐的办法并没有减轻计算负担。而得出的雅各比矩阵存在着2x2子矩阵的分块对称性，一方面是原本是复数的偏导数被分为实部和虚部，复功率被表达为有功和无功所以从两个复数变为4个实数，另一方面，是因为导纳矩阵为对称阵，因此分块对称。

下面，将以直角坐标系下的潮流计算为例，推导潮流计算雅各比矩阵的矩阵表达形式：

首先要介绍一些向量求导的知识:为了使向量 $A$ 和向量 $B$ 相乘得到相同维度的向量 $C$ （即 $c_i=a_i*b_i$ ),可以采用以下形式：
$C = [A]B = [B]A$
然后，向量的导数就可以用下面公式表达：
$\frac{dC}{dX} =C_{X}= [A]\frac{dB}{dX}+ [B]\frac{dA}{dX}$

" $[]$ "是取向量的对角化矩阵。
设 $X$ 是潮流计算所要求的实数解向量，即:

$X=[e,f]^{T}$

为了方便起见，先用一个全部为PQ节点的系统举例。

设 $U$ 是一个 $n_{B}\times 1$ 的向量,表示各节点电压，其中每个元素 $u_{i}=e_{i}+f_{i}j$ ，其对向量 $e,f$ 的偏导为：
$\frac{\partial U}{\partial e}=[E]$
$\frac{\partial U}{\partial f}=j[E]$

其中 $[E]$ 为 $n_{B}\times n_{B}$ 的单位矩阵。之所以这个向量对向量的偏导是一个矩阵，是因为其中每个元素 $u_{i}$ 都要对 $e$ 这个向量求偏导，这样才符合向量求导的定义。

潮流计算中，复功率可表示为

$S_{bus}=[U]I^{*}_{bus}$

其中 $I_{bus}=Y_{bus}U$

$\frac{\partial I_{bus}}{\partial e}=Y_{bus} \frac{\partial U}{\partial e}=Y_{bus} [E]$

$\frac{\partial I_{bus}}{\partial f}=Y_{bus} \frac{\partial U}{\partial f}=jY_{bus} [E]$

雅各比矩阵其实就是

$J=\frac{dS_{bus}}{dX}=[\frac{\partial S}{\partial e},\frac{\partial S}{\partial f}]$
其中
$\frac{\partial S}{\partial e}=[V]\frac{\partial I_{bus}^{*}}{\partial e}+[I_{bus}^{*}]\frac{\partial U}{\partial e}$

$=[V]Y_{bus}^{*}+I_{bus}^{*}$

$\frac{\partial S}{\partial f}=[V]\frac{\partial I_{bus}^{*}}{\partial f}+[I_{bus}^{*}]\frac{\partial U}{\partial f}$
$=j([I_{bus}^{*}]-Y_{bus}^{*}[V])$

下面是程序验证，用书上的例子进行对比：

import numpy as np

# 陈珩《电力系统稳态分析》教材例4-3
Y=np.array([[6.25-18.75j,-5+15j,-1.25+3.75j,0,0],
           [-5+15j,10.834-32.5j,-1.667+5j,-1.667+5j,-2.5+7.5j],
           [-1.25+3.75j,-1.667+5j,12.917-38.75j,-10+30j,0],
           [0,-1.667+5j,-10+30j,12.917-38.75j,-1.25+3.75j],
           [0,-2.5+7.5j,0,-1.25+3.75j,3.75-11.25j]])

U=np.array([1.06,1,1,1,1])
#注入功率计算式
Ibus=Y.dot(U)
P0=U*np.conj(Ibus)
P1=np.diag(np.conj(U))*Ibus

#N,L
J1=np.diag(U).dot((Y.conj()))+np.diag(Ibus.conj())
#H,J
J2=1j*(np.conj(np.diag(Ibus))-np.conj(Y).dot(np.diag(U)))

L=(len(J1)-1)*2
J=np.zeros((L,L))

J[0::2,1::2]=J1[1:,1:].real
J[::2,::2]=J2[1:,1:].real
J[1::2,1::2]=J1[1:,1:].imag
J[1::2,::2]=J2[1:,1:].imag

print(J)

[[ 33.4    10.534  -5.     -1.667  -5.     -1.667  -7.5    -2.5  ]
 [-11.134  31.6     1.667  -5.      1.667  -5.      2.5    -7.5  ]
 [ -5.     -1.667  38.975  12.842 -30.    -10.      0.      0.   ]
 [  1.667  -5.    -12.992  38.525  10.    -30.      0.      0.   ]
 [ -5.     -1.667 -30.    -10.     38.75   12.917  -3.75   -1.25 ]
 [  1.667  -5.     10.    -30.    -12.917  38.75    1.25   -3.75 ]
 [ -7.5    -2.5     0.      0.     -3.75   -1.25   11.25    3.75 ]
 [  2.5    -7.5     0.      0.      1.25   -3.75   -3.75   11.25 ]]

可见结果与教材完全一致。

显然对于PV节点由于方程形式不同，应变量由单纯的复功率变为了有功功率和电压的平方，其偏导数也不同，但我们可以知道，有功功率的偏导和上面的结论一致，而电压平方的偏导很容易求出，即:

$U_{i}^{2}=e_{i}^{2}+f_{i}^{2}$

$\frac{\partial U^{2}}{\partial e}=2[e]$

$\frac{\partial U^{2}}{\partial f}=2j[f]$

与教材结论完全一致。

参考文献MatPower Maunal