简介面向风力发电出力随机性与间歇性带来的调度难题这份基于Matlab的参考代码以分布鲁棒优化为方法对机组组合问题进行建模与仿真适合电气工程、运筹优化方向的本硕学生及科研人员开展教研复现。相比随机优化与经典鲁棒优化分布鲁棒优化能在给定分布集合内提供更稳健的决策策略代码为此给出了可运行的实践起点。代码围绕分布鲁棒优化核心流程展开依次处理风电不确定集构建、模糊集参数设置、鲁棒对应模型求解以及机组启停与出力决策等关键环节并内含可直接运行的仿真结果便于结合论文图表逐段验证。压缩包共6个文件主体为Matlab源文件另附说明文档、PDF理论参考以及仿真效果示意图整体仅1.1MB结构清晰而轻量尤其适合在Matlab2014/2019a中快速部署。目前已有293人学习使用既可直接打开工程复现亦可作为理解风电不确定性与分布鲁棒优化的入门案例为后续扩展研究提供算法骨架与调参基础。1. 风力发电不确定性下的机组组合为什么用分布鲁棒优化风电在出力预测上留给调度员的余量往往比实际需要的多得多。把预测误差当成某个精确分布来建模一旦真实现场与假设不符备用容量和爬坡约束要么过于紧张、要么浪费严重。分布鲁棒优化处理的是“分布本身不确定”的情形用模糊集圈住所有可能的风电出力分布再在这些分布中寻找最坏期望成本。标题里这套“考虑风力发电不确定性的机组分布鲁棒优化matlab参考代码”就是把机组启停拼成一个混合整数规划把风速分布不确定性转成可计算的凸规划最后用求解器跑通的两阶段决策方案。对做电力系统调度的研究者和写论文的硕士研究生来说这套思路最值得借鉴的就是“模糊集半径怎么设、主问题和子问题怎么拆、场景数怎么验”三个动作。下文按建模、实现、分解、验证的顺序展开代码骨架可以直接套到自己的IEEE节点系统上。2. 不确定性建模分布鲁棒优化的核心前提与模糊集选择2.1 风电场出力不确定性数据的“形状”风功率预测误差不是正态分布实际中偏度和厚尾都很明显。风速服从威布尔分布但功率曲线把风速映射成功率之后误差分布会在低出力区和高出力区表现出完全不同的方差。传统做法是用场景scenario来描述这种不确定性从历史数据中采样一批可能的出力序列每个序列代表一种未来状态。场景数量的选择直接影响计算量1000个场景跑机组组合很容易把求解器拖到超时。分布鲁棒优化里的模糊集就建立在样本之上。样本的质量决定了模糊集的有效半径样本之间的“距离度量”决定了最坏情况长什么样。我在实际项目中一般先把风电历史出力归一化到[0,1]区间然后统计均值向量μ和协方差矩阵Σ。这两者是构造矩模糊集的最小输入也是后续Wasserstein模糊集标定的对照基准。归一化这一步容易被忽略但模糊集半径的数值大小完全取决于出力是否做了标幺化处理量纲不统一时半径本身没有任何物理含义。对场景数据的预处理建议保留三类指标各时段出力的均值、各时段出力的方差、相邻时段出力的协方差。这三类指标在模糊集对偶变换后都会出现在最坏期望的表达式中提前算好可以避免在每一次迭代中重复统计。2.2 从随机规划到分布鲁棒三个模型的递进随机规划Stochastic Programming假设风功率分布P是已知的目标函数写成期望形式 $ \min \mathbb{E}P[f(x,\xi)] $。实际操作中这个分布是用历史场景拟合出来的一旦真实分布与假设不符优化结果就会在边缘场景上失效。鲁棒优化Robust Optimization走向另一个极端假设出力只落在集合Ξ内目标写成 $ \min \max{\xi \in \Xi} f(x,\xi) $完全不使用概率信息结果通常过度保守。分布鲁棒优化Distributionally Robust Optimization, DRO把这两个极端做了折中。它假设真实分布P属于一个模糊集F目标函数写成$$ \min_{x \in X} \max_{P \in \mathcal{F}} \mathbb{E}_P[f(x,\xi)] $$三个模型的差别在“你知道什么”这个问题上随机规划知道分布的形状鲁棒优化只知道上下界分布鲁棒知道分布落在某个邻域内。电力系统场景里风电预测误差的均值、协方差相对容易获得但分布形状难以精确刻画这正好落在DRO的适用区间内。表三种不确定性建模方式的对比模型输入数据保守程度计算复杂度随机规划精确分布低低单层期望鲁棒优化不确定集合边界高中min-max两层分布鲁棒优化模糊集半径样本可调高两阶段三层模糊集半径是DRO中控制保守程度的关键参数。半径设为零DRO退化为随机规划半径趋于无穷大DRO趋近于鲁棒优化。这个性质的工程意义在最后一章会详细展开。2.3 矩模糊集与Wasserstein模糊集的建模纵深矩模糊集通过约束分布的矩信息来定义模糊集。最常用的一阶矩模糊集和二阶矩模糊集可以写成$$ \mathcal{F}_{mom} \left{ P : (\mathbb{E}_P[\xi] - \mu_0)^T \Sigma_0^{-1} (\mathbb{E}_P[\xi] - \mu_0) \leq \gamma_1, ; \mathbb{E}_P[(\xi - \mu_0)(\xi - \mu_0)^T] \preceq \gamma_2 \Sigma_0 \right} $$其中μ0和Σ0是名义均值与名义协方差γ1控制均值偏离程度γ2控制协方差放缩系数。实际计算中二阶矩约束会引入一个半定矩阵变量W所以矩模糊集DRO的最终形式是一个SDP问题。对小规模算例YALMIP加SDPT3可以处理规模稍大时半定约束会成为性能瓶颈。Wasserstein模糊集是另一个方向。它直接使用经验分布 $ \hat{P}N \frac{1}{N} \sum{n1}^{N} \delta_{\hat{\xi}_n} $用Wasserstein距离圈出一个球$$ \mathcal{F}_{W} \left{ P : W(P, \hat{P}_N) \leq \varepsilon \right} $$工程上Wasserstein模糊集的好处是不需要做矩估计直接使用历史样本代价是半径ε的标定与样本维度相关维度越高同样半径所覆盖的分布越多。凸性方面Wasserstein球内的最坏期望有一致的对偶表达形式缩放性好这也是近年电力系统文献中Wasserstein词频远高于矩模糊集的原因之一。2.4 最坏期望的对偶变换把min-max变成可计算的凸优化两阶段DRO中第二阶段问题写成$$ Q(x, \xi) \min_{y \in \mathcal{Y}(x, \xi)} ; c^T y $$第一阶段决策x影响可行域第二阶段决策y吸收不确定性的冲击。最坏期望 $ \max_{P \in \mathcal{F}} \mathbb{E}_P[Q(x,\xi)] $ 不能直接计算因为Q是ξ的隐函数。Wasserstein模糊集下利用对偶原理可以得到$$ \max_{P \in \mathcal{F}W} \mathbb{E}P[Q(x,\xi)] \inf{\lambda \geq 0} \left{ \lambda \varepsilon \frac{1}{N} \sum{n1}^{N} \max_{\xi \in \Xi} \left[ Q(x,\xi) - \lambda | \xi - \hat{\xi}_n | \right] \right} $$这个变换的关键在于对偶后外层是尺度为1的变量λ内层是对每个样本点做一个以样本为中心的“惩罚项”。λ的物理意义是模糊集半径的边际价格——λ越大模糊集半径对目标函数的影响越敏感。矩模糊集的对偶形式不同最坏期望会转成$$ \mu^T v \langle \Sigma, W \rangle \text{常数项} $$其中μ和Σ是名义矩v和W是拉格朗日乘子。两者代码实现上的差距很明显Wasserstein模糊集只需要一个标量λ和N个辅助变量矩模糊集需要一个向量和一个矩阵乘子。对机组组合这种本身已经包含大量二值变量的模型来说Wasserstein模糊集在求解器友好度上更优。提示对偶变换成立的前提是第二阶段问题为线性规划或可以等价改写成LP。机组组合的二阶段子问题如果不包含整数变量强对偶成立变换无误如果子问题里也出现了0/1变量则要先做整数变量的枚举或决策变量仿射重构否则对偶间隙会破坏最优性。3. 机组组合模型目标、约束与分布鲁棒等价变换3.1 系统结构与机组组合的混合整数线性模型机组组合Unit Commitment是要确定每台火电机组在每个调度时段内的启停状态与出力水平。系统结构上一个典型算例包含多台火电机组、一个风电场和一个负荷。调度周期通常取24小时时段间隔为1小时状态变量为0/1启停变量连续变量为机组出力。目标函数包含两部分确定性部分是第一阶段的启停成本与空载成本不确定性相关部分是第二阶段的可变燃料成本。全部写在一起的完整目标函数为$$ \min ; \sum_{t1}^{T} \sum_{i1}^{I} \left[ SU_{i,t} SD_{i,t} C_i^0 u_{i,t} \max_{P \in \mathcal{F}} \mathbb{E}P \left( \sum{t1}^{T} \sum_{i1}^{I} a_i p_{i,t} b_i u_{i,t} \right) \right] $$其中SU和SD是启停成本C_i^0是空载成本a_i和b_i是燃料成本线性化系数。实际项目中燃料成本常用分段线性函数逼近段数设置在3到5段之间超过5段时MIP求解时间会指数增长。3.2 功率平衡、爬坡约束与最小启停时间约束功率平衡约束是整个模型的骨架$$ \sum_{i1}^{I} p_{i,t} w_t D_t, \quad \forall t $$其中w_t是风电场实际出力D_t是负荷。这条约束必须严格满足否则调度方案在物理上不可行。爬坡约束描述机组出力的时间耦合$$ p_{i,t} - p_{i,t-1} \leq RU_i \cdot u_{i,t-1} P^{max}i \cdot (1 - u{i,t-1}) $$$$ p_{i,t-1} - p_{i,t} \leq RD_i \cdot u_{i,t} P^{max}i \cdot (1 - u{i,t}) $$第一式表示机组在上一时段停机时本时段启动后可以直接跳到最大出力第二式对称。RU_i和RD_i是上下爬坡速率单位MW/h。很多初稿直接写成 $p_{i,t} - p_{i,t-1} \leq RU_i$这在连续运行时段没问题但在启停边界上会让模型过紧导致无解。最小启停时间约束用线性不等式表达$$ \sum_{kt}^{tUT_i-1} u_{i,k} \geq UT_i \cdot (u_{i,t} - u_{i,t-1}) $$$$ \sum_{kt}^{tDT_i-1} (1 - u_{i,k}) \geq DT_i \cdot (u_{i,t-1} - u_{i,t}) $$这两条约束的物理含义不必多解释值得留意的是它们的展开方式如果把每个时段都写成独立不等式约束总数会超过I×T×max(UT,DT)规模求解器预处理阶段会自动消去冗余但不建议手工省略。3.3 不确定变量的纳入仿射可调策略与最坏期望第二阶段中风功率w_t与预测值 $ \hat{w}_t $ 之间的偏差用随机变量ξ_t表示。分布式鲁棒优化的核心动作是把w_t的表达写成一个弹性形式$$ w_t \hat{w}_t - \alpha_t (\xi_t - \mu_t) $$其中α_t是调整系数表示风电功率被削减或由火电机组补偿的比例。α_t0时风功率完全刚性输出α_t1时风功率全额参与调节。将上式代入功率平衡约束后火电机组出力被分成基准点与调节量两部分。这个仿射可调策略是在计算复杂度和表达力之间的实际妥协。完全不做仿射近似第二阶段每个场景都要重新求解一个线性规划数据规模无法接受做了仿射近似子问题变成一个纯LP可以安全地做对偶变换。第二阶段子问题的LP形式如下$$ Q(x, \xi) \min_{\Delta p} ; \sum_{i,t} a_i \Delta p_{i,t} $$$$ \text{s.t.} \quad \sum_i \Delta p_{i,t} \alpha_t (\xi_t - \mu_t), \quad \forall t $$$$ P^{min}i u{i,t} \leq p_{i,t} \Delta p_{i,t} \leq P^{max}i u{i,t} $$$$ -RD_i \leq \Delta p_{i,t} - \Delta p_{i,t-1} \leq RU_i $$这里Δp是火电机组的调节量是ξ的线性函数。LP的强对偶性保证了与Wasserstein模糊集对偶后没有对偶间隙。3.4 对偶后的确定性MIP形式与参数表将2.4节的对偶变换代入3.3节的子问题目标函数从min-max双重结构改写为单层确定性MIP$$ \min ; \sum_{t,i} C^{start}{i,t} \lambda \varepsilon \frac{1}{N} \sum{n1}^{N} s_n $$其中s_n是每个场景对应的最坏情况惩罚变量服从约束$$ s_n \geq Q(x, \hat{\xi}_n) - \lambda | \hat{\xi}_n - \xi^*_n |_p, \quad \forall n $$这里 $ | \cdot |_p $ 是范数通常取1-范数或无穷范数。取1-范数时约束可以线性化为一组绝对值不等式取无穷范数时需要引入对角矩阵变量。1-范数在电力系统场景中经验表现最好推荐优先采用。表一个3机组小型系统的样本参数符号含义取值I火电机组数3T调度时段24P^max最大出力[200, 150, 100] MWRU/RD爬坡速率50/50, 40/40, 30/30 MW/hUT/DT最小启停时间4/4, 3/3, 2/2 hN历史场景数200εWasserstein半径0.05 至 0.5关于N和ε的比例关系经验是N越大ε标定值越小。样本量从50增加到500时ε的合理取值大约从0.8下降到0.05具体数值需要在留出数据上做网格搜索。4. 用Matlab实现YALMIP建模与可复现的参考代码结构4.1 场景生成与缩减从历史数据到模糊集输入风速分布按威布尔分布采样再经过功率曲线转换成功率场景。功率曲线的分段线性化可以用Matlab自带的interp1实现。场景缩减直接调用kmeans把500个原始场景聚类到50个代表性场景聚类中心作为 $ \hat{\xi}_n $每个类内样本数量作为权重。%% 生成风电出力场景并做kmeans缩减 % 输入w_history 为 N0 x T 的历史出力矩阵N0500, T24 % 输出xi_hat 为 N x T 的缩减场景N50 N0 500; N 50; T 24; w_history normalize(w_history, 2, range); % 按行归一化到[0,1] % 用kmeans按欧氏距离聚类到N个代表场景 [idx, C] kmeans(w_history, N, MaxIter, 500, Replicates, 5); xi_hat C; % 聚类中心50 x 24 % 样本权重类内样本数占总样本数比例 weights accumarray(idx, 1) / N0; % 50 x 1逻辑说明kmeans的Replicates参数设为5因为聚类算法对初始中心敏感多次重复取最优结果。归一化必须放在聚类之前否则量纲差异会主导欧氏距离。weights数组在后续对偶变换中作为样本分布权重使用4.2节你会看到它们出现在s_n的组合中。参数说明N从原始500缩减到50保留的目的是让第二阶段对偶子问题只含50个场景约束。缩减比例控制在10到50倍之间缩减过多会丢失极端场景让模糊集半径失去意义缩减过少则对偶后的问题规模仍然过大。4.2 YALMIP建模主问题与Wasserstein对偶子问题的写法%% 两阶段DRO机组组合的核心建模片段 % 前置变量xi_hat (N x T), weights (N x 1), epsilon % 系统常量load (T x 1), Pmax, Pmin, RU, RD (I x 1) % u binvar(I, T, full); % 启停状态 p sdpvar(I, T, full); % 基准出力 lambda sdpvar(1, 1); % Wasserstein对偶变量 lambda 0 s sdpvar(N, 1); % 每个场景的最坏惩罚 Constraints []; % ---- 第一阶段的确定性约束 ---- for t 1:T Constraints [Constraints, sum(p(:,t)) w_base(t) load(t)]; for i 1:I Constraints [Constraints, Pmin(i)*u(i,t) p(i,t) Pmax(i)*u(i,t)]; end end % 爬坡约束、最小启停时间约束按第3章公式逐条追加 % ---- 第二阶段对偶变量的约束 ---- % s_n Q(x, xi_n) - lambda * ||xi_n - xi_j||_1 线性化 for n 1:N % Q(x, xi_n) 是LP的有限值省略内部变量推导直接加到约束里 Constraints [Constraints, ... s(n) Q_value(xi_hat(n,:), p, u) - lambda * norm(xi_hat(n,:) - xi_hat(1,:), 1)]; end Constraints [Constraints, lambda 0]; % ---- 目标函数一阶段成本 lambda*epsilon 场景惩罚加权 ---- obj sum(sum(SU.*max(0,u(:,2:end)-u(:,1:end-1)))) ... sum(sum(a.*p b.*u)) ... lambda * epsilon sum(weights .* s); optimize(Constraints, obj, sdpsettings(solver,gurobi,verbose,1));逻辑说明Q_value是一个子函数内部用linprog或对偶条件计算固定x和ξ下的第二阶段最优值。YALMIP的约束拼接用中括号加逗号逐条追加每次往Constraints里加一行而不是一次性写全部约束——后者在约束数量上千时会拖慢建模速度。参数说明SU是启停成本矩阵启动成本只在从0变1的时段计取上式用max(0, ...)实现这一点。weights与s做点积得到的是样本分布下的期望惩罚这正是Wasserstein对偶中 $ \frac{1}{N} \sum s_n $ 的加权版本。sdpsettings中的solver字段设置成gurobi如果本机没有Gurobi授权可以用matlab优化工具箱自带的intlinprog兜底小算例完全能跑。提示YALMIP对约束的构建方式很敏感。如果model检测结果为infeasible优先检查功率平衡约束中w_base的维度是否和xi_hat转置后一致。维度不匹配在YALMIP里不会报错但会产生一个空约束非常隐蔽。4.3 求解器配置Gurobi与Matlab优化工具箱的关键参数表常见求解器参数及推荐值参数推荐值说明MIPGap0.01求解器认为最优的可行间隙TimeLimit600 秒超过时限返回当前最好可行解IntFeasTol1e-5整数可行性容忍度LazyConstraints1CCG中启用惰性约束OutputFlag0/1关闭/开启求解过程输出NumericFocus1数值困难时提高范围0-3MIPGap设成0.01对论文级别的计算够用再小会明显增加求解时间。TimeLimit设成600秒可以保证离线分析时最多等10分钟在线调度场景一般设到60秒以内。Gurobi的一处使用细节如果模型里用到了YALMIP的value函数回读解就把OutputFlag设为0避免反复刷屏影响日志检查。4.4 Matlab环境下常见的三类报错与排查路径第一类报错是检测不到求解器典型提示是No suitable solver。检查顺序是先跑yalmiptest看YALMIP能否识别已安装的求解器再确认Gurobi的安装路径是否在Matlab的setenv中YALMIP通过mex接口调用Gurobi时需要在启动Matlab前把gurobi的路径加入环境变量。第二类报错是求解过程一直无解且模型的问题分类显示为Infeasible or unbounded。这时候把目标函数改成零先只做可行性检查。如果可行再逐步恢复目标函数定位是哪个二次项或非线性项让模型失去凸性。第三类报错来自稀缺的inf表达式代码中如果出现1/norm(...)这类分母含变量的表达式YALMIP会抛出次凸性错误解法是把范数约束转成一组线性不等式再写。5. 大规模求解Benders分解与列约束生成CCG的Matlab实现5.1 为什么要分解MIP规模与求解器瓶颈3机组的算例YALMIP生成的约束大约几千条直接用MIP求解器可以在几秒内出结果。把机组数提升到IEEE 118节点系统约50台机组、调度周期96时段之后变量数会超过5万其中二值变量接近5000直接求解的MIP间隙在30秒内很难压到1%以下。分解算法把原问题拆成主问题Master Problem和子问题Subproblem。主问题只保留第一阶段决策变量和二值变量子问题在固定第一阶段决策后求解第二阶段的线性规划。子问题的对偶信息被转换成割cut回传给主问题循环迭代直到上下界收敛。这样做的好处是每次只求解一个不含或少含整数变量的线性规划求解器压力显著下降。5.2 Benders分解主-子问题分工与割的形式Benders分解中主问题是O.A.$$ \min ; c^T x \eta $$其中η是最坏期望成本的逼近值。每次迭代子问题求出的对偶乘子π会生成一个Benders割$$ \eta \geq \pi_k^T (h - T x) $$这个割是原问题的下半估计。Benders与CCG的关键区别在于Benders每次从子问题取对偶乘子割与主问题的决策变量是线性关系CCG则是把子问题发现的最坏场景直接作为新列加入主问题生成的是强割而非弱割。强割的收紧速度远快于弱割所以工程实践里CCG在同迭代次数下几乎总是优于Benders。5.3 列约束生成CCG的实现要点%% CCG 双层循环伪代码结构Matlab风格 % 初始化 LB -inf; UB inf; k 0; x_hat ones(I, T); % 第一阶段初值 y_k {}; % 存储加入主问题的场景索引 tol 1e-3; while (UB - LB) / UB tol k k 1; % 步骤1求解主问题包含前k-1轮添加的所有场景 [x_current, obj_mp] solve_master_problem(y_k); LB obj_mp; % 主问题给出下界 % 步骤2固定x_current求解子问题 [worst_cost, xi_worst] solve_subproblem(x_current); UB min(UB, c * x_current worst_cost); % 步骤3子问题无界时添加可行性割有界时添加最优性割 if isinf(worst_cost) add_feasibility_cut(); else y_k{end1} xi_worst; % 新场景加入主问题 add_optimality_cut(xi_worst); end % 步骤4显示迭代信息 fprintf(Iter%d LB%.4f UB%.4f gap%.4f\n, k, LB, UB, (UB-LB)/UB); end逻辑说明solve_problem返回的最坏场景xi_worst是子问题中让第二阶段成本最大的那个场景索引来自对偶问题的取值点。CCG把xi_worst作为新列补充到主问题的场景集合中主问题规模一次比一次大但每次新增的列携带了子问题最关键的局部信息收敛速度远快于Benders。参数说明tol取1e-3在实际运行中兼具速度和精度。CCG的LB单调递增但UB可能波动所以UB取历史值的最小值而LB直接取主问题目标。输出fprintf的每一行都保留方便结束后绘制收敛曲线。5.4 收敛性、容忍度与停机判据停机判据写成(UB-LB)/UB小于容差这里注意分母是UB而非LB防止LB为负时误判。子问题求解后的可行性割和最优性割需要分开存储YALMIP实现时用cell数组分别管理。收敛速度上3机组算例到5次迭代内收敛到1%间隙50机组算例到12次左右。如果迭代到20次仍未收敛基本可以判定模糊集半径ε设得过小导致最坏期望值在场景间剧烈跳变。表Benders与CCG在UC问题中的实现差异方面BendersCCG割的类型弱割线性割强割场景列子问题信息对偶乘子π最坏场景索引主问题规模增长线性约束增加场景变量增加典型收敛迭代数15-305-15求解器友好度中高6. 验证与调试模糊集半径标定、场景数交互与成本-鲁棒性权衡曲线用留出数据验证DRO模型是否真的比随机规划更稳健。把历史数据按时间分成训练集和验证集训练集用来构造Wasserstein模糊集验证集用于模拟真实分布。在这个验证集上计算两种策略的实际总成本DRO的该方法设计目的是在偏离训练分布时仍然能维持可接受的性能。先关注“模糊集半径—总成本”曲线。ε从0.05增长到0.5目标函数中的最坏期望成本单调上升。这条曲线的斜率在某个点发生明显转折转折点左侧是数据驱动区域右侧是过度保守区域。经验上ε在该点的值为最优值它恰好覆盖验证集上观测到的分布偏移又不包含太多物理上不可能存在的分布。其次是场景数N的交互影响。N从20增加到200时相同ε下的期望成本会下降因为更多样本让模糊集更紧致。但N超过200后对偶后场景数量通过变量数的增长开始侵蚀求解效率。一个有效的组合策略是固定ε之后把N从50逐步提升到200观察最优成本的变化幅度变化超过2%就说明N太小需要增加。最后一处调试技巧是绘制ε-成本灵敏度图后叠加一条“零时刻爬坡违约率”曲线。做法是在验证集上固定调度方案模拟各个时序下的爬坡约束是否被击穿。违约率低于2%是工程可接受阈值。参考代码里预留了违约率统计脚本的接口把求解结果导出为矩阵后对每个验证场景逐时段回放约束即可。验证集上的分布偏移如果过大先检查历史数据的时间跨度是否覆盖了极端气象事件。风电大发和枯风期混在一起的样本集会让模糊集半径偏向高风险状态。合理的做法是把不同季节的数据分开建模针对夏季和冬季分别标定ε再在跨季节的验证集上做整体评估。本文还有配套的精品资源点击获取