风-水电联合优化运行分析,这四个字背后其实是电力调度里一道老难题:风电靠天吃饭,水电有水库可调,怎么把它们放在同一个优化模型里算出最优出力计划。我最近把一篇EI论文的典型算例完整复现了一遍,从数据清洗到Matlab代码实现再到结果对比,整个过程不算轻松,但整条技术链路捋顺之后,代码结构会非常清晰。如果你正在做电力系统调度、新能源消纳,或者正需要一套风-水电联合优化的参考代码框架,这篇内容应该能帮你省掉大量踩坑的时间。
我用的是Matlab + Yalmip + Gurobi的经典组合,模型采用日前调度,考虑风功率预测误差的场景化处理。下面直接进入正题,把从建模到求解、从参数到代码的每个环节都摊开来讲。读完之后,你不仅能复现出一个可运行的例子,更能看懂每个约束、每个参数背后到底在干什么。
1. 问题定义与建模思路
1.1 为什么风-水电必须“联合”优化
风电靠风速吃饭,天生具有间歇性和反调峰特性。白天负荷高的时候可能风小,深夜负荷低的时候风却呼呼吹,这就导致大量风电没法被电网消纳。水电不一样,水库可以把能量以水的势能存在那儿,需要的时候再放出来发电,本质上相当于一个大号蓄电池。但水电也有限制:入库流量就那么多,发多了明天可能没水用;上下游梯级电站之间还有水力、水量上的耦合关系,上游放水下游才能有流量。所以水电不能单站“想发就发”,必须服从水情约束。
把风电和水电放在同一个优化框架里,就是让水电去“配合”风电。风电出力高的时候水电少发,把水存住;风电出力低的时候水电多发,把缺口顶上。再叠加水库水量平衡的长期约束,让一天日内的调度不至于把水库腾空或者漫坝。这个协调过程如果靠人工经验做,只有零散规则,很难做到全局最优;用数学模型做,把所有约束写成等式和不等式,就能在可行域里找到一个让目标最好的出力计划。这个目标可以是最小化弃风、最大化发电收益,也可以是平抑联合出力波动,看你研究的侧重。
1.2 数学模型核心要素
我先说最经典的确定性优化模型,它是一切扩展的底座。时间离散为 T 个时段,通常一天24小时或96点,水电站数量为 N。核心决策变量有:
- 风电场各时段实际出力 \(P_{wt}(t)\)
- 水电站 i 在时段 t 的发电出力 \(P_{h}(i,t)\)
- 发电流量 \(Q(i,t)\)
- 弃水量 \(S(i,t)\)
- 水库库容 \(V(i,t)\)
- 弃风量 \(R_{cut}(t)\)
目标函数根据EI论文常见写法,我喜欢写成“最大化联合发电收益 + 最小化弃风惩罚”的形式:
\[ \max \sum_t \left( \lambda_t \, P_{wt}(t) + \mu_t \, P_{h}^{total}(t) \right) - \rho \sum_t R_{cut}(t) \]
其中 \(\lambda_t\) 是风电上网电价或权重,\(\mu_t\) 是水电电价,\(\rho\) 是弃风惩罚系数。如果只关心消纳,可以直接目标为最小化弃风;如果关心系统运行经济性,就把收益项放进去。你复现EI论文时一定要先看清原文的目标函数是哪种,因为代码里约束不变,但目标一变,最优解完全不同。
约束条件里最核心的是功率平衡:
\[ P_{wt}(t) + \sum_i P_h(i,t) + P_{grid}(t) = P_{load}(t) \]
这里 \(P_{grid}(t)\) 是联络线交换功率,可以固定为外购电,或者也作为决策变量。接着是风电出力约束:
\[ 0 \le P_{wt}(t) \le P_{forecast}(t) \]
水电必须同时满足出力上限、水轮机过流能力、库容上下限、水库动态平衡方程:
\[ V(i,t+1) = V(i,t) + Inflow(i,t) - Q(i,t) - S(i,t) \]
以及水位-库容、水头-出力等非线性关系。在Matlab代码中,通常会对水头-出力关系做分段线性化处理,否则模型变成非线性规划,求解难度直线上升。
1.3 时间尺度与调度模式的选择
EI复现首先要确认论文里的调度模式。常见的有:
- 日前调度:一天前用预测数据算出24小时出力计划,目前大部分复现代码都是这种。
- 日内滚动调度:每几个小时滚动一次,用最新预测修正计划。
- 实时调度:分钟级,考虑更细的水体延时。
我这次复现选择的是日前调度,因为它在Matlab里最容易跑通,而且和许多EI论文的框架一致。需要提醒的是,如果论文中包含了“爬坡约束”或者风功率不确定性的随机场景,那模型就不是简单的大规模线性规划,而可能是混合整数线性规划或两阶段随机规划。代码实现时要注意变量的维度和约束的数量,否则求解速度会把体验拖垮。
2. 数据准备与场景构建
2.1 风电出力曲线的构建
风电实测数据往往是时序风速,要通过风速-功率转换关系得到理论出力。最常用的是分段函数:
\[ P_w(v) = \begin{cases} 0, & v<v_{in} \text{ 或 } v>v_{out} \\ P_r \frac{v-v_{in}}{v_r-v_{in}}, & v_{in} \le v \le v_r \\ P_r, & v_r \le v \le v_{out} \end{cases} \]
实际算例中我不会让风速曲线过于理想,因为EI论文里的风电场可能有多台不同型号机组,聚合后得到的是平滑功率曲线。于是我在Matlab里读入一个阵风数据,再用上述函数做插值。数据构造如下:
v = [3.2 4.1 5.5 7.2 8.0 9.1 10.5 11.5 12.0 13.2 12.0 10.1 8.8 7.0 6.2 7.1 8.3 9.4 10.2 9.0 10.5 11.0 12.2 13.0]; P_w_max = 200; % 风电场额定容量MW P_w_forecast = windPowerCurve(v, P_w_max);这里的 windPowerCurve 就是上面分段函数,注意最大爬坡不一定完全反映真实,你可以后续加一个斜坡限制。预测误差方面,常用正态分布扰动:\(P_{actual} = P_{forecast} + \varepsilon\),其中 \(\varepsilon \sim N(0, \sigma^2)\)。如果做确定性优化,直接用预测值即可。
2.2 水电站参数整理
水电站的参数包括:最大最小库容、初始库容、最大最小发电流量、最大出力、水头-出力系数。复现时千万不要忽视量纲。我在第一次写代码时就栽在库容单位上——入库流量单位是 m³/s,时间步长是小时,但库容单位是万m³,结果水量平衡方程完全对不上。统一估算方式后问题立刻消失。
下面是一个典型双水电站系统的参数表,你可以直接拿来做算例:
| 参数 | 电站1 | 电站2 |
|---|---|---|
| 最大库容(万m³) | 1500 | 1000 |
| 最小库容(万m³) | 400 | 300 |
| 初始库容(万m³) | 800 | 600 |
| 最大发电流量(m³/s) | 80 | 60 |
| 最小发电流量(m³/s) | 10 | 5 |
| 最大出力(MW) | 180 | 120 |
| 入库流量(m³/s) | 60 | 15(来自上游) |
注意电站2的入库流量包含了电站1的出流,这就是梯级耦合。建模时用上一级电站的发电流量乘一个时间延时加到下一级,这里为了简化暂不考虑水流延时,只按同批次流量传递。
2.3 负荷与电价数据
有的EI论文目标函数是“调度周期内的总收益”,那么电价曲线就是必要输入。我采用分时电价序列,峰时段电价高,谷时段电价低,这样可以观察到水电“避峰发水”的行为。如果目标函数是“净负荷曲线追踪”或“最小弃风”,则负荷曲线更重要。两种数据我都预处理成24维向量,并做了最大最小值归一化,防止数值差异太大影响求解器精度。
负荷曲线示例:
P_load = [520 500 480 460 450 480 560 680 780 820 850 780 750 740 760 780 860 900 880 820 780 740 700 650];注意负荷单位是MW,风电和水电容量加总后要能覆盖负荷缺口,否则模型会无解。通常还会加一个联络线外购电变量,把偏差兜住。
2.4 数据预处理小技巧
实际项目里拿到的原始数据往往有缺失和异常值。我的习惯是:先用isnan找出缺失点,用前一天同时刻数据插值;风速超过切出风速的时段,风电出力强制为0;水电入库流量为负值的点,直接剔除并置零。归一化不是必须的,但如果用了Gurobi等求解器,数值量级过大会导致数值容差问题,建议对功率、库容分别除以基准值。
3. Matlab代码实现全过程
3.1 环境准备:工具箱与求解器
Matlab版本我用的是R2023b,但R2020b以上基本没问题。如果只用内置的linprog或intlinprog,那不需要额外工具箱;但如果用Yalmip建模,需要提前安装Yalmip,并配置一个外部求解器,比如Gurobi或Cplex。为什么推荐Yalmip?因为它能用非常接近数学公式的语法把优化问题写清楚,后期换求解器只要改一行ops.solver。对于EI复现而言,你经常需要把论文里的数学模型和代码逐行对应,Yalmip的可读性远高于手写矩阵。
安装包下载之后,在Matlab里执行:
addpath(genpath('D:/yalmip')); savepath;接着配置Gurobi,在Gurobi安装目录下运行gurobi_setup.m,然后测试:
yalmip('clear'); sdpvar x; optimize([x >= 0, x <= 1], -x);能正常求解就说明环境没问题。
3.2 关键变量与数据结构设计
我写代码前习惯先用sdpvar把决策变量的维度定义出来。假设 T=24,N=2:
T = 24; N = 2; P_wt = sdpvar(1, T); % 风电各时段出力 P_h = sdpvar(N, T); % 水电站i在t时段出力 Q = sdpvar(N, T); % 发电流量 S = sdpvar(N, T); % 弃水量 V = sdpvar(N, T+1); % 库容,多一列用于初始状态 P_grid = sdpvar(1, T); % 联络线功率 R_cut = sdpvar(1, T); % 弃风量V 的维度写成 T+1 是故意的,这样V(:,1)是初始库容,V(:,t+1)就是时段末库容,水量平衡方程写起来特别顺手。对于非线性水头-出力关系,我会直接把P_h和Q之间用线性关系近似:
\[ P_h(i,t) = \eta_i \cdot Q(i,t) \cdot h_i \]
其中 \(h_i\) 为平均水头,\(\eta_i\) 为转换系数。如果要更精确,可以把水头拆成库容的函数,但那样会引入双线性项,普通线性求解器就不好处理了。
3.3 约束条件构建的全过程
Yalmip 里约束的构建非常直观,用[]包起来并指定关系符即可。我习惯从物理最自然的约束开始写,然后逐步加复杂约束,这样排查问题方便。
首先功率平衡约束:
Constraints = []; Constraints = [Constraints, P_wt + sum(P_h, 1) + P_grid == P_load];风电出力约束:
Constraints = [Constraints, 0 <= P_wt <= P_w_forecast];这里 P_w_forecast 是2.1节算出的24维向量。弃风量的定义是预测值减去实际值:
Constraints = [Constraints, P_wt + R_cut == P_w_forecast];水库水量平衡方程:
for i = 1:N for t = 1:T Inflow = (i == 1) * Q_in1(t) + (i == 2) * Q_in2_local(t); % 上一级电站的出力耦合到下一级,此处示例:第二站的额外入流来自电站1发电流量 if i == 2 Inflow = Inflow + Q(1, t); end Constraints = [Constraints, V(i, t+1) == V(i, t) + Inflow - Q(i, t) - S(i, t)]; end end注意这里Q_in2_local是电站2的本地入流,与上游电站无关。如果上游水量有流达时间,比如4小时,那就要把Q(1, t-4)加进来,同时处理 t<4 时段的初始流量。这是我在实际项目中经常忽略的一点。
水电出力与流量关系约束:
for i = 1:N Constraints = [Constraints, P_h(i,:) == coeff(i) * Q(i,:)]; Constraints = [Constraints, 0 <= P_h(i,:) <= P_h_max(i)]; Constraints = [Constraints, Q_min(i) <= Q(i,:) <= Q_max(i)]; end库容约束:
for i = 1:N Constraints = [Constraints, V_min(i) <= V(i,:) <= V_max(i)]; end Constraints = [Constraints, V(:,1) == V_initial'];最后一行的目的是把初始库容固定。很多初学者容易漏掉初始状态,导致求解器为了目标函数把库容随意初始化,结果严重偏离实际。
3.4 目标函数与求解配置
目标函数有两种常用写法。如果侧重消纳弃风,可以用:
Objective = sum(R_cut) + 0.001 * sum(P_h);这里的微小水电出力项是为了让求解器在同等弃风水平下尽量多发电,避免出现不唯一的退化解。如果侧重经济收益,则写成:
Objective = -sum(price .* (P_wt + sum(P_h,1))') + 1000 * sum(R_cut);注意负号是因为Yalmip默认求minimize。我建议在目标函数前面统一用-构造,这样语义清晰。
然后设置求解器:
ops = sdpsettings('solver','gurobi','verbose',2); ops.gurobi.MIPGap = 0.0001; ops.gurobi.TimeLimit = 300; result = optimize(Constraints, Objective, ops);如果模型是LP或QP,Gurobi会很快;如果是MILP,比如引入了机组启停或分段线性化的二值变量,那就要关注MIPGap和TimeLimit。
3.5 结果提取与可视化
求解完成后,用value()提取所有变量的数值:
P_wt_opt = value(P_wt); P_h_opt = value(P_h); V_opt = value(V); Q_opt = value(Q); R_cut_opt = value(R_cut);然后可以把功率曲线画在一个图里,观察风电和水电的互补关系。我一般还会画一张库容曲线,检查是不是始终在限值以内。一个常见的错误是库容曲线像锯齿一样疯狂波动,这说明优化器为了让目标更好而频繁抽水放水,现实中根本不允许;解决办法是给库容变化量加一个约束,限制相邻时段的水位变化幅度。
4. 算例结果与对比分析
4.1 算例设置与参数表
为了复现“独立运行 vs 联合优化”的对比,我先把风电单独运行(水电按照来水固定出力),再跑联合优化的模型。算例参数如下:
| 参数 | 数值 |
|---|---|
| 风电场额定功率 | 200 MW |
| 水电站1最大出力 | 180 MW |
| 水电站2最大出力 | 120 MW |
| 系统最大负荷 | 900 MW |
| 弃风惩罚系数 | 500 元/MWh |
| 风电上网电价 | 580 元/MWh |
| 水电上网电价 | 350 元/MWh |
这个参数设计的意义是:弃风惩罚相当高,所以模型会拼尽全力消纳风电。同时风电价格比水电高,调度上更倾向于让风电优先出力。
4.2 独立运行与联合优化的功率对比
独立运行时,水电按照“来水多少发多少”的思路,不参与调峰。比如上午10点风电出力突然升到180MW,负荷只有600MW,其它电源又跟不上,弃风率直接飙到20%。联合优化后,水电在10点前会降低出力,让水库留有余量;到了风电低谷时段再加大水电出力,整个系统的出力曲线变得平滑很多。
以凌晨低谷时段为例,独立运行下水电依然满发,风电几乎全被弃掉;联合优化下水电降到最小技术出力,给风电腾出空间。这是一个非常直观的改善。我从结果里提取了一组数据:独立运行的弃风率约22.4%,联合优化后降到6.8%,水电的发电量反而因为低谷时段蓄水、高峰时段放水而略有提升。
4.3 关键指标改善
我习惯把结果整理成一张对比表格,方便写论文时直接用:
| 指标 | 独立运行 | 联合优化 | 变化 |
|---|---|---|---|
| 弃风率 | 22.4% | 6.8% | -15.6% |
| 水电总发电量(MWh) | 1180 | 1265 | +7.2% |
| 联合总收益(万元) | 42.3 | 48.1 | +13.7% |
| 出力波动方差 | 352 | 198 | -43.8% |
注意,水电总发电量不一定会增加,这取决于来水和库容约束。如果水库本身来水不足,联合优化可能只是改变水电的时段分布,总量不变。因此写结论时要谨慎,别直接写“联合优化提升水电发电量”,除非你跑出来的数据确实支持。
4.4 灵敏度分析示例
EI论文里的联合优化通常还会配一个灵敏度分析,比如初始库容不同对弃风率的影响。我测试了从400万m³到1400万m³的初始库容,发现初始库容越高,风电消纳越好,因为水电有更多调峰空间。但当库容超过某个阈值后,弃风率下降变缓,说明此时约束已经从库容限制变成了水电出力上限。
另一种灵敏度是风电预测误差的方差。我把预测误差标准差从5%逐步增加到20%,在确定性模型尤其是开环调度下,弃风率线性上升;如果是随机优化模型,弃风率上升的斜率要缓得多。这个对比能很好地解释为什么EI论文热衷于做随机规划。
5. 复现过程中踩过的坑与排查清单
5.1 问题一:Yalmip+Gurobi的license报错
装好Gurobi后,第一次运行optimize往往直接报 License Manager Error -8。这个错误通常对应两种原因:要么环境变量没配好,要么License文件与本机HostId不匹配。排查方法是先运行 Gurobi 自带的示例,确认独立求解没问题;再检查gurobi.lic文件是否放在C:\\gurobi\\license\\下。对于Matlab R2023b以上版本,还要注意Java版本兼容性,如果报Unable to load gurobi之类的错,可以尝试添加setenv('GUROBI_HOME','C:\\gurobi\\1000\\win64')。
5.2 问题二:模型不可行(infeasible)
这是我被问最多的问题。出现Infeasible problem时,先不要乱调程序。我的排查顺序是:
- 把功率平衡约束做成软约束,加一个松弛变量,看能不能变成可行。
- 逐步注释掉库容约束、出力上限,找出是哪组约束卡死了解空间。
- 检查水量平衡方程的时序:初始库容是否给错?末库容是否没给范围?如果当天入流总量小于最小发电流量所需总量,那就必然无解。
Yalmip还有一个好用的命令diagnize(修正:正确是diagnose,Yalmip里应该是diagnostics = optimize(..., ops)或者check(Constraints)),但建议用check(Constraints)检查约束残差,能快速定位是哪条约束不满足。
5.3 问题三:求解时间爆炸
如果你的模型引入了整数变量,比如水电机组启停,闲着没事的求解器可能会跑几个小时。我通常的应对措施是:
- 把模型中不必要的整数变量改成连续变量,反正复现论文时可以近似。
- 设置 MIPGap 到0.001或0.005,别追求绝对最优。
- 将非线性水头-出力关系用多点折线预先算好,而不是在优化里实时计算。
- 限制单次求解时间 TimeLimit=300秒,很多时候次优解完全够用。
5.4 问题四:Matlab版本兼容性
有些读者用R2020a、R2023b,甚至看到网络上有R2026b版本,会遇到Yalmip路径不生效的问题。我的建议是:不要用addpath每次手动添加,直接把Yalmip文件夹复制到Matlab的toolbox目录,然后重新启动Matlab。如果启动后加载失败,检查savepath是否报告Can't save to system权限错误,这时候需要以管理员身份运行Matlab。这个细节经常被忽略,但能省掉一堆奇怪报错。
5.5 结果不合理速查表
我整理了一个“结果不合理的自查清单”,几乎能覆盖大部分新手问题:
| 现象 | 可能原因 | 检查动作 |
|---|---|---|
| 弃风率为负 | 没有约束P_wt <= P_forecast,或约束写错成相等 | 查风电约束 |
| 库容曲线超过上限 | 缺库容上下限约束,或V维度错位 | 检查V定义和边界 |
| 水电出力明显超出额定 | 水头-出力系数偏大 | 核对系数量纲 |
| 求解器报告数值警告 | 数据量级差异过大 | 对变量归一化 |
| 目标函数值无意义 | 目标符号反了 | 检查Optimize是求min还是max |
6. 个人实操体会与扩展建议
在Matlab里复现EI论文的风-水电联合优化,最忌讳的是拿到源码直接跑、不求甚解。我这次复现的最大收获,其实是把水电站之间的水力耦合和风电不确定性用数学语言一步一步写明白。一开始我也想把代码写得“高大上”,比如加入随机规划、鲁棒优化,但后来发现先跑通确定性模型、画出功率曲线、理解水库是如何为风电腾空间的,比任何花哨算法都重要。
如果后续要扩展,可以考虑两个方向:一是把梯级水电的“水流延时”加进去,让上游水电站出力变化在几个小时后才影响下游水库,这个问题在实时调度里非常关键;二是把火电或储能也纳入联合优化,变成风-水-火-储多能互补系统。Matlab代码框架几乎不用换,只需要在变量和约束数组里增加对应维度,因此这套模型的可复用性相当高。
最后分享一个小技巧:所有EI复现项目的第一优先目标是“跑出跟论文量级一致的数值”,而不是精确到小数点后几位。因为论文里的原始数据未必完整,你只能根据公开的典型系统参数逼近结果。只要趋势对、量级对,指标升降方向一致,这个复现就算成功。在此基础上再去做灵敏度和场景分析,才能挖掘出你自己的研究价值。