企业官网建设流程全解析

1. 项目概述与核心价值

1.1 这个项目到底解决什么问题

先说个大白话版本。综合能源系统这个词听起来挺吓人，但本质上就是把电、热、气、冷这些不同形式的能源放在一张网里协调着用。比如一个园区里既有燃气轮机发电，又有余热锅炉供热，还有储能电池和光伏板，各设备单独运行往往浪费严重——燃气轮机发电时排出的高温烟气直接扔掉太可惜，但你要让热力和电力部门自己商量怎么配合，又会出现利益纠纷。

这就是我做这个项目的初衷：用合作博弈论把经济利益分配这件事做成一个可计算的数学模型和优化调度程序。合作博弈处理的不是设备怎么运转，而是“谁该拿多少钱”的问题。很多做技术的人容易忽略一个事实：系统优化调度做得再好，如果收益分得不公平，联盟第二天就散了。所以这个项目把“怎么调度”和“怎么分钱”两条线拧在了一起，用Matlab完整实现。

从实际价值上看，这套方法主要面向三类人：第一类是搞综合能源系统规划的工程师，需要给多主体联合运行方案算经济账；第二类是研究博弈论应用的研究生，需要一个能改参数、能跑出图表、能对比不同分配方法的Matlab框架；第三类是园区能源管理系统的设计者，虽然不一定直接在代码层级用合作博弈，但需要理解多主体利益协调的底层逻辑。这篇文章我就把项目的建模思路、Matlab代码结构、核心算法实现细节和调试心得全部展开讲一讲。

1.2 为什么利益分配比调度本身更棘手

如果你做过能源系统优化，大概率熟悉这样的套路：建立目标函数、写约束条件、调一个求解器、给出最优出力曲线。这道流程走下来，得到的是一个“总成本最小”或“总收益最大”的方案。问题是，当这个系统里同时存在多个利益主体——比如燃气轮机属于A公司，储能属于B运营商，光伏属于C业主——你算出的那个“最优总收益”并不能直接落地，因为谁也不愿意在自己吃亏的情况下去成就“整体最优”。

举个更直观的类比。三个人合伙开一家餐馆，厨师负责做菜，服务员负责接待，采购负责买食材。三个人单独干都能赚钱，但合伙干效益更高。现在问题来了：合伙干多赚的那部分钱，到底怎么分才合理？是按工作时长？按贡献大小？还是按“如果没你，别人能赚多少”的差额来定？合作博弈里的Shapley值就是回答这个问题的经典方法，它的核心思想是：每个人拿到的收益，等于他对所有可能的合作组合带来的边际贡献的平均值。

2. 综合能源系统的数学建模与场景设计

2.1 典型系统结构与参与者划分

我做这个项目时，设计了一个典型的三主体综合能源系统，结构如下：

• 主体A：热电联产机组（CHP）。燃气轮机发电，同时通过余热回收装置供热，是系统的主要能量来源。
• 主体B：风电运营商。拥有风机，出力具有随机性，但边际成本极低。
• 主体C：储能运营商。拥有电储能设备，可以在电价低谷时充电、高峰时放电，同时也能配合风电消纳。

三者独立运行时各自为政，合作运行时的潜力很明显：CHP可以根据风电出力灵活调整自身发电量，避免弃风；储能可以在风电多发时充电，在负荷高峰时段放电，替代CHP的高成本出力；CHP的余热直接供给热负荷，而储能可以通过电转热设备辅助供热。

为什么要选这三个主体？因为它们具备典型的互补性特征。如果主体之间完全没有互补性，合作博弈算出来的分配结果会退化到“各干各的”，这个项目就失去意义了。风电的不确定性需要储能来平抑，CHP的发电余热需要热负荷来消纳，三者正好构成一个闭环。

2.2 目标函数与约束条件的标准化表述

项目采用日前调度模式，即提前一天根据预测数据制定24小时调度计划，时间分辨率取1小时。决策变量包括：CHP的发电功率和产热功率、储能的充放电功率、购电功率以及与电网的交互功率。

目标函数是合作联盟的总运行收益最大化，表达式如下：

max F = Σ [ p_e(t) * P_grid(t) + p_h(t) * Q_h(t) - c_gas * V_gas(t) - c_om * P_chp(t) ] - c_battery * P_dis(t) + p_elec(t) * P_wind(t)

其中，每个变量代表的含义是：p_e(t)和p_h(t)是电网电价和热价，P_grid(t)为与电网交互的电功率，Q_h(t)为CHP的供热量，c_gas为天然气单位热值成本，V_gas(t)为天然气耗量，c_om为机组运行维护成本系数，P_chp(t)为CHP发电功率，c_battery为储能单位放电成本，P_dis(t)为储能放电功率。

约束条件涵盖以下几类：

• 功率平衡约束：发电量 + 风电出力 + 储能放电 + 购电量 = 电负荷 + 储能充电 + 售电量。
• CHP运行约束：发电功率上下限、爬坡率限制、热电比运行区间（这是CHP建模里最容易写错的地方，后面重点讲）。
• 储能约束：充放电功率限制、SOC（荷电状态）连续性约束、容量上下限。
• 交互功率约束：与电网的最大购售电功率限制。

这一类模型用Yalmip写起来非常方便，核心代码大致如下：

% 定义决策变量 P_chp = sdpvar(T, 1); % CHP发电功率 H_chp = sdpvar(T, 1); % CHP供热功率 P_dis = sdpvar(T, 1); % 储能放电 P_ch = sdpvar(T, 1); % 储能充电 SOC = sdpvar(T+1, 1); % 荷电状态 P_grid = sdpvar(T, 1); % 电网交互功率（正为购电，负为售电） Constraints = [SOS_balance: P_chp + P_wind + P_dis + P_grid == P_load + P_ch]; Constraints = [Constraints, P_chp_min <= P_chp <= P_chp_max]; Constraints = [Constraints, -ramp_rate <= P_chp(2:end) - P_chp(1:end-1) <= ramp_rate]; Constraints = [Constraints, SOC(2:end) == SOC(1:end-1) + P_ch * eff_ch - P_dis / eff_dis]; Constraints = [Constraints, SOC_min <= SOC(2:end) <= SOC_max]; Objective = -sum(price_profile .* P_grid - gas_cost .* V_gas); optimize(Constraints, Objective, sdpsettings('solver', 'gurobi'));

注意CHP的热电比约束。很多初写CHP模型的人直接把电出力和热出力分别限定了上下限，但实际机组的热电比是存在一个可行域的——出力大的时候热电比偏向低值，出力小的时候热电比反而高。不加这个约束，优化结果可能让CHP在低电出力下产生超物理限制的热量，分钱时就会出问题。

2.3 为什么独立运行和合作运行的收益差这么大

做合作博弈之前，先要算出三种情况的数据：

1. 独立运行：A、B、C各自以自己的目标函数和约束求解，互不通信。
2. 部分联盟：AB合作、AC合作、BC合作，各自内部优化但整体独立于第三方。
3. 大联盟：三个主体共同优化，即上面那个统一目标函数。

计算独立运行收益是必须的一步，因为它是合作博弈中做“边际贡献”对比的基准线。我实际测算时发现，独立运行模式下CHP由于没有储能的配合，为了满足夜间低谷时段的电负荷，不得不在低效工况区运行，而且为了满足热负荷需求，电出力被热出力“绑架”，无法灵活调整——它的单位发电成本明显高于大联盟模式。风电独立运行时遇到大发时段只能弃风，因为电网购售电协议限制了反送功率的极限。储能独立运行时则是“低买高卖”赚电费差价，完全不知道还有“帮风电消纳、帮CHP调峰”这种更高价值的玩法。

合并系统后，CHP的灵活性大幅提升：白天负荷高时满发，晚间负荷低时降出力，不足的部分由风电和储能顶上，而热负荷需求可以由余热锅炉和储热协同满足。风电大发时段储能充电消纳，小发时段储能放电支撑负荷。三者配合的结果就是，大联盟总收益比三个独立运行收益之和高出大约18%到25%（具体数字取决于机组容量配比和负荷曲线，但方向是一致的）。多出来的这部分就是合作剩余，是博弈分配的对象。

3. 合作博弈分配方法的原理与实现

3.1 Shapley值分配法：按边际贡献分钱的经典方案

Shapley值的物理含义很直观：先计算每种联盟组合的收益，再通过排列组合计算每个参与者的边际贡献期望。公式长这样：

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

理解这个公式需要一点排列组合基础。|S|是子联盟S中的人数，N是总人数。v(S∪{i}) - v(S)表示当i加入联盟S时，联盟总收益的增加量。为了让所有可能的加入顺序都被考虑，前面乘以一个权重系数——它代表在所有排列中，参与者i恰好排在S后面加入的概率。

在Matlab中实现Shapley值计算时，我采用了矩阵化的写法而不是三重嵌套循环。效率提升非常明显，尤其当参与者数量从3扩展到5或6时：

function phi = shapley_value(v, N) % v: 向量，长度为 2^N，包含所有联盟的收益值 % 其中 v(1) 是空联盟收益，v(end) 是大联盟收益 % v(2) 代表参与者1单独收益，v(3) 代表参与者2单独收益... phi = zeros(N, 1); for i = 1:N total = 0; for S = 0:(2^N - 1) % 检查参与者i是否在联盟S中 if bitget(S, i) continue; end prev_v = v(S + 1); % 当前联盟S加上参与者i后的编号 new_v = v(bitor(S, bitshift(1, i-1)) + 1); marginal = new_v - prev_v; s_size = sum(bitget(S, 1:N)); weight = factorial(s_size) * factorial(N - s_size - 1) / factorial(N); total = total + weight * marginal; end phi(i) = total; end end

这段代码的关键点：用bitget和bitor做位运算来表示联盟集合。经常有人问为什么用位运算不用数组——因为3主体时还有穷举的价值，但到6主体时联盟数已经达到2的6次方64个，如果用数组存联盟状态，代码会变得笨重而且很难扩展维度。

Shapley值的最大优点是有唯一性和公平性公理支撑：对称参与者分配相等、有效分配、匿名性、可加性。它在理论上无懈可击，但计算代价是O(N!)，参与者超过7个时基本没法全枚举，需要用抽样估计算法来处理。

3.2 核仁法分配方案：最小化不满情绪的博弈解法

核仁法走的是另一条路线。它的目标不是“按边际贡献平均”，而是最小化联盟中最大“不满程度”。定义联盟S的自己，即S的成员认为自己在分配中应该拿到的最低收益为v(S)，而当前分配方案给S的总收益是Σφ_i(S)，那么联盟S的剩余为：

e(S, x) = Σ φ_i(S) - v(S)

核仁法的逻辑是：在所有满足整体均衡的分配方案中，找到一个方案，使得最不满的联盟的剩余尽量大——也就是尽量缩小每个联盟“应得”和“实得”之间的差距。

实现核仁法要用到线性规划的分层求解。第一轮求解目标函数是最大化最小的e(S, x)，得到第一个β1；第二轮在保证第一轮目标不被破坏的前提下，最大化次小剩余，以此类推。这个过程本质上是求解一系列双层线性规划，在Matlab里可以直接调用linprog完成，但中间需要维护一系列约束的累积。

有人会说核仁法太“平均主义”了，但在某些场景下它确实比Shapley值更贴近实际。比如当一个主体对联盟贡献不大但退盟损失惨重时，Shapley值会分给它很少的钱，这就有可能导致该主体拒绝合作——核仁法会多给它一些补偿，保证联盟稳定。换句话说，Shapley值按贡献分配，核仁法按稳定性分配。

3.3 纳什谈判解：一种基于议价能力的替代方案

第三种方法是纳什谈判解。它不再从“边际贡献平均值”的角度出发，而是从谈判起点出发。假设各主体独立运行时的收益为x0_i，也就是谈判破裂时的保底收益。纳什谈判模型求解的是以下优化问题：

max Π (x_i - x0_i) s.t. Σ x_i = V(N) x_i >= x0_i

其中乘积最大化代表帕累托效率与个体理性的叠加。

这个问题的几何意义很有趣：可行域是分配多面体，纳什谈判解是处于“帕累托前沿”和“个体理性约束”交界处的点，它在保持合作稳定性的同时尽量扩大了整体收益的乘积。在Matlab中实现也不难，用fmincon配合对数变换即可：

% 目标函数变换成 sum(log(x_i - x0_i)) obj_fun = @(x) -sum(log(x(1:3) - x0)); % 约束条件 Aeq = [1 1 1]; beq = V_total; % 总收益平衡 lb = x0 + eps; ub = [V_total, V_total, V_total]; x_nash = fmincon(obj_fun, init, [], [], Aeq, beq, lb, ub);

表格式对比三种分配方法效果如下：

3.4 三种分配方法的Matlab实现框架

项目中我把三种方法封装到一个统一的分配模块中，输入是各联盟的收益值，输出是分配结果和相关指标。核心代码结构如下：

function allocation = allocation_module(v, method, options) % v: 所有联盟收益向量 % method: 'shapley', 'kernel', 'nash' N = options.N; switch method case 'shapley' allocation = shapley_value(v, N); case 'kernel' allocation = kernel_nucleolus(v, N); case 'nash' allocation = nash_bargaining(v, N, options.x0); end % 后处理：检查个体理性、整体理性、联盟理性 if sum(allocation) ~= v(end) warning('分配结果与总收益不一致，请检查约束'); end if any(allocation < options.x0) error('分配结果不满足个体理性约束'); end end

写这个模块时要注意一个细节：浮点数精度问题。位数运算生成的联盟编号和权重乘法很容易引入微小误差，导致分配结果加起来不等于大联盟收益。建议分配结果输出前做一次归一化修正：

allocation = allocation * v(end) / sum(allocation);

不然到了后续展示阶段，柱状图里各个模块收益加总比总收益少几毛钱，检查数据的时候很容易让人误以为代码有bug。

4. 优化调度与利益分配的联动实现

4.1 双层优化架构：先调度、后分配

处理“调度+分配”问题，我采用的是顺序式、非迭代的双阶段架构。第一个阶段求解大联盟的收益最大化调度方案，得到机组出力和能源流动数据；第二个阶段把第一阶段的调度结果代入联盟收益向量，基于合作博弈方法完成利益分配。

为什么不把调度和分配一次性建模成一个双层规划？原因很简单：合作博弈的收益函数v(S)本身就是从优化问题里求出来的，把它嵌套进分配公式会让整个模型变成复杂的非线性混合整数规划，求解极其困难。实际中我采用的这种“算完再分”的做法已经是学界和工程界的主流选择。单层求解器如果设定不好很容易陷入局部最优，且调试难度加倍。

4.2 联盟收益向量v(S)的自动化计算

整个项目里最消耗时间但最机械化的部分，就是计算所有子联盟的收益。3个主体需要计算2的3次方=8个联盟的收益（包括空联盟和大联盟）。每个联盟都要运行一次优化调度，也就是8次求解。推到4个主体就是16次，5个主体32次——每多一个参与者，计算量翻倍。

代码层面我用一个循环来处理：

alliances = {}; for k = 1:N combos = nchoosek(1:N, k); for j = 1:size(combos, 1) active = combos(j, :); % 调用优化调度函数，只激活该联盟内的设备和负荷 [obj, result] = run_scheduling(active, system_data); alliance_id = sum(2.^(active - 1)) + 1; v_arr(alliance_id) = obj; end end

这里特别注意run_scheduling函数内部要对不参与联盟的设备做“代码级关闭”处理。如果设备不参与某个联盟，它的出力要被固定为0或强制退出运行状态，而不是靠约束边界去压制——因为混合整数模型里约束边界压制经常导致求解缓慢，还可能出现退化解。

4.3 调度结果的可视化与敏感性分析

Matlab做可视化是强项，这个项目我生成了三张经典图表：第一张是24小时各机组出力堆叠图，能直观看到合作前后CHP出力曲线的形状变化；第二张是各联盟收益柱状图，用来展示三种分配方法下各主体收益的差异；第三张是风速和电价波动条件下三种分配方法的对比曲线，做敏感性分析用。

做敏感性分析时我固定了系统参数不变，把风速最大出力从5MW调到10MW步长1MW，看看各分配比例怎么变。结论很有意思：Shapley值分配的份额不随风速变化而剧烈波动，而纳什谈判解的分配结果对谈判保底收益设定极其敏感。如果保底收益设得过高，纳什谈判解会直接退化成平均分配；设得太低，又会主导联盟内的话语权偏差。这个现象在论文里可以当亮点讲，在工程上则需要警惕——参数计算的鲁棒性，直接影响分配的稳定性。

5. 常见问题与调试心得

5.1 求解器选择：Gurobi还是Cplex还是内置的linprog

做混合整数线性规划时我一直用Gurobi，一方面是因为它在电力系统优化领域确实快，另一方面是Yalmip接口支持得最好。Cplex在部分分配类问题上表现也好，但遇到大规模混合整数规划时，Gurobi的MIP gap下降速度明显更快。

如果只有Matlab基础工具箱，即使用linprog也能走通，但求解速度会慢很多，尤其在联盟收益计算阶段——8次求解的累计时间会到分钟级别。如果坚持使用内置求解器，建议将模型改成纯线性规划，不要引入二进制设备启停变量，否则每天调度任务根本跑不完。

5.2 Shapley值计算组合爆炸：N变大后怎么办

当参与者数量突破7个，Shapley值的全枚举计算会卡到没法用。项目里我做了抽样近似版本：

function phi_approx = shapley_approx(v, payoff_func, N, sample_count) phi_approx = zeros(N, 1); for s = 1:sample_count perm = randperm(N); cum_v = 0; for i = 1:N current_player = perm(i); S = perm(1:i-1); val_with = payoff_func(sort([S, current_player])); marginal = val_with - cum_v; phi_approx(current_player) = phi_approx(current_player) + marginal; end end phi_approx = phi_approx / sample_count; end

这么做误差大概在2%到5%之间，但对于工程决策足够了。关键是，每次抽样调用payoff_func时，要向调度模块实时读取对应联盟的收益值。

5.3 热电比约束的“物理陷阱”

这是我实际写代码踩过的最深的一个坑。CHP热电比约束的正确写法是：

H_chp(t) <= r_max * P_chp(t) + b H_chp(t) >= r_min * P_chp(t)

不写b偏置量会出现什么问题？理想化的热电比约束会导致机组在小出力状态下被迫产热，优化求解器为了满足热负荷，会把电出力抬到不经济的水平。这样算出来的调度不可执行，随后算出来的联盟收益自然全都失真。模型校验后加一步可行性检验——随机给一组出力点看约束是否全满足，这类简单的小测试永远不要跳过。

5.4 联盟收益矩阵的对称性与重复计算问题

三个主体构成的联盟中，AB联盟和BA联盟本质上是同一个联盟，代码里用位运算自动会合并，但新手容易用元胞数组或字符串拼联盟ID，造成重复计算。解决方式：统一用整数位掩码表示联盟，1表示参与者A，2表示B，4表示C，那么AB就是3，AC是5，BC是6，大联盟是7。这样既方便存储，也方便调用时做成员判断。

5.5 数据驱动：负荷数据和电价序列的准备

项目需要用真实数据才能体现价值。我建议用以下来源：欧洲公开的能源市场电价数据、中国某园区典型日负荷曲线。这些数据目前都能在公开学术数据集和项目中找到，不过要注意时区、节假日处理。如果只用合成数据，算出来的结果很难说服审稿人或者领导。第一次跑通代码时可以用合成数据验证逻辑，正式实验时换成公开真实数据采集，两类数据跑出来的图表最好都保留，展示时可互补。

另外一个细节：电价序列必须按正负符号处理好购售电方向。很多初学者用正电价表示购电成本，用负电价表示售电收入，但在目标函数里容易符号写反，导致模型拼命购电而不是合理购电。我的建议是在数据预处理阶段就把电价拆成两个正变量——购电价和售电价，分别乘对应的功率项，这样逻辑清晰不容易错。

6. 项目扩展方向与个人实操体会

这套框架做完后，我最大的感觉是：博弈论方法落地的难度不在理论推导，而在建模尺度与数据粒度之间的匹配。Shapley值、核仁法这些概念说起来很漂亮，但真正要让它算出一个让各方信服的分配方案，中间隔着的是成千上万行代码、数不清的参数调整和无数次“这个约束怎么又退化”的抓狂时刻。

如果你打算把这个项目往后继续扩展，有几条路参考：

• 把多时段耦合度提高，从日前调度扩展到实时滚动调度，分配周期也随之变短。这种情况下，今天基于目前预测做的分配方案到了下午就可能需要动态调整，计算性能要求就上来了。
• 引入碳排放配额交易，把环境成本纳入收益函数。这样合作博弈不仅要分钱，还得分碳指标，分配对象变成了二维需求，求解难度和实用价值同时增加。
• 加不确定性建模，风电和负荷不只是用期望值，而是用场景树或机会约束，看看不确定性下的收益分配比例有何变化——这个方向学术发文的认可度很高，但代码复杂度会翻倍。

最后分享一个小技巧：项目里所有分配方法的结果最后都用一张横向柱状图汇总对比，每个主体三种颜色分别代表Shapley值、核仁法、纳什谈判解的分配金额。每次拿出这张图，无论对方是工程师还是管理层，沟通效率都远高于扔出几十页纯代码。优化结果的说服力，一半在计算精度，一半在可视化表达。这也是我做优化调度项目的一条总原则——让数据自己把话讲清楚。