含风电的电力系统动态经济调度：随机场景与MILP建模详解

1. 从静态到动态：风电随机性到底难在哪

接到这个题目的时候，我第一反应是：很多人把“动态经济调度”和“含风电的经济调度”当成两件事来做。实际上，这两个难点叠在一起，才是这个模型的真正核心——既要处理常规机组跨时段的启停、爬坡约束，又要处理风电出力不可精确预知带来的随机性。两个维度一旦耦合，问题性质立刻就从简单的线性规划变成了混合整数规划，而且在某些时段约束下还可能无解或产生病态解。

先说清楚“动态”是什么意思。传统静态经济调度只关心某个单独时段内，怎么分配各机组出力让总煤耗成本最低，它不关心机组上一个小时在不在运行、能不能在十分钟内把出力提上去。但实际电网调度是按小时甚至按一刻钟滚动执行的，机组不可能瞬间从50MW跳到200MW，也不可能刚停机就立刻开机。动态经济调度就是在时间轴上把这些物理限制全部加进去：机组启停状态是0/1整数变量，出力上下限和爬坡速率是连续约束，两个相邻时段之间的状态转移会被严格限制。这样一来，模型规模立刻增大，求解难度也上了一个台阶。

再说“随机性”这一层。风电出力取决于风速，而风速预测不可能绝对准确。一个常见的处理方式是：用预测误差的分布函数生成大量可能的风电出力场景，每个场景附带一个概率，然后让调度模型在所有场景下都满足安全约束，同时期望成本最小。这就是随机规划里的“场景法”，也是本文这套代码的核心思路。与之相对的还有一种鲁棒优化思路，只考虑最坏场景，但这会过于保守，实际工程中更常用场景法搭配概率约束。

这套Matlab代码解决的就是这样一个组合问题：常规机组动态经济调度 + 风电随机性多场景建模。它能给出每个时段各机组的启停计划、出力安排、旋转备用容量，以及整个调度周期内的期望总成本。适合电力系统方向的研究生做课题、工程师做初步方案论证，也适合想从静态模型过渡到动态模型的入门者作为模板去改。

我在下面会把模型怎么建、场景怎么生成、代码怎么组织、求解时哪些地方容易翻车，全部拆开讲一遍。代码本身是完整可跑的，但我更建议你先跟着文章把逻辑捋顺了再跑，不然出错的时候你根本分不清是数据问题、约束问题还是求解器参数问题。

2. 目标函数和约束条件：模型的骨架这样搭

2.1 目标函数：三部分成本各自怎么算

这套模型的目标函数包含三块：燃料成本、启停成本、弃风惩罚成本。缺一不可，但很多人容易漏掉第三块。

燃料成本用二次函数表示，经典形式是：

C_i(P_it) = a_i * P_it^2 + b_i * P_it + c_i

其中下标i是机组编号，t是时段编号。实际用求解器时，二次函数经常会拉低求解速度，特别是在大规模场景下。我的做法是把它做分段线性化处理，把连续出力范围切成若干段，每段用直线近似，精度损失很小，但intlinprog求解会顺畅很多。具体切段数可以按精度需求选，一般5到7段就够用。

启停成本相对直白：机组启动一次要付固定启动费用，停机一般假设不花钱（当然也有模型加停机费用的，本套代码没加，你可以自行扩展）。这里注意一个细节——启动成本在目标函数里要乘以一个“启动动作标志”，它是由机组启停状态变量推导出来的，后面代码部分会说清楚怎么推导。

弃风惩罚成本是我建议一定保留的。为什么不直接允许弃风？因为如果弃风没有代价，求解器会倾向于在每个时段都安排一定弃风来放松备用约束，得出一个看似成本很低、实际上浪费清洁能源的方案。加一个单位弃风惩罚系数（通常取较高值，比如300到500元/MWh），就能在“保安全”和“用风电”之间取得合理平衡。

2.2 关键约束：哪些必须写，哪些可以视情况简化

第一类约束是功率平衡约束，这是经济调度模型的基石：每个时段所有机组出力加上风电实际出力，必须等于该时段负荷。注意，因为是多场景建模，功率平衡不是对“期望场景”成立一次就够，而是对每个风电场景都要成立。这就是场景法模型规模膨胀的根源——场景数量越多，平衡约束的行数就越多。

第二类是机组出力上下限约束，简单但容易写错。写出力上限时，不能写P_i <= P_max，而是要写成P_i <= P_max * u_i，u_i是0/1启停变量。否则停机状态的机组也会被分配一个正出力，结果完全失真。下限同理。

第三类是爬坡约束，这是“动态”的核心。常规机组相邻时段出力变化量受爬坡速率限制，包括向上爬坡和向下爬坡。这里有个隐蔽的坑：如果机组当前时段停机（u_it=0），它下一时段开机，出力可能直接从0跳到某个大值，这个“启动过程”的出力变化是否受爬坡约束限制？严格来说，启动过程也有一个启动爬坡速率限制，而且通常比正常运行时的爬坡速率更严格。有些论文为了建模简洁会忽略启动爬坡约束，但工程上不建议忽略。本套代码里已经把启动/停机过程的爬坡约束一并写进去了。

第四类是旋转备用约束。这是应对风电随机性的关键安全网。系统需要在每个时段保留足够的向上备用容量，以便在风电出力低于预期时，常规机组能快速顶上。向上备用约束的写法是：

sum_i (P_i_max * u_i - P_i) >= 风电预测出力 * 备用系数 + 负荷备用需求

注意左边是“还能往上加多少出力”的总量，这是动态的，随机组实际出力变化。很多初学者把它写成固定值备用，那就完全丧失了备用约束的意义。

2.3 风电随机性如何进入约束：场景耦合与置信度的取舍

场景法里，风电出力是一个在多个场景间变化的量。每个场景s都有对应的风电出力序列W_s(t)。那么功率平衡、备用约束就要写成带场景下标的形式。问题来了：机组启停和出力变量要不要也带场景下标？

这里有两种建模流派。一种称为“非预期性约束”，要求机组决策只能基于当前可用信息，不能预知未来场景，所以所有决策变量对所有场景必须相同，共担一个解。这种做法的模型更符合实际调度逻辑，但约束形式复杂，需要额外加非预期性约束。另一种做法是让决策变量随场景变化，每个场景单独求解，最后按场景概率加权求期望成本——这其实就是多场景分别求解，模型简单，但结果在物理上不一定可执行，因为调度员在实时运行时并不知道未来哪个场景会真实发生。

本套代码采用了一个折中方案：机组启停变量在所有场景下保持一致（因为启停需要提前确定，不可能等风电场景揭晓后再决定），但机组出力变量允许随场景调整，因为实际运行中，调度员可以根据风电实时出力对常规机组做再调度。这个假设在实际工程里是合理的，也是大多数随机调度论文采用的处理方式。这样约束规模不会过度膨胀，同时保留了随机性对调度方案的实质性影响。

2.4 为什么选MILP而非智能算法

这个问题我几乎每隔一段时间就会被人问到。很多做课题的同学第一反应是用遗传算法、粒子群算法去求解，理由是“MILP建模太难了”。我的观点很明确：只要问题能写成MILP，就优先用MILP求解器。

原因有三。第一，MILP有全局最优性保证，intlinprog返回的解要么是最优解，要么给出了最优性间隙，你能知道自己离最优解有多远。智能算法跑出来的结果没有这个保证，很多时候你根本不知道结果靠不靠谱。第二，MILP求解的稳定性强，同样的模型和参数，不管跑多少次，结果都一样。粒子群算法每次结果都可能不同，写论文时审稿人问一句重复性，你很难答。第三，现代MILP求解器（包括Matlab内置的intlinprog）在分支定界和割平面技术上已经非常成熟，几千个变量、几千个约束的问题，一般分钟内能解出来。对学术研究和方案验证来说，这种速度和确定性远比“算法听起来高级”更重要。

当然也有例外：如果你的问题规模极大（比如上百台机组、上千个节点、上千个场景），intlinprog会比较吃力，这时可以转用商业求解器如Gurobi或CPLEX，再不行才考虑启发式算法。但那是另一个量级的工程优化问题了，本套代码的规模完全不需要走到那一步。

3. 风电随机性怎么建模：场景生成是重头戏

3.1 风速分布与出力转换：从Weibull到风电功率曲线

风电随机性建模的第一步是模拟风速。工程上最常用的风速分布是两参数Weibull分布，其概率密度函数为：

f(v) = (k/λ) * (v/λ)^(k-1) * exp(-(v/λ)^k)

k是形状参数，λ是尺度参数。参数怎么取？最简单的方法是根据风电场历史数据的平均风速和方差反推，近似公式是λ ≈ 平均风速 / Γ(1+1/k)，k通常在1.8到2.3之间。比如平均风速6.5m/s、k取2.0，λ大约在7.3左右。如果手头没有实际数据，用这些典型值也完全够跑通模型。

有了风速场景后，需要通过风机功率曲线把风速转换为出力。标准的转换规则是三段式：

• 风速小于切入风速Vin或大于切出风速Vout时，出力为0；
• 风速在额定风速Vr和切出风速Vout之间时，出力等于额定功率Pr；
• 风速在切入风速和额定风速之间时，出力通常按三次方关系近似：P = Pr * (v - Vin) / (Vr - Vin) 的三次方再乘以额定功率，也可以按风机厂家给定的功率曲线线性插值。

本套代码里用的是三次方近似，并对出力做了0到Pr的截断。这一部分很快，但很关键，因为后续所有的场景削减和约束构建都依赖这个转换结果。

3.2 蒙特卡洛抽样：生成原始场景集

风速场景的抽样分两步。第一步，用Weibull分布抽样得到基础风速序列，代表预测值；第二步，在基础风速上叠加预测误差。预测误差一般假设服从正态分布，标准差取预测风速的一定比例，例如10%到15%，并且误差可以在时间上做相关性处理——相邻时段的风速误差往往正相关，通过一个简单的一阶自回归模型可以体现。

我建议抽样的场景数初始设为500到1000个。有人可能会问，500个场景对应的约束行数太多怎么办？别急，场景削减就是用来解决这个问题的。

蒙特卡洛抽样的具体实现，在Matlab里无非是两类代码：用makedist和random函数从Weibull分布抽样，再加一个normrnd叠加误差。这里有一个容易踩的坑：抽样得到的风速可能是负值，必须在转功率前做max(0, v)处理，否则功率曲线公式里会出现负数的三次方，结果直接出错。

3.3 场景削减：500个场景变成10个，代价是什么

场景削减的核心思想是：把相似的场景合并，用少量代表性场景近似整个场景集的概率分布。削减后每个场景赋予一个概率值，所有概率之和等于1。这个近似是随机的代价，也是工程上权衡计算量和精度之后的必然选择。

常用的削减方法有两种，一是基于同步回代消除的快速前向选择法，二是K-means聚类法。我的经验是：小规模问题（10个场景以内）用K-means足够，速度快，实现简单；大规模问题用同步回代消除法更优，因为它考虑到了场景两两之间的距离和概率，削减结果在分布意义上更接近原始集合。

本套代码用的是同步回代消除法，具体逻辑如下：

1. 计算所有场景两两之间的欧氏距离；
2. 找到距离最近的一对场景，删除其中概率较小的那个，把它的概率累加到保留的场景上；
3. 重复上述过程，直到场景数达到预设目标数（比如10个）。

这段逻辑用Matlab写大概是二三十行，但要注意距离矩阵的计算要基于整个时间序列的出力曲线，而不是单点。很多初学者在这一点上做错，只算同一时段的风电出力距离，忽略了时间维度上的相关性，削减出来的场景集质量会很差。

削减到什么程度合适？经验值是：10到20个场景能保留原始分布90%以上的信息量，继续增加场景数对结果的改善有限，但求解时间会成倍增长。本套代码默认削减到10个场景，后续计算验证也证明了这组参数在精度和效率之间是比较合适的折中。

3.4 场景概率归一化：一个容易忽略的最后一步

削减完成后，每个保留场景都有一个概率，但这些概率之和未必正好等于1，因为削减过程中存在击穿（某些场景概率很小被直接删除）。最后一定要做一次归一化：每个场景的概率除以所有场景概率之和。

这一步不做的话，目标函数里期望成本的计算会系统性偏低或偏高，而且约束中的概率约束（比如旋转备用满足概率不小于某个阈值）也会失真。代码里我专门留了一行注释提醒这一点，因为我自己第一次写这块函数时也漏过——当时结果差得不多，但论文里的数字对不上，排查了大半天才发现是概率没归一化。

4. Matlab代码实现：数据怎么组织，求解器怎么调

4.1 数据结构设计：参数、变量、约束矩阵的分区管理

写Matlab代码求解MILP模型，最大的难点不是数学公式本身，而是把数学模型逐行翻译成intlinprog能接受的大矩阵Aeq、beq、A、b。我的习惯是把整个代码按功能分区：

• 数据准备区：负荷序列、机组参数、风速参数、场景生成与削减，全部在这里完成，输出是净负荷序列和每个场景的风电出力矩阵；
• 变量定义区：明确每个决策变量在总决策向量里的起始位置和长度，用一组索引变量（如idxP, idxU, idxStart等）去管理；
• 约束构建区：逐类构建约束行，每构建一类就立即写入Aeq/A矩阵，同时把对应的右侧项写入beq/b；
• 求解与结果输出区：调用intlinprog，然后把最优解拆分回各个变量，重新组织成可读的表格和图表。

这样的分区方式有几个直接好处。第一，调试定位快——某类约束出错了，直接跳到对应分区检查。第二，扩展容易——想增加一个约束类型，只需在约束构建区加一段代码，不用动其他部分。第三，把决策变量的索引管理集中在一起，后面对解的拆分完全不需要再猜。

决策变量的维度设计也要提前想清楚。本套代码的变量个数是：

机组数N_G × 时段数T × (出力变量1个 + 启停变量1个 + 启动动作变量1个) 再加上备用相关的松弛变量。

以10台机组、24时段为例，总变量数约720个，其中整数变量约480个。这个规模对intlinprog来说相当轻松，一般几秒到几十秒就能收敛到最优性间隙1%以内。

4.2 模型参数设定与常规机组数据

我推荐用两种典型算例来验证模型：3机6节点系统（适合快速验证逻辑）和10机39节点系统（适合正式课题分析）。套代码里两组数据都内置了，切换方式就是改一个N_G参数和对应的机组参数矩阵。

以一个典型的10机系统为例，机组参数含义如下：

• 机组编号
• 出力下限（MW）
• 出力上限（MW）
• 爬坡速率（MW/h）
• 二次燃料成本系数a、b、c
• 启动成本（元/次）
• 初始状态（0或1，表示第一个时段前是否已在线）

这里提醒一个细节：初始状态非常重要。如果第一时段某机组处于运行状态，那么它的最小出力下界约束在第一时段必须满足P_i1 >= P_i_min；如果处于停机状态，那么该时段出力必须为0而且不允许启动动作（否则启停变量逻辑就冲突了）。代码里的机组数据表第一行专门用来标注初始状态，很多人替换数据时容易漏掉这一行，导致第一时段的解完全不正常。

负荷数据建议用典型负荷曲线，早晚高峰和深夜低谷要体现出来——如果没有明显波动，爬坡约束和启停决策根本不会被激活，模型的“动态”特征也就体现不出来。24时段负荷曲线在本套代码里预设了一组数据，也可以读取外部Excel替换。

4.3 intlinprog的参数设置：求解质量与速度的平衡

intlinprog的调用并不复杂，核心是Options设置。我实际测试下来，这几个参数影响最大：

• IntegerTolerance：整数容忍度，默认1e-5，一般不用改；
• ConstraintTolerance：约束容忍度，默认1e-6，求解前可以放宽到1e-4，能显著减少数值病态问题，尤其是当约束矩阵里出现量级差异很大的系数时；
• RelativeGapTolerance：相对最优性间隙，默认是0（必须达到最优才停止），建议设为0.01（1%间隙）。对调度问题来说，1%的间隙代价误差可能只有几百元，完全可接受，但求解时间能降一个量级；
• MaxTime：限制最大求解时间，防止某些病态实例把脚本卡死。

我的建议是：先按默认参数跑一遍，确认模型逻辑无误后，再把RelativeGapTolerance设置为0.01进行后续批量实验。追求严格的0间隙在工程上没必要，还会拖长仿真时间。

还有一个非常实用的小技巧：给intlinprog传入一个合理的初始可行解。Matlab的intlinprog支持通过x0参数提供一个启动点，能极大加快分支定界早期的上下界收敛速度。在本套代码里，我是先按确定性模型（取风电预测期望值）跑一遍，把解作为随机模型的初始解传入，实测整体求解时间能下降30%到50%。

4.4 约束矩阵的构建细节：稀疏矩阵是必选项

约束矩阵必须用稀疏矩阵存储，这是我一直强调的。原因很直接：3机系统可能无所谓，但10机24时段10场景的模型，约束行数和变量列数都是千级别，如果用全矩阵存储，内存占用会达到几十MB甚至更多，而实际非零元素占比可能不到5%。Matlab的稀疏矩阵可以压缩存储这部分内存，更关键的是，intlinprog内部算法对稀疏结构有专门优化，求解速度快得多。

构建稀疏矩阵时，我的写法是先用三个数组（rowIdx, colIdx, valueIdx）收集所有非零元素的位置和值，最后一次性调用sparse函数生成矩阵。这样避免了在循环里不断拼接稀疏矩阵——那是Matlab性能的大忌。

约束构建的具体顺序也有讲究。我会按这个顺序来：

1. 功率平衡约束（行数 = T × S，S为场景数）；
2. 机组出力上下限约束（T × S × N_G个有效行）；
3. 爬坡约束（T-1 × S × N_G）；
4. 旋转备用约束（T × S）；
5. 启停状态逻辑约束（T × N_G）；
6. 启动动作推导约束（T × N_G）。

约束顺序不影响求解结果，但分组构建方便调试——哪一类约束出了问题，检查对应数据块的维度是否匹配就一目了然。

4.5 结果还原与输出：从解向量回到可读的调度表

intlinprog返回的解x是一个一维向量，必须按照之前定义好的变量索引拆分。我的做法是定义一个结构体solution，把每个变量段单独放好：

solution.P：N_G × T 的连续出力矩阵； solution.U：N_G × T 的0/1启停矩阵； solution.V：N_G × T 的启动动作标志位； solution.Wcurtail：T × S 的弃风量矩阵。

拆分后输出四类结果：每个时段的总发电成本、各机组出力曲线、机组启停计划表、风电消纳情况。我建议至少输出这三张图：

• 各机组出力叠加图 + 负荷曲线对照，能直观看出功率平衡是否满足；
• 机组启停状态甘特图，用stair函数画就能实现；
• 实际风电消纳量对比预测值的曲线，能看出弃风发生的时间和幅度。

有个结果检查的小技巧：把P矩阵每行求和加上风电消纳，对照负荷曲线，如果有个别时段的平衡误差超过0.1MW，多半是约束构建时分场景下标写串了。95%的模型bug都能通过这个“平衡校验”检查出来。

5. 算例验证：确定性模型和随机模型的结果对比

5.1 算例设置：2种规模、3组对照实验

我用两套算例来验证模型。第一套是3机6节点系统，T=24时段，场景数从5到20做收敛性测试；第二套是10机39节点系统，T=24时段，场景数固定为10。对照实验设置三组：

• 算例A：不考虑风电随机性，用风电预测期望值做确定性调度；
• 算例B：考虑随机性，但旋转备用需求按固定值设定；
• 算例C：完整的随机动态经济调度模型，含动态旋转备用约束和多场景功率平衡。

这三组对照能清楚地区分出“随机性”和“动态备用”各自对调度方案和成本的影响。很多论文只给一个孤立的结果，读者很难判断模型到底贡献在哪——三组对照才是真正有说服力的做法。

5.2 结果观察1：常规机组出力计划的变化

先看3机系统下三种算例的调度结果。确定性模型（算例A）给出的调度方案在负荷低谷时段会把部分机组压到最低出力，在傍晚高峰时段启用额外机组。把风电预测值替换成实际场景后，会发现一个关键问题：确定性模型在风电出力低于预测的时段出现了备用容量不足，最严重的一个时段，系统向上备用缺口达到45MW——这意味着如果风电真按那个低场景发展，系统将无法保证调频能力。

而随机模型（算例C）在同样的时段里预留了更充足的向上备用，代价是某些机组出力整体上调。成本对比如下：

看到随机模型总成本比确定性模型高约1.4万元，涨幅2.3%左右。这1.4万就是“为了应对不确定性而付出的安全成本”。如果你算出来随机模型成本反而更低，那几乎可以断定模型有问题——大概率是可再生能源出力被不合理放大，或者备用约束写得太松。

5.3 结果观察2：备用容量的时空分布差异

动态备用约束下，系统向上备用容量不再是一个常数，而是随时间变化。打开结果看每个时段的备用值，明显能看到：夜间负荷低谷、风电预测出力高的时段，备用值相对较小；白天负荷爬升、风电出力不确定度大的时段，备用值显著增大。这说明模型自动在“经济性”和“安全性”之间做了时空维度的权衡——这正是动态备用相对固定备用的核心区别。

对比固定备用的算例B，它的备用约束在每个时段都要求同样的裕度，结果是在风电出力确定高的深夜时段浪费了备用容量，而在风电预测误差大的白天时段，备用却不够。这就是固定备用模型“顾此失彼”的典型表现。场景法与动态备用结合后，备用需求实际上由“所有场景下最恶劣的备用缺口”决定，相当于系统自动识别关键场景、关键时段，并针对性地配置安全裕度。

5.4 结果观察3：场景数量对解的影响

我在3机系统上把场景数从5逐步加到50，观察总成本和求解时间的变化趋势：

• 5个场景：总成本62.21万元，求解时间约3秒；
• 10个场景：总成本62.55万元，求解时间约12秒；
• 20个场景：总成本62.67万元，求解时间约45秒；
• 40个场景：总成本62.71万元，求解时间约180秒；

可以看到，从10个场景增加到40个场景，总成本变化不到0.3%，但求解时间增加了15倍。这说明对这套数据而言，10个场景基本已经足够。继续增加场景数属于典型的边际效益递减。如果你在做自己的课题，建议做一个类似的场景数-成本收敛曲线，一方面能找到精度和速度的平衡点，另一方面这本身也是论文中一个很有说服力的数据支撑。

6. 实测中的五个坑：代码跑通之后才算真正的挑战

6.1 爬坡约束和启停状态之间的逻辑冲突

这是我在调试过程中遇到最多的一类问题，几乎每次改数据都会碰到。场景是这样的：某台机组在第一时段处于停机状态（u_i1=0），第二时段要启动并出力到200MW。这时爬坡约束里有一项P_i2 - P_i1 <= 爬坡上限，而P_i1=0，所以理论上只要P_i2不超过爬坡上限就满足。但前提是P_i1必须真的为0——如果机组初始状态设为运行但出力下限没设定，或者约束矩阵里对应的行没有乘启停变量，P_i1就可能被分配一个非零值，爬坡约束就形同虚设。

我的处理方式是：在约束构建区单独做一次自检，检查所有“启停状态为0”的时段，对应出力变量的解向量是否严格为0。如果出现非零出力，基本可以断定上下限约束漏乘了u_i变量，排查方向就很明确了。

6.2 备用约束里的“隐式耦合”很容易被忽视

旋转备用约束表面上只涉及机组出力上限和实际出力的差值，但“实际出力”本身受爬坡约束限制。也就是说，如果某个机组当前时段被要求出力接近上限，它下一时段想再往上爬就爬不动了，备用能力实际上被削弱了。这在约束构建时很容易被忽视——备用约束和爬坡约束是分开写的，但物理上它们通过“出力值”这一共同变量耦合在一起。

解决方案是：在旋转备用约束里，把“可上调能力”计算为min(机组出力上限, 上一时段出力 + 爬坡上限) 再减去当前出力，而不是简单用上限减出力。这样模型自动考虑了爬坡对备用能力的影响。否则，备用裕度会被高估，在风电突然下跌时可能出现实际响应能力不足的隐患。我在这套代码里已经按这个修正方式实现，你后续扩展时也要记得加上这一项。

6.3 intlinprog求解时间的突然恶化

有时候你会遇到这种情况：改了一个约束，求解时间从10秒暴涨到300秒还停不下来。常见原因是某个约束写得太“硬”，把可行域推到一个非常狭窄的角落，分支定界需要大量探索才能证明最优。

我的排查方法有两个。第一，把RelativeGapTolerance临时设成0.2，先快速拿到一个可行解，看看这个解是不是合理的；如果连可行解都拿不到，基本是约束之间矛盾导致无解。第二，用MaxTime限制求解时间，同时把求解过程日志打开（Display选项设为iter），观察MIP gap的下降趋势——如果gap一直卡在某个值不肯下降，通常是某个约束的非零系数存在数值问题，尝试把该约束的系数做归一化处理（比如把单位从MW改成GW），有时能立刻见效。

6.4 场景削减结果不稳定

同步回代消除法在实现时有个容易踩的细节：距离矩阵的计算必须基于“归一化后的出力序列”，否则风电场容量大的场景在距离计算中天然会“权重过高”，导致削减出的场景集合总是偏向高出力场景，削弱了代表低出力场景的能力。

另一个细节是随机种子的影响。蒙特卡洛抽样和削减过程如果涉及随机数生成，不同种子可能导致最终削减结果不同。这不是bug，但会让你的实验结果难以复现。建议在代码开头固定随机数种子，例如rng(2025)，论文里加一句“所有实验采用固定随机种子以保持可复现性”，既规范又省事。

6.5 结果表格里的数字单位混淆

最后这个坑很基础，但杀伤力非常大：单位不统一导致整张结果表报废。Matlab里如果机组容量单位是MW，成本系数单位是元/MW²和元/MW，目标函数里的量纲混起来后就很容易出现千、百万级别的偏差。我当时遇到过的情况是：缺了一项系数除以1000的转换，燃料成本多算了三倍，结论完全错误。

所以无论用什么数据，第一步就是把所有参数的物理单位列成一张表，逐项确认。建议在代码里统一使用“MW”、“元”作为基本单位，所有输入参数在数据准备区完成一次性的单位换算，后续模型里绝不再出现换算因子。

7. 这套代码还能往哪个方向扩展

模型本身改完之后，我觉得最有价值的事，是沿着“随机动态经济调度”这条线继续往三个方向延伸。第一个方向是加入储能。储能系统本质上是一个跨时段耦合元件，它的充放电决策天然适合并入动态经济调度框架——只需增加储能SOC状态变量和充放电功率变量，再补充SOC递推约束和容量约束。场景法框架完全兼容储能建模，改动的代码量不大，但能显著提升系统对风电不确定性的消纳能力。

第二个方向是改为滚动时域优化。把24时段一次求解改成每4小时滚动一次、每次求解未来8小时的窗口，这样可以把最新风电预测数据持续引入模型，决策更贴近实时变化。实现上只需要把数据准备区的负荷序列和风电场景按窗口切片，循环调用求解器即可。滚动优化的效果在风电预测误差随时间变化时体现尤为明显。

第三个方向是多目标化。目前模型只优化经济成本，你可以引入碳排放量作为第二个目标函数，用加权和法或ε-约束法求Pareto前沿。这个方向在“双碳”背景下论文价值很高，代码层面只需增加一个碳排放在目标函数里的线性项，再对不同权重做几次参数扫描就能出结果。

如果你只是临时跑通一个课程作业，那么本文的模型和代码已经足够；如果你想拿这套模型支撑一篇高质量论文或者做实际项目预研，建议花时间把上面三个扩展方向至少做通一个。这个模型框架的终点不在“跑通”，而在“怎么把它变成一个可以回答真实工程问题的工具”。我自己当初就是在这套模型上一步步加上储能和滚动优化，最后做出了一套可以适应不同风电场数据、不同负荷曲线、不同机组参数的通用调度分析工具——那种感觉比单纯调通一个函数要爽得多。