尧图网站建设 
1. 为什么微网调度要从确定性走向两阶段鲁棒1.1 一次光伏骤降让确定性调度直接翻车去年年底我在做一个园区级微网的日前经济调度项目碰到的第一个问题就非常典型光伏和风电的预测曲线在早上看起来非常漂亮预测中午光伏能到120kW结果上午十点一片云飘过来光伏出力从75kW直接掉到35kW。当时用的是最传统的确定性调度所有风光出力都按预测值当作已知参数喂给优化模型燃气轮机出力和储能充放电曲线也是按这个预测排好的。实际运行中光伏一掉燃气轮机因为爬坡约束没法立刻补上储能又已经在低谷时段放得差不多了最后只能从主网买高价电那一天的综合运行成本比调度计划超支了12%还要多。这个经历让我意识到微网优化调度如果只盯着预测值做经济调度本质上是在赌天气。后来我把求解框架换成了两阶段鲁棒优化配合CCG列与约束生成算法在MATLABYalmip环境里重新搭建了整个日前经济调度模型。这不是一个简单的多套几个场景的思路而是整个优化问题的结构都变了把决策分成今天就要定下来的事和等不确定性揭晓后再调整的事两层然后专门针对最坏的不确定性做应对。1.2 确定性优化、随机优化、鲁棒优化三种思路的取舍在做微网调度之前得先分清三种主流建模思路的适用边界不然很容易选错方向。第一种是确定性优化模型长这样[ \min_{x} ; c^Tx \quad s.t. ; Ax \ge b ]这里的负荷、光伏、风电全用预测值。优点是模型简单LP或者MILP一求解就完事速度快。缺点也明摆着预测和实际的偏差直接变成电网的安全风险或者成本风险。实际工程中靠人工预留旋转备用去对冲误差本质上还是经验主义。第二种是随机规划Stochastic Programming思路是把不确定性描述成若干离散场景每个场景有一个概率目标函数变成期望成本最小。比如用蒙特卡洛采样生成1000个光伏出力场景然后求期望。随机规划的强项是有概率意义理论上很优雅但它的痛点也很实际场景多了模型爆炸场景少了概率分布又不准。而且随机规划对分布假设敏感如果你把分布形状搞错了解出来的期望最优可能还不如确定性调度。第三种就是鲁棒优化Robust Optimization它不问概率分布只给一个不确定集合把问题写成[ \min_{x} \max_{u \in U} ; c^Tx f(x,u) ]要求的是在集合内任意一个不确定性实现下方案都可行且代价可控。鲁棒优化不追求期望最优它追求的是最坏情况下的保证。代价是方案必然偏保守但换来的是安全性。两阶段鲁棒优化是鲁棒优化在决策有先后场景下的升级版一部分决策必须在不确性实现之前敲定另一部分决策可以等看到实际风光出力之后再调整。1.3 两阶段鲁棒在微网场景下的适配性为什么微网经济调度特别适合两阶段鲁棒因为微网的运行逻辑本身就是分层的。日前阶段燃气轮机要不要开、开几台、是否和主网签订购电协议这些涉及整数变量的决策必须提前确定执行层面没法等到实时再改。而在实时阶段光伏实际出力出来了燃气轮机出力、储能充放电、主网交换功率这些连续量是可以快速调整的。两个阶段决策的自由度完全不同恰好对应两阶段鲁棒优化的结构。用一句话概括这个思路第一阶段花钱买选项第二阶段根据最坏情况动用选项把成本压住。这套东西做出来之后我最大的感受是两阶段鲁棒不是让你去准确预测不确定性而是让你在不确定性真正发生的时候有一个事先算好的、能兜底的调整策略。对微网这种源荷双重波动、设备响应速度又各不相同的系统来说这个特性太重要了。2. 两阶段鲁棒经济调度的数学模型变量、目标与不确定集2.1 微网设备、变量划分与运行假设先把我搭建的微网系统讲清楚。为了兼顾论文级严谨和代码可实现我选了包含微型燃气轮机、储能、光伏、风电、主网连接和本地负荷的经典拓扑。微型燃气轮机MT最小出力20kW最大出力80kW爬坡约束30kW/h燃料成本用二次函数近似储能系统ESS容量200kWh最大充放电功率50kW充放电效率0.95SOC限制在0.1~0.9初始SOC为0.5光伏风电不可控出力落在不确定集合内主网交互允许购电和售电购售电价分时且有交互功率上限负荷本文先按确定值处理实际工程中可以方便地把负荷波动也塞进不确定集合。变量划分是两阶段鲁棒模型的核心。第一阶段决策变量是所有必须在不确定性实现前定下来的整数量和部分合同量燃气轮机启停状态 (y_{on,t})与主网的购/售电状态 (y_{grid,t})整数储能是否处于充电状态 (y_{ch,t})防止同时充放电的整数约束。第二阶段决策变量则是在看到实际风光出力后可以调整的连续量燃气轮机出力 (P_{MT,t})储能充放电功率 (P_{ch,t}, P_{dis,t})主网购/售电功率 (P_{buy,t}, P_{sell,t})弃光/弃风量 (P_{curt,t})、切负荷量 (P_{cut,t})。这里有两条运行假设需要说明。第一储能同时充放电约束是用整数变量表达的在日前调度里基本是标配不然会出现充电又放电这种物理上不成立但数学上省钱的操作。第二我把切负荷和弃风弃光放进了第二阶段并赋予很高的惩罚系数目的是保证模型在极端场景下有退路同时也让经济性指标更有解释力。2.2 目标函数与约束的完整建模两阶段鲁棒优化目标的标准写法是[ \min_{y} ; c^T y \max_{u \in U} ; \min_{z \in F(y,u)} ; d^T z ]其中 (y) 是第一阶段决策(z) 是第二阶段调整变量(u) 是不确定变量。第一阶段的成本 (c^Ty) 包括燃气轮机启停成本我按单台单次启停300元计。第二阶段成本 (d^Tz) 则包含燃气轮机燃料成本用分段线性化逼近二次项向主网购电成本减去售电收入弃风弃光惩罚和切负荷惩罚后者单价设为1000元/MWh远高于正常发电成本。功率平衡约束是模型的主心骨对每个时段 (t)[ P_{MT,t} P_{dis,t} - P_{ch,t} P_{PV,t} P_{WT,t} P_{buy,t} - P_{sell,t} P_{cut,t} P_{load,t} P_{curt,t} ]这里我把弃风弃光量放在等式右侧表示多余电量可以被主动放弃切负荷量放在左侧表示需求不足时允许切除部分负荷。燃气轮机约束包括出力上下限、爬坡约束[ P_{MT}^{min} y_{on,t} \le P_{MT,t} \le P_{MT}^{max} y_{on,t} ][ -P_{ramp} \le P_{MT,t} - P_{MT,t-1} \le P_{ramp} ]储能约束包括功率限制、SOC递推、SOC上下限[ SOC_{t1} SOC_t \frac{\eta_{ch} P_{ch,t} - P_{dis,t}/\eta_{dis}}{E_{ESS}} ][ 0 \le P_{ch,t} \le P_{ch}^{max} y_{ch,t}, \quad 0 \le P_{dis,t} \le P_{dis}^{max} (1-y_{ch,t}) ]主网交互约束比较简单购售电功率各有一个上限且通过整数变量保证不会同时买卖[ 0 \le P_{buy,t} \le P_{grid}^{max} y_{grid,t}, \quad 0 \le P_{sell,t} \le P_{grid}^{max} (1-y_{grid,t}) ]模型到这里还是比较标准的。从代码实现角度看这些约束绝大多数只涉及第二阶段变量只有启停状态、购售电状态这些整数变量来自第一阶段这正是CCG算法适用的典型结构。2.3 不确定集设计盒式集合与预算约束两阶段鲁棒优化的鲁棒程度完全由不确定集合 (U) 决定。集合越大方案越保守成本越高集合越小方案越激进风险越大。这一步如果处理不好后面所有计算都是白搭。我最开始用的是最朴素的盒式不确定集[ U \left{ u ;|; u_j^{min} \le u_j \le u_j^{max}, \forall j \right} ]每个不确定量独立地在区间内波动。光伏预测120kW允许±15%那光伏出力就在102kW到138kW之间波动风电预测60kW允许±25%那就在45kW到75kW之间波动。但盒式集合的麻烦在于它会让模型默认所有不确定量同时取最坏值这是极端保守的。真实场景中光伏和风电同时达到最坏边界的概率极低用盒式集合做出来的调度成本会明显偏高。因此我在模型中引入了预算约束Budget Constraint[ \sum_j \frac{|u_j - u_j^{pred}|}{\Delta u_j} \le \Gamma ]这个 (\Gamma) 就是总偏差预算它限制的是所有不确定量偏离预测值的总幅度。比如24个时段(\Gamma4) 意味着最多等效于4个时段的不确定变量同时取到边界其他时段只能在预测值附近小幅波动。(\Gamma) 越大越保守(\Gamma0) 等价于确定性调度(\Gamma24) 等价于盒式集合。预算约束的好处非常多。第一它让不确定集合的体积可调给了调度员一个直观的经济鲁棒性权衡旋钮第二它在数学上保持了线性结构配合CCG迭代可以采用MILP求解第三实践中最坏场景一定会落在边界点组合上这让后面的对偶和线性化变得有迹可循。3. CCG算法拆解主问题-子问题迭代如何求解min-max-min3.1 为什么Benders不够CCG更合适第一次看到两阶段鲁棒模型的时候我的第一反应是能不能用Benders分解毕竟Benders分解在随机规划和确定性混合整数规划里是很成熟的技术。后来我仔细对比之后发现CCG在处理两阶段鲁棒优化上明显更顺手。Benders分解的思路是主问题求解第一阶段决策子问题求解第二阶段对偶然后把子问题的最优值函数用一条对偶割最优割或可行性割加回主问题迭代逼近原问题。问题在于子问题返回的信息是一个割平面本质上是目标值的下界逼近如果第二阶段存在整数变量Benders切割的处理会非常麻烦实际实现中经常出现振荡或收敛慢。CCG的核心区别是子问题每迭代一次生成一个具体的最坏场景 (u^*)并把该场景对应的第二阶段变量和全部约束原封不动地加入主问题。这相当于在主问题里新增一列和一组约束所以叫列与约束生成。因为每一轮都引入了完整场景的信息而不是一个压缩的割平面CCG的收敛速度通常远快于Benders。对于微网调度这种规模不太大、但场景约束很多的模型CCG往往在3~5次迭代内就能给出满意的间隙。3.2 主问题MP的构造与子问题SP的对偶化CCG迭代的第一步是把原始的两阶段模型拆成主问题MP和子问题SP。主问题MP在每一轮迭代时包含当前已发现的全部最坏场景。假设已经迭代到第 (k) 轮有 (k) 个已知场景 (u_1^, u_2^, ..., u_k^*)MP可以写为[ \min_{y, z_k, \eta} ; c^T y \eta ][ s.t. ; \eta \ge d^T z_k, \quad \forall k ][ z_k \in F(y, u_k^*), \quad \forall k ][ y \in Y, ; z_k \in Z ]辅助变量 (\eta) 的作用是逼近第二阶段的最坏成本下界。随着场景增多MP被约束逐步收紧目标函数值从下往上逼近真实最优值这就是下界LB。子问题SP是在给定第一阶段解 (y^*) 后找最坏场景和对应的最小第二阶段成本[ Q(y^*) \max_{u \in U} ; \min_{z} ; d^T z ][ s.t. ; Gz \ge h - Ey^* - Mu ]内层是关于连续变量 (z) 的LP因此可以用强对偶定理把它转成单层最大化问题。设 (\pi) 为内层约束对应的对偶变量SP等价为[ \max_{u, \pi} ; \pi^T (h - Ey^* - Mu) ][ s.t. ; G^T \pi \le d, \quad \pi \ge 0, \quad u \in U ]这里会出现一个双线性项 (\pi^T M u)这是整个CCG实现中第一个真正的技术难点。如果 (M) 矩阵的每个不确定变量只出现在一个约束里比如光伏只出现在功率平衡约束那么这项可以写成 (u)的线性组合系数由 (\pi) 决定的形态。当不确定集合是盒式时最优 (u) 一定取边界方向由 (\pi^T M) 的符号决定当加上预算约束时需要引入0-1变量选择哪些变量取边界SP就变成一个MILP。3.3 迭代流程、上下界更新与收敛条件CCG的整体流程我总结成下面这个循环这也是我在MATLAB里实际实现的骨架初始化给定初始场景 (u_1)我习惯取预测值设置 (LB -\infty)(UB \infty)迭代次数 (k1)求解MP得到最优解 (y^, \eta^)更新 (LB c^Ty^* \eta^*)固定 (y^)求解SP得到 (Q(y^)) 和新的最坏场景 (u_{k1})更新 (UB \min(UB, c^Ty^* Q(y^*)))计算相对间隙 (\frac{UB-LB}{|UB|})如果小于设定容差比如1e-3停止迭代否则把场景 (u_{k1}) 对应的第二阶段变量 (z_{k1}) 和约束 (F(y, z_{k1}, u_{k1})) 加入MP(k k1)回到第2步。这里有个细节值得注意MP目标值是松弛解所以LB一定小于等于原问题最优值SP给出的 (c^Ty^* Q(y^*)) 是某个可行方案的完整成本所以UB一定大于等于原问题最优值。两个界从两边夹逼收敛判据是有严格意义的。我在第一次跑通时看到LB和UB逐步靠近那种感觉很直观比Benders那种割平面越来越多但界还是慢吞吞动的体验好太多。4. MATLABYalmip环境下的CCG实现全流程4.1 环境准备Yalmip与求解器配置MATLAB这边我用的版本是R2023b建模工具是Yalmip求解器用Gurobi 10.x。这里多啰嗦一句两阶段鲁棒模型的内层和外层都涉及MILPYalmip自带的默认求解器是做不了大规模MILP的必须配置好商业求解器。Gurobi和CPLEX都行我个人偏Gurobi原因后面踩坑部分会说。Yalmip的安装很简单把文件夹加进MATLAB路径就行然后用yalmiptest验证一下求解器是否被识别。求解器配置代码ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 1e-4; ops.gurobi.TimeLimit 600;我的习惯是外层CCG的迭代容差设1e-3内层Gurobi的MIPGap设1e-4。内层比外层紧一个量级可以避免外层还没收敛、内层精度不够导致的UB抖动。4.2 主问题MP的代码框架在Yalmip里建MP核心是按场景循环累加约束。因为MP会随着迭代不断长大我建议一开始就用cell数组或者多维sdpvar把第二阶段变量组织好避免在循环里反复创建变量。我采用的变量组织方式是T 24; % 时段数 nMT 1; % 燃机台数 nCont 4; % 第二阶段连续变量个数P_MT, P_ch, P_dis, P_buy, P_sell... K 0; % 当前已入库场景数 y_on binvar(nMT, T, full); % 燃机启停第一阶段 y_grid binvar(1, T, full); % 主网购售电状态第一阶段 z cell(K, 1); % 每个已发现场景对应的第二阶段变量 eta sdpvar(1);主问题约束的构建遵循一个套路先写第一阶段公共约束再对每个已发现场景 (k)把第二阶段约束完整写一遍。功率平衡约束的代码示例C []; for k 1:K for t 1:T u_PV u_set{k}(t, 1); u_WT u_set{k}(t, 2); C [C, z{k}(1,t) z{k}(2,t) - z{k}(3,t) ... % P_MT P_dis - P_ch z{k}(4,t) - z{k}(5,t) u_PV u_WT - z{k}(6,t) P_load(t)]; % 其中 z(k,:) [P_MT, P_dis, P_ch, P_buy, P_sell, P_cut] end end目标函数里辅助变量(\eta)对应所有场景第二阶段成本的最大值下界所以MP目标就是启停成本加上(\eta)obj sum(sum(c_on .* y_on)) 1e-6 * sum(y_on(:)) eta;这里1e-6的微调是为了防止整数变量出现不希望的退化解属于经验性的小技巧。MP是MILP第一轮没有场景时可以用一个空场景或者直接设MP只有第一步触发实际操作中我会先手动把初始场景(u_1)加进MP保证模型非退化。4.3 子问题SP与双线性项线性化代码SP的实现是CCG里最考验功底的部分。我先把一般形式列出来再给代码。SP内层LP写成标准形式[ \min_z d^T z \quad s.t. ; Gz \ge h - Ey^* - Mu ]对偶后是[ \max_{\pi \ge 0} ; \pi^T(h - Ey^* - Mu) \quad s.t. ; G^T\pi \le d ]在Yalmip里直接建模这个对偶问题pi_var sdpvar(n_cons, 1); % 对偶变量 v binvar(n_u, 1, full); % 不确定变量是否取上界的0-1指示 u u_lb (u_ub - u_lb) .* v; % 不确定量 C [G * pi_var d, pi_var 0]; C [C, sum(v) Gamma]; % 目标中的双线性项处理需要先定义 w(i,j) pi_var(i) * v(j) M_big 1000; w sdpvar(n_cons, n_u, full); for i 1:n_cons for j 1:n_u C [C, 0 w(i,j) M_big * v(j)]; C [C, pi_var(i) - M_big*(1-v(j)) w(i,j) pi_var(i)]; end end obj_sp pi_var * (h - E * y_star) - sum(sum(M .* w));这里的目标函数用了一个数学等价(\pi^T M u \sum_i \sum_j \pi_i M_{ij} u_j)把 (u_j) 展开成 (u_{lb,j} (u_{ub,j} - u_{lb,j})v_j) 后常数项直接保留变量部分用大M线性化结果(w_{i,j})来替代 (\pi_i v_j)。这样SP就是一个标准MILPGurobi可以直接求解。不过说实话大M法线性化不是唯一选择也未必是最好选择。我在迭代了几轮后发现如果模型规模不大还有一种更省事的做法因为(\pi)的最优值在LP极点处取得而(u_j)在最优时必取边界可以先用符号规则求解一个线性目标再通过分支定界选择预算。这个技巧比较复杂新手阶段建议还是稳妥地走大M法先把流程跑通。4.4 CCG主循环与LB/UB收敛过程输出主循环代码写出来非常短但逻辑要非常小心。伪代码如下LB -1e6; UB 1e6; K 1; u_set{K} u_pred; % 初始场景取预测值 gap 1; while gap 1e-3 K 10 % 1. 建立并求解MP [MP, obj_MP] build_MP(u_set, K, ...); optimize(MP, obj_MP, ops); y_star value(y_on); LB max(LB, value(obj_MP)); % 2. 固定y_star求解SP [SP, obj_SP_val, u_new] solve_SP(y_star, ...); UB min(UB, value(obj_SP_val)); % 3. 计算gap gap abs((UB - LB) / abs(UB)); % 4. 入库新场景 K K 1; u_set{K} u_new; end每轮迭代我都会记录LB、UB、gap、MP求解时间、SP求解时间最后用plot画出收敛曲线。这里有一个工程细节SP求解后拿到的u_new是连续值考虑到最坏场景应该取边界我拿到结果后会做一个归边处理把接近边界的值直接顶到边界上避免数值误差导致后续MP出现轻微的病态。5. 算例验证收敛曲线、成本对比与鲁棒性分析5.1 24时段微网算例参数设定为了验证CCG算法在MATLAB里到底能不能用、效果如何我搭了一个24时段的小微网算例。参数如下参数数值负荷峰值300kW晚高峰光伏预测峰值120kW波动±15%风电预测峰值60kW波动±25%预算参数Γ4燃气轮机容量20~80kW爬坡30kW/h燃机燃料成本系数0.5元/kWh储能容量/功率200kWh / 50kW购电价/售电价峰时0.9/0.7元谷时0.4/0.2元切负荷惩罚1000元/kWh弃风弃光惩罚50元/kWh这个算例的大小大概是24时段第二阶段连续变量约240个第一阶段整数变量48个不确定变量48个光伏风电各24个时段SP的MILP规模也就在几百个变量量级Gurobi求解速度飞快。5.2 CCG迭代收敛过程解读我记录的迭代过程如下迭代次数LB元UB元相对间隙/%15623.406054.187.1125867.325902.750.6035878.025884.610.1145879.105880.030.016第1轮LB和UB之间差7%左右第2轮直接压到0.6%第4轮就已经到0.016%。这个收敛速度验证了CCG在小规模微网调度中的高效性。从曲线形态看UB是单调下降的LB是单调上升的。UB来自给定第一阶段方案后算出来的完整成本随着迭代找到更好的第一阶段方案UB只会越来越低LB随着MP中场景增多、松弛被收紧只会越来越高。两条线一夹最优值就被限定在一个很小的区间里。5.3 确定性调度与鲁棒调度的成本、切负荷对比作为对比我还跑了一个确定性方案风光出力直接取预测值目标函数里的(\max_{u\in U})消失。两种方案的结果差别很有意思对比项确定性调度鲁棒调度Γ4日前计划成本元5726.405798.21最坏场景下的实际总成本元6385.905879.10最坏场景切负荷量kWh1260最坏场景购电成本占比高中等确定性调度在预测场景下成本便宜了约1.3%但如果光伏、风电真的按最坏情况发生它需要切126kWh负荷按惩罚单价算总成本直接飙到6385.9元。鲁棒调度在实际最坏场景下的成本只有5879.1元比确定性方案在最坏场景下的成本低了约8%。为什么会有这种差距看调度方案就能明白。鲁棒方案在最坏场景的边界条件下会提前把储能充得更满燃气轮机的基荷安排得更保守主网购电计划也预留了更多后手。这些在预测场景下都是多余成本但在最坏场景下就是救命稻草。这就是经济调度里最核心的权衡用一点正常场景下的经济性牺牲换取极端场景下的成本可控和供电可靠性。实际项目里到底选多大的Γ不是拍脑袋定的而是要结合当地天气规律和风险偏好画一条Γ-成本曲线让管理层拿这个曲线去决策。6. 踩过的坑与调参心得6.1 对偶转换中的符号与方向问题我在搭SP模型时遇到的第一个拦路虎就是对偶方向。很多教材写对偶都是标准形式max问题对偶成min但实际建模时约束方向经常是反的。功率平衡是等式约束它的对偶变量是自由变量不是非负的。我在第一版代码里把所有对偶变量都设成了sdpvar(n,1)并加了0约束结果SP目标值算出来完全不对而且UB和LB根本不可能收敛。正确的做法是等式约束对应的对偶变量用自由sdpvar不等式约束对应的对偶变量取非负对应的取非正。如果你不想每次都对冲方向有个土办法把模型写在纸上先把SP内层写成min dz s.t. Az b的统一形式再写成对偶最后再敲进Yalmip。手推一遍永远比直接凭感觉写代码快。6.2 大M值的选取宁可大一点不能病态大M法线性化(\pi_i v_j)的时候M的取值直接决定Gurobi的数值稳定性。M太小的后果是约束起不到惩罚效果最坏场景被错误地限制住UB算出来偏低M太大又会让MILP的LP松弛变得很差分支定界效率下降甚至出现数值警告。我的经验做法是先跑一次纯LP版本的SP不引入0-1变量和预算观察对偶变量(\pi)的量级然后取M为最大对偶变量量级的10~20倍。比如这里(\pi)量级在几百到一两千之间我取了M10000。如果模型规模变大建议用gurobi的数值诊断工具看condition number。另一个备选方案是使用Gurobi的指示约束indicator constraintsYalmip里可以通过implies实现避免大M带来的数值问题缺点是部分老版本求解器不支持移植性差一些。6.3 第二阶段含整数变量时对偶失效怎么办这是我踩过最深的坑。两阶段鲁棒优化的标准CCG推导有一个隐含前提第二阶段内层问题必须是连续LP才能用强对偶转成单层MILP。如果第二阶段变量里有整数变量内层就不是LP强对偶定理不成立对偶模型求出来的目标值就不是原问题的可行成本整个迭代逻辑就崩了。微网调度里最典型的诱导变量就是储能的防同时充放整数变量。我第一次建模时把储能充放电状态放在第二阶段结果SP对偶完之后的目标值和原内层目标值对不上检查了很久才发现是这个原因。解决办法有三个我都试过把所有整数变量全部挪到第一阶段。这是最推荐的做法也是很多文献的默认假设。储能是否充放电、燃机是否启停、主网是否购电这些状态确实应该在日前就定下来放到第一阶段从物理意义上也说得通使用基于KKT条件的单层化方法把内层原问题的最优性条件包含互补松弛写进约束然后引入大M虽然能处理整数变量但模型规模和数值难度都上了一个台阶用迭代启发式方法比如SP内部再做一层小循环把整数变量枚举出来但这种方法只能针对极小规模实用性有限。经过这次教训我现在建模时的原则是第一阶段尽量把所有离散状态都吃掉第二阶段只保留连续调整量。这不只是为了CCG代码好写也符合微网实际运行中模式切换提前定、出力调整实时做的操作习惯。6.4 从盒式到预算先跑通再降保守性最后一个经验是关于调试顺序的。我第一次跑CCG时不确定性直接上预算约束和0-1选择变量SP是MILP大M线性化、符号规则、归边处理各种因素搅在一起出了问题根本不知道是算法逻辑错还是数值实现错排查成本极高。后来我换了一个更聪明的路径先把不确定集退化成纯盒式集合也就是去掉预算约束相当于(\Gamma \infty)这时SP的对偶问题里没有0-1变量不需要大M线性化最优(u_j)直接由对偶变量系数的符号决定整个SP就是一个纯LP跑得飞快。在这个简化版本上我先确认CCG的整体迭代逻辑、LB/UB收敛行为全部正确然后再把预算约束和0-1变量加回去解决大M数值问题。这样每一步出了问题定位范围就小得多。这个调试顺序其实反映了这类算法问题的通用方法论先在一个结构更简单但核心逻辑完整的版本上把骨架跑通再逐步增加复杂度。两阶段鲁棒CCG的坑大多数不在算法不懂而在逻辑被复杂的模型细节淹没了。最后再分享一个我在实际运行中验证过的小技巧CCG收敛后不要把最后那个(y^)直接拿去用建议把这个(y^)带回到SP重新求解一次确认UB对应的是同一个第一阶段方案。因为MP求解时(\eta)只是约束下界理论上它等于SP值但数值上可能差一两个数量级的尾数。这一步校验成本很低却能避免很多莫名其妙的收敛了但方案不可用问题。 微网两阶段鲁棒经济调度与CCG算法MATLAB实现
1. 为什么微网调度要从确定性走向两阶段鲁棒1.1 一次光伏骤降让确定性调度直接翻车去年年底我在做一个园区级微网的日前经济调度项目碰到的第一个问题就非常典型光伏和风电的预测曲线在早上看起来非常漂亮预测中午光伏能到120kW结果上午十点一片云飘过来光伏出力从75kW直接掉到35kW。当时用的是最传统的确定性调度所有风光出力都按预测值当作已知参数喂给优化模型燃气轮机出力和储能充放电曲线也是按这个预测排好的。实际运行中光伏一掉燃气轮机因为爬坡约束没法立刻补上储能又已经在低谷时段放得差不多了最后只能从主网买高价电那一天的综合运行成本比调度计划超支了12%还要多。这个经历让我意识到微网优化调度如果只盯着预测值做经济调度本质上是在赌天气。后来我把求解框架换成了两阶段鲁棒优化配合CCG列与约束生成算法在MATLABYalmip环境里重新搭建了整个日前经济调度模型。这不是一个简单的多套几个场景的思路而是整个优化问题的结构都变了把决策分成今天就要定下来的事和等不确定性揭晓后再调整的事两层然后专门针对最坏的不确定性做应对。1.2 确定性优化、随机优化、鲁棒优化三种思路的取舍在做微网调度之前得先分清三种主流建模思路的适用边界不然很容易选错方向。第一种是确定性优化模型长这样[ \min_{x} ; c^Tx \quad s.t. ; Ax \ge b ]这里的负荷、光伏、风电全用预测值。优点是模型简单LP或者MILP一求解就完事速度快。缺点也明摆着预测和实际的偏差直接变成电网的安全风险或者成本风险。实际工程中靠人工预留旋转备用去对冲误差本质上还是经验主义。第二种是随机规划Stochastic Programming思路是把不确定性描述成若干离散场景每个场景有一个概率目标函数变成期望成本最小。比如用蒙特卡洛采样生成1000个光伏出力场景然后求期望。随机规划的强项是有概率意义理论上很优雅但它的痛点也很实际场景多了模型爆炸场景少了概率分布又不准。而且随机规划对分布假设敏感如果你把分布形状搞错了解出来的期望最优可能还不如确定性调度。第三种就是鲁棒优化Robust Optimization它不问概率分布只给一个不确定集合把问题写成[ \min_{x} \max_{u \in U} ; c^Tx f(x,u) ]要求的是在集合内任意一个不确定性实现下方案都可行且代价可控。鲁棒优化不追求期望最优它追求的是最坏情况下的保证。代价是方案必然偏保守但换来的是安全性。两阶段鲁棒优化是鲁棒优化在决策有先后场景下的升级版一部分决策必须在不确性实现之前敲定另一部分决策可以等看到实际风光出力之后再调整。1.3 两阶段鲁棒在微网场景下的适配性为什么微网经济调度特别适合两阶段鲁棒因为微网的运行逻辑本身就是分层的。日前阶段燃气轮机要不要开、开几台、是否和主网签订购电协议这些涉及整数变量的决策必须提前确定执行层面没法等到实时再改。而在实时阶段光伏实际出力出来了燃气轮机出力、储能充放电、主网交换功率这些连续量是可以快速调整的。两个阶段决策的自由度完全不同恰好对应两阶段鲁棒优化的结构。用一句话概括这个思路第一阶段花钱买选项第二阶段根据最坏情况动用选项把成本压住。这套东西做出来之后我最大的感受是两阶段鲁棒不是让你去准确预测不确定性而是让你在不确定性真正发生的时候有一个事先算好的、能兜底的调整策略。对微网这种源荷双重波动、设备响应速度又各不相同的系统来说这个特性太重要了。2. 两阶段鲁棒经济调度的数学模型变量、目标与不确定集2.1 微网设备、变量划分与运行假设先把我搭建的微网系统讲清楚。为了兼顾论文级严谨和代码可实现我选了包含微型燃气轮机、储能、光伏、风电、主网连接和本地负荷的经典拓扑。微型燃气轮机MT最小出力20kW最大出力80kW爬坡约束30kW/h燃料成本用二次函数近似储能系统ESS容量200kWh最大充放电功率50kW充放电效率0.95SOC限制在0.1~0.9初始SOC为0.5光伏风电不可控出力落在不确定集合内主网交互允许购电和售电购售电价分时且有交互功率上限负荷本文先按确定值处理实际工程中可以方便地把负荷波动也塞进不确定集合。变量划分是两阶段鲁棒模型的核心。第一阶段决策变量是所有必须在不确定性实现前定下来的整数量和部分合同量燃气轮机启停状态 (y_{on,t})与主网的购/售电状态 (y_{grid,t})整数储能是否处于充电状态 (y_{ch,t})防止同时充放电的整数约束。第二阶段决策变量则是在看到实际风光出力后可以调整的连续量燃气轮机出力 (P_{MT,t})储能充放电功率 (P_{ch,t}, P_{dis,t})主网购/售电功率 (P_{buy,t}, P_{sell,t})弃光/弃风量 (P_{curt,t})、切负荷量 (P_{cut,t})。这里有两条运行假设需要说明。第一储能同时充放电约束是用整数变量表达的在日前调度里基本是标配不然会出现充电又放电这种物理上不成立但数学上省钱的操作。第二我把切负荷和弃风弃光放进了第二阶段并赋予很高的惩罚系数目的是保证模型在极端场景下有退路同时也让经济性指标更有解释力。2.2 目标函数与约束的完整建模两阶段鲁棒优化目标的标准写法是[ \min_{y} ; c^T y \max_{u \in U} ; \min_{z \in F(y,u)} ; d^T z ]其中 (y) 是第一阶段决策(z) 是第二阶段调整变量(u) 是不确定变量。第一阶段的成本 (c^Ty) 包括燃气轮机启停成本我按单台单次启停300元计。第二阶段成本 (d^Tz) 则包含燃气轮机燃料成本用分段线性化逼近二次项向主网购电成本减去售电收入弃风弃光惩罚和切负荷惩罚后者单价设为1000元/MWh远高于正常发电成本。功率平衡约束是模型的主心骨对每个时段 (t)[ P_{MT,t} P_{dis,t} - P_{ch,t} P_{PV,t} P_{WT,t} P_{buy,t} - P_{sell,t} P_{cut,t} P_{load,t} P_{curt,t} ]这里我把弃风弃光量放在等式右侧表示多余电量可以被主动放弃切负荷量放在左侧表示需求不足时允许切除部分负荷。燃气轮机约束包括出力上下限、爬坡约束[ P_{MT}^{min} y_{on,t} \le P_{MT,t} \le P_{MT}^{max} y_{on,t} ][ -P_{ramp} \le P_{MT,t} - P_{MT,t-1} \le P_{ramp} ]储能约束包括功率限制、SOC递推、SOC上下限[ SOC_{t1} SOC_t \frac{\eta_{ch} P_{ch,t} - P_{dis,t}/\eta_{dis}}{E_{ESS}} ][ 0 \le P_{ch,t} \le P_{ch}^{max} y_{ch,t}, \quad 0 \le P_{dis,t} \le P_{dis}^{max} (1-y_{ch,t}) ]主网交互约束比较简单购售电功率各有一个上限且通过整数变量保证不会同时买卖[ 0 \le P_{buy,t} \le P_{grid}^{max} y_{grid,t}, \quad 0 \le P_{sell,t} \le P_{grid}^{max} (1-y_{grid,t}) ]模型到这里还是比较标准的。从代码实现角度看这些约束绝大多数只涉及第二阶段变量只有启停状态、购售电状态这些整数变量来自第一阶段这正是CCG算法适用的典型结构。2.3 不确定集设计盒式集合与预算约束两阶段鲁棒优化的鲁棒程度完全由不确定集合 (U) 决定。集合越大方案越保守成本越高集合越小方案越激进风险越大。这一步如果处理不好后面所有计算都是白搭。我最开始用的是最朴素的盒式不确定集[ U \left{ u ;|; u_j^{min} \le u_j \le u_j^{max}, \forall j \right} ]每个不确定量独立地在区间内波动。光伏预测120kW允许±15%那光伏出力就在102kW到138kW之间波动风电预测60kW允许±25%那就在45kW到75kW之间波动。但盒式集合的麻烦在于它会让模型默认所有不确定量同时取最坏值这是极端保守的。真实场景中光伏和风电同时达到最坏边界的概率极低用盒式集合做出来的调度成本会明显偏高。因此我在模型中引入了预算约束Budget Constraint[ \sum_j \frac{|u_j - u_j^{pred}|}{\Delta u_j} \le \Gamma ]这个 (\Gamma) 就是总偏差预算它限制的是所有不确定量偏离预测值的总幅度。比如24个时段(\Gamma4) 意味着最多等效于4个时段的不确定变量同时取到边界其他时段只能在预测值附近小幅波动。(\Gamma) 越大越保守(\Gamma0) 等价于确定性调度(\Gamma24) 等价于盒式集合。预算约束的好处非常多。第一它让不确定集合的体积可调给了调度员一个直观的经济鲁棒性权衡旋钮第二它在数学上保持了线性结构配合CCG迭代可以采用MILP求解第三实践中最坏场景一定会落在边界点组合上这让后面的对偶和线性化变得有迹可循。3. CCG算法拆解主问题-子问题迭代如何求解min-max-min3.1 为什么Benders不够CCG更合适第一次看到两阶段鲁棒模型的时候我的第一反应是能不能用Benders分解毕竟Benders分解在随机规划和确定性混合整数规划里是很成熟的技术。后来我仔细对比之后发现CCG在处理两阶段鲁棒优化上明显更顺手。Benders分解的思路是主问题求解第一阶段决策子问题求解第二阶段对偶然后把子问题的最优值函数用一条对偶割最优割或可行性割加回主问题迭代逼近原问题。问题在于子问题返回的信息是一个割平面本质上是目标值的下界逼近如果第二阶段存在整数变量Benders切割的处理会非常麻烦实际实现中经常出现振荡或收敛慢。CCG的核心区别是子问题每迭代一次生成一个具体的最坏场景 (u^*)并把该场景对应的第二阶段变量和全部约束原封不动地加入主问题。这相当于在主问题里新增一列和一组约束所以叫列与约束生成。因为每一轮都引入了完整场景的信息而不是一个压缩的割平面CCG的收敛速度通常远快于Benders。对于微网调度这种规模不太大、但场景约束很多的模型CCG往往在3~5次迭代内就能给出满意的间隙。3.2 主问题MP的构造与子问题SP的对偶化CCG迭代的第一步是把原始的两阶段模型拆成主问题MP和子问题SP。主问题MP在每一轮迭代时包含当前已发现的全部最坏场景。假设已经迭代到第 (k) 轮有 (k) 个已知场景 (u_1^, u_2^, ..., u_k^*)MP可以写为[ \min_{y, z_k, \eta} ; c^T y \eta ][ s.t. ; \eta \ge d^T z_k, \quad \forall k ][ z_k \in F(y, u_k^*), \quad \forall k ][ y \in Y, ; z_k \in Z ]辅助变量 (\eta) 的作用是逼近第二阶段的最坏成本下界。随着场景增多MP被约束逐步收紧目标函数值从下往上逼近真实最优值这就是下界LB。子问题SP是在给定第一阶段解 (y^*) 后找最坏场景和对应的最小第二阶段成本[ Q(y^*) \max_{u \in U} ; \min_{z} ; d^T z ][ s.t. ; Gz \ge h - Ey^* - Mu ]内层是关于连续变量 (z) 的LP因此可以用强对偶定理把它转成单层最大化问题。设 (\pi) 为内层约束对应的对偶变量SP等价为[ \max_{u, \pi} ; \pi^T (h - Ey^* - Mu) ][ s.t. ; G^T \pi \le d, \quad \pi \ge 0, \quad u \in U ]这里会出现一个双线性项 (\pi^T M u)这是整个CCG实现中第一个真正的技术难点。如果 (M) 矩阵的每个不确定变量只出现在一个约束里比如光伏只出现在功率平衡约束那么这项可以写成 (u)的线性组合系数由 (\pi) 决定的形态。当不确定集合是盒式时最优 (u) 一定取边界方向由 (\pi^T M) 的符号决定当加上预算约束时需要引入0-1变量选择哪些变量取边界SP就变成一个MILP。3.3 迭代流程、上下界更新与收敛条件CCG的整体流程我总结成下面这个循环这也是我在MATLAB里实际实现的骨架初始化给定初始场景 (u_1)我习惯取预测值设置 (LB -\infty)(UB \infty)迭代次数 (k1)求解MP得到最优解 (y^, \eta^)更新 (LB c^Ty^* \eta^*)固定 (y^)求解SP得到 (Q(y^)) 和新的最坏场景 (u_{k1})更新 (UB \min(UB, c^Ty^* Q(y^*)))计算相对间隙 (\frac{UB-LB}{|UB|})如果小于设定容差比如1e-3停止迭代否则把场景 (u_{k1}) 对应的第二阶段变量 (z_{k1}) 和约束 (F(y, z_{k1}, u_{k1})) 加入MP(k k1)回到第2步。这里有个细节值得注意MP目标值是松弛解所以LB一定小于等于原问题最优值SP给出的 (c^Ty^* Q(y^*)) 是某个可行方案的完整成本所以UB一定大于等于原问题最优值。两个界从两边夹逼收敛判据是有严格意义的。我在第一次跑通时看到LB和UB逐步靠近那种感觉很直观比Benders那种割平面越来越多但界还是慢吞吞动的体验好太多。4. MATLABYalmip环境下的CCG实现全流程4.1 环境准备Yalmip与求解器配置MATLAB这边我用的版本是R2023b建模工具是Yalmip求解器用Gurobi 10.x。这里多啰嗦一句两阶段鲁棒模型的内层和外层都涉及MILPYalmip自带的默认求解器是做不了大规模MILP的必须配置好商业求解器。Gurobi和CPLEX都行我个人偏Gurobi原因后面踩坑部分会说。Yalmip的安装很简单把文件夹加进MATLAB路径就行然后用yalmiptest验证一下求解器是否被识别。求解器配置代码ops sdpsettings(solver, gurobi, verbose, 2); ops.gurobi.MIPGap 1e-4; ops.gurobi.TimeLimit 600;我的习惯是外层CCG的迭代容差设1e-3内层Gurobi的MIPGap设1e-4。内层比外层紧一个量级可以避免外层还没收敛、内层精度不够导致的UB抖动。4.2 主问题MP的代码框架在Yalmip里建MP核心是按场景循环累加约束。因为MP会随着迭代不断长大我建议一开始就用cell数组或者多维sdpvar把第二阶段变量组织好避免在循环里反复创建变量。我采用的变量组织方式是T 24; % 时段数 nMT 1; % 燃机台数 nCont 4; % 第二阶段连续变量个数P_MT, P_ch, P_dis, P_buy, P_sell... K 0; % 当前已入库场景数 y_on binvar(nMT, T, full); % 燃机启停第一阶段 y_grid binvar(1, T, full); % 主网购售电状态第一阶段 z cell(K, 1); % 每个已发现场景对应的第二阶段变量 eta sdpvar(1);主问题约束的构建遵循一个套路先写第一阶段公共约束再对每个已发现场景 (k)把第二阶段约束完整写一遍。功率平衡约束的代码示例C []; for k 1:K for t 1:T u_PV u_set{k}(t, 1); u_WT u_set{k}(t, 2); C [C, z{k}(1,t) z{k}(2,t) - z{k}(3,t) ... % P_MT P_dis - P_ch z{k}(4,t) - z{k}(5,t) u_PV u_WT - z{k}(6,t) P_load(t)]; % 其中 z(k,:) [P_MT, P_dis, P_ch, P_buy, P_sell, P_cut] end end目标函数里辅助变量(\eta)对应所有场景第二阶段成本的最大值下界所以MP目标就是启停成本加上(\eta)obj sum(sum(c_on .* y_on)) 1e-6 * sum(y_on(:)) eta;这里1e-6的微调是为了防止整数变量出现不希望的退化解属于经验性的小技巧。MP是MILP第一轮没有场景时可以用一个空场景或者直接设MP只有第一步触发实际操作中我会先手动把初始场景(u_1)加进MP保证模型非退化。4.3 子问题SP与双线性项线性化代码SP的实现是CCG里最考验功底的部分。我先把一般形式列出来再给代码。SP内层LP写成标准形式[ \min_z d^T z \quad s.t. ; Gz \ge h - Ey^* - Mu ]对偶后是[ \max_{\pi \ge 0} ; \pi^T(h - Ey^* - Mu) \quad s.t. ; G^T\pi \le d ]在Yalmip里直接建模这个对偶问题pi_var sdpvar(n_cons, 1); % 对偶变量 v binvar(n_u, 1, full); % 不确定变量是否取上界的0-1指示 u u_lb (u_ub - u_lb) .* v; % 不确定量 C [G * pi_var d, pi_var 0]; C [C, sum(v) Gamma]; % 目标中的双线性项处理需要先定义 w(i,j) pi_var(i) * v(j) M_big 1000; w sdpvar(n_cons, n_u, full); for i 1:n_cons for j 1:n_u C [C, 0 w(i,j) M_big * v(j)]; C [C, pi_var(i) - M_big*(1-v(j)) w(i,j) pi_var(i)]; end end obj_sp pi_var * (h - E * y_star) - sum(sum(M .* w));这里的目标函数用了一个数学等价(\pi^T M u \sum_i \sum_j \pi_i M_{ij} u_j)把 (u_j) 展开成 (u_{lb,j} (u_{ub,j} - u_{lb,j})v_j) 后常数项直接保留变量部分用大M线性化结果(w_{i,j})来替代 (\pi_i v_j)。这样SP就是一个标准MILPGurobi可以直接求解。不过说实话大M法线性化不是唯一选择也未必是最好选择。我在迭代了几轮后发现如果模型规模不大还有一种更省事的做法因为(\pi)的最优值在LP极点处取得而(u_j)在最优时必取边界可以先用符号规则求解一个线性目标再通过分支定界选择预算。这个技巧比较复杂新手阶段建议还是稳妥地走大M法先把流程跑通。4.4 CCG主循环与LB/UB收敛过程输出主循环代码写出来非常短但逻辑要非常小心。伪代码如下LB -1e6; UB 1e6; K 1; u_set{K} u_pred; % 初始场景取预测值 gap 1; while gap 1e-3 K 10 % 1. 建立并求解MP [MP, obj_MP] build_MP(u_set, K, ...); optimize(MP, obj_MP, ops); y_star value(y_on); LB max(LB, value(obj_MP)); % 2. 固定y_star求解SP [SP, obj_SP_val, u_new] solve_SP(y_star, ...); UB min(UB, value(obj_SP_val)); % 3. 计算gap gap abs((UB - LB) / abs(UB)); % 4. 入库新场景 K K 1; u_set{K} u_new; end每轮迭代我都会记录LB、UB、gap、MP求解时间、SP求解时间最后用plot画出收敛曲线。这里有一个工程细节SP求解后拿到的u_new是连续值考虑到最坏场景应该取边界我拿到结果后会做一个归边处理把接近边界的值直接顶到边界上避免数值误差导致后续MP出现轻微的病态。5. 算例验证收敛曲线、成本对比与鲁棒性分析5.1 24时段微网算例参数设定为了验证CCG算法在MATLAB里到底能不能用、效果如何我搭了一个24时段的小微网算例。参数如下参数数值负荷峰值300kW晚高峰光伏预测峰值120kW波动±15%风电预测峰值60kW波动±25%预算参数Γ4燃气轮机容量20~80kW爬坡30kW/h燃机燃料成本系数0.5元/kWh储能容量/功率200kWh / 50kW购电价/售电价峰时0.9/0.7元谷时0.4/0.2元切负荷惩罚1000元/kWh弃风弃光惩罚50元/kWh这个算例的大小大概是24时段第二阶段连续变量约240个第一阶段整数变量48个不确定变量48个光伏风电各24个时段SP的MILP规模也就在几百个变量量级Gurobi求解速度飞快。5.2 CCG迭代收敛过程解读我记录的迭代过程如下迭代次数LB元UB元相对间隙/%15623.406054.187.1125867.325902.750.6035878.025884.610.1145879.105880.030.016第1轮LB和UB之间差7%左右第2轮直接压到0.6%第4轮就已经到0.016%。这个收敛速度验证了CCG在小规模微网调度中的高效性。从曲线形态看UB是单调下降的LB是单调上升的。UB来自给定第一阶段方案后算出来的完整成本随着迭代找到更好的第一阶段方案UB只会越来越低LB随着MP中场景增多、松弛被收紧只会越来越高。两条线一夹最优值就被限定在一个很小的区间里。5.3 确定性调度与鲁棒调度的成本、切负荷对比作为对比我还跑了一个确定性方案风光出力直接取预测值目标函数里的(\max_{u\in U})消失。两种方案的结果差别很有意思对比项确定性调度鲁棒调度Γ4日前计划成本元5726.405798.21最坏场景下的实际总成本元6385.905879.10最坏场景切负荷量kWh1260最坏场景购电成本占比高中等确定性调度在预测场景下成本便宜了约1.3%但如果光伏、风电真的按最坏情况发生它需要切126kWh负荷按惩罚单价算总成本直接飙到6385.9元。鲁棒调度在实际最坏场景下的成本只有5879.1元比确定性方案在最坏场景下的成本低了约8%。为什么会有这种差距看调度方案就能明白。鲁棒方案在最坏场景的边界条件下会提前把储能充得更满燃气轮机的基荷安排得更保守主网购电计划也预留了更多后手。这些在预测场景下都是多余成本但在最坏场景下就是救命稻草。这就是经济调度里最核心的权衡用一点正常场景下的经济性牺牲换取极端场景下的成本可控和供电可靠性。实际项目里到底选多大的Γ不是拍脑袋定的而是要结合当地天气规律和风险偏好画一条Γ-成本曲线让管理层拿这个曲线去决策。6. 踩过的坑与调参心得6.1 对偶转换中的符号与方向问题我在搭SP模型时遇到的第一个拦路虎就是对偶方向。很多教材写对偶都是标准形式max问题对偶成min但实际建模时约束方向经常是反的。功率平衡是等式约束它的对偶变量是自由变量不是非负的。我在第一版代码里把所有对偶变量都设成了sdpvar(n,1)并加了0约束结果SP目标值算出来完全不对而且UB和LB根本不可能收敛。正确的做法是等式约束对应的对偶变量用自由sdpvar不等式约束对应的对偶变量取非负对应的取非正。如果你不想每次都对冲方向有个土办法把模型写在纸上先把SP内层写成min dz s.t. Az b的统一形式再写成对偶最后再敲进Yalmip。手推一遍永远比直接凭感觉写代码快。6.2 大M值的选取宁可大一点不能病态大M法线性化(\pi_i v_j)的时候M的取值直接决定Gurobi的数值稳定性。M太小的后果是约束起不到惩罚效果最坏场景被错误地限制住UB算出来偏低M太大又会让MILP的LP松弛变得很差分支定界效率下降甚至出现数值警告。我的经验做法是先跑一次纯LP版本的SP不引入0-1变量和预算观察对偶变量(\pi)的量级然后取M为最大对偶变量量级的10~20倍。比如这里(\pi)量级在几百到一两千之间我取了M10000。如果模型规模变大建议用gurobi的数值诊断工具看condition number。另一个备选方案是使用Gurobi的指示约束indicator constraintsYalmip里可以通过implies实现避免大M带来的数值问题缺点是部分老版本求解器不支持移植性差一些。6.3 第二阶段含整数变量时对偶失效怎么办这是我踩过最深的坑。两阶段鲁棒优化的标准CCG推导有一个隐含前提第二阶段内层问题必须是连续LP才能用强对偶转成单层MILP。如果第二阶段变量里有整数变量内层就不是LP强对偶定理不成立对偶模型求出来的目标值就不是原问题的可行成本整个迭代逻辑就崩了。微网调度里最典型的诱导变量就是储能的防同时充放整数变量。我第一次建模时把储能充放电状态放在第二阶段结果SP对偶完之后的目标值和原内层目标值对不上检查了很久才发现是这个原因。解决办法有三个我都试过把所有整数变量全部挪到第一阶段。这是最推荐的做法也是很多文献的默认假设。储能是否充放电、燃机是否启停、主网是否购电这些状态确实应该在日前就定下来放到第一阶段从物理意义上也说得通使用基于KKT条件的单层化方法把内层原问题的最优性条件包含互补松弛写进约束然后引入大M虽然能处理整数变量但模型规模和数值难度都上了一个台阶用迭代启发式方法比如SP内部再做一层小循环把整数变量枚举出来但这种方法只能针对极小规模实用性有限。经过这次教训我现在建模时的原则是第一阶段尽量把所有离散状态都吃掉第二阶段只保留连续调整量。这不只是为了CCG代码好写也符合微网实际运行中模式切换提前定、出力调整实时做的操作习惯。6.4 从盒式到预算先跑通再降保守性最后一个经验是关于调试顺序的。我第一次跑CCG时不确定性直接上预算约束和0-1选择变量SP是MILP大M线性化、符号规则、归边处理各种因素搅在一起出了问题根本不知道是算法逻辑错还是数值实现错排查成本极高。后来我换了一个更聪明的路径先把不确定集退化成纯盒式集合也就是去掉预算约束相当于(\Gamma \infty)这时SP的对偶问题里没有0-1变量不需要大M线性化最优(u_j)直接由对偶变量系数的符号决定整个SP就是一个纯LP跑得飞快。在这个简化版本上我先确认CCG的整体迭代逻辑、LB/UB收敛行为全部正确然后再把预算约束和0-1变量加回去解决大M数值问题。这样每一步出了问题定位范围就小得多。这个调试顺序其实反映了这类算法问题的通用方法论先在一个结构更简单但核心逻辑完整的版本上把骨架跑通再逐步增加复杂度。两阶段鲁棒CCG的坑大多数不在算法不懂而在逻辑被复杂的模型细节淹没了。最后再分享一个我在实际运行中验证过的小技巧CCG收敛后不要把最后那个(y^)直接拿去用建议把这个(y^)带回到SP重新求解一次确认UB对应的是同一个第一阶段方案。因为MP求解时(\eta)只是约束下界理论上它等于SP值但数值上可能差一两个数量级的尾数。这一步校验成本很低却能避免很多莫名其妙的收敛了但方案不可用问题。