最近在搭一个计及需求响应和电能交互的多主体综合能源系统调度算例时,又被主从博弈折腾了一轮。这套东西看着高深,说白了就是解决一个现实矛盾:综合能源系统运营商想让收益最大,用户想让用能成本最低,两边都不想被对方完全拿捏。集中式优化在这里行不通,因为用户没有义务服从你制定的用能计划。所以要用主从博弈,让运营商作为领导者先定能源价格,用户作为跟随者再根据价格调整负荷,最终收敛到一个双方都能接受的均衡点。这篇文章我打算用Matlab代码把整个思路落地,把建模、KKT转化、Yalmip实现和调参踩坑完整梳理一遍。正在做综合能源系统优化调度、微电网/园区能源管理相关课题,或者想用博弈论做论文创新点的同学,可以参考。
1. 多主体综合能源系统调度为何需要主从博弈
1.1 从集中式优化到主从博弈:一个实际场景
先说一个我实际遇到的场景。一个园区综合能源系统,里面有风电、光伏、燃气轮机、电锅炉和储能,下面带着商业楼、工厂、居民区好几个用户主体。最开始我按传统思路做集中式优化:采集所有用户负荷数据,全局求解一个最小化总成本或最大化总收益的问题,然后给每个用户下达“最优用电指令”。模型跑得很顺,结果也漂亮,但导师一句话问住我:用户凭什么都听你的?
现实中的用户有独立利益,商业楼希望电费最低但舒适度不能降,工厂希望生产计划不被破坏,居民用户则更在意电价波动和用能习惯。强行统一调度,即使全局最优,也落不了地。这就是多主体综合能源系统调度和传统单主体调度的本质区别:信息私有、目标冲突、决策分散。
主从博弈(Stackelberg game)恰好是刻画这种层级决策关系的工具。它天然有一个“领导者”和多个“跟随者”:领导者先行动,公布策略;跟随者观察到领导者的策略后,再做出最优反应;领导者又会根据跟随者的反应调整策略,直到没人愿意单方面改变。在综合能源系统里,最自然的角色分配就是:综合能源服务商(Integrated Energy Service Operator,简称IESO)当领导者,制定售电和售热价格;多个负荷聚合商(Load Aggregator,简称LA)当跟随者,根据价格决定用能计划和需求响应量。最后得到的均衡解,就是双方在利益博弈下都能接受的结果。
1.2 需求响应与电能交互在博弈中的角色
有了博弈框架,需求响应和电能交互这两个词就不再是孤立的技术名词,而是博弈中真正的“耦合变量”。
需求响应,本质上是给跟随者(用户)松绑,让用户对价格信号产生弹性。如果用户负荷完全刚性,运营商定价定多少用户都得买,那博弈就没有意义。可一旦用户具备可转移负荷、可削减负荷、蓄热式电采暖等柔性资源,价格升高时用户会主动削减或转移用电,价格降低时用户会增加用电。运营商想高价卖电,用户会“用脚投票”,这就是供需双方的相互制约。
电能交互,体现在两个层面。一是运营商与上级电网之间的购售电交互:运营商既可以从上级电网买电,也可以在本地富余时向上级电网卖电,这个交互功率本身就是运营商的重要决策变量。二是运营商与多个用户主体之间的电能交互:用户向运营商购电,在不同价格信号下改变购电量,本质上是一种“主从之间”的功率交互。在我的模型里,所有用户主体的购电需求汇聚到一起,再由运营商统一与上级电网交互,形成一个电能交换的枢纽。多主体的博弈结果,最终都会反映到这些交互功率曲线里。
1.3 本文采用的系统结构与博弈关系
我搭建的系统结构比较典型,也方便用Matlab实现。考虑一个IESO作为上层领导者,服务N个负荷聚合商作为下层跟随者。IESO拥有分布式电源(风电、光伏、燃气轮机)、电锅炉、储能装置,向下层用户出售电和热;IESO同时与上级电网进行电能交互,电价按分时电价执行。每个LA内部包含柔性负荷(可转移、可削减)和蓄热式电采暖负荷,LA可以自主决定从IESO购买多少电、多少热,以及参与需求响应的功率调整量。
博弈流程如下:IESO首先公布每个时段的售电价格和售热价格;LA们看到价格后,各自以用能成本最小化为目标,优化购电计划、购热计划和需求响应量;IESO根据LA的购能计划,核算自己的收益和运行约束,再调整价格;重复这个过程,直到价格和用能计划都不再变化。这个均衡点就是Stackelberg均衡。
如果读者想把这个模型扩展到多个IESO之间通过联络线交互电能,可以在上层增加相邻系统交互功率变量,但求解复杂度会明显上升。下面我讲的求解方法,本质上仍可以扩展到多领导者场景,只是收敛性验证会更麻烦,建议大家先把单领导者多跟随者跑通再说。
2. 上层运营商优化模型:定价、购售电与运行约束
2.1 目标函数:收益最大化
上层IESO的目标是最大化整个调度周期内的净收益。这里有个容易混淆的地方:很多教材喜欢把上层写成“成本最小化”,但成本最小化看不到定价带来的收益弹性。在主从博弈里,领导者要决定价格,所以目标函数写成收益最大化更自然。
上层目标函数包含这几部分:
- 售电收入:所有LA在t时段的购电功率之和乘以IESO制定的售电价格;
- 售热收入:所有LA在t时段的购热功率之和乘以售热价格;
- 向上级电网售电收入:如果IESO本地发电量超过负荷需求,可以把多余的电卖给上级电网;
- 向上级电网购电成本:本地发电不足时,从上级电网买电的费用;
- 燃气轮机的燃料成本、储能和新能源设备的运维成本;
- 对用户需求响应的补偿成本:用户削减或转移负荷,需要给予一定补偿,否则用户没有动力参与。
为了简洁,上层目标可以写成下面的形式:
max Σ_t [ p_t^e * (Σ_i L_i,t^e) + p_t^h * (Σ_i H_i,t^h) + p_t^grid_out * P_t^out - p_t^grid_in * P_t^in - C_gt(P_gt,t) - C_om(...) - C_dr(...) ]
其中 L_i,t^e 和 H_i,t^h 不是上层直接控制的变量,而是下层LA对上层公布的 p_t^e、p_t^h 的最优响应。也就是说,上层目标里的售电收入项 p_t^e * L_i,t^e 同时包含上层变量(价格)和下层变量(购电量),这是一个双线性项,也是后面KKT转化时要重点处理的难点。
2.2 约束条件:功率平衡、交互容量与设备运行范围
上层模型需要满足一系列运行约束,这些约束在Matlab里用Yalmip写起来很直接,但容易被忽略。我列一下我在代码里实际用到的约束:
- 电功率平衡约束:燃气轮机出力 + 风电出力 + 光伏出力 + 储能放电功率 + 上级电网购电功率 = LA总购电功率 + 上级电网售电功率 + 储能充电功率 + 电锅炉耗电功率。这里每个时段的等式平衡都是强约束,求解不可行时第一个检查它。
- 热功率平衡约束:燃气轮机余热回收 + 电锅炉产热 + 蓄热装置放热 = LA总购热功率 + 蓄热装置充热。热网动态可以简化成静态平衡,如果考虑管道延时,模型会变成复杂的时空耦合,建议初学者先忽略。
- 上级电网交互约束:购电和售电不能同时为真,需要引入0-1变量约束 P_t^in * P_t^out = 0,或者用一个大M约束限制。实际计算中,如果不加互斥约束,求解器经常给出既买又卖的荒唐结果。
- 燃气轮机约束:出力上下限、爬坡约束。注意爬坡约束要按分钟或小时折算。
- 储能约束:荷电状态SOC动态、充放电功率上下限、以及始末SOC相等。最后一个约束很多论文不提,但实际上为了保证调度方案的可持续性,必须加。
- 电锅炉约束:输入电功率与输出热功率之间有转换效率,还有容量上下限。
- 价格约束:售电/售热价格有上下限,这既来自监管,也避免领导者随意定价导致下层问题无界。
这些约束都相对成熟,不再展开推导。重点是:上层问题的可行域里,价格变量 p_t^e 和 p_t^h 是决策变量,而不是固定参数,这会让整个问题变成带均衡约束的数学规划,也就是MPEC。MPEC在Matlab里直接写给求解器是解不动的,必须处理下层反应函数。
2.3 为什么说“定价+购售电”是上层博弈的核心抓手
我在调试过程中最大的感受是:上层模型里真正体现“博弈能力”的决策变量,就是价格和与电网交互的功率。机组出力、储能充放电这些,说到底是在响应价格信号和负荷需求,属于执行层。如果只把储能和机组出力当作上层变量、把价格固定,那模型就退化成了单主体的运行优化,主从博弈就名存实亡。
所以大家在复现论文时,一定要确认清楚:上层有没有把价格作为决策变量?如果没有,那只是“考虑需求响应的多主体优化”,而不是严格意义上的主从博弈。审稿人或者答辩老师其实很看重这个逻辑。
3. 下层用户需求响应模型:柔性负荷与用能成本最小化
3.1 可转移负荷、可削减负荷与蓄热式电采暖
下层用户不是铁板一块,我按负荷性质把每个LA内部的用电需求拆成三部分:刚性负荷、可转移负荷、可削减负荷。刚性负荷没有调节空间,就是必须消费的;可转移负荷允许在时间轴上平移,比如工厂的流水线、洗衣房的洗涤任务;可削减负荷允许在高峰时段削减一部分,比如空调温度上调一度、照明亮度下调一档,代价是要给用户补偿。
可转移负荷的建模要点:转移前后总用电量保持守恒。比如某个用户一天需要10 MWh的转移负荷,不能转移到最后总电量变成12 MWh。这个守恒约束用数学表达很简单:
Σ_t P_shift,i,t = E_shift,i
同时每个时段的转移功率要满足上下限,并且只能转移到允许的时间窗内。很多初学者容易漏掉“同一个转移任务不能被拆散”的约束,如果忽略离散档位,纯连续模型也是常用简化。
可削减负荷的建模要点:削减量有上限,通常表示成基线负荷的一定比例;同时要限制削减次数或连续削减时长,否则求解器会让用户在某个时段削减到零,这在现实中是 unacceptable 的。我习惯加一个0-1变量,表示用户在某时段是否参与削减,然后令削减功率小于等于该变量乘以上限值。
蓄热式电采暖是北方园区常见的柔性热负荷。它本质是一个小型储热装置,SOC状态方程和电储能类似,但输入是热功率,并且热负荷有舒适度区间约束,允许室温在某个范围内波动。如果你觉得同时建模三类负荷太复杂,最低限度也要包含可转移和可削减两类,这样需求响应才算真正“计及”了。
3.2 用户用能成本最小化目标
每个LA的决策目标是自己在调度周期内的净用能成本最小。净用能成本=购电成本+购热成本-需求响应补偿收益。注意,需求响应补偿是IESO付给用户的,所以对用户来说是负成本。
min Σ_t [ p_t^e * (P_base,i,t^e + P_shift,i,t + P_cut,i,t) + p_t^h * H_i,t + λ_shift * P_shift,i,t^abs + λ_cut * P_cut,i,t - ρ_t^dr * (P_shift,i,t + P_cut,i,t) ]
这里有个关键点:当IESO公布价格 p_t^e、p_t^h 和补偿价格 ρ_t^dr 后,用户的目标函数是线性的,约束也是线性的,所以下层是一个线性规划问题。线性规划非常重要,因为它保证了后续KKT转化的正确性。如果下层是非线性凸规划,KKT条件也能用,但需要求解器支持二次约束;如果是非凸问题,KKT条件只是必要条件,Stackelberg均衡的存在性就无法保证。所以我在构建下层模型时,会刻意避免引入非线性项,比如避免用“目标函数里含二次成本”这类导致下层非线性的设计。
3.3 需求响应约束的凸性分析与KKT可行性
为什么我反复强调下层要是线性规划?因为主从博弈求解有一种非常漂亮的方式:用下层问题的KKT条件替代下层优化问题。KKT条件成立的前提是下层问题满足约束规格(比如Slater条件)并且是凸的。线性规划天然满足。如果下层模型里出现了双线性项,比如用户购电量乘以用户自己的某个决策变量,那就把凸性破坏了,后面所有KKT推导都失去依据。
在Matlab里验证下层问题凸性很简单:把所有约束打出来看一遍,如果有变量相乘、绝对值、分段函数等,就需要先线性化。我遇到过一个案例,用户侧电动汽车的充电费用包含了“充电量乘以充电效率”的变量相乘,必须把效率归到系数里,避免出现变量相乘。这种细节,写论文时很多人不注意,但是到了编程阶段就一定会暴露出来。
KKT可行性的另一个注意点是:下层变量的非负约束也要写进KKT的互补松弛条件里。Yalmip里定义变量时可以声明为nonnegative,这样在手动写KKT时,要记得为每个非负变量生成对应的互补条件。漏掉一个,就会导致求解出的“均衡”实际上只是单方面最优,下层并没有真正达到最优。
4. 主从博弈求解:KKT转化与线性化处理
4.1 用KKT条件将双层问题化为单层
主从博弈直接求解很难,问题在于上层目标里包含下层的最优反应函数。好在当下层是连续线性规划时,我们可以用它的KKT条件把“最优反应”完整地刻画成一组约束和互补条件。这样双层问题就变成单层MPEC问题。
具体做法:对每个LA的优化模型,写出拉格朗日函数;对每个变量求偏导得到驻点条件;列原始可行性约束;列对偶可行性约束;再叠加互补松弛条件。这些条件一起构成了“用户最优反应”的数学描述。把这组KKT条件替代到上层问题中,就得到一个单层的约束优化问题。
这一步说说容易,手推真的很费时间。以4个LA、24个时段为例,每个LA有约10个决策变量、15个不等式约束,KKT条件会多出来上百条互补约束。全部手推容易出错,我曾经因为漏了一个对偶变量,结果算出来的电价高得离谱,排查了整整两天。
4.2 互补松弛条件的线性化
KKT条件里最麻烦的是互补松弛条件,形式是 λ * g(x) = 0,其中 λ≥0 是对偶变量,g(x)≥0 是原始不等式约束。非线性乘积,没法直接交给MILP求解器。标准做法是用大M法引入二进制变量,把互补条件拆成一堆线性逻辑约束。
以条件 λ * g(x) = 0 为例,等价于:
当 λ > 0 时,g(x) = 0;当 g(x) > 0 时,λ = 0。用大M法写:
g(x) ≤ M * z λ ≤ M * (1 - z) λ ≥ 0, z ∈ {0, 1}
这里z是二进制变量,M是一个足够大的正数。这样就把乘积关系变成了线性约束。需要注意的是,一个不等式约束对应一组互补约束,每个互补约束需要一个二进制变量。如果下层约束多,二进制变量数量会涨得很厉害,所以实际求解时M的值和初始点对求解速度影响巨大。
4.3 强对偶定理消除双线性项
上层目标函数里的 p_t^e * L_t^e 仍然是个双线性项。因为L_t^e是下层变量,它和上层价格相乘,属于两个不同层级变量相乘的项,MILP求解器不接受。
我的处理方法是利用强对偶定理:当下层LP达到最优时,它的目标函数值等于其对偶问题的目标函数值。于是可以把下层目标函数中用价格乘以购电量产生的线性项,替换成对偶问题里对应的对偶变量乘以常数系数的项。这样原始目标中的双线性项就消失了,换来的是对偶变量和常数相乘的线性项。具体推导需要根据下层模型的每个约束逐一配平,网上很多论文里叫“强对偶转化”或“dual rewrite”。
这一步是整个算法里最容易出错的环节。我推荐一个检查技巧:在完成强对偶替换后,把求解结果代入,分别计算下层原问题和对偶问题的目标值,看是否相等。如果差了很多,说明替换公式有误。我在调试时加了一个校准函数,专门做这个检查。
经过上述三步——KKT条件、互补松弛线性化、强对偶消双线性——原双层问题就转化成了一个混合整数线性规划问题,可以直接用Gurobi或Cplex求解。
5. Matlab代码实现框架与核心片段
5.1 主程序流程与数据准备
Matlab环境建议R2020b以上,配Yalmip工具箱和一个商业求解器。我个人用Gurobi 10,Cplex也可以,但Gurobi在MILP求解上通常速度更快。
主程序流程我分成五步:
- 写入基础数据:24个时段、机组参数、储能参数、上级电网分时电价、用户基线负荷、可转移负荷量、可削减比例、价格上下限、需求响应补偿价格。
- 用Yalmip定义所有变量,区分上层变量、下层变量、对偶变量、二进制变量。
- 添加上层约束、下层KKT约束、互补松弛线性化约束和强对偶目标替换。
- 调用求解器求解MILP。
- 提取结果,绘制电价曲线、负荷响应曲线、储能SOC和收益柱状图。
数据准备部分我直接写脚本生成,没有用Excel读取。理由很简单:调参时改数据太频繁,每次改Excel还要重新读表,不如在Matlab脚本里用结构体统一管理。比如:
% 基础时段 T = 24; % 上级电网分时购电价,单位元/kWh grid_buy_price = [0.35*ones(1,8), 0.75*ones(1,7), 1.2*ones(1,5), 0.75*ones(1,4)]; % 用户个数 N = 4; % 可转移负荷总量,kWh,每个用户一个值 E_shift_total = [1200, 2400, 800, 1500]; % 可削减比例上限 cut_ratio_max = [0.15, 0.20, 0.10, 0.25];5.2 Yalmip建模核心代码
定义变量阶段最需要注意的是:价格变量是上层决策变量,但在下层KKT条件中又作为固定参数出现。实际处理时,价格变量仍然用sdpvar声明,因为在KKT中它要参与拉格朗日函数的偏导。用Matlab的大致写法如下:
% 上层变量 pe = sdpvar(1, T); % 售电价 ph = sdpvar(1, T); % 售热价 Pgrid_in = sdpvar(1, T); % 上级电网购电 Pgrid_out = sdpvar(1, T); % 向上级电网售电 Pgt = sdpvar(1, T); % 燃气轮机出力 Pes_c = sdpvar(1, T); % 储能充电 Pes_d = sdpvar(1, T); % 储能放电 SOC = sdpvar(1, T+1); % 储能SOC Heb = sdpvar(1, T); % 电锅炉产热 % 下层变量 Pshift = sdpvar(N, T); % 可转移负荷 Pcut = sdpvar(N, T); % 可削减负荷 Hbuy = sdpvar(N, T); % 用户购热 Heat_store = sdpvar(N, T); % 蓄热式电采暖充放热 % 对偶变量 lambda_eq = sdpvar(N, T); % 用于下层等式约束 mu_ineq = sdpvar(N, T); % 用于下层不等式约束 ...写约束和KKT时,可以封装成函数,把对偶变量和原始变量传进去,返回约束数组。这样主函数清爽,排查问题时也方便。
5.3 求解器配置与结果输出
求解前用Yalmip的ops设置求解器:
ops = sdpsettings('solver', 'gurobi', 'verbose', 2); ops.gurobi.MIPGap = 0.001; ops.gurobi.TimeLimit = 600; ops.gurobi.NumericFocus = 3;NumericFocus = 3这个选项非常重要,因为经过KKT和大M转化后,模型数值尺度差异很大,不同数量级变量混在一起,默认参数容易让求解器在单纯形迭代里出现数值困难。我一开始没设置这个,Gurobi经常报“Numerical trouble”警告,设成3之后明显稳定。
结果输出部分,我一般会把决策结果保存到结构体里,然后用plot画四张图:电价与负荷响应曲线、各类电源出力柱状堆叠图、储能SOC和交互功率、各用户成本占比。如果有对比场景,就画柱状图对比收益和成本,方便写报告。
6. 算例验证与结果分析
6.1 参数设置与场景设计
为了让结果能说明问题,我设计了三个对比场景:
- 场景A:无需求响应、无电能交互(用户负荷固定,IESO只能向上级电网购电满足刚性负荷)。
- 场景B:仅考虑需求响应(用户具备可转移/可削减负荷,但不允许IESO向上级电网售电)。
- 场景C:同时考虑需求响应和电能交互(即本文的完整主从博弈模型)。
基础参数如下:24个时段,4个LA;上级电网购电价采用峰谷平三段,峰段1.2元/kWh、平段0.75元/kWh、谷段0.35元/kWh;IESO售电价上下限为0.3~1.5元/kWh;售热价上下限为0.2~0.8元/kWh;燃气轮机容量2 MW,效率0.4,天然气热值按9.7 kWh/m³计算;储能容量1 MWh,充放电效率0.95,最大充放电功率0.25 MW;电锅炉效率0.9,容量1.5 MW;用户可转移负荷比例约15%~25%,可削减比例约10%~25%。
6.2 博弈均衡结果分析
以场景C为例,求解结束后得到的售电价曲线很有意思:高峰时段电价明显抬升,但并不是直接顶到上限1.5元/kWh,而是停留在1.32元/kWh附近。原因是如果电价定得太高,用户的可转移负荷会大规模转移出高峰时段,虽然单价高,但销售量下降,总收益反而变小;如果定得太低,高峰时段上级电网购电成本高,IESO利润空间不够。均衡电价正好是两者博弈后的折中。低谷时段电价下探到0.38元/kWh,吸引用户把可转移负荷转移到夜间,既提高了用户侧的用能经济性,也提高了IESO低谷时段风电的消纳能力。
用户侧的需求响应行为也比较直观:高峰时段,4个LA的可削减负荷平均削减了基线负荷的18%左右;可转移负荷中约70%转移到了22:00到次日6:00的谷段。蓄热式电采暖在谷段自动充热,早高峰时段放热,热负荷的峰谷差显著减小。
IESO与上级电网的电能交互上,购电功率主要出现在16:00到21:00的峰段,而凌晨2:00到6:00之间,由于燃气轮机和风电出力充裕,IESO会选择向上级电网售电。这里的“交互”不是固定的,而是根据价格信号动态调整的,充分体现了“电能交互”在博弈中的作用。
6.3 需求响应和电能交互带来的效益对比
三个场景的优化结果如下表:
| 指标 | 场景A(无需求响应/无交互) | 场景B(仅需求响应) | 场景C(需求响应+电能交互) |
|---|---|---|---|
| IESO总收益(元) | 25600 | 29700 | 34150 |
| 用户总用能成本(元) | 20100 | 17600 | 15700 |
| 系统峰时购电功率(MW) | 2.41 | 1.93 | 1.67 |
| 系统总购电成本(元) | 12600 | 11400 | 9800 |
从表里能看出,场景C相比场景A,IESO收益提高了33%,用户成本降低了22%,峰时购电功率下降了30%。需求响应让用户把负荷从高峰挪到低谷,电能交互让IESO能利用地区间的价差灵活买卖电,两个机制叠加,实现了双赢。
有一点必须提醒:这些数字是特定参数下的结果,不代表所有场景下都有这么大的提升。如果把上级电网的峰谷价差调小,电能交互的收益就会变小;如果把用户可转移负荷比例调低,需求响应的作用也会降低。所以论文里做敏感性分析是有必要的,这能验证模型在不同环境下的适应性。
7. 实际代码调参中的踩坑记录
7.1 迭代法 vs KKT单层法:为什么我放弃迭代
一开始我图省事,用“交替迭代”实现主从博弈:上层先给定初始价格,求解下层得到负荷响应;再把负荷响应代回上层,用启发式规则更新价格;重复直到收敛。结果发现价格经常在上下限之间来回振荡,怎么调步长都压不住,而且更新规则没有理论保证,审稿人一问“怎么证明收敛”我就没法回答。
后来老老实实改用KKT单层转化。虽然推导和编码量大了不少,但求解过程是确定性的,Gurobi给出的是全局最优解(在给定离散变量取值下),不需要担心振荡问题。如果你的目标只是快速验证想法,迭代法可以跑着看趋势,但要做严谨的科研或者实际项目,KKT转化是更稳妥的路线。
7.2 大M值选取和数值病态
互补松弛线性化里的M值选择,是影响求解稳定性的最直接因素。M取得太小,可能把本该可行的解错误截断;M取得太大,会让约束矩阵条件数变差,求解器陷入数值困难。我踩过的最夸张一次,M取了1e8,结果Gurobi直接警告“Matrix is ill-conditioned”,求解出的电价出现很多肉眼可见的不合理尖峰。
我的经验是:M值根据对应约束中变量的数量级来取,通常取该约束可能达到的最大值的2~3倍。比如功率平衡约束右侧最大负荷接近3 MW,M取10就够;对偶变量的M值可以按对偶变量理论上限取,比如1000。不要图省事全模型统一用一个很大的M。必要的话,可以先解一个RELAXATION版本,统计原始约束的最大mismatch,再据此定M。
7.3 不可行问题的系统排查经验
加了KKT和互补条件后,模型规模激增,不可行几乎是必然遇到的情况。我排查不可行时有一套固定流程:
先去掉所有互补松弛约束,求解一个松弛问题,看看原始约束和目标函数是否可行。如果这时候还不可行,说明问题出在功率平衡、设备容量这些基础约束上,检查和上层、下层约束定义是否冲突。如果松弛问题可行,再逐一加入互补约束,定位是哪一组互补约束导致不可行。通常是大M值太小,或者是某个变量忘记定义非负。
另外,Yalmip默认的求解信息里会给出“infeasible”提示,但不会直接告诉你哪条约束坏了。我习惯用yalmip('clear')之后,把约束数组拆分成小块逐一验证,配合check(Constraint)查看残差。这个方法比较笨,但确实能救命。
7.4 关于Yalmip的kkt命令
Yalmip自带一个kkt函数,可以自动生成大规模优化问题的KKT条件。听起来很方便,但实测有几个限制:它要求问题必须是凸的,而且对变量类型和约束形式有比较严格的要求;遇到含二进制变量的上层问题时,直接把kkt用在目标里会报错;生成的KKT条件会比较冗余,导致求解性能下降。另外,kkt生成的互补条件也是非线性,仍需自己处理。所以我的建议是:模型简单时可以用kkt快速验证思路,但正式项目还是手动推导更可控。手动推导的过程虽然痛苦,但能让你对模型结构理解更深,后面扩展模型时也不至于一脸懵。
以上这些经验,都是我从最初能跑通简单算例到最终实现完整主从博弈模型的过程中一步步攒下来的。代码不是一次写成的,模型也不是一步到位的。如果你现在正在折腾双层优化、多主体博弈,建议先把一个最小的两主体、两时段模型跑通,再慢慢扩展到24时段、多用户。小模型上的KKT推导正确了,大模型只是“复制粘贴+改索引”的事,但前提是你真的理解了每一行KKT条件是从哪条约束来的。