简介面向电力系统、电子信息和计算机等专业的学生及研究人员此MATLAB代码包聚焦节点电压灵敏度系数的解析计算可服务于课程设计、期末大作业与毕业设计中的系统稳定性评估和参数影响分析。包内共11个文件包含5个.m脚本、3个.mat数据文件、2个.png结果示意图及1个.xlsx表格覆盖主程序、潮流计算和IEEE34节点算例等模块整体仅117KB。代码采用参数化编程运行环境兼容MATLAB2014、2019a、2024a并附赠可直接运行的案例数据便于快速换参复现。代码思路清晰、注释详细特别适合初学者理解灵敏度计算的完整逻辑。目前已有56人学习使用。通过这套代码使用者能掌握节点电压对系统参数变化的敏感程度计算为电力系统优化设计、故障分析及可靠性评估提供直接可用的工具与参考。1. 节点电压灵敏度系数的解析计算从潮流雅可比矩阵中直接取数做无功电压优化时我经常要回答这样一个问题在某条母线补一组电容器周围电压到底能抬起来多少。最直接的办法是把负荷功率加一点再重新跑一次潮流但那样要反复调用潮流收敛计算量线性增长更稳的办法是在牛顿-拉夫逊潮流已经收敛的工况点上把雅可比矩阵求逆直接读出电压随注入功率变化的系数。这种做法就叫节点电压灵敏度系数的解析计算。它和数值扰动法在数学上是同一阶近似但只需要一次矩阵分解就能得到全部节点的电压-有功、电压-无功对应关系。对新稳态分析、薄弱点筛选、优化模型里的灵敏度约束都很适合。下面把这套计算从公式推到 MATLAB 代码再给出验证参数写代码的人能直接照着搭。2. 节点电压灵敏度系数的解析计算数学基础V-Q 与 V-P 灵敏度从哪来2.1 灵敏度不是一张能沿用很久的固定表格有人会把电压灵敏度理解成一张“节点阻抗表”离电源越远灵敏度越高。这个印象在轻载时成立但运行点一变灵敏度矩阵变化很大。重载、无功不足、电压偏低时雅可比矩阵接近奇异个别节点的灵敏度系数会突然放大好几倍。换句话说它是潮流解在当前工况下的一阶偏导不在通用表格里。解析计算存在的意义就是每次拿到新的潮流断面马上在同一个稳态点上算出一组与工况一致的系数供调度或优化使用。2.2 从潮流方程推到 ΔV K⁻¹ΔQ极坐标潮流方程中对第 i 条母线的注入功率可以写成P_i V_i Σ_j V_j (G_ij cosθ_ij B_ij sinθ_ij)Q_i V_i Σ_j V_j (G_ij sinθ_ij - B_ij cosθ_ij)在潮流收敛点附近做一阶线性化得到ΔP J_Pθ Δθ J_PV ΔVΔQ J_Qθ Δθ J_QV ΔV这里 ΔP、ΔQ 是节点注入功率的扰动量Δθ 是相角摄动ΔV 是电压幅值摄动。J 的下标表示对哪个变量求偏导它在内容上就是牛顿-拉夫逊潮流最后一次迭代使用的雅可比矩阵。如果只关心无功注入变化对电压的影响最自然的约束是认为有功注入在本次扰动中不变即 ΔP0。由第一个方程解出 ΔθΔθ -J_Pθ⁻¹ J_PV ΔV代入第二式ΔQ ( J_QV - J_Qθ J_Pθ⁻¹ J_PV ) ΔV定义 K J_QV - J_Qθ J_Pθ⁻¹ J_PV这就是通常说的 V-Q 灵敏度矩阵。于是ΔV K⁻¹ ΔQ节点电压灵敏度系数就是 K⁻¹ 中的元素。K⁻¹ 的第 i 行第 j 列物理含义是第 j 条 PQ 母线的无功注入变化一个单位时第 i 条母线电压幅值的变化量单位是 p.u./p.u.。日常说“母线自电压灵敏度”指的是对角元“互电压灵敏度”指的是非对角元。2.3 解析法和数值扰动法同一阶近似两种路径实际编程时多数 MATLAB 代码并不先把 J_Pθ⁻¹ 显式解出来再构造 K。更常见的做法是把完整潮流雅可比矩阵求逆取右下角与电压幅值对应的分块。这两种方法在数学上完全相同但完整求逆在 MATPOWER 里更容易写因为四个子块的索引可以直接从母线类型里得到。下面这张表列出了两种路径的区别对比项解析法矩阵逆/分解数值扰动法差分计算量一次矩阵分解多个右端项复用每个扰动点都要重算一次潮流精度由潮流收敛精度决定受差分步长和迭代精度双重影响非线性误差只有一阶项天然线性化步长大会混入二阶项适用面批量灵敏度、优化模型求导单点校验、小规模交叉验证理解这一点之后代码实现的关键就不再是“要不要求逆”而是怎样把雅可比矩阵的行列顺序和母线类型对齐。顺序错一位取出来的块就不知道是在描述哪两个节点之间的关系了。3. MATLAB 中节点电压灵敏度系数的解析计算核心程序实现3.1 从潮流结果拿到雅可比矩阵的三个准备步骤在 MATPOWER 体系里节点数据用 bus 矩阵表达第 2 列是母线类型1 表示 PQ 节点2 表示 PV 节点3 表示平衡节点。计算灵敏度的第一步是跑一次潮流第二步是根据潮流结果重建导纳矩阵第三步是用 MATPOWER 提供的接口计算复功率对相角和电压幅值的偏导。下面代码用 IEEE 9 节点算例说明流程。% 节点电压灵敏度系数解析计算主入口 mpc loadcase(case9); % 读入 IEEE 9 节点算例 res runpf(mpc); % 先求潮流收敛解 baseMVA mpc.baseMVA; bus res.bus; n size(bus, 1); % 按母线类型分类pv 是电压控制节点pq 是负荷节点 pv find(bus(:, 2) 2); pq find(bus(:, 2) 1); pvpq [pv; pq]; % 相角变量除平衡节点外全都要这一段的关键是把 pv 和 pq 按固定顺序拼好。因为最后雅可比矩阵的行列顺序完全由这里决定后面取灵敏度子块时也用同一组索引顺序不一致就会张冠李戴。3.2 雅可比矩阵组装与右下角灵敏度子块有了分类索引再从潮流结果重建电压相量用 MATPOWER 的dSbus_dV接口得到四个子块最后按牛顿法的经典排法拼成完整雅可比矩阵。% 从潮流结果重建导纳矩阵和电压相量 [Ybus, ~, ~] makeYbus(mpc.bus, mpc.branch); Vm bus(:, 8); Va bus(:, 9) * pi / 180; V Vm .* exp(1j * Va); % dSbus_dV 返回对相角和电压幅度的复功率偏导 [dSbus_dVa, dSbus_dVm] dSbus_dV(Ybus, V); % 组装牛顿法雅可比矩阵 J11 real(dSbus_dVa(pvpq, pvpq)); % dP/dVa J12 real(dSbus_dVm(pvpq, pq)); % dP/dVm J21 imag(dSbus_dVa(pq, pvpq)); % dQ/dVa J22 imag(dSbus_dVm(pq, pq)); % dQ/dVm J [J11 J12; J21 J22]; % 完整求逆小算例可以直接做 Ainv inv(J); npq length(pq); row_base length(pvpq); Vrows row_base (1:npq); % 左下角开始的 PQ 电压行 % 右下角子块即 dVm/dQ dVm_dQ Ainv(Vrows, Vrows);这段代码的逻辑顺序很重要J 的左半部分是相角列右半部分是 PQ 电压列上半行是 PVPQ 节点的有功方程下半行是 PQ 节点的无功方程。因此逆矩阵的右下角块正好对应 PQ 节点电压对 PQ 节点无功注入的灵敏度。Vrows的起点是length(pvpq)也就是把相角变量全部数完之后才轮到电压变量。inv(J)只适合几百节点以内的小算例。真实电网规模超过几千节点时完整求逆既慢又浪费内存应该改成稀疏 LU 分解后逐列回代。% 大型电网改用它对每个右端项只做两次三角回代 [Lm, Um, pm, qm] lu(J, vector); % 以第 i 个 PQ 母线无功注入为扰动源 i 1; b zeros(size(J, 1), 1); b(Vrows(i)) 1; % 只在第 i 列加单位注入 y Lm \ b(pm); % 先解 L x zeros(size(J, 1), 1); x(qm) Um \ y; % 再解 U dVdq_i x(Vrows); % 第 i 列解析灵敏度b(pm)先按行置换取索引求解后再把解放回qm对应的位置。对多个 PQ 母线同时仿真时只需循环不同的b向量分解只做一次这套写法在节点规模大时优势非常明显。3.3 参数说明和常见顺序错误这里的dSbus_dV是 MATPOWER 的内部函数调用它时电压相量必须是复数形式不能只传幅值。.mat文件、Excel 导入的母线编号也建议先统一转换成 Matpower 的 bus 格式再走同样流程。常见顺序错误是pv 和 pq 用find取出来后没有拼接直接按各自顺序参与索引。这样 J11 等子块的行列坐标就不是同一个坐标空间最后dVm_dQ的每一行无法对应到具体母线编号。写代码时记住一点所有子块共用的行、列索引都必须来自同一份pvpq和pq向量不要新造变量。4. 节点电压灵敏度系数的应用场景与 MATLAB 参数设置4.1 用 V-Q 灵敏度定位电压薄弱母线拿到了 dVm_dQ 矩阵第一件能做的事是筛选电压薄弱节点。做法是统计每条 PQ 母线受其他 PQ 节点无功扰动时可能产生的最大电压变化。变化越大的母线通常在系统中离无功电源越远或处在重载通道末端优先把它列为无功补偿候选点。% 对每行取绝对值最大值作为母线电压强度指标 row_max max(abs(dVm_dQ), [], 2); [~, order] sort(row_max, descend); for k 1:min(5, length(order)) j order(k); fprintf(PQ母线 %d: 最大灵敏度 %.5f p.u./p.u.\n, ... pq(j), row_max(j)); end取“行最大值”而不只取对角元是因为某条母线电压可能更多受邻近其他节点的无功影响仅看自灵敏度会低估问题。输出结果里灵敏度单位是 p.u./p.u.意思是给节点注入 1 p.u. 无功电压变化多少 p.u.。这里 1 p.u. 在 100 MVA 系统里是 100 Mvar所以 0.02 p.u. 的灵敏度已经算比较高工程上要认真对待。4.2 由目标电压反推最小无功补偿容量灵敏度系数最直接的用途是估算补偿容量。假设第 k 条母线当前电压比目标值低 ΔV_target想靠本地无功补偿把电压抬起来所需的注入无功增量可以近似为ΔQ ΔV_target / dVm_dQ(k, k)这里的 dVm_dQ(k,k) 是本节点的自灵敏度。若目标抬升 0.02 p.u.自灵敏度是 0.1 p.u./p.u.则 ΔQ 0.02 / 0.1 0.2 p.u.乘上 baseMVA 后得到有名值容量。dV_target 0.02; % 目标电压抬升量p.u. idx order(1); % 选灵敏度最高的母线 dq_pu dV_target / dVm_dQ(idx, idx); fprintf(建议无功注入 %.2f Mvar\n, dq_pu * baseMVA);这个公式是单节点近似没有考虑其他节点无功同时变化时的相互作用。实际多节点补偿时应该把灵敏度子矩阵拿出来做最小二乘解比如给三个候选母线同时分配无功dq_vector pinv(dVm_dQ(candidate, candidate)) * dV_target * ones(length(candidate), 1);pinv给出的解满足最小二乘意义下最平稳的无功分配适合作为优化模型初值。4.3 影响计算结果的四个参数解析计算的准确性不完全来自矩阵求逆代码更多取决于输入参数的设置。下面是四个最常影响结果的地方参数建议取值影响baseMVA与算例内所有功率数据一致灵敏度换算成 Mvar/MW 时必须乘 baseMVA取错则容量计算差一个数量级母线类型严格区分 1PQ、2PV、3平衡类型分类错雅可比矩阵维度和子块位置全乱潮流收敛精度Matpower 默认 1e-8收敛门槛太松灵敏度会带着潮流求解误差节点编号不要求连续但顺序必须与 bus 矩阵一致编号重复或乱序会破坏局部索引定位调试时建议先打印size(J)和length(pvpq)length(pq)两者必须相等。如果调试对象是老版本 MATPOWER 生成的 .mat 文件还要留意 bus 矩阵里的单位是不是已经在标幺体系下避免把有名值混进灵敏度计算。5. 用数值摄动法验证解析电压灵敏度系数的 3 个检查点5.1 检查点一矩阵维数与变量排序自检验证从维度开始。解析灵敏度矩阵 dVm_dQ 的行数必须等于 PQ 节点总数列数也等于 PQ 节点总数。如果行数列数对不上先回去看 pv 和 pq 是否分类完整。assert(isequal(size(dVm_dQ), [npq npq]), 灵敏度矩阵维度错误);这一条能挡住绝大多数索引错位问题。维度正确后再看对角线自灵敏度应当普遍大于互灵敏度的平均量级若发现某列全为零多半是 J 子块拼错或母线类型识别有误。5.2 检查点二与潮流重解的差商做交叉验证数值摄动法是最直观的验证手段。给某一个 PQ 母线的无功负荷减小一点相当于注入无功增加重新跑潮流比较电压差与解析估计。i 1; dq 1e-4; % 扰动步长p.u. mpc2 mpc; mpc2.bus(pq(i), 5) mpc2.bus(pq(i), 5) - dq; % Qd 减小代表注入增加 res2 runpf(mpc2); numeric res2.bus(pq, 8) - res.bus(pq, 8); analytic dVm_dQ(:, i) * dq; max_dev max(abs(numeric - analytic)); fprintf(最大偏差 %.3e\n, max_dev);因为解析法是一阶近似扰动步长不能取得太大。1e-4 p.u. 是兼顾差分噪声和线性化误差的常用值。偏差量级如果小于 1e-6说明解析矩阵和潮流实现使用相同的雅可比模型如果偏差在 1e-4 左右要检查是不是扰动步长踩到了非线性区或者潮流重新收敛后运行点偏移太远。5.3 检查点三用多步长扫描排除非线性干扰单个步长验证通过还不够。把 dq 分别取 1e-3、1e-4、1e-5 各跑一遍观察偏差是否随步长缩小。如果步长减小时偏差方向相反或呈跳跃变化说明存在单位换算错误如果步长很小时偏差不再下降说明潮流本身的收敛误差成为主导这是正常现象。多步长扫描还可以帮你确认矩阵分解没有发生在错误的潮流断面上。运行状态变化后只需要重新调用runpf和makeYbus数值验证代码可以原样复用。把这套校验写进回归测试用例以后每次改母线参数后都能自动盯住灵敏度矩阵的变化。本文还有配套的精品资源点击获取