★ FEATURED ARTICLE 简介MATLAB平台上的社区综合能源系统双层优化代码面向研究微网、综合能源、需求响应与动态定价的科研人员与研究生。程序以综合能源系统整体收益为上层目标将电价、热价等作为决策变量下层运营商与负荷聚合商作为跟随者构建Stackelberg主从博弈模型并计入电/热功率平衡约束上层采用DE优化算法下层调用CPLEX求解器实现上下层嵌套寻优解决多主体交互决策问题。压缩包共11个文件全部为m脚本大小约13KB包含主函数与多个子函数注释清晰运行即可输出全部图支持在原始数据基础上修改与扩展。目前已有1805人学习下载。代码创新性较强适合理解主从博弈、双层优化及综合能源动态定价的读者上手学习。1. 考虑需求响应微网与社区综合能源系统优化代码先搞清它优化的是什么拿到一套带注释、能跑通、还能扩展的微网/社区综合能源系统优化代码第一件事不是点运行而是先问一句这套代码到底在算什么在综合能源系统领域所谓“优化”主体永远是设备出力与能量平衡核心诉求是“在满足电、热、冷负荷的前提下把一天24小时的运行成本压到最低”。需求响应则是这个框架里的一层柔性手段——它把一部分用户负荷当作可调节资源让负荷曲线跟着电价或补偿信号走。这套代码的价值在于把需求响应、冷热电联供、储电储热储冷放进同一个混合整数线性规划MILP模型里用求解器在几十秒内给出全局最优的调度方案并把全部结果画成可用于EI论文的曲线。适合谁正在做微网、社区综合能源系统优化调度方向的研究生需要一套可复现、可扩展的基线代码来支撑对比实验的人。2. 模型先立住目标函数与约束的写法决定代码能不能扩展优化代码的骨架是数学模型。代码可以换语言、换求解器但只要目标函数和约束的表达方式不变结果就基本一致。这章节先讲清楚模型层面的三种核心选择——如何选目标、如何写能量平衡、如何处理设备非线性——再去对照代码看实现。2.1 目标函数运行成本最小化但成本项不是只有购电综合能源系统的运行成本最常见的表达式是min C C_grid C_gas C_om C_dr C_curtail其中C_grid 是向电网购电的费用购电单价随24小时分时电价变化C_gas 是燃气轮机和燃气锅炉消耗天然气折合成的费用按天然气热值和单价计算C_om 是各设备的运行维护成本一般按出力量乘以运维单价比如燃料电池每发1kWh电对应0.02元的维护费用C_dr 是需求响应补偿成本C_curtail 是弃风弃光的惩罚项设成很高的单价比如5元/kWh让求解器主动避免切出力。这套写法在EI期刊里很常见也是代码默认的目标。做扩展时你要注意一点大多数模型不会在目标函数里直接写“碳排放最小”而是把碳排放转成碳税或配额成本加进去。如果你的研究方向是双碳建议保留这个目标结构在C_gas里增加一个碳税系数即可不需要推翻重写。2.2 能量平衡约束电、热、冷三条母线的写法微网/社区综合能源系统的核心约束是能量平衡代码里通常会把它写成三条母线方程电平衡购电 光伏出力 燃气轮机发电 储能放电 需求响应削减量 电负荷 电锅炉耗电 电制冷机耗电 储能充电热平衡燃气轮机余热回收 燃气锅炉供热 储热罐放热 热负荷 吸收式制冷机耗热冷平衡电制冷机供冷 吸收式制冷机供冷 储冷罐放冷 冷负荷注意热平衡里的吸收式制冷机耗热这一项是冷热电联供系统特有的耦合约束——燃气轮机的排烟余热既能直接供热也能驱动溴化锂机组制冷。代码里会用一个系数把吸收式制冷机的出力折算成耗热量这个折算系数通常在1.3到1.5之间取决于COP制冷效能比。调试时如果热平衡总是对不上先检查这个系数是否写反。2.3 设备模型燃气轮机、储能与需求响应的线性化燃气轮机的发电效率和出力区间是非线性的代码常用的做法是把出力区间分成几段用分段线性化逼近效率曲线或者按工程惯例让效率在某个负荷率区间内近似恒定。另一种更常见的做法是给燃气轮机加一个最小出力约束比如“出力的连续变量取值落在[Pg_min, Pg_max]区间内”用半连续变量实现——这也是YALMIP或Gurobi里最容易处理的方式。储能部分的核心是SOC荷电状态递推方程SOC(t) SOC(t-1) eta_ch * P_ch(t) * dt / Cap - P_dis(t) * dt / (eta_dis * Cap)约束里还要限制SOC上下限一般是10%~90%、充放电功率上限以及“不能同时充放电”的逻辑。这个逻辑在MILP里就是两个二进制变量之和小于等于1。代码里会写成ch_on dis_on 1求解器才能正确处理。需求响应在这层模型里的本质是给负荷方程增加一个可变量。后面单独章节细讲。3. 需求响应建模价格型与激励型怎么把柔性负荷塞进MILP需求响应是这套代码最大的亮点也是新手最容易写烂的地方。社区综合能源系统里的需求响应一般分两种玩法价格型需求响应和激励型需求响应。代码里通常同时建模但实现思路完全不同。3.1 价格型需求响应可平移负荷的时段挪移价格型需求响应的核心思想是用户会自发地把洗衣、充电等可推迟负荷从高电价时段挪到低电价时段。代码里用一个可平移负荷矩阵L_shift(t, d)来表达——t是开始时段d是持续时长。举个例子一台持续2小时的洗衣机用户原本打算在18点开始用系统给它算了一笔账如果把开始时间挪到凌晨2点电费能省2.3元于是求解器把L_shift(18,2)置1L_shift(2,2)置118点对应的2小时负荷被平移走。代码里对这种负荷的关键约束有三个每个可平移负荷只能选择一个开始时段被平移后的负荷要保证连续供电不能拆成两半同一时刻开始的可平移负荷总数不能超过线路容量限制。价格型需求响应不需要在目标函数里加补偿成本因为用户省下的电费就是激励本身。3.2 激励型需求响应可削减负荷的0-1变量激励型需求响应则完全不同用户签了协议允许运营方在高峰时段切掉一部分负荷比如空调温度上调2度、工业设备降载作为交换运营方按削减量给用户补偿。代码里用Load_shed(t)表示每个时段的削减量并定义一个0-1变量u_shed(t)u_shed(t)1表示该时段执行了削减Load_shed(t)的范围是[0, Load_max_shed]一天内累计削减次数有限制避免把用户折腾得太过分。补偿成本在目标函数里单独成项C_dr sum( Load_shed(t) * price_dr(t) )price_dr(t)是单位削减补偿单价通常比正常电价高出20%~50%具体数值要看需求响应项目的合同定法。3.3 需求响应适度参与为什么负荷平移要设上限代码里如果允许负荷完全自由平移求解器会给出一个极端答案把白天所有负荷都挪到凌晨电费最小了但用户舒适度归零了。这在工程上完全不可行。所以代码里定了两个硬约束平移负荷占总负荷的比例不超过15%~20%具体比例可以在参数区修改价格型需求响应部分每个时段的最大平移量受限于该时段可平移负荷总量和线路容量。类似地激励型需求响应也有“每日最大削减时长”“连续削减间隔”等约束。这些约束不是论文里堆字数的装饰它们是决定模型能不能用于实际项目、能不能过审稿人那关的关键细节。4. 代码落地从参数区到绘图脚本逐段拆解运行流程模型搞清楚之后再看代码就轻松了。以MATLAB YALMIP Gurobi/Cplex这套最常用的技术栈为例把这套代码从头到尾的典型文件结构和核心片段拆开讲。4.1 输入参数区把24小时负荷、电价、天然气价格放在一张表里代码的第一步永远是载入基础数据。社区综合能源系统的典型做法是直接在脚本里定义一个结构体或用Excel读入24小时数据。从代码维护角度看我建议把负荷、电价、天然气价格、光照强度这四组数据单独存成CSV或Excel主程序用readtable读入。% 读取24小时基础数据 data readtable(load_price.xlsx); P_load data.P_load; % 电负荷单位kW24x1 H_load data.H_load; % 热负荷单位kW24x1 C_load data.C_load; % 冷负荷单位kW24x1 price_e data.price_e; % 分时电价单位元/kWh, 24x1 price_g 0.35; % 天然气单价元/kWh按热值折算这里的price_e是分时电价代码里常见的是峰谷平三段费率比如高峰1.2元、平段0.75元、低谷0.4元。注意天然气价格折算到元/kWh时要用天然气热值除以锅炉效率。比如天然气低位热值9.7元/Nm³折合每kWh约0.35元不同地区差异很大代码里这个参数要按当地实际气价改。4.2 决策变量定义sdpvar与binvar怎么混用YALMIP里连续变量用sdpvar定义0-1整数变量用binvar定义。MILP模型里两类变量同时存在求解器才能正确处理逻辑约束。% 决策变量定义 Pg sdpvar(24,1); % 燃气轮机发电出力 Hg sdpvar(24,1); % 燃气锅炉供热出力 Pchp sdpvar(24,1); % 热电联产余热回收功率 Pgrid sdpvar(24,1); % 从电网购电功率 SOC sdpvar(24,1); % 储能荷电状态 Pch sdpvar(24,1); % 储能充电功率 Pdis sdpvar(24,1); % 储能放电功率 % 需求响应变量 load_shift binvar(24,1); % 每个时段是否执行可平移负荷 load_shed sdpvar(24,1); % 每个时段的负荷削减量 u_shed binvar(24,1); % 负荷削减状态标志 % 储能同时充放电标志 u_ch binvar(24,1); u_dis binvar(24,1);binvar数量会直接影响求解速度。一个典型社区能源系统只有24个时段、几条母线约束binvar通常在200个以内Gurobi在几秒内就能求到MIP gap小于0.1%的最优解。如果你的代码里binvar超过500个可以先把24小时改成8个时段来验证模型正确性再逐步加密。4.3 约束装配把能量平衡与设备模型写成YALMIP约束数组约束装配是代码里最长、也最容易出错的部分。每个设备的约束独立写最后拼成一个约束数组F交给optimize求解。% 电平衡约束含需求响应削减量 F [F, Pgrid Pg Pdis load_shed P_load - load_shift P_eb P_ec]; % 热平衡约束 F [F, Pchp Hg H_dis H_load H_ac]; % 冷平衡约束 F [F, C_ec C_ac C_dis C_load]; % 燃气轮机约束出力上下限与爬坡约束 F [F, 0 Pg Pg_max]; F [F, -ramp_limit Pg(2:24) - Pg(1:23) ramp_limit]; % 储能SOC递推与充放电互斥 F [F, SOC [SOC_0; SOC(1:23)] eta_ch * Pch / Cap - Pdis / (eta_dis * Cap)]; F [F, 0 SOC 0.95]; F [F, u_ch u_dis 1]; F [F, Pch u_ch * Pch_max, Pdis u_dis * Pdis_max];储能SOC这行用了一个技巧SOC(1:23)表示将SOC的前23个元素作为t-1时刻的值向后推移在YALMIP里这样可以直接写出一整条时序递推约束不用for循环代码更紧凑。很多人第一次接触这段会看不懂手动展开成24行递推式会更直观。4.4 求解与结果校验不要直接信最优解ops sdpsettings(solver, gurobi, showprogress, 1, verbose, 2); optimize(F, objective, ops); % 校验求解状态 if sol.problem 0 Pg_opt value(Pg); Pgrid_opt value(Pgrid); % 计算各项成本 cost_grid sum(Pgrid_opt .* price_e); cost_gas sum(Pg_opt ./ eta_g) * price_g; else error(Solver failed with code %d, sol.problem); end求解完成之后代码里的统计输出会把购电成本、天然气成本、需求响应补偿成本、储能收益拆成清单打印出来。这个成本对账表非常重要——论文里的经济性分析数据全是从这里来的。建议每个实验跑完都先看这张表对不上账就说明某个约束写错了。绘图部分截图代码figure(Position, [100 100 1200 700]); subplot(3,2,1); bar(1:24, [Pg_opt, Pgrid_opt, Pdis_opt, Pch_opt], stacked); legend(燃气轮机, 购电, 储能放电, 储能充电); xlabel(时刻(h)); ylabel(功率(kW)); title(电功率平衡);这一组图最后会拼成论文里的“系统优化调度结果图”要保证图例顺序、坐标单位、字体大小都一致投稿时才不会被打回。5. 避坑与常见问题为什么代码跑不出论文里的图这套代码能跑通和能跑出正确结果是两回事。实际使用中出现的问题大多是模型层和参数层的不是代码层的。以下是五个高频问题按“现象→原因→解决”写照方抓药即可。5.1 现象求解器提示infeasible但模型看着没问题这是MILP模型最常见的问题。原因大概率是某个约束相互矛盾比如储能的初始SOC、最小SOC和最大充放电功率之间在24小时的窗口内根本不可能满足。另一种典型原因是SOC递推约束里使用了SOC_00.5但储能容量Cap设得过小导致某个时段必须同时满足充电需求和放电需求。解决方法是把SOC_0改成0.5把Cap的值加大或者把SOC初始值从固定值改成约束SOC(1)SOC_0且SOC(24)SOC_0给求解器留余地。也可以用YALMIP的assign和check命令逐一检查每条约束的残差先定位在哪一行。5.2 现象热平衡“差一口气”吸收式制冷机的耗热老是溢界现象是优化结果里热母线的总有功供应与热负荷总是不平差的那部分刚好等于吸收式制冷机的耗热量。原因很常见热平衡约束里忘了把“吸收式制冷机消耗的热量H_ac”从热负荷侧扣掉或者H_ac本身是用COP折算的但折算方向写反了。解决方法是回到公式吸收式制冷机供出C_ac的冷量需要消耗H_ac C_ac / COP_abs的热量这个量要同时出现在热平衡方程右侧作为耗热项和冷平衡方程左侧作为供冷项。这两行如果分开写在代码的不同位置漏写一个就出这种“幽灵差量”。5.3 现象SOC曲线乱跳或者充放电功率同时非零SOC曲线出现阶梯状跳变通常是部分SOC值越过上下界导致的。如果SOC曲线一天之内反复冲到0.95又跌回0.1说明目标函数里缺少对储能“寿命损耗”的考量求解器让电池满充满放。工程上可以给SOC增加变化幅度惩罚或者在目标函数里加一个很小的单位容量折旧成本。如果充放电功率同时非零则是充放电互斥约束没生效检查u_ch u_dis 1是否真的进入了约束数组F。一个容易犯的错是在YALMIP里约束条件用分号结尾直接就丢掉不装配了所以这种逻辑约束要单独一行并用逗号或赋值号连接。5.4 现象换一台电脑脚本报“未定义函数或变量”这通常是路径问题。代码里用到了自定义函数比如读取数据、画图但主程序没有把自定义函数所在的文件夹加入MATLAB路径。解决方法是把整个项目目录一次性加进路径或者把主程序改成用绝对路径引用数据文件。如果代码里用了gurobi和yalmip还要确认这两个工具箱在新机器上已安装并配置好。建议多看求解器的verbose输出——如果Gurobi的license有问题YALMIP往往会在求解前直接报错一查便知。5.5 现象画出图来图例漂移、中文乱码、导出PDF字号太小中文乱码是老问题。MATLAB在Windows下用SimHei或Microsoft YaHei字体能显示但exportgraphics导出时经常把中文变成方框。我的常用做法是所有图例和坐标轴标签用英文只在标题里放中文论文中需要中文说明的在PPT里后补。字号方面MATLAB默认的FontSize是10期刊图一般要求最小字号不小于6号字投稿前用exportgraphics(fig, result.pdf, ContentType, vector)导出矢量图确保所有线条和文字清晰。6. 在基础上扩展从单日示范到多场景对比的验证技巧拿到这套代码后不要急着改模型先做一次“论文复现验证”把代码默认算例的结果和你目标论文里的数据对上。对不上的话先把目标函数的成本单价、设备容量参数逐项对齐再往下走。三个值得优先做的扩展方向第一改成多日连续优化。单日优化的边界问题是储能初始SOC需要猜测多日优化则让储能连续跨天运行更贴近工程实态。修改方法是把维度从24扩展为24*N天循环约束里注意日与日之间的SOC衔接。第二加入多场景对比。把光伏出力从单一日照曲线改成晴、多云、阴雨三组数据分别优化并比较成本和需求响应效果。这种对比是EI论文最吃香的图——一张三维堆叠图就能直观展示不同光照条件下的设备出力与成本差异。第三把确定性优化换成鲁棒优化。在现有约束里加入光伏出力的不确定集目标函数改成“最坏情况下成本最小化”。这一步改动相对独立——在YALMIP里本质上是用bilinear项替换原线性项但要注意Gurobi对bilinear的处理较慢建议优先用Cplex或改用线性化公式。我这几年带学生做这类复现遇到最多的问题不是模型写不出来而是“能画图但不敢确认结果是对的”。对付这个问题的土办法就一个手算一个最简单算例——把设备数量缩减到1台、时段缩减到4个把结果和代码输出对一遍确认逻辑通了再放回去跑24小时全模型。这个习惯救过我很多次也希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
