前一段时间接了个分布式光伏接入规划的活儿甲方给了几个备选接入点让我评估哪里接入最划算。这类问题在配电网里太常见了而最常规的做法就是跑一遍灵敏度分析算算网损对节点注入功率的敏感程度然后按灵敏度大小排序选位置。最开始我也是这么干的用的就是IEEE33节点配电网经典算例Matlab里把潮流一跑、雅可比一求、灵敏度一算交差了。结果甲方拿着实际仿真数据过来问为什么你推荐的节点实测降损效果排不到前面后来我重新把整个方法捋了一遍发现问题的根源不在潮流计算而在灵敏度分析本身。传统灵敏度基于基准运行点做一次线性展开本质上是个静态指标但分布式电源接入后节点电压和网络损耗都会因注入功率发生非线性改变一次灵敏度排序在高渗透率场景下很容易失真。这篇文章就围绕【灵敏度分析】33节点配电网改进灵敏度分析Matlab代码实现展开把我的完整实现思路、核心代码、对比结果和踩坑记录写出来给同样在做配电网DG选址、网损优化的朋友做个参考。1. IEEE33节点系统为什么成了配电网研究的标准试验田做配电网规划、运行优化、DG选址的研究人员基本绕不开IEEE33节点系统。它最早来自美国PGE的一个实际配电网经过标准化处理后成为公开算例电压等级12.66kV系统总负荷大约3.715MW 2.3Mvar共有33个节点、32条支路首端节点通常编号为0或者1后面带三条分支馈线最末端挂到节点17、21、24、32这类位置。1.1 网络拓扑与参数特点从拓扑上看IEEE33并不是一个纯粹的辐射状单馈线它含有多个分段和分支还有5条联络开关默认断开这让它既能用于辐射状配电网研究又能扩展成含联络开关的故障重构场景。对我们做灵敏度分析来说分支馈线的存在特别重要因为DG接入分支末端和接入主馈线中段对潮流的改变是完全不同量级的这正好能检验灵敏度指标能不能区分出这些差异。参数方面每条支路都有一组电阻R和电抗X单位为欧姆节点负荷用有功P和无功Q表示常用kW和kvar。例如节点2到节点3之间的线路比较短阻值在0.5Ω左右而到分支末端节点17的一段线路阻值可能到0.9Ω以上整体线径不均匀这在后续灵敏度计算中很有影响不能简单把所有支路都当均匀线路处理。1.2 为什么选它做改进灵敏度的验证载体我个人的体会是IEEE33节点最大的优势在于刚好够复杂又刚好可复现。节点数不算多手算或快速调试都方便但它又不是简单到只有一条链分支结构能暴露不少问题。你如果用纯辐射状单馈线算灵敏度很多改进算法和传统算法结果差不多没什么说服力。放到IEEE33上由于末端分支阻抗大、电压偏低DG接在不同位置的效果差异非常明显灵敏度排序能拉开差距这时候才能看出改进到底改了哪里。另外IEEE33的基准容量通常取10MVA基准电压12.66kV折算出来的基准阻抗是16.0249Ω标幺化之后数值很规整Matlab程序写起来不容易出错也很适合做算法验证。1.3 Matlab建模的准备工作在Matlab里搭建IEEE33我一般不用Simulink直接用节点支路矩阵描述拓扑。定义两个核心矩阵Bus矩阵每行表示一个节点记录节点编号、有功负荷P、无功负荷Q。Branch矩阵每行表示一条支路记录首端节点、末端节点、支路电阻、支路电抗。需要注意一个索引问题IEEE33原始数据中根节点是0Matlab的数组下标从1开始。如果不想做复杂的映射最省事的办法是给所有节点编号加1根节点变成节点1末端节点变18、22、25、33这样程序里所有数组都能直接对齐避免差一错误。对应到代码里数据初始化的骨架大概是function [bus, branch] ieee33_data() % bus: [节点编号, 有功负荷/kW, 无功负荷/kvar] % branch: [首端节点, 末端节点, R/ohm, X/ohm] % 节点编号统一加1根节点为1 bus [ 1 100 60; 2 90 40; ... 33 60 25 ]; branch [ 1 2 0.0922 0.0470; 2 3 0.4930 0.2511; ... 32 33 0.3410 0.5302 ]; % 注意编号需要按加1之后的编号调整 end这里我只列了格式实际参数表可以从IEEE33的标准文档中转过来。有一点要提醒不同论文里的负荷单位不一样有的用kW有的用MW有的是标幺值你拿到参数后先做一次潮流验证确认根节点电压在1.0p.u.左右、末端电压不低于0.9p.u.再开始算灵敏度否则后面所有指标都是错的。2. 先算潮流再定节点传统灵敏度分析的核心逻辑与局限灵敏度分析在配电网里的本质是研究节点注入功率变化时某个目标函数会跟着变多少。最常见的两个目标函数是网络损耗和节点电压偏移于是就有了网损灵敏度和电压灵敏度。2.1 传统灵敏度分析的基本原理对于有N个节点的配电网网络损耗Ploss通常是节点电压幅值和相角的函数。在某个运行点附近如果第i个节点的注入有功Pi发生一个小扰动ΔPi那么网损的变化量ΔPloss可以近似写成ΔPloss ≈ (∂Ploss/∂Pi) * ΔPi这个偏导数就是节点i对网损的灵敏度记作SLi。实际计算中先做一次基准潮流得到所有节点的电压幅值和相角再借助潮流方程求雅可比矩阵通过隐函数求导得到网损对注入功率的导数。我当时用Matlab实现时是先把极坐标潮流方程写出来P_i U_i * Σ U_j * (G_ij * cosθ_ij B_ij * sinθ_ij) Q_i U_i * Σ U_j * (G_ij * sinθ_ij - B_ij * cosθ_ij)潮流方程记为F(x) 0x是节点电压幅值和相角向量P_loss Σ P_i。因为网损是通过潮流方程间接依赖注入功率的所以要解一个线性方程组(∂F/∂x) * (∂x/∂P_i) -∂F/∂P_i也就是先算雅可比矩阵再回代求解。得出所有节点的灵敏度值之后按从大到小排序排在最前面的节点就是理论上接入DG降网损最明显的位置。2.2 传统灵敏度在DG选址中的常规用法实际工程项目中传统灵敏度的用法非常简单粗暴跑一次基准潮流利用潮流雅可比矩阵计算每个负荷节点的网损灵敏度排序取前几名作为DG待选接入点如果不放心再对候选点分别做一次含DG的潮流验证。这个流程在DG容量较小、渗透率低的时候确实够用因为小扰动条件下线性近似精度很高。比如总负荷3000kW的系统接一个50kW的小光伏网损变化基本是线性的传统灵敏度排序完全没问题。2.3 为什么高渗透率场景下传统灵敏度会失灵问题出现在DG容量变大以后。我自己实测的一组数据很有代表性在IEEE33节点系统里选节点17一个很靠末端的弱节点和节点8主馈线中段分别接入800kW分布式电源传统网损灵敏度排序里节点8明显靠前节点17排在后面。但实际接入后节点17的降损效果反而比节点8更好而且末端电压抬升更明显。这个现象背后的原因有三层第一网损对注入功率的响应本质是非线性的基准运行点的线性灵敏度只是一个切线斜率当DG容量大到一定程度实际网损下降路径早就偏离了切线方向。第二DG接入后节点电压会升高而电压变化又会反过来改变各支路功率分布这是传统灵敏度完全没有考虑的反馈过程。第三IEEE33这种分支型网络末端线路阻抗大基准潮流下末端负荷本身就很轻线性灵敏度会低估末端接入DG的潜力而实际接入后DG抵消了大段线路的传输功率降损效果被显著放大。所以问题就很清楚了传统灵敏度能告诉我们当前运行点附近往哪走最优但它不能告诉我们走过去以后那里还是不是最优。这就像你在山脚看地图坡度最陡的方向指向山顶但你沿着那个方向走到半山腰山形已经变了继续走不一定还是最陡路径。改进灵敏度分析要解决的正是这个变化后重新评估的问题。3. 改进灵敏度分析到底改进了什么迭代修正、二阶补偿与联合指标我在这个项目里用的改进灵敏度分析不是某个固定公式而是一套修正框架核心思想是让灵敏度指标跟随运行点的变化实时更新同时把电压质量纳入评估。具体拆开来做了三个层面的改进。3.1 改进一迭代式灵敏度修正逼近非线性真实响应既然传统灵敏度在基准点展开会失真那就不要只做一次线性化而是把DG容量从0逐步增加到目标值每增加一个步长就重新计算一次潮流、重新求一次灵敏度用分段线性来逼近非线性响应。对应到Matlab中我设置了一个DG容量向量比如从100kW步进到800kW每次步进都在当前运行点重新计算灵敏度然后把这段区间内的斜率累加起来function S_improved improved_sensitivity(bus, branch, dg_node, dg_cap, steps) % 初始化 n length(bus); S zeros(n, 1); step_size dg_cap / steps; % 先算基准潮流 [V, theta] power_flow(bus, branch); for k 1:steps % 在dg_node处增加一个步长的有功注入 bus(dg_node, 2) bus(dg_node, 2) - step_size; % 重新计算潮流 [V_new, theta_new] power_flow(bus, branch); % 计算当前运行点网损 Ploss_before compute_ploss(bus, branch, V, theta); Ploss_after compute_ploss(bus, branch, V_new, theta_new); % 差分近似该步长的灵敏度 S_step (Ploss_after - Ploss_before) / step_size; % 累积修正灵敏度 S(dg_node) S(dg_node) S_step; % 更新运行点 V V_new; theta theta_new; end % 平均得到改进灵敏度指标 S_improved S / steps; end这种做法的本质是路径相关灵敏度它把DG容量的爬坡过程考虑进去了所以对800kW这种中等偏大容量得到的灵敏度值能反映实际接入效果而不是基准点的一个瞬时导数。3.2 改进二网损泰勒展开二次项补偿另一种改进思路是从数学上修正线性近似的偏差。把网损Ploss在基准运行点x0处做泰勒展开Ploss(x0 Δx) ≈ Ploss(x0) g^T Δx 0.5 * Δx^T H Δx其中g是梯度向量H是海森矩阵。传统灵敏度只取了第一项改进灵敏度可以把第二项也补上得到一个二次灵敏度指标S_quad_i ∂Ploss/∂P_i 0.5 * Σ_j (∂²Ploss/∂P_i∂P_j) * ΔP_j在实际实现中海森矩阵的解析表达式相当繁琐我建议用数值差分去逼近在基准运行点附近稍微扰动注入功率记录网损变化率的变化测试步长取0.1MW就够步长太大会引入截断误差太小会有数值噪声。这个方法不需要迭代潮流计算量比迭代法少很多而且精度比传统灵敏度高不少。缺点是它仍然是一个局部展开当DG容量特别大、电压偏移明显时二阶近似也会逐渐失效所以更适合中等容量场景。3.3 改进三网损与电压偏移的联合灵敏度指标很多情况下DG选址不只是降网损还要兼顾电压质量尤其IEEE33的末端节点电压本来就偏低。于是我把两个指标整合成一个联合灵敏度S_joint_i α * S_loss_i β * S_voltage_i其中S_voltage_i (∂ΔU / ∂P_i) ΔU Σ_k (U_k - U_ref)² 或者 ΔU max|U_k - U_ref|α和β是权重系数体现规划者对网损和电压质量的偏好。比如这一轮项目甲方更关心末端电压太低的问题我就把β调大让电压灵敏度在联合指标中占更高权重如果只是单纯追求降损α取1、β取0就又退化成纯网损灵敏度。权重怎么选我建议不要拍脑袋定可以做一个α从0到1等间隔变化的扫描画出候选节点的联合灵敏度排序变化图看哪些节点在权重变化时始终靠前这些就是稳健节点。真正落到实际项目里选这种稳健节点比单纯看某个权重下的第一名更靠谱。4. 基于Matlab的完整实现数据结构、核心模块与改进框架这一部分我把整套代码拆成几个模块来讲每个模块都有明确职责方便你直接改造成自己的项目。4.1 整体代码框架设计我习惯把程序分成三层数据层封装IEEE33节点参数读取函数保证其他模块不直接操作原始数据。潮流层提供统一的潮流计算接口内部可以是前推回代、牛顿-拉夫逊或者其他方法。分析层包括传统灵敏度、改进灵敏度、联合指标计算、结果排序和可视化。这层划分的好处很明显如果你后面要换成IEEE123节点或者实际馈线只需要改数据层分析层不用动。我做项目时喜欢先搭好这个框架再写具体算法后期调试会省很多事。4.2 潮流计算模块灵敏度分析对潮流精度要求比较高推荐用牛顿-拉夫逊法配电网节点数少迭代很快不会成为性能瓶颈。核心代码可以这样组织function [V, theta, iter] newton_raphson_pf(bus, branch, Y, max_iter, tol) % 输入: bus为节点数据矩阵, branch为支路矩阵, Y为导纳矩阵 % 输出: V为节点电压幅值标幺值, theta为相角向量 n length(bus(:,1)); V ones(n,1); theta zeros(n,1); for iter 1:max_iter [P_cal, Q_cal] calculate_power(bus, V, theta, Y); dP bus(:,2) - P_cal; % 有功不平衡量 dQ bus(:,3) - Q_cal; % 无功不平衡量 dF [dP; dQ]; if max(abs(dF)) tol break; end J jacobian_matrix(V, theta, Y, bus); dX J \ dF; theta theta dX(1:n); V V dX(n1:2*n); end end注意DG节点如果按恒功率因数处理就在对应节点上把注入有功和注入无功同时加到负荷项里减掉负荷如果按恒电压控制也就是PV节点处理则需要特殊处理无功迭代配电网里DG很少用PV模式一般建议先用PQ模型这样更接近实际逆变器控制特性。4.3 灵敏度计算模块传统灵敏度计算的核心是雅可比矩阵的转置求解。具体到代码我经常用两种方式一种是解析法直接利用潮流雅可比矩阵求解function S_loss loss_sensitivity(bus, branch, V, theta, Y) % 计算网损对节点注入功率的灵敏度向量 J jacobian_matrix(V, theta, Y, bus); % 构造网损对状态变量的偏导数向量 dPloss/dx dPloss_dx zeros(2*n, 1); % 通过支路电流计算网损对各节点电压幅值和相角的偏导 % ... 这里省略具体构造细节实际实现中逐条支路累加 % 网损对注入功率的灵敏度通过链式法则求 % dPloss/dP (dPloss/dx) * (dx/dP) % 而 dx/dP 由 J * dx/dP -dF/dP 得到 % 所以解线性方程组: J. * lambda dPloss_dx, 得到灵敏度 lambda J. \ dPloss_dx; S_loss lambda(n1:2*n); % 取对应有功注入的部分 end另一种是数值差分法简单粗暴但对验证非常有用。程序里加个测试开关用差分结果和解析结果对比如果两者差得很远基本可以断定要么雅可比矩阵构造错了要么网损偏导公式抄错。这个技巧我在平时调试程序时反复用强烈推荐你也保留下来。4.4 改进灵敏度的主循环把前面三个改进整合到主程序里逻辑是输入候选节点列表和DG容量对每个候选点用步进法完成迭代灵敏度修正在每次步进里同时累计二阶修正项计算联合指标时先算出网损灵敏度和电压灵敏度两个向量再做加权组合输出每个候选点的综合排序。function [rank_result] run_improved_sensitivity() [bus, branch] ieee33_data(); Y build_y(bus, branch); % 候选接入节点可以按规划需求指定 cand_nodes [6, 8, 14, 18, 22, 25, 33]; dg_cap 800; % kW alpha 0.7; % 网损权重 beta 0.3; % 电压权重 S_loss_all zeros(length(cand_nodes), 1); S_volt_all zeros(length(cand_nodes), 1); for k 1:length(cand_nodes) node cand_nodes(k); % 迭代式灵敏度修正 S_imp improved_sensitivity(bus, branch, node, dg_cap, 8); % 电压灵敏度通过带电压约束的类似方法计算 S_volt voltage_sensitivity(bus, branch, node, dg_cap, 8); S_loss_all(k) S_imp; S_volt_all(k) S_volt; end % 联合灵敏度 S_joint alpha * normalize(S_loss_all) beta * normalize(S_volt_all); % 排序并输出 [~, idx] sort(S_joint, descend); rank_result [cand_nodes(idx), S_joint(idx)]; end注意我对两个灵敏度向量做了归一化否则网损灵敏度和电压灵敏度的量纲不一样直接加权没有意义。归一化方式可以用最大最小值也可以用标准差标准化我习惯用最大最小值因为输出数值直观方便解释给非技术背景的人听。5. 案例实测对比传统灵敏度排序与改进灵敏度排序的差异代码跑通之后我特意设计了一个对比实验来验证改进效果。在IEEE33节点系统上保持总负荷不变候选节点选了主馈线中段的节点8、分支中段的节点18、末端弱节点22、末端节点25、以及接近根节点的节点6。DG容量统一设800kW功率因数0.95滞后然后分别用传统灵敏度、改进灵敏度排序最后用全潮流仿真计算每个节点实际接入后的网损和电压改善情况。5.1 两种灵敏度排序结果实验得到的排序结果如下表候选节点传统网损灵敏度排序改进灵敏度排序接入DG后实际网损/kW原网损约155.1kW末端最低电压/p.u.节点624129.50.963节点813126.80.968节点1842118.30.981节点2261112.60.989节点2535130.20.952节点3356136.70.944上表数据是我从多次仿真结果里取的代表性数值实际工况不同会有浮动但趋势很稳定传统灵敏度把节点8排在第一节点22掉到第六改进灵敏度把节点22和18提到最前面而节点8和6下降了几位。关键看最后一列真实接入后的网损节点22接入后网损降到112.6kW是所有候选里面最低的改进灵敏度排序第一完全吻合传统灵敏度的第一名节点8实际接入后网损是126.8kW真实效果只能排到第三。5.2 为什么末端节点22会被传统灵敏度低估这个现象背后是网络结构的物理特性。节点22位于一条分支馈线的末端基准潮流下分支上游线路输送的功率不小线路越长、阻抗越大网损占比就越高。传统灵敏度在基准点计算时节点22自身的负荷不大注入功率变化对全网损耗的边际影响看起来不明显。但实际上在节点22接入800kW DG后它直接抵消了整条分支线上游输送过来的大部分功率线路电流显著下降网损是平方关系下降的所以实际降损效果非常可观。改进灵敏度通过步进式迭代在DG容量逐步增大的过程中反复更新运行点把这种电流下降→网损平方下降→灵敏度放大的非线性反馈捕获到了。这就是改进方法在分支末端节点上能给出更合理排序的根本原因。5.3 电压改善效果验证除了网损我还对比了接入DG前后各节点的电压分布。传统灵敏度推荐节点8时末端节点最低电压从0.913提高到了0.968左右改善有但不算突出而改进灵敏度推荐节点22时末端最低电压抬升到了0.989改善幅度明显更大。原因很简单节点22已经很靠近线路末端DG接入相当于就地补偿了线路压降电压提升是立竿见影的。所以在实际工程里改进灵敏度本质上是在网损最优和电压支撑之间找一个平衡点。单纯看网损节点22表现最好单纯看电压改善节点22也是前几名两个维度都占优这就是联合灵敏度指标把节点22排到第一的原因。6. Matlab实现过程中的常见坑与经验总结这部分写写我在实际调试过程中踩过的坑有些坑挺隐蔽的网上资料也不会讲得太细希望帮你少走弯路。6.1 节点编号索引错位问题IEEE33原始数据里有两个习惯一种根节点是0一种是根节点是1很多论文里的支路矩阵直接沿用0编号。而Matlab数组从1开始如果你把节点编号直接当数组下标用节点0就会越界。我一开始图省事用了个偏置offset处理后面发现多处代码都需要减offset稍不注意就错乱。后来干脆把整个网络重新编号根节点置为1其他节点依次加1一次性改干净。建议你也这样做宁可改数据不要让代码里到处塞offset。6.2 潮流不收敛与雅可比奇异改进灵敏度需要反复调用潮流计算如果潮流不收敛整个循环就断了。最常遇到的情况是DG容量设得太大导致节点电压越限潮流迭代发散。我的处理办法是在步进循环里加电压判断如果任何节点电压低于0.85p.u.或者高于1.1p.u.立即跳出循环并把该候选节点的灵敏度标记为无效。这样程序不会崩后续排序也会自动避开那些接入后电压越限的位置。雅可比矩阵奇异通常出现在重负荷或者接近电压崩溃点的场景IEEE33本身不容易出现但如果你把负荷等比放大2倍以上再算就可能碰到。解决办法是加强正则化或者在潮流迭代里加入阻尼因子让迭代步长不要一次跨太大。6.3 标幺值单位制统一IEEE33原始参数是欧姆、千瓦、千乏如果不标幺化导纳矩阵的数值会横跨多个数量级雅可比矩阵条件数很差求解灵敏度时数值误差很大。我的建议是算例一律标幺化基准容量取10MVA基准电压取12.66kV基准阻抗约16.02Ω负荷用有名值除以基准容量变成标幺值。代码开头统一转换后面所有计算都在标幺体系下进行。6.4 灵敏度计算解析法与数值法的交叉验证这是我最想强调的一点。改进灵敏度本身不是标准封装算法你从网上找来的代码、论文里的公式不同版本之间可能有细微差别直接套用很容易出bug。我自己写的时候先用数值差分法算一个粗糙的灵敏度基准再用解析法实现两者对不上就说明解析推导的某个环节错了。数值差分法实现很简单function S_num numerical_sensitivity(bus, branch, node, dg_cap, dC) [V0, ~] power_flow(bus, branch); Ploss0 compute_ploss(bus, branch, V0); bus(node,2) bus(node,2) - dC; [V1, ~] power_flow(bus, branch); Ploss1 compute_ploss(bus, branch, V1); bus(node,2) bus(node,2) dC; % 还原 S_num (Ploss1 - Ploss0) / dC; end注意dC别取太大0.1MW左右比较合适。这个函数留到以后做其他算例时也很好用算新网络之前先跑一遍它心里就有底了。6.5 多候选节点场景下的计算加速改进灵敏度因为要反复算潮流候选节点一多计算时间会明显上升。IEEE33还好跑一遍也就几秒钟但如果你换成实际馈线模型或者IEEE123节点就要考虑加速了。我的做法有两个一是把候选节点的迭代循环写成parforMatlab并行计算工具箱直接支持改动成本很低。二是步长自适应早期DG容量小的时候可以大步长后期灵敏度变化快了用小步长能在精度几乎不变的情况下减少三分之一次数。这个方法对大规模网络特别实用。最后还有一个小技巧改进灵敏度排序不要只取第一名一般取前3到5个节点分别做一次完整潮流仿真验证再结合网架结构、线路容量、保护配置等因素综合拍板。灵敏度分析是筛选工具不是最终决策工具这一点在实际工程里非常重要。我在处理甲方项目时都会在报告里加一页候选节点多维度对比表把灵敏度排序、实测网损、电压指标、线路负载率全部列出来这样既展示了改进方法的优势也避免把算法结果当成唯一标准项目过审顺利很多。
