尧图网页设计 
1. 微电网两阶段鲁棒优化到底解决什么问题1.1 从确定性调度到不确定性调度做微电网优化调度的人最早接触的模型基本都是确定性优化光伏出力取预测值、负荷取预测值、电价取已知的分时电价然后丢给求解器一把梭得到一个“最优”的机组出力和储能充放电计划。这个做法最大的问题在于——预测永远不可能完全准确。光伏在晴天中午可能预报偏差只有5%到了多云、雷阵雨天气15分钟级滚动预测的误差轻轻松松超过20%。等到了实际运行那天你拿着昨天算好的调度计划发现光伏根本没发到预测值储能该充电的时候在放电该放电的时候又没电可放整个计划直接失效。两阶段鲁棒优化Two-Stage Robust Optimization的核心思路就是不再把不确定性参数当成一个确定的数而是给它们划一个范围——一个不确定集Uncertainty Set。所有可能出现的风光出力、负荷波动都被约束在这个集合内然后优化的目标变成在最恶劣的不确定性实现下依然能让系统运行成本最低、且所有约束都满足。换句话说它求的不是“预测场景下的最优”而是“最坏情况下的最优”。这个思路在实际工程里非常重要。微电网往往接在配电网末端系统容量小、惯性小光伏和负荷的小幅波动就可能造成电压越限或者功率倒送。你不可能指望运行人员每天盯着天气云图手动调整机组出力必须让调度策略自带“抗扰动能力”。两阶段鲁棒优化恰恰就是干这个的它能在不确定性到来之前先做好防御性的决策安排。1.2 为什么用两阶段而不是单阶段或随机规划单阶段鲁棒优化是最保守的模型——所有决策都在不确定性实现之前敲定并且要求对不确定集内的所有场景都可行。它的问题是过于悲观算出来的运行成本往往高得离谱储能基本不敢用柴发一直开在最高出力实际运行中根本没人愿意接受这种方案。随机规划是另一种思路它给不确定参数赋概率分布然后优化期望成本。这个方向理论上很美但工程落地有个尴尬问题光伏出力的概率分布你很难精确刻画尤其是短时波动往往不是标准正态分布用场景法做的话需要生成大量场景才能收敛计算量爆炸式增长而且场景缩减Scenario Reduction操作不当还可能改变问题的可行域结构。两阶段鲁棒优化把决策拆成“现在做”和“等不确定实现后再做”两步。第一阶段决策Here-and-Now是在不确定性没有揭开时就必须拍板的量比如机组启停状态、与主网的购售电协议、储能的日前充放电计划第二阶段决策Wait-and-See是在不确定性实现后可以调整的量比如具体每台机组的出力微调、储能的实时充放电功率、弃光弃风量。这样做的好处是既保留了鲁棒优化“保证所有场景可行”的安全性又不像单阶段那样僵硬——第二阶段可以根据实际情况灵活调整成本自然不会那么离谱。1.3 一个直观的类比打个比方这就像你计划明天出门旅行。确定性优化相当于你完全相信天气预报说明天不下雨于是只带了一顶遮阳帽单阶段鲁棒优化相当于你假设明天可能下暴雨、刮台风、下冰雹甚至地震于是把整个家都搬上车后备箱都盖不上两阶段鲁棒优化则是你先带上雨伞、冲锋衣这种能应付大多数恶劣天气的装备第一阶段决策等第二天早上真正看到天气情况后再决定穿哪件衣服、走哪条路第二阶段决策——既保证了不会淋雨又不会出门带太多累赘。2. 模型构建的核心逻辑2.1 第一阶段与第二阶段的决策变量划分建模之前最需要想清楚的问题就是哪些变量放第一阶段哪些变量放第二阶段。这个划分直接决定了模型的保守程度和求解难度。在我做的这个微电网系统里典型配置是光伏阵列、风力发电机、蓄电池储能、柴油发电机或者微型燃气轮机、与外部电网的联络线以及本地负荷。第一阶段决策变量选的是柴发的启停状态0-1整数变量、储能是否处于充电/放电状态的二进制变量、从外部电网购电/售电的状态变量以及与主网的交换功率计划值。这些变量对应的是日前时间尺度需要在不确定性揭晓之前就确定下来。第二阶段决策变量则是每台机组的实际出力调整量、储能的实际充电/放电功率、弃光弃风率、负荷削减量、与主网实际交换功率的偏差修正值。这些量可以在实时运行中根据实际光伏出力和负荷情况做调整。需要特别提醒的是储能是否充放电这个0-1变量最好不要放在第一阶段除非你有特殊需求。原因是这会让模型变成混合整数两阶段鲁棒优化MIRU求解复杂度直接上一个台阶。很多时候可以采用一个简化技巧允许储能在第二阶段连续调整充放电功率但通过一个足够大的M值约束防止它同时充放电这样主问题Master Problem依然是个MILP但子问题Subproblem变成了纯LP处理起来省事得多。2.2 不确定集的设计与参数调节不确定集的选择是两阶段鲁棒优化建模的灵魂。常用的有盒式不确定集Box、椭球不确定集Ellipsoidal和多面体不确定集Polyhedral。盒式不确定集形式最简单每个不确定参数独立波动在其上下界之间即 $u \in [\hat{u} - \hat{u}\Delta, \hat{u} \hat{u}\Delta]$。但是这种集合的保守度很高因为它允许所有不确定参数同时取到最坏值而现实中这种情况概率极低。椭球不确定集考虑了参数之间的相关性形状更贴近实际但问题会变成二阶锥规划SOCP甚至更复杂的形式工程实现门槛高。多面体不确定集也叫预算不确定集Budget Uncertainty Set在盒式集合的基础上增加一个总偏差预算约束例如 $\sum_i |u_i - \hat{u}_i| / \hat{u}_i \leq \Gamma$这个 $\Gamma$ 就是鲁棒预算Budget of Uncertainty。我在代码里采用的正是多面体不确定集因为它在工程实用性和保守度之间取得了很好的平衡。$\Gamma$ 的物理含义是“最多同时允许多少个不确定参数偏离预测值”取0时就是确定性模型取值越大就越保守。实际操作中我习惯先跑一次确定性模型得到基准成本然后逐步增大 $\Gamma$观察总成本的变化曲线。如果成本增长斜率在某个点之后突然变陡那说明再增加保守度就不划算了——这个拐点附近就是 $\Gamma$ 的最佳取值区间。这里分享一个我踩过的坑盒式集合在某几个子问题上解出来的目标值偏大一开始我还以为是求解器数值问题后来排查了半天才发现是不确定参数之间的相关性没有建模导致最坏场景同时乐观估计了充电电价又悲观估计了放电电价双重叠加把成本推高了。改成多面体集合加上预算约束后结果就正常了。2.3 目标函数与约束的矩阵化表达微电网经济调度的目标函数一般包含以下几项柴油发电机燃料成本、启停成本、从主网购电成本减去售电收益、储能充放电循环老化成本可选以及弃光弃风惩罚项可选如果允许弃电的话。两阶段鲁棒优化的标准形式是$$\min_{x} \ c^T x \max_{u \in \mathcal{U}} \min_{y \in \Omega(x, u)} d^T y$$其中 $x$ 是第一阶段变量$y$ 是第二阶段变量$u$ 是不确定参数。内层 $\max \min$ 表示在给定第一阶段决策后寻找最恶劣的不确定实现然后在最恶劣情况下最小化运行成本。对应地约束分为几个层次第一阶段固有约束机组启停逻辑约束、最小启停时间约束、储能充放电状态互斥约束、购售电状态互斥约束。耦合约束功率平衡约束中既含第一阶段变量也含第二阶段变量还有储能SOC的状态转移方程。第二阶段约束机组出力上下限、爬坡约束、储能充放电功率限制、SOC上下限、与主网交换功率限制、弃电量约束。写Matlab代码的时候我强烈建议把所有约束整理成矩阵形式而不是标量循环。原因有两个一是YALMIP对矩阵化约束的预处理效率远高于逐条添加二是调试的时候你能直接看约束矩阵的维度是否匹配问题定位快得多。比如功率平衡约束如果你有24个时段的优化直接构造一个24行的等式约束矩阵而不是写24次constraint [constraint; ...]。3. Matlab代码的实现思路与关键细节3.1 整体代码结构与主循环框架拿到一个两阶段鲁棒优化问题最常用的求解方法是列与约束生成算法Column-and-Constraint GenerationCCG也常称为Benders分解的变体。核心思想是把原问题分解成一个主问题MP和一个子问题SP迭代求解。主问题是在已知的有限个最恶劣场景下做优化得到一个下界对最小化问题而言子问题是在给定主问题解的基础上寻找新的最恶劣不确定场景并给主问题上界。整个代码框架我分成了几个模块main.m主程序负责设置系统参数、调用求解循环、输出结果。data_define.m所有基础数据的定义——负荷曲线、光伏预测出力、电价曲线、储能参数、柴发参数、不确定集参数。build_MP.m构建主问题的YALMIP模型。solve_SP.m求解子问题找出当前迭代的最恶劣场景。plot_results.m将迭代过程中的调度结果可视化。CCG主循环的逻辑如下% 初始化 LB -inf; UB inf; iter 1; max_iter 10; % 或者用收敛判据 best_scenarios {}; % 存储所有已发现的最恶劣场景 while (UB - LB) / abs(UB) 1e-4 iter max_iter % 1. 求解主问题加入已知的最恶劣场景 [x_opt, obj_MP] solve_MP(best_scenarios); LB obj_MP; % 2. 代入主问题解求解子问题 [u_new, obj_SP, y_opt] solve_SP(x_opt); UB min(UB, obj_MP - obj_SP); % 注意这里根据对偶形式调整 % 3. 如果子问题目标值对应的场景尚未覆盖则加入场景并继续 if abs(obj_SP) 1e-6 best_scenarios{end1} u_new; end % 4. 检查收敛 fprintf(Iteration %d: LB%.4f, UB%.4f, gap%.4f%%\n, iter, LB, UB, ... abs(UB - LB) / abs(UB) * 100); iter iter 1; end有一点需要注意在标准的两阶段鲁棒优化中子问题通常是“给定第一阶段决策 $x^*$寻找最大化的最恶劣场景下的最小运行成本”。如果子问题是LP可以通过KKT条件或者对偶变换把它转成单层优化如果子问题含有整数变量就得用分解算法或者大M法线性化。在我这个模型中第二阶段只含连续变量所以子问题通过强对偶转成单层的max问题然后用YALMIPGurobi一步求出最恶劣场景和对应的成本。3.2 YALMIP建模与求解器设置代码实现层面我用的是 YALMIP 作为建模语言求解器选的是 Gurobi实测比CPLEX在MIQP上更快也可以用SCIP或者CBC做开源替代。一个完整的MP构建代码示意如下function [Constraints, Obj, x_vars] build_MP(best_scenarios) % 第一阶段变量 x_vars.u_on binvar(24, 1); % 柴发启停 x_vars.x_charge binvar(24, 1); % 储能充电状态 x_vars.x_discharge binvar(24, 1); % 储能放电状态 x_vars.x_buy binvar(24, 1); % 购电状态 x_vars.x_sell binvar(24, 1); % 售电状态 x_vars.p_buy sdpvar(24, 1); % 购电功率 x_vars.p_sell sdpvar(24, 1); % 售电功率 x_vars.soc sdpvar(24, 1); % 储能荷电状态 x_vars.p_sto_out sdpvar(24, 1); % 储能放电功率 x_vars.p_sto_in sdpvar(24, 1); % 储能充电功率 x_vars.p_dg sdpvar(24, 1); % 柴发出力 Constraints []; % 储能状态互斥 Constraints [Constraints, x_vars.x_charge x_vars.x_discharge 1]; Constraints [Constraints, x_vars.x_buy x_vars.x_sell 1]; % SOC递推 for t 2:24 Constraints [Constraints, ... x_vars.soc(t) x_vars.soc(t-1) ... x_vars.p_sto_in(t) * eta_ch / Cap - ... x_vars.p_sto_out(t) / eta_dis / Cap]; end % 购售电功率上限大M形式 Constraints [Constraints, x_vars.p_buy x_vars.x_buy * P_grid_max]; Constraints [Constraints, x_vars.p_sell x_vars.x_sell * P_grid_max]; % 柴发出力范围 Constraints [Constraints, x_vars.p_dg x_vars.u_on * P_dg_max]; Constraints [Constraints, x_vars.p_dg x_vars.u_on * P_dg_min]; % 对每一个已知的最恶劣场景建立第二阶段约束 Obj sum(fuel_cost(x_vars.p_dg)) sum(start_cost .* max(0, diff([0; x_vars.u_on]))) ... sum(price_buy .* x_vars.p_buy) - sum(price_sell .* x_vars.p_sell); for k 1:length(best_scenarios) u_k best_scenarios{k}; % 第k个最恶劣场景 % 对应场景下的第二阶段变量 y_pv_curtail sdpvar(24, 1); y_load_curtail sdpvar(24, 1); y_dg_adjust sdpvar(24, 1); % 功率平衡约束 Constraints [Constraints, ... u_k.p_pv u_k.p_wind x_vars.p_dg y_dg_adjust ... x_vars.p_sto_out - x_vars.p_sto_in x_vars.p_buy - x_vars.p_sell ... y_load_curtail u_k.p_load - y_pv_curtail]; % 弃光限制 Constraints [Constraints, 0 y_pv_curtail u_k.p_pv]; Constraints [Constraints, y_load_curtail 0]; end end这段代码只展示了核心思路实际项目中还需要考虑柴发最小启停时间、储能SOC初值设定等细节。但基本框架就是上面这样——主问题每轮迭代都会重新构建一次把新发现的场景追加到约束里。子问题的求解稍微复杂一点。给定第一阶段变量 $x^*$对偶问题可以写成function [u_new, obj_SP] solve_SP(x_vars) % 不确定变量 u_pv sdpvar(24, 1); u_load sdpvar(24, 1); % 多面体不确定集约束 Constraints [Constraints, ... u_pv pv_forecast * (1 - delta_pv), ... u_pv pv_forecast * (1 delta_pv), ... u_load load_forecast * (1 - delta_load), ... u_load load_forecast * (1 delta_load)]; % 预算约束 Constraints [Constraints, ... sum(abs(u_pv - pv_forecast) / (pv_forecast * delta_pv 1e-6)) ... sum(abs(u_load - load_forecast) / (load_forecast * delta_load 1e-6)) Gamma]; % 第二阶段成本最小化 obj_SP sum(c_buy .* max(p_buy - p_buy_ref, 0)) ... sum(c_curtail .* u_pv_curtail); optimize(Constraints, -obj_SP, sdpsettings(solver, gurobi)); u_new.p_pv value(u_pv); u_new.p_load value(u_load); obj_SP value(obj_SP); end这里的-obj_SP是因为我们要找的是最大化的最恶劣成本而YALMIP默认是最小化。3.3 几个容易出错的实现细节第一个细节是预算约束中分母的数值问题。如果预测值恰好是0比如夜间光伏预测出力为0直接做除法会出现NaN。我在实际代码里加了1e-6的平滑项这个技巧虽然看起来不太优雅但在数值计算里非常实用。也可以用max(pv_forecast * delta_pv, 0.01)这种方式效果类似。第二个细节是购售电价不一致。微电网从外部电网购电的电价通常高于向外部电网售电的电价上网电价这意味着你在目标函数里不能简单用一个净交换功率乘以一个价格必须显式区分购电变量和售电变量并且加互斥约束。如果忽略这个价格差优化器可能会出现同时购电和售电的循环套利模式——虽然这在数学上仍然满足约束但在物理世界中毫无意义而且会让储能SOC曲线出现奇怪的“充满-放空-再充满”震动。第三个细节是储能SOC的初值和终值约束。通常我们会设置SOC初值为某个水平比如0.5并且要求一天结束后SOC回到初值否则模型会倾向于把储能电量全部放光来省成本。但是如果你在多面体不确定集下同时施加“最恶劣场景下也能回到初值”的约束模型的可行域会变得非常紧甚至可能无解。我的做法是终值SOC约束只施加在预测场景下鲁棒场景下允许SOC终值偏离初值但在目标函数里加入SOC偏差惩罚项。这样既保证了调度的可持续性又不会因为过度约束导致无解。4. 常见问题与排查技巧实录4.1 收敛慢、迭代次数过多的原因两阶段鲁棒优化最常见的尴尬情况就是明明问题规模不大但CCG循环了很多轮就是不收敛上下界一直震荡。根据我的经验95%的情况出在以下几个原因。第一种情况是不确定集定义得太“宽”。比如你把光伏的波动范围设成±50%同时把负荷波动也设成±30%预算约束又给得很大这时最恶劣场景的组合空间非常大子问题每次都能找到一个新的极端场景主问题就要不断加入新变量和新约束自然收敛慢。解决方法是缩小不确定集范围或者用历史数据的5%-95%分位数来标定边界而不是拍脑袋定系数。第二种情况是子问题没有求到真正的全局最优解。如果第二阶段的LP存在数值病态Gurobi可能给出了次优解导致你找到的“最恶劣场景”并不是最恶劣的于是主问题解出来的上界比真实值低上下界差距虚高循环就停不下来。排查方法很简单——手工构造一个极端场景代入子问题对比一下目标值是否比算法给出的还大。如果答案是肯定的基本就是数值问题需要检查约束矩阵的条件数。第三种情况是目标函数中缺少必要的惩罚项。比如你没有给弃光惩罚或切负荷惩罚那么子问题在最恶劣场景下可能通过“无限弃光”来维持功率平衡这样最恶劣场景的辨识度就降低了子问题每次返回的 $\lambda$ 和 $u$ 可能都在换方向。加上合理的惩罚系数后最恶劣场景就有了明确的物理指向收敛速度会有明显改善。4.2 数值病态问题如何处理YALMIP Gurobi的组合虽然强大但数值病态问题依然防不胜防。最典型的场景是储能容量300kWh、柴发最大出力200kW、电价单位是元/kWh大概是0.5到1.2这个量级、SOC是0到1的小数。这些变量之间的量级差距可能有几百倍对于线性规划来说还可以扛一扛但到了混合整数规划求解器的容差设置就变得非常敏感。我推荐的统一做法是把所有物理量归一化到“0.1到100”的范围。具体来说功率用“MW”而不是“kW”储能容量用“MWh”价格用“元/MWh”。这样一来所有决策变量基本都在0到100这个量级Gurobi的数值容差设置就不需要额外调整。如果你始终坚持用kW和元的组合建议至少把sdpsettings(gurobi.NumericFocus, 2)打开让求解器启用更费时但更稳定的数值处理策略。另外如果约束中出现类似a * x b * y c但 $c$ 的量级在1e-6以下比如某些惩罚项阈值强烈建议先把约束整体乘以1000再写入模型。这个操作在数学上完全等价但对求解器的数值鲁棒性帮助巨大。4.3 关于代码“升级优化版本”的几点说明这个标题里提到的“升级优化版本”在我看来主要体现在三个方面。第一是求解效率的升级。旧版本里有人喜欢用枚举法处理所有不确定场景——把光伏、负荷各离散成3个水平组合出9个场景一起丢给MILP求解。这种方法在24时段下问题尚可处理但一旦扩展成96时段15分钟粒度就会直接内存爆炸。升级版的CCG算法把计算复杂度从“场景数量的指数级”降到了“迭代次数的线性级”实际工程意义非常明显。第二是模型表达的升级。旧版本往往用一个统一的确定性模型加一个“最恶劣场景”约束来近似鲁棒这其实是启发式方法严格来说不保证最优性。升级版严格对偶推导了第二阶段子问题的KKT条件保证每次迭代找到的都是真正的全局最恶劣场景。第三是代码结构上的升级。旧版本往往是几百行全部堆在一个脚本里改一个系数要找半天。升级版按模块拆分主循环、MP、SP、数据文件各自独立跑参数扫描的时候只需要在外层套一个for循环改 $\Gamma$不需要动任何模型代码。4.4 调试建议与排查速查表调试两阶段鲁棒优化程序时我建议按以下顺序排查现象可能原因检查方法解决思路一上来就无解第一阶段约束过紧多种运行状态互斥条件矛盾先去掉鲁棒场景约束只保留预测场景看能否求解检查状态变量的M值是否设置过小迭代1次就收敛但解明显不合理不确定集参数设置过小退化为确定性模型打印 $\Gamma$ 和不确定参数的实际边界对比预测值增大 $\Delta$ 或 $\Gamma$ 重新运行上界持续下降但不触发收敛子问题没有正确对偶UB计算有误人工构造一个场景手算目标值与代码输出对照重新检查子问题的对偶推导过程求解时间超过预期不确定场景数量过多或MILP分支复杂度过高查看Gurobi输出日志确认MIP Gap是否长时间停滞增加求解器MIPFocus设置或减少场景数结果曲线出现锯齿状震荡目标函数缺少平滑项或惩罚项画出SOC与购售电功率曲线观察是否频繁切换状态增加储能容量损耗成本项或启停惩罚还有一个很实用的小技巧当你不确定自己的鲁棒模型是否正确时可以把不确定集范围收敛到几乎为0$\Delta 0.001$$\Gamma 0$这时代码退化成确定性优化。拿这个结果和一个独立的确定性模型做对比如果两者输出一致说明你的鲁棒化改造没有引入模型错误。5. 一些个人实操心得做了这么多微电网调度的项目我最深的体会是两阶段鲁棒优化模型的复杂度其实不在数学推导而在工程参数的选择。不确定集设大了模型很“安全”但经济性极差储能基本成了摆设柴发永远在运行不确定集设小了模型倒是经济但一遇到极端天气就抓瞎。比较好的做法是用历史数据驱动的方式标定不确定参数波动范围比如取过去三年同期数据的90%置信区间作为上下界再用业务人员对天气风险的判断来调节 $\Gamma$ 值而不是拍脑袋给定一个数。最后再分享一个实用技巧如果你们实际的微电网项目使用的是15分钟级调度周期我建议把代码里的时段数从24扩展到96同时把不确定参数的相关性也改成相邻时段关联的形式。这种“滚动时域两阶段鲁棒”的组合才是工程落地的完全体。具体实现上只需要在不确定集约束里加入类似 $|u_{t} - u_{t-1}|_\infty \leq \eta$ 的时域平滑约束其余代码结构完全不用动。实测下来这种方式的调度计划在跟实际光伏出力曲线的贴合度上比独立时段的多面体集合要好得多成本平均能再省3%到5%。 微电网两阶段鲁棒优化：原理、Matlab实现与工程实践
1. 微电网两阶段鲁棒优化到底解决什么问题1.1 从确定性调度到不确定性调度做微电网优化调度的人最早接触的模型基本都是确定性优化光伏出力取预测值、负荷取预测值、电价取已知的分时电价然后丢给求解器一把梭得到一个“最优”的机组出力和储能充放电计划。这个做法最大的问题在于——预测永远不可能完全准确。光伏在晴天中午可能预报偏差只有5%到了多云、雷阵雨天气15分钟级滚动预测的误差轻轻松松超过20%。等到了实际运行那天你拿着昨天算好的调度计划发现光伏根本没发到预测值储能该充电的时候在放电该放电的时候又没电可放整个计划直接失效。两阶段鲁棒优化Two-Stage Robust Optimization的核心思路就是不再把不确定性参数当成一个确定的数而是给它们划一个范围——一个不确定集Uncertainty Set。所有可能出现的风光出力、负荷波动都被约束在这个集合内然后优化的目标变成在最恶劣的不确定性实现下依然能让系统运行成本最低、且所有约束都满足。换句话说它求的不是“预测场景下的最优”而是“最坏情况下的最优”。这个思路在实际工程里非常重要。微电网往往接在配电网末端系统容量小、惯性小光伏和负荷的小幅波动就可能造成电压越限或者功率倒送。你不可能指望运行人员每天盯着天气云图手动调整机组出力必须让调度策略自带“抗扰动能力”。两阶段鲁棒优化恰恰就是干这个的它能在不确定性到来之前先做好防御性的决策安排。1.2 为什么用两阶段而不是单阶段或随机规划单阶段鲁棒优化是最保守的模型——所有决策都在不确定性实现之前敲定并且要求对不确定集内的所有场景都可行。它的问题是过于悲观算出来的运行成本往往高得离谱储能基本不敢用柴发一直开在最高出力实际运行中根本没人愿意接受这种方案。随机规划是另一种思路它给不确定参数赋概率分布然后优化期望成本。这个方向理论上很美但工程落地有个尴尬问题光伏出力的概率分布你很难精确刻画尤其是短时波动往往不是标准正态分布用场景法做的话需要生成大量场景才能收敛计算量爆炸式增长而且场景缩减Scenario Reduction操作不当还可能改变问题的可行域结构。两阶段鲁棒优化把决策拆成“现在做”和“等不确定实现后再做”两步。第一阶段决策Here-and-Now是在不确定性没有揭开时就必须拍板的量比如机组启停状态、与主网的购售电协议、储能的日前充放电计划第二阶段决策Wait-and-See是在不确定性实现后可以调整的量比如具体每台机组的出力微调、储能的实时充放电功率、弃光弃风量。这样做的好处是既保留了鲁棒优化“保证所有场景可行”的安全性又不像单阶段那样僵硬——第二阶段可以根据实际情况灵活调整成本自然不会那么离谱。1.3 一个直观的类比打个比方这就像你计划明天出门旅行。确定性优化相当于你完全相信天气预报说明天不下雨于是只带了一顶遮阳帽单阶段鲁棒优化相当于你假设明天可能下暴雨、刮台风、下冰雹甚至地震于是把整个家都搬上车后备箱都盖不上两阶段鲁棒优化则是你先带上雨伞、冲锋衣这种能应付大多数恶劣天气的装备第一阶段决策等第二天早上真正看到天气情况后再决定穿哪件衣服、走哪条路第二阶段决策——既保证了不会淋雨又不会出门带太多累赘。2. 模型构建的核心逻辑2.1 第一阶段与第二阶段的决策变量划分建模之前最需要想清楚的问题就是哪些变量放第一阶段哪些变量放第二阶段。这个划分直接决定了模型的保守程度和求解难度。在我做的这个微电网系统里典型配置是光伏阵列、风力发电机、蓄电池储能、柴油发电机或者微型燃气轮机、与外部电网的联络线以及本地负荷。第一阶段决策变量选的是柴发的启停状态0-1整数变量、储能是否处于充电/放电状态的二进制变量、从外部电网购电/售电的状态变量以及与主网的交换功率计划值。这些变量对应的是日前时间尺度需要在不确定性揭晓之前就确定下来。第二阶段决策变量则是每台机组的实际出力调整量、储能的实际充电/放电功率、弃光弃风率、负荷削减量、与主网实际交换功率的偏差修正值。这些量可以在实时运行中根据实际光伏出力和负荷情况做调整。需要特别提醒的是储能是否充放电这个0-1变量最好不要放在第一阶段除非你有特殊需求。原因是这会让模型变成混合整数两阶段鲁棒优化MIRU求解复杂度直接上一个台阶。很多时候可以采用一个简化技巧允许储能在第二阶段连续调整充放电功率但通过一个足够大的M值约束防止它同时充放电这样主问题Master Problem依然是个MILP但子问题Subproblem变成了纯LP处理起来省事得多。2.2 不确定集的设计与参数调节不确定集的选择是两阶段鲁棒优化建模的灵魂。常用的有盒式不确定集Box、椭球不确定集Ellipsoidal和多面体不确定集Polyhedral。盒式不确定集形式最简单每个不确定参数独立波动在其上下界之间即 $u \in [\hat{u} - \hat{u}\Delta, \hat{u} \hat{u}\Delta]$。但是这种集合的保守度很高因为它允许所有不确定参数同时取到最坏值而现实中这种情况概率极低。椭球不确定集考虑了参数之间的相关性形状更贴近实际但问题会变成二阶锥规划SOCP甚至更复杂的形式工程实现门槛高。多面体不确定集也叫预算不确定集Budget Uncertainty Set在盒式集合的基础上增加一个总偏差预算约束例如 $\sum_i |u_i - \hat{u}_i| / \hat{u}_i \leq \Gamma$这个 $\Gamma$ 就是鲁棒预算Budget of Uncertainty。我在代码里采用的正是多面体不确定集因为它在工程实用性和保守度之间取得了很好的平衡。$\Gamma$ 的物理含义是“最多同时允许多少个不确定参数偏离预测值”取0时就是确定性模型取值越大就越保守。实际操作中我习惯先跑一次确定性模型得到基准成本然后逐步增大 $\Gamma$观察总成本的变化曲线。如果成本增长斜率在某个点之后突然变陡那说明再增加保守度就不划算了——这个拐点附近就是 $\Gamma$ 的最佳取值区间。这里分享一个我踩过的坑盒式集合在某几个子问题上解出来的目标值偏大一开始我还以为是求解器数值问题后来排查了半天才发现是不确定参数之间的相关性没有建模导致最坏场景同时乐观估计了充电电价又悲观估计了放电电价双重叠加把成本推高了。改成多面体集合加上预算约束后结果就正常了。2.3 目标函数与约束的矩阵化表达微电网经济调度的目标函数一般包含以下几项柴油发电机燃料成本、启停成本、从主网购电成本减去售电收益、储能充放电循环老化成本可选以及弃光弃风惩罚项可选如果允许弃电的话。两阶段鲁棒优化的标准形式是$$\min_{x} \ c^T x \max_{u \in \mathcal{U}} \min_{y \in \Omega(x, u)} d^T y$$其中 $x$ 是第一阶段变量$y$ 是第二阶段变量$u$ 是不确定参数。内层 $\max \min$ 表示在给定第一阶段决策后寻找最恶劣的不确定实现然后在最恶劣情况下最小化运行成本。对应地约束分为几个层次第一阶段固有约束机组启停逻辑约束、最小启停时间约束、储能充放电状态互斥约束、购售电状态互斥约束。耦合约束功率平衡约束中既含第一阶段变量也含第二阶段变量还有储能SOC的状态转移方程。第二阶段约束机组出力上下限、爬坡约束、储能充放电功率限制、SOC上下限、与主网交换功率限制、弃电量约束。写Matlab代码的时候我强烈建议把所有约束整理成矩阵形式而不是标量循环。原因有两个一是YALMIP对矩阵化约束的预处理效率远高于逐条添加二是调试的时候你能直接看约束矩阵的维度是否匹配问题定位快得多。比如功率平衡约束如果你有24个时段的优化直接构造一个24行的等式约束矩阵而不是写24次constraint [constraint; ...]。3. Matlab代码的实现思路与关键细节3.1 整体代码结构与主循环框架拿到一个两阶段鲁棒优化问题最常用的求解方法是列与约束生成算法Column-and-Constraint GenerationCCG也常称为Benders分解的变体。核心思想是把原问题分解成一个主问题MP和一个子问题SP迭代求解。主问题是在已知的有限个最恶劣场景下做优化得到一个下界对最小化问题而言子问题是在给定主问题解的基础上寻找新的最恶劣不确定场景并给主问题上界。整个代码框架我分成了几个模块main.m主程序负责设置系统参数、调用求解循环、输出结果。data_define.m所有基础数据的定义——负荷曲线、光伏预测出力、电价曲线、储能参数、柴发参数、不确定集参数。build_MP.m构建主问题的YALMIP模型。solve_SP.m求解子问题找出当前迭代的最恶劣场景。plot_results.m将迭代过程中的调度结果可视化。CCG主循环的逻辑如下% 初始化 LB -inf; UB inf; iter 1; max_iter 10; % 或者用收敛判据 best_scenarios {}; % 存储所有已发现的最恶劣场景 while (UB - LB) / abs(UB) 1e-4 iter max_iter % 1. 求解主问题加入已知的最恶劣场景 [x_opt, obj_MP] solve_MP(best_scenarios); LB obj_MP; % 2. 代入主问题解求解子问题 [u_new, obj_SP, y_opt] solve_SP(x_opt); UB min(UB, obj_MP - obj_SP); % 注意这里根据对偶形式调整 % 3. 如果子问题目标值对应的场景尚未覆盖则加入场景并继续 if abs(obj_SP) 1e-6 best_scenarios{end1} u_new; end % 4. 检查收敛 fprintf(Iteration %d: LB%.4f, UB%.4f, gap%.4f%%\n, iter, LB, UB, ... abs(UB - LB) / abs(UB) * 100); iter iter 1; end有一点需要注意在标准的两阶段鲁棒优化中子问题通常是“给定第一阶段决策 $x^*$寻找最大化的最恶劣场景下的最小运行成本”。如果子问题是LP可以通过KKT条件或者对偶变换把它转成单层优化如果子问题含有整数变量就得用分解算法或者大M法线性化。在我这个模型中第二阶段只含连续变量所以子问题通过强对偶转成单层的max问题然后用YALMIPGurobi一步求出最恶劣场景和对应的成本。3.2 YALMIP建模与求解器设置代码实现层面我用的是 YALMIP 作为建模语言求解器选的是 Gurobi实测比CPLEX在MIQP上更快也可以用SCIP或者CBC做开源替代。一个完整的MP构建代码示意如下function [Constraints, Obj, x_vars] build_MP(best_scenarios) % 第一阶段变量 x_vars.u_on binvar(24, 1); % 柴发启停 x_vars.x_charge binvar(24, 1); % 储能充电状态 x_vars.x_discharge binvar(24, 1); % 储能放电状态 x_vars.x_buy binvar(24, 1); % 购电状态 x_vars.x_sell binvar(24, 1); % 售电状态 x_vars.p_buy sdpvar(24, 1); % 购电功率 x_vars.p_sell sdpvar(24, 1); % 售电功率 x_vars.soc sdpvar(24, 1); % 储能荷电状态 x_vars.p_sto_out sdpvar(24, 1); % 储能放电功率 x_vars.p_sto_in sdpvar(24, 1); % 储能充电功率 x_vars.p_dg sdpvar(24, 1); % 柴发出力 Constraints []; % 储能状态互斥 Constraints [Constraints, x_vars.x_charge x_vars.x_discharge 1]; Constraints [Constraints, x_vars.x_buy x_vars.x_sell 1]; % SOC递推 for t 2:24 Constraints [Constraints, ... x_vars.soc(t) x_vars.soc(t-1) ... x_vars.p_sto_in(t) * eta_ch / Cap - ... x_vars.p_sto_out(t) / eta_dis / Cap]; end % 购售电功率上限大M形式 Constraints [Constraints, x_vars.p_buy x_vars.x_buy * P_grid_max]; Constraints [Constraints, x_vars.p_sell x_vars.x_sell * P_grid_max]; % 柴发出力范围 Constraints [Constraints, x_vars.p_dg x_vars.u_on * P_dg_max]; Constraints [Constraints, x_vars.p_dg x_vars.u_on * P_dg_min]; % 对每一个已知的最恶劣场景建立第二阶段约束 Obj sum(fuel_cost(x_vars.p_dg)) sum(start_cost .* max(0, diff([0; x_vars.u_on]))) ... sum(price_buy .* x_vars.p_buy) - sum(price_sell .* x_vars.p_sell); for k 1:length(best_scenarios) u_k best_scenarios{k}; % 第k个最恶劣场景 % 对应场景下的第二阶段变量 y_pv_curtail sdpvar(24, 1); y_load_curtail sdpvar(24, 1); y_dg_adjust sdpvar(24, 1); % 功率平衡约束 Constraints [Constraints, ... u_k.p_pv u_k.p_wind x_vars.p_dg y_dg_adjust ... x_vars.p_sto_out - x_vars.p_sto_in x_vars.p_buy - x_vars.p_sell ... y_load_curtail u_k.p_load - y_pv_curtail]; % 弃光限制 Constraints [Constraints, 0 y_pv_curtail u_k.p_pv]; Constraints [Constraints, y_load_curtail 0]; end end这段代码只展示了核心思路实际项目中还需要考虑柴发最小启停时间、储能SOC初值设定等细节。但基本框架就是上面这样——主问题每轮迭代都会重新构建一次把新发现的场景追加到约束里。子问题的求解稍微复杂一点。给定第一阶段变量 $x^*$对偶问题可以写成function [u_new, obj_SP] solve_SP(x_vars) % 不确定变量 u_pv sdpvar(24, 1); u_load sdpvar(24, 1); % 多面体不确定集约束 Constraints [Constraints, ... u_pv pv_forecast * (1 - delta_pv), ... u_pv pv_forecast * (1 delta_pv), ... u_load load_forecast * (1 - delta_load), ... u_load load_forecast * (1 delta_load)]; % 预算约束 Constraints [Constraints, ... sum(abs(u_pv - pv_forecast) / (pv_forecast * delta_pv 1e-6)) ... sum(abs(u_load - load_forecast) / (load_forecast * delta_load 1e-6)) Gamma]; % 第二阶段成本最小化 obj_SP sum(c_buy .* max(p_buy - p_buy_ref, 0)) ... sum(c_curtail .* u_pv_curtail); optimize(Constraints, -obj_SP, sdpsettings(solver, gurobi)); u_new.p_pv value(u_pv); u_new.p_load value(u_load); obj_SP value(obj_SP); end这里的-obj_SP是因为我们要找的是最大化的最恶劣成本而YALMIP默认是最小化。3.3 几个容易出错的实现细节第一个细节是预算约束中分母的数值问题。如果预测值恰好是0比如夜间光伏预测出力为0直接做除法会出现NaN。我在实际代码里加了1e-6的平滑项这个技巧虽然看起来不太优雅但在数值计算里非常实用。也可以用max(pv_forecast * delta_pv, 0.01)这种方式效果类似。第二个细节是购售电价不一致。微电网从外部电网购电的电价通常高于向外部电网售电的电价上网电价这意味着你在目标函数里不能简单用一个净交换功率乘以一个价格必须显式区分购电变量和售电变量并且加互斥约束。如果忽略这个价格差优化器可能会出现同时购电和售电的循环套利模式——虽然这在数学上仍然满足约束但在物理世界中毫无意义而且会让储能SOC曲线出现奇怪的“充满-放空-再充满”震动。第三个细节是储能SOC的初值和终值约束。通常我们会设置SOC初值为某个水平比如0.5并且要求一天结束后SOC回到初值否则模型会倾向于把储能电量全部放光来省成本。但是如果你在多面体不确定集下同时施加“最恶劣场景下也能回到初值”的约束模型的可行域会变得非常紧甚至可能无解。我的做法是终值SOC约束只施加在预测场景下鲁棒场景下允许SOC终值偏离初值但在目标函数里加入SOC偏差惩罚项。这样既保证了调度的可持续性又不会因为过度约束导致无解。4. 常见问题与排查技巧实录4.1 收敛慢、迭代次数过多的原因两阶段鲁棒优化最常见的尴尬情况就是明明问题规模不大但CCG循环了很多轮就是不收敛上下界一直震荡。根据我的经验95%的情况出在以下几个原因。第一种情况是不确定集定义得太“宽”。比如你把光伏的波动范围设成±50%同时把负荷波动也设成±30%预算约束又给得很大这时最恶劣场景的组合空间非常大子问题每次都能找到一个新的极端场景主问题就要不断加入新变量和新约束自然收敛慢。解决方法是缩小不确定集范围或者用历史数据的5%-95%分位数来标定边界而不是拍脑袋定系数。第二种情况是子问题没有求到真正的全局最优解。如果第二阶段的LP存在数值病态Gurobi可能给出了次优解导致你找到的“最恶劣场景”并不是最恶劣的于是主问题解出来的上界比真实值低上下界差距虚高循环就停不下来。排查方法很简单——手工构造一个极端场景代入子问题对比一下目标值是否比算法给出的还大。如果答案是肯定的基本就是数值问题需要检查约束矩阵的条件数。第三种情况是目标函数中缺少必要的惩罚项。比如你没有给弃光惩罚或切负荷惩罚那么子问题在最恶劣场景下可能通过“无限弃光”来维持功率平衡这样最恶劣场景的辨识度就降低了子问题每次返回的 $\lambda$ 和 $u$ 可能都在换方向。加上合理的惩罚系数后最恶劣场景就有了明确的物理指向收敛速度会有明显改善。4.2 数值病态问题如何处理YALMIP Gurobi的组合虽然强大但数值病态问题依然防不胜防。最典型的场景是储能容量300kWh、柴发最大出力200kW、电价单位是元/kWh大概是0.5到1.2这个量级、SOC是0到1的小数。这些变量之间的量级差距可能有几百倍对于线性规划来说还可以扛一扛但到了混合整数规划求解器的容差设置就变得非常敏感。我推荐的统一做法是把所有物理量归一化到“0.1到100”的范围。具体来说功率用“MW”而不是“kW”储能容量用“MWh”价格用“元/MWh”。这样一来所有决策变量基本都在0到100这个量级Gurobi的数值容差设置就不需要额外调整。如果你始终坚持用kW和元的组合建议至少把sdpsettings(gurobi.NumericFocus, 2)打开让求解器启用更费时但更稳定的数值处理策略。另外如果约束中出现类似a * x b * y c但 $c$ 的量级在1e-6以下比如某些惩罚项阈值强烈建议先把约束整体乘以1000再写入模型。这个操作在数学上完全等价但对求解器的数值鲁棒性帮助巨大。4.3 关于代码“升级优化版本”的几点说明这个标题里提到的“升级优化版本”在我看来主要体现在三个方面。第一是求解效率的升级。旧版本里有人喜欢用枚举法处理所有不确定场景——把光伏、负荷各离散成3个水平组合出9个场景一起丢给MILP求解。这种方法在24时段下问题尚可处理但一旦扩展成96时段15分钟粒度就会直接内存爆炸。升级版的CCG算法把计算复杂度从“场景数量的指数级”降到了“迭代次数的线性级”实际工程意义非常明显。第二是模型表达的升级。旧版本往往用一个统一的确定性模型加一个“最恶劣场景”约束来近似鲁棒这其实是启发式方法严格来说不保证最优性。升级版严格对偶推导了第二阶段子问题的KKT条件保证每次迭代找到的都是真正的全局最恶劣场景。第三是代码结构上的升级。旧版本往往是几百行全部堆在一个脚本里改一个系数要找半天。升级版按模块拆分主循环、MP、SP、数据文件各自独立跑参数扫描的时候只需要在外层套一个for循环改 $\Gamma$不需要动任何模型代码。4.4 调试建议与排查速查表调试两阶段鲁棒优化程序时我建议按以下顺序排查现象可能原因检查方法解决思路一上来就无解第一阶段约束过紧多种运行状态互斥条件矛盾先去掉鲁棒场景约束只保留预测场景看能否求解检查状态变量的M值是否设置过小迭代1次就收敛但解明显不合理不确定集参数设置过小退化为确定性模型打印 $\Gamma$ 和不确定参数的实际边界对比预测值增大 $\Delta$ 或 $\Gamma$ 重新运行上界持续下降但不触发收敛子问题没有正确对偶UB计算有误人工构造一个场景手算目标值与代码输出对照重新检查子问题的对偶推导过程求解时间超过预期不确定场景数量过多或MILP分支复杂度过高查看Gurobi输出日志确认MIP Gap是否长时间停滞增加求解器MIPFocus设置或减少场景数结果曲线出现锯齿状震荡目标函数缺少平滑项或惩罚项画出SOC与购售电功率曲线观察是否频繁切换状态增加储能容量损耗成本项或启停惩罚还有一个很实用的小技巧当你不确定自己的鲁棒模型是否正确时可以把不确定集范围收敛到几乎为0$\Delta 0.001$$\Gamma 0$这时代码退化成确定性优化。拿这个结果和一个独立的确定性模型做对比如果两者输出一致说明你的鲁棒化改造没有引入模型错误。5. 一些个人实操心得做了这么多微电网调度的项目我最深的体会是两阶段鲁棒优化模型的复杂度其实不在数学推导而在工程参数的选择。不确定集设大了模型很“安全”但经济性极差储能基本成了摆设柴发永远在运行不确定集设小了模型倒是经济但一遇到极端天气就抓瞎。比较好的做法是用历史数据驱动的方式标定不确定参数波动范围比如取过去三年同期数据的90%置信区间作为上下界再用业务人员对天气风险的判断来调节 $\Gamma$ 值而不是拍脑袋给定一个数。最后再分享一个实用技巧如果你们实际的微电网项目使用的是15分钟级调度周期我建议把代码里的时段数从24扩展到96同时把不确定参数的相关性也改成相邻时段关联的形式。这种“滚动时域两阶段鲁棒”的组合才是工程落地的完全体。具体实现上只需要在不确定集约束里加入类似 $|u_{t} - u_{t-1}|_\infty \leq \eta$ 的时域平滑约束其余代码结构完全不用动。实测下来这种方式的调度计划在跟实际光伏出力曲线的贴合度上比独立时段的多面体集合要好得多成本平均能再省3%到5%。