两阶段鲁棒优化在微网调度中的Matlab实现与CCG算法详解

很久之前就想写一篇关于微网调度鲁棒优化的文章因为这一块理论公式多、实现门槛高网上大多是PPT级科普真正能跑通的Matlab代码却少得可怜。这个项目正是冲着这个痛点去的用两阶段鲁棒优化处理风电、光伏和负荷的不确定性再配合关键场景辨别算法也就是CCG列与约束生成的思想逐步找出最恶劣场景并迭代求解最后在Matlab里完整落地。核心关键词就四个两阶段鲁棒、微网调度、不确定性、关键场景辨别。这篇文章会把模型怎么搭、算法怎么迭代、代码怎么写、坑怎么踩一次讲透适合正在做微网优化调度课题的研究生、刚接触鲁棒优化的工程师以及想从确定性调度转鲁棒调度的朋友。读完不是只懂概念而是能直接照着改写自己的模型。1. 项目整体设计与思路拆解1.1 微网调度为什么必须考虑不确定性微网调度和传统大电网调度的最大区别在于“底子薄”。风电、光伏出力随天气波动负荷曲线也受气温、生产计划影响而微网内部可调机组少、储能容量有限、与主网的联络线功率又常常受限。如果只按预测值做确定性调度实际运行中很容易出现两个问题一是功率不平衡导致切负荷或弃风弃光二是被迫在实时市场高价购电运行成本大幅上升。举个例子某微网按光伏预测出力80kW安排了日前计划但实际光照受云层影响只有30kW缺口50kW不可能瞬间由燃气轮机补上储能也可能已经充满或放空。此时调度员只能切负荷或者以惩罚性电价从主网购电。一次两次还能接受长期下来可靠性和经济性都失控。更麻烦的是微网中的负荷预测也不是精确的用户侧行为、电价响应都会带来偏差。所以调度模型必须把“不确定性”本身建模进去而不是简单地把预测误差当成事后修正。确定性调度之所以在学术和工程中都被认为不够用是因为它只处理一条预测曲线得到的解往往是“在预测场景下最优但稍微偏离一点就不可行”。两阶段鲁棒优化正是为了处理这类问题在日前计划阶段就考虑最坏情况并在日内通过第二阶段变量进行最经济的调整。这样得到的计划可能不是预测场景下成本最低的但一定是面对不确定性时仍然安全、且综合代价可控的方案。1.2 两阶段鲁棒的核心逻辑min-max-min两阶段鲁棒的数学形式看起来很唬人其实思想非常朴素。以出行规划来类比我个人觉得最容易理解周末去一个陌生城市你能提前确定的是机票、酒店、主要景点这是第一阶段决策当地天气可能下雨、可能暴晒、可能交通管制这是不确定性真到了当地你可以临时调整打车路线、改签门票这是第二阶段调整。两阶段鲁棒要求的就是“提前计划不能太死要保证在每一种天气下都有办法应对并且最坏那天的总花费也可以接受”。对应到微网调度第一阶段是日前计划包括机组启停状态、储能的充电/放电安排、与主网购售电的日前计划等。这些变量在知道实际风光和负荷之前就要定下来特点是“定了很难改”。第二阶段是日内再调度在实际风光、负荷场景暴露后根据不平衡量调整机组出力、修正储能功率、必要时切负荷。这样整个优化目标可以写成$$\min_{x} \left( c^T x \max_{u \in U} \min_{y \in \Omega(x,u)} d^T y \right)$$最外层的min对应日前计划成本中间的max对应寻找最坏不确定性场景最内层的min对应在给定场景下做最经济的内日调整。三个层级嵌套就是常说的min-max-min结构。使用鲁棒而不是随机优化的原因也很实际随机优化需要知道不确定量的概率分布可微网历史数据往往不够分布估计误差很大鲁棒优化只需要给出波动范围或者一个不确定集合数据需求低得多保守性也可以通过参数调节。1.3 关键场景辨别算法在整体方案中的位置两阶段鲁棒模型理论上面对的不确定集合是连续的比如光伏出力在一个区间内连续变化。如果直接枚举所有可能的场景计算量无法接受如果不枚举又无法准确评估最坏情况。关键场景辨别算法就是来解决这个矛盾的它不追求一次性考虑全部场景而是在“一个有限场景集合”上反复迭代每次找出当前计划下最恶劣的几个场景把它们加入主问题重新优化直到结果收敛。这里说的“最恶劣”可以从两个角度看一是让系统运行成本最高的场景二是让约束最难满足的场景。在工程实现中通常是先生成大量候选场景然后通过第二阶段子问题逐个评估选出当前调度方案下成本最高的候选场景作为关键场景。整个过程和CCG列与约束生成非常相似主问题在已有关键场景集合上求最优计划子问题负责“辨别”下一个需要加入的关键场景。因为关键场景数量通常很少往往几个到十几个就够所以求解速度远快于枚举所有场景。这也是为什么这个项目明明面对复杂的不确定优化却能在Matlab上用通用求解器跑通的核心原因。2. 数学模型把调度问题写成两阶段鲁棒优化2.1 输入侧的不确定性建模与预算参数微网中的不确定性主要来自风电出力、光伏出力和负荷。实际建模时不能只给一个预测值还要给出预测偏差的上下界。以风电为例t时段出力可以写成$$P_{w}(t) \in [P_{w}^{fore}(t) - \delta_{w}(t), ; P_{w}^{fore}(t) \delta_{w}(t)]$$光伏和负荷同理。把所有不确定量放在一起就得到一个盒式不确定集合U。盒式集合虽然简单但直接使用会带来一个麻烦如果每个时段、每个不确定量都取到边界模型会非常保守实际运营者往往觉得成本高得离谱。为了解决这个问题可以引入不确定预算Γ限制所有时段中最多只能有Γ个时段的偏差达到边界其余时段只能在较小范围内波动。这个思想在鲁棒电力系统调度里很常见本质上是在最坏情况和实际可能之间找一个折中点。Γ值的选择很有讲究。Γ0时所有不确定量都取预测值模型退化为确定性调度ΓT时每个时段都可以达到最坏边界模型最保守。实际项目中我会把Γ当成一个可调的旋钮从0到T逐步扫描结合微网的风险承受能力确定。一般来说微网内部容量小安全裕度要求高我会取Γ在0.6T到0.8T之间既不会把成本抬得太高又能覆盖大部分极端情况。2.2 目标函数与两阶段变量划分整个模型的目标函数是总成本最小化。第一阶段成本包括燃气轮机的燃料成本、启停成本以及日前与主网购电的成本。第二阶段成本是在具体场景下的再调度成本和惩罚成本包括机组出力调整成本、储能充放电的运维折算、切负荷惩罚、弃风弃光惩罚。第一阶段决策变量具体包括燃气轮机启停状态、机组日前出力计划、储能日前充电计划和放电计划、日前购售电计划。这些变量的特点是必须在知道实际风光和负荷之前确定。第二阶段决策变量包括每个场景下机组出力的调整量、储能实际充放电功率的修正量、切负荷量、弃风弃光量。第二阶段变量的自由度更高相当于给系统一个“场景发生后的补救手段”。为什么切负荷和弃风弃光要进第二阶段而不是直接在约束里硬性禁止因为如果不允许切负荷最坏场景下模型可能直接不可行而且硬性禁止会导致第一阶段计划过度保守。把切负荷作为带高惩罚的松弛变量放进目标函数可以让模型在最坏情况下仍然可行同时通过惩罚系数的大小让优化器尽量不去用它。惩罚系数一般设置为正常购电成本的5到10倍具体可以根据项目对供电可靠性的要求调整。2.3 约束条件与松弛惩罚模型的约束大致可以分为四类。第一类是第一阶段的运行约束。燃气轮机出力有上下限和爬坡约束启停状态与出力之间需要逻辑约束储能的充电功率和放电功率不能同时为正SOC要满足递推关系并且限制在安全上下限之间与主网的购售电功率不能超过联络线容量。第二类是功率平衡约束。在每个时段所有电源出力、储能放电、购电之和等于负荷、储能充电、售电之和。在确定性调度里这是一个等式约束在两阶段鲁棒模型里需要分场景写第一阶段不考虑不确定量时有一个“日前预测平衡约束”第二阶段在每个关键场景下都写一个实时功率平衡约束。第三类是第二阶段场景约束。给定一个场景第二阶段调整后的机组出力、储能功率、弃风和切负荷量必须满足功率平衡和所有运行上下限。也就是说无论出现哪个关键场景都存在一套可行的调整方案。这是鲁棒解“可行”的核心保证。第四类是松弛变量的取值范围约束。切负荷量不能超过该时段负荷值弃风弃光量不能超过该时段风光预测出力松弛变量均为非负。在实际代码里这些变量会让模型从“可行或不可行”变成“可行但昂贵”大大提高了求解稳定性。2.4 紧凑形式与求解架构把上面这些约束写成矩阵形式两阶段鲁棒模型可以整理成下面这样的标准形式$$\min_{x} ; c^T x \max_{u \in U} \min_{y} ; d^T y$$$$\text{s.t.} ; A x \le b, ; x \in X$$$$B x C y \le g - D u, ; y \in Y$$主问题是在一个有限关键场景集合K上求解混合整数线性规划MILP。给定N个关键场景后主问题可以写成$$\min_{x, y_1,...,y_N, \eta} ; c^T x \eta$$$$\text{s.t.} ; A x \le b$$$$\eta \ge d^T y_k, ; \forall k$$$$B x C y_k \le g - D u_k, ; \forall k$$这里的η用来逼近所有已选场景中最大的第二阶段成本。由于主问题只考虑了部分场景它的最优目标值只是原问题的一个下界。子问题则在给定第一阶段决策x的情况下搜索最恶劣场景并计算对应的第二阶段成本这个成本又构成了原问题的一个上界。上下界不断靠近就完成了鲁棒迭代求解。3. Matlab实现主问题、子问题与关键场景迭代3.1 代码结构与环境准备在Matlab中实现这套模型我推荐使用YALMIP作为建模层求解器可以选择Gurobi或CPLEX。YALMIP的优势是可以用接近数学表达式的方式写约束对项目前期快速验证特别友好求解器负责真正解MILP和LP。环境准备很简单安装Matlab R2020b以上版本下载YALMIP源码并加入路径再安装Gurobi并配置yalmip的solver选项。装好后运行yalmiptest看到所有测试通过即可。整个项目代码按功能拆成几个文件不建议把所有逻辑写在一个大脚本里。我习惯的文件结构是这样的文件名作用main_two_stage_robust.m主脚本控制迭代流程case_data.m微网参数、预测曲线、不确定度配置generate_candidate_scenarios.m生成候选场景池master_problem.m构建并求解主问题MILPsub_problem.m对给定阶段一方案求第二阶段成本select_key_scenarios.m在候选场景中辨别关键场景固定随机种子这一点非常关键。候选场景生成如果用了随机抽样不固定种子会导致每次运行结果不一样调试的时候会怀疑人生。项目里我会在脚本开头写rng(1234)保证可复现。3.2 主问题建模有限场景集下的MILP主问题的核心是用YALMIP定义变量、约束和目标函数。先定义第一阶段变量以燃气轮机为例启停是二进制变量出力是连续变量P_gt sdpvar(T, 1, full); C_gt binvar(T, 1, full); % 1表示开机 E_s sdpvar(T1, 1, full); % 储能SOC P_ch sdpvar(T, 1, full); % 储能充电功率 P_dis sdpvar(T, 1, full); % 储能放电功率 P_buy sdpvar(T, 1, full); % 从主网购电 P_sell sdpvar(T, 1, full); % 向主网售电 eta sdpvar(1, 1); % 近似最坏场景二阶段成本主问题中还需要为每个关键场景定义一套第二阶段变量。为了不把变量定义得杂乱我通常用一个cell数组来装每个场景的二阶段变量。场景循环里添加约束给定场景u_k第二阶段运行变量必须满足功率平衡和上下限约束并且用目标变量obj_k记录该场景的二阶段成本然后加约束eta obj_k。功率平衡约束示例Constraints []; Constraints [Constraints, P_gt(t) P_dis(t) P_buy(t) P_pv_fore(t) P_wind_fore(t) ... P_load_fore(t) P_ch(t) P_sell(t)];注意第一阶段功率平衡用的是预测值第二阶段功率平衡再用实际场景值形成“预测计划场景调整”的组合。目标函数写成Objective sum(C_fuel .* P_gt) sum(C_start .* max(0, diff(C_gt))) ... sum(C_buy .* P_buy - C_sell .* P_sell) eta;调用optimize(Constraints, Objective, ops)即可。如果Gurobi装好了ops sdpsettings(solver,gurobi,verbose,0)。3.3 子问题与关键场景辨别如何找到最恶劣场景子问题要回答的问题是给定主问题求出来的第一阶段方案x_star在哪个场景下系统的第二阶段成本最高这个场景就是当前的关键场景。在严格的两阶段鲁棒求解中子问题需要在连续不确定集合上寻找最坏场景通常需要通过强对偶转化为一个带双线性项的单层优化再用大M法线性化。这种做法严谨但代码量和调试成本都不小。我在实际工程中更常用一种简化但足够有效的实现预先通过蒙特卡洛或拉丁超立方抽样生成一个较大的候选场景池然后在每次迭代中对候选池里的每个场景都计算一次第二阶段成本取成本最高的那个作为关键场景加入主问题。第二阶段子问题在YALMIP中的建模非常直接。把主问题得到的第一阶段变量value()之后当作常数传入子问题函数function cost eval_scenario(x_star, u_k) y sdpvar(...); % 定义第二阶段调整变量 Constraints [ ... ]; % 固定x_star写入u_k的实际值 Objective ...; % 第二阶段成本加惩罚 optimize(Constraints, Objective, ops); cost value(Objective); end在select_key_scenarios.m里遍历候选池记录最大成本对应的场景编号。如果发现连续两轮选出的都是同一个场景说明关键场景集合已经稳定可以直接结束迭代。这样做的好处是逻辑简单、不容易出数值问题而且候选场景池本身如果构建得好逼近效果已经相当接近严格鲁棒解。严格连续鲁棒的CCG路径可以作为后续扩展方向整体框架并不冲突。3.4 迭代主循环与收敛判据整个程序的核心迭代逻辑如下key_scenarios []; LB -inf; UB inf; for iter 1:max_iter % 1. 求解主问题 [x_star, obj_master] solve_master_problem(key_scenarios); LB max(LB, obj_master); % 2. 遍历候选场景辨别最恶劣场景 [worst_scene, worst_cost] find_worst_scenario(x_star, candidate_scenarios); total_ub obj_master - value(eta) worst_cost; % 更严谨用cx_star worst_cost UB min(UB, total_ub); % 3. 收敛判断 gap (UB - LB) / abs(UB) * 100; if gap tol break; end % 4. 把关键场景加入主问题场景集 key_scenarios(end1) worst_scene; end这里有个容易踩坑的地方UB不能直接用主问题目标函数值因为主问题里eta只是一个下界不是真实最坏成本。必须用当前x_star的“第一阶段实际成本 子问题找到的最坏二阶段成本”来更新UB。我第一次写的时候直接拿master的obj当UB结果上下界越迭代越乱折腾了半天才反应过来。收敛判据我采用相对gap阈值一般设为0.5%。如果候选场景池足够大通常迭代三到八次就能收敛耗时从几十秒到几分钟不等完全在可接受范围内。4. 算例结果与分析4.1 测试微网结构与参数为了验证代码我搭了一个典型小型微网算例。系统包含一台燃气轮机、一组风电、一组光伏、一套储能以及与主网的联络线。调度周期取24小时时间间隔1小时。主要参数如下参数数值燃气轮机容量50 ~ 100 kW燃气轮机燃料成本0.42 元/kWh启动成本2.5 元/次风电容量80 kW光伏容量60 kW风电/光伏预测偏差±20%储能容量60 kWh储能最大充放功率15 kW储能SOC范围0.2 ~ 0.9储能初始SOC0.5购电价格0.6 元/kWh售电价格0.4 元/kWh联络线功率上限50 kW切负荷惩罚5 元/kWh不确定预算Γ14候选场景池生成1000个场景随机种子固定为1234。每个场景包含24时段的风电、光伏和负荷偏差。生成后先用简单指标做一次预筛选比如按净负荷极端程度排序留下500个候选场景再进入鲁棒迭代这样既保留多样性又减少子问题评估次数。4.2 迭代过程与调度结果用上面的参数跑程序迭代过程如下迭代次数下界LB上界UBgap11136.21258.79.7%21219.51250.42.5%31242.81249.60.54%41247.31249.10.14%第4次迭代后gap降到0.14%满足收敛条件程序停止。可以看到第一次迭代的主问题没有加入任何关键场景等价于只按预测场景调度成本很低但子问题一评估就发现最坏场景下总成本要高出120多元说明原计划在极端场景下非常脆弱。随着关键场景逐个加入下界快速上升最终稳定在1249元左右。从关键场景本身来看程序辨识出的第一个关键场景是一个“晚间负荷高峰但风电远低于预测”的组合场景第二个关键场景是“中午光伏骤降且负荷较高”的场景第三个是凌晨风电过剩、必须大量弃风同时储能充满的场景。这几个场景几乎覆盖了微网最典型的三种风险高峰缺电、光伏骤跌、风电反调峰。算法自动找出来的场景和人工经验判断的极端情况高度一致这是让我比较放心的一点。对比确定性调度确定性模型优化结果总成本是1122.8元比鲁棒方案低约10%。但如果把这个确定性计划放到最坏场景下实际上运行会发现功率无法平衡需要额外切负荷或高价购电真实成本会达到1280元以上而且伴随可靠性风险。鲁棒方案用约10%的成本上升换来了所有候选场景下的绝对可行这笔账在微网工程里是划算的。4.3 鲁棒参数灵敏度与工程平衡不确定预算Γ是一个非常有用的调节旋钮。我分别用Γ0、7、14、21、24跑了一遍Γ值鲁棒调度总成本关键场景迭代次数01122.8171181.52141249.14211287.65241301.36Γ从0到24成本从1122元涨到1301元涨幅约16%。这说明鲁棒性不是免费的越保守成本越高。实际项目中我会建议根据微网运营方的风险偏好选择一个“成本-可靠性”的平衡点。如果切负荷惩罚很高或者有重要负荷就选较大的Γ如果成本敏感就选小一些的Γ同时通过储能配置和机组组合来弥补安全裕度不足。储能容量对鲁棒调度结果的影响也值得注意。在同样条件下把储能容量从60kWh提升到90kWh鲁棒最优成本下降了约3.5%关键场景迭代次数也减少到3次。储能相当于给系统提供了更多第二阶段调整手段让最坏场景下的功率平衡更容易实现从而缓解了第一阶段计划的保守程度。这个结论在写报告或者做方案论证时非常有用提升鲁棒性不一定只靠增加不确定预算也可以通过合理配置灵活性资源实现。5. 常见问题与排查技巧实录5.1 主问题求解慢或内存不足主问题是MILP随着关键场景数量增加变量和约束会线性增长。第一次跑的时候我加入了20个关键场景主问题包含超过800个二进制变量和几万个约束Gurobi默认参数下求解时间迅速上升最后甚至出现内存不足。解决办法有几个。第一给Gurobi设置合理的MIP gap比如ops sdpsettings(solver,gurobi,mipgap,0.01)牺牲少量最优性换取速度。第二检查是否有无效变量比如某些场景下弃风量变量恒为0就不要定义。第三如果场景数量真的很大可以一次迭代同时加入多个关键场景而不是每次只加一个虽然单次主问题变重但总迭代次数减少整体反而更快。还有一个容易被忽视的点YALMIP建模时如果循环里拼接约束用了[Constraints, Constraints, new]这种写法会在变量数量大时产生很多中间矩阵拖慢求解器预处理。建议先用空约束初始化然后在循环里用Constraints [Constraints, new_constraint]尽量保持结构紧凑。5.2 子问题不可行或目标值异常子问题不可行多半是第一阶段决策x_star没有给第二阶段留出调整空间。比如燃气轮机已经满发、储能SOC也到了上限但最坏场景下负荷仍然高于所有可用电源之和这时第二阶段无论如何都无法平衡功率。我的处理方法是确保模型中始终保留切负荷和弃风弃光松弛变量。这两个变量几乎能让子问题从“不可行”变成“可行但高成本”所以子问题函数里不要省掉它们。如果加入松弛变量后仍然报infeasible则需要检查第二阶段约束里是否不小心把x_star也当成变量了。在YALMIP中第一阶段变量传入子问题时必须用value(x_star)或直接赋值成数值不能把sdpvar对象传入后再使用的用错否则会建立一个本不存在的联合优化问题。另外子问题目标出现NaN或负无穷通常是因为目标函数中某个成本系数写成了0或者整数变量参与计算。逐项检查value()结果尤其注意储能充放电同时不为0的约束是否写对。如果储能充电和放电可以同时为正功率平衡中的“充放电抵消”会让目标函数出现奇怪的退化现象。5.3 迭代振荡或不收敛上下界不收敛的原因很大概率是LB或UB更新逻辑写反了。主问题的obj是下界这一点很多人会记成上界。实际上主问题只考虑了部分场景相当于一个松弛问题目标值只会偏低所以应该是LB。相反用x_star和子问题最坏场景算出来的总成本才是这个可行方案在“至少一个场景”下的真实成本可以作为UB。如果LB和UB都正常但gap始终停在5%左右不动可能是候选场景池里根本没有真正恶劣的场景。比如你只用了预测值附近小范围抽样或者随机种子刚好抽出一批温和场景。这时候程序会很快收敛到一个“看似鲁棒”但实际不够安全的解。我的经验是先生成一批较大的候选池然后故意加入几个边界场景比如光伏出力全天偏低、负荷全天偏高用来“逼”鲁棒解更稳妥。关键场景辨别算法的价值恰恰在于即使候选池很大真正进入主问题的也只会有那么几个所以没必要担心边界场景让求解变慢。5.4 求解器环境与YALMIP相关坑YALMIPGurobi的组合在学术用户里很流行但配置过程仍然有坑。最常见的是Gurobi许可证没有正确激活导致optimize直接报错。建议安装后用命令行执行gurobi_cl --license确认许可证状态。另一个常见问题是Matlab版本和YALMIP版本不匹配表现为有些函数找不到或者求解器识别失败。遇到这类问题升级YALMIP到最新版基本能解决。求解器选择也需要注意。主问题MILP用Gurobi或CPLEX都没问题子问题是LP也可以直接用同一套求解器。如果机器上没有安装商业求解器可以用YALMIP内置的sedumi或glpk做小规模测试但大算例会很慢不建议生产环境依赖。我在项目里通常都是先装好Gurobi再开始建模避免中途切换求解器导致约束写法不兼容。5.5 一定要记住的四个实操细节第一单位统一。Matlab代码里很可能把kW、kWh、元混在一起一旦出现1000倍误差结果会完全失真。我习惯在所有数据文件里统一采用“kW-元-小时”单位并在参数表注释里写清楚。第二固定随机种子。候选场景生成如果不设rng每次运行结果都不一样很难复现论文数据。正式跑批量实验时我会专门写一个run_experiment.m循环里每次设置不同种子但通过随机种子列表记录下来保证实验可复现。第三不要一开始就上大规模。先用24小时、一台机组、5个候选场景的小例子验证主程序和子程序逻辑确认LB/UB变化趋势正确后再扩大到完整算例。否则变量一多模型出问题根本无从排查。第四保留每个迭代轮次的关键场景编号和场景数据。优化完成后把这些场景的负荷曲线和风电光伏曲线画出来对比你会发现它们都是非常有物理意义的极端情形。这些图既可以用在论文里也能帮助解释为什么鲁棒解比确定性解成本更高让结果更有说服力。最后再分享一个小技巧调试两阶段鲁棒程序时最有效的办法是打印每一轮的关键场景编号、LB和UB。如果看到LB持续单调上升、UB持续单调下降哪怕gap还大程序大概率是对的如果出现LB下降或者UB上升先别急着改模型回到主问题和子问题的目标函数定义去检查绝大多数问题都出在上下界更新或者变量固定上。这个思路帮我省了非常多排查时间也推荐给你。