前阵子帮某园区做多能互补方案时遇到一件挺有意思的事几个供能主体原本各干各的分开算成本加上电网购电费、燃气费、设备折旧一个月总支出大概是10.5万。后来我们把热电联产机组、燃气锅炉和储能拉到一起联合调度总成本直接降到9万。账面上多了1.5万“合作红利”但问题也来了——这笔钱怎么分按投资比例分储能方觉得吃亏按用电量分热力方不答应。谈了两轮差点回到各自为政的状态。这个场景就是综合能源系统里典型的利益分配问题。做优化调度的人都知道综合能源系统里的电、热、冷、气不是各走各的管道而是通过热电联产、电锅炉、热泵这些设备深度耦合。合作运行能带来实实在在的红利但“蛋糕做大以后怎么切”才是项目落地的真正瓶颈。这篇文章我想把合作博弈在这类问题里的建模思路、分配算法、Matlab实现和实战中容易踩的坑完整梳理一遍适合正在做综合能源优化调度、被“怎么分钱”卡住的研究生和工程师参考。1. 综合能源系统的“集体红利”与分配难题合作博弈为什么在场1.1 多主体协同优化的成本与收益综合能源系统之所以“综合”核心在于设备之间的互补性。一个带热电联产机组的系统发电的同时产生余热这本是天然的资源但如果热力系统归A主体管、电力系统归B主体管余热就只能在A自己的管网里消化消化不了就白白散掉。B这边可能还在用燃气锅炉专门烧热等于一边浪费热、一边又花钱买热。把多个主体拉进同一个优化框架后情况完全不同。用电低谷时段燃气锅炉的出力可以压下来由热电联产机组的余热接管供热电网电价高企时储能在低谷充的电可以放出来顶上电网电价低时又可以用电锅炉替代部分燃气供应。这些操作每一项单独看都不大合在一起就是5%到30%的总成本下降空间。这片下降空间就是“集体红利”。从数学上看红利来自联盟成员的资源互补和负荷叠加效应。负荷曲线叠加后峰值削减、设备利用小时数提高、购能曲线更平缓这些都是单打独斗时拿不到的好处。1.2 没有公平分配合作随时可能破裂问题在于联盟总成本的下降不等于每个主体都能接受最后的分账结果。某主体可能只出了1台闲置锅炉却因为名义资产规模大而分走了大头真正承担调峰压力的储能方反而分得少。这种方案就算算得再精确也落不了地。合作博弈在这里的价值是把“公平”从概念变成可计算的量。它不回答“谁道德上应该多拿”而是回答“按什么规则计算每个人都愿意继续留在联盟里”。一套完整的利益分配方案必须满足两个基本条件个体理性即每个主体分到的收益不低于自己单独干时的收益集体理性即所有主体分到的收益加总等于联盟总收益。这两个条件保证联盟“散不了伙”而Shapley值、核仁法这些分配规则是在这两个条件之上进一步给出唯一且可解释的“切法”。2. 合作博弈建模特征函数与联盟成本的计算逻辑2.1 特征函数如何定义联盟的价值要把合作博弈套到综合能源系统上第一步是定义特征函数v(S)。这里的S是任意一个参与人子集v(S)表示该联盟通过内部优化调度能获得的价值。如果做成本型博弈v(S)可以定义为联盟S内所有设备协同运行时的最小总运行成本如果做收益型博弈则定义为联盟相对各成员单独运行时的成本节省总额。我自己的习惯是采用“成本节省”口径先算出每个成员单独运行的成本c({i})再算出任意联盟S的最低运行成本c(S)联盟价值就是 v(S) Σ_{i∈S} c({i}) - c(S)。这个口径的好处是直观分配时大家讨论的是“节省下来的钱”而不是抽象的“联盟总效用”。实际谈判中决策者听的是一张张账单这种口径容易对齐。这里有个容易忽略的细节v(S)不是简单地把各成员成本相加而是要重新优化计算。因为联盟里多台设备协同运行后最优运行策略会完全改变。比如单独运行时热电联产机组可能只在白天满发加入储能联盟后它可能改成夜间多发电、白天停机让储能放电整体策略都重排了。所以每个联盟S都要单独求解一次调度优化问题这个成本计算量会在后面成为主要瓶颈。2.2 能源枢纽模型与联盟成本计算联盟内部的设备耦合关系最常用的建模方式是能源枢纽模型。核心思想是把联盟看作一个输入输出端口输入侧购买电网电力、天然气输出侧满足电负荷和热负荷中间的设备用耦合矩阵串联起来。以一个小型三方联盟为例A主体拥有电储能B主体拥有热电联产机组C主体拥有燃气锅炉和热储能。联盟S优化调度的目标函数可以写成min Σ_{t1}^{T} [ c_e(t)·P_grid(t) c_g(t)·F_gas(t) ]约束条件包括电力平衡、热力平衡、设备出力上下限、爬坡约束、储能SOC递推方程等。热电联产机组的电出力P_chp和热出力H_chp之间满足热电比约束通常写成 H_chp(t) δ_c·P_chp(t)。这里的δ_c是热电比不同机组差异很大实际代码里得逐个设备标定。特征函数求解过程可以抽象成黑盒输入联盟成员集合S、设备参数、负荷曲线、能源价格输出该联盟的最优成本c(S)。这一步是把合作博弈和经典优化调度衔接起来的枢纽也是整个代码框架的核心函数。3. 分配方法选择Shapley值、核仁法与纳什谈判解的取舍3.1 Shapley值每个成员的边际贡献期望Shapley值是合作博弈里名气最大、学术论文里用得最多的分配规则。它的思想很简单一个成员的价值等于它在所有可能的“加入联盟顺序”中对联盟价值贡献的平均值。以三人联盟{N A,B,C}为例所有可能的加入顺序有6种。A的边际贡献在每种顺序下各不相同Shapley值就是这6个边际贡献的算术平均φ_A (1/6) × [ v(A) - v(∅) v(AB) - v(B) v(AC) - v(C) v(ABC) - v(BC) ... ]推广到n人Shapley值的计算公式为φ_i Σ_{S ⊆ N{i}} [ |S|! · (n - |S| - 1)! / n! ] · [ v(S∪{i}) - v(S) ]Shapley值是唯一满足对称性、有效性、可加性和哑元性质这四个公理的分配规则所以它被视为“最公平”的基准。但也正因为要遍历所有可能的联盟排列计算复杂度是指数级的联盟数超过10个之后就很难直接枚举。3.2 核仁法照顾最不满意的成员核仁算法走的不是“平均贡献”路线而是“稳定性优先”路线。对任何一个分配向量x可以计算每个联盟S的“超额”e(S,x)v(S)-x(S)也就是联盟成员对分到手的收益的不满程度。核仁法通过线性规划逐步压缩最大的超额值最终找到一组分配让最不满意的联盟也尽可能满意。从数学性质上说核仁法得到的分配一定在核core内只要核非空。核内的分配满足所有联盟的联盟理性对任意SΣ_{i∈S} x_i ≥ v(S)。换句话说任何一个联盟拿到手的总收益都不低于它单干的收益任何子联盟都没有理由脱离大联盟。核仁的优点是稳定性极强实际谈判中各方很难挑出毛病缺点是计算复杂且分配结果有时看起来“不公平”。在我前面提到的算例里Shapley值给贡献大的成员分了600核仁法却可能只给500因为核仁优先喂饱“最不满意的二人联盟”反而压制了强势成员的讨价还价能力。3.3 不同分配方法的适用场景我在实际项目里总结出一个经验先想清楚项目追求的是“说服力”还是“稳定性”。分配方法核心逻辑计算复杂度典型适用场景Shapley值按边际贡献平均指数级需要枚举联盟学术研究、多方博弈仲裁、追求公平时核仁法最小化最大超额中高需迭代线性规划联盟稳定性要求高、担心合作破裂时纳什谈判解最大化各方效用乘积中常需迭代求解双方谈判、收益弹性差异大的场景按投资比例分摊按资产占比分配低项目初期快速过账、各方出资额明确时如果项目只有两三个主体Shapley值计算起来很轻松直接用它就好。如果主体数量很多或者联盟内有明显的弱势联盟我会建议用核仁法兜底再用Shapley值做一个参考对照。两种方法结果差异不大时方案的真实可接受度通常很高。4. Matlab环境下的调度建模与求解流程4.1 代码框架数据、联盟成本函数、分配计算Matlab代码实现的核心不是“写公式”而是设计好调用关系。我一般把代码分成三个层次。第一层是数据层负责读入设备参数、负荷曲线、能源价格所有参数都打包进一个结构体params。第二层是联盟成本函数层输入一个联盟成员标识数组求解该联盟的最优调度成本。第三层是分配计算层枚举所有联盟、调用成本函数、计算Shapley值或核仁法结果。% 主程序示例计算所有联盟成本 n 3; % 参与主体数量 v zeros(1, 2^n); for k 1:(2^n-1) members logical(bitget(k-1, 1:n)); v(k) coalitionCost(members, params); end % 后期用 v 计算 Shapley值或核仁这种结构的好处是后期可以随意替换分配算法不用动最底层的调度模型。调度的核心我用YALMIP搭求解器可以接CPLEX或Gurobi。没装商用求解器的话小算例用linprog也能对付就是电池储能的0-1状态变量会麻烦一些。4.2 联盟成本函数的核心实现联盟成本函数是整个程序的心脏。函数内部先根据联盟成员列表启用对应设备变量再写目标函数和约束。设备启用与否我会用一段索引逻辑控制而不是在约束里硬塞大M。function cost coalitionCost(members, params) % members: [1 0 1] 表示A和C加入联盟B不加入 T params.T; pv sdpvar(1, T); % 向上级购电功率 gb sdpvar(1, T); % 燃气锅炉热出力 chp_p sdpvar(1, T); % 热电联产电出力 chp_h sdpvar(1, T); % 热电联产热出力 es_p sdpvar(1, T); % 储能充放电功率, 正放负充 es_soc sdpvar(1, T1); % 储能荷电状态 obj sum(params.price_e .* pv) sum(params.price_g .* gb / params.gas_eff); cons []; if members(1) % A主体参与接入储能设备 cons [cons, es_soc(1) params.es_init]; cons [cons, es_soc(2:T1) es_soc(1:T) - es_p * params.dt / params.es_cap]; cons [cons, es_p params.es_pmax, es_p -params.es_pmax]; cons [cons, es_soc params.es_socmax, es_soc params.es_socmin]; cons [cons, params.load_e chp_p pv es_p]; else cons [cons, params.load_e chp_p pv]; end if members(2) % B主体参与接入热电联产机组 cons [cons, chp_h params.chp_rh * chp_p]; cons [cons, chp_p params.chp_pmin, chp_p params.chp_pmax]; else cons [cons, chp_p 0, chp_h 0]; end if members(3) % C主体参与接入燃气锅炉 cons [cons, gb params.gb_min, gb params.gb_max]; else cons [cons, gb 0]; end cons [cons, chp_h gb params.load_h]; optimize(cons, obj, sdpsettings(solver, gurobi)); cost value(obj); end上面这段代码我刻意做了简化储能功率统一用正负号表示充放没加0-1变量防止同一时刻既充又放。严格做项目时这块一定要补上二元变量约束否则可能出现“充放同时进行”的伪优化结果。4.3 Shapley值的计算实现与优化小规模联盟直接枚举所有排列即可代码非常简洁。n人联盟共有n!种排列每个排列里依次累计边际贡献。n3时只有6种完全没压力。function phi shapleyValue(v, n) phi zeros(1, n); perms_list perms(1:n); for r 1:size(perms_list, 1) order perms_list(r, :); prev_set 0; for i 1:n curr_set bitxor(prev_set, 2^(order(i)-1)); phi(order(i)) phi(order(i)) v(curr_set 1) - v(prev_set 1); prev_set curr_set; end end phi phi / factorial(n); end这里的v数组下标要按联盟位编码对齐空集对应下标1。整体逻辑就是把每个排列的边际贡献累计到对应成员头上最后除以排列总数。代码本身没什么技巧真正要花时间的反而是前面联盟成本函数能不能快速、稳定地收敛。5. 算例设计三个能源主体合作前后的成本对比5.1 算例参数设定与仿真场景我用一个虚拟的三主体系统来演示整个流程。A主体拥有一套电储能B主体拥有一台热电联产机组C主体拥有一台燃气锅炉。负荷数据采用某工业园区的典型日曲线时间分辨率为1小时调度周期24小时。电网采用分时电价峰时段1.2元/kWh谷时段0.4元/kWh天然气价格折算为2.5元/m³。初始单独运行方式很简单A只能从电网购电B用热电联产同时满足部分电和热负荷C用燃气锅炉供热。各自独立优化后得到三个主体单独运行成本如下。主体设备配置单独运行成本/元A电储能上级购电4000B热电联产机组3500C燃气锅炉热储能3000单独运行总成本10500元。接下来把所有组合都跑一遍联盟成本函数得到的结果非常有意思。5.2 合作博弈结果解读从联盟成本到Shapley值分配各联盟的最优运行成本如下表所示。两两联盟产生的成本节省分别是AB联盟省1000元AC联盟省1000元BC联盟省700元三个主体全部合作后总成本9000元节省1500元。联盟最优运行成本/元成本节省/元A40000B35000C30000AB65001000AC60001000BC5800700ABC90001500按Shapley值计算A分到600元B分到450元C分到450元。这个结果跟直觉基本一致A的储能设备和电网交互最有调度弹性在多个联盟里都能创造价值所以分得最多B和C的热力设备虽然也重要但互补性略弱贡献相当。三个人分到的都是正数加总正好等于1500元没人会拿到比自己单干更差的结果。5.3 怎么判断分配结果站得住脚分配方案能不能让各方点头不能只看个人理性和集体理性还得验证联盟理性。把分到的收益代入所有子联盟检查每个二人联盟分到的总和都必须不低于该联盟自己的成本节省。比如A和B合计分到1050元大于AB联盟的节省1000元A和C合计1050元大于AC联盟的1000元B和C合计900元也大于BC联盟的700元。三个条件全部满足说明这个分配向量在核内任何一对主体都没有动力脱离大联盟单独合作。如果任何一条不满足方案就得打回重算通常要先检查是不是某个联盟的成本函数被算错了再考虑换核仁法重新分配。这种“先算数、后验核”的做法在实际项目汇报里几乎每次都能帮我把技术方案从会议室里顺利带走。6. 代码实现中的常见坑与工程化建议6.1 Shapley值计算组合爆炸超过8个主体就麻烦Scen前面代码里用了perms枚举排列这个写法在3到5个主体时非常爽快。但主体一旦到8个排列数是40320种每个排列又要顺序累计边际贡献算起来勉强能忍到10个主体排列数变成362万种每个排列还对应一组联盟成本求解彻底算不动。常用的处理思路有这么几个。一是用采样法近似Shapley值随机抽一部分排列算均值误差可控但对随机数种子敏感。二是按联盟价值函数结构压缩发现某类设备同质化严重时可以把同质成员聚合成一个参与人。三是改用双层优化或者先聚类负荷曲线、再算分配降低联盟成本求解的时间常数。我自己的建议是做方案演示3到5个主体用精确解就行做到实际工程项目时提前想清楚参与人粒度别把一个主体的多台设备拆成多个参与人只会给自己找麻烦。6.2 建模中的线性化与不可行问题综合能源系统里到处是非线性关系最典型的是热电联产机组的热电比不是常数机组在部分负荷下热电比会漂移储能充放效率不是同一个值天然气价格可能是阶梯式的。这些直接塞进YALMIP里Gurobi会当场报错或者求解极慢。我处理这些问题的优先级是先想办法保持MILP结构能线性化就线性化不能线性化就用分段线性逼近。充放电效率的0-1变量必须写完整不然会出“充放同时”这种鬼畜结果燃气锅炉的效率可以近似常数但热电联产的热电比灵活性必须保留否则机组完全失去运行弹性联盟红利会被严重低估。联盟成本函数一旦出现不可行首先检查负荷平衡约束里有没有漏掉某台设备的出力项这类问题十次里有八次是最基础的能量平衡写错了。6.3 从学术代码到实际项目交付的几个建议代码能跑通Shapley值只是第一步。做实际项目时还要注意这么几件事参数标定要比论文严谨得多设备效率、出力上下限最好有厂家数据或者历史运行数据支撑不能全靠文献估值分配结果要做敏感性分析电价、气价波动区间内方案是否依然稳定这直接影响谈判时各方对“未来风险”的接受度边界条件必须写清楚比如输配网络阻塞、电压约束这些如果在模型里被忽略分配结果在真实系统里会失真。可视化也是个容易被低估的环节。把各联盟成本、Shapley值分配结果、各主体逐时段出力画成图谈判时比单扔一张表格管用得多。我在项目里习惯做一张“成本节省瀑布图”从单独运行总成本一路降到联盟总成本每一步标注是哪项协同操作带来的节省再用饼图展示各方分配对方一眼就知道钱从哪来、到哪去。最后再分享一个我自己用过很多次的小技巧给联盟打标签。别用“联盟A”“联盟B”这种抽象编号而是用“储能热电联盟”“热电燃气联盟”这种方式命名。合作博弈的数学规则本身是抽象的但落到具体设备和负荷上时项目参与人真正在意的是“我的锅炉跟谁的储能在合作”。把模型语言翻译成业务语言方案才推得动。