基于Matlab的高比例新能源调峰成本量化与分摊模型

这两年做新能源并网分析，遇到最多的需求不是算新能源发了多少电，而是“电网为了消纳风电光伏，到底多花了多少钱、这笔钱该由谁承担”。调峰成本这个词，以前基本只出现在学术论文里，现在却成了电网调度、电力交易部门天天要面对的现实问题。我最近把高比例可再生能源电力系统的调峰成本量化与分摊模型在 Matlab 里完整实现了一遍，从机组组合求解、成本增量计算到基于合作博弈的分摊，整条链路都跑通了。这篇就把模型思路、代码框架、参数设置和踩过的坑一次性整理出来，给正在做电力系统经济调度、电力市场课题的同学做个参考。

1. 项目背景与整体设计思路

1.1 高比例新能源系统为什么绕不开调峰成本

风电和光伏的出力由天气决定，不具备传统火电那样的可调度性。风大的夜晚、光照强烈的午间，新能源大发往往压过负荷需求，逼迫火电机组从原本高效的基荷运行状态变成频繁升降出力、甚至深度压负荷运行。而负荷高峰时段新能源出力骤减，又需要机组快速爬坡顶上。这个“追着净负荷跑”的过程，本质上是在消耗常规机组的调节能力，而这种消耗直接体现为煤耗增加、机组启停次数增多、设备寿命折损，这些加起来就是调峰成本。

高比例可再生能源系统的净负荷曲线呈现两个典型特征：一是峰谷差被显著拉大，二是净负荷变化率（爬坡速度）变得很陡。前者考验系统的调峰深度，后者考验系统的爬坡速率。火电在深度调峰区间运行时会偏离设计工况，锅炉燃烧效率下降，单位电量的煤耗明显上升；频繁启停则带来大量的启动燃料消耗和汽轮机热应力损耗。这些成本在低比例新能源时代不明显，因为机组大部分时间都稳定在额定出力附近运行，调峰需求有限。当新能源装机占比超过一定阈值后，调峰成本开始指数级增长，成了必须量化、必须分摊的显性费用。

1.2 模型整体框架与设计取舍

我搭建的这套模型分三层：数据层、优化层和分摊层。数据层负责输入负荷曲线、风光出力序列、机组技术参数与成本参数；优化层通过求解机组组合（Unit Commitment，UC）问题，得到系统在不同场景下的最小调度成本；分摊层则基于优化层输出的成本增量，把调峰成本公平合理地在各责任主体之间分配。

设计时我做了三个关键取舍：

第一，量化调峰成本采用“增量成本法”，而不是直接核算绝对成本。调峰成本本身不是一种独立可观测的费用，它混杂在燃料成本、启停成本、检修成本当中。增量成本法的思路是：先求一个不含新能源接入的基准场景调度成本，再求含新能源场景的调度成本，两者之差就是新能源并网带来的调峰成本增量。这样做的好处是基准清晰、物理意义明确，也方便后续分摊。

第二，机组组合模型用混合整数线性规划（MILP）求解。UC本质上是含启停状态变量的非凸优化问题，整数变量和连续变量并存，目前工程界标准做法就是MILP。Matlab 自带的intlinprog可以直接解，或者用 YALMIP 调 Gurobi、CPLEX 等商用求解器。我在模型里对煤耗曲线做了分段线性化处理，既保留了精度，又避开了二次规划带来的求解效率问题。

第三，分摊阶段同时实现了两种路线：一种是指标分摊法，按净负荷峰谷差增量的贡献比例分配；另一种是合作博弈法，用 Shapley 值衡量每个新能源场站对调峰成本增量的边际贡献。两种方法互相校验，工程报告里一般用指标法讲清楚“谁的责任大”，学术论文里则用 Shapley 值论证公平性。

2. 调峰成本的构成与量化建模

2.1 常规煤耗成本与出力特性建模

火电机组的煤耗特性通常用二次函数拟合：

F_i(P_i) = a_i · P_i² + b_i · P_i + c_i

其中 F_i 是单位时间燃料成本，P_i 是出力，a_i、b_i、c_i 是由机组热力试验得到的煤耗系数。二次函数能反映机组偏离最优工况后煤耗率上升的特性，但把它直接放进 MILP 会引入非线性项。工程上常用的处理办法是分段线性化：把出力范围切成若干段，每段用一条直线近似，再引入连续变量和 SOS2 约束（或者直接用分段线性约束）。

我在代码里用的是更简单的一种近似：把二次项在一个运行点附近做一阶展开，或者干脆把煤耗曲线拆成三段直线拟合。实际算例显示，分段数取 3~4 段时，煤耗成本误差能控制在 0.5% 以内，对调峰成本这种本身就带有估算性质的指标来说完全够用。如果你用的是 YALMIP，其实可以保留二次目标函数直接甩给 Gurobi 处理非凸 MIQP，但求解时间会长不少，而且非凸问题容易陷入局部最优，不如线性化稳当。

2.2 深度调峰与启停带来的额外成本

这是调峰成本量化的核心难点，也是最容易被人忽略的部分。常规机组在设计上有最小稳定出力限制，一般为其额定容量的 40%~50%。当净负荷低谷迫使机组压到最小技术出力以下运行时，就进入了深度调峰区间。

深度调峰阶段，锅炉燃烧稳定性变差，可能需要投油助燃；汽轮机末级叶片面临湿蒸汽侵蚀加剧；金属部件热应力循环加速疲劳。这些损耗在会计上不直接体现在燃料账单里，但对机组寿命的影响是实打实的，所以在量化模型中必须折算成费用。我采用的模型是：

C_dp,i = k_dp,i · (P_min,i − P_dp,i) · Δt

其中 P_min,i 是常规最小稳定出力，P_dp,i 是深度调峰阶段的实际出力，k_dp,i 是深度调峰补偿系数，单位为元/MWh，取值一般参考当地调峰辅助服务市场的补偿标准，大致在 100~300 元/MWh 区间。如果你手头有特定机组的寿命损耗试验数据，也可以把 k_dp 换成关于出力深度和持续时间的函数，精度更高但参数获取难度也更大。

启停成本同样不可忽略。一次冷态启动的燃料消耗可能相当于机组满负荷运行好几个小时，而且热态、温态、冷态启动成本差异很大，划分依据一般是停机时长。我在模型中用分段函数处理：停机 8 小时内为热启动，8~48 小时为温启动，超过 48 小时为冷启动，对应不同的启动成本系数。这个细节对结果影响不小，因为高比例新能源场景下机组的启停频次明显增加，如果统一按热启动成本算，会把调峰成本明显低估。

2.3 以基准场景差值为核心的量化逻辑

有了上面各类成本的分项模型，量化过程就变成两次 UC 求解：

第一次求基准场景：假设没有新能源接入，系统只靠常规火电满足负荷曲线，得到一个总运行成本 C_base。这个场景里机组基本处于稳定的调峰状态，成本构成相对简单。

第二次求含新能源场景：把风电、光伏出力作为负的负荷叠加进功率平衡约束，重新做 UC，得到总运行成本 C_re。由于新能源出力具有反调峰特性，这个场景的成本会显著高于基准场景。

调峰成本增量 ΔC = C_re − C_base。

这里有个容易被质疑的地方：基准场景选择不同，量化结果会差很多。如果把基准场景设为“机组全部按最低出力持续运行”，那调峰成本增量会非常大；如果设为“机组完全自由优化”，则增量会偏小。我在项目中采用的是后者，也就是从调度优化的角度，衡量新能源接入后系统最优运行方式变化带来的成本上升。这样做逻辑自洽，也方便后续用 UC 模型直接复现。报告中要写清楚基准场景的定义方式，否则评审专家一定会问。

3. 调峰成本分摊机制设计

3.1 按净负荷增量责任的分摊法

成本量化解决的是“总共花了多少”的问题，分摊解决的是“谁该承担多少”的问题。我用的第一种分摊方法基于净负荷增量责任，核心逻辑很简单：每个新能源场站对系统净负荷峰谷差和爬坡需求的贡献不同，谁造成的调峰压力大，谁就多承担调峰成本。

具体计算时，先对各场站的出力序列做基线修正。设场站 i 的日前预测出力为 w_i(t)，整个系统的原始负荷为 L(t)，净负荷为：

L_net(t) = L(t) − Σ w_i(t)

然后分别计算接入单一场站 i 后的净负荷曲线，与无该场站时的净负荷曲线比较，提取峰谷差增量 ΔD_i 和最大爬坡增量 ΔR_i。将两者归一化加权后得到每个场站的调峰责任系数：

ρ_i = α · ΔD_i / Σ ΔD_i + (1−α) · ΔR_i / Σ ΔR_i

α 是权重系数，反映调峰深度和爬坡速率在成本形成中的相对重要性，一般取 0.5~0.7。最后，场站 i 分摊的调峰成本为 ρ_i · ΔC。

这个方法计算简单、直观、可解释性强，适合工程落地。但它的缺陷也明显：没有考虑场站之间的交互作用。两个出力曲线相似的场站一起接入时，它们对净负荷峰谷差的影响不是简单叠加，而是存在“1+1<2”或者“1+1>2”的非线性效应，按比例法会高估或低估某些场站的责任。这个缺陷促成了第二种方法。

3.2 合作博弈与 Shapley 值分摊

Shapley 值来自合作博弈论，解决的就是上述“交互作用”问题。把每个新能源场站视为博弈参与者，任意一个参与者集合 S 构成一个联盟，联盟对应的调峰成本 v(S) 定义为“只有该联盟内场站接入系统时所产生的调峰成本增量”。这样，v(S) 就包含了联盟内部的协同效应。场站 i 的 Shapley 值定义为：

φ_i(v) = Σ_{S⊆N{i}} [ |S|!·(n−|S|−1)! / n! ] · [ v(S∪{i}) − v(S) ]

这个公式的物理解释是：在所有可能的接入顺序中，把场站 i 加入某个联盟时引起的调峰成本增量取平均。由于每一个接入顺序都被等概率考虑，Shapley 值具备唯一性和公平性，且满足有效性（所有场站的 Shapley 值之和等于总增量成本 ΔC）。

我在 Matlab 里实现了一个完整的 Shapley 值计算模块。前提是写一个函数v(S)，输入联盟 S 的成员集合，内部构造对应的 UC 模型并求解，返回该联盟场景下的调峰成本增量。然后枚举所有不含 i 的子集，按公式累加。核心循环代码如下：

function phi = shapley_cost(v_func, N) % v_func: 联盟成本函数句柄，输入为0/1掩码向量，返回该联盟的调峰成本增量 phi = zeros(N, 1); for i = 1:N for mask = 0:2^N - 1 % 跳过包含参与者i的联盟 if bitand(mask, 2^(i-1)) ~= 0 continue; end nS = sum(bitget(mask, 1:N)); % 联盟S的大小 w = factorial(nS) * factorial(N - nS - 1) / factorial(N); cost_with = v_func(bitor(mask, 2^(i-1))); cost_without = v_func(mask); phi(i) = phi(i) + w * (cost_with - cost_without); end end end

这段代码的瓶颈在于v_func每次都要调用一次 UC 求解器。当参与分摊的新能源场站数量为 N 时，需要求解 2^N 次 UC 问题。N=8 时是 256 次，N=12 时是 4096 次，再往上就非常吃力了。所以该方法更适合场站数量较少、或者先把同质场站聚合成几个集群再计算的场景。

3.3 核仁法与工程简化

Shapley 值处理的是“公平性”问题，但对部分联盟的约束可能引发争议。核仁法（Nucleolus）的思路是寻找一组分摊结果，使得所有联盟中“对分得不满意程度”最大的那个联盟的不满意度尽可能小，也就是最小化最大超额。核仁法的计算结果通常需要通过线性规划迭代求解，Matlab 里可以用linprog实现。

从项目实际交付角度，我建议不要一上来就上核仁法。原因很简单：核仁法需要解一个规模不小的 LP 序列，而且结果对成本函数 v(S) 的数值误差很敏感，一旦 UC 求解精度不够，核仁法很容易算出数值上不稳定的小数。工程报告中更常见的折中方案是：用 Shapley 值计算结果作为参考，用指标法计算结果作为结算依据，两者偏差控制在合理范围内即可。

我做过一个 6 个新能源场站的算例对比，指标法和 Shapley 值的分摊结果偏差大约在 3%~8% 之间，在可接受范围内。所以如果你不是做纯理论研究，建议优先把指标法做扎实，再补一个 Shapley 值做交叉验证。

4. Matlab 实现的前期准备

4.1 数据结构与算例设计

开始写代码前，先把数据结构设计好。我建议用结构体或表格把输入数据组织清晰，不要到处散落变量。推荐建一个input_data结构体，包含以下字段：

机组参数表gen_param的每一行对应一台火电机组，六列分别是额定容量、最大出力、最小稳态出力、深度调峰极限出力、爬坡速率、初始状态。成本参数表则包含煤耗系数 a、b、c，启停成本，深度调峰补偿系数。这些数据可以从 IEEE RTS 算例系统获取，也可以根据实际电厂参数整理。算例的典型配置是 24 小时调度周期、时间分辨率为 1 小时，机组数量取 6~10 台，新能源场站取 3~6 个。

时间尺度上注意一点：如果研究的是调峰成本和系统爬坡压力，1 小时分辨率会低估爬坡约束的严格性。建议在关键场景下把时间分辨率细化到 15 分钟，代价是 UC 求解时间翻几倍。我一般先在 1 小时尺度上跑通逻辑，确认模型无误后再加密到 15 分钟做灵敏度分析。

4.2 建模方式选择：intlinprog 还是 YALMIP

Matlab 里实现 UC 模型有两条主流路线：纯intlinprog和YALMIP + 外部求解器。

纯intlinprog的好处是不依赖第三方工具箱，Matlab 自带优化工具箱就能跑。坏处是所有约束都需要手工写成矩阵不等式形式，变量索引关系稍不留神就错。尤其是最小启停时间约束和爬坡约束涉及相邻时段耦合，矩阵构造非常繁琐，调试体验比较痛苦。

YALMIP 的好处是建模语言接近数学表达式，binvar、sdpvar定义变量，约束直接写在Constraints集合中，目标函数也很直观。求解时只要机器上装了 Gurobi 或 CPLEX，通过sdpsettings('solver', 'gurobi')指定求解器即可。缺点是 YALMIP 本身是第三方工具箱，需要额外安装，而且它对高版本 Matlab 的适配偶尔会有小问题。

我的选择是：项目原型用 YALMIP 建模，快速验证逻辑；交付版本如果客户环境没有 YALMIP，再改写为intlinprog。下面给出的核心代码以 YALMIP 为主，同时我会标注intlinprog的等价实现思路。

5. 核心代码实现与结果解析

5.1 主程序与机组组合求解

主程序分四步：加载数据、定义决策变量、构造约束与目标、调用求解器。先看整体框架：

% 主程序 main_uc.m clear; clc; close all; load('input_data.mat'); T = 24; % 调度时段数 ng = size(gen_param, 1); % 火电机组数 % 决策变量 u = binvar(T, ng, 'full'); % 0/1启停状态 P = sdpvar(T, ng, 'full'); % 机组出力 SU = binvar(T, ng, 'full'); % 启动动作 SD = binvar(T, ng, 'full'); % 停机动作 Constraints = []; % 功率平衡约束 for t = 1:T wind_total = sum(wind_profile(t, :)); pv_total = sum(pv_profile(t, :)); Constraints = [Constraints, sum(P(t, :)) + wind_total + pv_total == load_profile(t)]; end % 出力上下限约束（含深度调峰区间） for t = 1:T for i = 1:ng % Pmin_cv是常规最小稳定出力，P_deep是深度调峰极限 Constraints = [Constraints, u(t,i) * gen_param(i,4) <= P(t,i) <= u(t,i) * gen_param(i,2)]; end end % 爬坡约束 for t = 2:T for i = 1:ng Constraints = [Constraints, P(t,i) - P(t-1,i) <= gen_param(i,5)]; Constraints = [Constraints, P(t-1,i) - P(t,i) <= gen_param(i,5)]; end end % 启动/停机动作约束 for t = 2:T Constraints = [Constraints, SU(t,:) >= u(t,:) - u(t-1,:)]; Constraints = [Constraints, SD(t,:) >= u(t-1,:) - u(t,:)]; end % 目标函数：煤耗 + 启停 + 深度调峰损耗 Objective = 0; for t = 1:T for i = 1:ng Objective = Objective + cost_param(i,1) * P(t,i)^2 ... + cost_param(i,2) * P(t,i) ... + cost_param(i,3) * u(t,i); Objective = Objective + cost_param(i,4) * SU(t,i) ... + cost_param(i,4) * SD(t,i); % 深度调峰损耗：当出力低于常规最小稳定出力时按差值计费 Objective = Objective + cost_param(i,5) * max(0, gen_param(i,3) - P(t,i)); end end ops = sdpsettings('solver', 'gurobi', 'verbose', 0); optimize(Constraints, Objective, ops);

这里目标函数里出现了平方项P(t,i)^2，YALMIP 会把问题转成 MIQP 丢给 Gurobi。如果求解器不支持 MIQP，或者想用intlinprog，就必须先做分段线性化。我常用的是把二次函数在几个运行点处线性插值，然后用二进制变量激活对应分段，这部分代码比较长，放在函数linearized_cost.m里统一处理。实际项目中，我强烈建议默认走线性化路线，虽然系数矩阵变大，但求解稳定性好得多。

深度调峰损耗项里的max(0, gen_param(i,3) - P(t,i))本质上是 max 函数，非光滑。YALMIP 在 MIQP 框架下能处理这种形式，但为了稳妥，我一般引入辅助变量dp(t,i)代替 max，并加约束：

dp(t,i) >= gen_param(i,3) - P(t,i); dp(t,i) >= 0;

然后用 dp 的加权和作为该部分成本。这个技巧建议收藏，很多需要分段计费的成本（比如阶梯电价、阶梯碳价）都可以照这个套路建模。

5.2 成本量化模块

UC 求解完成后，需要把结果拆解成成本结构清单。我写了一个后处理函数，输出各类成本的数值：

% 成本量化模块 cost_report.m function report = cost_report(P, u, SU, SD, gen_param, cost_param, T, ng) report.fuel_cost = 0; report.startup_cost = 0; report.shutdown_cost = 0; report.deep_peakshaving_cost = 0; for t = 1:T for i = 1:ng report.fuel_cost = report.fuel_cost + ... cost_param(i,1)*P(t,i)^2 + cost_param(i,2)*P(t,i) + cost_param(i,3)*u(t,i); report.startup_cost = report.startup_cost + cost_param(i,4)*SU(t,i); report.shutdown_cost = report.shutdown_cost + cost_param(i,4)*SD(t,i); if P(t,i) < gen_param(i,3) report.deep_peakshaving_cost = report.deep_peakshaving_cost + ... cost_param(i,5) * (gen_param(i,3) - P(t,i)); end end end end

跑完基准场景和含新能源场景后，ΔC = C_re − C_base就得到调峰成本增量。我习惯把这个增量再拆成“煤耗增量”“启停增量”“深度调峰增量”三部分，画一张堆积柱状图。这样能直观看出新能源并网后成本涨在哪里：是频繁启停带来的，还是深度调峰工况变多导致的。

5.3 Shapley 值分摊模块

在 3.2 节已经给出了 Shapley 值计算的核心循环。实际使用时，v_func函数需要接收联盟掩码向量，返回该联盟对应的调峰成本增量。实现如下：

% 联盟成本函数 function deltaC = coalition_cost(mask, Nw, base_cost, ... wind_profile, load_profile, gen_param, cost_param) % 根据mask确定接入哪些新能源场站 active = find(mask(1:Nw) == 1); wind_active = sum(wind_profile(:, active), 2); T = length(load_profile); % 构造UC模型并求解，返回总成本 total_cost = solve_uc(wind_active, 0, gen_param, cost_param); deltaC = total_cost - base_cost; end

注意每次调用solve_uc都是独立求解一个 UC 问题，中间不共享信息。当联盟数量较多时，这是很自然的并行化切入点。Matlab 的parfor可以直接替换外层循环，把 N 个 Shapley 值计算并行化，实测 6 个场站、24 个时段的问题，并行后耗时能从几分钟降到几十秒。

另一个工程技巧是对称性缩减。如果两个场站的出力序列完全一致（比如同一风电基地的两期项目），它们在 Shapley 值计算中地位等同，可以合并为一个参与者，权重按场站容量折算，能大幅减少联盟数量。同时，对于出力曲线相似度超过阈值的场站，可以先用 K-means 聚类，把同质场站聚合成一个集群，再对集群做 Shapley 值分摊，最后按组内容量比例二次分摊。这种方法在参与主体多、数据差异小的场景下特别好用。

5.4 结果可视化

结果可视化我用三张图搞定：

第一张是净负荷曲线图，把原始负荷、新能源出力、净负荷画在同一张图上，直观展示峰谷差变化。

第二张是机组组合甘特图，横轴是时间，纵轴是机组编号，用色块表示启停状态。这张图能快速暴露模型问题，比如某台机组频繁启停、或者出力长时间贴着上下限。

第三张是分摊结果堆积柱状图，每个新能源场站一根柱子，柱子里按常规煤耗、启停、深度调峰三类成本分层填色。这张图是交付报告的核心，直接回答“谁该承担多少”。

figure; bar_handle = bar(cost_share_matrix, 'stacked'); legend({'煤耗成本增量', '启停成本增量', '深度调峰成本增量'}, 'Location', 'northwest'); set(gca, 'XTickLabel', {'风电场1', '风电场2', '光伏站1', '光伏站2'}); ylabel('分摊成本（万元）'); title('各新能源场站调峰成本分摊结果'); grid on;

6. 工程落地中的常见坑与调试心得

6.1 求解不出可行解的排查思路

这是做 UC 最常遇到的坑。模型报“infeasible”的时候，先不要怀疑求解器，按下面的顺序排查：

第一，检查功率平衡约束是否和机组总容量矛盾。比如负荷低谷时段，所有机组都压到最小技术出力，仍然大于净负荷，那就说明系统中常规机组的最小出力总和过高，必须允许某些机组停机。这时要检查最小启停时间约束是否阻止了这些机组停机。

第二，检查爬坡约束和出力上下限是否冲突。有些机组爬坡速率太慢，在 1 小时尺度下，从最小出力爬到最大出力需要好几个时段，如果在相邻时段强制要求它大幅度调整，就会无解。解决办法是把爬坡约束放到机组启停时间窗口中处理，或者细化时间分辨率。

第三，用松弛法定位不可行约束。在约束后面加一个小的非负松弛变量，把约束写成“变量 + slack >= 需求”的形式，求解后看哪些 slack 不为零，不可行的根源就暴露出来了。这个方法虽然会让目标函数略微失真，但作为调试工具非常有效。

6.2 Shapley 值计算组合爆炸的简化处理

前面提过，N 个参与者的 Shapley 值需要 2^N 次 UC 求解，N 超过 12 基本就受不了。我实际项目里用过两个替代方案：

方案一是抽样估计。从所有联盟中等概率抽取 M 个子集，用这些子集的 v(S) 值通过回归或者随机采样近似 Shapley 值，Matlab 里可以直接调fitrgp做高斯过程代理模型，M 取 1000~3000 时精度就比较可观了。

方案二是 Douglas 分解或者按时间窗口拆解。调峰成本本质上与净负荷曲线形态强相关，可以把 24 小时按峰、平、谷分成三个典型时段，先分别计算各时段的调峰成本增量，再在各时段内部用 Shapley 值分摊，最后按时段加权汇总成全天结果。时间窗口拆解后每个窗口内的 UC 问题规模变小，求解速度显著提升，而且不会损失太多精度。

如果只是做工程结算，我更推荐前面的指标分摊法，配合 3~5 个典型场景做灵敏度分析，基本就能满足报告要求了。

6.3 结果合理性校验技巧

模型跑完不等于结果可靠，拿到数字先做三件事：

第一，检查总成本与分项成本比例。正常情况下煤耗成本占大头，启停成本次之，深度调峰损耗最小。如果算出来深度调峰成本占了一半以上，先检查是不是把深度调峰区间设置得太宽，或者补偿系数 k_dp 取值过大。

第二，检查 Shapley 值之和是否严格等于调峰成本增量 ΔC。这是 Shapley 值有效性的数学性质，任何代码实现都不能违背。如果对不上，优先检查联盟成本函数 v(S) 是否一致，特别是基准成本 base_cost 有没有在每次调用时都被正确扣除。

第三，把结果和物理直觉对照。比如某风电场夜间的出力高峰正好对应负荷低谷，那么它的深度调峰责任一定大；某光伏站在午间强出力时段反而帮助削峰，它的分摊额应相对较小。如果计算结果与这些直觉明显矛盾，大概率是数据预处理出了问题，比如时间对齐错误、出力序列单位不统一。

我在实际项目里曾经因为把小时数据和分钟数据混在一起，导致净负荷曲线出现严重跳变，UC 结果完全失真。查了一天最后发现是xlsread读数据时多读了一列。所以强烈建议在做任何计算前，先用plot把每条输入曲线画出来看一遍，这个习惯能省掉大量排查时间。

最后再说一个小技巧。整套模型跑完后，一定要把机组组合结果和成本结构导出成 Excel 或 CSV，方便在报告中直接引用。Matlab 的writetable可以一行代码把表格数据写出去：

result_table = table(load_profile, sum(wind_profile,2), sum(P,2), ... 'VariableNames', {'Load', 'Wind', 'Thermal'}); writetable(result_table, 'dispatch_result.xlsx');

这个项目做完后，我最大的体会是：调峰成本建模的难点往往不在数学公式，而在“基准场景定义”和“成本边界划分”这两个看似简单的问题上。基准场景稍微变一下，成本增量就能差出百分之三四十；深度调峰成本算不算启停损失，各家理解也不一致。做这类课题之前先把边界条件和各方利益诉求摸清楚，想清楚分摊结果拿来做什么用——是作为市场结算依据，还是作为政策评估参考——再决定用多复杂的模型。模型不是越复杂越好，但代码框架一定得留足扩展空间，因为等新能源占比再上一个台阶，储能、需求响应、抽蓄这些调节资源都会被纳入进来，现在的分摊模型很可能要在这个框架上继续长出新东西。