简介：本资源面向电力系统优化调度方向的毕业设计与科研入门者，围绕《计及调峰主动性的风光水火储多能系统互补协调优化调度》一文提供完整MATLAB复现代码。针对大规模风光并网后系统调峰能力不足的问题，代码实现了含调峰补偿与分摊、火电调峰主动性约束的分层优化调度方案：上层以净负荷波动最小和储能收益最大为目标，下层以火电运行成本最小和弃电量最小为目标，并基于分解协调思想交替迭代求解，算例采用改进IEEE30节点系统进行多场景仿真。压缩包共21个文件，约509KB，以19个m脚本为核心，涵盖约束构建、目标函数、导纳与关联矩阵生成、深度调峰及主程序等模块，另含1个png结果图与1个lp数据文件，便于直接运行与二次修改。已有227人学习下载，适合需要快速理解模型结构、复现仿真结果并撰写论文的读者参考。

1. 从一篇知网论文到可运行代码：风光水火储互补调度到底在算什么

如果你在知网搜到《计及调峰主动性的风光水火储多能系统互补协调优化调度》，下载了 PDF，看完公式却不知道怎么落地成能跑的 MATLAB 程序，那这篇就是写给你的。它解决的核心问题很具体：给定一个包含风电、光伏、水电、火电和储能的多能系统，在一天 96 个调度时段里，怎么分配各机组的出力，使得总运行成本最低，同时让火电机组主动参与调峰而不是被动被调度。适合两类人：一类是电力系统优化调度的研究生，需要复现论文里的模型跑出仿真曲线；另一类是做新能源场站功率管理的工程师，想搞清楚多能互补的优化逻辑到底怎么写成代码。

调峰主动性这个词听起来玄学，说白了就是给火电机组一个额外的收益信号，让它愿意在负荷低谷时压低出力、在高峰时快速顶上去，而不是只按调度指令被动执行。传统经济调度模型里火电只算煤耗成本，调峰主动性模型会额外引入调峰补偿系数和深度调峰成本项。风光水火储五类电源的互补协调，本质是一个带约束的大规模非线性优化问题，MATLAB 里通常用 YALMIP 或 CVX 建模，再调 Gurobi 或 CPLEX 求解。下面从模型结构、参数设置、代码实现到避坑，一步步拆开讲。

2. 多能系统互补调度的数学模型与 MATLAB 建模路径

2.1 目标函数里到底有哪几项成本

论文里的目标函数通常写成总运行成本最小化，展开后包含五块：火电燃料成本、火电调峰补偿成本、水电运行成本、储能充放电损耗成本、弃风弃光惩罚成本。火电燃料成本一般用二次函数拟合，形式是 a·P² + b·P + c，其中 P 是火电出力。调峰补偿成本是这篇论文的核心创新点，通常写成调峰深度乘以补偿单价，调峰深度定义为火电额定容量与当前出力的差值再除以额定容量。

风光水火储的互补体现在约束的耦合上。功率平衡约束要求每个时段五类电源出力之和等于负荷需求。火电有爬坡约束和出力上下限，水电有水量平衡和库容约束，储能有 SOC 递推和充放电功率限制，风光出力在论文里一般作为预测值直接给定，不参与优化变量。这些约束在 MATLAB 里用矩阵形式表达，YALMIP 的语法最直观。

我一般会先把目标函数和约束在纸上写成标准形式，再翻译成代码。直接上手写代码容易漏约束，尤其是储能的 SOC 递推和火电的爬坡耦合，漏一个约束求解结果就会离谱。

2.2 用 YALMIP 搭建优化模型的完整代码框架

下面是一个可运行的最小框架，以 24 个调度时段为例，实际论文用 96 点只需改 T 的值。代码里用 YALMIP 定义决策变量和约束，求解器指定为 Gurobi，如果没有 Gurobi 可以换成 intlinprog 或 fmincon，但速度会慢很多。

% 多能系统互补协调优化调度 - 最小可运行框架 clear; clc; close all; T = 24; % 调度时段数，论文常用96 dt = 1; % 每时段小时数 % ---- 负荷与风光预测数据（此处用典型日曲线替代）---- P_load = [220 210 200 195 200 230 280 320 350 360 370 380 ... 375 365 355 340 330 340 360 370 350 310 270 240]; % MW P_wind = [80 75 70 65 60 55 50 45 50 60 70 85 ... 90 95 100 105 110 115 120 115 105 95 90 85]; % MW P_pv = [0 0 0 0 0 5 20 45 70 95 115 130 ... 135 130 120 100 75 45 15 0 0 0 0 0]; % MW % ---- 火电参数 ---- P_th_max = 300; P_th_min = 120; % 火电出力上下限 MW a_th = 0.0005; b_th = 0.3; c_th = 20; % 燃料成本系数 ramp_up = 60; ramp_down = 60; % 爬坡速率 MW/h k_peak = 50; % 调峰补偿单价 元/MW % ---- 水电参数 ---- P_h_max = 100; P_h_min = 0; W_in = 200 * ones(1, T); % 来水量 m3/s eta_h = 8.5; % 水电出力系数 % ---- 储能参数 ---- E_max = 200; E_min = 20; % 储能容量上下限 MWh P_ch_max = 50; P_dis_max = 50; % 充放电功率上限 MW eta_ch = 0.95; eta_dis = 0.95; % 充放电效率 SOC0 = 100; % 初始SOC MWh % ---- 决策变量 ---- P_th = sdpvar(1, T); % 火电出力 P_h = sdpvar(1, T); % 水电出力 P_ch = sdpvar(1, T); % 储能充电功率 P_dis = sdpvar(1, T); % 储能放电功率 E_soc = sdpvar(1, T); % 储能SOC u_ch = binvar(1, T); % 充电状态0-1 u_dis = binvar(1, T); % 放电状态0-1 % ---- 目标函数 ---- Cost_fuel = sum(a_th * P_th.^2 + b_th * P_th + c_th); Cost_peak = sum(k_peak * max(0, (P_th_max - P_th) / P_th_max)); Cost_hydro = sum(0.1 * P_h); Cost_ess = sum(0.05 * (P_ch + P_dis)); Cost_curtail = sum(100 * max(0, P_wind + P_pv - (P_load - P_th_min))); Objective = Cost_fuel + Cost_peak + Cost_hydro + Cost_ess + Cost_curtail; % ---- 约束条件 ---- Constraints = []; % 功率平衡 for t = 1:T Constraints = [Constraints, ... P_th(t) + P_h(t) + P_wind(t) + P_pv(t) + P_dis(t) - P_ch(t) == P_load(t)]; end % 火电约束 Constraints = [Constraints, P_th_min <= P_th <= P_th_max]; for t = 2:T Constraints = [Constraints, ... -ramp_down <= P_th(t) - P_th(t-1) <= ramp_up]; end % 水电约束 Constraints = [Constraints, P_h_min <= P_h <= P_h_max]; for t = 1:T Constraints = [Constraints, P_h(t) <= eta_h * W_in(t)]; end % 储能约束 Constraints = [Constraints, E_min <= E_soc <= E_max]; Constraints = [Constraints, 0 <= P_ch <= P_ch_max * u_ch]; Constraints = [Constraints, 0 <= P_dis <= P_dis_max * u_dis]; Constraints = [Constraints, u_ch + u_dis <= 1]; % 充放电互斥 for t = 2:T Constraints = [Constraints, ... E_soc(t) == E_soc(t-1) + eta_ch * P_ch(t) - P_dis(t) / eta_dis]; end Constraints = [Constraints, E_soc(1) == SOC0 + eta_ch * P_ch(1) - P_dis(1) / eta_dis]; % ---- 求解 ---- ops = sdpsettings('solver', 'gurobi', 'verbose', 1); sol = optimize(Constraints, Objective, ops); if sol.problem == 0 P_th_val = value(P_th); P_h_val = value(P_h); P_ch_val = value(P_ch); P_dis_val = value(P_dis); E_soc_val = value(E_soc); fprintf('总运行成本: %.2f 元\n', value(Objective)); else disp('求解失败，检查约束是否可行'); disp(sol.info); end

这段代码的逻辑分四层：第一层定义物理参数和预测数据，第二层声明决策变量，第三层构造目标函数和约束，第四层调求解器。参数说明几个关键点：k_peak是调峰补偿单价，论文里一般取 30 到 80 元/MW 之间，取值越大火电越愿意压低出力；eta_h是水电出力系数，等于水头乘以流量转换效率再乘以重力加速度，典型值在 8 到 9 之间；储能充放电互斥约束u_ch + u_dis <= 1必须加，否则求解器会让储能同时充放电来套利，结果物理上不成立。

2.3 从 24 时段扩展到 96 时段要注意什么

论文里通常用 96 点，每 15 分钟一个时段。直接把 T 改成 96 后，爬坡约束要按 15 分钟折算，原来 60 MW/h 的爬坡速率变成 15 MW/15min。储能 SOC 递推里的效率系数也要注意，如果充放电功率单位是 MW，时段长度是 0.25 小时，那么 SOC 变化量要乘以 dt。很多复现翻车就翻在这里，SOC 曲线跑出来是一条直线或者直接越限，就是因为时段长度没折算。

另外 96 时段下二进制变量有 192 个，加上连续变量总共接近 500 个，Gurobi 求解一般 10 秒内能出结果。如果换成 fmincon 会非常慢，因为 fmincon 不擅长处理混合整数问题。建议至少装一个 Gurobi 学术版，或者用 YALMIP 自带的bnb求解器但接受更长的求解时间。

3. 调峰主动性怎么在代码里体现：补偿系数与深度调峰约束

3.1 调峰补偿成本的三种建模方式对比

调峰主动性在论文里有三种常见写法，选哪种直接影响求解难度和结果合理性。第一种是线性补偿，成本项写成 k·(P_max - P)，P 越小补偿越高，优点是线性好求解，缺点是火电会无限制压低出力直到下限。第二种是分段补偿，把调峰深度分成常规调峰和深度调峰两段，深度调峰给更高补偿单价，更贴近实际电力市场规则，但需要引入额外的 0-1 变量来标记分段。第三种是阶梯补偿，用阶梯函数拟合，非线性最强，一般用 MILP 线性化后求解。

我一般先用线性补偿跑通，确认功率平衡和储能逻辑没问题，再换成分段补偿。分段补偿的代码改动集中在目标函数和约束上，下面给出分段补偿的核心片段。

% 分段调峰补偿：常规调峰深度0-40%，深度调峰40%-60% d_peak = sdpvar(1, T); % 调峰深度 d_reg = sdpvar(1, T); % 常规调峰深度 d_deep = sdpvar(1, T); % 深度调峰深度 y_deep = binvar(1, T); % 是否进入深度调峰 for t = 1:T Constraints = [Constraints, ... d_peak(t) == (P_th_max - P_th(t)) / P_th_max, ... d_reg(t) + d_deep(t) == d_peak(t), ... d_reg(t) <= 0.4, ... d_deep(t) <= 0.2 * y_deep(t), ... d_deep(t) >= 0.001 * y_deep(t)]; end Cost_peak = sum(30 * d_reg * P_th_max + 80 * d_deep * P_th_max);

这里y_deep是二进制变量，当调峰深度超过 40% 时置 1，深度调峰部分才会计入高补偿。d_deep的下界设 0.001 是为了避免求解器把y_deep置 1 但d_deep为 0 的退化情况。补偿单价 30 和 80 是示例值，实际论文里会根据当地调峰辅助服务市场规则设定。

3.2 深度调峰对火电最小出力的影响

火电在深度调峰时最小出力可以降到额定容量的 30% 甚至更低，但常规调峰时最小出力一般不低于 50%。这个差异要在约束里体现，否则优化结果会给出一个物理上不可行的出力曲线。常见做法是引入两个最小出力参数：P_th_min_reg和P_th_min_deep，用y_deep来切换。

P_th_min_reg = 0.5 * P_th_max; % 常规调峰最小出力 P_th_min_deep = 0.3 * P_th_max; % 深度调峰最小出力 for t = 1:T Constraints = [Constraints, ... P_th(t) >= P_th_min_reg - (P_th_min_reg - P_th_min_deep) * y_deep(t), ... P_th(t) <= P_th_max]; end

这个约束的逻辑是：当y_deep为 0 时，火电出力下限是P_th_min_reg；当y_deep为 1 时，下限放宽到P_th_min_deep。这样火电只有在进入深度调峰状态时才能压到更低出力，符合实际机组的运行特性。

3.3 求解结果怎么验证合理性

跑出结果后不要只看总成本，要检查四条曲线：火电出力曲线是否平滑、储能 SOC 是否在合理区间内波动、弃风弃光量是否集中在风光大发时段、调峰深度是否在补偿单价变化时发生跳变。我一般会画四张子图放在一个 figure 里，横轴都是时段，纵轴分别是各电源出力、SOC、调峰深度和弃风弃光功率。

如果火电出力曲线出现剧烈震荡，说明爬坡约束没起作用或者时段长度没折算。如果 SOC 曲线一直贴在上限或下限，说明储能容量设置太小或者充放电功率限制太松。如果弃风弃光量在负荷高峰时段反而很大，说明功率平衡约束写反了符号。这些检查比看总成本数字有用得多。

4. 风光水火储五类电源的约束细节与参数取值

4.1 风电和光伏在优化模型里怎么处理

论文里风电和光伏一般不作为决策变量，而是作为预测值直接代入功率平衡约束。但有些论文会引入弃风弃光变量，把实际消纳的风光功率作为决策变量，上限是预测值。这两种处理方式的区别在于：前者假设风光全额消纳，后者允许弃风弃光但计入惩罚成本。

我倾向于后者，因为实际系统里弃风弃光不可避免，而且引入弃风弃光变量后模型更灵活。代码改动很小，把P_wind和P_pv换成决策变量P_wind_use和P_pv_use，加上上限约束等于预测值，目标函数里加惩罚项。

P_wind_use = sdpvar(1, T); P_pv_use = sdpvar(1, T); Constraints = [Constraints, 0 <= P_wind_use <= P_wind, 0 <= P_pv_use <= P_pv]; Cost_curtail = sum(200 * (P_wind - P_wind_use) + 200 * (P_pv - P_pv_use));

惩罚系数取 200 元/MWh 是常见值，取值越大越倾向于全额消纳。如果取值太小，求解器会大量弃风弃光来降低火电调峰压力，结果不符合实际调度逻辑。

4.2 水电的水量平衡约束怎么写

水电约束比火电复杂，因为涉及水量平衡和库容约束。简化模型里只考虑发电用水，水量平衡写成：上时段库容 + 来水量 - 发电用水量 = 本时段库容。发电用水量和出力之间的关系是 P_h = eta_h * Q_h，其中 Q_h 是发电流量。

V_max = 5000; V_min = 1000; % 库容上下限 万m3 V0 = 3000; % 初始库容 Q_h = sdpvar(1, T); % 发电流量 V = sdpvar(1, T); % 库容 for t = 1:T Constraints = [Constraints, ... V_min <= V(t) <= V_max, ... 0 <= Q_h(t) <= W_in(t), ... P_h(t) == eta_h * Q_h(t)]; end for t = 2:T Constraints = [Constraints, ... V(t) == V(t-1) + (W_in(t) - Q_h(t)) * 3600 / 1e4]; end Constraints = [Constraints, V(1) == V0 + (W_in(1) - Q_h(1)) * 3600 / 1e4];

这里单位换算是坑点：来水量单位是 m³/s，时段长度是小时，库容单位是万 m³，所以乘以 3600 再除以 1e4。如果单位不统一，库容曲线会跑出离谱的值。我一般会在代码里加注释标明每个变量的单位，避免后面改参数时搞混。

4.3 储能 SOC 递推的效率和自放电怎么处理

储能模型的核心是 SOC 递推，前面代码里已经写了基本形式。实际论文里还会考虑自放电率，每月自放电 1% 到 3%，折算到每 15 分钟大约是 0.001% 到 0.003%。自放电项写成E_soc(t) = (1 - sigma) * E_soc(t-1) + ...，其中 sigma 是自放电率。

自放电率对结果影响很小，但加了之后模型更完整。如果审稿人要求考虑自放电，就加上；如果只跑仿真看趋势，不加也行。另外储能初始 SOC 和终止 SOC 是否相等也是一个选择，有些论文要求一天结束后 SOC 回到初始值，有些允许自由结束。要求回到初始值的话加一条约束E_soc(T) == SOC0，但可能导致求解不可行，因为一天内充放电量不一定能完全平衡。

4.4 五类电源参数汇总表

下面这张表汇总了常见参数取值范围，方便复现时对照检查。参数取值来自公开论文的典型值，实际项目要根据具体机组和场站数据调整。

5. 复现时最容易翻车的五个地方

5.1 求解器返回 infeasible 但不知道哪条约束冲突

现象：optimize返回 problem 非 0，提示 infeasible，但约束有几十条，不知道哪条导致不可行。原因：功率平衡约束和机组出力上下限冲突，比如负荷太低而火电最小出力太高，风光又必须全额消纳，导致没有可行解。解决：先把二进制变量和分段补偿去掉，用线性模型跑一遍，确认基本约束可行后再逐步加复杂约束。YALMIP 可以用optimize(Constraints, [], ops)只检查可行性，不优化目标，这样能快速定位是约束冲突还是求解器问题。

5.2 储能 SOC 曲线跑成一条直线

现象：SOC 始终等于初始值，充放电功率都是 0。原因：储能充放电损耗成本设得太高，或者充放电互斥约束写成了u_ch + u_dis == 1，导致每个时段必须充或放，求解器选择不充不放时约束冲突。解决：检查互斥约束应该是<= 1而不是== 1，损耗成本系数不要超过 0.1 元/MWh，否则储能没有套利空间。

5.3 火电出力在爬坡约束下仍然跳变

现象：相邻时段火电出力变化超过爬坡速率限制。原因：爬坡约束的时段长度没折算，96 点调度时爬坡速率还是按小时给的。解决：把爬坡约束写成-ramp_down * dt <= P_th(t) - P_th(t-1) <= ramp_up * dt，其中 dt 是时段长度，96 点时 dt=0.25。

5.4 调峰补偿导致火电出力贴下限不参与调峰

现象：火电出力一直等于最小出力，调峰深度拉满，补偿成本很高但总成本反而更低。原因：补偿单价设得太高，火电压出力拿补偿比多发电省燃料更划算。解决：补偿单价不能超过燃料成本的边际值，一般控制在 50 元/MW 以内，同时加一条约束限制调峰深度不超过 60%，避免火电无限制压低出力。

5.5 96 点调度求解时间超过 10 分钟

现象：T=96 时 Gurobi 跑很久不出结果。原因：二进制变量太多，分段补偿和储能互斥各引入 96 个 binvar，加上其他 0-1 变量总共超过 300 个。解决：把分段补偿的二进制变量用 SOS2 约束替代，或者把储能互斥约束松弛成连续约束，求解后再取整。另外可以设置 Gurobi 的 MIPGap 为 0.01，允许 1% 的次优解，求解时间能缩短一半以上。

6. 把论文模型改成自己的场景：参数替换与结果对比技巧

复现论文模型只是第一步，真正有用的是把它改成你自己的场景。我一般会做三组对比：基准场景用论文原始参数，场景 A 把调峰补偿单价提高 50%，场景 B 把储能容量扩大一倍，然后对比三组的总成本、弃风弃光率和火电调峰深度分布。这样能看出哪个参数对结果影响最大，也方便写论文的灵敏度分析章节。

参数替换时最容易忽略的是单位一致性。论文里功率用 MW，电量用 MWh，时间用小时，但有些论文用 kW 和 kWh，替换时要把所有参数统一乘或除 1000。我习惯在代码开头写一个单位换算区块，把所有参数统一到 MW 和 MWh，后面不再改。

结果对比我一般用三张图：第一张是各场景下的火电出力曲线叠在一起，看调峰深度差异；第二张是储能 SOC 曲线对比，看储能利用率；第三张是弃风弃光功率柱状图，看消纳改善。三张图放在一个 figure 里用 subplot 排列，导出 EPS 格式方便插入论文。MATLAB 导出 EPS 用print -depsc2 filename.eps，如果中文乱码就把字体设成 SimSun 或 SimHei。

还有一个技巧是把优化结果存成 mat 文件，后面画图不用重新求解。save('result.mat', 'P_th_val', 'P_h_val', 'E_soc_val', 'P_ch_val', 'P_dis_val')，下次直接load就行。这样调图的时候不用等求解器，效率高很多。

最后说一个我踩过的坑：论文里的参数表看起来很清楚，但有些参数是标幺值，有些是实际值，混在一起用会导致结果差几个数量级。我现在的习惯是拿到论文先列一张参数表，标清楚每个参数的单位和来源，确认全部统一后再写代码。这个习惯帮我省了很多后悔药。希望帮到你。

本文还有配套的精品资源，点击获取