做风-水电联合优化的人应该都有同感：单独看风电，随机性、间歇性大到让人头疼；单独看水电，又受来水、水库调节能力和生态流量约束。但把两者放在同一个调度框架里做联合优化运行分析，往往能在几乎不新增硬件投资的情况下，把系统的灵活性和可再生能源消纳空间拉高一大截。用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(t+1) = 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。这里给你一张可以直接抄作业的参数表模板：
| 参数 | 符号 | 数值 | 单位 |
|---|---|---|---|
| 调度时段数 | T | 24 | h |
| 风电装机 | cap_w | 300 | MW |
| 风电预测出力均值 | Pw_avg | 120 | MW |
| 风电预测误差比例 | sigma_w | 0.15 | % |
| 水库初始库容 | V_init | 1.2 | 亿m3 |
| 水库最大库容 | V_max | 2.0 | 亿m3 |
| 水库最小库容 | V_min | 0.5 | 亿m3 |
| 天然入库流量 | Q_in | 40（随季节变化） | m3/s |
| 水电最大出力 | Ph_max | 120 | MW |
| 水电最小出力 | Ph_min | 10 | MW |
| 综合出力系数 | K | 8.5 | - |
| 参考水头 | H_ref | 30 | m |
| 弃风惩罚系数 | lambda_curt | 300 | 元/MWh |
| 风电电价 | rho_w | 350 | 元/MWh |
| 水电电价 | rho_h | 280 | 元/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_h+n_w+1 : 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_h+t) = 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(t+1) - 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(T+1:2*T); P_curt_opt = x_opt(2*T+1: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，如果直接乘dt=3600，得到的是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这套环境的优势就在这儿：建模灵活，代码复用率高，换个场景就是一篇新的算例分析。
