前一阵子接到“风-水电联合优化运行”这个EI复现任务原本以为最费时间的会是代码实现真正跑起来以后才发现从论文公式到一份能稳定复现出结果的Matlab工程中间隔着大量没有写在纸上的建模细节和数据假设。这篇文章就围绕我完整走通的一条复现链路展开包括目标函数怎么定、水电出力特性怎么简化、风电不确定性怎么进模型、Matlab代码怎么搭、最后怎么验证“复现成功”给打算做EI复现或做电力系统优化调度的同学一个可以直接参考的模板。这套模型要解决的核心问题很直观风电出力有天然波动火电机组调节慢且成本高水电站依托水库具备快速响应能力。把三者放进同一个优化框架在满足负荷需求和水力约束的前提下让系统总运行成本最低、同时尽量少弃风少弃水这就是联合优化运行的全部内核。我下面的所有内容都会围绕这条主线展开不管是复现论文、准备毕业设计还是想评估水电对新能源的调节价值应该都能从中找到可以直接落地的思路。1. 复现前我先把这三个问题想清楚了写代码之前我花了大概两天时间反复读原文最后总结出三个必须先回答的问题。这三个问题如果没想清楚后面写多少代码都是白费。1.1 先搞清楚目标函数到底在优化什么很多复现失败不是代码问题而是连“优化什么”都没定义清楚。风-水电联合调度的目标函数在同类论文里有好几种写法有的最小化系统总运行成本有的最大化新能源消纳量还有的做成双目标去权衡经济性与碳排放。我复现的这篇以经济调度为基调目标函数取的是系统总运行成本最小化具体包含三项火电机组煤耗成本通常写成二次函数 a·P² b·P c弃风惩罚即预测风电没有被消纳的部分乘以一个惩罚单价弃水惩罚即水电站可以发但被迫放弃的水量乘以惩罚单价。对应的目标函数表达式就是min Σ时段 Σ机组 [ a_i · P_Ti² b_i · P_Ti c_i ] λ_w · Σ P_弃风 λ_h · Σ Q_弃水这里最容易被忽略的是惩罚系数λ_w和λ_h的取值。它们不是越大越好也不是越小越好。λ_w太小模型会为了省煤耗成本主动弃风违背“尽量消纳新能源”的调度倾向λ_w太大则会强迫系统在所有场景下全额消纳风电遇上极端风况可能直接无解或者把水电调度搞得过度激进。我的经验是先从火电边际成本的0.6~1.2倍取初值然后跑几组敏感性测试观察弃风量和总成本的变化拐点。如果某个λ取值让弃风量从10%骤降到接近0同时总成本跳升一大截那说明这个点附近就是合理的惩罚量级。1.2 约束条件的“硬”与“软”要分开处理优化模型里不是所有约束都处于同等地位。我的做法是把约束分成两类处理硬约束功率平衡、火电上下限、火电爬坡速率、水电站出力上下限、水库库容上下限。这些是物理上必须满足的违反任何一条求解器给出的解都不能用。软约束弃风、弃水。它们在目标函数里以惩罚项出现等于给“尽量消纳”一个优先级而不是必须做到。这样做的好处是模型在极端场景下依然存在可行解只是成本高一些如果非把弃风弃水都写成等式约束“必须为零”那随机场景耦合下经常得到infeasible的结果。功率平衡是等式约束必须严格满足这一点在代码里我会用容差判断去校验后面专门有一节讲。火电爬坡约束和水电调节速率约束是一对矛盾火电爬坡慢、水电爬坡快如果把火电爬坡设得太紧水电就不得不频繁调整出力去追负荷和风电波动。复现时我建议先给火电一个较宽松的爬坡范围比如每分钟1%~3%额定出力跑通后再收紧看影响这样能帮你判断论文里“水电发挥了调节作用”到底是模型内生出来的结论还是仅仅因为火电被约束绑死了。1.3 论文没给的数据怎么补全而不失真EI论文最让人头疼的就是关键数据经常不给你完整表格机组参数、负荷曲线、来水过程全在图片里只能靠肉眼读数。我当时为了凑齐一套能跑的系统数据用了三步处理系统规模采用3台火电、2座梯级水电站、1个风电场的简化系统。负荷曲线取典型日负荷峰值大致在600MW~800MW这个区间。水电参数补全每座水电站设置装机容量、库容上下限、最大发电流量、综合出力系数K。这些参数不完全来自论文部分参照同类型水电站的公开运行数据做了合理预估。敏感性检验对每个补全的参数做±20%的扰动重新求解观察最终调度结论是否翻轉。如果参数变化20%结论就变了说明模型对那个参数过于敏感要么参数估计有问题要么模型结构需要调整。这一步花的时间不少但非常值得。数据补全的过程本身就是对模型理解深度的检验——你不知道哪些参数重要就说明还没把论文吃透。2. 风-水电联合调度模型数学上到底堆了几层模型本身并不复杂但每一层都有需要注意的细节。我按从整体到局部的顺序拆一下。2.1 功率平衡与目标函数合在一起看任意时段t系统必须满足Σ 火电出力 Σ 水电出力 实际消纳风电 负荷这里“实际消纳风电”是一个变量它等于预测风电减去弃风P_w_use(t) P_w_forecast(t) - P_curtail(t)所以功率平衡可以写成Σ_i P_Ti(t) Σ_j P_Hj(t) P_w_use(t) P_load(t)火电煤耗成本按二次函数求和水电和风电的变动成本在目标函数里可以视为零。这种结构下优化器的逻辑就变成只要水电和风电还能出力就优先让它们顶上去火电只填补剩下的净负荷。2.2 水电站出力计算水头-流量非线性的处理水电站的出力表达式是P_H g · η · ρ · Q · H其中Q是发电流量H是净水头。麻烦在于H本身又是库容的函数水库水位随调度过程变化导致这个约束本质上是非凸的直接放进MILP框架很难处理。多数EI论文并不会展开这一层非线性而是采用简化表达。我复现时用的是固定水头近似调度周期只有24小时如果水库不是特别小库容变化引起的水头变化通常不超过几个百分点。把H近似成常数后出力表达式就变成P_H K · Q_flowK为综合出力系数单位是MW/(m³/s)。这样做有两个直接好处一是约束变成线性求解稳定性大幅提升二是物理意义依然清晰发电流量直接决定出力水量平衡依然把库容、入流、出流绑在一起。如果对精度要求更高可以改成分段线性化把库容区间分成三段每一段取不同的K系数再用0-1变量选择当前时段处于哪个区间。本质上就是在用混合整数线性规划拟合一族非线性曲线。我测试过分段线性化对最终成本的影响大约在1%以内但求解时间会成倍增加所以复现阶段取固定水头近似足矣。2.3 风电不确定性的场景化表达确定性模型直接用预测曲线等于假设预测完全准确。这个假设在真实调度里站不住脚也是很多论文复现后结果对不上原图的原因——原文明明在讨论不确定性你拿确定性模型去对结论当然对不上。我的方案是两层处理场景生成对预测风电曲线叠加正态分布误差误差标准差取预测值的10%~20%。用蒙特卡洛方法生成200个随机场景。场景削减用同步回代法backward reduction把200个场景削减到10个有代表性的场景每个场景附带一个概率权重π_s。随机规划下目标函数变成各个场景成本的加权期望min Σ_s π_s · cost_s关键是非预期性约束所有场景共享第一阶段的火电开机决策或者至少保证同一时段、同样信息状态下的决策一致。水电出力可以随场景调整因为水电的调节价值恰恰体现在“知道风电不准之后还能补救”。如果你把所有阶段的决策都按场景独立优化结果会过于乐观也违背了“决策必须基于当前可用信息”的调度逻辑。3. Matlab代码架构从脚本堆砌到能跑的工程这部分是我花时间最多的地方也是复现价值的直接体现。3.1 建模环境选型Yalmip Cplex我没有从零手写求解算法那不现实也没必要。我的技术栈是Yalmip做建模层Cplex或Gurobi做求解层。Yalmip在Matlab下声明优化变量、写约束、调用求解器都非常顺手比直接调Cplex的C API舒服太多。如果不用商业求解器Matlab自带的intlinprog也能跑但处理多场景随机规划时会明显变慢。10个场景、24个时段、几十个变量的规模下Cplex通常在几十秒内收敛intlinprog可能要跑到数分钟甚至更久。所以有条件的尽量上Cplex或Gurobi学生版授权也不难申请。3.2 核心代码段与逐段解释决策变量定义部分我用Yalmip写成这样%% 系统参数 T 24; % 调度时段数 nT 3; % 火电机组台数 nH 2; % 水电站数量 nW 1; % 风电场数量 tau 1; % 时段长度单位小时 %% 决策变量 P_T sdpvar(nT, T, full); % 火电出力 P_H sdpvar(nH, T, full); % 水电出力 Q_flow sdpvar(nH, T, full); % 发电流量 Q_spill sdpvar(nH, T, full); % 弃水流量 V sdpvar(nH, T1, full); % 库容注意多出一列 P_w_use sdpvar(nW, T, full); % 实际消纳风电 P_curtail sdpvar(nW, T, full); % 弃风功率这里有一个我一开始没注意的细节库容V定义成T1列因为水量平衡要从0时段推进到T时段最后需要T1个状态点。约束构建部分功率平衡和水量平衡是最核心的两块Constraints []; %% 功率平衡约束任意时段火电水电风电消纳负荷 for t 1:T Constraints [Constraints, sum(P_T(:,t)) sum(P_H(:,t)) sum(P_w_use(:,t)) P_load(t)]; end %% 水量平衡约束V(t1) V(t) (入流 - 发电流量 - 弃水) * tau for h 1:nH for t 1:T Constraints [Constraints, V(h,t1) V(h,t) (Q_in(h,t) - Q_flow(h,t) - Q_spill(h,t)) * tau]; end end %% 风电消纳关系消纳量 预测值 - 弃风 for t 1:T Constraints [Constraints, P_w_use(:,t) P_w_forecast(:,t) - P_curtail(:,t)]; end目标函数部分火电成本加上惩罚项Cost_thermal 0; for i 1:nT Cost_thermal Cost_thermal a(i) * P_T(i,:) * P_T(i,:) b(i) * sum(P_T(i,:)) c(i) * T; end Cost_curtail_wind lambda_w * sum(P_curtail(:)); Cost_spill lambda_h * sum(Q_spill(:)); Objective Cost_thermal Cost_curtail_wind Cost_spill;求解设置ops sdpsettings(solver, cplex, verbose, 2, savesolveroutput, 1); ops.cplex.mip.tolerances.mipgap 0.001; % 设置1e-3相对gap result optimize(Constraints, Objective, ops);这里有一个经验值mipgap设到1e-3足够了再往下压除了让求解器多跑几分钟对结果的实际影响很小。如果你的目标函数里存在非常接近的数量级差异可以适当放宽到5e-3换来求解时间的大幅下降。3.3 求解效率优化从一小时压到五分钟第一次把200个场景全部塞进去跑求解器直接卡住不动。后面我做了几件事把求解时间从一小时量级压到了五分钟左右。第一尽量矩阵化减少for循环。Yalmip里变量本身支持向量化求和像sum(P_T, 2)这样一次算完所有时段的某一台机组出力能少写循环就能少给求解器增加不必要的内存占用。第二做量纲统一。所有机组出力和负荷统一用MW所有水流量统一用万m³/h。如果火电煤耗系数是10⁻⁴量级水库库容却是10⁷量级目标函数和约束里的数值会差12个数量级这对Cplex/Gurobi的数值稳定性非常不友好会出现ill-conditioned的警告甚至莫名其妙的不可行解下一章节细说。第三把二次煤耗成本分段线性化。如果直接用二次函数模型变成MIQP虽然也能解但速度比纯MILP慢很多。我试验过后直接在论文中采用了分段线性化处理比如把出力区间切成四段每段线性化成本斜率整体变成纯MILP。这样既能保持精度又能让求解速度大幅提升。优化前后的对比大概是这样方案求解时间目标函数值说明直接MIQP10场景约40分钟仅作参考太慢调参困难MILP分段线性10场景约4分钟相对误差0.5%推荐MILP场景削减到10个约2分钟相对误差1%日常快速试验用4. 复现结果到底算不算“复现成功”复现不是“跑出一个好看的图”就结束了。我给自己定了三条判断标准。4.1 逐时段的物理约束校验输出方案以后第一件事不是看图而是写一段校验脚本把每个时段的功率平衡残差、水量平衡残差、库容上下限、爬坡约束逐个检查一遍。功率平衡的残差应该小于1e-6 MW量级水量平衡残差应该小于1e-6万m³量级。如果残差偏大基本可以判定建模或求解环节有问题。我习惯的做法是res_power sum(P_T, 1) sum(P_H, 1) P_w_use - P_load; max_abs_res max(abs(res_power)); fprintf(最大功率平衡残差: %.2e\n, max_abs_res);如果这里出现明显残差先回头检查约束有没有漏写再检查求解器设置是不是把容差放得太宽了。Cplex默认容差通常够用但如果你手动改过eprhs或mipgap要留意它们对解的影响。4.2 与论文图表对比的量化方法论文里的曲线是图像没有原始数据所以对比只能做数值化处理。我的对比方式分三层目标函数值对比同规模、同时段下论文的总成本和我的目标函数值应在正负10%以内。超过这个范围说明某些参数或者约束条件设置存在较大偏差。曲线形状对比计算火电、水电、风电消纳曲线与论文图的RMSE和Pearson相关系数。相关系数在0.9以上基本可以认为趋势吻合。典型模式验证看水电出力是否在负荷高峰时段显著增加、在低谷时段减小火电则跟随净负荷趋势变化。如果水电出力曲线是平的说明水库的调节潜力没有被模型发挥出来约束设置或者目标函数大概率有问题。4.3 压力测试边界条件下的稳定性最后一步也是我认为最容易被复现者跳过的一步压力测试。调整输入场景观察模型是否稳定给出合理的调度方案。极端风电场景把所有时段的风电预测设为零模型应该自动增加火电和水电出力度过“无风期”把所有时段设为额定出力模型应该能靠水库调节能力消纳多余风电实在消纳不了时才弃风。枯水场景把来水大幅降低水电出力受限火电应当顶上库容约束不出现越界。丰水场景来水很足模型应当优先让水电多发同时警惕弃水惩罚设置是否合理。如果模型在边界场景下经常报infeasible或者给出明显不合常理的调度结果问题往往不出在求解器而出在约束条件本身存在互相矛盾的地方。5. 我这一趟踩过的坑按重现顺序梳理踩坑是复现过程中最有价值的部分写出来帮大家省时间。5.1 库容时间索引的错位导致水电站“凭空抽水”第一次跑通模型后我检查水量平衡发现第一个时段的库容变化和出入流对不上库容凭空跳了一大截。排查了半天最后发现是索引错位我把水量平衡写成了V(h,t) V(h,t-1) ...然后在t1时引用了V(h,0)Matlab里索引从1开始导致这个约束被默认忽略或者读取了错误的初始值。正确写法是V(h,t1) V(h,t) ...让时间索引从初始库容逐步向前推进。排查这个问题的办法就是在校验脚本里逐时段打印水量平衡各项不要只打印残差总量。残差为零可能是约束没生效打印每一项才能暴露索引问题。5.2 场景削减过猛导致无可行解最初我把200个场景削减到5个结果求解器报了infeasible。原因在于场景数量太少时爬坡约束和功率平衡约束在少数几个场景里互相牵扯解空间就被锁死了。200个场景虽然解大但冗余场景提供了更多的缓冲余地。削减到5个后每个场景都必须被满足约束过紧。解决方法是削减后的场景数不要低于10个并且削减完成后对权重做归一化重新检查场景矩阵的条件数。如果条件数异常说明保留的场景之间相关性太高代表性不足。5.3 求解器的数值病态量纲不统一害死人有一段时间Cplex老是警告ill-conditioned解出来的结果时好时坏有时候弃风量突然变成负数。后来我把所有量纲拉出来对比发现火电煤耗系数在10⁻⁴量级水库库容在10⁶量级惩罚系数λ又取了比较大的数整个约束矩阵的条件数已经差到十分恶劣的程度。解决办法就是前面提到的量纲归一化把出力统一到MW库容统一到万m³所有惩罚系数和煤耗成本保持在相近数量级。改造之后求解器警告彻底消失结果也稳定了。这事给我的教训是代码里优先写好量纲注释一劳永逸。写在最后的体会我整个复现过程做完的体会是好的复现不是把公式抄一遍而是要在“模型保真度”和“可求解性”之间找平衡。固定水头近似、煤耗分段线性化、场景削减到10个这些看上去都是让模型“变简单”的操作却能让求解稳定性和结果可用性大幅上升反而不是坏事。另一个比较深的感受是数据准备和结果验证几乎占掉了六成以上时间真正写模型框架反而是快的那部分。建议后来者从确定性模型先跑通再加场景和不确定性一步到位出问题的可能性会高很多。如果你之后想把这套模型往更实用的方向扩展可以考虑把调度周期从24小时拉长到几天把末库容从固定值改成优化边界让水库在更长的时间尺度上做跨日调节。改动量不大但结论的丰富程度会明显不一样。