☰
含风电的低碳经济调度:源荷不确定性建模与Matlab实现
2026/9/28 15:45:07 网站建设 项目流程

1. 源荷双侧不确定性建模:这是调度决策的地基

搞含风电的低碳调度,第一个绕不开的问题就是不确定性怎么描述。很多人一上来就急着写目标函数、列约束,结果把风电出力和负荷预测误差随随便便用个期望值替代,最后仿真结果漂亮得不行,但拿到实际系统里根本不敢用。这就是典型的“模型失真”。

我个人的经验是:不确定性建模的精细程度,直接决定了调度方案能不能落地。所谓“源荷两侧不确定性”,拆开看就是电源侧的随机性和负荷侧的随机性。电源侧主要是风电出力的波动性和间歇性,跟风速的随机特性强相关;负荷侧则是用户用电行为的不可精确预知性,哪怕有短期负荷预测,误差也永远存在。两侧不确定性叠加起来,如果调度模型不考虑,那备用容量配置、碳排放量计算、机组启停安排全部都会跑偏。

处理源荷不确定性,主流做法分两派:随机规划(Stochastic Programming)和鲁棒优化(Robust Optimization)。随机规划的思路是给不确定性变量赋予概率分布,通过蒙特卡洛抽样生成大量场景,然后用场景期望来逼近真实情况。鲁棒优化则是构造一个不确定集,要求调度方案在不确定集内的所有可能实现下都可行,牺牲经济性换鲁棒性。

1.1 风电出力场景:从风速分布到功率曲线的完整链路

风电出力的不确定性,根子上来自风速的随机性。做场景模拟的时候,最常见的手段是用**威布尔分布(Weibull Distribution)**去拟合历史风速数据。形状参数和尺度参数可以用极大似然估计从实测风速序列里标定出来。

风速模拟出来之后,还要经过风电功率转换关系才能得到出力场景。工程上通常用分段函数来描述:

[ P_{WT}(v) = \begin{cases} 0 & v < v_{in} \text{ 或 } v \geq v_{out} \ \dfrac{v - v_{in}}{v_r - v_{in}} P_r & v_{in} \leq v < v_r \ P_r & v_r \leq v < v_{out} \end{cases} ]

其中 (v_{in}) 是切入风速,(v_{out}) 是切出风速,(v_r) 是额定风速,(P_r) 是风电机组额定功率。这套公式做仿真的人应该都很熟了,但这里有个细节经常被忽略:实际风速到功率的转换并不是完全确定的,还存在尾流效应、风机老化、偏航误差等因素带来的转换不确定性。如果想把模型做细,可以在功率计算后面再加一个正态分布扰动项,模拟转换过程中的随机偏差。

1.2 负荷预测误差:不是简单加个正态分布就完事了

负荷侧的不确定性相对温和一些,毕竟负荷的时序规律比风速强得多。但“温和”不代表可以马虎处理。负荷预测误差一般假设服从均值为零的正态分布,标准差跟预测时域长度成正相关,预测时段越远,误差越大。这个特性在做多时段调度的时候特别重要——你不能用同一个误差方差去描述第1小时和第24小时的负荷预测。

实际做算例的时候,很多论文直接在原始负荷曲线上叠加一个固定比例的正态扰动就完事了。这当然可以,但更合理的做法是让扰动幅度随预测时段递增:

[ P_{L,t}^{actual} = P_{L,t}^{forecast} + \varepsilon_t, \quad \varepsilon_t \sim \mathcal{N}(0, \sigma_t^2) ]

[ \sigma_t = \alpha \cdot P_{L,t}^{forecast} \cdot \sqrt{t/24} ]

如果做的是日前调度,24小时的预测误差方差拉成一个递增序列会更贴近实际。别小看这个细节,它会影响机组备用容量在各个时段的分配结果。

1.3 场景削减:不削减的随机优化就是灾难

MC抽样出来的场景数量直接决定了求解规模。假设你抽了1000组源荷联合场景,每个场景都要满足功率平衡、机组出力、爬坡等一堆约束,模型求解的时间会膨胀到不可接受。

我自己跑过的算例里面,IEEE 30节点系统,5台火电机组加1个风电场,如果把调度时段展开成24小时,不做场景削减直接丢给CPLEX去解,经常要跑几个小时甚至解不出来。所以场景削减不是可选项,是必选项。

主流的削减手段是K-means聚类或者同步回代消除法(Simultaneous Backward Reduction)。以K-means为例,流程就是:

  1. 对MC生成的 (S) 个源荷联合场景进行聚类,把相似场景合并成一个典型场景。
  2. 聚类的特征向量可以是“各时段风电出力+各时段负荷”拼接成的高维向量。
  3. 用轮廓系数(Silhouette Coefficient)或者肘部法则选择最优聚类数 (K)。
  4. 聚类的中心作为代表场景,每个类内场景数量占比作为该场景的权重。

这样处理后,几千个场景压缩成5~10个典型场景,求解效率直接质变。我在代码里用的是K-means,实现简单,效果也稳定。同步回代的效果其实更好,但代码实现复杂一些,适合做进阶研究。

提示:场景削减的时候要特别注意,削减后场景的总体统计特性(均值、方差、相关性)会发生变化。做完削减后最好做个检验,对比一下削减前后的风电出力和负荷的期望值曲线,偏差超过5%就要考虑增加聚类数。

2. 低碳调度模型:碳排放约束是怎么“嵌”进去的

低碳调度的“低碳”二字,不是把碳排放当成事后统计指标,而是要把排放约束引入优化模型,让调度决策朝着低碳方向主动调整。这里涉及碳交易机制、碳捕集设备、低碳机组组合等多个维度的建模选择。

2.1 目标函数:经济成本 + 碳成本的双层结构

传统经济调度的目标函数是发电成本最小化。低碳调度通常在此基础上叠加碳排放成本,变成一个双目标、或者单目标加权的问题。最常用的做法是用碳交易价格把碳排放量折算成成本项,从而合并成一个单目标:

[ \min ; C = C_{fuel} + C_{carbon} + C_{curtail} ]

其中 (C_{fuel}) 是机组燃料成本,用二次函数拟合:

[ C_{fuel} = \sum_{t=1}^T \sum_{i=1}^N \left( a_i P_{i,t}^2 + b_i P_{i,t} + c_i \right) ]

(C_{carbon}) 则是碳排放成本:

[ C_{carbon} = \lambda_{CO2} \sum_{t=1}^T \sum_{i=1}^N \left( \delta_i P_{i,t} - E_{free,i} \right) ]

这里 (\lambda_{CO2}) 是碳交易价格(元/吨),(\delta_i) 是机组 i 的碳排放强度(吨/MWh),(E_{free,i}) 是免费碳排放配额。当实际排放超过配额时需要购买碳排放权,低于配额则可以出售获利。这个机制做出来的效果就是:低碳机组因为碳成本低,在调度中会更有竞争力;高碳机组如果碳价足够高,就可能被压到最小出力甚至停机。

如果考虑阶梯碳价(碳价随排放量分段递增),效果更贴近实际,但模型会引入整数变量,复杂度上一个台阶。初学者建议先做线性碳价,跑通了再升级。

2.2 约束条件:不是列完就完事了

约束条件通常包括系统功率平衡、机组出力上下限、爬坡约束、旋转备用约束、碳排放上限约束等。功率平衡约束是所有调度模型的核心:

[ \sum_{i=1}^N P_{i,t} + P_{W,t} = P_{L,t} ]

这里 (P_{W,t}) 是风电并网功率。注意,在含风电的模型中,风电通常有两种处理方式:一种是“全额消纳”,即 (P_{W,t}) 等于预测值;另一种是“允许弃风”,即可以在功率平衡约束中引入弃风变量,让模型自行决定风电出力,代价是弃风惩罚项进入目标函数。第二种更贴近实际,因为极端场景下全额消纳会迫使火电机组深度调峰甚至停机,经济性和安全性都受不了。

爬坡约束是多时段调度里的关键:

[ -Ramp_i^{down} \leq P_{i,t} - P_{i,t-1} \leq Ramp_i^{up} ]

这个约束最坑的地方在于:它会引入时间耦合,让各时段的决策变量不再是独立的,模型求解规模和难度都会上升。如果你用YALMIP建模,这个约束用循环逐时段添加就能搞定,但循环速度慢;用矩阵形式批量添加会快很多。

旋转备用约束是不确定性建模的核心体现:

[ \sum_{i=1}^N (P_{i}^{max} - P_{i,t}) \geq R_{up} + \xi_{W} \cdot P_{W,t} ]

[ \sum_{i=1}^N (P_{i,t} - P_{i}^{min}) \geq R_{down} + \xi_{L} \cdot P_{L,t} ]

典型做法是按“负荷的5%~10% + 风电出力的15%~20%”来设定备用需求。这样风电出力越大的场景,需要预留的上调备用越多,确保风电突然掉出力时系统不会失负荷。

碳排放约束一般用系统碳排放总量上限来控制:

[ \sum_{t=1}^T \sum_{i=1}^N \delta_i P_{i,t} \leq E_{total}^{max} ]

这个约束的作用是在碳交易成本之外再叠加一个刚性排放上限。两种机制同时启用时,碳价兜底、总量封顶,模型会同时受经济和环保双轮驱动,结果更接近真实碳约束环境下的调度行为。

2.3 典型场景集下的运行约束处理

既然场景削减后得到的是一个典型场景集合,调度模型就需要保证方案在这些场景下都可行。最常用的方式是把功率平衡约束扩展成带松弛变量的形式,或者对每个场景单独列约束。

比如,针对每个典型场景 (s),功率平衡写成:

[ \sum_{i=1}^N P_{i,t} + P_{W,s,t} - P_{cur,s,t} = P_{L,s,t} ]

(P_{cur,s,t}) 是场景 s 时段 t 的弃风功率,在目标函数中用惩罚系数约束它尽量小。

这种做法的好处是:机组出力的决策变量是第一阶段的“here-and-now”决策,不随场景变化;而弃风、切负荷是第二阶段的“wait-and-see”决策,可以随场景调整。这个其实是两阶段随机规划的雏形,代码实现上就是让一部分变量带场景下标,一部分不带。YALMIP里建这个非常方便,就是分别定义两套变量的问题。

2.4 低碳机制:碳捕集与电转气

如果题目里的“低碳”只是碳排放目标加一个总额约束,那其实不够“碳达峰、碳中和”时代的技术感。更进一步的低碳调度会引入**碳捕集与封存(CCS)设备和电转气(P2G)**设备。

碳捕集设备运行时有个典型特性:捕碳需要消耗额外电能,这部分能耗会降低机组的净出力。所以带CCS的机组约束变成:

[ P_{i,t}^{net} = P_{i,t}^{gross} - P_{i,t}^{capture} ]

净出力等于毛出力减去捕集能耗。捕集的CO₂量跟机组排放量成正比,被捕获的碳量还可以卖给P2G或封存。这样一来,机组出力约束、碳排放约束都跟捕集状态耦合在一起。

说实话,CCS和P2G的建模细节非常多,初学者不建议一上来就全模型铺开。先做好常规的碳成本模型和排放约束,把基础调度结果做扎实,再逐步添加上这些低碳单元。

3. Matlab实现:从数学模型到可运行代码的全流程

学电力系统优化调度,最大的坎不是数学本身,而是怎么把一个有约束的最优化模型在Matlab里正确表达出来。很多人看着模型公式觉得都懂了,一写代码就报错。分享下我的完整实现流程,照这个思路走能少踩很多坑。

3.1 求解框架选型:YALMIP+CPLEX是黄金组合

Matlab里求解带约束的非线性优化问题,工具箱选择很关键。我强烈推荐用YALMIP做建模层,配CPLEX或Gurobi做求解器。下面的MATLAB代码示例展示了如何配置YALMIP与CPLEX求解器。

YALMIP的核心理念是你只需要声明“决策变量”和“约束条件”,然后调用optimize函数求解,它自动处理变量类型检测、线性化、求解器切换等脏活累活。相比直接用CPLEX的底层API,YALMIP的代码可读性和调试性都高一个档次,后续要改模型复杂度也方便得多。我在论文复现里基本都是YALMIP+CPLEX,很少翻车。

注意:YALMIP自带默认求解器,如果你装了CPLEX但没配置,YALMIP会在第一次运行时自动检测。确认安装成功后,可以用yalmiptest命令检查求解器是否正常挂载。我第一次用的时候CPLEX权限没配好,求解器一直走的默认分支定界,速度慢到怀疑人生,排查半天才找到原因。

3.2 代码整体架构:主函数+数据文件+建模函数三段式

我不建议把所有代码塞进一个大脚本文件里。模块化是工程世界的水电煤,科研代码也不例外。推荐这样一个结构:

  • main.m:主入口,负责数据加载、模型调用、结果绘图。
  • load_data.m/case_data.m:数据文件,存放机组参数、负荷曲线、风速参数、碳价等。
  • generate_scenarios.m:生成源荷联合场景并削减。
  • build_model.m:用YALMIP表达式构建决策变量、目标函数和约束。
  • solve_model.m:调用求解器,处理求解结果和输出。

各模块用函数封装,参数通过结构体传递。这样换一组成数据测试跟你迁移代码改三个地方的约束又不小心引进奇怪的bug之间,只隔了一层窗户纸。

3.3 核心变量定义与约束添加

决策变量分两类:不随场景变化的机组出力和随场景变化的弃风切负荷量。第一类用矩阵 (N_unit \times T),第二类用三维矩阵 (N_scenario \times T \times N_unit) 的形式定义。MATLAB代码示例如下:

P = sdpvar(n_gen, T, 'full'); % 机组出力,第二阶段Isn't决策 P_w = sdpvar(n_scenario, T, 'full'); % 风电接纳功率(随场景变化) P_cur = sdpvar(n_scenario, T, 'full'); % 弃风功率 slack_load = sdpvar(n_scenario, T, 'full'); % 切负荷松弛变量

约束添加用的是Constraints = [Constraints; ...]累积语法。下面是一个功率平衡和备用约束的片段:

Constraints = []; for t = 1:T % 功率平衡(对所有场景) for s = 1:n_scenario Constraints = [Constraints, ... sum(P(:, t)) + P_w(s, t) + slack_load(s, t) == ... P_load(s, t) + P_cur(s, t)]; end % 旋转备用约束 Constraints = [Constraints, ... sum(Pmax - P(:, t)) >= reserve_up(t) + 0.15 * P_w_forecast(t)]; end

这里有个我在实际调试中养成的习惯:功率平衡约束不要用“等于”,用“大于等于”或“小于等于”加松弛变量。因为只要有微小数值误差,严格等式约束会导致求解器判定不可行,而引入一个非负的松弛变量后,模型会自动在可接受范围内找到可行解,代价是目标函数里加个极大的惩罚系数。这个技巧在工程中用得非常普遍。

3.4 目标函数构建与求解器调用

目标函数分三块:燃料成本、碳成本、弃风惩罚。燃料成本的二次项可以通过YALMIP自动处理,CPLEX会自行判断问题类型并选择适当的求解算法,以下代码展示了关键的建模部分:

% 燃料成本(二次型) fuel_cost = 0; for i = 1:n_gen for t = 1:T fuel_cost = fuel_cost + a(i) * P(i,t)^2 + b(i) * P(i,t) + c(i); end end % 碳排放成本(碳排放量按线性强度折算) carbon_cost = 0; for i = 1:n_gen for t = 1:T carbon_cost = carbon_cost + lambda_co2 * (delta(i) * P(i,t)); end end % 弃风惩罚 penalty = 100 * sum(P_cur(:)); % 总目标 Objective = fuel_cost + carbon_cost + penalty;

然后一行调用求解器:

ops = sdpsettings('solver', 'cplex', 'verbose', 1, 'solver', 'cplex', 'showprogress', 1); optimize(Constraints, Objective, ops);

求解完成后用value()提取结果:

P_opt = value(P); P_w_opt = value(P_w); P_cur_opt = value(P_cur);

3.5 数据处理:风电场景生成与削减的代码实现

场景生成我用的是Weibull分布抽样加功率转换,代码如下:

% 风速抽样(Weibull分布) v = wblrnd(A_shape, B_scale, n_scenario * T, 1); v = reshape(v, n_scenario, T); % 风速到功率转换 P_w_raw = zeros(n_scenario, T); for s = 1:n_scenario for t = 1:T if v(s,t) < v_in || v(s,t) >= v_out P_w_raw(s,t) = 0; elseif v(s,t) >= v_in && v(s,t) < v_r P_w_raw(s,t) = P_r * (v(s,t) - v_in) / (v_r - v_in); else P_w_raw(s,t) = P_r; end end end

这样得到的是原始场景集,接下来用K-means聚类削减:

% K-means聚类削减场景 [idx, C] = kmeans(P_w_raw, K, 'Distance', 'sqeuclidean', 'MaxIter', 500); % C: 聚类中心(代表场景),每个类的大小作为权重 scenario_prob = histcounts(idx, K) / n_scenario;

注意这里只对风电场景做了聚类。工程上更正规的做法是把风电和负荷联合场景拼成一个特征矩阵一起聚类,比如 ( [P_w(s,:), P_load(s,:)]) 作为一行的维度。因为源荷有一定的负相关性(风大的时候往往负荷低),分开聚类会丢失这种相关性信息,削减后的联合场景就不够保真。我在实践中是把两部分拼在一起做聚类的,效果确实更好。

3.6 结果呈现:绘图与对比分析

仿真的最后是画图。必备的三张图:调度结果堆叠图、风电计划与实际出力对比图、各场景下的碳排量柱状图。调度结果图建议用area图,把每台机组的出力垂直堆叠,能直观看到火电和风电各自承担多少负荷。

我之前踩过的坑是画图时坐标系没调好,风功率曲线和火电曲线叠在一起根本分不清谁是谁。加个legend并把LineWidth设成1.5以上,配色用matlab的经典色序就够了,别整花活。

4. 场景分析与结果解读:调度决策背后的逻辑

模型跑通之后,最关键的环节是对结果做系统性剖析。很多人仿真出图后就算完事了,这是很可惜的——因为你根本没发挥出模型的分析价值。

4.1 不同碳价下的调度结果对比

如果保持模型其余参数不变,单独把碳交易价格从50元/吨逐步提高到300元/吨,你会看到三个显著变化:

第一,高碳机组的出力占比明显下降。碳价上涨增加了高碳机组的边际发电成本,调度模型自然会用低碳机组或者风电去替代它。第二,弃风率下降。高价碳环境下风电的替代效益更突出,模型更愿意承担风电不确定性带来的备用成本去接纳风电。第三,系统总碳排放量下降,但系统总成本可能先降后升——因为碳价涨到一定水平后,低成本发电组合已经挖掘殆尽,新增的碳成本没有足够的低碳改造空间吸收,只能净增加成本。

4.2 场景数与求解精度的权衡

我在实际测试中分别用2、5、10、20个典型场景跑过同一组数据,结果是:从2个场景增加到5个场景时,目标函数值变化明显;从5个到10个也有变化,但幅度缩小;再往上,结果趋于平稳。

所以实际做研究时,场景数不是越多越好。取5~10个典型场景是一个性价比很高的区间,既保留了源荷不确定性的大部分统计特性,又不会让求解时间失控。具体数字根据系统规模微调。

4.3 确定性模型 vs 随机模型的对比

这个对比是审稿人和答辩老师最爱问的问题之一。做法很简单:跑一个用预测值做确定性调度的模型(也就是把不确定性和场景全部去掉),再用同一组数据跑场景模型。你会看到:确定性模型由于免费占用了“完美预测”的红利,机组出力曲线更激进——火电出力更低、风电替换量更大、碳排放更少。但它的“好结果”经不起推敲,因为一旦实际风速或负荷偏离预测值,它的备用不足就会暴露,实际损失远大于优化节省的成本。

两者之差,本质上就是信息价值(Value of Stochastic Solution),这个指标在研究随机规划类问题时一定要提出来。审稿人看到这个分析,你稿子的学术分量会大幅提升。

5. 常见问题与调试经验速查

下面是实际运行这套模型时容易遇到的问题,比看十遍理论都有用。

5.1 YALMIP求解不可行(Infeasible Problem)

这类问题的核心原因往往集中在约束上。我常用的排查步骤如下:

  1. 看Power Balance约束。这是低频翻车点。尤其是加入了备用约束后,系统的可调功率上限可能真的撑不住,这时候不是模型写错,而是物理上无解。把备用约束调低,或者给功率平衡加松弛变量试试。
  2. 检查变量维度是否一致。YALMIP对维度比较宽容,但矩阵不匹配时照样报错。比如P是 (N \times T),而P_load是 (T \times 1) 的列向量,相加的时候维度对不上。
  3. 检查约束列表是否为空。我自己就干过忘记把某条约束append进Constraint集合,模型出于求解器的指令找不到约束而误报不可行的蠢事。
  4. 用诊断工具。YALMIP自带了Check命令,可以在求解前手动验证约束与决策变量是否一致。

5.2 CPLEX错误:No variable with name ... exists

这类问题多是YALMIP与CPLEX之间的接口解析问题。通常原因是:你在同一个Matlab工作区里多次调用YALMIP,导致CPLEX内部变量名称映射错乱。最简单的解决方式是每次求解前clear all,或者把模型构建代码放在独立的函数内,避免全局变量污染。另外,CPLEX版本与YALMIP版本不匹配也会导致该问题。

5.3 场景削减后风电置信区间失真

如果削减后的风电出力范围明显比原始场景集窄,说明聚类数K选太小了。解决方案:一方面增加K值,另一方面在做聚类特征时加入各时段风电出力与负荷的上下界统计量,保证聚类的“极端场景”也能被捕获到。

5.4 Matlab求解速度太慢

如果模型规模不大但求解很慢,先检查模型类型。如果存在二次目标函数导致模型变为QCQP/二次约束,CPLEX会动用更耗时的内点法。若目标是仅含二次项的单目标优化,可以试着将二次项分段线性化,转换成线性规划后再求解,速度会有数量级的提升。

5.5 考虑冷启动时的机组组合误停机

如果原始调度结果里机组频繁启停(比如某一时段出力压到下限附近的机组,下一时段又换成了另一台机组),很可能是爬坡约束太松,或者没有启停成本约束。加入机组启停成本和最小开关机时间约束后,能有效压低机组启停次数。但这个优化会让模型变成混合整数规划,CPLEX求解速度会大幅牺牲,需要权衡。

5.6 风电场上调备用系数取值过大导致无解

如果你把备用系数(风电出力的15%~20%)调高后无解,大概率是系统热备用容量本身不足。两种方案:一是引入机组附加备用容量决策变量(也就是允许备用不足但附带惩罚成本),二是在目标函数中增加切负荷惩罚项作为软约束,优先保可行。

5.7 碳价高时目标函数中出现负成本

碳价升高后,低碳机组因为配额盈余可卖碳而获得收益,目标函数可能为负。这不是bug,这是经济信号:系统正在从“高碳发电组合”切换到“低碳发电组合”。结果画图时注意纵轴负值的表示方式,避免读者误读。

6. 扩展方向与最后的经验

模型本身已经能出完整结果之后,还可以往几个方向扩展:考虑碳捕集设备与P2G的联合调度、引入需求响应资源、把确定性场景调度升级为两阶段分布鲁棒优化等。我个人的建议是先把基础模型吃透,场景调度和不确定性建模搞扎实,再逐步进阶。

最后分享一个我在调试代码时养成的小习惯:每次修改模型参数后,先跑一个简化的小规模算例验证逻辑,再跑全规模。比如先让 (K=3)、(T=6)、机组3台,跑通跑快,再逐步恢复完整配置。这样能极大减少虫子(bug)的排查时间,把精力放在优化本身。

这套Matlab代码如果能配合一组合理的参数,其实复现性很强。整个项目的完整流程——从不确定性建模、场景削减、低碳目标构建到求解分析——如果你能独立走通一遍,你对电力系统调度优化的理解会上一个台阶。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询