搞居民负荷分层调度那会儿，我一度被非合作博弈和双层优化搞得头大。电网侧想通过分时电价引导用户错峰，用户侧又不愿意被强行控制用电习惯，双方目标冲突，这种关系用单层优化根本说不清。后来我把模型拆成上下两层，再用双层鲸鱼算法（Bi-level WOA）去求解，整套东西在Matlab里跑通了。这篇文章就记录一下我在建模、编码、调参过程中摸出来的经验，给同样在做智能用电、虚拟电厂、负荷聚合的朋友做个参考。内容会尽量讲清“为什么这么做”，而不是只给一段能跑的代码。
1. 问题背景与模型思路
1.1 居民负荷调度为什么要“分层”
居民负荷最大的特点是分散、随机、单个设备可控性弱。空调、热水器、电动汽车、洗衣机这些东西，你没法像工厂负荷那样直接下发调度指令。强行控制用户开关，体验很差，也不符合电力市场的逻辑。所以主流做法是用价格信号来“诱导”——你电价高的时候用户自然少用，电价低的时候用户愿意把可转移负荷挪过去。这个过程里，电价就是上层决策，用户响应就是下层行为。
分层结构是天然存在的。电网公司或负荷聚合商站在全局视角，考虑批发购电成本、零售收入、峰谷差、网络容量，他们是领导者；居民用户站在自家电表和舒适度视角，考虑用电费用和生活方式，他们是跟随者。两个层面的目标函数、决策变量、约束条件完全不一样。如果硬塞进一个单层优化问题里，等于强迫中央控制器同时决定每一户的每一个设备应该怎么开，这在现实中根本执行不下去。而且用户设备信息属于隐私，不可能全部上报给调度中心。
分层模型的好处就在这：上层只下达价格、激励等公共信号，下层根据信号自主优化，最后只把“总用电曲线”反馈给上层。既保护隐私，又保留了用户的选择权。数学上，这个结构就是典型的双层优化：上层问题以“下层问题的最优解”作为自己的约束或输入，下层问题在上层给定参数下独立求解。
1.2 非合作博弈与Stackelberg均衡
非合作博弈听起来很学术，拆开看就是每个参与者只干对自己有利的事，不搞合作“协议”。在这个模型里，上层和下层没有结盟关系：上层想多赚钱少风险，下层想少花钱保舒适，目标天然打架。先动的是上层，它宣布电价策略；后动的是下层，它看到电价后再决定怎么用电。这种一先一后的博弈叫Stackelberg博弈，也叫“领导者-跟随者”博弈。
为什么最后能稳定下来？因为存在一个Stackelberg均衡：上层在预测到下层会如何最优响应的前提下，选择自己的最优策略；下层在给定上层策略后，也确实做到了自己的最优反应。在这个均衡点上，任何单方面改变自己策略的一方都不会变好。把这种博弈写成数学形式，就是一个双层优化。所以你可以这样记：非合作博弈是经济学的描述，双层优化是数学的框架，而双层鲸鱼算法就是我们用来解这个框架的工具。
有的文献会把多个用户也建模成互相博弈的多个参与者，用户之间也存在竞争关系，比如都在同一时段充电会导致负荷过载。但大多数时候，如果把每个用户都当成独立博弈主角，模型规模和求解复杂度会爆炸。比较实用的做法是把居民用户看作“价格接受者”——他们不对电价产生策略性影响，只对给定的电价做最优响应。这样做既能保留上下层的非合作博弈特征，又不至于让模型复杂到没法跑。
2. 双层鲸鱼算法的设计思路与Matlab实现骨架
2.1 鲸鱼优化算法的核心机制回顾
鲸鱼优化算法（WOA）是2016年提出的一种元启发式算法，灵感来自座头鲸用“气泡网”捕食。想象一头鲸鱼围绕鱼群收缩包围圈，同时螺旋向上吐泡泡把鱼群赶到中间，时不时还会随机冲向某个方向寻找新鱼群。把这三个行为翻译成数学更新公式，就成了WOA。
包围猎物：
X_new = X* − A · D
其中 D = | C · X* − X |，A = 2a·r − a，C = 2r。a从2线性降到0，r是[0,1]随机数。当|A| < 1时收缩包围，当|A| > 1时跳出当前包围圈去做随机探索。
螺旋气泡攻击：
X_new = D′ · e^(bl) · cos(2πl) + X*
D′是当前个体到最优个体之间的距离，b是螺旋形状常数，l是[-1,1]随机数。
标准WOA通过一个随机概率p来选择收缩包围还是螺旋更新，以50%概率切换。这算法结构简单、参数少——本质上就一个a需要线性递减，其余都是随机数。它的全局搜索能力不错，尤其在变量维度不高（比如24维分时电价）、目标函数非凸的情况下，比传统梯度法好用很多，也容易嵌套成双层结构。
2.2 双层嵌套的实现流程与伪代码
双层WOA的思想非常直观：外层一个WOA搜上层决策变量（比如24小时分时电价），内层再跑一个WOA负责求解给定电价下用户的最优负荷曲线。外层每评估一次某个电价方案，就要完整调用一次内层WOA，拿到下层最优响应后再计算外层的适应度。
写成伪代码是这样的：
% 双层鲸鱼算法伪代码 初始化上层种群 P_up (N_up 个个体，维度=电价时段数) 评价上层适应度: for i = 1:N_up price = 解码(P_up(i)) [L_best, f_best] = lower_woa(price) % 内层WOA求用户最优响应 Fit_up(i) = cal_upper_fitness(price, L_best) end for t = 1:T_up for i = 1:N_up % 按WOA三种机制更新上层个体位置 % 边界处理 % 重新调用 lower_woa 计算新的 Fit_up(i) end 更新全局最优 end注意这里有一个关键点：下层WOA不是只运行一次，而是上层每个个体每次评价都要运行一次。假设上层种群30、迭代100次，下层种群20、迭代200次，总的目标函数评估次数就是30×100×20×200 = 1200万次（还要算上上下层目标函数的复杂度）。所以下层模型能简化就尽量简化，能用解析表达式的就别让用户内部再嵌套一层复杂的模拟。
2.3 初始化和参数设置的讲究
参数设置直接影响求解速度。我的经验值是：上层种群N_up取20~50，上层迭代T_up取50~200；下层种群N_low取20~40，下层迭代T_low取100~500。如果你只是做小规模概念验证，可以全部往小了设，先把链路跑通，再逐步加大。对于24时段的电价优化问题，N_up=30、T_up=100，N_low=20、T_low=200这个组合在普通台式机上大概要跑十几分钟，还能接受。
初始化时所有位置要在边界内均匀随机生成。电价不能是负数，用户负荷功率不能超上限。对于离散决策变量（比如洗衣机开关），我一般不直接用0/1去套WOA，而是把变量编码成连续量，比如表示“在某时段开启的功率”，再用约束保证整数性或开关逻辑。如果实在要用离散变量，可以用四舍五入+罚函数，但要注意梯度不连续带来的搜索停滞问题。
还有一个容易被忽略的点：固定随机种子。Matlab里直接跑rand每次结果都不同，上层的最优解可能每次都变，这对调试非常痛苦。写代码第一行就加rng(2024)，保证结果可复现。后面要测试不同随机种子对结果稳健性的影响时，再把这个种子列表循环跑。
3. 居民负荷分层调度模型的详细建模
3.1 上层运营商的优化目标与约束
上层是负荷聚合商或者售电公司，它从批发市场买电，再以零售价格卖给居民用户。它的成本包括购电成本和可能存在的峰谷差惩罚，收入是卖电所得。所以上层目标可以写成最小化“购电成本 − 售电收入 + 峰谷差惩罚”，如果收益为正，目标值整体为负，这不影响优化。
一个典型的24小时模型：
min F = Σ_t (C_wholesale_t × L_t) − Σ_t (price_t × L_t) + α × Σ_t (L_t − L_avg)²
其中：
- L_t = P_base_t + Σ_j x_jt，即基础负荷加上可控负荷；
- price_t 是上层要决策的分时电价；
- C_wholesale_t 是已知的批发市场价格曲线；
- α 是峰谷差/波动惩罚系数；
- L_avg 是日平均负荷。
约束条件要包括：
- 电价上下限：price_min ≤ price_t ≤ price_max
- 总负荷上限：max(L_t) ≤ 变压器容量
- 价格平滑约束：|price_t − price_{t−1}| ≤ Δ，防止电价突变导致用户无所适从
为什么加上峰谷差惩罚？因为如果只让售电公司最大化卖电收入，它完全可以一直定高价，但用户响应后负荷会急剧萎缩，系统峰谷差也会很难看。把波动项放进目标函数，相当于让上层在“多卖电”和“平稳用电”之间找平衡。
3.2 下层用户的响应行为建模
下层用户不是毫无约束地追求电费最低，还要考虑舒适度。你让用户半夜爬起来用洗衣机，虽然电费便宜，但舒适度损失太大，用户不会干。所以下层目标函数通常是“电费 + 舒适度损失”最小化。为了简化，舒适度损失可以用“与原始用电习惯的偏移量平方和”来刻画。
一个单用户的下层模型：
min f = Σ_t price_t × x_t + ω × Σ_t (x_t − x_ref_t)²
约束：
- Σ_t x_t = E_total，总用电量必须完成；
- 0 ≤ x_t ≤ x_max，每时段功率上限；
- 如果是电动汽车，可能还要加充电时间窗约束。
这里 x_t 是该用户可控负荷在t时段的功率，x_ref_t 是用户原本计划中的功率曲线，E_total是该负荷一天的总耗电量。ω越大，说明用户越不愿意改变习惯；ω越小，用户越愿意为了省钱而转移负荷。实际算例中，ω要根据用户调查或历史数据来标定，没有数据时可以从0.01到10之间扫描几组值。
多用户的情况就是每个用户单独调用一次下层WOA，最后把各自的 x_t 累加得到总负荷曲线，返回给上层。每个用户之间互不通信，各自优化，这也符合非合作博弈的设定。
3.3 双层目标函数怎么衔接迭代
上下层的衔接靠两个向量：从上层传到下层的是电价向量 price，从下层传回上层的是最优总负荷曲线 L_best。在每次上层个体评价时，流程是：
- 上层WOA产生一个电价向量 price；
- 把这个 price 作为参数传给下层；
- 下层WOA在当前 price 下求解用户模型，得到最优负荷；
- 上层用它自己的目标函数，把 price 和 L_best 代进去，计算适应度；
- 上层WOA根据这个适应度更新种群位置。
这里要注意“评价阶段”和“更新阶段”不能混。同一代的上层个体应该都基于“同一批新位置完成后”才能更新，不能边更新边评价。原因很简单：如果个体A在迭代t代评价，而个体B在第t代还没更新完就开始评价，那么两者比较的不是同一状态下的结果，适应度没有可比性。
下层返回的内容不需要是用户所有内部变量的详细快照，只需要总负荷曲线。这也正好符合隐私保护的需求。
4. 完整代码实现与运行调试
4.1 主程序框架与文件清单
Matlab代码建议按功能拆成几个文件，别把500行堆在一个脚本里：
main_double_woa.m：设置参数，固定随机种子，调用主循环，输出结果。load_scenario.m：定义批发价格、基础负荷、用户数量、可控负荷参数。upper_woa.m：实现上层WOA循环。lower_woa.m：实现下层WOA，输入price，输出最优负荷曲线。cal_upper_fitness.m：计算上层目标函数。cal_lower_fitness.m：计算用户目标函数和罚函数。
我先写一个最精简的main框架：
% main_double_woa.m clear; close all; clc; rng(2024); scn = load_scenario(); params_up.N = 30; params_up.T = 100; params_up.dim = 24; % 分时电价维度 params_up.lb = 0.3 * ones(1,24); params_up.ub = 1.2 * ones(1,24); params_low.N = 20; params_low.T = 200; params_low.dim = 24; % 可控负荷功率曲线维度 params_low.lb = zeros(1,24); params_low.ub = 20 * ones(1,24); [best_price, best_fit, history] = upper_woa(params_up, params_low, scn); plot(history);这个小框架能让你一眼看清有哪些输入输出。后面所有细节都围绕这个框架填充。
4.2 上下层代码如何安全传递数据
建议用参数结构体传递，不要用全局变量。全局变量在Matlab里虽然写着方便，但一旦函数嵌套多了，很容易出现莫名其妙的数据覆盖，改bug时想死。
下层函数签名：
function [L_best, f_best] = lower_woa(price, user_params) % user_params 包含基础负荷、可控负荷能量、功率上限、omega等 % 返回 L_best: 24维总负荷曲线 end上层个体评价代码：
for i = 1:N_up price = positions(i, :); [L_best, ~] = lower_woa(price, user_params); upper_fitness(i) = cal_upper_fitness(price, L_best, grid_params); end这里grid_params是批发价格、变压器容量、惩罚系数等。用struct一次性打包好，调用非常清晰。
另一个细节：下层WOA内部迭代时，如果某次更新后个体越界，千万不要直接丢弃，要把它拉回边界附近。边界吸收策略：x(x<lb) = lb(x<lb); x(x>ub) = ub(x>ub)。对于需要用罚函数处理的约束，可以在适应度末尾加上一个极大惩罚项，但别让惩罚项大到把所有正常目标淹没，否则搜索会偏向“满足约束但不优”的区域。
4.3 五类高频报错与解决记录
我整理了一些实际调试中遇到的高频问题，做成速查表：
| 问题 | 可能原因 | 解决方法 |
|---|---|---|
| 下标越界 | 索引非整数或超维度 | 检查round()与size()，断点打印索引 |
| 适应度全是NaN | 下层返回空/罚函数溢出 | 加小对数保护，检查下层是否收敛且返回值非空 |
| 跑得太慢 | 下层评价次数过多 | 减小下层种群/迭代，或启用parfor并行 |
| 结果上下震荡 | 上层电价变化太剧烈，下层尚未收敛 | 增加价格变化平滑约束，提高下层迭代次数 |
| 不可复现 | 未固定随机种子 | 在main开头使用rng(seed)；并行时用RandStream |
对于NaN问题，尤其要检查下层返回的负荷向量是否包含空数组。可以在下层函数末尾加一句判断：
if isempty(L_best) L_best = zeros(1,24); % 给一个默认值，至少不会让上层崩掉 end不建议直接让上层目标函数变成NaN，因为WOA比较适应度时NaN会导致排序失效。宁可用大数惩罚。
5. 仿真算例与结果分析
5.1 算例场景与参数表
设计一个可复现的小算例：一个小区有100户居民，每一户有基础负荷和可控负荷，可控负荷总能量按小区汇总。上层决策24小时分时电价，批发价格曲线模拟典型的“双峰”形态：早高峰8-11点，晚高峰18-22点，其余时段较低。这样上层为了压低峰谷差，可能会在低谷时段降低零售价，鼓励用户把可转移负荷挪到半夜。
参数表如下：
| 参数 | 数值 | 说明 |
|---|---|---|
| 上层种群 N_up | 30 | 电价个体数量 |
| 上层迭代 T_up | 100 | 主循环代数 |
| 下层种群 N_low | 20 | 用户响应个体数量 |
| 下层迭代 T_low | 200 | 内层优化代数 |
| 电价下界 | 0.3 元/kWh | 防止恶意低价 |
| 电价上界 | 1.2 元/kWh | 防止恶意高价 |
| 用户数 | 100 | 简单等比聚合 |
| 可控负荷总能量 | 1500 kWh | 所有可转移设备汇总 |
| 每时段小区功率上限 | 500 kW | 变压器容量约束 |
| 舒适度权重 ω | 0.5 | 用户敏感度权重 |
这些参数是我调得比较稳的一组，你可以按自己需求改。注意用户数多但下层模型如果是线性凸的，可以用统一等效用户，这里因为用WOA解，直接聚合求解省时间。
5.2 收敛性能与优化效果对比
用上述参数跑出来的典型结果是：上层适应度在前30代下降很快，后面逐渐平缓。原因很简单，WOA前期以随机探索为主，很快找到好的电价区间，后期收敛到局部精细搜索，变化就小了。如果你发现适应度后面还在剧烈波动，大概率是下层响应没有稳定下来，或者电价平滑约束没加。
优化前后对比可以用一张示例表：
| 指标 | 优化前 | 优化后 | 变化幅度 |
|---|---|---|---|
| 峰时电价 | 0.9 元/kWh | 1.08 元/kWh | 上调20% |
| 谷时电价 | 0.5 元/kWh | 0.35 元/kWh | 下调30% |
| 用户总电费 | 100% | 约92% | 下降8%左右 |
| 峰谷差 | 1200 kW | 约1050 kW | 下降12%左右 |
| 上层利润 | 基准值 | 约基准的110% | 提升约10% |
不要盯着具体数字，关键是趋势：电价确实向“峰高谷低”演化，用户费用略降，峰谷差收窄，上层利润略升。这基本符合Stackelberg均衡的预期——领导者虽然提峰价，但用户会转移负荷，总体成本结构变优；跟随者用舒适度换电费下降，双方都获得好处。
5.3 负荷转移效果分析
观察优化后的典型负荷曲线，能明显看到原本18-22点的负荷尖峰被削掉了一部分，转移到了凌晨1-5点。可转移负荷里，洗衣机、洗碗机这类对时间不敏感的设备贡献最大；电动汽车因为充电时间长，如果时间窗允许，也会大量挪到谷时。但空调等温控负荷受用户舒适度约束较强，挪动不明显。
6. 调参经验与避坑指南
6.1 让双层WOA稳定收敛的调参顺序
很多人一上来就调上层参数，结果发现下层响应乱七八糟。我建议的顺序是：先固定一组合理的初始电价，只调下层WOA。让下层重复运行10次，如果10次得到的最优负荷曲线基本一致，说明下层够稳；如果每次都不一样，说明下层种群太小或迭代太少，请把下层的N_low和T_low调大，直到稳定。
然后再调上层。一开始上层种群和迭代用较小值，比如N_up=15、T_up=30，只要能跑通主循环就行。跑通后再增大参数。这样做的好处是避免上层在不断变化的噪声中找方向，浪费时间还找不到。
还有一个技巧：给电价变化加“惯性”或“平滑限制”。比如每代允许电价相对上一代整体移动最多20%。这相当于给上层搜索加了阻尼，防止因为用户响应过于敏感而出现价格来回震荡的“死循环”。
6.2 我在实操中踩过的三个坑
第一个坑：直接用标准WOA解带约束的下层问题，结果大量个体越界，尤其是洗衣机这种有开关次数约束的设备，越界后罚函数计算出的适应度全是异常值。后来我调整策略：变量采用连续功率表示，约束用“总能量=固定值”这种线性约束，通过罚函数和边界吸收配合，才稳定下来。
第二个坑：下层WOA每次都是从同一个随机初始位置开始，但每次迭代的随机数序列不同，导致同一个电价下两次调用下层，返回的最优负荷曲线差异很大。上层适应度在这个噪声下根本没法比较。后来我在上个层的每个个体开始内层优化之前，用rng(1000+i)固定一个子随机流，保证同价格返回相同的下层结果，上层更新才变得可比较。
第三个坑：为了提速用parfor并行评价上层个体，但没有控制每个worker的随机数流。结果并行跑完发现所有worker结果一样，等于没并行。后来我在parfor循环内部使用RandStream.create('mt19937ar','Seed', seed_per_worker)来独立分配随机流，或者干脆在并行前把所有需要调用的下层任务按“个体编号”分配好种子，这样才能真正并行出不同结果。
这些小问题看起来不起眼，但每一个都能让你白跑一晚上。希望我的记录能帮你少走几步弯路。双层鲸鱼算法本身不强求“最优收敛曲线漂亮”，它更像一个灵活的求解骨架——上面挂的是博弈模型，下面挂的是用户响应模型。你用Matlab跑通之后，再换别的算法、改目标函数、加新的负荷类型，都会顺手很多。
