看到“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}^2x_{ij}^2)l_{ij} \end{cases} ]其中 (l_{ij} (P_{ij}^2Q_{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_{t1}^{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_{k1}^T e_{tk} \frac{\Delta \rho_k}{\rho_k} ]其中 (e_{tk}) 是弹性系数矩阵的元素。当 (tk) 时是自弹性描述本时段电价对负荷的影响一般为负值当 (t\neq k) 时是交叉弹性描述其他时段电价变化对当前时段负荷的转移效应一般为正值。这段代码里建模的关键是这个 ( \Delta P_t ) 不是直接算出来的数值而是要和电价变量、负荷变量一起放进优化模型里求解。也就是说优化器可以同时决定最优电价激励和用户响应后的负荷水平。我在第一次复现时就犯了个错误先把DR后的负荷曲线在外部算好当成常数喂给重构模型结果等于完全没把需求响应优化进来。3.2 可平移、可削减的可行域约束弹性系数模型描述了负荷随电价变化的“方向”但还不足以构成数学约束因为负荷不能无限平移。论文里通常还会给两个辅助约束一是调节幅度约束。需求响应后的负荷在基准负荷的一定比例范围内比如 (0.7P_t^0 \leq P_t^{\mathrm{DR}} \leq 1.3P_t^0)。这个约束防止优化器把负荷调到零或者翻倍等不现实的值。二是总用电量守恒约束。对可平移负荷来说一天内用电总量近似不变所以会有 (\sum_{t1}^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 决策变量与约束构建从表达式到YALMIPYALMIP的优势是能把数学表达式几乎原样翻译成代码。定义完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_fbj),t)) - sum(Pij(find(br_tbj),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 结果和论文对不上时先调什么如果整体结构对了但数值和论文差了一截我建议按以下顺序排查权重系数。检查网损、弃电、开关操作的量纲是否一致。很多论文里系数是统一到“万元”级别的你在代码里如果用“元”或者“kW”为单位不折算结果会差好几个数量级。初始开关状态。配电网重构的结果高度依赖初始拓扑。论文用了不同的初始开关状态最终优化结果就会不同。这是复现中最容易被忽略的变量。分布式电源的接入位置和出力曲线。同一套IEEE 33节点系统论文可能在节点12、22、29分别接了光伏、风电、储能如果你接错位置结果必然不对。DR弹性矩阵数值。很多人喜欢用对称矩阵但实际弹性系数往往不对称。比如工作时间电价上涨对早间负荷影响大夜间交叉弹性很小。6.3 参数灵敏度让复现结果更有说服力的小技巧最后一件事我不建议只复现论文的主算例。主算例跑通后做一组灵敏度分析会让你的复现更有价值。常见做法是将分布式电源渗透率从20%、30%调到40%观察网损变化和重构次数变化或者把DR弹性系数整体乘0.8、1.0、1.2观察最优开关组合是否发生改变。这类敏感性分析不需要大改代码只需要把关键参数抽成一个变量循环重算即可。做完之后你会发现论文里那些结论往往不是某一个具体数值而是一组“随着参数变化而变化的趋势”。能复现出这个趋势才算真正吃透了文章。我最后补了一个“DR弹性系数→重构开关次数”的曲线和论文图对比时数值可能有小差距但曲线形状完全一致那一刻比单点数据完全吻合要让人踏实得多。