做微电网优化的人应该都有同感多能互补热电联供型微网优化在论文里看起来非常顺理成章无非是建立目标函数、列约束、求最优解但真正动手把它变成一套能跑出曲线的MATLAB代码事情就完全不一样了。变量维度、单位、求解器配置、约束写法每一步都可能让程序报错更别提还要保证算出来的结果跟文献趋势一致。这套代码是我照着相关文献思路重新整理出来的MATLAB完整实现注释很细从数据初始化到最终画图都有说明。不管是刚接触微网调度优化的学生还是想验证模型思路的工程师都可以直接用这套代码做底子。1. 我为什么花两周时间把多能互补热电联供微网优化模型写成MATLAB代码先说点背景。多能互补热电联供型微网听起来复杂核心其实就是把电、热两类能量放在一个系统里统一调度让不同能源设备各司其职在满足负荷的前提下把运行成本压到最低。传统微网调度通常只盯电功率平衡加入热负荷和热储能之后问题规模翻倍变量之间的耦合关系也多了代码实现难度明显上升。我最初接触这类优化模型时最头疼的不是数学模型本身而是怎么把论文里那些数学符号准确变成可执行的代码。论文里写一个P_t^{CHP}就完事了但代码里你得定义它是一个1×24的变量向量还是24×1的列向量你要算一天的运行成本就要考虑电价曲线怎么读、单位怎么统一、热负荷和电负荷的数据从哪来。这些细节论文不会告诉你只有自己动手写一遍才能真正理解。这套代码解决的问题很直接给出一套可复现的框架系统内包含光伏、风电、燃气轮机CHP机组、燃气锅炉、电储能、热储能目标函数是运行成本最小化约束涵盖能量平衡、设备出力上下限、爬坡约束、储能SOC约束以及电网交互约束。通过MATLAB建模并调用求解器一次性算出全天24小时各设备的出力计划。如果你只需要看结果曲线代码跑完就能出图如果你想深入改造模型代码里每个约束都标记了来源和物理含义可以直接替换设备参数或增加新约束。这也是我把它称作“完美复现”的原因——不是指结果跟某篇文献一模一样而是指建模逻辑完整、注释覆盖到位拿到手不靠猜。2. 模型侧先搞清楚热电联供微网优化到底在优化什么写代码之前必须先把优化模型本身梳理清楚。这一部分如果模糊后面所有代码都站不住脚。2.1 系统拓扑与能量流这个微网系统里涉及的主要设备按能量类型可以分成三类能量类型设备作用电源光伏、风电可再生出力优先消纳通常不作为优化决策电/热耦合燃气轮机CHP机组同时产生电功率和热功率是多能互补的关键热源燃气锅炉、热储能补充供热平衡热负荷与CHP产热之间的差额储能电储能、热储能削峰填谷平滑新能源波动外部交互配电网允许从电网购电电价按分时电价设定能量流的核心逻辑是电负荷由光伏、风电、CHP发电、电储能放电以及电网购电共同满足热负荷由CHP余热、燃气锅炉供热、热储能放热共同满足。注意CHP机组在这里处于“电热耦合”位置发电多的同时产热也多所以调度时必须同时权衡电平衡和热平衡不能单独看某一边。2.2 目标函数最小化全天运行成本目标函数是优化问题的“指挥棒”它决定了系统会往哪个方向搜索最优解。这套代码选择最小化一天内的总运行成本公式分三块第一外购电成本也就是从配电网买电的费用等于各时段购电功率乘以对应时段的电价再累加。分时电价结构直接影响储能充放电策略和CHP机组的运行方式。第二燃料成本包括燃气轮机CHP机组消耗的天然气和燃气锅炉消耗的天然气。燃料成本按热值折算注意CHP机组同时输出电和热所以它的发电效率和热电比都要写进计算。第三如果模型考虑设备启停还可以加入启停成本项但我这个版本里没有加因为很多文献的基础模型也不考虑启停简化处理更利于初学者上手。目标函数写成数学形式就是[ \min \sum_{t1}^{24} \left[ \text{grid_price}(t) \times P_{grid}(t) C_{gas} \times \left( \frac{P_{chp}(t)}{\eta_{chp,e}} \frac{H_{gb}(t)}{\eta_{gb}} \right) \right] ]这里grid_price(t)是分时电价C_gas是天然气单价折算系数P_chp(t)、H_gb(t)分别是CHP电功率和燃气锅炉热功率η_chp,e和η_gb是设备效率。目标函数每项都要统一单位我用的是kW和元/kWh读者拿到别的项目时第一件事也应该是确认单位。2.3 约束条件把物理极限翻译成数学不等式没有约束的优化问题没有意义。这套代码里的约束分四层第一层是电功率平衡。每一时刻光伏出力、风电出力、CHP电出力、电储能放电功率、电网购电功率之和必须等于电负荷加上电储能充电功率。这里有个新手容易漏掉的点储能充电是负荷不是电源所以它出现在等号右侧。第二层是热功率平衡。CHP余热回收功率、燃气锅炉热出力、热储能放热功率之和等于热负荷加上热储能充热功率。第三层是设备运行约束。每台设备有出力上下限CHP机组还要考虑爬坡约束储能需要同时设置充放电功率上限和SOC上下限。第四层是电网交互约束。与配电网的购电功率不能超过联络线容量如果需要还可以限制不能向电网反送电根据所复现的文献要求来定。这些约束在代码里对应的是Constraints [Constraints, ...]这样的YALMIP约束累加写法。后续我会把这部分拆开解释。3. MATLAB代码结构让“注释详细”不只是口号很多复现代码的毛病是脚本一锅烩所有数据、变量、约束、求解全部堆在一个文件里看的人很难下手。这套代码在结构上做了拆分逻辑清晰每部分职责单一注释详细到每个变量都标注了物理含义和单位。3.1 文件组织方式代码包含以下几个文件和脚本各司其职文件职责main.m主脚本负责初始化、建模、求解、调用绘图data_input.m定义所有设备参数、负荷曲线、电价曲线build_model.m定义决策变量、目标函数和约束条件solve_model.m配置求解器并求解输出结果结构体plot_results.m绘制各设备出力曲线、储能SOC曲线、成本对比图实际运行时main.m依次调用后面几个模块整个过程大概几秒钟就能跑完。把代码拆成模块的好处是你想改参数就在data_input.m里改想换模型就在build_model.m里改动量或约束不需要在几百行代码里反复搜索。3.2 参数定义与数据准备data_input.m里定义了所有基础数据。以24小时为调度周期时间分辨率取1小时因为绝大多数热电联供类文献都采用这个粒度既能反映日出荷特性和分时电价变化又不至于因为时间尺度过细导致混合整数规划求解过慢。负荷数据方面电负荷和热负荷曲线是直接从典型日数据里提取的画出来是典型的双峰形状电负荷早晚各有一个峰值热负荷在夜间偏高。如果不方便抄文献数据也可以用MATLAB自带的随机函数生成平滑曲线但要注意加个随机种子保证每次运行结果可复现。设备参数方面CHP机组的额定发电功率、发电效率、热电比燃气锅炉的额定热功率和效率电储能和热储能的容量、最大充放电功率、初始SOC、SOC上下限都要在这里一次性定义。单位统一用kW和kWh电价用元/kWh天然气费用折算到元/kWh。3.3 求解器与建模工具的选择代码采用YALMIP作为建模层用CPLEX求解底层优化问题。为什么这么选YALMIP的变量定义方式和约束累加方式非常接近数学描述写出来的代码可读性极高适合复现论文模型。CPLEX作为商业求解器处理线性规划和混合整数线性规划的速度很快特别是当模型加入储能二进制状态变量后整数变量一多普通MATLAB内置求解器跑起来会很吃力。如果你的环境里没有CPLEX也可以把求解器换成Gurobi或MATLAB自带的linprog。除非模型加入了储能状态二进制变量必须用intlinprog否则linprog也能运行。我在代码里加了一段求解器自动检测逻辑如果检测到CPLEX就用CPLEX检测不到就回退到optimize默认求解器这样在没装CPLEX的电脑上也不会一上来就报错。4. 关键代码片段逐段拆解这一节是整套代码的核心我挑几个最关键的部分逐段说明你会看到数学公式是怎么变成MATLAB代码的。4.1 决策变量定义在YALMIP里变量用sdpvar定义连续变量用binvar定义二进制变量。这里定义如下% 决策变量定义24小时 P_pv sdpvar(1, 24); % 光伏出力单位kW实际上通常作为已知参数 P_wt sdpvar(1, 24); % 风电出力单位kW P_chp sdpvar(1, 24); % CHP机组发电功率单位kW H_chp sdpvar(1, 24); % CHP机组余热回收功率单位kW H_gb sdpvar(1, 24); % 燃气锅炉热出力单位kW P_cha sdpvar(1, 24); % 电储能充电功率单位kW P_dis sdpvar(1, 24); % 电储能放电功率单位kW H_cha sdpvar(1, 24); % 热储能充热功率单位kW H_dis sdpvar(1, 24); % 热储能放热功率单位kW P_grid sdpvar(1, 24); % 电网购电功率单位kW SOC_e sdpvar(1, 24); % 电储能荷电状态单位% SOC_h sdpvar(1, 24); % 热储能荷电状态单位%这里说明一下光伏和风电在大多数确定性调度模型里作为已知的预测出力输入不参与优化决策所以代码里虽然把它们写成了sdpvar实际使用时应该直接赋值为预测曲线数据。这个写法是为了让能量平衡约束写得统一不需要区分变量和参数。4.2 能量平衡约束能量平衡是微网优化里最不能出错的约束。电功率平衡约束如下% 电功率平衡发电 负荷 储能充电 Constraints [Constraints, ... P_pv P_wt P_chp P_dis P_grid P_load P_cha];热功率平衡约束类似% 热功率平衡产热 热负荷 热储能充热 Constraints [Constraints, ... H_chp H_gb H_dis H_load H_cha];这两行代码是YALMIP的精髓所在。在YALMIP中不是执行逻辑判断而是生成一个约束对象这个对象被添加进约束数组后求解时由求解器统一处理。如果你只是想在MATLAB里做逻辑判断绝不能这样写但在YALMIP建模里这就是标准写法。4.3 CHP机组运行约束与热电耦合关系CHP机组最关键的约束是热电耦合关系。燃气轮机产生电功率的同时回收余热两者满足以下关系% H_chp 与 P_chp 通过热电比耦合 H_chp CHP_rh * P_chp; % CHP_rh 为热电比这个约束把热平衡和电平衡联系在了一起。如果热电比是固定的模型是线性约束如果采用变热电比模型约束会变成非线性需要引入额外的分段线性化处理。代码基于固定热电比的简化模型这也是绝大多数文献采用的简化方式。CHP机组的出力上下限和爬坡约束% CHP机组出力上下限 Constraints [Constraints, ... P_chp_min P_chp P_chp_max]; % 爬坡约束相邻时段出力变化幅度不能超过限值 for t 2:24 Constraints [Constraints, ... -delta_Chp P_chp(t) - P_chp(t-1) delta_Chp]; end为什么用循环来写爬坡约束因为爬坡约束是相邻时段之间的耦合关系无法像容量约束那样用一个向量表达式整体描述。当然也可以把约束写成矩阵形式但对可读性要求高的复现代码来说循环写法更直观而且24个时段的循环开销非常小完全不影响求解性能。4.4 储能SOC约束储能的SOC动态方程是时序耦合的必须从第2个时段开始递推% 电储能SOC动态方程 Constraints [Constraints, SOC_e(1) SOC_e_initial ... (P_cha(1) * eta_cha_e - P_dis(1) / eta_dis_e) / Cap_e]; for t 2:24 Constraints [Constraints, SOC_e(t) SOC_e(t-1) ... (P_cha(t) * eta_cha_e - P_dis(t) / eta_dis_e) / Cap_e]; end % SOC上下限约束 Constraints [Constraints, SOC_e_min SOC_e SOC_e_max];这里有个容易出错的细节充电效率η_cha和放电效率η_dis所在的位置不同。充电时储能吸收的功率需要乘以效率才能转化为SOC增量放点时储能释放到电网的功率除以效率才是SOC减少量。很多复现代码跑出来的SOC曲线越界或最终SOC不为0多半是这个效率位置写反了。为了避免储能同时充放电这种不物理的情况可以加二进制变量约束也可以用一组线性不等式强制充放电状态互斥% 防止同时充放电使用大M法降低复杂度 M 1e3; u_cha binvar(1, 24); % 充电状态二进制变量 u_dis binvar(1, 24); % 放电状态二进制变量 Constraints [Constraints, u_cha u_dis 1]; Constraints [Constraints, P_cha u_cha * P_cha_max]; Constraints [Constraints, P_dis u_dis * P_dis_max];引入二进制变量后问题从线性规划变成混合整数线性规划求解时间会明显增加但它保证了结果的物理合理性这一步不该省。4.5 目标函数与求解调用目标函数定义上面已经给出了公式代码里对应如下% 成本项 Cost_grid sum(grid_price .* P_grid); % 购电成本 Cost_gas sum(C_gas * (P_chp / eta_chp_e H_gb / eta_gb)); % 燃气成本 Objective Cost_grid Cost_gas; % 求解 options sdpsettings(verbose, 1, solver, cplex, ... showprogress, 1, debug, 0); sol optimize(Constraints, Objective, options);求解完成后最好加一段结果检查判断是否真的收敛到了最优解if sol.problem 0 disp(求解成功); else disp(求解出现问题); sol.info endsol.problem 0是YALMIP返回成功状态的标志非零值分别对应不可行、数值问题、求解器未安装等不同错误排查时先看这个值能省不少时间。5. 实测中的坑变量维度、求解器精度、SOC曲线异常复现代码最花时间的不是写代码而是排查那些莫名其妙的问题。我把自己踩过的坑总结出来希望对你有帮助。5.1 维度不一致导致的约束崩溃YALMIP对维度非常敏感。如果你定义的P_pv是1×24行向量但P_grid是24×1列向量两者相加时MATLAB会自动广播但YALMIP生成的约束可能是24×1也可能报维度不匹配。这个坑的典型表现是solve之前一切正常一求解就报Constraint is not satisfied或维度错误。排查思路是先检查所有变量定义是否一致。保险起见所有决策变量我统一用sdpvar(1, 24)行向量所有外部参数也统一定义为1×24行向量。如果数据从Excel或CSV读入读完后强制转成行向量P_load P_load(:); % 强制转成行向量这个小技巧能在很大概率上避免维度导致的奇怪报错。5.2 求解器返回“不可行”的定位方法模型不可行说明约束之间互相矛盾但问题不会告诉你矛盾在哪个约束。我的做法是把约束分成几组逐组放开测试。先只保留能量平衡约束求解成功再叠加上下限约束成功后再加SOC约束一层层往上加直到出现不可行为止。这样能快速定位是哪一类约束导致的问题。最常见的不可行原因是储能SOC初始值和终止值设置不合理。比如说SOC初始是0.5但负荷曲线和出力曲线决定了储能白天大量放电到晚上SOC被压到0.1以下而SOC下限是0.2那就必然不可行。解决方法是放松SOC下限或者改变初始SOC或者允许储能在一天结束后SOC不等于初始值具体取决于所复现文献的设定。如果文献要求SOC日始日终相等即SOC(24) SOC_initial那要特别注意储能容量是否够用否则无解。5.3 储能SOC曲线异常振荡有一次跑完代码电储能SOC曲线在0.3到0.9之间锯齿状振荡非常不自然。检查后发现是电价在峰谷交替时段目标函数倾向于让储能频繁切换充放电状态来套利而我在约束里没有加充放电状态互斥的二进制变量导致SOC出现高频抖动。解决方法是上面提过的增加二进制状态变量强制u_cha u_dis 1。加了之后振荡消失结果也更有物理意义。代价是求解时间从不到1秒涨到3到5秒但对24小时规模的模型来说完全可接受。5.4 求解器相关的问题使用CPLEX时如果提示没有许可证会直接报错但这属于环境问题。另一类问题是数值尺度差异过大——储能容量可能是100 kWh而购电功率峰值只有几十kW两者相差不大但如果电价单位是元/MWh购电成本可能达到10的4次方量级目标函数数值过大可能让GBD通用分支定界算法收敛变慢。解决办法是把目标函数各项统一为元/kWh并检查所有变量的数值范围尽量控制在1e-3到1e3之间。6. 结果怎么看用三条曲线判断优化是否“复现成功”跑完代码出了图到底对不对这个判断标准比大多数人想的更清楚。我会按以下几步验证结果。6.1 检查能量平衡约束是否闭合既然是等式约束理论上最优解肯定严格满足平衡但因为求解器数值误差实际结果可能存在极小的偏差。可以把约束左侧和右侧的值分别算出来比较差值residual (P_pv P_wt P_chp P_dis P_grid) - (P_load P_cha); disp(max(abs(residual)));如果这个残差小于1e-6说明求解结果数值可靠。6.2 看典型日出力曲线plot_results.m里用area画电功率平衡的堆叠面积图从上到下依次是光伏、风电、CHP、储能放电、电网购电最上面叠加电负荷曲线。一个健康的调度结果应该具备这些特征光伏出力的时段电网购电明显减少储能充电大多发生在电价低谷时段放电发生在电价高峰时段CHP机组在热负荷高的时段出力更大发挥热电联产优势。如果看到储能逆着电价充放电——电价高时充电、电价低时放电那说明目标函数或约束写反了需要回到代码检查。热功率曲线方面主要看燃气锅炉是否只在CHP余热不足时启动。如果燃气锅炉全天满发那说明CHP装机太小或热负荷太大系统的“热电联供”优势没有体现出来。6.3 对比文献数值所谓“完美复现”最直接的验证是把运行成本总值、各设备累计出力和文献给的结果做对比。因为参数和电价可能不完全一致数值不完全相等很正常但趋势应该一致购电成本占比、燃气成本占比应该在合理区间CHP机组的利用小时数也应该在工程常识范围内。我通常会把最终总运行成本打印出来再手动估算一个大致区间如果代码算出来几万块钱而估算值只有几千那目标函数单位多半出了问题。7. 代码扩展方向从“能跑”到“能用”最后分享一些个人经验这套基础代码可以往几个方向扩展改造成本都不高。第一个方向是加入碳交易成本。多能互补微网优化经常和低碳经济挂钩此时目标函数里除了运行成本还要增加碳排放成本项等于在购电和燃气成本之外按碳排放量乘碳价。因为碳排放量和出力是线性关系这种扩展只需要在目标函数里加一行对模型结构和求解难度几乎没有影响。第二个方向是改成两阶段鲁棒优化或随机优化。具体做法是把光伏、风电出力从确定值改成不确定集合第一阶段做日前计划第二阶段做实时调整。这种扩展要把约束结构改成MPEC形式工作量大得多但这套确定性代码作为基础框架仍然适用。第三个方向是增加需求响应负荷。把一部分可转移负荷作为决策变量并入电功率平衡约束并加上转移量和转移时段的限制。这样目标函数不变约束增加几条线性不等式改动也不算复杂。我自己现在还在另一个项目里把这套代码往多微网互联方向扩展核心改动是把每个微网当成一个独立节点微网之间增加联络线功率变量和平衡约束思路是一样的。上面的经验就是先让模型跑通再考虑算法复杂度先保证结果合理再考虑是否引入不确定性。在具体操作中有一点体会很深注释写得越详细后期改代码越轻松。当时为了赶进度写的简化代码过了一个月再看连自己都要重新推演一遍。这套代码我把每个变量的单位、每个约束的物理含义、每个公式对应文献的哪一处都写进了注释整个过程确实更繁琐但后续调整参数、替换数据时省下的时间远超过当时写注释的投入。