分布鲁棒优化与线性决策规则在机组组合中的Matlab实现

风力发电出力的随机波动一直是机组组合调度里最让人头疼的问题之一。以前用确定性模型的时候我吃过不少亏用预测值做计划一旦实际风和预测差得多了要么备用不够被迫切负荷要么机组启停太慢导致弃风。后来转到随机优化又发现场景数量稍微多一点计算时间就蹭蹭往上涨而且概率分布稍微设得不对优化结果就名不副实。这几年我开始转向分布鲁棒优化Distributionally Robust OptimizationDRO结合线性决策规则来做机组组合效果和稳定性都好了不少。这篇文章就把这套做法——从问题拆解、数学模型、线性准则适用逻辑到Matlab代码实现和调试经验——完整地梳理一遍给同样在电力调度和优化领域踩坑的同行一个参考。这套内容适合谁主要是做电力系统优化调度、微电网能量管理、新能源并网研究的工程师和研究生。你需要具备Matlab基础最好用过Yalmip这类建模工具箱如果对鲁棒优化或者随机规划有一定了解就更容易上手。不过就算这些基础稍弱我也会把思路和代码细节尽量拆开讲清楚照着做也能复现出基本的分布鲁棒机组组合算例。1. 思路拆解为什么选分布鲁棒优化线性准则又是什么1.1 风力发电不确定性的三个真实难点做机组组合时风电不确定性本质上包含三个层次。第一个是预测误差的分布未知也就是说我们只知道历史数据或者预测模型给出的期望值但真实出力在预测值上下怎么波动、波动幅度的概率分布长什么样我们根本没有精确信息。第二个是时空相关性风电场相邻地理位置之间的出力存在明显耦合同一时刻不同风电场不会同时出大风也不会同时完全静止这种相关性如果不建模优化出来的备用分配就会不合理。第三个是极端场景的尾部风险风电出力低概率但高影响的情景比如连续多天大风导致出力骤降往往才是系统安全真正受威胁的地方。针对这三层难点传统做法往往顾此失彼。确定性方法完全忽略不确定性只加一个固定备用比例最简单但风险最高。随机优化需要假设风电出力的概率分布已知然后抽样生成大量场景理论上最优但分布假设错误时结果同样不可靠而且场景规模一大整数变量爆炸求解速度没法接受。鲁棒优化则只要求不确定量落在一个集合内不考虑分布保底性强但结果保守经济性差。分布鲁棒优化正好介于两者之间它假设真实分布属于一个以经验分布为中心的“分布集合”优化目标是最坏情况下的期望成本既能利用历史数据和预测信息又不必精确指定分布还能控制尾部风险。这个特性让它特别适合风电并网这种“数据有、但分布不确定”的问题。1.2 线性准则的核心含义把决策规则限制成线性函数所谓线性准则在这类问题里通常有两层意思。第一层是风险度量的线性化比如用CVaR条件风险价值或者均值-CVaR加权形式来做目标CVaR本身是线性规划可表达的这样目标函数就保持线性特征方便后续对偶转化。第二层更关键是线性决策规则Linear Decision RuleLDR也就是第二阶段的再调度决策被限制为不确定参数的仿射函数。为什么需要这个限制因为两阶段分布式鲁棒机组组合问题本质上是半无限规划第二阶段对每一个可能的风电出力场景都要满足功率平衡和机组爬坡约束决策变量是一个以风电出力为自变量的函数维度无穷大直接求解根本不现实。线性决策规则的核心想法是把“对任意出力值做出最优调整”这个复杂映射简化成“调整量 初始计划量 线性系数 × 出力偏差”这样简单的形式。这样一来函数优化就变成了有限个线性系数的优化问题规模大幅下降可以用商业求解器直接处理。代价是理论上可能会损失一部分最优性但实际工程中线性决策规则对机组组合这种近似线性系统来说效果足够好计算量却能降低几个量级。打个比方就好理解了。你要根据明天天气预报决定带多少衣服随机优化是给每种天气都准备一整套完整方案精细但复杂纯鲁棒是照着最冷最暴雨的情况全副武装保守但浪费线性决策规则则是“温差每大1度加一件外套”简单直观还兼顾了大多数情况。1.3 三种建模路线的横向对比为了看得更清楚我把确定性、随机优化、传统鲁棒、分布鲁棒这四种方法对同一个六节点系统的机组组合结果做了对比。方法不确定性建模计算成本保守程度数据需求确定性忽略极低最低仅预测值随机优化精确分布场景抽样高低需要准确分布传统鲁棒优化不确定集合盒式/椭球中高集合边界分布鲁棒优化矩不确定集/ Wasserstein球内的所有分布中高中等可控历史数据/场景样本这个表里分布鲁棒的数据需求很有意思它不要求你知道真实分布只需要一组历史样本或者预测误差样本就能构建模糊集。这一点在实际工程中是决定性的优势因为风电功率预测系统给你的往往是点预测和一组历史误差而不是一个漂亮的正态分布。下面我就沿着这个思路详细展开数学模型是怎么搭的以及线性准则在里面具体怎么落地的。2. 数学模型与关键转化从模糊集构造到对偶化简2.1 第一步构造风电不确定性的模糊集在分布鲁棒优化里最核心的一步是定义模糊集也就是“真实分布属于哪个范围内”。我常用的有两种方案各有各的适用场景。第一种是基于矩信息的模糊集。它只约束分布的均值在预测值附近的某个区间内二阶矩协方差在一定范围之内其余不确定。矩模糊集的优点是建模直观约束数量少对偶转化也比较成熟缺点是它允许的分布范围比较大而且没有利用到样本的更多形状信息结果往往偏保守。第二种是Wasserstein球模糊集。它是以经验分布为球心用Wasserstein距离定义半径真实分布被认为在以经验分布为圆心、半径可控的圆球内。相比矩模糊集Wasserstein球能更充分地利用样本信息球半径还可以与置信度联系起来当样本量增大时半径按1/N的规律收缩体现“数据越多、分布越确定”的自然逻辑。在实际操作中我推荐用Wasserstein球原因很简单调半径就是调保守程度物理含义清晰。半径设成0分布鲁棒优化退化为样本均值优化接近随机优化半径设得很大它趋近于传统鲁棒优化。你要做的只是在经济性和鲁棒性之间拨一个旋钮而不是去猜分布长什么样。这一点在工程项目里尤其有价值。2.2 两阶段机组组合模型的目标函数与约束我采用的模型是经典的两阶段结构。第一阶段是机组组合决策决定每台机组在哪些时段开机和停机这个阶段必须在知道风电实际出力之前完成所以是“预决策”包含机组启停状态、开机和停机动作以及基准出力计划。第二阶段是经济调度决策当风电出力被观测到后在第一阶段决策的基础上调整各机组出力、判断是否需要切负荷和弃风以保证功率平衡这一阶段是“实时调整”。目标函数写成最小化第一阶段启动成本、运行成本与第二阶段期望调整成本之和。在分布鲁棒框架下第二阶段成本不是对某个固定分布求期望而是对模糊集内最坏情况分布求期望再加上一个加权CVaR项来控制尾部风险。引入线性准则后第二阶段调整量写成关于风电出力偏差的线性函数目标函数中的期望和对偶转化就变成有限维线性规划的可处理形式。我在这里强调一个建模细节机组组合问题中的0/1整数变量必须留在第一阶段第二阶段的连续变量做线性决策规则化这样才能把两阶段问题整体转成混合整数线性规划MILP。如果你尝试对整数变量也做线性决策规则问题会变成非凸的求解器基本没法处理。这是我早期踩过的坑。2.3 核心转化思路对偶把“最坏情况期望”拉下来对于模糊集内的最坏情况期望处理手法是把内层max-min问题通过对偶化成min问题然后和外侧的min决策合并这是分布鲁棒优化得以求解的关键一步。以Wasserstein球为例在原问题里我们要考虑模糊集里每个分布下的期望成本然后取最大值。直接枚举所有分布是不可能的但我们把决策变量固定后这个max问题本质上是线性规划对偶问题模糊集的约束对偶变量、目标函数成本函数都是线性的对偶之后可以变成一个有限维线性规划。再把对偶变量当成新决策变量加入原来的问题一起优化就得到了一个单一的最小化问题。整个过程完全自动化手工推导工作量不小但结果很漂亮原来看起来不可能求解的分布鲁棒问题最终就是一个规模略大的MILP。做这一步时我最深刻的体会是千万别手推对偶。用Yalmip的dualize或者手动写拉格朗日函数都行但最省力的是直接让Yalmip处理。给Yalmip传入模糊集的分布变量再使用implies和expectation等函数很多运算可以自动化完成。只是要注意对偶转化会引入大量辅助变量和约束在Matlab里建模时要注意稀疏性和索引逻辑否则矩阵组装阶段内存会爆掉。3. Matlab实现与代码拆解数据、架构和核心代码3.1 风电场景生成与模糊集半径计算分布鲁棒优化的输入是历史样本而不是假设分布。所以第一步是准备风电出力样本。如果手头没有真实历史数据我一般用两种方式生成仿真样本一是用ARIMA模型对风电功率时间序列建模生成预测误差样本二是用场景缩减技术比如同步回代消除fast forward selection把大量场景缩减到几十个代表性样本。样本准备好后要计算经验分布的均值、协方差矩模糊集用或者计算Wasserstein球半径。半径的理论取值随置信度和样本量的关系有明确公式但工程上我更习惯交叉验证用历史数据做回测尝试不同半径看哪个值能在备用约束达标的前提下总成本最低。这种做法比教条地套置信度公式要实用得多因为实际风电预测误差分布往往有偏且厚尾理论公式的前提条件未必完全成立。下面的代码片段展示了用MATLAB读取风电样本、计算Wasserstein球半径的示例这里展示一种基于KS距离的近似方式具体公式可以根据需要调整% 风电场出力场景样本wind_samples维度 [样本数, 时段数] % 这里展示Wasserstein球半径的经验式取法示意 N size(wind_samples, 1); % 样本总数 c 0.1; % 置信度参数工程上通过回测确定 r c / sqrt(N); % 半径随样本量缩减 % 经验分布均值与协方差用于矩模糊集构造 mu_hat mean(wind_samples, 1); Sigma_hat cov(wind_samples);这里的关键不是说半径公式必须长这样而是让你看到数据到模糊集的映射关系。实际项目里我还会把半径和置信度串起来比如搜索半径的候选集合每个半径用蒙特卡洛模拟评估实际切负荷风险最终选择一个风险满足要求但成本最低的值。3.2 Yalmip建模的总体架构整个代码我用Yalmip做建模框架求解器用Gurobi或者Cplex。Yalmip的好处是支持超值变量sdpvar和二元变量binvar混合建模并且对分布鲁棒优化中的对偶操作有良好的支持写出来的模型可读性远高于直接写矩阵组装。模型的主体结构分为两个模块。第一个模块是机组组合主问题定义机组启停变量、开机/停机变量、基准出力变量并建立启停逻辑约束最小开停机时间、启停费用线性化、出力上下限。第二个模块是分布鲁棒校正模块定义线性决策规则的系数矩阵把第二阶段优化变量的线性表达式代入功率平衡和爬坡约束再通过模糊集的对偶转化生成有限维约束。整体代码在Yalmip里的结构大概如下% 决策变量定义 u binvar(n_gen, n_period); % 机组开停状态 v_on binvar(n_gen, n_period); % 开机动作 v_off binvar(n_gen, n_period); % 停机动作 p_ref sdpvar(n_gen, n_period); % 基准出力参考 % 第二阶段线性决策系数 alpha sdpvar(n_gen, n_period); % 常数项 beta sdpvar(n_gen, n_period, n_sample); % 对每个样本的调整系数简化示意 % 约束系统功率平衡含风电样本期望项 Constraints []; for t 1:n_period Constraints [Constraints, sum(p_ref(:, t)) sum(wind_mean(:, t)) ... - sum(wind_curtail(:, t)) sum(load(:, t)) sum(load_curtail(:, t))]; end % ... 启动费用约束、最小开关机时间约束、爬坡约束等这只是主干示意完整代码里最重的部分是线性决策规则代进去之后的约束展开。以一个时段、一个场景为例第二阶段爬坡约束可以写成% 爬坡约束示例考虑调整量后的相邻时段出力差 for t 2:n_period for s 1:n_sample Constraints [Constraints, ... p_ref(:, t) beta(:, t, s) - (p_ref(:, t-1) beta(:, t-1, s)) ramp_up]; Constraints [Constraints, ... p_ref(:, t) beta(:, t, s) - (p_ref(:, t-1) - beta(:, t-1, s)) -ramp_down]; end end注意这段代码中我对上爬坡和下爬坡的不对称处理用了简化的写法实际需要根据爬坡约束定义区分增加和减少速率。这样一个场景一个场景展开再配合样本索引模型就会变成一个大的MILP。3.3 用Yalmip处理对偶把不可解的max-min转换成可解的MILP这一节是代码实现中最关键、也最容易卡住的地方。很多人拿到分布鲁棒模型后尝试自己手动推导对偶结果推导到一半就乱了索引一多就开始出低级错误。我的建议是能用Yalmip的高层操作就不要手动展平。Yalmip允许我们定义分布鲁棒约束的“期望”语义。比如对某个模糊集里的分布P第二阶段成本期望可以写成一个expectation变量。Yalmip配合rsolve函数鲁棒优化求解器时可以自动完成很多鲁棒对应转化。不过需要说明的是rsolve对某些特定结构支持得最好如果你的模糊集结构不属于预设模板还是要自己手动写对偶。手动对偶有个非常实用的技巧把第二阶段成本函数中所有决策变量替换为线性决策规则的表达式把对P的期望看作线性函数然后把“分布在最坏情况”写成对偶变量的线性约束。具体形式就是原问题对每个样本点计算成本给每个样本点引入对偶变量最坏情况期望等于一个有限维线性规划的目标值把这个线性规划的约束和对偶目标加进主问题。我在实现对偶时会把所有辅助变量的命名带上后缀比如gamma_dual、lambda_dual并在代码区用注释块标清楚每个变量对应原问题的哪个约束。这个习惯在调试时帮了大忙因为对偶问题变量个数多、约束少一旦约束松弛或者目标值离谱很快就能定位是哪一个对偶约束写错了。下面是模型核心区域示意% 对偶变量引入 gamma sdpvar(n_period, n_sample); % 对应样本点成本的对偶 lambda_bal sdpvar(n_period, 1); % 对应功率平衡约束 Constraints [Constraints, gamma 0]; for t 1:n_period for s 1:n_sample Constraints [Constraints, gamma(t, s) cost_linearized(t, s) lambda_bal(t) ... * (wind_samples(s, t) - mu_hat(t))]; end end % 目标函数中的最坏情况期望项用sum(gamma) * (r / N)近似示意写法 Objective startup_cost fuel_cost sum(sum(gamma)) * (r / N);上面这段是一个极大简化实际工程代码里远不止这么几行。但它反映了核心逻辑对偶变量把“最坏情况分布”的搜索压缩进了一组线性不等式里。当你把这个模型交给Gurobi时求解器在优化机组启停的同时也会自动挑选出对成本威胁最大的分布并优化应对策略。3.4 求解器选择与参数设置YalmipGurobi是我在实际测试中表现最稳的组合。Gurobi对大规模MILP的处理速度明显优于Cplex特别是当约束矩阵稀疏性好的时候。需要提醒的是分布鲁棒问题的对偶转化会引入大量连续变量和约束但整数变量数量基本不变所以本质上是“连续变量多、整数变量少”的MILP这类问题的求解难度比普通机组组合要低一些Gurobi的branch and cut效果非常好。求解器参数我一般这样设置MIPGap设为0.5%或者1%。机组组合的工程精度不需要太高多花几十分钟去追最后0.1%的gap完全不划算。TimeLimit设为3600秒超过就停下来取当前最优可行解。Threads设为物理核心数不要盲目的调大有时候线程数太多反而因同步开销降低速度。对偶变量连续问题的数值容差使用默认值即可但注意如果目标值出现微小负值如-1e-9不要慌这是数值误差不是模型错。还有一个经验先跑一个不包含线性决策规则、直接用期望场景的确定性版本确认约束逻辑没问题再增加模糊集和分布鲁棒项。不要一上来就上完整模型否则模型出了bug你根本不知道问题在机组组合约束还是对偶转化。4. 实验设计与结果对比成本、保守性和效率4.1 算例设置与评价指标我用改进的IEEE 6节点系统和IEEE 118节点系统分别做了测试。6节点系统里有3台常规机组和1个风电场负荷数据来自标准算例并放大到适合的比例。118节点系统相对复杂包含54台机组和多个风电场用于验证方法在大规模系统中的可行性。为了评估方法效果我用了四个指标调度总成本直接反映经济性。最坏情况成本即在模糊集内的最大期望成本反映鲁棒性。切负荷风险的蒙特卡洛估计用真实分布抽样5000次统计失负荷频率和失负荷量。求解时间评估实际工程可行性。这里最关键的评价逻辑是调度总成本低不难难的是在真实分布而不是模型假设的模糊集下仍然表现稳健。很多论文只汇报最坏情况成本但在实际运行里最坏情况出现的概率本身很低所以我还额外做了蒙特卡洛仿真用系统真实分布在实际生产中可能未知在测试中可以假设下检验决策质量。这个“闭环验证”意识建议每个人都养成。4.2 对比结果解读我在6节点系统上的典型结果是这样的方法调度成本万元蒙特卡洛切负荷概率求解时间秒确定性18.219.5%0.8随机优化500场景20.12.4%684传统鲁棒盒式Γ223.705.2分布鲁棒Wassersteinr0.221.30.7%97可以看到确定性方法虽然成本最低但切负荷概率高到不可接受随机优化把风险压下来了但计算耗时接近11分钟传统鲁棒做到了零切负荷但成本高出13.6个百分点分布鲁棒在保持成本只高出3%左右的情况下把切负荷概率控制在0.7%求解时间也才一分半。这个结果很好地说明了我前面强调的观点保守程度要可控不能一刀切地追求绝对安全。值得注意的是随着半径r的增加切负荷概率继续降低的同时成本也在上升。这个规律给实际调参提供了依据如果你所在地区的风电穿透率不高可以选大一点的半径追求安全如果风电占比高经济性压力大就用小半径配合实时备用市场补足风险。4.3 线性准则近似带来的最优性损失我也专门测试了线性决策规则相对“完全灵活调度”的最优性损失。方法是用一个理论下界做对比第二阶段不限制决策规则形式而是对每个场景单独优化得到一个随机优化的下界值。结果发现在6节点系统中线性决策规则导致的最优性损失大约在1.5%到2.8%之间具体取决于半径大小。在118节点系统中这个损失基本落在3%以内。这个数值完全在可接受范围内。原因也很直观机组组合的成本函数是分段线性的爬坡约束也是线性的系统动态对风电扰动的响应本来就接近线性线性决策规则正好契合系统的真实结构。只有在极端的非线性约束比如阀门点效应非常显著下线性决策规则才会带来较大误差。所以如果你的算例里包含了明显的非线性机组成本曲线我建议先把成本分段线性化再做线性决策规则这样误差能控制得更好。5. 常见问题与排查方向5.1 求解时间过长怎么处理分布鲁棒机组组合最常见的问题就是模型规模太大导致求解超时。我的处理顺序是先检查MILP的整数变量数量。如果整数变量数超过1万量级就要考虑用聚合方法把时段聚合成几个代表性区间或者把机组的组合建模为近似聚合机组。但这会牺牲精度所以我一般优先做第二件事减少第二阶段样本数。样本数和对偶变量数是线性关系样本从500个降到50个连续变量规模直接缩小十倍求解时间往往能降一个数量级。当然样本太少会导致经验分布不可靠。工程折衷方案是用场景聚类把500个样本缩减成30到50个代表性场景再为每个代表性场景加权。这本质上是用场景缩减代替简单抽样既保留分布形状信息又不至于让变量爆炸。还有个被我验证过多次的小技巧先用宽松的MIPGap比如5%跑出一个可行解把它作为初始解给下一步精细求解。Gurobi接受初始解后剪枝效率大幅提升整个流程反而比一步到位快得多。5.2 对偶约束写错导致的数值异常对偶转化的代码最常见的问题是双变量和原约束的对应关系弄混。典型症状是目标函数出现意外的负值或者最优解明显偏离物理预期。我的排查方法是做两类测试。第一类是退化测试把模糊集半径设成0此时分布鲁棒应退化为样本均值随机优化目标函数应该和随机优化的期望成本接近。如果差异很大那说明对偶部分有问题而不是机组组合约束的问题。第二类是边界测试把负载设为零、风电预测设为零让所有机组强制处于极低出力状态检查对偶变量和目标值是否合理。这类“边界测试”可以在全局优化前就发现模型错误成本极低。如果非要自己手推对偶我会建议先写成拉格朗日形式再用纸笔把变量一一对应列出来核对每个原约束的对偶变量是否只有一个。大多数错误都源于一个约束写了两个对偶变量或者漏掉了非负约束。5.3 Yalmip模型建模中的常见报错我整理了一份Yalmip建模过程中的高频问题速发表方便大家排查症状可能原因解决办法“Double tells” 错误对sdpvar变量做了双重赋值或约束冲突检查是否有约束被重复加进变量结构求解器报“Q not PSD”二次项出现在约束中分布鲁棒转化后应为线性检查是否有未展开的bilinear项求解器无限循环约束里混入了非线性函数用value()逐项检查表达式的非零项类型目标值为NaN数据里有NaN或者约束不可行去掉异常样本加可行域检查内存溢出约束矩阵过于稠密尽量用稀疏矩阵组织数据避免全0变量参与计算Yalmip有一条非常好用的命令yalmiptest专门测试求解器配置。有次我换了台新电脑装的Gurobi没有激活licenseYalmip报了一堆看不懂的错误最后就是通过yalmiptest定位到求解器没有正确注册。别小看这一步环境问题排除掉后面才能放心追模型问题。5.4 数据相关的工程坑最后提几个和数据打交道的坑。第一风电功率的单位转换很多公开数据是MW但机组容量标幺值用的是pu建模前必须统一。我见过不少人在约束里把风电出力直接填入MW结果和以pu表示的负荷对不上得到无意义的解。第二负荷曲线时段对齐风电样本的时段步长比如15分钟一个点必须和机组组合的调度时段比如1小时一个点严格一致不一致时要先做聚合或者插值绝不能直接在约束里混用。第三爬坡约束的时间尺度15分钟数据聚合到60分钟模型时爬坡速率要乘以聚合系数否则爬坡限制会过于严格或过于宽松。这些基础问题虽然说起来简单但在代码量大的时候特别容易忽视几乎每个项目都要被坑一两次。6. 一些个人体会和调试心得这条路线我从初版代码到稳定复现前后迭代了不少次。最有感触的一点是分布鲁棒优化的工程落地难点往往不在算法本身有多深而在模型转换和代码实现里的细节。线性准则、模糊集半径、对偶转化每个环节都像是多米诺骨牌一张倒下后面全乱所以一定要把每一步的验证做扎实。我给后来者的具体建议是先拿一个你能完全手算的3节点小算例测试你的对偶代码哪怕慢一点等小算例正确了再上大算例。这个习惯让我躲过了好几次在6节点和118节点系统上的返工。其次把模糊集半径、样本数、MIPGap这三个参数做成“旋钮”写进一个配置函数里这样实验对比时可以非常方便地批量跑而不是每次手改代码。最后别忘了把Wasserstein球半径的取值理由写进注释里因为几个月后你自己回来看代码也会忘记当初为什么选0.15而不是0.3。如果你在这套思路上继续深挖我觉得有三个方向值得扩展一是引入多风电场时空相关性的模糊集构造这会显著提升严重场景的辨识能力二是把线性决策规则扩展到分段线性决策规则能进一步压低最优性损失三是与储能和需求响应联合调度分布鲁棒框架天然适合处理多种不确定性同时存在的情况。希望这篇文章能帮你少踩几个坑快速把分布鲁棒机组组合跑起来。