从接手第一个综合能源系统协调调度项目开始我就被“利益分配”这四个字卡了很长时间。一开始想得简单多能互补、联合优化把总成本压下来大家都省钱不就完事了但真正做下来才发现最大的难点根本不在物理层面的潮流平衡或者储能充放电策略而在于一个特别现实的问题——“省下来的钱该怎么分”。分得合理各方才有动力配合调度分得不合理哪怕总收益再高第二天联盟就散了。这也是为什么合作博弈在综合能源系统调度里会被频繁提及因为它本身就是研究“多人合作产生的超额收益如何公平分配”的一套数学框架。这篇文章就围绕我自己在Matlab里完整实现的一套方案来聊。先梳理清楚利益分配和优化调度是怎么耦合在一起的再拆解合作博弈里几种主流分配方法核心、Shapley值、核仁到底该选谁然后给出代码层面的完整实现链路和调试过程中踩过的坑。适合正在做园区级综合能源系统优化、多主体协调调度或者研究电力市场P2P交易、微网能量管理的朋友参考。无论你是把这篇当课程设计的起点还是作为工程项目的前期验证都应该能从里面找到可以直接拿去用的东西。1. 一张图看懂利益分配问题的核心矛盾为什么调度结果会反过来影响分配结果先说清楚一个问题本质综合能源系统里的“利益分配优化调度”不是一个纯粹的“先求出最优调度、再按比例分钱”的两段式问题。它比这个复杂原因是调度方案本身和分配方案之间存在相互影响。1.1 利益主体之间的耦合多个能流、多个主体、一套共享物理网络综合能源系统的典型形态是电、气、热、冷多种能源在同一个物理区域里耦合运行。比如热电联产机组CHP既发电又供热电锅炉可以消纳新能源富余电力来产热储电设备和储热设备在时间维度上提供灵活性。如果这套系统的业主只有一个那问题很干净——单目标优化最小化总运行成本即可。可一旦系统里有多个参与主体情况就变了。某园区可能有这样三类主体投资CHP机组的热力公司A、运营储能电站的储能服务商B、负责园区用户售电和负荷聚合的综合能源服务商C。A、B、C各自有独立的收益核算却又共享同一套电网接入点、同一套热力管网和同一个调度中心。任何一个主体的运行策略都会影响其他主体的收益空间。举个例子CHP机组的启停策略会改变园区关口处的购电功率曲线而这条曲线恰好决定储能服务商的峰谷套利空间储能的充放电行为又会影响热力公司余热利用的时间窗口。这就形成了一种相互咬合的状态每个主体的“独立最优”并不是全局最优只有协调调度才能释放整体的合作剩余。1.2 合作剩余从哪来协同运行释放的三类红利那所谓“合作剩余”具体是怎么产生的我做了不少仿真之后把它归纳成三类来源。第一类是时间互补。储能设备和热储能设备在时间轴上平移能量。白天光伏大发时多充电、夜间低谷时蓄热把低价时段能量搬到高价时段用这是最直接的额外收益来源。第二类是能质匹配。电锅炉和燃气锅炉之间的替换不是简单等价而是根据实时电价和气价的比值动态取舍。气价低、电价高时多烧气少购电电价低、气价高时反过来。单主体各自运行时通常没有这种跨能流比较的全局视角。第三类是容量共享。比如园区配变容量固定储能服务商和售电服务商如果各自独立报装容量会浪费协调运行后可以错峰使用同一个关口容量容量费用和需量电费都能降下来。这三类红利在仿真中的体现通常是总运行成本下降8%到15%具体数字取决于系统规模和设备配置。但这里有个关键点这些红利不是平均从天上掉下来的不同调度方案下不同主体对总收益的边际贡献是完全不同的。1.3 分配方案反作用于调度方案一个容易被忽视的循环这是最让我在产品设计上头疼的地方也是我觉得很多人做这类课题时没有想透的地方。如果我们先求出一个全系统最优调度方案算出总收益然后分钱这套流程其实隐含一个前提——各主体接受这个调度方案。但现实中每个主体都有“退出合作单独运行”的选项。如果某个主体发现自己在合作中的收益还不如单干时它的理性选择就是退出这时候所谓的“全局最优调度”就变成了一张废纸。反过来推理每个主体的单干收益决定了它参与合作的最低心理门槛。这个门槛在博弈论里叫“个体理性约束”会反过来限制调度方案的可行域——调度中心必须保证最终方案让每个主体分到的收益不低于其单干收益否则调度根本落不了地。所以“利益分配”不是一个调度结束以后才开始的事而是从建模开始就要考虑的约束条件。这也是标题里把“利益分配”和“优化调度”并列的原因。2. 合作博弈分配算法选型为什么我更常用Shapley值而不是核仁既然分配和调度要联动那就得选一套合理的分配法则。合作博弈的分配方法不少我在实际项目里主要对比过三类核心Core、Shapley值、核仁Nucleolus。2.1 核心Core理想但常常不存在核心的定义是所有能够同时满足个体理性和联盟理性的分配方案集合。联盟理性要求的是任何一组主体形成子联盟后其总收益不能低于这个子联盟单独运行时的收益这样就没人有动力拉帮结派搞“内部小团体”推翻大联盟方案。这个概念在数学上非常漂亮但实际应用中有个大麻烦核心经常是空集尤其是主体数量比较多、设备耦合复杂的时候。算出来是空集意味着这个合作结构本身就无法维持你得回到系统设计去改方案而不是在博弈论里找解。而且在Matlab实现时核心需要用线性规划求解一组不等式解出来的还可能是一条线段多个解不好直接落地。2.2 Shapley值边际贡献加权平均工程上最容易解释Shapley值的思路是每个主体应得的收益等于它对所有可能的加入顺序下“边际贡献”的平均值。通俗点说参加合作的元老们逐一加入大联盟每个人加入时给联盟带来的收益增量取所有加入顺序的平均就是这个人应得的分配。这个分配方法有几个非常好的性质效率性所有主体分配之和等于大联盟总收益、对称性贡献相同的主体分配相同、可加性两个独立博弈合并后分配等于分别分配之和、以及哑元性质不贡献任何增量的主体分到零。这些公理化的性质使得Shapley值的计算结果天然带着“公平”的标签任何主体想挑战分配结果都很难从数学上找出毛病——这在多主体谈判现场真的是个巨大的优势。2.3 为什么在实际Matlab实现里我优先选Shapley值从实现层面讲Shapley值的计算框架最规整。特征函数任何一个子集S合作后能产生的总收益只需要对每个子集求解一次优化调度问题然后按公式加权求和即可。计算模式高度统一适合写循环调用求解器。核仁虽然在某些状态下更稳健但需要多次迭代求解线性规划每一步都要更新约束条件代码复杂得多。我实测下来同样是五主体系统Shapley值的代码量大约是核仁的三分之一。当然Shapley值也有它的不足——当主体数量超过8个时子集数量按指数增长2的N次方计算代价会涨得很快。如果真遇到这种场景可以考虑用抽样估计的近似Shapley值这个我在后面代码解析部分会展开讲。3. Matlab代码框架全解析从问题矩阵到Shapley计算的完整链路接下来是重点部分我把自己在项目里跑通的代码框架按模块拆开讲。整个代码的核心思路是“外壳处理博弈分配、内核处理调度优化”——用Shapley值做分配外壳用混合整数线性规划MILP做优化内核。3.1 整体架构与文件组织工程文件建议这样组织main_IES_cooperative.m主入口读取数据、调用各模块、输出结果build_problem.m构建优化模型的目标函数系数和约束矩阵solve_subset.m对给定子集S求解优化调度问题返回该子集的最大总收益特征函数值compute_shapley.m根据特征函数值按Shapley公式计算各主体分配结果check_stability.m校验个体理性和联盟理性约束是否满足data/config.xlsx存放电价、气价、负荷曲线、设备参数等输入数据主流程分四个阶段。第一阶段加载数据第二阶段对所有非空子集循环调用solve_subset得到全部特征函数值第三阶段把特征函数值输入compute_shapley得到分配结果第四阶段把分配结果代回check_stability做校验如果不满足约束修改调度约束后重新算。3.2 优化调度内核如何把一个子集S的运行问题构建成MILPsolve_subset是整个代码里最核心的函数。它的任务是对任意给定的主体集合S求这个集合作为一个独立联盟运行时能得到的最大总收益。以一天24小时、每时段1小时为例决策变量包括x_gt(t)燃气轮机发电功率连续变量x_gb(t)燃气锅炉产热功率连续变量x_eb(t)电锅炉产热功率连续变量x_chp(t)CHP机组热出力连续变量y_gt(t)燃气轮机启停状态0-1变量x_es_c(t)、x_es_d(t)储能充电、放电功率连续变量soc(t)储电设备荷电状态连续变量目标函数是联盟总净收益最大化。净收益等于售电收入加售热收入减去购电成本、燃气成本和设备运维成本f -(sum(price_sell_e .* x_load_supply) ... % 售电收入取负转为最小化 sum(price_sell_h .* x_heat_supply) ... % 售热收入 - sum(price_buy_e .* x_grid_buy) ... % 购电成本 - sum(price_gas .* (x_gt_fuel x_gb_fuel)) ... % 燃气成本 - sum(om_cost .* x_power_total)); % 运维成本注意Matlab的intlinprog是求解最小化问题的所以收益项全部加负号。关键约束用矩阵形式写电功率平衡园区总用电负荷加储能充电功率等于光伏发电加电网购入加燃气轮机发电加储能放电。热功率平衡热负荷由CHP热出力、燃气锅炉、电锅炉和储热设备共同满足。储能约束SOC递推方程加充放电状态互斥约束。这个互斥约束需要引入另一个0-1变量y_es_c保证同一时刻储能不会同时充放电。设备出力上限约束每台设备出力在0到额定容量之间。CHP热电耦合约束CHP的发电出力和热出力之间满足热电比约束这是综合能源系统建模里最容易出错的一条约束。电网交互约束关口购电功率不能超过配变容量限制。构建约束矩阵时一个小技巧是把变量索引提前做成映射表用结构体存起来。比如idx.x_gt 1:24; idx.x_gb 25:48;这样后续写约束时直接用索引赋值避免反复寻找变量位置。3.3 特征函数计算不是每个子集都调用完整的MILP所有子集的集合S遍历时我会做一个剪枝直接计算单主体子集与大联盟子集中间的子集按参与主体的类型特征判断是否需要求解。例如两个燃气相关设备主体组成的子集如果其没有任何储能或新能源资源该子集的收益就是两个主体独立运行收益之和不需要再算一次联合优化。这能让计算量下降不少。每个子集求解得到的最优目标值v(S)就是特征函数值也是后续Shapley计算的原料。3.4 Shapley值的Matlab实现公式落地与验证Shapley值的标准公式是phi_i sum_{S包含i} (|S| - 1)! * (N - |S|)! / N! * [v(S) - v(S\{i})]在Matlab里实现时我会先枚举所有包含主体i的子集然后按公式累加。为了便于理解这里给出一个主体数为3的简化版本% N 3, 计算主体1的Shapley值 N 3; phi_1 0; for k 1:N % 枚举包含主体1且大小为k的所有子集 subsets nchoosek(2:N, k-1); for j 1:size(subsets, 1) S [1, subsets(j, :)]; S_minus_i subsets(j, :); % 当前子集收益和去掉主体1后的子集收益 v_S v_func(S); v_Si v_func(S_minus_i); % 权重因子 weight factorial(k-1) * factorial(N-k) / factorial(N); phi_1 phi_1 weight * (v_S - v_Si); end end这段代码的逻辑和数学公式完全对应。在实际封装时我把v_func从函数句柄改为查表索引——因为预先算好的特征函数值都存到了数组v_table里查表比重复调用求解器快得多。验证Shapley值计算是否正确有一个很简单的办法算完所有主体的分配后求和检查是否等于大联盟的总收益v(全集)。如果不等一定是你枚举子集时漏了或者权重算错了。这个效率公理性检验每次跑完结果我都必做算是自检的最好习惯。4. 实际运行中的四个坑数值振荡、维度错配与个体理性失效代码框架搭好后真正调试的过程才是消耗时间的大头。下面这几个坑我都是亲手踩过并逐一定位解决的分享出来帮你们绕开。4.1 坑一intlinprog数值振荡导致结果不稳定第一次跑通全流程时我发现每次运行得到的总收益会差几千块钱而且连续两次运行结果不一致。查了代码逻辑没发现问题后来才意识到这是求解器的数值容差问题。因为目标函数里不同项的单位不一致——购电成本是千瓦时乘以电价几百到上千运维成本却是千瓦时乘以一个很小的系数百分之几两者相差好几个数量级。intlinprog在求解含有0-1变量的混合整数规划时整数预处理的容差过小会导致某些变量在最优解附近振荡从而给出不同结果。解决方式很粗暴也有效把所有量纲统一。收益、成本统一换算成“万元/天”为单位保证目标函数系数落在0.01到100之间。这样改动后结果稳定了求解速度也快了不少。4.2 坑二约束矩阵维度错配报错信息却指向不明有一次报错提示Index exceeds matrix dimensions我花了一个多小时排查最后发现是电功率平衡约束里储能充电变量用了1:24的索引但储能设备在子集S不包含时根本没定义这个变量——变量数量变了约束矩阵的列数没跟着变维度自然就对不上了。这个坑的根因是变量集合需要随子集S动态变化。子集S不包含储能设备时储能相关变量直接置零并从决策变量集合里删除。在build_problem函数里我专门加了一段判断逻辑先检查子集S包含哪些设备类型再按类型拼接约束需要的变量索引段。从那以后我再没遇到过维度错配。4.3 坑三Shapley值计算量随主体数量指数膨胀前面说了三主体的Shapley计算很直观可当系统增加到五主体时就有点恼火了。五主体需要枚举31个子集每个子集都要跑一次MILP。如果MILP里又有24个时段的0-1变量每次求解耗时几十秒整个计算流程需要十几分钟。而如果系统有八到十个主体枚举所有子集基本不可能在可接受时间内跑完。我的方法是引入抽样子集估计。随机抽取一定数量的子集比如100个用这些子集的边际贡献平均值作为近似Shapley值。实测下来当样本量到150个左右时误差能控制在3%以内而计算时间只是精确计算的十分之一。当然用抽样前必须先算一次大联盟收益用来标定结果保证分配值之和仍然等于大联盟总收益。4.4 坑四个体理性约束不满足导致的联盟破裂这是整个项目里最刺头的一个问题。某次仿真结果里热力公司A分到的收益比它单干时还低了8000块这意味着该主体有充足的理由退出合作。问题出在特征函数值计算时单独运行部分的负荷参数没设对——单干时A可以按自己的方式报价卖热而合作后热价被内部结算价替代导致A的收入核算逻辑出现偏差。解决方案分两步。第一步修改单干场景的参数设定单干时允许主体按市场价卖电卖热合作时按内部结算价两者分开核算。第二步在分配环节加入个体理性约束检查一旦某主体分配值低于单干收益就把该主体的“最低参与条件”作为附加约束加入调度模型重新求解。这样虽然总收益会略降一点但至少保证分配方案是可行的——大家可以坐下来继续谈合作而不是直接掀桌子。5. 一个三主体算例的全过程拆解从特征函数到分配数值为了让整套方法更落地我拿一个简化三主体算例完整走一遍流程。这三个主体分别是热力公司A持有CHP机组、储能服务商B持有电化学储能、售电服务商C聚合园区用户负荷。系统基础参数是CHP额定容量2MW储电容量4MWh园区峰时负荷2.5MW谷时负荷1.2MW。5.1 特征函数计算五个关键子集收益先算单主体独立运行收益。A单干只靠CHP发电自用加卖热日收益约1.82万元B单干靠峰谷套利日收益约0.65万元C单干做售电代理赚购销价差日收益约0.93万元。三者单干收益之和为3.40万元。再算两两合作子集。AB合作CHP和储能协同后A可以把谷时多发电存起来峰时卖B获得更稳定的充放电策略两者合计3.45万元。AC合作售电契约让CHP发电可以更贴合园区负荷曲线合计3.12万元。BC合作储能直接服务用户负荷优化合计2.23万元。最后算大联盟ABC合作收益三类设备充分协同——CHP发电、储能调峰、售电聚合共同优化日总收益3.84万元。5.2 合作剩余的来源解构大联盟总收益3.84万元减去各单干收益之和3.40万元合作剩余是0.44万元。这4400块就是这个合作博弈的“总价值增量”。用12小时时段的理论锋谷差来说明更直观原来自各支峰的收益空间因为协调调度被释放出来储能增加了每天约两次深度循环CHP提高了约15%的发电利用小时数——这些隐含的价值最终都折算进了那4400块。5.3 Shapley值分配结果与稳定性校验手算三个主体Shapley值的过程不展开公式只看结论A分配到1.94万元B分配到0.82万元C分配到1.08万元。三项之和恰好等于3.84万元通过效率公理化检验。对比各自单干收益A增加了0.12万元B增加了0.17万元C增加了0.15万元。增量最大的B符合直觉——储能是整个系统灵活性的关键提供者也是让CHP能“削峰填谷”、让售电代理能“错峰购电”的核心推手。C的增量居中A的增量最少但也都是正收益个体理性约束全部满足。这个结果在项目讨论时很有说服力——没有任何一个主体会觉得自己吃亏了所有人都有动力维持大联盟方案。这也是合作博弈分配相比“按投资比例分收益”或者“按用能量分摊”在公平性方面更强的地方。6. 实用经验扩展从三主体到十主体这套方案怎么演进如果你做的系统比这个算例更复杂有几点扩展建议你可以直接用。6.1 多主体时不能只靠Shapley一种方法兜底主体数量上去之后如果把Shapley值结果做稳定性校验发现有个别主体不满足个体理性我的建议是引入“核仁作为修正层”——先用Shapley值得到公平基准再用核仁方法在满足个体理性的前提下修正分配量。这种“Shapley起步、核仁修正”的组合兼顾了公平性和可行性比单一算法更抗造。但请注意核仁求解需要用到线性规划工具箱代码量会增加所以只在确实需要时再上。6.2 考虑碳约束时的扩展方式现在越来越多的项目会叠加碳排放配额约束。在这种情况下特征函数计算里就要把碳排放约束加进去子集S的收益是在允许的排放配额内能获得的最大净收益。碳价可以作为影子价格的参数写进目标函数这样分配结果也能反映碳成本在各主体之间的分摊机制。实测下来加入碳约束后Shapley分配结果会更偏向新能源占比高的主体——因为它们的边际贡献在碳约束下更值钱。6.3 多时段耦合与随机性处理我的代码是以24小时为周期、以小时为步长做确定性优化的。如果项目要做日前调度反映光伏出力和负荷的不确定性建议把场景法加进来对几个典型场景分别求解特征函数再按场景概率加权求Shapley值。这个改动只涉及外层数据生成逻辑核心的MILP和Shapley代码不用动。我做这套方案时前前后后大概花了三周其中一半时间消耗在调数值稳定性和验证分配公平性上。建议刚开始做的朋友不要一上来就上最复杂的模型先把三到五个主体的确定性模型跑通再逐步加约束、加场景。先把分配机制看清楚再谈调度复杂度这个顺序能帮你省掉大量走弯路的时间。