1. 从单层优化到双层优化:多出来的那层到底在哪
如果你在微网调度这个圈子里待过一阵子,大概率先接触的是单层优化:给定负荷曲线,给定从上级电网买电的分时电价,然后让储能、光伏、可控机组老老实实地算出一个成本最低的调度方案。这套模型的好处是干净、好解、结果稳定,几乎是所有入门教材的标配。但真到了实际项目里,只要涉及"动态电价"或者"需求响应",单层模型就会露出一个非常尴尬的破绽——电价和负荷是互相咬合的,你没法把一头固定住去优化另一头。
我做的这个微网场景里,运营商手里握着一套储能,还带着一个可以跟上级电网交易的联络线。如果我只按单层模型去优化,那就等于默认用户在任何电价下都会保持同一个用电行为。可现实是用户对电价有反应:电价高了,他们削减用电或者把负荷挪到别的时段;电价低了,他们可能多开空调、多充电。你上层把电价当作外生参数,下层又把负荷当作外生参数,这个循环根本闭环不了。真正的问题是:运营商的定价决定用户的用电量,用户的用电量又反过来决定运营商的最优调度,这是一个典型的博弈过程,而不是一个单纯的最优化过程。
这就是双层优化模型存在的理由。学术界一般叫它领导者-跟随者结构,也就是斯坦克尔伯格博弈。微网运营商作为领导者,先抛出自己确定的售电价;用户作为跟随者,在看到电价之后,按照自己的用电效用去决定买多少电。运营商在决策时不能强迫用户接受任何用电量,只能通过调整电价来间接"遥控"用户的反应。这个关系写成数学形式是这样的:
min F(x, y) s.t. G(x, y) ≤ 0 y ∈ argmin { f(x, y) : g(x, y) ≤ 0 }
看起来只是在上层约束里多了一个 y 的 argmin,但这个 argmin 才是整个问题的灵魂。它代表用户拥有自己的优化目标,上层不能替用户做主。我第一次接触这个式子时很不以为然,觉得无非就是把下层问题的解用一个函数带进去,结果查了半天文献才发现,这件事远比想象中麻烦。下层对不同的 x 会有不同的最优反应 y(x),这个反应函数往往是非线性的,你根本没办法直接把 y 显式解出来塞回上层。必须换一种思路处理"下层最优性"这个约束。
我建议任何一个准备做微网双层优化的同学,先别急着写代码,第一步就是反复确认自己到底能不能讲清楚:上层有哪些决策变量、下层有哪些决策变量、哪几个变量是两端共用的耦合变量。如果这三样说不清楚,后面所有的MATLAB代码都只是在堆砌错误的约束。
2. 上下层模型怎么落笔:变量、目标、约束的边界划分
我的算例规模刻意控制得比较小:24个时段,3个用户聚合体,1套储能,外加一条上级电网联络线。为什么要这么小?因为双层优化模型从建立到验证需要反复修改,算例太大每次调试都在等求解器跑半天,心态很容易崩。等模型验证可信了,再把用户数从3个扩到30个、时段从24变成96都不难。小算例不是偷懒,是为了留出快速试错的余地。
2.1 下层用户问题:购电费与用电效用的权衡
下层用户的目标说白了就是一句话:在电价确定的前提下,找到让自己最舒服的用电量。但"最舒服"不能拍脑袋定义,我采用的是二次凹效用函数:
U_i(Q) = a_i * Q - 0.5 * b_i * Q²
其中 Q 是用户购电量,a_i 反映用户对电力的基础依赖程度,b_i 反映边际效用递减的速度。用户实际面对的问题是:
min π * Q - U_i(Q) s.t. Qmin_i ≤ Q ≤ Qmax_i
这里 π 是运营商设定的售电价,Qmin 和 Qmax 分别是需求响应协议允许的最小、最大用电量。减去效用函数是因为用户要同时权衡"用电带来的满足感"和"需要付出的电费",只有二者的差值最小化时才是用户的利益最大点。
这个模型最重要的数学性质是目标函数为严格凸二次函数,约束条件只有简单的区间约束,因此下层是一个严格凸二次规划。严格凸意味着KKT条件不只是必要条件,它同时也是充分条件,KKT点就是全局唯一最优解。这个性质在后面做KKT转换时起了决定性作用,如果下层问题不够凸,后面所有推导都会站不住脚。
2.2 上层运营商问题:定价、储能调度与配网购电
上层运营商的决策变量有三组:各时段的售电价 π(k)、储能净放电功率 P_ess(k)、向上一级电网的购电功率 P_grid(k)。P_ess 为正时表示放电,为负时表示充电。运营商的日净收益表达式为:
max Σ [ π(k) * Q_total(k) - c_buy(k) * P_grid(k) - c_ess * |P_ess(k)| ]
Q_total(k) 是当前时段所有用户购电量之和,也就是微网总负荷;c_buy(k) 是运营商向配电网购电的分时电价;c_ess 是储能单位充放电损耗成本。加这个很小的损耗项,目的不只是经济性,更重要的是避免储能模型在目标函数中"无成本"地反复充放电,生成一些震荡解。
上层约束更多是常规项:功率平衡约束 P_grid + P_ess + P_pv = Q_total,储能SOC递推式 SOC(k+1) = SOC(k) - P_ess(k)*Δt/E_cap,SOC上下限约束,储能功率限值,售电价上下限,联络线购电功率上限。这些约束在单层调度里都见过,真正的难点在第2.3节的耦合关系上。
2.3 上下层耦合点与参数量级预算
上下层的第一个耦合点是售电价 π(k) 出现在下层用户的目标函数里。第二个耦合点是用户总购电量 Q_total(k) 同时出现在上层的功率平衡约束和目标函数中。这两个耦合点会让整个模型形成闭环:上层定价 → 用户响应 → 用户响应结果反哺上层收益与功率平衡。
在设计模型阶段,还有一个容易被忽略的工作是提前估算KKT线性化需要的M值。电价量级大约在0.3到1.5元/千瓦时,用户购电量通常在几十千瓦的量级,效用系数 a_i 为了让用户目标函数在量级上合理,也要取到几元/千瓦时,b_i 则一般在10的负二次方到负三次方之间。粗略估算下来,KKT条件中拉格朗日乘子的量级通常在几十以内,预留M取500左右是安全的。这个估算看似随手,实际帮我省掉了大量后期数值调试的时间,我见过太多人上来就M取一万,结果模型数值病态到求解器根本没法收敛。
3. 求解方案选型:为什么KKT转换最终成了我的选择
双层优化不是只有一个解法。我在动手之前梳理了三条最主流的路线,也分别试跑过,最终才确定KKT转换这条路线最稳妥。
3.1 三条常见路线对比
| 求解路线 | 基本思路 | 优点 | 缺点 | 适用场景 |
|---|---|---|---|---|
| 枚举/离散化 | 把上层变量按网格离散,对每个组合求解一次下层 | 实现简单,不需要推导 | 组合爆炸,只能得到近似解 | 上层变量少、精度要求低的演示算例 |
| 智能优化算法 | 外层用遗传/粒子群搜上层变量,内层每次求解下层 | 几乎不受模型非线性限制 | 目标函数不光滑,容易陷入局部最优,且无最优性保证 | 非线性强、精度要求不高的探索 |
| KKT/强对偶转换 | 把下层问题的KKT条件作为约束并向上层合并 | 可转成严格MILP,求解器给出全局最优 | 推导工作量大,互补松弛线性化易出错 | 下层为凸、约束连续的典型微网双层模型 |
我自己最初用枚举法跑过一个1时段1用户的极简算例,效果不错,甚至能画出"电价-收益"曲线。但我知道一旦把上层电价做成24时段连续变量,枚举的复杂度直接指数爆炸。智能算法我也跑过,粒子群搜电价、内层用quadprog求用户响应,问题在于上层目标函数不是光滑函数,粒子群很容易在一个局部鞍点附近反复横跳,且完全无法证明解的最优性。
最终我选择了KKT转换路线。理由有三:第一,下层问题是严格凸二次规划,KKT条件完全等价于下层原问题;第二,互补松弛经过大M线性化后可以转化为MILP,成熟的商业求解器能给出带最优性保证的全局解;第三,只要推导一次,后续更换场景、增加约束都只需要在模型里追加条件,不需要换求解框架。
3.2 KKT转换的完整推导与线性化逻辑
对每个下层用户 i、每个时段 k,写下拉格朗日函数并对 Q 求导,得到稳定性条件:
a_i - b_i * Q_i(k) - π(k) - λ_i_ub(k) + λ_i_lb(k) = 0
λ_i_ub 对应上限约束 Q_i ≤ Qmax_i 的乘子,λ_i_lb 对应下限约束 Qmin_i ≤ Q_i 的乘子。与此同时还要满足互补松弛条件:
λ_i_ub * (Qmax_i - Q_i(k)) = 0 λ_i_lb * (Q_i(k) - Qmin_i) = 0
这两个等式里带乘积项,是非线性约束,不能直接扔给MILP求解器。需要用大M法引入0-1变量 z,把它拆成两组线性约束:
λ_i_ub(k) ≤ M * z_ub,Qmax_i - Q_i(k) ≤ M * (1 - z_ub) λ_i_lb(k) ≤ M * z_lb,Q_i(k) - Qmin_i ≤ M * (1 - z_lb)
这组约束等价于两句话:要么乘子非零且约束取等号,要么约束放松且乘子为零。换句话说,z 是"该约束是否被激活"的开关。很多初学者把第二行写成 Q 减去上下限的方向,结果方向一错,模型不报错但解出来的互补松弛完全不成立。这种事真的需要一个一个小算例去验证,不能指望求解器帮你发现逻辑错误。
3.3 为什么我不建议直接用现成的bilevel求解器
很多人看到YALMIP自带bilevel solver接口就很兴奋,我一开始也一样。但实测发现这个接口限制太多:下层必须能写成特定的二次规划形式,约束结构稍有变化就报错,而且报错信息对调试没有任何帮助。更关键的是,使用现成bilevel求解器意味着你把模型的推导过程外包给了黑箱,一旦结果不合理,你甚至不知道是模型错了还是求解器内部出了问题。
自己手动做KKT转换虽然前期费些时间,但每一步都能掌握清楚,后续构建不同场景时也能灵活修改约束条件。如果你实在想验证自己推导是否正确,可以先用一个极小规模的算例,把用KKT转化后的MILP解和暴力枚举解做对比,一致了再放心放大规模。
4. Matlab代码骨架:从零把一个双层模型搭成可求解的MILP
到代码部分,我不打算把整个工程几百行全部贴出来,那样太长且针对性太强。我更想把骨架结构讲清楚,尤其是那些容易写错的细节。整体文件组织可以分成三块。
4.1 三个文件的结构与YALMIP建模优势
- case_data.m:集中定义所有参数,包括时段数、用户数、储能容量与功率限值、分时购电价、光伏出力序列、用户效用系数、用电上下限,以及上一步估算好的M值。
- build_model.m:核心建模函数,输入参数结构体 params,输出YALMIP的决策变量、约束集合和优化目标。
- main.m:主脚本,依次执行载入参数、构建模型、配置求解器、求解、画图。
我之所以坚持用YALMIP而不是直接用原生求解器接口,是因为YALMIP允许你以矩阵形式声明变量和约束,建模效率高,而且以后想从Gurobi切换到CPLEX或者其他求解器,只需要改一处求解器设置。唯一要付出的代价是YALMIP在建模预处理时有一点性能开销,对中小规模微网算例来说完全无感。
4.2 核心代码:KKT条件、储能SOC与目标函数
build_model.m 里最核心的变量声明部分是这样的:
function m = build_model(p) N = p.N; I = p.I; pi_sell = sdpvar(N, 1); % 上层变量:各时段售电价 P_ess = sdpvar(N, 1); % 上层变量:储能净放电功率 P_grid = sdpvar(N, 1); % 上层变量:向上级购电功率 Q = sdpvar(N, I); % 下层用户购电量,也是上下层耦合变量 lam_u = sdpvar(N, I); % 上限约束乘子 lam_l = sdpvar(N, I); % 下限约束乘子 z_u = binvar(N, I); % 上限互补松弛的0-1指示 z_l = binvar(N, I); % 下限互补松弛的0-1指示然后写KKT稳定性条件和互补松弛约束:
cons = []; for k = 1:N for i = 1:I cons = [cons, p.a(i) - p.b(i)*Q(k,i) - pi_sell(k) ... - lam_u(k,i) + lam_l(k,i) == 0]; cons = [cons, lam_u(k,i) >= 0, p.qmax(i) - Q(k,i) >= 0, ... lam_l(k,i) >= 0, Q(k,i) - p.qmin(i) >= 0]; % 互补松弛线性化:λ_u>0 当且仅当 Q取到上界 cons = [cons, lam_u(k,i) <= p.M(k,i) * z_u(k,i), ... p.qmax(i) - Q(k,i) <= p.M(k,i) * (1 - z_u(k,i))]; % 互补松弛线性化:λ_l>0 当且仅当 Q取到下界 cons = [cons, lam_l(k,i) <= p.M(k,i) * z_l(k,i), ... Q(k,i) - p.qmin(i) <= p.M(k,i) * (1 - z_l(k,i))]; end end注意上面 p.qmax(i) - Q(k,i) >= 0 这一行约束是用户问题本身的上限约束,不能省略。有些简化版本只写KKT稳定性条件和互补松弛约束,把原始可行性约束漏掉了,模型照样能出解,但解出的Q可能根本不在用户允许范围内。
接着写储能SOC递推和功率平衡:
SOC = sdpvar(N, 1); cons = [cons, SOC(1) == p.soc0 - P_ess(1)*p.dt/p.ess_cap]; for k = 2:N cons = [cons, SOC(k) == SOC(k-1) - P_ess(k)*p.dt/p.ess_cap]; end cons = [cons, p.soc_min <= SOC <= p.soc_max, ... -p.ess_pmax <= P_ess <= p.ess_pmax]; Qsum = sum(Q, 2); cons = [cons, P_grid + P_ess + p.p_pv == Qsum]; % 功率平衡 cons = [cons, p.pi_lb <= pi_sell <= p.pi_ub]; % 电价限值 cons = [cons, 0 <= P_grid <= p.grid_pmax]; % 联络线功率限值储能递推式这个方向很多人会搞混:P_ess 为正表示放电,SOC 应该下降,所以是减去 P_ess*Δt/E_cap;如果写成加号,储能就会一边放电一边SOC上升,模型结果看起来风平浪静,实际上物理过程完全拧了。
目标函数方面,YALMIP默认最小化,所以上层收益最大化要写成最小化负收益:
ess_abs = sdpvar(N, 1); % 储能损耗的绝对值线性化 cons = [cons, ess_abs >= P_ess, ess_abs >= -P_ess]; m.cons = cons; m.obj = -( sum(pi_sell .* Qsum) ... - sum(p.c_buy .* P_grid) ... - p.c_ess * sum(ess_abs) ); m.pi_sell = pi_sell; m.Q = Q; m.P_ess = P_ess; m.P_grid = P_grid; end4.3 主脚本与求解器参数设置
main.m 的调用逻辑非常直接:
case_data; m = build_model(params); ops = sdpsettings('solver', 'gurobi', 'mipgap', 0.01, 'verbose', 2); result = optimize(m.cons, m.obj, ops); if result.problem ~= 0 error('求解失败: %s', result.info); end这里有一个很关键的求解器经验:MILP的默认MIPGap通常设到1e-4,但对带大M约束的模型,这个值会让求解器花费大量时间去证明一个本身就带数值噪声的整数解达到极小gap。我实际测下来,MIPGap放宽到1e-2,求解时间可以下降一个量级,而结果质量几乎不受影响。如果你机器上没装Gurobi,YALMIP也能调用MATLAB自带的intlinprog,只是算例规模稍大时速度会明显变慢,建议还是装一个免费的学术版Gurobi或CPLEX。
5. 调试期中让我反复折腾的五个细节
从模型第一次跑通到结果真正可信,我花了将近三周。这中间踩过的坑,每一个都值得单独写一段。
5.1 大M的选择与数值病态
我第一次图省事,把所有M统一设成10000。结果是Gurobi报numerical difficulties,解出来的值看起来合理,但对账时发现互补松弛条件被数值误差彻底穿透,用户的最优性条件实际上不再成立。后来我把M改成按用户和时段分别估算,总体控制在500的量级,模型立刻稳定下来。M取值没有通解,核心思路是先通过线性松弛求一次乘子范围的估计值,再留出3到5倍的裕量。这个方法虽然笨,但非常有效。
5.2 互补松弛线性化的"或逻辑"不能写反
线性化时很多初学者会写成 λ ≤ Mz 且 g ≤ Mz,认为"g不为0时z必须为1",但这是错的。正确的配对必须是 λ ≤ Mz 与 h ≤ M(1-z),规定z=1时h被压到0,z=0时λ被压到0。这个方向问题模型不会报错,但会静悄悄地破坏KKT条件。我排查这个问题的方式是搭一个两时段的微算例,把MILP结果和暴力枚举结果逐步对比,最后才定位到某一行约束方向写反。
5.3 下层目标不够凸,KKT转换就失去严谨性
一开始我把用户效用写成线性函数 U_i = a_iQ,下层问题退化成线性规划。线性规划的最优解不一定唯一,这就意味着上层在KKT点处无法确定下层用户到底选取了哪一个最优响应,模型的整体结果不可复现。我的解决办法是在用户目标函数里加一个很小的正定二次项 0.5ε*Q²,ε取1e-4。它不影响实际最优解位置,但能让下层问题重新变回严格凸优化,KKT转换的充要性也就恢复了。
5.4 量纲换算会悄无声息地毁掉结果
我有一版模型用的是96时段、每时段15分钟,结果收益数值比预期大出4倍。排查到最后才发现,所有费用项都忘了乘时间段长度Δt。在24时段、Δt=1小时的情况下,这个bug完全隐身;一改成15分钟时段,收益和SOC递推式同时出错。从那以后,我把Δt单独定义在参数文件里,所有费用项统一写成功率×Δt×电价,绝不在公式里写裸的功率乘电价。这个习惯直接杜绝了同类问题。
5.5 用暴力网格搜索做小型对照实验
每个双层优化模型我都会保留一个"对照组":取1个时段、1个用户,把电价从0.3到1.5按步长0.01穷举一遍,每个电价用quadprog求用户的最优Q,再算上层收益。画出"电价-收益"曲线后,把MILP解出来的结果点投到曲线上看是否吻合。这个方法成本极低,却能在十分钟内发现KKT推导、方向或M值中的大部分问题。我强烈建议大家都保留这个验证步骤,它比看任何求解器日志都直观。
6. 结果曲线里的博弈特征,以及我总结的三层进阶
第一批结果图出来后,最有意思的现象是售电价曲线并不是简单地贴着购电价上边界走。高峰时段的售电价会略高于购电价,低谷时段的售电价会被压得很低甚至靠近下限;储能则顺着电价差在低谷充电、高峰放电。部分时段的定价恰好落在"再提一点价,用户就会明显削减用电"的临界点附近,这就是双层博弈均衡在曲线上的具体表现。
对比固定售电价的单层调度,双层模型算出的运营商总收益要高一些,而且高出的部分不是靠单纯抬价,而是来自储能套利和负荷削峰填谷带来的整体效率提升。这个结果本身也印证了做双层优化的价值:不考虑用户响应侧的调度,容易高估定价的自由度,运营商以为价格越贵越好,实际在需求弹性面前,最优价格往往是那个"让用户愿意买又不至于买太多"的平衡点。
回看整条实现路线,我把自己的经验归纳成三层进阶。第一层是能说清楚上下层各自的决策变量、收益边界和耦合变量;第二层是能独立完成KKT推导,并理解每个线性化约束背后的互补语义;第三层是能设计小规模对照实验,验证模型结果在逻辑上和数值上都是可信的。这三层只要有一层含糊,后面的Matlab实现就很容易变成"看起来很能跑、实际结果经不起推敲"的假仿真。
最后再分享一个让我受益很久的工程习惯:所有参数从第一天起就用结构体绑在一起传参,而不是在脚本里散落一堆全局变量。我用结构体之后,从3个用户改成30个用户、从24时段改成96时段,只需要重写参数文件,build_model函数本身一行都不用动。如果你现在刚开始写自己的双层优化模型,这一个习惯能帮你省下后续无数次重构的功夫。