最近要做一篇EI论文的复现，方向是风-水电联合优化运行分析，平台选Matlab。论文前前后后读了三遍，代码写了五个晚上，中间推翻过一次建模思路，最后总算把输出曲线和论文图表对上了。这篇就把完整的复现过程摊开讲一讲，包含模型怎么从论文里提炼出来、Matlab代码怎么组织、求解器怎么选、结果怎么验证，以及那些只有真跑过代码才会踩到的问题。如果你正打算复现一篇EI论文，或者刚开始接触电力系统优化调度，这篇内容应该用得上。
这类模型解决的现实问题很直接。风电场出力受风速影响，波动大且预测难，半夜风大时电网常常消纳不完，只能弃风；水电站响应速度快，可以在十几分钟内调整出力，正好用来和风电打配合。把风电和水电放进同一个优化框架，让它们按照负荷需求协同安排发电计划，就是风-水电联合优化运行分析要解决的核心问题。它适合拿来练习复现，是因为模型本身不算大，却包含了时序耦合、不确定性和非线性处理三类难点，麻雀虽小五脏俱全。
下面按我实际操作的顺序来讲：先讲拿到论文后怎么读，再讲风、水两个子系统的建模，然后是优化模型和求解器的组合，最后是验证方法和踩坑记录。
1. 先从论文里读出目标函数、约束条件和默认假设
1.1 决策变量和目标函数要画成一张关系图
拿到论文后最忌讳的事，就是急着打开编辑器开始写代码。我自己的习惯是先把数学模型部分出现的全部符号摘出来，整理成一张变量表：哪些是已知参数，哪些是决策变量，哪些是中间计算量。联合优化模型里，决策变量一般是水电站的发电流量、弃水量、库容状态、外购电力，目标函数则通常是最小化系统总运行成本，成本项包括购电费用、弃风惩罚、缺电惩罚，有时还叠加碳排放成本。
搞清目标函数的时候，可以顺手列一张成本项对应表，写代码时直接把各项写进目标表达式，不容易漏：
| 成本项 | 对应变量 | 典型量纲 |
|---|---|---|
| 外购电力成本 | 每个时段从电网购电的功率 | 元/MWh |
| 弃风惩罚 | 每个时段弃风功率 | 元/MWh |
| 缺电惩罚 | 每个时段失负荷功率 | 元/MWh |
| 碳排放成本 | 外购电对应排放量 | 元/t |
这一步千万别跳过。EI论文的数学模型通常是压缩过的，符号复用很常见，同一个字母在不同章节里可能含义完全不同。把变量关系图和数学表达式对应上，后面写约束才不会返工。我第一次复现时就因为漏标了一组中间变量，导致第一版代码跑出来的出力曲线乱成一团，白白浪费了半个晚上。
1.2 约束条件里真正难处理的是时序耦合约束
风电的约束很好写：每个时段出力在0和可用功率之间。外购电力和负荷平衡也很容易表达成线性等式。真正的骨头在水电部分——水量平衡约束把24个时段像锁链一样串起来：
V(t+1) = V(t) + (I(t) - Q(t) - S(t)) × Δt
V是库容，I是天然来水，Q是发电流量，S是弃水流量，Δt是时段长度，一般取1小时。这个约束意味着t时段的水电站决策会直接影响后面所有时段的库容状态，属于典型的时序耦合。再加上库容上下限、发电流量上限、出力上下限，以及调度期末库容约束——很多论文要求末库容等于初库容，体现水库的可持续运行——水力子系统的约束数量一下子就上去了。
处理方式上用循环逐时段追加约束即可，不用手写复杂矩阵。但要注意变量维度和索引对齐：库容变量最好定义成T+1维，循环里t和t+1别写反。我见过不少人在这类索引问题上卡很久，其实只要打印一次维度信息就能立刻发现。
1.3 论文没写出来的简化假设，恰恰是复现成败的分水岭
这是整个复现过程中最微妙的环节。EI论文篇幅有限，不会把所有假设都写进正文，但复现者必须把它们找出来，否则结果永远对不上。据我观察，风-水电联合优化方向的常见隐藏假设至少有这几项：
- 风速预测用的是点预测还是场景集。如果是随机规划论文，会有一组典型场景和对应概率，模型结构完全不同。
- 水电站水头取常数还是随库容变化。这直接决定出力方程是线性还是非线性。
- 是否考虑梯级电站联动。单库模型和梯级模型的约束结构差异很大。
- 调度周期的颗粒度是1小时还是15分钟。科研复现常取1小时，工程应用通常需要更细。
我建议在代码文件头部写一段注释，把上述假设逐一列出来。这样不仅自己心里有数，代码给别人看的时候也能直接理解建模边界。很多复现代码没法用，问题就出在读者完全搞不清作者在什么前提下建模。
2. 风电不确定性与水电调节模型在Matlab里的落地方式
2.1 从风速到风功率：一条曲线解决的事别写一串if-else
风电模型第一步是生成或读取风速序列。风速一般用两参数Weibull分布描述，Matlab里用wblrnd函数可以直接生成。功率曲线是标准的分段函数：低于切入风速出力为0，在切入风速和额定风速之间近似线性上升，额定风速以上保持满发，超过切出风速停机保护。
v_ci = 3; v_r = 12; v_co = 25; P_r = 50; wind_speed = wblrnd(8, 2, [T, 1]); % 威布尔分布随机风速 Pw_avail = zeros(T, 1); idx_linear = (wind_speed >= v_ci) & (wind_speed < v_r); idx_rated = (wind_speed >= v_r) & (wind_speed < v_co); Pw_avail(idx_linear) = P_r .* (wind_speed(idx_linear) - v_ci) / (v_r - v_ci); Pw_avail(idx_rated) = P_r;向量化写法比在for循环里堆if-else干净得多，速度也快。Pw_avail就是不确定性模型和优化模型之间的接口：在确定性模型里，它作为风电出力上限；在随机规划模型里，它变成多维数组。要注意这里的P_r是风电场装机容量，不是单机容量，量纲上不要搞混。
2.2 水量平衡和出力方程：线性还是非线性取决于水头假设
水电站的出力和发电流量、有效水头直接相关：P_h = η × ρ × g × Q × H。如果论文假设水头恒定，那么η × ρ × g × H可以合并成一个常数系数K，出力方程变成P_h = K × Q，这是个漂亮的线性项。如果水头随库容变化，H就是关于库容V的函数，表达式会带上非线性，求解难度和模型表达的信息量完全不同。
复现时优先按照论文的原始假设来。如果论文把水头简化成常数，但你想做得更细，可以在灵敏度分析阶段再加回随库容变化的水头模型，对比两次结果差异。这样既保证了复现的忠实度，又做出一小点自己的扩展。我建议在代码里把K做成可配置参数，默认值按论文给出的数据倒推，之后调整也方便。
水量平衡方程在Matlab里就是一个for循环逐时段追加约束：
% YALMIP约束写法 Constraints = []; for t = 1:T Constraints = [Constraints, V(t+1) == V(t) + (I(t) - Q_h(t) - S_spill(t)) * dt]; Constraints = [Constraints, V_min <= V(t+1) <= V_max]; Constraints = [Constraints, 0 <= Q_h(t) <= Q_max]; Constraints = [Constraints, P_h(t) == K * Q_h(t)]; Constraints = [Constraints, P_h_min <= P_h(t) <= P_h_max]; end Constraints = [Constraints, V(1) == V_start, V(T+1) == V_start];用YALMIP写这种双端不等式非常直观，可读性好，后面如果换求解器也不需要改约束声明。
2.3 从确定性调度升级到多场景调度
EI论文里一旦讨论风电不确定性，常用做法是构建多个风速场景，对每个场景求解同一调度问题，以期望成本最小化为目标。Matlab实现的步骤是：先生成大量风速样本，再做场景缩减。拉丁超立方采样配合kmeans聚类是主流做法，代码量很小：
n_scen = 500; wind_speed_all = wblrnd(8, 2, [T, n_scen]); idx_k = kmeans(wind_speed_all', 5); % 聚成5个典型场景 scen_prob = accumarray(idx_k, 1) / n_scen; % 场景概率 for s = 1:5 avg_speed = mean(wind_speed_all(:, idx_k == s), 2); Pw_avail(:, s) = power_curve(avg_speed); end多场景模型的变量维度会多出一维：每个时段、每个场景都要有对应的风电出力和水电出力决策，而库容决策可以按场景独立，也可以共享，取决于论文用的是两阶段随机规划还是多阶段模型。这个选择直接决定YALMIP的变量声明方式，建议一开始就把维度写清楚，不要写一半再改。
3. 联合优化模型矩阵化：为什么优选MILP和YALMIP/CPLEX
3.1 线性化带来的收益：全局最优和求解效率
复现这类模型时，很多人纠结要不要用非线性规划。我的经验是：凡是可以线性化的，就线性化。原因有两点。第一，MILP求解器能保证全局最优，而非线性规划很容易陷入局部最优，拿一个局部最优解去对比论文结果，完全没有说服力。第二，MILP求解器经过了几十年的工程优化，处理几千个变量、几万个约束非常成熟，求解速度和稳定性都比通用非线性算法好得多。
线性化的对象通常是水电出力方程中水头和流量的乘积项。可以按水头区间做分段线性近似，也可以用big-M法把逻辑条件变成线性不等式。做完这一步，模型整体就是线性约束加整数变量的标准MILP，可以直接交给商业求解器。
3.2 YALMIP建模代码骨架
YALMIP是Matlab生态下最顺手的优化建模工具箱，语法直观，核心价值在于不需要手写A矩阵。下面给一个最小可跑的确定性风-水联合调度骨架：
T = 24; dt = 1; % 参数 Pw_avail = wind_power_curve(wind_speed); % 风电可用功率序列 Load = load_profile(); % 负荷序列 I_t = inflow_profile(); % 来水序列 K = 0.85; % 水电出力系数 % 变量 P_w = sdpvar(T,1); % 风电出力 Q_h = sdpvar(T,1); % 发电流量 V = sdpvar(T+1,1); % 库容 P_h = sdpvar(T,1); % 水电出力 P_buy = sdpvar(T,1); % 外购电力 P_cur = sdpvar(T,1); % 弃风量 Constraints = []; Constraints = [Constraints, P_w + P_h + P_buy == Load]; % 功率平衡 Constraints = [Constraints, P_w + P_cur == Pw_avail]; % 风电实际出力+弃风=可用 Constraints = [Constraints, 0 <= P_w <= 50, 0 <= P_cur <= 50]; for t = 1:T Constraints = [Constraints, V(t+1) == V(t) + (I_t(t) - Q_h(t)) * dt]; Constraints = [Constraints, V_min <= V(t+1) <= V_max]; Constraints = [Constraints, Q_min <= Q_h(t) <= Q_max]; Constraints = [Constraints, P_h(t) == K * Q_h(t)]; Constraints = [Constraints, 0 <= P_h(t) <= 40]; Constraints = [Constraints, 0 <= P_buy(t) <= 100]; end Constraints = [Constraints, V(1) == V_start, V(T+1) == V_start]; Objective = sum(Price_buy .* P_buy) + 200 * sum(P_cur); ops = sdpsettings('solver', 'cplex', 'verbose', 1); optimize(Constraints, Objective, ops);注意这里弃风量P_cur是通过等式约束P_w + P_cur == Pw_avail推导出来的，目标函数用惩罚项抑制弃风。如果负荷不足或者外购电价便宜，求解器可能选择弃风而不是硬买高价电。这种"有取舍"的结果，正是分析联合调度价值的好素材。
3.3 求解器怎么选：一个对比表
| 求解器 | 类型 | 适用场景 | 许可证情况 |
|---|---|---|---|
| CPLEX | 商业MILP/LP | 大规模调度，科研标准配置 | 学术免费，需注册 |
| Gurobi | 商业MILP/LP | 与CPLEX同级，性能相近 | 学术免费 |
| CBC | 开源MILP | 小规模验证、避免商业依赖 | 完全免费 |
| linprog/intlinprog | Matlab自带 | 中小规模，无需额外安装 | 随Matlab |
| fmincon | 连续非线性 | 小规模非线性模型 | 随Matlab |
复现EI论文建议直接用CPLEX或Gurobi。论文里的算例规模通常上百个变量、上千条约束，CBC偶尔会出现数值问题，fmincon又处理不了整数变量。没有商业求解器的读者，先用intlinprog把逻辑跑通，再换YALMIP加CPLEX看性能差距，这样也能走通。
4. 跑通之后别急着收工：结果校验和灵敏度分析
4.1 设计三组基准算例：纯风、纯水、联合
模型第一次跑通之后，不要直接看结果就宣布成功。先设计对照组，否则你根本不知道结果合不合理。我通常会跑三组算例：纯风电加外购电力、纯水电加外购电力、风-水联合调度。三组用同一份负荷曲线和同一组参数，对比才有意义。
联合调度的目标函数值应该优于或至少不差于前两者——弃风量更小、购电成本更低，或者两者兼得。如果联合调度的目标函数值比纯风电还差，模型里大概率有约束写错了，常见的是水电调节能力被某个错误的参数限制住，比如库容上下限设得太窄、发电流量上限设得太小，导致水电完全没法配合风电出力。
4.2 从结果逆向验证模型正确性
拿到优化结果后，我会做四件事来验证。第一，检查各时段功率平衡约束是否满足，虽然优化过程已经强制满足，但数据导出时要复核一遍；第二，画出库容曲线，看它是否在上下限之间平滑变化，有没有不自然的跳变；第三，反向计算弃风比例，和论文报告的数值做对比；第四，把目标函数的成本分项拆出来看，检查每个分项的数量级是否合理。
有一个真实教训：有次我复现时，水电出力连续几个时段顶在上限，库容却还在上限附近。回头查发现是水量平衡方程里漏掉了蒸发项，论文正文里其实留了一行小字说明，只是没有放进数学模型里。这种问题只能靠逐条核对物理关系来找出来，哪怕求解器给的是可行解，依然可能违反实际物理逻辑。
4.3 灵敏度分析怎么设计，能让讨论部分直接有素材
灵敏度分析是复现工作的增值项，也是把"抄论文"变成"理解论文"的关键一步。常见做法有两个方向。一是对参数做变化：比如把弃风惩罚系数从50元/MWh改到200元/MWh，观察弃风量和购电成本的变化，画出关系曲线；二是对结构做变化：比如对比固定水头模型和变水头模型下的调度结果，量化简化假设带来的误差范围。
这些分析做完，不仅验证了模型行为符合预期，也让复现不再是简单照搬论文。如果后续要写自己的报告或文章，这些图表可以直接用作讨论部分的素材。我个人的习惯是每次只改一个参数，保存一份结果，最后汇总成一张大表。工作量不大，但结论会非常有说服力。
5. 复现过程中最容易翻车的五个坑
5.1 论文没给参数，怎么"科学地猜"
EI论文普遍存在参数不全的问题。来水序列、负荷曲线、风速数据这些算例参数，经常只画在图上而不放进表格里。我的处理方式是三步：先看论文算例描述里有没有标注数据来源，很多会注明来自某个标准测试系统或某地区实际数据；没有的话，参考同领域经典文献设定数量级一致的数据；最后在代码注释里标明"该参数为基于常见算例的合理估计"。
实测下来，数据大致符合即可，关键是量纲和数据形态要对。负荷要有早晚高峰，来水要有季节性趋势，风速要有波动特征。如果曲线形态都不对，后面所有结果分析都会失去意义。
5.2 水头常数的诱惑：别为了线性把物理特征丢掉
把水头简化成常数确实能让模型变成纯线性、求解飞快，但代价是放弃了水库调节的核心物理逻辑。水头会随库容变化而变化，库容低了，发电流量再大也可能发不出额定功率。如果论文明确用了常数水头，那按论文来即可；如果论文模型没有写清水头处理方式，建议至少做一次变水头对比，看看简化假设到底引进了多大误差。就我的经验来说，这一步往往比调求解器参数更能体现复现者的水平。
5.3 big-M的M取值：小模型耗时长的大半原因
只要模型里用了逻辑约束或者分段线性化，就一定会碰到big-M。M取得太大，数值稳定性差、求解时间暴涨；M取得太小，可能把可行域切掉一块，得到次优解甚至无解。经验做法是给每个约束单独设置与物理量纲匹配的M值，比如功率相关的M取装机容量的两倍，而不是全文共用一个10000。在YALMIP里尽量给变量设置合理边界，也能帮助求解器预处理，效果比加大M值好得多。
5.4 日前调度结果和实时出力对不上
风电预测本来就有误差，日前调度给的是基于预测风速的计划，真实风速出来以后还要做实时修正。复现时如果发现计划出力和实际可用出力差异很大，不要急着说模型错了。正确做法是明确区分两层：日前层做优化决策，日内层做偏差调整。如果论文只讨论日前调度，那复现到计划层就够了，不要强行混入实时修正逻辑，否则结果会不伦不类。
5.5 求解器报Infeasible：先检查这四处
模型不可行时，95%的问题出在四个地方：库容初值和末值设置矛盾、功率平衡约束少了一项、水电出力上下限和流量上下限换算不一致、某个big-M参数过小而切掉了可行域。我的排查顺序是：先去掉目标函数只看约束可行性，再逐个放宽功率平衡、库容边界、水量平衡，看哪一步把模型从不可行变成可行。这套流程下来，几乎都能在半个小时以内定位根因，比盲目改参数高效得多。
6. 从论文到代码的整体复盘与操作心得
整套流程走完，我最深的体会是：复现的价值不在跑通代码那一刻，而在跑通之后你不得不去理解论文里每一处"为什么"。比如为什么目标函数里弃风惩罚要设置成某个量级，为什么水库调度期末要回落到初库容，为什么风速场景要缩减成五个而不是二十个——这些问题在纯读论文时很容易滑过去，写代码时却绕不开。
如果你也要复现一篇EI论文，有件事值得尽早动手：建一个变量字典，把论文里的隐含假设全部记录到代码注释里，同时准备好三组基准算例，最后把关键参数做一轮灵敏度分析。这四条做完，复现的成功率和理解深度都会有明显提升。
最后再分享一个小技巧：EI论文的数学符号经常不统一，同一个符号在不同章节可能含义不同。遇到这种情况不要靠猜，优先以前文最早的定义为准，并且把定义写进代码注释。完整代码和测试数据我已经整理打包，需要的朋友可以到我的资源区获取，欢迎在评论区交流建模和求解过程中遇到的具体问题。
