企业官网建设流程全解析

前一阵子接到“风-水电联合优化运行”这个EI复现任务，原本以为最费时间的会是代码实现，真正跑起来以后才发现，从论文公式到一份能稳定复现出结果的Matlab工程，中间隔着大量没有写在纸上的建模细节和数据假设。这篇文章就围绕我完整走通的一条复现链路展开，包括目标函数怎么定、水电出力特性怎么简化、风电不确定性怎么进模型、Matlab代码怎么搭、最后怎么验证“复现成功”，给打算做EI复现或做电力系统优化调度的同学一个可以直接参考的模板。

这套模型要解决的核心问题很直观：风电出力有天然波动，火电机组调节慢且成本高，水电站依托水库具备快速响应能力。把三者放进同一个优化框架，在满足负荷需求和水力约束的前提下，让系统总运行成本最低、同时尽量少弃风少弃水，这就是联合优化运行的全部内核。我下面的所有内容都会围绕这条主线展开，不管是复现论文、准备毕业设计，还是想评估水电对新能源的调节价值，应该都能从中找到可以直接落地的思路。

1. 复现前我先把这三个问题想清楚了

写代码之前，我花了大概两天时间反复读原文，最后总结出三个必须先回答的问题。这三个问题如果没想清楚，后面写多少代码都是白费。

1.1 先搞清楚目标函数到底在优化什么

很多复现失败不是代码问题，而是连“优化什么”都没定义清楚。风-水电联合调度的目标函数在同类论文里有好几种写法，有的最小化系统总运行成本，有的最大化新能源消纳量，还有的做成双目标去权衡经济性与碳排放。我复现的这篇以经济调度为基调，目标函数取的是系统总运行成本最小化，具体包含三项：

• 火电机组煤耗成本，通常写成二次函数 a·P² + b·P + c；
• 弃风惩罚，即预测风电没有被消纳的部分乘以一个惩罚单价；
• 弃水惩罚，即水电站可以发但被迫放弃的水量乘以惩罚单价。

对应的目标函数表达式就是：

min Σ时段 Σ机组 [ a_i · P_Ti² + b_i · P_Ti + c_i ] + λ_w · Σ P_弃风 + λ_h · Σ Q_弃水

这里最容易被忽略的是惩罚系数λ_w和λ_h的取值。它们不是越大越好，也不是越小越好。λ_w太小，模型会为了省煤耗成本主动弃风，违背“尽量消纳新能源”的调度倾向；λ_w太大，则会强迫系统在所有场景下全额消纳风电，遇上极端风况可能直接无解，或者把水电调度搞得过度激进。我的经验是先从火电边际成本的0.6~1.2倍取初值，然后跑几组敏感性测试，观察弃风量和总成本的变化拐点。如果某个λ取值让弃风量从10%骤降到接近0，同时总成本跳升一大截，那说明这个点附近就是合理的惩罚量级。

1.2 约束条件的“硬”与“软”要分开处理

优化模型里不是所有约束都处于同等地位。我的做法是把约束分成两类处理：

• 硬约束：功率平衡、火电上下限、火电爬坡速率、水电站出力上下限、水库库容上下限。这些是物理上必须满足的，违反任何一条，求解器给出的解都不能用。
• 软约束：弃风、弃水。它们在目标函数里以惩罚项出现，等于给“尽量消纳”一个优先级，而不是必须做到。这样做的好处是模型在极端场景下依然存在可行解，只是成本高一些；如果非把弃风弃水都写成等式约束“必须为零”，那随机场景耦合下经常得到infeasible的结果。

功率平衡是等式约束，必须严格满足，这一点在代码里我会用容差判断去校验，后面专门有一节讲。火电爬坡约束和水电调节速率约束是一对矛盾：火电爬坡慢、水电爬坡快，如果把火电爬坡设得太紧，水电就不得不频繁调整出力去追负荷和风电波动。复现时我建议先给火电一个较宽松的爬坡范围（比如每分钟1%~3%额定出力），跑通后再收紧看影响，这样能帮你判断论文里“水电发挥了调节作用”到底是模型内生出来的结论，还是仅仅因为火电被约束绑死了。

1.3 论文没给的数据，怎么补全而不失真

EI论文最让人头疼的就是关键数据经常不给你完整表格，机组参数、负荷曲线、来水过程全在图片里，只能靠肉眼读数。我当时为了凑齐一套能跑的系统数据，用了三步处理：

• 系统规模：采用3台火电、2座梯级水电站、1个风电场的简化系统。负荷曲线取典型日负荷，峰值大致在600MW~800MW这个区间。
• 水电参数补全：每座水电站设置装机容量、库容上下限、最大发电流量、综合出力系数K。这些参数不完全来自论文，部分参照同类型水电站的公开运行数据做了合理预估。
• 敏感性检验：对每个补全的参数做±20%的扰动，重新求解，观察最终调度结论是否翻轉。如果参数变化20%结论就变了，说明模型对那个参数过于敏感，要么参数估计有问题，要么模型结构需要调整。

这一步花的时间不少，但非常值得。数据补全的过程本身就是对模型理解深度的检验——你不知道哪些参数重要，就说明还没把论文吃透。

2. 风-水电联合调度模型，数学上到底堆了几层

模型本身并不复杂，但每一层都有需要注意的细节。我按从整体到局部的顺序拆一下。

2.1 功率平衡与目标函数合在一起看

任意时段t，系统必须满足：

Σ 火电出力 + Σ 水电出力 + 实际消纳风电 = 负荷

这里“实际消纳风电”是一个变量，它等于预测风电减去弃风：

P_w_use(t) = P_w_forecast(t) - P_curtail(t)

所以功率平衡可以写成：

Σ_i P_Ti(t) + Σ_j P_Hj(t) + P_w_use(t) = P_load(t)

火电煤耗成本按二次函数求和，水电和风电的变动成本在目标函数里可以视为零。这种结构下，优化器的逻辑就变成：只要水电和风电还能出力，就优先让它们顶上去，火电只填补剩下的净负荷。

2.2 水电站出力计算：水头-流量非线性的处理

水电站的出力表达式是：

P_H = g · η · ρ · Q · H

其中Q是发电流量，H是净水头。麻烦在于H本身又是库容的函数，水库水位随调度过程变化，导致这个约束本质上是非凸的，直接放进MILP框架很难处理。

多数EI论文并不会展开这一层非线性，而是采用简化表达。我复现时用的是固定水头近似：调度周期只有24小时，如果水库不是特别小，库容变化引起的水头变化通常不超过几个百分点。把H近似成常数后，出力表达式就变成：

P_H = K · Q_flow

K为综合出力系数，单位是MW/(m³/s)。这样做有两个直接好处，一是约束变成线性，求解稳定性大幅提升；二是物理意义依然清晰，发电流量直接决定出力，水量平衡依然把库容、入流、出流绑在一起。

如果对精度要求更高，可以改成分段线性化：把库容区间分成三段，每一段取不同的K系数，再用0-1变量选择当前时段处于哪个区间。本质上就是在用混合整数线性规划拟合一族非线性曲线。我测试过，分段线性化对最终成本的影响大约在1%以内，但求解时间会成倍增加，所以复现阶段取固定水头近似足矣。

2.3 风电不确定性的场景化表达

确定性模型直接用预测曲线，等于假设预测完全准确。这个假设在真实调度里站不住脚，也是很多论文复现后结果对不上原图的原因——原文明明在讨论不确定性，你拿确定性模型去对结论，当然对不上。

我的方案是两层处理：

• 场景生成：对预测风电曲线叠加正态分布误差，误差标准差取预测值的10%~20%。用蒙特卡洛方法生成200个随机场景。
• 场景削减：用同步回代法（backward reduction）把200个场景削减到10个有代表性的场景，每个场景附带一个概率权重π_s。

随机规划下目标函数变成各个场景成本的加权期望：

min Σ_s π_s · cost_s

关键是非预期性约束：所有场景共享第一阶段的火电开机决策，或者至少保证同一时段、同样信息状态下的决策一致。水电出力可以随场景调整，因为水电的调节价值恰恰体现在“知道风电不准之后还能补救”。如果你把所有阶段的决策都按场景独立优化，结果会过于乐观，也违背了“决策必须基于当前可用信息”的调度逻辑。

3. Matlab代码架构：从脚本堆砌到能跑的工程

这部分是我花时间最多的地方，也是复现价值的直接体现。

3.1 建模环境选型：Yalmip + Cplex

我没有从零手写求解算法，那不现实也没必要。我的技术栈是Yalmip做建模层，Cplex或Gurobi做求解层。Yalmip在Matlab下声明优化变量、写约束、调用求解器都非常顺手，比直接调Cplex的C API舒服太多。

如果不用商业求解器，Matlab自带的intlinprog也能跑，但处理多场景随机规划时会明显变慢。10个场景、24个时段、几十个变量的规模下，Cplex通常在几十秒内收敛，intlinprog可能要跑到数分钟甚至更久。所以有条件的尽量上Cplex或Gurobi，学生版授权也不难申请。

3.2 核心代码段与逐段解释

决策变量定义部分，我用Yalmip写成这样：

%% 系统参数 T = 24; % 调度时段数 nT = 3; % 火电机组台数 nH = 2; % 水电站数量 nW = 1; % 风电场数量 tau = 1; % 时段长度，单位小时 %% 决策变量 P_T = sdpvar(nT, T, 'full'); % 火电出力 P_H = sdpvar(nH, T, 'full'); % 水电出力 Q_flow = sdpvar(nH, T, 'full'); % 发电流量 Q_spill = sdpvar(nH, T, 'full'); % 弃水流量 V = sdpvar(nH, T+1, 'full'); % 库容，注意多出一列 P_w_use = sdpvar(nW, T, 'full'); % 实际消纳风电 P_curtail = sdpvar(nW, T, 'full'); % 弃风功率

这里有一个我一开始没注意的细节：库容V定义成T+1列，因为水量平衡要从0时段推进到T时段，最后需要T+1个状态点。

约束构建部分，功率平衡和水量平衡是最核心的两块：

Constraints = []; %% 功率平衡约束：任意时段，火电+水电+风电消纳=负荷 for t = 1:T Constraints = [Constraints, sum(P_T(:,t)) + sum(P_H(:,t)) + sum(P_w_use(:,t)) == P_load(t)]; end %% 水量平衡约束：V(t+1) = V(t) + (入流 - 发电流量 - 弃水) * tau for h = 1:nH for t = 1:T Constraints = [Constraints, V(h,t+1) == V(h,t) + (Q_in(h,t) - Q_flow(h,t) - Q_spill(h,t)) * tau]; end end %% 风电消纳关系：消纳量 = 预测值 - 弃风 for t = 1:T Constraints = [Constraints, P_w_use(:,t) == P_w_forecast(:,t) - P_curtail(:,t)]; end

目标函数部分，火电成本加上惩罚项：

Cost_thermal = 0; for i = 1:nT Cost_thermal = Cost_thermal + a(i) * P_T(i,:) * P_T(i,:)' + b(i) * sum(P_T(i,:)) + c(i) * T; end Cost_curtail_wind = lambda_w * sum(P_curtail(:)); Cost_spill = lambda_h * sum(Q_spill(:)); Objective = Cost_thermal + Cost_curtail_wind + Cost_spill;

求解设置：

ops = sdpsettings('solver', 'cplex', 'verbose', 2, 'savesolveroutput', 1); ops.cplex.mip.tolerances.mipgap = 0.001; % 设置1e-3相对gap result = optimize(Constraints, Objective, ops);

这里有一个经验值：mipgap设到1e-3足够了，再往下压除了让求解器多跑几分钟，对结果的实际影响很小。如果你的目标函数里存在非常接近的数量级差异，可以适当放宽到5e-3，换来求解时间的大幅下降。

3.3 求解效率优化：从一小时压到五分钟

第一次把200个场景全部塞进去跑，求解器直接卡住不动。后面我做了几件事，把求解时间从一小时量级压到了五分钟左右。

第一，尽量矩阵化，减少for循环。Yalmip里变量本身支持向量化求和，像sum(P_T, 2)这样一次算完所有时段的某一台机组出力，能少写循环就能少给求解器增加不必要的内存占用。

第二，做量纲统一。所有机组出力和负荷统一用MW，所有水流量统一用万m³/h。如果火电煤耗系数是10⁻⁴量级，水库库容却是10⁷量级，目标函数和约束里的数值会差12个数量级，这对Cplex/Gurobi的数值稳定性非常不友好，会出现ill-conditioned的警告，甚至莫名其妙的不可行解，下一章节细说。

第三，把二次煤耗成本分段线性化。如果直接用二次函数，模型变成MIQP，虽然也能解，但速度比纯MILP慢很多。我试验过后直接在论文中采用了分段线性化处理，比如把出力区间切成四段，每段线性化成本斜率，整体变成纯MILP。这样既能保持精度，又能让求解速度大幅提升。

优化前后的对比大概是这样：

4. 复现结果到底算不算“复现成功”

复现不是“跑出一个好看的图”就结束了。我给自己定了三条判断标准。

4.1 逐时段的物理约束校验

输出方案以后，第一件事不是看图，而是写一段校验脚本，把每个时段的功率平衡残差、水量平衡残差、库容上下限、爬坡约束逐个检查一遍。功率平衡的残差应该小于1e-6 MW量级，水量平衡残差应该小于1e-6万m³量级。如果残差偏大，基本可以判定建模或求解环节有问题。

我习惯的做法是：

res_power = sum(P_T, 1) + sum(P_H, 1) + P_w_use - P_load; max_abs_res = max(abs(res_power)); fprintf('最大功率平衡残差: %.2e\n', max_abs_res);

如果这里出现明显残差，先回头检查约束有没有漏写，再检查求解器设置是不是把容差放得太宽了。Cplex默认容差通常够用，但如果你手动改过eprhs或mipgap，要留意它们对解的影响。

4.2 与论文图表对比的量化方法

论文里的曲线是图像，没有原始数据，所以对比只能做数值化处理。我的对比方式分三层：

• 目标函数值对比：同规模、同时段下，论文的总成本和我的目标函数值应在正负10%以内。超过这个范围，说明某些参数或者约束条件设置存在较大偏差。
• 曲线形状对比：计算火电、水电、风电消纳曲线与论文图的RMSE和Pearson相关系数。相关系数在0.9以上，基本可以认为趋势吻合。
• 典型模式验证：看水电出力是否在负荷高峰时段显著增加、在低谷时段减小，火电则跟随净负荷趋势变化。如果水电出力曲线是平的，说明水库的调节潜力没有被模型发挥出来，约束设置或者目标函数大概率有问题。

4.3 压力测试：边界条件下的稳定性

最后一步，也是我认为最容易被复现者跳过的一步：压力测试。调整输入场景，观察模型是否稳定给出合理的调度方案。

• 极端风电场景：把所有时段的风电预测设为零，模型应该自动增加火电和水电出力度过“无风期”；把所有时段设为额定出力，模型应该能靠水库调节能力消纳多余风电，实在消纳不了时才弃风。
• 枯水场景：把来水大幅降低，水电出力受限，火电应当顶上，库容约束不出现越界。
• 丰水场景：来水很足，模型应当优先让水电多发，同时警惕弃水惩罚设置是否合理。

如果模型在边界场景下经常报infeasible或者给出明显不合常理的调度结果，问题往往不出在求解器，而出在约束条件本身存在互相矛盾的地方。

5. 我这一趟踩过的坑，按重现顺序梳理

踩坑是复现过程中最有价值的部分，写出来帮大家省时间。

5.1 库容时间索引的错位导致水电站“凭空抽水”

第一次跑通模型后，我检查水量平衡，发现第一个时段的库容变化和出入流对不上，库容凭空跳了一大截。排查了半天，最后发现是索引错位：我把水量平衡写成了V(h,t) == V(h,t-1) + ...，然后在t=1时引用了V(h,0)，Matlab里索引从1开始，导致这个约束被默认忽略或者读取了错误的初始值。正确写法是V(h,t+1) == V(h,t) + ...，让时间索引从初始库容逐步向前推进。

排查这个问题的办法就是在校验脚本里逐时段打印水量平衡各项，不要只打印残差总量。残差为零可能是约束没生效，打印每一项才能暴露索引问题。

5.2 场景削减过猛导致无可行解

最初我把200个场景削减到5个，结果求解器报了infeasible。原因在于场景数量太少时，爬坡约束和功率平衡约束在少数几个场景里互相牵扯，解空间就被锁死了。200个场景虽然解大，但冗余场景提供了更多的缓冲余地。削减到5个后，每个场景都必须被满足，约束过紧。

解决方法是：削减后的场景数不要低于10个，并且削减完成后对权重做归一化，重新检查场景矩阵的条件数。如果条件数异常，说明保留的场景之间相关性太高，代表性不足。

5.3 求解器的数值病态：量纲不统一害死人

有一段时间Cplex老是警告ill-conditioned，解出来的结果时好时坏，有时候弃风量突然变成负数。后来我把所有量纲拉出来对比，发现火电煤耗系数在10⁻⁴量级，水库库容在10⁶量级，惩罚系数λ又取了比较大的数，整个约束矩阵的条件数已经差到十分恶劣的程度。

解决办法就是前面提到的量纲归一化：把出力统一到MW，库容统一到万m³，所有惩罚系数和煤耗成本保持在相近数量级。改造之后，求解器警告彻底消失，结果也稳定了。这事给我的教训是：代码里优先写好量纲注释，一劳永逸。

写在最后的体会

我整个复现过程做完的体会是：好的复现不是把公式抄一遍，而是要在“模型保真度”和“可求解性”之间找平衡。固定水头近似、煤耗分段线性化、场景削减到10个，这些看上去都是让模型“变简单”的操作，却能让求解稳定性和结果可用性大幅上升，反而不是坏事。另一个比较深的感受是数据准备和结果验证几乎占掉了六成以上时间，真正写模型框架反而是快的那部分。建议后来者从确定性模型先跑通，再加场景和不确定性，一步到位出问题的可能性会高很多。

如果你之后想把这套模型往更实用的方向扩展，可以考虑把调度周期从24小时拉长到几天，把末库容从固定值改成优化边界，让水库在更长的时间尺度上做跨日调节。改动量不大，但结论的丰富程度会明显不一样。