MATLAB代码电力系统低碳调度源荷不确定性与风电接入——从建模思路到yalmip落地去年我做了一个含风电接入的电力系统低碳调度项目代码用MATLAByalmip实现。接手之前我以为最大的麻烦会是碳约束和风电随机性怎么建模真正动手之后才发现模型原理其实都能查得到真正让人头疼的是yalmip的变量设计、约束组装和求解器之间的配合。标题里这个组合——低碳调度、源荷不确定、风电——几乎涵盖了当前调度类代码最常被问到的三个点这篇文章我把我实际写代码的过程、踩过的坑和最后跑通的思路完整梳理一遍给正准备做同类仿真的朋友一个参考。整个模型解决的核心问题可以概括成一句话在风电出力和负荷都有预测误差的前提下安排常规机组出力、风电消纳和碳交易策略让系统总运行成本最低同时满足功率平衡、机组爬坡、备用容量和碳排放等约束。代码本身不算复杂但里面有几个边界条件特别容易翻车比如爬坡约束与启停状态的耦合、备用容量的分配逻辑以及碳交易项的线性化处理。下文我按建模顺序展开代码片段都是我实测能跑通的版本。1. 低碳调度问题在写什么从碳排放约束到碳交易机制1.1 碳排放从哪儿来常规机组的煤耗特性与碳排模型低碳调度和传统经济调度的根本区别在于目标函数和约束中多了“碳排放”这条线。电力系统里碳排放的主要来源是火电机组燃气机组虽然也有排放但强度低一个量级风电和光伏在运行阶段基本可以按零碳处理。所以建模的第一步就是把常规机组的碳排放量和出力之间的关系写清楚。实际工程中火电机组的煤耗特性常用二次函数表示C_fuel a * P^2 b * P c其中P是机组出力a、b、c是煤耗系数。碳排放量和燃料消耗量基本成正比所以碳排量也可以写成类似的形式。为了方便优化求解通常是取二次函数因为yalmip配合gurobi/cplex可以直接处理凸二次项。如果你用的是线性求解器或者想加快求解速度也可以把二次函数分段线性化但代价是需要引入额外的连续变量和约束代码量会增加不少。我在项目中采用的碳排模型更简单直接——认为碳排量和出力近似线性关系加上一个空载排放系数E_i(t) alpha_i * P_i(t) beta_i * u_i(t)这里u_i(t)是机组i在时段t的启停状态0/1变量beta_i表示机组空转时的基础排放。虽然精度比二次模型略低但在低碳调度场景里调度的核心是“让谁多发电、让谁少发电”线性模型已经能反映机组之间的碳排强度差异而且求解稳定性好很多。如果你想追求更精确的结果改成二次模型在代码上也只是加一个quadratic term的事不影响整体框架。1.2 碳交易机制如何进优化模型低碳调度里最常被问到的就是“碳交易怎么建模”。这个其实不复杂核心就一个公式碳交易成本 碳价 * (实际碳排放量 - 免费碳配额)如果实际排放量低于配额右边是负的相当于系统通过减排赚到了卖碳的收益如果超过配额就需要花钱买碳配额。把这个项放进目标函数调度模型就会自动权衡——“多花点煤钱让低碳机组多发”和“多买碳配额让高碳机组多发”之间求解器会选总成本更低的那条路。关键参数的设置需要注意参数作用典型取值参考碳价决定碳约束的惩罚强度50-100元/吨随政策场景变化免费配额决定系统的减排压力按历史排放量或装机容量比例分配碳排系数影响各机组被调用的优先级燃煤0.8-1.0 kg/kWh燃气0.4-0.5 kg/kWh我一开始把碳配额设置得很宽裕结果发现碳交易项对整个调度结果几乎没影响模型退化成普通经济调度了。后来把配额收紧到比无损调度的碳排放量低10%左右碳价项才开始真正起作用——这是低碳调度代码调试中很值得注意的一点。1.3 低碳调度和传统经济调度的本质区别用大白话说传统经济调度是“怎么发电最便宜”低碳调度是“怎么发电既便宜又少排碳”。但这个区别不是在目标函数里简单加一个碳价系数就行它会连锁改变系统的运行方式实际表现有三点第一高碳机组的出力会被压低。即使某些煤电厂的发电成本很低只要碳价够高调度结果会倾向让燃气机组甚至部分弃风来替代它。第二机组的启停策略会变化为了减少空转排放模型可能让部分机组在低谷时段直接停机而不是低负荷运行。第三风电的消纳优先级会提高因为风电不仅边际成本为零碳排也为零。这个区别其实提醒了我们低碳调度代码里的关键不是碳排约束本身有多复杂而是它在整个目标函数中的权重设置是否合理是否真的影响了调度决策。这也是我下面代码里目标函数设计的一个核心考量。2. 源荷不确定性的数学表达与场景生成2.1 风电出力不确定性的两个层面风电接入带来最本质的问题就是“不确定性”你预测的风电曲线和实际曲线永远对不上。从建模角度说这种不确定性体现在两个层面一是预测误差的分布特性二是极端场景下风电大幅波动对系统备用容量的冲击。经典的表达方式是假设风电预测误差服从正态分布以预测值为期望标准差取预测值的10%-20%P_wind_real P_wind_forecast epsilon_wind epsilon_wind ~ N(0, sigma_w^2)但正态分布有一个问题它的尾部太薄很难刻画风电“突然从满发掉到零”的极端情况。实际工程里我见过用β分布、拉普拉斯分布甚至直接用历史误差统计的。在yalmip框架里这些分布本身不影响建模影响的是你怎么把分布信息转化为约束或场景。对于这个项目我采用了场景法。所谓场景法就是通过蒙特卡洛抽样生成大量可能的风电出力曲线再对场景进行削减用少数有代表性的场景代替连续分布。这样做的好处是模型可以直接写成确定性的混合整数规划用gurobi求解即可。2.2 负荷预测误差怎么处理负荷不确定性经常和风电不确定性放在一起处理所以叫“源荷不确定”——源是电源荷是负荷。负荷的预测误差同样可以用正态分布刻画P_load_real P_load_forecast epsilon_load epsilon_load ~ N(0, sigma_l^2)相对风电而言负荷预测误差的波动要小一些典型标准差是预测值的3%-5%。但它的影响不能忽略因为功率平衡约束必须时刻满足如果负荷实际值高于预测值就需要有额外的向上备用。因此在调度代码里负荷不确定性一方面通过随机场景影响功率平衡另一方面通过备用需求约束来保证系统安全。我实际跑下来的经验是负荷误差和风电误差可以在同一个场景里耦合抽样它们的相关系数设为0就可以因为从实际数据看负荷高峰期往往伴随着风电出力的某些统计特征但短期预测误差层面的相关性很弱设为0不会带来明显偏差。2.3 场景生成与削减从1000个场景到10个场景场景生成的流程我在代码里分三步第一步用蒙特卡洛抽样生成N个风电出力场景和N个负荷场景N通常取500到1000。第二步用同步回代削减法fast backward reduction把N个场景削减到S个代表性场景S根据你期望的计算时间取5到20。第三步为每个场景计算一个概率权重所有权重之和为1。场景削减的逻辑不复杂核心思想就是“保留概率大、空间上差异大的场景去掉概率小或者跟其他场景很像的场景”。MATLAB实现时可以调用函数或者自己写一个简单的贪婪算法每次计算所有场景对之间的距离删除一个与已保留场景最近且概率最小的场景然后把它的概率加到最近邻场景上。这里想强调一下场景数量的选择。很多初学者一上来就保留20个甚至30个场景结果模型规模爆炸求解器几分钟都算不完。我测试下来对一个3台火电机组加一座风电场的测试系统10个场景和30个场景的调度结果差异很小目标函数值差异一般不超过2%但求解时间可能差好几倍。所以先跑10个场景调通逻辑再根据精度需求增加是比较务实的路径。2.4 随机优化与鲁棒优化的取舍写代码之前还要做一个路线选择用随机优化还是鲁棒优化。这两种方法的数学表达不同代码结构也不同。随机优化就是我上面说的场景法的目标是最小化所有场景下总成本的期望值它给出的是一个“平均意义”上的最优方案。鲁棒优化则是在最恶劣场景下寻找最优方案保证无论风电和负荷怎么波动系统都能安全运行代价是运行成本偏高。如果Title里明确写了“源荷不确定”两种方法都算关键看你要解决什么问题。我的建议是做数据分析或学术对比用随机优化做工程保守设计用鲁棒优化。后面代码部分我会以随机优化为主线因为它更直观而且能直接给出“风电消纳率”“碳排期望值”这种大家关心的指标。鲁棒优化的代码我最后也补充了一些思路。3. yalmip建模的核心代码骨架与变量设计思路3.1 环境准备与求解器配置先说一下环境。MATLAB yalmip 只是建模环境真正干活的是底层求解器。低碳调度问题本质是混合整数线性规划MILP或混合整数二次规划MIQP所以你需要一个能处理整数变量的求解器。我实测推荐组合是MATLAB R2019b及以上 yalmip gurobi或cplex。gurobi对整数规划的求解速度明显优于cplex而且学术许可证免费。装好之后在MATLAB里验证是否可用% 设置求解器路径 addpath(D:\gurobi\gurobi952\matlab); % 改成你自己的路径 savepath % 验证yalmip能找到gurobi yalmiptestyalmiptest这个命令会检查yalmip已安装的求解器列表如果gurobi出现在列表里就说明配置成功。如果是用cplex同理检查cplex是否被识别。3.2 决策变量的定义连续变量与二进制变量的分工这是yalmip建模中最核心的一步。变量设计得好不好直接决定约束能不能写清楚、模型能不能快速求解。低碳调度模型中的变量大致分三类调度变量连续变量P(i,t)机组i在时段t的出力P_wind(t)风电场在时段t的并网功率注意不一定是可用风电功率P_load_real(t)实际负荷需求场景相关状态变量二进制变量u(i,t)机组i在时段t的启停状态1为开机0为停机辅助变量场景相关每个随机场景下的功率偏差、备用容量等用yalmip定义变量的标准写法如下%% 定义变量 % 基础参数 I 3; % 火电机组数 T 24; % 调度时段24小时 S 10; % 场景数 % 决策变量 P sdpvar(I, T, full); % 机组出力 u binvar(I, T, full); % 机组启停状态 P_wind sdpvar(1, T, full); % 风电上网功率 P_ls sdpvar(1, T, full); % 弃风功率 % 场景相关变量每个场景下各机组出力 风电实际消纳 P_scn sdpvar(I, T, S, full); % 场景s下机组i的出力 P_wind_scn sdpvar(1, T, S, full); % 场景s下的风电消纳这里有一组容易混淆的变量P、P_scn、P_wind、P_wind_scn。我的设计逻辑是P是基准场景预测场景的调度计划P_scn是各随机场景下的功率平衡修正。在随机优化模型中通常只对第一阶段的决策变量P、u做所有场景共同约束第二阶段变量P_scn允许随场景变化。这就是两阶段随机规划的基本思想。3.3 目标函数怎么拼装低碳调度的目标函数由四部分组成经典发电成本二次煤耗成本 sum(aP^2 bP c)机组启停成本碳交易成本正或负弃风惩罚项为了保证风电优先消纳代码中我用sdpvar和sum拼接目标函数格式如下%% 目标函数 Cost 0; % 1) 常规机组发电成本二次项可被gurobi处理 for i 1:I Cost Cost sum(a(i) * P(i,:).^2 b(i) * P(i,:) c(i) * u(i,:)); end % 2) 启停成本 Cost Cost sum(sum(C_start * max(0, diff([u(:,1), u], 1, 2)))); % 注意启停差分写法 % 3) 碳交易成本 E_total sum(sum(alpha * P beta .* u)); % 总碳排放量 E_quota quota_total; % 免费碳配额 Cost_Co2 carbon_price * (E_total - E_quota); Cost Cost Cost_Co2; % 4) 弃风惩罚项 Cost Cost wind_penalty * sum(P_wind_max - P_wind); % 设置目标为最小化 optimize(Constraints, Cost, sdpsettings(solver,gurobi,verbose,2));其中启停成本那一行需要单独解释一下。diff([u(:,1), u], 1, 2)的实际效果是对每一台机组取相邻时段的启停状态差值如果从0变1差值为1表示冷启动计一次启动成本。这种写法比用循环更简洁而且yalmip可以直接处理max函数。3.4 约束组装与求解指令yalmip的约束组装就是把所有等式和不等式用方括号拼接最终传给optimize%% 约束 Constraints []; % 功率平衡基准场景 Constraints [Constraints, sum(P, 1) P_wind - P_ls P_load_forecast]; % 机组出力上下限与启停耦合 for i 1:I Constraints [Constraints, P_min(i) * u(i,:) P(i,:) P_max(i) * u(i,:)]; end % 风电有功上限 Constraints [Constraints, 0 P_wind P_wind_max]; % 备用约束向上备用 Constraints [Constraints, sum(P_max .* u - P, 1) (P_wind_max - P_wind) reserve_req]; % 求解 ops sdpsettings(solver,gurobi,verbose,2); result optimize(Constraints, Cost, ops);注意批量约束的写法MATLAB的矩阵运算在这里很方便不需要写循环。但有一个坑当约束里出现P_max(i) * u(i,:)这类混合矩阵时一定要确保矩阵维度匹配否则yalmip会报维度错误。4. 约束条件构建的工程细节爬坡、备用、网络约束怎么落地4.1 机组出力上下限与启停逻辑约束condition最容易出错的地方就是变量之间的耦合逻辑。机组出力上下限这一条如果写成P_min(i) * u(i,t) P(i,t) P_max(i) * u(i,t)等于同时实现了两个功能机组开机时出力在[P_min, P_max]内、机组停机时出力被强制为0。这个写法比分开写P(i,t) 0和P(i,t) P_max * u(i,t)更严谨因为它同时约束了出力下限防止模型把一台开着机的机组出力压到接近0来节省燃料成本——现实中机组有最小技术出力低于这个值机组会不稳。我在第一次写代码时忽略了P_min的耦合结果模型给出的解里有一台机组出力只有5MW明显偏离实际。后来加上下限约束后才收敛到合理结果。4.2 爬坡约束的正确写法爬坡约束是调度模型里最容易把模型写“死”的地方。大部分教材里的爬坡约束长这样P(i,t) - P(i,t-1) R_up(i) % 向上爬坡限制 P(i,t-1) - P(i,t) R_down(i) % 向下爬坡限制但这个写法有个大问题如果t-1时段机组是停机状态P(i,t-1) 0那么t时段开机时出力的变化量可以直接达到P_max因为0到P_max的“爬坡”超过了R_up会被第一个约束卡死。也就是说这种写法会阻止机组从停机状态开机或者从开机状态停机。正确的做法是给爬坡约束加上启停状态的修正项for i 1:I for t 2:T % 向上爬坡考虑从停机状态启动的情况 Constraints [Constraints, P(i,t) - P(i,t-1) R_up(i) * u(i,t-1) P_max(i) * (1 - u(i,t-1))]; % 等价地如果上一时段停机右端项为P_max(i)允许直接启动 end end这个约束的逻辑是如果u(i,t-1)1上一时段开机右端项等于R_up(i)限制正常爬坡如果u(i,t-1)0上一时段停机右端项等于P_max(i)等于把爬坡限制放开了允许机组从0直接启动到任意出力。类似的向下爬坡也要加一项P_max * (1 - u(i,t))考虑本时段停机的情况。这个细节不处理好模型很容易报“infeasible problem”。4.3 系统功率平衡与备用约束功率平衡约束是整个模型的基本骨架。在确定性模型中写sum(P) P_wind P_load就够了但在随机优化框架下要更细致基准场景功率平衡这个约束用预测值决定主要的调度方案。场景功率平衡每个随机场景下要允许第二阶段变量进行调整。备用约束保证系统在风电/负荷波动时能够通过爬坡等方式应对。备用约束是相对容易被忽略的一条。它的物理含义是系统在某一时刻能够“多发的功率”要大于预测误差可能带来的缺口。我用向上备用为例% 向上备用约束 Constraints [Constraints, sum(P_max .* u - P, 1) (P_wind_max - P_wind) reserve_require];P_max .* u - P的含义是每台机组还可以多发的空间P_wind_max - P_wind是风电还可以多发的空间。备用的需求量通常取系统最大负荷的5%-10%再叠加一个风电出力预测误差的百分比比如风电装机容量的10%。如果备用需求设得过高成本会明显上升设得过低极端场景下可能出现功率失衡。建议先跑一版不加备用约束的模型看系统最大波动的量级再反推备用需求值。4.4 网络安全约束直流潮流法的简化应用如果你做的是单节点系统所有机组和负荷都挂在一个母线上功率平衡就够了。但实际电网是有网络结构的线路传输容量会限制调度结果。这时候需要加直流潮流约束。直流潮流法是把交流潮流做线性化处理标准的表达式是PL(k) sum over all buses of PTDF(k, n) * (P_gen(n) - P_load(n))其中PTDF是功率传输分配因子矩阵。在yalmip里写成矩阵约束很容易% 假设已知线路潮流转移因子矩阵PTDF P_line PTDF * (P_gen_vector - P_load_vector); % 各线路潮流 Constraints [Constraints, P_line_min P_line P_line_max];不过要注意如果模型节点数比较小比如IEEE 30节点PTDF矩阵通常是稠密的直接相乘会增加不少约束表达时间。我建议把PTDF计算放在模型搭建之前作为常量矩阵传入不要每次迭代都重新计算。另外如果需要考虑网损直流潮流法就有些不够用了可以考虑直流潮流加网损修正但那会让模型变成二次约束规划求解复杂度上升不少。5. 求解效率与结果分析谁卡住了你的求解器5.1 模型求解速度的优化技巧跑通模型不难跑快是另一回事。我实测的调参经验如下第一个影响求解速度的因素是场景数。场景从5个加到20个求解时间可能从几秒涨到几分钟而且这个增长不是线性的。如果你只做趋势分析5个场景就够了如果要做论文级别的精度对比控制在10-15个场景是效率平衡点。第二个因素是大M参数的选择。在机组启停、备用约束里经常需要用到逻辑约束写的时候会用到M这个常数。M取得太大比如1e6会让求解器的数值稳定性变差M取得太小又会把可行域错误地截断。我给个经验值M一般是相关变量取值范围的上界比如机组出力上限P_max最多再乘1.1的安全系数不要拍脑袋取一个超大数。yalmip中有一个implies函数可以直接处理逻辑约束可以避免手写大M极力推荐使用Constraints [Constraints, implies(u(i,t), P(i,t) P_min(i))];第三点是尽量少引入额外的整数变量。整数变量是MIP求解慢的根源。能用连续变量表达的逻辑就不要用整数能用一个0/1变量表达的不要用两个。比如“机组启停”这个状态用一个u变量就够了不需要额外引入y_on和y_off两个变量。5.2 结果分析思路模型求解完之后光拿到目标函数值是不够的需要从调度结果中提取信息判断模型是否正确。我通常按下面几步检查第一步看功率平衡是否满足。把各时段的机组出力风电出力弃风功率和负荷预测值进行比较如果偏差超过1e-6说明约束有问题。第二步看碳排放总量和碳配额的差距这个比值反映了碳价在多大程度上影响调度决策。第三步看风电消纳率如果弃风率异常高比如超过30%很可能是备用约束设得过于严格或者然后机组的最小出力总和太高。第四步看每台机组的出力曲线检查是否存在频繁启停或不正常的出力跳变。结果可视化我常用两个图一个是24小时机组出力堆叠图一个是各场景下总成本概率分布直方图。堆叠图能直观看到火电、风电在每个时段的分工直方图能反映不确定性对总成本的影响范围。5.3 常见错误与调试经验这一节专门说说我跑这类模型时遇到的高频报错。最常见的错误是Infeasible problem。可见约束冲突我会先检查爬坡约束是否忘了加启停修正项见4.2再检查备用约束是否过于苛刻最后检查碳配额相关的目标函数是否导致模型无界。第二个常见错误是求解器返回NaN或Inf。这种情况大多是因为某个变量没有被约束完全比如机组启停变量u和出力变量P之间只有上限耦合、没有下限耦合导致某些机组的P可以取负值在sdpvar默认没有非负约束时。解决方法是给所有物理量显式加上界约束不要依赖隐式的非负性。第三个经典坑是目标函数里出现了非线性函数而求解器不支持。比如用yalmip写max(P, 0)或abs(P)时yalmip会自动引入辅助变量大多数情况下没问题但如果用了exp、log这类非线性函数就会变成NLP问题gurobi无法处理报错提示通常很隐晦。低碳调度里凡是涉及逻辑函数的地方尽量用线性改写实在不行用分段线性化。关于这个模型往后还能怎么扩展——我自己下一步可能会加上储能设备因为储能的加入会引入“时序耦合”的约束让模型更复杂也更有实际价值。另外碳价如果做成随时间变化的波动序列结果也会比固定碳价更贴近现实。这套MATLAByalmip的框架都能比较自然地延展过去变量多加一组约束多加几行整体逻辑不需要推倒重来。最后分享一个实操细节跑正式的算例之前一定要先跑一个2机、24时段、2场景的最小模型验证逻辑再逐渐扩大规模。我吃过一次亏直接从10机组、IEEE 39节点、20场景起步报错之后定位到问题是爬坡约束的启停耦合但在大模型里排查这个问题的成本远高于在迷你模型里的成本。先小后大是这类调度代码最稳的调试路径。