做微电网日前调度的人迟早要跟“不确定性”这三个字正面硬碰。风光出力说变就变负荷预测也不可能永远精确但调度方案又必须在前一天定下来。最怕的情况就是你按某个“确定值”排好了机组启停和购电计划第二天实际来风一吹安全约束直接崩掉。我最近在做的微电网两阶段鲁棒优化matlab代码就是专门解决这类问题的整套实现基于matlabyalmipcplex三件套用列与约束生成CCG算法把min-max-min结构拆成主问题和子问题交替求解。这篇文章适合正在做微电网经济调度、选修过最优化基础、或者想从确定性优化转向鲁棒优化的同学参考我会把建模思路、代码骨架、以及我在调试时踩过的坑全部写出来。1. 两阶段鲁棒优化到底在解决什么问题1.1 先从一个普通调度模型的不确定性说起常规的微电网经济调度模型长什么样目标函数一般是购电成本加机组燃料成本加弃风弃光惩罚约束有功率平衡、机组爬坡约束、储能充放电约束、联络线传输功率约束等。这本质上是一个混合整数线性规划MILP问题用yalmipcplex几分钟就能解完。但这套东西有个致命前提所有参数都是“给定”的。风电预测曲线是什么就用什么可实际运行中预测误差是必然存在的短时间尺度内偏差20%甚至更多都正常。如果调度方案建立在单点预测之上一旦实际出力偏低系统就得临时高价购电或者甩负荷如果实际偏高又得弃风弃光。确定性模型不是不能用而是它把风险成本完全忽略掉了。于是有两条进阶路线随机优化和鲁棒优化。随机优化要假设不确定性服从某个概率分布然后对大量场景求期望最优。它的缺点是场景数据获取成本和求解成本都很高而且分布假设一旦错了结果反而误导决策。鲁棒优化则走了另一个极端不关心概率只保证“所有可能情况都安全”但传统的单阶段鲁棒优化把所有决策都压在“最坏场景”上做结果往往保守到没法用机组全开、储能全备、购电拉满成本高得离谱调度员看一眼就想把它删掉。两阶段鲁棒优化是这两者之间的实用折中第一阶段做“看到不确定性之前必须提前定下来”的决策比如机组启停、日前购电计划第二阶段做“等不确定性暴露后再调节”的决策比如实际出力调整、储能充放电微调。它的思想是所有可调节的柔性留给第二阶段去消化第一阶段的决策只需要保证“无论最坏情况怎么来第二阶段的调整都能兜住”。这样既保证了安全性又不会把柔性提前浪费掉。1.2 min-max-min结构与CCG的核心逻辑两阶段鲁棒优化的数学形式核心就是一个min-max-min问题。外层的min对应第一阶段决策目标包含第一阶段成本以及一个“最坏情况下二期成本”的贡献。中间的max对应“大自然在不确定性集合里挑一个最不利的场景”让系统付出最大代价。内层的min对应第二阶段在固定了第一阶段决策和具体场景之后做出的最优调整。如果直接把这个三层问题丢给求解器它会直接崩溃内层和外层都是优化嵌在一起没法直接求解。标准解法是列与约束生成Column-and-Constraint GenerationCCG核心思路是用迭代逼近主问题MP把“最坏场景”替换成一组已经枚举出来的具体场景在保证这些场景下第二阶段都有可行调整的前提下求第一阶段成本和对应场景成本的总和得到当前的第一阶段决策和下界。子问题SP固定第一阶段决策去求解最坏场景下的调整成本得到一个新的“最坏场景”和对应的最优调整方案把这个场景对应的第二阶段约束和变量“列”添加到主问题里。两个问题交替求解上下界收敛到同一个值迭代停止。这个思路理解起来就像“谈判”第一阶段先出一个方案子问题找到这个方案最致命的场景主问题根据这个场景修补方案然后再被找茬再修补直到方案在任何场景下都站稳为止。需要注意CCG和Benders分解长得很像但区别在于CCG每次往主问题里加的是“新场景对应的完整第二阶段变量和约束”而Benders加的是一堆Benders割。CCG的收敛性通常明显更快尤其适合这种场景数不多但约束结构复杂的两阶段模型。2. 为什么选matlabyalmipcplex这套技术栈2.1 三件套的分工与优势先说说我为什么一直在用这套组合因为这是目前工程和学术界做两阶段鲁棒优化最顺手的配置。matlab负责数据和矩阵层面的组织调度。微电网调度问题天然就是矩阵问题节点注入功率、时间断面、机组状态、储能SOC这些数据用matlab的数组和结构体组织起来写起来和调试起来都比C或者Python的原生数据结构顺手。当然你可以说Pythonnumpy也不差但matlab在处理带时序的调度模型时矩阵思维上的顺畅程度还是明显更舒服。yalmip是建模层。它的意义在于把“优化问题”用接近数学表达式的语言写出来比如约束直接写 P_balance 0变量用 sdpvar 或 binvar 声明目标函数直接求和。这意味着从论文推导到代码实现之间几乎不需要翻译出错的概率自然就低。用yalmip建过模的人都懂一旦用久了你会非常嫌弃一个个手写Aeq、beq矩阵的方式。cplex是求解层。它是个商业MILP/QP求解器对大规模混合整数线性规划的求解速度和数值稳定性都相当能打。两阶段鲁棒优化的主问题会被求解很多轮次每一轮都是完整的MILP如果求解器不够强整个算法会被拖垮。cplex在这种高频次、重负载的MILP求解场景下稳定性确实比开源求解器要让人放心得多。2.2 为啥不用其他求解器和优化框架我并不是说只有这一套能用。Yalmip其实也支持gurobi、scip、osqp等求解器如果你非得用纯Pythonpyomo开源求解器也能凑合跑cvxpygurobi也是一种选择甚至有些人在用小众的鲁棒优化工具箱ROME但一是文档少二是出问题排查资料非常有限。我的具体建议是如果你是先用matlab做研究的就直接上yalmipcplex别犹豫如果你本身就在Python生态里可以试试gurobi的Python API建模自由度其实更高如果只是学生做课程设计用开源求解器scip顶替cplex也足够了但一旦模型规模扩大求解时间会肉眼可见地拉长。做两阶段鲁棒优化这种“求解器调用次数多到令人发指”的问题求解器的效率直接决定了你的实验能做多大、论文数据能填多满。cplex的强项在于MILP的预求解、cut生成和并行搜索这些能力在这种迭代算法里几乎每轮都会吃满。对比项matlabyalmipcplexPythonpyomogurobi开源求解器scip等建模便利度接近数学表达式也不错结构化更强一般MILP求解速度快快慢规模一大明显吃力调试资料量多多少建议人群matlab已有基础的科研人员Python生态用户课程设计、模型验证2.3 环境准备与版本注意事项基本环境是matlab2016以上的大多数版本都能跑、yalmip更新到最新版、cplex12.8或更高版本。这里有几个非常容易踩的坑我先列出来yalmip的版本和matlab版本兼容性问题通常表现为 “Undefined function or variable yalmiptest” 或者sdpvar创建后变量显示不出。解决方案删掉旧的yalmip路径重新添加新下载的路径确保只用一套。cplex在matlab中的API不是所有版本都能无缝对接yalmip。装完cplex后最好跑一次cplex.optimize测试或者直接用yalmiptest看看有没有识别到cplex。不要迷信最新版。我有一次换了某新版matlab之后旧版cplex动态库直接加载失败花了一个下午才搞清楚是版本兼容矩阵问题。稳定组合一句话matlab 2021a yalmip 2021 cplex 12.10这个组合我长期使用没有问题。3. 代码实现的关键环节与建模要点3.1 不确定性集合到底怎么建这是整个模型最关键的地方也是最容易出错的地方。不确定性集合U描述“大自然到底能在这个集合里怎么折腾”。集合的大小直接决定解的保守程度和计算负担。最常用的是盒式集合配合预算约束。例如风电预测功率为 w(t)^pre实际出力 w(t) 满足w(t) ∈ [w(t)^pre - w(t)^down, w(t)^pre w(t)^up]同时加上一个预算约束所有时段的偏差总和不超过预算Γ防止每个时刻都走到极端。这个Γ的本质是对不确定性“总变异量”的约束工程上常用 sqrt(T)T为时段数来作为经验取值例如24时段可取5左右。为什么要用预算约束因为如果没有它最坏场景必然是“每个时段的风电都取最小”这太极端了实际几乎不会发生会让结果过于保守。预算约束的本质是允许每个时段有偏差但限制总偏差的大小让最坏场景落在“虽然糟糕但还算有可能发生”的范围内。实现上在matlab里可以先定义偏差变量和上下限w_dev sdpvar(T, 1, full); % 每个时段的偏差量 w_act sdpvar(T, 1, full); % 实际风电出力也是子问题里的u变量 cons_uncer []; for t 1:T % 盒式约束delta_up和delta_down是各时段上下偏差范围 cons_uncer [cons_uncer, w_act(t) w_pre(t) - w_dev(t) * delta_down(t)]; cons_uncer [cons_uncer, w_act(t) w_pre(t) w_dev(t) * delta_up(t)]; end % 预算约束 cons_uncer [cons_uncer, sum(w_dev) Gamma, w_dev 0];注意不同时段允许的向上、向下偏差大小最好从历史预测误差统计得出不要拍脑袋定。另外Γ的取值会直接影响保守程度调试时可以多跑几组Γ做敏感性分析。3.2 主问题的建模与约束生成主问题在第K次迭代时包含第一阶段决策变量x以及前面迭代生成的K组场景对应的第二阶段变量y_k。模型可以写为min cx tau s.t. tau dy_k, k 1,...,K A x a B x C y_k E u_k g, k 1,...,K y_k 0用tau作为第二阶段最大成本的替代变量配合约束tau dy_k进行线性化。这是标准CCG主问题的写法好处是避免把K个场景的目标全部求和导致目标函数规模膨胀。用yalmip实现时关键是按迭代次数动态增加约束和变量。yalmip不像AMPL有集合成员机制一个实用的办法是先预分配一个足够大的变量数组然后每个迭代只给第k组变量赋新值Kmax 50; % 预留最大迭代次数 y2 sdpvar(2*T, Kmax, full); % 预分配第二阶段的决策变量矩阵 % 在第k次迭代中只对y2(:,k)添加约束和变量 MP_cons [MP_cons, B * x C * y2(:,k) g - E * u_k(:,k)]; MP_cons [MP_cons, y2(:,k) 0];这个预分配逐列激活的办法避免了在迭代中频繁重建整个模型求解速度会明显快很多。我第一次写的时候每次迭代都用[MP_cons, MP_cons_new]这种拼接方式跑24时段的算例时能明显感到模型一变大就卡顿改成预分配后流畅很多。3.3 子问题的对偶化实现固定第一阶段决策x之后子问题长这样Q(x) max_{u in U} min_{y} dy s.t. C y g - B x - E u y 0内层min问题是一个线性规划可以直接取对偶把min换成maxmax_{lambda 0} lambda (g - B x - E u) s.t. C lambda d于是原max-min问题变成max_{u in U, lambda 0} lambda (g - B x - E u) s.t. C lambda d这是一个单层优化但目标函数里出现了 lambda E u 的双线性项。如果u是连续变量双线性项处理起来比较麻烦通常用大M法引入辅助变量线性化。我的建议对中小规模算例直接对双线性项用大M处理或者用McCormick松弛这是最实用的做法虽然松弛会有微小gap但工程上完全可接受。如果u本身可以离散化为有限个预算组合也可以枚举场景再对每个离散u求解一个LP取最大目标值但枚举数量可能爆炸一般只在验证小算例时用。大M取值是个玄学。太多会让cplex在预求解阶段磨蹭很久太少又会把最优解砍掉。我一般先抛一个50~100的数量级试跑看最优解是否贴着边界再逐步减小。3.4 迭代收敛判据与代码骨架整个算法的迭代流程可以写成如下流程初始化下界LB -inf上界UB inf迭代次数k 0设置初始场景u_1为预测值的期望场景。求解主问题MP得到第一阶段决策x_k和主问题目标f_MP更新LB max(LB, f_MP)。固定x_k求解子问题SP得到最坏场景u_{k1}和对应目标Q(x_k)更新UB min(UB, cx_k Q(x_k))。如果(UB - LB) / abs(UB) epsilon比如0.01停止否则向主问题添加场景u_{k1}对应的变量和约束k k 1回到步骤2。这里有两个细节坑一是上下界更新方向别搞反LB从主问题来只增不小UB从子问题来只减不增。二是如果第二阶段问题规模较大子问题每次求解都比较费时可以考虑把多个历史场景同时加入主问题来加速收敛这是CCG的一种变体。核心代码骨架如下% 初始化 LB -inf; UB inf; k 1; u_k w_pre; % 初始最坏场景取预测值 while (UB - LB) / abs(UB) 0.01 % 求解主问题 MP_cons [MP_cons, B*x C*y2(:,k) g - E*u_k]; optimize(MP_cons, obj, opt_settings); LB max(LB, value(obj)); x_star value(x); % 固定x求解子问题 q_sub build_SP(x_star); % 对偶后得到的单层max问题 optimize(q_sub_cons, -q_sub_obj, opt_settings); q_val value(q_sub_obj); u_k1 value(w_act); UB min(UB, value(c*x) q_val); % 生成新场景列并迭代 u_k u_k1; k k 1; end4. 实操过程从零跑通两阶段鲁棒优化代码4.1 小算例先行我第一次写这类代码时上来就套了24时段、十几个节点、多台机组的完整模型结果调试时每一步都很难判断是模型错还是算法错。后来我就学乖了先做一个小算例。具体可以构造一个3节点微电网、4个时段、单台柴发、一组风电的场景甚至更极端把时段数缩减到2小时验证算法能收敛、上下界gap能下来、最坏场景不出意料再改回24时段。这种做法的价值在于小模型的变量和约束少任何一个约束写反了都能快速看出结果异常。同时你还可以手推一遍小规模的min-max问题用穷举法验证子问题的目标值是否合理。我经常用穷举法算2时段、风电两档取值的情况如果CCG给出的最坏场景和目标值跟穷举一致那对偶推导基本就是对的。4.2 核心参数的初始设置初始参数直接影响能不能跑通。我一般这样设置风电预测曲线和偏差范围用历史数据的误差分位数来定。预算Γ5对24时段先跑通再调。cplex的MIP gap给到0.01%~0.1%。默认值太紧会导致某些大规模场景长时间不收敛给的太松又会让每次主问题求解的结果偏差过大影响整个CCG的收敛。外循环epsilon取0.01或0.02研究用途取0.01足够。另外所有成本系数和功率单位最好统一。我吃过一次亏某个参数用的是MW另一个用的是kW导致约束严重违反结果看起来像模型错误排查了一整天才发现是单位问题。4.3 迭代过程中的数据记录跑CCG时如果不记录中间量基本等于瞎跑。我在每次迭代里都会把主问题目标、子问题目标、当前x的启停状态、最坏场景曲线存下来画在一张图上一眼就能看出收敛过程是否健康。可以画三张图第一张是gap收敛曲线正常应该是平滑单减如果大幅振荡说明模型或代码有问题第二张是最坏场景下的风电曲线它应该是有“针对性”的——哪里机组的爬坡约束最紧它往往就在哪里压低出力第三张是机组启停状态的演化能直观看到方案为何保守。这些图的价值在调试时怎么说都不为过。有一回我发现gap在0.15附近卡住不动看场景曲线才发现子问题每次返回的最坏场景几乎一样而主问题在该场景下的成本预估和子问题的目标对不上一查果然是第二阶段目标函数里的时段权重写错了。5. 常见问题与排查技巧实录5.1 求解器配置与yalmip报错最常见的坑是明明装了cplex但yalmip solver显示还是默认用sedumi或gurobi。此时在求解前可以调用solvesdp(MP_cons, obj, sdpsettings(solver,cplex))或者跑yalmiptest(cplex)看看返回状态。另一个坑是求解报“Infeasible problem”但不是模型真不可行而是yalmip把某类约束写成了非线性。检查方法用check(MP_cons)逐条检查约束残差看哪条不满足然后把连续变量转成双线性或添加辅助变量。也有可能是初始场景u_1选得太极端导致主问题在第一次迭代就不可行。解决方法是先把u_1设为预测值保证主问题至少有一个可行解然后再进入迭代。5.2 收敛困难与振荡如果gap震荡不降通常有三类原因一是子问题的对偶写错了导致返回的“最坏场景”并不是真正最坏主问题基于错误场景修正方案自然来回摆动二是主问题的第二阶段目标在目标函数中的表达方式不一致比如MP里用平均值、子问题用总和的混搭三是上下界更新逻辑有bug导致中间结果时大时小。排查方法在小算例上人工穷举几个典型的u固定当前x后分别算子问题目标和子问题的返回值对比如果对不上就对偶错了。另外一个笨办法是手动打印每轮迭代的最坏场景u_k看它是否随时间在几个固定值之间蹦跳如果是多半是场景生成逻辑有问题。5.3 结果不合常理的排查路径我遇到最搞笑的一次是模型结果里弃风量巨大、但购电量也巨大两相矛盾。后来发现是新能源出力与负荷之间的功率平衡约束写成了“可以放行任何功率”也就是不等式方向写反了导致可行域极其宽松。排查这种问题最快的方法是固定所有变量为已知合理值代入约束检查是否违反违反的那条就是bug所在。还有一类情况是结果“过度保守”所有时段都开机组、都买高价电。这种通常是不确定性集合设得太大尤其是预算Γ设得过高。可以先跑一组Γ0的算例如果结果回归正常就说明问题出在不确定性集合的参数上而不是模型错误。5.4 性能调优的几个土办法如果算例规模偏大有几个实用优化手段一是预分配变量数组避免每次迭代重建模型前面已经提过这是最有效的一招。二是把约束矩阵稀疏化用sparse构造cplex对这种格式处理更快。三是第二阶段目标函数单纯求和的项很多合并同类项后模型体积能小不少。四是如果第一阶段变量就是0-1的启停变量可以先用LP松弛的MP求一个近似解作为热启动节省不少时间。另外cplex有很多隐藏参数可以调。比如sdpsettings(cplex.timelimit, 300)可以限制每次主问题求解时间sdpsettings(cplex.mip.tolerances.mipgap, 1e-4)可以放宽MIP gap。对CCG这种需要反复求解的场景这些参数合在一起往往能让总时间缩短一半以上。最后聊点实在的。做微电网两阶段鲁棒优化很多人刚上手时容易陷进纯粹的数学推导里忽略了代码层面的大量细节但反过来只贪图“把模型跑通”而不理解min-max-min结构和CCG的迭代逻辑那换一个场景你立刻就不会改。我个人体会最深的一件事是这类代码真正值钱的地方不在于某一次能跑出多漂亮的结果而在于你把不确定性集合从盒式改成椭球式、加入储能寿命约束、或者把单微电网扩成多微电网后还能不能快速修改并稳定出结果。这个扩展能力才是这套matlabyalmipcplex代码长期的价值所在。如果你也正在做类似的工作建议一定从小算例起步先把对偶推理和调试链路打通再考虑大规模场景的事。
