搞微网双层优化的人看到“双层”两个字心里要先有个预期这不是一个普通单层MPC能解决的事。我这个项目踩过不少坑今天把完整思路、建模方式、基于Yalmip的代码骨架、以及调试过程中的问题全梳理出来希望能帮你少走弯路。微网双层优化的本质就是主从博弈上层决策者——微网运营商制定电价策略或调度指令下层决策者——用户、负荷聚合商、分布式储能业主在接收到上层的信号后做自己利益最大化的响应。上下层各自有自己的目标函数和约束最终收敛到某种均衡。这种结构在需求响应、虚拟电厂、含多主体的微网交易中非常常见。实现这件事Matlab加Yalmip组合是很多人入门时的选择。Yalmip的最大价值是让你不需要关心不同求解器的低层语法用统一的建模语言描述约束然后一键切换Gurobi、CPLEX、Mosek或者国内常用的求解器。但我也要提醒一句Yalmip只负责建模双层模型里最难的互补条件线性化、大M法的边界问题都需要自己处理Yalmip本身不会替你完成。1. 微网双层优化到底要解决什么问题1.1 单层模型的局限与双层模型的必然很多刚接触微网调度的朋友会有个疑问经典的经济调度也是求成本最小为什么要整双层我先用一个例子说明单层的尴尬。假设一个微网内既有储能又有若干可调节的商业用户。如果只建一个混合整数线性规划把储能投资、用户电费、上级电网购电费全部叠加成一个目标函数得到的确实是一个数学意义上的全局最优解。但关键是这个“最优解”里用户可能被强制削减了用电或者储能频繁充放来实现总成本最低而这些牺牲没有“人”出来反对。实际问题中用户不是任人摆布的资产他们会根据电价做需求响应储能业主也不是免费工具他们要回收投资和保证收益。这就需要一个能反映不同主体之间决策顺序和博弈关系的模型。双层模型的价值不是让结果变得好看而是让结果更加符合现实中“你有你的策略、我有我的对策”的运行逻辑。1.2 微网中典型的上下层利益格局我这次大致按下面这个结构来搭建实例上层微网运营商。它管理着配电网关口处的购售电、储能系统的充放电策略同时向下层用户发布一个内部零售电价。下层用户/负荷聚合商。用户收到电价后在一定舒适度约束和总用电量约束下调整用电曲线。目标是最大化自己的效用或最小化自己的费用。上下层之间的连接点有两个第一个是电价信号上层给定电价下层以此做出响应第二个是用电功率曲线下层响应出的负荷曲线反馈给上层参与微网内的功率平衡约束。这里要注意不是随便传几个参数过去就叫双层关键是上下层各自的优化目标独立、决策顺序明显、交互通道明确。如果在模型里把用户成本也并入运营商的成本函数就又退化回单层了。1.3 为什么选择Matlab加Yalmip这套组合在确定用Matlab和Yalmip之前我也比较过其他方案用Python的Pyomo或CVXPY、用纯GAMS、或者直接在商业求解器里写模型。最终我还是选择Matlab加Yalmip理由有三点第一Matlab在数据前处理和结果可视化上效率高读入风光出力数据、负荷曲线、电价序列非常顺手画结果图时调整样式也很灵活。第二Yalmip写约束的直观度很高尤其对于带双层结构、大量上下标变量的模型可以直接用们与编辑距离。第三Yalmip对求解器的适配做得足够好从开源的求解器到商业求解器几乎无缝切换这在对比不同求解器性能时非常有用。2. 数学模型上下层怎么定义、怎么耦合2.1 上层优化模型的目标与约束我以上的项目实验配置为例一个含光伏PV、储能ESS、常规基础负荷和可转移负荷的微网向外通过公共连接点与配电网交换功率。调度周期 24 小时步长 1 小时。上层决策变量关口处购电功率与售电功率单位 kW储能充电功率和放电功率单位 kW储能荷电状态SOC单位面向用户的动态零售电价单位元/kWh。上层目标函数考虑四部分成本min ∑ [ C_buy_t * P_grid_buy_t - C_sell_t * P_grid_sell_t ] ∑ [ C_ess * (P_ch_t P_dis_t) ] - ∑ [ Price_t * P_user_t ]其中第一部分是从配网购电的成本第二部分是储能充放电的损耗与运行成本第三部分是微网向用户售电的收入。在表达式里售电收入前面是负号因为对运营商来说是收益。上层约束部分比较常规列出几个容易出问题的约束微网内功率平衡约束P_PV_t P_grid_buy_t P_dis_t P_load_base_t P_user_t P_ch_t P_grid_sell_t这里注意购电和售电不能同时为正否则模型会出现“从电网买高价电再低价卖回”的套利漏洞。有几种处理方式最简单的是用二元变量互斥或者用下面这种更稳的线性化0 ≤ P_grid_buy_t ≤ M * u_buy_t 0 ≤ P_grid_sell_t ≤ M * u_sell_t u_buy_t u_sell_t ≤ 1还有一种常见做法是不引入互斥变量而在目标函数里让购电价一直大于售电价这样模型本身就倾向于不同时买卖。但我实测发现如果电价数据有异常或边界情况还是要用互斥变量更保险。储能动态与SOC约束SOC_t SOC_{t-1} eta_ch * P_ch_t - P_dis_t / eta_dis SOC_min ≤ SOC_t ≤ SOC_max 0 ≤ P_ch_t ≤ P_chmax 0 ≤ P_dis_t ≤ P_dismaxSOC初值设为 0.5 左右末状态是否约束回 0.5 视需求而定。如果想要结果具备周期性就加上等值约束。电价调整范围约束Price_t Price_base_t d_Price_t Price_min ≤ Price_t ≤ Price_max固定基础电价加可调偏移这是很多双层电价优化论文里常用的简化。你也可以让电价直接连续可调但要注意下层用户反应对过大的电价波动会有抵触现实场景中电价调整一般有上限。2.2 下层优化模型的目标与约束下层是用户的用电成本最小化问题。用户面对上层下发的电价序列决定自己的用电功率曲线。目标函数如下min ∑ [ Price_t * P_user_t Miu * |P_user_t - P_ref_t| ]第一项是电费第二项是舒适度惩罚项。P_ref_t 是用户最舒服的自然用电功率一旦被要求转移到其他时段用户会产生不满。Miu 越大用户越不愿意响应价格信号。这个绝对值项看起来简单但处理不好会破坏模型的线性性。用辅助变量 z_t 做等价线性化z_t ≥ P_user_t - P_ref_t z_t ≥ P_ref_t - P_user_t z_t ≥ 0目标函数换成min ∑ [ Price_t * P_user_t Miu * z_t ]下层约束包括可调负荷上下限P_user,min ≤ P_user_t ≤ P_user,max一天总用电量约束∑ P_user_t E_total用户最大转移速率约束|P_user_t - P_user_{t-1}| ≤ ΔP_user_max这个约束描述用户用电习惯不允许剧烈变化。2.3 上下层耦合机制与求解思路上下层之间就是通过 Price_t 和 P_user_t 互相牵制。上层定电价下层返回用电量下层用电量又进入上层功率平衡约束和收益项。这种循环关系就是最典型的“主从决策”。求解这类双层模型有两条主流路径。第一条是KKT 单层化方法。把下层的线性规划问题用 Karush-Kuhn-Tucker 条件替换下层优化就变成了上层模型中的一组约束和一组互补条件。这样双层问题被转换为一个带互补约束的数学规划问题再通过大M法引入二进制变量将互补约束线性化最终变成一个混合整数线性规划。这个方法的优点是能求出全局最优解缺点是当下层问题不是线性规划或者含整数变量时KKT 条件不再适用。第二条是迭代启发式方法。上层采用粒子群算法、遗传算法等启发式算法生成电价策略然后在每代中调用 Yalmip 求解下层问题下层最优解返回后计算上层的适应度。这种思路可以处理更复杂的非线性下层模型但只能获得次优解而且计算量比较大。我在这次实现中优先采用了 KKT 单层化方法原因是规模不大、下层是 LP具备转化为 MILP 的条件求解稳定性较好而且用 Yalmip 写起来还算顺手。下面具体展开。3. 基于Yalmip的完整实现过程3.1 环境准备与求解器配置我用的环境是 Matlab 较新版本Yalmip 建议用最新版本求解器用的 Gurobi。需要注意几点第一Yalmip 不是一个求解器它只负责把模型转换成求解器能理解的格式。所以你的电脑上必须先装好求解器并且要让 Mexican 文件在路径中可被 Yalmip 调用。第二求解器的 License 也要检查。Gurobi 有免费学术版也有可以申请试用的评估版。装好之后在 Matlab 命令行运行optimize前建议先跑一遍check函数来验证求解器是否工作。简单验证代码optimize(sdpvar(1)0, -sdpvar(1), sdpsettings(solver,gurobi));如果求解成功并返回 infeasible 或 unbounded就说明链路通了。3.2 用Yalmip定义变量与约束顶层变量定义我习惯这样组织分块定义T 24; P_grid_buy sdpvar(1, T, full); P_grid_sell sdpvar(1, T, full); P_ch sdpvar(1, T, full); P_dis sdpvar(1, T, full); SOC sdpvar(1, T, full); d_price sdpvar(1, T, full); Price Price_base d_price; u_buy binvar(1, T); u_sell binvar(1, T); u_ch binvar(1, T); u_dis binvar(1, T);注意如果你是做连续变量默认 sdpvar 是 full 矩阵不必加第三个参数只有需要结构约束时才用模式参数。但明确写 full 语义更清晰。下层变量P_user sdpvar(1, T, full); z sdpvar(1, T, full); % 舒适度偏差辅助变量 lambda_soc sdpvar(1, T1, full); % 动态约束乘子3.3 下层KKT条件如何写这是这个项目里最关键的部分。先写出下层拉格朗日函数。下层目标分三块电费项、舒适度项、总用电量均衡约束的拉格朗日乘子部分、上下限约束的部分。对于每个约束要对应引入乘子。以标准 LP 形式的求最小问题为例若约束形式为 Ax ≤ b则乘子非负。互补松弛条件为lambda ≥ 0 Ax - b ≤ 0 lambda_i * (Ax - b)_i 0展开后你会得到每个时段 P_user_t 的驻点条件、每个约束的互补条件。这套互补条件里有“乘积等于0”的非线性项Yalmip 不能直接处理必须通过大M法线性化。互补条件的形式0 ≤ a ≤ M * bin 0 ≤ b ≤ M * (1 - bin)这里 a 代表乘子或松弛量b 代表约束左侧松弛量bin 是辅助二进制变量。核心在于大M值的选择这个后面我会重点讲。3.4 主程序代码一个可扩展的范本下面给出我实际整理过的核心代码框架。为便于理解我去掉了海量数据读入只保留模型构造逻辑。%% 参数设置 T 24; eta_ch 0.95; eta_dis 0.95; SOC_init 0.5; SOC_min 0.2; SOC_max 0.9; P_chmax 60; % 储能最大充电功率 kW P_dismax 60; % 储能最大放电功率 kW P_grid_limit 100; E_sto_cap 200; PV [0 0 0 0 0 2 6 12 18 25 30 33 31 27 22 15 8 4 2 0 0 0 0 0]; Load_base [35 34 32 31 30 38 45 52 58 60 55 53 50 46 42 45 48 52 50 44 42 40 38 36]; Price_base [0.6 0.55 0.5 0.5 0.5 0.65 0.8 0.9 1.0 1.0 0.95 0.9 0.85 0.8 0.85 0.9 0.95 1.0 1.05 1.0 0.9 0.7 0.6 0.55]; c_buy [0.65 0.6 0.55 0.55 0.55 0.7 0.9 1.05 1.1 1.05 1.0 0.95 0.9 0.85 0.9 0.95 1.0 1.1 1.15 1.1 0.95 0.75 0.65 0.6]; c_sell 0.45 * ones(T,1); M_cost 0.02; % 储能运行损耗系数 元/kWh Miu 15; % 舒适度惩罚权重 %% 决策变量 P_buy sdpvar(T,1); P_sell sdpvar(T,1); P_ch sdpvar(T,1); P_dis sdpvar(T,1); SOC sdpvar(T1,1); d_price sdpvar(T,1); Price Price_base d_price; u_buy binvar(T,1); u_sell binvar(T,1); P_user sdpvar(T,1); z_abs sdpvar(T,1); %% 约束集合 Constraints []; % 上层功率平衡 for t 1:T Constraints [Constraints, PV(t) P_buy(t) P_dis(t) Load_base(t) P_user(t) P_ch(t) P_sell(t)]; end % 关口购售电互斥 Constraints [Constraints, P_buy P_grid_limit * u_buy, P_sell P_grid_limit * u_sell, u_buy u_sell 1, P_buy 0, P_sell 0]; % 储能约束 Constraints [Constraints, P_ch 0, P_ch P_chmax * (P_ch 0), P_ch P_chmax]; Constraints [Constraints, P_dis 0, P_dis P_dismax]; Constraints [Constraints, SOC(1) SOC_init]; for t 1:T Constraints [Constraints, SOC(t1) SOC(t) eta_ch * P_ch(t) - P_dis(t)/eta_dis]; end Constraints [Constraints, SOC_min SOC(2:end) SOC_max]; % 电价区间 Constraints [Constraints, Price 0.4, Price 1.5]; %% 下层约束与KKT条件 % 下层可行域 P_user_ref Load_base; % 近似认为自然用电量 Constraints [Constraints, P_user 0.6*Load_base, P_user 1.2*Load_base]; Constraints [Constraints, sum(P_user) sum(Load_base)]; Constraints [Constraints, z_abs P_user - P_user_ref, z_abs P_user_ref - P_user, z_abs 0]; % 拉格朗日乘子使用辅助变量替代直接写KKT的乘积项 % 简化处理直接利用驻点条件与互补条件线性化 % 为了阅读清晰这里用一个统一的线性互补实现函数需自行实现或展开 % 实际建模时把每个下层约束对应的乘子写出后套用下节大M法模板 %% 下层驻点条件示例手动展开形式 % 对每个时段 t: 2*Miu*z_abs(t) Price(t) - mu_equ - mu_min mu_max 0 % 以及 z_abs 和互补条件的线性化 %% 目标函数 Objective sum(c_buy .* P_buy - c_sell .* P_sell M_cost * (P_ch P_dis)) - sum(Price .* P_user); %% 求解 ops sdpsettings(solver, gurobi, verbose, 2, debug, 0); sol optimize(Constraints, Objective, ops); if sol.problem 0 fprintf(求解成功\n); else sol.info end这只是框架代码。实际跑通还需要把每个互补条件的大M形式完整写出来。下面单列一节专门讲KKT和大M因为这里最容易被卡死。3.5 KKT单层化的完整展开与数值实现我把下层模型简写为标准形式下层变量 x[P_user;z_abs]下层约束包括不等式约束集 G(x)≤0、等式约束 H(x)0目标 f(x)。则有拉格朗日函数L f(x) λ^T * G(x) μ^T * H(x)KKT条件梯度条件∇f(x) J_G^T λ J_H^T μ 0原始可行性G(x) ≤ 0H(x) 0对偶可行性λ ≥ 0互补松弛λ_i * G_i(x) 0核心难点在互补松弛条件。转成线性需要G_i(x) ≤ 0 λ_i ≥ 0 λ_i ≤ M * bin_i -G_i(x) ≤ M * (1 - bin_i)这里 G_i(x) 是某个下层约束左侧的表达量。注意G_i(x) 本来就是 ≤0 的取 -G_i(x) 为松弛量所以 ≥0。由此每个互补条件只需要两个约束和一个二进制变量。大M取值要足够大但不能过大。过大可能导致数值病态、求解变慢甚至精度丢失。我建议根据每个约束的物理量纲先估算上界。例如 P_user 不超过 1.2 倍基础负荷那么对 P_user ≤ P_usermax 这个约束松弛量最大不超过 P_usermax - 0取 MP_usermax 就够如果取 1000 倍则没必要。当然实际找 M 常量时我一般先设一个保守的上界再通过求解后的对偶变量检查有没有乘子异常接近 M。还有一个容易被忽略的点如果下层存在两个相邻时段之间的耦合约束对应的乘子要成对写KKT中全局梯度条件要对整段时间的拉格朗日函数逐变量求偏导别漏项。我当时就漏过 SOC 的乘子项导致结果长期不收敛。3.6 求解器调用与结果输出用 Yalmip 的 optimize 函数选择 solver 为 gurobi。实际运行中我建议把 verbose 开到 2这样能看到 MIP gap 的下降过程如果长时间不动就可以判断是互补条件线性化有问题。求解完成后用 value() 提取各变量绘制电价曲线、功率平衡曲线、SOC曲线。结果呈现上我一般会把上层购售电、下层用户响应前后负荷、储能充放电放在同一张图里对比既能检查功率平衡也能直观看到需求响应的削峰填谷效果。4. 实现过程里的坑与排查技巧4.1 模型无解时先别急着调求解器这是我最想强调的一点。很多人一遇到 infeasible 就去看求解器日志或者改 MIP gap这是方向错了。绝大多数无解问题都出在约束本身自相矛盾或者边界条件不对尤其在双层模型单层化之后KKT条件本身就可能引入虚假不可行。建议的排查顺序第一先把下层KKT条件整体注释掉只保留上层自己的一套约束和变量用上层目标跑一遍确认上层自身可行。第二再把下层自身模型单独扔进 Yalmip 里固定一组 Price 值验证下层模型在单独求解时是否可行。这一步能把上下层问题拆开定位到具体那一侧出了问题。第三最后把互补条件加回去逐步放行。可以先只放一部分乘子试跑。我实践中最常见的问题是SOC 初始约束写到了 SOC(1)但在变量定义时又错误地把 SOC 当成了 T 维导致下标越界后的终端 SOC 约束根本没生效看起来模型一堆约束“成立”实际上却把路径堵死了。所以变量长度和下标一一对应检查时真的要一个萝卜一个坑。4.2 大M值取太大求解慢取太小会锁错解大M法的本质就是把一个非线性乘积条件变成线性条件代价是引入额外的二进制变量。而 M 的取值范围直接决定了解的可行域边界松弛程度。原则一M 必须大于松弛量的所有可能最大值。比如 P_dis 最大是 60那么某个互补条件中涉及 P_dis 的松弛量M 至少取 60 或略大。原则二不要用同一个 M 套所有约束。不同物理量纲差很多用一个巨大的通用 M 会导致矩阵条件数恶化Gurobi 也会给出数值警告。原则三如果求完解发现某个互补条件对应乘子奇异这里“奇异”不是报错而是求解结果上有不合理的尖峰就返回去降低 M。一个常用做法是先用模型里各个变量的自然上限公式化计算 M例如M_P_user 1.2 * max(Load_base) M_P_grid P_grid_limit M_soc SOC_max - SOC_min然后分别在每个互补约束里指定不同的 M。Yalmip 里这就是一个参数或常量不存在额外难度。4.3 求解结果里出现同时购售电怎么办如果你用了互斥二进制变量但没有把目标函数中购售电单价处理好仍可能出现购电和售电同时大于零。原因就是购电和售电都带上了很小的合理值且功率平衡约束刚好能通过比如 5 kW 购电和 0.1 kW 售电同时存在物理上不是明显违法但这个结果本质上不合理。虽然加了 u_buy u_sell ≤ 1理论上能保证互斥但如果某个时段购电量和售电量都大于零是因为松弛变量或互补条件放宽问题还是在对应时段的电价成本差。简单一句话关注每个时段购电价是否小于售电价目标里购电是成本、售电是收益如果购电价大于售电价模型自动会避免同时进出。如有异常多半是大M参数把购售电上限放开了。4.4 模型收敛慢怎么定位瓶颈MILP 规模一大分支定界会很吃力。Yalmip 在输出时有一个字段叫sol.solveroutput.infoGurobi 的日志会打印 BestBd 和 Incumbent如果 gap 一直卡在比如 1e-2 不动可能就是约束里存在数值对称性或重复约束。我一般会做三件事第一打印模型信息model export(Constraints, Objective, sdpsettings(solver,gurobi));看看二进制变量总数、约束总数、非零元素数。如果二进制变量超过预期检查辅助变量定义是不是冗余。第二给二元变量加合理的初始热启动。Yalmip 里可以直接赋值assign(u_buy, zeros(T,1)); assign(u_sell, ones(T,1)); % 根据时段电价特点给经验初值这通常会大大加快启动阶段。第三把 M 尽量收紧减少无效整数节点的探索空间。在一次测试里我只是把某个 M 从 1000 降到 120求解时间从 480 秒降到了 60 秒左右效果极其明显。4.5 别忘了数值单位的一致性很多初学的人容易犯的低级错误光伏功率用的是 kW储能容量用的是 MWh电价用的是元/kWh最后在功率平衡约束里数值差了一千倍还不自知。我的建议是建模前把所有数据统一到同一单位体系。这个项目里统一为 kW 和 kWh电价统一为元/kWh。数据加载后先做一道“极值审查”看看 PV 最大值、负荷最大值、购电功率上限之间有没有离谱差异。这个习惯能避免后面花几小时排查无解问题。5. 从教学Demo到工程化扩展的路径5.1 从确定性模型到考虑不确定性基础双层模型假设光伏出力和负荷是已知的确定值。实际微网中光伏预测误差和负荷波动都很大。扩展方向主要有两个一是鲁棒优化将 PV 出力的不确定集合引入上层约束目标函数改为在最坏情况下依然可行且成本尽量低。这会让模型变成典型的 min-max 结构跟双层博弈叠加后层级更复杂。二是机会约束规划把某些约束以概率形式给出比如“功率不平衡的概率小于 5%”。Yalmip 对机会约束的转化支持不如专门工具那么自动需要手动线性化确定等价式。5.2 从单微网到多微网互联多微网场景下每个微网既可以看作是上层运营者的受控对象也可以看作是其他微网的竞争对手。此时双层会变成多层或均衡约束规划EPEC用KKT单层化法依然可以建模但模型中会同时出现多套互补条件二进制变量数量大幅增加规模会膨胀。这种情况下我强烈建议先用单微网验证所有逻辑再逐步增加微网数量。一次加太多节点调试复杂度会指数上升。5.3 把模型嵌入滚动优化框架工程上更常见的用法是“双层滚动调度”——比如以 15 分钟为一个时间步每次重新求解未来 4 小时的优化问题。模型本身还是那些约束但需要在 Matlab 循环里反复调用 optimize。这里有个经验滚动框架下初始 SOC 不再是固定值而是上一轮解出的 SOC 末端值这个值必须从求解结果中提取并赋给下一轮变量初值。如果每次滚动都把 SOC 初始化成固定 0.5整个调度曲线会严重失真相当于每次都重新从满荷开始规划。5.4 求解速度优化的一些实测经验我在这个项目里实测过几种提速手段排序如下收紧大M收益最大几乎零成本给二进制变量热启动效果好但初始赋值要经验减少不必要的时间耦合变量如果多时段储能 SOC 不涉及阶梯电价可适当解耦但这会牺牲准确性换求解器参数把 Gurobi 的 MIPFocus 调到 1适合可行解容易找但证明最优困难的问题调到 3 则适合需要严格下界的问题最后才考虑用启发式替代精确求解当规模确实大到无法求解时用上层启发式套下层Yalmip这是次优解但至少能给你一个可用的方案。6. 我对这类项目的一点体会玩过几个这样的项目后我明显感觉双层优化的门槛不在目标函数如何复杂而在KKT单层化的工程细节和大M参数的数值敏感性。网上很多代码只给了框架一到关键地方就用省略号带过导致新手复现时处处碰壁。如果谁要上手我真心建议从最小规模做起先把时间周期压到4小时先不考虑储能SOC的严格等约束先把单一互补条件和大M调通过再逐步放大。这个小样本调试策略比一开始直接跑24小时要省太多时间。最后再分享一个小技巧调试时可以把目标函数的各项成本分别输出逐项对比一旦某个成本分量数值出现明显尖峰通常就能顺藤摸瓜找到互补条件或电价约束的问题。模型是死的但数值特征会向你透露出很多解题线索。
