简介:这份基于双层优化的微电网系统规划设计MATLAB程序,面向电气工程与能源系统方向的研究生、科研人员及工程师,重点解决微电网中电源容量配置、运行调度策略协同优化等问题。压缩包共6个文件,以MATLAB脚本(.m)为主,并包含一个Excel典型日数据文件(.xlsx),资源整体仅38KB,涵盖上层规划模型、下层运行优化、参数设置模块及主程序入口,模块边界清晰,便于按需拆解与复用。目前已有165人学习下载,是理解双层优化建模思路与算法实现的紧凑型参考。借助该程序,可快速复现微电网规划设计的完整流程,掌握上下层模型交互逻辑、约束条件处理以及典型日数据的组织方式;同时可以在现有代码基础上调整参数、替换数据场景,开展多情形算例验证,为后续改进优化模型、扩展算法功能提供了可运行的代码基础。
1. 微电网规划为什么要用双层优化,而不是一次优化算到底?
看“matlab程序-基于双层优化的微电网系统规划设计方法”这个标题,核心矛盾在于:规划者要做的是二十年长期投资决策,运行者却必须设法应对每一个小时的功率平衡。把这两个时间尺度硬塞进一个单层优化模型,要么用典型日比例系数糊弄运行过程,要么把调度策略简化成基载加峰谷叠加,算完的容量不是偏大就是偏小。双层优化的思路是把设备选型与容量配置放上层,把多场景逐时调度放下层,上层算年化总成本,下层按典型日回传运行费用,两层在迭代中收敛到同一组容量解。这套方法在配电网扩容、独立微网和并网型微网规划中都是主流做法,适合做选址定容、储能容量测算,或者写论文需要方法对比的工程师与研究生。
2. 双层优化在微电网里算什么——上层投资决策与下层运行调度的变量交接
2.1 单层模型在微电网场景里为什么会失真
常见的单层规划做法是“年费用加约束一并优化”,将运行约束用最大负荷、最小负荷、年利用小时数等统计量近似。这样算出来的储能容量往往偏大,原因在于储能收益来自电价差和削峰填谷,而在单层模型里全天8760小时的时序信息被压成几个峰谷比,储能的调度价值被严重低估。反过来,如果只按最恶劣场景规划,设备又会过度配置。微电网规划要同时回答“装多少”和“怎么用”,两个问题时间尺度相差巨大——一个以年计,一个以小时甚至分钟计,因此工程实践和文献里都会把它构造成双层优化。
2.2 上层模型:全寿命周期投资成本最小化
上层规划层的目标函数一般写作:
min C_inv + C_OM + C_rep + C_buy + C_fuel
其中 C_inv 为光伏、风电、储能、PCS 的初始投资折现值,C_OM 为年运行维护费用折现,C_rep 为电池寿命期内的更换成本,C_buy 为从电网购电费用,C_fuel 为柴油发电机燃料费用。决策变量是设备安装容量,即 x = [Ppv_cap, Pwt_cap, Ebatt_cap, Pbatt_pcs]。约束包括可用安装面积、最大投资预算、单类设备容量上下限等。
这里的关键在于 C_buy 和 C_fuel 无法由上层直接求出,它们取决于逐时运行策略。因此上层目标函数必须嵌套一个下层优化回传的“运行总成本”,这就构成了两层之间的第一次变量交接。
2.3 下层模型:多场景逐时运行成本最小化
下层运行调度层按季节或聚类后的典型日展开,每个典型日取24个时段。给定设备容量后,对每个典型日求解:
min sum(购电成本 + 柴油燃料成本 + 切负荷惩罚)
约束至少包括:功率平衡约束、储能荷电状态递推公式、充放电功率上下限、爬坡约束、联络线功率上限、备用容量约束。SOC 递推的简化形式为:
SOC(t+1) = SOC(t) + eta_ch * Pch * dt / E - Pdch * dt / (eta_dch * E)
下层的解输出一组逐时调度文件,代回上层计算真实运行成本,再让上层调整容量。两层之间的交接变量可归纳如下:
| 变量类别 | 上层(规划层) | 下层(运行层) | 传递方向 |
|---|---|---|---|
| 设备容量 | 决策变量 | 固定参数 | 上→下 |
| 逐时出力 | —— | 决策变量 | —— |
| 运行成本 | 目标函数组成项 | 目标函数值 | 下→上 |
| SOC 状态 | —— | 决策变量 | —— |
2.4 双层优化的常见数学形式与解法分类
| 上下层形式 | 典型特征 | 常用解法 |
|---|---|---|
| 混合整数线性双层 | 两层均为 MILP | KKT 条件加大 M 法转单层,再用 Gurobi/CPLEX |
| 上层 MILP、下层 LP | 下层线性可微 | 强对偶或 KKT 转单层 |
| 非凸双层 | 下层含整数变量或非线性 | 智能算法外层,内层求解器 |
| 迭代枚举法 | 容量变量离散化 | 遍历容量组合,内层逐次求解 |
在 MATLAB 里做规划,我会先判断下层的连续性和凸性。若下层是纯线性规划,KKT 条件替换是效率最高的路径,配合 YALMIP 可以省掉手写对偶矩阵的推导;若下层包含 0/1 调度变量,则 4.2 节的智能算法方案更现实。
3. MATLAB 里的实现骨架——双层模型的代码落点与参数传递
3.1 数据准备:把全年 8760 小时聚合成典型日场景
写优化模型之前,先把负荷、风速、光照的全年序列读进来,按工作与非工作日加季节特征进行 k-means 聚类,挑出 12 个典型日。每个典型日必须带权重系数,即该类场景在全年的天数占比,这样下层结果按权重加权后才是真实的全年运行成本。
% 典型日聚类与权重计算 load hourly_data.mat; % 含 load_pu, wind_pu, pv_pu, price 四列 rng(42); % 固定随机种子,保证结果可复现 [idx, C] = kmeans([load_pu, wind_pu, pv_pu, price], 12, 'MaxIter', 500); weights = accumarray(idx, ones(size(idx)), [12, 1])' / length(idx); % C 的第 k 行对应第 k 个典型日:负荷、风、光、电价 % idx 是每天所属的簇编号kmeans 的第一个参数是特征矩阵,第三组参数 MaxIter 控制最大迭代次数,一般取 300 到 500 够用。rng(42) 固定随机种子,确保每次聚类结果一致,否则权重数组在两次运行间会变化,下层成本不可比。C 矩阵用于从全年数据中抽取典型日曲线,也可以直接用簇中心作为调度场景输入。
3.2 上层规划的 YALMIP 建模骨架
上层模型适合用 YALMIP 的声明式写法,因为模型在迭代中经常要增加面积约束、预算上限或可靠性约束,声明式改动成本比手工拼矩阵低得多。也可以用 MATLAB 自带优化工具箱的 linprog 或 intlinprog,但可读性差一些。
% 上层决策变量定义 Ppv_cap = sdpvar(1,1); % 光伏安装容量 kW Pwt_cap = sdpvar(1,1); % 风电安装容量 kW Ebatt_cap = sdpvar(1,1); % 储能额定容量 kWh Pbatt_pcs = sdpvar(1,1); % 储能 PCS 功率 kW unit_cost = [6500, 9000, 1800, 1200]; % 元/kW 或 元/kWh inv_cost = unit_cost * [Ppv_cap; Pwt_cap; Ebatt_cap; Pbatt_pcs]; Constraints = [Ppv_cap >= 0, Pwt_cap >= 0, ... Ebatt_cap >= 0, Pbatt_pcs >= 0]; % 可根据项目再补充面积上限、预算上限等约束unit_cost 行向量与容量列向量做内积得到初始投资。注意储能能量与功率是两个独立变量,对应单价不同,常见错误是只建一个电池变量,导致 PCS 容量与电池容量比例失衡。下层运行成本的接口通过下一小节的函数回传,上层求解器不能直接计算购电与燃料费用。
3.3 下层运行调度的函数化封装
下层模型我习惯独立写成函数,输入设备容量,输出加权后的全年运行成本。这样上层无论用遗传算法、粒子群还是枚举法,都只需要调用同一个入口,不会重复建模。
function total_cost = lower_level(Ppv_cap, Pwt_cap, Ebatt_cap, Pbatt_pcs, scen, weights) % 输入: 设备容量, 场景结构体, 典型日权重 % 输出: 12 个典型日按权重加权后的年运行成本 T = 24; total_cost = 0; for k = 1:length(weights) Pch = sdpvar(T,1); % 储能充电功率 Pdch = sdpvar(T,1); % 储能放电功率 SOC = sdpvar(T+1,1); % SOC 序列 Pgrid = sdpvar(T,1); % 向电网购电功率 Pdie = sdpvar(T,1); % 柴油机出力 PV = scen.pv_pu(:,k) .* Ppv_cap; % 光伏按典型日系数折算 WT = scen.wt_pu(:,k) .* Pwt_cap; % 风电按典型日系数折算 Pload = scen.load(:,k); % 功率平衡: 光伏+风电+柴油+放电+购电 = 负荷+充电 Cons = [PV + WT + Pdie + Pdch + Pgrid == Pload + Pch, SOC(2:end) == SOC(1:end-1) + ... eta_ch*Pch/Ebatt_cap - Pdch/(eta_dch*Ebatt_cap), SOC(1) == SOC(end), % 末端回初值 SOC >= 0.2*Ebatt_cap, SOC <= Ebatt_cap, Pch >= 0, Pch <= Pbatt_pcs, Pdch >= 0, Pdch <= Pbatt_pcs, Pgrid >= 0, Pgrid <= Pgrid_max]; Obj = sum(scen.price(:,k).*Pgrid + cost_diesel*Pdie); optimize(Cons, Obj, sdpsettings('solver','cplex')); total_cost = total_cost + weights(k) * value(Obj); end end这段代码的关键有两处。一是 SOC 维度设为 T+1,通过 SOC(1) == SOC(end) 实现运行周期内储能电量闭环,避免储能“只放不充”的白嫖行为;二是 eta_ch 与 eta_dch 分别表示充放电效率,实际项目中充电效率一般在 0.9 到 0.95,放电效率略低,二者不能合并。若想严谨区分充放电互斥,还需要引入二元变量或 Pch.Pdch == 2... 的非线性约束,示例中省略该行以便阅读,工程上建议用二元变量加大 M 法处理。
3.4 上下层之间的迭代接口
采用智能算法做外层时,上下联耦合可以收敛成单个适应度函数。上层每轮生成一组容量解,传给 lower_level,回传运行成本后与投资成本合并,得到年化总成本:
function obj = fitness(x) Ppv_cap = x(1); Pwt_cap = x(2); Ebatt_cap = x(3); Pbatt_pcs = x(4); run_cost = lower_level(Ppv_cap, Pwt_cap, Ebatt_cap, Pbatt_pcs, scen, weights); inv_cost = unit_cost * x'; % 示意写法: 实际需将投资按资本回收系数折算到年 obj = inv_cost + run_cost; end这里最大的坑是时间尺度不统一。inv_cost 是建设期一次投资,run_cost 是运营期一年费用,两者直接相加会高估投资占比。常见做法是用资本回收因子 CRF 把投资折算成等年值,即 obj = CRF * inv_cost + run_cost,CRF = r*(1+r)^n / ((1+r)^n - 1),r 为折现率,n 为项目寿命。
4. 求解实践与参数配置——KKT 替换、智能算法与收敛判据
4.1 KKT 条件把双层转成单层的 MATLAB 路径
当上层连续或混合整数线性、下层是 LP 时,标准做法是对下层拉格朗日函数写 KKT 条件,将互补松弛条件用大 M 法线性化,整成单层 MILP。YALMIP 可以逐层建模后再合并,但更快的路径是直接借助求解器的二次约束支持,把下层用强对偶替换。互补松弛条件无法直接用 optimize 表达,常见处理是引入二元变量:
% 互补松弛条件线性化示例 % dual_var 对应某条不等式约束的 KKT 乘子 % slack 是原约束的松弛量 M = 1e5; % 需比最大可行松弛量放大一个量级 z = binvar(size(dual_var)); Constraints = Constraints + ... [dual_var >= 0, slack >= 0, ... dual_var <= M*z, slack <= M*(1-z)];M 的取值决定求解成功与否。取得过大会让整数规划数值病态,取得过小可能剪掉最优解。工程上我会先跑一次去掉互补松弛条件的松弛模型,统计各约束松弛量的最大量级,再取该量级的 50 到 100 倍作为 M,这样既保证数值稳定又不至于剪枝。
4.2 用粒子群或遗传算法处理非凸双层优化
如果下层带 0/1 变量或储能寿命衰减非线性,双层问题不再等价于单层 MILP。此时常见的做法是把上层容量决策交给智能算法,下层照常用求解器求最优。MATLAB 的 Global Optimization Toolbox(国内资料里常叫 matlab 优化工具箱)自带 particleswarm 和 ga,中等规模的微电网规划够用:
lb = [0, 0, 0, 0]; % 容量下限 ub = [500, 300, 1000, 500]; % 容量上限,单位同上 options = optimoptions('particleswarm', ... 'SwarmSize', 40, ... 'MaxIterations', 60, ... 'FunctionTolerance', 1e-4, ... 'Display', 'iter'); [x_opt, fval] = particleswarm(@fitness, 4, lb, ub, options);SwarmSize 取 40 是微电网场景的经验值,典型日数量不多时单轮适应度计算只需几秒,40 个粒子可以在合理时间内覆盖容量搜索空间。FunctionTolerance 表示种群最优改进小于 1e-4 时提前停机,这个值设太大容易卡局部最优,建议先放松到 1e-6 观察收敛曲线形态。
注意:particleswarm 对适应度函数的连续性没有要求,因此在处理带 0/1 变量的下层模型时比 KKT 方法更宽容。但它的随机性意味着同一组参数每次运行结果略有不同,正式投产或发表结果前务必固定随机种子并在多组起始点上重复验证。
4.3 求解器选型与三个必调参数
| 场景特征 | 推荐求解器 | 三个关键参数 |
|---|---|---|
| 线性下层,规模小于 10 万变量 | Gurobi / CPLEX | MipGap=0.01, TimeLimit=3600, Threads=8 |
| 线性下层但无商业求解器 | intlinprog | MaxTime, IntegerTolerance, LPPreprocess |
| 非线性上层,反复迭代 | particleswarm + 内层 Gurobi | SwarmSize, StallIterLimit, FunctionTolerance |
| 储能非线性寿命模型 | ga + fmincon 混合 | PopulationSize, MaxGenerations, FunctionTolerance |
有一个容易踩的坑:分时电价曲线变化很陡时,默认 MipGap 会让求解器过早停在一个次优解上。规划项目里建议把 MipGap 调到 0.001 甚至更低,运行时间从几分钟涨到几十分钟,但方案总成本往往能再压下来 3% 到 5%,对于二十年项目寿命来说这笔计算成本完全值得。
5. 验证技巧——如何判断两层优化结果是可信方案
5.1 用单层模型做下界对照
写完双层模型先别急着收工。我要做的第一件事是用同一数据集跑一遍单层等效模型:把典型日按最大负荷比例叠加,忽略逐时调度约束,只算总容量与总费用。双层结果的总成本应当高于单层结果,因为前者更贴近真实运行约束。如果高了超过 15% 到 20%,大概率是下层运行策略没给足自由度,比如储能被爬坡约束或联络线限制锁死。把双层与单层的容量差列成对比表,哪个设备容量差异最大,就去检查对应约束的边界值是否过于保守。
5.2 场景数量与权重敏感性验证
典型日数量从 4 个增加到 12 个,规划容量变化通常应在 5% 以内。如果容量跳变超过 10%,说明聚类场景没有覆盖极端光伏日或极端负荷日。这时可以手动把“夏季工作大负荷加夜间无光”和“冬季连续阴天”两个极端日作为强制场景加入,再在 5% 范围内扰动全年权重,观察目标函数变化率。实现上只需在 scen 结构体追加两组手工场景,并把 weights 对应位置手动赋值,不需要改动模型主体。
5.3 收敛性观察与一个实用技巧
用智能算法迭代时,把每代最优值画成收敛曲线,正常形态是前十代快速下降后趋于平缓。如果曲线呈锯齿状且后期大范围跳跃,说明上层容量与下层成本之间的映射存在剧烈非凸,常见于电池容量恰好越过某条功率约束的场景。我常用的加速技巧是前三轮迭代只用四个典型日粗筛容量可行域,确认搜索范围后再切换为 12 个典型日精算,同一套代码只改 weights 和 scen 的数据来源,能省掉 30% 到 50% 的总计算时间。
最后提一个直接可用的小动作:优化完成后把典型日内储能 SOC 曲线画出来,检查是否在日前电价谷段充电、峰段放电。如果 SOC 曲线几乎不变,说明下层目标函数权重配比有问题,或者 SOC 末端闭合条件写反了,这两行 plot 加在 lower_level 函数返回值之后就能完成可视化校核。
本文还有配套的精品资源,点击获取