冷热电多微网双层优化配置：基于储能电站服务的MATLAB仿真-平芜编程栈

从事综合能源系统仿真这一行的人，大概都逃不过“冷热电多微网”和“双层优化配置”这两座大山：一个是把电、热、冷三种能源形式捆在一起搞协同，另一个是两层优化模型相互嵌套、来回迭代。我本人因为项目需要，用MATLAB完整做过一套“基于储能电站服务的冷热电多微网系统双层优化配置仿真”，从数学模型搭建、YALMIP建模、CPLEX求解到结果分析，踩了不少坑，也积累了一些能直接复用的经验。今天这篇就把整个思路、代码结构和调试技巧一次说透，给正在写论文、做课题或者搞工程仿真的朋友一个可参考的路线。

这套仿真解决的核心问题很明确：多个冷热电联供型微网，共用一座储能电站，储能既做电量的“蓄水池”又能降低各微网的容量投资；上层优化决定各微网的燃气轮机、吸收式制冷机、电制冷机以及储能电站的容量配置，下层优化在每个典型日内做冷热电联合调度，两层之间通过运行成本来回反馈，最终找到总投资和总运行成本综合最优的方案。适合谁看呢？准备用MATLAB做微电网/综合能源系统优化的研究生，做园区多能互补方案设计的工程师，以及对双层规划代码实现有兴趣的读者都可以参考。

1. 项目背景与整体构架

1.1 为什么做冷热电多微网+共享储能

单个微网做冷热电联供其实已经很常见：燃气轮机发电，余热用来供暖或者驱动吸收式制冷机，实现能量梯级利用。但单个微网的系统容量小、负荷波动大，设备配置往往要么偏大浪费、要么偏小不够用。多个微网放在一起，负荷曲线可以错峰互补，再配上共享的储能电站，就能用更小的总容量满足所有微网的峰值需求。这就是我当初选择“多微网+共享储能”而不是“单微网独立配置”的根本原因。

从系统架构来说，每个微网内部都有燃气轮机（GT）、余热锅炉或换热装置、吸收式制冷机（AC）、电制冷机（EC），微网之间通过低压交流母线相连，并统一接入一个共享储能电站（ESS）。储能电站既可以从某个微网购电，也可以向另一个微网放电，起到跨微网的电量转移作用。这样设计的好处是节省投资，代价是调度模型更复杂，因为储能的服务对象是多个微网，充放电功率分配必须和微网的实时运行状态耦合。

1.2 双层优化的核心思路

双层优化（Bi-level Optimization）乍一听很抽象，用大白话说就是：上层定“买多大设备”，下层决定“设备怎么运行”。下层运行成本是上层容量决策的“打分器”，容量配置得好，下层调度就灵活、运行成本低；容量配置得差，下层再怎么优化也兜不住运行成本。两层之间形成一种领导者和跟随者的博弈关系，数学上可以用带上下层嵌套的优化问题描述。

在MATLAB仿真中，我采取的是“主从递阶迭代”思路：先给定一组初始容量方案，交给下层做典型日运行优化，得到运行成本和调度策略；然后把这个运行成本反馈给上层，上层在容量投资和运行成本的加权和最小化目标下更新容量方案；再带着新容量方案重新做下层优化，反复迭代直到容量方案和运行成本收敛。这个过程说起来简单，代码实现却有各种细节要注意，后面专门讲。

2. 数学模型：从目标函数到约束条件

2.1 上层容量配置模型

上层优化的决策变量是各微网内主要设备的安装容量，以及共享储能电站的额定功率和额定容量：

第i个微网：燃气轮机容量$P_{GT,i}^{cap}$、吸收式制冷机容量$C_{AC,i}^{cap}$、电制冷机容量$C_{EC,i}^{cap}$；
共享储能：额定功率$P_{ESS}^{cap}$、额定容量$E_{ESS}^{cap}$。

上层目标函数是年综合成本最小化，包括设备投资的等年值成本、年运行成本（由下层优化结果折算）和维护成本。计算公式大致是：

$$\min C_{total} = C_{inv} + C_{ope}$$

其中投资成本：

$$C_{inv} = \sum_{i=1}^{N} \left( \alpha_{GT} P_{GT,i}^{cap} + \alpha_{AC} C_{AC,i}^{cap} + \alpha_{EC} C_{EC,i}^{cap} \right) + \alpha_{ESS} \left( P_{ESS}^{cap} + \beta \cdot E_{ESS}^{cap} \right)$$

这里的$\alpha_{GT}$、$\alpha_{AC}$、$\alpha_{EC}$分别是燃气轮机、吸收式制冷机、电制冷机的单位容量投资成本，$\alpha_{ESS}$是储能单位功率成本，$\beta$是储能单位容量成本的折算系数，具体按常见工程数据取值。等年值时把设备寿命和折现率算进去，比如燃气轮机寿命20年、折现率8%，用等年值系数折算。上层约束主要是设备容量上下限、微网可用安装面积折算的最大容量、储能电站总占地约束等。

2.2 下层运行调度模型

下层的任务是在给定设备容量下，对每个典型日做逐时段（我取24小时，步长1小时）的冷热电联合调度。决策变量包括各微网燃气轮机的电出力$P_{GT,i}(t)$、吸收式制冷机从余热中制取的冷量$C_{AC,i}(t)$、电制冷机耗电制取的冷量$C_{EC,i}(t)$、与上级电网的购售电功率$P_{grid,i}(t)$、储能充放电功率$P_{ch,i}(t)/P_{dis,i}(t)$，以及各微网间的交互功率。

下层目标函数是典型日运行成本最小化：

$$\min C_{day} = \sum_{t=1}^{T} \left[ \sum_{i=1}^{N} \left( c_{buy}(t) P_{grid,i}(t) + c_{gas} G_{fuel,i}(t) \right) + c_{OM} \right]$$

这里$c_{buy}(t)$是分时电价，$G_{fuel,i}(t)$是燃气轮机耗气量，$c_{gas}$是天然气价格，$c_{OM}$是运维成本。约束条件要仔细列，少一条都会让仿真结果失真：

电力平衡约束：任意时段任意微网，电源出力加购电、储能放电，等于电负荷加电制冷机耗电、储能充电。

$$P_{GT,i}(t) + P_{grid,i}(t) + P_{dis,i}(t) = P_{load,i}(t) + P_{EC,i}(t) + P_{ch,i}(t)$$

热平衡约束：燃气轮机余热回收的热量，一部分供给热负荷，一部分驱动吸收式制冷机。

$$H_{rec,i}(t) = H_{load,i}(t) + H_{AC,i}(t)$$

这里余热回收量和燃气轮机发电量之间存在热电比关系。燃气轮机效率$\eta_e$、余热回收效率$\eta_{rec}$满足：

$$H_{rec,i}(t) = \frac{\eta_{rec} (1-\eta_e)}{\eta_e} P_{GT,i}(t)$$

冷平衡约束：吸收式制冷机和电制冷机的制冷量之和等于冷负荷。

$$C_{AC,i}(t) + C_{EC,i}(t) = C_{load,i}(t)$$

设备出力上下限、爬坡约束、储能SOC约束。储能约束是下层模型里最容易出问题的地方，尤其SOC递推公式：

$$SOC(t+1) = SOC(t) + \eta_{ch} P_{ch}(t) \Delta t - \frac{P_{dis}(t)}{\eta_{dis}} \Delta t$$

同时要保证SOC在$[SOC_{min}, SOC_{max}]$内，避免过充过放。

2.3 双层问题的耦合关系

上层和下层通过“容量方案→运行成本”这条线耦合在一起。上层一个容量决策，对应下层一个最优运行成本，下层的最优值函数本身是上层的约束或目标的一部分。严格意义上，这是一个带均衡约束的数学规划问题（MPEC），典型解法是用KKT条件把下层替换成上层约束，但KKT推导和编程对初学者都偏难。我的实现采用迭代反馈法，工程精度够用，且代码容易调试，后面会详细说。

3. 求解策略与MATLAB代码实现

3.1 求解框架选型

付完模型，先别急着写代码。求解框架主要有三条路可选：

KKT条件法：把下层问题用KKT条件替换，形成单层混合整数非线性问题，用求解器一把梭。优点是一次求解，理论最优，缺点是KKT推导复杂，大规模问题求解慢。
启发式迭代法：外层用粒子群或遗传算法搜容量方案，内层用确定性优化算运行成本。优点是实现直观，缺点是要反复调粒子群参数，容易陷入局部最优。
交替迭代法：固定下层、解上层，固定上层、解下层，循环几次直到收敛。速度快，但收敛性没有严格保证。

我的项目选择第三种：外层用整数/连续混合优化（MATLAB自带的fmincon或YALMIP），内层用YALMIP+CPLEX做线性规划。原因是冷热电联供系统的下层模型很大程度是线性或可线性化的，CPLEX求解快且稳；上层变量不多，fmincon足够。如果你追求严谨，可以在迭代完成后用KKT条件做一次校验，至少验证一下最优性gap。

3.2 核心代码结构逐段解析

整个MATLAB工程分五个文件：主程序main.m、参数配置文件config_data.m、上层优化函数upper_optimization.m、下层优化函数lower_optimization.m和结果绘图plot_results.m。这里给出关键代码段的思路和可复用片段。

主程序框架如下：

%% main.m 双层优化主程序 clear; clc; close all; config_data; % 载入所有基础数据，包括负荷曲线、电价、设备参数、初始容量 max_iter = 30; % 最大迭代次数 tol = 1e-3; % 收敛阈值：相邻两次运行成本相对偏差 Capex_cur = init_cap; % 初始容量方案 Cday_prev = 1e6; for iter = 1:max_iter % 下层运行优化：输入容量，输出典型日运行成本 Cday_cur = lower_optimization(Capex_cur, data); % 上层容量优化：输入运行成本，更新容量方案 Capex_new = upper_optimization(Cday_cur, Capex_cur, data); % 收敛判断 if abs(Cday_cur - Cday_prev) / Cday_prev < tol break; end Cday_prev = Cday_cur; Capex_cur = Capex_new; end

注意上层更新容量后，一定要回到下层重新算运行成本，而不是直接拿上一次的运行成本来算综合成本。因为容量变了，下层最优调度和运行成本都会变，若沿用旧运行成本会导致结果失真。

下层优化函数用YALMIP建模，核心如下：

function Cday = lower_optimization(cap, data) % cap 包含各微网GT容量、AC容量、EC容量、ESS功率和容量 n_mg = data.n_mg; T = 24; % 定义变量：各微网各时段的出力、购电、储能功率 P_GT = sdpvar(n_mg, T, 'full'); P_grid = sdpvar(n_mg, T, 'full'); C_AC = sdpvar(n_mg, T, 'full'); P_EC = sdpvar(n_mg, T, 'full'); P_ch = sdpvar(n_mg, T, 'full'); P_dis = sdpvar(n_mg, T, 'full'); SOC = sdpvar(n_mg, T+1, 'full'); % 约束集合 Constraints = []; Objective = 0; for i = 1:n_mg % 冷热电平衡、设备限值、储能SOC递推等约束 Constraints = [Constraints, P_GT(i,:) >= 0, P_GT(i,:) <= cap.P_GT_cap(i)]; Constraints = [Constraints, C_AC(i,:) >= 0, C_AC(i,:) <= cap.C_AC_cap(i)]; Constraints = [Constraints, P_EC(i,:) >= 0, P_EC(i,:) <= cap.P_EC_cap(i)]; Constraints = [Constraints, P_dis(i,:) - P_ch(i,:) <= cap.P_ESS_cap]; Constraints = [Constraints, P_dis(i,:) >= 0, P_ch(i,:) >= 0]; % SOC递推，充放电效率 Constraints = [Constraints, SOC(i,1) == 0.2 * cap.E_ESS_cap]; Constraints = [Constraints, SOC(i,2:T+1) == SOC(i,1:T) ... + data.eta_ch * P_ch(i,:) * data.dt ... - P_dis(i,:) / data.eta_dis * data.dt]; Constraints = [Constraints, SOC(i,:) >= 0.1 * cap.E_ESS_cap]; Constraints = [Constraints, SOC(i,:) <= 0.9 * cap.E_ESS_cap]; end % 目标：购电成本 + 燃气成本 + 运维成本 for t = 1:T for i = 1:n_mg G_fuel(i,t) = P_GT(i,t) / data.eta_e(i) / data.LHV; % 燃气耗量 Objective = Objective + data.price(t) * P_grid(i,t) ... + data.gas_price * G_fuel(i,t) ... + data.c_om * (P_GT(i,t) + P_EC(i,t)); end end ops = sdpsettings('solver', 'cplex', 'verbose', 0); optimize(Constraints, Objective, ops); Cday = value(Objective); end

这段代码要说明几个关键点。第一，储能容量E_ESS_cap在SOC初始化时参与了计算，SOC初值取20%的额定容量，避免初始时刻储能未投入或过满。第二，充放电功率用了P_dis - P_ch合并表达，实际表示净放电功率，这样写的好处是避免同时充放电的死锁问题，但约束里要单独限制两者非负且不超过额定功率。第三，燃气耗量计算依赖燃气轮机效率eta_e和天然气低位热值LHV，单位要注意统一。

G_fuel变量在代码里是临时计算的，但我在真实项目里是把它定义为sdpvar，因为如果有热电比耦合约束，燃气耗量和余热回收量都需要作为变量参与约束。上面代码为了清晰简化了，实操建议把G_fuel定义为显式变量，并在约束中关联。

上层优化函数相对简单：

function cap_new = upper_optimization(Cday, cap_old, data) % 上层目标：投资等年值 + 年运行成本 % 年运行成本 = Cday * data.day_num（典型日天数加权） x0 = [cap_old.P_GT_cap, cap_old.C_AC_cap, cap_old.P_EC_cap, ... cap_old.P_ESS_cap, cap_old.E_ESS_cap]; lb = [data.P_GT_min, data.C_AC_min, data.P_EC_min, data.P_ESS_min, data.E_ESS_min]; ub = [data.P_GT_max, data.C_AC_max, data.P_EC_max, data.P_ESS_max, data.E_ESS_max]; obj = @(x) inv_cost(x, data) + data.day_num * Cday; options = optimoptions('fmincon', 'Display', 'off', ... 'MaxIterations', 500, 'OptimalityTolerance', 1e-4); [x_opt, ~] = fmincon(obj, x0, [], [], [], [], lb, ub, [], options); cap_new.P_GT_cap = x_opt(1:data.n_mg); % ... 同理赋值其他容量 end

这里最大的坑是：fmincon是连续变量优化，如果你给的是连续容量，那么下层模型的整数变量（比如GT启停状态）就完全没考虑。严格讲，微网内GT应该考虑启停机状态，0/1变量会让下层变成混合整数规划，上层容量也是离散候选集。我的处理方式是把容量连续化，下层也不考虑启停，把GT最小出力比例限制在30%，用连续出力近似，工程上说得通，论文里加一句“本文忽略机组启停成本，按连续可调处理”即可。

3.3 数据输入与结果输出模块

config_data.m处理所有原始数据，包括：

典型日选取：一般取春秋、夏、冬三个典型日，或者按季节性负荷聚类。我按夏季制冷、冬季采暖、过渡季通风三类场景设置权重，日数按30天、90天、240天近似，或者用实际气象天数加权。这个权重会直接影响年运行成本，是个关键输入参数。
分时电价曲线：峰段1.2元/kWh、平段0.8元/kWh、谷段0.4元/kWh，具体峰谷时段按当地政策设置。
负荷曲线：每个微网都有一套独立的电、热、冷日负荷曲线，我用的是实测数据生成，如无实测可以用典型负荷比例加随机波动生成。
设备参数：燃气轮机效率、热电比、吸收式制冷机性能系数COP、电制冷机COP、储能充放电效率等。这些数据来源要可靠，我一般查厂家手册和文献，不要随意取。

结果输出模块核心是绘制三类图：各微网典型日电负荷平衡图（面积图）、储能SOC变化曲线、上层迭代收敛曲线。我以前还加画一个“不同容量方案成本散点图”，可以直观看出容量增大导致投资上升、运行成本下降的trade-off关系。

4. 仿真算例与结果分析

4.1 算例场景与参数设定

算例设定三个微网，每个微网最大电负荷分别为200kW、150kW、180kW，最大冷负荷分别为180kW、120kW、160kW，最大热负荷分别为150kW、100kW、120kW。各设备候选容量范围：GT 0~300kW，AC 0~250kW，EC 0~200kW。共享储能额定功率候选范围10~100kW，额定容量50~300kWh。单位投资成本按经验赋值：GT约8000元/kW，AC约1200元/kW，EC约800元/kW，储能按2000元/kWh计算容量成本、500元/kW计算功率成本。折现率8%，项目寿命20年。

不同时段电价和负荷波动设计成夏季中午光伏反调峰效应明显，晚间冷负荷高峰，储能正好在午间低谷充电、晚间高峰放电。

4.2 典型日运行结果解读

跑完仿真后，典型的夏季日调度结果呈现出“燃气轮机全天基荷+电制冷机跟随峰谷”的形态：白天电价平段和峰段，GT满发，余热驱动AC承担大部分冷负荷；晚间电价峰段，GT仍满发，但电负荷下降，富裕电量给储能充电；深夜谷段，GT降低出力，不足电量从电网购电，储能放电补足。电制冷机EC白天几乎不出力，只在AC容量不够或余热不足时才启动。这说明系统确实实现了“热跟随电、冷跟随热”的梯级利用逻辑。

储能SOC曲线是判断逻辑正确性的重要指标。正常结果里，SOC应该呈现“夜间谷段充电→白天峰段放电→午间光伏/低价时段再充电→傍晚高峰再放电”的波动。如果SOC曲线整天贴着上限或下限走，说明储能容量设置不合理或电价引导失效。

4.3 与单微网独立配置方案的成本对比

为了验证共享储能+多微网协同的价值，我单独跑了一个对照算例：三个微网各自独立配置储能，储能容量之和与共享方案总容量相同，彼此之间无功率交互。结果通常共享储能方案的投资成本能降低15%左右，年运行成本降低5%~8%，综合年成本下降10%以上。原因很好解释：多微网的负荷峰值时间错开，共享储能总容量不用按三个微网峰值叠加，且储能充放电策略可以在全网范围内优化，单微网独立运作时只能各自为政，峰谷利用不充分。

这个对比实验强烈建议保留，因为论文评审或项目汇报里这是最直观的“价值证明”。你需要计算清楚各方案的投资差和运行差，甚至画一个堆叠柱状图展示成本构成。

5. 常见问题与调试经验

5.1 求解器选择与配置坑

YALMIP后端我用的是CPLEX，但CPLEX需要单独安装并配置许可证，这是新手最容易卡住的点。有几个替代方案：

MATLAB R2024a以后自带intlinprog、linprog，可以直接用solvesdp或者手动转换成MATLAB优化工具箱格式，不用额外装CPLEX；
Gurobi也可以，许可证申请流程和CPLEX类似，但学术许可更方便；
如果纯线性模型，linprog完全够用；有整数变量，用intlinprog。

如果使用CPLEX，一个常见报错是“CPLEX not found”，这时候在sdpsettings里指定'solver','cplex'之前，先运行yalmiptest检查YALMIP是否识别到求解器。还有个精细问题是，CPLEX对MILP问题默认启用并行，如果模型小反而变慢，可以在solveroptions里设置'cplex.mip.strategy.fp'调整搜索策略。

5.2 双层迭代不收敛的排查

我实际调试时最常碰到的就是迭代震荡——容量方案在两三个值之间来回跳，运行成本也上下波动不收敛。原因主要有三个：

上层目标函数对运行成本的反馈过于敏感，即运行成本在目标函数中占比过高。解决方法是给运行成本加一个平滑因子或者限制每次迭代容量变化的步长，比如cap_new = cap_old + alpha * (cap_solution - cap_old)，alpha取0.5~0.8，相当于阻尼更新。
下层模型有多重最优解，同样的运行成本对应不同调度策略，而上层计算投资时又依赖容量值，造成波动。这时候需要给下层目标函数加一个很小的正则项，比如加上储能充放电次数惩罚，打破对称性。
初始容量方案离最优解太远。建议先用启发式估算初始值：按各微网最大电负荷的40%作为GT初值、最大冷负荷的50%作为AC初值、最大冷负荷的20%作为EC初值，储能初值可以设成微网平均日用电量的20%~30%。