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为例流程就是对MC生成的 (S) 个源荷联合场景进行聚类把相似场景合并成一个典型场景。聚类的特征向量可以是“各时段风电出力各时段负荷”拼接成的高维向量。用轮廓系数Silhouette Coefficient或者肘部法则选择最优聚类数 (K)。聚类的中心作为代表场景每个类内场景数量占比作为该场景的权重。这样处理后几千个场景压缩成5~10个典型场景求解效率直接质变。我在代码里用的是K-means实现简单效果也稳定。同步回代的效果其实更好但代码实现复杂一些适合做进阶研究。提示场景削减的时候要特别注意削减后场景的总体统计特性均值、方差、相关性会发生变化。做完削减后最好做个检验对比一下削减前后的风电出力和负荷的期望值曲线偏差超过5%就要考虑增加聚类数。2. 低碳调度模型碳排放约束是怎么“嵌”进去的低碳调度的“低碳”二字不是把碳排放当成事后统计指标而是要把排放约束引入优化模型让调度决策朝着低碳方向主动调整。这里涉及碳交易机制、碳捕集设备、低碳机组组合等多个维度的建模选择。2.1 目标函数经济成本 碳成本的双层结构传统经济调度的目标函数是发电成本最小化。低碳调度通常在此基础上叠加碳排放成本变成一个双目标、或者单目标加权的问题。最常用的做法是用碳交易价格把碳排放量折算成成本项从而合并成一个单目标[ \min ; C C_{fuel} C_{carbon} C_{curtail} ]其中 (C_{fuel}) 是机组燃料成本用二次函数拟合[ C_{fuel} \sum_{t1}^T \sum_{i1}^N \left( a_i P_{i,t}^2 b_i P_{i,t} c_i \right) ](C_{carbon}) 则是碳排放成本[ C_{carbon} \lambda_{CO2} \sum_{t1}^T \sum_{i1}^N \left( \delta_i P_{i,t} - E_{free,i} \right) ]这里 (\lambda_{CO2}) 是碳交易价格元/吨(\delta_i) 是机组 i 的碳排放强度吨/MWh(E_{free,i}) 是免费碳排放配额。当实际排放超过配额时需要购买碳排放权低于配额则可以出售获利。这个机制做出来的效果就是低碳机组因为碳成本低在调度中会更有竞争力高碳机组如果碳价足够高就可能被压到最小出力甚至停机。如果考虑阶梯碳价碳价随排放量分段递增效果更贴近实际但模型会引入整数变量复杂度上一个台阶。初学者建议先做线性碳价跑通了再升级。2.2 约束条件不是列完就完事了约束条件通常包括系统功率平衡、机组出力上下限、爬坡约束、旋转备用约束、碳排放上限约束等。功率平衡约束是所有调度模型的核心[ \sum_{i1}^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_{i1}^N (P_{i}^{max} - P_{i,t}) \geq R_{up} \xi_{W} \cdot P_{W,t} ][ \sum_{i1}^N (P_{i,t} - P_{i}^{min}) \geq R_{down} \xi_{L} \cdot P_{L,t} ]典型做法是按“负荷的5%~10% 风电出力的15%~20%”来设定备用需求。这样风电出力越大的场景需要预留的上调备用越多确保风电突然掉出力时系统不会失负荷。碳排放约束一般用系统碳排放总量上限来控制[ \sum_{t1}^T \sum_{i1}^N \delta_i P_{i,t} \leq E_{total}^{max} ]这个约束的作用是在碳交易成本之外再叠加一个刚性排放上限。两种机制同时启用时碳价兜底、总量封顶模型会同时受经济和环保双轮驱动结果更接近真实碳约束环境下的调度行为。2.3 典型场景集下的运行约束处理既然场景削减后得到的是一个典型场景集合调度模型就需要保证方案在这些场景下都可行。最常用的方式是把功率平衡约束扩展成带松弛变量的形式或者对每个场景单独列约束。比如针对每个典型场景 (s)功率平衡写成[ \sum_{i1}^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 求解框架选型YALMIPCPLEX是黄金组合Matlab里求解带约束的非线性优化问题工具箱选择很关键。我强烈推荐用YALMIP做建模层配CPLEX或Gurobi做求解器。下面的MATLAB代码示例展示了如何配置YALMIP与CPLEX求解器。YALMIP的核心理念是你只需要声明“决策变量”和“约束条件”然后调用optimize函数求解它自动处理变量类型检测、线性化、求解器切换等脏活累活。相比直接用CPLEX的底层APIYALMIP的代码可读性和调试性都高一个档次后续要改模型复杂度也方便得多。我在论文复现里基本都是YALMIPCPLEX很少翻车。注意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); % 机组出力第二阶段Isnt决策 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这类问题的核心原因往往集中在约束上。我常用的排查步骤如下看Power Balance约束。这是低频翻车点。尤其是加入了备用约束后系统的可调功率上限可能真的撑不住这时候不是模型写错而是物理上无解。把备用约束调低或者给功率平衡加松弛变量试试。检查变量维度是否一致。YALMIP对维度比较宽容但矩阵不匹配时照样报错。比如P是 (N \times T)而P_load是 (T \times 1) 的列向量相加的时候维度对不上。检查约束列表是否为空。我自己就干过忘记把某条约束append进Constraint集合模型出于求解器的指令找不到约束而误报不可行的蠢事。用诊断工具。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的联合调度、引入需求响应资源、把确定性场景调度升级为两阶段分布鲁棒优化等。我个人的建议是先把基础模型吃透场景调度和不确定性建模搞扎实再逐步进阶。最后分享一个我在调试代码时养成的小习惯每次修改模型参数后先跑一个简化的小规模算例验证逻辑再跑全规模。比如先让 (K3)、(T6)、机组3台跑通跑快再逐步恢复完整配置。这样能极大减少虫子bug的排查时间把精力放在优化本身。这套Matlab代码如果能配合一组合理的参数其实复现性很强。整个项目的完整流程——从不确定性建模、场景削减、低碳目标构建到求解分析——如果你能独立走通一遍你对电力系统调度优化的理解会上一个台阶。
含风电的低碳经济调度：源荷不确定性建模与Matlab实现
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为例流程就是对MC生成的 (S) 个源荷联合场景进行聚类把相似场景合并成一个典型场景。聚类的特征向量可以是“各时段风电出力各时段负荷”拼接成的高维向量。用轮廓系数Silhouette Coefficient或者肘部法则选择最优聚类数 (K)。聚类的中心作为代表场景每个类内场景数量占比作为该场景的权重。这样处理后几千个场景压缩成5~10个典型场景求解效率直接质变。我在代码里用的是K-means实现简单效果也稳定。同步回代的效果其实更好但代码实现复杂一些适合做进阶研究。提示场景削减的时候要特别注意削减后场景的总体统计特性均值、方差、相关性会发生变化。做完削减后最好做个检验对比一下削减前后的风电出力和负荷的期望值曲线偏差超过5%就要考虑增加聚类数。2. 低碳调度模型碳排放约束是怎么“嵌”进去的低碳调度的“低碳”二字不是把碳排放当成事后统计指标而是要把排放约束引入优化模型让调度决策朝着低碳方向主动调整。这里涉及碳交易机制、碳捕集设备、低碳机组组合等多个维度的建模选择。2.1 目标函数经济成本 碳成本的双层结构传统经济调度的目标函数是发电成本最小化。低碳调度通常在此基础上叠加碳排放成本变成一个双目标、或者单目标加权的问题。最常用的做法是用碳交易价格把碳排放量折算成成本项从而合并成一个单目标[ \min ; C C_{fuel} C_{carbon} C_{curtail} ]其中 (C_{fuel}) 是机组燃料成本用二次函数拟合[ C_{fuel} \sum_{t1}^T \sum_{i1}^N \left( a_i P_{i,t}^2 b_i P_{i,t} c_i \right) ](C_{carbon}) 则是碳排放成本[ C_{carbon} \lambda_{CO2} \sum_{t1}^T \sum_{i1}^N \left( \delta_i P_{i,t} - E_{free,i} \right) ]这里 (\lambda_{CO2}) 是碳交易价格元/吨(\delta_i) 是机组 i 的碳排放强度吨/MWh(E_{free,i}) 是免费碳排放配额。当实际排放超过配额时需要购买碳排放权低于配额则可以出售获利。这个机制做出来的效果就是低碳机组因为碳成本低在调度中会更有竞争力高碳机组如果碳价足够高就可能被压到最小出力甚至停机。如果考虑阶梯碳价碳价随排放量分段递增效果更贴近实际但模型会引入整数变量复杂度上一个台阶。初学者建议先做线性碳价跑通了再升级。2.2 约束条件不是列完就完事了约束条件通常包括系统功率平衡、机组出力上下限、爬坡约束、旋转备用约束、碳排放上限约束等。功率平衡约束是所有调度模型的核心[ \sum_{i1}^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_{i1}^N (P_{i}^{max} - P_{i,t}) \geq R_{up} \xi_{W} \cdot P_{W,t} ][ \sum_{i1}^N (P_{i,t} - P_{i}^{min}) \geq R_{down} \xi_{L} \cdot P_{L,t} ]典型做法是按“负荷的5%~10% 风电出力的15%~20%”来设定备用需求。这样风电出力越大的场景需要预留的上调备用越多确保风电突然掉出力时系统不会失负荷。碳排放约束一般用系统碳排放总量上限来控制[ \sum_{t1}^T \sum_{i1}^N \delta_i P_{i,t} \leq E_{total}^{max} ]这个约束的作用是在碳交易成本之外再叠加一个刚性排放上限。两种机制同时启用时碳价兜底、总量封顶模型会同时受经济和环保双轮驱动结果更接近真实碳约束环境下的调度行为。2.3 典型场景集下的运行约束处理既然场景削减后得到的是一个典型场景集合调度模型就需要保证方案在这些场景下都可行。最常用的方式是把功率平衡约束扩展成带松弛变量的形式或者对每个场景单独列约束。比如针对每个典型场景 (s)功率平衡写成[ \sum_{i1}^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 求解框架选型YALMIPCPLEX是黄金组合Matlab里求解带约束的非线性优化问题工具箱选择很关键。我强烈推荐用YALMIP做建模层配CPLEX或Gurobi做求解器。下面的MATLAB代码示例展示了如何配置YALMIP与CPLEX求解器。YALMIP的核心理念是你只需要声明“决策变量”和“约束条件”然后调用optimize函数求解它自动处理变量类型检测、线性化、求解器切换等脏活累活。相比直接用CPLEX的底层APIYALMIP的代码可读性和调试性都高一个档次后续要改模型复杂度也方便得多。我在论文复现里基本都是YALMIPCPLEX很少翻车。注意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); % 机组出力第二阶段Isnt决策 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这类问题的核心原因往往集中在约束上。我常用的排查步骤如下看Power Balance约束。这是低频翻车点。尤其是加入了备用约束后系统的可调功率上限可能真的撑不住这时候不是模型写错而是物理上无解。把备用约束调低或者给功率平衡加松弛变量试试。检查变量维度是否一致。YALMIP对维度比较宽容但矩阵不匹配时照样报错。比如P是 (N \times T)而P_load是 (T \times 1) 的列向量相加的时候维度对不上。检查约束列表是否为空。我自己就干过忘记把某条约束append进Constraint集合模型出于求解器的指令找不到约束而误报不可行的蠢事。用诊断工具。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的联合调度、引入需求响应资源、把确定性场景调度升级为两阶段分布鲁棒优化等。我个人的建议是先把基础模型吃透场景调度和不确定性建模搞扎实再逐步进阶。最后分享一个我在调试代码时养成的小习惯每次修改模型参数后先跑一个简化的小规模算例验证逻辑再跑全规模。比如先让 (K3)、(T6)、机组3台跑通跑快再逐步恢复完整配置。这样能极大减少虫子bug的排查时间把精力放在优化本身。这套Matlab代码如果能配合一组合理的参数其实复现性很强。整个项目的完整流程——从不确定性建模、场景削减、低碳目标构建到求解分析——如果你能独立走通一遍你对电力系统调度优化的理解会上一个台阶。
