做微电网调度的复现我印象最深的一课是在MATLAB里把确定性经济调度跑通之后信心满满地把风光预测误差从±5%改到±30%结果同一套计划直接让储能越限、功率平衡失配。从那之后我就明白只要系统里装了可再生能源和储能调度模型里不考虑不确定性的最优大概率只是某个理想场景下的最优。这篇内容围绕的就是一套含可再生能源和储能的区域微电网多阶段鲁棒调度模型基于MATLAB完整复现核心思路是建立日前计划日内调整的两阶段鲁棒框架用区间集合刻画风光和负荷的预测偏差再通过列和约束生成算法CCG迭代求解最坏场景下的最优运行方案。文章会从模型为什么这样设计、不确定性集合怎么构建、CCG算法如何落地到MATLABYALMIPGurobi的代码骨架和复现时容易踩的坑完整过一遍适合正在做微电网优化调度、储能配置或鲁棒优化应用研究的朋友参考。1. 可再生渗透率攀升后传统确定性调度为什么越来越不好用1.1 确定性调度的假设与失效场景传统的微电网经济调度模型无论是目标函数还是约束条件都默认风电、光伏和负荷在某个时段的出力是一个已知的确定值。模型看起来非常干净目标函数是最小化运行成本包括向上级电网购电费用、微型燃气轮机的燃料费用、储能损耗成本约束包括功率平衡、储能SOC上下限、充放电功率限制、机组出力限制、联络线功率限制。这种模型在可再生能源渗透率比较低、系统备用充足的时候还能凑合用。但一旦风机和光伏在微电网里的装机占比上来了问题就暴露了。风电预测误差通常在15%到25%之间波动光伏在天气骤变时误差可以超过30%负荷预测虽然相对准但在极端天气下也会出现明显偏离。确定性调度把所有预测值当成真实值来处理等于把全部风险都压在了预测完全准确这个脆弱的假设上。我做过一个对比测试在一个24节点微电网测试系统里用确定性模型求出一组日前计划然后把实际风光出力设定为预测值的70%和130%两种偏差场景重新跑功率平衡校验结果有多组约束直接违反储能SOC甚至越过了下限。这就是典型的计划赶不上变化。1.2 随机优化、机会约束与鲁棒优化三条技术路线的实际取舍面对不确定性学术界和工程界常用的解决方案有三类很多初学者一上来容易混在一起方法基本思路优点主要问题随机规划用大量场景近似不确定性对场景求期望最优结果贴合概率分布经济性较好需要精确的概率分布场景数量大时计算量爆炸分布估计不准时结果不可靠机会约束规划允许约束以较高概率满足如95%概率满足功率平衡可以在安全性和经济性之间折中概率参数如何标定比较主观分布信息不足时难以求解鲁棒优化假设不确定性落在某个集合内保证最坏情况下约束也满足不需要精确概率分布结果对极端场景稳健计算效率相对可控结果偏保守成本一般高于随机规划在微电网这种规模不大、但安全约束要求高的场景下鲁棒优化是比较务实的选择。原因也很直接你通常拿不到足够多的历史数据来精确拟合风电和光伏误差的分布函数就算拟合出来分布参数本身也有估计误差。与其在概率上赌一个精度不如直接说我知道误差不会超过某个范围在这个范围内怎么折腾我的调度方案都兜得住。1.3 标题里鲁棒性和不确定性到底指什么标题里的鲁棒性英文对应Robustness放在调度模型里指的就是当实际的风电、光伏出力偏离预测值时系统依然能通过储能充放电调整、机组出力再调度、与主网功率交换等手段维持功率平衡并满足所有安全约束。衡量鲁棒性的标准不是平均工况下的成本而是最坏工况下的约束满足能力。不确定性在这个问题里主要来自三个方向源端不确定性风机出力和光伏出力受气象条件影响预测模型再准也有偏差荷端不确定性区域负荷受生产生活规律、天气温度等因素影响同样存在预测误差设备状态不确定性比如某台机组临时故障、某条线路检修但在调度主问题里通常把这类视为N-1预想故障单独处理不会混进源荷不确定性里。在工程复现时通常把源荷不确定性统一建模为一个不确定集合U集合里的每一个元素代表一组可能的源荷出力偏差向量。鲁棒优化要求对U里的任意一个偏差向量系统都能通过第二阶段调整找到可行解。这就是后面整个模型的核心出发点。2. 多阶段鲁棒调度模型的数学框架与不确定性集合设计2.1 两阶段模型的前后配合日前决策与日内调整标题里写的是多阶段鲁棒调度模型实际复现时最常见的落地形式是两阶段鲁棒优化模型也就是经典的三层结构min-max-min。之所以叫多阶段是因为调度决策本身是分层分时段的第一阶段日前决策阶段在获知次日风光和负荷预测值之后提前决定机组的启停状态、储能充放电基准计划、向上级电网购电的日前计划。这些决策的特点是必须在看到实际偏差之前做出所以也叫here-and-now决策。第二阶段日内调整阶段当实际风光出力、负荷与预测出现偏差后在日前计划的基础上通过调节储能出力、机组出力微调、可削减负荷等手段消除功率不平衡。这个阶段的决策是在偏差揭晓之后才能做的叫wait-and-see决策。两阶段之间通过功率平衡约束和设备运行约束耦合在一起。模型的目标函数也相应地拆成两块min (第一阶段运行成本) max_{u∈U} min_{y∈Ω(x,u)} (第二阶段调整成本)外层min决定第一阶段决策x内层max寻找使第二阶段成本最高的最坏偏差场景u中间层min则是在给定x和u之后寻找最优的调整策略y。整个模型求的是在所有可能出现在的偏差场景下系统总成本的最坏情况最优值。2.2 盒式、预算约束与椭球式不确定集合的适用场景不确定集合怎么建直接决定了模型的保守程度和求解难度。常见的有三种形式盒式集合Box直接把每个时段的不确定量限制在一个区间内如u(t)∈[u_min(t), u_max(t)]。这种集合最简单但要求所有时段同时达到最大偏差结果是最坏中的最坏过于保守。预算约束集合Budget在盒式集合基础上增加一个总偏差比例限制形如Σ|u(t)-u_hat(t)|/δ(t) ≤ Γ。这里的Γ就是预算参数表示在调度周期内最多允许多少个时段的偏差同时达到上限。Γ越小鲁棒约束越宽松成本越低Γ越大方案越保守。这个设置非常贴合实际因为真实的风光出力偏差不太可能24个小时每个小时都顶在最大值上。椭球式集合Ellipsoidal用2-范数约束替代1-范数额外限制形如‖W(u-u_hat)‖₂ ≤ Ω。椭球集合考虑了变量之间的相关性数学上更精细但会引入二次约束计算复杂度明显上升。在微电网调度里我最推荐的是盒式预算约束的组合。原因有三个物理意义清晰Γ可以理解为最坏情况下同时偏离预测的时间段数方便向工程人员解释模型保持线性结构经过对偶变换后仍然是LP或MILP可以用YALMIPGurobi稳定求解通过调节Γ的取值可以画出从激进到保守的完整成本曲线这是论文和报告里很有说服力的一张图。2.3 目标函数与约束的完整数学表达以含风电、光伏、储能和微型燃气轮机的区域微电网为例完整的鲁棒调度模型可以写成决策变量第一阶段燃气轮机开机状态z(t)、储能充电功率P_ess_c(t)、储能放电功率P_ess_d(t)、日前购电功率P_grid(t)、储能SOC(t)。第二阶段针对偏差场景u的储能调整量ΔP_ess(t)、机组调整量ΔP_g(t)、可削减负荷量ΔP_load(t)。目标函数min Σ_t [ C_grid * P_grid(t) C_g * P_g(t) C_ess * (P_ess_c(t)P_ess_d(t)) ] max_{u∈U} min_Δ Σ_t [ C_penalty * (ΔP_load(t) 弃风弃光惩罚项) ]约束条件功率平衡约束 P_grid(t) P_wind(t) P_pv(t) P_g(t) P_ess_d(t) P_load(t) P_ess_c(t) ΔP_load(t)其中P_wind(t)和P_pv(t)在第二阶段表示为P_wind_hat(t) u_w(t)、P_pv_hat(t) u_pv(t)。储能约束 SOC(t1) SOC(t) η_c * P_ess_c(t) * Δt - P_ess_d(t) * Δt / η_d SOC_min ≤ SOC(t) ≤ SOC_max 0 ≤ P_ess_c(t) ≤ P_ess_c_max * z_ess(t) 0 ≤ P_ess_d(t) ≤ P_ess_d_max * (1 - z_ess(t))这里需要注意的是充放电状态变量z_ess(t)的引入它保证储能不会同时充电和放电这在后面代码实现时是个关键细节。机组约束 z(t) * P_g_min ≤ P_g(t) ≤ z(t) * P_g_max 爬坡约束|P_g(t) - P_g(t-1)| ≤ Ramp_g与主网交互约束 P_grid_min ≤ P_grid(t) ≤ P_grid_max可削减负荷约束 0 ≤ ΔP_load(t) ≤ α * P_load(t)2.4 为什么多阶段不等于多周期很多初学者会把多阶段鲁棒理解成把24小时的调度周期分成很多段每段做一个鲁棒优化这是一个明显的误区。多周期优化只是把时间维度拉长决策变量在整个周期内一次性规划本质上仍是开环调度。多阶段鲁棒优化的关键不在于时间段的多少而在于决策的可调整性——第一阶段决策做出后可以看到不确定性的实现然后第二阶段再做出调整决策。这个特征对应的是闭环调度也就是日前计划日内实时调整的真实运行流程。用数学语言说真正的多阶段鲁棒模型必须具有min-max-min的嵌套结构而不是简单的min-max或者包含很多时段的一次性min。在复现时如果你的模型只是每个时段都预留备用容量然后用鲁棒约束包住所有时段这虽然也算鲁棒但不是多阶段鲁棒。要体现多阶段就必须显式区分可调整变量如储能的实时充放电调整量和不可调整变量如机组启停、日前购电计划并在目标函数中体现第二阶段对最坏场景的应对成本。3. 列和约束生成CCG算法迭代框架的收敛逻辑3.1 主问题的最小费用加最坏场景切割线两阶段鲁棒模型直接求解极其困难因为内层嵌套了max-min优化。工程上最常用的解法是列和约束生成算法核心思路是把原问题拆成一个主问题MP和一个子问题SP通过迭代逐步逼近最优解。主问题的形式如下min c * x η s.t. 第一阶段约束 η ≥ d * y_l, 对于已知的场景u_l x、y_l之间的耦合约束这里的关键是引入了辅助变量η它表示第二阶段成本的上界估计。每轮迭代中子问题会找到一个最坏场景u_l主问题就把这个场景对应的第二阶段变量y_l和约束作为切割加入模型让η不断收紧。主问题因为只考虑了已发现的最坏场景所以实际上是原问题的松弛版本它的最优值给出原问题的下界LB。3.2 子问题的max-min对偶转化子问题是在给定第一阶段解x_star之后寻找使第二阶段成本最大的场景u以及对应的最优调整成本形式是f(x_star) max_{u∈U} min_{y∈Ω(x_star, u)} d * y这个max-min结构不能直接丢给求解器。标准做法是对内层的min问题取强对偶。由于内层是关于y的线性规划只要可行域非空且有界强对偶成立就可以把min问题等价转化为其对偶max问题然后和外层的max合并变成一个单层max问题f(x_star) max_{u∈U, λ∈Λ} λ * (b - A * x_star - C * u)其中λ是对偶变量它的符号由原约束方向决定。如果原约束是≤型λ非负如果是等式约束λ自由变量。合并之后子问题变成带二次项的非线性问题因为λ和u相乘。在实际求解中这个双线性项通常可以用大M法或直接交给支持MIQP的非线性求解器处理Gurobi可以直接求解这类带双线性项的小规模问题。需要强调一点内层min问题必须满足强对偶条件。因此第二阶段调整问题必须建模为线性规划不能随意加入整数变量。如果第二阶段需要决策是否切除某条线路这类0-1变量标准CCG就不适用了需要改成嵌套CCG或者引入其他分解方法。3.3 收敛判据与最大迭代次数CCG算法的主循环流程如下初始化UBinf, LB-inf, k0设置初始场景u_0通常取预测均值 while k ≤ K_max 且 UB - LB tol: 1. 求解主问题得到x_k*和η_k*更新LB c*x_k* η_k* 2. 将x_k*代入子问题求解最坏场景u_k*和最坏场景成本f(x_k*) 3. 更新UB min(UB, c*x_k* f(x_k*)) 4. if UB - LB ≤ tol: 收敛退出 5. 否则将u_k*作为新场景加入主问题添加对应的第二阶段变量y_k和割约束 k k 1收敛判据一般取相对间隙或绝对间隙。我在复现时常用的组合是绝对误差取1e-4最大迭代次数取20轮。在24时段微电网模型里通常3到8轮就能收敛到很紧的间隙很少会真的跑到20轮。3.4 CCG与Benders分解的选择建议很多文献会同时提到Benders分解和CCG两者都是求解两阶段问题和分解主-子问题的经典算法但适用场景有明显差异对比维度Benders分解CCG切割信息通过对偶函数添加割平面直接添加最坏场景对应的完整约束迭代收敛对小规模问题收敛较快通常需要更少迭代次数尤其有连续变量时对整数变量的适应性一般需要处理整数对偶gap主问题中可以自然处理整数变量实现难度较低中等在微电网鲁棒调度里第一阶段通常包含机组启停这类整数变量而CCG的切割是把具体场景u_l下的约束完整塞进主问题对整数变量天然友好所以实际复现时我基本都用CCG。只有当子问题规模特别大、每次求解SP耗时很高时才会回头考虑Benders或者近似方法。4. MATLABYALMIPGurobi的关键代码实现4.1 环境准备与求解器适配复现前先确认环境MATLAB版本R2020b及以上即可建议R2022b以上对sdpvar和optimizer的兼容性更好YALMIP直接从GitHub拉取最新release不要用十年前的老版本求解器Gurobi 10.x或11.x学术许可免费申请Mosek 10.x也可以但Gurobi在MILP上的速度通常更优。在MATLAB里先跑一句验证yalmiptest能列出Gurobi的LP、QP、MILP、MIQP支持状态就说明环境没问题。如果用了旧版MATLAB配新版Gurobi经常会出现java相关报错那是MATLAB版本太老、Java版本不匹配导致的属于环境问题先升级MATLAB再排查。4.2 主问题的YALMIP建模代码下面是主问题的核心建模骨架采用24时段、1小时为步长的设定。nT 24; % 第一阶段变量 z_g binvar(1, nT); % 燃气轮机启停 P_g sdpvar(1, nT); % 燃气轮机出力 P_ess_c sdpvar(1, nT); % 储能充电功率 P_ess_d sdpvar(1, nT); % 储能放电功率 z_ess binvar(1, nT); % 储能充放电状态1充电0放电 SOC sdpvar(1, nT1); % 荷电状态 P_grid sdpvar(1, nT); % 购电功率 eta sdpvar(1); % 最坏场景成本上界 % 目标函数 objective_mp sum(grid_price .* P_grid) ... sum(g_fuel_cost * P_g) ... sum(ess_cost * (P_ess_c P_ess_d)) ... eta; % 约束 constraints_mp []; constraints_mp [constraints_mp, SOC(1) SOC_init]; for t 1:nT constraints_mp [constraints_mp, ... SOC(t1) SOC(t) eta_charge * P_ess_c(t) - P_ess_d(t) / eta_discharge, ... SOC_min SOC(t1) SOC_max, ... 0 P_ess_c(t) P_ess_c_max * z_ess(t), ... 0 P_ess_d(t) P_ess_d_max * (1 - z_ess(t)), ... z_g(t) * P_g_min P_g(t) z_g(t) * P_g_max, ... P_grid_min P_grid(t) P_grid_max, ... P_wind_forecast(t) P_pv_forecast(t) P_g(t) P_grid(t) P_ess_d(t) ... P_load_forecast(t) P_ess_c(t)]; end ops sdpsettings(solver, gurobi, verbose, 0, gurobi.NumericFocus, 3); optimize(constraints_mp, objective_mp, ops);这里有个细节给Gurobi设置NumericFocus为3可以明显降低大M约束带来的数值病态问题。如果模型里用到了大M这个设置几乎是必加的。4.3 子问题的对偶化与最坏场景求解代码第二阶段的最坏场景求解是核心代码分三步走第一步固定第一阶段变量第二步构造第二阶段LP并取对偶第三步求解单层max问题。% 固定第一阶段变量 x_star value([P_g; P_grid; P_ess_c; P_ess_d; SOC]); % 不确定变量 u_w sdpvar(1, nT); % 风电偏差 u_pv sdpvar(1, nT); % 光伏偏差 u_l sdpvar(1, nT); % 负荷偏差 % 第二阶段调整变量 delta_ess sdpvar(1, nT); % 储能功率调整 delta_g sdpvar(1, nT); % 机组出力调整 delta_load sdpvar(1, nT); % 可削减负荷 % 不确定集合盒式预算 U_con [-u_w_max u_w u_w_max, ... -u_pv_max u_pv u_pv_max, ... -u_l_max u_l u_l_max, ... sum(abs(u_w) ./ u_w_max) sum(abs(u_pv) ./ u_pv_max) Gamma]; % 第二阶段约束仅对应给定x_star SP_con [...]; SP_obj sum(penalty_load * delta_load penalty_curtail * (u_w u_pv) ... C_delta_g * abs(delta_g) C_delta_ess * abs(delta_ess)); % 通过YALMIP取对偶或直接用求解器求解最大最小问题 % 更通用的做法是手动构造对偶变量把max-min转为单层max dual_lambda sdpvar(size(A_sp, 1), 1); dual_mu sdpvar(size(A_eq, 1), 1); % 对偶目标函数注意线性化双线性项 dual_obj dual_lambda * (b_sp - A_sp * x_star_vec - B_u * [u_w; u_pv; u_l]) ... dual_mu * (b_eq - A_eq * x_star_vec); % 合并后的最坏场景问题 SP_merged_con [dual_lambda 0, dual_mu free, U_con]; result_sp optimize(SP_merged_con, -dual_obj, ops_sp);实际上当双线性项存在时Gurobi会以MIQP方式求解变量数不大时速度可以接受。如果模型规模较大建议用big-M法把双线性项线性化把问题转成MILP。4.4 迭代主循环与结果保存主循环代码框架如下UB inf; LB -inf; k 0; tol 1e-4; Kmax 20; cut_results {}; while k Kmax (UB - LB) tol k k 1; % 1. 求解主问题 optimize(constraints_mp, objective_mp, ops); LB value(objective_mp); x_cur value([P_g; P_grid; P_ess_c; P_ess_d; SOC]); % 2. 求解子问题 optimize(SP_merged_con, -dual_obj, ops_sp); sp_cost value(dual_obj); u_cur value([u_w; u_pv; u_l]); day_ahead_cost value(sum(grid_price .* P_grid) ... sum(g_fuel_cost * P_g) ... sum(ess_cost * (P_ess_c P_ess_d))); UB min(UB, day_ahead_cost sp_cost); fprintf(迭代%d: LB%.4f, UB%.4f, gap%.6f\n, k, LB, UB, UB-LB); % 3. 添加新场景切割到主问题 y_new sdpvar(size_y, 1); constraints_mp [constraints_mp, ... y_new对应的第二阶段约束把u_cur作为常数代入, ... eta 第二阶段成本表达式(y_new)]; cut_results{k}.u u_cur; cut_results{k}.cost sp_cost; end有几个实现上的细节值得说不要把主问题的sdpvar变量每轮重新创建。变量结构要保持不变只是追加约束。否则MATLAB每轮都要重新解析整个符号表达式速度直接掉一个数量级。主问题的optimize结果如果返回problemsolver 信息先打印optimize的第二个输出看看是primal infeasible还是numerical issues再决定调参数还是调模型。保存每轮的最坏场景u_cur最后画图时可以直接把最坏场景下的风电光伏曲线和基准预测曲线对比很直观。4.5 代码性能优化的几个细节复现完跑第一版10轮迭代可能要跑几分钟如果觉得慢按下面这几个方向优化通常能把时间压到原来的三分之一用矩阵约束代替for循环。上面为了可读性用了for循环实际写代码时尽量把所有时段的约束拼接成大矩阵一次性传给constraints。YALMIP在拼接大规模符号矩阵时效率一般减少循环能明显提升解析速度。子问题求解时不要用默认的optimize改用optimizer对象把x_star作为参数传入省去每轮重建约束的开销。子问题的对偶变量数量通常和约束数量一致如果第二阶段约束很多先做约束化简把明显冗余的约束删掉。Gurobi的线程数设为4到8ops.gurobi.Threads默认是0自动在共享计算环境里可以显式设置。模型里大量使用abs()时YALMIP会引入辅助变量和不等式变量规模膨胀很快尽量在建模阶段就用正负分解来避开abs。5. 完全复现中的参数设定、结果校核与常见坑5.1 推荐参考输入参数参数设定直接影响收敛性和结果合理性下面这组参数是我在多个复现案例里校准过的参考值可以直接作为起点参数推荐值说明调度周期24时段1h/时段多阶段鲁棒最基本的设定风电装机500 kW预测误差取±15%光伏装机300 kW预测误差取±10%负荷峰值600 kW预测误差取±5%储能容量400 kWhSOC范围20%~95%储能最大充放电功率100 kW充放电效率均取0.95燃气轮机容量200 kW最小出力50 kW爬坡率50kW/h购电价格上限1.2 元/kWh分时电价峰谷差越大越能体现储能价值预算参数Γ3~5表示最多3~5个时段同时达到最大偏差收敛阈值1e-4绝对间隙最大迭代次数20实际通常8轮内收敛5.2 怎么判断你复现的模型确实鲁棒跑通模型只是第一步关键是要会验证。我复现完之后一般做三组校验任何一组不过都说明模型有问题第一组最坏场景重放。把CCG最后一轮得到的最坏场景u_cur代入第二阶段模型在固定第一阶段决策的情况下重新求解第二阶段调整问题。如果找不到可行解说明子问题建模或者切割添加有问题。第二组蒙特卡洛模拟。从风电、光伏、负荷的预测误差分布中随机抽样生成1000个场景对每个场景调用第二阶段模型检验是否可行。如果可行率低于98%说明不确定性集合设计偏激进需要调大Γ或扩大盒式范围。第三组Γ对比分析。分别取Γ0、1、3、6、12观察总成本和最坏场景成本的变化。Γ0时模型退化为确定性问题总成本最低随着Γ增大总成本单调上升这是符合预期的。如果曲线出现非单调或者剧烈跳变大概率是不确定集合的约束写错了。5.3 我实际踩过的几个大坑第一个坑是储能的SOC约束没有处理好充放电同时性。模型里P_ess_c和P_ess_d是两个独立变量如果不加z_ess互斥约束求解器会出现一个时段既充电又放电的情况利用效率空转降低成本看起来成本很漂亮实际上物理不可行。解决办法就是代码里用z_ess做互斥或者加入一个很小的同时充放电惩罚项。第二个坑是对偶变量符号搞错。子问题在对偶化时如果原约束是Ax Bu ≤ b对偶变量必须是非负的如果是等式约束对偶变量自由。符号一错对偶目标的值会偏离真实最坏场景成本导致UB和LB无法收敛甚至发散。第三个坑是无脑扩大M。在线性化双线性项时M取值太小会剪掉有效解M取值太大会导致数值病态。解决方法是尽量用不确定变量的区间上下界来推导M的理论值而不是随便取一个1e6。Gurobi的NumericFocus开启后有一定缓解但不解决根本问题。第四个坑是YALMIP的符号展开。模型复杂时YALMIP默认会把所有约束展开成完整系数矩阵变量多的时候内存占用会爆炸。设置sdpsettings(expand, 0)可以减轻这个问题但代价是求解器内部的presolve工作量增加需要实测权衡。5.4 复现后如何扩展这个模型跑通基础版之后有几个很自然的扩展方向都能在不推翻框架的前提下增加价值分布鲁棒优化把固定盒式集合换成Wasserstein球利用历史数据构建经验分布再以球内最坏分布为目标。这个扩展只需要改不确定集合的建模和对偶约束CCG框架可以复用。多能互补在电储能的基础上加入蓄热罐、燃气锅炉把区域微电网扩展成综合能源系统。多能耦合会引入热功率平衡约束但两阶段鲁棒的结构完全一致。碳约束在目标函数中加入碳排放惩罚或碳配额约束分析碳交易价格对鲁棒调度策略的影响。这类扩展对接双碳主题在实际汇报中更有说服力。储能容量优化把储能的额定容量和功率作为第一阶段决策变量目标函数中加入储能投资年化成本模型就变成了容量规划鲁棒调度联合优化问题CCG依然适用。我个人在多次复现这个模型后的体会是两阶段鲁棒框架本身并不算特别难难的是把每个环节的物理意义和数学细节对齐。比如Γ怎么解释给非优化背景的同事听比如储能SOC方程里为什么充电效率在乘法位置、放电效率在除法位置这些细节如果不理清楚代码复制过来改了参数结果不对劲都不知道去哪查。建议第一次复现的朋友先用小系统比如6节点、4时段把CCG迭代曲线画出来确认LB和UB的收敛趋势没问题了再放大到24时段完整系统。小系统里能一眼看出的逻辑错误放大以后往往会被数值噪音掩盖排查难度翻倍。这个模型跑通之后你再去读其他鲁棒调度相关的论文思路会清晰很多。 微电网多阶段鲁棒调度模型MATLAB复现与CCG算法实践
做微电网调度的复现我印象最深的一课是在MATLAB里把确定性经济调度跑通之后信心满满地把风光预测误差从±5%改到±30%结果同一套计划直接让储能越限、功率平衡失配。从那之后我就明白只要系统里装了可再生能源和储能调度模型里不考虑不确定性的最优大概率只是某个理想场景下的最优。这篇内容围绕的就是一套含可再生能源和储能的区域微电网多阶段鲁棒调度模型基于MATLAB完整复现核心思路是建立日前计划日内调整的两阶段鲁棒框架用区间集合刻画风光和负荷的预测偏差再通过列和约束生成算法CCG迭代求解最坏场景下的最优运行方案。文章会从模型为什么这样设计、不确定性集合怎么构建、CCG算法如何落地到MATLABYALMIPGurobi的代码骨架和复现时容易踩的坑完整过一遍适合正在做微电网优化调度、储能配置或鲁棒优化应用研究的朋友参考。1. 可再生渗透率攀升后传统确定性调度为什么越来越不好用1.1 确定性调度的假设与失效场景传统的微电网经济调度模型无论是目标函数还是约束条件都默认风电、光伏和负荷在某个时段的出力是一个已知的确定值。模型看起来非常干净目标函数是最小化运行成本包括向上级电网购电费用、微型燃气轮机的燃料费用、储能损耗成本约束包括功率平衡、储能SOC上下限、充放电功率限制、机组出力限制、联络线功率限制。这种模型在可再生能源渗透率比较低、系统备用充足的时候还能凑合用。但一旦风机和光伏在微电网里的装机占比上来了问题就暴露了。风电预测误差通常在15%到25%之间波动光伏在天气骤变时误差可以超过30%负荷预测虽然相对准但在极端天气下也会出现明显偏离。确定性调度把所有预测值当成真实值来处理等于把全部风险都压在了预测完全准确这个脆弱的假设上。我做过一个对比测试在一个24节点微电网测试系统里用确定性模型求出一组日前计划然后把实际风光出力设定为预测值的70%和130%两种偏差场景重新跑功率平衡校验结果有多组约束直接违反储能SOC甚至越过了下限。这就是典型的计划赶不上变化。1.2 随机优化、机会约束与鲁棒优化三条技术路线的实际取舍面对不确定性学术界和工程界常用的解决方案有三类很多初学者一上来容易混在一起方法基本思路优点主要问题随机规划用大量场景近似不确定性对场景求期望最优结果贴合概率分布经济性较好需要精确的概率分布场景数量大时计算量爆炸分布估计不准时结果不可靠机会约束规划允许约束以较高概率满足如95%概率满足功率平衡可以在安全性和经济性之间折中概率参数如何标定比较主观分布信息不足时难以求解鲁棒优化假设不确定性落在某个集合内保证最坏情况下约束也满足不需要精确概率分布结果对极端场景稳健计算效率相对可控结果偏保守成本一般高于随机规划在微电网这种规模不大、但安全约束要求高的场景下鲁棒优化是比较务实的选择。原因也很直接你通常拿不到足够多的历史数据来精确拟合风电和光伏误差的分布函数就算拟合出来分布参数本身也有估计误差。与其在概率上赌一个精度不如直接说我知道误差不会超过某个范围在这个范围内怎么折腾我的调度方案都兜得住。1.3 标题里鲁棒性和不确定性到底指什么标题里的鲁棒性英文对应Robustness放在调度模型里指的就是当实际的风电、光伏出力偏离预测值时系统依然能通过储能充放电调整、机组出力再调度、与主网功率交换等手段维持功率平衡并满足所有安全约束。衡量鲁棒性的标准不是平均工况下的成本而是最坏工况下的约束满足能力。不确定性在这个问题里主要来自三个方向源端不确定性风机出力和光伏出力受气象条件影响预测模型再准也有偏差荷端不确定性区域负荷受生产生活规律、天气温度等因素影响同样存在预测误差设备状态不确定性比如某台机组临时故障、某条线路检修但在调度主问题里通常把这类视为N-1预想故障单独处理不会混进源荷不确定性里。在工程复现时通常把源荷不确定性统一建模为一个不确定集合U集合里的每一个元素代表一组可能的源荷出力偏差向量。鲁棒优化要求对U里的任意一个偏差向量系统都能通过第二阶段调整找到可行解。这就是后面整个模型的核心出发点。2. 多阶段鲁棒调度模型的数学框架与不确定性集合设计2.1 两阶段模型的前后配合日前决策与日内调整标题里写的是多阶段鲁棒调度模型实际复现时最常见的落地形式是两阶段鲁棒优化模型也就是经典的三层结构min-max-min。之所以叫多阶段是因为调度决策本身是分层分时段的第一阶段日前决策阶段在获知次日风光和负荷预测值之后提前决定机组的启停状态、储能充放电基准计划、向上级电网购电的日前计划。这些决策的特点是必须在看到实际偏差之前做出所以也叫here-and-now决策。第二阶段日内调整阶段当实际风光出力、负荷与预测出现偏差后在日前计划的基础上通过调节储能出力、机组出力微调、可削减负荷等手段消除功率不平衡。这个阶段的决策是在偏差揭晓之后才能做的叫wait-and-see决策。两阶段之间通过功率平衡约束和设备运行约束耦合在一起。模型的目标函数也相应地拆成两块min (第一阶段运行成本) max_{u∈U} min_{y∈Ω(x,u)} (第二阶段调整成本)外层min决定第一阶段决策x内层max寻找使第二阶段成本最高的最坏偏差场景u中间层min则是在给定x和u之后寻找最优的调整策略y。整个模型求的是在所有可能出现在的偏差场景下系统总成本的最坏情况最优值。2.2 盒式、预算约束与椭球式不确定集合的适用场景不确定集合怎么建直接决定了模型的保守程度和求解难度。常见的有三种形式盒式集合Box直接把每个时段的不确定量限制在一个区间内如u(t)∈[u_min(t), u_max(t)]。这种集合最简单但要求所有时段同时达到最大偏差结果是最坏中的最坏过于保守。预算约束集合Budget在盒式集合基础上增加一个总偏差比例限制形如Σ|u(t)-u_hat(t)|/δ(t) ≤ Γ。这里的Γ就是预算参数表示在调度周期内最多允许多少个时段的偏差同时达到上限。Γ越小鲁棒约束越宽松成本越低Γ越大方案越保守。这个设置非常贴合实际因为真实的风光出力偏差不太可能24个小时每个小时都顶在最大值上。椭球式集合Ellipsoidal用2-范数约束替代1-范数额外限制形如‖W(u-u_hat)‖₂ ≤ Ω。椭球集合考虑了变量之间的相关性数学上更精细但会引入二次约束计算复杂度明显上升。在微电网调度里我最推荐的是盒式预算约束的组合。原因有三个物理意义清晰Γ可以理解为最坏情况下同时偏离预测的时间段数方便向工程人员解释模型保持线性结构经过对偶变换后仍然是LP或MILP可以用YALMIPGurobi稳定求解通过调节Γ的取值可以画出从激进到保守的完整成本曲线这是论文和报告里很有说服力的一张图。2.3 目标函数与约束的完整数学表达以含风电、光伏、储能和微型燃气轮机的区域微电网为例完整的鲁棒调度模型可以写成决策变量第一阶段燃气轮机开机状态z(t)、储能充电功率P_ess_c(t)、储能放电功率P_ess_d(t)、日前购电功率P_grid(t)、储能SOC(t)。第二阶段针对偏差场景u的储能调整量ΔP_ess(t)、机组调整量ΔP_g(t)、可削减负荷量ΔP_load(t)。目标函数min Σ_t [ C_grid * P_grid(t) C_g * P_g(t) C_ess * (P_ess_c(t)P_ess_d(t)) ] max_{u∈U} min_Δ Σ_t [ C_penalty * (ΔP_load(t) 弃风弃光惩罚项) ]约束条件功率平衡约束 P_grid(t) P_wind(t) P_pv(t) P_g(t) P_ess_d(t) P_load(t) P_ess_c(t) ΔP_load(t)其中P_wind(t)和P_pv(t)在第二阶段表示为P_wind_hat(t) u_w(t)、P_pv_hat(t) u_pv(t)。储能约束 SOC(t1) SOC(t) η_c * P_ess_c(t) * Δt - P_ess_d(t) * Δt / η_d SOC_min ≤ SOC(t) ≤ SOC_max 0 ≤ P_ess_c(t) ≤ P_ess_c_max * z_ess(t) 0 ≤ P_ess_d(t) ≤ P_ess_d_max * (1 - z_ess(t))这里需要注意的是充放电状态变量z_ess(t)的引入它保证储能不会同时充电和放电这在后面代码实现时是个关键细节。机组约束 z(t) * P_g_min ≤ P_g(t) ≤ z(t) * P_g_max 爬坡约束|P_g(t) - P_g(t-1)| ≤ Ramp_g与主网交互约束 P_grid_min ≤ P_grid(t) ≤ P_grid_max可削减负荷约束 0 ≤ ΔP_load(t) ≤ α * P_load(t)2.4 为什么多阶段不等于多周期很多初学者会把多阶段鲁棒理解成把24小时的调度周期分成很多段每段做一个鲁棒优化这是一个明显的误区。多周期优化只是把时间维度拉长决策变量在整个周期内一次性规划本质上仍是开环调度。多阶段鲁棒优化的关键不在于时间段的多少而在于决策的可调整性——第一阶段决策做出后可以看到不确定性的实现然后第二阶段再做出调整决策。这个特征对应的是闭环调度也就是日前计划日内实时调整的真实运行流程。用数学语言说真正的多阶段鲁棒模型必须具有min-max-min的嵌套结构而不是简单的min-max或者包含很多时段的一次性min。在复现时如果你的模型只是每个时段都预留备用容量然后用鲁棒约束包住所有时段这虽然也算鲁棒但不是多阶段鲁棒。要体现多阶段就必须显式区分可调整变量如储能的实时充放电调整量和不可调整变量如机组启停、日前购电计划并在目标函数中体现第二阶段对最坏场景的应对成本。3. 列和约束生成CCG算法迭代框架的收敛逻辑3.1 主问题的最小费用加最坏场景切割线两阶段鲁棒模型直接求解极其困难因为内层嵌套了max-min优化。工程上最常用的解法是列和约束生成算法核心思路是把原问题拆成一个主问题MP和一个子问题SP通过迭代逐步逼近最优解。主问题的形式如下min c * x η s.t. 第一阶段约束 η ≥ d * y_l, 对于已知的场景u_l x、y_l之间的耦合约束这里的关键是引入了辅助变量η它表示第二阶段成本的上界估计。每轮迭代中子问题会找到一个最坏场景u_l主问题就把这个场景对应的第二阶段变量y_l和约束作为切割加入模型让η不断收紧。主问题因为只考虑了已发现的最坏场景所以实际上是原问题的松弛版本它的最优值给出原问题的下界LB。3.2 子问题的max-min对偶转化子问题是在给定第一阶段解x_star之后寻找使第二阶段成本最大的场景u以及对应的最优调整成本形式是f(x_star) max_{u∈U} min_{y∈Ω(x_star, u)} d * y这个max-min结构不能直接丢给求解器。标准做法是对内层的min问题取强对偶。由于内层是关于y的线性规划只要可行域非空且有界强对偶成立就可以把min问题等价转化为其对偶max问题然后和外层的max合并变成一个单层max问题f(x_star) max_{u∈U, λ∈Λ} λ * (b - A * x_star - C * u)其中λ是对偶变量它的符号由原约束方向决定。如果原约束是≤型λ非负如果是等式约束λ自由变量。合并之后子问题变成带二次项的非线性问题因为λ和u相乘。在实际求解中这个双线性项通常可以用大M法或直接交给支持MIQP的非线性求解器处理Gurobi可以直接求解这类带双线性项的小规模问题。需要强调一点内层min问题必须满足强对偶条件。因此第二阶段调整问题必须建模为线性规划不能随意加入整数变量。如果第二阶段需要决策是否切除某条线路这类0-1变量标准CCG就不适用了需要改成嵌套CCG或者引入其他分解方法。3.3 收敛判据与最大迭代次数CCG算法的主循环流程如下初始化UBinf, LB-inf, k0设置初始场景u_0通常取预测均值 while k ≤ K_max 且 UB - LB tol: 1. 求解主问题得到x_k*和η_k*更新LB c*x_k* η_k* 2. 将x_k*代入子问题求解最坏场景u_k*和最坏场景成本f(x_k*) 3. 更新UB min(UB, c*x_k* f(x_k*)) 4. if UB - LB ≤ tol: 收敛退出 5. 否则将u_k*作为新场景加入主问题添加对应的第二阶段变量y_k和割约束 k k 1收敛判据一般取相对间隙或绝对间隙。我在复现时常用的组合是绝对误差取1e-4最大迭代次数取20轮。在24时段微电网模型里通常3到8轮就能收敛到很紧的间隙很少会真的跑到20轮。3.4 CCG与Benders分解的选择建议很多文献会同时提到Benders分解和CCG两者都是求解两阶段问题和分解主-子问题的经典算法但适用场景有明显差异对比维度Benders分解CCG切割信息通过对偶函数添加割平面直接添加最坏场景对应的完整约束迭代收敛对小规模问题收敛较快通常需要更少迭代次数尤其有连续变量时对整数变量的适应性一般需要处理整数对偶gap主问题中可以自然处理整数变量实现难度较低中等在微电网鲁棒调度里第一阶段通常包含机组启停这类整数变量而CCG的切割是把具体场景u_l下的约束完整塞进主问题对整数变量天然友好所以实际复现时我基本都用CCG。只有当子问题规模特别大、每次求解SP耗时很高时才会回头考虑Benders或者近似方法。4. MATLABYALMIPGurobi的关键代码实现4.1 环境准备与求解器适配复现前先确认环境MATLAB版本R2020b及以上即可建议R2022b以上对sdpvar和optimizer的兼容性更好YALMIP直接从GitHub拉取最新release不要用十年前的老版本求解器Gurobi 10.x或11.x学术许可免费申请Mosek 10.x也可以但Gurobi在MILP上的速度通常更优。在MATLAB里先跑一句验证yalmiptest能列出Gurobi的LP、QP、MILP、MIQP支持状态就说明环境没问题。如果用了旧版MATLAB配新版Gurobi经常会出现java相关报错那是MATLAB版本太老、Java版本不匹配导致的属于环境问题先升级MATLAB再排查。4.2 主问题的YALMIP建模代码下面是主问题的核心建模骨架采用24时段、1小时为步长的设定。nT 24; % 第一阶段变量 z_g binvar(1, nT); % 燃气轮机启停 P_g sdpvar(1, nT); % 燃气轮机出力 P_ess_c sdpvar(1, nT); % 储能充电功率 P_ess_d sdpvar(1, nT); % 储能放电功率 z_ess binvar(1, nT); % 储能充放电状态1充电0放电 SOC sdpvar(1, nT1); % 荷电状态 P_grid sdpvar(1, nT); % 购电功率 eta sdpvar(1); % 最坏场景成本上界 % 目标函数 objective_mp sum(grid_price .* P_grid) ... sum(g_fuel_cost * P_g) ... sum(ess_cost * (P_ess_c P_ess_d)) ... eta; % 约束 constraints_mp []; constraints_mp [constraints_mp, SOC(1) SOC_init]; for t 1:nT constraints_mp [constraints_mp, ... SOC(t1) SOC(t) eta_charge * P_ess_c(t) - P_ess_d(t) / eta_discharge, ... SOC_min SOC(t1) SOC_max, ... 0 P_ess_c(t) P_ess_c_max * z_ess(t), ... 0 P_ess_d(t) P_ess_d_max * (1 - z_ess(t)), ... z_g(t) * P_g_min P_g(t) z_g(t) * P_g_max, ... P_grid_min P_grid(t) P_grid_max, ... P_wind_forecast(t) P_pv_forecast(t) P_g(t) P_grid(t) P_ess_d(t) ... P_load_forecast(t) P_ess_c(t)]; end ops sdpsettings(solver, gurobi, verbose, 0, gurobi.NumericFocus, 3); optimize(constraints_mp, objective_mp, ops);这里有个细节给Gurobi设置NumericFocus为3可以明显降低大M约束带来的数值病态问题。如果模型里用到了大M这个设置几乎是必加的。4.3 子问题的对偶化与最坏场景求解代码第二阶段的最坏场景求解是核心代码分三步走第一步固定第一阶段变量第二步构造第二阶段LP并取对偶第三步求解单层max问题。% 固定第一阶段变量 x_star value([P_g; P_grid; P_ess_c; P_ess_d; SOC]); % 不确定变量 u_w sdpvar(1, nT); % 风电偏差 u_pv sdpvar(1, nT); % 光伏偏差 u_l sdpvar(1, nT); % 负荷偏差 % 第二阶段调整变量 delta_ess sdpvar(1, nT); % 储能功率调整 delta_g sdpvar(1, nT); % 机组出力调整 delta_load sdpvar(1, nT); % 可削减负荷 % 不确定集合盒式预算 U_con [-u_w_max u_w u_w_max, ... -u_pv_max u_pv u_pv_max, ... -u_l_max u_l u_l_max, ... sum(abs(u_w) ./ u_w_max) sum(abs(u_pv) ./ u_pv_max) Gamma]; % 第二阶段约束仅对应给定x_star SP_con [...]; SP_obj sum(penalty_load * delta_load penalty_curtail * (u_w u_pv) ... C_delta_g * abs(delta_g) C_delta_ess * abs(delta_ess)); % 通过YALMIP取对偶或直接用求解器求解最大最小问题 % 更通用的做法是手动构造对偶变量把max-min转为单层max dual_lambda sdpvar(size(A_sp, 1), 1); dual_mu sdpvar(size(A_eq, 1), 1); % 对偶目标函数注意线性化双线性项 dual_obj dual_lambda * (b_sp - A_sp * x_star_vec - B_u * [u_w; u_pv; u_l]) ... dual_mu * (b_eq - A_eq * x_star_vec); % 合并后的最坏场景问题 SP_merged_con [dual_lambda 0, dual_mu free, U_con]; result_sp optimize(SP_merged_con, -dual_obj, ops_sp);实际上当双线性项存在时Gurobi会以MIQP方式求解变量数不大时速度可以接受。如果模型规模较大建议用big-M法把双线性项线性化把问题转成MILP。4.4 迭代主循环与结果保存主循环代码框架如下UB inf; LB -inf; k 0; tol 1e-4; Kmax 20; cut_results {}; while k Kmax (UB - LB) tol k k 1; % 1. 求解主问题 optimize(constraints_mp, objective_mp, ops); LB value(objective_mp); x_cur value([P_g; P_grid; P_ess_c; P_ess_d; SOC]); % 2. 求解子问题 optimize(SP_merged_con, -dual_obj, ops_sp); sp_cost value(dual_obj); u_cur value([u_w; u_pv; u_l]); day_ahead_cost value(sum(grid_price .* P_grid) ... sum(g_fuel_cost * P_g) ... sum(ess_cost * (P_ess_c P_ess_d))); UB min(UB, day_ahead_cost sp_cost); fprintf(迭代%d: LB%.4f, UB%.4f, gap%.6f\n, k, LB, UB, UB-LB); % 3. 添加新场景切割到主问题 y_new sdpvar(size_y, 1); constraints_mp [constraints_mp, ... y_new对应的第二阶段约束把u_cur作为常数代入, ... eta 第二阶段成本表达式(y_new)]; cut_results{k}.u u_cur; cut_results{k}.cost sp_cost; end有几个实现上的细节值得说不要把主问题的sdpvar变量每轮重新创建。变量结构要保持不变只是追加约束。否则MATLAB每轮都要重新解析整个符号表达式速度直接掉一个数量级。主问题的optimize结果如果返回problemsolver 信息先打印optimize的第二个输出看看是primal infeasible还是numerical issues再决定调参数还是调模型。保存每轮的最坏场景u_cur最后画图时可以直接把最坏场景下的风电光伏曲线和基准预测曲线对比很直观。4.5 代码性能优化的几个细节复现完跑第一版10轮迭代可能要跑几分钟如果觉得慢按下面这几个方向优化通常能把时间压到原来的三分之一用矩阵约束代替for循环。上面为了可读性用了for循环实际写代码时尽量把所有时段的约束拼接成大矩阵一次性传给constraints。YALMIP在拼接大规模符号矩阵时效率一般减少循环能明显提升解析速度。子问题求解时不要用默认的optimize改用optimizer对象把x_star作为参数传入省去每轮重建约束的开销。子问题的对偶变量数量通常和约束数量一致如果第二阶段约束很多先做约束化简把明显冗余的约束删掉。Gurobi的线程数设为4到8ops.gurobi.Threads默认是0自动在共享计算环境里可以显式设置。模型里大量使用abs()时YALMIP会引入辅助变量和不等式变量规模膨胀很快尽量在建模阶段就用正负分解来避开abs。5. 完全复现中的参数设定、结果校核与常见坑5.1 推荐参考输入参数参数设定直接影响收敛性和结果合理性下面这组参数是我在多个复现案例里校准过的参考值可以直接作为起点参数推荐值说明调度周期24时段1h/时段多阶段鲁棒最基本的设定风电装机500 kW预测误差取±15%光伏装机300 kW预测误差取±10%负荷峰值600 kW预测误差取±5%储能容量400 kWhSOC范围20%~95%储能最大充放电功率100 kW充放电效率均取0.95燃气轮机容量200 kW最小出力50 kW爬坡率50kW/h购电价格上限1.2 元/kWh分时电价峰谷差越大越能体现储能价值预算参数Γ3~5表示最多3~5个时段同时达到最大偏差收敛阈值1e-4绝对间隙最大迭代次数20实际通常8轮内收敛5.2 怎么判断你复现的模型确实鲁棒跑通模型只是第一步关键是要会验证。我复现完之后一般做三组校验任何一组不过都说明模型有问题第一组最坏场景重放。把CCG最后一轮得到的最坏场景u_cur代入第二阶段模型在固定第一阶段决策的情况下重新求解第二阶段调整问题。如果找不到可行解说明子问题建模或者切割添加有问题。第二组蒙特卡洛模拟。从风电、光伏、负荷的预测误差分布中随机抽样生成1000个场景对每个场景调用第二阶段模型检验是否可行。如果可行率低于98%说明不确定性集合设计偏激进需要调大Γ或扩大盒式范围。第三组Γ对比分析。分别取Γ0、1、3、6、12观察总成本和最坏场景成本的变化。Γ0时模型退化为确定性问题总成本最低随着Γ增大总成本单调上升这是符合预期的。如果曲线出现非单调或者剧烈跳变大概率是不确定集合的约束写错了。5.3 我实际踩过的几个大坑第一个坑是储能的SOC约束没有处理好充放电同时性。模型里P_ess_c和P_ess_d是两个独立变量如果不加z_ess互斥约束求解器会出现一个时段既充电又放电的情况利用效率空转降低成本看起来成本很漂亮实际上物理不可行。解决办法就是代码里用z_ess做互斥或者加入一个很小的同时充放电惩罚项。第二个坑是对偶变量符号搞错。子问题在对偶化时如果原约束是Ax Bu ≤ b对偶变量必须是非负的如果是等式约束对偶变量自由。符号一错对偶目标的值会偏离真实最坏场景成本导致UB和LB无法收敛甚至发散。第三个坑是无脑扩大M。在线性化双线性项时M取值太小会剪掉有效解M取值太大会导致数值病态。解决方法是尽量用不确定变量的区间上下界来推导M的理论值而不是随便取一个1e6。Gurobi的NumericFocus开启后有一定缓解但不解决根本问题。第四个坑是YALMIP的符号展开。模型复杂时YALMIP默认会把所有约束展开成完整系数矩阵变量多的时候内存占用会爆炸。设置sdpsettings(expand, 0)可以减轻这个问题但代价是求解器内部的presolve工作量增加需要实测权衡。5.4 复现后如何扩展这个模型跑通基础版之后有几个很自然的扩展方向都能在不推翻框架的前提下增加价值分布鲁棒优化把固定盒式集合换成Wasserstein球利用历史数据构建经验分布再以球内最坏分布为目标。这个扩展只需要改不确定集合的建模和对偶约束CCG框架可以复用。多能互补在电储能的基础上加入蓄热罐、燃气锅炉把区域微电网扩展成综合能源系统。多能耦合会引入热功率平衡约束但两阶段鲁棒的结构完全一致。碳约束在目标函数中加入碳排放惩罚或碳配额约束分析碳交易价格对鲁棒调度策略的影响。这类扩展对接双碳主题在实际汇报中更有说服力。储能容量优化把储能的额定容量和功率作为第一阶段决策变量目标函数中加入储能投资年化成本模型就变成了容量规划鲁棒调度联合优化问题CCG依然适用。我个人在多次复现这个模型后的体会是两阶段鲁棒框架本身并不算特别难难的是把每个环节的物理意义和数学细节对齐。比如Γ怎么解释给非优化背景的同事听比如储能SOC方程里为什么充电效率在乘法位置、放电效率在除法位置这些细节如果不理清楚代码复制过来改了参数结果不对劲都不知道去哪查。建议第一次复现的朋友先用小系统比如6节点、4时段把CCG迭代曲线画出来确认LB和UB的收敛趋势没问题了再放大到24时段完整系统。小系统里能一眼看出的逻辑错误放大以后往往会被数值噪音掩盖排查难度翻倍。这个模型跑通之后你再去读其他鲁棒调度相关的论文思路会清晰很多。