合作博弈在园区综合能源系统优化调度与利益分配中的Matlab实现

去年我在复现一个综合能源系统的优化调度项目，正好卡在一个特别现实的问题上：园区里风电场、光伏电站、燃气轮机和储能电池分属不同的运营主体，各自单独算账时，燃气机组要为风光波动频繁调峰，成本高得吓人；光伏正午被限电，白花花的出力只能扔掉；储能更是委屈，充放循环多了电池寿命掉得飞快。可一旦把这些主体拉进同一个调度平台联合优化，总运行成本唰地一下降下来了。但问题也随之而来——省下来的钱到底怎么分？分得不公平，联盟第二天就散伙。这里用到的正是合作博弈理论，而把博弈分配和优化调度在Matlab里完整落地，才是整套方案的真正难点。这篇就把我的思路、模型、代码和踩过的坑一并聊清楚，给打算复现同类项目的朋友一个可以直接参考的样本。

1. 单打独斗的困局：为什么综合能源系统一定要“结盟”

1.1 独立优化时到底亏在哪里

先明确一下“主体”的概念。在这个项目里，我假设园区内有三类独立运营的能源实体：主体1是风电场，拥有60MW风机；主体2是光伏电站，装机40MW；主体3是“灵活调节方”，手里握着30MW燃气轮机和20MWh/10MW的储能电池，同时掌握着向大电网购售电的合同权限。三者各自独立运行时，谁的日子都不好过。

风电的问题是出力随机性大。夜间大风多发时段，本地负荷根本吃不掉，又没有储能配合，只能被迫弃风，或者低价卖给电网。光伏更惨，午间出力高峰和负荷低谷错位，限电是家常便饭。主体3虽然调节能力强，但风光强的时候它得压低燃气出力，风光弱的时候又要紧急拉满，机组长期在非设计工况下运行，燃料耗率和运维成本都在涨。

把三者独立运行的成本简单加起来，会得到一个很高的数字。但这不意味着它们分开跑是“各凭本事”，而是整个园区的资源配置根本没有打通。风电光伏弃掉的电量，本可以用来替代燃气发电；储能的闲置容量，本可以在风光大发时吸收盈余、在负荷尖峰时释放。这种“互补价值”不通过合作是拿不到的。

1.2 合作博弈建模的适用性判断

合作博弈研究的核心问题是：一群参与者通过形成联盟能获得比单独行动更多的收益，那么这个多出来的收益该如何分配，才能让联盟在自愿参与的前提下稳定维持下去。

在综合能源系统的语境下，参与者集合N就是上述三个主体。特征函数v(S)表示任意子联盟S能实现的“最优收益”。如果从成本视角建模，v(S)就是联盟S在最优联合调度下的最小总运行成本。这里我需要强调一个关键前提：只有满足“成本次可加性”（即v(S∪T) ≤ v(S) + v(T)）时，合作才有经济动力。

实际项目中这个性质并不总是成立。比如两个主体之间输电通道容量受限，硬要合作反而导致潮流拥塞，联合运行成本可能高于独立之和。遇到这种情况，要么调整联盟结构，要么把网络约束从模型里明确拿出来，先判断可行域再说。我在测试中也发现，风光资源天然互补的场景（夜间风大、白天光强）最容易满足合作条件，这也是为什么大多数论文都喜欢拿风光储组合做例子——它有现实的经济逻辑做支撑。

2. 利益分配才是真正的难点：Shapley值、核仁与纳什谈判的取舍

2.1 先理解合作博弈的分配公理

合作博弈的分配向量x=(x1, x2, ..., xn)必须满足两条基本约束。第一条是“个体理性”：每个主体分配后的成本不能高于它单独运行的成本，否则它没有理由加入联盟。第二条是“整体理性”（也叫帕累托最优）：所有主体的分配之和必须等于大联盟的总成本，即∑xi = v(N)，也就是说合作省下的钱必须全部分完，不能有剩余。

这两条只是底线。实际分配方法还应该满足更高的公平性要求，这就引出了Shapley值、核仁这些工具。很多人第一次接触时觉得Shapley值就是个数学公式，其实它背后的思想很朴素：一个主体应该分到多少收益，取决于它对联盟整体的“边际贡献”有多大。一个主体如果走到哪里都能带来显著的成本下降，那它理应拿到更大的收益份额。

2.2 Shapley值的计算逻辑与直观理解

Shapley值的统一公式是：

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

其中v(S∪{i}) - v(S)是主体i加入联盟S时带来的边际收益。所有不含i的子联盟S都要参与求和，权重|S|!(n-|S|-1)! / n!反映的是“随机排列中，S恰好是i之前那些参与者”这个事件出现的概率。

我在给同行讲解时常用这样一个类比：一场宴席有三道菜，分别由三个人带来，现在要算每个人对“整桌丰盛程度”的贡献。Shapley值就是让三个人以所有可能的顺序依次入场，每人都记录自己入场前后整桌菜的变化，然后取所有顺序下的平均值。对三个主体的情况，可以手工列一张表来算。比如主体1的特征值φ1考虑了四种情况：空集时入场（边际值就是v({1})）、在主体2之后入场（边际值是v({1,2}) - v({2})）、在主体3之后入场（边际值是v({1,3}) - v({3})）、在主体2和3之后入场（边际值是v({1,2,3}) - v({2,3})），权重分别是1/3、1/6、1/6、1/3。这张表在后面的实例分析里我会直接算给大家看。

2.3 核仁：让最不满的联盟也挑不出毛病

Shapley值的缺点是它把所有子联盟的边际贡献一视同仁。实际谈判中，某个具体的子联盟（比如风电和光伏联合体）可能觉得自己的贡献被低估了。核仁的出发点完全不同：它先计算每个子联盟S对分配结果的不满程度，定义为“超额”e(S, x) = v(S) - ∑_{i∈S} xi。如果某个子联盟分配到的成本之和大于它独立运行的成本，这个超额就是正的，表示该子联盟有脱离联盟单干的动机。

核仁就是寻找一个分配，使得所有子联盟中最大超额尽可能小，并且按“从大到小逐级最小化”的顺序来确定唯一解。计算核仁通常要解一系列线性规划，第一轮是：

min ε s.t. ∑_{i∈N} xi = v(N) xi ≥ v({i})，∀i v(S) - ∑_{i∈S} xi ≤ ε，∀S ⊂ N

其中变量是x和ε。解出最优的ε后，把所有使超额恰好等于ε的约束固定为等式，再继续最小化次大的超额。对3~4个主体，手算和代码都能搞定；主体一多，约束规模指数增长，就得用约束生成法按需添加最紧的约束。

2.4 三种分配方法怎么选

我整理了在实践中对这几个方法的对比感受：

Shapley值的最大问题是它未必在核内，也就是说按Shapley值分配后，某个子联盟拿到的总收益可能反而不如自己单干，这会让联盟在实际落地时面临拆伙风险。核仁虽然计算繁琐，但天然把稳定性放第一位，在需要长期合作的综合能源项目里，工程价值更高。我最终的代码里两个都实现了，方便对照。

3. 优化调度模型：用数学语言描述“联盟总成本最小化”

3.1 目标函数的组成与物理意义

特征函数v(S)并不是凭空算出来的，它来自一个线性规划——联盟S在满足全部设备约束、负荷平衡约束和电网交互约束的前提下，最小化总运行成本。对每个联盟都要单独解一次这个优化，这正是整个代码框架的核心循环。

目标函数按规划时段T（通常取24小时）求和，包含五个部分。第一部分是燃气轮机的燃料成本，用出力除以发电效率再乘天然气热值与气价得到。第二部分是向大电网购电的成本，按分时电价计；如果联盟向电网反送电，则按较低的上网电价计入负成本。第三部分是各设备的运行维护成本，风电、光伏、燃气轮机、储能按各自的单位运维单价乘出力或充放电功率。第四部分是弃风弃光惩罚项，这个很容易被初学者忽略，但必须加：如果不给弃能惩罚，优化模型出于省钱考虑会“故意”弃掉新能源，得出一个不真实的低成本；加了惩罚后，模型才会在设备约束允许的范围内尽量消纳风光。第五部分是储能循环退化成本，用充放电功率乘一个等价老化损失系数，防止模型为了套利频繁强充强放。

目标函数的完整形式可以写成：

min Σ_{t=1..T} [ c_gas·p_gt(t)/η_gt + c_buy(t)·p_buy(t) - c_sell(t)·p_sell(t) + c_om_gt·p_gt(t) + c_om_wt·p_wt(t) + c_om_pv·p_pv(t) + c_om_ess·(p_ch(t)+p_dis(t)) + c_pen_wt·(p_wt_pred(t)-p_wt(t)) + c_pen_pv·(p_pv_pred(t)-p_pv(t)) ]

3.2 约束条件：一个都不能少

约束条件分五组，每组都有实际的工程含义。

第一组是功率平衡约束，联盟内所有电源出力、储能放电和购电量之和，必须等于负荷、储能充电和售电量之和。这条等式反映的是“同一时刻电能必须实时平衡”的物理规律，在Matlab里用等式约束直接写入优化模型。

第二组是设备出力上下限约束。风电、光伏的实际出力必须介于0和预测出力之间，燃气轮机不超过装机容量，储能充放电功率不超最大功率，从电网购售电也不超过联络线容量。

第三组是燃气轮机的爬坡约束。机组在两个相邻时段的出力变化不能超过爬坡速率上限。这个约束特别重要，因为联合调度之所以能省钱，很大程度体现在燃气机组可以更平稳地运行，而非跟着风光剧烈波动。

第四组是储能运行约束。储能荷电状态（SOC）的动态递推关系、上下限约束，以及调度周期初末SOC相等的循环约束。这条约束在模型中让储能“只搬运电量、不凭空生电”，否则优化器会白捡能量。

第五组是风光出力与预测上限的关系，同时把弃风弃光量定义为预测值与实际出力的差。这里有一点要特别注意：不同联盟里可用的设备集合不同，约束条件也随之裁剪。比如联盟只有风电时，模型里根本没有储能和燃机的变量；联盟加入储能后，功率平衡方式立刻丰富起来。在代码实现时，我用一个members列表动态生成变量和约束，而不是把所有设备的约束一次性写死。

3.3 求解框架：YALMIP建模加线性规划求解器

优化模型本质上是线性规划（目标函数和约束都是线性的），所以求解非常成熟。建模阶段我推荐用YALMIP工具箱，它最大的好处是变量声明和约束集合的写法贴近数学公式，代码可读性好，适合快速迭代。底层求解器可以用Gurobi、Cplex，也可以用Matlab自带的linprog。对于24×几小时的连续线性规划，求解时间通常不到一秒，整个合作博弈流程的瓶颈根本不在求解器，而在联盟数量。

在Matlab中，核心求解代码用YALMIP可以这样组织：

function [cost, result] = coalition_optimization(members, data) % members: 联盟成员索引，如 [1] 或 [1 3] 或 [1 2 3] % data: 负荷、预测出力、设备参数、电价等结构体 T = data.T; % 声明决策变量 p_gt = sdpvar(1, T); p_wt = sdpvar(1, T); p_pv = sdpvar(1, T); p_ch = sdpvar(1, T); p_dis = sdpvar(1, T); soc = sdpvar(1, T); p_buy = sdpvar(1, T); p_sell = sdpvar(1, T); variables_gt = 0; % 标识是否有燃机 % ...（根据members动态生成变量集合） Constraints = []; if ismember(3, members) Constraints = [Constraints, 0 <= p_gt <= data.P_gt_max]; Constraints = [Constraints, -data.R_gt_down <= diff(p_gt) <= data.R_gt_up]; variables_gt = 1; end % ...（其余设备同理） % 功率平衡 Constraints = [Constraints, p_wt + p_pv + variables_gt*p_gt + p_dis + p_buy == ... data.L_load + p_ch + p_sell]; % 储能SOC soc(1) = data.soc_init + data.eta_ch * p_ch(1) / data.E_cap - p_dis(1) / (data.eta_dis * data.E_cap); % ...（递推约束循环写入） % 目标函数 Objective = sum(data.c_gas * p_gt / data.eta_gt + data.c_buy .* p_buy ... - data.c_sell .* p_sell + ...各运维项... + ... 数据.惩罚项); ops = sdpsettings('solver', data.solver, 'verbose', 0); optimize(Constraints, Objective, ops); cost = value(Objective); result = value([p_gt; p_wt; p_pv; p_ch; p_dis; soc; p_buy; p_sell]); end

这段代码里我用ismember判断联盟成员，从而决定要不要加入燃机、储能、风电或光伏的约束。注意，当某个联盟没有燃机时，p_gt变量虽然仍被声明，但对应的约束没有加入，目标函数里它的系数为0，优化器自然会把它置为0，不会造成逻辑混乱。这个小技巧让代码简洁了不少，但前提是目标函数和约束中所有涉及该变量的项都必须同步“屏蔽”，否则变量会被白白优化出非零值。

4. Matlab代码实现：从联盟特征值计算到Shapley值的一站式流程

4.1 代码总架构与关键数据结构

整个项目的代码我拆成了五个模块：数据准备模块负责读入负荷预测、风光预测、设备参数、分时电价，封装成data结构体；联盟枚举模块负责生成所有非空子联盟并建立索引表；联盟优化模块就是上一节的coalition_optimization函数，对任意联盟求解最小运行成本；分配计算模块实现Shapley值和核仁；结果分析模块负责验证理性约束、输出对比图表。

联盟索引是代码的地基。n个主体共有2^n-1个非空子集，我用一个整数mask来表示联盟：mask的第i位为1表示主体i在联盟中。比如n=3时，mask=3的二进制是011，代表联盟{1,2}。特征函数v可以存成一个长度为2^n的向量，v(0)=0对应空集，v(mask)就是联盟mask的最小总成本。这个映射关系在后面的Shapley值计算里会频繁用到，务必保持清晰。

数据准备阶段的负荷曲线、风光预测曲线我用历史数据加上随机扰动生成，时间分辨率为1小时，周期24小时。设备参数和电价分别放在data结构体里。为了后续分析方便，我把独立运行的成本（各单元素联盟）和大联盟成本都单独存下来，作为个体理性和整体理性的校验基准。

4.2 特征函数的并行计算与缓存

计算特征函数是整个流程里最耗时的一环。对n个主体，需要求解2^n-1个优化问题。n=3时只要解7个，毫秒级完成；n=7时要解127个，仍然可以接受；n=10时数量突破1023个，求解时间开始变得明显；再往上就必须考虑并行和近似方法了。

我代码里直接用parfor并行求解，每个联盟调用一次coalition_optimization。要注意parfor里不能访问需要动态更新的工作区变量，所以我把所有参数都封装到data中，确保每个worker独立运算。另外，我加了一层缓存机制：如果某个联盟的子集结构和大联盟高度相似（比如只差一个风机），可以用大联盟的优化结果做热启动，把初始解传给求解器，能减少相当多迭代时间。不过Matlab的linprog和YALMIP对初始解的接口处理不太一样，热启动的价值取决于求解器，测试下来Gurobi的热启动收益最明显。

4.3 Shapley值的代码实现

Shapley值的代码实现最核心的环节是把联盟mask映射到特征值索引。我直接用mask本身作为v数组的下标，省去了字典映射的开销。完整的核心代码如下：

function phi = shapley_value(v, N) % v: 长度为2^N的向量，v(0)=0, v(mask)=联盟mask的最小运行成本 phi = zeros(1, N); for i = 1:N bit_i = 2^(i-1); phi_i = 0; for mask = 0:(2^N - 1) if bitget(mask, i) == 0 S_card = nnz(bitget(mask, 1:N)); weight = factorial(S_card) * factorial(N - S_card - 1) / factorial(N); mask_new = bitor(mask, bit_i); v_old = 0; if mask > 0 v_old = v(mask); end phi_i = phi_i + weight * (v(mask_new) - v_old); end end phi(i) = phi_i; end end

这里最容易被忽略的是mask=0的情况。空集的成本是0，v(0)不是从特征函数数组里取出来的，而是在代码里显式赋0。如果你直接在数组下标上用mask，Matlab不允许下标为0，所以要么把v数组整体偏移一位（v(1)对应空集，v(2)对应联盟{1}），要么像上面这样在访问前判空。我选择了后者，代码直观且不容易错。

4.4 核仁求解的线性规划实现

核仁的第一层线性规划用YALMIP写成这样：

function x = nucleolus(v, N, data) % v: 特征函数向量，v(0)=0 x_var = sdpvar(N, 1); eps_var = sdpvar(1); Objective = eps_var; Constraints = [sum(x_var) == v(2^N - 1)]; % 个体理性：每个主体分配后成本不高于独立运行成本 for i = 1:N Constraints = [Constraints, x_var(i) >= v(2^(i-1))]; end % 所有非空真子联盟：超额不大于eps for mask = 1:(2^N - 2) S = find(bitget(mask, 1:N)); Constraints = [Constraints, v(mask) - sum(x_var(S)) <= eps_var]; end ops = sdpsettings('solver', data.solver, 'verbose', 0); optimize(Constraints, Objective, ops); x = value(x_var); end

这一层求出的是“最大超额最小化”的分配，但还不是严格意义的核仁——严格的核仁在最大超额被压缩到最低后，还要继续最小化第二大的超额，直到唯一确定。对于大多数工程场景，第一层解已经足够稳定。我在项目里也加了后续迭代：先记录使超额等于当前最优ε的联盟，把它们转为等式约束，再求解次小超额的最小化，循环直到分配唯一。代码篇幅会明显增加，但对n=3~5的情况运行依然很快。

5. 算例验证：三主体系统的完整结果与理性校验

5.1 场景参数设置

为了把结果讲明白，我给出一个具体的参数表。负荷曲线设为典型工业园区日负荷，峰值为100MW；风电预测曲线在凌晨达到高点，光伏预测曲线在午间达到高点，两者错峰分布；分时电价按峰谷平三段设置。设备参数如下：

基于这组参数，我执行了全部7个非空联盟的优化调度，得到特征函数表：

5.2 合作红利与Shapley值的直接计算

从表中可以验证成本次可加性：三个主体独立运行总成本是58.2 + 41.6 + 77.3 = 177.1万元，大联盟成本146.9万元，合作节约30.2万元，相当于整体成本下降17.1%。任意两主体联盟也满足次可加性，比如风电和光伏单独运行是99.8万，联合运行是92.5万，节约7.3万。这说明在这个场景里合作带来的是普遍增益，不是个别组合的偶然结果。

接下来算Shapley值。以主体1（风电）为例，我按四个不包含它的联盟逐项计算：

• S=∅，权重1/3，边际贡献v({1}) - 0 = 58.2，贡献19.4
• S={2}，权重1/6，边际贡献v({1,2}) - v({2}) = 92.5 - 41.6 = 50.9，贡献8.483
• S={3}，权重1/6，边际贡献v({1,3}) - v({3}) = 118.7 - 77.3 = 41.4，贡献6.9
• S={2,3}，权重1/3，边际贡献v({1,2,3}) - v({2,3}) = 146.9 - 103.1 = 43.8，贡献14.6

四者相加得到φ1 = 49.383万元。同理算出φ2 = 33.283万元，φ3 = 64.233万元，三者之和恰好等于146.9万元，满足整体理性。

5.3 分配结果的理性校验与直觉印证

校验个体理性：风电分配后成本49.383 < 58.2，光伏33.283 < 41.6，灵活资源64.233 < 77.3，三个主体都从合作中受益，谁也不吃亏。从节约额看，风电节约8.817万，光伏节约8.317万，灵活资源节约13.067万。灵活资源分得的红利最大，完全符合直觉——燃气轮机和储能承担了平衡风光波动的核心角色，边际贡献最高，自然应获得更大的合作收益。风电和光伏由于出力曲线错峰，互为补充，也各自拿到了可观的节约。

核仁计算得到的分配值是49.52、33.43、63.95万元，和Shapley值相差不大。这说明在这个数据条件下，分配问题本身不敏感，两种方法给出的结果都在核内，联盟稳定性都能得到保证。如果你改用一组负荷与风光曲线特别极端的数据，两者的差异会更明显，那时核仁的相对优势就会体现出来。

6. 复现这套代码时容易踩的坑与我的改进思路

6.1 特征函数不满足次可加性：先查数据再查求解器

我最开始跑完整流程时遇到过一个问题：联盟{1,2}的成本竟然高于v({1}) + v({2})，合作反而亏本。当时第一反应是模型约束写错了，排查半天发现是求解器收敛容差设置太宽松，导致多个联盟的优化结果精度不一致，成本误差被放大。把Gurobi的MIPGap和最优性容差收紧后，次可加性就正常了。

另一个常见原因纯粹是数据层面的：如果两个主体的出力曲线峰值完全重叠，又缺乏储能配合，它们联合消纳的空间确实有限，合作的收益可能趋近于零甚至为负。这不是代码bug，而是物理上就不适合结盟。遇到这种情况，应该从联盟结构入手，把第三个调节资源拉进来，而不是强行让这两个主体单独合作。

6.2 Shapley值计算量爆炸的降维思路

主体数量从3涨到8时，2^n-1就从7涨到255；涨到12时，4095个联盟每个都要解一个大型优化问题，已经不太现实。我的处理有三个方向。

第一，并行化。parfor能线性加速，8~12个联盟规模时效果非常明显。第二，近似抽样。随机抽取一部分有代表性的联盟计算特征值，用蒙特卡洛方法估计Shapley值，工程上误差控制在2%以内。第三，结构简化。把功能相近的主体合并，比如多台同型号燃气机组合成一个聚合体，减少参与者数量，从根源上降低联盟规模。

6.3 面对不确定性：从确定性调度走向场景鲁棒

这个项目做的是日前确定性调度，但实际风光预测误差无法避免。合作博弈的分配结果对预测误差很敏感——如果某天风光实际出力远低于预测，原计划的联合调度方案就会失真，分配的基础也随之瓦解。

我尝试过的改进是引入场景法：生成多个风光出力场景，每个场景下都执行一次合作博弈优化，最终的分配结果是各场景分配值的概率加权平均。这样算出来的Shapley值更稳健，但计算量也成倍增长。更前沿的思路是做分布鲁棒合作博弈，把极端场景纳入约束，但实现复杂度较高，适合作为后续研究的方向，不适合第一次复现时就纳入。

6.4 落地时别忘了可解释性

最后分享一点非技术的心得。整套模型在学术上很漂亮，但真正要拿去和园区运营方谈落地时，对方最关心的不是Shapley值和核仁的数学定义，而是“我到底能省多少钱、凭什么我拿这么多”。我后来在项目报告里做了一张饼图，把30.2万元合作红利的构成拆解开——多少来自减少弃风弃光、多少来自提高燃气机组运行效率、多少来自储能峰谷套利——再和每个主体的分配额对应起来。这种“红利溯源”的做法极大增强了各方对分配结果的接受度，比任何数学证明都管用。合作博弈的最终目的是让联盟稳定、让系统更高效地运行，如果分配结果不能让每个成员信服，再精确的算法也只是纸上谈兵。