做水电调度优化的课题通常都有同一个痛点仿真模型搭起来了但论文里的结果就是复现不出来。梯级水光互补系统的短期优化调度尤其如此——上游连着几个水电站中间嵌着大片光伏基地调度目标不是单纯的“发满”而是在光伏出力的不确定性之下系统还能可靠外送多少电。这篇文章记录的是一个基于Matlab的EI论文复现过程模型的核心是梯级水光互补系统短期优化调度模型目标函数为最大化可消纳电量期望。代码在Matlab R2023a环境下用Yalmip建模、Gurobi求解完整跑通适合正在做新能源消纳、水电优化调度方向的硕博生也适合刚接触随机优化建模的读者当一份可参考的样例。1. 这个模型到底在优化什么问题拆解与场景设定1.1 梯级水光互补调度的现实约束把物理场景先梳理清楚。所谓梯级电站就是同一流域上下游串了多个水电站上游电站的出库流量经过一段时间延迟变成下游电站的入库流量。水光互补的意思是光伏电站和水电站共用一条外送通道光伏出力波动时由水电来“托底”和“补差”。举个例子比如某流域有A、B、C三个梯级电站总装机分别约为800MW、500MW、600MW光伏装机1200MW。白天中午光伏出力逼近峰值就要求水电压低出力给光伏让出通道傍晚光伏骤减水电必须快速顶上。如果调度只按光伏预测曲线做确定性安排实际出力一旦偏离预测就会造成外送通道闲置或者弃光。更极端的情况云层遮挡导致光伏瞬时出力大幅下降水电调整不及时就会破坏功率平衡。这些现实约束决定了模型不能是简单的确定性优化必须把光伏随机性纳入决策框架。1.2 “可消纳电量期望”的准确含义标题里最容易被误解的词是“可消纳电量期望”。很多人第一反应以为“期望”就是把预测出力曲线对时间积分其实不是。在短期调度中“可消纳电量”指的是在满足外送通道容量、水电站运行约束、光伏实际出力等条件之后系统真正送出去的电量。由于光伏出力存在随机性这个电量是一个随机变量所以目标函数取它在多个场景下的期望值。数学上写成max E[∑ P_net(ω,t)·Δt]其中ω表示光伏出力场景P_net是场景ω下t时刻的外送功率。这样做的现实意义是调度方案不针对某一条光伏曲线优化而是对所有可能场景的平均表现都尽量好避免“预测准了就赚、预测不准就亏”的极端情况。代码实现时期望运算用场景法落地生成N个光伏场景给每个场景一个概率把期望展开为场景概率加权求和的线性表达式再交给求解器。这也是这类EI论文里最主流的建模思路。1.3 短期优化调度的时间尺度与决策变量短期优化调度通常指日前调度时间窗口T24小时步长Δt1小时。如果做日内滚动可能会用15分钟分辨率但原理完全一样。决策变量分两类第一类是水电相关的各电站的出力、发电流量、弃水流量、库容或水位第二类是外送与弃电相关的各时段外送功率、弃光功率。这里有一个非常关键、也是很多复现翻车的点在期望最大化模型里水电出力通常作为第一阶段决策也就是在光伏场景实现之前就要确定下来而外送、弃电是第二阶段决策等看到光伏实际出力后做调整。在随机规划里这就是“here-and-now”和“wait-and-see”的区别。代码层面体现为水电出力变量不标场景索引外送和弃电变量标场景索引。如果把所有变量都带场景下标模型就退化成一堆独立子问题的拼接优化结果没有任何实际调度意义。这一点后面我还会反复强调。2. 数学模型目标函数与约束条件的完整拆解2.1 目标函数如何把期望值写成可求解的形式模型的目标函数可以写成max ∑_{s1}^{N_s} ρ_s · ∑_{t1}^{T} P_deliver(s,t)·Δt其中P_deliver(s,t)是场景s下t时刻实际外送功率ρ_s是场景s的概率N_s是场景总数。为什么要引入P_deliver而不直接用P_hP_pv因为外送功率要同时受制于发电能力和外送通道上限也就是P_deliver(s,t) min(P_h_total(t)P_pv(s,t), P_limit(t))。min函数不可导直接放进目标函数会让求解器无从下手。标准做法是引入辅助变量P_deliver用两个约束限定它P_deliver(s,t) ≤ P_h_total(t) P_pv(s,t) P_deliver(s,t) ≤ P_limit(t) P_deliver(s,t) ≥ 0由于目标函数是最大化P_deliver之和求解器会自动把P_deliver推到两个上限的较小值等价实现了min运算。这组约束连接了第一阶段水电变量和第二阶段场景变量是整个模型“期望可消纳电量最大化”的关键桥梁。注意目标函数里的Δt。如果Δt1小时且P_deliver单位是MW那么∑P_deliver·Δt就是电量MWh。很多人复现时忘了乘Δt导致目标值和论文对不上差了一个小时量级。2.2 梯级水电约束水量平衡与水头耦合梯级水电约束是模型里工程量最大的部分核心是四组约束第一组是水量平衡V_i(t1) V_i(t) [Q_nat_i(t) Q_up_i(t-τ) - Q_turb_i(t) - Q_spill_i(t)]·Δt其中Q_nat是天然来水Q_up是上游电站总出库发电流量加弃水流量τ是水流时滞。如果上游有支流汇入Q_up还要按拓扑结构拆分。第二组是库容上下限V_min_i ≤ V_i(t) ≤ V_max_i第三组是出库流量约束Q_turb_min_i ≤ Q_turb_i(t) ≤ Q_turb_max_i第四组是出力上限P_h_i(t) ≤ 9.81·η_i·Q_turb_i(t)·H_i(t)H_i(t)是净水头等于坝前水位减尾水位再减水头损失。坝前水位是库容的非线性函数净水头又和发电流量有耦合所以水电出力本质上是“发电流量×净水头”的双线性项。在Matlab复现里最常见的有两种处理一是忽略水头变化用恒定水头近似P_h_i(t) ≈ K_i·Q_turb_i(t)模型变成线性规划跑起来非常快二是预先算好“出力-流量-库容”关系表用分段线性约束拟合精度高但会引入0-1变量变成混合整数规划。论文里如果明确说“忽略水头变化”直接线性化就行如果没提建议先用恒定水头跑通再升级成查表法。2.3 光伏出力约束与不确定性场景生成光伏出力在模型里不是决策变量而是已知场景参数。场景生成的常用做法是基于预测曲线加误差扰动P_pv(s,t) max(0, P_pv_forecast(t)·(1 ε(s,t)))其中ε(s,t)是预测误差采样值可以假设服从正态分布或Beta分布。Matlab里的生成代码大致如下N_scen 15; T 24; sigma 0.1; epsilon normrnd(0, sigma, [N_scen, T]); P_pv_scen max(0, repmat(P_pv_forecast, N_scen, 1) .* (1 epsilon)); rho ones(N_scen, 1) / N_scen;如果论文给的场景数据是历史出力曲线那就要先做场景削减K-means聚类或者同步回代削减都行削减完再把概率归一化。场景削减的核心是保“典型性”让剩下场景里同时包含高、中、低出力的代表。2.4 电网消纳通道与弃电惩罚外送通道约束是决定“可消纳”的另一只手P_deliver(s,t) ≤ P_limit(t)P_limit是外送断面容量或者电网需求上限可以是24小时的时序曲线。弃电量定义为P_curtail(s,t) ≥ P_h_total(t) P_pv(s,t) - P_limit(t) P_curtail(s,t) ≥ 0如果论文里区分弃光和弃水就把P_curtail拆成两个变量分别给惩罚权重。通常光伏弃电代价高于弃水因为光伏没有燃料成本而且白天发电窗口错过就没了。还有一个容易漏的约束是功率平衡。并网场景下应该是P_h_total(t) P_pv(s,t) - P_curtail(s,t) P_load(t)也就是系统净负荷功率平衡。如果是孤网或外送场景则替换成外送通道约束。复现时先对照论文把所有等式写全再逐个检查哪些约束是有效的、哪些是松弛的。我见过太多复现失败案例都是漏写了一个不起眼但关键的等式约束。3. Matlab实现细节从公式到可运行代码的三层结构3.1 数据准备参数表与时间序列的组织方式代码实现我先说数据层。强烈建议把所有参数封装成结构体避免脚本里散落一堆魔法数字。下面是我习惯的组织方式% hydro_data.m hydro(1).name A水库; hydro(1).Vmin 1.2e8; % 死库容m^3 hydro(1).Vmax 4.5e8; % 正常蓄水位对应库容m^3 hydro(1).V0 3.0e8; % 初始库容m^3 hydro(1).Qmin 0; hydro(1).Qmax 250; % 最大发电流量m^3/s hydro(1).eta 0.88; % 综合效率系数 hydro(1).maxPower 800; % 装机容量MW hydro(1).K 8.5; % 出力系数恒定水头近似时使用 % 光伏与电网 pv.capacity 1200; % 光伏装机MW pv.forecast load(pv_forecast.mat); % 24小时预测出力MW grid.Plimit 1800; % 外送通道上限MW grid.load load(load_curve.mat); % 净负荷曲线MW这里有个小细节库容单位是m³发电流量单位是m³/s功率单位是MW。约束里把这些不同量纲的数值混在一起数量级可能差10^8求解器数值稳定性会很差。我的处理方式是把所有流量乘以Δt换算成时段内水量单位统一成“百万立方米”功率统一为MW目标函数的Δt用小时。这样各类数值都在几十到几千的量级Gurobi和CPLEX跑起来都舒服很多。3.2 Yalmip建模变量定义、约束装配与求解器调用建模层是代码的核心用Yalmip把第二章的数学模型逐句翻译成约束。核心框架如下% 变量定义 P_h sdpvar(N_hydro, T, full); % 第一阶段水电出力 Q_turb sdpvar(N_hydro, T, full); % 发电流量 Q_spill sdpvar(N_hydro, T, full); % 弃水流量 V sdpvar(N_hydro, T1, full); % 库容T1因为要覆盖首末时段 P_deliver sdpvar(N_scen, T, full); % 第二阶段外送功率 P_curtail sdpvar(N_scen, T, full); % 弃电功率 % 目标函数Yalmip默认求最小化所以取负号 objective -sum(sum(repmat(rho, 1, T) .* P_deliver)); % 约束装配 C []; % 水量平衡约束、库容约束、出力约束... C [C, V(:, t1) V(:, t) (Q_nat Q_up - Q_turb - Q_spill) * dt]; C [C, V(:, t) Vmin, V(:, t) Vmax]; % P_deliver与发电能力、外送上限的关系 C [C, P_deliver(s, t) sum(P_h(:, t)) P_pv_scen(s, t)]; C [C, P_deliver(s, t) Plimit(t)]; C [C, P_deliver(s, t) 0]; % 求解 ops sdpsettings(solver, gurobi, verbose, 2, mipgap, 1e-4); result optimize(C, objective, ops);这里有三个容易踩的坑第一sdpvar的维度方向。sdpvar(N,T,full)生成的是N行T列矩阵而sdpvar(N,T)默认是稀疏结构维度虽然对但有些运算会慢。更要命的是sdpvar(1,T)和sdpvar(T,1)方向完全不同约束里向量维度和矩阵维度对不上Yalmip会报错或者静默地做广播结果完全不对。第二如果用了分段线性化0-1变量要显式声明binary。Yalmip根据约束类型有时能自动识别但保险起见用binvar单独定义。第三约束拼接不要用cell数组逐个append再转矩阵直接在变量上累加即可。约束量大的时候Yalmip会每行单独存储用cell转矩阵反而拖慢速度。3.3 后处理出图、指标计算与结果校验求解完成后用value()函数把变量解析出来P_h_opt value(P_h); V_opt value(V); P_deliver_opt value(P_deliver); P_curtail_opt value(P_curtail);然后计算几个关键指标expect_energy mean(sum(P_deliver_opt, 2)) * 1; % MWhΔt1h curtail_rate mean(sum(P_curtail_opt, 2)) / mean(sum(P_pv_scen repmat(sum(P_h_opt,1), N_scen, 1), 2));出图建议把24小时的光伏出力、水电出力和外送功率画在同一张图里一眼就能看出光水互补效果。我还会额外画一张“所有场景外送功率带”的图把各场景下的P_deliver曲线叠成半透明线条直观展示调度方案在不确定性面前的稳定性。如果带宽很宽说明方案对光伏波动敏感需要进一步考虑备用或储能。4. 复现过程中踩过的坑与对应解法4.1 水电出力非线性函数的线性化处理第一次跑模型我直接写P_h 9.81etaQ_turb*H_netH_net又用库容线性函数近似结果求解器直接报nonconvex QP错误。问题本质是发电流量和净水头相乘形成了双线性项在优化框架里是非凸的。我的解决思路是分两步走。第一步把所有电站都改成恒定水头近似P_h K·Q_turb跑通全流程第二步对调节能力弱、水头变化明显的电站升级成“出力-流量-库容”二维查表法用分段线性拟合。这样既能保证精度又不会让模型一开始就陷进非凸泥潭。很多论文给的结果之所以好看恰恰是因为他们做了类似的线性化处理复现时不要看到原公式是非线性的就照抄。4.2 场景数与求解速度的权衡光伏场景数直接决定模型规模。我试过20个场景约束数膨胀20倍Gurobi求解时间从几十秒膨胀到十几分钟MIP gap还迟迟不收敛。后来把场景削减到15个并且检查削减后场景是否覆盖了高、中、低出力三类典型情况速度立刻回到可接受范围。这里有个实践技巧把场景共享的约束水电水量平衡、库容上下限从场景循环里提出来只写一次。这样约束数从“场景数×约束组”变成“共享约束组场景数×场景特有约束组”对求解器加速非常明显。另外目标函数里场景概率ρ_s要归一化。如果场景削减后忘了重新归一化期望电量会被系统性高估或低估和论文结果对不上时很容易让人怀疑是自己的算法错了其实只是概率没除干净。4.3 初始水位与末时段库容约束的敏感性短期调度里初始水位是已知输入末水位通常要求落在合理范围内。末水位约束设太紧会直接导致模型无解设太松模型会“聪明”地把水库在期末放空来多发电结果失去实际运行意义。我的做法是把末水位约束从单点等式改成区间区间宽度参考实际调度规程比如在死水位之上预留5%到10%库容。然后做一组敏感性扫描把末水位下限从低到高依次代入画出目标函数与末水位约束的关系曲线。这条曲线能直接告诉你“保住下一周期库容”的代价是多少。顺便说一句起始水位也必须和论文明细一致。如果论文用的是丰水期初水位3.0×10^8 m³你用了枯水期初水位1.5×10^8 m³那不管目标函数怎么调结果都不可能对得上。4.4 求解器数值问题单位、数量级与容差这个坑藏得很深。库容是10^8 m³量级发电流量是几十到几百m³/s出力是几百MW光伏是上千MW。这些数字混在同一个约束矩阵条件数会非常差Gurobi和CPLEX都会出现“numerical trouble”警告甚至误报不可行。我最终采用的量纲方案是流量统一乘Δt变成时段水量单位换算为百万立方米功率统一为MW目标函数电量用MWh。这样所有约束的系数都在0.1到1000之间求解器基本不会闹脾气。如果还是遇到数值警告再把Gurobi的feasibility tolerance调紧到1e-6同时把MIP gap设到1e-4通常就能稳定收敛。5. 怎么判断复现是否成功验证方法与实践体会5.1 与论文图表对拍复现的第一步验证是把论文里最核心的几张图——水位过程线、出力曲线、弃电率柱状图——提取出来和你的复现结果叠在一起对比。注意对比时检查三件事横轴时间单位是小时还是15分钟功率是MW还是p.u.用的典型日数据是否一致。如果对不上优先排查目标函数的Δt、流量水量是否统一、场景概率是否归一化、拓扑约束是否遗漏。很多时候不是算法问题是细节问题。5.2 极端场景下的合理性检验把光伏场景削减到只剩两种极端情况零出力场景和满发场景。此时模型退化成确定性优化结果可以手工验算零出力场景下水电应当尽可能维持高效运行外送不超过P_limit满发场景下中午时段水电主动让路光伏尽量满发弃电集中在光伏高发而通道受限的时段。我靠这个极端场景测试抓出过一个很隐蔽的bug水量平衡约束里上游来水的符号写反了导致水库入流被当成出流扣减物理上完全不成立但完整场景下很难察觉因为求解器会通过调整其他变量把结果“圆”回来。所以做任何优化模型极端场景测试都是必须的第一步不是可选项。5.3 敏感性分析验证模型稳健性最后一个层次是敏感性分析。让光伏预测误差标准差从5%逐步增加到20%观察期望可消纳电量目标值的变化。合理的结果是误差越大期望可消纳电量越低因为光伏极端场景增多弃电风险变大水电被迫预留更多调节空间。再扫一组外送通道容量通道容量逐渐增加期望可消纳电量应当单调上升然后饱和。饱和点对应的物理含义是“水电全出力光伏满发”的系统极限。如果这个趋势出不来说明不确定性约束没有真正参与决策你的模型很可能被简化过度了。这几组敏感性测试做完才能比较有底气地说这个复现结果是可信的。跑了几轮扫描之后我对这套模型最深的感受是真正决定复现成败的往往不是公式推导而是数值细节和数据口径。如果你在新数据集上第一版就跑通了大概率是运气好——别急着高兴先做一次极端场景检验再做一次敏感性扫描确认模型行为符合物理直觉这比什么技巧都重要。
基于Matlab的梯级水光互补系统短期优化调度模型复现与解析
做水电调度优化的课题通常都有同一个痛点仿真模型搭起来了但论文里的结果就是复现不出来。梯级水光互补系统的短期优化调度尤其如此——上游连着几个水电站中间嵌着大片光伏基地调度目标不是单纯的“发满”而是在光伏出力的不确定性之下系统还能可靠外送多少电。这篇文章记录的是一个基于Matlab的EI论文复现过程模型的核心是梯级水光互补系统短期优化调度模型目标函数为最大化可消纳电量期望。代码在Matlab R2023a环境下用Yalmip建模、Gurobi求解完整跑通适合正在做新能源消纳、水电优化调度方向的硕博生也适合刚接触随机优化建模的读者当一份可参考的样例。1. 这个模型到底在优化什么问题拆解与场景设定1.1 梯级水光互补调度的现实约束把物理场景先梳理清楚。所谓梯级电站就是同一流域上下游串了多个水电站上游电站的出库流量经过一段时间延迟变成下游电站的入库流量。水光互补的意思是光伏电站和水电站共用一条外送通道光伏出力波动时由水电来“托底”和“补差”。举个例子比如某流域有A、B、C三个梯级电站总装机分别约为800MW、500MW、600MW光伏装机1200MW。白天中午光伏出力逼近峰值就要求水电压低出力给光伏让出通道傍晚光伏骤减水电必须快速顶上。如果调度只按光伏预测曲线做确定性安排实际出力一旦偏离预测就会造成外送通道闲置或者弃光。更极端的情况云层遮挡导致光伏瞬时出力大幅下降水电调整不及时就会破坏功率平衡。这些现实约束决定了模型不能是简单的确定性优化必须把光伏随机性纳入决策框架。1.2 “可消纳电量期望”的准确含义标题里最容易被误解的词是“可消纳电量期望”。很多人第一反应以为“期望”就是把预测出力曲线对时间积分其实不是。在短期调度中“可消纳电量”指的是在满足外送通道容量、水电站运行约束、光伏实际出力等条件之后系统真正送出去的电量。由于光伏出力存在随机性这个电量是一个随机变量所以目标函数取它在多个场景下的期望值。数学上写成max E[∑ P_net(ω,t)·Δt]其中ω表示光伏出力场景P_net是场景ω下t时刻的外送功率。这样做的现实意义是调度方案不针对某一条光伏曲线优化而是对所有可能场景的平均表现都尽量好避免“预测准了就赚、预测不准就亏”的极端情况。代码实现时期望运算用场景法落地生成N个光伏场景给每个场景一个概率把期望展开为场景概率加权求和的线性表达式再交给求解器。这也是这类EI论文里最主流的建模思路。1.3 短期优化调度的时间尺度与决策变量短期优化调度通常指日前调度时间窗口T24小时步长Δt1小时。如果做日内滚动可能会用15分钟分辨率但原理完全一样。决策变量分两类第一类是水电相关的各电站的出力、发电流量、弃水流量、库容或水位第二类是外送与弃电相关的各时段外送功率、弃光功率。这里有一个非常关键、也是很多复现翻车的点在期望最大化模型里水电出力通常作为第一阶段决策也就是在光伏场景实现之前就要确定下来而外送、弃电是第二阶段决策等看到光伏实际出力后做调整。在随机规划里这就是“here-and-now”和“wait-and-see”的区别。代码层面体现为水电出力变量不标场景索引外送和弃电变量标场景索引。如果把所有变量都带场景下标模型就退化成一堆独立子问题的拼接优化结果没有任何实际调度意义。这一点后面我还会反复强调。2. 数学模型目标函数与约束条件的完整拆解2.1 目标函数如何把期望值写成可求解的形式模型的目标函数可以写成max ∑_{s1}^{N_s} ρ_s · ∑_{t1}^{T} P_deliver(s,t)·Δt其中P_deliver(s,t)是场景s下t时刻实际外送功率ρ_s是场景s的概率N_s是场景总数。为什么要引入P_deliver而不直接用P_hP_pv因为外送功率要同时受制于发电能力和外送通道上限也就是P_deliver(s,t) min(P_h_total(t)P_pv(s,t), P_limit(t))。min函数不可导直接放进目标函数会让求解器无从下手。标准做法是引入辅助变量P_deliver用两个约束限定它P_deliver(s,t) ≤ P_h_total(t) P_pv(s,t) P_deliver(s,t) ≤ P_limit(t) P_deliver(s,t) ≥ 0由于目标函数是最大化P_deliver之和求解器会自动把P_deliver推到两个上限的较小值等价实现了min运算。这组约束连接了第一阶段水电变量和第二阶段场景变量是整个模型“期望可消纳电量最大化”的关键桥梁。注意目标函数里的Δt。如果Δt1小时且P_deliver单位是MW那么∑P_deliver·Δt就是电量MWh。很多人复现时忘了乘Δt导致目标值和论文对不上差了一个小时量级。2.2 梯级水电约束水量平衡与水头耦合梯级水电约束是模型里工程量最大的部分核心是四组约束第一组是水量平衡V_i(t1) V_i(t) [Q_nat_i(t) Q_up_i(t-τ) - Q_turb_i(t) - Q_spill_i(t)]·Δt其中Q_nat是天然来水Q_up是上游电站总出库发电流量加弃水流量τ是水流时滞。如果上游有支流汇入Q_up还要按拓扑结构拆分。第二组是库容上下限V_min_i ≤ V_i(t) ≤ V_max_i第三组是出库流量约束Q_turb_min_i ≤ Q_turb_i(t) ≤ Q_turb_max_i第四组是出力上限P_h_i(t) ≤ 9.81·η_i·Q_turb_i(t)·H_i(t)H_i(t)是净水头等于坝前水位减尾水位再减水头损失。坝前水位是库容的非线性函数净水头又和发电流量有耦合所以水电出力本质上是“发电流量×净水头”的双线性项。在Matlab复现里最常见的有两种处理一是忽略水头变化用恒定水头近似P_h_i(t) ≈ K_i·Q_turb_i(t)模型变成线性规划跑起来非常快二是预先算好“出力-流量-库容”关系表用分段线性约束拟合精度高但会引入0-1变量变成混合整数规划。论文里如果明确说“忽略水头变化”直接线性化就行如果没提建议先用恒定水头跑通再升级成查表法。2.3 光伏出力约束与不确定性场景生成光伏出力在模型里不是决策变量而是已知场景参数。场景生成的常用做法是基于预测曲线加误差扰动P_pv(s,t) max(0, P_pv_forecast(t)·(1 ε(s,t)))其中ε(s,t)是预测误差采样值可以假设服从正态分布或Beta分布。Matlab里的生成代码大致如下N_scen 15; T 24; sigma 0.1; epsilon normrnd(0, sigma, [N_scen, T]); P_pv_scen max(0, repmat(P_pv_forecast, N_scen, 1) .* (1 epsilon)); rho ones(N_scen, 1) / N_scen;如果论文给的场景数据是历史出力曲线那就要先做场景削减K-means聚类或者同步回代削减都行削减完再把概率归一化。场景削减的核心是保“典型性”让剩下场景里同时包含高、中、低出力的代表。2.4 电网消纳通道与弃电惩罚外送通道约束是决定“可消纳”的另一只手P_deliver(s,t) ≤ P_limit(t)P_limit是外送断面容量或者电网需求上限可以是24小时的时序曲线。弃电量定义为P_curtail(s,t) ≥ P_h_total(t) P_pv(s,t) - P_limit(t) P_curtail(s,t) ≥ 0如果论文里区分弃光和弃水就把P_curtail拆成两个变量分别给惩罚权重。通常光伏弃电代价高于弃水因为光伏没有燃料成本而且白天发电窗口错过就没了。还有一个容易漏的约束是功率平衡。并网场景下应该是P_h_total(t) P_pv(s,t) - P_curtail(s,t) P_load(t)也就是系统净负荷功率平衡。如果是孤网或外送场景则替换成外送通道约束。复现时先对照论文把所有等式写全再逐个检查哪些约束是有效的、哪些是松弛的。我见过太多复现失败案例都是漏写了一个不起眼但关键的等式约束。3. Matlab实现细节从公式到可运行代码的三层结构3.1 数据准备参数表与时间序列的组织方式代码实现我先说数据层。强烈建议把所有参数封装成结构体避免脚本里散落一堆魔法数字。下面是我习惯的组织方式% hydro_data.m hydro(1).name A水库; hydro(1).Vmin 1.2e8; % 死库容m^3 hydro(1).Vmax 4.5e8; % 正常蓄水位对应库容m^3 hydro(1).V0 3.0e8; % 初始库容m^3 hydro(1).Qmin 0; hydro(1).Qmax 250; % 最大发电流量m^3/s hydro(1).eta 0.88; % 综合效率系数 hydro(1).maxPower 800; % 装机容量MW hydro(1).K 8.5; % 出力系数恒定水头近似时使用 % 光伏与电网 pv.capacity 1200; % 光伏装机MW pv.forecast load(pv_forecast.mat); % 24小时预测出力MW grid.Plimit 1800; % 外送通道上限MW grid.load load(load_curve.mat); % 净负荷曲线MW这里有个小细节库容单位是m³发电流量单位是m³/s功率单位是MW。约束里把这些不同量纲的数值混在一起数量级可能差10^8求解器数值稳定性会很差。我的处理方式是把所有流量乘以Δt换算成时段内水量单位统一成“百万立方米”功率统一为MW目标函数的Δt用小时。这样各类数值都在几十到几千的量级Gurobi和CPLEX跑起来都舒服很多。3.2 Yalmip建模变量定义、约束装配与求解器调用建模层是代码的核心用Yalmip把第二章的数学模型逐句翻译成约束。核心框架如下% 变量定义 P_h sdpvar(N_hydro, T, full); % 第一阶段水电出力 Q_turb sdpvar(N_hydro, T, full); % 发电流量 Q_spill sdpvar(N_hydro, T, full); % 弃水流量 V sdpvar(N_hydro, T1, full); % 库容T1因为要覆盖首末时段 P_deliver sdpvar(N_scen, T, full); % 第二阶段外送功率 P_curtail sdpvar(N_scen, T, full); % 弃电功率 % 目标函数Yalmip默认求最小化所以取负号 objective -sum(sum(repmat(rho, 1, T) .* P_deliver)); % 约束装配 C []; % 水量平衡约束、库容约束、出力约束... C [C, V(:, t1) V(:, t) (Q_nat Q_up - Q_turb - Q_spill) * dt]; C [C, V(:, t) Vmin, V(:, t) Vmax]; % P_deliver与发电能力、外送上限的关系 C [C, P_deliver(s, t) sum(P_h(:, t)) P_pv_scen(s, t)]; C [C, P_deliver(s, t) Plimit(t)]; C [C, P_deliver(s, t) 0]; % 求解 ops sdpsettings(solver, gurobi, verbose, 2, mipgap, 1e-4); result optimize(C, objective, ops);这里有三个容易踩的坑第一sdpvar的维度方向。sdpvar(N,T,full)生成的是N行T列矩阵而sdpvar(N,T)默认是稀疏结构维度虽然对但有些运算会慢。更要命的是sdpvar(1,T)和sdpvar(T,1)方向完全不同约束里向量维度和矩阵维度对不上Yalmip会报错或者静默地做广播结果完全不对。第二如果用了分段线性化0-1变量要显式声明binary。Yalmip根据约束类型有时能自动识别但保险起见用binvar单独定义。第三约束拼接不要用cell数组逐个append再转矩阵直接在变量上累加即可。约束量大的时候Yalmip会每行单独存储用cell转矩阵反而拖慢速度。3.3 后处理出图、指标计算与结果校验求解完成后用value()函数把变量解析出来P_h_opt value(P_h); V_opt value(V); P_deliver_opt value(P_deliver); P_curtail_opt value(P_curtail);然后计算几个关键指标expect_energy mean(sum(P_deliver_opt, 2)) * 1; % MWhΔt1h curtail_rate mean(sum(P_curtail_opt, 2)) / mean(sum(P_pv_scen repmat(sum(P_h_opt,1), N_scen, 1), 2));出图建议把24小时的光伏出力、水电出力和外送功率画在同一张图里一眼就能看出光水互补效果。我还会额外画一张“所有场景外送功率带”的图把各场景下的P_deliver曲线叠成半透明线条直观展示调度方案在不确定性面前的稳定性。如果带宽很宽说明方案对光伏波动敏感需要进一步考虑备用或储能。4. 复现过程中踩过的坑与对应解法4.1 水电出力非线性函数的线性化处理第一次跑模型我直接写P_h 9.81etaQ_turb*H_netH_net又用库容线性函数近似结果求解器直接报nonconvex QP错误。问题本质是发电流量和净水头相乘形成了双线性项在优化框架里是非凸的。我的解决思路是分两步走。第一步把所有电站都改成恒定水头近似P_h K·Q_turb跑通全流程第二步对调节能力弱、水头变化明显的电站升级成“出力-流量-库容”二维查表法用分段线性拟合。这样既能保证精度又不会让模型一开始就陷进非凸泥潭。很多论文给的结果之所以好看恰恰是因为他们做了类似的线性化处理复现时不要看到原公式是非线性的就照抄。4.2 场景数与求解速度的权衡光伏场景数直接决定模型规模。我试过20个场景约束数膨胀20倍Gurobi求解时间从几十秒膨胀到十几分钟MIP gap还迟迟不收敛。后来把场景削减到15个并且检查削减后场景是否覆盖了高、中、低出力三类典型情况速度立刻回到可接受范围。这里有个实践技巧把场景共享的约束水电水量平衡、库容上下限从场景循环里提出来只写一次。这样约束数从“场景数×约束组”变成“共享约束组场景数×场景特有约束组”对求解器加速非常明显。另外目标函数里场景概率ρ_s要归一化。如果场景削减后忘了重新归一化期望电量会被系统性高估或低估和论文结果对不上时很容易让人怀疑是自己的算法错了其实只是概率没除干净。4.3 初始水位与末时段库容约束的敏感性短期调度里初始水位是已知输入末水位通常要求落在合理范围内。末水位约束设太紧会直接导致模型无解设太松模型会“聪明”地把水库在期末放空来多发电结果失去实际运行意义。我的做法是把末水位约束从单点等式改成区间区间宽度参考实际调度规程比如在死水位之上预留5%到10%库容。然后做一组敏感性扫描把末水位下限从低到高依次代入画出目标函数与末水位约束的关系曲线。这条曲线能直接告诉你“保住下一周期库容”的代价是多少。顺便说一句起始水位也必须和论文明细一致。如果论文用的是丰水期初水位3.0×10^8 m³你用了枯水期初水位1.5×10^8 m³那不管目标函数怎么调结果都不可能对得上。4.4 求解器数值问题单位、数量级与容差这个坑藏得很深。库容是10^8 m³量级发电流量是几十到几百m³/s出力是几百MW光伏是上千MW。这些数字混在同一个约束矩阵条件数会非常差Gurobi和CPLEX都会出现“numerical trouble”警告甚至误报不可行。我最终采用的量纲方案是流量统一乘Δt变成时段水量单位换算为百万立方米功率统一为MW目标函数电量用MWh。这样所有约束的系数都在0.1到1000之间求解器基本不会闹脾气。如果还是遇到数值警告再把Gurobi的feasibility tolerance调紧到1e-6同时把MIP gap设到1e-4通常就能稳定收敛。5. 怎么判断复现是否成功验证方法与实践体会5.1 与论文图表对拍复现的第一步验证是把论文里最核心的几张图——水位过程线、出力曲线、弃电率柱状图——提取出来和你的复现结果叠在一起对比。注意对比时检查三件事横轴时间单位是小时还是15分钟功率是MW还是p.u.用的典型日数据是否一致。如果对不上优先排查目标函数的Δt、流量水量是否统一、场景概率是否归一化、拓扑约束是否遗漏。很多时候不是算法问题是细节问题。5.2 极端场景下的合理性检验把光伏场景削减到只剩两种极端情况零出力场景和满发场景。此时模型退化成确定性优化结果可以手工验算零出力场景下水电应当尽可能维持高效运行外送不超过P_limit满发场景下中午时段水电主动让路光伏尽量满发弃电集中在光伏高发而通道受限的时段。我靠这个极端场景测试抓出过一个很隐蔽的bug水量平衡约束里上游来水的符号写反了导致水库入流被当成出流扣减物理上完全不成立但完整场景下很难察觉因为求解器会通过调整其他变量把结果“圆”回来。所以做任何优化模型极端场景测试都是必须的第一步不是可选项。5.3 敏感性分析验证模型稳健性最后一个层次是敏感性分析。让光伏预测误差标准差从5%逐步增加到20%观察期望可消纳电量目标值的变化。合理的结果是误差越大期望可消纳电量越低因为光伏极端场景增多弃电风险变大水电被迫预留更多调节空间。再扫一组外送通道容量通道容量逐渐增加期望可消纳电量应当单调上升然后饱和。饱和点对应的物理含义是“水电全出力光伏满发”的系统极限。如果这个趋势出不来说明不确定性约束没有真正参与决策你的模型很可能被简化过度了。这几组敏感性测试做完才能比较有底气地说这个复现结果是可信的。跑了几轮扫描之后我对这套模型最深的感受是真正决定复现成败的往往不是公式推导而是数值细节和数据口径。如果你在新数据集上第一版就跑通了大概率是运气好——别急着高兴先做一次极端场景检验再做一次敏感性扫描确认模型行为符合物理直觉这比什么技巧都重要。
