搞电力系统优化调度的同学十有八九都绕不开“含风电电力系统低碳调度”这个方向。如果再叠加“源荷两侧不确定性”还要用Matlab把代码完整落地那确实是个看着简单、实际坑不少的事。这篇文章就把我这边的模型思路、代码框架和踩坑记录完整捋一遍适合正在做毕业设计、课题研究或者刚接触调度优化的工程师参考。你会看到一个从问题定义、数学建模、Matlab编码到结果分析的全过程所有关键决策我都会说清楚为什么这么做。1. 先把问题定义说清楚什么叫“源荷两侧不确定性”下的低碳调度1.1 源侧不确定性风电出力不是拍脑袋定的风电出力的随机性来自风速。风速本身是一个强随机过程加上预测模型、气象数据精度的限制风电功率预测很难做到准确。做过实际风电数据拟合的人都知道预测误差会随着预测时长的增加明显变大而且不是简单的正态分布特别是在大风过程切换和早晚时段误差波动往往更剧烈。如果在日前调度里我们把风电出力固定成一个预测值那么实际运行中一旦风电比预测低火电机组可能来不及爬坡系统就得面临切负荷风险反过来如果风电比预测高又可能因为消纳能力不足而弃风。所以源侧不确定性的本质是让风电出力在优化模型里不是一个点而是一组“可能发生的出力过程”每个过程对应一个概率。我们把这组过程叫场景集后面所有优化决策都要在这些场景下都可行才算真正考虑了不确定性。1.2 荷侧不确定性负荷预测误差同样不能忽略很多人一做不确定性调度眼睛就只盯着风电把负荷当成完全已知的参数处理。这在实际工程里是有隐患的。负荷预测同样存在误差尤其节假日、极端天气、大型赛事活动期间预测偏差可能达到5%以上。如果只考虑风电不确定性而忽略负荷侧那么调度方案在源荷双重波动下依然可能不够可靠。更合理的做法是把源侧和荷侧一起建模。最简单的方式是对风电预测误差和负荷预测误差分别采样生成初始场景然后把风电场景和负荷场景组合成一个“源荷联合场景”。联合场景的概率等于两者概率的乘积最后再用场景缩减算法砍到目标数量。这样做的好处是功率平衡约束、备用约束被放在同一个随机环境下考虑不会出现“风电场景对应高负荷、负荷场景却取固定值”这种错位。实际跑下来你会发现加入负荷不确定性之后备用的需求量会明显上升这比单纯增加风电场景更贴近真实运行。1.3 低碳调度和传统经济调度的本质区别传统经济调度的目标是最小化发电成本主要就是煤耗成本、启停成本和运维成本。低碳调度则在这个基础上引入“碳排放”相关成本或约束把环境成本变成调度决策的一部分。常见做法有三种碳税按实际排放量直接罚款简单粗暴但不够贴近市场机制碳排放上限约束直接限制系统总排放量属于硬约束必须满足碳排放配额与碳交易给系统一个免费配额排放低于配额可以卖出碳权获利超过配额则必须购买。目前在学术复现和工程方案里用得最多的是碳交易机制尤其是阶梯碳交易。所谓阶梯就是购买的碳权越多单价越高。这比单一碳价更真实因为碳市场中的大量购入必然抬高成交价。低碳调度和传统经济调度还有一个关键区别传统调度只按煤耗成本排序机组低碳调度还要看边际碳排放成本。如果高碳机组的煤耗低但碳排放强度高在碳价足够大的时候调度结果可能反而降低它的出力换成煤耗稍高但更清洁的机组。这个“排序变化”是低碳调度最值得分析的现象。2. 模型设计从目标函数到约束条件的完整建模过程2.1 目标函数运行成本、碳排放成本、弃风惩罚怎么平衡我常用的目标函数写成期望成本形式min E[ 燃料成本 碳交易成本 弃风惩罚 ] 启停成本其中E[...]表示对所有场景取概率加权平均。燃料成本使用二次函数F_i(P_i) a_i * P_i^2 b_i * P_i c_i这个二次项会让模型变成MIQP或MIQCQP。如果你担心非线性模型求解慢可以把煤耗曲线分段线性化转成混合整数线性规划。我在复现课题时一般直接用二次函数Yalmip能很好地传给Cplex或Gurobi最新版本求解器处理这类问题并不吃力。碳交易成本我建议写成CarbonCost carbon_price * (总排放量 - 免费配额)这个表达式是线性的不包含max函数。也就是说允许系统“卖碳”体现的是净碳交易成本。如果你不希望出现卖碳收入那就要引入辅助变量把超过配额的部分单独表示写法稍微复杂一点。我个人建议先保留卖碳机制这样目标函数更平滑后期如果想改阶梯碳价再叠加分段函数也不迟。弃风惩罚必须加。如果不加模型为了降低碳成本可能会让风电出力尽可能大但现实中风电超过负荷需求时只能弃掉。惩罚系数一般取100到300元/MWh具体取决于你对弃风的容忍度。这个值不能设太小否则结果里会莫名其妙出现大量弃风也不能设太大否则会掩盖其他成本信号。2.2 常规约束功率平衡、机组出力、爬坡、备用约束条件按场景逐一列写。核心约束包括功率平衡每个场景每个时段下所有常规机组出力、风电出力和负荷之间必须相等。如果需要考虑切负荷可以增加负荷失负荷变量但目标函数要加大惩罚。机组出力上下限P_i_min P_i P_i_max同时要乘以启停状态。爬坡约束机组在两个相邻时段间的出力变化不能超过爬坡速率。风电出力上限实际风电出力不能超过该场景下的风电预测值。注意这里允许“低于预测值”也就是允许弃风。旋转备用约束正备用容量要能应对负荷突增和风电突减负备用容量要能应对风电突增。对于备用约束我采用工程上常见的简化正备用不小于系统负荷的5%加上风电预测出力的15%。这个经验公式不一定最优但很容易嵌进模型而且能明显提升结果可靠性。如果你想做得更严格可以把备用容量也作为决策变量或者引入机会约束表示“备用不足的概率不超过某个阈值”。但这会让模型复杂度上一个台阶新手不建议一上来就碰。还有一点很多人容易搞错机组启停状态是第一阶段决策不随场景变化出力是第二阶段决策可以随场景变化。对应到日前调度场景机组开停机计划在知道实际风电和负荷之前就要确定而出力可以在运行中调整。建模时u不能定义成(nGen, T, Ns)而应该定义成(nGen, T)否则每个场景都有一套独立开停机计划既不现实也会让整数变量爆炸求解极慢。2.3 低碳约束碳排放配额与碳交易机制的引入方式免费碳排放配额的计算方式会直接影响调度结果。最稳妥的方式是把配额设为一个常数比如根据系统总容量和运行时段折算FreeQuota quota_rate * sum(P_max) * Tquota_rate的单位是t/MWh可以取0.3到0.5之间的值。这样配额是一个固定数碳交易成本只与机组实际排放量线性相关模型好解。另一种常见方式是按机组实际发电量给配额也就是配额等于某个比例乘以出力。看起来公平但配额率会进入目标函数的系数实际上等于改变了每个机组的边际碳成本模型依然是线性的。这种做法的缺点是“出力越多配额越多”相当于给高排放机组打了折扣低碳激励效果会减弱。阶梯碳价是很多论文里会用到的进阶方式。把碳排放购买量分成几段每段对应不同单价。在Yalmip里实现分段线性函数最简单的方式是用implies和二进制变量或者直接调用求解器内置的分段约束。如果你只是想要一个大概结论建议先用单一碳价跑通再改成阶梯。阶梯碳价的主要作用是观察碳价对调度结果的非线性影响比如超过某个排放量后系统会明显调整机组组合而这个现象在单一碳价下可能不明显。2.4 不确定性处理场景法 vs 鲁棒优化的选择处理源荷两侧不确定性主流路线就两条随机场景优化和鲁棒优化。场景优化的核心思想是把不确定参数离散成若干场景求期望成本最小的调度策略。它的优点是信息利用充分结果不保守代码也直观缺点是需要合理生成场景且场景数量影响计算精度和速度。鲁棒优化则不需要概率分布只需要一个不确定集合比如“风电出力在预测值的±20%范围内波动”“负荷在预测值的±10%内波动”。求出的解必须对集合内所有情况都可行因此保守性更强但能提供最恶劣情形下的保障。鲁棒模型的求解通常要转化成对偶问题或使用列约束生成算法实现门槛更高。如果问我的个人偏好做科研复现和课设先用场景法做实际工程安全校核再用鲁棒优化或分布鲁棒优化。场景法能更直观地展示“不确定性对成本的影响”而且非常适合和Matlab代码结合。你完全可以在同一个框架里先跑场景法再替换成鲁棒模型做对比这样论文里的分析维度会更丰富。3. Matlab代码实现核心模块拆分与关键代码片段3.1 整体文件结构与数据准备写Matlab代码不是简单堆脚本。我自己的项目结构一般是这样main.m % 主入口控制流程 case_data.m % 火电、风电、负荷、碳交易参数 gen_scenario.m % 源荷场景生成 reduce_scenario.m % 场景缩减 build_model.m % 构建Yalmip模型 solve_model.m % 求解并返回结果 plot_results.m % 结果可视化用结构体存数据比如para.gen、para.load、para.carbon这样后面写约束时不用一堆全局变量。以6机系统为例火电参数通常包括各机组最大/最小出力爬坡速率启停成本煤耗二次系数a、b、c碳排放强度系数风电预测曲线和负荷预测曲线。数据可以直接写在case_data.m里也可以用readmatrix读取Excel。我建议把原始数据放Excel参数放.m文件因为调度模型要改参数测试不同工况Excel改起来更方便。3.2 5分钟生成并缩减源荷场景蒙特卡洛K-means场景生成我一般分两步。第一步用蒙特卡洛采样出大量初始场景第二步用K-means聚类缩减到目标数量。虽然严格意义上的最优场景缩减是基于Wasserstein距离的同步回代法但K-means实现简单、结果稳定作为工程近似完全够用。下面是一个简化版的函数输入风电预测和负荷预测曲线输入目标场景数输出缩减后的联合场景和概率function [scen, prob] gen_scenario(wind_pred, load_pred, Ns) % wind_pred: T x 1 风电预测 % load_pred: T x 1 负荷预测 % Ns: 目标场景数 T length(load_pred); N0 500; % 初始采样数 % 预测误差标准差单位MW wind_std 0.1 * wind_pred 5; load_std 0.03 * load_pred 10; scen_all zeros(2*T, N0); for k 1:N0 wind_cur max(0, wind_pred wind_std .* randn(T, 1)); load_cur max(0, load_pred load_std .* randn(T, 1)); scen_all(:, k) [wind_cur; load_cur]; end % K-means缩减 [idx, C] kmeans(scen_all, Ns, Distance, sqeuclidean, Replicates, 5); prob accumarray(idx, 1) / N0; % 每个簇的概率 scen C; % Ns个代表场景 end这里注意两点。第一风电场景必须做非负截断否则会出现负出力这种物理上不存在的场景。第二初始场景数量N0不能太小我一般取300到1000太小代表分布不够充分缩减后场景会比较集中丢失尾部风险。目标场景数Ns的选取我放在后面的常见问题部分展开。3.3 优化模型求解YalmipCplex/Gurobi的配置模型构建用Yalmip求解器用Cplex或Gurobi这一套在Matlab里非常顺。Yalmip的好处是建模语言接近数学表达约束怎么写代码就怎么写不用自己拼矩阵。以下是核心建模片段%% 决策变量 % 第一阶段启停状态不随场景变化 u binvar(nGen, T, full); % 第二阶段常规机组出力、风电出力、弃风量随场景变化 P sdpvar(nGen, T, Ns, full); Pw sdpvar(nWind, T, Ns, full); Wsp sdpvar(nWind, T, Ns, full); % 弃风量 Constraints []; objective 0; for s 1:Ns for t 1:T % 功率平衡 Constraints [Constraints, sum(P(:,t,s)) sum(Pw(:,t,s)) ... Load_scen(t, s)]; % 机组出力与启停耦合 Constraints [Constraints, P(:,t,s) u(:,t) .* Pmax]; Constraints [Constraints, P(:,t,s) u(:,t) .* Pmin]; % 风电出力和弃风量 Constraints [Constraints, Pw(:,t,s) Wind_scen(t, s)]; Constraints [Constraints, Pw(:,t,s) Wsp(:,t,s) Wind_scen(t, s)]; % 爬坡约束注意跨时段 if t 1 Constraints [Constraints, P(:,t,s) - P(:,t-1,s) Ramp_up]; Constraints [Constraints, P(:,t,s) - P(:,t-1,s) -Ramp_down]; end % 目标函数期望成本碳交易成本包含在这里 fuel_cost sum(a .* P(:,t,s).^2 b .* P(:,t,s) c); carbon_cost carbon_price * (sum(co2_intensity .* P(:,t,s)) ... - FreeQuotaPerHour(t)); wind_cost wind_penalty * sum(Wsp(:,t,s)); objective objective prob(s) * (fuel_cost carbon_cost wind_cost); end end % 启停成本属于第一阶段只计算一次 start_cost sum(sum(StartCost .* max(u(:,2:end) - u(:,1:end-1), 0))); objective objective start_cost; %% 求解 ops sdpsettings(solver, gurobi, showprogress, 1, ... gurobi.MIPgap, 0.01); sol optimize(Constraints, objective, ops);这段代码里我把启停变量独立于场景这是两阶段随机规划的关键。如果你把u也定义成含场景的变量问题规模立刻膨胀好几倍而且结果解释起来很别扭。爬坡约束这里只做了同一场景内的时间耦合如果你的场景是“路径式场景”也就是每个场景包含完整24小时时序这样没问题如果你的场景是“快照式场景”那爬坡跨时段约束需要换一种写法这点要特别注意。3.4 结果分析与可视化碳价灵敏度曲线怎么出求解之后第一件事是检查sol.problem是否为0。如果非0说明求解失败或不收敛这时候不要急着看结果先回到约束排查。成功求解后我会用value()取出各个决策变量然后画三张图调度结果堆叠图横轴为24个时段纵轴为火电、风电出力弃风情况图风电预测、风电实际消纳、弃风量碳排放与碳成本图各机组碳排放量、系统总排放、碳交易成本。碳价灵敏度分析是最能体现低碳调度价值的部分。我在做论文时通常让碳价从0变化到200元/t间隔20元跑11次优化记录每次的总排放量和总成本然后画曲线。代码思路很简单price_list 0:20:200; emission_list zeros(size(price_list)); cost_list zeros(size(price_list)); for k 1:length(price_list) yalmip(clear); % 清除上一轮模型防止变量堆积 carbon_price price_list(k); build_model; % 重新构建目标函数和约束 optimize(Constraints, objective, ops); emission_list(k) value(sum(co2_intensity .* P, all)); cost_list(k) value(objective); end plot(price_list, emission_list);注意循环里最好调用一个独立的模型构建函数而不是在脚本里反复修改全局变量。yalmip(clear)这行很关键不然内存里会累积大量sdpvar对象越跑越慢。另外把build_model写成一个函数每次输入碳价输出Constraints、objective可以避免变量覆盖踩坑。4. 常见问题与排查技巧我在调试中踩过的坑4.1 求解器不收敛或内存溢出怎么办这个坑几乎百分之百会遇到。第一个原因就是场景太多或二进制变量太多。如果场景数上了50机组数上了10时段24个二进制变量就是10×24240个看起来不多但二次约束加随机场景会让MIP搜索树迅速膨胀。解决办法是把启停变量独立于场景先降到20个场景再把MIP gap从默认的1e-4放宽到0.01。个人经验是MIP gap取0.01在调度问题上完全够用目标函数误差在1%以内但求解速度能快好几倍。如果还是慢检查是不是用了binvar加上二次目标二次项会让MIP求解更难。这时可以把煤耗成本做分段线性化转成纯MILPGurobi处理纯MILP比MIQP快不少。还有一个容易忽略的问题Yalmip默认会对所有约束做预处理如果约束里出现了像u(:,t) .* P(:,t,s) xxx这种双线性约束问题性质会变成非凸求解器直接崩溃或死循环。所以我上文的示例里专门用了分解写法先让P u .* Pmax再用P u .* Pmin而不是写u * P这类乘积项。内存溢出一般出现在场景生成阶段。如果你一次性生成5000个初始场景每个场景长度48维矩阵也就5000×48不大但如果你把5000个场景全部塞进优化模型那模型会巨大。优化模型里的场景必须在缩减之后使用初始场景只用于缩减训练不能直接进约束。这个顺序别搞反。4.2 场景数量怎么选才合理场景数不是越大越好。我用同一套模型分别测试过5、10、20、30、50个场景结果很有意思5个场景时目标函数波动很大10个场景勉强稳定20个场景以后目标函数和碳排放几乎不再变化50个场景只是增加计算时间。所以我自己一般先取20个场景跑基准结果再取50个场景做一次验证。如果你的风电和负荷波动特别大比如日风电预测误差超过30%建议增大到30到40个场景。场景缩减还有一个容易踩的坑K-means聚类后每个簇的中心可能不代表一个真实可能的时序而是一个“平均场景”。这种平均场景会平滑尖锐波动导致风电的极端情况被抹掉。如果你要重点考虑高风电、低负荷这类极端场景建议改用同步回代消除法或者保留K-means结果中那些概率虽小但出力特别大或特别小的场景。说白了场景缩减不只是减少数量还要保留分布的形状。4.3 碳交易价格灵敏度分析怎么快速做如果只是手动改碳价然后点运行你会累死而且容易出错。写循环是标准做法但有三个注意点。第一碳价不是定义成一个固定常量然后让你在循环里给它赋值就完事因为目标函数里的碳价参与构建sdpvar表达式一旦构建了后面改变量本身不会自动更新。正确做法是把碳价作为模型构建函数的输入参数每次都重新构造目标函数。第二循环次数多的时候必须用yalmip(clear)清理上一轮模型否则变量越积越多Matlab越来越卡。第三如果你用了Gurobi不同碳价下同一个问题结构类似求解器可以用warm start加速但Yalmip默认不会自动传上一个解需要手动assign和sdpsettings(usex0,1)。对代价不敏感的课设场景没必要折腾warm start直接跑就行。4.4 常见问题速查表现象可能原因解决办法求解器报“Gurobi not found”或“Cplex not found”求解器未安装、路径未添加或license无效运行yalmiptest检查求解器状态确认gurobi或cplex可被Yalmip识别结果中风电出力始终等于预测值没有弃风未加风电出力上限或弃风惩罚太大/太小检查约束Pw Wind_scen调整弃风惩罚系数碳交易成本恒为0配额设置过大系统排放始终低于配额调低quota_rate或先画出系统总排放曲线再回填配额场景生成出现负风电用正态分布直接采样改用截断正态或用max(0, ...)处理最好换Beta分布启用阶梯碳价后模型求解时间暴增分段线性约束引入大量二进制变量先用单一碳价跑通确认逻辑后再加阶梯或减少分段数多次循环后Matlab内存不断增长未清理Yalmip变量在循环开头执行yalmip(clear)把模型构建封装成函数这些坑我都踩过最后留几句实在话我在把这个项目从零跑通的过程中最大的体会是不要把模型和代码混在一起调试。我第一次上手时直接就把场景、碳交易、两阶段全塞进一个脚本结果一报错根本分不清是模型逻辑错了还是代码写错了。后来改成先写一个确定性版本确认功率平衡、爬坡、机组组合都合理再逐步加入场景、加入碳交易整个过程顺畅很多。另一个体会是源荷两侧不确定性的核心不是场景堆得越多越好而是每个场景都能讲出物理意义风电场景要符合风速波动规律负荷场景不能出现离谱的尖峰否则优化结果只是数学上可行工程上没法用。如果你也想快速复现这个课题建议按我上文的顺序走先定模型再写数据再写场景生成最后写优化求解。这样哪怕中间卡住了你也知道该回头改哪里。