最近在梳理综合能源系统方向的一个项目主题是“基于合作博弈的综合能源系统利益分配优化调度”。做多能互补或者园区能量管理的朋友应该都有体会单主体的优化调度很好写但一旦涉及多个主体共享同一个配电网、同一级天然气管网比如光储微网、燃气热电联产站、售电公司聚在一块儿各自只管自己的账整个系统的运行成本往往压不下来设备利用率也上不去。这个课题的核心思路就是用合作博弈来协调多个主体之间的关系把“各自为政”变成“先合作、再分利”最终落到一个可复现的Matlab求解框架上。这套方法论能解决的问题很具体多主体联合调度时合作带来的总成本节约怎么计算这个节约又该如何在各参与主体之间公平分配才能保证没有人觉得自己吃亏、联盟关系稳定可持续。适合两类人来参考一类是做园区综合能源、多微网协调控制方向的研究生和工程师另一类是刚接触博弈论在能源系统中应用、想把理论模型落地成代码的初学者。下面的内容都是我在实际算例上跑过的思路和代码按从建模到求解再到分配的顺序来讲。1. 项目整体定位与技术思路拆解1.1 这个课题到底在解决什么问题区域综合能源系统的典型场景是一个园区里有几个相对独立的能源主体比如A主体拥有光伏和蓄电池白天光伏出力充足但本地负荷小B主体是一套燃气轮机热电联产装置能同时产电产热但发电边际成本较高C主体是负荷聚合商手里有电负荷和热负荷对能量有刚性需求。如果把三者分开独立调度A会为了减少弃光在自己的小系统里强行给电池充满B会按自己的气价和电价判断开不开机C只能按实时市场价购电购热。问题就出来了晚间光伏归零A要高价从电网买电实际上B的热电联产如果开机可以更便宜地供电白天光伏大发的时段B可能因为电价低不愿意开机发电但C恰恰需要热导致B要用低效率的燃气锅炉单独供热储能设备的充放电策略只基于本地电价判断没有考虑相邻主体的需求曲线导致整个系统的峰谷调节能力被浪费。这些问题本质上是“信息割裂”和“利益冲突”造成的局部最优。合作博弈的思路则是先计算各个主体及任意子集形成联盟后的联合运行最优成本再把“单独运行成本之和”和“联盟运行成本”之间的差值——也就是合作带来的成本节约——通过某种规则分给每个参与者。1.2 为什么选合作博弈不选主从博弈或集中式优化很多人在做多主体能量管理时会先想到两类方法一类是集中式优化把所有设备、负荷、网络约束放进一个大模型里求全局最优另一类是主从博弈比如园区运营商做领导者用户或微网做跟随者。集中式优化解决不了隐私和公平性问题——各主体的内部数据、设备参数都要上交而且分配结果没法服众。主从博弈更适合层级不对等的场景比如售电商和终端用户。当参与者地位基本平等、大家都要出钱出力共同组成一个虚拟聚合体时合作博弈是最自然的选择。合作博弈和另外两者的核心差异在于它把问题拆成了“怎么联合运行更优”和“怎么分配收益更公平”两个先后阶段。这种先后逻辑符合园区实际落地流程先签合作协议、确立联合调度的运行方式再定期算账、按贡献分配收益。相比Nash非合作博弈中寻找稳定均衡点的做法合作博弈中的Shapley值、核仁等方法有明确的分配公式和线性规划求解路径代码实现难度低得多结果也容易向非专业背景的管理层解释。1.3 总体框架两层模型四个步骤我在这套代码里采用的框架不复杂分四步第一步定义参与者集合N{1,2,...,n}初始化各主体的设备参数、负荷曲线、风光出力序列、分时电价、气价数据第二步枚举所有非空子联盟S⊆N对每个联盟分别求解一个多主体联合经济调度优化模型得到联盟最优运行成本C(S)再据此计算其特征函数v(S)第三步基于特征函数v(S)用矩阵方式组织联盟值计算Shapley值或用线性规划求解核仁Nucleolus、判断核心Core非空性第四步用分配结果倒推各主体实际应分摊的成本结合各主体的独立运行成本计算净收益并校验个体理性、联盟理性、集体理性三个条件。值得强调特征函数不是随便定义的。常见做法是令特征函数等于单独成本减去联盟成本体现联盟创造的价值。我在后面的代码里用的就是这种定义即v(S) Σ_{i∈S} C_i^ind - C^coal(S)其中C_i^ind是主体i独立运行时的最优成本C^coal(S)是联盟S联合调度的最优成本。这样v(S)为正值时代表联盟创造的成本节约数值越大说明这个联盟协同带来的收益越明显。2. 多主体综合能源系统建模与特征函数计算2.1 典型三主体的设备配置与关键约束项目算例我设计为三主体已经能完整展现合作博弈的分配逻辑又不会因为子联盟数量太多导致枚举复杂度过高主体A光储微网含光伏机组、蓄电池储能、基础电负荷主体B热电联产站含燃气轮机、余热回收装置、燃气锅炉可同时输出电力和热力主体C综合负荷聚合商含电负荷、热负荷自身无发电设备只能通过电网购电和购买天然气。每个主体的关键设备模型如下光伏出力PPV,t PPV_rated × ηPV,t其中ηPV,t为归一化出力系数由当地光照强度决定短期调度中可当作已知参数。蓄电池储能模型SOC_t SOC_{t-1} ηch × Pch,t × Δt / E_cap - Pdis,t × Δt / (ηdis × E_cap)约束条件0 ≤ SOC_t ≤ 1Pch,t、Pdis,t不同时为正且有上下限限制。燃气轮机CHP机组模型P_gt,min ≤ P_gt,t ≤ P_gt,maxH_gt,t η_h/η_e × P_gt,t 热电比固定其中η_e为发电效率η_h为余热回收效率。如果机组是可调热电比的抽凝机组热电比会在一个区间内变化模型变为混合整数规划这里先按背压式固定热电比处理。燃气锅炉模型0 ≤ H_gb,t ≤ H_gb,maxH_gb,t η_gb × F_gb,t × LHV其中F_gb,t为耗气量m³LHV为天然气低热值。电锅炉模型0 ≤ P_eb,t ≤ P_eb_maxH_eb,t η_eb × P_eb,t电锅炉是电转热的枢纽也是B和C之间热力耦合的关键。2.2 电力与热力平衡约束每个联盟内部的联合调度核心是满足整个联盟的总电平衡和总热平衡约束电平衡P_buy,t Σ P_pv,t Σ P_gt,t P_dis,t Σ L_e,t P_ch,t P_eb,t热平衡Σ H_gt,t Σ H_gb,t Σ H_eb,t Σ L_h,t H_chs,t其中P_buy,t为从上级电网购电功率L_e,t为电负荷L_h,t为热负荷。这里有个容易踩坑的点多主体联盟和单主体孤岛运行时平衡方程的形式完全一致只是把各主体的变量聚合到一起。这也是为什么我建议把调度模型的求解写成函数形式输入是“参与者集合S内部的设备参数、负荷、价格”输出是该联盟的最优成本。这样枚举所有子联盟时就不用重复建模直接改参数调函数即可。2.3 联盟枚举与特征函数缓存策略对于n个参与者非空子联盟总数是2^n - 1。n3时只有7个联盟n4时是15个n5时是31个直接枚举完全没有压力。但如果参与主体数量到了7、8个以上联盟数量就变成127、255个每个联盟都要单独求解一个MILP模型计算时间会迅速膨胀。我在代码里做了两个优化用二进制掩码0-1组合表示子集避免重复计算公共子联盟将已经算出的联盟成本C(S)存入Map对象中作为全局缓存。后续计算Shapley值时公式中的v(S)和v(S∪{i})直接查缓存不再重复求解。联盟枚举的Matlab代码片段如下N 3; % 参与者数量 v zeros(1, 2^N - 1); % 特征函数存储编号对应二进制掩码 for mask 1:(2^N - 1) S find(bitget(mask, 1:N)); % 将掩码转换为成员列表 C_ind_sum 0; for i S C_ind_sum C_ind_sum C_ind(i); % 单主体独立最优成本 end C_coal solve_coalition(S); % 求解联盟联合调度成本 v(mask) C_ind_sum - C_coal; % 特征函数 成本节约量 end各位不用纠结这个solve_coalition到底怎么实现下一部分会展开。这里先把特征函数的结构说清楚v的维度是2^N - 1索引方式用mask的十进制编号这对后续Shapley值的实现特别舒服。3. 利益分配核心方法及其Matlab实现3.1 “分什么”和“怎么分”的基本逻辑很多刚接触这个课题的人会误以为合作博弈分配的是“总收益”其实在我这个框架里更准确的说法是“分成本节约”。各主体独立运行的总成本是C_ind_total Σ C_i^ind合作后的联合运行总成本是C(N)。由于联合运行可以用燃气轮机代替高价的电网购电、可以用储能协调不同主体的峰谷互补通常C(N) C_ind_total中间的差值ΔC C_ind_total - C(N)就是可供分配的合作剩余。分配的目标是给每个主体i一个成本分摊结果x_i使得Σ x_i C(N) 预算平衡或者等价地Σ φ_i ΔC 把节约分光x_i ≤ C_i^ind 个体理性即合作后谁也不会比单干亏对任意子联盟SΣ_{i∈S} x_i ≤ C(S) 联盟理性维持联盟稳定性。如果用φ_i表示主体i分到的节约量则x_i C_i_ind - φ_i。下面所有分配代码都围绕φ_i进行。3.2 Shapley值的Matlab实现Shapley值的定义大家应该很熟φ_i(v) Σ_{S⊆N{i}} w(|S|) × [v(S∪{i}) - v(S)]其中w(|S|) |S|! × (n - |S| - 1)! / n!。这个公式的直观理解是主体i加入联盟S时带来的边际贡献在所有可能加入顺序中取平均。公平性很好但计算量是指数级的。好在三到五主体场景下几乎瞬间完成。我的实现方式是按mask循环避免递归更符合Matlab的矩阵思维function phi shapley_value(v, N) phi zeros(1, N); n N; for mask 1:(2^n - 1) S find(bitget(mask, 1:n)); s length(S); for i 1:n if ~ismember(i, S) % mask_with_i 是 S ∪ {i} 的掩码 mask_with_i mask 2^(i-1); v_S v(mask); v_Si v(mask_with_i); w factorial(s) * factorial(n - s - 1) / factorial(n); phi(i) phi(i) w * (v_Si - v_S); end end end end这段代码有两点值得说明。第一mask和2^(i-1)的加法运算在Matlab里非常快比用子集字符串、表格组织数据安全得多。第二v数组索引从1开始mask1对应的二进制数是001表示成员1mask2对应010表示成员2mask3对应011表示子联盟{1,2}以此类推。这也是上面特征函数存储方式能无缝衔接的原因。这个函数我在实际运行中测过N5时循环次数为31×5155次耗时几乎可以忽略。N7时是127×7889次依然很快。真正慢的是每个联盟的调度求解而不是Shapley值本身。3.3 核仁Nucleolus与核心Core约束Shapley值有唯一性、对称性、可加性等公理但在联盟理性上并不保证分配结果一定处于核心之内。也就是说可能存在某个子联盟S使得Σ_{i∈S} φ_i v(S)这样的分配对这个子联盟没有吸引力联盟稳定性存疑。所以我在工程实践中还会算一下核仁用线性规划的方式找一组分配让最大“不满意量”达到最小。核仁的求解思路如下。定义不满意量ε(S) Σ_{i∈S} φ_i - v(S)希望最小化最大的ε(S)。加上预算平衡约束Σ φ_i v(N)得到的线性规划问题是minimize ε subject to: Σ_{i∈S} φ_i - ε ≤ v(S)对所有S⊂N Σ_{i∈N} φ_i v(N) φ_i ≥ 0Matlab可以用Yalmip调用Gurobi或Cplex来解phi sdpvar(1, N); eps sdpvar(1, 1); cons [sum(phi) v(end)]; % 注意v(end)对应全联盟2^N-1 for mask 1:(2^N - 2) S find(bitget(mask, 1:N)); cons [cons, sum(phi(S)) - eps v(mask)]; end obj eps; optimize(cons, obj); phi_nucleolus value(phi);核仁的优势在于只要核心非空核仁给出的分配必定位于核心内部而且结果唯一。缺点是当联盟数量较大时约束数太多线性规划规模迅速膨胀。不过对三五个主体的小算例这个方案非常闭眼可用。我在实际项目中的做法是先算Shapley值然后检查它是否满足所有联盟理性约束如果满足就直接用Shapley值出报告因为它解释性好、计算简单如果不满足再切换到核仁。两种算法都实现一遍其实没多少代码量但对最终结果的鲁棒性提升非常明显。4. 优化调度模型与Yalmip求解关键细节4.1 目标函数购电成本、燃料成本与运行维护成本联盟联合调度的目标函数是整个联盟在调度周期内的总运行成本最小。我采用的成本项包括三类从上级电网购电成本Σ_t price_e(t) × P_buy(t)天然气购买成本Σ_t price_g × (F_gt(t) F_gb(t))设备运行维护成本Σ_t Σ_k cost_om(k) × P_k(t)其中F_gt、F_gb分别是燃气轮机和燃气锅炉的耗气量通过效率折算得到。具体为F_gt(t) P_gt(t) / (η_e × LHV)F_gb(t) H_gb(t) / (η_gb × LHV)。如果考虑碳排放还可以再加一项碳税成本但为了先跑通利益分配主流程我建议第一版模型先不加碳约束把设备模型、平衡约束、分配流程跑通之后再扩展。目标函数写成Yalmip形式cost sum(price_e * P_buy) price_g * sum(F_gt F_gb) ... sum(om_cost_gt * P_gt) sum(om_cost_gb * H_gb) ... sum(om_cost_ch * P_ch) sum(om_cost_dis * P_dis) ... sum(om_cost_eb * P_eb) sum(om_cost_pv * P_pv);注意Yalmip中的*对向量和标量混合时要小心先把所有时间序列变量定义成1×T的sdpvar向量再用总和汇总不容易出错。4.2 调度模型的Yalmip关键约束代码下面我给出联盟调度的核心代码框架。输入是联盟成员列表S以及各主体对应的参数结构体输出是该联盟的最优成本。function C_coal solve_coalition(S, params) T params.T; % 决策变量 P_buy sdpvar(1, T); P_gt sdpvar(1, T); H_gb sdpvar(1, T); P_eb sdpvar(1, T); P_ch sdpvar(1, T); P_dis sdpvar(1, T); SOC sdpvar(1, T); F_gt sdpvar(1, T); F_gb sdpvar(1, T); cons []; % 根据联盟成员计算集合内总负荷与总光伏出力 [L_e, L_h, P_pv] aggregate_loads(S, params); % 电平衡 cons [cons, P_buy P_pv P_gt P_dis L_e P_ch P_eb]; % 热平衡 cons [cons, params.eta_chp_h / params.eta_chp_e * P_gt H_gb params.eta_eb * P_eb L_h]; % 燃气轮机上下限与爬坡约束 cons [cons, P_gt params.P_gt_min, P_gt params.P_gt_max]; cons [cons, diff(P_gt) -params.ramp_gt, diff(P_gt) params.ramp_gt]; cons [cons, F_gt P_gt / (params.eta_gt * params.LHV)]; % 燃气锅炉约束 cons [cons, H_gb 0, H_gb params.H_gb_max, F_gb H_gb / (params.eta_gb * params.LHV)]; % 电锅炉约束 cons [cons, P_eb 0, P_eb params.P_eb_max]; % 储能约束 cons [cons, P_ch 0, P_ch params.P_ch_max, P_dis 0, P_dis params.P_dis_max]; cons [cons, SOC(2:T) SOC(1:T-1) params.eta_ch * P_ch(1:T-1) / params.E_cap ... - P_dis(1:T-1) / (params.eta_dis * params.E_cap)]; cons [cons, SOC params.SOC_min, SOC params.SOC_max, SOC(1) params.SOC_init, SOC(T) params.SOC_end]; % 购电上限 cons [cons, P_buy 0, P_buy params.P_line_max]; % 目标函数 cost sum(params.price_e .* P_buy) params.price_g * sum(F_gt F_gb) ... params.om_gt * sum(P_gt) params.om_gb * sum(H_gb) ... params.om_eb * sum(P_eb) params.om_ch * sum(P_ch) params.om_dis * sum(P_dis); ops sdpsettings(solver, gurobi, verbose, 0); optimize(cons, cost, ops); C_coal value(cost); end有几个细节要注意。燃气轮机的爬坡约束diff(P_gt)在T24时长度是23项要注意首末时段的对齐。储能最后时段的SOC_end约束不能省否则模型会把储能电量放到零来省成本结果看似省钱实际是透支电池。所有带效率参数的变量在计算耗气量时要严格注意LHV单位是kWh/m³还是MJ/m³单位不一致会导致结果偏差很大。4.3 求解器选择与求解速度优化这个模型本质上是一个混合整数线性规划问题——没错里面有连续变量本来没有0-1变量时是纯LP但由于电锅炉、储能充放电同时性、机组启停逻辑在实际中往往需要引入二进制变量就会变成MILP。在这种场景下我建议优先使用Gurobi和Cplex其次用CBC或SCIP兜底。Yalmip在solvesdp或optimize之前可以通过sdpsettings(solver,gurobi)指定求解器。如果嫌MILP慢第一版可以先不引入启停变量把燃气轮机出力下限设为一个正的小值储能充放电通过SOC自然避免同时充放——纯LP问题求解飞快。这样做出来的分配结论不会有质的差别但调试效率高一个量级。等流程全部跑通需要做更精细的分析时再升级成MILP。另外一个提速技巧是不同联盟的调度模型中大量参数和约束结构是共享的。可以先用常量参数一次性生成Yalmip模型模板然后循环里只修改负荷和光伏数据利用assign重置相关参数。这样比每个联盟都重新建模、重新构造sdpvar对象要快很多。小算例上感受不明显但对6个以上主体的场景建模开销甚至会超过求解开销这个优化非常值得做。5. 一个完整算例三主体24小时典型日5.1 算例基础参数这里我给出一个典型日算例的原始参数方便大家按同样数据复现出接近的结果。调度周期为24小时步长为1小时。天然气价格设为2.5元/m³天然气低热值LHV取9.7 kWh/m³。分时电价采用三段式峰谷平电价峰段10点到15点、18点到21点为1.2元/kWh谷段23点到次日7点为0.4元/kWh其余时段0.75元/kWh。主体A的配置光伏装机容量800 kW蓄电池容量1000 kWh充电效率0.95放电效率0.95最大充放电功率200 kW基础电负荷在300到600 kW之间波动。主体B的配置燃气轮机额定出力500 kW发电效率0.35热电比1.2出力下限100 kW爬坡速率100 kW/h燃气锅炉容量800 kW效率0.9同时带200到400 kW的热负荷。主体C的配置电负荷400到900 kW热负荷300到500 kW无自发电设备有独立的电锅炉500 kW和蓄热罐。光伏出力曲线我用了一个典型日归一化系数序列6点开始出力12点左右达到峰值0.8518点后归零。这里就不贴完整24个数字了大家手头有当地典型日数据可以直接替换。5.2 独立运行与联合运行成本对比三主体独立运行时各自求解本主体的最小成本。A主体的最优成本约为5200元其中光伏大发时段储能充电、晚间放电供应一部分负荷B主体成本约为6100元热电联产根据气价电价决定开停机盈亏平衡点附近可能有若干时段停机C主体成本约为7300元完全没有发电能力只能按市场价购电购热。三主体独立总成本C_ind_total为18600元。全联盟联合运行时总成本出现了明显下降。在这个算例里联合运行成本C(N)约为16200元成本节约ΔC为2400元约占独立总成本的13%。这个节约主要来自两处一是燃气轮机在电价高峰时段多发替换掉C主体从电网高价购电的功率二是蓄热电锅炉在谷电时段蓄热、峰电时段放热替代了部分燃气锅炉的天然气消耗。5.3 Shapley值分配结果示例用上一部分的Shapley函数计算后得到三个主体的节约分配为主体A分到约760元主体B分到约900元主体C分到约740元。三者加总正好等于2400元预算平衡成立。各主体的实际成本分摊为主体独立运行成本分得节约量实际分摊成本A5200元760元4440元B6100元900元5200元C7300元740元6560元从这个结果能明显看出B主体作为联盟中的容量提供者获得了最高的分配比例因为它承担了联合调度中最关键的出力调整角色A主体通过储能灵活性获得了可观收益C主体虽然自身没有太多灵活性但通过购买联盟内部便宜能源享受了成本下降因此也有参与合作的意愿。这种结果在面向管理层的汇报里很容易讲清楚没有人因为合作而吃亏而且分摊后的每项成本都低于自己单干时的成本这是合作博弈能落地的核心前提。5.4 核仁校验与结果对照我再把核仁求解跑一遍得到的分配结果是A约745元、B约915元、C约740元整体分配格局与Shapley值非常接近核仁结果也落入了核心内部。差异主要出现在A和B之间A在Shapley值中的贡献被更平均地分摊给所有联盟而核仁更关注最不满意子联盟的边界会削弱“边际贡献最大者”的优势。两者结果接近也说明这个算例的博弈结构相对均衡不依赖某一个人形成垄断性联盟。在向客户交付时我通常把Shapley值和核仁的分配结果同时附上让各方自己选择分配规则。实践中经常有这种情况某个主体觉得Shapley方法中自己贡献计算复杂不透明但核仁方法只需要确认自己没有吃亏就行另一个主体则认为Shapley值的公理性更强。规则本身总要选一个但两个都摆出来谈判空间就大了很多。6. 常见问题与排查技巧实录6.1 特征函数不满足超可加性怎么办合作博弈的一个基础假设是联盟创造的价值大于成员单独创造的价值之和也就是特征函数满足超可加性对任意两个不相交的联盟S和Tv(S∪T) ≥ v(S) v(T)。但在实际算例里我遇到过联合优化成本反而比独立运行更高的情况尤其是当联盟内部没有互补设备、反而增加了系统复杂度时。原因通常是设备运行约束过强比如燃气轮机最小出力限制高、储能强制末时段SOC回满、爬坡约束过紧导致联合调度里灵活性被锁死。排查方式先拿单个联盟的调度结果和成员独立成本之和做对比看是哪一个约束导致了成本上升然后判断该约束是否合理。如果在现实中确实无法放松那就说明这些主体本身就不适合结盟硬算出来的负特征函数会让Shapley值出现负分配工程上不可接受。6.2 联盟枚举太多运行时间爆炸解决思路有三层。第一层先确认是否真的需要枚举所有联盟如果参与主体数为5以下不用优化直接枚举即可。第二层把联盟数量和最终利益分配的关系做敏感性分析很多情况下某些主体是否合作对组分配结果影响很小可以用剪枝思想剔除边际贡献过小的联盟。第三层用抽样Shapley值也就是随机抽取联盟子集估计平均边际贡献。我在实际项目里遇到过8个参与者的场景全枚举255个联盟每个联盟一个MILP半天跑不完。后来改用随机抽样的Shapley值估算配合并行计算工具箱在误差可控的范围内把时间压缩到半小时以内。虽然理论上的精确Shapley值被替代成了估计值但对实际利益分配谈判来说精度足够。6.3 Yalmip模型提示Infeasible模型不可行的排查路径建议按顺序来先查平衡约束。最常见的问题是联盟内所有主体的总负荷、总光伏出力、总设备容量没有正确聚合并集比如C主体的热负荷写成了kW而B主体设备容量用的是kWh单位不一致导致热平衡约束量纲错误。再查储能约束。SOC递推公式里首末时段索引容易错位导致矩阵维度不匹配或初值约束与递推方程冲突。最后才是检查求解器选项比如Gurobi的数值容差设置过于严格。一个小技巧在求解前把模型里所有变量的上下限列表打出来逐项核对容量数据是否超出合理范围。很多时候模型不可行是因为某个负荷数据填错了一个数量级找了个零成本“虚拟机组”硬平衡了结果其他主体成本异常Yalmip直接报不可行。这种数据问题单靠看模型代码很难发现必须回到输入数据表里去查。6.4 分配结果不满足核心或出现负分配如果Shapley值结果存在某个子联盟权重异常多半是该子联盟的特征函数计算发生了错误比如求解器没有真正收敛到最优解或者某个联盟的调度模型里约束缺失导致成本被低估。此时v(S)偏大分配给该联盟成员的收益被高估。解决方法是把每个联盟的成本数据单独核实先确认特征函数的非负性和联盟理性成立再去做分配。另一种情况是合作剩余本身太小同时某个主体对联盟贡献的边际增量在特殊价格下偏低出现负的φ_i。我在一个电价差异很大的算例里碰到过类似情况某主体的设备参数和联盟目标完全不一致加入联盟反而要为了配合其他主体而牺牲自身最优运行点Shapley给它的分配接近于零。此时的工程处理是先检查是否出现了负边际贡献如果有说明该主体的设备配置和联盟整体方向不匹配要么优化其设备参数要么把它排除在联盟之外。6.5 常见问题速查表问题现象最可能的原因排查与解决建议联盟成本反而高于独立成本之和约束过严或联盟互补性不足检查燃气轮机最小出力、SOC末值约束重新评估结盟组合Shapley值出现负值某个特征函数v(S∪{i})-v(S)为负是设备不匹配问题建议调整主体设备参数Yalmip报Infeasible平衡约束或SOC约束索引错位逐个打印变量的上下限核对单位与维度求解时间过长联盟数量多且每个联盟是MILP引入二进制变量前先跑LP版本大规模时用随机Shapley估计分配结果对某个主体不公平单一分配规则不满足主观期待Shapley、核仁、核心约束一起算把多规则结果交给谈判最后再分享一个我个人在实际操作中的体会。做合作博弈的利益分配最容易忽略的其实不是博弈算法本身而是底层调度模型的可信度。v(N)再完美如果每个联盟的联合优化模型里藏着没写全的约束算出来的分配结果就没有任何落地的说服力。所以我建议在写这个课题的Matlab代码时一定把“单主体独立优化”当作一个基准测试先跑通然后在这个基础上逐步增加联盟成员。每次增加一个主体都手动验证一次联合成本是否低于成员独立成本之和。这条路径虽然笨但它能让你在写复杂特征函数和Shapley函数之前就对调度模型的正确性有足够信心。如果后续想扩展这个项目可以从三个方向入手一是把单一典型日扩展到多场景随机规划用场景概率加权构造期望成本二是加入碳交易成本让利益分配同时体现环境效益三是把当前集中式联合调度改成分布式的ADMM求解让各主体不交换内部数据也能完成协同。这些扩展的方向都建立在今天这个框架之上核心的建模、特征函数、分配算法部分完全不用动。 合作博弈驱动的综合能源系统利益分配优化调度：Matlab建模与求解实践
最近在梳理综合能源系统方向的一个项目主题是“基于合作博弈的综合能源系统利益分配优化调度”。做多能互补或者园区能量管理的朋友应该都有体会单主体的优化调度很好写但一旦涉及多个主体共享同一个配电网、同一级天然气管网比如光储微网、燃气热电联产站、售电公司聚在一块儿各自只管自己的账整个系统的运行成本往往压不下来设备利用率也上不去。这个课题的核心思路就是用合作博弈来协调多个主体之间的关系把“各自为政”变成“先合作、再分利”最终落到一个可复现的Matlab求解框架上。这套方法论能解决的问题很具体多主体联合调度时合作带来的总成本节约怎么计算这个节约又该如何在各参与主体之间公平分配才能保证没有人觉得自己吃亏、联盟关系稳定可持续。适合两类人来参考一类是做园区综合能源、多微网协调控制方向的研究生和工程师另一类是刚接触博弈论在能源系统中应用、想把理论模型落地成代码的初学者。下面的内容都是我在实际算例上跑过的思路和代码按从建模到求解再到分配的顺序来讲。1. 项目整体定位与技术思路拆解1.1 这个课题到底在解决什么问题区域综合能源系统的典型场景是一个园区里有几个相对独立的能源主体比如A主体拥有光伏和蓄电池白天光伏出力充足但本地负荷小B主体是一套燃气轮机热电联产装置能同时产电产热但发电边际成本较高C主体是负荷聚合商手里有电负荷和热负荷对能量有刚性需求。如果把三者分开独立调度A会为了减少弃光在自己的小系统里强行给电池充满B会按自己的气价和电价判断开不开机C只能按实时市场价购电购热。问题就出来了晚间光伏归零A要高价从电网买电实际上B的热电联产如果开机可以更便宜地供电白天光伏大发的时段B可能因为电价低不愿意开机发电但C恰恰需要热导致B要用低效率的燃气锅炉单独供热储能设备的充放电策略只基于本地电价判断没有考虑相邻主体的需求曲线导致整个系统的峰谷调节能力被浪费。这些问题本质上是“信息割裂”和“利益冲突”造成的局部最优。合作博弈的思路则是先计算各个主体及任意子集形成联盟后的联合运行最优成本再把“单独运行成本之和”和“联盟运行成本”之间的差值——也就是合作带来的成本节约——通过某种规则分给每个参与者。1.2 为什么选合作博弈不选主从博弈或集中式优化很多人在做多主体能量管理时会先想到两类方法一类是集中式优化把所有设备、负荷、网络约束放进一个大模型里求全局最优另一类是主从博弈比如园区运营商做领导者用户或微网做跟随者。集中式优化解决不了隐私和公平性问题——各主体的内部数据、设备参数都要上交而且分配结果没法服众。主从博弈更适合层级不对等的场景比如售电商和终端用户。当参与者地位基本平等、大家都要出钱出力共同组成一个虚拟聚合体时合作博弈是最自然的选择。合作博弈和另外两者的核心差异在于它把问题拆成了“怎么联合运行更优”和“怎么分配收益更公平”两个先后阶段。这种先后逻辑符合园区实际落地流程先签合作协议、确立联合调度的运行方式再定期算账、按贡献分配收益。相比Nash非合作博弈中寻找稳定均衡点的做法合作博弈中的Shapley值、核仁等方法有明确的分配公式和线性规划求解路径代码实现难度低得多结果也容易向非专业背景的管理层解释。1.3 总体框架两层模型四个步骤我在这套代码里采用的框架不复杂分四步第一步定义参与者集合N{1,2,...,n}初始化各主体的设备参数、负荷曲线、风光出力序列、分时电价、气价数据第二步枚举所有非空子联盟S⊆N对每个联盟分别求解一个多主体联合经济调度优化模型得到联盟最优运行成本C(S)再据此计算其特征函数v(S)第三步基于特征函数v(S)用矩阵方式组织联盟值计算Shapley值或用线性规划求解核仁Nucleolus、判断核心Core非空性第四步用分配结果倒推各主体实际应分摊的成本结合各主体的独立运行成本计算净收益并校验个体理性、联盟理性、集体理性三个条件。值得强调特征函数不是随便定义的。常见做法是令特征函数等于单独成本减去联盟成本体现联盟创造的价值。我在后面的代码里用的就是这种定义即v(S) Σ_{i∈S} C_i^ind - C^coal(S)其中C_i^ind是主体i独立运行时的最优成本C^coal(S)是联盟S联合调度的最优成本。这样v(S)为正值时代表联盟创造的成本节约数值越大说明这个联盟协同带来的收益越明显。2. 多主体综合能源系统建模与特征函数计算2.1 典型三主体的设备配置与关键约束项目算例我设计为三主体已经能完整展现合作博弈的分配逻辑又不会因为子联盟数量太多导致枚举复杂度过高主体A光储微网含光伏机组、蓄电池储能、基础电负荷主体B热电联产站含燃气轮机、余热回收装置、燃气锅炉可同时输出电力和热力主体C综合负荷聚合商含电负荷、热负荷自身无发电设备只能通过电网购电和购买天然气。每个主体的关键设备模型如下光伏出力PPV,t PPV_rated × ηPV,t其中ηPV,t为归一化出力系数由当地光照强度决定短期调度中可当作已知参数。蓄电池储能模型SOC_t SOC_{t-1} ηch × Pch,t × Δt / E_cap - Pdis,t × Δt / (ηdis × E_cap)约束条件0 ≤ SOC_t ≤ 1Pch,t、Pdis,t不同时为正且有上下限限制。燃气轮机CHP机组模型P_gt,min ≤ P_gt,t ≤ P_gt,maxH_gt,t η_h/η_e × P_gt,t 热电比固定其中η_e为发电效率η_h为余热回收效率。如果机组是可调热电比的抽凝机组热电比会在一个区间内变化模型变为混合整数规划这里先按背压式固定热电比处理。燃气锅炉模型0 ≤ H_gb,t ≤ H_gb,maxH_gb,t η_gb × F_gb,t × LHV其中F_gb,t为耗气量m³LHV为天然气低热值。电锅炉模型0 ≤ P_eb,t ≤ P_eb_maxH_eb,t η_eb × P_eb,t电锅炉是电转热的枢纽也是B和C之间热力耦合的关键。2.2 电力与热力平衡约束每个联盟内部的联合调度核心是满足整个联盟的总电平衡和总热平衡约束电平衡P_buy,t Σ P_pv,t Σ P_gt,t P_dis,t Σ L_e,t P_ch,t P_eb,t热平衡Σ H_gt,t Σ H_gb,t Σ H_eb,t Σ L_h,t H_chs,t其中P_buy,t为从上级电网购电功率L_e,t为电负荷L_h,t为热负荷。这里有个容易踩坑的点多主体联盟和单主体孤岛运行时平衡方程的形式完全一致只是把各主体的变量聚合到一起。这也是为什么我建议把调度模型的求解写成函数形式输入是“参与者集合S内部的设备参数、负荷、价格”输出是该联盟的最优成本。这样枚举所有子联盟时就不用重复建模直接改参数调函数即可。2.3 联盟枚举与特征函数缓存策略对于n个参与者非空子联盟总数是2^n - 1。n3时只有7个联盟n4时是15个n5时是31个直接枚举完全没有压力。但如果参与主体数量到了7、8个以上联盟数量就变成127、255个每个联盟都要单独求解一个MILP模型计算时间会迅速膨胀。我在代码里做了两个优化用二进制掩码0-1组合表示子集避免重复计算公共子联盟将已经算出的联盟成本C(S)存入Map对象中作为全局缓存。后续计算Shapley值时公式中的v(S)和v(S∪{i})直接查缓存不再重复求解。联盟枚举的Matlab代码片段如下N 3; % 参与者数量 v zeros(1, 2^N - 1); % 特征函数存储编号对应二进制掩码 for mask 1:(2^N - 1) S find(bitget(mask, 1:N)); % 将掩码转换为成员列表 C_ind_sum 0; for i S C_ind_sum C_ind_sum C_ind(i); % 单主体独立最优成本 end C_coal solve_coalition(S); % 求解联盟联合调度成本 v(mask) C_ind_sum - C_coal; % 特征函数 成本节约量 end各位不用纠结这个solve_coalition到底怎么实现下一部分会展开。这里先把特征函数的结构说清楚v的维度是2^N - 1索引方式用mask的十进制编号这对后续Shapley值的实现特别舒服。3. 利益分配核心方法及其Matlab实现3.1 “分什么”和“怎么分”的基本逻辑很多刚接触这个课题的人会误以为合作博弈分配的是“总收益”其实在我这个框架里更准确的说法是“分成本节约”。各主体独立运行的总成本是C_ind_total Σ C_i^ind合作后的联合运行总成本是C(N)。由于联合运行可以用燃气轮机代替高价的电网购电、可以用储能协调不同主体的峰谷互补通常C(N) C_ind_total中间的差值ΔC C_ind_total - C(N)就是可供分配的合作剩余。分配的目标是给每个主体i一个成本分摊结果x_i使得Σ x_i C(N) 预算平衡或者等价地Σ φ_i ΔC 把节约分光x_i ≤ C_i^ind 个体理性即合作后谁也不会比单干亏对任意子联盟SΣ_{i∈S} x_i ≤ C(S) 联盟理性维持联盟稳定性。如果用φ_i表示主体i分到的节约量则x_i C_i_ind - φ_i。下面所有分配代码都围绕φ_i进行。3.2 Shapley值的Matlab实现Shapley值的定义大家应该很熟φ_i(v) Σ_{S⊆N{i}} w(|S|) × [v(S∪{i}) - v(S)]其中w(|S|) |S|! × (n - |S| - 1)! / n!。这个公式的直观理解是主体i加入联盟S时带来的边际贡献在所有可能加入顺序中取平均。公平性很好但计算量是指数级的。好在三到五主体场景下几乎瞬间完成。我的实现方式是按mask循环避免递归更符合Matlab的矩阵思维function phi shapley_value(v, N) phi zeros(1, N); n N; for mask 1:(2^n - 1) S find(bitget(mask, 1:n)); s length(S); for i 1:n if ~ismember(i, S) % mask_with_i 是 S ∪ {i} 的掩码 mask_with_i mask 2^(i-1); v_S v(mask); v_Si v(mask_with_i); w factorial(s) * factorial(n - s - 1) / factorial(n); phi(i) phi(i) w * (v_Si - v_S); end end end end这段代码有两点值得说明。第一mask和2^(i-1)的加法运算在Matlab里非常快比用子集字符串、表格组织数据安全得多。第二v数组索引从1开始mask1对应的二进制数是001表示成员1mask2对应010表示成员2mask3对应011表示子联盟{1,2}以此类推。这也是上面特征函数存储方式能无缝衔接的原因。这个函数我在实际运行中测过N5时循环次数为31×5155次耗时几乎可以忽略。N7时是127×7889次依然很快。真正慢的是每个联盟的调度求解而不是Shapley值本身。3.3 核仁Nucleolus与核心Core约束Shapley值有唯一性、对称性、可加性等公理但在联盟理性上并不保证分配结果一定处于核心之内。也就是说可能存在某个子联盟S使得Σ_{i∈S} φ_i v(S)这样的分配对这个子联盟没有吸引力联盟稳定性存疑。所以我在工程实践中还会算一下核仁用线性规划的方式找一组分配让最大“不满意量”达到最小。核仁的求解思路如下。定义不满意量ε(S) Σ_{i∈S} φ_i - v(S)希望最小化最大的ε(S)。加上预算平衡约束Σ φ_i v(N)得到的线性规划问题是minimize ε subject to: Σ_{i∈S} φ_i - ε ≤ v(S)对所有S⊂N Σ_{i∈N} φ_i v(N) φ_i ≥ 0Matlab可以用Yalmip调用Gurobi或Cplex来解phi sdpvar(1, N); eps sdpvar(1, 1); cons [sum(phi) v(end)]; % 注意v(end)对应全联盟2^N-1 for mask 1:(2^N - 2) S find(bitget(mask, 1:N)); cons [cons, sum(phi(S)) - eps v(mask)]; end obj eps; optimize(cons, obj); phi_nucleolus value(phi);核仁的优势在于只要核心非空核仁给出的分配必定位于核心内部而且结果唯一。缺点是当联盟数量较大时约束数太多线性规划规模迅速膨胀。不过对三五个主体的小算例这个方案非常闭眼可用。我在实际项目中的做法是先算Shapley值然后检查它是否满足所有联盟理性约束如果满足就直接用Shapley值出报告因为它解释性好、计算简单如果不满足再切换到核仁。两种算法都实现一遍其实没多少代码量但对最终结果的鲁棒性提升非常明显。4. 优化调度模型与Yalmip求解关键细节4.1 目标函数购电成本、燃料成本与运行维护成本联盟联合调度的目标函数是整个联盟在调度周期内的总运行成本最小。我采用的成本项包括三类从上级电网购电成本Σ_t price_e(t) × P_buy(t)天然气购买成本Σ_t price_g × (F_gt(t) F_gb(t))设备运行维护成本Σ_t Σ_k cost_om(k) × P_k(t)其中F_gt、F_gb分别是燃气轮机和燃气锅炉的耗气量通过效率折算得到。具体为F_gt(t) P_gt(t) / (η_e × LHV)F_gb(t) H_gb(t) / (η_gb × LHV)。如果考虑碳排放还可以再加一项碳税成本但为了先跑通利益分配主流程我建议第一版模型先不加碳约束把设备模型、平衡约束、分配流程跑通之后再扩展。目标函数写成Yalmip形式cost sum(price_e * P_buy) price_g * sum(F_gt F_gb) ... sum(om_cost_gt * P_gt) sum(om_cost_gb * H_gb) ... sum(om_cost_ch * P_ch) sum(om_cost_dis * P_dis) ... sum(om_cost_eb * P_eb) sum(om_cost_pv * P_pv);注意Yalmip中的*对向量和标量混合时要小心先把所有时间序列变量定义成1×T的sdpvar向量再用总和汇总不容易出错。4.2 调度模型的Yalmip关键约束代码下面我给出联盟调度的核心代码框架。输入是联盟成员列表S以及各主体对应的参数结构体输出是该联盟的最优成本。function C_coal solve_coalition(S, params) T params.T; % 决策变量 P_buy sdpvar(1, T); P_gt sdpvar(1, T); H_gb sdpvar(1, T); P_eb sdpvar(1, T); P_ch sdpvar(1, T); P_dis sdpvar(1, T); SOC sdpvar(1, T); F_gt sdpvar(1, T); F_gb sdpvar(1, T); cons []; % 根据联盟成员计算集合内总负荷与总光伏出力 [L_e, L_h, P_pv] aggregate_loads(S, params); % 电平衡 cons [cons, P_buy P_pv P_gt P_dis L_e P_ch P_eb]; % 热平衡 cons [cons, params.eta_chp_h / params.eta_chp_e * P_gt H_gb params.eta_eb * P_eb L_h]; % 燃气轮机上下限与爬坡约束 cons [cons, P_gt params.P_gt_min, P_gt params.P_gt_max]; cons [cons, diff(P_gt) -params.ramp_gt, diff(P_gt) params.ramp_gt]; cons [cons, F_gt P_gt / (params.eta_gt * params.LHV)]; % 燃气锅炉约束 cons [cons, H_gb 0, H_gb params.H_gb_max, F_gb H_gb / (params.eta_gb * params.LHV)]; % 电锅炉约束 cons [cons, P_eb 0, P_eb params.P_eb_max]; % 储能约束 cons [cons, P_ch 0, P_ch params.P_ch_max, P_dis 0, P_dis params.P_dis_max]; cons [cons, SOC(2:T) SOC(1:T-1) params.eta_ch * P_ch(1:T-1) / params.E_cap ... - P_dis(1:T-1) / (params.eta_dis * params.E_cap)]; cons [cons, SOC params.SOC_min, SOC params.SOC_max, SOC(1) params.SOC_init, SOC(T) params.SOC_end]; % 购电上限 cons [cons, P_buy 0, P_buy params.P_line_max]; % 目标函数 cost sum(params.price_e .* P_buy) params.price_g * sum(F_gt F_gb) ... params.om_gt * sum(P_gt) params.om_gb * sum(H_gb) ... params.om_eb * sum(P_eb) params.om_ch * sum(P_ch) params.om_dis * sum(P_dis); ops sdpsettings(solver, gurobi, verbose, 0); optimize(cons, cost, ops); C_coal value(cost); end有几个细节要注意。燃气轮机的爬坡约束diff(P_gt)在T24时长度是23项要注意首末时段的对齐。储能最后时段的SOC_end约束不能省否则模型会把储能电量放到零来省成本结果看似省钱实际是透支电池。所有带效率参数的变量在计算耗气量时要严格注意LHV单位是kWh/m³还是MJ/m³单位不一致会导致结果偏差很大。4.3 求解器选择与求解速度优化这个模型本质上是一个混合整数线性规划问题——没错里面有连续变量本来没有0-1变量时是纯LP但由于电锅炉、储能充放电同时性、机组启停逻辑在实际中往往需要引入二进制变量就会变成MILP。在这种场景下我建议优先使用Gurobi和Cplex其次用CBC或SCIP兜底。Yalmip在solvesdp或optimize之前可以通过sdpsettings(solver,gurobi)指定求解器。如果嫌MILP慢第一版可以先不引入启停变量把燃气轮机出力下限设为一个正的小值储能充放电通过SOC自然避免同时充放——纯LP问题求解飞快。这样做出来的分配结论不会有质的差别但调试效率高一个量级。等流程全部跑通需要做更精细的分析时再升级成MILP。另外一个提速技巧是不同联盟的调度模型中大量参数和约束结构是共享的。可以先用常量参数一次性生成Yalmip模型模板然后循环里只修改负荷和光伏数据利用assign重置相关参数。这样比每个联盟都重新建模、重新构造sdpvar对象要快很多。小算例上感受不明显但对6个以上主体的场景建模开销甚至会超过求解开销这个优化非常值得做。5. 一个完整算例三主体24小时典型日5.1 算例基础参数这里我给出一个典型日算例的原始参数方便大家按同样数据复现出接近的结果。调度周期为24小时步长为1小时。天然气价格设为2.5元/m³天然气低热值LHV取9.7 kWh/m³。分时电价采用三段式峰谷平电价峰段10点到15点、18点到21点为1.2元/kWh谷段23点到次日7点为0.4元/kWh其余时段0.75元/kWh。主体A的配置光伏装机容量800 kW蓄电池容量1000 kWh充电效率0.95放电效率0.95最大充放电功率200 kW基础电负荷在300到600 kW之间波动。主体B的配置燃气轮机额定出力500 kW发电效率0.35热电比1.2出力下限100 kW爬坡速率100 kW/h燃气锅炉容量800 kW效率0.9同时带200到400 kW的热负荷。主体C的配置电负荷400到900 kW热负荷300到500 kW无自发电设备有独立的电锅炉500 kW和蓄热罐。光伏出力曲线我用了一个典型日归一化系数序列6点开始出力12点左右达到峰值0.8518点后归零。这里就不贴完整24个数字了大家手头有当地典型日数据可以直接替换。5.2 独立运行与联合运行成本对比三主体独立运行时各自求解本主体的最小成本。A主体的最优成本约为5200元其中光伏大发时段储能充电、晚间放电供应一部分负荷B主体成本约为6100元热电联产根据气价电价决定开停机盈亏平衡点附近可能有若干时段停机C主体成本约为7300元完全没有发电能力只能按市场价购电购热。三主体独立总成本C_ind_total为18600元。全联盟联合运行时总成本出现了明显下降。在这个算例里联合运行成本C(N)约为16200元成本节约ΔC为2400元约占独立总成本的13%。这个节约主要来自两处一是燃气轮机在电价高峰时段多发替换掉C主体从电网高价购电的功率二是蓄热电锅炉在谷电时段蓄热、峰电时段放热替代了部分燃气锅炉的天然气消耗。5.3 Shapley值分配结果示例用上一部分的Shapley函数计算后得到三个主体的节约分配为主体A分到约760元主体B分到约900元主体C分到约740元。三者加总正好等于2400元预算平衡成立。各主体的实际成本分摊为主体独立运行成本分得节约量实际分摊成本A5200元760元4440元B6100元900元5200元C7300元740元6560元从这个结果能明显看出B主体作为联盟中的容量提供者获得了最高的分配比例因为它承担了联合调度中最关键的出力调整角色A主体通过储能灵活性获得了可观收益C主体虽然自身没有太多灵活性但通过购买联盟内部便宜能源享受了成本下降因此也有参与合作的意愿。这种结果在面向管理层的汇报里很容易讲清楚没有人因为合作而吃亏而且分摊后的每项成本都低于自己单干时的成本这是合作博弈能落地的核心前提。5.4 核仁校验与结果对照我再把核仁求解跑一遍得到的分配结果是A约745元、B约915元、C约740元整体分配格局与Shapley值非常接近核仁结果也落入了核心内部。差异主要出现在A和B之间A在Shapley值中的贡献被更平均地分摊给所有联盟而核仁更关注最不满意子联盟的边界会削弱“边际贡献最大者”的优势。两者结果接近也说明这个算例的博弈结构相对均衡不依赖某一个人形成垄断性联盟。在向客户交付时我通常把Shapley值和核仁的分配结果同时附上让各方自己选择分配规则。实践中经常有这种情况某个主体觉得Shapley方法中自己贡献计算复杂不透明但核仁方法只需要确认自己没有吃亏就行另一个主体则认为Shapley值的公理性更强。规则本身总要选一个但两个都摆出来谈判空间就大了很多。6. 常见问题与排查技巧实录6.1 特征函数不满足超可加性怎么办合作博弈的一个基础假设是联盟创造的价值大于成员单独创造的价值之和也就是特征函数满足超可加性对任意两个不相交的联盟S和Tv(S∪T) ≥ v(S) v(T)。但在实际算例里我遇到过联合优化成本反而比独立运行更高的情况尤其是当联盟内部没有互补设备、反而增加了系统复杂度时。原因通常是设备运行约束过强比如燃气轮机最小出力限制高、储能强制末时段SOC回满、爬坡约束过紧导致联合调度里灵活性被锁死。排查方式先拿单个联盟的调度结果和成员独立成本之和做对比看是哪一个约束导致了成本上升然后判断该约束是否合理。如果在现实中确实无法放松那就说明这些主体本身就不适合结盟硬算出来的负特征函数会让Shapley值出现负分配工程上不可接受。6.2 联盟枚举太多运行时间爆炸解决思路有三层。第一层先确认是否真的需要枚举所有联盟如果参与主体数为5以下不用优化直接枚举即可。第二层把联盟数量和最终利益分配的关系做敏感性分析很多情况下某些主体是否合作对组分配结果影响很小可以用剪枝思想剔除边际贡献过小的联盟。第三层用抽样Shapley值也就是随机抽取联盟子集估计平均边际贡献。我在实际项目里遇到过8个参与者的场景全枚举255个联盟每个联盟一个MILP半天跑不完。后来改用随机抽样的Shapley值估算配合并行计算工具箱在误差可控的范围内把时间压缩到半小时以内。虽然理论上的精确Shapley值被替代成了估计值但对实际利益分配谈判来说精度足够。6.3 Yalmip模型提示Infeasible模型不可行的排查路径建议按顺序来先查平衡约束。最常见的问题是联盟内所有主体的总负荷、总光伏出力、总设备容量没有正确聚合并集比如C主体的热负荷写成了kW而B主体设备容量用的是kWh单位不一致导致热平衡约束量纲错误。再查储能约束。SOC递推公式里首末时段索引容易错位导致矩阵维度不匹配或初值约束与递推方程冲突。最后才是检查求解器选项比如Gurobi的数值容差设置过于严格。一个小技巧在求解前把模型里所有变量的上下限列表打出来逐项核对容量数据是否超出合理范围。很多时候模型不可行是因为某个负荷数据填错了一个数量级找了个零成本“虚拟机组”硬平衡了结果其他主体成本异常Yalmip直接报不可行。这种数据问题单靠看模型代码很难发现必须回到输入数据表里去查。6.4 分配结果不满足核心或出现负分配如果Shapley值结果存在某个子联盟权重异常多半是该子联盟的特征函数计算发生了错误比如求解器没有真正收敛到最优解或者某个联盟的调度模型里约束缺失导致成本被低估。此时v(S)偏大分配给该联盟成员的收益被高估。解决方法是把每个联盟的成本数据单独核实先确认特征函数的非负性和联盟理性成立再去做分配。另一种情况是合作剩余本身太小同时某个主体对联盟贡献的边际增量在特殊价格下偏低出现负的φ_i。我在一个电价差异很大的算例里碰到过类似情况某主体的设备参数和联盟目标完全不一致加入联盟反而要为了配合其他主体而牺牲自身最优运行点Shapley给它的分配接近于零。此时的工程处理是先检查是否出现了负边际贡献如果有说明该主体的设备配置和联盟整体方向不匹配要么优化其设备参数要么把它排除在联盟之外。6.5 常见问题速查表问题现象最可能的原因排查与解决建议联盟成本反而高于独立成本之和约束过严或联盟互补性不足检查燃气轮机最小出力、SOC末值约束重新评估结盟组合Shapley值出现负值某个特征函数v(S∪{i})-v(S)为负是设备不匹配问题建议调整主体设备参数Yalmip报Infeasible平衡约束或SOC约束索引错位逐个打印变量的上下限核对单位与维度求解时间过长联盟数量多且每个联盟是MILP引入二进制变量前先跑LP版本大规模时用随机Shapley估计分配结果对某个主体不公平单一分配规则不满足主观期待Shapley、核仁、核心约束一起算把多规则结果交给谈判最后再分享一个我个人在实际操作中的体会。做合作博弈的利益分配最容易忽略的其实不是博弈算法本身而是底层调度模型的可信度。v(N)再完美如果每个联盟的联合优化模型里藏着没写全的约束算出来的分配结果就没有任何落地的说服力。所以我建议在写这个课题的Matlab代码时一定把“单主体独立优化”当作一个基准测试先跑通然后在这个基础上逐步增加联盟成员。每次增加一个主体都手动验证一次联合成本是否低于成员独立成本之和。这条路径虽然笨但它能让你在写复杂特征函数和Shapley函数之前就对调度模型的正确性有足够信心。如果后续想扩展这个项目可以从三个方向入手一是把单一典型日扩展到多场景随机规划用场景概率加权构造期望成本二是加入碳交易成本让利益分配同时体现环境效益三是把当前集中式联合调度改成分布式的ADMM求解让各主体不交换内部数据也能完成协同。这些扩展的方向都建立在今天这个框架之上核心的建模、特征函数、分配算法部分完全不用动。