调度台前最怕的不是风电突然来一阵大波动而是我们根本不知道误差到底服从什么分布。第二天风电出力预测值是350兆瓦实际可能落在180到420兆瓦之间这种偏差的“分布形状”往往只有几十个历史样本谁也说不准。所以当“基于线性准则的考虑风力发电不确定性的分布鲁棒优化机组组合”这个思路进入视野时我最关心的是分布鲁棒优化能不能从论文变成能跑的代码。传统机组组合把风电当确定数字用随机优化则硬给误差套一个概率分布分布鲁棒优化走的是中间路线——只要真实分布藏在我用历史数据圈出的模糊集里调度方案就保证可行。本文就用Matlab把这个流程完整跑了一遍把建模思路、线性准则的作用、求解时踩过的坑都摆出来。1. 为什么机组组合必须正视风电误差的“分布未知”1.1 确定性调度误差不是噪声而是系统性风险很多教材里的机组组合模型输入侧只有一条风电预测曲线启停计划、出力基点、备用容量都围着这条单一曲线转。这种做法在风电占比小、误差不算离谱的时候还能接受但风电一多问题就藏不住了预测误差会在某些时段集中爆发比如夜间预测偏高而实际出力很低系统不得不用高价机组紧急顶上。更关键的是确定性模型的“备用设置”通常是拍脑袋的比如固定加10%备用或按最大单机容量留备用。它没有把“误差有多大、误差分布长什么样”和成本关联起来导致两种极端要么备用给少了负荷平衡被打破要么备用给多了经济性被白白牺牲。我见过一个很典型的案例某地区风电渗透率到30%之后固定比例备用策略在强风过程天气下连续三天出现备用不足最后是靠临时切负荷才兜住。问题根源不是预测算法太差而是调度模型根本没有量化“预测偏差的分布”对可行性的影响。1.2 随机优化与经典鲁棒优化的两难既然把风电当确定值不行自然会想到随机优化。随机优化要求先给出风电预测误差的确切概率分布比如已知误差服从均值为0、方差为σ的正态分布然后对大量抽样场景求期望最优。但实际中“分布是谁”本身就是未知数。文献里经常说“用一个正态分布拟合误差”就够了可你会发现同一风电场在不同季节、不同天气过程下的误差形态差别很大轻尾、重尾、偏态都可能出现。样本量不够时用指定分布做出来的解在样本外测试里频繁违约原因很简单你假设错了分布。反过来经典鲁棒优化把误差限定在一个确定的“盒式不确定集合”里只要误差落在这个盒子内约束必须全部满足。这个思路很稳但它等价于假设误差的每一个极端值同时发生。实际中不同风电场间误差存在明显抵消效应全网总偏差很少同时向最坏方向拉满。结果是鲁棒解的成本比确定性方案高出一大截启停机组频繁运行人员看了报价只想摇头。1.3 分布鲁棒优化给出的第三条路分布鲁棒优化的想法很自然别去精确指定分布也别用纯集合最坏情况而是在历史数据附近构造一个“模糊集”让真实分布以大概率落在这个集合里。调度方案只需要对该集合内所有可能的分布都可行或足够优。这个思路同时避开了两个坑第一不需要知道真实分布的具体形式只用样本第二通过模糊集的半径控制了保守程度半径取0退化成随机优化半径取无穷大退化成经典鲁棒优化。所以标题里“考虑风力发电不确定性”的关键不是把某个分布参数写得天花乱坠而是把“分布未知”这件事本身作为建模对象。我从一开始就确定这篇实现的核心难点在模糊集的构造、以及后续怎么把分布鲁棒问题转成Matlab能求解的混合整数线性规划。2. 模糊集与线性准则模型里真正值钱的数学设计2.1 把“分布未知”翻译成约束两阶段机组组合的一般形式机组组合天然是两阶段决策。第一阶段是“日前”决策决定每台机组在哪些时段开机、关机给出出力基点第二阶段是“实时/再调度”决策看到风电实际出力偏差之后通过调整部分机组出力、调用备用等手段保证负荷平衡和网络安全。用变量来写大致是这个结构第一阶段变量机组开停状态 (u_{i,t})、启动/停机变量、基点出力 (p^0_{i,t})。不确定参数各风电场在时段的出力偏差 (\xi_{w,t})可以按时段堆成一个向量 (\xi)。第二阶段变量实际调整量 (p_{i,t}(\xi))以及弃风、切负荷、备用调用等。问题的一般形式是[ \min_{u,p^0} \left{ c^T u d^T p^0 \sup_{\mathbb{P}\in \mathcal{F}} \mathbb{E}_{\mathbb{P}}\left[ Q(p^0,\xi) \right] \right} ]其中 (Q(p^0,\xi)) 是给定第一阶段方案和误差 (\xi) 后的再调度成本函数。难点就在 (\sup_{\mathbb{P}\in \mathcal{F}} \mathbb{E}_{\mathbb{P}}[\cdot]) 这一项它要求我们枚举模糊集内所有可能分布下的最坏期望成本直接算几乎不可能。2.2 用Wasserstein球圈出真实分布模糊集有很多种构造方式比如矩约束集、统计距离球、机器学习里的对抗样本集。本文采用了一个在工程上接受度很高的选择以历史经验分布为球心、以Wasserstein距离为半径构造模糊球。Wasserstein距离衡量的是“把一个分布搬运成另一个分布所需的最小成本”因为考虑到误差向量之间的坐标尺度问题实际计算时常选用1-范数或无穷范数等。它的好处是模糊集里不仅包含与经验分布接近的点分布还允许支撑集发生平移不会把某些真实可能出现的极端误差值从根上排除掉。设历史偏差样本为 (\hat{\xi}_1,\dots,\hat{\xi}_N)经验分布为 (\hat{\mathbb{P}}_N)则模糊集写为[ \mathcal{F}_\varepsilon{\mathbb{P}: W(\mathbb{P},\hat{\mathbb{P}}_N)\le \varepsilon} ]这里的 (\varepsilon) 是模糊集半径集中体现了模型对“分布不确定”的容忍程度。半径太小模型把历史样本当成金科玉律半径太大模型又宁可信最坏情况不信任任何数据。在这个框架下最坏期望问题可以通过对偶理论转变为一个增广的有限维优化问题。对于线性准则配合下的目标函数和约束最后得到的是混合整数线性规划这为Matlab下的求解扫清了最大障碍。2.3 线性决策规则把“看风下单”写成线性函数第二阶段最理想的调整策略是 (\xi) 的任意函数 (p_{i,t}(\xi))因为再调度完全跟随误差走。但任意函数是不可求解的工程中通行做法是限制函数形式其中最常用的是仿射决策规则也就是标题里说的“线性准则”。线性准则的含义非常直接机组出力对风电偏差的反应是线性响应写成[ p_{i,t}(\xi)p^0_{i,t}\sum_{w\in \mathcal{W}} \alpha_{i,t,w}\xi_w ]其中系数 (\alpha_{i,t,w}) 表示机组i在时段t对风电场w偏差的响应斜率。正值表示风电出力不足时多带出力负值表示风电出力过剩时少发或弃风。为什么敢用线性准则一方面是因为机组爬坡约束、成本函数在线性化后整个问题的最优调整策略在一定条件下本来就接近分段线性线性近似已经能捕获大部分经济性收益另一方面是只有把第二阶段策略设成线性函数对偶后的模型才能保持线性结构否则就要引入非线性规划求解难度会指数级上升。我当时做完第一版非线性场景测试后又用线性决策规则对比发现两者的期望成本差距在1%到3%之间但求解时间从几小时降到几分钟。做工程调度这个代价完全可以接受。3. Matlab实现的完整链路从历史样本到MILP落地3.1 数据准备从历史出力记录到偏差样本实现的第一步是整理风电预测偏差数据。每个风电场需要同时具备历史预测值和实际值确保数据在时区上对齐。时间颗粒度常见的是15分钟或1小时我建议用1小时起步先把模型跑通再细化。偏差样本的计算很简单% 假设 pred 是预测值矩阵(nHistorical x nWind) % actual 是实际出力矩阵(nHistorical x nWind) xi actual - pred; % 舍入到保留两位小数减少求解器数值压力 xi round(xi, 2);每个历史时刻对应一个偏差向量 (\xi_k)N个历史时刻就得到N个样本点。如果风电场数量很多建议先做相关性分析和主成分降维避免偏差向量维度过高导致对偶后的辅助变量爆炸式增长。我在实验中用了两个风电场、6个节点的小系统偏差样本取了最近90天的数据。注意样本量不必贪多但必须覆盖不同的天气过程和季节形态否则模糊集球心本身就有偏。3.2 上下层变量组织与目标函数构造模型变量主要包括三层第一层是整数变量描述机组启停状态例如6台机组、24个时段就有144个0-1变量。Matlab中可以用optimvar定义二进制变量。第二层是连续变量包括各机组各时段的基点出力、启动/停机成本相关辅助变量以及第二阶段调整策略的仿射系数 (\alpha_{i,t,w})。第三层是分布鲁棒对偶后的辅助变量这一步是模型能否落入MILP的关键通常包括对偶变量、范数约束里的辅助标量等。粗算下来6节点小系统大约有几千个连续变量和两百多个整数变量规模不算大。目标函数我分三块写% 燃料成本二次成本线性化后的分段系数 % 启停成本启动成本停机成本 % 分布鲁棒最坏期望再调度成本 objective sum(sum(bidcost .* p0)) ... sum(sum(startcost .* u_start)) ... dro_expected_adjust_cost;dro_expected_adjust_cost由对偶问题中的辅助变量构成求解器看不见“最坏分布”这些概念只能看见一堆线性变量和约束。3.3 对偶转换、线性化与intlinprog求解Wasserstein分布鲁棒问题对偶后指数项会被一组有限的变量和约束替换。对线性准则情形核心是把[ \sup_{\mathbb{P}\in \mathcal{F}\varepsilon} \mathbb{E}{\mathbb{P}}[Q(p^0,\xi)] ]转化为关于辅助变量的线性目标并在约束中加入关于每个历史样本 (\xi_k) 的不等式。在实际建模时我使用了Matlab的problem-based优化工具箱把约束一行行读进去再调用内置的混合整数线性规划求解器。关键代码示意如下prob optimproblem; u optimvar(u, nGen, T, Type, integer, LowerBound, 0, UpperBound, 1); p0 optimvar(p0, nGen, T, LowerBound, 0); alpha optimvar(alpha, nGen, T, nWind, LowerBound, -1, UpperBound, 1); s optimvar(s, nSample, 1, LowerBound, 0); % 对偶辅助变量 q optimvar(q, 1, 1); % 对偶标量变量 prob.Constraints.load_balance ...; prob.Constraints.dro_dual ...; prob.Objective ...;这里最容易出错的是范数约束的线性化。如果历史偏差按2-范数计算Wasserstein距离对偶后会出现二阶锥约束MATLAB自带intlinprog并不能直接处理。解决方式有两种要么把范数改成1-范数或无穷范数得到纯线性约束要么引入额外变量做锥规划近似。我建议第一版代码直接选1-范数干净、快速和2-范数在最坏成本上的差异很小。3.4 用一个小系统算例验证流程我用包含6台机组、24时段、2个风电场和90个历史偏差样本的模拟系统做验证。机组参数来自常见算例的公开形式风电预测曲线和实际出力由历史数据驱动。跑通之后第一步检查的是负荷平衡约束是否满足。因为分布鲁棒模型允许误差在一定范围内变化所以不能只查预测场景下的平衡还要查所有历史样本对应的平衡是否满足。把每个样本代入线性决策规则后逐点检查节点注入功率发现有一半时段因为风电功率离散化精度问题出现微小越限把偏差样本保留两位小数之后恢复正常。最终得到的启停计划相比确定性方案多启动了半台到一台机组但切负荷概率从确定性方案的10%以上降到3%以内总成本增加大约6%。在调度场景里这个代价换来的可靠性提升是值得的。4. 参数敏感性、求解性能与那些容易踩的坑4.1 模糊集半径ε怎么取过小等于自欺过大等于躺平分布鲁棒优化最难以解释也最需要调参的就是模糊集半径 ε。它本质上是“样本容量和真实置信度之间的扳手”。理论上可以通过统计方法给出ε关于样本数量N、置信度β的估计公式但工程上我更推荐直接用样本外测试定半径。做法是留出一部分历史数据不参与建模把不同ε下求得的调度方案代入样本外数据统计约束违反率和总成本。我实测的一组数据大致如下注意这里不是某个系统的标准答案只用于说明趋势ε取值总成本万元求解时间秒样本外约束违反率0.0110254218.7%0.051089596.5%0.101153872.1%0.3012481240.3%从表里能清楚看到ε0.01时模型约等于把历史样本当成精确分布样本外违反率飙升到接近两成这个方案根本不敢用ε0.30时样本外表现很稳但成本涨了20%多调度员也会抱怨太保守。我的选择习惯是先根据安全标准确定可接受违反率上限比如切负荷概率低于5%然后选能满足该约束的最小ε。这比纠结统计置信度公式更直观也更贴近调度要求。4.2 大M与辅助变量带来的数值陷阱分布鲁棒对偶转换后会引入很多形如 (\eta_k \ge \ldots) 的线性约束其中经常出现表示“样本k最坏成本”的辅助变量和二进制选择逻辑。为了把逻辑关系写进整数规划我用了大M法。大M法本身不复杂但M取多大会直接影响求解质量。M太小会把可行域错误压缩导致明明应该有解的时节电计划直接无解M太大则会让lp松弛很差整数求解时间暴涨。一个真实的排查经历我把启动成本的线性化M值设成1e6结果求解器在连续节点上反复探索无用区域跑2000秒还收敛不了。把M值缩小到和该时段最大可能成本同量级后问题在几十秒内解决。经验是M值不是越大越好最好通过历史数据推算出成本上界再放大1.5到2倍留一点余量即可。数值上还要注意向量维度假象。Matlab的problem-based模型会把每个约束自动向量化当你写sum(p, 1) demand时求解器内部生成的约束数量可能远大于你预期导致内存暴涨。我建议把24时段按循环逐条写清楚或者用索引生成约束矩阵不要依赖过度花哨的向量化技巧。4.3 求解时间与精度权衡的实测明细分布鲁棒机组组合的求解时间主要由整数变量规模和对偶辅助变量规模共同决定。我实测发现同等规模下把再调度策略从线性准则放松成普通连续变量再做场景枚举求解时间会呈指数级增长而线性准则将每个机组的调整策略压缩成少量仿射系数整数变量数量基本不变瓶颈从“如何枚举分布”转移到“如何切分整数节点”对商用求解器比较友好。另一个容易忽略的性能因素是风电偏差向量的维度。两个风电场时仿射系数矩阵是 机组数×时段数×风电场数还算可控如果系统变成20个风电场系数数量立刻膨胀到原来的10倍单纯增加风电场数量对模型负担非常大。针对这个情况我采用了两种缓解手段一是对高度相关的风电场进行聚类合并成若干等效风电场二是对偏差样本做稀疏编码只保留对负荷平衡影响最大的分量。这样可以在损失极小精度的情况下显著压减变量规模。5. 往储能与滚动调度扩展时我踩出的几条经验模型跑通之后我自然想把它往含储能的联合调度方向推。第一个直观做法是把储能充放电功率也放进第二阶段线性策略里让储能随风电偏差自动调整。试算后发现虽然储能提升了灵活性但储能自身的充放电效率、容量衰减等参数也存在不确定性如果只把风电偏差放进模糊集储能参数误差带来的风险会被低估。一个可行的改进是把储能效率波动也放入随机向量模糊集半径相应增大。代价是求解规模进一步变大但结果明显更扎实。另一个经验来自滚动调度场景。如果每15分钟就重新估计一次模糊集ε会随着最新样本进进出出发生抖动导致相邻两个时段的启停方案互相冲突。我最终选择了一个折中短期内不调整ε只在每天重新滚动时更新一次历史样本并重新标定半径。如果你也想复现这套流程别急着把模糊集半径、范数类型、辅助变量全部调成“理论上最优”。先把一个最小系统的完整链路跑通把负荷平衡、爬坡约束、直流潮流约束都检查一遍再逐步增加风电场数和样本量。数据驱动的调度模型最关键的不是公式有多漂亮而是每一个约束在Matlab里都真实、可复现、敢拿到明天的调度台上运行。