企业官网建设流程全解析

从接手第一个综合能源系统协调调度项目开始，我就被“利益分配”这四个字卡了很长时间。一开始想得简单：多能互补、联合优化，把总成本压下来，大家都省钱不就完事了？但真正做下来才发现，最大的难点根本不在物理层面的潮流平衡或者储能充放电策略，而在于一个特别现实的问题——“省下来的钱该怎么分”。分得合理，各方才有动力配合调度；分得不合理，哪怕总收益再高，第二天联盟就散了。这也是为什么合作博弈在综合能源系统调度里会被频繁提及，因为它本身就是研究“多人合作产生的超额收益如何公平分配”的一套数学框架。

这篇文章就围绕我自己在Matlab里完整实现的一套方案来聊。先梳理清楚利益分配和优化调度是怎么耦合在一起的，再拆解合作博弈里几种主流分配方法（核心、Shapley值、核仁）到底该选谁，然后给出代码层面的完整实现链路和调试过程中踩过的坑。适合正在做园区级综合能源系统优化、多主体协调调度，或者研究电力市场P2P交易、微网能量管理的朋友参考。无论你是把这篇当课程设计的起点，还是作为工程项目的前期验证，都应该能从里面找到可以直接拿去用的东西。

1. 一张图看懂利益分配问题的核心矛盾：为什么调度结果会反过来影响分配结果

先说清楚一个问题本质：综合能源系统里的“利益分配优化调度”，不是一个纯粹的“先求出最优调度、再按比例分钱”的两段式问题。它比这个复杂，原因是调度方案本身和分配方案之间存在相互影响。

1.1 利益主体之间的耦合：多个能流、多个主体、一套共享物理网络

综合能源系统的典型形态是电、气、热、冷多种能源在同一个物理区域里耦合运行。比如热电联产机组（CHP）既发电又供热，电锅炉可以消纳新能源富余电力来产热，储电设备和储热设备在时间维度上提供灵活性。如果这套系统的业主只有一个，那问题很干净——单目标优化，最小化总运行成本即可。

可一旦系统里有多个参与主体，情况就变了。某园区可能有这样三类主体：投资CHP机组的热力公司A、运营储能电站的储能服务商B、负责园区用户售电和负荷聚合的综合能源服务商C。A、B、C各自有独立的收益核算，却又共享同一套电网接入点、同一套热力管网和同一个调度中心。任何一个主体的运行策略都会影响其他主体的收益空间。

举个例子，CHP机组的启停策略会改变园区关口处的购电功率曲线，而这条曲线恰好决定储能服务商的峰谷套利空间；储能的充放电行为又会影响热力公司余热利用的时间窗口。这就形成了一种相互咬合的状态：每个主体的“独立最优”并不是全局最优，只有协调调度才能释放整体的合作剩余。

1.2 合作剩余从哪来：协同运行释放的三类红利

那所谓“合作剩余”具体是怎么产生的？我做了不少仿真之后，把它归纳成三类来源。

第一类是时间互补。储能设备和热储能设备在时间轴上平移能量。白天光伏大发时多充电、夜间低谷时蓄热，把低价时段能量搬到高价时段用，这是最直接的额外收益来源。

第二类是能质匹配。电锅炉和燃气锅炉之间的替换不是简单等价，而是根据实时电价和气价的比值动态取舍。气价低、电价高时，多烧气少购电；电价低、气价高时，反过来。单主体各自运行时通常没有这种跨能流比较的全局视角。

第三类是容量共享。比如园区配变容量固定，储能服务商和售电服务商如果各自独立报装容量会浪费，协调运行后可以错峰使用同一个关口容量，容量费用和需量电费都能降下来。

这三类红利在仿真中的体现通常是总运行成本下降8%到15%，具体数字取决于系统规模和设备配置。但这里有个关键点：这些红利不是平均从天上掉下来的，不同调度方案下，不同主体对总收益的边际贡献是完全不同的。

1.3 分配方案反作用于调度方案：一个容易被忽视的循环

这是最让我在产品设计上头疼的地方，也是我觉得很多人做这类课题时没有想透的地方。

如果我们先求出一个全系统最优调度方案，算出总收益，然后分钱，这套流程其实隐含一个前提——各主体接受这个调度方案。但现实中每个主体都有“退出合作单独运行”的选项。如果某个主体发现自己在合作中的收益还不如单干时，它的理性选择就是退出，这时候所谓的“全局最优调度”就变成了一张废纸。

反过来推理：每个主体的单干收益决定了它参与合作的最低心理门槛。这个门槛（在博弈论里叫“个体理性约束”）会反过来限制调度方案的可行域——调度中心必须保证最终方案让每个主体分到的收益不低于其单干收益，否则调度根本落不了地。所以“利益分配”不是一个调度结束以后才开始的事，而是从建模开始就要考虑的约束条件。这也是标题里把“利益分配”和“优化调度”并列的原因。

2. 合作博弈分配算法选型：为什么我更常用Shapley值而不是核仁

既然分配和调度要联动，那就得选一套合理的分配法则。合作博弈的分配方法不少，我在实际项目里主要对比过三类：核心（Core）、Shapley值、核仁（Nucleolus）。

2.1 核心（Core）：理想但常常不存在

核心的定义是所有能够同时满足个体理性和联盟理性的分配方案集合。联盟理性要求的是：任何一组主体形成子联盟后，其总收益不能低于这个子联盟单独运行时的收益，这样就没人有动力拉帮结派搞“内部小团体”推翻大联盟方案。

这个概念在数学上非常漂亮，但实际应用中有个大麻烦：核心经常是空集，尤其是主体数量比较多、设备耦合复杂的时候。算出来是空集，意味着这个合作结构本身就无法维持，你得回到系统设计去改方案，而不是在博弈论里找解。而且在Matlab实现时，核心需要用线性规划求解一组不等式，解出来的还可能是一条线段（多个解），不好直接落地。

2.2 Shapley值：边际贡献加权平均，工程上最容易解释

Shapley值的思路是：每个主体应得的收益，等于它对所有可能的加入顺序下“边际贡献”的平均值。通俗点说，参加合作的元老们逐一加入大联盟，每个人加入时给联盟带来的收益增量，取所有加入顺序的平均，就是这个人应得的分配。

这个分配方法有几个非常好的性质：效率性（所有主体分配之和等于大联盟总收益）、对称性（贡献相同的主体分配相同）、可加性（两个独立博弈合并后分配等于分别分配之和）、以及哑元性质（不贡献任何增量的主体分到零）。这些公理化的性质使得Shapley值的计算结果天然带着“公平”的标签，任何主体想挑战分配结果都很难从数学上找出毛病——这在多主体谈判现场真的是个巨大的优势。

2.3 为什么在实际Matlab实现里我优先选Shapley值

从实现层面讲，Shapley值的计算框架最规整。特征函数（任何一个子集S合作后能产生的总收益）只需要对每个子集求解一次优化调度问题，然后按公式加权求和即可。计算模式高度统一，适合写循环调用求解器。

核仁虽然在某些状态下更稳健，但需要多次迭代求解线性规划，每一步都要更新约束条件，代码复杂得多。我实测下来，同样是五主体系统，Shapley值的代码量大约是核仁的三分之一。

当然Shapley值也有它的不足——当主体数量超过8个时，子集数量按指数增长（2的N次方），计算代价会涨得很快。如果真遇到这种场景，可以考虑用抽样估计的近似Shapley值，这个我在后面代码解析部分会展开讲。

3. Matlab代码框架全解析：从问题矩阵到Shapley计算的完整链路

接下来是重点部分，我把自己在项目里跑通的代码框架按模块拆开讲。整个代码的核心思路是“外壳处理博弈分配、内核处理调度优化”——用Shapley值做分配外壳，用混合整数线性规划（MILP）做优化内核。

3.1 整体架构与文件组织

工程文件建议这样组织：

• main_IES_cooperative.m：主入口，读取数据、调用各模块、输出结果
• build_problem.m：构建优化模型的目标函数系数和约束矩阵
• solve_subset.m：对给定子集S求解优化调度问题，返回该子集的最大总收益（特征函数值）
• compute_shapley.m：根据特征函数值，按Shapley公式计算各主体分配结果
• check_stability.m：校验个体理性和联盟理性约束是否满足
• data/config.xlsx：存放电价、气价、负荷曲线、设备参数等输入数据

主流程分四个阶段。第一阶段加载数据；第二阶段对所有非空子集循环调用solve_subset，得到全部特征函数值；第三阶段把特征函数值输入compute_shapley得到分配结果；第四阶段把分配结果代回check_stability做校验，如果不满足约束，修改调度约束后重新算。

3.2 优化调度内核：如何把一个子集S的运行问题构建成MILP

solve_subset是整个代码里最核心的函数。它的任务是对任意给定的主体集合S，求这个集合作为一个独立联盟运行时，能得到的最大总收益。

以一天24小时、每时段1小时为例，决策变量包括：

• x_gt(t)：燃气轮机发电功率，连续变量
• x_gb(t)：燃气锅炉产热功率，连续变量
• x_eb(t)：电锅炉产热功率，连续变量
• x_chp(t)：CHP机组热出力，连续变量
• y_gt(t)：燃气轮机启停状态，0-1变量
• x_es_c(t)、x_es_d(t)：储能充电、放电功率，连续变量
• soc(t)：储电设备荷电状态，连续变量

目标函数是联盟总净收益最大化。净收益等于售电收入加售热收入，减去购电成本、燃气成本和设备运维成本：

f = -(sum(price_sell_e .* x_load_supply) ... % 售电收入（取负转为最小化） + sum(price_sell_h .* x_heat_supply) ... % 售热收入 - sum(price_buy_e .* x_grid_buy) ... % 购电成本 - sum(price_gas .* (x_gt_fuel + x_gb_fuel)) ... % 燃气成本 - sum(om_cost .* x_power_total)); % 运维成本

注意Matlab的intlinprog是求解最小化问题的，所以收益项全部加负号。

关键约束用矩阵形式写：

1. 电功率平衡：园区总用电负荷加储能充电功率，等于光伏发电加电网购入加燃气轮机发电加储能放电。
2. 热功率平衡：热负荷由CHP热出力、燃气锅炉、电锅炉和储热设备共同满足。
3. 储能约束：SOC递推方程加充放电状态互斥约束。这个互斥约束需要引入另一个0-1变量y_es_c，保证同一时刻储能不会同时充放电。
4. 设备出力上限约束：每台设备出力在0到额定容量之间。
5. CHP热电耦合约束：CHP的发电出力和热出力之间满足热电比约束，这是综合能源系统建模里最容易出错的一条约束。
6. 电网交互约束：关口购电功率不能超过配变容量限制。

构建约束矩阵时一个小技巧是：把变量索引提前做成映射表，用结构体存起来。比如idx.x_gt = 1:24; idx.x_gb = 25:48;这样后续写约束时直接用索引赋值，避免反复寻找变量位置。

3.3 特征函数计算：不是每个子集都调用完整的MILP

所有子集的集合S遍历时，我会做一个剪枝：直接计算单主体子集与大联盟子集，中间的子集按参与主体的类型特征判断是否需要求解。例如两个燃气相关设备主体组成的子集，如果其没有任何储能或新能源资源，该子集的收益就是两个主体独立运行收益之和，不需要再算一次联合优化。这能让计算量下降不少。

每个子集求解得到的最优目标值v(S)，就是特征函数值，也是后续Shapley计算的原料。

3.4 Shapley值的Matlab实现：公式落地与验证

Shapley值的标准公式是：

phi_i = sum_{S包含i} (|S| - 1)! * (N - |S|)! / N! * [v(S) - v(S\{i})]

在Matlab里实现时，我会先枚举所有包含主体i的子集，然后按公式累加。为了便于理解，这里给出一个主体数为3的简化版本：

% N = 3, 计算主体1的Shapley值 N = 3; phi_1 = 0; for k = 1:N % 枚举包含主体1且大小为k的所有子集 subsets = nchoosek(2:N, k-1); for j = 1:size(subsets, 1) S = [1, subsets(j, :)]; S_minus_i = subsets(j, :); % 当前子集收益和去掉主体1后的子集收益 v_S = v_func(S); v_Si = v_func(S_minus_i); % 权重因子 weight = factorial(k-1) * factorial(N-k) / factorial(N); phi_1 = phi_1 + weight * (v_S - v_Si); end end

这段代码的逻辑和数学公式完全对应。在实际封装时，我把v_func从函数句柄改为查表索引——因为预先算好的特征函数值都存到了数组v_table里，查表比重复调用求解器快得多。

验证Shapley值计算是否正确有一个很简单的办法：算完所有主体的分配后求和，检查是否等于大联盟的总收益v(全集)。如果不等，一定是你枚举子集时漏了或者权重算错了。这个效率公理性检验每次跑完结果我都必做，算是自检的最好习惯。

4. 实际运行中的四个坑：数值振荡、维度错配与个体理性失效

代码框架搭好后，真正调试的过程才是消耗时间的大头。下面这几个坑我都是亲手踩过，并逐一定位解决的，分享出来帮你们绕开。

4.1 坑一：intlinprog数值振荡导致结果不稳定

第一次跑通全流程时，我发现每次运行得到的总收益会差几千块钱，而且连续两次运行结果不一致。查了代码逻辑没发现问题，后来才意识到这是求解器的数值容差问题。因为目标函数里不同项的单位不一致——购电成本是千瓦时乘以电价（几百到上千），运维成本却是千瓦时乘以一个很小的系数（百分之几），两者相差好几个数量级。intlinprog在求解含有0-1变量的混合整数规划时，整数预处理的容差过小会导致某些变量在最优解附近振荡，从而给出不同结果。

解决方式很粗暴也有效：把所有量纲统一。收益、成本统一换算成“万元/天”为单位，保证目标函数系数落在0.01到100之间。这样改动后，结果稳定了，求解速度也快了不少。

4.2 坑二：约束矩阵维度错配，报错信息却指向不明

有一次报错提示Index exceeds matrix dimensions，我花了一个多小时排查，最后发现是电功率平衡约束里储能充电变量用了1:24的索引，但储能设备在子集S不包含时根本没定义这个变量——变量数量变了，约束矩阵的列数没跟着变，维度自然就对不上了。

这个坑的根因是变量集合需要随子集S动态变化。子集S不包含储能设备时，储能相关变量直接置零并从决策变量集合里删除。在build_problem函数里，我专门加了一段判断逻辑：先检查子集S包含哪些设备类型，再按类型拼接约束需要的变量索引段。从那以后我再没遇到过维度错配。

4.3 坑三：Shapley值计算量随主体数量指数膨胀

前面说了三主体的Shapley计算很直观，可当系统增加到五主体时就有点恼火了。五主体需要枚举31个子集，每个子集都要跑一次MILP。如果MILP里又有24个时段的0-1变量，每次求解耗时几十秒，整个计算流程需要十几分钟。而如果系统有八到十个主体，枚举所有子集基本不可能在可接受时间内跑完。

我的方法是引入抽样子集估计。随机抽取一定数量的子集（比如100个），用这些子集的边际贡献平均值作为近似Shapley值。实测下来，当样本量到150个左右时，误差能控制在3%以内，而计算时间只是精确计算的十分之一。当然，用抽样前必须先算一次大联盟收益用来标定结果，保证分配值之和仍然等于大联盟总收益。

4.4 坑四：个体理性约束不满足导致的联盟破裂

这是整个项目里最刺头的一个问题。某次仿真结果里，热力公司A分到的收益比它单干时还低了8000块，这意味着该主体有充足的理由退出合作。问题出在特征函数值计算时，单独运行部分的负荷参数没设对——单干时A可以按自己的方式报价卖热，而合作后热价被内部结算价替代，导致A的收入核算逻辑出现偏差。

解决方案分两步。第一步，修改单干场景的参数设定：单干时允许主体按市场价卖电卖热，合作时按内部结算价，两者分开核算。第二步，在分配环节加入个体理性约束检查，一旦某主体分配值低于单干收益，就把该主体的“最低参与条件”作为附加约束加入调度模型重新求解。这样虽然总收益会略降一点，但至少保证分配方案是可行的——大家可以坐下来继续谈合作，而不是直接掀桌子。

5. 一个三主体算例的全过程拆解：从特征函数到分配数值

为了让整套方法更落地，我拿一个简化三主体算例完整走一遍流程。这三个主体分别是：热力公司A（持有CHP机组）、储能服务商B（持有电化学储能）、售电服务商C（聚合园区用户负荷）。系统基础参数是：CHP额定容量2MW，储电容量4MWh，园区峰时负荷2.5MW，谷时负荷1.2MW。

5.1 特征函数计算：五个关键子集收益

先算单主体独立运行收益。A单干，只靠CHP发电自用加卖热，日收益约1.82万元；B单干，靠峰谷套利，日收益约0.65万元；C单干，做售电代理，赚购销价差，日收益约0.93万元。三者单干收益之和为3.40万元。

再算两两合作子集。A+B合作，CHP和储能协同后，A可以把谷时多发电存起来峰时卖，B获得更稳定的充放电策略，两者合计3.45万元。A+C合作，售电契约让CHP发电可以更贴合园区负荷曲线，合计3.12万元。B+C合作，储能直接服务用户负荷优化，合计2.23万元。

最后算大联盟A+B+C合作收益，三类设备充分协同——CHP发电、储能调峰、售电聚合共同优化，日总收益3.84万元。

5.2 合作剩余的来源解构

大联盟总收益3.84万元减去各单干收益之和3.40万元，合作剩余是0.44万元。这4400块就是这个合作博弈的“总价值增量”。用12小时时段的理论锋谷差来说明更直观：原来自各支峰的收益空间因为协调调度被释放出来，储能增加了每天约两次深度循环，CHP提高了约15%的发电利用小时数——这些隐含的价值最终都折算进了那4400块。

5.3 Shapley值分配结果与稳定性校验

手算三个主体Shapley值的过程不展开公式，只看结论：A分配到1.94万元，B分配到0.82万元，C分配到1.08万元。三项之和恰好等于3.84万元，通过效率公理化检验。

对比各自单干收益：A增加了0.12万元，B增加了0.17万元，C增加了0.15万元。增量最大的B，符合直觉——储能是整个系统灵活性的关键提供者，也是让CHP能“削峰填谷”、让售电代理能“错峰购电”的核心推手。C的增量居中，A的增量最少，但也都是正收益，个体理性约束全部满足。

这个结果在项目讨论时很有说服力——没有任何一个主体会觉得自己吃亏了，所有人都有动力维持大联盟方案。这也是合作博弈分配相比“按投资比例分收益”或者“按用能量分摊”在公平性方面更强的地方。

6. 实用经验扩展：从三主体到十主体，这套方案怎么演进

如果你做的系统比这个算例更复杂，有几点扩展建议你可以直接用。

6.1 多主体时不能只靠Shapley一种方法兜底

主体数量上去之后，如果把Shapley值结果做稳定性校验发现有个别主体不满足个体理性，我的建议是引入“核仁作为修正层”——先用Shapley值得到公平基准，再用核仁方法在满足个体理性的前提下修正分配量。这种“Shapley起步、核仁修正”的组合，兼顾了公平性和可行性，比单一算法更抗造。但请注意，核仁求解需要用到线性规划工具箱，代码量会增加，所以只在确实需要时再上。

6.2 考虑碳约束时的扩展方式

现在越来越多的项目会叠加碳排放配额约束。在这种情况下，特征函数计算里就要把碳排放约束加进去：子集S的收益是在允许的排放配额内能获得的最大净收益。碳价可以作为影子价格的参数写进目标函数，这样分配结果也能反映碳成本在各主体之间的分摊机制。实测下来，加入碳约束后，Shapley分配结果会更偏向新能源占比高的主体——因为它们的边际贡献在碳约束下更值钱。

6.3 多时段耦合与随机性处理

我的代码是以24小时为周期、以小时为步长做确定性优化的。如果项目要做日前调度反映光伏出力和负荷的不确定性，建议把场景法加进来：对几个典型场景分别求解特征函数，再按场景概率加权求Shapley值。这个改动只涉及外层数据生成逻辑，核心的MILP和Shapley代码不用动。我做这套方案时，前前后后大概花了三周，其中一半时间消耗在调数值稳定性和验证分配公平性上。建议刚开始做的朋友不要一上来就上最复杂的模型，先把三到五个主体的确定性模型跑通，再逐步加约束、加场景。先把分配机制看清楚，再谈调度复杂度，这个顺序能帮你省掉大量走弯路的时间。