简介一套MATLAB电力系统机组组合优化实践资料面向电力系统调度、运行规划研究人员及相关专业学生重点覆盖机组启停决策、最优出力分配与多类约束处理。压缩包共10个文件以Excel结果表格、Visio图表、MATLAB脚本和Word说明文档为主整体312KB内容包括热备用系数0.05与0.2两种场景下的机组组合求解结果、直流潮流下节点导纳矩阵及其逆矩阵、核心优化程序与基本要求说明图表和数据相互对应便于完整复盘建模与求解流程。已有304人学习下载。借助该包可理解机组组合的数学建模思路运行MATLAB脚本复现优化过程对比不同热备用水平对机组出力的影响数据与代码衔接完整主程序可直接运行并输出结果表格适合用于课程设计、毕业设计或实际调度分析。1. 电力系统机组组合优化是什么一份zip包背后的调度难题这份编号为034的zip包里装的是什么其实不用解压也能猜到八成——电力系统机组组合优化Unit Commitment简称UC。它解决的是每个调度员每天都要做两次的问题未来24小时哪些机组该开、哪些该关、每台发多少电才能在满足负荷和备用要求的前提下把总成本压到最低。别小看这个“开开关关”写成混合整数规划后规模可以从十几台机组膨胀到上千台求解时间从秒级变成小时级。这篇笔记适合三类人电力系统专业学生想找一份能跑通、能改参数的学习案例售电公司或园区微网的运行工程师打算把调度方案从拍脑袋换成优化求解刚入优化工具、想拿真实业务练手的算法工程师。我会从建模、代码、调参一路讲到踩坑把这类zip包背后最常用的落地路径一次说透。2. 机组组合建模与MILP选型从经济调度到启停决策2.1 目标函数里最容易算漏的三项成本机组组合和经济调度最大的区别是多了“启停”这个整数决策。经济调度默认机组已经在线上只分配出力UC要把机组开起来、关掉还得决定什么时刻操作。所以目标函数至少包含三块燃料成本、启动成本和停机成本。燃料成本通常用出力的二次函数拟合工程上为了保持MILP线性结构会做成三段分段线性函数。启动成本一般细分为冷启动、温启动、热启动三种价格相差很大取决于停机时长。停机成本虽然比启动成本小但一些大机组关停时要吹扫、要消耗辅助燃料算总账时漏掉它结果就偏向“频繁启停”的假便宜方案。很多入门案例只写燃料成本求解器为了省启动费就会疯狂启停机组最后给出的总成本很低但调度员一看启动次数就知道不能执行。这也是新手对不上实际成本最常见的原因——目标函数不完整优化出来的方案天然偏门。2.2 五类约束功率平衡、备用、爬坡、最小启停、初始状态UC的标准约束可以归成五组。第一功率平衡约束每个时段所有在线机组出力之和必须等于负荷这是硬等式。第二备用约束系统要预留旋转备用一般写成在线机组最大出力之和大于等于“负荷备用容量”。第三爬坡约束相邻时段之间出力变化不能超过爬坡率。第四最小运行和最小停机时间约束开机后至少要维持若干小时停机后也要冷却若干小时不能频繁启停。第五初始状态约束调度起始时刻之前各机组是开着还是关着、已连续运行多久停机多久都要作为参数传进去。前四类好理解初始状态是新手最容易爆雷的地方。比如某台机组在调度起跑点之前已经连续运行了20小时最小运行时间是8小时那它在第一个时段完全不用额外限制但如果这台机组前一天刚停机2小时而最小停机时间是8小时那么从调度开始的前6个小时里它都不能重新开机。如果没有这些初始信息模型会把第一时段当成“起点”算出违反物理约束的开机计划。这类错误在结果里很难排查因为总成本看着正常拉出来做时序回放才发现前面的机组状态全是错的。2.3 为什么选MILP启发式算法的两个硬伤很多人问“这类问题用遗传算法不也行吗”但电力系统行业从二十年前就开始转向MILP到现在商业求解器和开源CBC、SCIP都能处理几十万变量的问题原因有两点。第一启发式算法给不出最优性差距跑完一组结果你只知道“还行”不知道离最优解差百分之几MILP求解器明确告诉你gap调度员可以直接拍板“gap 1%以内就用”。第二UC带大量等式约束和逻辑约束启发式算法处理约束只能靠罚函数罚系数调起来是玄学一不小心功率平衡就根本不满足。MILP把约束写进模型可行域由求解器保证解出来就能直接执行。以pulp为例目标函数写起来很直观import pulp as lp # prob是问题对象u[i,t]是0/1开机状态p[i,t]是连续出力 # start[i,t]、shutdown[i,t]分别表示t时段是否发生启动和停机 prob lp.lpSum( fuel_cost[i] * p[i, t] start_cost[i] * start[i, t] shutdown_cost[i] * shutdown[i, t] for i in units for t in periods )这段代码有三个关键点需要说明。第一start[i,t]和shutdown[i,t]是二进制变量它们不是独立自由变量必须通过状态转移约束与u[i,t]关联否则求解器会“白嫖”启动成本——把每个时段都定义成启动但实际机组状态根本没变化。第二燃料成本系数这里用的是线性系数如果原始数据是二次曲线要提前做分段线性化。第三shutdown_cost容易漏我见过好几个模型只有启动成本没有停机成本优化结果里机组启停次数明显偏高。为什么强调MILP而不是动态规划动态规划做小规模UC很精确但状态爆炸太快机组超过50台基本没法用拉格朗日松弛可以解决更大规模但对偶间隙恢复可行解的代码复杂度高一般团队很难短时间搞定。因此从标准做法看MILP是成熟的工业选择后续换商业求解器也只需改一行参数。3. 把UC跑起来用python和开源求解器求解最小算例3.1 准备数据四台机组参数表与负荷曲线讨论理论够了先做一个能跑出结果的最小算例。四台机组24个时段负荷曲线直接写在脚本里。机组参数如下表燃料成本按一次线性近似单位是元/MWh方便计算。机组Pmin (MW)Pmax (MW)爬坡率 (MW/h)最小运行 (h)最小停机 (h)启动成本 (元)燃料成本 (元/MWh)G110040012042450060G28030010032300075G3502008022180095G4301005011800120为了简单不单独列停机成本在启动成本里按50%折算。负荷曲线用一段Python列表从凌晨低谷的350MW升到晚间高峰的850MW模拟典型的双峰形状。3.2 用pulp建模机组组合核心代码与逻辑说明先定义决策变量u[i,t]表示机组i在t时段是否在线p[i,t]是出力start[i,t]和shutdown[i,t]分别记录启动和停机动作。pulp里二进制变量直接指定catBinary。特别注意p[i,t]不能只设上下界还要与u[i,t]耦合机组离线时出力必须为0。import pulp as lp T 24 units { 1: {pmin: 100, pmax: 400, ramp: 120, minup: 4, mindown: 2, sc: 4500, fc: 60}, 2: {pmin: 80, pmax: 300, ramp: 100, minup: 3, mindown: 2, sc: 3000, fc: 75}, 3: {pmin: 50, pmax: 200, ramp: 80, minup: 2, mindown: 2, sc: 1800, fc: 95}, 4: {pmin: 30, pmax: 100, ramp: 50, minup: 1, mindown: 1, sc: 800, fc: 120}, } load [ 350, 340, 330, 350, 400, 500, 650, 750, 780, 790, 780, 730, 700, 680, 670, 690, 730, 760, 800, 850, 820, 700, 520, 380, ] prob lp.LpProblem(UC_Demo, lp.LpMinimize) u {}, p {}, start {}, shutdown {} for i in units: for t in range(T): u[i, t] lp.LpVariable(fu_{i}_{t}, catBinary) p[i, t] lp.LpVariable(fp_{i}_{t}, lowBoundunits[i][pmin], upBoundunits[i][pmax]) start[i, t] lp.LpVariable(fs_{i}_{t}, catBinary) shutdown[i, t] lp.LpVariable(fd_{i}_{t}, catBinary) # 目标函数燃料成本 启动成本 0.5*启动成本作为停机成本近似 prob lp.lpSum( units[i][fc] * p[i, t] 1.5 * units[i][sc] * start[i, t] for i in units for t in range(T) ) # 功率平衡 for t in range(T): prob lp.lpSum(p[i, t] for i in units) load[t] # 启停与出力耦合离线时出力为0在线时出力在Pmin~Pmax之间 for i in units: for t in range(T): prob p[i, t] units[i][pmax] * u[i, t] prob p[i, t] units[i][pmin] * u[i, t] # 状态转移定义启动和停机动作 for i in units: for t in range(1, T): prob start[i, t] u[i, t] - u[i, t-1] prob shutdown[i, t] u[i, t-1] - u[i, t] # 最小运行时间一旦启动至少连续运行minup个时段 for i in units: mu units[i][minup] for t in range(T - mu 1): prob lp.lpSum(u[i, tt] for tt in range(t, t mu)) mu * start[i, t] # 最小停机时间一旦停机至少连续停机mindown个时段 for i in units: md units[i][mindown] for t in range(T - md 1): prob lp.lpSum(1 - u[i, tt] for tt in range(t, t md)) md * shutdown[i, t] # 爬坡约束 for i in units: for t in range(1, T): prob p[i, t] - p[i, t-1] units[i][ramp] prob p[i, t-1] - p[i, t] units[i][ramp] status prob.solve(lp.PULP_CBC_CMD(msg1, timeLimit30)) print(求解状态:, lp.LpStatus[status]) print(总成本:, lp.value(prob.objective))核心逻辑说明三点。一是逆变电器耦合约束p[i,t] pmax * u[i,t]保证停机时机组出力被锁到0p[i,t] pmin * u[i,t]保证在线机组不能低于技术出力下限。少了这条求解器可能让一台离线机组按Pmin出力功率平衡虽然满足但状态和出力自相矛盾。二是最小运行时间的标准写法累计开机时长必须大于等于minup * start[i,t]这是一个经典的MILP建模技巧。三是启动和停机动作与状态的关系上面的写法是单侧松约束配合目标函数最小化求解器不会故意多报启动动作。跑这个算例时如果CBC在30秒内没达到gap 0通常也能返回一个可行解。四台机组规模很小一般几秒内就能找到最优解。3.3 求解参数设置与结果解读把timeLimit设成30秒是为了演示参数入口。真正跑的时候四台机组即使不用时间限制也很快但等到后面做几十台机组时这行参数就是救命稻草。结果解读建议额外打印每个时段的机组组合和出力明细for t in range(T): gen {i: lp.value(p[i, t]) for i in units if lp.value(u[i, t]) 0.5} total sum(gen.values()) print(ft{t1:2d} 负荷{load[t]:4.0f} 出力{total:4.0f} f在线{list(gen.keys())})判断结果好坏先看两件事第一每个时段出力总和必须等于负荷差值超过0.001就说明模型约束写错了第二启动次数是否合理比如四台机组24小时启动总数一般在3~8次如果出现十几二十次就要回头检查目标函数里是否漏了启动成本。4. 从玩具算例到真调度求解性能与参数调优4.1 MIP gap不是越小越好给求解器留余地很多人在MILP里只会设gapRel0.0001然后等几个小时。实际调度场景根本没必要。负荷预测本身有±3%的误差即使给你一个gap为0.001%的完美解执行时也会被预测误差吃掉。行业内做日前计划gap在0.5%到1%就足够日内滚动计划可以放宽到2%甚至5%。gapRel参数本质是“相对最优性间隙”公式是(best_bound - best_solution) / best_bound。设置1%不是“允许误差1%”而是告诉求解器如果当前解距离最优界在1%以内立刻停止。这一点新手上手时容易误读以为设了1%就等于1%的误差不可控。事实上求解器给出的解是可行解约束全部满足只是不保证继续改进一定超过1%。如果想看求解器当前最优界可以用solver lp.PULP_CBC_CMD(msg0, timeLimit300, options[gapRel 0.005]) prob.solve(solver) # bestBound在pulp中通常存在但某些版本可能取不到 print(目标值:, lp.value(prob.objective))这是对时间和精度的主动取舍。做滚动调度时每轮只有几分钟甚至几十秒窗口与其让求解器跑到gap0再出结果不如先拿一个1%以内的可行解把因预测更新产生的收益吃回来。4.2 时间限制与MIP start热启动实际业务中UC很少只跑一次。典型做法是每15分钟滚动一次用最新负荷预测重算未来4小时。如果每次都从零开始求解时间会累加得很可观。这时候热启动warm start也称MIP start就是最大的免费午餐。热启动的本质就是把你认为“可能比较优”的机组启停状态事先告诉求解器让分支定界从一开始就贴着这个解走剪枝效率高很多。在Gurobi里可以这样传递上一轮的最优解from gurobipy import Model, GRB m Model(UC) # 这里假设u_vars已经是上一轮求出的0/1矩阵 for i in unit_list: for t in T_list: u_var[i, t].start last_solution[i, t] # start是Gurobi为变量预留的MIP start属性 m.optimize()上一轮解出来的机组状态和这一轮差别通常不大热启动能让求解时间缩短30%~70%。但要注意如果负荷突变特别大比如某台机组故障导致负荷转移热启动解可能不在可行域内求解器会自动剔除不可行部分继续从头探索。所以热启动是“加速器”不是“保证书”。开源求解器CBC也支持MIP start但pulp里传入的API在不同版本间不一致使用前要先验证当前版本是否支持如果不支持可以放弃热启动改用“上一轮解作为初始值”手动构造可行解但意义不大。因此若项目规模上千台机组商业求解器的热启动价值会非常明显。4.3 大规模场景的分解与聚合思路机组数量上千时直接求解MILP的变量数轻松破百万时间和内存都受不了。业界常见做法有三个机组聚合、时间粒度和拉格朗日松弛。机组聚合是把同类型、同容量区间的机组合并成一台等效机组。比如某风电场或光伏电站有50台同型号逆变器聚合成一台整体参与优化再按等比例或边际成本排序分配出力。对于火电同容量机组也有效。代码如下def aggregate_units(unit_params, key_fn): groups {} for uid, param in unit_params.items(): key key_fn(param) groups.setdefault(key, []).append(uid) agg {} for key, ids in groups.items(): agg[key] { pmax: sum(unit_params[i][pmax] for i in ids), pmin: sum(unit_params[i][pmin] for i in ids), # 成本系数按额定容量加权避免简单平均带来的偏差 fc: sum(unit_params[i][fc] * unit_params[i][pmax] for i in ids) / sum(unit_params[i][pmax] for i in ids), } return agg注意加权平均的目的是让聚合后的成本特性尽量接近原群体。如果简单把燃料成本相加聚合机组成本会虚高优化结果会把电量偏向真正的低成本机组。时间粒度分解则是先跑每2小时一个大时段的计划确定大致启停方案再把第一个小时到第四个小时细化成15分钟粒度重优化。这样既保留了可执行性又降低了问题规模。拉格朗日松弛把系统功率平衡约束乘子化逐台机组独立求解再更新乘子迭代工程上常用它来生成初始可行解再交给MILP做最终优化。5. 踩坑记录从zip解压到约束漏加的五个典型问题5.1 zip解压报invalid zip archive文件没下完整现象拿到“034电力系统机组组合优化.zip”后双击解压直接弹“文件损坏”或“invalid zip archive: could not find eocd”也有人说成“failed to copy spatial iop zip”这类奇奇怪怪的报错。原因zip文件尾部有一个End of Central Directory记录Linux解压工具找不到这条记录就是典型的“文件没下完整”。可能是下载中断、网络缓存残留、或者从网盘拿到的分包没有全部放在同一目录就急着解压。解决先核对文件大小跟下载页标注的字节数一致再解压。第二步用命令行测试unzip -t 034电力系统机组组合优化.zip如果输出里有“End-of-central-directory signature not found”字样别犹豫重新下载。若文件特别大且下载工具支持断点续传直接续传后再次测试如果是分包压缩务必把分卷文件都放在同一个文件夹否则即使扩展名正确也解不出来。还有一个常见情况普通zip但下载时被浏览器改成“下载.zip.zip”文件大小正常却解不了这时改回原名即可。5.2 zip内中文乱码编码不匹配现象解压后文件名全是乱码类似“鎶楃伀”或“缃戠粶”这样的字符。原因zip包在Windows系统压缩时文件名一般按GBK编码写入在Linux或macOS上用默认UTF-8解压中文名就变成了乱码。这跟zip包本身无关纯粹是编码不一致。解决Linux环境用unzip -O GBK重新解压如果用的压缩软件不支持指定编码可以先用Python把zip里的文件名修正出来再根据修正后的名字手动解压到本地。代码如下import zipfile with zipfile.ZipFile(034电力系统机组组合优化.zip) as zf: for info in zf.infolist(): # 先把字节按cp437解码再按gbk重新编码还原中文 safe_name info.filename.encode(cp437).decode(gbk, errorsignore) print(safe_name)这个技巧对含中文的zip非常通用。如果zip本身加了密码编码问题反而是小事——没有密码只能在脚本里跑字典属于合规性以外的话题不做展开。5.3 结果违反功率平衡单位与约束两处都查现象求解完成后打印各个时段出力之和发现和负荷差了几十甚至几百MW或者某台机组状态显示为停机但出力却等于Pmin。原因先从约束排查最常见的是漏写了p[i,t] pmax * u[i,t]和p[i,t] pmin * u[i,t]停机机组照样被分配到出力。其次查数据单位负荷如果用的是MWh而机组出力是MW一个24小时的时段里差了24倍。解决先加一个结果校验模块每次求解完成后自动对账for t in range(T): gen_sum sum(lp.value(p[i, t]) for i in units) if abs(gen_sum - load[t]) 1e-3: print(f时段{t1}功率不平衡: 负荷{load[t]:.2f} 出力{gen_sum:.2f})这段校验不是可选项是必选项。尤其是换数据源时经常遇到“MW当成MWh”“某机组Pmax写错”之类的低级问题有了它对账翻车概率大大降低。5.4 最小启停约束导致无解初始状态写错现象模型报infeasible但把最小运行时间和最小停机时间约束全部注释掉后立刻有解说明问题就出在这组约束上。原因最小启停约束是按“调度起点”写的如果机组在前一天已经运行或停机了若干小时你必须在模型里把这段“历史状态”接过来。我见过一个典型案例某机组最小运行时间8小时初始状态写成了“在线”但实际它只连续运行了2小时模型要求它前6小时必须保持运行可功率平衡又被其他机组占满导致无解。根本原因是初始状态写错或根本没初始化。解决在建模前把每台机组的初始运行/停机时间单独建表并在约束中用分段判断弥补剩余时长。一个常见写法是如果机组初始在线时长已超过最小运行时间则首个时段不额外施加minup约束如果还没达到则让它在剩余不足的时段里保持在线。最小停机时间同理初始停机时长不足时前几个时段禁止开机。# init_status: 1在线 0停机 # init_time: 已经连续运行或停机的小时数 if init_status[i] 1 and init_time[i] units[i][minup]: need units[i][minup] - init_time[i] for t in range(min(need, T)): prob u[i, t] 1这就是前文说的“初始状态约束”。很多教程里不写这层逻辑是因为假设机组已满足历史状态但真实调度场景里凌晨1点出现无解往往就是这个根因。5.5 求解器报infeasible先做约束松弛现象模型无解求解器只返回“problem infeasible”不告诉你哪条约束冲突排错全靠猜。原因MILP求解器内部用的是分支定界它发现可行域为空时会直接停止不会提供一条人话解释。爬坡约束跨时段耦合一处无解会传导到后续所有时段排查特别费劲。解决先用“连续化诊断”找到冲突范围。把机组状态变量从二进制放松成0到1的连续变量再求解同一模型。如果连续模型有解说明是整数约束和状态转移之间的逻辑冲突如果连续模型也无解说明线性约束之间就矛盾了。连续模型求解后查看哪些u[i,t]落在0到1之间这些时段和机组往往就是约束冲突的高发区域for i in units: for t in range(T): u[i, t] lp.LpVariable(fu_{i}_{t}, lowBound0, upBound1) # 由Binary改成Continuous只用于诊断 prob.solve() for i in units: for t in range(T): val lp.value(u[i, t]) if 0.01 val 0.99: print(f诊断: 机组{i} 时段{t1} 状态值{val:.3f})中间值意味着这个时段的状态被两组约束拉扯重点检查该机组的爬坡、最小启停和备用约束。这种“连续松弛诊断法”在大型模型里比人肉对数小时管用得多是我现在排查无解的第一动作。6. 进阶考虑不确定性与回测验证的实用技巧6.1 从确定性UC到随机UC场景法怎么改确定性UC里负荷是单条曲线随机UC常见做法是场景法构造多个负荷预测场景每个场景设一个概率。状态变量u[i,t]是所有场景共用的非预期变量而出力p[i,t,s]允许随场景变化。目标函数改为所有场景成本按概率加权。代码层面的改动很小只需要给p和load加一个场景下标再把功率平衡约束对每个场景分别写。代价是变量规模乘以场景数所以工程上场景常用概率距离方法剪枝到5~20个。6.2 用历史负荷做回测验证结果是否可信不管模型写得多复杂最终要回答“结果能不能用于明天调度”。我的习惯是拿过去30天的历史负荷做回测每天用前一天的预测负荷跑SCUC再把当天实际负荷代入结果统计功率平衡偏差和爬坡越限次数。验证脚本核心就一段def validate(plan_result, actual_load, tolerance10.0): imbalance (actual_load - plan_result[total_gen]).abs() violations (imbalance tolerance).sum() return violations # 返回越限时段数越少说明方案越稳能通过回测的方案才敢排进调度计划否则模型再精巧也只是纸面好看。这类优化项目里我吃过的最大教训就是“迷信gap和最优性”却忽略了约束本身是否正确。有一版模型把爬坡约束方向写反求解器跑出完美的0.1% gap拿去现场用才发现无法执行。从那以后任何UC代码跑完结果我都强制自己先看启动次数、单位时段启停频繁度再做平衡校验。如果你手头正打开这份zip包希望这篇笔记能帮你少走同样的弯路希望帮到你。本文还有配套的精品资源点击获取