合作博弈与Shapley值：综合能源系统优化调度及利益分配的MATLAB实现

简介：面向综合能源系统与电力市场方向的毕业设计，这里提供一套基于合作博弈的综合能源系统利益分配优化调度MATLAB源程序。代码覆盖P2G电转气、碳捕集、燃气轮机、热储能等设备建模，并实现三方主体合作联盟构建，利用Shapley值法对合作剩余按贡献分配，可帮助读者复现内蒙古地区实例仿真中的运行成本、风电消纳与碳排放分析。压缩包共13个文件，整体约215KB，以10个.m源码文件为主，配套1个.lp优化问题定义、1个.xlsx算例数据和1个.png示意图，便于对照模型结构调试。目前已有231人学习浏览，适合需要深入理解综合能源系统低碳经济调度与合作博弈利益分配建模的本科生与研究生。代码注释与模块划分清晰，可快速定位设备建模、联盟收益计算和Shapley值求解环节，可作为论文复现与二次开发的起点。

1. 合作博弈的综合能源系统优化调度源程序：到底能帮你省多少事

做综合能源系统（IES）优化调度，最让人头疼的不是理论推导，而是从论文公式到一份能跑的“合作博弈源程序”之间的那段路。这套基于合作博弈的综合能源系统利益分配优化调度MATLAB源码，把电转气（P2G）、碳捕集、燃气轮机、热储能等设备建模，以及Shapley值分配算法打包成了可直接运行的项目。你拿到的不是一个孤零零的shapley.m，而是完整的联盟成本计算框架：sc1、sc2、sc3分别处理三个运营主体的独立运行，sc12、sc13、sc23处理两两联盟，main123.m负责总联盟的协同优化，最终用Shapley值把合作剩余按贡献分给各方。对正在做IES方向毕业设计的人，这相当于把论文里的算例变成了可以上手改的MATLAB工程；对想快速验证合作博弈策略的从业者，也能省掉从头搭模型的几天时间。

2. 先想明白：综合能源系统为什么需要“合作博弈 + Shapley值”

2.1 从“各主体独立运行”到“合作联盟”：成本下降的真实来源

综合能源系统里通常有电力、热力和燃气这三个运营主体，各自有独立的设备和使用场景。过去常见的做法是各算各账，电力系统只管电功率平衡，热力系统只管热负荷，燃气系统只管气源调度。但这种独立运行模式忽略了设备之间的耦合潜力，最典型的就是电转气（P2G）设备——在风电高出力时段，电网可能面临弃风，而气网恰好缺气，这时候如果两个主体合作，让多余风电去制氢、制甲烷，就能把弃风变成气网的有用入力，省下一部分购气费用。这就是“合作剩余”的来源。

合作博弈的核心思想不是要求每个主体都让利，而是设计一个“联盟”模式，使联盟的总运行成本小于各主体独立运行成本之和。从数学上看，设特征函数为联盟的负成本，合作的净收益等于独立成本之和减去联盟成本。这个收益不是凭空出现的，而是设备间能量互补所带来的“耦合红利”。P2G设备在电-气协同中扮演了柔性联结角色，燃气轮机则把气转电，热储能把热电耦合在一起，碳捕集装置再折算成碳排放成本降低。这些设备之间的双向流转，只有放在联盟框架下才能充分体现。

2.2 模型里有哪些设备和主体：P2G、碳捕集、燃气轮机、热储能

这份源程序的算例里，综合能源系统的框架包含四类关键设备：P2G、碳捕集、燃气轮机和热储能。P2G的建模重点是电转气的转换效率，通常在0.5到0.7之间，参数过高会放大风电消纳能力，过低则导致联盟收益不明显。碳捕集设备负责捕获燃气轮机和锅炉排放的CO2，捕集率一般设在0.85左右，它能直接降低系统碳排放成本，但也会增加额外的运行和维护电耗。燃气轮机承担热电联产功能，发电和产热同时进行，建模时要写清楚热电比和爬坡约束。热储能则像一个缓冲池，在热负荷低谷蓄热、高峰放热，让热力侧的运行曲线更灵活。

系统内的三个运营主体可以按能源类型划分：主体1可能是电力侧，管理风电场、燃气轮机和P2G设备；主体2是热力侧，管理热储能和供热锅炉；主体3是燃气侧，管理气源、储气罐和燃气负荷。这样划分的好处是每个主体的约束相对独立，联盟通过设备交互变量耦合在一起。资源中的sc1.m、sc2.m、sc3.m分别对应三个主体的独立运行优化，sc12.m、sc13.m、sc23.m则处理两两联盟，main123.m统筹三方联盟。这种模块化命名方式，让改动和排查都比较直观，不用在一个大文件里翻几百行代码。

2.3 利益分配不是AA制：为什么选Shapley值

联盟总成本降下来了，但怎么分这笔“省下来的钱”是所有合作博弈模型绕不开的问题。如果把收益平均分给三方，出力贡献大的主体会有意见，贡献小的主体则可能“搭便车”，合作就难以维持。Shapley值的价值在于它提供了一套满足公理化公平性的分配规则：每个主体获得的收益等于它对所有可能联盟组合的边际贡献的加权平均值。具体来说，对三个主体的情形，要分别计算主体在空集、单人联盟和两人联盟中的边际贡献，再按组合数赋予权重。

选择Shapley值而不是其他分配方案（如核仁、核心法）的原因主要有三点：一是计算逻辑清晰，三主体情况只需要7个联盟组合（空集、三个单人、三个两人、一个三人），手算也能核对；二是它满足可加性、对称性和有效性这三个公理，分配结果不会被“谁先提合作谁吃亏”这种博弈顺序影响；三是在IES这类多能互补场景中，各主体的投入相对可量化，边际贡献可以通过运行成本直接度量。资源里shapley.m和shapley1.m就是在做这件事，前者通常负责三个主体的完整Shapley计算，后者可能是针对某个主体或某种联盟结构的具体实现。

3. 拿到源程序后：文件结构、运行环境和第一遍跑通

3.1 文件清单与每个文件的作用

先把这份资源里出现的文件理一遍。主程序是main123.m，从命名就能看出它处理的是“主体1+主体2+主体3”的总联盟优化；sc1.m、sc2.m、sc3.m对应三个单主体独立运行；sc12.m、sc13.m、sc23.m两两联盟成本求解；shapley.m和shapley1.m做利益分配计算；eigen.m从这个名称猜测，可能是用于求解矩阵特征值或者验证联盟稳定性的辅助函数；ab.lp是一个线性规划模型文件，很可能是联盟“a-b”对应的LP模型描述；data.xlsx保存所有输入数据，包括负荷、风电预测、设备参数和碳交易价格等。具体项目可能还有少量辅助脚本，但核心骨架就是这样一张表。

3.2 MATLAB运行前的环境检查

跑通这份源程序并不需要特别高配的机器，但环境上必须先满足两个条件：第一，MATLAB版本建议在R2016a以上，太低不支持readtable等函数；第二，优化求解器必须可用。部分程序会直接使用MATLAB自带的linprog，这种情况只需要Optimization Toolbox；但如果主程序调用了ab.lp这个LP文件，可能依赖lpsolve或者CPLEX的MATLAB接口，这时候只装Optimization Toolbox会报错。我一般会在第一次运行前先把环境探一遍，避免被“未定义函数”这种错误耽误半小时。

% 环境检查脚本：确保必要文件都能被MATLAB找到 files = {'main123.m', 'sc1.m', 'sc2.m', 'sc3.m', ... 'sc12.m', 'sc13.m', 'sc23.m', 'shapley.m', 'data.xlsx'}; missing = {}; for i = 1:length(files) if exist(files{i}, 'file') fprintf('[OK] %s\n', files{i}); else fprintf('[缺失] %s\n', files{i}); missing{end+1} = files{i}; %#ok<SAGROW> end end if ~isempty(missing) fprintf('请先 cd 到程序解压目录，或把以上文件加入MATLAB路径\n'); end

这个脚本用exist检查每个文件是否在当前工作目录或搜索路径下，返回“file”表示是文件而不是变量。如果提示缺失，最常见的原因是当前工作目录不对——很多人下载压缩包之后没解压，或者解压后直接双击脚本，MATLAB当前目录还在“文档”里。下一步用绝对路径cd到文件夹就行。注意exist对文件名区分大小写，Windows下不严格但macOS和Linux下必须完整匹配。

3.3 main123.m 跑通的最小流程

第一次运行不要急着看结果，先理解main123.m的执行顺序。以这份资源按常规模块化写法推测，它的流程一般是：读入data.xlsx中的基础参数和负荷曲线，然后依次计算各单主体成本、两两联盟成本，再求解三方联盟总成本，最后调用shapley.m完成分配。具体代码不一定按这个顺序写，但数据流大致是这个方向。下面给一个按相同逻辑组织的示意结构，用来对照理解资源中的main123.m。

% main123.m —— 主流程示意框架 % 读取输入数据 T = readtable('data.xlsx', 'Sheet', '参数表'); % 若原文件无sheet名则去掉第二个参数 % 解析时段数和各主体数据 % 假设data.xlsx中第一列是时段，后面依次是电负荷、热负荷、风电出力等 T_horizon = height(T); % 调度时段数（24或96） data_struct = struct('elec_load', T{:,2}, 'heat_load', T{:,3}, ... 'wind_power', T{:,4}, 'gas_price', T{:,5}); % 独立运行成本 c1 = sc1(data_struct); % 主体1独立成本 c2 = sc2(data_struct); % 主体2独立成本 c3 = sc3(data_struct); % 主体3独立成本 % 两两联盟成本 c12 = sc12(data_struct); c13 = sc13(data_struct); c23 = sc23(data_struct); % 三方联盟总成本（一般在main内部或由sc23等组合得到） c123 = compute_total_union(data_struct); % 需查看main123实际实现 % Shapley值分配 phi = shapley([c1 c2 c3], [c12 c13 c23], c123); fprintf('独立成本：c1=%.2f c2=%.2f c3=%.2f\n', c1, c2, c3); fprintf('联盟成本：c12=%.2f c13=%.2f c23=%.2f c123=%.2f\n', c12, c13, c23, c123); fprintf('分配结果：phi1=%.2f phi2=%.2f phi3=%.2f\n', phi(1), phi(2), phi(3));

这里的compute_total_union不是现成函数，资源里更可能是通过sc23.m或main123.m内部的优化循环完成的。你拿到代码后，先找到计算总联盟成本那一段，确认它输入的是什么变量、输出有没有存到工作表里。跑通之后你会看到所有成本输出，其中分配结果能直接看出三方各分到了多少合作剩余。如果输出的phi很小甚至负数，先不要怀疑算法，大概率是传入的联盟成本顺序不对，这一点后面避坑章节会专门讲。

4. 看懂Shapley值分配代码：从联盟成本到边际贡献

4.1 shapley.m 与 shapley1.m 的分工

资源里同时存在shapley.m和shapley1.m，名字上都跟Shapley值有关。从工程经验看，shapley.m多半是通用实现，输入7个联盟成本（三个单人、三个两两、一个三人）并返回三个主体的分配值；shapley1.m则可能是单主体视角的边际计算函数，或者用来做图形化的对比，比如把某个主体在不同联盟中的边际贡献画出来。如果你修改了主体数量，比如从三个主体改成四个主体，shapley.m里的权重公式就得重写，shapley1.m则可能需要重新定义“主体1”的输入，否则维度对不上。

4.2 输入与输出：data.xlsx 和 ab.lp 的关键作用

data.xlsx是外界与核心算法的唯一数据接口，它的列顺序和单位直接决定优化结果。以三主体IES算例为例，通常包含电负荷（MW）、热负荷（MW）、风电预测出力（MW）、天然气价格（元/m³）、碳排放价格（元/t）、设备容量、效率等参数。ab.lp这个LP文件提供的是联盟优化问题的标准形式，很多旧版MATLAB程序习惯把约束写成lp格式，再用lpsolve或其他求解器读取。如果你只装了MATLAB自带优化工具箱，遇到lp文件可能会不知所措——其实可以绕过它，把边际成本和发电机爬坡等约束直接用linprog的密集矩阵重写。

% 读取data.xlsx示例：优先使用readtable，便于按列名操作 T = readtable('data.xlsx', 'PreserveVariableNames', true); % 查看前几行 head(T); % 提取数值与列名对照 col_names = T.Properties.VariableNames; fprintf('列名列表：\n'); disp(col_names);

readtable的主要优势是能直接读取数值和文本混合的表格，但要求工作表内没有合并单元格，表头是唯一的、不重复的行。如果原文件里有中文注释或单位行，readtable会把它们当作数据行读进来，导致后续矩阵运算报错。这时最简单的处理是另存一份data_clean.xlsx，只保留纯数值区域，并把第一行表头改成简单英文，比如time、elec_load、heat_load、wind、carbon_price。

4.3 一个三主体Shapley值的MATLAB对照实现

为了让你能亲手验算分配逻辑，这里直接给一个标准的三主体Shapley值计算函数。它的输入是7个成本值，输出为三个分配结果。注意合作博弈里的特征函数一般取收益为正，而这里我们手头是成本，所以要把成本取负再代入公式，否则数学关系就反了。

function phi = shapley_three(c1, c2, c3, c12, c13, c23, c123) % 输入：独立成本c1,c2,c3；两两联盟成本c12,c13,c23；三方联盟成本c123 % 输出：三个主体的Shapley分配值（分配的是合作收益，正值为收益） % 说明：成本越低越好，收益特征 v(S) = -cost(S) v1 = -c1; v2 = -c2; v3 = -c3; v12 = -c12; v13 = -c13; v23 = -c23; v123 = -c123; % 对n=3，权重分别为 1/3（单元素边际）、1/6（两元素边际）、1/3（三元素边际） % 主体1的贡献 = (1/3)*(v({1})-0) + (1/6)*(v({1,2})-v({2})) % + (1/6)*(v({1,3})-v({3})) + (1/3)*(v({1,2,3})-v({2,3})) phi1 = (1/3)*(v1) + (1/6)*(v12 - v2) + (1/6)*(v13 - v3) + (1/3)*(v123 - v23); phi2 = (1/3)*(v2) + (1/6)*(v12 - v1) + (1/6)*(v23 - v3) + (1/3)*(v123 - v13); phi3 = (1/3)*(v3) + (1/6)*(v13 - v1) + (1/6)*(v23 - v2) + (1/3)*(v123 - v12); % 收益分配值通常是正数；若需要转成成本分担，用总成本减去phi即可 end

代码里的权重为什么要分1/3和1/6？因为对三个主体的Shapley值，边际贡献需要按所有可能的排列情况取平均。共3!种排列（ABC、ACB、BAC、BCA、CAB、CBA），其中单主体出现在第一位的情况有2!种，占比2/6=1/3；单主体出现在第二位、插入到两元素联盟的情况有1!×1!/6=1/6；单主体出现在第三位、加入已成形的两元素联盟的情况2!/6=1/3。所以公式里你看到单元素和两元素联盟的权重全是1/6，只有“单独启动”和“最后收尾”这两种情况的权重才是1/3。这个手算逻辑跟资源中shapley.m的实现应当一致，你可以用相同数据对比验证。

5. 避坑指南：IES合作博弈程序最常见的五个坑

5.1 报错“未定义函数或变量sc1”，根源是工作目录不对

现象：运行main123.m，MATLAB报错“未定义函数或变量'sc1'”，但你的确在程序包里看到了sc1.m。

原因：脚本所在目录不在MATLAB搜索路径中。现在很多同学把压缩包解压到桌面某子文件夹，然后直接在MATLAB的命令窗口敲main123，因为当前目录是“当前文件夹”默认位置，根本找不到脚本。

解决：先使用cd命令切换到程序目录，再用addpath把该目录加入路径。具体操作：

% 切换到你的实际路径 cd('D:\IES_Shapley\182-基于合作博弈的综合能源系统利益分配优化调度'); addpath(genpath(pwd)); % 将当前目录及其子目录都加入搜索路径

用genpath(pwd)可以把所有子文件夹也加进来，避免后续再遇到子目录内的函数找不到。从那以后我每次拿到新资源第一件事就是跑一下上节的检查脚本，确认所有文件都在路径里再开始。

5.2 求解器报错：ab.lp文件无法被MATLAB直接读取

现象：运行到涉及ab.lp的部分时，报错“无法打开文件ab.lp”或“未找到求解器”。

原因：MATLAB原生函数并不直接支持读取.lp格式的线性规划文件。ab.lp是给lpsolve或CPLEX这类外部求解器准备的，如果你没装对应工具箱，MATLAB解析不了这个文件，自然卡死。

解决：先看程序中是怎么调用求解器的。如果代码里写着cplex.lp或cplex.readModel，则需要安装CPLEX；如果用的是lpsolve，则需要下载lpsolve的MATLAB接口并注册路径。急着跑通算例的话，更推荐的做法是把ab.lp中的目标函数和约束改写成linprog能接受的格式：

% 用linprog替代LP文件读取的示例 f = zeros(1, 24); % 把24时段的决策变量代入目标系数 A_ineq = []; b_ineq = []; A_eq = []; b_eq = []; lb = zeros(1, 24); ub = ones(1, 24) * 500; x0 = ones(size(f))*10; [x, fval] = linprog(f, A_ineq, b_ineq, A_eq, b_eq, lb, ub);

这里fval就是联盟的最小成本，对应sc12等输出的c12。你可能需要从.lp文件里把目标函数系数和约束矩阵人工抄出来，虽然麻烦，但比装求解器快得多。

5.3 data.xlsx读入后全是NaN或矩阵尺寸不匹配

现象：用readtable读取后，数据列全是NaN，或者程序报“矩阵维度必须一致”。

原因：Excel文件里常见合并单元格、表头占多行、单位行夹在数据中、数字被保存成文本带单位，比如“1000元”而不是1000。readtable读进来时，这些非数值行会被置为NaN，导致后续矩阵拼接错位。

解决：把数据区复制到新工作表，只保留纯数值。我一般会先新建一个“clean”工作表，然后选中数据列，使用“数据→分列→分隔符号”强制转换列格式为数字。最终检查sum(sum(isnan(T)))是否为零，再进入正式运行。这里有一个重要提示：修改data.xlsx时不要增加行数，否则24时段的调度模型的维度会不一致。

提示：如果原程序在数据读取部分用了xlsread而不是readtable，旧版xlsread遇到非数值区会直接报错，遇到这种情况优先把Excel另存为xlsx并删除所有表头。

5.4 联盟成本传入顺序错位，导致Shapley收益出现负值

现象：shapley计算后，某个主体的分配收益为负数，或者三个主体分配收益之和明显不等于合作总剩余。

原因：这通常不是算法问题，而是传入shapley.m的联盟成本顺序与主体编号不对应。例如sc12.m计算的可能是“主体1+主体2”，但你在调用时把结果放在了c13的位置，公式中v({1,3})-v({3})算出来的就不是主体1的真实边际贡献了。我自己拿到这个程序时就踩过这个坑，因为sc12和sc12.m的文件名重复，差点把变量名搞混。

解决：在运行前打印每个联盟的输出值，并手算核对。一个快速的验证方法是：让所有设备参数都设置为“独立运行等于联盟运行”的极端情况，即c12 = c1 + c2、c13 = c1 + c3、c23 = c2 + c3、c123 = c1 + c2 + c3，这时Shapley收益应该全为0。如果跑出来非0，说明输入顺序错了。

5.5 复现论文结果有偏差：关键参数不一致

现象：程序能跑通，但输出的总成本曲线和论文对不上，数值差20%以上，不过趋势相近。

原因：论文算例里的碳交易价格、P2G效率、碳捕集率、热储能容量等参数，很可能与原程序中的默认值不同。资源里data.xlsx给出的是一组参考值，未必是与你手头论文完全匹配的算例值。

解决：翻论文参数表，逐项核对data.xlsx中的列。尤其注意碳交易价格的单位，有的是元/t，有的是元/kg，换算差出1000倍。把参数改成论文值后再运行，曲线基本能重合。这条经验在后续章节会展开。

6. 把案例改成你自己的数据：从修改参数到验证结果

6.1 修改data.xlsx中的负荷与风电数据

实际工程或毕业设计中最常见的诉求不是跑通现成算例，而是代入自己所在区域或课题给定的负荷数据。原程序data.xlsx里的数据是内蒙古地区某算例，你可以直接覆盖对应列。操作步骤里最关键的是保持时间列的行数不变，例如原文件是24小时调度，你的数据也必须是24行；如果是96点，则需要把程序里跟时间维度相关的循环变量全部改成96，否则各种小优化子问题都会维度爆掉。

% 按列名修改数据，并另存为新文件 T = readtable('data.xlsx'); % 假设你的数据已放在myData.xlsx中，三列分别为电负荷、热负荷、风电 new_data = xlsread('myData.xlsx'); T.elec_load = new_data(:,1); T.heat_load = new_data(:,2); T.wind = new_data(:,3); writetable(T, 'data_modified.xlsx');

这里的xlsread要求myData.xlsx中只有数值，不带单位。写入之后，还要回到main123.m中把读取文件名改成data_modified.xlsx，否则程序还在用旧数据。如果你想保留原始文件作为对照，可以不改原文件，而是复制一份到工作目录下再修改。

6.2 调整碳交易价格、P2G效率和碳捕集率

碳交易价格和P2G效率是影响Shapley分配结果最敏感的五个参数中的两个。碳价提高时，碳捕集设备的收益增大，含碳捕集的系统联盟更愿意合作；P2G效率升高时，电转气的边际收益变大，风电消纳更加积极，电-气联盟的成本下降也更明显。以下是我常用的参数调整参考区间：

修改位置一般在main123.m或sc系列脚本的开头参数区。找一个名为params的结构体，里面应该有carbon_price、eta_p2g、capture_rate等字段。注意单位统一：如果data.xlsx中碳排放成本已经给出，就不要再在代码里重复乘以系数，否则会翻倍。

6.3 复现论文中的系统总成本曲线并验证Shapley值

改完参数后，最后一步是验证结果是否与论文趋势一致。通常在main123.m运行完成后，会有一个变量保存各时段的系统总成本，可能是cost_series或者total_cost。把它画出来与论文曲线对比，重点关注波峰波谷的时段和量级。只要趋势一致、数值差异在10%以内，就说明你的模型和参数复现是成功的。

% 绘图验证总成本曲线 t = 1:length(cost_series); figure('Position',[100 100 800 400]); plot(t, cost_series, 'b-o', 'LineWidth', 1.2); hold on; % 如果论文里有对比曲线，自己先存为数组 paper_cost，然后叠加 % plot(t, paper_cost, 'r--s', 'LineWidth', 1.2); xlabel('时段 / h'); ylabel('系统总成本 / 元'); legend('本文模型', '论文算例', 'Location', 'best'); grid on; title('综合能源系统合作博弈调度总成本');

这段代码里，cost_series从main123.m中读取，如果原程序没有公开这个变量，你可以在main123.m里找到cost累加的位置，手动加一行“cost_series = …;”。验证完总成本，还要验证Shapley分配之和等于合作剩余（独立成本总和减总联盟成本）。我用这个公式作为所有分配计算正确性的唯一判据：sum(phi) == c1 + c2 + c3 - c123，差值小于1e-6才算通过。从那以后，我每次拿到这类合作博弈源码都会先把data.xlsx里的主体编号和sc12/sc13/sc23的拼写对应关系写在便签上，再跑main123——这个动作至少帮我省掉半天排查成本。希望帮到你。

本文还有配套的精品资源，点击获取