前阵子帮某园区做多能互补方案时遇到一件挺有意思的事:几个供能主体原本各干各的,分开算成本,加上电网购电费、燃气费、设备折旧,一个月总支出大概是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 Σ_{t=1}^{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, T+1); % 储能荷电状态 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:T+1) == 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!种排列,每个排列里依次累计边际贡献。n=3时只有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 | 电储能+上级购电 | 4000 |
| B | 热电联产机组 | 3500 |
| C | 燃气锅炉+热储能 | 3000 |
单独运行总成本10500元。接下来把所有组合都跑一遍联盟成本函数,得到的结果非常有意思。
5.2 合作博弈结果解读:从联盟成本到Shapley值分配
各联盟的最优运行成本如下表所示。两两联盟产生的成本节省分别是:AB联盟省1000元,AC联盟省1000元,BC联盟省700元;三个主体全部合作后总成本9000元,节省1500元。
| 联盟 | 最优运行成本/元 | 成本节省/元 |
|---|---|---|
| A | 4000 | 0 |
| B | 3500 | 0 |
| C | 3000 | 0 |
| AB | 6500 | 1000 |
| AC | 6000 | 1000 |
| BC | 5800 | 700 |
| ABC | 9000 | 1500 |
按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”这种抽象编号,而是用“储能+热电联盟”“热电+燃气联盟”这种方式命名。合作博弈的数学规则本身是抽象的,但落到具体设备和负荷上时,项目参与人真正在意的是“我的锅炉跟谁的储能在合作”。把模型语言翻译成业务语言,方案才推得动。