基于场景法的含风电低碳调度源荷不确定性建模与求解

1. 为什么源荷两侧不确定性必须放在一个模型里

如果你正在做电力系统优化调度，尤其是含风电的低碳调度，那么“源荷两侧不确定性”这几个字一定会出现在开题报告或者项目需求里。这个题目看起来不大，但真正动手用Matlab实现一遍后你会发现，90%的工作量在建模和场景处理上，剩下那10%才是写求解器约束。这篇分享我就把整个实现过程拆开讲，包括模型怎么建、场景怎么生成和削减、YALMIP里怎么写约束，以及实际调参时踩过的一些坑。

先说这个题目的核心矛盾。传统调度里，机组出力可以按确定性的负荷预测来安排，顶多留一个固定旋转备用率。但一旦加入大规模风电，源头出力变得不确定，负荷预测也不可能完全准，两个不确定量叠在一起，原来的“确定”调度就变得非常脆弱。可能风功率预测值明明够用，实际风小了，系统就要拉闸限电；也可能负荷预测偏低，实际用电猛涨，机组爬坡跟不上，频率直接跌落。所以我一直建议做含风电调度项目的人，不要只在风电侧做文章，负荷侧的不确定性必须一起进模型，否则算出来的结果到实际运行里很容易“翻车”。

1.1 风电出力的不确定性到底难在哪

风电出力本质上是一个强随机过程，取决于风速。研究里最常用Weibull分布去拟合风速的统计特性，再用风机功率曲线把风速映射成出力。这个映射不是线性的，切入风速、额定风速、切出风速三个断点会让场景分布产生明显的截断效应。再加上风电还有反调峰特性：白天负荷高峰风速可能很低，深夜负荷低谷反而风大。你要是只用确定性预测曲线做调度，大概率会得到“低谷时段大量弃风、高峰时段备用不足”的结果，实际并网运行根本不敢这么执行。

还有一个大家容易忽略的问题：不同时间断面之间的风电出力具有很强的时序相关性。比如凌晨西北风系统过境，很可能连续十几个小时风电出力都偏高，而不是每个小时独立波动。如果场景生成时忽略了这种时间相关性，调度结果的备用配置会被严重低估。我早期用独立抽样生成风速场景，每个时段单独抽分布，结果模型算出来的成本很低，一加上时序检验就崩了，后来改成时序相关的抽样才正常。

1.2 负荷侧的波动不只是“加个误差”

负荷预测误差一般用正态分布描述，均值为预测值，标准差取预测值的2%到5%。但负荷侧不确定性不只是“数值上偏一点”这么简单。它同样存在日内的连续变化和突发突变，比如气温骤降、大型工业用户临时启停、节假日负荷特性偏移等。更关键的是，风电出力误差和负荷预测误差之间存在相关性：如果风速预测偏低，往往意味着天气系统判断错误，温度变化也会带动负荷预测偏差。所以做源荷联合场景的时候，不能简单地把两个独立随机变量拼在一起，最好用协方差矩阵描述它们之间的相关性，再通过Cholesky分解生成联合场景。这一步看起来只是数学处理，但对调度结果的保守程度影响很大。

用一个生活类比：出门前看天气预报说“下午有雨”，于是带伞；但又看到商场大促可能很多人排队，这是两件独立的事。如果你把“下雨”和“商场人少”误当成独立事件，可能带伞却人挤人，或者不带伞却被淋湿。源荷两侧不确定性也是同一个逻辑，必须当成一个联合随机过程来考虑，而不是各算各的。

1.3 低碳调度不只是“少烧煤”

低碳调度的本质是在传统经济调度的基础上，把碳排放变成一个可量化的成本项或硬约束。现在国内研究里最常用的是碳交易机制：给每台火电机组分配一个免费的碳排放配额，实际排放量超出配额要花钱买，低于配额可以把差额卖掉。这相当于给“低碳”标了个价格，调度目标从单纯省煤耗变成了“煤耗成本+碳交易成本+弃风惩罚”综合最小化。

如果只把碳排放当成一个固定上限约束，问题会简单很多，但实际效果并不好。因为碳上限设得太紧，系统可能无解或被迫大量上高成本机组；设得太松，又起不到减排作用。碳交易机制的好处是它给了系统一个弹性：在负荷顶峰，允许用高排放但低边际成本的机组顶上，多买碳配额；在负荷低谷，让高排放机组少发，省下来的配额卖出去赚钱。我在模型里把碳交易成本直接放进目标函数，同时保留一个系统总碳排放的软约束场景化处理，这样既体现低碳目标，又不会因为刚性约束导致无解。

2. 低碳调度模型怎么搭：目标、约束与不确定性方法

这个部分我说一下模型的主干。这里不贴完整推导，先把思路捋清楚，因为完整程序里最复杂的其实就是两件事：一是目标函数里怎么把成本、碳交易、弃风切负荷惩罚揉在一起；二是每个时段的功率平衡和备用约束怎么在不确定性场景下成立。

2.1 目标函数怎么设计

假设系统里有常规火电机组、风电场和负荷，调度周期为24小时。目标函数我通常写成三个部分之和。

第一是运行成本，包含煤耗成本和启停成本。煤耗成本一般用二次函数表示，但在混合整数规划里二次项会让求解变慢，实际程序里我会用分段线性近似。第二是碳交易成本。采用基准线法时，碳交易成本等于碳价乘以（实际碳排放量减去免费配额）。这里的实际碳排放量用机组出力乘排放强度累加得到，配额可以按机组容量或历史排放水平设定。第三是弃风惩罚和失负荷惩罚。风电消纳不是无限度的，场景化调度里，某些极端场景可能需要弃风或者切除少量负荷，否则功率平衡可能无解。给这些操作设一个较高的惩罚系数，让优化器只在极端场景下才允许使用。

目标函数的表达式不复杂，但要正确处理“期望成本”。因为场景法里每个场景有不同风电出力、不同负荷，优化目标应该是对所有场景概率加权后的总期望成本，而不是拿某一个确定场景算。这样算出来的调度方案，才是面对未来不确定性时的平均最优，而不是某个“预测值”下的最优。

2.2 约束条件清单

约束条件是整个模型里最容易漏项也最容易导致无解的地方。我按类别列一下：

• 功率平衡约束：每个时段，所有机组出力加上风电消纳量要等于负荷需求。在场景法中，这个约束需要针对每个场景单独成立。
• 机组出力上下限：常规机组的出力不能越过技术最小出力和最大出力，这个不难，但要注意启动状态变量对出力的限制。
• 爬坡约束：机组在相邻时段的出力变化不能超过爬坡速率，风电波动越大，这里越容易成为瓶颈。
• 旋转备用约束：为了保证可靠性，系统需要预留一定上备用和下备用。上备用用来应对风电突然减少或负荷突增，下备用用来应对风电突增或负荷突减。场景法中备用通常取场景偏差的一定比例或固定值。
• 碳排放约束：除了目标函数里的碳交易成本，我会再加一个系统总碳排放量的上限约束。这个上限不要设成硬性死值，可以设置成基准排放的百分比，否则极端场景下容易无解。
• 弃风和切负荷变量约束：弃风量和切负荷量必须非负，且不能超过该时段的可用风电或负荷值。这一步给求解器留了“应急出口”。

每条约束在Matlab里都对应一组矩阵或者YALMIP的表达式，真正写起来并不难。难的是场景数量和约束数量的增长关系：每增加一个场景，功率平衡和备用约束就扩展一份。如果不做场景削减，500个场景直接进MIP，内存和求解时间都会爆炸。

2.3 不确定性建模：场景法与鲁棒优化怎么选

含不确定性的调度建模主流方法大致分三类。第一类是随机规划，核心思路是生成大量场景，用期望值做目标，代表方法是两阶段随机规划。第二类是鲁棒优化，不知道概率分布时用区间或盒式不确定集，追求最坏情况下的安全，但结果通常偏保守。第三类是分布鲁棒优化，介于两者之间，用模糊集描述分布不确定，理论漂亮但实现复杂度高。

如果只是做Matlab代码实现，我个人推荐先用场景法。原因很实在：场景法概念直观，代码容易调试，而且能够充分利用风电出力和负荷预测误差的历史统计信息。你生成500个初始场景，削减到10个左右，求解结果不仅有调度计划，还能看到不同场景下的风电消纳情况和切负荷风险，这对写分析报告非常有帮助。鲁棒优化虽然不用生成场景，但是不确定集合边界怎么定、对偶约束怎么转化，都要花不少功夫，而且最终解的保守性很难向导师或甲方解释清楚。

场景法里还有一个细节：两阶段决策结构。第一阶段决策是机组启停和基本出力，这些必须在看到实际风电和负荷之前确定；第二阶段是弃风、切负荷和机组调整量，这些可以根据场景实现后确定。在代码中，第一阶段变量用binvar和sdpvar声明，第二阶段变量只在对应场景约束里使用。目标函数里的期望成本，就是对所有场景的第二阶段惩罚项做概率加权。这部分写代码时很容易犯错，我后面会详细说。

3. Matlab实现：从场景生成到求解出图

这一部分是实操量最大的地方，我把自己的程序框架和核心代码思路都放出来。由于完整代码有上千行，这里只保留最关键的逻辑片段。你完全可以按照这个框架自己搭一套，数据换成自己系统的就行。

3.1 程序架构与数据准备

我习惯把程序分成四个模块：参数设置、场景生成、模型求解、结果输出。参数设置模块包括机组参数、负荷预测曲线、风速/负荷不确定性参数、碳交易价格和配额等。机组参数建议放在表格里，读入后方便修改。为了测试方便，我用了6台火电机组加1个风电场的简化系统，24小时调度周期，初始场景数500，削减后保留10个场景。

机组参数表大致长这个样子：

风电场额定容量300MW，切入风速3m/s，额定风速12m/s，切出风速25m/s。负荷预测数据我直接取了一个典型日曲线，这里不逐小时列出。碳价设为40元/吨，免费配额按火电装机容量的一定比例分配，实际排放超过配额的部分按碳价计入成本。

3.2 用Matlab生成源荷联合场景的三个关键步骤

第一步是生成初始场景。风速用Weibull分布抽样，负荷误差用正态分布抽样。需要注意的是风速是时序相关的，所以我会用马尔可夫链或者简单的一阶自回归模型来生成风速序列。自回归系数根据历史风速数据拟合，通常设成0.85左右，表示当前时段风速对下一时段影响较强。代码如下：

% 关键变量 N = 500; % 初始场景数 T = 24; % 调度时段 v_shape = 2.0; % Weibull形状参数 v_scale = 6.0; % Weibull尺度参数 rho = 0.85; % 风速自回归系数 wind_scn = zeros(N, T); load_base = load_forecast; % 24x1负荷预测列向量 load_err_std = 0.03; % 负荷预测标准差比例 load_scn = zeros(N, T); % 生成时序相关风速场景 for i = 1:N v = zeros(1, T); for t = 1:T innovation = wblrnd(v_shape, v_scale); % Weibull噪声 if t == 1 v(t) = innovation; else v(t) = rho * v(t-1) + sqrt(1 - rho^2) * innovation; end end wind_scn(i, :) = v; load_scn(i, :) = load_base' .* (1 + normrnd(0, load_err_std, 1, T)); end

第二步是把风速转成风电出力。这一步对应风机功率曲线，可以用分段函数处理。切入风速以下或切出风速以上出力为0，切入风速到额定风速之间近似线性上升，额定风速到切出风速之间保持额定出力：

% 风机功率曲线 v_ci = 3; v_r = 12; v_co = 25; P_r = 300; wind_power = zeros(N, T); idx_linear = (wind_scn >= v_ci) & (wind_scn <= v_r); wind_power(idx_linear) = P_r .* (wind_scn(idx_linear) - v_ci) / (v_r - v_ci); idx_rated = (wind_scn > v_r) & (wind_scn < v_co); wind_power(idx_rated) = P_r;

第三步是考虑源荷相关性。严格的实现是两个随机变量通过协方差矩阵联合抽样，但工程上可以先独立抽样，再用Cholesky分解对已生成的风电出力和负荷误差做线性变换。要注意变换之后风电功率范围可能越界，所以最后需要再做一次截断或重新归一化。这一步我一般放在场景削减之前，避免把相关性破坏掉。

3.3 场景削减：同步回代消除法

初始场景往往有几百个，如果全塞进MIP，变量数量和约束数量会剧增。我实测过，500个场景、24时段、6台机组的目标函数展开之后，YALMIP可能要卡上几分钟甚至报内存错误。场景削减的标准做法很多，我比较推荐同步回代消除法，它实现简单且效果稳定。

核心思路是每次删掉一个场景，并把被删场景的概率累加到离它最近的场景上，使剩余场景集合与原始场景分布之间的概率距离增量最小。判断“最近”的距离可以用欧氏距离，也可以把风速场景和负荷场景合并成一个高维向量再计算。代码如下：

scn = [wind_power, load_scn]; % 合并源荷维度 probs = ones(N, 1) / N; K_keep = 10; % 削减后保留场景数 D = pdist2(scn, scn, 'euclidean'); D(1:N+1:end) = Inf; while size(scn, 1) > K_keep min_val = inf; for i = 1:size(scn, 1) [v_j, idx_j] = min(D(i, :)); inc = probs(i) * v_j; if inc < min_val min_val = inc; del_i = i; del_j = idx_j; end end probs(del_j) = probs(del_j) + probs(del_i); probs(del_i) = []; scn(del_i, :) = []; D = pdist2(scn, scn, 'euclidean'); D(1:size(scn,1)+1:end) = Inf; end

这个削减算法在实际运行时有个小坑：pdist2每次重算距离矩阵，如果场景数很大，循环次数很多，会非常慢。我处理的办法是先用K-means之类的粗聚类把场景缩到50个以下，再用同步回代精确削减到10个。这样精度损失很小，但速度能快一两个数量级。

3.4 在YALMIP中建立优化模型并求解

求解部分我用YALMIP做建模层，求解器用Gurobi或者Cplex都行。如果没装商业求解器，也能用SCS或OSQP跑部分模型，但混合整数问题最好还是用商业求解器，单纯形分支定界的性能差别很大。

模型变量分成三组：

• 机组启停变量：u，二值变量，24×6。
• 机组出力变量：P，连续变量，24×6。
• 弃风变量和切负荷变量：每个场景下24×1的连续变量。

注意第二阶段变量是按场景设置的，比如弃风变量就是24×6×场景数的大矩阵。YALMIP里我喜欢用细胞数组存，这样约束循环写起来更清晰。

核心建模代码如下：

% 变量定义 u = binvar(nG, T); % 启停状态 P = sdpvar(nG, T); % 机组出力 for s = 1:K_keep curt{s} = sdpvar(1, T); % 弃风 shed{s} = sdpvar(1, T); % 切负荷 end % 目标函数 fuel_cost = sum(sum(coe_a .* repmat(P.^2, [1 1]) + coe_b .* P)); % 实际代码中二次项替换为分段线性 carbon_cost = carbon_price * ... sum(sum(emission_factor .* P)) - free_quota; obj = fuel_cost + carbon_cost + ... sum(probs .* (sum(curt{s}, 2)' * penal_curt + sum(shed{s}, 2)' * penal_shed));

约束添加的时候，功率平衡约束写成每个场景下成立：

Constraints = []; for s = 1:K_keep Constraints = [Constraints, ... sum(P, 1) + wind_power_sce(s, :) - curt{s} == load_sce(s, :) - shed{s}]; Constraints = [Constraints, ... 0 <= curt{s} <= wind_power_sce(s, :)]; Constraints = [Constraints, ... 0 <= shed{s} <= load_sce(s, :)]; end

比较关键的一点是：备用约束不要写成“所有场景都要满足”的确定性约束，否则场景法就变成鲁棒优化了，结果会非常保守。更合理的做法是让备用需求与负荷水平绑定，并考虑风电出力预测误差的标准差：

reserve_up = 0.05 * load_forecast + 0.08 * (pred_wind_power - min_wind_power); Constraints = [Constraints, ... sum(P_max .* u, 1) >= load_forecast + reserve_up];

然后调用求解器：

ops = sdpsettings('solver', 'gurobi', 'verbose', 2, ... 'gurobi.MIPGap', 0.01, 'gurobi.TimeLimit', 300); optimize(Constraints, obj, ops);

MIPGap设成0.01就够了，没必要追求严格最优。很多场景下，相对误差从1%降到0.1%要花几倍时间，但对工程结果没有本质区别。我第一次跑的时候把MIPGap设成0.0001，结果一个小时都没跑完，后来改成0.01，两分钟内出解。

3.5 结果分析与可视化

求解之后我习惯画三张图。第一张是机组出力堆叠图，横轴是24小时，纵轴是各机组出力，风电部分也堆叠上去。第二张是不同场景下的风电消纳对比图，把削减后10个场景的消纳区间画出来，能直观看到预测误差对风电消纳的影响。第三张是碳排放成本随碳价变化的敏感性曲线，改变碳价从20到100元/吨，看系统总成本和碳排放量的变化趋势。

画图用Matlab自带bar和plot就行。要注意把YALMIP里sdpvar求解后的结果用value()取出来，直接画会出错。另外堆叠图需要把风电放在最上面或最下面，不然负值会导致图很难看。图里的坐标轴、图例都写好，论文直接用这个图基本不用再加工。

4. 跑代码时容易踩的坑和排查手册

做这类调度模型，最大的敌人不是数学推导，而是模型写好了之后求解器给你一个“无解”或者“结果明显不合理”的输出。我把自己调试过程中遇到的典型问题整理成一张速查表，希望能省掉你半天排查时间。

4.1 场景数选多少合适

场景数是个绕不开的权衡。初始场景太少，分布覆盖不足；削减后保留场景太少，极端风电情况可能被削掉。我的经验值是初始场景300到500个，削减后保留10到15个。如果你用了K-means粗聚类，可以保留20个再精削减。这样既能保留典型极端情况，又能让MIP在2分钟内收敛。如果场景里包含强相关性，不要用随机抽样，尽量用LHS或Sobol序列生成初始场景，分布覆盖率会好很多。

4.2 风电功率序列越界和负值问题

风速转功率时，最容易出现的就是功率曲线函数写得不够严谨，导致风速低于切入风速时功率非零，或者风速超过额定风速时功率超过额定容量。我建议转功率之后立即加一行截断：

wind_power = max(0, min(P_r, wind_power));

另外，如果做了场景相关性变换，变换后的功率值可能不在正常范围内，也必须在变换后重新截断。但要注意，截断会破坏一部分相关性，所以协方差矩阵设计时要预留一点裕度，不要让两个变量强相关到0.95以上，否则截断之后相关性和预期差很远。

4.3 碳交易参数的敏感性

碳价和配额会影响优化结果的方向。碳价太低，系统发现买碳配额比调整出力结构更便宜，就不会主动减排，风电消纳率可能不高。碳价太高，系统会过度弃火电，导致煤耗成本上升，但碳排放确实下降。我建议在一开始先用一个中等碳价跑通模型，然后做敏感性分析，观察碳排放总量和总成本的帕累托前沿。如果你的场景里碳配额是按历史排放设置的，记得每个场景的实际排放不同，而配额通常是一个固定值，不要不小心把配额写成场景相关变量。

4.4 备用系数怎么取值

备用需求可以用负荷比例加风电出力波动区间来表示。我常用的是上备用=5%负荷+8%风电装机容量×当前预测出力比例，下备用=5%负荷。这个系数没有统一标准，要根据你系统的实际情况调。原则上，备用系数越大，模型越稳健，但调度成本越高；系数太低，可能在极端场景下出现切负荷。可以先跑几组敏感性实验，找到成本和风险的拐点。

4.5 二次煤耗成本的处理

很多教科书上煤耗成本用二次函数，但直接放进MIP会带来凸二次项。Gurobi/Cplex能直接处理二次目标，但和线性约束混在一起时求解速度会慢。实际做工程时，我通常把机组出力区间分成三段，用分段线性函数逼近二次成本曲线。分段点取在技术出力范围的1/3和2/3附近，每段的斜率先算好，然后用一组连续变量表示各段出力，加上总和约束。这样目标函数变成线性，求解速度明显加快，精度损失不到1%。如果你的场景数很多，这一步非常值得做。

说了这么多，最后聊一点个人体会。我一开始也迷信复杂的不确定性建模方法，觉得鲁棒优化、分布鲁棒才是“高水平”。但真到了写代码、调参、出结果的时候，发现把场景法做得扎实，把源荷相关性处理好，已经把95%的工程问题解决了。这个题目后续其实还有很多可以扩展的方向，比如加入需求响应、碳捕集电厂、储能，或者把调度尺度从日前延伸到日内滚动，都是在这个框架上做加法。如果你也在做类似的东西，建议先跑通这个基本版，再一步一步加复杂度。模型不是越复杂越好，能稳定出结果、能解释清楚，才是眼下最要紧的。