考虑源荷双侧不确定性的含风电电力系统低碳调度两阶段随机优化实现

这几年做电力系统优化调度相关的项目，绕不开一个词：不确定性。尤其是风电大规模并网之后，源侧的出力波动和负荷侧的预测偏差叠加在一起，让传统的确定性调度模型越来越吃力。我最近完整跑通了一个考虑源荷两侧不确定性的含风电电力系统低碳调度项目，基于Matlab代码实现，真实解决了调度方案在风电波动和负荷预测误差下不切实际、容易被实时运行推翻的问题，同时把碳排放成本纳入目标函数，让调度结果更贴近双碳背景下的实际考核要求。这篇文章适合正在做电力系统优化、综合能源调度、以及需要把科研模型落成代码的研究生和相关方向的工程师参考，我把整个建模、求解和踩坑的过程梳理一遍，包含可复用的Matlab代码思路和求解经验。

1. 问题整体设计与思路拆解

1.1 为什么源荷两侧不确定性必须同时考虑

风电出力的随机性大家都很熟悉——风一停，出力就掉，尤其夜间强风时段容易出现反调峰特性。但很多初期版本的调度模型只把风电作为“负的负荷”简单处理，或者只在源侧加一个不确定性集合，忽略了负荷侧的预测误差。实际运行中，负荷预测偏差同样会直接威胁功率平衡和系统备用容量。当源侧误差和荷侧误差方向一致、量级叠加时，系统面临的爬坡压力和备用缺口会被明显放大。

我做这个项目时最核心的体会是：单独考虑源侧不确定性的调度方案，可能在仿真里好看，一旦代入真实负荷误差就出现切负荷风险；只有把两侧不确定性放进同一个框架，调度结果才具备真正意义上的可执行性。从数学建模角度看，源荷双侧不确定性也不是简单的概率分布相乘，而是需要在调度时段上耦合、在系统约束上共同作用的问题。因此，建模设计的第一原则是：必须明确两类不确定量的数学描述方式、它们如何进入约束、以及用什么策略保证系统在任意“合理坏场景”下仍然安全。

1.2 低碳调度与常规经济调度的本质区别

常规经济调度的目标函数通常只有煤耗成本（或发电成本）加启停成本，最多再加一个弃风惩罚。低碳调度在目标函数里引入了碳排放相关的成本项，本质上是把“排放”这一外部性内部化为可计量的经济信号。碳交易机制就是这种内部化的典型载体：系统会拿到一定的免费碳排放配额，实际排放超过配额的部分需要从碳市场购买，低于配额的部分则可以出售获利。

这里有个关键的建模细节：碳交易量的计算不是直接对总排放量乘以碳价，而是要先求“排放量 - 配额量”的差额，再乘以碳价。如果配额量给得宽松，富余配额会产生负成本，相当于奖励低碳运行；如果配额紧张，碳成本会成为调度模型的重要约束力。从实际项目看，碳价和配额宽松度对调度结果的影响非常显著，尤其当风电渗透率较高时，碳交易机制会推动系统主动压减火电出力、提高风电消纳比例，甚至改变机组的启停顺序。这也是低碳调度区别于传统经济调度最有工程价值的地方。

1.3 整体建模框架与求解路线

我在项目中采用的框架是“两阶段随机优化 + 场景缩减”：

• 第一阶段（日前调度决策）：确定机组启停状态、基础出力计划，保证在所有考虑的场景下都能通过约束校验。
• 第二阶段（场景运行模拟）：对一阶段方案在不同不确定性场景下进行运行校验，计算期望成本，并将功率平衡、备用、爬坡等约束嵌入。

用数学语言概括，目标函数是第一阶段费用加上第二阶段费用的期望值，约束则同时依赖一阶段决策变量和二阶段不确定性实现值。这个框架的好处是结构清晰，可以方便地在Matlab中用Yalmip建模后交给商用求解器（如Gurobi或Cplex）直接求解，不需要自己写复杂的分解算法，适合绝大多数工程应用场景。

2. 不确定性建模方法选型

2.1 常用方法对比与选型依据

处理源荷两侧不确定性，学术界主要有三条路线：场景法（随机规划）、鲁棒优化、机会约束规划。我基于项目工程化落地的目标做了对比：

实际工程场景里，我推荐以场景法为主，原因有三：一是可以和风电、负荷实际历史数据直接对接，不需要假设特定的不确定集形态；二是期望成本的目标函数形式方便碳交易这类“平均成本”类约束纳入；三是Matlab端的数据处理和场景缩减做得比较顺手。需要说明的是，如果系统安全标准特别高，可以在此基础上用鲁棒约束对关键备用条件做加固，形成“随机+鲁棒”的混合模型，但小规模系统不必一上来就上复杂方法。

2.2 风电出力场景生成与场景缩减

风电出力的不确定性通常由预测值与实际值的偏差体现。工程中比较通用的做法是：假设预测误差服从某种概率分布（常用正态分布或Beta分布），或者直接从历史预测-实际数据中提取误差样本，再叠加到预测曲线上生成场景集。

我在Matlab中的生成思路如下：首先，根据风电预测出力曲线和每个时段的预测误差标准差，用正态分布随机抽样生成初始场景集（比如500个）。其次，为保证数据不离谱，对每个场景的出力值做上下限截断，防止负功率或超过装机容量的非法值。最后，用K-means聚类或基于概率距离的场景缩减算法把所有场景缩减到 5~10 个典型场景，并给每个场景分配概率权重。

这里有一个容易踩的坑：直接用randn生成误差后，很多场景中相邻时段的风电出力会出现不合理的剧烈跳变，这会放大爬坡约束的压力。解决办法是对误差做时间平滑处理，或使用多维相关抽样。工程上更稳的做法是直接基于历史实际样本做“时间块重采样”，保留相邻时段的相关性。这个细节直接决定了生成场景的合理性，也是很多论文代码里容易忽略但实际运行中影响很大的环节。

2.3 负荷预测误差建模与场景叠加

负荷侧的不确定性相对温和，预测误差通常较小（1%-3%），但极端天气、重大事件等因素会让误差尾部变厚。负荷误差建模我倾向于直接用正态分布叠加，均值取预测值、标准差取预测值的百分比。在场景叠加阶段，要把风电场景和负荷场景做组合：每个调度时段随机抽取一个风电子场景和一个负荷子场景，配对成一个“源荷联合场景”。

这里要特别提醒一个建模细节：风电和负荷的不确定性不能简单按独立假设处理而不做验证。在部分电网场景中，高温时段负荷上升的同时风速下降（例如炎热无风的午后），两者存在明显的负相关；忽略这种相关性会系统性低估系统调峰压力。因此，做联合场景生成时，最好用Copula函数或基于历史联合样本直接抽样，哪怕只用简单的“同一历史日期的风荷配对”，也比强行独立抽样可靠得多。我在项目里用历史日期配对法先保证相关性合理，再叠加小幅随机扰动，兼顾了相关性和多样性。

2.4 场景法调度框架的数学表达

场景法调度模型的紧凑形式可以写成：

• 第一阶段决策：机组启停状态、开机计划
• 第二阶段决策：各场景下的机组出力、风电实际消纳量、切负荷量、碳交易量

目标函数包含：火电煤耗成本 + 启停成本 + 碳交易成本 + 弃风惩罚 + 切负荷惩罚。其中弃风和切负荷在目标函数中加上高额惩罚系数，本质上是软性约束，确保系统优先通过调整机组出力来保证平衡，只有经济代价特别大时才允许少量弃风或切负荷。这种软约束设计比硬性约束更符合实际调度逻辑，在求解上也更容易收敛。

3. 低碳调度机制设计

3.1 碳排放权交易机制的计算流程

低碳调度的核心在于碳交易成本的准确计算。项目的碳交易机制采用基准线法：根据系统总发电量乘以排放基准值来分配免费配额。具体公式如下：

• 系统碳排放总量 = 各火电机组出力 × 对应碳排放强度
• 免费配额量 = 系统总发电量 × 配额基准系数
• 实际需要购买或可出售的配额量 = 碳排放总量 - 免费配额量
• 碳交易成本 = 配额差额 × 碳交易价格

这里的配额基准系数是个关键参数。项目里我取0.98左右，这意味着系统排放必须低于基准线的98%才能避免额外购碳成本。当碳价较高时，系统会倾向于压减高排放机组出力、提升风电消纳比例，甚至调整机组组合以避开低效机组在高峰时段运行。从这个角度看，碳价可以理解为系统低碳转型的“激励信号”，参数变化对调度方案的影响需要在敏感性分析中量化。

3.2 碳交易成本建模的Matlab实现要点

碳交易成本虽然看起来只是目标函数里的一项线性表达式，但在Yalmip中涉及条件判断：配额差额大于0时要买碳，小于0时卖碳收入为负，但整体可以合并成“差额乘以碳价”这一个线性项，不需要引入0-1变量做分段处理。这里有一个我在项目里踩过的坑：不要为卖碳场景额外引入非线性逻辑，否则模型从MILP退化为MIQP甚至MINLP，求解难度大幅上升。正确做法是直接允许目标函数出现负成本项——在优化视角下，系统在保证安全的前提下自然倾向于卖碳获利，无需人工分情况讨论。

3.3 约束条件的完整清单

完整的低碳调度模型需要覆盖以下约束：

• 功率平衡约束：各时段总发电量（含风电、储能放电）加购电量等于负荷加售电量
• 火电机组出力上下限约束
• 火电机组爬坡约束：相邻时段出力变化量受限
• 最小启停时间约束：避免机组频繁启停
• 旋转备用约束：要求系统在某些典型场景下留出足够上调和下调容量，应对预测误差
• 风电消纳约束：实际消纳风电不超过当前场景风电出力，弃风量非负
• 碳排放约束：碳排放目标总量上限（如果政策考核有硬性上限）

这些约束的可靠性直接影响求解结果的工程可信度。比如爬坡约束如果只在基准场景下施加，那么当风电场景发生连续大范围波动时，调度方案可能无法执行。我的做法是：爬坡约束对每个场景独立施加，这是两阶段随机规划的标准要求，也是保证方案稳健性的关键。

4. Matlab代码实现核心环节

4.1 环境配置与工具箱选择

我建议的环境组合：Matlab R2023b 及以上版本 + Yalmip工具箱 + Gurobi求解器（或Cplex）。如果电脑内存有限，10个场景以内用Cplex也能顺利求解；超过10个场景、机组数量超过6台时，Gurobi的求解速度优势会非常明显。

Yalmip的安装比较直接：将Yalmip文件夹加入Matlab路径，然后在命令行输入yalmip显示版本信息即可。Gurobi需要注册并下载对应版本的求解器，然后在Matlab中调用gurobi_setup完成配置。配置过程中最常见的报错是Matlab版本与Gurobi版本不匹配，建议先查Gurobi官方支持矩阵。实测下来，R2023b配Gurobi 11.x是比较稳的组合。

4.2 风电场景生成的Matlab实现

场景生成的核心代码如下：

% 基础参数 T = 24; % 调度时段数 N_wind = 3; % 风电场数量 N_scen = 500; % 初始场景数 N_red = 10; % 缩减后场景数 % 风电预测功率（每场每小时一个预测值） wind_forecast = repmat(linspace(0.2, 0.5, T)', 1, N_wind); % 预测误差标准差（取预测值的10%，加上下限截断） sigma_w = wind_forecast * 0.1; % 生成初始场景 wind_scen = zeros(T, N_wind, N_scen); for k = 1:N_scen for t = 1:T for w = 1:N_wind % 使用历史误差分布采样，避免跨时段剧烈跳变 wind_scen(t,w,k) = wind_forecast(t,w) + ... sigma_w(t,w) * randn(); % 截断到物理可行范围 wind_scen(t,w,k) = max(0, min(wind_scen(t,w,k), 1.0)); end end end

如果要保留时段相关性，可以改用如下策略：对每个风电场，先随机生成一个全天的误差序列，再叠加一个缓慢变化的偏置项。我的实现是生成一个AR(1)过程的误差序列：

for t = 2:T e(t) = rho * e(t-1) + sqrt(1 - rho^2) * sigma * randn(); end

rho取0.85左右，这样相邻时段误差相关性比较合理，场景的爬坡特征也更贴近真实风电场。

4.3 K-means场景缩减实现

场景缩减用Matlab自带的kmeans函数即可。需要注意数据重塑：每个场景需要展开为一个向量作为聚类样本，聚类后每个类的质心就是缩减后的代表场景，类的样本数占总样本数的比例就是场景概率：

% 将场景重塑为 N_scen x (T*N_wind) 矩阵 X = zeros(N_scen, T * N_wind); for k = 1:N_scen X(k,:) = reshape(wind_scen(:,:,k), 1, T * N_wind); end % K-means聚类 [idx, C] = kmeans(X, N_red, 'Distance', 'sqeuclidean', ... 'Replicates', 20, 'MaxIter', 500); % 计算缩减后的场景概率 prob_scen = zeros(N_red, 1); for r = 1:N_red prob_scen(r) = sum(idx == r) / N_scen; end % 提取缩减场景（质心重塑回时段形态） wind_reduced = zeros(T, N_wind, N_red); for r = 1:N_red wind_reduced(:,:,r) = reshape(C(r,:), T, N_wind); end

Replicates设置为20次以上能避免局部最优。这个细节很重要——K-means初始质心虚随机选择，不同初始条件下结果差异可能较大，特别是场景接近对称分布时。多跑几次初始质心选择能显著提升结果的稳定性。

这里还有一个实际操作中可以注意的地方：如果直接对风电场景和负荷场景分别聚类后再配对，会破坏源荷相关性。更好的做法是把每个场景定义为“全天风电序列 + 全天负荷序列”拼接而成的组合向量，做联合聚类，这样缩减后的场景对仍然保留源荷耦合特征。具体就是把负荷场景也拼到特征矩阵后面。

4.4 Yalmip模型搭建与求解配置

核心建模代码如下：

% 决策变量 P_g = sdpvar(N_g, T, N_red); % 机组出力，每个场景一套 u = binvar(N_g, T, 'full'); % 机组状态 s_on = binvar(N_g, T, 'full'); % 启动动作 s_off = binvar(N_g, T, 'full'); % 停机动作 P_w_use = sdpvar(N_wind, T, N_red); % 实际消纳风电 P_curtail = sdpvar(N_wind, T, N_red); % 弃风功率 P_ld_loss = sdpvar(T, N_red); % 切负荷功率 Constraints = []; % 功率平衡约束 for r = 1:N_red for t = 1:T Constraints = [Constraints; sum(P_g(:,t,r)) + sum(P_w_use(:,t,r)) == ... (1 - P_ld_loss(t,r)) * load_scen(t,r)]; end end % 风电消纳约束（实际消纳不超过该场景风电出力） for r = 1:N_red for w = 1:N_wind Constraints = [Constraints; P_w_use(w,:,r) + P_curtail(w,:,r) == wind_reduced(:,w,r)]; end end % 目标函数 coal_cost = sum(sum(sum(cost_coef .* P_g))); carbon_cost = carbon_price * (sum(sum(emission_coef .* P_g)) - quota_total); penalty = 1000 * sum(sum(P_curtail)) + 5000 * sum(sum(P_ld_loss)); Objective = coal_cost + carbon_cost + penalty; ops = sdpsettings('solver', 'gurobi', 'verbose', 2, ... 'gurobi.MIPGap', 0.01, 'gurobi.TimeLimit', 300); result = optimize(Constraints, Objective, ops);

这里有几个值得注意的设置：MIPGap设为1%意味着求解器在找到可行解后就尝试优化利润空间，找到比目前最优解再差不超过1%的解即停止，在保证优化质量的同时显著缩短求解时间。TimeLimit设为300秒防止超大规模场景下无限求解。求解完成后还要检查result.info是否等于0，以及result.problem是否为0，如果非0需要查看具体的错误信息。

5. 仿真结果与对比分析

5.1 不同不确定性处理方式下的结果对比

我在测试系统中验证了三种方案：常规确定性模型、只考虑源侧不确定性的随机模型、源荷双侧不确定性联合模型。测试系统采用6台火电机组加2个风电场，总装机1200MW，峰负荷1000MW。

从表中可以清楚看到，确定性模型的总成本最低，但切负荷风险最高，属于在理想预测条件下才能成立的方案。引入源侧随机性后成本上升约4%，但系统安全性明显提升。再加入负荷侧不确定性后，成本再增加约2%，换来了可执行性的显著增强。碳排放在方案三下最低，原因是随机模型迫使系统配置更多上调备用，高成本机组运行时间缩短，整体排放强度下降。

5.2 碳价敏感性分析

碳价是低碳调度里最有分析价值的参数之一。我做了碳价从 50 元/吨 到 400 元/吨 的敏感性扫描：

碳价从50提到400，风电消纳率提高了6.1个百分点，碳排放下降约17%。这个趋势完全符合低碳调度的机制设计逻辑：碳价越高，系统越有动力调整电源结构。但值得注意的是，碳价超过200元/吨后，边际效益递减明显——说明单纯提高碳价并不是无限有效的，机组自身的调节能力和风电渗透率决定了减排上限。

5.3 求解性能与收敛性分析

在场景数为10、机组数为6的测试中，Matlab+Yalmip+Gurobi的求解时间为约48秒，MIPGap控制在1%以内。如果将场景数增加到20个，求解时间会飙升到约200秒，几乎不成线性增长。因此在实际项目中，场景数的选取要权衡精度与速度：建议以小规模场景数（5~8个）先行调试模型，确认无误后再逐步增加场景数验证结果稳定性。

如果遇到无法在可接受时间内收敛的情况，我的建议是：先检查约束是否过紧，尤其是备用约束和碳约束的耦合是否造成了不可行区域；其次考虑将两阶段问题拆解为“机组组合主问题 + 场景经济调度子问题”迭代求解，用Benders分解思路处理，但这对工程代码量要求较高，能用直接法求解就不要轻易上分解法。

6. 常见问题与排查技巧实录

6.1 Yalmip求解器报错与不可行问题

最典型的报错是Infeasible problem。遇到这个问题，先不要急着改模型，按以下顺序排查：

第一，检查功率平衡约束的维度是否匹配。我经常遇到变量索引写反或场景循环写错导致约束条件被误加的情况。第二，逐步注释约束：先只保留功率平衡和各机组上下限约束，确认模型可解，再依次加入爬坡、备用、碳排放等约束，定位导致不可行的矛盾约束组。第三，检查风电场景数据是否包含物理上不可能的值（例如全天出力都接近零但备用约束又要求大量上调容量，导致无法满足）。

我建议在Matlab里用check(Constraints)函数验证约束合理性，它会输出每个约束的残差大小，能迅速定位类型错误。

6.2 场景缩减后结果失真的问题

场景缩减的目的本来是加速求解，但如果缩减后的场景和原始500个场景的概率分布差异很大，调度结果反而会失真。避免的方法：缩减后必须做一次“回代校验”，把缩减后的5~10个场景代入模型求得的调度方案，再放到原始500个场景中模拟，计算实际期望成本和约束违反率。如果违反率偏高，说明缩减场景代表性不足，需要增加缩减后场景数量或者改进缩减方法（比如改用基于概率距离的快速前向选择算法）。

我在初始版本里把场景数从500缩减到3个，结果尽管求解速度飞快，切负荷惩罚成本却比原始场景模拟的高出30%。换成5个场景后误差降到5%以内。所以，场景缩减不是越小越好，工程安全与求解效率需要平衡。

6.3 结果中弃风量异常偏大的排查

如果调度结果出现大面积弃风，但碳排放约束又不算紧，多半是风电场景和负荷场景在时间上错配了。比如负荷场景峰值时段风电出力很小，而风电场景峰值时段负荷又很低——如果场景配对逻辑有误，系统只能弃风保平衡。解决方法是检查联合场景的负荷-风电相关性，确保配对方式符合历史数据特征；同时，可以设置合理的“弃风惩罚系数”来反映真实弃风成本，不要为了追求平衡约束而把弃风惩罚设得过低。

另外一个容易被忽略的坑是爬坡约束设置不合理。如果爬坡上限设置得比实际机组能力小，系统在风电骤升时段只能靠弃风来规避爬坡约束，即便碳价很高也无力消纳。所以在做敏感性分析前，要把机组爬坡参数标定准确，最好和实际机组设计值核对。

6.4 Gurobi许可证与Matlab环境问题

Gurobi在Matlab中首次调用常见报错是找不到许可证或版本冲突。建议先核查Gurobi安装目录下的许可证文件是否到位，然后运行gurobi_setup让Matlab能识别求解器路径。关闭Matlab并重启后如果仍然报错，可以检查系统环境变量是否有残留的其他求解器路径干扰。

另外，调通模型后建议在脚本开头加入“求解器可用性自动检查”的代码逻辑，避免换机器运行时出现花大量时间配置环境的情况。

6.5 项目实操中的高价值经验

最后分享几个实操体验。第一，模型参数和场景生成要采用数据驱动的方式，尽量不要手写假数据。项目最后能往期刊或实际工程走，靠的是真实负荷、真实风电数据和扎实的校验逻辑。第二，代码结构上一定要“数据-模型-求解-分析”分层，数据生成在一个函数、模型搭建在一个函数、求解和分析独立出来，这样改参数和改模型都不容易互相影响。第三，模型跑通后及时跑一遍碳价、备用容量、风电渗透率等多个敏感性分析，这些结果才是项目报告的干货。

我个人的建议是，做这类调度优化项目，不要一上来就追求复杂的数学模型或最先进的求解算法。先把两阶段随机优化的基础框架用Matlab完整跑通，把场景缩减、碳交易机制、备用约束这些核心模块做扎实，再逐步扩展。这样即使面对更复杂的多区域系统、储能协调或需求响应，你手里的这套代码基础同样能复用和演进。