做综合能源系统调度的朋友应该都有体会模型本身难的不是数学推导而是“一群利益不一致的主体为什么愿意坐到一起把能量湖打通”。这个标题里的“合作博弈”说到底就是在回答这个问题——先证明合作比单干划算再把合作产生的额外收益公平分掉。这篇博文我基于自己跑过的一个三主体园区算例把建模、求解、利益分配到Matlab实现串起来讲清楚文末还有我实际调试中踩过的几个坑希望帮你少走弯路。适合正在做IES方向研究、或想用合作博弈写调度/分配代码的同学参考。1. 先拆清楚合作博弈在综合能源系统里到底解决什么问题1.1 多能互补背后藏着的是“多主体的利益博弈”综合能源系统IES表面上是一个能量系统电、热、气耦合在一起有热电联产、热泵、储能、燃气锅炉这些设备。但你往运营层面看它其实是一个“多人利益场”。比如一个园区里热电联产机组可能属于能源服务公司储能站可能是第三方投资运营热泵和用户负荷则是用户侧资产各主体之间并没有一个“上帝视角”的统一指挥中心。传统优化调度默认有一个集中调度者可以全局安排每台设备的出力目标是系统总成本最小。这在理论上是漂亮的现实却经常推不动某个主体觉得自己出力了、设备被占用了、风险被拉高了最终却只分到很少的收益那下次它就不配合。合作博弈的价值正在于此它不是站在“上帝视角”做调度而是承认每个主体都是理性决策者通过公平的分配机制让各方自愿参与合作。1.2 什么时候该用合作博弈什么时候不该用需要先说清楚合作博弈不是替代优化调度而是覆盖在优化调度之上的“利益分配层”。底层调度问题依然要解只不过要对所有可能的联盟组合分别求解一遍得到每个联盟的收益然后拿这套收益表去算分配。什么情况下适合用合作博弈我的判断标准有三条系统里存在多个独立利益主体且它们有能力选择“不合作”主体之间联动后能产生明显的协同收益最常见的就是余热利用、互供电、共享储能需要一个让各方都认可的分配结果而不是简单的按比例分钱。反过来如果整个系统就是一个运营主体内部完全是成本中心管理那就没必要上合作博弈集中优化就够了。合作博弈的本质是为“分蛋糕”提供一套令人信服的规则。对比维度集中式优化调度合作博弈调度决策假设一个中心决策者多主体自愿协作目标函数系统总成本最小/收益最大联盟收益最大分配公平关键输出各设备出力计划分配给各主体的收益值落地痛点主体可能不配合需要枚举/筛选联盟典型场景一个电网公司管整个片区园区多投资方协同、微网群互济2. 建立这个模型需要哪几块“积木”2.1 多能流系统的物理模型要算联盟收益第一步是把每个主体的物理模型建出来。学术界最常用的框架是能量枢纽Energy Hub模型核心是一个输入-输出的转换矩阵。比如一个简单的电-热-气耦合系统输入是天然气、电网购电输出是电负荷、热负荷中间设备包括CHP、燃气锅炉、热泵、储能。我在实际项目中更喜欢把设备模型写成“通用单元”形式每个设备至少包含输入功率、输出功率、转换效率、爬坡约束、启停状态等结构这样后续算联盟时只需要在集合层面增删设备非常省事。2.2 调度模型的目标与约束调度模型是联盟收益的计算引擎。目标函数通常写为[ \min \sum_{t1}^{T} \left( C_{\text{grid},t} P_{\text{grid},t} C_{\text{gas}} F_{\text{gas},t} - C_{\text{sell}} P_{\text{sell},t} \right) ]其中第一项是购电成本第二项是天然气成本第三项是向电网售电的收益。常见约束包括电功率平衡发电购电放电 电负荷用电设备充电售电热功率平衡产热蓄热放热 热负荷设备出力上下限与爬坡约束储能SOC递推关系与容量约束联络线功率、燃气购买上下限。需要提一句对于多主体联盟每个联盟的最优调度是在“该联盟内部设备可用、且只能通过公共联络线对外交易”的假设下求解的。不同联盟的设备集合不同所以需要为每个联盟单独建一次调度优化问题。这也是特征函数计算成本高的根源。2.3 特征函数所有联盟收益的集合合作博弈的核心数据是特征函数 (v(S))表示联盟 (S \subseteq N) 内部成员合作时能获得的最大总收益。这里的 (v(S)) 不一定是“收益”也可以是负值即最小运行成本。实际操作中我会统一转成收益口径[ v(S) - \min \sum_{t} C ]这样所有联盟的 (v) 都是“越大越好”。算例中我取了三个主体CHP热电联产、HP热泵、BESS储能。独立运行时各自收益分别是 30、18、22 万元两两联盟的收益分别是 CHPHP 61 万元、CHPBESS 63 万元、HPBESS 58 万元三家全部合作收益是 101 万元。可以看到合作收益明显大于个体简单加总这超额部分就是协同红利。3. 利益分配不患寡而患不均3.1 Shapley值的核心思想和计算Shapley值是目前用得最多的分配规则之一。它的思想非常直观某个主体应得的分配等于它对所有可能联盟顺序的边际贡献的平均值。[ \phi_i(v) \sum_{S \subseteq N \setminus {i}} \frac{|S|!(n - |S| - 1)!}{n!} \left[ v(S \cup {i}) - v(S) \right] ]用生活类比解释一个团队里每个人的“功劳”应该通过他在不同入场顺序下给团队带来的增量收益来评估。你加入之前团队收益 50 万加入之后变成 61 万那你单次的边际贡献就是 11 万在所有可能的加入顺序下取平均就是你最终的分配值。继续用上面的数字计算三个主体的Shapley值。先看CHP的边际贡献加入空集(v({CHP}) - v(\emptyset) 30)加入HP所在联盟(v({CHP,HP}) - v({HP}) 61 - 18 43)加入BESS所在联盟(v({CHP,BESS}) - v({BESS}) 63 - 22 41)加入HP和BESS已结成的联盟(v({CHP,HP,BESS}) - v({HP,BESS}) 101 - 58 43)注意Shapley值需要对每一种可能的加入顺序取平均。三个主体共有 6 种顺序我这边直接按公式权重计算[ \phi_{CHP} \frac{2!}{6} \cdot 30 \frac{1!1!}{6}(4341) \frac{2!}{6} \cdot 43 10 14 14.33 38.33 ]同理可得 HP 约 30.67BESS 约 32。三者加起来为 101正好等于大联盟总收益。这个“总和等于大联盟收益”的性质非常重要在代码里可以作为校验条件。3.2 其他分配方法和选型依据除了Shapley值还有几种常见方法。核仁法Nucleolus最关注稳定性它通过求解一系列线性规划问题最小化联盟中“最不满的人”的不满程度。好处是分配结果一定处于核心Core里不会出现任何子联盟想脱离团队单干的情况坏处是实现复杂度高而且可能把超额收益全部分配给少数关键主体导致某些主体分到的收益反而比Shapley值低。按比例分配是最容易实现的。可以按初始投资比例、最大出力比例或独立运行收益比例来分。优点是简单透明缺点是没有反映真实的边际贡献合作中出力大、协同价值高的主体容易感到不公平。我的建议是工程项目要落地优先用Shapley值如果希望合作结构特别稳定、不能有“分裂倾向”可以考虑核仁法。在做对比分析时可以同时算两种把差异作为一个灵敏度分析维度呈现给决策者。3.3 大联盟收益是否一定大于子联盟之和答案是不一定。我带过的算例里就出现过储能与热泵合作后热泵选择在自己电制热最优的时段大量用电导致储能原本想做的低谷充电高峰放电空间被压缩两主体联盟收益反而低于各自独立收益之和。这种情况下合作博弈的正收益可能不成立。所以在代码流程里我建议先做一个“合作收益非负性”检查% 检查特征函数的超可加性 superadditive true; for i 1:n for j i1:n if v([i,j]) v(i) v(j) superadditive false; end end end如果发现某些组合不满足超可加性就需要检查模型参数是否合理比如是否缺少了热网、储热等互补设备或者中间产品计价方式是否存在扭曲。不是所有系统天生适合合作博弈这也是选型前要先算一次的原因。4. 调度与分配的联动Matlab代码实现的核心流程4.1 整体代码架构我习惯把整个项目拆成三层设备参数层、调度求解层、分配计算层。设备参数层用 struct 存好所有设备参数调度求解层针对每一个联盟分别建模型并调用求解器分配计算层负责枚举联盟、汇总收益、计算Shapley值。整体流程是这样的定义N个主体的设备参数枚举所有非空联盟子集共 (2^N -1) 个对每个联盟构建优化调度模型并求解得到联盟总收益把所有联盟收益填入特征函数表用Shapley公式或核仁法求解分配输出分配结果并校验“总和大联盟收益”。4.2 联盟枚举与特征函数计算联盟枚举我用循环生成二进制掩码实现N10以内都很轻松。N 3; v zeros(1, 2^N-1); % 存储所有非空联盟收益 for k 1:(2^N-1) mask bitget(k, 1:N); % 二进制掩码标记哪些主体参与 member find(mask 1); v(k) solve_alliance(member, params); endsolve_alliance会读取当前联盟包含的设备集合动态拼装约束矩阵并调用求解器。这里有个容易疏忽的地方联盟的收益必须扣掉联盟内部交易的结算成本吗不需要。特征函数定义的是联盟作为一个整体对外的净收益内部资金流转是分配阶段的事情在计算特征函数时不要做内部结算否则会重复计算。4.3 调度模型的Matlab实现我一般用Yalmip建模调用gurobi或cplex求解因为这类调度问题大多带0-1变量设备启停纯linprog不够用。一个典型的联盟调度代码如下function val solve_alliance(members, params) T params.T; yalmip(clear); % 决策变量定义 P_chp sdpvar(1, T); P_hp sdpvar(1, T); P_bess_in sdpvar(1, T); P_bess_out sdpvar(1, T); P_grid sdpvar(1, T); P_sell sdpvar(1, T); SOC sdpvar(1, T1); % 目标函数 Cost sum(params.c_buy .* P_grid) params.gas_price * sum(P_chp / params.chp.eta_p) ... - sum(params.c_sell .* P_sell); Constraints []; % 电平衡联盟内部更复杂省略 Constraints [Constraints, P_grid P_chp P_bess_out params.e_load P_hp P_bess_in P_sell]; % 热平衡 Constraints [Constraints, params.chp.eta_h * P_chp params.hp.cop * P_hp params.h_load]; % 储能SOC Constraints [Constraints, SOC(2:end) SOC(1:end-1) ... P_bess_in * params.bess.eta - P_bess_out / params.bess.eta]; optimize(Constraints, Cost, sdpsettings(solver, gurobi)); val -value(Cost); end这里我故意把电平衡写简化了实际工程里还要考虑CHP和热泵是否属于当前联盟。如果某个主体没参与对应设备的决策变量要置零或不存在而不是让它免费服务于联盟。这也是初学者最容易出错的地方。4.4 Shapley值的代码实现我提供两种实现。第一种是精确版N较小或联盟数可枚举时使用function phi shapley_exact(v, N) phi zeros(1, N); n N; for i 1:N sum_value 0; for S_all 0:(2^n - 1) mask_i bitget(S_all, i); if mask_i 1 continue; end member_ids find(bitget(S_all, 1:n) 1); S_size length(member_ids); S_index S_all 1; % 注意空集对应 index?需要修正 if S_size 0 vS 0; else vS get_v(S_all 1); % 特征函数查表 end S_with_i bitor(S_all, bitshift(1, i-1)); vS_i get_v(S_with_i 1); weight factorial(S_size) * factorial(n - S_size - 1) / factorial(n); sum_value sum_value weight * (vS_i - vS); end phi(i) sum_value; end end第二种是蒙特卡洛近似版N较大没办法枚举所有联盟时用。核心思想是随机生成大量加入顺序按顺序逐个计算边际贡献后取平均。实测下来联盟数在百万级时抽样 5000 次就能得到相当稳定的结果。function phi shapley_mc(v_func, N, n_samples) phi zeros(1, N); for r 1:n_samples order randperm(N); v_prev 0; for k 1:N member order(1:k); v_cur v_func(member); % 特征函数查询 phi(order(k)) phi(order(k)) (v_cur - v_prev); v_prev v_cur; end end phi phi / n_samples; end注意蒙特卡洛版必须保证每次抽样用的加入顺序是均匀随机的否则边际贡献估计会有偏。5. 常见问题与排查技巧实录5.1 明明合作了收益却变低了我之前调试一个包含储热系统的案例分析时就遇到过三主体联盟总收益低于两两联盟的情况。排查后发现是两个设备之间存在强耦合CHP需要在夜间低负荷运行来配合储能充电导致它的热出力严重偏离热负荷被迫消耗更多天然气。解决方案是在调度模型里加入“热负荷灵活区间”约束并让热储能承担解耦任务。做了之后大联盟收益明显回升。所以我建议如果发现特征函数不满足超可加性先看模型里是否缺少灵活性设备而不是直接怀疑合作博弈方法。5.2 枚举联盟组合爆炸怎么办主体数量到10个以后枚举 (2^{10}-11023) 个联盟还勉强能接受到15个主体时32767个联盟就需要一些计算时间了到20个主体时百万级联盟基本不可能逐一精确求解。工程上我用了两招第一招业务剪枝。很多联盟在物理上就不成立比如某联盟里只有负荷没有电源这种联盟特征函数收益很低且不参与Shapley值的主要贡献可以直接跳过不精确求解用一个固定惩罚估计值。第二招蒙特卡洛Shapley采样。先大量抽样联盟顺序统计哪些主体间边际贡献近似相同把高度对称的主体合并成等价类再用等价类减少求解次数。实测在16个主体场景下求解次数从65535次降到约4000次误差控制在2%以内。5.3 Shapley值出现负分配怎么办严格来说负分配的边界值并不违反Shapley性质的逻辑但在实际项目里没法向决策者交代。我在某印刷厂能源托管项目中就遇到类似情况某主体独立运转时能耗极低一旦合作接入共享母线它的用电曲线被优化到“高电价时段少用电”导致它在内部结算中被分到了负收益。处理方式是在调度阶段加入“个体理性约束”[ \phi_i \ge v({i}) ]即在分配层面强行保证每个主体分到的不少于独立运行收益。Shapley值本身不一定满足这个性质此时需要做修正分配比如做最小二乘调整或在分配约束下重新求解一个线性规划目标。从工程角度我更推荐后者因为决策者要的不是数学最优而是“比单干强”的心理底线。5.4 求解器选型与性能调优综合能源系统调度往往带整数变量我优先推荐Yalmipgurobi组合。没有gurobi授权也可以用cplex或scip但性能上稍有差距。如果只有Matlab自带工具纯LP问题用linprog含整数变量用intlinprog注意把变量范围收窄能显著加快求解速度。联盟求解的耗时大头往往在重新建约束矩阵。每次调用solve_alliance都要重新构造约束这是很浪费的。优化思路是预先构造一个“超集”模型包含所有主体和设备然后用0-1变量控制每个设备是否参与当前联盟。这样只需要求解一次就得到所有联盟的结果其实不是因为每变动一个设备集合最优解都会变仍需要重新优化但至少约束矩阵的结构可以复用省掉重复建模的时间。5.5 数值问题与敏感性分析特征函数收益如果量级相差太大比如一个几百块一个几十万在Shapley值公式里会出现大数吃小数的问题。我建议在分配前对特征函数做归一化处理或至少做中心化先减去每个主体的独立运行收益 (v({i}))把超额收益单独拿出来分配再回加独立收益。这样数字更稳定解释也更清楚。敏感性分析一定要做。项目汇报时决策者一定会问如果天然气价格涨10%分配结果怎么变如果储能容量降一半谁受影响最大这类分析其实不难把参数循环批量跑就行关键是最后要整理成一张“价格-分配”对照表比出一段文字更有说服力。6. 最后分享一点实操体会这个项目真正让我觉得值得分享的地方不在于Shapley值的数学有多深而是“把调度和分配串成一条线”的工程实现。我最早踩过的坑就是只做了调度优化拿着总成本下降了10%的结论去汇报结果被追问“这10%到底是谁创造的、谁来兑现”一时间答不上来。后来把合作博弈的分配层补上才发现一套完整的项目报告必须同时回答三个问题合作能多赚多少、谁贡献大、分配方案大家认不认。MATLAB代码实现本身不难难点在于把特征函数算准、把分配结果解读透。如果你也想在自己的项目中复现这套方法我建议从小规模三主体案例入手先按文中的流程跑通调度Shapley分配再逐步扩展到更复杂的多能互补场景。等你习惯了特征函数那套思维再回头看集中式调度就会觉得少了一层东西——那层东西就是「利益」。