做微网优化调度的人应该都遇到过这种情况光伏预测出力是50kW但云一飘过来实际只剩30kW。按预测值做好的日前调度计划立刻失效购电曲线、储能充放电安排全得推翻重来极端情况下还会触发切负荷。传统确定性模型把预测值当确定值用调度结果在理想场景下非常好看一落到实际运行就是“纸面最优、现场打脸”。这两年我一直在折腾两阶段鲁棒微网优化调度结合关键场景辨别算法在Matlab里做了一套可以跑通、可以分析结果的实现。今天这篇把整套方案的关键逻辑和实现细节梳理出来适合正在做微网调度、储能调度、新能源消纳的研究生和工程师参考尤其是已经接触过鲁棒优化但卡在代码实现层面的那批人。这套方案解决的核心问题就是一句话在风光和负荷都有很强不确定性的情况下怎么提前制定一个“最坏情况下也不至于太亏”的调度计划同时计算量还能控制在工程可接受范围内。1. 微网不确定性处理的必要性与方案选型1.1 确定性调度模型为什么不够用我刚接触微网调度时第一个模型就是把光伏预测出力、风电预测出力、负荷预测曲线直接塞进去配一个储能和燃气轮机用YALMIP一调很快出结果。当时觉得“优化调度也不过如此”直到拿历史数据回测才发现问题严重。用一个最简单的例子来说明某时刻光伏预测出力是100kW但实际可能只有70kW差额30kW。假设微网本身没有旋转备用储能又在按照预测值放电那这30kW的缺口只能由主网购电来补。如果日前计划里的购电上限已经卡死系统就会被迫切负荷或者买高价电。一次两次还好长期运行下来惩罚成本累积起来非常可观。这就是确定性模型的缺陷它把不确定参数当成固定常数省略了一个重要的风险维度。只要预测误差存在所谓最优解就只是“预测场景下的最优解”而不是“所有可能场景下都可靠的解”。对于微网这种小规模系统这个问题被放得更大。微网容量小、惯性小、调节手段有限不像大电网有强大的频率支撑和备用容量。风电、光伏的间歇性、负荷的随机波动叠加储能SOC限制、燃气轮机爬坡约束调度方案稍有偏差就可能违反运行约束。所以不确定性建模不是一个可选项而是必须做的功课。1.2 三种主流不确定性处理方案的对比做不确定性调度的主流方案粗分下来是三类随机规划也叫场景规划。先给不确定参数假设概率分布再抽样生成大量场景最后求期望成本最小化方案。优点是比较贴近真实缺点是需要准确的概率分布而且为了精度经常要几千个场景求解规模蹭蹭涨。传统鲁棒优化。只给定不确定参数的上下界区间不关心概率保证最坏情况下系统依然可行。优点是模型稳健、不需要概率信息缺点是太保守了所有变量都在最坏点取值真实成本往往远低于最优值。两阶段鲁棒优化。前两个方案的折中。第一阶段做预调度决策第二阶段在不确定场景揭示后再进行低成本调整目标是最小化“最坏场景下的总成本”。保守性比传统鲁棒低稳健性比随机规划高数学结构上正好对应微网“日前预调度日内调整”的实际操作模式。从这个对比可以看出两阶段鲁棒在微网场景里不是故意炫技而是物理过程与优化结构的高度匹配。微网的运行方式就是分层的日前定启停、定购电计划日内根据实际风光出力重新调整储能充放电和机组出力。两阶段鲁棒优化天然地把这两个层面拆开了。1.3 两阶段鲁棒为什么适配微网场景微网调度里有个典型现象储能充放电决策有“后效性”。今天多充一度电明天可能就有放电空间今天的决策会约束明天。这种跨时间耦合的决策特别适合两阶段结构。第一阶段做的是需要提前承诺、无法快速反悔的决策第二阶段做的是发生时可以灵活微调的决策。表两阶段决策变量划分逻辑对比决策类型划分到哪一阶段原因机组启停状态第一阶段日前启停有最小运行/停机时间约束不能随时改变与主网的日前购电计划第一阶段日前购电协议一般提前签订实时市场波动大储能充放电计划第一阶段定基准第二阶段修正储能响应快日内可以根据实际出力调整燃气轮机实时出力第二阶段日内分钟级响应可按偏差调整切负荷量第二阶段日内作为惩罚手段留到最后用弃风弃光量第二阶段日内应对大发时段的风光超发这样划分之后优化问题自然变成min-min-max结构第一阶段先做预调度第二阶段面对最恶劣的不确定场景做矫正最终目标是最小化这两部分成本之和的上界。这个结构在后面章节会展开。2. 两阶段鲁棒微网模型的数学结构与关键假设2.1 目标函数与min-min-max结构两阶段鲁棒优化的目标函数表达方式比较固定但在给代码之前先要把这个结构拆明白。它的基本形式是min_{x} { c*x max_{xi in U} min_{y in F(x, xi)} d*y }其中x是第一阶段的预调度决策向量xi是不确定性向量光伏出力、风电出力、负荷U是不确定集合y是第二阶段的调整决策向量F表示在给定x和xi下y的可行域。理解这个结构的一个日常类比是出门旅行装箱先确定带哪几件大件行李x这是出发前必须定下来的到了目的地发现天气变化xi再做微调把某个厚外套换成薄外套y。你不可能为了所有可能的天气都准备一箱子衣服但你要保证最恶劣天气下的损失在自己的容忍度内。在微网里第一阶段成本包括机组启停成本、日前购电成本、储能充放电的基础损耗第二阶段成本包括实时购电调整成本、切负荷惩罚成本、弃风弃光惩罚成本。max部分负责在不确定集合里找出“对系统最不利”的场景min部分负责在这个场景下找到成本最小的调整策略。2.2 不确定集合的建模方式与参数含义不确定集合的选取直接影响解的保守性和求解难度。我在实现中用过两种效果都不错盒式不确定集。所有不确定变量都在预测值附近一个区间内波动U { xi : xi_min xi xi_max }这种方式简单直接但保守性偏强等于允许所有变量同时取最坏值。实际中风光和负荷不会同时到极限全部按最坏情况准备成本往往偏高。多面体不确定集也叫budget不确定集。在盒式基础上增加一个总量约束限制偏离预测值的总量不能超过某个阈值U { xi : xi_min xi xi_max, sum(abs(xi - xi_pred)) Gamma }Gamma就是鲁棒调节参数。Gamma取0时所有变量固定为预测值模型退化成确定性模型Gamma取最大可能值时等于盒式不确定集中间取值则可以在保守性和经济性之间调节。我一般把Gamma设置为可用不确定性的60%到80%回测下来成本和鲁棒性平衡得比较好。这里有个很多人会忽略的点Gamma不是一个神秘参数它的实际含义是“整个调度周期内允许不确定量偏离预测值的总幅度”。实际系统中可以统计历史预测误差的分布取90%分位数作为Gamma的参考值而不是拍脑袋定。我在代码里留了Gamma这个参数跑案例时建议都做一次敏感性分析看看成本随Gamma的变化曲线这个曲线本身就能给决策者提供非常有价值的风险偏好参考。2.3 微网模型的主要约束模型里的约束分三类。第一类是第一阶段决策本身的约束包括机组启停时间约束、日前购电容量约束、储能SOC初值约束。第二类是耦合约束把x和y联系起来比如功率平衡约束、储能SOC递推关系、机组出力上下限。第三类是第二阶段调整决策的可行域约束包括实时购电上下限、储能充放电功率限制、切负荷不超过总负荷的比例约束。功率平衡方程是核心形式是P_pv(xi) P_wt(xi) P_g P_es P_buy P_shed P_load(xi) P_curtail这个等式左边是供电侧右边是负荷侧。注意P_pv和P_load都是xi的函数这意味着不确定性通过功率平衡约束直接影响可行性。这个约束在确定性模型里只是个简单等式在鲁棒模型里却成了最关键的约束来源因为不同场景下等式右边不同第一阶段决策必须对所有场景都能找到可行的第二阶段调整方案。储能SOC递推关系也值得留意。储能是第二阶段主要的调节资源它的SOC更新公式在两个阶段之间建立了时间耦合SOC(t1) SOC(t) eta_ch * P_ch(t) - P_dch(t) / eta_dchSOC初始值在日前给定中间时刻的SOC在第二阶段根据实际风光出力调整最终SOC一般要求回到一个期望范围内这个约束会显著影响储能的灵活度。我在调参时发现如果最终SOC约束给得太紧两阶段模型很容易无解因为第二阶段找不到同时满足功率平衡和SOC终值的调整方案。3. 关键场景辨别算法的原理与实现思路3.1 枚举所有场景为什么不现实你可能会想既然要保证最坏场景下可行那我干脆把不确定集合里所有可能的场景都列出来找到成本最高的那个不就行了问题在于场景数量爆炸。就拿一个包含3台光伏、2台风电、1个负荷节点的小微网来说每个不确定性源取10个离散点组合起来就有10的6次方个场景。每个场景求解一个优化问题就算每个问题只要0.1秒也要10万秒也就是将近28小时。这还只是一个微网做调度方案优化肯定不能接受这种计算量。而且很多场景对调度方案的影响高度相似把它们全部考虑一遍纯粹是浪费算力。关键场景辨别算法要解决的就是这个问题能不能只找出少数几个对调度成本影响最大的场景用它们来代替整个不确定集合使得最终决策的安全性和全场景枚举几乎一致但求解规模大幅缩小。3.2 关键场景的动态辨识CCG方法我采用的核心算法是列与约束生成Column-and-Constraint GenerationCCG。你看名字可能觉得陌生实际上它的工作方式很好理解。它不是一个固定场景集合里的人肉筛选而是在迭代过程中动态地“长”出关键场景。具体思路分成四步初始化把名义场景也就是预测值场景加入场景集合。 求解主问题主问题是在当前场景集合下找最小化预调度成本加最坏调整成本上限的方案。 固定第一阶段决策求解子问题子问题是在给定x的前提下搜索不确定集合内让第二阶段成本最大的场景。 判断收敛如果新发现的最坏场景导致的总成本和主问题给出的下限之差小于容差停止否则把这个场景加入主问题生成新的约束重新求解主问题。这里的主问题和子问题交替求解每一次子问题识别出来的场景就是所谓的“关键场景”。它之所以关键是因为它在当前决策下对系统威胁最大不把它纳入主问题当前的决策就不安全。把它加进去之后主问题会调整决策然后再搜索新的威胁最大场景循环往复。这和场景削减的思路完全不同场景削减是事先从历史数据里挑代表CCG是从不确定集合里主动搜索危险点针对性更强。3.3 子问题为什么需要专门处理子问题本身是个max-min问题不能直接扔给求解器。一般做法是把内层min问题用对偶理论转换成max问题然后与外层max合并成一个单层max问题。这样做之后子问题变成了一个以对偶变量和不确定变量为决策变量的优化问题。但这里出现了一个新的麻烦目标函数里会出现对偶变量和不确定性变量的乘积项。这是典型的双线性项不是线性规划能处理的。我在实现中用的办法是大M线性化把双线性乘积用辅助变量和一组约束展开整体变成一个混合整数线性规划交给Gurobi求解。具体来说如果某两个连续变量的乘积是w_i * xi_j引入一个辅助变量z_ij来表示乘积再通过以下形式约束把z的取值限制在合理范围内z 0 z M * w z xi_j M * (1 - w) z xi_j - M * (1 - w)这里是拿一个二进制指示变量配合大M约束实现的M要取一个足够大的数但也不能大得离谱否则求解器数值稳定性会变差。我一般取不确定变量量纲的10到100倍之间。如果M太小可行域被错误缩小结果不对如果M太大对偶变量取值接近M时会产生大量小量级误差甚至出现无界解。3.4 与场景削减方法的对比有朋友问过我关键场景辨别和K-means聚类削减历史场景有什么区别。这里整理一下对比维度K-means场景削减关键场景辨别CCG输入大量历史或抽样场景不确定集合的描述区间、预算选择标准场景之间的相似度场景对调度成本的威胁程度场景来源原有场景集的子集或中心点动态生成的新场景不一定来自原场景集保守性可能低估最坏情况理论上收敛到最坏场景适用场景随机规划预处理鲁棒优化求核心场景两者可以配合使用。比如随机规划里先用聚类削减场景数再用期望成本做优化鲁棒优化里直接把CCG识别出来的关键场景作为鲁棒调度的依据。我自己在实际项目中先用CCG跑出关键场景再把这些关键场景作为随机规划的输入能得到一个兼顾鲁棒性和经济性的混合方案这算是一个工程上的取巧做法。4. Matlab实现与代码结构4.1 工具选型与求解器配置我的实现环境是Matlab配合YALMIP工具箱求解器用的Gurobi。YALMIP负责建模Gurobi负责求解线性规划和混合整数线性规划。为什么不用Matlab内置的linprog和intlinprog因为两阶段鲁棒的迭代过程要反复求解主问题、子问题内置求解器在MILP上的性能比不上商业求解器特别是不确定集比较复杂、约束较多的时候差距非常明显。如果手头没有GurobiCPLEX也是可以的YALMIP对两者的支持都很完善。跑代码前先确认环境yalmip(clear) options sdpsettings(solver, gurobi, verbose, 2, gurobi.MIPGap, 0.01);verbose设为2可以看清每步求解过程调试阶段有必要开。正式批量跑的时候我一般关掉只留汇总信息。4.2 主问题建模骨架主问题的核心是给定一个已经识别的关键场景集合最小化“预调度成本 一个表达最坏情况成本的上限值theta”。每加入一个关键场景就增加一组约束。这里给出一个示意性代码骨架变量命名简化过实际项目里请按自己的系统参数扩展% 第一阶段决策变量 x sdpvar(nx, 1); % 代表启停、购电计划等 theta sdpvar(1, 1); % 第二阶段最坏成本的上限 % 主问题约束集合 mp_constr [basic_constraints]; % 第一阶段基本约束 % 对每个已识别的关键场景添加约束 for k 1:size(xi_set, 2) xi_k xi_set(:, k); % 第k个关键场景 y_k sdpvar(ny, 1); % 该场景下的第二阶段变量 mp_constr [mp_constr, ... theta second_stage_cost(y_k), ... operation_constraints(x, y_k, xi_k)]; end % 求解主问题 optimize(mp_constr, first_stage_cost(x) theta, options); x_star value(x);这段代码有几个关键点。第一第二阶段变量y_k是为每个关键场景单独创建的副本场景之间互不共享变量这体现了鲁棒优化的核心思想不管哪个场景出现第一阶段决策都能为它找到合理的调整方案。第二theta是全局共享的它被约束在大于等于每个场景的第二阶段成本。第三随着关键场景数量增加主问题规模线性增长这就是为什么场景数量不能太多的直接原因。4.3 子问题与对偶处理骨架子问题负责在给定x_star时寻找最恶劣场景。直接写max-min问题没法求解需要先对偶。对偶后的标准形式是max_{lam, xi} { (h - E*x_star - G*xi) * lam } subject to: D*lam b, lam 0, xi in U这里lam是对偶变量xi是不确定变量。前面说过目标函数里lam和xi是乘积关系属于双线性项需要用大M法处理。示意代码如下% 对偶变量与不确定变量 lam sdpvar(n_dual, 1); xi sdpvar(n_xi, 1); % 双线性项线性化辅助变量 z sdpvar(n_dual, n_xi, full); bigM 1e4; % 大M约束把 lam_i * xi_j 替换为 z(i,j) for i 1:n_dual for j 1:n_xi sp_constr [sp_constr, ... z(i,j) bigM * lam(i), ... z(i,j) xi(j) bigM * (1 - binvar)], ... z(i,j) xi(j) - bigM * (1 - binvar), ... z(i,j) 0]; end end % 子问题目标取负号转为min形式 sub_obj -(h - E * x_star) * lam sum(sum(z .* G)); optimize([sp_constr, dual_constraints, uncertain_set_constr], sub_obj, options); % 得到最恶劣场景和最坏成本上界 xi_worst value(xi); UB value(sub_obj) first_stage_cost(x_star);注意z(i,j)的引入是因为双线性项展开后G*xi会乘以lam形成多组乘积项。这段代码里的binvar是临时引入的二进制变量数量等于双线性项个数所以子问题的求解规模会比主问题大一些。如果双线性项实在太多也可以换一种思路当不确定集合是盒式且维度不高时最恶劣场景通常出现在不确定区间的端点组合上可以枚举端点每个端点固定xi后直接求min问题取最坏值。这个方法在维度低于10时非常高效我的经验是能枚举端点就优先枚举计算量反而更可控。4.4 主循环与收敛判定有了主问题和子问题的求解模块主循环就清晰了。伪代码如下xi_set []; UB inf; LB -inf; gap 0.01; while abs(UB - LB) / abs(LB) gap % 第一步求解主问题得到x_star和theta_star % 主问题的目标函数值是当前场景集下的成本下界 LB value(theta) first_stage_cost(x_star); % 第二步固定x_star求解子问题得到最恶劣场景xi_worst if sub_problem_feasible(x_star) UB min(UB, value(sub_obj) first_stage_cost(x_star)); else % 子问题无可行解说明不确定集合里有场景无法满足约束 % 需要把相应的可行性割约束加入主问题 add_feasibility_cut(x_star); end % 第三步收敛判断 if abs(UB - LB)/abs(LB) gap break; end % 第四步把最恶劣场景加入场景集下一轮主问题重新优化 xi_set [xi_set, xi_worst]; end这里有一个容易被忽略的细节主问题求得的下界LB并不严格单调但UB会随着迭代逐步下降最终收敛。收敛判据里的gap我建议设成1%到5%不用追求0.1%的高精度因为微网调度参数本身的精度远达不到这个水平把gap设得太小只会白白增加迭代次数。我的实际经验是从2个初始场景开始通常5到7次迭代就能收敛到2%以内的gap求解时间在几十秒到几分钟视系统规模而定。4.5 参数设置与结果输出模型里的参数建议单独放在初始化脚本里不要散落在各段代码中。我习惯把系统参数分成三类物理参数机组容量、储能容量、爬坡率、线路容量、成本参数购电电价、燃气价格、切负荷惩罚、弃风弃光惩罚、鲁棒参数预测值、预测误差区间、Gamma值。每次跑新案例只改初始化脚本不动算法主体。结果输出方面至少要输出四样东西第一第一阶段决策变量的取值第二识别出来的关键场景列表第三每个关键场景对应的第二阶段调整策略第四迭代收敛曲线。第四样特别重要通过收敛曲线可以直观看出gap下降趋势如果曲线震荡剧烈说明代码里大概率存在bug需要回到子问题检查对偶推导或者大M约束的合理性。5. 常见问题与排查技巧实录5.1 求解器报错与连接问题跑代码时遇到最多的求解器问题是YALMIP找不到Gurobi。检查方法是在Matlab命令行输入yalmiptest看输出中Gurobi的状态是不是OK返回OK就不必担心。如果显示找不到Gurobi通常是安装路径没被Matlab识别或者license过期。还有一类问题是运行时报“License manager error -8”这个和Matlab自身license绑定有关检查license环境变量是否正确指向许可证文件即可。日常排查单上求解器问题的优先级最高因为问题出在这一层后面所有优化结果都是无效的。5.2 子问题对偶推导实现不一致对偶推导是最容易出错的地方。常见症状是子问题求出的目标值特别大甚至无穷大导致UB直接爆掉。排查时我一般先把xi固定为名义预测值只求解内层min问题并对比手算值如果这一步就错了问题十有八九出在约束漏写或对偶方向搞反。再进一步固定对偶变量lam检查目标函数里xi相关项的符号。大M线性化实现时二进制变量的指示意义必须和约束一一对应很容易出现写反的情况。顺着这个思路排查再麻烦的bug也能定位到。最大体会是不要试图一步把所有数学推导直接写成代码。老老实实先写一个干净的小规模测试用例比如2个节点、3个时段手算验证正确后再扩展可以节省大量时间。5.3 收敛慢或不收敛收敛慢的原因大体有三种。第一种是起始关键场景选得不好。我刚开始做的时候只拿预测场景作为初始化结果前面几轮迭代gap一直不降因为最恶劣场景离预测场景太远。后来把盒式端点组合也加进初始场景集收敛速度快了不少。第二种是gap设置太严格。鲁棒优化求解的是带近似的不确定性模型精确值本身就有建模误差gap设成0.001只会让循环多跑好几轮而结果几乎不变。合理做法是设置在0.01到0.03区间。第三种是大M参数不合理。大M取太小会截断可行域子问题永远找不到真正的恶劣场景取太大会带来数值问题Gurobi求解MILP时可能出现虚假的最优解。排查大M是否合适一个经验法是观察识别出来的最恶劣场景是否集中在不确定集边界。如果大量关键场景出现在不确定集合的内部说明大M约束很可能把边界截断了需要增大M值或者调整线性化约束形式。5.4 结果异常的合理性检查最后说说结果合理性检查。鲁棒优化模型很容易跑出一个“看起来合理但实际很怪”的解。我通常用一张核对表逐项检查检查项正常表现异常表现切负荷量只在少数极端场景出现量不大几乎所有场景都切负荷储能SOC轨迹峰谷变化明显两端约束留有余量频繁顶到上下限SOC波动很小弃光弃风量与光伏大发时段基本相关与天气完全不相关主购电出力跟随净负荷变化偏离净负荷曲线很远Gamma敏感性Gamma增大总成本缓慢上升Gamma微调成本剧烈跳变如果切负荷量在所有场景中都很大我会首先检查储能容量和爬坡约束是否设置过紧其次检查惩罚成本系数是不是远低于购电成本导致模型主动选择切负荷而不去购电。如果Gamma敏感性特别强烈说明模型对不确定集参数特别敏感这时候不能直接采信结果要先回归数据把预测误差区间和Gamma的真实分布校准一遍再重新跑优化。做微网调度的核心不只是把数学公式换成代码更重要的是理解每一个参数调节背后对应的物理含义。比如Gamma变大意味着系统要应对更严峻的源荷偏差成本自然上升储能容量加大两阶段调整空间变大鲁棒性变好但初始投资也变大。这些都是优化结果之外真正用于工程决策的信息。6. 实操经验与后续扩展思考整个项目跑下来我的感受是两阶段鲁棒优化并没有想象中那么遥远。它本质上就是把微网调度固有的“先承诺、后调整”过程数学化再通过关键场景辨别算法把计算复杂度控制住。比起把大部分场景都塞进一个模型追求理论最优我更认可这种“识别少数威胁场景、确保最坏情况不失控”的思路它对实际运行的指导意义更强。在实际工程中还有一个小技巧把优化得到的预调度方案和关键场景对应的第二阶段调整策略一起输出给运行人员。这样运行人员看到的不只是一条调度曲线而是一套完整的“如果光伏低于预期怎么办、如果负荷突然升高怎么办”的应对手册。这在小型微网运维场景里特别实用因为现场工作人员未必懂鲁棒优化但都看得懂一套场景化的预案表。如果后续想深化我建议顺着两条线走。一条是把分布鲁棒优化引入到模型里它用历史数据的矩信息构造不确定集比盒式不确定集更精确保守性更低但求解难度更高需要进一步研究线性化和分解算法。另一条是加入电动汽车充电负荷、空调温控负荷这类柔性负荷的调节能力它们作为第二阶段的弹性资源效果比单纯依赖储能更明显。等这两块跑通了这套框架基本就能覆盖实际微网运行的大部分不确定性来源了。