企业官网建设流程全解析

做综合能源系统调度这两年，最让我头疼的不是设备模型多复杂，而是系统里电、气、热设备根本不属于同一个运营方。光伏可能在用户侧屋顶，燃气机组归能源公司，蓄热罐又是热力站的资产。如果硬按单一运营商做全局优化，物理上能算出最低成本，但现实中根本执行不下去——因为没人愿意无条件牺牲自己的利益去成全系统。后来我把思路切到合作博弈框架，用 Matlab 写了一套“联盟优化 + 利益分配 + 分配结果再反馈到调度约束”的完整实现，才把账算明白。这篇就聊聊这套实现里最容易被忽略、也最能决定论文和工程方案能不能站住脚的细节，主要面向正在用 Matlab 做综合能源系统优化调度、想把合作博弈落下去的同学们。

1. 为什么综合能源系统的调度要引入合作博弈

1.1 同一张网，多个账本

传统调度模型通常假设一个“上帝视角”的调度中心，所有设备、负荷、储能在同一面旗下运行，目标函数是总成本最小或总收益最大。但在实际项目里，多能互补系统经常是多个主体各自投资、各自运营的：电网侧关心峰谷套利，热力公司关心供暖保障，燃气公司关心气价波动，用户侧聚合商关心可调负荷的补偿。这种场景下你按总成本最小去调度，就可能出现某个主体的设备被压到极低出力、收益下降的情况，协调自然失败。

合作博弈解决的正是“多个独立决策主体如何行动能整体最好，同时每个人又愿意参与”的问题。它把各主体看作博弈玩家，先把大家绑定成一个联盟去优化联盟总收益，再把联盟产生的合作收益按某种兼顾公平和稳定的方式分回去。引入合作博弈后，调度问题变成一个“先算锅里的肉，再分锅里的肉”的过程。肉不够大，方案没有吸引力；肉分不匀，方案没法落地。

1.2 合作节省的钱从哪来

用一句话概括：合作收益来自各主体单独运行时无法利用的互补性。举个我在某区域项目里简化的例子。光伏富余时段，用户侧用不完的电如果只馈到电网，收益很低；但隔壁有台电锅炉可以吃电产热给热力管网，这就是互补。燃气机组余热可以进余热锅炉供给热负荷，减少电锅炉耗电，也是互补。储能、蓄热罐还能把时段间的价差套出来。

这些互补性意味着，如果每个主体各自优化后再简单并网交易，会漏掉很大一部分协调收益。联盟优化则把所有资源共享到一个统一优化问题里，让燃气轮机出力、电锅炉功率、储能充放策略、热罐蓄放热策略同时决策，从而拿到更低的总运行成本。系统单主体优化得出的最低成本，和合作联盟优化得出的最低成本，差额就是合作博弈要分的合作节省额。

1.3 调度与分配的耦合才是难点

有一点必须提醒：很多文章把利益分配当成“事后算账”，跑完联盟优化用 Shapley 值分一下账就结束。但在工程上，分配结果会影响下一步的调度决策——如果某个主体分到的收益低于他单独经营的收益，他会退出联盟，然后联盟结构变了，调度方案也得跟着变。所以严格的模型应该是“调度—分配—校验参与约束—必要时重新调度”的迭代闭环。

我在模拟项目 X 里就是用双层结构实现的：内层做各主体可形成的联盟运行优化，外层做收益分配与稳定性校验，如果校验不通过就把参与约束以惩罚项的形式加回优化目标再求解，循环若干轮直到所有主体都愿意留在联盟内。

2. 从物理模型到特征函数：IES合作博弈建模的关键细节

2.1 能源集线器：把多能耦合写进约束

综合能源系统的物理模型我习惯用能源集线器（Energy Hub）形式描述。它把输入端（电网购电、天然气站购气、光伏、风电）和输出端（电负荷、热负荷、冷负荷）通过设备转换矩阵连接起来。建模时我通常把每个参与主体视为拥有若干设备的集合，设备类型主要有这几类：

• 燃气轮机：天然气转电力，附带热电比约束
• 余热锅炉：吸收燃气轮机余热生产热力
• 电锅炉：电力转热，效率一般取 0.9 到 0.98
• 热泵：用电从低温热源提热，COP 通常 2.5 到 4
• 蓄电/蓄热：储能类，需要 SOC 递推和容量约束
• 光伏/风电：出力曲线给定，可弃风弃光

电网购电、天然气购气和设备出力构成决策变量，约束至少包括电功率平衡、热功率平衡、设备出力上下限、爬坡约束、储能 SOC 连续性。时间粒度我常用 1 小时，一天 24 个时段总量不算大；如果你做高比例可再生场景，可以细化到 15 分钟，但求解时间会明显上涨。

2.2 特征函数怎么算：每个联盟都要重新优化

合作博弈里最关键的概念是特征函数 v(S)，表示联盟 S 合作能带来的额外收益。在我这套模型里，联盟收益统一定义为“独立运行总成本 - 联盟最优运行成本”，也就是成本节约额。因此必须先对每个非空联盟 S 都做一次优化调度，得到对应最小总成本 C(S)，再算 v(S)=ΣC({i}) - C(S)。

这里有个特别容易踩的坑：计算 C({i}) 的口径必须和 C(S) 完全一致。单独运行时主体 i 的设备只有它自己的，购电、购气按它自己的负荷；联盟运行时，原属不同主体的设备可以协同，但固定成本、惩罚项、违约惩罚这些非运行成本项目不要左边算进去右边没算。口径不一致会导致分配结果出现负数或超额，后面校验根本过不了。

所有非空联盟的数量是 2^N-1。N=3 时只有 7 个子问题，Matlab 跑起来很轻松；N=4 时 15 个；N=5 时 31 个，如果每个子问题都是 MILP，单核可能要跑几十分钟，这时候就该考虑并行或近似方法，具体在第 5 节展开。

2.3 三主体算例：手工验证特征函数与 Shapley 值

为了确认程序没有算错，我在搭代码之前先手推了一个 3 主体简化案例。假设三个主体独立运行成本分别是 c1=12、c2=10、c3=8，联盟后运行成本及成本节约额如下：

这个表翻译成程序逻辑就是：对每个联盟独立调一次度，取成本，求差。程序跑出来的 v 如果跟手算不一致，一定先查是不是联盟约束写漏了。

3. Shapley值、核仁法与最小核心：分配方法怎么选，怎么用

3.1 Shapley值：公平性最强的默认选项

Shapley 值的基本思想是：一个主体得到的分配等于它所有可能加入顺序中对联盟边际贡献的平均值。拿上面的例子手算：主体 1 的 Shapley 值为 1.5，主体 2 为 1.5，主体 3 为 2，合计 5，正好等于大联盟的总节约额。

计算上，每个主体 i 的 Shapley 值公式是：

φ_i = Σ_{S⊆N{i}} |S|!(n-|S|-1)! / n! × [v(S∪{i}) - v(S)]

代码实现时直接枚举不含 i 的联盟子集，乘上权重即可。当 N 不超过 5 时，穷举完全没问题。Shapley 值有很好的性质：总和等于大联盟总收益、对称、哑元为 0、可加性。工程上大家默认它最公平，不容易引起争议。

但 Shapley 值没有保证分配结果一定落在核心内，也就是说可能出现某个子联盟分到的总收益少于它独立合作能拿到的收益，这时子联盟成员有动机退出。我在项目里遇到过：某个含储热主体的小联盟分到的收益偏低，而它单独跟热力公司抱团能拿更多，于是靠 Shapley 分完的联盟直接不稳定。

3.2 核仁与最小核心：稳定性优先的替代方案

如果稳定性是第一诉求，可以用核仁法或最小核心法。核仁的思路是寻找一个分配，使“最不满意”的联盟的不满程度尽量小。不满程度通常用超额（excess）衡量：e(S)=v(S)-x(S)，其中 x(S) 是联盟 S 所有成员分到的收益之和。核仁算法要逐步最小化排序后的超额向量，本质是一个系列线性规划。

最小核心法是更简单的松弛版本：引入一个变量 ε，约束所有联盟的 x(S)+ε≥v(S)，然后最小化 ε。它的含义是允许所有联盟的整体分配都往公平性上稍微偏离一点，但保证没有一个联盟会退出太多。好处是线性规划即可求解，坏处是解不唯一，两个不同的解都是最小核心点，所以有时还需要再叠一个利润均等的目标来选点。

我之前做过三种方法的数值对比，见下表：

3.3 稳定性校验：分配结果能不能让所有人留下

不管用哪种方法，我最终都会做一项校验：对每个非空联盟 S，检查 Σx_i ≥ v(S)。如果某个联盟不满足，就意味着那批成员合起来没拿到他们独自合作能拿到的收益，联盟有解体风险。我的处理方法是先尝试调整分配规则，比如改用核仁；如果还是不行，就说明特征函数本身有问题，需要回查 v(S) 的计算口径。

校验通过之后，再把分配结果回填到调度模型里。具体做法是：把各主体的参与约束（自己的分配收益不低于单独经营收益）写成软约束，加进联盟优化目标函数里，稳定解作为下一步调度的边界条件。

4. Matlab实现框架与核心代码拆解

4.1 总框架：枚举联盟、求解特征值、计算分配

我用 Matlab 实现时把整个流程拆成五步，每步一个函数，方便单独调试：

1. 数据准备：读入负荷曲线、设备参数、分时电价、气价
2. 联盟枚举：生成 2^N-1 个联盟的成员组合
3. 特征值求解：对每个联盟调用统一优化函数，得到 C(S)，保存到向量
4. 分配计算：根据特征值向量分别算 Shapley 值、核仁或最小核心
5. 结果后处理：输出各主体分配收益、运行成本、稳定性检查表

下面是联盟枚举与特征值计算的骨架：

N = 3; % 参与主体数量 numCoal = 2^N - 1; c = zeros(numCoal, 1); v = zeros(numCoal, 1); for idx = 1:numCoal members = find(bitget(idx, 1:N)); c(idx) = solve_union_cost(members, data, params); end cSingle = zeros(N, 1); for i = 1:N cSingle(i) = c(find(bitget(1, i) & (1:numCoal))); % 简化示意 end totalStandalone = sum(cSingle); for idx = 1:numCoal members = find(bitget(idx, 1:N)); singleCost = sum(cSingle(members)); v(idx) = singleCost - c(idx); end

实际工程里我不会用find反复匹配单主体索引，而是直接把单主体成本单独存一份，避免逻辑混乱。核心是solve_union_cost这个函数，它接受成员编号数组，内部构建对应设备的优化模型并调用求解器。

4.2 统一优化函数：YALMIP 建模与求解器配置

我习惯用 YALMIP 做建模界面，底层求解器用 Gurobi。如果你没有 Gurobi 授权，也可以用 Matlab 自带的 intlinprog 处理 MILP，线性规化问题用 linprog。快速安装配置参考如下：

% 求解器路径与测试（以Gurobi为例） addpath('D:/gurobi/win64/matlab'); gurobi_setup(); currentSolver = yalmp;

模拟项目 X 中，单联盟 24 时段 MILP 模型大概有几百个变量、上千个约束，Gurobi 默认参数就能一两分钟内收敛。YALMIP 里核心建模片段大致是这样：

x = sdpvar(repmat(nT, 1, nGen), repmat(1, 1, nGen)); % 简写示意 F = [x(1) >= lower(1), x(1) <= upper(1)]; % 电平衡约束 F = [F, sum(pGen) + buyPower + pv - loadElec == chargeP - dischargeP]; % 热平衡约束 F = [F, sum(hGen) + boilerHeat + hpHeat - loadHeat == heatIn - heatOut]; % SOC递推 for t = 2:nT F = [F, soc(t) == soc(t-1) + etaC * chargeP(t) - dischargeP(t) / etaD]; end optimize(F, sum(operationCost), options);

写约束时一定注意电、热、气单位的统一。热负荷和电负荷经常一个是 GJ，一个是 MWh；燃气轮机热电比决定余热量，如果单位换错，特征值会差出好几倍。我吃过这个亏，后来所有热功率统一折算成 kW，燃气的单位热值也换算成以电等价的 kW 口径，问题才消失。

4.3 Shapley 值和核仁法的Matlab实现

Shapley 值实现按公式穷举即可：

phi = zeros(1, N); for i = 1:N for S = 0:(2^N - 1) if bitget(S, i) == 0 % S不包含i s = nnz(bitget(S, 1:N)); idxS = S + 1; idxSUi = bitor(S, bitshift(1, i-1)) + 1; if idxS < 1 || idxS > numCoal || idxSUi < 1 || idxSUi > numCoal continue; end phi(i) = phi(i) + factorial(s) * factorial(N - s - 1) / factorial(N) ... * (v(idxSUi) - v(idxS)); end end end

核仁的 LP 可以用 YALMIP 快速搭出来，核心是epsilon和分配向量x的优化：

x = sdpvar(1, N); epsilon = sdpvar(1, 1); F = [sum(x) == v(end), x >= 0]; for idx = 1:numCoal members = find(bitget(idx, 1:N)); F = [F, sum(x(members)) + epsilon >= v(idx)]; end optimize(F, epsilon, options);

这个模型求解非常快，几行代码就能替代多轮 LP 的繁琐实现。如果你要标准核仁，还需要对超额向量做词典序最小化，我一般用迭代加约束的方式实现，代码会稍复杂一点，但论文里通常只需给出核心法结果即可。

4.4 并行计算和结果缓存

特征值求解是整个流程最耗时的部分。N 达到 4 或 5 时，我会用 Parallel Computing Toolbox 的 parfor 替代 for：

parfor idx = 1:numCoal c(idx) = solve_union_cost(find(bitget(idx, 1:N)), data, params); end

还要把算过的特征值存成 mat 文件。因为做参数敏感性分析时，电价或气价一变，某些联盟的最优成本变化很小，缓存后能省掉大量重复求解。这是我在处理某高比例可再生场景时总结出来的高效经验，看似不起眼，实际能节省一个数量级的时间。

5. 从集中式到分布式：大规模场景的可扩展路径

5.1 为什么必须考虑分布式化

当 N 超过 6，穷举 2^N-1 个联盟就变得不现实；更麻烦的是，集成系统内部各主体未必愿意把自己的设备参数和负荷数据全部交给一个中心求解器。这种情况下需要引入分布式优化框架。我在扩展模拟项目 X 时采用了 ADMM 思路：每个主体内部保留自己的设备模型和隐私数据，只通过公共耦合变量和价格信号跟其他主体交换信息。

ADMM 的耦合变量通常是联络线功率、热网交换热量、储能共享电量。每个主体求解自己的子问题，然后更新拉格朗日惩罚项，循环迭代。这个框架的优点是能保留 Shapley 值计算所需的各联盟成本估计，缺点是收敛速度受惩罚参数影响很大，需要手动调参。我的建议是先用集中式结果作为基线，验证分布式优化得到的系统总成本和集中式相差不大，再继续后面的分配环节。

5.2 Shapley 值的近似计算

大规模下精确 Shapley 值必须近似。我试过两种比较有效：

一是基于采样的边际贡献估计。随机抽取大量排列，统计主体 i 在每个排列中的边际贡献，取平均值。排列数量从几千到几万，看你对误差的容忍度。

二是基于多线性扩展的近似公式。Shapley 值本质上是值函数在某个点上的离散梯度平均，当 N 很大但设备规模相似时，可以用连续化后的解析近似。

更简单但不那么严谨的做法是直接用比例分摊或最小核心做粗分配，然后检验稳定性。我在某些快速方案里会用最小核心替代 Shapley，作为初步筛选；如果最小的几个联盟都没有退出动机，就不再跑采样近似了。

5.3 鲁棒和不确定性场景的扩展方向

合作收益很大一部分来自可调节设备对价格波动的响应。如果光伏、风电出力不确定，后面还可以把模型扩展成两阶段鲁棒或随机规划，让特征函数变成期望值或最坏情况下的收益。比如天然气价格上浮 10% 时，某个联盟的收益可能由正转负，这种情况在分配时必须考虑进去。我在项目中把分时电价和负荷预测的不确定性用若干典型场景替代，发现各联盟收益排序会变化，Shapley 值也跟着变，所以论文里做参数敏感性时一定要写明“分配结果基于预测场景集”。

6. 调试与验证：几个最常见的翻车现场和排查方法

6.1 手工算例先行

我在 2.3 节用了三主体算例。这个算例小到什么程度？小到可以在纸上验算。程序跑完后，先检查 v 是否满足超可加性：

v({1,2}) = 2，v({3}) = 0，v({1,2,3}) = 5，满足；再看 v({1,3}) + v({2}) = 3，也满足。

如果发现超可加性不满足，比如 v({1,2})+v({3}) > v({1,2,3})，那说明联盟模型里存在“负耦合”，比如某个设备重复投资、容量互相挤占。这通常不是算法问题，而是物理模型建模问题。

6.2 分配结果出现负数或超额不合理

我遇到过用 Shapley 值分完，某个主体分到的收益为负。查到最后原因很简单：算单主体成本时把固定电价的容量费也算了，算联盟成本时却忘了加总固定费用。所以我的建议是跑任何验证前先检查 C({i}) 和 C(S) 是否在同一口径下。把所有非运行成本统一抽出来，单独做成常数项，而不是揉进每个子问题的目标函数，否则很容易重复计算。

6.3 MILP 求解时间爆炸

如果联盟优化模型里启停变量太多，24 时段也会卡住。第一次尝试可以把启停约束改成连续变量，用线性松弛判断成本差异，这样能快速定位是不是整数变量导致的复杂度。另一种办法是给 Gurobi 提供启发式初始解：先按不含启停的 LP 求解，然后用结果去设定启停变量的初始值，能大幅缩短分支定界时间。

6.4 YALMIP 与求解器接口报错

这类问题大多出在求解器路径没有配置好，或者版本不匹配。报错信息只要出现“Solver not found”或“No appropriate solver”，先检查 addpath 是否真的生效，再用yalmiptest跑一下诊断。我在某次升级后 Gurobi 版本变了，YALMIP 提示找不到求解器，重装接口之后恢复。这里不建议在模型里混用多个求解器，因为不同求解器的数值容差不一样，会对特征值微小的差值产生影响。

6.5 单位不统一导致特征值巨变

热负荷和电负荷单位不同、燃气的热值与电力换算系数不一致，是隐藏最深的错误。我在某次计算中把热负荷误用成 MWh，实际题目给的是 GJ，导致单主体成本和联盟成本同时放大，虽然 v 看起来还算正常，但分配比例失真。遇到数值异常，先检查单位矩阵，再检查效率乘以输入功率的方向是否一致。

最后分享几个我自己的实操习惯

纸质手算三主体算例永远比先跑 N=5 的大算例高效得多。每次改完模型，我都先把 N 临时设成 3，跑一组手工可验证的数据，通过后再放大规模。分配方法我一般同时算 Shapley 值和最小核心，论文里两张表放一起对比，工程汇报里只给 Shapley 值加稳定性校验结果。

如果你正要把这套 Matlab 代码写进研究或项目报告，我的建议是：不要把分配结果当成最终输出，一定要把它作为调度方案的一个校验环节，另外把特征值求解结果做成 Excel 或 CSV 表导出，方便回溯哪一次求解用的哪一套参数。这套流程我用下来最大的体会是，它最磨人的不是优化算法本身，而是把合作博弈的数学定义翻译成准确、口径一致的运行优化模型。一旦这一步做扎实，后面的 Shapley、核仁、分布式扩展就都有了稳定基础。