企业官网建设流程全解析

简介：这是一份面向综合能源系统（IES）优化调度方向的毕业设计源程序包，对应知网论文《基于合作博弈的综合能源系统利益分配优化调度》。内容围绕“双碳”背景下IES低碳经济运行问题，构建包含电转气、碳捕集、燃气轮机与热储能等设备的合作博弈模型，并通过Shapley值法完成合作剩余分配，适合电力系统、新能源或能源经济方向的本科生与研究生参考复现。压缩包共13个文件，以10个MATLAB脚本（.m）为主，分别承担模型构建、优化计算与结果分析等功能；另含1个图片文件用于展示仿真结果、1个Excel数据文件存放算例输入参数，以及1个.lp文件描述优化问题，整体大小约215KB。资源附带数据文件与结果图，可直接在MATLAB中运行，便于对照论文验证调度结果与利益分配逻辑。目前已有231人学习下载，对理解多主体协同优化与碳交易机制下的IES运行策略具有较高参考价值。

1. 综合能源系统的合作博弈调度源程序：先别急着写代码，想清楚“分钱逻辑”再动手

做综合能源系统优化调度的人，多半会撞上同一个墙：模型建出来了、Gurobi 也跑通了、联合调度的总成本确实降了，但问一句“省下来的钱怎么分”，整篇论文就卡住了。很多做这个方向的源程序，其实是把“合作博弈利益分配”和“优化调度”两件事焊在了一起——先让多个能源主体组成联盟一起运行，再用 Shapley 值把联盟总收益切成每个主体的贡献。你看到的标题《基于合作博弈的综合能源系统利益分配优化调度》，说的正是这套完整闭环：论文全文在知网按标题直接能检索到，本人博客里还有源码解读。这篇我按自己复现这类代码的路线，把它拆成建模层、调度层、分配层三层来讲，适合正在复现论文、做毕设或者准备把博弈调度写成期刊论文的人。

先说清楚一个基本前提：合作博弈在这里不是耍花架子。综合能源系统里通常有好几个独立运营商（燃气轮机方、电锅炉方、储能方），各自独立运行也能活，但一旦联合起来共用设备、互济功率，总运行成本就会下降。这个下降空间来自热电联产余热利用、储能共享、购买力的折扣与峰谷套利。联盟总收益是实打实的，但“谁贡献大、谁该多得”没有天然答案。源程序的价值就是要让这个答案可计算、可复现、可论证——也正是在这一步，大量仿真代码会翻车。因为收益分配不只看优化调度的结果，还要看你有没有把特征函数、联盟枚举和 Shapley 权重写对。下面我就按这三个层次，把每一步的落地细节和参数盘清楚。

2. 建模层：先把“不合作时各自的钱”算明白

2.1 综合能源系统的模型结构与设备参数表

复现这种源代码，第一步不是直接上博弈，而是先把物理模型坐实。常见做法是设置三个能源主体共用一套园区电热负荷：主体 A 是燃气轮机（CHP）运营商，主体 B 是电锅炉运营商，主体 C 是储能运营商。这样设置的好处是：主体数量少，联盟组合只有 7 个，Shapley 值的幂集展开能手动验算；同时每个主体的设备特征差异明显，合作收益来源清晰。

模型参数一般包括三类：设备容量与效率、能源价格、负荷曲线。我在自己搭模型时常用下面这组参数做基准，单位统一成 MW 和 万元/MWh，避免后面换算收益时对不上账：

需要说明，这里的效率影响的是热电联产约束系数。燃气轮机发 1 MW 电，会有 0.9 MW 左右的热回收进热网（具体按热效率/电效率折算），电锅炉则直接吃电产热；储能只做充放电套利，不参与热平衡。给这三个主体配上典型日的电负荷和热负荷曲线（24 个时点），就能开始算独立运行。

负荷数据没有公开统一标准时，我一般用某园区典型日负荷的差值曲线，电负荷峰谷差控制在 40% 左右，热负荷集中在早晚两个峰。注意不要让热负荷高到电锅炉单独顶不住，否则独立运行模式下某个主体完全无法满足自身热负荷，那“不合作的底线收益”就成负数了，后面的谈判破裂点会失真。

2.2 独立运行模式：把“谈判破裂点”算出来

合作博弈里有个概念叫谈判破裂点（disagreement point），在综合能源系统里它就是每个主体不加入联盟、独自运行时的最优成本。Shapley 值分配出来的结果必须保证每个主体分到的收益不低于这个破裂点，否则联盟没有存在意义。源码里的第一步，就是把 M个 主体各自独立跑一次经济调度：

import gurobipy as gp from gurobipy import GRB def independent_operation(agent, load_e, load_h, price_e, price_g): """ agent: 主体字典，包含该主体拥有的设备参数 load_e / load_h: 该主体负责的电负荷、热负荷曲线（24h） price_e / price_g: 网购电价、购气价（万元/MWh） """ m = gp.Model("independent_%s" % agent["id"]) # 决策变量（以 CHP 主体为例，所有设备出力 + 外部购能） p_gt = m.addVars(24, lb=0, ub=agent.get("gt_cap", 0), name="gt_power") h_gt = m.addVars(24, lb=0, ub=agent.get("gt_heat_cap", 0), name="gt_heat") p_buy = m.addVars(24, lb=0, ub=200, name="buy_grid") g_buy = m.addVars(24, lb=0, name="buy_gas") # 目标：24h总运行成本最小 = 购电成本 + 购气成本 m.setObjective( gp.quicksum(price_e[t] * p_buy[t] + price_g[t] * g_buy[t] for t in range(24)), GRB.MINIMIZE ) # 电平衡约束：燃气轮机发电 + 网购电 = 电负荷 m.addConstrs((p_gt[t] + p_buy[t] == load_e[t] for t in range(24)), "elec_balance") # 热平衡约束：燃气轮机余热 = 热负荷 m.addConstrs((h_gt[t] == load_h[t] for t in range(24)), "heat_balance") # 热电联产耦合：余热产出与发电量成比例（这里是简化的定热电比） m.addConstrs((h_gt[t] == 0.9 * p_gt[t] for t in range(24)), "chp_coupling") m.optimize() return m.objVal

注意几个关键参数：决策变量的上界ub直接由主体拥有的设备容量决定——储能不能发电，电锅炉不能产热，源程序里每个主体独立运行时“可调用的设备”要和自身属性严格匹配。还要注意燃气轮机的热电耦合约束，这是 CHP 主体区别于其他主体的灵魂约束。独立运行模式下，如果你的代码里不限制 CHP 的热出力比例，它会变成一台既能卖电又能卖热的“完美机器”，那独立运行成本会低得离谱，联盟收益会被高估，分配结果立刻失真。

运行完独立模式后，保存每个主体的最优成本，这就是后续特征函数计算中的基准。通常你会发现储能主体的独立收益不高甚至为负——它只能靠峰谷价差套利，如果峰谷电价差不够，单纯储能的破裂点可能不理想。论文源码里这部分的处理办法是保留储能主体但给它一定的峰谷套利空间，或者让储能和光伏绑定成一个主体。我习惯直接用峰谷价差 0.3 元/kWh 的场景，这样储能主体勉强能活，后续 Shapley 值也不会因为某个成员“天生是负担”而出现极端分配。

3. 合作调度层：从“各自维护平衡”到“共享设备池”

3.1 联盟运行的目标函数与约束怎么改

让多个主体组成联盟后，要做的第一件事是把“各主体必须自己满足自己的负荷”这一条拆掉。独立运行时有 N 个电平衡方程、N 个热平衡方程，联盟运行时整个联盟只保留一个联合电平衡和一个联合热平衡——这就是合作调度层的核心特征来源。

目标函数从“单个主体的购电购气成本”变成“联盟的总购能成本”。比如 CHP 主体和电锅炉主体组成联盟，原本电锅炉要自己买电来产热，CHP 余热却能直接供给电锅炉的热负荷，联盟内部出现了一次能源替代：热负荷由燃气轮机的余热承担，电锅炉可以少出力甚至不出力，省下的电费就是联盟收益。如果不拆掉“各自平衡”的约束，这种互补根本不会发生。这是我在复现时最容易踩的坑：合作调度和独立调度的差别，不是目标函数加一项，而是平衡约束从“主体级”改成“联盟级”。

设备参数还是原来那套，变量含义不变。储能主体参与联盟时，充放电功率可以在联盟内任意时点使用，不再受“只能服务于自己负荷”的限制——这会让储能的利用价值大幅提升，也直接反映在 Shapley 的边际贡献上。编写联盟运行函数时，我在源代码里通常用一组member_ids列表作为入参，然后从全局设备池里筛出该联盟可用的设备：

def alliance_operation(member_ids, common_data): """ member_ids: 联盟成员编号列表，如 [0, 1] 表示主体 A 和 B 合作 common_data: 包含全部设备参数、负荷、电价的全局字典 """ m = gp.Model("alliance_" + "_".join(map(str, member_ids))) # 联盟可用设备由成员列表决定：成员越多，可用设备池越大 has_gt = 0 in member_ids # 燃气轮机是否可用 has_eb = 1 in member_ids # 电锅炉是否可用 has_storage = 2 in member_ids # 储能是否可用 p_buy = m.addVars(24, lb=0, ub=200, name="buy_grid") g_buy = m.addVars(24, lb=0, name="buy_gas") # 按成员身份条件声明设备变量 p_gt = m.addVars(24, lb=0, ub=120, name="gt_power") if has_gt else None p_eb = m.addVars(24, lb=0, ub=60, name="eb_power") if has_eb else None ... # 联盟总电负荷 = 各成员电负荷求和 load_e_total = [sum(common_data["load_e"][k][t] for k in member_ids) for t in range(24)] load_h_total = [sum(common_data["load_h"][k][t] for k in member_ids) for t in range(24)] # 目标：联盟总购能成本 m.setObjective(gp.quicksum(price_e[t] * p_buy[t] + price_g[t] * g_buy[t] for t in range(24)), GRB.MINIMIZE) # 联合电平衡 m.addConstrs(( (p_gt[t] if has_gt else 0) + (p_dis[t] if has_storage else 0) + p_buy[t] == load_e_total[t] + (p_char[t] if has_storage else 0) + (p_eb[t] if has_eb else 0) for t in range(24)), "joint_elec_balance") m.optimize() return m.objVal

这段代码的要点是：联盟收益不来自设备本身的效率提升，而是来自设备之间的互补。变量声明用if has_gt做条件控制，是为了模拟“该成员没加入时设备不存在”的场景。空联盟（只有一个主体）调用这个函数，结果应当和独立运行一致——这是验证合作调度代码是否正确的第一步测试。

3.2 联盟调度的求解器调用与参数设定

求解器方面，源码大多用 Gurobi 或 CPLEX 的 Python 接口写 MILP。如果你的机器没有 Gurobi 授权，用开源的CBC求解器也能跑，但要注意两件事：一是 MIP 求解速度会明显偏慢，联盟数量一多，CBC 可能要跑几分钟才出一个可行解；二是 CBC 对某些数值尺度敏感，建议把所有成本系数都放大到“元”而不是“万元”，避免求解器在容差范围内无法收敛。

求解器参数有两个是关键：MIPGap和TimeLimit。Shapley 值要求每个联盟的解都是全局最优的，因为特征函数值一旦有偏差，后续边际贡献的权重就会算错，最终分配结果的对账会对不上。所以我一般把 MIPGap 设为 0.001（即 0.1%），TimeLimit 设为 300 秒，不收敛就直接跳过错报，不让后续的分配步骤用次优解去算。

另外要做一次“可加性”测试：把两个主体分别独立求解的成本相加，和它们组成联盟后的总成本比较。如果联盟总成本不低于两个独立之和，说明这两个成员之间没有互补性，联盟收益为负，那这个联盟不该出现。在做 3 主体完整 Shapley 值时，你总会发现某些二元联盟的收益很小甚至接近零，这是正常的——正是这些低边际贡献的联盟，让你能看出谁才是真正的“关键角色”。

4. 利益分配层：用合作博弈的 Shapley 值算“谁贡献大”

4.1 特征函数先从调度结果里怎么切出来

合作博弈的核心是特征函数 v(S)，它定义每一个联盟 S 能创造多少“价值”。但我们在调度层得到的是成本，不是收益，所以要在这一层做一次关键换算：每个联盟的收益 = 联盟成员的独立运行成本之和 − 联盟联合运行成本。

这个换算不做，Shapley 值算出来就是负数，符号完全反了。这是整个源程序里最隐蔽的一个逻辑转换，很多复现者把代码跑通后，发现分配结果怎么是负的，就是因为少了这一步。空联盟的收益定义为 0，单成员联盟的收益也必须是 0——因为单成员“联盟”就是独立运行，成本和独立成本相等，收益为零。用一个 Python 字典来缓存所有联盟的收益，是避免后续 Shapley 计算重复调用求解器的关键：

def compute_profit_for_alliance(member_ids, independent_cost, common_data): """ 返回联盟 member_ids 的收益 independent_cost: 字典，key 是成员编号，value 是该成员独立运行成本 """ # 联盟运行成本 alliance_cost = alliance_operation(member_ids, common_data) # 独立成本之和 base_cost = sum(independent_cost[k] for k in member_ids) # 收益 = 省下来的钱 profit = base_cost - alliance_cost return profit # 缓存示例，避免同一个联盟反复求解 profit_cache = {} def get_profit(member_ids, independent_cost, common_data): key = tuple(sorted(member_ids)) if key not in profit_cache: profit_cache[key] = compute_profit_for_alliance(list(key), independent_cost, common_data) return profit_cache[key]

注意independent_cost必须是每个主体独立运行时的最优成本，不能拿联盟运行结果里的某个数值去填充。profit_cache这个字典用联盟成员的排序元组做键，在后面的幂集枚举中非常实用，否则三主体还好，六主体以上联盟组合几十个，每个都重新调求解器会非常慢。

4.2 Shapley 值的幂集展开与权重公式

Shapley 值的本质是边际贡献的加权平均。某个主体的分配值等于：它在所有不包含自己的联盟中，加入前后联盟收益变化量的加权和。权重公式看起来绕，但写成代码就几行：

from itertools import combinations from math import factorial def shapley_value(players, independent_cost, common_data): """ players: 总成员列表，如 [0, 1, 2] """ n = len(players) phi = {i: 0.0 for i in players} # 枚举所有不包含主体 i 的联盟 S for i in players: others = [p for p in players if p != i] for r in range(0, n): # S 的大小从 0 到 n-1 for S in combinations(others, r): S_key = tuple(sorted(S)) # 联盟 S 的收益 v_S = get_profit(list(S_key), independent_cost, common_data) if S_key else 0.0 # 加入 i 后的联盟收益 Si_key = tuple(sorted(S + (i,))) v_Si = get_profit(list(Si_key), independent_cost, common_data) marginal = v_Si - v_S # Shapley 权重：|S|! * (n - |S| - 1)! / n! s = len(S) weight = factorial(s) * factorial(n - s - 1) / factorial(n) phi[i] += weight * marginal return phi

整个计算的复杂度和 n 的关系是 O(n·2^n)。n=3 时有 7 个联盟，n=5 时有 31 个联盟，n=7 时有 127 个，每次都要调用 MILP 求解器。因此很多论文源码只做到 4~5 个主体是合理的工程取舍，不是模型能力不足，而是 Shapley 值枚举天然有这个指数瓶颈。

如果你要扩展到 10 个以上主体，常见做法是用蒙特卡洛抽样 Shapley 值——随机采样大量联盟，用边际贡献的平均值逼近精确 Shapley 值。但这个方案需要注意收敛性，至少采样几万次才能稳定，而且每次采样仍然要跑 MILP，实际工程中并不可观。我自己的处理方式是：正文示例用 3 主体精确计算，扩展实验单独写一个抽样函数做趋势分析，不追求全部联盟穷举。

4.3 三主体示例：一张表看穿整个分配结果

用前面三个主体跑一遍，假设独立运行成本分别是 CHP 12.3 万元、电锅炉 8.7 万元、储能 5.2 万元。各种联盟的合作调度成本算出后，总收益如下：

按 Shapley 权重公式计算，CHP 分得 1.83 万元，电锅炉分得 1.63 万元，储能分得 0.73 万元，加总约 4.19 万元（微差来自四舍五入）。数字背后是一条规则：CHP 在任意联盟中都能提供便宜的余热，所以边际贡献最高；储能必须和别的成员组队才有价值，单独存在时几乎无法创造收益，所以分得最少。

这张表非常值得对着源码去复现，因为绝大多数论文里的结果表就是这样生成的。只要你能独立复现出这张表，说明你从建模、调度到 Shapley 的整条链路已经跑通，后面改数据、改主体数量只是工作量问题。

5. 避坑指南：复现和改编时最容易翻车的五个点

5.1 现象：Shapley 算出的“收益”出现负值，甚至比独立运行的成本还高

这是我见过最多的复现失败案例。原因几乎只有一个：把联盟运行成本直接当成特征函数 v(S) 代入 Shapley 公式。因为成本是正的，边际贡献差值算出来当然是正负混乱。解决方法是回到第 4.1 节，用“独立成本之和 − 联盟运行成本”做一次收益换算，确保空的、单成员的联盟收益为 0，所有多成员联盟收益为正。换算之后再算一遍，符号就正常了。

5.2 现象：主体数量加到 6 个以后，程序跑几个小时都没结果

原因：Shapley 值的全联盟枚举是 2^n 级别，每个联盟又需要求解一次 MILP，复杂度乘起来直接爆表。解决：要么把示例规模控制在 5 个以内，要么改用蒙特卡洛采样 Shapley，把联盟枚举改成随机抽样。我做扩展实验时会先把求解器的 TimeLimit 降下来，让每个联盟最多花 20 秒求解，然后单独记录哪些联盟超时并调整参数，避免整个程序卡死。

5.3 现象：Gurobi 报错Model is infeasible或License expired

模型不可行通常不是求解器的问题，而是你给的负荷曲线和设备容量不匹配。例如热负荷峰值 100 MW，但联盟里只有一台 60 MW 的电锅炉，热平衡约束永远无法满足，模型自然跑不出可行解。解决：先跑一个“所有外部购能不设上限”的松弛版本，确认问题本身有解；再把外部购能上限定得足够高，让求解器能找到可行域。授权报错则是环境变量问题，检查GRB_LICENSE_FILE是否指向了正确的证书路径，学术版要确认 IP 是否在学校许可范围内。

5.4 现象：Shapley 分配结果满足公式，却不被某个主体接受（低于独立收益）

Shapley 值的特点是公平性（对称性、可加性、虚拟博弈者性质），但它不保证个体理性——也就是分配结果不一定落在博弈的核心（Core）里。如果某主体分到的钱比它单干还少，联盟在实际中就无法成立。解决：在分配层加一个“让利”机制，或者改用核仁（Nucleolus）计算分配结果。很多源码只算 Shapley 不给这个约束，复现时别直接照搬，要先检查分配是否满足个人理性。

5.5 现象：和博客里的表格数值对不上

原因通常是数据口径不一致：博客里的成本单位可能用的是“元”，你用的是“万元”；或者典型日负荷曲线取的时点密度不同（1 小时间隔和 15 分钟间隔算出的总成本差很多）。解决：先把所有单位统一到MW·h和单一货币单位，再核对负荷曲线和价格曲线的峰值位置。如果数值差距依然很大，建议逐条对比热平衡约束里的热效率折算系数——这是综合能源系统里最容易差异化处理的地方，多一点少一点都会完全改变联盟收益。

6. 验证和进阶：分完钱只是开始，走上稳定分配还得靠这三步

Shapley 值算出结果之后，我每回都会继续做三道验证，不通过就不敢拿去写结论。

第一道是个体与集体理性检验：每个主体的分配值必须大于等于独立运行收益，所有分配值之和必须等于大联盟总收益。这个检验用两行代码就能做，但能立刻暴露模型里的逻辑漏洞。第二道是用核心法核验：把分配结果代入所有联盟的约束不等式 v(S) ≤ Σ_{i∈S} φ_i 中，一旦有某个子联盟被违反，说明该联盟有动机退出大联盟独立运营。遇到这种情况，我一般会调整收益口径，比如把网损节约额也计入联盟收益，或者引入让利因子，让核心约束重新满足。第三道是灵敏度测试：把天然气价格上浮 10%、峰谷电价差拉大，重算整个流程，看各主体的 Shapley 分配比例是否会剧烈变化。如果某个主体的分配比例从 30% 跳变到 60%，说明模型对某个参数过于敏感，论文结论里应该主动讨论这个边界，而不是藏起来。

做完这三步，这套“合作博弈 + 优化调度”的方案才算真正立住了。实际写论文或做汇报时，我还会额外做一张“联盟收益来源分解图”——把总收益拆成余热利用收益、储能套利收益、购能折扣收益三类，这样审稿人或导师能一眼看出收益不是从天上掉下来的。如果你拿到标题里的源程序，我建议先别急着改主体数量，而是用三主体复现一张我前面写的收益表，确认没有符号问题后再加复杂度。这一步走稳了，后面任何扩展都是体力活，不会再有玄学问题。希望这些调试路径能帮到你。

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