做风-水电联合优化的人应该都有同感单独看风电随机性、间歇性大到让人头疼单独看水电又受来水、水库调节能力和生态流量约束。但把两者放在同一个调度框架里做联合优化运行分析往往能在几乎不新增硬件投资的情况下把系统的灵活性和可再生能源消纳空间拉高一大截。用Matlab来做这块的代码复现是我个人觉得性价比最高的路径。这篇文章就围绕“风-水电联合优化运行分析”这个标题把EI论文里常出现的建模思路、目标函数与约束条件怎么落地成Matlab代码、求解器怎么选、结果怎么验证、以及复现过程中最容易踩的坑全部过一遍。适合正在复现论文、做课程设计、或者准备写小论文的电力方向研究生也适合想快速上手调度优化模型的工程师参考。1. 先说清楚风-水电联合优化运行到底在解决什么问题1.1 为什么风电场要跟水电站“绑”在一起风电出力靠天吃饭风速曲线长什么样出力就长什么样基本没有太多人工干预的空间。水电则不同只要水库有足够调节库容机组出力可以在一定范围内快速调节。把这两个电源放在一个调度模型里核心思路是用可控性好、响应速度快的水电去平抑风电的波动和预测误差。风电大发的时候水库少放水、多蓄水风电小发的时候水库多放水、补足缺口。这样一来电网调度看到的是一条更平稳的外送功率曲线原先不得不弃掉的风电也能被吸收进系统。这在术语里叫“互补运行”但你要真去复现论文就会发现难点不在概念而在模型。水文过程有水量平衡风电有预测误差两者耦合在一起之后约束条件呈指数级膨胀调度时段一旦取到96点甚至288点普通Excel根本算不动必须依靠Matlab这样的数值计算环境。1.2 EI论文里最常见的三类模型我在复现不同文献时发现风-水电联合优化模型的差别主要在三处时间尺度、目标函数、不确定性处理方式。给你一个速查表模型类型典型时间尺度目标函数常见求解方法长期调度月/旬发电量最大、蓄能最大DP、逐次逼近短期调度日/96点运行成本最小、弃风最小MILP、NLP实时/滚动15min级跟踪计划偏差最小MPC、鲁棒优化真正的EI复现绝大多数集中在“短期调度”这一档。再往上文献会在这个框架里叠加风电不确定性产生两类变体场景法用蒙特卡洛抽样生成大量风电出力场景然后用场景削减技术浓缩成几个典型场景模型变成随机规划。鲁棒法不关心具体概率分布只给风电出力一个波动区间模型追求最坏情况下的调度可行性。这两条路线在Matlab里的实现路径差别很大后文我会分别讲怎么处理。2. 建模环节从物理对象到数学表达2.1 目标函数怎么选选目标函数首要判断依据是你复现的论文里系统调度主体是谁追求的经济指标是什么。如果是独立发电企业视角目标函数通常是“系统运行收益最大”即水电上网收益加上风电上网收益减去运行维护成本、弃风惩罚成本。公式上写作max F sum( sum( P_h(i,t) * price_h(t) * dt ) sum( P_w(i,t) * price_w(t) * dt ) ) - sum( curtail_w(i,t) * penalty * dt )如果是电网调度视角目标函数更多写成“系统运行成本最小”这时水电出力成本用可变成本函数近似风电成本近似为0或者说优先消纳目标则变成减少火电出力和弃风惩罚min C sum( sum( a_i * P_g(i,t)^2 b_i * P_g(i,t) c_i ) ) sum( lambda_cur * curt_w(t) )复现时要留意EI论文常把目标函数写成带时间权重的求和形式你拿到的数据、价格要是以小时为单位而论文用的是15分钟一个时段换算错了整个结果曲线都不对。2.2 水电调度约束组水电建模里最核心、也最容易写错的就是水量平衡约束。单库单机组的形式V(t1) V(t) ( Q_in(t) - Q_release(t) ) * dt - S(t)V是库容单位立方米Q_in是天然入库流量单位立方米每秒Q_release是出库流量含发电流量S是弃水流量。乘上dt换算成时段水量。这里出库流量和出力之间还有一层转换关系P_h(t) K * Q_turbine(t) * H_eff(t)K是综合出力系数H_eff是净水头。问题在于H_eff本身又跟库容V相关。很多初学者第一次写约束直接把P写成Q的线性函数拿去跑模型。对于有些论文这种简化是允许的但如果审稿人较真水头变化就必须当作变量处理这时候就得引入库容-水位曲线把水头变成一个随库容变化的非线性函数模型立刻从LP变成NLP/MILP。除了水量平衡还有一组常用约束我列全给你库容上下限V_min ≤ V(t) ≤ V_max防止出现调度方案让水库干涸或漫坝。出库流量限值Q_min ≤ Q_release(t) ≤ Q_max兼顾下游生态和泄洪能力。水电出力限值P_h_min ≤ P_h(t) ≤ P_h_max受机组技术出力限制。爬坡约束| P_h(t) - P_h(t-1) | ≤ ramp_h水电机组出力调整速率不是无限的。末库容约束V(T) V_target这是调度周期结束时的控制目标。梯级水电站还要加“水流滞时”约束即上游电站的出库要经过一个时间延迟才能进入下游电站的入库Q_in_down(t) Q_release_up(t - tau) Q_local(t)这个滞时一加约束矩阵就不是简单对角结构了写Matlab时要注意时序索引别越界。2.3 风电不确定性的两种处理套路风电出力最简单的表达是乘以装机容量和预测出力系数P_w(t) ≤ cap_w * u_w(t)再叠加预测误差。场景法的处理步骤是对风速预测值加随机扰动生成N条风电出力样本常用拉丁超立方采样比普通的蒙特卡洛效率高得多。用场景削减算法把N条聚合成S条典型场景快进快出的话可以直接用Matlab自带的kmeans但要跟论文对上的话还是推荐写同步回代削减算法。目标函数变成各场景期望F_total sum( prob_s * F_s )决策变量要满足所有场景下的约束也就是常说的时间耦合约束。鲁棒法的处理则不同它给风电设一个波动范围把约束改写成P_w(t) ∈ [ P_w_pred(t) - delta_down(t), P_w_pred(t) delta_up(t) ]模型求解时寻找对任意区间内风电出力都可行的调度方案。这种模型在Matlab里常用“对偶重构”或“列与约束生成”来做代码量比场景法大但结果更保守、更贴近电网调度需求。3. Matlab复现的关键设计与求解器选型3.1 我为什么推荐“建模层求解层”分离很多人第一次做优化调度喜欢把目标函数和约束全部写进一个巨长的主程序里然后调fmincon一把梭。代码短的时候还行一旦约束数量上百、变量上千排查维度不匹配错误会让人怀疑人生。我自己复现这类EI论文用的是分层结构数据层所有物理参数、预测曲线、价格曲线集中在一个脚本里运行一次生成结构体保存。模型层目标函数、约束条件按逻辑分模块编写每个模块就是一个函数输入决策变量、输出目标值或约束残差。求解层统一调用求解器并做好求解状态校验solver status、退出标志、残差分析。后处理层画曲线、计算指标、输出表格。这样做的最大好处是论文改了目标函数里的一个惩罚系数你只需要改模型层一个函数求解器从linprog换成Gurobi也不用动其他代码。我实测下来这个结构在复现“多场景随机优化”时尤其关键因为要跑的算例多可复用性直接决定你一个晚上能出几个实验。3.2 自带linprog与外部求解器的取舍Matlab自带优化工具箱里最常用的是linprog线性规划纯LP问题直接用它。intlinprog混合整数线性规划处理机组启停、0-1状态变量时用它。fmincon非线性约束优化水头变化、非线性运行成本等场景用它。先说我的结论能线性化就线性化优先用linprog/intlinprog没有把握就先用fmincon把模型跑通再考虑升级求解器。为什么这么建议因为linprog求解速度快、结果稳定、几乎不会出现数值病态的问题而且Matlab自带的求解器不需要额外安装工具复现门槛低。fmincon对初值敏感初值给不好结果可能落在局部最优。现在很多EI论文里的风-水电联合调度已经写成MILP形式我自己拿到代码的第一步就是检查它的约束里有没有整数变量有就上intlinprog。如果你有商业求解器比如Gurobi、CPLEX建议配合YALMIP一起用。YALMIP是一个Matlab建模工具箱它把LP、MILP、SDP统一封装你只需要写符号化的约束模型层代码能缩短一半。不过用YALMIP要注意版本匹配和路径配置问题见后文的坑。3.3 代码模块设计与数据结构我通常把24小时或96时段的调度问题组织成这样的数据结构%% 数据层示例params.m params.T 24; % 时段数1h步长 params.dt 3600; % 时段长度单位秒 params.rho_w 3.5; % 风电上网电价元/kWh params.rho_h 2.8; % 水电上网电价元/kWh params.lambda_curt 8.0; % 弃风惩罚系数元/kWh params.P_w_forecast ...; % 1xT 风电预测出力MW params.cap_w 300; % 风电场装机MW params.V_init 1.2e8; % 水库初始库容m3 params.V_max 2.0e8; params.V_min 5.0e7; params.V_target 1.2e8; % 末库容目标 params.Q_in ...; % 入库流量序列m3/s params.P_h_max 120; % 水电最大出力MW params.P_h_min 10; params.K 8.5; % 出力系数约 MW/(m3/s * m) params.H_ref 30; % 参考水头m把参数集中放的结构后面调试时改起来方便。再强调一下单位水头单位米流量单位立方米每秒出力单位兆瓦电价按千瓦时算这三个量纲之间的系数如果不在一个量级会导致求解器数值上出现巨大差异后文专门说。4. 核心代码模块实现与参数配置实操4.1 基础数据准备含参数表我自己复现时最喜欢先把论文里的系统参数表格化再对照填入Matlab。这里给你一张可以直接抄作业的参数表模板参数符号数值单位调度时段数T24h风电装机cap_w300MW风电预测出力均值Pw_avg120MW风电预测误差比例sigma_w0.15%水库初始库容V_init1.2亿m3水库最大库容V_max2.0亿m3水库最小库容V_min0.5亿m3天然入库流量Q_in40随季节变化m3/s水电最大出力Ph_max120MW水电最小出力Ph_min10MW综合出力系数K8.5-参考水头H_ref30m弃风惩罚系数lambda_curt300元/MWh风电电价rho_w350元/MWh水电电价rho_h280元/MWh这里有个细节弃风惩罚系数如果设得太低模型会倾向弃风而不调用水电调节能力结果看着“最优”实际完全违背联合调度的初衷。我一般先把惩罚系数设成风电电价的2~3倍确保解出的方案优先消纳风电。4.2 目标函数与约束的矩阵化写法假定你已经把问题简化为LP水头恒定、不引入整数变量那么目标函数可以写成标准形min f * x其中决策变量x包含水电出力序列 P_h(T)、风电实际出力 P_w(T)、弃风量 P_curt(T)。%% 模型层示例build_model.m % 决策变量排序x [P_h(1..T), P_w(1..T), P_curt(1..T)] n_h T; n_w T; n_c T; xlen n_h n_w n_c; % 目标最小化运行成本 % 简化水电成本忽略目标主要是弃风惩罚 f zeros(xlen,1); f(n_hn_w1 : end) params.lambda_curt .* ones(T,1);功率平衡约束写成等号约束P_h(t) P_w(t) load(t)其中负荷场景按论文设置。线性约束统一放到A*x b和Aeq*x beq里% 功率平衡约束 Aeq Aeq zeros(T, xlen); beq load_profile(:); % T x 1 for t 1:T Aeq(t, t) 1; % P_h(t) Aeq(t, n_ht) 1; % P_w(t) % P_curt不参与平衡 end % 出力上下限 lb ub lb [params.P_h_min * ones(T,1); zeros(T,1); zeros(T,1)]; ub [params.P_h_max * ones(T,1); params.cap_w * ones(T,1); params.cap_w * ones(T,1)];这里有个很隐蔽的问题风电实际出力 P_w 被纳入决策变量但没有限制它与预测值的关系。要加一条约束实际出力不超过预测值即 P_w(t) ≤ P_w_forecast(t)。如果不加模型可能会“凭空”让风电多出力来满足负荷平衡优化结果完全失真这是新人最容易犯的错误之一。水库水量平衡约束如果只靠矩阵拼你会觉得写起来很痛苦所以我推荐一个小技巧用稀疏矩阵构造时间耦合项。% 水量平衡V(t1) - V(t) - Q_rel(t)*dt Q_in(t)*dt % V是辅助变量时会增加变量维度简化版本直接把V消去只保留Q_rel约束 % 上面简化模型已经假设水头恒定,所以 P_h(t) K*H_ref*Q_turb(t) % 反解得发电流量Q_turb P_h / (K*H_ref)单位换算后代入平衡约束其实在做纯LP快速复现时我不建议直接引入库容变量而是把水量平衡约束转换为“水电日发电水量总数限制”sum(P_h(t)) / (K * H_ref) * dt some_total_energy这样就将水量平衡变成了对水电全天发电量的整体约束相当于约束发电用水总量等于可用水量适合调度周期内水位变化不大的场景。当然这只是一种近似完整论文里你要保留逐时段库容约束那就得把V(t)也作为决策变量加进去并配合V_MIN、V_MAX上下限。具体哪种合适取决于原论文有没有考虑库容动态变化。稳妥起见我强烈建议第一次复现先跑通整体水量约束出结果后再升级到逐时段库容模型。4.3 主程序循环与结果输出主程序负责把上述模块串起来并做求解状态判断。我的主程序往往这样写%% 主程序run_scheduling.m params data_init(); % 构造模型矩阵 [f, A, b, Aeq, beq, lb, ub] build_model(params); % 求解 options optimoptions(linprog, Algorithm, dual-simplex, Display, iter); [x_opt, fval, exitflag] linprog(f, A, b, Aeq, beq, lb, ub, options); if exitflag ~ 1 error(求解失败exitflag %d, exitflag); end % 提取结果 P_h_opt x_opt(1:T); P_w_opt x_opt(T1:2*T); P_curt_opt x_opt(2*T1:3*T); % 后处理 plot_schedule(P_h_opt, P_w_opt, P_curt_opt); compute_metrics(P_h_opt, P_w_opt, P_curt_opt);如果报错“Linprog stopped because it exceeded its iteration limit”加一行options optimoptions(linprog, MaxIterations, 10000, OptimalityTolerance, 1e-8);求解器默认的迭代限制是针对标准规模问题的你加了288个时段、几千个变量之后迭代上限经常不够用。5. 结果解读与验证方法5.1 先跑一个确定性基线拿到代码第一步不是马上加随机性而是先把风电预测设为一条确定性曲线跑一遍纯确定性调度。这一步的目的是验证模型本身对不对。验证点有三个功率平衡是否全时段满足输出每个时段的误差应该小于1e-6。水电出力是否在可行域内不会出现负值或者越过装机上限。水库水量是否守恒把各时段的发电流量对时间积分对比可用水量误差应该在可接受范围内。我通常会画一张调度图横轴是24小时纵轴是功率把水电出力、风电出力、负荷和弃风量画在同一张图里。如果看到弃风集中在风电大发而负荷低谷时段说明模型逻辑正确如果弃风遍地都是那大概率是约束写错了比如忘了加“实际出力不超预测”那条限制。5.2 随机场景下的对比指标确定性模型通过后再上随机性。场景法的结果处理有固定的套路生成500~1000条风电场景用同步回代削减到10~20个代表场景。求解多场景随机优化模型得到第一阶段的调度计划。用这个计划去回代测试所有原始场景统计平均弃风率、失负荷率等指标。注意这里有个“阶段”的概念。如果模型假设所有决策都在看到风电出力之前就定好属于完全前瞻性的调度过于乐观如果允许水电在每个时段根据风电实时出力滚动调整则结果更贴近实际。论文里叫“非预期性约束”复现时务必看清原文属于哪一种。计算指标时我建议输出四个核心指标指标公式意义弃风率弃风量 / 风电可发总量衡量风电消纳水平水电利用率水电发电量 / 水电可发电量衡量水量利用程度运行成本目标函数值经济性爬坡次数统计水电出力方向变化的次数衡量调节压力这四个指标放在对比表格里基本能满足课程作业、学位论文里“算例分析”一节的需要也能支撑你写完一篇数据扎实的复现笔记。6. 复现中高频踩坑与排查实录6.1 求解器数值异常先说我最常遇到的linprog报“The dual-simplex algorithm exhausted its options”或者“Unable to find a feasible solution”。这种问题的原因比例大概是约束本身相矛盾占30%数值尺度差异过大占50%模型数据错误占20%。数值尺度差异是最隐蔽的。库容是亿立方米电价是几百元每兆瓦时出力是几十到几百兆瓦三个量级放在一起约束矩阵的条件数会非常差。解法很简单把所有物理量折算到统一量纲。比如库容改用百万立方米流量改用百万立方米每小时电价改成元每兆瓦时你会发现收敛速度提升非常明显。再者约束矛盾常见于水量平衡和末库容约束之间。比如天然来水很少但你又要求调度期末库容回到一个很高的目标值同时还要保证水电全天大出力这组约束在物理上就是无解的。排查方法先注释末库容约束跑一遍看水电出力上限是否被顶格再决定调整来水数据还是放宽末库容。fmincon场景下最常见报错是“Objective function returned NaN”。十有八九是决策变量超出可行域导致计算了负数的平方根或者除零。我会在目标函数最前面加一行if any(x 0) fval 1e10; return; end用罚函数快速把搜索拉回可行域。6.2 约束维度与单位换算另一个高频坑是矩阵维度对不上。linprog的约束矩阵要求是A * x b其中A的行数等于约束条数列数等于变量数。一旦你中途改了时段数、加了变量忘记同步扩容A就会出现“Matrix dimensions must agree”或者“Linprog received the incorrect number and/or types of parameters.”。我的习惯是每次修改模型之后先跑一段自动检查assert(size(A,2) xlen, A列数%d不等于变量数%d, size(A,2), xlen); assert(size(Aeq,2) xlen, Aeq列数与变量数不一致); assert(length(lb) xlen length(ub) xlen, 边界维度错误);这段代码能省下大量低级排错时间。单位换算的问题前面反复提过这里再给一个具体例子。入库流量40 m³/s如果直接乘dt3600得到的是144000 m³这个数跟库容一亿立方米放在一起做水量平衡精度上没问题。但如果你用linprog默认的容差很容易把这个小量直接忽略掉导致水量平衡约束形同虚设。解决办法是给水量平衡约束单独配置容差或者在建模时把水量单位从立方米放大到百万立方米。6.3 YALMIP与求解器配置问题如果你选择用YALMIP加外部求解器最常见的坑是YALMIP安装后找不到求解器。原因大多是没有把安装目录添加到Matlab路径或者是Gurobi的许可证没有正确设置。启动时先跑yalmiptest它会自动检测可用求解器看到输出的汇总表格里出现了你要用的求解器就说明配置成功。如果提示“No suitable solver for this problem type”则是你建的模型里有变量类型或约束类型当前求解器不支持。排查思路先确认问题是LP还是MILP再确认是否出现了非线性项。很多时候YALMIP会把约束判断为非线性实际只是你在写约束时用了sdpvar变量的二次项比如写成了x(1)^2其实应该是线性化后的变量。6.4 结果合理性检查最后分享几个我复现完以后必做的合理性检查项也是写论文时审稿人最爱问的点全时段功率平衡误差小于1e-5否则说明约束漏写。任何弃风时段前面应该堆积了水库蓄水或者水电出力下限顶格的情况弃风不会无缘无故出现。末库容约束是否满足如果水库蓄水大量用光但未来几天没有来水数据调度结果就很危险。水电爬坡出力要符合机组实际调节速率不要出现15分钟内从10 MW跳到120 MW的“神仙调度”。我个人体会是这条检查清单比我第一次看论文时以为的更重要。很多论文复现失败往往不是算法不对而是这些基础物理约束没满足导致结果曲线看上去“很漂亮、但没法解释”。写在最后的实操建议如果只能给你一条建议那就是先把只含水电的单库优化跑通再叠加风电最后才处理不确定性。我见过太多人一上来就直奔随机规划结果模型线性化、场景削减、非预期性约束叠在一起出问题了根本不知道是哪一层的锅。以我自己的经验一篇典型的EI论文风电-水电联合调度复现正常节奏是第一天搭数据层加确定性模型第二天加场景生成和削减第三天补后处理和算例对比。三天能完成大部分工作剩下时间都花在排查约束和调参上。等你这套代码跑顺之后后面扩展空间很大。把水电换成抽水蓄能模型结构几乎不用大改把目标函数从发电收益换成系统备用容量最大化也只是多几个约束的事。Matlab这套环境的优势就在这儿建模灵活代码复用率高换个场景就是一篇新的算例分析。
做风-水电联合优化的人应该都有同感单独看风电随机性、间歇性大到让人头疼单独看水电又受来水、水库调节能力和生态流量约束。但把两者放在同一个调度框架里做联合优化运行分析往往能在几乎不新增硬件投资的情况下把系统的灵活性和可再生能源消纳空间拉高一大截。用Matlab来做这块的代码复现是我个人觉得性价比最高的路径。这篇文章就围绕“风-水电联合优化运行分析”这个标题把EI论文里常出现的建模思路、目标函数与约束条件怎么落地成Matlab代码、求解器怎么选、结果怎么验证、以及复现过程中最容易踩的坑全部过一遍。适合正在复现论文、做课程设计、或者准备写小论文的电力方向研究生也适合想快速上手调度优化模型的工程师参考。1. 先说清楚风-水电联合优化运行到底在解决什么问题1.1 为什么风电场要跟水电站“绑”在一起风电出力靠天吃饭风速曲线长什么样出力就长什么样基本没有太多人工干预的空间。水电则不同只要水库有足够调节库容机组出力可以在一定范围内快速调节。把这两个电源放在一个调度模型里核心思路是用可控性好、响应速度快的水电去平抑风电的波动和预测误差。风电大发的时候水库少放水、多蓄水风电小发的时候水库多放水、补足缺口。这样一来电网调度看到的是一条更平稳的外送功率曲线原先不得不弃掉的风电也能被吸收进系统。这在术语里叫“互补运行”但你要真去复现论文就会发现难点不在概念而在模型。水文过程有水量平衡风电有预测误差两者耦合在一起之后约束条件呈指数级膨胀调度时段一旦取到96点甚至288点普通Excel根本算不动必须依靠Matlab这样的数值计算环境。1.2 EI论文里最常见的三类模型我在复现不同文献时发现风-水电联合优化模型的差别主要在三处时间尺度、目标函数、不确定性处理方式。给你一个速查表模型类型典型时间尺度目标函数常见求解方法长期调度月/旬发电量最大、蓄能最大DP、逐次逼近短期调度日/96点运行成本最小、弃风最小MILP、NLP实时/滚动15min级跟踪计划偏差最小MPC、鲁棒优化真正的EI复现绝大多数集中在“短期调度”这一档。再往上文献会在这个框架里叠加风电不确定性产生两类变体场景法用蒙特卡洛抽样生成大量风电出力场景然后用场景削减技术浓缩成几个典型场景模型变成随机规划。鲁棒法不关心具体概率分布只给风电出力一个波动区间模型追求最坏情况下的调度可行性。这两条路线在Matlab里的实现路径差别很大后文我会分别讲怎么处理。2. 建模环节从物理对象到数学表达2.1 目标函数怎么选选目标函数首要判断依据是你复现的论文里系统调度主体是谁追求的经济指标是什么。如果是独立发电企业视角目标函数通常是“系统运行收益最大”即水电上网收益加上风电上网收益减去运行维护成本、弃风惩罚成本。公式上写作max F sum( sum( P_h(i,t) * price_h(t) * dt ) sum( P_w(i,t) * price_w(t) * dt ) ) - sum( curtail_w(i,t) * penalty * dt )如果是电网调度视角目标函数更多写成“系统运行成本最小”这时水电出力成本用可变成本函数近似风电成本近似为0或者说优先消纳目标则变成减少火电出力和弃风惩罚min C sum( sum( a_i * P_g(i,t)^2 b_i * P_g(i,t) c_i ) ) sum( lambda_cur * curt_w(t) )复现时要留意EI论文常把目标函数写成带时间权重的求和形式你拿到的数据、价格要是以小时为单位而论文用的是15分钟一个时段换算错了整个结果曲线都不对。2.2 水电调度约束组水电建模里最核心、也最容易写错的就是水量平衡约束。单库单机组的形式V(t1) V(t) ( Q_in(t) - Q_release(t) ) * dt - S(t)V是库容单位立方米Q_in是天然入库流量单位立方米每秒Q_release是出库流量含发电流量S是弃水流量。乘上dt换算成时段水量。这里出库流量和出力之间还有一层转换关系P_h(t) K * Q_turbine(t) * H_eff(t)K是综合出力系数H_eff是净水头。问题在于H_eff本身又跟库容V相关。很多初学者第一次写约束直接把P写成Q的线性函数拿去跑模型。对于有些论文这种简化是允许的但如果审稿人较真水头变化就必须当作变量处理这时候就得引入库容-水位曲线把水头变成一个随库容变化的非线性函数模型立刻从LP变成NLP/MILP。除了水量平衡还有一组常用约束我列全给你库容上下限V_min ≤ V(t) ≤ V_max防止出现调度方案让水库干涸或漫坝。出库流量限值Q_min ≤ Q_release(t) ≤ Q_max兼顾下游生态和泄洪能力。水电出力限值P_h_min ≤ P_h(t) ≤ P_h_max受机组技术出力限制。爬坡约束| P_h(t) - P_h(t-1) | ≤ ramp_h水电机组出力调整速率不是无限的。末库容约束V(T) V_target这是调度周期结束时的控制目标。梯级水电站还要加“水流滞时”约束即上游电站的出库要经过一个时间延迟才能进入下游电站的入库Q_in_down(t) Q_release_up(t - tau) Q_local(t)这个滞时一加约束矩阵就不是简单对角结构了写Matlab时要注意时序索引别越界。2.3 风电不确定性的两种处理套路风电出力最简单的表达是乘以装机容量和预测出力系数P_w(t) ≤ cap_w * u_w(t)再叠加预测误差。场景法的处理步骤是对风速预测值加随机扰动生成N条风电出力样本常用拉丁超立方采样比普通的蒙特卡洛效率高得多。用场景削减算法把N条聚合成S条典型场景快进快出的话可以直接用Matlab自带的kmeans但要跟论文对上的话还是推荐写同步回代削减算法。目标函数变成各场景期望F_total sum( prob_s * F_s )决策变量要满足所有场景下的约束也就是常说的时间耦合约束。鲁棒法的处理则不同它给风电设一个波动范围把约束改写成P_w(t) ∈ [ P_w_pred(t) - delta_down(t), P_w_pred(t) delta_up(t) ]模型求解时寻找对任意区间内风电出力都可行的调度方案。这种模型在Matlab里常用“对偶重构”或“列与约束生成”来做代码量比场景法大但结果更保守、更贴近电网调度需求。3. Matlab复现的关键设计与求解器选型3.1 我为什么推荐“建模层求解层”分离很多人第一次做优化调度喜欢把目标函数和约束全部写进一个巨长的主程序里然后调fmincon一把梭。代码短的时候还行一旦约束数量上百、变量上千排查维度不匹配错误会让人怀疑人生。我自己复现这类EI论文用的是分层结构数据层所有物理参数、预测曲线、价格曲线集中在一个脚本里运行一次生成结构体保存。模型层目标函数、约束条件按逻辑分模块编写每个模块就是一个函数输入决策变量、输出目标值或约束残差。求解层统一调用求解器并做好求解状态校验solver status、退出标志、残差分析。后处理层画曲线、计算指标、输出表格。这样做的最大好处是论文改了目标函数里的一个惩罚系数你只需要改模型层一个函数求解器从linprog换成Gurobi也不用动其他代码。我实测下来这个结构在复现“多场景随机优化”时尤其关键因为要跑的算例多可复用性直接决定你一个晚上能出几个实验。3.2 自带linprog与外部求解器的取舍Matlab自带优化工具箱里最常用的是linprog线性规划纯LP问题直接用它。intlinprog混合整数线性规划处理机组启停、0-1状态变量时用它。fmincon非线性约束优化水头变化、非线性运行成本等场景用它。先说我的结论能线性化就线性化优先用linprog/intlinprog没有把握就先用fmincon把模型跑通再考虑升级求解器。为什么这么建议因为linprog求解速度快、结果稳定、几乎不会出现数值病态的问题而且Matlab自带的求解器不需要额外安装工具复现门槛低。fmincon对初值敏感初值给不好结果可能落在局部最优。现在很多EI论文里的风-水电联合调度已经写成MILP形式我自己拿到代码的第一步就是检查它的约束里有没有整数变量有就上intlinprog。如果你有商业求解器比如Gurobi、CPLEX建议配合YALMIP一起用。YALMIP是一个Matlab建模工具箱它把LP、MILP、SDP统一封装你只需要写符号化的约束模型层代码能缩短一半。不过用YALMIP要注意版本匹配和路径配置问题见后文的坑。3.3 代码模块设计与数据结构我通常把24小时或96时段的调度问题组织成这样的数据结构%% 数据层示例params.m params.T 24; % 时段数1h步长 params.dt 3600; % 时段长度单位秒 params.rho_w 3.5; % 风电上网电价元/kWh params.rho_h 2.8; % 水电上网电价元/kWh params.lambda_curt 8.0; % 弃风惩罚系数元/kWh params.P_w_forecast ...; % 1xT 风电预测出力MW params.cap_w 300; % 风电场装机MW params.V_init 1.2e8; % 水库初始库容m3 params.V_max 2.0e8; params.V_min 5.0e7; params.V_target 1.2e8; % 末库容目标 params.Q_in ...; % 入库流量序列m3/s params.P_h_max 120; % 水电最大出力MW params.P_h_min 10; params.K 8.5; % 出力系数约 MW/(m3/s * m) params.H_ref 30; % 参考水头m把参数集中放的结构后面调试时改起来方便。再强调一下单位水头单位米流量单位立方米每秒出力单位兆瓦电价按千瓦时算这三个量纲之间的系数如果不在一个量级会导致求解器数值上出现巨大差异后文专门说。4. 核心代码模块实现与参数配置实操4.1 基础数据准备含参数表我自己复现时最喜欢先把论文里的系统参数表格化再对照填入Matlab。这里给你一张可以直接抄作业的参数表模板参数符号数值单位调度时段数T24h风电装机cap_w300MW风电预测出力均值Pw_avg120MW风电预测误差比例sigma_w0.15%水库初始库容V_init1.2亿m3水库最大库容V_max2.0亿m3水库最小库容V_min0.5亿m3天然入库流量Q_in40随季节变化m3/s水电最大出力Ph_max120MW水电最小出力Ph_min10MW综合出力系数K8.5-参考水头H_ref30m弃风惩罚系数lambda_curt300元/MWh风电电价rho_w350元/MWh水电电价rho_h280元/MWh这里有个细节弃风惩罚系数如果设得太低模型会倾向弃风而不调用水电调节能力结果看着“最优”实际完全违背联合调度的初衷。我一般先把惩罚系数设成风电电价的2~3倍确保解出的方案优先消纳风电。4.2 目标函数与约束的矩阵化写法假定你已经把问题简化为LP水头恒定、不引入整数变量那么目标函数可以写成标准形min f * x其中决策变量x包含水电出力序列 P_h(T)、风电实际出力 P_w(T)、弃风量 P_curt(T)。%% 模型层示例build_model.m % 决策变量排序x [P_h(1..T), P_w(1..T), P_curt(1..T)] n_h T; n_w T; n_c T; xlen n_h n_w n_c; % 目标最小化运行成本 % 简化水电成本忽略目标主要是弃风惩罚 f zeros(xlen,1); f(n_hn_w1 : end) params.lambda_curt .* ones(T,1);功率平衡约束写成等号约束P_h(t) P_w(t) load(t)其中负荷场景按论文设置。线性约束统一放到A*x b和Aeq*x beq里% 功率平衡约束 Aeq Aeq zeros(T, xlen); beq load_profile(:); % T x 1 for t 1:T Aeq(t, t) 1; % P_h(t) Aeq(t, n_ht) 1; % P_w(t) % P_curt不参与平衡 end % 出力上下限 lb ub lb [params.P_h_min * ones(T,1); zeros(T,1); zeros(T,1)]; ub [params.P_h_max * ones(T,1); params.cap_w * ones(T,1); params.cap_w * ones(T,1)];这里有个很隐蔽的问题风电实际出力 P_w 被纳入决策变量但没有限制它与预测值的关系。要加一条约束实际出力不超过预测值即 P_w(t) ≤ P_w_forecast(t)。如果不加模型可能会“凭空”让风电多出力来满足负荷平衡优化结果完全失真这是新人最容易犯的错误之一。水库水量平衡约束如果只靠矩阵拼你会觉得写起来很痛苦所以我推荐一个小技巧用稀疏矩阵构造时间耦合项。% 水量平衡V(t1) - V(t) - Q_rel(t)*dt Q_in(t)*dt % V是辅助变量时会增加变量维度简化版本直接把V消去只保留Q_rel约束 % 上面简化模型已经假设水头恒定,所以 P_h(t) K*H_ref*Q_turb(t) % 反解得发电流量Q_turb P_h / (K*H_ref)单位换算后代入平衡约束其实在做纯LP快速复现时我不建议直接引入库容变量而是把水量平衡约束转换为“水电日发电水量总数限制”sum(P_h(t)) / (K * H_ref) * dt some_total_energy这样就将水量平衡变成了对水电全天发电量的整体约束相当于约束发电用水总量等于可用水量适合调度周期内水位变化不大的场景。当然这只是一种近似完整论文里你要保留逐时段库容约束那就得把V(t)也作为决策变量加进去并配合V_MIN、V_MAX上下限。具体哪种合适取决于原论文有没有考虑库容动态变化。稳妥起见我强烈建议第一次复现先跑通整体水量约束出结果后再升级到逐时段库容模型。4.3 主程序循环与结果输出主程序负责把上述模块串起来并做求解状态判断。我的主程序往往这样写%% 主程序run_scheduling.m params data_init(); % 构造模型矩阵 [f, A, b, Aeq, beq, lb, ub] build_model(params); % 求解 options optimoptions(linprog, Algorithm, dual-simplex, Display, iter); [x_opt, fval, exitflag] linprog(f, A, b, Aeq, beq, lb, ub, options); if exitflag ~ 1 error(求解失败exitflag %d, exitflag); end % 提取结果 P_h_opt x_opt(1:T); P_w_opt x_opt(T1:2*T); P_curt_opt x_opt(2*T1:3*T); % 后处理 plot_schedule(P_h_opt, P_w_opt, P_curt_opt); compute_metrics(P_h_opt, P_w_opt, P_curt_opt);如果报错“Linprog stopped because it exceeded its iteration limit”加一行options optimoptions(linprog, MaxIterations, 10000, OptimalityTolerance, 1e-8);求解器默认的迭代限制是针对标准规模问题的你加了288个时段、几千个变量之后迭代上限经常不够用。5. 结果解读与验证方法5.1 先跑一个确定性基线拿到代码第一步不是马上加随机性而是先把风电预测设为一条确定性曲线跑一遍纯确定性调度。这一步的目的是验证模型本身对不对。验证点有三个功率平衡是否全时段满足输出每个时段的误差应该小于1e-6。水电出力是否在可行域内不会出现负值或者越过装机上限。水库水量是否守恒把各时段的发电流量对时间积分对比可用水量误差应该在可接受范围内。我通常会画一张调度图横轴是24小时纵轴是功率把水电出力、风电出力、负荷和弃风量画在同一张图里。如果看到弃风集中在风电大发而负荷低谷时段说明模型逻辑正确如果弃风遍地都是那大概率是约束写错了比如忘了加“实际出力不超预测”那条限制。5.2 随机场景下的对比指标确定性模型通过后再上随机性。场景法的结果处理有固定的套路生成500~1000条风电场景用同步回代削减到10~20个代表场景。求解多场景随机优化模型得到第一阶段的调度计划。用这个计划去回代测试所有原始场景统计平均弃风率、失负荷率等指标。注意这里有个“阶段”的概念。如果模型假设所有决策都在看到风电出力之前就定好属于完全前瞻性的调度过于乐观如果允许水电在每个时段根据风电实时出力滚动调整则结果更贴近实际。论文里叫“非预期性约束”复现时务必看清原文属于哪一种。计算指标时我建议输出四个核心指标指标公式意义弃风率弃风量 / 风电可发总量衡量风电消纳水平水电利用率水电发电量 / 水电可发电量衡量水量利用程度运行成本目标函数值经济性爬坡次数统计水电出力方向变化的次数衡量调节压力这四个指标放在对比表格里基本能满足课程作业、学位论文里“算例分析”一节的需要也能支撑你写完一篇数据扎实的复现笔记。6. 复现中高频踩坑与排查实录6.1 求解器数值异常先说我最常遇到的linprog报“The dual-simplex algorithm exhausted its options”或者“Unable to find a feasible solution”。这种问题的原因比例大概是约束本身相矛盾占30%数值尺度差异过大占50%模型数据错误占20%。数值尺度差异是最隐蔽的。库容是亿立方米电价是几百元每兆瓦时出力是几十到几百兆瓦三个量级放在一起约束矩阵的条件数会非常差。解法很简单把所有物理量折算到统一量纲。比如库容改用百万立方米流量改用百万立方米每小时电价改成元每兆瓦时你会发现收敛速度提升非常明显。再者约束矛盾常见于水量平衡和末库容约束之间。比如天然来水很少但你又要求调度期末库容回到一个很高的目标值同时还要保证水电全天大出力这组约束在物理上就是无解的。排查方法先注释末库容约束跑一遍看水电出力上限是否被顶格再决定调整来水数据还是放宽末库容。fmincon场景下最常见报错是“Objective function returned NaN”。十有八九是决策变量超出可行域导致计算了负数的平方根或者除零。我会在目标函数最前面加一行if any(x 0) fval 1e10; return; end用罚函数快速把搜索拉回可行域。6.2 约束维度与单位换算另一个高频坑是矩阵维度对不上。linprog的约束矩阵要求是A * x b其中A的行数等于约束条数列数等于变量数。一旦你中途改了时段数、加了变量忘记同步扩容A就会出现“Matrix dimensions must agree”或者“Linprog received the incorrect number and/or types of parameters.”。我的习惯是每次修改模型之后先跑一段自动检查assert(size(A,2) xlen, A列数%d不等于变量数%d, size(A,2), xlen); assert(size(Aeq,2) xlen, Aeq列数与变量数不一致); assert(length(lb) xlen length(ub) xlen, 边界维度错误);这段代码能省下大量低级排错时间。单位换算的问题前面反复提过这里再给一个具体例子。入库流量40 m³/s如果直接乘dt3600得到的是144000 m³这个数跟库容一亿立方米放在一起做水量平衡精度上没问题。但如果你用linprog默认的容差很容易把这个小量直接忽略掉导致水量平衡约束形同虚设。解决办法是给水量平衡约束单独配置容差或者在建模时把水量单位从立方米放大到百万立方米。6.3 YALMIP与求解器配置问题如果你选择用YALMIP加外部求解器最常见的坑是YALMIP安装后找不到求解器。原因大多是没有把安装目录添加到Matlab路径或者是Gurobi的许可证没有正确设置。启动时先跑yalmiptest它会自动检测可用求解器看到输出的汇总表格里出现了你要用的求解器就说明配置成功。如果提示“No suitable solver for this problem type”则是你建的模型里有变量类型或约束类型当前求解器不支持。排查思路先确认问题是LP还是MILP再确认是否出现了非线性项。很多时候YALMIP会把约束判断为非线性实际只是你在写约束时用了sdpvar变量的二次项比如写成了x(1)^2其实应该是线性化后的变量。6.4 结果合理性检查最后分享几个我复现完以后必做的合理性检查项也是写论文时审稿人最爱问的点全时段功率平衡误差小于1e-5否则说明约束漏写。任何弃风时段前面应该堆积了水库蓄水或者水电出力下限顶格的情况弃风不会无缘无故出现。末库容约束是否满足如果水库蓄水大量用光但未来几天没有来水数据调度结果就很危险。水电爬坡出力要符合机组实际调节速率不要出现15分钟内从10 MW跳到120 MW的“神仙调度”。我个人体会是这条检查清单比我第一次看论文时以为的更重要。很多论文复现失败往往不是算法不对而是这些基础物理约束没满足导致结果曲线看上去“很漂亮、但没法解释”。写在最后的实操建议如果只能给你一条建议那就是先把只含水电的单库优化跑通再叠加风电最后才处理不确定性。我见过太多人一上来就直奔随机规划结果模型线性化、场景削减、非预期性约束叠在一起出问题了根本不知道是哪一层的锅。以我自己的经验一篇典型的EI论文风电-水电联合调度复现正常节奏是第一天搭数据层加确定性模型第二天加场景生成和削减第三天补后处理和算例对比。三天能完成大部分工作剩下时间都花在排查约束和调参上。等你这套代码跑顺之后后面扩展空间很大。把水电换成抽水蓄能模型结构几乎不用大改把目标函数从发电收益换成系统备用容量最大化也只是多几个约束的事。Matlab这套环境的优势就在这儿建模灵活代码复用率高换个场景就是一篇新的算例分析。
