企业官网建设流程全解析

看到“EI复现”这四个字，我就知道点进来的大概率是准备写电网方向论文的研究生，或者刚接触配电网重构的工程师。这个题目本身是近几年电力系统领域的典型组合：高比例清洁能源接入改变了传统配电网的单一功率流向，需求响应又给负荷侧加了灵活性，两者与配电网重构叠加，本质上是在回答一个问题——当发、用两侧都变得不确定时，网络开关拓扑该怎样调整才能既安全又经济。这篇内容把我从读论文、推导公式到在MATLAB里用YALMIP跑通算例的完整过程整理出来，可参考的不仅是代码本身，更是每一步选择背后的逻辑和踩过的坑。

为什么我觉得这类论文值得复现？因为分布式光伏和风电大量接入后，传统辐射状配电网的潮流方向不再是简单的“变电站单向往下送”，而是出现局部反送、电压越限等问题；同时需求响应的加入让负荷不再是一个固定数字，而是一组可以平移、削减的区间。把这三件事放在一个优化模型里，重构就不再是单纯的开关组合优化，而是一个多时段、多目标、带非线性潮流约束的混合整数规划问题。

1. 复现前，先把论文的“骨架”拆清楚

1.1 三个关键词在模型中分别承担什么角色

我在复现之前习惯先问一个问题：题目里的每个关键词，到底对应优化模型里的哪一块？

“高比例清洁能源接入”对应的是分布式电源模型。光伏、风电的出力有随机性，论文里通常简化为预测值加一定置信区间，或者直接给典型场景。优化模型里要做的是把分布式电源的出力当作可控变量参与功率平衡，同时为了防止算法“硬塞”分布式电源导致电压越限，目标函数里往往带一个弃电惩罚项。

“计及需求响应”对应的是负荷侧模型。需求响应分两种，一种叫价格型DR，用户根据分时电价调整用电量；一种叫激励型DR，用户直接响应调度指令削减或转移负荷。论文里多数用的是价格型需求响应，具体会落到一个弹性系数矩阵上。

“配电网重构”对应的是拓扑决策变量。重构的本质是改变联络开关和分段开关的开合状态，从而改变潮流分布、降低网损、消除过载或电压越限。这是一个离散决策问题，和潮流方程、需求响应模型咬合在一起，最终形成MISOCP（混合整数二阶锥规划）或者MILP框架。

1.2 整个模型的“输入-优化-输出”闭环

读完摘要和引言，我会先画出模型构成。一个典型的复现闭环长这样：输入是配电网原始拓扑、线路阻抗、各时段负荷曲线、分布式电源预测出力、分时电价；决策变量是各时段开关状态、各节点分布式电源实际出力、需求响应后的负荷曲线；约束条件是DistFlow潮流方程、径向拓扑约束、电压上下限、线路容量、DR调节范围；输出是最优开关组合、网损、电压分布、DR前后的负荷曲线对比。

我建议你也先画这么一张图，不用多精美，但一定要能回答一个问题：哪些量是已知的，哪些量是模型算出来的。很多论文复现失败，根源不是代码写错，而是边界没画清楚——把该固定的参数弄成了变量，或者反过来。

1.3 复现这篇论文要做到什么程度才算完成

先说结论：复现不等于把论文里的公式抄一遍，更不等于照扒代码。我的标准是三条：第一，能用自己的代码在标准算例（比如IEEE 33节点系统）上跑出和论文量级相近的结果；第二，能解释清每个约束和每个参数为什么这么设定；第三，能修改关键参数（分布式电源渗透率、弹性系数）后，现象变化的方向和论文一致。

举个例子，如果论文说“随着清洁能源渗透率提高，最优开关组合会向分布式电源集中区域靠拢”，那你的复现结果至少要能复现这个趋势。如果趋势反了，那大概率不是参数问题，而是模型结构理解错了。这个标准会指导你后面每一步调试，而不是只盯着一个网损数值看。

2. 配电网重构的数学核心：DistFlow约束与径向拓扑

2.1 DistFlow潮流方程与二阶锥松弛

配电网潮流计算几乎不会用常规牛拉法里的极坐标方程，因为配电网是辐射状、R/X比较大，牛顿法容易不收敛，而且优化模型需要潮流约束的梯度性质好看。主流做法是用DistFlow支路潮流方程，以支路有功、无功和节点电压幅值平方作为变量，推导过程在很多文献里都有，我不重复推了，直接写最终形式。

对一条支路(i,j)，假设功率从i流向j，那么有如下方程：

[ \begin{cases} P_{ij} - r_{ij}l_{ij} = \sum_{(j,k)\in \mathcal{E}} P_{jk} + P_j^{\mathrm{load}} - P_j^{\mathrm{DG}} \ Q_{ij} - x_{ij}l_{ij} = \sum_{(j,k)\in \mathcal{E}} Q_{jk} + Q_j^{\mathrm{load}} - Q_j^{\mathrm{DG}} \ V_j^2 = V_i^2 - 2(r_{ij}P_{ij} + x_{ij}Q_{ij}) + (r_{ij}^2+x_{ij}^2)l_{ij} \end{cases} ]

其中 (l_{ij} = (P_{ij}^2+Q_{ij}^2)/V_i^2)，本质就是支路电流的平方。这个方程是精确的，但非凸，因为 (l_{ij}) 的定义式是一个二次等式约束。

标准做法是把它松弛成二阶锥约束：

[ \left\lVert \begin{matrix} 2P_{ij} \ 2Q_{ij} \ l_{ij} - V_i^2 \end{matrix} \right\rVert_2 \leq l_{ij} + V_i^2 ]

为什么能松弛？从物理意义看，这样约束下 (l_{ij}) 只被要求“不小于”电流平方，而目标函数里网损是 (r_{ij}l_{ij})，为了压低网损，求解器会自动把 (l_{ij}) 往小了压。只要目标函数里网损权重为正，松弛在最优解处就是紧的。这意味着我们能安心用一个凸问题去近似原非凸问题。

2.2 径向拓扑约束：虚拟潮流法

配电网重构最大的难点不是潮流方程，而是“重构结果必须是辐射状”这个拓扑约束。辐射状的意思是：整个网络连通，且正好有 (N-1) 条闭合支路，不能成环。

直接写“不能成环”在数学上很麻烦，论文里常用一个技巧：虚拟潮流法。想象每个负荷节点都向根节点（变电站母线）发出单位虚拟流量，那么：

• 对于每个非根节点，流出该节点的虚拟流量总和等于流入减1；
• 对于根节点，总流出等于 (N-1)；
• 只有被选中的支路（开关闭合）才允许传输虚拟流量。

用变量 (f_{ij}) 表示支路虚拟流量，(s_{ij}) 表示开关状态，那么约束写出来是：

[ \sum_{(i,j)\in\mathcal{E}} f_{ij} - \sum_{(j,k)\in\mathcal{E}} f_{jk} = 1 \quad \forall j \neq \text{根节点} ]

同时辅助约束 ( -M s_{ij} \leq f_{ij} \leq M s_{ij})，确保断开支路不传导虚拟流。加上 (\sum s_{ij} = N-1)，就能把辐射状结构锁死。我最初在YALMIP里为了省事只加了“闭合支路总数等于N-1”，结果跑出来的方案经常出现“孤岛”和“环网”并存的结构，后来补上虚拟流约束才正常。

2.3 目标函数的层次：网损、弃电、开关操作代价

配电网重构的目标函数通常不是单一项。论文里最常见是三层相加。

第一层是网损，用全网所有支路的 (r_{ij}l_{ij}) 求时段和：

[ F_1 = \sum_{t=1}^{T}\sum_{(i,j)\in\mathcal{E}} r_{ij} l_{ij,t} ]

第二层是弃电惩罚。高比例清洁能源接入下，模型为了让电压不越限，允许削减部分分布式电源出力，但削减要付出代价，一般用削减量和单位惩罚成本的乘积。

第三层是开关操作代价。如果每个时段都重新重构，一天24小时开关来回切不现实，所以目标里带上相邻时段开关状态变化量的惩罚，系数一般取得很小，只起“温柔限制”作用。

这三项之间的权重系数，论文正文或者附录里一般都会给出。如果没给，复现时的默认做法是先只跑网损最小化，确定最优拓扑，再把其他项用较小的权重加进去，确保它们不改变主目标的最优解结构。这也是我和论文结果对不上时最常检查的地方。

3. 需求响应如何“长”进优化模型

3.1 弹性系数矩阵：最经典的价格型DR建模

需求响应建模有几十种，复现时先看论文用的是哪一种。最经典的是将电价变化映射到负荷变化。假设基准电价为 (\rho_0)，基准时段t的负荷为 (P_t^0)，电价变动量为 (\Delta \rho_k)，那么时段t的负荷变化量写作：

[ \Delta P_t = P_t^0 \sum_{k=1}^T e_{tk} \frac{\Delta \rho_k}{\rho_k} ]

其中 (e_{tk}) 是弹性系数矩阵的元素。当 (t=k) 时是自弹性，描述本时段电价对负荷的影响，一般为负值；当 (t\neq k) 时是交叉弹性，描述其他时段电价变化对当前时段负荷的转移效应，一般为正值。

这段代码里建模的关键是：这个 ( \Delta P_t ) 不是直接算出来的数值，而是要和电价变量、负荷变量一起放进优化模型里求解。也就是说，优化器可以同时决定最优电价激励和用户响应后的负荷水平。我在第一次复现时就犯了个错误，先把DR后的负荷曲线在外部算好，当成常数喂给重构模型，结果等于完全没把需求响应优化进来。

3.2 可平移、可削减的可行域约束

弹性系数模型描述了负荷随电价变化的“方向”，但还不足以构成数学约束，因为负荷不能无限平移。论文里通常还会给两个辅助约束：

一是调节幅度约束。需求响应后的负荷在基准负荷的一定比例范围内，比如 (0.7P_t^0 \leq P_t^{\mathrm{DR}} \leq 1.3P_t^0)。这个约束防止优化器把负荷调到零或者翻倍等不现实的值。

二是总用电量守恒约束。对可平移负荷来说，一天内用电总量近似不变，所以会有 (\sum_{t=1}^T \Delta P_t = 0)。如果论文里假设的是可削减负荷，则不需要这个硬等式，但会加上一个削减量上限，比如总削减不超过总负荷的5%。

这两个约束有没有写对，直接决定DR模型在求解器中是“削峰填谷”还是“瞎调”。我曾经遇到过一种情况：不加上电量守恒约束时，最优解把白天负荷全挪到半夜，网损确实降到很低，但完全不符合需求响应的物理语义。加上总电量守恒后，结果才合理。

3.3 DR在目标函数和约束中的“位置”决定求解复杂度

写模型时要清醒认识到：需求响应变量一旦引入，它不只在约束里出现，还经常和目标函数联动。比如需求响应后的负荷曲线改变了节点功率注入，进而改变DistFlow方程里的负荷项；同时，如果论文使用了激励型DR，那么还给参与DR的用户补偿成本，目标函数里就会多一项

[ F_{\mathrm{dr}} = \sum_t c_{\mathrm{dr}} (P_t^{\mathrm{DR}} - P_t^0) ]

复现时需要注意：这个成本和网损代价相比是很大还是很小？如果数量级不对，求解器要么让DR承担所有调节，要么完全不用DR。一般论文会给补偿单价，折合成标幺值后，我们也要跟着折算，不能直接拿原始市场电价代入。

4. 算法选型：为什么不用遗传算法和粒子群优化

4.1 启发式算法的痛点：组合爆炸与不可验证

配电网重构本质是选择哪些开关闭合、哪些断开，开关组合数是随网络规模指数增长的。局限在小网络时，遗传算法和粒子群优化看起来也能跑出结果，但遇到两个问题就难受了。

第一是每次迭代要对每个开关组合做一次潮流计算或牛顿法校验，计算成本很高。第二是启发式算法的结果依赖初始种群、变异概率、交叉概率，换一组随机种子可能得到完全不同的拓扑，很难说清解到底是不是最优的。对复现论文来说，如果参考的文献没有给出详细的算法参数，你很难把结果复现到同一水平。

4.2 从非凸到凸：二阶锥松弛带来的可求解性

把DistFlow的二次等式约束松弛成二阶锥约束之后，整个潮流区域是凸的。凸区域有很好的性质：局部最优就是全局最优，且可以用成熟的内点法高效求解。配合开关状态这个二进制变量后，模型变成混合整数二阶锥规划，YALMIP可以直接调用MOSEK这类商业求解器处理。

为什么很多论文选择这条路，而不是继续堆启发式算法？因为MISOCP模型有几个实际优点：结果可复现、有清晰的最优性间隙、软件生态成熟。你在跑同一个模型时，MOSEK给出间隙为0的解，换成Gurobi大概率也得到相同结果，这对学术复现非常有价值。

4.3 求解器选择与参数无关的小细节

这里必须提醒一个坑：MISOCP的求解时间对开关数量极其敏感。IEEE 33节点系统有32个分段开关加5个联络开关，如果每个时段都设一组二进制变量，24时段就是888个二进制变量，听起来不算多，但配上潮流连续变量和Big-M约束，MOSEK也要跑几分钟甚至更久。

我个人的习惯是：先做单时段优化，确认模型正确，再扩展到24时段。如果你一上来就多时段拉满，一旦模型写错，连定位问题的时间都被浪费了。另外，求解器参数要把相对间隙设到如 (10^{-4}) 级别，太大时开关状态基本不等，太小时个别网络容易卡死。

5. 从原理到代码：基于YALMIP和MOSEK的实现

5.1 数据准备：不用手打IEEE 33节点数据

IEEE 33节点系统的标准拓扑是33个节点、37条支路（含5条联络开关支路），基准电压12.66kV，基准功率10MVA。这个算例在各类配电网重构论文里出现频率极高，数据可以直接从MATPOWER的case33bw里读，不用手抄。

mpc = loadcase('case33bw'); branch = mpc.branch; bus = mpc.bus; % 折算到标幺值 Sbase = mpc.baseMVA * 1e6; % 通常10MVA Vbase = 12.66e3; % IEEE 33系统基准电压 Zbase = Vbase^2 / Sbase; br_r = branch(:,3) .* branch(:,1) ...

我这里写了个示意，实际使用时更简洁的做法是直接用p.u.表示，MATPOWER默认的R pu和X pu已经是标幺值。注意点只有一个：branch里开的联络开关在数据里是“0状态”，但初始拓扑的s0向量要按照论文给的初始状态来设置，别直接用MATPOWER里的通断状态当初始状态，这个是复现常见偏差来源之一。

5.2 决策变量与约束构建：从表达式到YALMIP

YALMIP的优势是能把数学表达式几乎原样翻译成代码。定义完sdpvar和binvar之后，核心工作是逐条写约束。

% 决策变量定义 u = binvar(E, NT); % 各时段开关状态 Pij = sdpvar(E, NT, 'full'); % 支路有功 Qij = sdpvar(E, NT, 'full'); % 支路无功 Ujm = sdpvar(N, NT); % 节点电压平方 I2 = sdpvar(E, NT); % 支路电流平方 fvf = sdpvar(E, NT, 'full'); % 虚拟潮流 % 约束集合 Constraints = []; for t = 1:NT for e = 1:E i = br_fb(e); j = br_tb(e); % DistFlow电压方程，开关断开时通过Big-M放松 Constraints = [Constraints, Ujm(i,t) - Ujm(j,t) - 2*(br_r(e)*Pij(e,t) + br_x(e)*Qij(e,t)) ... + (br_r(e)^2 + br_x(e)^2)*I2(e,t) == 0]; % 开关断开时，潮流/电压方程不成立，用大M处理 Constraints = [Constraints, Ujm(i,t) - Ujm(j,t) - 2*(br_r(e)*Pij(e,t) + br_x(e)*Qij(e,t)) ... + (br_r(e)^2 + br_x(e)^2)*I2(e,t) <= M*(1-u(e,t))]; Constraints = [Constraints, -(Ujm(i,t) - Ujm(j,t) - 2*(br_r(e)*Pij(e,t) + br_x(e)*Qij(e,t)) ... + (br_r(e)^2 + br_x(e)^2)*I2(e,t)) <= M*(1-u(e,t))]; % 二阶锥约束 Constraints = [Constraints, bounds([2*Pij(e,t); 2*Qij(e,t); I2(e,t)-Ujm(i,t)]) ... <= I2(e,t)+Ujm(i,t)]; end end

写约束时有两件事要特别留意。第一，Big-M的M值不能取太大，太大会造成数值病态，MOSEK警告信息会刷屏；一般取电压基准值的平方，比如电压上限1.05的平方乘上若干倍，大概100左右足够。第二，二阶锥约束在YALMIP里也可以用cone(2*Pij(e,t), 2*Qij(e,t), I2(e,t)-Ujm(i,t), I2(e,t)+Ujm(i,t))这种形式，效果一样，选一种自己顺手的。

5.3 需求响应代码块：把DR模型“插”进电力约束

需求响应变量接入的位置是节点功率平衡方程。以价格型DR为例，节点j的净注入要写成：

% 节点功率平衡，以节点j为例 % Pdr是DR优化后的有功负荷，Pd是原始负荷，Pdg是分布式电源出力 for t = 1:NT for j = 1:N Constraints = [Constraints, sum(Pij(find(br_fb==j),t)) - sum(Pij(find(br_tb==j),t)) ... == Pdr(j,t) - Pd(j,t) - Pdg(j,t)]; end end % 需求响应弹性约束，节点j为例 for t = 1:NT Constraints = [Constraints, Pdr(j,t) == Pd0(j,t) + Pd0(j,t)*sum(elasticMat(t,:).*dPrice./rho0)]; end % 调节范围和总电量约束 Constraints = [Constraints, 0.8*Pd0 <= Pdr <= 1.2*Pd0]; Constraints = [Constraints, sum(Pdr(j,:)) == sum(Pd0(j,:))];

这里有一个藏在细节里的麻烦：elasticMat的维度需要是 (T \times T)，(T) 是时段数。如果取24时段，这个矩阵就是24×24，每个元素都要有具体数值。论文通常会用表格给出自弹性、交叉弹性分区值，比如“峰时自弹性-0.15，交叉弹性0.03”。复现的时候千万不要自己瞎填，否则DR模型的行为会和论文完全对不上，表现为“负荷曲线被过度拉平”或者“几乎没反应”。

5.4 目标函数与求解设置：跑通只是一个开始

目标函数里加网损、弃电成本和开关动作惩罚。

obj_loss = sum(sum(repmat(br_r, 1, NT) .* I2)); obj_curtail = C_curtail * sum(sum(Pdgmax - Pdg)); obj_switch = C_switch * sum(sum(abs(u(:, 2:end) - u(:, 1:end-1)))); objective = obj_loss + obj_curtail + obj_switch; ops = sdpsettings('solver', 'mosek', 'verbose', 2); ops.mosek.MSK_DPAR_MIO_REL_GAP_TOL = 1e-4; optimize(Constraints, objective, ops);

跑通之后不要急着看网损数值，先看三件事：

• 开关状态是否满足“闭合线路总数等于N-1”；
• 每个节点电压是否落在合理区间；
• DR之后的负荷曲线是否出现了削峰填谷的形状。

如果这三样都对，再回头和论文对比数据。很多人复现时只盯着网损一个指标，结果开关组合是错的，网损数字碰巧很接近，就觉得代码没问题，这是最容易被审稿人抓包的情况。

6. 复现调试的几条真实心得

6.1 三种最常见的复现失败原因

我复现过程中踩过的坑，基本可以归为三类。

第一类是径向约束没写全。只加了“闭合线路总数等于N-1”而没有虚拟潮流约束，结果求解器为了降低网损，擅自把某个负荷孤立出去，让整个网络处于不连通状态。看结果时网损特别小，但只要画一下拓扑图就发现节点被甩在孤岛上。所以复现完一定要输出开关状态，自己画个图验证连通性。

第二类是Big-M参数设置不当。M取得太小会切断可行域，得到次优解；取得太大则数值不稳定。我的经验是，潮流约束相关的大M取100，虚拟潮流约束的大M取节点数N就足够，不要一刀切全用1000。

第三类是DR模型被“架空”。如果需求响应变量没有正确出现在节点功率平衡方程里，那么DR对负荷曲线的影响就完全不存在，但模型还能照常运行，表面看不出问题。这种情况下跑出来的开关组合和纯重构模型几乎一样，论文里那种“DR改变重构方案”的核心结论自然复现不出来。

6.2 结果和论文对不上时，先调什么

如果整体结构对了，但数值和论文差了一截，我建议按以下顺序排查：

1. 权重系数。检查网损、弃电、开关操作的量纲是否一致。很多论文里系数是统一到“万元”级别的，你在代码里如果用“元”或者“kW”为单位不折算，结果会差好几个数量级。

2. 初始开关状态。配电网重构的结果高度依赖初始拓扑。论文用了不同的初始开关状态，最终优化结果就会不同。这是复现中最容易被忽略的变量。

3. 分布式电源的接入位置和出力曲线。同一套IEEE 33节点系统，论文可能在节点12、22、29分别接了光伏、风电、储能，如果你接错位置，结果必然不对。

4. DR弹性矩阵数值。很多人喜欢用对称矩阵，但实际弹性系数往往不对称。比如工作时间电价上涨对早间负荷影响大，夜间交叉弹性很小。

6.3 参数灵敏度：让复现结果更有说服力的小技巧

最后一件事，我不建议只复现论文的主算例。主算例跑通后，做一组灵敏度分析会让你的复现更有价值。常见做法是：

将分布式电源渗透率从20%、30%调到40%，观察网损变化和重构次数变化；或者把DR弹性系数整体乘0.8、1.0、1.2，观察最优开关组合是否发生改变。

这类敏感性分析不需要大改代码，只需要把关键参数抽成一个变量，循环重算即可。做完之后你会发现，论文里那些结论往往不是某一个具体数值，而是一组“随着参数变化而变化的趋势”。能复现出这个趋势，才算真正吃透了文章。我最后补了一个“DR弹性系数→重构开关次数”的曲线，和论文图对比时，数值可能有小差距，但曲线形状完全一致，那一刻比单点数据完全吻合要让人踏实得多。