打开PyCharm新建一个空白工程合作博弈的数学工具箱哐当一声砸在桌面上——这个画面我现在回想起来还挺有仪式感的。综合能源系统的利益分配问题确实像块硬骨头单靠传统比例分摊或者拍脑袋谈判最后大概率谈崩。咱们今天要做的就是拿Shapley值这把牙口把这块骨头啃碎了、咽下去还得把每个参与方的账算得明明白白。这篇文章是我自己从问题建模、数学推导到PyCharm编码实测的全过程记录。适合三类人看一是做综合能源系统规划、微电网运营的工程师天天被钱怎么分逼到头大二是学运筹优化、博弈论的研究生想看看Shapley值在工程场景里到底怎么落地三是刚接触PyCharm和Python的初学者想找一个不那么玩具级的实战案例练手。文章里不会堆教科书公式我会把每个关键步骤都摊开讲包括我踩过的坑和返工的经历。1. 先把桌子和工具收拾好PyCharm环境与工程搭建1.1 为什么选PyCharm来做这个计算我见过太多人用记事本写Python跑出结果就算交差了。但一旦数据量上来、调试次数变多没有一个像样的IDE真的会让人抓狂。PyCharm是我最常用的Python开发环境选它来做Shapley值计算有几个很实际的理由。首先PyCharm对工程结构的组织能力很强。咱们这个项目看着只是算一个公式实际上涉及特征函数定义、子集生成、边际贡献计算、结果可视化好几个模块。直接把所有代码堆在单个.py文件里能跑但等到你要改某一个联盟的收益数值、加一个新的参与方、或者换个场景复用时就会知道模块化工程有多重要。PyCharm的Project工具窗口能让你一眼看清整个项目的文件结构这在迭代调试时省太多时间了。其次PyCharm的调试器堪称神器。Shapley值的计算过程有大量循环嵌套和字典查询一旦结果不对你得知道是哪一步出了问题。PyCharm里打上断点逐行查看每个子集的价值、每个边际贡献的计算结果问题几分钟就能定位。这一点对理解Shapley值的内在逻辑尤其有帮助——你不再是拿着黑盒跑结果而是亲眼看着每个数值是怎么一步步累加出来的。最后PyCharm对科学计算生态的支持很完备。numpy、pandas、matplotlib这些库的安装、导入、智能提示都做得顺滑虚拟环境的管理也直观。咱们这个项目虽然不需要重型科学计算但用pandas整理特征函数表、用matplotlib画分配结果对比图体验比在命令行里折腾好太多。1.2 新建工程与配置解释器打开PyCharm后选择New Project建议先别急着把文件建在默认路径而是建一个语义清晰的目录比如shapley_ies。工程名和文件命名尽量用英文避免某些库在中文路径下出现莫名其妙的编码问题这个教训我在后文会细说。新建工程时有几个关键选项需要注意Location工程存放路径建议专门建一个projects文件夹别散落在桌面和下载目录里。Interpreter解释器类型。如果你之前装了Anaconda可以直接选择已有的Conda环境如果没有选择Virtualenv或者直接用系统Python都行。对咱们这种中等规模的计算三者没有本质区别关键是别选错环境否则后面装库全装到别的环境里去了运行时又报ModuleNotFoundError能让人怀疑人生。Create a main.py welcome script默认勾选可以去掉咱们不需要那个示例脚本。工程创建完成后第一时间检查解释器是否正常显示。在PyCharm右下角状态栏可以看到当前解释器路径点击可以切换。如果你用的是Conda环境建议顺手装一下conda插件并在Settings里把Conda环境配好这样PyCharm能正确识别已安装的包代码里写import numpy时不会有红色波浪线。1.3 装依赖库numpy、pandas、matplotlib一个都不能少在PyCharm中安装第三方库有两种方式我个人的习惯是在Terminal里直接操作因为能顺便看到安装日志出了问题也好排查。pip install numpy pandas matplotlib如果你的网络环境较差导致下载超时可以临时换用国内镜像源pip install numpy pandas matplotlib -i https://pypi.tuna.tsinghua.edu.cn/simple安装完成后在代码里依次导入确认无报错即可import numpy as np import pandas as pd import matplotlib.pyplot as plt这里有几个小细节值得一提。第一matplotlib在中文Windows上默认字体显示不了中文画图时坐标轴和标签全是小方框。解决办法是显式指定支持中文的字体比如SimHei或Microsoft YaHei。第二如果你装的PyCharm版本较新创建工程时可能默认启用了PEP 8检查代码风格问题会以黄色波浪线提示不影响运行但建议顺手规范一下毕竟代码是要给别人看的。2. 拆解问题综合能源系统为什么需要合作博弈2.1 多主体综合能源系统的利益纠葛综合能源系统说人话就是光伏、风电、燃气轮机、储能、负荷等多类能源设备和用户组合在一起通过协调运行实现多能互补、提质增效的系统。理想很丰满现实很骨感——因为参与方往往不是同一个业主。我做过的一个实际测算场景里园区综合能源系统里有四类典型主体光伏电站、风电场、燃气轮机和储能电站。它们分属不同投资方有各自独立的运营账本。光伏电站关心年利用小时数和补贴收益风电场担心弃风限电燃气轮机盯着燃料成本和启停损耗储能则纠结于峰谷价差能不能覆盖电池衰减成本。关键问题在于这四类主体单独运营时效益都算不上好。光伏存在弃光问题白天出力过剩送不出去晚上又彻底歇菜风电的出力波动性强电网调度经常要求限出力燃气轮机单独支撑负荷时得频繁调整出力效率低、燃料成本高储能最惨单纯靠低买高卖峰谷价差赚的那点钱连电池维护费都悬。但是如果它们联合起来情况就完全不同了。光伏白天多发电储能把多余电量存下来晚上储能放电支撑晚高峰燃气轮机只在光伏和风电都出力不足时启动发挥快速调节和兜底作用。风电出力波动时储能和平稳运行的燃气轮机共同补偿。整条链条一配合弃光弃风率大幅下降燃气轮机不用频繁变负荷储能利用率也上去了。问题来了总收益确实增加了但这增加的收益该如何分配如果分得不公平光伏觉得我贡献了最多的电量储能觉得我帮大家消纳了那么多弃电燃气轮机说我兜底出力关键时刻全靠我风电说没有我你们晚上那缺口你补吗——各方都觉得自己亏了合作第二天就告吹。这种多方投入、联合产出、收益需要科学分摊的场景博弈论中的合作博弈理论就是专门解决这个问题的。2.2 为什么不能各干各的互补性分析在做正式建模之前有必要把合作到底增加了多少价值这件事量化清楚。还是用上面的四主体场景假设独立运营时各主体一个月能获得的净收益分别是单位万元主体独立运营净收益光伏40风电35燃气轮机60储能15单独加起来一共是150万元。但如果两两合作或者三者合作经过协调运行收益会比各自单干之和更大。比如光伏加储能联合运营储能能消纳部分弃光发电量光伏不必低价甚至免费送电储能也能获得更多充放电循环一个月总收益可以做到62万元比两个主体独立运营之和55万元多出7万元。这种112的性质在博弈论里叫超可加性。综合能源系统天然具备这种特性因为不同能源品种在时间尺度和出力特性上存在很强的互补性。光伏和风电可以在一定程度上互补出力曲线燃气轮机和储能可以配合进行调峰调频几方联合后整体对外供电的可靠性和经济性都提升了。这正是合作博弈能够成立的前提——如果合作不能带来增量价值那么任何分配方案都是空谈。2.3 为什么选Shapley值而不是简单按比例分摊很多人的第一反应是既然合作之后总收益增加了那就按投资比例分呗或者按独立收益占比分。这样做问题很大。按投资比例分储能往往会吃大亏。储能设备单位投资高、回收周期长独立运营时看似收益最低但在联合运营中它承担的削峰填谷、消纳弃电的角色极其重要。如果仅按投资比例分配储能方案成了冤大头下次没人愿意再投资储能了。按独立收益占比分结果类似本质是存量分配而非增量分配完全没有体现某个主体加入联盟后带来多少额外价值。光伏独立收益高可能分到很大一块但真正让总体收益从150涨到180的边际贡献有一部分是储能提供的。忽略了边际贡献的分配方案对参与合作的积极性和工程项目的长期可靠性都是极大的伤害。Shapley值之所以在合作博弈里地位极高是因为它严格满足了四个公理有效性、对称性、虚拟参与人性质、可加性。简单理解就是分配结果的总和必须等于合作总收益编号互换不影响分配对合作没有任何贡献的人分不到钱多个独立博弈合并后分配自动对应相加。这四个公理保证了分配结果在数学上是公平的让每一个参与方拿到的份额精确反映它对各联盟产出的边际贡献。这在多主体能源系统里太重要了——因为没人愿意接受一个你说公平但我说不公平的方案。3. Shapley值的数学原理与实践建模3.1 Shapley值到底在算什么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]$$看着吓人拆开就清楚了。$N$是所有参与人的集合$n$是参与人总数$S$是不包含参与人$i$的任意子联盟$v(S)$是子联盟$S$独自合作能创造的价值$v(S \cup {i})$是$i$加入$S$之后创造的价值两者的差就是$i$加入$S$时的边际贡献。前面的系数$\frac{|S|! (n - |S| - 1)!}{n!}$本质上是所有加入顺序中能形成$S$在$i$之前加入这种情形的排列数占总排列数的比例。换成生活化类比几个合伙人开公司不同合伙人加入公司的先后顺序不一样每个人加入前后公司估值的差值也不一样。Shapley值就是把这所有128种如果有7个合伙人就是7!种可能的加入顺序下的估值差取平均算出一个最能反映此人真实价值的股份数。在n4的例子里有4! 24种不同的排列顺序。逐一去枚举排列再算平均是理解这一步最直观的方式代码实现也很顺手。后面我用代码实现的也是这条路径。3.2 特征函数设计联盟收益怎么来特征函数$v(S)$是整个Shapley值计算的核心输入它规定了任意一个玩家集合$S$通过内部合作能够创造的总收益。这个函数不是凭空拍出来的在综合能源系统项目里它通常由两部分工作得到第一步是场景仿真。针对每一个可能的联盟组合用电力系统优化调度模型比如混合整数线性规划模拟该联盟成员联合运行时的最优经营策略求解出该联盟一个周期内的最优总收益。这里面要包含设备出力约束、储能SOC约束、供需平衡约束、电价曲线等细节。我用的这套数值虽然是我为了演示手算过程构造的示例数据但真实项目中每一个$v(S)$都应当来自这类仿真计算。第二步是对数值做合理化校验。正常的特征函数必须满足两条性质第一空联盟收益为0也就是$v(\varnothing)0$第二超可加性即对任意两个不相交的联盟$S$和$T$必须保证$v(S \cup T) \geq v(S) v(T)$。否则合作反而亏钱这个博弈就没有维持合作的必要了。我在文章里构造了一套满足上述性质的演示数据覆盖全部16个非空子集四类主体两两组合一共6个三三组合4个四主体联盟1个单主体4个加上空集。这套数据对应表里列出的就是特征函数字典的原始依据。联盟联盟收益万元{光伏}40{风电}35{燃气轮机}60{储能}15{光伏, 风电}80{光伏, 燃气轮机}105{光伏, 储能}62{风电, 燃气轮机}100{风电, 储能}58{燃气轮机, 储能}82{光伏, 风电, 燃气轮机}150{光伏, 风电, 储能}102{光伏, 燃气轮机, 储能}132{风电, 燃气轮机, 储能}125{光伏, 风电, 燃气轮机, 储能}180你可以看到总联盟收益180万元比四个主体独立收益之和150万元多出30万元这30万元就是协同运营创造的增量价值也是Shapley值要分配的红利。3.3 手算一遍Shapley值数据让我心服口服在做代码实现之前我先手动把光伏的Shapley值算了一遍这一步对理解算法逻辑非常有帮助也便于后续验证代码结果是否正确。光伏的边际贡献需要考察它加入各种不包含它的联盟时的增量。结果如下联盟Sv(S)v(S∪{光伏})边际贡献权重加权值∅040401/410.00{风电}3580451/123.75{燃气轮机}60105451/123.75{储能}1562471/123.92{风电, 燃气轮机}100150501/124.17{风电, 储能}58102441/123.67{燃气轮机, 储能}82132501/124.17{风电, 燃气轮机, 储能}125180551/413.75四类主体时$|S|0$和$|S|3$的权重都是$1/4$$|S|1$和$|S|2$的权重都是$1/12$逐项累加后光伏的Shapley值约为47.17万元。对比独立运营的40万元光伏通过合作多拿了7.17万元。按同样的方法把其余三个主体的Shapley值也算出来最终结果主体独立收益Shapley分配增收增幅光伏4047.177.1717.9%风电3541.676.6719.0%燃气轮机6068.178.1713.6%储能1523.008.0053.3%合计150180.0030.00—储能分配结果的增幅最大达到了53.3%。这个结果很耐人寻味储能独立运营时最不赚钱但在合作体系里它通过消纳弃电、削峰填谷带来的边际贡献很大理应在增量收益中分到大头。Shapley值把这个隐性价值显性化了。四者分配总和正好等于180万元一分钟验算无误。4. 用PyCharm写代码从特征函数到分配结果4.1 数据准备与特征函数字典在PyCharm工程里新建一个shapley_ies.py文件第一步是把特征函数写成Python字典的形式。为了让代码可读性高我用玩家名称字符串拼接后的排序结果作为字典键。比如联盟光伏、风电记作PVWT联盟风电、储能记作WTBESS。from itertools import permutations from collections import defaultdict PLAYERS [PV, WT, GT, BESS] V { PV: 40, WT: 35, GT: 60, BESS: 15, PVWT: 80, PVGT: 105, PVBESS: 62, WTGT: 100, WTBESS: 58, GTBESS: 82, PVWTGT: 150, PVWTBESS: 102, PVGTBESS: 132, WTGTBESS: 125, PVWTGTBESS: 180, } def coalition_value(coalition): if not coalition: return 0 key .join(sorted(coalition)) return V[key]这里的coalition_value函数是特征函数的统一入口空联盟返回0非空联盟通过拼接排序后的成员名去字典里取值。这种写法在后面计算边际贡献时非常顺手不会因为集合的哈希性问题引发麻烦。4.2 核心算法实现用全排列把逻辑讲到透Shapley值实现方式有两种一种是带权重公式计算另一种是全排列统计平均。公式方法执行效率高但对新手不太友好尤其是在理解权重系数出处时容易卡壳。全排列方法更贴近定义本身——把$n!$种玩家加入顺序全部枚举出来对某个玩家来说每种顺序下它加入时的边际贡献求和再除以$n!$就是它的Shapley值。我在这篇文章里选用全排列实现逻辑清楚、不易出错而且n4时只有24种排列计算量可以忽略不计。def shapley_by_permutations(players, value_func): n len(players) phi defaultdict(float) count 0 for perm in permutations(players): coalition set() for p in perm: marginal value_func(coalition | {p}) - value_func(coalition) phi[p] marginal coalition.add(p) count 1 for p in players: phi[p] / count return dict(phi)这段代码里最核心的操作是value_func(coalition | {p}) - value_func(coalition)它计算的就是玩家p加入当前联盟前后的收益差。跑完24种排列后每个玩家把每一次的边际贡献加起来再除以24得到的就是Shapley值。如果玩家数量增多比如到了8个、10个全排列$n!$会爆炸式增长。那时候可以改用基于子集的加权公式或者使用蒙特卡洛抽样近似计算。我对代码做了个扩展同样用公式法实现了一遍from itertools import combinations from math import factorial def shapley_by_formula(players, value_func): n len(players) phi defaultdict(float) for i in players: others set(players) - {i} for k in range(n): for S_set in combinations(others, k): S_set set(S_set) weight factorial(k) * factorial(n - k - 1) / factorial(n) marginal value_func(S_set | {i}) - value_func(S_set) phi[i] weight * marginal return dict(phi)两种方法的结果应当完全一致可以交叉验证。我在实际跑代码时验证过输出结果是47.17、41.67、68.17、23.00和手算结果一分不差。4.3 可视化与结果分析算完Shapley值如果只打印一串数字既不好看也不利于向各方汇报。用matplotlib画一张柱状图把独立收益和Shapley分配值放在一起对比增量部分用明显的颜色标出。别忘记配置中文字体不然图上全是方框。import matplotlib.pyplot as plt plt.rcParams[font.sans-serif] [SimHei] plt.rcParams[axes.unicode_minus] False results shapley_by_permutations(PLAYERS, coalition_value) names [光伏, 风电, 燃气轮机, 储能] ind [40, 35, 60, 15] shap [results[p] for p in PLAYERS] x range(len(names)) plt.figure(figsize(10, 6)) bars1 plt.bar(x, ind, width0.35, label独立运营收益, color#8ecae6) bars2 plt.bar([i 0.35 for i in x], shap, width0.35, labelShapley分配, color#ffb703) for bar in bars1: plt.text(bar.get_x() bar.get_width() / 2, bar.get_height() 0.8, f{bar.get_height():.1f}, hacenter, fontsize10) for bar in bars2: plt.text(bar.get_x() bar.get_width() / 2, bar.get_height() 0.8, f{bar.get_height():.2f}, hacenter, fontsize10) plt.xticks([i 0.175 for i in x], names, fontsize12) plt.ylabel(收益万元) plt.title(综合能源系统合作收益分配独立运营 vs Shapley值) plt.legend() plt.grid(axisy, linestyle--, alpha0.3) plt.tight_layout() plt.savefig(shapley_result.png, dpi200) plt.show()画完图你会有更直观的感受储能主体的两栏柱子高度差肉眼可见它从15万涨到23万涨幅接近一半燃气轮机虽然拿到的绝对值最高但独立运营时它本来就是60万的高基数涨幅反而最小。把这段可视化代码在PyCharm里运行后输出文件会保存在工程根目录下。需要提醒的是如果电脑上装了多个Python环境PyCharm的Run窗口和环境变量偶尔会出现编码或路径问题报错No module named matplotlib时先检查右下角的解释器是不是你安装matplotlib的那一个别急着重装。5. 实操中的几个坑与排查建议5.1 特征函数设计的常见错误这个坑是所有做Shapley值的人必踩的。特征函数定义错了后面一切计算都白费。最常见的错误有两种。第一种是忽略空联盟的收益必须为0。有些人在特征函数里把某个主体的独立收益写成了包含固定成本扣除后的净利润导致空集处理不当计算结果出现系统性偏差。处理方式很简单——在我的coalition_value函数里第一行就判断了空集并返回0。第二种是数据不满足超可加性。比如某两个主体合作后总收益反而比它们单独运营之和还低这种情况下Shapley值仍然能算出结果但这个结果没有实际指导意义——因为理性主体根本不乐意合作。检查方法很粗暴但有效把所有特征函数值打印出来逐一验证$v(S \cup T) \geq v(S) v(T)$对所有不相交的S和T成立。我在项目里写了一个自动检查的小函数很实用。def check_superadditivity(players, value_func): n len(players) for mask in range(1, 1 n): S {players[i] for i in range(n) if mask (1 i)} rest set(players) - S for sub_mask in range(1, 1 len(rest)): rest_list list(rest) T {rest_list[i] for i in range(len(rest)) if sub_mask (1 i)} if value_func(S | T) value_func(S) value_func(T): print(f不满足超可加性: S{S}, T{T})5.2 计算复杂度问题玩家一多就跑不动Shapley值最大的软肋是计算复杂度。$n$个玩家需要遍历$2^n$个子集用全排列实现则需要$n!$次操作。我项目里只有4个玩家随便跑但如果你要给一个包含10个以上主体的综合能源系统算分配直接暴力枚举就崩溃了。我给两个工程中切实可用的优化手段。第一改用子集加权公式把全排列的$n!$降到$2^n$——还是指数级但n10时至少有实用价值。第二采用蒙特卡洛近似随机抽取大量排列计算边际贡献的平均值样本数量越大越接近精确Shapley值。工程实践中通常抽几千到几万次就能得到足够精度。import random def shapley_monte_carlo(players, value_func, n_samples10000): phi defaultdict(float) n len(players) for _ in range(n_samples): perm players[:] random.shuffle(perm) coalition set() for p in perm: marginal value_func(coalition | {p}) - value_func(coalition) phi[p] marginal coalition.add(p) return {p: v / n_samples for p, v in phi.items()}5.3 稳定性检查分配结果是否被接受Shapley值算出结果不意味着各方一定心服口服。合作博弈理论里有个核的概念——如果某种分配方案下某个子联盟能获得的收益比它在总联盟分配中拿到的还多这个子联盟就有动机另起炉灶所谓分配方案就是不稳定的。Shapley值并不总是落在核里面这是它的一个天然局限。我在实际项目中会额外做一个稳定性验证。针对每个真子集$S$检查分配给$S$中成员的Shapley值之和是否不少于$v(S)$。如果某个子集不满足说明存在分裂动机。这时候通常的补救措施有两种一是微调分配方案使其落入最小核least core范围内二是设计补偿机制。比如前面的结果中如果我检查出风电加储能的联盟收益58万元高于两方Shapley值之和64.67万元我算一下风电41.67 储能23.00 64.67大于58所以稳定。所有子联盟我都验证过没有发现问题。但这个检查步骤不能省尤其是数据规模增大以后。5.4 PyCharm运行时的几个常见报错最后说说代码运行时容易遇到的几个问题。ModuleNotFoundError是最常见的。原因无非两种解释器选错环境、包没装进当前环境。解决方法是在PyCharm的Terminal里执行pip list确认包是否存在再看右下角解释器路径是否属于同一个环境。SyntaxError或者奇怪的缩进报错通常是代码从网页或文档拷贝时混入了中文全角空格。遇到这类问题先在PyCharm的Settings里开启显式空白符显示然后把出错的代码行前后空格全部删除重打一遍。matplotlib中文标题显示为方框的问题我在前文已经提到在代码开头加上两行中文字体配置即可。如果你用的不是Windows而是macOS或者Linux中文字体路径可能不同可以用plt.rcParams[font.sans-serif] [PingFang SC]或者查看系统可用字体列表后选择合适的字体名称替换。数据精度问题也值得提一句。Shapley值计算涉及浮点数累加多次循环后可能出现微小误差我在输出结果时保留两位小数即可满足工程精度。如果对精度要求极高可以改用Python的decimal模块。结尾一点个人体会这套PyCharm加Shapley值的方案我现在回头复盘觉得最值得推荐的倒不是某个具体算法或者某个IDE功能而是那种把问题建模清楚、把逻辑拆解到底的工作方式。Shapley值的好处在于它的可解释性极强——每个人都拿到了自己那份边际贡献的加权平均谁也挑不出毛病。放在综合能源系统的商业谈判里这种能讲清楚道明的分配机制比拍脑袋分账靠谱得多。而PyCharm在这件事里提供的工程化支持让从仿真数据到最终可视化的整条链路都能顺畅跑通。如果你之后要处理比这四个主体更复杂的工程场景比如把电、热、气等多种能源耦合在一起考虑或者涉及数十个分布式主体的合作运营我建议从特征函数构建就开始考虑模块化设计最好把特征函数的计算和Shapley值求解封装成独立的函数方便复用。另外合作博弈的工具箱里其实还有核仁、班茨哈夫指数等分配方法不同方法在公平性定义上各有侧重工程上可以并行算几套方案供决策层选。至少从我这几次实操来看Shapley值作为首选的公平基准线在综合能源系统的多主体利益分配问题里是经得起推敲的。