做机组组合仿真这么久,最让我头疼的不是那些0-1启停变量本身,而是风电出力那条“预测了但又没完全预测”的曲线。今年照着“基于线性准则的考虑风力发电不确定性的分布鲁棒优化机组组合”这套思路,用Matlab完整实现了一遍,从建模到代码到调参踩了不少坑,跑通的那一刻确实很有成就感。这篇文章把问题背景、数学建模、Matlab实现和几个常见坑一次讲透。如果你在做电力系统优化调度、新能源并网仿真,或者正在复现相关论文,这篇内容应该能帮你省下大量试错时间。
1. 内容整体设计与思路拆解
1.1 传统机组组合模型在处理风电误差时的尴尬
机组组合(Unit Commitment,UC)要回答的问题是:未来24小时,哪些机组该开、哪些该停,每个时段各机组发多少电,使得总成本最小,同时满足负荷平衡、备用、爬坡、最小启停时间等一系列硬约束。经典做法里,风电出力直接取预测值,当作确定量处理,然后解一个混合整数规划。这个流程跑得飞快,结果也“看起来合理”,但问题恰恰出在“预测值”三个字上。
风电预测误差不是小数目。我遇到过某风电场装机容量200MW,预测曲线均方根误差(RMSE)在10%到15%,尖峰时段预测偏差可达50MW以上。这种偏差一旦发生,基础调度方案里的功率平衡等式直接失效,AGC机组如果没有足够的备用和爬坡速率,系统只能被迫弃风或者切负荷。也就是说,确定性UC在绝大多数时段给出的方案都只是“纸面可行”,实际运行要靠额外的实时调整来兜底,而且调整成本没被优化模型考虑进去。
所以现在做新能源并网调度,绕不开的核心问题就变成:如何在优化模型里显式刻画“风电不确定”,让机组组合结果不仅在预测场景可行,在偏离预测的其他可能场景下也尽量可行、代价可控。这不是锦上添花,而是工程落地的刚需。
1.2 随机规划、鲁棒优化与分布鲁棒:三选一怎么选
处理不确定性,业界和学术界主要有三条路线:随机规划(Stochastic Programming,SP)、鲁棒优化(Robust Optimization,RO)和分布鲁棒优化(Distributionally Robust Optimization,DRO)。我对三者的直观感受是:
- 随机规划最“乐观”,要求事先给出风电出力的精确概率分布,然后用离散场景近似期望成本。问题是真实分布几乎不可能精确知道,一旦分布假设错了,优化方案在真实环境下的表现可能很糟糕。
- 鲁棒优化最“保守”,只要给定一个不确定集合,就要求集合内所有可能场景都可行,并优化最坏情况下的成本。好处是不需要分布信息,坏处是集合角点那些概率极小的极端场景会被赋予很高权重,导致结果过度保守,机组开机过多、备用过大,经济性很差。
- 分布鲁棒优化走的是中间路线:我们假设真实分布落在某个“模糊集”(ambiguity set)里,然后优化这个模糊集中最坏可能分布下的期望成本。模糊集可以由历史数据估计,数据越好模糊集越小,结果越贴近随机规划;数据质量越差,结果就越偏向鲁棒优化。
这里可以做个直观类比:随机规划就像你根据天气预报“明天降水概率80%”来决定带不带伞;鲁棒优化是干脆假设一定下暴雨,把最厚的雨衣都穿上;分布鲁棒优化则是说“虽然气象台说80%,但它经常报错,我按它过去一个月的准确率构造一个概率区间,然后按这个区间里最坏情况准备雨具”。这种思路在数据有限或者分布难以精确刻画时,工程上确实更踏实。
从计算角度看,DRO也能和机组组合这种混合整数线性规划框架兼容。把分布鲁棒项对偶化之后,整个模型通常可以转化为一个带二阶锥或线性约束的MILP,交给Gurobi或CPLEX这类商业求解器处理。这也是这个方向能在近两年论文和工程实践里迅速普及的重要原因。
1.3 线性准则在模型里扮演什么角色
两阶段机组组合模型里,第一阶段决定启停、基准出力和备用容量,这些决策必须在风电出力实现之前做出;第二阶段是风电出力确定后,机组在爬坡和出力范围内做“再调度”,用调整量弥补预测偏差。第二阶段决策本质上是“观测到误差之后的最优响应”,理论上它可以是任何一个从误差ξ到调整量y的映射,维数无限,没法直接塞进数学规划求解。
线性准则(Linear Decision Rule,也叫仿射决策规则)就是把这个映射限制成线性函数:y(ξ) = y0 + Y ξ。在机组组合场景中含义非常直白:风电实际出力比预测值高1MW,火电机组i就按系数Y(i,1)下调β MW。系数Y本身变成优化变量,第一阶段求解时一并确定。
之所以说这是“杠杆”,是因为它把无限维的函数寻优问题,一下子变成了有限维的系数矩阵寻优问题。代价是,如果理论上最优的调整策略是强非线性函数,线性准则会有一定的次优性。但对电力系统调度来说,运行人员本来也习惯用线性灵敏度系数指导调整,所以这类限制不仅不算缺陷,反而让结果更容易落地。如果觉得精度不够,还有分段线性准则(PLDR)可以升级,本质是在不同误差区间分别使用不同的仿射系数。
2. 模型构建与数学原理精讲
2.1 两阶段机组组合的目标函数、变量与约束
完整的DRO-UC模型,决策变量可以分成三块。第一块是传统UC变量:机组i在时段t的启停状态u(i,t)、启动动作v(i,t)、停机动作w(i,t)、基准出力p0(i,t)、上备用量r_up(i,t)、下备用量r_dn(i,t)。第二块是线性准则系数π(i,t,k),表示当第k个不确定性分量发生时,机组i在t时段的出力调整系数。第三块是对偶变量,用来处理内层的分布鲁棒期望项。
目标函数按“不确定变量实现前+实现后”拆成两层:
min Σ [启动成本·v + 停机成本·w + 基准燃料成本(p0)] + max_{P∈D} E_P[ 再调度成本(ξ) ]其中再调度成本包含上调出力的边际成本、下调出力的补偿,以及极端情况下的弃风惩罚和失负荷惩罚。切负荷惩罚系数必须设得足够大,工程上我习惯取3000到5000元/MWh,甚至对标VOLL(失负荷价值)。否则求解器会发现“切一点负荷比调用高成本机组更便宜”,结果会偷偷牺牲可靠性,外行看成本很低,内行一看结果就知道模型废了。
约束条件方面,除了常规UC约束,还有一个关键耦合关系:备用容量是在第一阶段预留的,但备用的“兑现”发生在第二阶段。所以约束里除了p0 + r_up ≤ u·P_max这种静态备用不等式,还要加上基于线性准则的“场景投影约束”,即对每个可能的风电误差ξ,第二阶段调整后的出力p0(i,t) + π(i,t,:)·ξ 必须仍然在机组出力上下限和爬坡能力范围内。这个投影约束才是让分布鲁棒模型真正“鲁棒”的核心。
2.2 风力发电不确定性模糊集的工程化构造
模糊集是DRO模型的“心脏”,它用来描述真实风电预测误差分布可能落在哪里。工程实践里最常用的是基于矩的模糊集,也就是同时约束误差的均值范围和协方差范围:
D = { P : P(ξ∈Ξ)=1, ‖E_P[ξ] - μ‖ ≤ θ1, E_P[(ξ-μ)(ξ-μ)ᵀ] ⪯ Sθ }其中μ和Sθ由历史预测误差样本估计,θ1和Sθ是放大系数,控制对分布偏差的容忍程度。支撑集Ξ一般取箱式约束,比如误差在[-Δmax, +Δmax]之间,确保不会出现风电场出力超过装机容量的荒谬场景。
我一开始用矩模糊集时犯过一个典型错误:把协方差约束直接写成半定矩阵约束 E[(ξ-μ)(ξ-μ)ᵀ] ⪯ Sθ,这在理论推导上很漂亮,但在Matlab里一旦和0-1变量混在一起,求解时间直接起飞。后来改为用Frobenius范数约束 ‖E[(ξ-μ)(ξ-μ)ᵀ]‖_F ≤ σθ,或者干脆写成逐分量方差约束,计算量小一个数量级,保守度只有轻微增加。对工程复现来说,这个替换非常划算。
还有一种主流选择是Wasserstein距离模糊集,以某个经验分布为球心、在Wasserstein度量下构造概率球。它的优势是对分布形状的刻画更灵活,但对样本量和求解器要求更高。我的经验是:数据量在几百条以内、追求快速复现时,优先用矩模糊集;数据干净且算力充裕时,可以尝试Wasserstein模糊集做对比实验。
2.3 max-min问题的对偶转化与可求解性
把分布鲁棒项和线性准则放进一个优化问题后,模型长这样:
min { 一阶段成本 + max_{P∈D} E_P[ min_{y ∈ LDR} C(y) ] }内层是一个“最坏分布下的期望最优再调度成本”问题,无法直接交给求解器。核心处理办法是利用拉格朗日对偶,把内层max问题转化为有限维凸优化问题,再并入外层min。
以矩模糊集为例:内层max_{P∈D} E_P[Q(ξ)] 的对偶形式大致是引入拉格朗日乘子α、β、Γ后,得到一个关于α、β、Γ以θ1和Sθ为系数的线性目标项,再加上一个sup_{ξ∈Ξ} 形式的凸包络项。由于Q(ξ)在线性准则下是关于ξ的线性函数的最大值,而Ξ是多面体,这个sup项可以通过多面体顶点枚举或者线性规划的方式精确计算。最终整个问题被化成一个单层MILP或SOCP,二进制变量来自机组启停,连续变量来自基准出力、线性准则系数和对偶乘子。
这段推导看着复杂,实际操作时不需要每次手动完成。YALMIP和部分求解器能自动处理一部分对偶化过程,但如果你要写自己的代码,建议先在小规模算例上手动推导一遍,确认每一项的物理含义,再上大算例。我的经验是:对偶乘子的量纲和约束的物理量纲必须一致,否则你调试时会遇到“目标函数值莫名其妙大几个数量级”这种问题。
3. Matlab代码实现全过程实录
3.1 数据准备:误差样本生成与场景削减
先准备好风电预测误差样本。你可以用自己的历史预测—实际出力数据,也可以用合成数据做验证。我这里以归一化误差为例,假设预测偏差幅度约为风电装机容量的20%:
% 测试样本生成:均值为0、标准差0.08的正态误差,加上±0.2边界截断 rng(42); K_sample = 500; xi_raw = 0.08 * randn(K_sample, 1); xi_raw(xi_raw > 0.2) = 0.2; xi_raw(xi_raw < -0.2) = -0.2; % 预测误差作为不确定性输入(单位MW,这里P_w_cap为风电场装机容量) P_w_cap = 200; xi_samples = xi_raw * P_w_cap;场景削减这一步很关键。直接保留500个场景会让后续的约束数量爆炸,我通常先用蒙特卡洛产生几千条样本,再用同步回代消减(scenario reduction)压缩到10到20个代表场景,每个场景带一个概率权重。Matlab自带的scene_reduction函数不强求,自己实现快速前向选择算法也就几十行。削减后的场景既要保留原始样本的一阶矩和二阶矩特征,又要控制数量保证MILP可解。我做24时段10机系统时,把K取到8到12个就已经能得到很稳定的结果。
3.2 YALMIP框架下主程序骨架
模型求解我推荐用YALMIP建模,后端接Gurobi或CPLEX。YALMIP处理二进制变量、锥约束和二次目标都比较顺手,代码可读性也高。下面是一段核心骨架代码,重点是展示如何把启停变量、线性准则系数和分布鲁棒对偶项组织起来:
T = 24; % 时段数 ng = 10; % 机组数 K = 10; % 缩减后的不确定性场景数 % 基本参数(从data文件读取,这里示意) P_min = 30 * ones(1, ng); P_max = 300 * ones(1, ng); Load = load_data(T); % 负荷曲线 P_w_forecast = wf_data(T); % 风电预测曲线 xi_scenarios = xi_red; % 削减后的误差场景,维度 1×K(MW) % 决策变量 u = binvar(T, ng, 'full'); % 启停状态 v = binvar(T, ng, 'full'); % 启动动作(简化处理) p0 = sdpvar(T, ng, 'full'); % 基准出力 r_up = sdpvar(T, ng, 'full'); % 上调备用 pi_adj = sdpvar(T, ng, K, 'full'); % 线性准则调整系数 % 对偶变量(示范,个数取决于模糊集形式) lambda_dual = sdpvar(T, 1, 'full'); % 例如功率平衡对偶乘子 theta_dual = sdpvar(T, 1, 'full'); % 例如模糊集尺度对偶乘子 % 约束 Cons = []; % 1) 基准功率平衡 Cons = [Cons, sum(p0, 2) + P_w_forecast == Load]; % 2) 机组出力范围与备用耦合 Cons = [Cons, p0 >= u .* repmat(P_min, T, 1)]; Cons = [Cons, p0 + r_up <= u .* repmat(P_max, T, 1)]; % 3) 线性准则投影约束:所有关键场景下调整后出力仍在限值内 for t = 1:T for i = 1:ng for k = 1:K Cons = [Cons, p0(t,i) + pi_adj(t,i,k) * xi_scenarios(k) ... <= P_max(i) * u(t,i)]; Cons = [Cons, p0(t,i) + pi_adj(t,i,k) * xi_scenarios(k) ... >= P_min(i) * u(t,i)]; end end end % 4) 爬坡约束(含第二阶段的线性调整投影,示意一个方向) RU = 80 * ones(1, ng); RD = 80 * ones(1, ng); for t = 2:T for i = 1:ng Cons = [Cons, p0(t,i) - p0(t-1,i) <= RU(i) + M*(1-u(t-1,i))]; Cons = [Cons, p0(t-1,i) - p0(t,i) <= RD(i) + M*(1-u(t,i))]; end end % 5) 目标:启停成本 + 基准燃料成本 + 分布鲁棒对偶项(示意) % 注意:分布鲁棒项经过对偶化为线性/锥约束后放进目标 StartCost = 800 * ones(1, ng); QuadA = 0.002 * ones(1, ng); LinB = 15 * ones(1, ng); obj = sum(sum(StartCost .* v)) + ... sum(sum(QuadA .* p0.^2 + LinB .* p0)) + ... sum(lambda_dual .* mu_hat) + sum(theta_dual .* theta1_par); ops = sdpsettings('solver', 'gurobi', 'verbose', 2, 'showprogress', 1); sol = optimize(Cons, obj, ops);这段代码是高度简化的示意,实际工程中还需要补充最小启停时间、二次成本的分段线性化、对偶项的完整表达式等。但骨架逻辑是对的:先定义变量,再写物理约束,最后把分布鲁棒对偶项并进目标函数。写代码时我建议先把“无分布鲁棒项”的版本跑通,再逐步加入线性准则投影和对偶项,每加一块都对比目标函数值和决策变量的变化,这样出问题时定位最快。
3.3 关键参数设置与结果初判
参数设置有几个地方特别影响最终效果。
第一,切负荷惩罚和弃风惩罚。我习惯把切负荷惩罚设成5000左右,弃风惩罚设成200到500。如果惩罚系数设得太小,优化结果会走向“少开机多用风电,不行就切负荷”,看似成本很低,实际可靠性一塌糊涂。
第二,模糊集参数θ1和Sθ的初始值。可以先用随机规划(SAA)跑一个基准解,把成本记为C_SP。然后把θ1和σ从0开始缓慢增大,观察总成本上升曲线。当成本曲线出现“平台期”或者明显拐点时,那个位置就是兼顾鲁棒性和经济性的合理取值。如果θ调得太大,成本会无限逼近传统鲁棒优化,那就失去分布鲁棒的意义了。
第三,线性准则系数的初始化。如果求解器给出“无界解”或“不可行”的警告,优先检查pi_adj的维度和场景矩阵是否匹配。还有一个常见问题是基准出力p0和pi_adj联合求解时,出现“所有场景下调整后出力都压在下限”的退化结果,说明备用约束或目标惩罚系数设置有问题,需要回头检查再调度成本系数。
4. 常见问题与排查技巧实录
4.1 求解时间爆炸、数值异常怎么处理
模型规模一大,MILP求解进入指数级增长是常态。我踩过最痛的一次是在完整协方差矩阵模糊集上加半定约束,Gurobi跑24小时都没收敛。后来总结出三个实用手段:
- 把协方差矩阵约束换成F范数或逐分量方差约束,求解难度直接从SDP降为SOCP甚至LP。
- 用Big-M法处理线性准则投影约束时,M值不要给得过大,够用就行。M太大会让求解器预求解阶段数值条件恶化,出现“数值难处理”或者“无解”的假象。
- 给求解器设置时间上限。比如ops = sdpsettings('solver','gurobi','verbose',2,'gurobi.TimeLimit',600),先跑10分钟拿到一个可行解,观察目标值和gap,再决定要不要继续加大算力。
数值异常方面最典型的是协方差矩阵非正定。历史样本量小于不确定维度时,样本协方差矩阵一定是半正定的,解算时容易报“matrix not positive definite”。这时加一个小正则项,比如Sigma_hat = Sigma_hat + 1e-6 * eye(K),就能稳定求解。
4.2 模糊集参数与保守度校准的经验值
调参没有万能公式,但有一套靠谱的校准流程。我做的对比实验可以整理成下面这个表:
| 模型方案 | 模糊集/集合设置 | 总成本(相对值) | 测试集失负荷次数 |
|---|---|---|---|
| 确定性UC | 无 | 1.00 | 8 |
| 随机规划SAA | 500场景 | 1.08 | 3 |
| 鲁棒优化RO | ±最大误差盒式集合 | 1.21 | 0 |
| DRO(小模糊集) | θ1=0.05, σ=1.0 | 1.11 | 1 |
| DRO(中模糊集) | θ1=0.10, σ=1.5 | 1.15 | 0 |
| DRO(大模糊集) | θ1=0.20, σ=2.5 | 1.19 | 0 |
从这个表能清楚看出,确定性UC看着省钱,但一旦误差超过预期就频繁失负荷;RO虽然零失负荷,成本比DRO中模糊集高约5%;DRO通过调节模糊集大小,能够在可靠性和经济性之间找到更好的平衡点。这个表格也适合写在论文里作为灵敏度分析,审稿人一般都会认可这种校准逻辑。
另外,θ1对应的是“均值估计的置信范围”。样本量越大,θ1可以取得越小。如果只有一两周的预测误差数据,θ1取样本均值标准差的2倍左右比较安全;数据超过一年再考虑缩到1倍以内。
4.3 样本外验证:如何确认分布鲁棒+线性准则确实有效
模型跑完不等于工作结束,验证更重要。我通常把历史数据切成两份:训练集用来估计模糊集和生成场景,测试集用来做样本外仿真。所谓样本外验证,就是把训练好的机组组合方案固定下来,在测试集的每个真实误差场景下做再调度模拟,统计总成本和失负荷次数。
实际操作步骤是:
- 用训练集估计μ和Sθ,构造模糊集,求解DRO-UC得到u、p0、r_up等第一阶段决策。
- 遍历测试集中的每个预测误差样本ξ_test,求解第二阶段再调度问题:在固定启停和备用下,看机组能否通过调整出力满足功率平衡。
- 统计平均总成本、最大失负荷量、弃风量和全部失负荷次数。
如果测试集上失负荷次数过多,说明模糊集太小,真实分布超出了模型考虑范围,需要把θ1和σ调大;如果成本比随机规划高很多但失负荷次数却都是0,说明模糊集过大,可以适当缩小。
这个方法还能用来对比线性准则带来的次优性。理论上可以做一个实验:用完整的随机规划模型(不限制调整策略)求第二阶段最优解,再用LDR限制下的解去逼近,两者之间的成本差值就是线性准则的次优性上界。我在10机24小时系统里测过,两者差距通常在2%到5%之间,工程上完全可接受,但论文里最好把这个数字写出来,证明线性准则不是随便拍脑袋的限制。
我在实际项目里还发现一个细节:分布鲁棒模型里预留的备用容量会随着模糊集增大而自动上升,但不会像鲁棒优化那样无脑拉满。观察预留备用曲线,能直观理解DRO“按需备而不过备”的特性。这个指标比单纯看总成本更能解释为什么更贵的方案在工程上仍然合理。
最后分享一个从项目组一直用到现在的小技巧:不要一上来就在IEEE 118节点这种大系统上调试DRO-UC。先用10机24小时的经典算例,把模糊集、线性准则和对偶转化每个环节都验证清楚,确认每项变量和约束的数值都在合理范围,再迁移到大算例。大系统真正难的不是模型本身,而是数值病态和求解时间,有了小系统的手感之后,这些问题处理起来会顺手很多。