做综合能源系统调度这两年最让我头疼的不是设备模型多复杂而是系统里电、气、热设备根本不属于同一个运营方。光伏可能在用户侧屋顶燃气机组归能源公司蓄热罐又是热力站的资产。如果硬按单一运营商做全局优化物理上能算出最低成本但现实中根本执行不下去——因为没人愿意无条件牺牲自己的利益去成全系统。后来我把思路切到合作博弈框架用 Matlab 写了一套“联盟优化 利益分配 分配结果再反馈到调度约束”的完整实现才把账算明白。这篇就聊聊这套实现里最容易被忽略、也最能决定论文和工程方案能不能站住脚的细节主要面向正在用 Matlab 做综合能源系统优化调度、想把合作博弈落下去的同学们。1. 为什么综合能源系统的调度要引入合作博弈1.1 同一张网多个账本传统调度模型通常假设一个“上帝视角”的调度中心所有设备、负荷、储能在同一面旗下运行目标函数是总成本最小或总收益最大。但在实际项目里多能互补系统经常是多个主体各自投资、各自运营的电网侧关心峰谷套利热力公司关心供暖保障燃气公司关心气价波动用户侧聚合商关心可调负荷的补偿。这种场景下你按总成本最小去调度就可能出现某个主体的设备被压到极低出力、收益下降的情况协调自然失败。合作博弈解决的正是“多个独立决策主体如何行动能整体最好同时每个人又愿意参与”的问题。它把各主体看作博弈玩家先把大家绑定成一个联盟去优化联盟总收益再把联盟产生的合作收益按某种兼顾公平和稳定的方式分回去。引入合作博弈后调度问题变成一个“先算锅里的肉再分锅里的肉”的过程。肉不够大方案没有吸引力肉分不匀方案没法落地。1.2 合作节省的钱从哪来用一句话概括合作收益来自各主体单独运行时无法利用的互补性。举个我在某区域项目里简化的例子。光伏富余时段用户侧用不完的电如果只馈到电网收益很低但隔壁有台电锅炉可以吃电产热给热力管网这就是互补。燃气机组余热可以进余热锅炉供给热负荷减少电锅炉耗电也是互补。储能、蓄热罐还能把时段间的价差套出来。这些互补性意味着如果每个主体各自优化后再简单并网交易会漏掉很大一部分协调收益。联盟优化则把所有资源共享到一个统一优化问题里让燃气轮机出力、电锅炉功率、储能充放策略、热罐蓄放热策略同时决策从而拿到更低的总运行成本。系统单主体优化得出的最低成本和合作联盟优化得出的最低成本差额就是合作博弈要分的合作节省额。1.3 调度与分配的耦合才是难点有一点必须提醒很多文章把利益分配当成“事后算账”跑完联盟优化用 Shapley 值分一下账就结束。但在工程上分配结果会影响下一步的调度决策——如果某个主体分到的收益低于他单独经营的收益他会退出联盟然后联盟结构变了调度方案也得跟着变。所以严格的模型应该是“调度—分配—校验参与约束—必要时重新调度”的迭代闭环。我在模拟项目 X 里就是用双层结构实现的内层做各主体可形成的联盟运行优化外层做收益分配与稳定性校验如果校验不通过就把参与约束以惩罚项的形式加回优化目标再求解循环若干轮直到所有主体都愿意留在联盟内。2. 从物理模型到特征函数IES合作博弈建模的关键细节2.1 能源集线器把多能耦合写进约束综合能源系统的物理模型我习惯用能源集线器Energy Hub形式描述。它把输入端电网购电、天然气站购气、光伏、风电和输出端电负荷、热负荷、冷负荷通过设备转换矩阵连接起来。建模时我通常把每个参与主体视为拥有若干设备的集合设备类型主要有这几类燃气轮机天然气转电力附带热电比约束余热锅炉吸收燃气轮机余热生产热力电锅炉电力转热效率一般取 0.9 到 0.98热泵用电从低温热源提热COP 通常 2.5 到 4蓄电/蓄热储能类需要 SOC 递推和容量约束光伏/风电出力曲线给定可弃风弃光电网购电、天然气购气和设备出力构成决策变量约束至少包括电功率平衡、热功率平衡、设备出力上下限、爬坡约束、储能 SOC 连续性。时间粒度我常用 1 小时一天 24 个时段总量不算大如果你做高比例可再生场景可以细化到 15 分钟但求解时间会明显上涨。2.2 特征函数怎么算每个联盟都要重新优化合作博弈里最关键的概念是特征函数 v(S)表示联盟 S 合作能带来的额外收益。在我这套模型里联盟收益统一定义为“独立运行总成本 - 联盟最优运行成本”也就是成本节约额。因此必须先对每个非空联盟 S 都做一次优化调度得到对应最小总成本 C(S)再算 v(S)ΣC({i}) - C(S)。这里有个特别容易踩的坑计算 C({i}) 的口径必须和 C(S) 完全一致。单独运行时主体 i 的设备只有它自己的购电、购气按它自己的负荷联盟运行时原属不同主体的设备可以协同但固定成本、惩罚项、违约惩罚这些非运行成本项目不要左边算进去右边没算。口径不一致会导致分配结果出现负数或超额后面校验根本过不了。所有非空联盟的数量是 2^N-1。N3 时只有 7 个子问题Matlab 跑起来很轻松N4 时 15 个N5 时 31 个如果每个子问题都是 MILP单核可能要跑几十分钟这时候就该考虑并行或近似方法具体在第 5 节展开。2.3 三主体算例手工验证特征函数与 Shapley 值为了确认程序没有算错我在搭代码之前先手推了一个 3 主体简化案例。假设三个主体独立运行成本分别是 c112、c210、c38联盟后运行成本及成本节约额如下联盟最优运行成本成本节约额 v(S){1}120{2}100{3}80{1,2}202{1,3}173{2,3}153{1,2,3}255这个表翻译成程序逻辑就是对每个联盟独立调一次度取成本求差。程序跑出来的 v 如果跟手算不一致一定先查是不是联盟约束写漏了。3. Shapley值、核仁法与最小核心分配方法怎么选怎么用3.1 Shapley值公平性最强的默认选项Shapley 值的基本思想是一个主体得到的分配等于它所有可能加入顺序中对联盟边际贡献的平均值。拿上面的例子手算主体 1 的 Shapley 值为 1.5主体 2 为 1.5主体 3 为 2合计 5正好等于大联盟的总节约额。计算上每个主体 i 的 Shapley 值公式是φ_i Σ_{S⊆N{i}} |S|!(n-|S|-1)! / n! × [v(S∪{i}) - v(S)]代码实现时直接枚举不含 i 的联盟子集乘上权重即可。当 N 不超过 5 时穷举完全没问题。Shapley 值有很好的性质总和等于大联盟总收益、对称、哑元为 0、可加性。工程上大家默认它最公平不容易引起争议。但 Shapley 值没有保证分配结果一定落在核心内也就是说可能出现某个子联盟分到的总收益少于它独立合作能拿到的收益这时子联盟成员有动机退出。我在项目里遇到过某个含储热主体的小联盟分到的收益偏低而它单独跟热力公司抱团能拿更多于是靠 Shapley 分完的联盟直接不稳定。3.2 核仁与最小核心稳定性优先的替代方案如果稳定性是第一诉求可以用核仁法或最小核心法。核仁的思路是寻找一个分配使“最不满意”的联盟的不满程度尽量小。不满程度通常用超额excess衡量e(S)v(S)-x(S)其中 x(S) 是联盟 S 所有成员分到的收益之和。核仁算法要逐步最小化排序后的超额向量本质是一个系列线性规划。最小核心法是更简单的松弛版本引入一个变量 ε约束所有联盟的 x(S)ε≥v(S)然后最小化 ε。它的含义是允许所有联盟的整体分配都往公平性上稍微偏离一点但保证没有一个联盟会退出太多。好处是线性规划即可求解坏处是解不唯一两个不同的解都是最小核心点所以有时还需要再叠一个利润均等的目标来选点。我之前做过三种方法的数值对比见下表指标Shapley值核仁法最小核心法唯一性唯一唯一不唯一分配公平性最强一般较弱联盟稳定性可能不稳定强制稳定近似稳定计算复杂度O(n·2^n)多次LP迭代一次LP适用场景论文基准、多方协商强调联盟牢固大规模近似3.3 稳定性校验分配结果能不能让所有人留下不管用哪种方法我最终都会做一项校验对每个非空联盟 S检查 Σx_i ≥ v(S)。如果某个联盟不满足就意味着那批成员合起来没拿到他们独自合作能拿到的收益联盟有解体风险。我的处理方法是先尝试调整分配规则比如改用核仁如果还是不行就说明特征函数本身有问题需要回查 v(S) 的计算口径。校验通过之后再把分配结果回填到调度模型里。具体做法是把各主体的参与约束自己的分配收益不低于单独经营收益写成软约束加进联盟优化目标函数里稳定解作为下一步调度的边界条件。4. Matlab实现框架与核心代码拆解4.1 总框架枚举联盟、求解特征值、计算分配我用 Matlab 实现时把整个流程拆成五步每步一个函数方便单独调试数据准备读入负荷曲线、设备参数、分时电价、气价联盟枚举生成 2^N-1 个联盟的成员组合特征值求解对每个联盟调用统一优化函数得到 C(S)保存到向量分配计算根据特征值向量分别算 Shapley 值、核仁或最小核心结果后处理输出各主体分配收益、运行成本、稳定性检查表下面是联盟枚举与特征值计算的骨架N 3; % 参与主体数量 numCoal 2^N - 1; c zeros(numCoal, 1); v zeros(numCoal, 1); for idx 1:numCoal members find(bitget(idx, 1:N)); c(idx) solve_union_cost(members, data, params); end cSingle zeros(N, 1); for i 1:N cSingle(i) c(find(bitget(1, i) (1:numCoal))); % 简化示意 end totalStandalone sum(cSingle); for idx 1:numCoal members find(bitget(idx, 1:N)); singleCost sum(cSingle(members)); v(idx) singleCost - c(idx); end实际工程里我不会用find反复匹配单主体索引而是直接把单主体成本单独存一份避免逻辑混乱。核心是solve_union_cost这个函数它接受成员编号数组内部构建对应设备的优化模型并调用求解器。4.2 统一优化函数YALMIP 建模与求解器配置我习惯用 YALMIP 做建模界面底层求解器用 Gurobi。如果你没有 Gurobi 授权也可以用 Matlab 自带的 intlinprog 处理 MILP线性规化问题用 linprog。快速安装配置参考如下% 求解器路径与测试以Gurobi为例 addpath(D:/gurobi/win64/matlab); gurobi_setup(); currentSolver yalmp;模拟项目 X 中单联盟 24 时段 MILP 模型大概有几百个变量、上千个约束Gurobi 默认参数就能一两分钟内收敛。YALMIP 里核心建模片段大致是这样x sdpvar(repmat(nT, 1, nGen), repmat(1, 1, nGen)); % 简写示意 F [x(1) lower(1), x(1) upper(1)]; % 电平衡约束 F [F, sum(pGen) buyPower pv - loadElec chargeP - dischargeP]; % 热平衡约束 F [F, sum(hGen) boilerHeat hpHeat - loadHeat heatIn - heatOut]; % SOC递推 for t 2:nT F [F, soc(t) soc(t-1) etaC * chargeP(t) - dischargeP(t) / etaD]; end optimize(F, sum(operationCost), options);写约束时一定注意电、热、气单位的统一。热负荷和电负荷经常一个是 GJ一个是 MWh燃气轮机热电比决定余热量如果单位换错特征值会差出好几倍。我吃过这个亏后来所有热功率统一折算成 kW燃气的单位热值也换算成以电等价的 kW 口径问题才消失。4.3 Shapley 值和核仁法的Matlab实现Shapley 值实现按公式穷举即可phi zeros(1, N); for i 1:N for S 0:(2^N - 1) if bitget(S, i) 0 % S不包含i s nnz(bitget(S, 1:N)); idxS S 1; idxSUi bitor(S, bitshift(1, i-1)) 1; if idxS 1 || idxS numCoal || idxSUi 1 || idxSUi numCoal continue; end phi(i) phi(i) factorial(s) * factorial(N - s - 1) / factorial(N) ... * (v(idxSUi) - v(idxS)); end end end核仁的 LP 可以用 YALMIP 快速搭出来核心是epsilon和分配向量x的优化x sdpvar(1, N); epsilon sdpvar(1, 1); F [sum(x) v(end), x 0]; for idx 1:numCoal members find(bitget(idx, 1:N)); F [F, sum(x(members)) epsilon v(idx)]; end optimize(F, epsilon, options);这个模型求解非常快几行代码就能替代多轮 LP 的繁琐实现。如果你要标准核仁还需要对超额向量做词典序最小化我一般用迭代加约束的方式实现代码会稍复杂一点但论文里通常只需给出核心法结果即可。4.4 并行计算和结果缓存特征值求解是整个流程最耗时的部分。N 达到 4 或 5 时我会用 Parallel Computing Toolbox 的 parfor 替代 forparfor idx 1:numCoal c(idx) solve_union_cost(find(bitget(idx, 1:N)), data, params); end还要把算过的特征值存成 mat 文件。因为做参数敏感性分析时电价或气价一变某些联盟的最优成本变化很小缓存后能省掉大量重复求解。这是我在处理某高比例可再生场景时总结出来的高效经验看似不起眼实际能节省一个数量级的时间。5. 从集中式到分布式大规模场景的可扩展路径5.1 为什么必须考虑分布式化当 N 超过 6穷举 2^N-1 个联盟就变得不现实更麻烦的是集成系统内部各主体未必愿意把自己的设备参数和负荷数据全部交给一个中心求解器。这种情况下需要引入分布式优化框架。我在扩展模拟项目 X 时采用了 ADMM 思路每个主体内部保留自己的设备模型和隐私数据只通过公共耦合变量和价格信号跟其他主体交换信息。ADMM 的耦合变量通常是联络线功率、热网交换热量、储能共享电量。每个主体求解自己的子问题然后更新拉格朗日惩罚项循环迭代。这个框架的优点是能保留 Shapley 值计算所需的各联盟成本估计缺点是收敛速度受惩罚参数影响很大需要手动调参。我的建议是先用集中式结果作为基线验证分布式优化得到的系统总成本和集中式相差不大再继续后面的分配环节。5.2 Shapley 值的近似计算大规模下精确 Shapley 值必须近似。我试过两种比较有效一是基于采样的边际贡献估计。随机抽取大量排列统计主体 i 在每个排列中的边际贡献取平均值。排列数量从几千到几万看你对误差的容忍度。二是基于多线性扩展的近似公式。Shapley 值本质上是值函数在某个点上的离散梯度平均当 N 很大但设备规模相似时可以用连续化后的解析近似。更简单但不那么严谨的做法是直接用比例分摊或最小核心做粗分配然后检验稳定性。我在某些快速方案里会用最小核心替代 Shapley作为初步筛选如果最小的几个联盟都没有退出动机就不再跑采样近似了。5.3 鲁棒和不确定性场景的扩展方向合作收益很大一部分来自可调节设备对价格波动的响应。如果光伏、风电出力不确定后面还可以把模型扩展成两阶段鲁棒或随机规划让特征函数变成期望值或最坏情况下的收益。比如天然气价格上浮 10% 时某个联盟的收益可能由正转负这种情况在分配时必须考虑进去。我在项目中把分时电价和负荷预测的不确定性用若干典型场景替代发现各联盟收益排序会变化Shapley 值也跟着变所以论文里做参数敏感性时一定要写明“分配结果基于预测场景集”。6. 调试与验证几个最常见的翻车现场和排查方法6.1 手工算例先行我在 2.3 节用了三主体算例。这个算例小到什么程度小到可以在纸上验算。程序跑完后先检查 v 是否满足超可加性v({1,2}) 2v({3}) 0v({1,2,3}) 5满足再看 v({1,3}) v({2}) 3也满足。如果发现超可加性不满足比如 v({1,2})v({3}) v({1,2,3})那说明联盟模型里存在“负耦合”比如某个设备重复投资、容量互相挤占。这通常不是算法问题而是物理模型建模问题。6.2 分配结果出现负数或超额不合理我遇到过用 Shapley 值分完某个主体分到的收益为负。查到最后原因很简单算单主体成本时把固定电价的容量费也算了算联盟成本时却忘了加总固定费用。所以我的建议是跑任何验证前先检查 C({i}) 和 C(S) 是否在同一口径下。把所有非运行成本统一抽出来单独做成常数项而不是揉进每个子问题的目标函数否则很容易重复计算。6.3 MILP 求解时间爆炸如果联盟优化模型里启停变量太多24 时段也会卡住。第一次尝试可以把启停约束改成连续变量用线性松弛判断成本差异这样能快速定位是不是整数变量导致的复杂度。另一种办法是给 Gurobi 提供启发式初始解先按不含启停的 LP 求解然后用结果去设定启停变量的初始值能大幅缩短分支定界时间。6.4 YALMIP 与求解器接口报错这类问题大多出在求解器路径没有配置好或者版本不匹配。报错信息只要出现“Solver not found”或“No appropriate solver”先检查 addpath 是否真的生效再用yalmiptest跑一下诊断。我在某次升级后 Gurobi 版本变了YALMIP 提示找不到求解器重装接口之后恢复。这里不建议在模型里混用多个求解器因为不同求解器的数值容差不一样会对特征值微小的差值产生影响。6.5 单位不统一导致特征值巨变热负荷和电负荷单位不同、燃气的热值与电力换算系数不一致是隐藏最深的错误。我在某次计算中把热负荷误用成 MWh实际题目给的是 GJ导致单主体成本和联盟成本同时放大虽然 v 看起来还算正常但分配比例失真。遇到数值异常先检查单位矩阵再检查效率乘以输入功率的方向是否一致。最后分享几个我自己的实操习惯纸质手算三主体算例永远比先跑 N5 的大算例高效得多。每次改完模型我都先把 N 临时设成 3跑一组手工可验证的数据通过后再放大规模。分配方法我一般同时算 Shapley 值和最小核心论文里两张表放一起对比工程汇报里只给 Shapley 值加稳定性校验结果。如果你正要把这套 Matlab 代码写进研究或项目报告我的建议是不要把分配结果当成最终输出一定要把它作为调度方案的一个校验环节另外把特征值求解结果做成 Excel 或 CSV 表导出方便回溯哪一次求解用的哪一套参数。这套流程我用下来最大的体会是它最磨人的不是优化算法本身而是把合作博弈的数学定义翻译成准确、口径一致的运行优化模型。一旦这一步做扎实后面的 Shapley、核仁、分布式扩展就都有了稳定基础。