1. 问题背景为什么风电光伏接入后传统优化调度不好用了先聊点实际的。做电力系统优化调度的朋友应该都有感触前几年做确定性优化还够用无非是负荷预测给一个固定值机组组合按这个值去算。但现在风、光渗透率一上来情况完全变了——风电场今天出力可能800兆瓦明天同一时刻可能只有80兆瓦光伏更夸张一片云飘过来五分钟内出力能掉一半。这时候你要还敢用单一场景去做决策调度员心里肯定发慌。不确定性带来的真正麻烦不只是“预测不准”而是“决策之后没有退路”。火电机组的启停决策今天做了明天凌晨才能变可风电出力明天中午到底是多少没人知道。这就引出了我这次要分享的核心内容计及风、光、负荷不确定性的两阶段鲁棒优化以及实现它必须掌握的两个关键工具——大M法和CCG算法。代码用Matlab实现完整思路和踩坑记录都在下面。这套方案解决什么问题简单说就是帮你在“最恶劣的出力场景”下依然能给出一个可行的、不过分保守的调度方案。它特别适合做微电网经济调度、主动配电网优化、机组组合、储能容量配置这些方向的研究生和工程师。你要是正在写这方面的论文或者导师让你把确定性模型改成鲁棒模型这篇文章应该能帮你少走不少弯路。2. 两阶段鲁棒优化的模型构建到底在优化什么2.1 两阶段模型的结构先决策后调整先说清楚“两阶段”是个什么意思。它不是时间上的两个时段而是决策逻辑上的两层结构第一阶段Here-and-Now决策在不确定性实现之前就要做出决策。典型的就是机组启停状态、储能充放电状态这类0-1变量一旦定了就不能改。第二阶段Wait-and-See决策等不确定性场景暴露出来之后在第一阶段决策的基础上做经济调度调整。典型的就是各机组的实际出力、弃风弃光量、切负荷量。第一阶段相当于“定盘子”第二阶段相当于“填细节”。两层决策串起来目标函数就是最小化第一阶段成本加上最恶劣场景下的第二阶段成本。写成标准形式就是min_{x} ( c^T x max_{u∈U} min_{y∈Ω(x,u)} d^T y )其中x是第一阶段的二进制/连续变量u是不确定参数风速、光照、负荷y是第二阶段调整变量U是不确定集合Ω(x,u)是给定x和u后的可行域。我第一次看这个式子的时候最困惑的地方是那个“max-min”嵌套结构。后来我找了个生活化的类比你开一家餐馆菜单x今天早上就得印好但今天到底来多少客人u你不知道。你需要确保的是就算今天来了最刁钻的一桌客人最恶劣场景你也能靠后厨调整y把菜做出来不至于砸招牌。鲁棒优化就是让你在“最坏情况”下依然有办法。2.2 不确定性集合怎么构建盒式区间为主流选择既然要处理不确定性第一步就是先把“不确定性”数学化也就是构建不确定集合U。常用的有三种盒式Box、椭球式Ellipsoidal、多面体式Polyhedral。实践里最常用的是盒式集合U { u : |u - u_hat| ≤ Δu, 且满足上下限约束 }也就是每个不确定参数都在预测值附近一个区间内波动。比如风速预测值是10m/s允许偏差2m/s那风速就在[8,12]之间取值。这里有个关键概念叫做“预算约束”Budget of Uncertainty公式是sum |u_i - u_hat_i| / Δu_i ≤ ΓΓ的取值在0到N之间N为不确定参数个数它控制着“最坏情况”的保守程度。Γ0时退化为确定性模型ΓN时就是所有参数同时取极端值模型最保守。我在实际项目里通常把Γ取在N/3到N/2之间。比如我有10个不确定参数Γ取4就差不多了。这样既不会因为太保守导致成本虚高太多又能保证一定的鲁棒性。这个取值没有绝对标准我一般会跑几组Γ的灵敏度分析画一条“成本-鲁棒性”曲线然后把拐点处的Γ作为推荐值论文里也好看。还有一点要提醒风电、光伏、负荷的不确定性是相互独立的建模时把它们拼到一个大向量里就行但它们的波动幅度Δu是不同的——风电波动大光伏次之负荷相对平稳。这个幅度会直接影响不确定集合的“体积”进而影响最终结果的保守程度。2.3 目标函数与约束条件具体问题里的数学表达拿一个典型的微电网系统来说假设里面有火电机组、风电场、光伏电站、储能系统和本地负荷。两阶段鲁棒模型可以写成第一阶段机组启停min sum( C_startup * v_i C_shutdown * w_i ) 第二阶段成本约束包括机组启停逻辑约束、最小启停时间约束等。第二阶段经济调度min sum( C_fuel * P_g ) C_curtail * P_curtail C_shed * P_shed约束包括功率平衡P_g P_w P_pv - P_curtail P_dis - P_ch P_load - P_shed机组出力上下限P_g_min ≤ P_g ≤ P_g_max储能约束SOC动态方程、充放电功率限值弃风弃光约束0 ≤ P_curtail ≤ P_w P_pv切负荷约束0 ≤ P_shed ≤ 0.1 * P_load给个5%-10%的允许范围就行第二阶段决策变量是连续的所以内层min是个线性规划这对后面用对偶变换或者CCG算法都是很重要的前提——如果内层混入整数变量处理难度会陡增。3. 求解方法拆解大M法怎么用CCG怎么迭代3.1 大M法把“或”逻辑变成线性约束大M法在鲁棒优化里几乎是绕不开的尤其是处理机组启停状态和约束激活条件的时候。它的核心思想是把一个逻辑条件如果某个状态成立则某条约束才生效转换成一个带大M的线性不等式。举个例子机组i如果处于停机状态u_i0那它的出力必须是0如果开机u_i1出力在[P_min, P_max]之间。这个逻辑可以写成P_i ≤ P_max * u_i P_i ≥ P_min * u_i当u_i0时约束强制P_i0当u_i1时约束退化为正常的出力上下限。这里的逻辑很直接但另一个更常见的场景是“条件约束”如果储能处于充电状态z1充电功率≤P_ch_max 如果储能处于放电状态z0放电功率≤P_dis_max。写成大M形式P_ch ≤ P_ch_max * z P_dis ≤ P_dis_max * (1 - z)还有支路潮流约束、网络安全约束里也经常出现类似的“条件激活”逻辑原理一样。大M怎么取值是这门手艺的关键。M取得太小会错误地砍掉可行域导致找不到解M取得太大数值稳定性会变差求解器很容易报“数值困难”或者给出错误的可行解。我的经验是M取目标变量物理上限的2-5倍就行不要动辄取到1e6这种吓人的数。比如机组最大出力是300MW那M取600到1000就足够了。如果你不确定可以先用一个较大的M跑一遍然后逐步缩小M观察最优解是否变化如果缩到某个区间内结果稳定了就取那个区间的中点。3.2 CCG算法核心列与约束生成主问题与子问题的攻防CCGColumn-and-Constraint Generation算法是解决两阶段鲁棒优化最主流的精确算法之一思路比Benders分解更简单粗暴也更高效。它把原问题拆成主问题MP和子问题SP然后反复迭代主问题一个松弛后的两阶段问题只考虑有限个“已发现”的恶劣场景。求出的最优解x*是第一阶段的候选决策。子问题固定x x*在不确定集合U里寻找让第二阶段成本最大的那个u*。这个u*就是当前最恶劣的场景。关键步骤来了找到u之后不是像Benders那样加一条割平面约束而是**直接把u作为一个新场景往主问题里加一组新的第二阶段变量和对应的约束**。这相当于“扩列”——每迭代一次主问题里就多出一个完整的场景。CCG的好处是迭代次数比Benders少得多实践中一般5-10轮就收敛了。子问题本身是一个两层优化结构如下max_{u∈U} min_{y∈Ω(x*,u)} d^T y内层min是线性规划可以用KKT条件或者对偶变换把它变成一个单层max问题。最常用的做法是对内层取对偶与外层max合并得到一个单层的max问题。如果内层min是LP而且满足强对偶条件那么这个操作是完全等价的。3.3 对偶变换的具体操作别在符号上翻车这是实现里最容易出错的地方我单独拎出来详细讲。假设内层问题写成标准形式min d^T y s.t. A y ≤ b B u C y e F u y ≥ 0对偶之后变成max -(bBu)^T λ (eFu)^T μ s.t. A^T λ C^T μ ≤ d 对应y≥0的对偶约束 λ ≥ 0 μ 自由变量实操提醒两个点第一等式约束的对偶变量是自由变量第二不等式方向要仔细核对差一个符号结果就完全不对。我一开始经常犯的错误是把λ的非负约束写成自由变量导致对偶目标函数在max时直接无穷大MATLAB里YALMIP报错或者求解器返回Inf。如果你不想手动推导对偶也可以用YALMIP的dualize函数或者直接让求解器内部处理但作为研究者我强烈建议你至少手推一次简单案例这样你对模型结构的理解会深很多出了问题也知道往哪查。3.4 整体求解流程CCG迭代步骤完整梳理完整的CCG求解流程我整理成下面这个步骤表照着走基本不会乱步骤操作内容输出结果1初始化设定UB∞LB-∞迭代次数k1定义初始不确定场景u^(1)通常取预测值初始场景2求解主问题含已知的k个场景得到x^(k)和第一阶段成本更新LBmax(LB, 当前主问题最优值)x^(k)3固定xx^(k)求解子问题得到最优目标值Q(x^(k))和对应最恶劣场景u^(k1)u^(k1)4更新UB min(UB, c^T x^(k) Q(x^(k)))UB5判断gap (UB - LB) / abs(LB) ≤ ε收敛则停止否则向主问题添加新变量和约束对应场景u^(k1)kk1返回步骤2收敛判断注意上下界的更新逻辑主问题因为场景在增加是逐渐“收紧”的所以下界LB是单调上升的子问题在最坏场景下求解目标值一般偏大所以上界UB是单调下降的。两者逐渐靠拢gap趋于0。收敛判据我一般取ε0.01即1%就够了。如果你追求精度可以取0.001但迭代次数会明显增加而且对于实际工程问题1%的偏差完全可接受。4. Matlab代码实现从建模到求解的完整拆解4.1 工具选型YALMIP求解器还是纯YALMIPMatlab里实现鲁棒优化的主流方案是YALMIP建模 商业求解器求解。YALMIP对鲁棒优化的支持相当好它内部可以直接处理一部分uncertain类型的变量但我个人的建议是不要依赖YALMIP内置的鲁棒优化模块而是自己手动实现CCG迭代。原因有两个第一手动实现能让你对算法流程有完全的控制模型稍微变一变比如加入整数变量、非线性约束也能应对第二YALMIP内置鲁棒优化模块的处理方式像个黑箱出了问题你很难排查。求解器方面线性规划子问题用linprog或者gurobi都行主问题里有二进制变量所以需要混合整数线性规划求解器首选gurobi没有的话cplex也可以再不行用Matlab自带的intlinprog对付着跑小规模算例。提示如果是学生或者没有商业求解器许可intlinprog在规模不大几百个变量的情况下完全够用但超过这个规模之后求解速度会明显下降。有条件还是建议装Gurobi学术许可免费速度能快一个数量级。4.2 核心代码框架主循环就是这么写的为了不写成占篇幅的流水账我核心代码只展示关键逻辑完整代码可以按这个框架自行补全。第一阶段主问题代码用YALMIP建模% 主问题定义 x binvar(n_g, 1); % 机组启停状态 P_g sdpvar(n_g, K); % 每个场景下的机组出力 % 目标函数第一阶段启停成本 第二阶段运行成本所有场景的加权这里不weighted是因为要取max obj c_start * x; for k 1:K % 当前场景u_k下的约束 Constraints [Constraints, P_g(:,k) P_min .* x]; Constraints [Constraints, P_g(:,k) P_max .* x]; Constraints [Constraints, sum(P_g(:,k)) P_w(:,k) P_pv(:,k) P_load(:,k) - P_shed(:,k)]; obj obj d_cost * P_g(:,k) shed_cost * P_shed(:,k); end optimize(Constraints, obj, sdpsettings(solver,gurobi));第二阶段子问题固定x后求解最恶劣场景% 子问题固定x求最恶劣u x_fixed value(x); % 从主问题拿到的解 % 内层变量 P_g2 sdpvar(n_g, 1); P_shed2 sdpvar(n_load, 1); % 内层LP对偶前形式 obj_inner d_cost * P_g2 shed_cost * P_shed2; Constraints_inner [P_g2 P_min .* x_fixed]; Constraints_inner [Constraints_inner, P_g2 P_max .* x_fixed]; % ... 其他约束 % 手动对偶或用YALMIP自动对偶 % 这里展示直接构造对偶问题的方式省去中间变量 lambda1 sdpvar(n_g, 1); % 对应下界约束的对偶变量 lambda2 sdpvar(n_g, 1); % 对应上界约束的对偶变量 obj_dual ... % 按对偶公式填写 Constraints_dual [lambda1 0, lambda2 0, ...]; % 外层maxu是不确定变量 u_w sdpvar(n_w, 1); u_pv sdpvar(n_pv, 1); u_load sdpvar(n_load, 1); % 不确定集合约束盒式预算约束 Constraints_unc [abs(u_w - u_w_hat) du_w]; Constraints_unc [Constraints_unc, abs(u_pv - u_pv_hat) du_pv]; Constraints_unc [Constraints_unc, abs(u_load - u_load_hat) du_load]; % 预算约束 Constraints_unc [Constraints_unc, sum(abs(u_w - u_w_hat)./du_w) ... sum(abs(u_pv - u_pv_hat)./du_pv) ... sum(abs(u_load - u_load_hat)./du_load) Gamma]; optimize([Constraints_dual, Constraints_unc], -obj_dual, sdpsettings(solver,gurobi)); u_worst [value(u_w); value(u_pv); value(u_load)];这里有个细节容易踩坑对偶问题和不确定集合要写在一个sdpvar优化问题里联合求解不能分开求解。因为外层max和内层对偶后是一个整体问题分开解就得不到最恶劣场景的正确值。4.3 主问题扩展场景的细节变量索引怎么设计才不出错CCG算法里最让人头疼的就是主问题中变量和约束的“动态扩展”。每轮迭代都要新增一个场景意味着你要新增一组P_g变量、一组P_shed变量、一组对应的平衡约束和上下限约束。我的建议是从一开始就把P_g定义成二维矩阵行是机组列是场景。这样每轮迭代只需要增加一列已经存在的列完全不受影响。% 迭代过程中的扩展操作 k k 1; % 新场景编号 P_g [P_g, sdpvar(n_g, 1)]; % 扩展一列 P_shed [P_shed, sdpvar(n_load, 1)]; % 新增约束 Constraints [Constraints, P_g(:,k) P_min .* x]; Constraints [Constraints, P_g(:,k) P_max .* x]; Constraints [Constraints, sum(P_g(:,k)) u_w_found u_pv_found u_load_found - sum(P_shed(:,k))];这里有个性能优化的点新增变量时用sdpvar动态扩展在YALMIP里是可以的但如果迭代次数多建议用cell数组或者预先分配一个足够大的结构避免反复动态扩展导致的内存碎片。我在一个案例里跑了50轮迭代初始P_g是20x1最后变成20x51动态扩展方式在时间上还是可接受的但如果机组数量超过100还是提前分配好空间更稳。4.4 数据处理风电、光伏、负荷场景怎么生成真实数据要么来自实际历史数据要么用预测误差分布采样生成。处理方式% 风速预测序列和误差带 wind_forecast [5.2; 6.1; 7.3; 8.0; 7.5]; % 单位m/s wind_delta 0.2 * wind_forecast; % 20%波动 % 风速转出力简化模型 P_w_hat 0.5 * 1.2 * 1.2 * wind_forecast.^3 / 1000; % 1.2kg/m3空气密度扫风面积1.2m2 % 实际用必须查风机功率曲线这里只是演示在子问题里不确定变量就是u_w它落在区间[P_w_hat - du_w, P_w_hat du_w]内。当风速超过切入风速或低于切出风速时出力为0这个逻辑也要体现在约束里或者在实际数据预处理时把边界值截断。负荷数据同理取历史98%置信区间的上下界作为波动范围。4.5 求解结果怎么解读不只是看一个最优值CCG跑完之后除了记录UB、LB的收敛曲线我还建议额外输出几个东西最恶劣场景的构成看最恶劣场景下哪些不确定参数偏离预测值最大这能告诉你系统最怕什么样的天气/负荷组合。各阶段成本的拆分第一阶段启停成本多少第二阶段燃料成本多少弃风弃光成本多少切负荷成本多少。如果切负荷成本很高说明系统裕度不够如果弃风弃光很多说明储能或线路容量配置有问题。与确定性优化的成本对比同样的数据跑一遍确定性调度Γ0对比成本增量。这个增量就是“鲁棒代价”——你为抵抗不确定性额外付出的成本。论文里这个对比几乎是必备的。画收敛图的时候别只画目标值把UB和LB两条线都画上看它们怎么从两头往中间“夹逼”。如果两条线波动很大或者不收敛大概率是主问题没加够场景或者子问题对偶写错了。5. 常见问题与排查技巧实录5.1 子问题不可行的排查思路这是最常遇到的问题。固定x之后子问题报Infeasible。第一步先检查x是否真的是可行解——比如机组启停组合中开着的机组总容量低于最恶劣场景下的最大负荷那子问题必然不可行。第二步检查第二阶段约束是否束得太紧。常见原因是功率平衡约束要求严格相等而系统总出力区间和最恶劣负荷区间没有交集。解决办法是把功率平衡写成带松弛变量的形式比如加上一个切负荷变量允许在极端场景下少量切除负荷这样模型从“硬鲁棒”变成“软鲁棒”可行性大大提升也更符合工程实际。第三步检查不确定集合是否过大。如果把Γ取到N所有参数同时取极端很可能在数学上就没有可行解。这时候要么缩小Γ要么扩大系统可调资源的范围比如增加储能容量、允许切负荷比例。5.2 收敛过慢或震荡的原因及调整CCG正常收敛是单调的上界单调降、下界单调升。如果出现震荡最常见的原因是子问题求解不精确导致返回的场景u不对。排查步骤检查子问题里的对偶问题是不是写对了尤其是变量符号和约束方向。检查求解器精度设置Gurobi默认的MIPGap如果是1e-4子问题返回的是近似最优解这会让UB的更新产生微小波动一般不要紧但如果震荡明显把精度调到1e-6。还有一种情况子问题有多个最优解多个不同场景对应同一个最优值每次迭代返回不同的u导致主问题增长的场景“不聚焦”。解决方法是在子问题目标里加一个很小的正则项比如0.0001 * norm(u - u_hat)^2优先选偏离预测值最远的场景这样能加速收敛。5.3 大M取值不当导致的数值问题YALMIP/Gurobi报“Numerical trouble”九成和大M有关。典型表现是同样的模型跑两遍结果不一样或者明明有解却告诉你不可行。我的排查习惯是把所有系数矩阵的尺度打印出来看看如果某个约束的系数在1e-8到1e8之间横跳那数值稳定性肯定完蛋。处理方法所有物理量统一单位制MW、MWh这些别一会在kW一会又在MW大M统一取值别一处用1000另一处用1e6如果还有问题把Gurobi的NumericFocus参数调到2或3。5.4 快速速查表我踩过的坑和对应的措施症状可能原因解决措施子问题返回Inf对偶问题方向错误/未约束变量符号重写对偶等式约束对偶变量设为自由主问题变量爆炸迭代次数过多场景全被加入检查子问题是否收敛过慢加正则项UB-LB收敛gap不降子问题场景更新不“恶劣”检查不确定集合约束是否生效求解器数值警告大M取值过大或量纲不统一统一单位制缩小大M调NumericFocus结果对Γ不敏感不确定集合没被正确传递到子问题打印u_worst确认确实偏离预测值Matlab内存溢出动态扩展变量太多预分配结构或分块处理场景6. 案例实测与结果分析6.1 算例设置与数据我这边用的测试系统是改装的IEEE 6节点微电网3台火电机组50MW/100MW/200MW一个50MW风电场一个30MW光伏电站峰值负荷280MW。负荷预测曲线取24小时典型日数据风速用Weibull分布拟合历史数据后取预测值光照用Beta分布拟合。不确定集合参数风电波动20%光伏波动25%负荷波动10%预算参数Γ取6一共20个不确定参数节点相当于限制同时偏极端的最多6个变量。6.2 CCG收敛行为分析运行结果第1轮迭代时UB-LB差距非常大gap接近35%因为主问题里只有一个预测场景下界很低而子问题找到了一个恶劣场景上界很高。第2轮加完场景后gap骤降到8%第4轮降到1.5%第6轮gap0.8%收敛。这个收敛速度验证了CCG在实际问题中是高效的。相比之下我之前用Benders分解做过同样的问题到第15轮才收敛到2%以内。所以CCG在鲁棒优化里确实有性能优势这也是它成为主流方法的原因。6.3 与确定性优化的成本对比与敏感性分析确定性方案Γ0的24小时总成本是18.6万元鲁棒方案Γ6是21.4万元成本上升了15%。多出来的2.8万元就是为抵抗最恶劣场景付出的“保险金”。再看不同Γ下的成本Γ总成本万元切负荷量占比弃风弃光率018.60%8.2%319.80%10.5%621.41.2%12.8%1023.12.5%15.3%2025.84.8%18.9%从数据看Γ从0到10成本线性增长还可以接受但Γ超过10之后成本加速上升而切负荷量也在增加——这说明系统在极端场景下确实“力不从心”了。工程上这个曲线很有用你可以拿给导师或者甲方看说明你选的Γ是有依据的不是随手拍脑袋。6.4 鲁棒方案的调度策略解读观察最恶劣场景下的调度结果有几个有意思的发现最恶劣场景下风电出力几乎全部处于波动下界光伏正午时段也打了七折而负荷恰好处于区间上界。三个“坏消息”叠加在一起。储能系统在鲁棒方案里充电策略更保守确定性方案里储能会在夜间低谷充满、白天峰值放光鲁棒方案里储能始终保留20%的SOC作为“应急储备”。机组组合层面鲁棒方案会多开一台小机组做备用而不是像确定性方案那样只开两台大机组。虽然煤耗高了一点但换来了负荷跟踪的灵活性。这些结果其实反映了鲁棒优化的本质它不是在寻找一个“最省钱的方案”而是在寻找一个“最让系统安全的省钱方案”。理解了这一点你就能明白为什么论文里总爱说鲁棒优化“牺牲一定的经济性换取更高的安全性”。7. 实操心得与后续扩展建议代码这一轮做完我最大的体会是两阶段鲁棒优化的难点不在求解器而在模型构建的前几步。你如果把第二阶段的LP写对、对偶变换不出错、CCG迭代框架搭清楚了后面的求解就是水到渠成的事。多花时间检查对偶式子的符号多打印中间结果验证每个环节比闷着头调求解器参数有效得多。如果你想把这套方案继续深化我建议按这三个方向扩展第一个方向把约束从盒式集合升级为多面体集合。盒式集合虽然简单但有时候过于保守。多面体集合通过引入预算约束能更精细地控制不确定参数之间的关联性让结果更经济。第二个方向引入分布鲁棒优化。两阶段鲁棒优化的升级版它假设不确定参数服从某种分布但分布本身有不确定性通过矩约束或者Wasserstein距离来定义模糊集。这个方法这几年很火发表论文的价值更高。第三个方向把CCG算法加速。当系统规模大到几千个节点时CCG的收敛速度会变慢。你可以试试在子问题上引入Benders割的组合策略或者在主问题里加入对偶信息做“加速割”这些改进思路都是可以写进论文的亮点。最后分享一个实际调试中的小技巧先跑一个超小规模的算例比如2台机组、3个节点手算出每一步的期望结果再去验证代码。别一上来就甩给求解器。这样做虽然前期慢一点但能让你精准定位问题出在建模、对偶、还是迭代逻辑上后面大算例才会顺。我就是靠这个笨办法把很多看起来“玄学”的bug变成了“一分钟定位”的明确错误。希望这篇内容能让你少走点弯路。做优化的人都知道算法本身不复杂复杂的是把现实问题翻译成数学模型时那些“说不清道不明”的取舍。坚持做下去你会越做越有感觉。