解析计算节点电压灵敏度：雅可比矩阵方法及MATLAB实现

简介：面向电力系统分析与电气工程领域学生、研究者的MATLAB计算工具包，聚焦节点电压灵敏度系数的解析计算，可服务于课程设计、期末大作业及毕业设计等场景。压缩包共11个文件，包含5个.m主程序与函数脚本、3个.mat案例数据集、2个.png结果示意图及1个.xlsx数据表格，整体仅117KB，轻量易部署。代码支持MATLAB 2014/2019a/2024a，采用参数化编程，参数可灵活修改，注释详细，思路清晰；附赠案例数据可直接运行，免去额外预处理。通过该工具，使用者能快速掌握灵敏度系数的计算流程，并应用于电力系统稳定性评估、故障分析与优化调整。目前已有56人学习下载，适合需要结合理论解析与工程实现的高校学生及初级研究人员。

1. 为什么说解析计算是节点电压灵敏度的理想打开方式

节点电压灵敏度系数的解析计算最近在配电网和输电网分析里又翻红，原因在于它比传统扰动法快一个数量级，且不会因步长选择而引入截断误差。简单说，灵敏度系数回答的是“当某个节点的注入功率变化 1 个标幺值时，其他节点电压幅值会变多少”这一类问题，对电压稳定评估、无功补偿配置和分布式电源选址都有直接参考价值。这个 MATLAB 资源正是围绕这一需求设计的：它自带 IEEE34 三相算例，提供 main.m 、SC_Voltage.m 和 LoadFlow 模块，用户在 MATLAB2014/2019a/2024a 任一版本上都能直接跑。适合电力系统课程设计，也适合刚接触灵敏度分析的工程师快速对比“解析法”与“数值扰动法”的差异。

2. 灵敏度系数的解析推导：从雅可比矩阵到电压对注入功率的偏导

2.1 潮流方程与雅可比矩阵的关系

节点电压灵敏度并不是一个凭空定义的概念，它本质上是潮流方程在某运行点的一阶 Taylor 展开系数。先看最常用的极坐标形式潮流方程：

[ P_i = V_i \sum_{j \in i} V_j (G_{ij}\cos\theta_{ij} + B_{ij}\sin\theta_{ij}) ]

[ Q_i = V_i \sum_{j \in i} V_j (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij}) ]

将全部节点的有功、无功功率不平衡量写成向量，与节点电压相角和幅值偏差之间就形成线性化关系：

[ \begin{bmatrix} \Delta P \ \Delta Q \end{bmatrix} = J \begin{bmatrix} \Delta \theta \ \Delta V \end{bmatrix} ]

这里 (J) 就是牛顿-拉夫逊法潮流计算中反复构造的雅可比矩阵。习惯上把子矩阵记为：

[ J = \begin{bmatrix} H & N \ M & L \end{bmatrix} ]

其中 (H) 对应 (\partial P / \partial \theta)，(N) 对应 (\partial P / \partial V)，(M) 对应 (\partial Q / \partial \theta)，(L) 对应 (\partial Q / \partial V)。

我们要求的电压灵敏度系数，即 (\partial V / \partial P) 和 (\partial V / \partial Q)，正好藏在 (J^{-1}) 里。因为对两侧同时求逆可以得到：

[ \begin{bmatrix} \Delta \theta \ \Delta V \end{bmatrix} = J^{-1} \begin{bmatrix} \Delta P \ \Delta Q \end{bmatrix} ]

把 (J^{-1}) 写成对应维度的四块：

[ J^{-1} = \begin{bmatrix} S_{\theta P} & S_{\theta Q} \ S_{V P} & S_{V Q} \end{bmatrix} ]

这里的 (S_{V P}) 就是电压幅值对节点有功注入的灵敏度矩阵，(S_{V Q}) 是电压幅值对无功注入的灵敏度矩阵。也就是说，解析计算的核心就是“先形成潮流雅可比矩阵，再对其求逆并提取子块”，完全不需要对每个节点做扰动重算潮流。

2.2 节点类型如何影响灵敏度矩阵的维度

实际电网里节点不全是 PQ 节点，还有平衡节点和 PV 节点。PV 节点电压幅值给定，所以其对电压幅值的灵敏度没有物理意义；平衡节点相角给定，其相角灵敏度也不参与讨论。因此，标准做法是先将雅可比矩阵写成全维“未处理”形式，然后把 PV 节点对应的 (Q) 平衡方程和电压幅值列去掉，把平衡节点对应的 (P,\theta) 行和列也去掉，得到一个降维可逆矩阵。

这个资源包里的 SC_Voltage.m 处理的正是这个过程。常见做法是维护一个节点类型向量type，其中 1 表示 PQ，2 表示 PV，3 表示平衡节点。降维索引可以这样构造：

% 节点类型：1=PQ 2=PV 3=平衡 pq = find(type == 1); pv = find(type == 2); slack = find(type == 3); % 保留的功率平衡方程索引：PQ全部，PV只保留有功方程 p_idx = [pq; pv]; q_idx = pq; % 状态变量索引：相角仅保留PQ和PV，电压幅值仅保留PQ theta_idx = [pq; pv]; v_idx = pq; % 从全维雅可比J_full中抽取降维矩阵 J11 = J_full(theta_idx, theta_idx); J12 = J_full(theta_idx, v_idx); J21 = J_full(q_idx, theta_idx); J22 = J_full(q_idx, v_idx); J_reduced = [J11 J12; J21 J22];

这段代码的逻辑是：平衡节点不参与方程，PV 节点的无功方程被删掉，其电压幅值也不作为状态量。这样 (J_{reduced}) 才是方阵且通常可逆。代码里如果用inv(J_reduced)会直接得到灵敏度，但实际工程中更推荐用J_reduced \ eye(n)或S = inv(J_reduced)后的子块提取。

2.3 从逆矩阵提取电压灵敏度子块

一旦求出 (J_{reduced}^{-1})，按之前的索引映射回去即可得到原始节点的电压灵敏度。由于我们最终关心每个 PQ 节点的电压幅值对所有 PQ 节点有功、无功注入的偏导，直接取：

S_inv = inv(J_reduced); % S_inv 按行是[theta_idx中的节点; v_idx中的节点] % 按列是[p_idx注入; q_idx注入] n_pq = length(pq); S_VP = S_inv(n_pq+1 : n_pq+length(v_idx), 1 : length(p_idx)); S_VQ = S_inv(n_pq+1 : n_pq+length(v_idx), length(p_idx)+1 : end);

这里S_VP的第 (i) 行第 (j) 列表示 PQ 节点 (i) 的电压幅值对 PQ/PV 节点 (j) 的有功注入的灵敏度。S_VQ同理，但注入节点只限于 PQ 节点，因为 PV 节点无功是待定状态，不能作为独立扰动源。

值得注意的是，解析法求得的灵敏度矩阵是稠密的。如果你用的是 IEEE34 这类多节点系统，直接inv(J_reduced)在内存上没压力，但到几千节点规模时，最好用稀疏 LU 分解并求解一组单位向量，否则内存墙会先拦住你。后面第 5 章会专门讲这个优化。

2.4 SC_Voltage.m 核心代码段注释

我按项目里SC_Voltage.m的功能，把最关键的求逆和子块提取浓缩成下面这段带注释代码：

function [S_VP, S_VQ] = SC_Voltage(Jac, idx) % Jac : 牛顿潮流最后一次迭代的雅可比矩阵 % idx : 结构体，包含pq, pv, slack索引 p = [idx.pq; idx.pv]; q = idx.pq; t = [idx.pq; idx.pv]; v = idx.pq; J11 = Jac(t, t); J12 = Jac(t, v); J21 = Jac(q, t); J22 = Jac(q, v); Jred = [J11 J12; J21 J22]; % 用左除避免显式求逆，数值更稳定 S = Jred \ eye(size(Jred)); npq = length(idx.pq); npv = length(idx.pv); ncol_p = length(p); % 有功注入数量 ncol_q = length(q); % 无功注入数量 % 电压幅值灵敏度位于S的后v_idx块 SV = S(npq+npv+1 : npq+npv+length(v), :); S_VP = SV(:, 1:ncol_p); S_VQ = SV(:, ncol_p+1 : ncol_p+ncol_q); end

这里用Jred \ eye(size(Jred))等效于求逆，但底层走高斯消元，比inv(Jred)更稳定，尤其在矩阵接近奇异时能少一点数值误差。S_VP的行号顺序与idx.pq一致，列号顺序与[idx.pq; idx.pv]一致。我在实际项目中会在函数入口打印这几个维度的提示，避免后续和别人写的数据格式混淆。

3. 程序模块拆解：LoadFlow 如何支撑 SC_Voltage.m

3.1 项目文件结构与执行流程

下载解压后，核心文件包括main.m、SC_Voltage.m、LoadFlow（可能是子目录或脚本）、input examples下的 IEEE34 数据文件，以及两张参考输出图1.png、2.png。整体流程是：

这个流程把一个复杂的工程问题拆成了三个可独立替换的单元。想换电网模型，只改input examples；想换潮流算法，只要保证输出接口还是雅可比矩阵；想换灵敏度求解方式，只动SC_Voltage.m。这就是所谓“参数化编程”的典型思路。

3.2 LoadFlow 的接口约定

从代码衔接看，LoadFlow至少需要返回三个东西：节点电压向量V、节点功率注入向量S、以及牛顿法最后一步雅可比矩阵Jac。在 MATLAB 中，常见的调用语句是：

[V, S, Jac, type, base] = LoadFlow('ieee34', '3ph');

这里'ieee34'指定算例名称，'3ph'表示三相潮流模式。IEEE34 是一个三相不平衡配电网，节点编号达到 30 个以上，线路参数包含不对称线路电抗。LoadFlow 内部会先把三相线路转成节点导纳矩阵，再按三相分量分别建立功率方程。

雅可比矩阵的输出顺序必须和节点编号一致。我见过不少读者自己写灵敏度程序时踩坑：潮流算完但不输出雅可比，最后用有限差分去近似灵敏度，绕了一大圈。实际上牛顿法的最后一次雅可比迭代已经收敛，直接交给SC_Voltage.m就能用，没必要重算。

3.3 参数化编程：哪些参数可以改

代码注释里强调“参数可方便更改”，这一点对自主扩展很有价值。最常见的可调参数有：

% main.m 顶部参数区 baseMVA = 1; % 基准容量，单位 MVA tol = 1e-8; % 潮流收敛误差 maxIter = 30; % 最大迭代次数 pvNodes = [12, 18]; % 指定哪些节点作为 PV 节点（如有分布式电源） loadScale = 0.85; % 负荷整体缩放系数

这些参数改变后，潮流解和灵敏度矩阵都会随之变化。特别是loadScale，它能在不修改原始 IEEE34 数据的前提下，模拟系统重负荷工况。我一般会把它从 0.5 扫到 1.5，观察同一节点的S_VQ变化趋势，判断是否接近电压崩溃。

3.4 如何快速替换成自己的电网数据

如果你想算自己的网络，需要把节点、支路、变压器数据整理成与input examples相同的格式。IEEE34 示例里每一行代表一个支路段，带有始端节点、终端节点、电阻、电抗、电容等。对应到 MATLAB 代码中，LoadFlow内部通过稀疏矩阵组装导纳阵。替换数据时，注意编号必须从 1 开始连续，否则稀疏组装会丢节点。

下面是典型的支路数据格式（示例，非完整）：

% bus_i, bus_j, R_ohm_km, X_ohm_km, C_nF_km, 长度km branch_data = [ 1, 2, 0.1053, 0.1238, 5.3, 0.85; 2, 3, 0.1053, 0.1238, 5.3, 1.20; 3, 4, 0.0752, 0.0961, 4.1, 0.90; ];

替换后只需确保type向量正确设置每个节点类型，SC_Voltage.m不关心具体物理参数，只看雅可比矩阵和索引。这也是把灵敏度计算从潮流中解耦出来的价值所在。

4. 在 IEEE34 三相节点算例上的复现与结果验证

4.1 运行 main.m 的准备工作

解压后，先在 MATLAB 中将当前路径设为项目根目录。在命令窗口执行：

run('main.m')

如果 MATLAB 版本是 2014 或更早，请先确认代码里没有使用arguments块、string新语法等 2016b 之后才有的特性。项目标注支持 2014/2019a/2024a，说明作者已尽量用兼容写法，但为保险起见，运行前可以用checkcode('main.m')检查潜在语法警告。

main.m 内部一般会调用LoadFlow获得基准潮流，然后调用SC_Voltage.m。运行结束后，工作区会生成S_VP和S_VQ两个矩阵。同时会输出两张图：第一张是各节点电压幅值随某个注入节点有功变化的折线图，第二张是灵敏度矩阵的热力图。

4.2 结果矩阵的直观解读

以 IEEE34 的 PQ 节点为例，假设系统中有 26 个 PQ 节点、4 个 PV 节点。运行后S_VQ是 26×26 矩阵，第 (i) 行第 (j) 列表示节点 (j) 注入 1 MVar 无功时，节点 (i) 电压幅值的变化（标幺值）。正常情况下，对角元是正值，且明显大于非对角元，说明本地无功注入对本节点电压支撑最直接。

我运行后得到的前 5 个对角元大致如下（数值为标幺值，示意）：

这说明离电源端越远的节点，其电压对本地无功注入越敏感。注意，这里的值受基准容量影响，如果你把baseMVA从 1 改成 10，对角线会缩小约 10 倍，所以比较不同运行点时必须统一基准。

4.3 用灵敏度焦耳热力图排序

SC_Voltage.m返回的S_VQ矩阵可以直接用imagesc可视化：

figure; imagesc(S_VQ); colorbar; xlabel('注入节点编号'); ylabel('电压观察节点编号'); title('S_VQ，各节点无功注入对电压的灵敏度');

从热力图能看到明显的对角线亮带，如果某些非对角线元素也特别亮，说明该节点对远处节点的无功注入也敏感，这就是薄弱节点。我通常会把每列最大值所在的位置找出来，标记为“最影响该节点的注入源”，再结合loadScale参数做重复试验，基本能锁定系统中最需要无功补偿的 2 到 3 个节点。

4.4 与数值扰动法对比验证

解析法算出的灵敏度是否可信，最直接的验证是拿扰动法对比。做法是：对第 (j) 个 PQ 节点注入无功增加一个小量 (\Delta Q_j)（比如 0.01 MVar），重新运行一次潮流，记录所有 PQ 节点电压变化 (\Delta V_i)，得到：

dQ = 0.01; S_VQ_numeric = zeros(n_pq, n_pq); for j = 1:n_pq Q_new = Q0; Q_new(j) = Q0(j) + dQ; V_new = LoadFlow_with_injection(Q_new); S_VQ_numeric(:, j) = (abs(V_new(pq)) - abs(V0(pq))) / dQ; end

然后计算误差：

err = max(max(abs(S_VQ - S_VQ_numeric))); fprintf('最大误差：%.4e\n', err);

这一步很关键。如果采用的是收敛后的雅可比矩阵，那么解析法和扰动法的误差通常在 (10^{-6}) 量级，偏差完全来自潮流重新求解时的迭代精度。如果发现误差大到 (10^{-3})，最可能的原因是SC_Voltage.m里用了第一次迭代的雅可比矩阵，或者降维索引对不上。

5. 提升解析灵敏度实用性的三个排错技巧

5.1 用稀疏 LU 分解代替显式求逆

当节点数超过 500 时，inv(Jred)会变得又慢又占内存。我的习惯是改用lu分解后回代：

[L, U] = lu(Jred); S = U \ (L \ eye(size(Jred)));

更省内存的做法是不直接构造整个S，而是对每个需要输出的行 (e_i) 求解一次：

S_i = e_i' / Jred; % 等价于 (Jred' \ e_i)' ?

但 MATLAB 里更推荐S_i = e_i' / Jred;，它求解的是行向量的线性方程。计算全部节点时可以循环，或者直接一次性S = Jred \ speye(n)，后者在稀疏模式下会采用稀疏 LU，比inv快不少。

5.2 雅可比矩阵奇异时的处理

接近电压崩溃点时，(J_{reduced}) 的行列式趋近于零，求逆会出现巨大数值。此时矩阵的条件数condest(Jred)可能超过 (10^{12})。我一般先检查条件数，如果过大，就不做全矩阵求逆，而是改用 Tikhonov 正则化：

lambda = 1e-6 * max(abs(Jred(:))); Jreg = Jred + lambda * eye(size(Jred)); S_reg = Jreg \ eye(size(Jreg));

这种做法会带来一定误差，但在崩溃点附近，解析值本身已经失去意义，正则化后的结果反而能提示“最危险的灵敏度方向”。在工程报告中，我会把这种情况标注为“接近电压稳定极限，灵敏度值仅供参考”。

5.3 三相不平衡系统的灵敏度合并

IEEE34 是三相模型，LoadFlow输出的电压和雅可比矩阵实际上按三相分别列出。SC_Voltage.m 默认可能只处理正序分量，但用户如果想得到“三相综合灵敏度”，可以用电压幅值的两点间偏差比例合成。常见做法是先提取三相电压幅值 (V_a, V_b, V_c)，用平均电压 (V_{avg}) 计算灵敏度：

V_avg = (abs(V_a) + abs(V_b) + abs(V_c)) / 3; S_VQ_avg = S_VQ_a + S_VQ_b + S_VQ_c; % 按注入三相功率等量扰动

此时要确保潮流中注入功率也按三相分别赋值，否则平均灵敏度与单相灵敏度的量纲会不一致。我在处理有单相光伏接入的案例时，会进一步拆出每个节点的相别，否则容易把同一点的不同相灵敏度混成一个错误值。这个坑在配电网三相分析里几乎必踩，建模时把相别信息一并传入SC_Voltage.m就能规避。

本文还有配套的精品资源，点击获取