1. 先说清楚这项目到底要解决什么问题1.1 微电网、综合能源系统与优化调度三者到底怎么咬合拿到这个题目的时候我第一反应是“又一套经典组合拳”。微电网、综合能源系统、Matlab、CPLEX凡是做过电力系统或能源方向课题的人几乎都能在这几个词里找到自己熬夜调模型的影子。我最初接手这个项目时甲方给的需求很朴素一个混合了光伏、风机、柴油发电机、电储能、燃气热电联产机组和储热罐的小型园区微电网要求在规定运行约束下把一天24小时内的运行成本压到最低同时保证供电和供热都不能断。这个需求听起来简单但真正拆开之后就会发现它已经不只是“微电网调度”一个问题而是“综合能源系统优化”的典型场景。园区里既有电力负荷又有热力负荷电可以由光伏、风机、柴油机、电网和储电共同承担热可以由燃气热电联产机组、电锅炉和储热罐协调供给。电和热之间通过热电联产和电锅炉互相耦合牵一发而动全身。微电网想要经济高效地运行就不能再像传统方式那样“电是电、热是热”分开考虑而是要用一套模型同时描述两类能源的产、储、用全过程。这类问题的学名叫“电-热综合能源系统日前优化调度”落到工程上就是一个典型的混合整数线性规划问题。我在这类项目里常用的技术路线是用Matlab做建模和数据处理把物理系统的约束翻译成线性不等式然后用CPLEX做底层求解器。Matlab负责“讲人话”CPLEX负责“算得快”。这个组合之所以成为业界常规方案是因为Matlab生态里的数据分析和可视化工具实在太顺手而CPLEX在求解大规模混合整数线性规划时稳定性和速度都足够可靠尤其在整数变量较多、约束矩阵规模大的场景下比某些开源求解器要省心得多。这个项目适合谁来参考如果你正在做微电网调度、园区综合能源规划、储能容量配置或者单纯想搞明白“YALMIPCPLEX到底怎么建模”这篇文章都能给你一条可以直接抄作业的路线。我会把从数学建模、代码实现到调试心得的完整过程都摊开来讲包括遇到的坑和应对方式。1.2 为什么这个问题会变成Mixed Integer Linear Programming很多人第一次接触这类项目时会困惑为什么非得用“混合整数线性规划”直接做非线性优化不行吗其实关键在于设备运行特性的数学表达方式。以柴油发电机为例它有一个最低技术出力。如果一台机组的最小出力是30%额定功率那么它的输出功率只有两种状态要么是0要么是在30%~100%额定功率之间的连续值。这个“要么停机、要么在一定区间内运行”的特性天然就需要一个0/1整数变量来表征启停状态。燃气热电联产机组也有类似的启停逻辑甚至还有“开机最小持续时间”这类时间耦合约束这些都必须用整数变量配合线性不等式来描述。另一方面储能设备的电量递推关系、功率平衡方程、设备出力上下限这些本身都是线性关系。电池的充电效率、放电效率在固定工况下也可以近似为常数于是整个问题中连续变量的约束都能用线性不等式表达。最终模型里同时出现了连续变量各设备出力、储能SOC、购售电功率和整数变量机组启停、储能充放电状态目标函数是运行成本的最小化约束全部线性这就是标准的混合整数线性规划。我还在早期版本里尝试过把柴油机燃料成本曲线拟合成二次函数做成混合整数二次规划。结果发现在CPLEX里求解时二次项带来的计算负担明显增加而且工程上燃料成本曲线在运行区间内用分段线性逼近已经足够精确。所以在正式项目里我又回到了MILP框架用分段线性函数拟合非线性成本曲线既保留了精度又显著提升了求解速度。这一点后面在模型部分会详细讲。2. 数学模型把能量系统翻译成计算机能懂的语言2.1 目标函数一天的运行成本明细优化模型的第一步就是定目标函数。在这个项目里系统需要支付的成本大致有五块从上级电网购电的费用、柴油发电机燃料成本、燃气热电联产机组的燃料成本、储能系统充放电产生的退化损耗成本以及设备的启停惩罚成本。如果园区有富余电量还可以反向卖给电网这部分售电收益要从总成本里扣掉。写成数学形式就是min sum_{t1}^{24} [ C_elec(t)*P_buy(t) - S_elec(t)*P_sell(t) C_diesel*F_diesel(t) C_gas*F_chp(t) C_bat*(P_ch(t) P_dis(t)) C_start*delta_diesel(t) ]其中各个符号的含义如下P_buy(t)、P_sell(t)分别是t时段从电网购电和向电网售电的功率F_diesel(t)和F_chp(t)是柴油机和热电联产机组在t时段的燃料消耗量P_ch(t)、P_dis(t)是储能的充、放电功率delta_diesel(t)是柴油机的启动标志等于1时表示该时段机组启动。这里有个容易忽略的细节储能的退化成本。很多入门教程只考虑“买电便宜还是放电便宜”却忘了电池本身循环寿命有限频繁充放电会带来实际的寿命折损。我在目标函数里给充放电功率各加了一个很小的单位成本系数相当于人为地给“折腾电池”施加惩罚让优化结果不会出现设备疯狂来回充放的病态解。这个系数不需要很大只要大到能抑制无意义的充放电循环又小到不影响正常的削峰填谷策略就行。燃料成本的处理也需要注意。柴油机的耗量特性是出力的一次函数可以写成 F_diesel(t) a * P_diesel(t) b * u_diesel(t)其中u_diesel(t)是启停状态整数变量。这里的b是空载损耗系数直接乘在状态变量上天然就是线性的。热电联产机组类似但它的耗气量还要同时反映电出力和热出力的耦合关系要单独建模这个放到约束条件里一起说。2.2 约束条件功率平衡、设备出力上下限、储能电量递推、热电联产目标函数定完接下来是约束条件。约束决定了解的可行域一个约束写错整个结果就可能直接失真。我在这个项目里把约束分成了五组下面逐一说明。第一组是电功率平衡约束。无论任何时候系统内的电力供应和用电需求必须相等P_buy(t) - P_sell(t) P_pv(t) P_wt(t) P_diesel(t) P_chp_e(t) P_dis(t) P_load(t) P_ch(t) P_eb(t) P_hp(t)等式左边是各类电源的出力等式右边是电负荷加上电锅炉、热泵、储能充电等用电设备的消耗。这里的P_pv和P_wt在日前调度中是给定的预测值不作为决策变量因为光伏和风机在无弃光弃风情况下是“发多少算多少”的不可调度电源。第二组是各设备的出力上下限约束。柴油机、热电联产机组、储能充放电功率、电锅炉产热功率都要限制在额定范围内。特别要强调的是储能电池的充放电功率不能同时为正所以需要引入一个充电状态整数变量eta_bat(t)写成P_ch(t) P_ch_max * eta_bat(t) P_dis(t) P_dis_max * (1 - eta_bat(t))第三组是储能电量递推约束。电池的荷电状态按小时递推SOC(t1) SOC(t) eta_ch * P_ch(t) * delta_t - P_dis(t) / eta_dis * delta_t SOC_min SOC(t) SOC_max SOC(0) SOC_initial, SOC(24) SOC_final_mindelta_t取1小时。注意充放电效率在SOC递推式的不同位置充电时是“消了电、存了更多”所以乘效率放电时是“放出了电、减少了库存”所以除效率。很多人容易在这两个效率的位置上搞反导致SOC凭空多出来或者少下去。储热罐的递推公式和电池完全同构只是单位换成热量、效率换成储热损失系数不再赘述。第四组是热电联产机组的电热耦合约束。燃气热电联产机组的电出力和热出力遵循固定热电比同时受总产能限制P_chp_h(t) COP_chp * P_chp_e(t) P_chp_e(t) P_chp_e_max * u_chp(t) P_chp_e(t) P_chp_e_min * u_chp(t)这里的COP_chp是热电联产机组的热电比系数表示每发一度电能同时产生多少热量。这种线性等式约束写法前提是机组运行在热电比基本恒定的额定工况区。如果要做更精细的变工况建模就得把工作区间划分成若干区域做分段线性化这个项目里没有用到但值得知道它的存在。第五组是热功率平衡约束。热负荷由热电联产机组、电锅炉、储热罐共同满足P_chp_h(t) P_eb_h(t) P_tes_dis(t) H_load(t) P_tes_ch(t)别小看这一条热平衡约束。很多微电网项目做到最后发现“电调好了但热对不上”就是因为把电和热当成两个独立系统分开优化。综合能源系统的核心价值恰恰在于电热耦合这个等式是整套模型里最能体现“综合能源”四个字的地方。2.3 决策变量的选择与为什么要引入整数变量决策变量的选取直接决定模型的规模和求解难度。这个项目里的连续决策变量包括购电功率、售电功率、柴油机出力、热电联产机组电出力和热出力、电锅炉产热功率、储能充放电功率、储热罐充放热功率、各时段SOC。整数决策变量则包括柴油机启停状态、热电联产机组启停状态、储能充放电状态。为什么要引入整数变量前面已经提到了设备工作特性。这里补充一个经常被忽略的点设备启停状态变量还会引入“最小启停时间约束”。比如柴油机要求启动后至少连续运行2小时写成约束就是u_diesel(t) - u_diesel(t-1) u_diesel(tau)其中tau取t到t最小运行时间-1之间的所有时段。这个约束的工程意义是防止优化结果让机组每小时启停一次那在实际运行中既不安全也不经济。类似的时间耦合约束会让模型从“单时段独立优化”升级为“跨时段联合优化”计算复杂度会上一个台阶但也让结果真正具备可执行性。在构建模型时我还会额外注意一个细节购电和售电不能同时发生。这个可以通过价格信号自然避免——只要售电价始终低于购电价优化结果就不会出现“低价买高价卖”的套利循环。但为了保险起见我还是加了整数变量m_grid(t)来控制购售电互斥。这个变量在初期调试时可以去掉它会明显加快求解速度等到模型跑通后再加上做最终验证。3. Matlab里怎么把模型建起来并让CPLEX跑起来3.1 环境准备与工具箱配置我在这个项目里选择的软件组合是Matlab R2021a YALMIP CPLEX 12.10。YALMIP是一个免费的Matlab建模工具箱它最大的价值在于把优化模型从“人工拼矩阵”中解放出来——你只需要用接近数学表达式的语法把决策变量、约束和目标函数描述清楚它自动生成CPLEX需要的标准矩阵形式。这对工程技术人员来说简直太友好了相当于把“手写汇编”升级成“写高级语言”。安装配置的步骤不多但有几个细节值得注意。CPLEX解压后要确保把cplex\matlab\文件夹加入Matlab路径同时把对应的dll文件目录加入系统PATH。我遇到过好几次“明明装好了CPLEX但YALMIP调用时一直报找不到求解器”的情况最后发现都是路径没加全。调试方法很简单在Matlab里执行yalmiptest命令它会列出所有可用的求解器及其状态。如果CPLEX显示“not found”就去检查路径配置重点确认cplex和cplexlink这两个文件是否真的在Matlab路径里。3.2 数据初始化设备参数和负荷/新能源曲线模型建好之前先得有一份完整且靠谱的数据。我习惯在模型代码前专门放一个“数据区”把所有设备参数和预测曲线集中管理这样后续做参数敏感性分析时只需要改这一块不需要动后面的建模代码。设备参数包括柴油机额定功率、最小出力、燃料成本系数热电联产机组额定电功率、电热比、燃料成本系数储能电池容量、SOC上下限、充放电效率、最大充放电功率储热罐容量、热损失系数电网购售电价。负荷曲线和新能源出力曲线则是24个时点的数组可以是实测历史数据也可以是预测数据。我在项目里用的典型参数大致如下设备/参数数值单位光伏额定容量500kW风机额定容量300kW柴油机额定功率400kW热电联产机组电功率上限300kW热电联产热电比1.2-储能电池容量800kWh储能最大充/放电功率200kW峰时电价1.2元/kWh谷时电价0.3元/kWh天然气价格2.5元/m^3这套参数的合理性在于模拟了一个有一定规模但不算太大的园区微电网。光伏和风机加起来超过柴油机容量说明新能源渗透率很高储能在峰谷电价差1.2元对0.3元的情况下有充分的套利空间。热负荷曲线设定为日间低、夜间高贴合北方园区冬季供暖的实际形态。3.3 用YALMIP建模并调用CPLEX求解数据准备完毕接下来就是核心步骤用YALMIP把第二节的数学模型翻译成代码。我用一段简化但可运行的核心代码来说明整体思路。注意这里面省略了部分细节但结构完整可以直接迁移到自己的项目里%% 数据加载示例 T 24; dt 1; P_load [ ... 24个负荷值 ... ]; H_load [ ... 24个热负荷值 ... ]; P_pv [ ... 24个光伏预测值 ... ]; P_wt [ ... 24个风机预测值 ... ]; %% 决策变量 P_buy sdpvar(T,1); % 购电功率 P_sell sdpvar(T,1); % 售电功率 P_diesel sdpvar(T,1); % 柴油机出力 u_diesel binvar(T,1); % 柴油机启停状态 P_chp_e sdpvar(T,1); % 热电联产电出力 P_chp_h sdpvar(T,1); % 热电联产热出力 u_chp binvar(T,1); % 热电联产启停状态 P_ch sdpvar(T,1); % 储能充电功率 P_dis sdpvar(T,1); % 储能放电功率 eta_bat binvar(T,1); % 储能充/放状态 SOC sdpvar(T1,1); % 荷电状态 P_eb sdpvar(T,1); % 电锅炉电功率 P_tes_ch sdpvar(T,1); % 储热罐充热功率 P_tes_dis sdpvar(T,1); % 储热罐放热功率 SOC_tes sdpvar(T1,1); % 储热罐蓄热状态 P_start_diesel binvar(T,1); % 柴油机启动标志 %% 目标函数 C_buy 1.2 * ones(T,1); % 分时电价示意实际应分段 C_sell 0.3 * ones(T,1); C_diesel_fuel 0.6; % 柴油单位燃料成本系数 C_gas 2.5; % 天然气单价折算系数 C_bat 0.02; % 储能退化成本系数 Objective sum(C_buy.*P_buy - C_sell.*P_sell ... C_diesel_fuel.*P_diesel C_gas.*(P_chp_e*0.35) ... C_bat.*(P_ch P_dis) ... 50 * P_start_diesel); %% 约束条件 Constraints []; % 电功率平衡 Constraints [Constraints, ... P_buy - P_sell P_pv P_wt P_diesel P_chp_e P_dis ... P_load P_ch P_eb]; % 热功率平衡 Constraints [Constraints, ... P_chp_h 0.95*P_eb P_tes_dis H_load P_tes_ch]; % 柴油机出力限幅 Constraints [Constraints, ... 0.3*400*u_diesel P_diesel 400*u_diesel]; % 热电联产出力限幅与热电比 Constraints [Constraints, ... 0.3*300*u_chp P_chp_e 300*u_chp]; Constraints [Constraints, P_chp_h 1.2 * P_chp_e]; % 储能约束 Constraints [Constraints, ... 0 P_ch 200 * eta_bat]; Constraints [Constraints, ... 0 P_dis 200 * (1 - eta_bat)]; Constraints [Constraints, SOC(2:T1) SOC(1:T) ... 0.95*P_ch*dt - P_dis/0.95*dt]; Constraints [Constraints, 100 SOC 800]; Constraints [Constraints, SOC(1) 400, SOC(T1) 400]; % 储热罐约束结构类似储能略 % ... %% 求解 ops sdpsettings(solver,cplex,verbose,2); result optimize(Constraints, Objective, ops);注意热电联产的气耗量我用了电出力的线性折算系数0.35表示每千瓦时电出力对应的天然气耗量。实际项目中这个系数要由机组性能曲线拟合得出不同机型差异很大。启动惩罚系数50加在柴油机启动标志上是为了避免机组频繁启停——这个值需要根据实际运行成本标定设大了会抑制正常启停设小了起不到惩罚作用。求解完成后用value()函数提取各决策变量的数值结果再用Matlab自带的绘图函数把24小时的功率曲线画出来。可视化这一步千万别省因为只看一堆数字很难发现“某时段功率平衡被悄悄破坏”这种隐蔽问题。3.4 求解器设置和结果整理CPLEX的求解器参数直接通过YALMIP的sdpsettings传递。我常用的设置有ops sdpsettings(solver,cplex,verbose,1); ops.cplex.mip.tolerances.mipgap 0.001; % 1%的MIP间隙即可接受 ops.cplex.timelimit 120; % 限制最大求解时间为120秒 ops.cplex.mip.strategy.startalgorithm 2; % 使用对偶单纯形做根节点松弛其中mipgap是最关键的一个参数。对于日前调度这类需要每天重复执行100次以上的场景没必要追求0.0001的极端精度把间隙放宽到0.01左右CPLEX求解时间往往能缩短一个数量级。我在这个项目里做过对比同样一套模型mipgap设0.05时求解时间在10秒以内设0.0001时能拉到好几分钟而目标函数值的差别不到0.3%。工程上省下的时间完全可以覆盖这点精度损失。结果整理时我固定输出四张图各电源出力堆叠图、储能SOC变化曲线、购电功率与电价对照图、热功率平衡图。这四张图基本能把一天的运行逻辑一眼讲清楚无论是自己复盘还是向项目团队汇报都非常有用。4. 从仿真结果看系统优化与经济收益4.1 基准场景的调度结果长什么样先给出一组基准场景的典型结果。假设当天光照充足光伏出力在10点到14点达到峰值峰时电价同样落在9点到15点之间。优化结果通常会呈现一个很明显的规律白天光伏大发时系统主动压低购电功率甚至把多余电量卖给电网储能则在午间光伏高峰时段充电等到晚高峰电价高的时候再放电。具体到数字我在一组典型数据下得到的结果是全天总购电量约4200 kWh其中峰时购电不到500 kWh大部分购电集中在谷时和平时段柴油机只在夜间热负荷尖峰时段启动了4小时热电联产机组白天保持满发产生的热量优先供给热负荷多余热量存入储热罐夜间热负荷上来时再释放。储能在13点到15点充电近400 kWh在18点到21点放电约370 kWh完成了一个完整的“低储高发”循环。对系统的全天天运行成本来看相比完全不优化、按经验规则运行的基准方案日运行成本下降了约18%。这个数字很能说明优化的价值——不是靠砍负荷而是靠合理地“错峰”和“多能互补”。4.2 并网与孤岛运行模式对比我再把并网模式和孤岛模式各跑了一遍。并网模式下系统可以随时从电网买电因此当光伏不足、储能放空时就直接从电网购电兜底经济性最优。孤岛模式下系统与电网完全断开一切电力缺口都必须靠柴油机、热电联产和储能来填补。结果差异非常明显。孤岛模式下柴油机几乎全天都在运行因为仅靠光伏、风机和储能在晚间根本撑不住负荷。热电联产机组也始终保持电出力因为它的热出力也是刚需。全天燃料成本是并网模式的两倍还多。这个对比验证了一个工程常识微电网的经济性极度依赖电网互动脱离电网硬抗只适合应急场景。在做孤岛验证时我还特意检查了储能SOC是否在每个时段都满足最低容量约束发现傍晚负荷高峰时段SOC会降到接近下限说明储能容量在这种模式下的裕度非常紧张。4.3 储能容量和热电联产对总成本的影响最后做了一组敏感性分析。把储能容量从400 kWh逐步提高到2000 kWh观察总成本变化曲线。结果是储能容量从400增加到1000时成本下降非常显著因为更大的容量允许系统在谷电时段存更多电、在峰时释放但容量超过1000之后成本下降趋缓边际效益明显递减。原因在于峰谷套利的收益受一天充放电循环次数和功率限制约束容量再大一天也只能充放一次多出来的容量大多数时间闲置反而增加了初始投资成本虽然这里没建投资模型但从运行成本角度看已经能看到回报率下降的趋势。热电联产机组的作用更加直接。我把热电联产的电出力上限从200 kW扫到500 kW热负荷不变系统总运行成本几乎线性下降。这说明在同时存在电、热负荷的综合能源系统里热电联产的扩容能有效替代电锅炉和柴油机的运行成本经济性非常强劲。这也是为什么很多园区级项目愿意优先上热电联产机组而不是扩大储能的原因。5. 实操中踩过的坑与排查经验5.1 模型不可行时要怎么查约束我做这个项目过程中遇到最多的问题就是模型无解。CPLEX返回“infeasible”时很多初学者的第一反应是瞎改参数这完全错误。正确做法是启用CPLEX的不可行性分析功能在YALMIP里设置ops.cplex.lpmethod 0; ops.cplex.emphasis.mip 4; % 让CPLEX生成不可行子系统报告CPLEX会输出一个IISIrreducible Inconsistent Subsystem不可约不一致子系统它会指出是哪几条约束联立起来导致了矛盾。我遇到的一个典型情况是储能SOC的初始值和最终值都设定为400 kWh但某天负荷极高、储能放电量很大中间某时段即使满功率充电也来不及恢复到400。这种情况下IIS会指向储能递推约束和SOC上下限约束。解决办法很简单把SOC(final)的要求从“等于400”放宽为“大于等于300”问题立即解除。另一个血泪教训是热负荷为0的时段不要设置储热罐必须填充到某个值。热力系统的惯性没有电力系统那么强硬性约束很容易造成无解。建议把储热罐的最终蓄热状态也做成“大于等于某值”而不是“等于某值”。5.2 CPLEX求解速度慢和收敛抖动的处理求解速度慢是个令人头疼的问题特别是引入最小启停时间约束和启动变量之后整数变量数量会明显增加。我处理速度问题有三个办法。第一简化不需要的整数变量。如果储能充放电互斥通过价格机制已经天然成立可以把eta_bat这个整数变量去掉直接给P_ch和P_dis加上一个很小的互相限制甚至干脆让它们同时为正但在目标函数里惩罚充放电总量。实测下来去掉这个整数变量后求解速度提升20%以上。第二设置合理的MIP间隙。前面提过把mipgap从默认的1e-4放宽到0.01求解时间可以从十几分钟降到几秒而且成本差异极小。很多商业化的调度系统上线后跑的就是这组宽松参数。第三给CPLEX提供一个初始可行解。用上一时段的解作为本时段的冷启动起点YALMIP支持通过assign给变量赋初值CPLEX会在求解过程中利用这个初值加速整数分支收敛。这在做滚动优化时效果尤其明显。5.3 YALMIP报错的几种典型原因最后汇总一下YALMIP本身的常见报错。遇到过YALMIP Data in constraint not matching多半是某个变量的维度不对。比如P_load是24×1数组而P_buy被声明成1×24减号运算时光标对齐失败。此时用size()逐个检查所有参与约束的变量即可。遇到过No suitable solver found那是路径配置或求解器安装有问题yalmiptest一下就知道。还遇到过This is a non-convex problem这种情况几乎总是因为某个变量被乘了两次比如把SOC(t)和SOC(t1)写进了同一项或者用了norm()这类非线性函数而CPLEX只支持线性模型。写在最后的一点个人体会这个项目做完我对微电网综合能源优化的理解比之前深了不少。最核心的收获是模型一定要从物理过程出发别一开始就想着怎么用高级算法。很多问题卡在“约束没写对”而不是“求解器不够快”。CPLEX非常强大但它只能帮你在给定模型里找最优解模型本身错了他也无能为力。另一个体会是储能退化成本、机组启停惩罚这些“小参数”千万别省它们在约束解的质量上起的作用远超你的直觉。最后的最后多花时间把数据管线做扎实把数据格式统一、曲线可视化做得顺眼比花时间调求解器参数更能提升整个项目的产出效率。这个项目后续还可以扩展到电-气耦合、阶梯电价、需求侧响应等方向每一块都是另一场硬仗也是这条路上最有意思的部分。
1. 先说清楚这项目到底要解决什么问题1.1 微电网、综合能源系统与优化调度三者到底怎么咬合拿到这个题目的时候我第一反应是“又一套经典组合拳”。微电网、综合能源系统、Matlab、CPLEX凡是做过电力系统或能源方向课题的人几乎都能在这几个词里找到自己熬夜调模型的影子。我最初接手这个项目时甲方给的需求很朴素一个混合了光伏、风机、柴油发电机、电储能、燃气热电联产机组和储热罐的小型园区微电网要求在规定运行约束下把一天24小时内的运行成本压到最低同时保证供电和供热都不能断。这个需求听起来简单但真正拆开之后就会发现它已经不只是“微电网调度”一个问题而是“综合能源系统优化”的典型场景。园区里既有电力负荷又有热力负荷电可以由光伏、风机、柴油机、电网和储电共同承担热可以由燃气热电联产机组、电锅炉和储热罐协调供给。电和热之间通过热电联产和电锅炉互相耦合牵一发而动全身。微电网想要经济高效地运行就不能再像传统方式那样“电是电、热是热”分开考虑而是要用一套模型同时描述两类能源的产、储、用全过程。这类问题的学名叫“电-热综合能源系统日前优化调度”落到工程上就是一个典型的混合整数线性规划问题。我在这类项目里常用的技术路线是用Matlab做建模和数据处理把物理系统的约束翻译成线性不等式然后用CPLEX做底层求解器。Matlab负责“讲人话”CPLEX负责“算得快”。这个组合之所以成为业界常规方案是因为Matlab生态里的数据分析和可视化工具实在太顺手而CPLEX在求解大规模混合整数线性规划时稳定性和速度都足够可靠尤其在整数变量较多、约束矩阵规模大的场景下比某些开源求解器要省心得多。这个项目适合谁来参考如果你正在做微电网调度、园区综合能源规划、储能容量配置或者单纯想搞明白“YALMIPCPLEX到底怎么建模”这篇文章都能给你一条可以直接抄作业的路线。我会把从数学建模、代码实现到调试心得的完整过程都摊开来讲包括遇到的坑和应对方式。1.2 为什么这个问题会变成Mixed Integer Linear Programming很多人第一次接触这类项目时会困惑为什么非得用“混合整数线性规划”直接做非线性优化不行吗其实关键在于设备运行特性的数学表达方式。以柴油发电机为例它有一个最低技术出力。如果一台机组的最小出力是30%额定功率那么它的输出功率只有两种状态要么是0要么是在30%~100%额定功率之间的连续值。这个“要么停机、要么在一定区间内运行”的特性天然就需要一个0/1整数变量来表征启停状态。燃气热电联产机组也有类似的启停逻辑甚至还有“开机最小持续时间”这类时间耦合约束这些都必须用整数变量配合线性不等式来描述。另一方面储能设备的电量递推关系、功率平衡方程、设备出力上下限这些本身都是线性关系。电池的充电效率、放电效率在固定工况下也可以近似为常数于是整个问题中连续变量的约束都能用线性不等式表达。最终模型里同时出现了连续变量各设备出力、储能SOC、购售电功率和整数变量机组启停、储能充放电状态目标函数是运行成本的最小化约束全部线性这就是标准的混合整数线性规划。我还在早期版本里尝试过把柴油机燃料成本曲线拟合成二次函数做成混合整数二次规划。结果发现在CPLEX里求解时二次项带来的计算负担明显增加而且工程上燃料成本曲线在运行区间内用分段线性逼近已经足够精确。所以在正式项目里我又回到了MILP框架用分段线性函数拟合非线性成本曲线既保留了精度又显著提升了求解速度。这一点后面在模型部分会详细讲。2. 数学模型把能量系统翻译成计算机能懂的语言2.1 目标函数一天的运行成本明细优化模型的第一步就是定目标函数。在这个项目里系统需要支付的成本大致有五块从上级电网购电的费用、柴油发电机燃料成本、燃气热电联产机组的燃料成本、储能系统充放电产生的退化损耗成本以及设备的启停惩罚成本。如果园区有富余电量还可以反向卖给电网这部分售电收益要从总成本里扣掉。写成数学形式就是min sum_{t1}^{24} [ C_elec(t)*P_buy(t) - S_elec(t)*P_sell(t) C_diesel*F_diesel(t) C_gas*F_chp(t) C_bat*(P_ch(t) P_dis(t)) C_start*delta_diesel(t) ]其中各个符号的含义如下P_buy(t)、P_sell(t)分别是t时段从电网购电和向电网售电的功率F_diesel(t)和F_chp(t)是柴油机和热电联产机组在t时段的燃料消耗量P_ch(t)、P_dis(t)是储能的充、放电功率delta_diesel(t)是柴油机的启动标志等于1时表示该时段机组启动。这里有个容易忽略的细节储能的退化成本。很多入门教程只考虑“买电便宜还是放电便宜”却忘了电池本身循环寿命有限频繁充放电会带来实际的寿命折损。我在目标函数里给充放电功率各加了一个很小的单位成本系数相当于人为地给“折腾电池”施加惩罚让优化结果不会出现设备疯狂来回充放的病态解。这个系数不需要很大只要大到能抑制无意义的充放电循环又小到不影响正常的削峰填谷策略就行。燃料成本的处理也需要注意。柴油机的耗量特性是出力的一次函数可以写成 F_diesel(t) a * P_diesel(t) b * u_diesel(t)其中u_diesel(t)是启停状态整数变量。这里的b是空载损耗系数直接乘在状态变量上天然就是线性的。热电联产机组类似但它的耗气量还要同时反映电出力和热出力的耦合关系要单独建模这个放到约束条件里一起说。2.2 约束条件功率平衡、设备出力上下限、储能电量递推、热电联产目标函数定完接下来是约束条件。约束决定了解的可行域一个约束写错整个结果就可能直接失真。我在这个项目里把约束分成了五组下面逐一说明。第一组是电功率平衡约束。无论任何时候系统内的电力供应和用电需求必须相等P_buy(t) - P_sell(t) P_pv(t) P_wt(t) P_diesel(t) P_chp_e(t) P_dis(t) P_load(t) P_ch(t) P_eb(t) P_hp(t)等式左边是各类电源的出力等式右边是电负荷加上电锅炉、热泵、储能充电等用电设备的消耗。这里的P_pv和P_wt在日前调度中是给定的预测值不作为决策变量因为光伏和风机在无弃光弃风情况下是“发多少算多少”的不可调度电源。第二组是各设备的出力上下限约束。柴油机、热电联产机组、储能充放电功率、电锅炉产热功率都要限制在额定范围内。特别要强调的是储能电池的充放电功率不能同时为正所以需要引入一个充电状态整数变量eta_bat(t)写成P_ch(t) P_ch_max * eta_bat(t) P_dis(t) P_dis_max * (1 - eta_bat(t))第三组是储能电量递推约束。电池的荷电状态按小时递推SOC(t1) SOC(t) eta_ch * P_ch(t) * delta_t - P_dis(t) / eta_dis * delta_t SOC_min SOC(t) SOC_max SOC(0) SOC_initial, SOC(24) SOC_final_mindelta_t取1小时。注意充放电效率在SOC递推式的不同位置充电时是“消了电、存了更多”所以乘效率放电时是“放出了电、减少了库存”所以除效率。很多人容易在这两个效率的位置上搞反导致SOC凭空多出来或者少下去。储热罐的递推公式和电池完全同构只是单位换成热量、效率换成储热损失系数不再赘述。第四组是热电联产机组的电热耦合约束。燃气热电联产机组的电出力和热出力遵循固定热电比同时受总产能限制P_chp_h(t) COP_chp * P_chp_e(t) P_chp_e(t) P_chp_e_max * u_chp(t) P_chp_e(t) P_chp_e_min * u_chp(t)这里的COP_chp是热电联产机组的热电比系数表示每发一度电能同时产生多少热量。这种线性等式约束写法前提是机组运行在热电比基本恒定的额定工况区。如果要做更精细的变工况建模就得把工作区间划分成若干区域做分段线性化这个项目里没有用到但值得知道它的存在。第五组是热功率平衡约束。热负荷由热电联产机组、电锅炉、储热罐共同满足P_chp_h(t) P_eb_h(t) P_tes_dis(t) H_load(t) P_tes_ch(t)别小看这一条热平衡约束。很多微电网项目做到最后发现“电调好了但热对不上”就是因为把电和热当成两个独立系统分开优化。综合能源系统的核心价值恰恰在于电热耦合这个等式是整套模型里最能体现“综合能源”四个字的地方。2.3 决策变量的选择与为什么要引入整数变量决策变量的选取直接决定模型的规模和求解难度。这个项目里的连续决策变量包括购电功率、售电功率、柴油机出力、热电联产机组电出力和热出力、电锅炉产热功率、储能充放电功率、储热罐充放热功率、各时段SOC。整数决策变量则包括柴油机启停状态、热电联产机组启停状态、储能充放电状态。为什么要引入整数变量前面已经提到了设备工作特性。这里补充一个经常被忽略的点设备启停状态变量还会引入“最小启停时间约束”。比如柴油机要求启动后至少连续运行2小时写成约束就是u_diesel(t) - u_diesel(t-1) u_diesel(tau)其中tau取t到t最小运行时间-1之间的所有时段。这个约束的工程意义是防止优化结果让机组每小时启停一次那在实际运行中既不安全也不经济。类似的时间耦合约束会让模型从“单时段独立优化”升级为“跨时段联合优化”计算复杂度会上一个台阶但也让结果真正具备可执行性。在构建模型时我还会额外注意一个细节购电和售电不能同时发生。这个可以通过价格信号自然避免——只要售电价始终低于购电价优化结果就不会出现“低价买高价卖”的套利循环。但为了保险起见我还是加了整数变量m_grid(t)来控制购售电互斥。这个变量在初期调试时可以去掉它会明显加快求解速度等到模型跑通后再加上做最终验证。3. Matlab里怎么把模型建起来并让CPLEX跑起来3.1 环境准备与工具箱配置我在这个项目里选择的软件组合是Matlab R2021a YALMIP CPLEX 12.10。YALMIP是一个免费的Matlab建模工具箱它最大的价值在于把优化模型从“人工拼矩阵”中解放出来——你只需要用接近数学表达式的语法把决策变量、约束和目标函数描述清楚它自动生成CPLEX需要的标准矩阵形式。这对工程技术人员来说简直太友好了相当于把“手写汇编”升级成“写高级语言”。安装配置的步骤不多但有几个细节值得注意。CPLEX解压后要确保把cplex\matlab\文件夹加入Matlab路径同时把对应的dll文件目录加入系统PATH。我遇到过好几次“明明装好了CPLEX但YALMIP调用时一直报找不到求解器”的情况最后发现都是路径没加全。调试方法很简单在Matlab里执行yalmiptest命令它会列出所有可用的求解器及其状态。如果CPLEX显示“not found”就去检查路径配置重点确认cplex和cplexlink这两个文件是否真的在Matlab路径里。3.2 数据初始化设备参数和负荷/新能源曲线模型建好之前先得有一份完整且靠谱的数据。我习惯在模型代码前专门放一个“数据区”把所有设备参数和预测曲线集中管理这样后续做参数敏感性分析时只需要改这一块不需要动后面的建模代码。设备参数包括柴油机额定功率、最小出力、燃料成本系数热电联产机组额定电功率、电热比、燃料成本系数储能电池容量、SOC上下限、充放电效率、最大充放电功率储热罐容量、热损失系数电网购售电价。负荷曲线和新能源出力曲线则是24个时点的数组可以是实测历史数据也可以是预测数据。我在项目里用的典型参数大致如下设备/参数数值单位光伏额定容量500kW风机额定容量300kW柴油机额定功率400kW热电联产机组电功率上限300kW热电联产热电比1.2-储能电池容量800kWh储能最大充/放电功率200kW峰时电价1.2元/kWh谷时电价0.3元/kWh天然气价格2.5元/m^3这套参数的合理性在于模拟了一个有一定规模但不算太大的园区微电网。光伏和风机加起来超过柴油机容量说明新能源渗透率很高储能在峰谷电价差1.2元对0.3元的情况下有充分的套利空间。热负荷曲线设定为日间低、夜间高贴合北方园区冬季供暖的实际形态。3.3 用YALMIP建模并调用CPLEX求解数据准备完毕接下来就是核心步骤用YALMIP把第二节的数学模型翻译成代码。我用一段简化但可运行的核心代码来说明整体思路。注意这里面省略了部分细节但结构完整可以直接迁移到自己的项目里%% 数据加载示例 T 24; dt 1; P_load [ ... 24个负荷值 ... ]; H_load [ ... 24个热负荷值 ... ]; P_pv [ ... 24个光伏预测值 ... ]; P_wt [ ... 24个风机预测值 ... ]; %% 决策变量 P_buy sdpvar(T,1); % 购电功率 P_sell sdpvar(T,1); % 售电功率 P_diesel sdpvar(T,1); % 柴油机出力 u_diesel binvar(T,1); % 柴油机启停状态 P_chp_e sdpvar(T,1); % 热电联产电出力 P_chp_h sdpvar(T,1); % 热电联产热出力 u_chp binvar(T,1); % 热电联产启停状态 P_ch sdpvar(T,1); % 储能充电功率 P_dis sdpvar(T,1); % 储能放电功率 eta_bat binvar(T,1); % 储能充/放状态 SOC sdpvar(T1,1); % 荷电状态 P_eb sdpvar(T,1); % 电锅炉电功率 P_tes_ch sdpvar(T,1); % 储热罐充热功率 P_tes_dis sdpvar(T,1); % 储热罐放热功率 SOC_tes sdpvar(T1,1); % 储热罐蓄热状态 P_start_diesel binvar(T,1); % 柴油机启动标志 %% 目标函数 C_buy 1.2 * ones(T,1); % 分时电价示意实际应分段 C_sell 0.3 * ones(T,1); C_diesel_fuel 0.6; % 柴油单位燃料成本系数 C_gas 2.5; % 天然气单价折算系数 C_bat 0.02; % 储能退化成本系数 Objective sum(C_buy.*P_buy - C_sell.*P_sell ... C_diesel_fuel.*P_diesel C_gas.*(P_chp_e*0.35) ... C_bat.*(P_ch P_dis) ... 50 * P_start_diesel); %% 约束条件 Constraints []; % 电功率平衡 Constraints [Constraints, ... P_buy - P_sell P_pv P_wt P_diesel P_chp_e P_dis ... P_load P_ch P_eb]; % 热功率平衡 Constraints [Constraints, ... P_chp_h 0.95*P_eb P_tes_dis H_load P_tes_ch]; % 柴油机出力限幅 Constraints [Constraints, ... 0.3*400*u_diesel P_diesel 400*u_diesel]; % 热电联产出力限幅与热电比 Constraints [Constraints, ... 0.3*300*u_chp P_chp_e 300*u_chp]; Constraints [Constraints, P_chp_h 1.2 * P_chp_e]; % 储能约束 Constraints [Constraints, ... 0 P_ch 200 * eta_bat]; Constraints [Constraints, ... 0 P_dis 200 * (1 - eta_bat)]; Constraints [Constraints, SOC(2:T1) SOC(1:T) ... 0.95*P_ch*dt - P_dis/0.95*dt]; Constraints [Constraints, 100 SOC 800]; Constraints [Constraints, SOC(1) 400, SOC(T1) 400]; % 储热罐约束结构类似储能略 % ... %% 求解 ops sdpsettings(solver,cplex,verbose,2); result optimize(Constraints, Objective, ops);注意热电联产的气耗量我用了电出力的线性折算系数0.35表示每千瓦时电出力对应的天然气耗量。实际项目中这个系数要由机组性能曲线拟合得出不同机型差异很大。启动惩罚系数50加在柴油机启动标志上是为了避免机组频繁启停——这个值需要根据实际运行成本标定设大了会抑制正常启停设小了起不到惩罚作用。求解完成后用value()函数提取各决策变量的数值结果再用Matlab自带的绘图函数把24小时的功率曲线画出来。可视化这一步千万别省因为只看一堆数字很难发现“某时段功率平衡被悄悄破坏”这种隐蔽问题。3.4 求解器设置和结果整理CPLEX的求解器参数直接通过YALMIP的sdpsettings传递。我常用的设置有ops sdpsettings(solver,cplex,verbose,1); ops.cplex.mip.tolerances.mipgap 0.001; % 1%的MIP间隙即可接受 ops.cplex.timelimit 120; % 限制最大求解时间为120秒 ops.cplex.mip.strategy.startalgorithm 2; % 使用对偶单纯形做根节点松弛其中mipgap是最关键的一个参数。对于日前调度这类需要每天重复执行100次以上的场景没必要追求0.0001的极端精度把间隙放宽到0.01左右CPLEX求解时间往往能缩短一个数量级。我在这个项目里做过对比同样一套模型mipgap设0.05时求解时间在10秒以内设0.0001时能拉到好几分钟而目标函数值的差别不到0.3%。工程上省下的时间完全可以覆盖这点精度损失。结果整理时我固定输出四张图各电源出力堆叠图、储能SOC变化曲线、购电功率与电价对照图、热功率平衡图。这四张图基本能把一天的运行逻辑一眼讲清楚无论是自己复盘还是向项目团队汇报都非常有用。4. 从仿真结果看系统优化与经济收益4.1 基准场景的调度结果长什么样先给出一组基准场景的典型结果。假设当天光照充足光伏出力在10点到14点达到峰值峰时电价同样落在9点到15点之间。优化结果通常会呈现一个很明显的规律白天光伏大发时系统主动压低购电功率甚至把多余电量卖给电网储能则在午间光伏高峰时段充电等到晚高峰电价高的时候再放电。具体到数字我在一组典型数据下得到的结果是全天总购电量约4200 kWh其中峰时购电不到500 kWh大部分购电集中在谷时和平时段柴油机只在夜间热负荷尖峰时段启动了4小时热电联产机组白天保持满发产生的热量优先供给热负荷多余热量存入储热罐夜间热负荷上来时再释放。储能在13点到15点充电近400 kWh在18点到21点放电约370 kWh完成了一个完整的“低储高发”循环。对系统的全天天运行成本来看相比完全不优化、按经验规则运行的基准方案日运行成本下降了约18%。这个数字很能说明优化的价值——不是靠砍负荷而是靠合理地“错峰”和“多能互补”。4.2 并网与孤岛运行模式对比我再把并网模式和孤岛模式各跑了一遍。并网模式下系统可以随时从电网买电因此当光伏不足、储能放空时就直接从电网购电兜底经济性最优。孤岛模式下系统与电网完全断开一切电力缺口都必须靠柴油机、热电联产和储能来填补。结果差异非常明显。孤岛模式下柴油机几乎全天都在运行因为仅靠光伏、风机和储能在晚间根本撑不住负荷。热电联产机组也始终保持电出力因为它的热出力也是刚需。全天燃料成本是并网模式的两倍还多。这个对比验证了一个工程常识微电网的经济性极度依赖电网互动脱离电网硬抗只适合应急场景。在做孤岛验证时我还特意检查了储能SOC是否在每个时段都满足最低容量约束发现傍晚负荷高峰时段SOC会降到接近下限说明储能容量在这种模式下的裕度非常紧张。4.3 储能容量和热电联产对总成本的影响最后做了一组敏感性分析。把储能容量从400 kWh逐步提高到2000 kWh观察总成本变化曲线。结果是储能容量从400增加到1000时成本下降非常显著因为更大的容量允许系统在谷电时段存更多电、在峰时释放但容量超过1000之后成本下降趋缓边际效益明显递减。原因在于峰谷套利的收益受一天充放电循环次数和功率限制约束容量再大一天也只能充放一次多出来的容量大多数时间闲置反而增加了初始投资成本虽然这里没建投资模型但从运行成本角度看已经能看到回报率下降的趋势。热电联产机组的作用更加直接。我把热电联产的电出力上限从200 kW扫到500 kW热负荷不变系统总运行成本几乎线性下降。这说明在同时存在电、热负荷的综合能源系统里热电联产的扩容能有效替代电锅炉和柴油机的运行成本经济性非常强劲。这也是为什么很多园区级项目愿意优先上热电联产机组而不是扩大储能的原因。5. 实操中踩过的坑与排查经验5.1 模型不可行时要怎么查约束我做这个项目过程中遇到最多的问题就是模型无解。CPLEX返回“infeasible”时很多初学者的第一反应是瞎改参数这完全错误。正确做法是启用CPLEX的不可行性分析功能在YALMIP里设置ops.cplex.lpmethod 0; ops.cplex.emphasis.mip 4; % 让CPLEX生成不可行子系统报告CPLEX会输出一个IISIrreducible Inconsistent Subsystem不可约不一致子系统它会指出是哪几条约束联立起来导致了矛盾。我遇到的一个典型情况是储能SOC的初始值和最终值都设定为400 kWh但某天负荷极高、储能放电量很大中间某时段即使满功率充电也来不及恢复到400。这种情况下IIS会指向储能递推约束和SOC上下限约束。解决办法很简单把SOC(final)的要求从“等于400”放宽为“大于等于300”问题立即解除。另一个血泪教训是热负荷为0的时段不要设置储热罐必须填充到某个值。热力系统的惯性没有电力系统那么强硬性约束很容易造成无解。建议把储热罐的最终蓄热状态也做成“大于等于某值”而不是“等于某值”。5.2 CPLEX求解速度慢和收敛抖动的处理求解速度慢是个令人头疼的问题特别是引入最小启停时间约束和启动变量之后整数变量数量会明显增加。我处理速度问题有三个办法。第一简化不需要的整数变量。如果储能充放电互斥通过价格机制已经天然成立可以把eta_bat这个整数变量去掉直接给P_ch和P_dis加上一个很小的互相限制甚至干脆让它们同时为正但在目标函数里惩罚充放电总量。实测下来去掉这个整数变量后求解速度提升20%以上。第二设置合理的MIP间隙。前面提过把mipgap从默认的1e-4放宽到0.01求解时间可以从十几分钟降到几秒而且成本差异极小。很多商业化的调度系统上线后跑的就是这组宽松参数。第三给CPLEX提供一个初始可行解。用上一时段的解作为本时段的冷启动起点YALMIP支持通过assign给变量赋初值CPLEX会在求解过程中利用这个初值加速整数分支收敛。这在做滚动优化时效果尤其明显。5.3 YALMIP报错的几种典型原因最后汇总一下YALMIP本身的常见报错。遇到过YALMIP Data in constraint not matching多半是某个变量的维度不对。比如P_load是24×1数组而P_buy被声明成1×24减号运算时光标对齐失败。此时用size()逐个检查所有参与约束的变量即可。遇到过No suitable solver found那是路径配置或求解器安装有问题yalmiptest一下就知道。还遇到过This is a non-convex problem这种情况几乎总是因为某个变量被乘了两次比如把SOC(t)和SOC(t1)写进了同一项或者用了norm()这类非线性函数而CPLEX只支持线性模型。写在最后的一点个人体会这个项目做完我对微电网综合能源优化的理解比之前深了不少。最核心的收获是模型一定要从物理过程出发别一开始就想着怎么用高级算法。很多问题卡在“约束没写对”而不是“求解器不够快”。CPLEX非常强大但它只能帮你在给定模型里找最优解模型本身错了他也无能为力。另一个体会是储能退化成本、机组启停惩罚这些“小参数”千万别省它们在约束解的质量上起的作用远超你的直觉。最后的最后多花时间把数据管线做扎实把数据格式统一、曲线可视化做得顺眼比花时间调求解器参数更能提升整个项目的产出效率。这个项目后续还可以扩展到电-气耦合、阶梯电价、需求侧响应等方向每一块都是另一场硬仗也是这条路上最有意思的部分。
