企业官网建设流程全解析

前一段时间接了个分布式光伏接入规划的活儿，甲方给了几个备选接入点，让我评估"哪里接入最划算"。这类问题在配电网里太常见了，而最常规的做法就是跑一遍灵敏度分析，算算"网损对节点注入功率的敏感程度"，然后按灵敏度大小排序选位置。最开始我也是这么干的，用的就是IEEE33节点配电网经典算例，Matlab里把潮流一跑、雅可比一求、灵敏度一算，交差了。结果甲方拿着实际仿真数据过来问：为什么你推荐的节点，实测降损效果排不到前面？

后来我重新把整个方法捋了一遍，发现问题的根源不在潮流计算，而在"灵敏度分析"本身。传统灵敏度基于基准运行点做一次线性展开，本质上是个静态指标，但分布式电源接入后，节点电压和网络损耗都会因注入功率发生非线性改变，一次灵敏度排序在高渗透率场景下很容易失真。这篇文章就围绕【灵敏度分析】33节点配电网改进灵敏度分析（Matlab代码实现）展开，把我的完整实现思路、核心代码、对比结果和踩坑记录写出来，给同样在做配电网DG选址、网损优化的朋友做个参考。

1. IEEE33节点系统为什么成了配电网研究的"标准试验田"

做配电网规划、运行优化、DG选址的研究人员，基本绕不开IEEE33节点系统。它最早来自美国PG&E的一个实际配电网，经过标准化处理后成为公开算例，电压等级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原始数据中根节点是0，Matlab的数组下标从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) = 0，x是节点电压幅值和相角向量，P_loss = Σ P_i。因为网损是通过潮流方程间接依赖注入功率的，所以要解一个线性方程组：

(∂F/∂x) * (∂x/∂P_i) = -∂F/∂P_i

也就是先算雅可比矩阵，再回代求解。得出所有节点的灵敏度值之后，按从大到小排序，排在最前面的节点就是"理论上接入DG降网损最明显的位置"。

2.2 传统灵敏度在DG选址中的常规用法

实际工程项目中，传统灵敏度的用法非常简单粗暴：

1. 跑一次基准潮流；
2. 利用潮流雅可比矩阵计算每个负荷节点的网损灵敏度；
3. 排序，取前几名作为DG待选接入点；
4. 如果不放心，再对候选点分别做一次含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(n+1: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(n+1:2*n); % 取对应有功注入的部分 end

另一种是数值差分法，简单粗暴但对验证非常有用。程序里加个测试开关，用差分结果和解析结果对比，如果两者差得很远，基本可以断定要么雅可比矩阵构造错了，要么网损偏导公式抄错。这个技巧我在平时调试程序时反复用，强烈推荐你也保留下来。

4.4 改进灵敏度的主循环

把前面三个改进整合到主程序里，逻辑是：

1. 输入候选节点列表和DG容量；
2. 对每个候选点，用步进法完成迭代灵敏度修正；
3. 在每次步进里同时累计二阶修正项；
4. 计算联合指标时，先算出网损灵敏度和电压灵敏度两个向量，再做加权组合；
5. 输出每个候选点的综合排序。

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 两种灵敏度排序结果

实验得到的排序结果如下表：

上表数据是我从多次仿真结果里取的代表性数值，实际工况不同会有浮动，但趋势很稳定：传统灵敏度把节点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节点，就要考虑加速了。我的做法有两个：

一是把候选节点的迭代循环写成parfor，Matlab并行计算工具箱直接支持，改动成本很低。二是步长自适应，早期DG容量小的时候可以大步长，后期灵敏度变化快了用小步长，能在精度几乎不变的情况下减少三分之一次数。这个方法对大规模网络特别实用。

最后还有一个小技巧：改进灵敏度排序不要只取第一名，一般取前3到5个节点，分别做一次完整潮流仿真验证，再结合网架结构、线路容量、保护配置等因素综合拍板。灵敏度分析是筛选工具，不是最终决策工具，这一点在实际工程里非常重要。我在处理甲方项目时，都会在报告里加一页"候选节点多维度对比表"，把灵敏度排序、实测网损、电压指标、线路负载率全部列出来，这样既展示了改进方法的优势，也避免把算法结果当成唯一标准，项目过审顺利很多。