拓冰网站定制 
1. 为什么储能电站服务下的多微网系统天然需要双层优化配置先交代一下背景。我最近一直在做储能电站相关的项目客户那边给的课题是基于储能电站服务的冷热电多微网系统双层优化配置要求用 MATLAB 实现给出代码框架和计算结果。说实话刚拿到这个题目的时候我也愣了一下——双层优化配置这个名词听起来很高大上但拆开来看就是两件事第一给储能电站和各微网里的设备定容量第二在容量定好之后安排各个微网在一天的运行计划。前者是建多大后者是怎么用。这两件事一旦放在同一个模型里单层优化是做不到的必须分层来解决。很多人第一次接触这类问题上来就想把容量配置和运行调度写进同一个目标函数里一次性求解。这种做法不是不行而是会导致一个很尴尬的结果设备容量和运行策略在模型里会互相打架。比如说储能容量定得大运行阶段就可以多充多放、降低购电成本但储能本身的投资成本又很高如果容量定得小投资省了可运行阶段又买高价电。这两个目标属于不同的时间尺度容量配置是年尺度决策运行调度是小时尺度决策硬揉在一起模型要么不收敛要么算出来的容量毫无工程意义。所以我在这篇文章里想把我从建模到 MATLAB 代码落地的全过程梳理一遍包括双层模型怎么拆、上层和下层分别优化什么、两层之间怎么迭代耦合以及代码里最容易踩的坑。如果你正在做微电网、综合能源系统、共享储能或者CCHP 多微网方向的课题这篇内容基本可以当成一份可抄作业的路线图来用。先明确一个概念什么叫储能电站服务。在传统的微网系统里储能是每个微网自建自用的容量小、利用率低而且投资主体和运行主体是同一个。但随着共享储能模式出现储能电站变成了一个独立的服务商它向多个微网提供充电、放电服务容量决策在一个主体手里运行调度却在各个微网手里。这种投资运营分离的结构天然就是双层优化的用武之地——上层做储能电站和各微网设备的容量配置下层做多微网系统的日运行调度。共享储能的意义在于降低储能成本、提高设备利用率而双层优化配置就是把这个意义真正落到数学建模和代码里的关键手段。2. 系统拓扑与设备建模先把冷热电多微网画清楚2.1 网络拓扑怎么搭在动手写代码之前我习惯先把系统拓扑在草稿纸上画出来。这个项目对应的结构大致是这样的系统包含若干个微网常见的做法是取 3 个典型微网每个微网内部有光伏、微型燃气轮机、燃气锅炉、电制冷机、吸收式制冷机这几类设备承担自己的电负荷、热负荷和冷负荷。微网与上级电网之间有购电和售电通道微网与天然气网之间有购气通道。储能电站则独立于各个微网之外但每个微网都可以通过联络线向储能电站充电或放电。这个拓扑结构是绝大多数冷热电多微网文献里的标准设定画图很容易但建模时要特别注意储能电站是共享的因此会有一个总容量约束也就是说各微网同一时刻向储能电站充放的总功率不能超过储能电站的额定功率总能量不能超过额定容量。这个约束是体现共享二字的灵魂代码里很多容易出错的地方都出在这里。2.2 储能电站的建模细节储能电站的建模是整个系统的核心因为它跨两层容量上层决策和充放电策略下层决策。在 MATLAB 代码里上层给储能电站的变量是额定功率 P_ess 和额定容量 E_ess下层则在每个时段给出充电功率 P_ch(t) 和放电功率 P_dis(t)。储能电站最重要的约束除了功率上下限就是电量 SOC荷电状态的递推关系SOC(t1) SOC(t) (eta_ch * P_ch(t) - P_dis(t) / eta_dis) * dt / E_ess这里 etta_ch 是充电效率eta_dis 是放电效率dt 是单个时段长度单位小时E_ess 是储能额定容量。很多人第一次写这个递推关系时容易漏掉效率系数或者把充电效率直接乘在 SOC 上导致储能电站凭空多出能量。另外别忘了全天最后一个时段要回到初始 SOC也就是 SOC(T) SOC(0)否则储能电站一天下来私吞了一部分电量运行成本会被严重低估。还有一个细节充放电功率不能同时为正这个约束在优化里叫互补约束。在 YALMIP 里我习惯引入二元变量或者直接用 big-M 约束来实现具体写法后面讲代码时会给出。2.3 冷热电联供设备的建模冷热电多微网里的设备建模说穿了就是能量转换效率和容量约束两个东西的组合。我把常用设备整理成一张表方便你对照着写代码设备类型输入输出关键参数微型燃气轮机 MT天然气电 高温烟气余热发电效率 eta_mt热电比燃气锅炉 GB天然气热热效率 eta_gb电制冷机 EC电冷制冷效率 COP_ec吸收式制冷机 AC热余热冷制冷效率 COP_ac光伏 PV太阳辐照电额定容量、预测出力系数微燃机是冷热电联供系统的核心设备它同时产生电和热余热可以一部分直接供给热负荷另一部分驱动吸收式制冷机来制冷。代码里最常写的电功率平衡方程是P_pv(t) P_mt(t) P_ess_dis(t) - P_ess_ch(t) P_buy(t) P_elec_load(t) P_ec(t)热功率平衡是P_mt_recover(t) P_gb(t) P_heat_load(t) P_ac(t)冷功率平衡是P_ec(t) * COP_ec P_ac(t) * COP_ac P_cool_load(t)注意冷负荷那里我用的是大于等于因为实际工程中允许适当少供冷或者用柔性负荷来缓解供冷压力。如果你想严格满足冷负荷改成等号即可但费用会相应高一些。2.4 数据准备负荷曲线和典型场景这个项目里最容易让代码跑飞的就是数据。光伏出力和负荷曲线必须提前准备好常见做法是对一年的历史数据做 k-means 聚类聚出几个典型日场景一般取 4 个季或 3~6 个典型日每个典型日取 24 小时的时序数据。聚类逻辑很简单把 365 天的 24 维负荷向量或者光伏向量作为样本用 MATLAB 自带的 kmeans 函数聚类再对每类内部取平均作为典型日。这样既保留了季节特征又把全年计算量压缩到可以接受的范围。我项目里最终取了 4 个典型日场景分给春夏秋冬效果不错。对了单位统一非常关键。我曾经在代码里把光伏出力写成 kW储能容量写成 kWh结果约束检查时怎么都不对。建议所有功率统一用 MW能量统一用 MWh投资成本统一用万元或百万元不然数值尺度差了几个量级求解器精度会大打折扣。3. 双层模型的数学结构上层定容量下层定运行3.1 上层配置模型等年值投资成本最小化上层优化要解决的问题是储能电站建多大、每个微网的微燃机配多大、锅炉配多大、制冷机配多大。目标函数通常写成等年值总成本最小min C_inv C_run C_maintC_inv 是设备投资的等年值成本计算方法是把总建设成本乘以一个资本回收系数CRF r * (1r)^n / ((1r)^n - 1)其中 r 是贴现率常取 0.06n 是设备寿命储能一般取 10 年微燃机 15 年。然后把各自的单位容量投资成本乘以容量再乘以 CRF就得到每年的投资分摊成本。这里特别强调等年值是因为不能把一次性投资直接和目标函数里的单年运行成本相加否则决策结果会严重偏向少投资、多运行。上层决策变量就是各设备的容量包括储能电站额定功率 P_ess 和额定容量 E_ess每个微网的微燃机容量 P_mt(i)每个微网的燃气锅炉容量 P_gb(i)每个微网的电制冷机容量 P_ec(i)每个微网的吸收式制冷机容量 P_ac(i)这些都是连续变量少数文献里会写成整数容量档位但在 MATLAB 实现时先按连续变量处理更稳妥。3.2 下层运行调度模型运行成本最小化给定一套容量配置之后下层要做的是在每个典型日 24 小时尺度上安排各设备的出力目标函数是典型日的运行成本最小化min C_fuel C_buy - C_sellC_fuel 是购买天然气的费用C_buy 是从上级电网购电的费用C_sell 是向电网售电的收入。电价采用分时电价峰值时段电价高低谷时段电价低这样储能电站才有套利空间。燃气价格一般是常量乘上天然气低热值和设备耗气量即可得到燃料成本。下层决策变量是各个设备在每个时段的出力以及储能电站的充放电功率、购电售电功率。下层约束除了前面说的功率平衡约束、设备出力上下限、储能 SOC 约束之外还有一个很容易漏掉的东西——向上级电网购电的功率上限。在大多数项目设定里多微网与上级电网的交换功率是受限的通常限制为联络线容量常见值在 1~2 MW 之间。这个约束不写模型就会在极端负荷时段出现想买多少买多少的不合理结果。3.3 两层之间到底怎么耦合双层模型的核心难点在于两层不是独立的而是互相依赖的。上层给下层一套容量参数比如储能 E_ess 20 MWh下层拿着这个容量去算日运行成本下层算出来的运行成本反过来成为上层目标函数里 C_run 的输入。所以双层问题本质上是一个领导者和跟随者的博弈问题——上层决定策略下层在这个策略下做出最优反应。用代码术语来说上层每一次迭代都要调用一次下层求解器。我见过很多人把双层优化理解为先优化上层再优化下层各自算各自的这是不对的。正确的关系是下层优化问题作为约束嵌套在上层优化里面上层任何一个候选方案都必须通过下层优化来评估其真实成本。这就解释了为什么双层优化的计算量远比单层大——每评估一次 PSO 里的粒子位置就要完整跑一遍下层优化。3.4 关于典型场景权重的一点补充如果聚了 4 个典型日每个典型日代表的天数不同那么下层运行成本要按场景天数加权求和再与上层等年值投资成本相加。比如春季典型日代表 92 天那么该场景的运行成本要乘以 92/365 的权重。这个权重若不处理优化结果会被某个代表性弱的场景带偏。代码里可以在上层目标函数中用循环累加多个典型日的运行成本实现起来并不复杂。4. MATLAB 实现路径我选的是 PSO 嵌套 YALMIP/CPLEX4.1 两条主流技术路线的对比双层优化在 MATLAB 里的实现方案其实就两条主流路线。第一条是 KKT 条件法把下层优化问题用其 KKT 条件或强对偶条件替换把双层问题转成单层 MPEC 问题然后用求解器直接解。优点是理论上严谨不用迭代但缺点非常明显——下层问题里如果有整数变量比如储能充放电的 0/1 变量KKT 条件的推导会变得极其繁琐而且 MPEC 问题的互补约束对求解器数值稳定性要求很高YALMIP 处理起来很吃力。我在实际实验中发现即便把小问题转成单层求解效率也低得吓人。第二条是启发式嵌套法上层用种群算法PSO、GA 等每个个体代表一套容量配置交给下层求解器YALMIP 调用 CPLEX/Gurobi计算运行成本再把成本返回上层作为适应度值。这种方式实现简单、思路直白工程上很稳定缺点是每次迭代都要调用多次下层求解计算量大。但对于典型日数量不超过 5 个、微网数量不超过 5 个的场景这个方案完全能接受。我个人强烈推荐第二条路线因为可解释性强出了问题也好排查。下面是它的整体流程初始化 PSO 种群每个粒子 一套设备容量 for 迭代次数: for 每个粒子: 把容量参数传给下层 下层用 YALMIP 求解多场景日运行优化 计算运行成本并返回 上层按 PSO 规则计算适应度投资成本 运行成本 更新粒子速度和位置 输出最优容量配置4.2 PSO 参数的设置心得PSO 在这里是上层寻优器参数设置直接影响收敛效果。我的实测经验是种群规模 40~60迭代次数 60~80惯性权重从 0.9 线性递减到 0.4学习因子 c1 c2 2。这些参数不是我拍脑袋定的是做过对比实验的——惯性权重不减的话粒子容易在后期来回震荡种群太小容易陷入局部最优。另外上下层变量的缩放很重要。存储容量的取值范围我建议设置成 [5, 50] MWh功率范围 [1, 10] MW微燃机容量 [0.5, 5] MW。PSO 对取值范围很敏感范围设得太大粒子分布稀疏收敛慢范围设得太小可能漏掉最优解。如果你不确定最优容量大概在什么量级可以先跑一次单层粗略优化找感觉。4.3 下层 YALMIP 求解器的调用骨架下层是标准的 LP/MILP 问题用 YALMIP 建模非常高效。下面这段代码是下层函数的核心骨架我把注释写在代码里方便你直接改来用function [run_cost, detail] lower_level(cap, load_data, price) % cap 是上层传下来的容量向量 % load_data 是 24h 负荷与光伏数据 % price 是分时电价向量 T 24; E_ess cap(1); P_ess cap(2); % 储能容量和功率 P_mt_max cap(3); % 微燃机容量 yalmip(clear); % 决策变量 P_mt sdpvar(T,1); % 微燃机出力 P_ch sdpvar(T,1); % 储能充电功率 P_dis sdpvar(T,1); % 储能放电功率 P_buy sdpvar(T,1); % 购电 P_sell sdpvar(T,1); % 售电 P_gb sdpvar(T,1); % 锅炉出力 P_ec sdpvar(T,1); % 电制冷机输入电功率 P_ac sdpvar(T,1); % 吸收式制冷机输入热功率 SOC sdpvar(T1,1); % 荷电状态 u_ch binvar(T,1); % 充电状态标志 u_dis binvar(T,1); % 放电状态标志 Constraints []; % 储能SOC递推eta_ch0.95, eta_dis0.95 dt 1; for t 1:T Constraints [Constraints, SOC(t1) SOC(t) ... (0.95*P_ch(t) - P_dis(t)/0.95)*dt/E_ess]; end Constraints [Constraints, SOC(1) 0.2, SOC(T1) 0.2]; Constraints [Constraints, 0.1 SOC 0.9]; Constraints [Constraints, 0 P_ch P_ess*u_ch]; Constraints [Constraints, 0 P_dis P_ess*u_dis]; Constraints [Constraints, u_ch u_dis 1]; % 电功率平衡 P_pv load_data.pv; P_load load_data.elec; Constraints [Constraints, P_pv P_mt P_dis P_buy ... P_load P_ec P_ch P_sell]; % 热/冷平衡、设备上下限等略... % 目标函数购电费用 - 售电收入 燃气费用 gas_price 2.5; % 元/m3 gas_lhv 9.7; % kWh/m3 Objective sum(price.buy .* P_buy - price.sell .* P_sell ... gas_price/gas_lhv * (P_mt P_gb)); options sdpsettings(solver,cplex,verbose,0, ... showprogress,0); optimize(Constraints, Objective, options); run_cost value(Objective); detail.P_mt value(P_mt); detail.P_ch value(P_ch); detail.P_dis value(P_dis); end这段代码里有一个很微妙的点SOC 递推公式中我把效率放在了充电功率项上放电功率除以效率这在物理上是合理的——充进去 1 kWh 电量扣掉充电损耗后只有 0.95 kWh 进入电池放出来时为了输出 1 kWh电池内部要消耗 1/0.95 kWh。如果你把效率位置放反储能运行成本会算错优化结果直接偏离。4.4 上层 PSO 主循环代码骨架上层主循环我单独拿一个脚本写这样日志、断点续跑都比较方便load load_data.mat; load price.mat; nPop 50; maxIter 60; w (iter) 0.9 - 0.5*iter/maxIter; c1 2; c2 2; lb [5, 1, 0.5]; ub [50, 10, 5]; % 容量上下限 pos repmat(lb, nPop, 1) rand(nPop, length(lb)).*(repmat(ub-lb, nPop, 1)); vel zeros(nPop, length(lb)); pbest pos; pbest_fit inf(nPop,1); for iter 1:maxIter for i 1:nPop cap pos(i,:); [run_cost, ~] lower_level(cap, load_data, price); inv_cost cal_invest(cap); % 等年值投资成本 fitness inv_cost run_cost; if fitness pbest_fit(i) pbest_fit(i) fitness; pbest(i,:) cap; end end [gbest_fit, idx] min(pbest_fit); gbest pbest(idx,:); for i 1:nPop vel(i,:) w(iter)*vel(i,:) ... c1*rand*(pbest(i,:) - pos(i,:)) ... c2*rand*(gbest - pos(i,:)); pos(i,:) pos(i,:) vel(i,:); pos(i,:) min(max(pos(i,:), lb), ub); % 边界约束 end fprintf(iter %d, best cost %.2f\n, iter, gbest_fit); end这里有个细节值得注意cal_invest(cap)的计算要单独写成一个函数它内部需要按设备寿命分别计算 CRF再换算成等年值。如果投资成本和运行成本不是一个时间尺度直接相加会让 PSO 的适应度函数数值失真。5. 代码落地中我踩过的 9 个坑以及排查思路5.1 YALMIP 与求解器的安装问题很多新手第一次跑这类代码问题不在模型而在环境。MATLAB 本身要先装好网上教程很多但装完以后经常出现由于找不到 msvcp140.dll 无法继续执行代码这类系统级报错这是 MATLAB 依赖的 Visual C 运行库缺失导致的去微软官网下载 vc_redist.x64.exe 装上就好了。这属于环境配置不解决的话连 MATLAB 都启动不了。然后是 YALMIP 和求解器。YALMIP 直接去 GitHub 下载并添加路径即可不用安装。求解器方面CPLEX 和 Gurobi 都有学术授权如果没有也可以用开源的 GLPK 或 MATLAB 自带的 linprog。注意 YALMIP 里指定求解器的方法是在sdpsettings中设置solver,cplex但前提是 MATLAB 能通过hasName(cplex)检测到它。我踩过最深的一个坑是装完 CPLEX 之后YALMIP 始终找不到 cplex后来发现问题出在 MATALB 的 java 路径没有指向 CPLEX 的 jar 包。这个问题在网上搜cplex yalmip not found基本能找到解决方法。5.2 数值单位不一致导致的假收敛前面提到过单位统一的问题。这里我再举个具体例子如果光伏出力数据是 kW 而储能容量是 MWh约束写出来就有 1000 倍的差异CPLEX 在数值求解时容易因为矩阵条件数过大而报numerical difficulties或者干脆给出一个很离谱的解。我在项目初期就吃过一次亏算出来的储能容量只有 0.02 MWh明显不合理排查半天发现是光伏数据用了 kW储能单位用了 MWh。所以做这类项目一定要在数据读入之后先做一次单位归一化检查。5.3 SOC 循环约束导致的下层问题不可行最容易让下层模型无解的坑就是 SOC 循环约束。如果设定 SOC(T) SOC(0) 0.2但在 PV 出力很小、电价很高的连续时段里储能既要满足负荷又要回到初始值很可能导致下层无可行解。应对方法有两个一是允许 SOC 终端值与初始值在一定范围内波动比如 SOC(T) SOC(0) - 0.05二是在目标函数里加一个很小的终端 SOC 偏差惩罚项。我自己偏好第二种因为代码改动小物理含义也明确。5.4 充放电互补约束的写法在 4.3 代码骨架里我用了两个二元变量 u_ch 和 u_dis 来保证不同时充放电。这个约束如果不加CPLEX 算出来的解可能出现既充电又放电的套利怪象——虽然目标函数不会因此变得更小但结果没有工程意义。还有另一种写法是在约束里直接加上 P_ch * P_dis 0但那会让问题变成非线性YALMIP 不太好处理。用二元变量是最常用的标准做法。5.5 PSO 迭代震荡我遇到过 PSO 跑了 40 代还不收敛适应度曲线像心电图一样上下波动。排查之后发现原因有两个一是惯性权重没有递减导致粒子后期速度依然很大二是下层求解器在极少数情况下会返回一个病态值比如某天场景权重很小但运行成本异常高把整体适应度曲线拉出尖峰。后来我做了两个改进惯性权重改为随迭代次数线性递减对种群中每个粒子的越界位置做随机重置而不是简单截断。修改之后收敛曲线明显平稳大概在 25 代左右就进入稳定状态。5.6 场景权重的处理疏漏如果你的模型聚了多个典型日每个场景的天数权重不同下层目标函数必须乘上对应权重后再累加。我见过有人把 4 个场景的运行成本直接相加没有加权导致夏季数据占比过低的场景被忽略储能容量整体偏小。这里的权重本质上代表该典型日在一年中出现的概率写代码时只需在外面套一层循环乘上系数即可但这个循环的位置很容易放错。我的做法是把每个场景的运行成本在上层函数内部单独算一遍乘以天数权重后再累加确保逻辑清晰。5.7 联络线功率上限缺失导致的异常购电有些论文在建模时会把购电上限写在约束里但很多人抄代码时容易忽略。实际系统中多微网与电网的连接容量是有限的一般不会超过几 MW。如果不在模型里加P_buy P_buy_max优化器会让系统在高峰时段大量购电甚至超过电网实际供电能力算出来的运行成本偏低但没有工程可行性。更合理的做法是同时限制购电和售电不能同时发生道理同储能互补约束。5.8 计算时间失控PSO 嵌套 YALMIP 最大的痛点是计算时间。我最初把所有约束写成逐时段的循环10 个粒子跑 30 代就已经消耗了快一个小时。后来我把能向量化的约束改成矩阵形式比如储能 SOC 递推关系用矩阵直接表示避免在约束集合里写 24 次 append计算时间缩短了大概 70%。如果你觉得代码还是慢优先检查约束构建部分是否用了大量 for 循环YALMIP 的 constraint append 是非常耗时的操作。5.9 结果检查的三看原则代码跑完之后一定要做三件事一看容量结果是否符合工程直觉储能太大还是太小二看典型日的 SOC 曲线是否平滑、是否频繁触顶触底三看每个微网的电功率平衡是否真的闭合。我通常会在输出函数里额外打印各时段的功率平衡残差如果残差超过 1e-6说明约束写错了或者求解器精度不够需要回调。6. 从运行结果到研究报告补充三个维度的对比分析6.1 基础结果表的输出优化完成后你需要把结果整理成表格方便在论文或报告里直接引用。我的习惯是至少输出三张表配置结果表各设备的容量、储能电站额定功率和容量、等年值投资成本运行结果表各典型日的购电量、购气量、运行成本、峰值负荷对比表有共享储能 vs 无共享储能的总成本双层优化 vs 单层优化的结果差异对比表尤其重要它能证明双层优化模型的必要性。我实验中得到的一组典型结果是共享储能模式下储能利用效率比自建储能高约 35%系统年总成本下降约 12.8%。当然不同项目参数不同数值仅供参考但方向是一致的。6.2 敏感性分析怎么做审稿人或导师几乎必定会问如果储能单位成本变了结果还成立吗。所以你需要补一组敏感性分析把储能的单位容量成本从 800 元/kWh 依次调到 1000、1200、1400 元/kWh再看最优配置怎么变化。实现方式很简单外层 PSO 不变把投资成本函数里的单位成本参数换成数组重新跑一遍即可。这样一来你得到的不再是单个最优解而是一组成本-容量的对应关系说服力强很多。6.3 三个有意义的扩展方向如果这个项目还有后续我建议从三个方向延伸需求响应把冷热负荷从刚性改成弹性加入可中断负荷和价格型需求响应下层模型会变成更复杂的混合整数二次规划碳交易机制在目标函数里加入阶梯式碳交易成本关注二氧化碳排放约束对容量配置的影响多主体博弈把各微网从受系统统一调度变为独立利益主体用纳什均衡或主从博弈框架表征微网之间的策略互动双层结构会升级为多层或 MPEC 问题这三个方向其实是综合能源系统领域目前最热的课题如果你做的是硕士或博士课题完全可以在当前代码基础上继续迭代。就我个人经验来说这类双层优化配置项目最难的从来不是某一段代码怎么写而是你能不能把哪个变量放上层、哪个变量放下层、两层之间靠什么耦合想清楚。这个问题的答案直接决定了模型行为的合理性也决定了论文评审老师会不会追问。我一开始就是急着写代码结果被导师一句话问住你凭什么把储能容量放在上层如果运行方同时也是投资方这个双层结构还成立吗这个问号逼着我回去把系统物理背景又捋了三遍。最终落在纸面上的模型是所有思考过程的自然结果而不是代码本身。最后补一个小技巧不管用什么求解器强烈建议在代码里固定随机种子rng(2024)因为 PSO 是随机算法不固定种子的话每次运行结果都有波动写报告时数据对不上会很痛苦。我用这个办法保证了所有对比实验都在同一基准下进行也让后续的敏感性分析结果更可信。 储能电站服务下冷热电多微网系统双层优化配置的MATLAB实现
1. 为什么储能电站服务下的多微网系统天然需要双层优化配置先交代一下背景。我最近一直在做储能电站相关的项目客户那边给的课题是基于储能电站服务的冷热电多微网系统双层优化配置要求用 MATLAB 实现给出代码框架和计算结果。说实话刚拿到这个题目的时候我也愣了一下——双层优化配置这个名词听起来很高大上但拆开来看就是两件事第一给储能电站和各微网里的设备定容量第二在容量定好之后安排各个微网在一天的运行计划。前者是建多大后者是怎么用。这两件事一旦放在同一个模型里单层优化是做不到的必须分层来解决。很多人第一次接触这类问题上来就想把容量配置和运行调度写进同一个目标函数里一次性求解。这种做法不是不行而是会导致一个很尴尬的结果设备容量和运行策略在模型里会互相打架。比如说储能容量定得大运行阶段就可以多充多放、降低购电成本但储能本身的投资成本又很高如果容量定得小投资省了可运行阶段又买高价电。这两个目标属于不同的时间尺度容量配置是年尺度决策运行调度是小时尺度决策硬揉在一起模型要么不收敛要么算出来的容量毫无工程意义。所以我在这篇文章里想把我从建模到 MATLAB 代码落地的全过程梳理一遍包括双层模型怎么拆、上层和下层分别优化什么、两层之间怎么迭代耦合以及代码里最容易踩的坑。如果你正在做微电网、综合能源系统、共享储能或者CCHP 多微网方向的课题这篇内容基本可以当成一份可抄作业的路线图来用。先明确一个概念什么叫储能电站服务。在传统的微网系统里储能是每个微网自建自用的容量小、利用率低而且投资主体和运行主体是同一个。但随着共享储能模式出现储能电站变成了一个独立的服务商它向多个微网提供充电、放电服务容量决策在一个主体手里运行调度却在各个微网手里。这种投资运营分离的结构天然就是双层优化的用武之地——上层做储能电站和各微网设备的容量配置下层做多微网系统的日运行调度。共享储能的意义在于降低储能成本、提高设备利用率而双层优化配置就是把这个意义真正落到数学建模和代码里的关键手段。2. 系统拓扑与设备建模先把冷热电多微网画清楚2.1 网络拓扑怎么搭在动手写代码之前我习惯先把系统拓扑在草稿纸上画出来。这个项目对应的结构大致是这样的系统包含若干个微网常见的做法是取 3 个典型微网每个微网内部有光伏、微型燃气轮机、燃气锅炉、电制冷机、吸收式制冷机这几类设备承担自己的电负荷、热负荷和冷负荷。微网与上级电网之间有购电和售电通道微网与天然气网之间有购气通道。储能电站则独立于各个微网之外但每个微网都可以通过联络线向储能电站充电或放电。这个拓扑结构是绝大多数冷热电多微网文献里的标准设定画图很容易但建模时要特别注意储能电站是共享的因此会有一个总容量约束也就是说各微网同一时刻向储能电站充放的总功率不能超过储能电站的额定功率总能量不能超过额定容量。这个约束是体现共享二字的灵魂代码里很多容易出错的地方都出在这里。2.2 储能电站的建模细节储能电站的建模是整个系统的核心因为它跨两层容量上层决策和充放电策略下层决策。在 MATLAB 代码里上层给储能电站的变量是额定功率 P_ess 和额定容量 E_ess下层则在每个时段给出充电功率 P_ch(t) 和放电功率 P_dis(t)。储能电站最重要的约束除了功率上下限就是电量 SOC荷电状态的递推关系SOC(t1) SOC(t) (eta_ch * P_ch(t) - P_dis(t) / eta_dis) * dt / E_ess这里 etta_ch 是充电效率eta_dis 是放电效率dt 是单个时段长度单位小时E_ess 是储能额定容量。很多人第一次写这个递推关系时容易漏掉效率系数或者把充电效率直接乘在 SOC 上导致储能电站凭空多出能量。另外别忘了全天最后一个时段要回到初始 SOC也就是 SOC(T) SOC(0)否则储能电站一天下来私吞了一部分电量运行成本会被严重低估。还有一个细节充放电功率不能同时为正这个约束在优化里叫互补约束。在 YALMIP 里我习惯引入二元变量或者直接用 big-M 约束来实现具体写法后面讲代码时会给出。2.3 冷热电联供设备的建模冷热电多微网里的设备建模说穿了就是能量转换效率和容量约束两个东西的组合。我把常用设备整理成一张表方便你对照着写代码设备类型输入输出关键参数微型燃气轮机 MT天然气电 高温烟气余热发电效率 eta_mt热电比燃气锅炉 GB天然气热热效率 eta_gb电制冷机 EC电冷制冷效率 COP_ec吸收式制冷机 AC热余热冷制冷效率 COP_ac光伏 PV太阳辐照电额定容量、预测出力系数微燃机是冷热电联供系统的核心设备它同时产生电和热余热可以一部分直接供给热负荷另一部分驱动吸收式制冷机来制冷。代码里最常写的电功率平衡方程是P_pv(t) P_mt(t) P_ess_dis(t) - P_ess_ch(t) P_buy(t) P_elec_load(t) P_ec(t)热功率平衡是P_mt_recover(t) P_gb(t) P_heat_load(t) P_ac(t)冷功率平衡是P_ec(t) * COP_ec P_ac(t) * COP_ac P_cool_load(t)注意冷负荷那里我用的是大于等于因为实际工程中允许适当少供冷或者用柔性负荷来缓解供冷压力。如果你想严格满足冷负荷改成等号即可但费用会相应高一些。2.4 数据准备负荷曲线和典型场景这个项目里最容易让代码跑飞的就是数据。光伏出力和负荷曲线必须提前准备好常见做法是对一年的历史数据做 k-means 聚类聚出几个典型日场景一般取 4 个季或 3~6 个典型日每个典型日取 24 小时的时序数据。聚类逻辑很简单把 365 天的 24 维负荷向量或者光伏向量作为样本用 MATLAB 自带的 kmeans 函数聚类再对每类内部取平均作为典型日。这样既保留了季节特征又把全年计算量压缩到可以接受的范围。我项目里最终取了 4 个典型日场景分给春夏秋冬效果不错。对了单位统一非常关键。我曾经在代码里把光伏出力写成 kW储能容量写成 kWh结果约束检查时怎么都不对。建议所有功率统一用 MW能量统一用 MWh投资成本统一用万元或百万元不然数值尺度差了几个量级求解器精度会大打折扣。3. 双层模型的数学结构上层定容量下层定运行3.1 上层配置模型等年值投资成本最小化上层优化要解决的问题是储能电站建多大、每个微网的微燃机配多大、锅炉配多大、制冷机配多大。目标函数通常写成等年值总成本最小min C_inv C_run C_maintC_inv 是设备投资的等年值成本计算方法是把总建设成本乘以一个资本回收系数CRF r * (1r)^n / ((1r)^n - 1)其中 r 是贴现率常取 0.06n 是设备寿命储能一般取 10 年微燃机 15 年。然后把各自的单位容量投资成本乘以容量再乘以 CRF就得到每年的投资分摊成本。这里特别强调等年值是因为不能把一次性投资直接和目标函数里的单年运行成本相加否则决策结果会严重偏向少投资、多运行。上层决策变量就是各设备的容量包括储能电站额定功率 P_ess 和额定容量 E_ess每个微网的微燃机容量 P_mt(i)每个微网的燃气锅炉容量 P_gb(i)每个微网的电制冷机容量 P_ec(i)每个微网的吸收式制冷机容量 P_ac(i)这些都是连续变量少数文献里会写成整数容量档位但在 MATLAB 实现时先按连续变量处理更稳妥。3.2 下层运行调度模型运行成本最小化给定一套容量配置之后下层要做的是在每个典型日 24 小时尺度上安排各设备的出力目标函数是典型日的运行成本最小化min C_fuel C_buy - C_sellC_fuel 是购买天然气的费用C_buy 是从上级电网购电的费用C_sell 是向电网售电的收入。电价采用分时电价峰值时段电价高低谷时段电价低这样储能电站才有套利空间。燃气价格一般是常量乘上天然气低热值和设备耗气量即可得到燃料成本。下层决策变量是各个设备在每个时段的出力以及储能电站的充放电功率、购电售电功率。下层约束除了前面说的功率平衡约束、设备出力上下限、储能 SOC 约束之外还有一个很容易漏掉的东西——向上级电网购电的功率上限。在大多数项目设定里多微网与上级电网的交换功率是受限的通常限制为联络线容量常见值在 1~2 MW 之间。这个约束不写模型就会在极端负荷时段出现想买多少买多少的不合理结果。3.3 两层之间到底怎么耦合双层模型的核心难点在于两层不是独立的而是互相依赖的。上层给下层一套容量参数比如储能 E_ess 20 MWh下层拿着这个容量去算日运行成本下层算出来的运行成本反过来成为上层目标函数里 C_run 的输入。所以双层问题本质上是一个领导者和跟随者的博弈问题——上层决定策略下层在这个策略下做出最优反应。用代码术语来说上层每一次迭代都要调用一次下层求解器。我见过很多人把双层优化理解为先优化上层再优化下层各自算各自的这是不对的。正确的关系是下层优化问题作为约束嵌套在上层优化里面上层任何一个候选方案都必须通过下层优化来评估其真实成本。这就解释了为什么双层优化的计算量远比单层大——每评估一次 PSO 里的粒子位置就要完整跑一遍下层优化。3.4 关于典型场景权重的一点补充如果聚了 4 个典型日每个典型日代表的天数不同那么下层运行成本要按场景天数加权求和再与上层等年值投资成本相加。比如春季典型日代表 92 天那么该场景的运行成本要乘以 92/365 的权重。这个权重若不处理优化结果会被某个代表性弱的场景带偏。代码里可以在上层目标函数中用循环累加多个典型日的运行成本实现起来并不复杂。4. MATLAB 实现路径我选的是 PSO 嵌套 YALMIP/CPLEX4.1 两条主流技术路线的对比双层优化在 MATLAB 里的实现方案其实就两条主流路线。第一条是 KKT 条件法把下层优化问题用其 KKT 条件或强对偶条件替换把双层问题转成单层 MPEC 问题然后用求解器直接解。优点是理论上严谨不用迭代但缺点非常明显——下层问题里如果有整数变量比如储能充放电的 0/1 变量KKT 条件的推导会变得极其繁琐而且 MPEC 问题的互补约束对求解器数值稳定性要求很高YALMIP 处理起来很吃力。我在实际实验中发现即便把小问题转成单层求解效率也低得吓人。第二条是启发式嵌套法上层用种群算法PSO、GA 等每个个体代表一套容量配置交给下层求解器YALMIP 调用 CPLEX/Gurobi计算运行成本再把成本返回上层作为适应度值。这种方式实现简单、思路直白工程上很稳定缺点是每次迭代都要调用多次下层求解计算量大。但对于典型日数量不超过 5 个、微网数量不超过 5 个的场景这个方案完全能接受。我个人强烈推荐第二条路线因为可解释性强出了问题也好排查。下面是它的整体流程初始化 PSO 种群每个粒子 一套设备容量 for 迭代次数: for 每个粒子: 把容量参数传给下层 下层用 YALMIP 求解多场景日运行优化 计算运行成本并返回 上层按 PSO 规则计算适应度投资成本 运行成本 更新粒子速度和位置 输出最优容量配置4.2 PSO 参数的设置心得PSO 在这里是上层寻优器参数设置直接影响收敛效果。我的实测经验是种群规模 40~60迭代次数 60~80惯性权重从 0.9 线性递减到 0.4学习因子 c1 c2 2。这些参数不是我拍脑袋定的是做过对比实验的——惯性权重不减的话粒子容易在后期来回震荡种群太小容易陷入局部最优。另外上下层变量的缩放很重要。存储容量的取值范围我建议设置成 [5, 50] MWh功率范围 [1, 10] MW微燃机容量 [0.5, 5] MW。PSO 对取值范围很敏感范围设得太大粒子分布稀疏收敛慢范围设得太小可能漏掉最优解。如果你不确定最优容量大概在什么量级可以先跑一次单层粗略优化找感觉。4.3 下层 YALMIP 求解器的调用骨架下层是标准的 LP/MILP 问题用 YALMIP 建模非常高效。下面这段代码是下层函数的核心骨架我把注释写在代码里方便你直接改来用function [run_cost, detail] lower_level(cap, load_data, price) % cap 是上层传下来的容量向量 % load_data 是 24h 负荷与光伏数据 % price 是分时电价向量 T 24; E_ess cap(1); P_ess cap(2); % 储能容量和功率 P_mt_max cap(3); % 微燃机容量 yalmip(clear); % 决策变量 P_mt sdpvar(T,1); % 微燃机出力 P_ch sdpvar(T,1); % 储能充电功率 P_dis sdpvar(T,1); % 储能放电功率 P_buy sdpvar(T,1); % 购电 P_sell sdpvar(T,1); % 售电 P_gb sdpvar(T,1); % 锅炉出力 P_ec sdpvar(T,1); % 电制冷机输入电功率 P_ac sdpvar(T,1); % 吸收式制冷机输入热功率 SOC sdpvar(T1,1); % 荷电状态 u_ch binvar(T,1); % 充电状态标志 u_dis binvar(T,1); % 放电状态标志 Constraints []; % 储能SOC递推eta_ch0.95, eta_dis0.95 dt 1; for t 1:T Constraints [Constraints, SOC(t1) SOC(t) ... (0.95*P_ch(t) - P_dis(t)/0.95)*dt/E_ess]; end Constraints [Constraints, SOC(1) 0.2, SOC(T1) 0.2]; Constraints [Constraints, 0.1 SOC 0.9]; Constraints [Constraints, 0 P_ch P_ess*u_ch]; Constraints [Constraints, 0 P_dis P_ess*u_dis]; Constraints [Constraints, u_ch u_dis 1]; % 电功率平衡 P_pv load_data.pv; P_load load_data.elec; Constraints [Constraints, P_pv P_mt P_dis P_buy ... P_load P_ec P_ch P_sell]; % 热/冷平衡、设备上下限等略... % 目标函数购电费用 - 售电收入 燃气费用 gas_price 2.5; % 元/m3 gas_lhv 9.7; % kWh/m3 Objective sum(price.buy .* P_buy - price.sell .* P_sell ... gas_price/gas_lhv * (P_mt P_gb)); options sdpsettings(solver,cplex,verbose,0, ... showprogress,0); optimize(Constraints, Objective, options); run_cost value(Objective); detail.P_mt value(P_mt); detail.P_ch value(P_ch); detail.P_dis value(P_dis); end这段代码里有一个很微妙的点SOC 递推公式中我把效率放在了充电功率项上放电功率除以效率这在物理上是合理的——充进去 1 kWh 电量扣掉充电损耗后只有 0.95 kWh 进入电池放出来时为了输出 1 kWh电池内部要消耗 1/0.95 kWh。如果你把效率位置放反储能运行成本会算错优化结果直接偏离。4.4 上层 PSO 主循环代码骨架上层主循环我单独拿一个脚本写这样日志、断点续跑都比较方便load load_data.mat; load price.mat; nPop 50; maxIter 60; w (iter) 0.9 - 0.5*iter/maxIter; c1 2; c2 2; lb [5, 1, 0.5]; ub [50, 10, 5]; % 容量上下限 pos repmat(lb, nPop, 1) rand(nPop, length(lb)).*(repmat(ub-lb, nPop, 1)); vel zeros(nPop, length(lb)); pbest pos; pbest_fit inf(nPop,1); for iter 1:maxIter for i 1:nPop cap pos(i,:); [run_cost, ~] lower_level(cap, load_data, price); inv_cost cal_invest(cap); % 等年值投资成本 fitness inv_cost run_cost; if fitness pbest_fit(i) pbest_fit(i) fitness; pbest(i,:) cap; end end [gbest_fit, idx] min(pbest_fit); gbest pbest(idx,:); for i 1:nPop vel(i,:) w(iter)*vel(i,:) ... c1*rand*(pbest(i,:) - pos(i,:)) ... c2*rand*(gbest - pos(i,:)); pos(i,:) pos(i,:) vel(i,:); pos(i,:) min(max(pos(i,:), lb), ub); % 边界约束 end fprintf(iter %d, best cost %.2f\n, iter, gbest_fit); end这里有个细节值得注意cal_invest(cap)的计算要单独写成一个函数它内部需要按设备寿命分别计算 CRF再换算成等年值。如果投资成本和运行成本不是一个时间尺度直接相加会让 PSO 的适应度函数数值失真。5. 代码落地中我踩过的 9 个坑以及排查思路5.1 YALMIP 与求解器的安装问题很多新手第一次跑这类代码问题不在模型而在环境。MATLAB 本身要先装好网上教程很多但装完以后经常出现由于找不到 msvcp140.dll 无法继续执行代码这类系统级报错这是 MATLAB 依赖的 Visual C 运行库缺失导致的去微软官网下载 vc_redist.x64.exe 装上就好了。这属于环境配置不解决的话连 MATLAB 都启动不了。然后是 YALMIP 和求解器。YALMIP 直接去 GitHub 下载并添加路径即可不用安装。求解器方面CPLEX 和 Gurobi 都有学术授权如果没有也可以用开源的 GLPK 或 MATLAB 自带的 linprog。注意 YALMIP 里指定求解器的方法是在sdpsettings中设置solver,cplex但前提是 MATLAB 能通过hasName(cplex)检测到它。我踩过最深的一个坑是装完 CPLEX 之后YALMIP 始终找不到 cplex后来发现问题出在 MATALB 的 java 路径没有指向 CPLEX 的 jar 包。这个问题在网上搜cplex yalmip not found基本能找到解决方法。5.2 数值单位不一致导致的假收敛前面提到过单位统一的问题。这里我再举个具体例子如果光伏出力数据是 kW 而储能容量是 MWh约束写出来就有 1000 倍的差异CPLEX 在数值求解时容易因为矩阵条件数过大而报numerical difficulties或者干脆给出一个很离谱的解。我在项目初期就吃过一次亏算出来的储能容量只有 0.02 MWh明显不合理排查半天发现是光伏数据用了 kW储能单位用了 MWh。所以做这类项目一定要在数据读入之后先做一次单位归一化检查。5.3 SOC 循环约束导致的下层问题不可行最容易让下层模型无解的坑就是 SOC 循环约束。如果设定 SOC(T) SOC(0) 0.2但在 PV 出力很小、电价很高的连续时段里储能既要满足负荷又要回到初始值很可能导致下层无可行解。应对方法有两个一是允许 SOC 终端值与初始值在一定范围内波动比如 SOC(T) SOC(0) - 0.05二是在目标函数里加一个很小的终端 SOC 偏差惩罚项。我自己偏好第二种因为代码改动小物理含义也明确。5.4 充放电互补约束的写法在 4.3 代码骨架里我用了两个二元变量 u_ch 和 u_dis 来保证不同时充放电。这个约束如果不加CPLEX 算出来的解可能出现既充电又放电的套利怪象——虽然目标函数不会因此变得更小但结果没有工程意义。还有另一种写法是在约束里直接加上 P_ch * P_dis 0但那会让问题变成非线性YALMIP 不太好处理。用二元变量是最常用的标准做法。5.5 PSO 迭代震荡我遇到过 PSO 跑了 40 代还不收敛适应度曲线像心电图一样上下波动。排查之后发现原因有两个一是惯性权重没有递减导致粒子后期速度依然很大二是下层求解器在极少数情况下会返回一个病态值比如某天场景权重很小但运行成本异常高把整体适应度曲线拉出尖峰。后来我做了两个改进惯性权重改为随迭代次数线性递减对种群中每个粒子的越界位置做随机重置而不是简单截断。修改之后收敛曲线明显平稳大概在 25 代左右就进入稳定状态。5.6 场景权重的处理疏漏如果你的模型聚了多个典型日每个场景的天数权重不同下层目标函数必须乘上对应权重后再累加。我见过有人把 4 个场景的运行成本直接相加没有加权导致夏季数据占比过低的场景被忽略储能容量整体偏小。这里的权重本质上代表该典型日在一年中出现的概率写代码时只需在外面套一层循环乘上系数即可但这个循环的位置很容易放错。我的做法是把每个场景的运行成本在上层函数内部单独算一遍乘以天数权重后再累加确保逻辑清晰。5.7 联络线功率上限缺失导致的异常购电有些论文在建模时会把购电上限写在约束里但很多人抄代码时容易忽略。实际系统中多微网与电网的连接容量是有限的一般不会超过几 MW。如果不在模型里加P_buy P_buy_max优化器会让系统在高峰时段大量购电甚至超过电网实际供电能力算出来的运行成本偏低但没有工程可行性。更合理的做法是同时限制购电和售电不能同时发生道理同储能互补约束。5.8 计算时间失控PSO 嵌套 YALMIP 最大的痛点是计算时间。我最初把所有约束写成逐时段的循环10 个粒子跑 30 代就已经消耗了快一个小时。后来我把能向量化的约束改成矩阵形式比如储能 SOC 递推关系用矩阵直接表示避免在约束集合里写 24 次 append计算时间缩短了大概 70%。如果你觉得代码还是慢优先检查约束构建部分是否用了大量 for 循环YALMIP 的 constraint append 是非常耗时的操作。5.9 结果检查的三看原则代码跑完之后一定要做三件事一看容量结果是否符合工程直觉储能太大还是太小二看典型日的 SOC 曲线是否平滑、是否频繁触顶触底三看每个微网的电功率平衡是否真的闭合。我通常会在输出函数里额外打印各时段的功率平衡残差如果残差超过 1e-6说明约束写错了或者求解器精度不够需要回调。6. 从运行结果到研究报告补充三个维度的对比分析6.1 基础结果表的输出优化完成后你需要把结果整理成表格方便在论文或报告里直接引用。我的习惯是至少输出三张表配置结果表各设备的容量、储能电站额定功率和容量、等年值投资成本运行结果表各典型日的购电量、购气量、运行成本、峰值负荷对比表有共享储能 vs 无共享储能的总成本双层优化 vs 单层优化的结果差异对比表尤其重要它能证明双层优化模型的必要性。我实验中得到的一组典型结果是共享储能模式下储能利用效率比自建储能高约 35%系统年总成本下降约 12.8%。当然不同项目参数不同数值仅供参考但方向是一致的。6.2 敏感性分析怎么做审稿人或导师几乎必定会问如果储能单位成本变了结果还成立吗。所以你需要补一组敏感性分析把储能的单位容量成本从 800 元/kWh 依次调到 1000、1200、1400 元/kWh再看最优配置怎么变化。实现方式很简单外层 PSO 不变把投资成本函数里的单位成本参数换成数组重新跑一遍即可。这样一来你得到的不再是单个最优解而是一组成本-容量的对应关系说服力强很多。6.3 三个有意义的扩展方向如果这个项目还有后续我建议从三个方向延伸需求响应把冷热负荷从刚性改成弹性加入可中断负荷和价格型需求响应下层模型会变成更复杂的混合整数二次规划碳交易机制在目标函数里加入阶梯式碳交易成本关注二氧化碳排放约束对容量配置的影响多主体博弈把各微网从受系统统一调度变为独立利益主体用纳什均衡或主从博弈框架表征微网之间的策略互动双层结构会升级为多层或 MPEC 问题这三个方向其实是综合能源系统领域目前最热的课题如果你做的是硕士或博士课题完全可以在当前代码基础上继续迭代。就我个人经验来说这类双层优化配置项目最难的从来不是某一段代码怎么写而是你能不能把哪个变量放上层、哪个变量放下层、两层之间靠什么耦合想清楚。这个问题的答案直接决定了模型行为的合理性也决定了论文评审老师会不会追问。我一开始就是急着写代码结果被导师一句话问住你凭什么把储能容量放在上层如果运行方同时也是投资方这个双层结构还成立吗这个问号逼着我回去把系统物理背景又捋了三遍。最终落在纸面上的模型是所有思考过程的自然结果而不是代码本身。最后补一个小技巧不管用什么求解器强烈建议在代码里固定随机种子rng(2024)因为 PSO 是随机算法不固定种子的话每次运行结果都有波动写报告时数据对不上会很痛苦。我用这个办法保证了所有对比实验都在同一基准下进行也让后续的敏感性分析结果更可信。