1. 为什么配电网调度需要“两阶段”而不是一步到位我在做配电网优化调度时最早用的是单阶段确定性模型——给定负荷预测和分布式电源出力预测直接求解一个最优潮流问题得出各机组和联络线的功率计划。这种思路在分布式电源渗透率低的场景下够用但一旦光伏、风电接入容量上去问题就来了预测误差导致的偏差会在实际运行中被放大甚至让日前计划完全不可执行。两阶段优化调度就是把“决策”拆成两个时间尺度第一阶段是在日前前一天基于预测信息确定机组启停、储能充放电计划、联络线交换功率等需要提前锁定的决策量第二阶段是在日内或者实时根据更准确的预测甚至实测数据对可控分布式电源、储能、需求响应等资源做修正调整消除预测误差带来的功率偏差和电压越限。说白了第一阶段管“明天大概怎么跑”第二阶段管“实际跑偏了怎么拉回来”。这种“先计划、后修正”的思路本质上是对不确定性的一种结构化处理。配电网里的不确定性来源很多光伏出力受云层遮挡影响分钟级波动能达到装机容量的30%到50%风电更不用说爬坡事件经常让调度员措手不及负荷侧的电动汽车、空调等随机性也在增加。把这些不确定性全部塞进一个确定性模型里结果必然保守或者冒进。两阶段模型的价值在于让“必须提前定的”和“可以等一等再定的”分开处理既保证了日前计划的可行性又给日内调整留了足够的自由度。从Matlab代码实现的角度看两阶段模型的难点不在于优化理论本身而在于如何把两阶段决策变量、约束条件和目标函数组织成求解器能接受的形式。很多同学拿到这个问题第一反应是直接写一个大规模混合整数线性规划扔给求解器结果不是内存爆炸就是求解时间不可接受。这篇文章我会把模型拆开讲清楚再给出一套在Matlab里可落地的代码框架包括数据准备、Yalmip建模、求解器配置和结果分析最后聊几个我在实际调试中踩过的坑。2. 日前两阶段优化调度模型的数学描述目标函数和约束是怎么定的2.1 目标函数经济性和安全性的权衡两阶段模型的目标函数通常写成第一阶段成本加上第二阶段期望成本的形式。第一阶段成本包括常规机组或者配电网里的微型燃气轮机、柴油发电机的启停成本和燃料成本联络线购电成本以及储能充放电的折旧成本。第二阶段成本则是在不确定性 realization 出来之后对调整量施加的惩罚成本。用数学语言写就是[ \min \quad \sum_{t \in T} \left[ C_{start} u_t C_{fuel}(P_{g,t}) C_{grid} P_{grid,t} \right] \sum_{t \in T} E_{\xi} \left[ Q(P_{g,t}, P_{ess,t}, \xi_t) \right] ]其中 ( C_{start} ) 是启停成本( u_t ) 是机组启停状态( P_{g,t} ) 是机组出力( P_{grid,t} ) 是联络线功率( \xi_t ) 代表不确定性比如光伏实际出力与预测的偏差( Q(\cdot) ) 是第二阶段最优调整成本函数。在实际代码里目标函数不会写这么抽象。通常的做法是第一阶段决策变量包括机组启停 ( u_t )、机组出力 ( P_{g,t} )、储能充放电功率、联络线功率第二阶段决策变量包括针对每个场景或者最坏情况的调整量比如切负荷量、弃光量、储能修正充放电量。目标函数里把这些调整量乘以一个较大的惩罚系数迫使优化结果尽量避免切负荷和弃光。这里有一个关键点惩罚系数的量级要设置合理。如果罚因子太小优化结果可能会出现弃光、切负荷来“省钱”的不合理现象如果太大又会数值病态影响求解器收敛。我一般把切负荷罚因子设为电价的10到20倍弃光罚因子设为分布式电源上网电价的5到10倍具体值需要根据算例调试。2.2 第一阶段约束网络拓扑和运行边界第一阶段的约束主要是配电网的稳态潮流约束以及各类设备的运行边界。潮流约束是配电网优化里最费时间的部分。对于辐射状配电网最常用的是DistFlow方程[ P_{i,t} \sum_{j \in child(i)} P_{j,t} r_{ij} l_{ij,t} P_{load,i,t} - P_{dg,i,t} ] [ Q_{i,t} \sum_{j \in child(i)} Q_{j,t} x_{ij} l_{ij,t} Q_{load,i,t} - Q_{dg,i,t} ] [ V_{i,t}^2 V_{j,t}^2 - 2(r_{ij} P_{ij,t} x_{ij} Q_{ij,t}) (r_{ij}^2 x_{ij}^2) l_{ij,t} ]其中 ( P_{ij,t} )、( Q_{ij,t} ) 是支路有功和无功功率( r_{ij} )、( x_{ij} ) 是支路电阻和电抗( l_{ij,t} ) 是电流幅值平方( V_{i,t} ) 是节点电压幅值。这些约束是非线性的因为 ( l_{ij,t} ) 和 ( V_{i,t} ) 是平方项在Matlab里通常会用二阶锥松弛SOCP把它转成凸约束[ l_{ij,t} \geq \frac{P_{ij,t}^2 Q_{ij,t}^2}{V_{i,t}^2} ]这样处理后模型就变成了混合整数二阶锥规划MISOCP可以用Yalmip配合Gurobi、CPLEX或者Mosek求解。设备约束方面常规机组有出力上下限和爬坡约束[ u_t P_{g,min} \leq P_{g,t} \leq u_t P_{g,max} ] [R_{down} \leq P_{g,t} - P_{g,t-1} \leq R_{up} ]储能约束包括充放电功率限制、容量状态SOC递推和SOC上下限。这里容易犯的错是把充电和放电用两个独立变量建模导致求解器同时充放电来“套利”。正确做法是引入一个二进制变量 ( b_{ch,t} ) 表示充电状态约束 ( P_{ch,t} \leq b_{ch,t} P_{ch,max} ) 和 ( P_{dis,t} \leq (1-b_{ch,t}) P_{dis,max} )从模型层面杜绝同时充放电。联络线功率约束则要反映配电网与上级电网的交换容量限制必要时还要考虑峰谷电价下的购电策略。2.3 第二阶段约束调整量和安全边界第二阶段的决策变量是各场景下的调整量。常见的建模方式有两种一种是随机规划里的多场景方法对每个场景都求解一个第二阶段子问题另一种是鲁棒优化里的盒式不确定集合求解最坏情况下的调整量。多场景方法相对容易实现代码上就是给第二阶段变量增加一个维度——场景维度。比如储能修正充放电功率写成 ( P_{ess,ch,s,t} )其中 ( s ) 是场景索引。这样第二阶段约束就是每个场景独立满足功率平衡、储能SOC递推和安全约束。关键的地方在于第一阶段变量对所有场景是共享的不能随场景变化这就是所谓的“非预期约束”non-anticipativity constraint。在代码实现里通过把第一阶段变量定义成不含场景维度的变量即可Yalmip会自动处理这种耦合关系。但要注意如果场景数量太多比如超过50个模型规模会急剧膨胀求解时间指数级上升。我在实际项目中一般会先用K-means聚类把原始场景缩减到10到20个代表性场景再带入优化模型效果和精度损失非常小。安全边界方面第二阶段要保证在最坏情况下节点电压不越限、支路电流不越限、不切负荷或者切负荷量在可接受范围内。这些约束在代码中表现为每个场景下都要满足的潮流约束和电压上下限约束。3. Matlab代码架构从数据准备到Yalmip建模全流程3.1 数据准备负荷、分布式电源和网络参数在写优化模型之前要把基础数据整理干净。我习惯用结构体或MATLAB表格统一管理避免散落一堆变量名。核心数据包括24小时或96个时段的负荷预测曲线光伏和风电的预测出力曲线标幺值或实际值网络拓扑参数节点数、支路数、支路阻抗、节点负荷比例储能参数容量、最大充放电功率、初始SOC、效率常规机组参数出力上下限、爬坡率、燃料成本系数分时电价以及分布式电源的接入位置和容量。这里要给新手一个建议数据格式统一用列向量长度等于时段数网络参数用矩阵存储第一列是支路首端节点第二列是支路末端节点第三、四列是电阻和电抗。这样在用Yalmip定义变量和约束时可以直接用矩阵索引代码简洁且不容易出错。3.2 用Yalmip构建两阶段模型的代码骨架Yalmip是Matlab里做优化建模最顺手的工具没有之一。它的好处是语法接近数学表达式不需要手动把约束展开成大矩阵。下面给出一个两阶段模型的精简代码骨架完整版需要根据算例数据补充。%% 定义时间尺度和场景数 T 24; % 日前调度时段数 N_s 10; % 场景数缩减后 N_bus 33; % 节点数 N_branch 32; % 支路数 %% 第一阶段决策变量 u_g binvar(N_g, T, full); % 机组启停状态 P_g sdpvar(N_g, T, full); % 机组出力 P_grid sdpvar(1, T, full); % 联络线功率 P_ch sdpvar(N_ess, T, full); % 储能充电功率 P_dis sdpvar(N_ess, T, full); % 储能放电功率 b_ch binvar(N_ess, T, full); % 储能充电状态 SOC sdpvar(N_ess, T1, full); % 储能荷电状态 %% 第二阶段决策变量带场景维度 P_curtail sdpvar(N_dg, T, N_s, full); % 分布式电源弃电功率 P_load_shed sdpvar(N_load, T, N_s, full); % 切负荷功率 P_ch_s sdpvar(N_ess, T, N_s, full); % 场景下储能充电修正 P_dis_s sdpvar(N_ess, T, N_s, full); % 场景下储能放电修正 %% 目标函数 objective 0; for t 1:T objective objective C_start * u_g(:,t) fuel_cost(P_g(:,t)) ... price(t) * P_grid(t) ... % 购电成本 lambda_curtail * sum(P_curtail(:,t,:), all) ... % 弃光/弃风惩罚 lambda_shed * sum(P_load_shed(:,t,:), all); % 切负荷惩罚 end目标的求和里第二阶段惩罚项用sum(..., all)把所有场景的惩罚加总。由于场景是有概率权重的严格来说应该加权重 ( \pi_s )这里为了简化展示没写完整代码里记得把场景概率乘进去。约束的写法类似核心是把第一阶段约束和第二阶段约束分开定义constraints []; %% 第一阶段约束 for t 1:T % 机组出力上下限 constraints [constraints, P_g(:,t) u_g(:,t) * P_g_min]; constraints [constraints, P_g(:,t) u_g(:,t) * P_g_max]; % 储能同时充放电约束 constraints [constraints, P_ch(:,t) b_ch(:,t) * P_ch_max]; constraints [constraints, P_dis(:,t) (1 - b_ch(:,t)) * P_dis_max]; % SOC递推 constraints [constraints, SOC(:,t1) SOC(:,t) eta_ch * P_ch(:,t) ... - P_dis(:,t) / eta_dis]; end %% 第二阶段约束每个场景 for s 1:N_s for t 1:T % 节点功率平衡DistFlow线性化版本——这里以根节点功率平衡示意 constraints [constraints, P_grid(t) sum(P_g(:,t)) sum(P_dis(:,t)) ... sum(P_ch(:,t)) P_load_total(t) - sum(P_curtail(:,t,s))]; end end %% 求解 ops sdpsettings(solver, gurobi, verbose, 2); optimize(constraints, objective, ops);这个骨架省略了潮流约束的完整展开实际项目里需要用DistFlow的SOCP松弛形式逐支路写约束。Yalmip里二阶锥约束可以写成cone([2*P_ij; 2*Q_ij; V_i - l_ij], V_i l_ij)的形式非常方便。3.3 求解器选择Gurobi还是CPLEX两阶段模型最终求解的是MISOCP问题。Matlab自带的linprog、intlinprog只能处理线性和混合整数线性问题无法直接处理二阶锥约束。所以必须外接求解器。我试过Gurobi、CPLEX和Mosek三个求解器结论是Gurobi在二阶锥规划和混合整数问题上的速度综合最优尤其在大规模场景下优势明显。CPLEX的稳定性好Mosek在纯凸优化上很强但在混合整数二阶锥上不如前两者。如果是学生做学术研究Gurobi有免费学术许可申请流程简单推荐优先尝试。Yalmip配置Gurobi很简单安装好Gurobi后在Matlab里运行yalmiptest检测一下然后求解时指定sdpsettings(solver, gurobi)即可。注意Gurobi的Matlab接口需要把Gurobi安装目录下的matlab文件夹加入Matlab路径否则会报找不到gurobi_mex的错误。3.4 结果可视化调度曲线和电压分布求解完成后第一件事是检查求解状态和目标函数值确认模型没有infeasible或unbounded。然后提取各决策变量画出调度曲线。我习惯画三张图第一张是24小时的功率平衡图把负荷、光伏出力、风电出力、机组出力、联络线功率、储能充放电画在同一个坐标系里一眼就能看出功率是否平衡第二张是各节点电压幅值分布图检查是否有电压越限第三张是储能SOC曲线确认储能没有出现频繁启停或者深度充放的异常行为。画图代码没什么特别的但有一点要提醒Yalmip求解完成后变量是sdpvar对象要用value()函数提取数值直接使用会得到符号对象画图时报错。4. 实际调试中的四个大坑求解失败、收敛慢和结果不合理4.1 二阶锥松弛不紧电压偏低但目标函数还在下降这是我在做33节点配电网算例时遇到的第一个大坑。加了SOCP松弛后模型求解速度确实快但画出来的电压分布明显不合理——某些节点电压掉到0.85标幺值以下而实际配电网不会出现这种情况。问题出在二阶锥松弛不紧。DistFlow的线性化版本——把 ( V_i^2 ) 替换成新的变量 ( U_i )把 ( l_{ij} ) 也作为变量——只有当目标函数对 ( l_{ij} ) 和 ( U_i ) 的系数满足一定条件时SOCP松弛才是紧的即松弛后的解等于原非线性问题的解。在我那个算例里目标函数中不含线损项导致锥松弛在某些支路上不收紧出现了“虚假”的低电压解。解决方法是目标函数中加入一个极小的线损惩罚项或者对松弛约束加上一个小的罚因子。我当时的做法是在目标函数中加上 ( \alpha \sum r_{ij} l_{ij,t} )其中 ( \alpha ) 取1e-4。加了之后锥松弛被拉紧电压恢复到合理范围。这个技巧在很多文献里叫“penalty convex-concave procedure”但在简单场景下一个小的线损惩罚项就够了。4.2 储能同时充放电一个被忽略的二进制变量储能同时充放电是配电网优化里频率极高的bug。我最早做的时候用P_ess一个变量表示储能功率正值放电、负值充电约束只有 ( -P_{ch,max} \leq P_ess \leq P_{dis,max} )。结果是每次优化完储能都会出现在电价低谷时“充电”到SOC上限、电价高峰时“放电”到SOC下限这没问题但某些时段会出现 ( P_ess ) 在正负之间来回跳变看起来像是在同时充放电。仔细检查后发现这是因为SOC递推约束里效率系数 ( \eta_{ch} ) 和 ( 1/\eta_{dis} ) 的存在使得同一时刻的充电和放电在SOC层面无法完全抵消——但目标函数里电价又鼓励这种行为所以求解器会在离散点上制造同时充放电的伪现象。解决方法是引入二进制变量 ( b_{ch} )并加上互斥约束。代价是模型复杂度上升因为每个时段每个储能多了一个整数变量。对于大规模系统如果整数变量太多导致求解时间不可接受可以用Big-M方法或者线性化互补约束来替代但互斥二进制变量是最稳妥的做法。4.3 场景缩减别拿100个场景直接扔给求解器多场景随机规划里场景数量直接决定模型规模。100个场景、96个时段、33个节点的配电网变量数轻松超过10万Gurobi再快也要跑到猴年马月。我建议的做法是先用场景削减算法比如K-means聚类、快速前向选择算法把原始场景缩到10到20个。K-means虽然简单但需要先把每个场景的时序曲线展平成高维向量然后用kmeans()函数聚类再用聚类中心作为代表场景场景概率就是该类场景的数量占比。这样做之后优化结果和保留全部场景的解差异非常小一般目标函数偏差在1%以内但求解时间从几十分钟降到几分钟性价比极高。还有一个小技巧如果在Yalmip里变量的场景维度太大导致内存不足可以尝试把第二阶段问题的约束写成稀疏矩阵形式而不是直接用sdpvar的3维数组。3维数组在Yalmip里展开成约束矩阵时会有大量的稀疏填充开销。4.4 求解器报“Numerical issues”的处理思路Gurobi报numerical trouble是最让人头疼的通常表现为求解过程反复振荡、收敛极慢甚至出现错误的infeasible判定。我的排查顺序是第一步检查量纲。配电网模型里功率单位用MW电压用kV阻抗用Ω数量级通常还在可控范围。如果用了kW和kV混用目标函数里的成本项和约束系数可能差好几个数量级求解器数值稳定性会大幅下降。我的习惯是全部用标幺值pu建模这样所有变量数量级都在0.01到1之间最稳定。第二步检查惩罚系数。切负荷罚因子如果设到1e6而其他成本都是100以内会出现严重的数值病态。我会把目标函数的各分量先打印出来看看量级确保最大项和最小项之比不超过1e4。第三步检查约束是否有冗余或者冲突。比如同时约束了联络线功率等于某值和大于某值这类矛盾约束在数值上会让求解器无所适从。我一般会在构建完约束集合后用Yalmip的check(constraints)命令检查每个约束的残差——虽然求解前的check主要是看有没有语法错误但至少能发现明显的维度不匹配。5. 两阶段模型的效果验证从目标函数到调度策略的三个对比实验5.1 对比基线单阶段确定性模型验证两阶段模型的价值最直接的方式是和单阶段确定性模型做对比。单阶段模型直接用预测曲线求解得到一组机组出力和储能计划两阶段模型则额外考虑了场景偏差的调整能力。对比的指标包括总运行成本、弃电量、切负荷量和电压越限次数。我在IEEE 33节点系统上做过这个对比结果很有代表性单阶段模型的名义成本比两阶段模型低3%左右但把预测误差对应的实际场景回代进去时单阶段方案出现了明显的功率不平衡需要额外切负荷或者弃光实际总成本反而比两阶段模型高6%到8%。这说明两阶段模型虽然名义成本高一点但面对不确定性时更稳健。5.2 场景数量对结果的影响我自己跑过一组实验保留5、10、20、50个场景对比求解时间和目标函数值。结果如下表场景数求解时间秒目标函数值万元相对偏差51252.31.8%103551.60.4%2010851.50.2%5062051.4基准可以看到从5个场景增加到10个场景目标函数改善了1.4%但从10增加到20只改善了0.2%。考虑到求解时间从35秒涨到108秒10个场景是性价比最高的选择。当然不同系统结论可能不同但“先跑5个场景看趋势再逐步增加”的调参思路是通用的。5.3 储能容量与调度策略的联动分析两阶段模型还有一个好处它能更真实地评估储能在不确定性环境下的价值。因为在单阶段确定性模型里储能的作用主要靠峰谷价差套利而在两阶段模型里储能额外承担了平抑预测误差的任务所以储能容量增大带来的收益不仅体现在电价套利上还体现在减少切负荷惩罚上。我在某实际工业园区配网数据上测算过储能容量从1MWh增加到3MWh时套利收益增幅放缓但切负荷惩罚成本显著下降整体收益依然可观。如果只算峰谷价差可能得不出“储能值得投资”的结论。这其实就是两阶段模型对规划问题的指导意义——运行优化和投资决策本来就不该用确定性逻辑来简化。6. 从算例到实际工程模型扩展的三个方向6.1 多时段耦合约束把冷热电联供纳入调度配电网里经常不只是电负荷还有热负荷、气负荷。如果园区里有冷热电联供CCHP系统两阶段模型需要扩展热电联产约束——发电量和产热量之间存在可行域约束不能独立调节。Matlab实现时可以在第一阶段变量里增加热出力变量 ( H_{cchp,t} )并加上可行域约束[ P_{cchp,t} \in [P_{min}(H_{cchp,t}), P_{max}(H_{cchp,t})] ]这个可行域通常是一个多边形可以用一组线性约束近似。我实际项目中是把厂商提供的运行曲线离散成多个顶点然后用凸组合形式表述Yalmip里可以直接定义。6.2 需求响应资源的柔性约束需求响应DR是配电网调度里越来越重要的灵活性资源。两阶段模型中DR适合放在第二阶段作为调整手段——第一天签订基线负荷和响应容量日内根据实际偏差调用。建模时负荷变成一个区间变量 ( P_{load,i,t} \in [P_{load,i,t}^{base} - \Delta P_{DR,i,t}, P_{load,i,t}^{base} \Delta P_{DR,i,t}] )同时约束全天的响应电量不超过合同电量。这个变量如果放在第一阶段会因为负荷不确定性而变得过于保守放在第二阶段则可以灵活地“吸收”分布式电源预测误差。代码实现上只需在第二阶段约束中把负荷变量改成区间变量并在目标函数中加入DR调用成本。6.3 与实时调度的闭环衔接两阶段模型解决的是日前计划问题但它输出的第二阶段调整量实际上可以给实时调度提供初始条件。我在工程中做过一个方案每天都跑日前两阶段模型生成机组启停计划和储能SOC参考轨迹实时调度层每5分钟一个周期把这些计划作为硬约束只优化储能和可控负荷的微小调整量。这样既保证了日前计划的全局最优性又满足了实时控制的快速性要求。这个架构的好处是日前模型里第二阶段场景对应的调整量实际上为实时层提供了“预案”信息。比如某场景下预测到午后光伏出力骤减日前模型已经算好了对应的储能放电方案实时层只需按预案执行不需要从零开始优化计算时间可以控制在几秒以内。7. 最后聊几点个人体会这套两阶段模型从理论到Matlab落地我前后折腾了不少时间总结几条经验供参考。第一别在一开始就追求模型的完美复杂度。先把单阶段确定性问题跑通再加需求响应、加多场景、加SOCP松弛每一步都验证约束是否合理。直接堆一个巨复杂模型出bug了根本无从排查。第二Yalmip对这个领域帮助极大但它的调试信息有限。如果模型infeasibleYalmip不会告诉你哪条约束错了。我常用的方法是二分法删约束——先把所有约束注释掉然后每次加一半约束看模型什么时候从feasible变成infeasible就能锁定问题约束所在区域。第三数据质量对优化结果的影响远超模型本身的改进。光伏预测和负荷预测的数据如果不做清洗有异常尖峰或者缺失优化结果会非常离谱。我在跑任何模型前都会把数据画出来看一遍这比任何算法调优都重要。第四两阶段模型的代码框架其实可以复用到很多其他问题上——微电网优化调度、综合能源系统运行优化、电动汽车有序充电本质都是“先计划、后修正”的结构。把这套框架吃透换一套数据就能迁移过去性价比很高。希望这篇文章能帮你把“含分布式电源的配电网日前两阶段优化调度模型”从论文里的公式真正变成能跑的Matlab代码。如果在实现过程中遇到问题沿着前面几个调试方向排查大概率能解决问题。