矩阵思维下的潮流计算
在知识星球“深入理解电力系统”,前段时间有同仁提问:”如何看待矩阵理论在电力系统中的应用?是让问题形式更简洁、更干净?还是说确实能从矩阵的角度发现不同的解决思路?”。这个问题非常有趣,我当时也做了较为详细的回答。但考虑到矩阵思维在电力系统课程和工程实践中仍未充分普及,本文以潮流计算为例,进一步展示基于矩阵思维的潮流计算推导和软件实现思路。
1. 导读
首次看本文的推导,许多人可能会有一种明显的不适应:熟悉的
\[P_i=\sum_j m_im_j\left(G_{ij}\cos\theta_{ij}+B_{ij}\sin\theta_{ij}\right)\]和
\[Q_i=\sum_j m_im_j\left(G_{ij}\sin\theta_{ij}-B_{ij}\cos\theta_{ij}\right)\]似乎突然退到了幕后,取而代之的是关联矩阵、节点导纳矩阵、复电压向量以及
\[I=YV, \qquad S=[V](YV)^*.\]随后,传统教材里需要分节点、分对角元与非对角元、分有功与无功反复计算的雅可比导数,在这里被压缩为极短的复数矩阵微分公式。表面上看,这像是在学习另一种潮流算法。但实际上,这并不是新的潮流算法,物理模型、未知量和最终求解结果都没有改变。真正改变的是:我们整体思考问题的思路,以及这种思路可用于编写现代高性能潮流计算软件。
2. 总体思路
传统教材通常先选定第 $i$ 个母线,将节点电流写成各相邻母线电压的加权和,再把复电压展开成幅值与相角,最后分别取实部和虚部,得到 $P_i$ 与 $Q_i$。为了构造牛顿法,还要进一步计算
\[\frac{\partial P_i}{\partial\theta_j}, \quad \frac{\partial P_i}{\partial m_j}, \quad \frac{\partial Q_i}{\partial\theta_j}, \quad \frac{\partial Q_i}{\partial m_j}.\]由于 $i=j$ 与 $i\neq j$ 时公式形式不同,教材往往还要分别列出对角元素和非对角元素。学生最终会看到许多带有 $G_{ij}$、$B_{ij}$、$\sin\theta_{ij}$ 和 $\cos\theta_{ij}$ 的表达式,并需要反复核对正负号。
这种路线并没有错误。它的优点是每一个式子都能落实到单个母线,适合小系统手算,也便于直接观察某条支路参数如何进入某个节点方程。但它有可能存在教学上的副作用:当所有公式都被完全展开以后,网络的共同结构被隐藏了。学生容易误以为潮流计算的本质就是三角函数求和,甚至把雅可比矩阵理解成四张需要单独背诵的公式表。
事实上,那些三角函数并不是潮流模型最原始的物理内容。它们只是复电压采用极坐标表示后,复数乘法被写成实数坐标时自然出现的结果。
矩阵方法把观察顺序倒了过来。它首先看的是问题的整体结构,例如:
- 支路怎样连接母线?
- 支路电压怎样产生支路电流?
- 支路电流怎样通过节点电流定律汇集为母线电流?
- 母线电压和母线电流怎样组成复功率?
- 当电压幅值和相角发生微小变化时,上述映射怎样整体变化?
回答这些问题,只需要依次使用四类关系:
\[\boxed{ \text{拓扑关系} \longrightarrow \text{元件伏安关系} \longrightarrow \text{节点电流关系} \longrightarrow \text{复功率关系} }\]3. 从支路拓扑得到节点导纳矩阵
3.1 简单串联支路
先考虑只有串联导纳、没有变压器和线路对地电纳的简单网络。为每条支路任意指定一个方向,定义有向关联矩阵
\[D\in\mathbb R^{\ell\times n}.\]若第 $k$ 条支路从母线 $i$ 指向母线 $j$,则该行满足
\[D_{ki}=1, \qquad D_{kj}=-1,\]其余元素为零。
于是支路两端电压差为
\[v_{\mathrm{br}}=DV.\]设各支路串联导纳组成向量
\[y=(y_1,y_2,\ldots,y_\ell)^T,\]则按所选方向定义的支路电流为
\[i_{\mathrm{br}}=[y]DV.\]利用节点电流定律,将支路电流汇集回母线:
\[I=D^T i_{\mathrm{br}}.\]代入支路电流可得
\[I=D^T[y]DV.\]因此
\[\boxed{Y=D^T[y]D.}\]若再考虑母线并联导纳 $Y_{\mathrm{sh}}$,则
\[\boxed{Y=D^T[y]D+Y_{\mathrm{sh}}.}\]支路方向是人为选择的。若把某条支路方向反转,$D$ 的相应行会整体变号,但 $D^T[y]D$ 不变,因此物理结果与方向选择无关。
3.2 实际线路的两端口形式
实际交流支路常含线路充电电纳、非标准变比和相移。此时,比 $D^T[y]D$ 更通用的做法是先写每条支路的两端口关系:
\[\begin{bmatrix} I_{f,k}\\ I_{t,k} \end{bmatrix} = \begin{bmatrix} y_{ff,k}&y_{ft,k}\\ y_{tf,k}&y_{tt,k} \end{bmatrix} \begin{bmatrix} V_{f,k}\\ V_{t,k} \end{bmatrix}.\]其中 $f$ 表示支路首端,$t$ 表示支路末端。
对普通 $\pi$ 型线路,若串联导纳为 $y_s$,总充电电纳为 $jb_c$,则
\[y_{ff}=y_s+j\frac{b_c}{2}, \qquad y_{ft}=-y_s,\] \[y_{tf}=-y_s, \qquad y_{tt}=y_s+j\frac{b_c}{2}.\]若首端还有复变比
\[t=\tau e^{j\phi},\]则常用原始导纳系数为
\[y_{ff}=\frac{y_s+jb_c/2}{|t|^2}, \qquad y_{ft}=-\frac{y_s}{t^*},\] \[y_{tf}=-\frac{y_s}{t}, \qquad y_{tt}=y_s+jb_c/2.\]相移变压器可能使 $y_{ft}\neq y_{tf}$,从而 $Y$ 不再是普通对称矩阵。后文推导并不要求 $Y=Y^T$。
3.3 全网矩阵形成
定义首端、末端连接矩阵
\[C_f,C_t\in\{0,1\}^{\ell\times n}.\]它们把母线电压抽取到支路两端:
\[V_f=C_fV, \qquad V_t=C_tV.\]将所有支路原始导纳系数排列成向量后,可构造
\[Y_f=[y_{ff}]C_f+[y_{ft}]C_t,\] \[Y_t=[y_{tf}]C_f+[y_{tt}]C_t.\]于是
\[I_f=Y_fV, \qquad I_t=Y_tV.\]节点电流由支路两端电流汇集得到:
\[\boxed{ Y=C_f^TY_f+C_t^TY_t+Y_{\mathrm{sh}}. }\]这说明,节点导纳矩阵不是一组孤立的节点公式,而是“网络拓扑”和“元件端口关系”装配后的全局线性算子。
4. 交流潮流的三个基本公式
统一矩阵推导只需要反复使用三个公式。
4.1 母线电流
\[\boxed{I=YV.}\]这是线性网络方程。给定母线电压后,所有母线注入电流一次矩阵乘法即可得到。
4.2 母线复功率
单个母线的复功率定义为
\[S_i=V_iI_i^*.\]将所有母线同时写成向量形式(本文的星号表示逐元素复共轭,不是共轭转置):
\[\boxed{S=[V]I^*=[V](YV)^*.}\]其中
\[S=P+jQ.\]因此
\[P=\operatorname{Re}(S), \qquad Q=\operatorname{Im}(S).\]网络方程 $I=YV$ 是线性的,潮流非线性来自功率定义中的乘积 $V\odot I^*$。
4.3 极坐标电压
\[\boxed{V=m\odot e^{j\theta}.}\]第 $i$ 个元素为
\[V_i=m_ie^{j\theta_i}.\]潮流求解通常以相角和幅值作为实状态变量,但在计算过程中仍将电压作为一个复向量整体保存。
5. 电压对相角和幅值的微分
这是整个推导中最关键、也最简单的一步。
5.1 单个母线
由
\[V_i=m_ie^{j\theta_i}\]可得
\[\frac{\partial V_i}{\partial\theta_i} =jm_ie^{j\theta_i} =jV_i,\]以及
\[\frac{\partial V_i}{\partial m_i} =e^{j\theta_i}.\]定义单位相量
\[u=V\oslash m=e^{j\theta}.\]于是
\[\frac{\partial V_i}{\partial m_i}=u_i.\]5.2 全部母线的矩阵形式
由于每个 $V_i$ 只直接依赖自己的 $\theta_i$ 和 $m_i$,全部母线的一阶微分为
\[\boxed{ dV=j[V]d\theta+[u]dm. }\]取共轭得到
\[\boxed{ dV^*=-j[V^*]d\theta+[u^*]dm. }\]这两个公式有直观的几何含义:
- 改变 $\theta_i$,会使相量沿圆周切向移动,方向为 $jV_i$;
- 改变 $m_i$,会使相量沿径向移动,方向为单位相量 $u_i$。
6. 节点复功率雅可比的整体推导
6.1 对复功率使用乘积法则
从
\[S=[V]I^*\]出发,一阶微分为
\[dS=[dV]I^*+[V]dI^*.\]利用对角矩阵恒等式
\[[dV]I^*=[I^*]dV,\]得到
\[dS=[I^*]dV+[V]dI^*.\]当 $Y$ 为固定网络参数时,
\[dI=YdV,\]因此
\[dI^*=Y^*dV^*.\]于是
\[\boxed{ dS=[I^*]dV+[V]Y^*dV^*. }\]6.2 代入电压微分
代入
\[dV=j[V]d\theta+[u]dm,\] \[dV^*=-j[V^*]d\theta+[u^*]dm,\]可得
\[\begin{aligned} dS ={}&[I^*]\left(j[V]d\theta+[u]dm\right)\\ &+[V]Y^*\left(-j[V^*]d\theta+[u^*]dm\right). \end{aligned}\]将 $d\theta$ 和 $dm$ 的系数分别收集:
\[dS=S_\theta d\theta+S_m dm.\]其中
\[S_\theta =j[I^*][V]-j[V]Y^*[V^*].\]由于两个对角矩阵可以交换次序,$[I^][V]=[V][I^]$,故
\[\boxed{ S_\theta =j[V]\left([I^*]-Y^*[V^*]\right). }\]同理,
\[\boxed{ S_m =[I^*][u]+[V]Y^*[u^*]. }\]这里
\[S_\theta=\frac{\partial S}{\partial\theta}, \qquad S_m=\frac{\partial S}{\partial m},\]都是 $n\times n$ 的复矩阵。
6.3 四个传统雅可比块
因为
\[S=P+jQ,\]所以
\[dP=\operatorname{Re}(dS), \qquad dQ=\operatorname{Im}(dS).\]因此关于完整状态 $(\theta,m)$ 的实数雅可比为
\[\boxed{ J_{PQ,\mathrm{full}} = \begin{bmatrix} \operatorname{Re}(S_\theta)&\operatorname{Re}(S_m)\\ \operatorname{Im}(S_\theta)&\operatorname{Im}(S_m) \end{bmatrix}. }\]这正是传统教材中的四个块:
\[H=\frac{\partial P}{\partial\theta}, \qquad N=\frac{\partial P}{\partial m},\] \[M=\frac{\partial Q}{\partial\theta}, \qquad L=\frac{\partial Q}{\partial m}.\]矩阵方法的关键优势是:四个块不是分别推导,而是从两个复矩阵 $S_\theta$、$S_m$ 直接取实部和虚部得到。
附注:实际潮流实现时,通常并不直接使用完整雅可比矩阵,而是删除参考母线角度列、删除 PV/参考母线电压幅值列。
7. 支路潮流的统一矩阵表达
潮流收敛后,还需要计算线路两端功率和线路损耗。
7.1 支路电压、电流与功率
支路两端电压:
\[V_f=C_fV, \qquad V_t=C_tV.\]支路两端电流:
\[I_f=Y_fV, \qquad I_t=Y_tV.\]支路两端复功率:
\[\boxed{ S_f=[V_f]I_f^*, \qquad S_t=[V_t]I_t^*. }\]按该方向约定,$S_f$、$S_t$ 都表示对应母线向支路注入的功率,因此支路复损耗为
\[\boxed{S_{\mathrm{loss}}=S_f+S_t.}\]7.2 支路功率导数
如果需要支路潮流灵敏度、最优潮流或状态估计,可以继续使用同一个乘积法则:
\[dS_f=[I_f^*]C_fdV+[V_f]Y_f^*dV^*.\]代入电压微分后:
\[\boxed{ S_{f,\theta} =j\left( [I_f^*]C_f[V]-[V_f]Y_f^*[V^*] \right), }\] \[\boxed{ S_{f,m} =[I_f^*]C_f[u]+[V_f]Y_f^*[u^*]. }\]末端公式只需将下标 $f$ 替换为 $t$。可见,节点功率和支路功率的导数具有完全相同的结构:
\[\boxed{ \text{功率微分} = \text{电压变化引起的项} + \text{电流变化引起的项}. }\]8. 软件实现的关键点
矩阵推导与现代软件工程实现天然一致。
一个最小计算内核只需完成:
V = m .* exp(j * theta)
Vnorm = V ./ abs(V)
I = Y * V
S = V .* conj(I)
# MATPOWER/RustPower 风格的完整复导数:
# dS_dVa = ∂S/∂theta, dS_dVm = ∂S/∂m
dS_dVa = 1j * diag(V) * conj(diag(I) - Y * diag(V))
dS_dVm = diag(V) * conj(Y * diag(Vnorm)) + diag(conj(I)) * diag(Vnorm)
随后根据 PV、PQ 索引抽取雅可比行列并求解线性方程。
工程实现应注意:
- $Y$、$Y_f$、$Y_t$ 和雅可比通常是稀疏矩阵;
- 不必真正生成大量稠密对角矩阵,可用稀疏对角算子或逐元素缩放;
- 相角必须使用弧度;
- 网络矩阵改变后才需要重建 $Y$,迭代中主要更新 $V,I,S,J$;
- 新实现应使用中心有限差分或实变量自动微分校验解析导数。
软件层面最重要的原则不是代码技巧,而是保持功率正号、支路方向、共轭符号和节点索引的一致性。
拓展专题1:RustPower代码示例
用 RustPower 的开源源码作为代码参照,说明“统一矩阵化交流潮流推导”如何落到一个高性能 Rust 实现中。
对于简单支路,教学推导写作
\[Y=D^T[y]D.\]RustPower 的 create_y_bus 先构造 incidence matrix,再执行稀疏乘法形成 y_bus。代码中 incidence_matrix 的每一列对应一条支路,首端填 $+1$、末端填 $-1$。
if topo.0[0] >= 0 {
incidence_matrix.push(topo.0[0] as usize, idx, Complex64::one());
}
if topo.0[1] >= 0 {
incidence_matrix.push(topo.0[1] as usize, idx, -Complex64::one());
}
let y_bus = &incidence_matrix * (diag_admit * incidence_matrix.transpose());
若有变压器或其他两端口补丁,RustPower 再向 trans_patch_matrix 中加入对应的四个端口导纳项:
trans_patch_matrix.push(from.0 as usize, from.0 as usize, p[(0, 0)]);
trans_patch_matrix.push(to.0 as usize, to.0 as usize, p[(1, 1)]);
trans_patch_matrix.push(from.0 as usize, to.0 as usize, p[(0, 1)]);
trans_patch_matrix.push(to.0 as usize, from.0 as usize, p[(1, 0)]);
这对应高阶推导中的
\[Y=C_f^TY_f+C_t^TY_t+Y_{\mathrm{sh}},\]或更一般的局部两端口矩阵装配。注意:实际工程中,变压器的变比方向、基准电压和标幺换算必须与数据接口一致;数学公式本身只给出结构。
Newton 主循环中最核心的一行是计算
\[S(V)=[V](YV)^*.\]RustPower 中的对应写法如下:
let mut mis = &v.component_mul(&(Ybus * &v).conjugate()) - Sbus;
let mut F = DVector::zeros(n_state);
assemble_f_v2(&mut F, n_bus, &mis, n_state, npq);
对于雅可比矩阵处理,教学公式为:
\[S_m=[I^*][u]+[V]Y^*[u^*],\] \[S_\theta=j[V]\left([I^*]-Y^*[V^*]\right).\]RustPower 早期实现 dSbus_dV_old 直接按 MATPOWER 记号构造稀疏对角矩阵。摘录如下:
let dS_dVm = &diagV * (Ybus * &diagVnorm).conjugate()
+ diagIbus.conjugate() * &diagVnorm;
let dS_dVa = &diagV * (diagIbus - Ybus * &diagV).conjugate()
* Complex::<f64>::i();
这两行与矩阵公式完全对应:
\[\texttt{diagV * conj(Ybus * diagVnorm)} \Longleftrightarrow [V]Y^*[u^*],\] \[\texttt{conj(diagIbus) * diagVnorm} \Longleftrightarrow [I^*][u],\] \[\texttt{diagV * conj(diagIbus - Ybus * diagV) * j} \Longleftrightarrow j[V]\left([I^*]-Y^*[V^*]\right).\]矩阵公式适合推导,但如果每次都构造对角矩阵并做稀疏矩阵乘法,会产生额外开销。RustPower 的优化实现预先进行符号分析,直接填充矩阵的数值,这里就不再展开。
最终结论是:矩阵推导不仅有助于教学上的统一理解,也更容易用于实现高性能稀疏矩阵计算和现代工程软件。
拓展专题2:把 PQ、PV 与平衡节点还原为设备控制模式
经典潮流通常把母线直接分成 PQ、PV 和平衡节点。这种分类对教学很方便,但容易产生一个误解:似乎节点类型属于网络拓扑本身。实际上,母线只承担连接作用;类型来自连接在该母线上的注入设备采用了何种控制。
1. 全空间功率平衡
设母线 $i$ 上聚合发电注入为
\[s_{g,i}=p_{g,i}+jq_{g,i},\]负荷消耗为
\[s_{\ell,i}=p_{\ell,i}+jq_{\ell,i},\]其他受控设备向交流网注入的功率为 $s_{c,i}$。按本文正号,完整节点平衡为
\[\boxed{ S_i(V)-s_{g,i}+s_{\ell,i}-s_{c,i}=0. }\]把所有母线堆叠,可写成
\[g_a(V,s_g,s_c)=S(V)-C_gs_g+S_\ell-C_cs_c=0.\]这一层对所有母线完全相同。所谓 PQ、PV 或平衡模式,只决定发电或可控注入设备还要附加哪两条控制方程。
2. 三种经典控制模式
对一个具有两个可调稳态自由度 $(p_g,q_g)$ 的聚合发电设备,可写出三种典型模式。
PQ 模式:有功和无功均给定,
\[g_P=p_g-p_g^\star=0, \qquad g_Q=q_g-q_g^\star=0.\]网络平衡决定该母线的 $\theta_i,m_i$。传统“PQ 母线”正是消去 $p_g,q_g$ 后的结果。固定功率负荷也是 PQ 注入,只是其净注入通常为负。
PV 模式:有功和电压幅值给定,
\[g_P=p_g-p_g^\star=0, \qquad g_V=m_i-m_i^\star=0.\]此时 $q_g$ 成为代数输出,用于实现电压调节;网络平衡决定 $\theta_i$ 和所需无功。传统 PV 潮流把 $m_i$ 直接从状态中消去,并在收敛后反算 $q_g$。
参考/平衡模式:电压相角和幅值给定,
\[g_\theta=\theta_i-\theta_i^\star=0, \qquad g_V=m_i-m_i^\star=0.\]此时 $p_g,q_g$ 均为代数输出,用来补偿系统总有功、无功不平衡及损耗。每个同步交流孤岛至少需要一个相角参考;但“承担相角参考”与“承担全部有功平衡”在数学上可以分开,后者可通过参与因子实现分布式平衡。
由此可以看出,三类节点的差别只是控制方程的选择:
| 模式 | 被控制量 | 自由设备输出 | 主要网络未知量 |
|---|---|---|---|
| PQ | $p_g,q_g$ | 无 | $\theta_i,m_i$ |
| PV | $p_g,m_i$ | $q_g$ | $\theta_i$ |
| 参考 | $\theta_i,m_i$ | $p_g,q_g$ | 无 |
在全空间形式中,节点 $P,Q$ 平衡始终保留;在经典约化形式中,则通过消元得到“PV 只保留 $P$ 方程、PQ 保留 $P,Q$ 方程”的结构。
3. 多台设备接在同一母线时
母线类型应由控制资源的聚合关系决定,而不是简单地给每台机组单独标注 PV。若同一母线有多台发电机,母线电压只有一个,因此只能有一个独立的本地电压控制方程。各机组无功可以通过参与因子、无功下垂或优先级分配,例如
\[q_{g,k}=q_{g,k}^{0}+\alpha_k q_{\mathrm{reg}}, \qquad \sum_k\alpha_k=1.\]某台机组达到限值后,应把它从自由调节集合移出,再由剩余机组重新归一化参与因子。若所有调节资源均达到限值,母线整体才失去电压控制并转入 PQ 限值模式。
远方电压控制也应按设备端口建模。若发电机位于母线 $i$,却调节母线 $r$ 的电压,则控制方程是
\[m_r-m_r^\star=0,\]而不是把母线 $r$ 粗略标记为普通 PV。此时控制设备所在节点、被控节点和无功平衡节点可能并不相同。
4. PV 到 PQ 的无功限值转换
PV 模式隐含假设
\[q_g^{\min}<q_g<q_g^{\max}.\]若求得
\[q_g>q_g^{\max},\]说明维持 $m_i=m_i^\star$ 所需的无功超过设备能力。正确处理不是简单裁剪输出后继续保留电压方程,而是用限值方程替换电压控制:
\[q_g-q_g^{\max}=0,\]同时释放
\[m_i-m_i^\star=0.\]设备于是从 PV 模式进入 $PQ^{\max}$ 模式,$m_i$ 重新成为未知量。下限同理:
\[q_g-q_g^{\min}=0.\]在经典约化雅可比中,这一转换表现为增加一个电压幅值未知量和一个无功平衡方程;在全空间形式中,变量维度可以保持不变,只需替换一条控制残差。
5. 状态机表达
可将单个电压调节设备的离散模式定义为
\[\sigma_g\in\{PV,PQ^{\max},PQ^{\min}\}.\]典型转移为:
| 当前模式 | 触发条件 | 新模式 | 被移除方程 | 被加入方程 |
|---|---|---|---|---|
| $PV$ | $q_g>q_g^{\max}+\varepsilon_{on}$ | $PQ^{\max}$ | $m_i-m_i^\star=0$ | $q_g-q_g^{\max}=0$ |
| $PV$ | $q_g<q_g^{\min}-\varepsilon_{on}$ | $PQ^{\min}$ | $m_i-m_i^\star=0$ | $q_g-q_g^{\min}=0$ |
| $PQ^{\max}$ | 恢复电压控制所需 $q_g<q_g^{\max}-\varepsilon_{off}$ | $PV$ | $q_g-q_g^{\max}=0$ | $m_i-m_i^\star=0$ |
| $PQ^{\min}$ | 恢复电压控制所需 $q_g>q_g^{\min}+\varepsilon_{off}$ | $PV$ | $q_g-q_g^{\min}=0$ | $m_i-m_i^\star=0$ |
其中通常取
\[\varepsilon_{off}>\varepsilon_{on}>0\]形成滞回,避免模式在边界附近反复振荡。所谓“恢复电压控制所需的无功”不能只读取当前限值模式中的 $q_g$,因为它已被固定在边界;应通过试探 PV 求解、灵敏度预测或外层活动集检验判断。
工程中可采用的默认策略之一,是在一次基态求解内只允许 $PV\rightarrow PQ$ 单向转换;若确需恢复,则在内层牛顿收敛后进行试探,而不是在每个未收敛迭代点来回切换。
6. 平衡节点达到限值时
参考角必须保留,但承担系统有功或无功平衡的设备可以更换。若参考母线发电机达到 $P$ 或 $Q$ 限值,不应删除孤岛的相角参考;应把“角度参考”与“功率平衡”拆开:
- 保持某个母线的 $\theta=\theta^\star$;
- 将有功失配变量分配给其他具有调节能力的机组;
- 无功限值按 PV/PQ 逻辑处理;
- 若所有平衡资源均耗尽,应报告该控制模式不可行,而不是继续输出一个越限解。
这一处理对含 HVDC 和储能的系统尤其重要,因为直流换流器也可能参与交流区域的有功平衡,但它并不能提供两个异步交流岛之间的相角参考。
附录:常用矩阵微分与复向量微分公式
本附录只列一般数学规则,不再把速查表写成电力系统专用公式。正文中的节点功率、支路功率和雅可比公式,本质上都是把这些通用规则代入具体对象后得到的结果。
1. 基本符号与微分约定
设
\[a,b,z\in\mathbb C^n,\qquad x\in\mathbb R^p,\qquad A\in\mathbb C^{m\times n}.\]本文使用
\[[a]=\operatorname{diag}(a)\]表示以向量 $a$ 为对角线的对角矩阵。逐元素乘法与对角矩阵乘法等价:
\[a\odot b=[a]b=[b]a.\]逐元素除法记为
\[a\oslash b.\]复共轭、转置和共轭转置必须严格区分:
\[A^*=\text{逐元素复共轭},\qquad A^T=\text{普通转置},\qquad A^H=(A^*)^T=\text{共轭转置}.\]若未明确写出转置符号,星号 $*$ 只表示逐元素复共轭。本文把复向量函数视为实状态变量的复值函数,即
\[f:\mathbb R^p\rightarrow\mathbb C^q, \qquad df=J_f\,dx, \qquad J_f=\frac{\partial f}{\partial x}\in\mathbb C^{q\times p}.\]这类微分不要求 $f$ 对复变量解析;凡含有共轭、模值、实部、虚部或相角的函数,都应按实变量微分理解。
2. 公式速查表
以下公式是本文反复使用的通用规则(包括状态估计也用到,本文未详细展开)。表中 $A,B$ 可为矩阵,$a,b,z$ 可为向量,$x$ 为实状态变量;维度均假定相容。
| 类型 | 一般公式 | 说明 |
|---|---|---|
| 线性映射 | 若 $y=Ax$ 且 $A$ 固定,则 $dy=A\,dx$ | 最基本的链式构件 |
| 仿射映射 | 若 $y=Ax+b$ 且 $A,b$ 固定,则 $dy=A\,dx$ | 常数项微分为零 |
| 矩阵乘积 | $d(AB)=(dA)B+A(dB)$ | 普通乘积法则 |
| 矩阵—向量乘积 | 若 $y=A(x)b(x)$,则 $dy=(dA)b+A\,db$ | 当矩阵依赖状态时必须保留 $(dA)b$ |
| 逆矩阵 | $d(A^{-1})=-A^{-1}(dA)A^{-1}$ | 要求 $A$ 非奇异 |
| 转置与共轭 | $d(A^T)=(dA)^T$,$d(A^)=(dA)^$,$d(A^H)=(dA)^H$ | 对实状态微分成立 |
| 对角算子 | $d([a])=[da]$ | 对角矩阵逐元素微分 |
| 对角乘向量 | $d([a]b)=[b]da+[a]db$ | 常用来处理逐元素乘积 |
| Hadamard 乘积 | $d(a\odot b)=[b]da+[a]db$ | 与上一行等价 |
| 实部投影 | 若 $dz=J_zdx$,则 $d\operatorname{Re}(z)=\operatorname{Re}(J_z)dx$ | 把复导数转为实雅可比 |
| 虚部投影 | 若 $dz=J_zdx$,则 $d\operatorname{Im}(z)=\operatorname{Im}(J_z)dx$ | 同上 |
| 复模平方 | $d(z^\odot z)=2\operatorname{Re}([z^]dz)$ | 对零点仍可微 |
| 复模 | $d\lvert z\rvert=\operatorname{Re}!\left(\left[z^*\oslash \lvert z\rvert\right]dz\right)$ | 要求每个 $z_i\ne0$ |
| 复相角 | $d\arg z=\operatorname{Im}([z]^{-1}dz)$ | 要求每个 $z_i\ne0$,且局部不跨越分支切口 |
| 链式法则 | 若 $y=f(x)$、$h=g(y)$,则 $dh=J_gJ_fdx$ | 所有组合映射的核心 |
| 行选择 | 若 $y_{\mathcal I}=R_{\mathcal I}y$,则 $dy_{\mathcal I}=R_{\mathcal I}dy$ | 用于抽取子向量或残差行 |
| 列嵌入 | 若 $dx=C_{\mathcal I}d\xi$,则 $df=J_fC_{\mathcal I}d\xi$ | 用于删除参考变量或固定变量 |
| 块堆叠 | 若 $h=\operatorname{col}(f,g)$,则 $dh=\operatorname{col}(df,dg)$ | 用于装配残差向量 |
| 二次型 | 若 $\phi=x^TAx$ 且 $A$ 固定,则 $d\phi=x^T(A+A^T)dx$ | 若 $A=A^T$,则 $d\phi=2x^TA dx$ |
| 加权平方残差 | 设$W=W^T$,若 $\phi=\frac12 r^TWr$ 且 $W$ 固定,则 $d\phi=r^TWdr$ | 最小二乘正规方程的来源 |
| 隐函数 | 若 $g(x,y)=0$ 且 $g_y$ 非奇异,则 $dy/dx=-g_y^{-1}g_x$ | 用于消去内部变量 |
3. 如何从通用公式回到正文推导
正文中的核心步骤只是上述规则的代入。例如,对任意复向量 $a,b$,由
\[d([a]b)=[b]da+[a]db\]可知,若某个量写成“一个向量的对角矩阵乘另一个向量”,则其微分一定由两项组成:第一项来自前一个向量的变化,第二项来自后一个向量的变化。本文中节点功率和支路功率的微分,正是这一规则在具体电力系统变量上的应用。
再例如,若一个复值母函数满足
\[da=A_xdx,\]那么它的实部、虚部、模值、模平方和相角量测,都可分别由上表的实部投影、虚部投影、复模、复模平方和复相角公式得到。这样做的意义是:不需要为每一种量测重新手算元素级三角函数导数,只需先得到复母函数的一阶微分,再按量测类型投影即可。