最近刚把一个EI论文复现项目收尾主题是风-水电联合优化运行的Matlab代码实现。这个方向在可再生能源调度这块太常见了论文里通常构建一个24时段的经济调度模型把风电场和水电站放在同一个优化框架里求解用水库的蓄放水去平抑风电波动最终得到每台机组的出力计划和水库运行策略。如果你也接到类似的复现任务或者想把论文里的数学模型变成真正能跑出结果的代码这篇文章应该能帮你少走不少弯路。拿到这种复现任务很多人第一反应是找“论文对应的源代码”但多数EI论文是不会公开代码的。你能拿到的通常只有数学模型、参数表和几张结果图。所以复现的本质是根据论文中的模型定义自己从头搭一套可运行的代码再用论文里的结果做对标。这个过程有几个核心节点比较考验人模型翻译、变量编码、约束矩阵搭建、求解器选择。下面我按自己做项目时的顺序把每个环节都拆开讲一遍。1. 为什么风、水要放在同一个优化框架里风电出力的主要特点是不可控和随机性。风速的波动直接反映在出力曲线上小时级变化非常明显而且风电场本身没有“存储”能力风来了电就得发出去负荷侧用不完就只能弃掉。你可以把风电想象成一场暴雨雨量大小完全看天你只能被动接收没法让它按你的需求定时定量下。水电则完全是另一种性格。水电站可以通过调节闸门和发电流量在几分钟内改变出力大小。更重要的是水库本身就是天然的能量储存装置——水量攒在那里相当于电量也存在那里。大风时段让水电机组少发、把水蓄起来风小的时候再把水放出来发电本质上就是一套大容量储能系统。联合优化要解决的核心问题就是在满足负荷需求的前提下把水电的出力轨迹和风电的随机波动匹配起来使系统总运行成本最低或者弃风电量最小。举个具体例子某时段风电场可用出力是180MW但系统负荷只有150MW水电机组最小出力20MW那么功率平衡下风电只能发130MW必须弃掉50MW。如果没有优化模型拍板这个决策可能很随意但有了约束和目标函数系统能自动算出每个时段该发多少、该弃多少、水库该蓄还是该放。在EI论文里这类模型按目标函数可以分成三类复现前必须看清楚你手里这篇属于哪一种目标类型数学思路典型场景经济调度型最小化水电/火电运行成本与弃风惩罚之和电网公司或发电集团的运行决策消纳最大化型最小化弃风电量或最大化风电容纳量高比例风电地区消纳评估综合效益最大化型计及上网电价后最大化总收益电力市场环境下的发电商决策不同目标对应不同的约束组合和量纲。如果拿到一个经济调度论文却照着消纳论文的模型去改目标函数量纲都对不上后面肯定白折腾。我见过有人在这个环节卡了一整天最后才发现论文的目标函数里根本没有弃风惩罚项。2. 先把EI论文里的模型翻译成可编程语言进入到具体建模之前我建议先在纸上把所有变量和约束列清楚这一步做好了后面写矩阵会非常顺畅。2.1 决策变量拆解一个典型的风-水联合日调度模型决策变量通常包括下面五个序列变量含义维度单位P_w(t)风电场实际出力1×24MWP_h(t)水电站出力1×24MWQ(t)发电流量1×24万m³/hV(t)水库蓄水量1×25万m³P_cur(t)弃风功率1×24MW注意V的维度比时段数多1因为水量动态方程会牵扯到调度开始和结束两个时刻的库容状态。2.2 目标函数与约束条件列表以最基本的“经济调度型”为例目标函数取min Σ_t [ c_h * P_h(t) λ * P_cur(t) ]其中 c_h 是水电运行成本系数λ 是弃风惩罚系数风电边际成本视为0。λ通常取得比水电成本大很多用来表达“尽量别弃风”的倾向。约束条件我用一张表总结每一条都要对应到Matlab矩阵里的一行或几行编号约束名称数学表达式物理含义1系统功率平衡P_w(t) P_h(t) P_load(t)每个时刻发电与负荷相等2风电出力与弃风关系P_w(t) P_cur(t) P_w_avail(t)用不完的风电才叫弃风3风电出力上下限0 ≤ P_w(t) ≤ P_w_avail(t)实际出力不超过可用功率4水电出力上下限P_h_min ≤ P_h(t) ≤ P_h_max水电机组物理出力范围5水量-出力转换P_h(t) k * Q(t)发电流量越大出力越大6库容水量平衡V(t1) V(t) I_t - Q(t)水库蓄水变量动态关系7库容上下限V_min ≤ V(t) ≤ V_max水库安全运行范围8调度期末库容V(25) V0模拟日不透支水资源其中第5条里的k是把水头、重力加速度、机组效率综合在一起折算出的常数单位是MW/(万m³/h)。在简化模型里通常按平均水头取一个近似值比如2.3左右。严格来说水头变化会让k变成非线性的但很多EI论文为了可解性都做了线性化处理复现时也跟着论文走就行。还有一点要特别注意第6条水量平衡是全局耦合约束它通过V把24个时段串联在一起。和风电、水电出力上下限这种“单时段独立约束”相比这类跨时段约束最容易写错。3. Matlab实现从变量编号到求解器的完整流程模型写清楚之后就到了最核心的代码环节。这里没有任何取巧的空间所有约束必须变成linprog认识的矩阵形式。3.1 整体思路先定变量再做矩阵linprog只认识一个列向量x以及一组等式约束 Aeq * x beq 和一组上下界。所以第一步就是把所有时段变量打平成一个大向量并记住每个子序列在x中的起始位置。我的做法是T 24; p_w_idx 1:T; % 风电出力 P_w p_h_idx T1:2*T; % 水电出力 P_h p_cur_idx 2*T1:3*T; % 弃风功率 P_cur q_idx 3*T1:4*T; % 发电流量 Q v_idx 4*T1:5*T; % 库容 V(2)~V(25) nvars 5*T;V(1)是初始库容这里作为常数处理不放到x里。这样x的长度就是5×24120。后面的Aeq每一行代表一条约束每一列对应一个变量千万别把行列号搞反。3.2 数据准备数据是整个复现的基础。负荷曲线、风电可用功率曲线、来水序列、水电站参数都需要在进入优化前读进来。下面是我构造的一组示例数据先跑通结构再替换成你自己论文里的数据% 24小时负荷MW P_load [180 175 170 168 170 180 190 200 205 210 205 200 ... 195 190 185 180 185 195 200 205 200 195 190 185]; % 24小时风电可用功率MW P_w_avail [100 95 90 80 75 70 80 95 110 115 105 90 ... 80 70 65 60 70 85 100 110 105 90 80 70]; % 24小时来水万m3/h I [42 42 43 44 45 46 46 47 47 46 45 44 44 43 43 44 ... 45 46 47 46 45 44 43 42]; % 水电站参数 P_h_max 180; % 最大出力 MW P_h_min 20; % 最小出力 MW k 2.31; % 水量-出力转换系数 MW/(万m3/h) Q_max P_h_max / k; % 最大流量 万m3/h Q_min P_h_min / k; % 最小流量 万m3/h V_min 200; % 最小库容 万m3 V_max 1200; % 最大库容 万m3 V0 600; % 初始库容 万m3 V_end V0; % 期末库容要求 % 成本系数 c_h 50; % 水电运行成本 元/MWh lambda 500; % 弃风惩罚 元/MWh这里有个容易被忽视的问题发电流量上下限其实是由水电出力上下限和k换算出来的如果论文给的是水轮机最大过机流量记得换算成一致的单位。单位混用的坑后面单独讲。3.3 构建Aeq和beq目标函数系数向量f是最简单的f zeros(nvars, 1); f(p_h_idx) c_h; f(p_cur_idx) lambda;然后是关键部分把8条约束转成Aeq和beq。我习惯用稀疏矩阵来存虽然120变量规模不大但如果你以后扩展到8760时段稀疏矩阵和稠密矩阵的性能差距会非常明显。% 等式约束总数功率平衡T 弃风定义T 水量平衡T 出力-流量关系T 期末库容1 n_eq 4*T 1; Aeq sparse([], [], [], n_eq, nvars); beq zeros(n_eq, 1); row 0; % 约束1功率平衡 P_w(t) P_h(t) P_load(t) for t 1:T row row 1; Aeq(row, p_w_idx(t)) 1; Aeq(row, p_h_idx(t)) 1; beq(row) P_load(t); end % 约束2弃风定义 P_w(t) P_cur(t) P_w_avail(t) for t 1:T row row 1; Aeq(row, p_w_idx(t)) 1; Aeq(row, p_cur_idx(t)) 1; beq(row) P_w_avail(t); end % 约束3水量-出力转换 P_h(t) k * Q(t) for t 1:T row row 1; Aeq(row, p_h_idx(t)) 1; Aeq(row, q_idx(t)) -k; beq(row) 0; end % 约束4水量平衡 V(t1) V(t) I(t) - Q(t) for t 1:T row row 1; if t 1 % V(2) - Q(1) V(1) I(1) Aeq(row, v_idx(t)) 1; Aeq(row, q_idx(t)) 1; beq(row) V0 I(t); else % V(t1) - V(t) - Q(t) I(t) - V(t) 移项后 Aeq(row, v_idx(t)) 1; Aeq(row, v_idx(t-1)) -1; Aeq(row, q_idx(t)) 1; beq(row) I(t); end end % 约束5期末库容 V(25) V_end row row 1; Aeq(row, v_idx(T)) 1; beq(row) V_end;第4条水量平衡里t1时V(1)是已知常数所以要把它挪到等式右侧这点特别容易漏。我之前第一次写的时候就忘了V0这一项结果模型跑出来的库容轨迹整体偏移了600万m³。变量上下界放在lb和ub中比写成不等式约束更高效lb zeros(nvars, 1); ub zeros(nvars, 1); lb(p_w_idx) 0; ub(p_w_idx) P_w_avail; lb(p_h_idx) P_h_min; ub(p_h_idx) P_h_max; lb(p_cur_idx) 0; ub(p_cur_idx) P_w_avail; lb(q_idx) Q_min; ub(q_idx) Q_max; lb(v_idx) V_min; ub(v_idx) V_max;3.4 求解与结果提取到这里就可以调用linprog了。我习惯用dual-simplex算法因为这类带大量上下界的稀疏线性规划问题dual-simplex通常比interior-point更快而且迭代过程更透明options optimoptions(linprog, Algorithm, dual-simplex, Display, iter); [x, fval, exitflag] linprog(f, [], [], Aeq, beq, lb, ub, options); if exitflag 0 error(求解失败exitflag %d, exitflag); end % 结果的索引还原 P_w x(p_w_idx); P_h x(p_h_idx); P_cur x(p_cur_idx); Q x(q_idx); V [V0, x(v_idx)]; % 加上初始库容我通常会把P_w、P_h、P_cur、Q、V这五个序列画在同一张图上检查有没有明显违背物理直觉的地方比如库容突然跳变、弃风和风电同时为正这类情况。4. 复现EI论文时最容易踩到的四个坑这部分都是我自己真实踩过的坑挨个排过来希望你能绕开。4.1 单位混用一个乘以36就崩掉整个模型复现论文时数据表格里的单位常常五花八门发电流量用m³/s径流量用亿m³库容用万m³出力用MW。如果不先统一单位水量平衡方程直接乱套。典型错误是论文里给的最大发电流量是40m³/s你在代码里直接当成40万m³/h用。1m³/s 3600m³/h 0.36万m³/h也就是说40m³/s等于14.4万m³/h量级差了快三倍。第一次跑出库容轨迹冲顶之后我逐项打印了每个变量的量级才定位到是Q统一成万m³/h后应该除以2.78而不是直接抄。建议在数据准备区把所有输入统一成“MW、万m³、万m³/h、小时”这套单位并在代码开头加一行注释标明单位。跑完后把日均流量、库容变化量打印出来和输入数据验算一下能瞬间暴露这类问题。4.2 约束写成了不可行域linprog直接无解linprog返回exitflag-2说明找不到可行解。我遇到这种情况第一反应先检查功率平衡能不能满足如果某个时段P_load大于P_h_max加上P_w_avail那无论怎么调都缺功率模型一定无解。这时候只有两个方向降低负荷或者把水电装机上限调大。如果负荷没问题就去查等式约束的符号。水量平衡里V(t1) V(t) I(t) - Q(t)移项写进Aeq时每一项的正负号都要对。一个很实用的调试手段把Aeq和beq打印出来手算一个已知可行解代入比如假设所有Q都取最小值、库容保持不动看beq残差是不是0。这样能很快定位到具体哪一行错了。还有一种情况是期末库容约束太紧。日调度模型如果强制V(25)V0而来水不足以支撑全天负荷模型就会无解。这种时候先去掉期末约束跑一遍如果库容掉到了下限以下再回去调整来水数据或者放宽期末约束而不是硬调约束逻辑。4.3 求解器选不对非线性模型硬套linprog我见过不少复现项目卡在“为什么我的linprog会报错”上最后发现模型里出现了两个变量相乘的非线性项。典型的就是计及水头变化的出力表达式P_h η * ρ * g * Q * H这里面Q和H都是变量乘积项是二次的linprog根本处理不了。正确的做法有两种一是按照论文的线性化方法把水头当作常数或者分段常数把乘积项拆成线性项二是改用fmincon做非线性规划。但fmincon没有全局收敛保证容易陷入局部解。所以如果EI论文本身是线性规划复现时务必严格保持线性不要自己引入非线性环节。如果模型里有0-1变量比如水电机组启停状态那么linprog也要让位于intlinprog。这个内置函数可以处理混合整数线性规划对中小规模问题够用。真要上大规模或者更快的求解器再考虑YALMIP配合Gurobi或者CPLEX。4.4 结果和论文对不上不是代码错了是场景数据不同复现论文最让人沮丧的是自己的结果和论文结果对不上。但这里有个现实问题EI论文里的负荷曲线、风电预测曲线、来水序列往往属于“测试场景”论文可能只给了汇总数据完整序列根本不会公开。你用自己的数据天然不可能得到和论文一模一样的数字。我验证模型是否“复现成功”看的不是绝对数值而是以下几点目标函数值是否在合理范围弃风率变化趋势是否一致水电出力轨迹和库容变化量是否在边界内如果论文做了敏感性分析改变λ时弃风率的趋势是否一致。把这些对齐了模型框架基本就是正确的。有一个我常用的交叉验证方法构造一个论文明确提到的极端场景。比如论文说“当风电为0时全部负荷由水电承担”那我把P_w_avail全设为0看模型是不是让水电满出力。这种人工构造的边界测试能快速判定约束写没写对。5. 用三个指标验证联合优化是不是真的有效模型跑出来不能只看一眼目标函数值就完事。我会用三个指标来判断联合运行方案是不是真的比单独运行更好。5.1 弃风率弃风率是衡量风电消纳水平最直接的指标total_cur sum(P_cur); total_avail sum(P_w_avail); curtail_rate total_cur / total_avail * 100;单独运行风电时只要负荷小于可用风电就只能弃风弃风率完全由负荷曲线决定联合运行后水电在风小的时段弥补缺口在风大的时段压低出力甚至停发相当于给风电让出一条路。对比两个场景的弃风率差异非常直观。我一般会写个基准模型只跑风电约束不加入水电的水量动态做一个对照。5.2 水电电量占比与库容轨迹合理性水电日发电量可以算E_h sum(P_h); E_load sum(P_load); share_h E_h / E_load * 100;占比太高说明水电承担了主要供电任务这时更要警惕库容能不能支撑。我的习惯是把V序列画出来和V_min、V_max放在同一张图里检查有没有触顶或触底的时段。如果库容在某个时段压到下限说明水电机组那时已经没有调节能力了模型的“联合效果”要打个折扣。5.3 出力平滑度与负荷匹配联合运行另一个好处是让系统总出力更加平滑地跟踪负荷。可以看水电出力的小时级变化幅度smooth_h mean(abs(diff(P_h)));smooth_h越小说明水电出力变化越平缓机组机械磨损和调度压力都更小。但要注意这个指标不是越小越好水电的核心价值就是快速调节有时候剧烈的出力变化恰恰是为了补偿风电波动。综合来看我会给出一张类似这样的对比表虽然数据是演示用的但形式上很直观指标风电单独运行风-水联合运行总发电量MWh18402070弃风率%16.88.1负荷满足率%93.2100水电日均出力MW-96联合优化之后负荷满足率提升是最关键的收益它意味着系统在极端场景下不再缺电。6. 做完这个复现后我想再说的几件小事这个项目做完我对“EI复现”这件事有了更实在的理解。第一模型框架的通用性比代码本身更重要。一旦把变量编号、Aeq构建这套思路跑通往里面加火电、加储能、加抽水蓄能都只是增加新的变量和约束行核心逻辑不用动。第二如果要把规模做大比如做全年8760小时的调度内置linprog还是会吃力这时候推荐YALMIP加Gurobi的组合建模更直观求解速度快一个量级。但作为复现和教学用linprog完全够用而且不依赖外部许可证。第三联合优化里水库的“时间搬移”效应值得多想一步。水库本质上是把电量从某一时段搬运到另一时段的工具白天负荷高放水发电夜间风电大发蓄水这和抽水蓄能的工作逻辑非常像。理解了这一点后面再做随机优化、鲁棒优化处理风电预测误差时思路会顺很多。最后想说复现论文不是较劲“必须一模一样”而是把论文里的建模逻辑真正吃到自己脑子里。当你看到自己的代码跑出一条合理的库容轨迹、一组成熟的弃风和出力数据时那种感觉比抄一份现成代码要踏实得多。希望这篇分享能帮你在复现路上省下几个周末。