做过微电网优化调度的朋友应该都有体会项目做到一半最让你头疼的往往不是风光预测、不是储能SOC约束而是“电、热、气、储”这几个主体坐在一起谁都不愿意先让步。设备模型都能建优化目标都能写但最后总被一个问题卡住——联合运行确实省钱了省下的钱该怎么分分得不公平合作当场破裂。这也是我最初接触“基于合作博弈的综合能源系统利益分配优化调度”这个课题时的真实感受。这个方向的实质就是用合作博弈把多主体联合运行的总收益拆解成每个参与者的合理回报再配合优化调度模型在Matlab里用一段可复现的代码把“先合作、再分账”整个流程跑通。它适合正在做微电网、综合能源系统方向的研究生也适合想在企业多能互补项目里落地利益分配机制的工程师。今天这篇文章我就按自己实际做过的项目流程把建模思路、Shapley值计算、Matlab代码结构、求解器选型和排查技巧一次讲透。1. 项目整体设计与博弈论选型先想清楚“分蛋糕”再“做蛋糕”1.1 为什么综合能源系统一定要做“利益分配”先看一个常见场景一个园区里同时有光伏、风机、储能、燃气轮机和燃气锅炉还有一个负荷聚合商。单独运行时每家各干各的光伏白天大发但可能弃电燃气轮机按自己的热电比硬顶着供热储能只给自己调峰。联合运行之后燃气轮机可以先发电、余热供热光伏多出来的电拿去给电锅炉制热储能低谷充电高峰放电整体买电量明显下降系统总运行成本比各干各的更低。这本是个双赢的局面但问题恰恰出在“双赢”上——总收益提升了谁多分、谁少分没有一个大家都能接受的规则合作就进行不下去。这就是综合能源系统里的“利益分配”问题也是合作博弈发挥作用的地方。说白了合作博弈研究的就是一群人组成的联盟可以获得超额收益这个收益如何按照某种公平性准则分配给联盟成员使得每个成员分到的收益不低于单独行动时的收益从而维持联盟稳定。没有这一步你的优化调度结果再漂亮也落不了地。1.2 Shapley值、核仁、纳什议价三个分配方案怎么选利益分配不是只有一种算法。我在实际项目里一般从三个方案里选Shapley值、核仁、纳什议价。它们的思路完全不同选错了后面会非常难受。Shapley值是一个算边际贡献的经典方法每个成员分到的收益等于他加入所有可能联盟时带来的收益增量的平均值。它的优点是有明确的公理化基础对称性、有效性、可加性、哑元性质都满足大家普遍觉得“公平”也最容易写论文、最容易向业主解释。缺点是需要计算2的n次方个子联盟的收益主体超过15个计算量就开始吃紧。核仁的思想则更“硬核”它追求的是让最不满意的联盟的“不满意程度”降到最小保证任何子联盟都没有动力脱离大联盟。稳定性极好但求解要套线性规划且结果对模型误差比较敏感。纳什议价是另一个思路它融入了谈判系数可以体现不同主体的议价能力差异在实际商业谈判中更好用但参数怎么定是个主观活写进论文里容易被评审追问“你怎么标定的”。我的建议是如果目标是学术论证和可解释性优先选Shapley值如果项目方特别强调联盟绝对不能破裂、要追求博弈论意义上的“核”就用核仁做验证如果涉及真实交易谈判可以在Shapley值基础上引入权重系数演变成加权Shapley或Nash议价。下面这张表是我项目里常用的选型对照分配方案核心思想优点缺点适用场景Shapley值按边际贡献分配公理化公平可解释性强子联盟数量指数增长学术研究、多主体收益分摊核仁最小化最大不满意程度联盟稳定性最强求解复杂度高强调联盟稳固的工程场景纳什议价带谈判系数的讨价还价贴近商业谈判参数主观性强有明确谈判主体的商业项目1.3 为什么我还是用Matlab而不是Python这几年Python在优化领域声势很大但做微电网优化调度我依然推荐Matlab。这不是情怀是实际效率问题。电力系统领域的老牌工具箱和示例代码绝大多数都是Matlab写的Yalmip这个建模层的生态更是绕不开你定义一个sdpvar决策变量、写一组约束、丢给Gurobi求解整个流程比Python里拼PuLP或者自己手写约束矩阵要快得多。再一个Matlab的调试体验对科研场景极其友好变量工作区随时可以检查约束矩阵的维度和数值遇到“维度不匹配”这类问题一眼就能定位。那些说Python更高效的人多半是没被密密麻麻的索引报错折磨过。当然如果你的场景是深度学习和海量数据处理Python确实有优势。文章开头热门搜索词里的bilstm、transformer、ppo这些Matlab代码也说明现在在Matlab里做数据驱动方法已经不稀奇了。我的习惯是优化调度的核心用Matlab数据预处理和预测部分偶尔用Python两边通过文件接口打通。但为了让你能最快复现下面全文给的方案都基于Matlab Yalmip Gurobi的组合。2. 数学模型搭建从独立运行到联盟运行的三阶段框架2.1 三阶段框架把“合作收益”先算出来再分出去利益分配的前置条件是“有增量可以分”所以整个项目必须走一个三阶段框架这也是我对接项目时最常用的流程。第一阶段是独立运行优化。每个主体各自运行追求自身收益最大化或成本最小化算出每个人单独行动时的收益记为初始收益。第二阶段是联盟运行优化。把所有主体打包成一个整体统一调度以系统总成本最小为目标算出联盟总收益。第三阶段是利益分配。联盟收益减去独立收益之和得到的就是合作带来的增量收益用Shapley值或核仁把增量分给每个主体再叠加各主体独立收益得到最终分配结果。我踩过的最大的坑是很多人会跳过第一阶段直接做联盟优化然后分钱结果分完发现某个主体分到的钱还不如他单独干联盟立刻“破裂”。没有第一阶段的独立收益做基线你根本没法验证“个体理性”约束——也就是说每个主体参与联盟后分到的收益必须不低于他单干的收益。三阶段缺一不可顺序也不能乱。2.2 联盟运行阶段的目标函数与关键约束联盟运行阶段是核心它的目标函数我一般写成系统总运行成本最小化因为成本口径在工程里最直观也最容易和电网公司的结算表对上。总成本一般包含四块从上级电网购电费用、购气费用、新能源弃电惩罚成本以及储能充放电老化成本。[ \min \sum_{t1}^{T} \left( c_{e,t} P_{buy,t} c_{g,t} F_{g,t} c_{curtail} P_{curtail,t} c_{bat} (P_{ch,t} P_{dis,t}) \right) ]其中T是调度周期内的时段数一般取24小时。约束方面有几个必须写全的功率平衡约束、燃气轮机热电比约束、储能SOC递推约束、蓄水/蓄热罐的容量约束、联络线功率上下限约束。储能SOC递推是这里最容易漏项的地方[ SOC_{t1} SOC_t \eta_{ch} P_{ch,t} / E_{rate} - P_{dis,t} / (\eta_{dis} E_{rate}) ]这个式子表示储能状态在相邻时段之间的递推关系充放电效率要分开标定而且要防止同一个时段内储能又充又放那在物理上是严重违背实际的。处理办法是引入两个二进制变量表示充放电状态并加互斥约束。2.3 Shapley值计算子联盟收益怎么算才不出错Shapley值的计算是整个项目里最容易被误解的部分先说公式。设有一个包含n个主体的集合Nv(S)表示子联盟S的收益或者成本节约量那么主体i分到的Shapley值可以写为[ \phi_i(v) \sum_{S \subseteq N \setminus {i}} \frac{|S|!(n-|S|-1)!}{n!} \big( v(S \cup {i}) - v(S) \big) ]这个式子的含义是把主体i随机插入到某个子联盟S中观察它加入前后联盟收益的变化值v(S∪{i})−v(S)这就是i对这个联盟的边际贡献。把所有可能的S都算一遍再按排列组合数加权平均就是主体i的Shapley值。加权因子|S|!(n−|S|−1)!/n!的来历是在随机排列中S中的成员恰好排在i前面且N\S{i}中的成员恰好排在i后面的概率。实际操作里每个子联盟S都要重新运行一次联盟优化模型求出v(S)。这就是Shapley值计算量的主要来源。n3时只要算6次n8时已经是128次到了n12就是2048次每一个都要完整调用求解器。我在项目里发现一个非常容易错的地方有些子联盟里根本没有某个主体的关键设备比如S里只有光伏和储能那燃气轮机的燃气成本自然不该被算进去反过来如果S不包含负荷聚合商那S对外购电的功率边界也要跟着变。也就是说每个子联盟的优化模型参数都要随之调整绝对不是换个人名重复算同一个模型。3. Matlab代码实现与核心环节调试3.1 工程化代码结构从main函数开始拆分很多初学者喜欢把全部代码堆在一个脚本里跑通了还行一旦要改一个约束或者换一种分配方案整个人就疯了。我建议模仿工程项目的方式拆分文件一个典型结构是main.m % 主程序定义主体、循环子联盟、调用优化、分配、画图 case_data.m % 参数集中管理气价、电价、负荷曲线、设备参数 independent_opt.m % 第一阶段独立运行优化输出独立收益 coalition_opt.m % 第二阶段给定成员集合建立并求解联盟优化模型 subsets_gen.m % 生成所有子联盟并映射到成员集合 shapley_value.m % 第三阶段基于各子联盟收益计算Shapley值 verify_stability.m % 验证个体理性与联盟稳定性 plot_results.m % 可视化main函数不要一上来就写公式先搭一个清晰的执行流加载参数、枚举子联盟、for循环调coalition_opt、存下每个v(S)、计算Shapley值、写结果表格。这样后面加新主体、改模型都只需要改对应的局部文件。3.2 Yalmip建模核心片段与注释联盟优化阶段用Yalmip建模非常顺手这里我贴一段我项目里精简过的核心代码帮你建立直观印象。假设联盟里有燃气轮机、储能、光伏和负荷约束包含电功率平衡、SOC递推、购电上限和光伏出力上限function [cost, detail] coalition_opt(members, data) % members: 1xN逻辑向量标记哪些主体参与本子联盟 T data.T; % 24小时 yalmip(clear); % 决策变量 P_buy sdpvar(1, T); % 从上级电网购电功率 P_chp sdpvar(1, T); % 燃气轮机发电功率 P_pv sdpvar(1, T); % 光伏实际出力 P_ch sdpvar(1, T); % 储能充电功率 P_dis sdpvar(1, T); % 储能放电功率 SOC sdpvar(1, T1); % 储能SOC状态 u_ch binvar(1, T); % 充电状态标志 u_dis binvar(1, T); % 放电状态标志 % 约束条件 cons []; % 电功率平衡购电 光伏 燃气轮机放电 储能放电 负荷 充电 cons [cons, P_buy P_pv P_chp P_dis data.P_load P_ch]; % 光伏出力限制 cons [cons, 0 P_pv data.P_pv_max]; % 购电功率限制 cons [cons, 0 P_buy data.P_buy_max]; % 燃气轮机出力上下限 cons [cons, data.P_chp_min P_chp data.P_chp_max]; % 储能SOC递推 for t 1:T cons [cons, SOC(t1) SOC(t) data.eta_ch * P_ch(t)/data.E_rate ... - P_dis(t)/(data.eta_dis * data.E_rate)]; end % 充放电互斥与功率限制 cons [cons, 0 P_ch u_ch * data.P_bat_max]; cons [cons, 0 P_dis u_dis * data.P_bat_max]; cons [cons, u_ch u_dis 1]; cons [cons, data.SOC_min SOC data.SOC_max]; cons [cons, SOC(1) data.SOC_init, SOC(T1) data.SOC_init]; % 目标函数按上面讲的四部分成本累加 obj sum(data.c_e .* P_buy) sum(data.c_g .* (P_chp / data.eta_chp)) ... sum(data.c_curtail .* (data.P_pv_max - P_pv)) ... sum(data.c_bat .* (P_ch P_dis)); ops sdpsettings(solver, gurobi, verbose, 0, mipgaptol, 1e-4); optimize(cons, obj, ops); cost value(obj); detail.P_buy value(P_buy); detail.P_chp value(P_chp); % ... 其他结果 end这段代码有两点值得注意第一充放电互斥用二进制变量处理比把充放电合写成功率变量再分段建模要稳定得多第二SOC递推约束里充放电效率分别乘在对应项上如果项目允许储能同时充放效率模型会失真物理意义也不对。如果你没有Gurobi把solver改成cplex或quadprog也能跑但整数变量多的时候内置求解器会非常吃力。3.3 求解器配置与非线性处理的实战经验求解器配置这件事值得单独说。Yalmip只是一个建模语言真正求解的是Gurobi、Cplex这类后端求解器。我实测下来在同一组约束和变量规模下Gurobi比Matlab内置的linprog/intlinprog快3到5倍是很正常的遇到混合整数线性规划问题差距更大。安装好Gurobi之后第一次使用前记得在Matlab里执行yalmip(clear)然后重新执行optimize否则Yalmip可能还在用旧的求解器路径报“No suitable solver”的错。非线性处理是另一个高频翻车点。燃气轮机的热电比约束、储能效率、电网潮流里的电压乘积项很多地方天然是非线性的。我强烈建议把模型尽量线化热负荷平衡用线性等式设备效率作为常数代入功率平衡用线性等式。如果真的躲不开乘积项比如电压与功率相乘那就用大M法引入辅助二进制变量线性化。能建MILP就不建MINLP这会让你的求解时间从几个钟头骤降到几十秒而且Gurobi对MILP的求解器远成熟于MINLP。4. 仿真结果分析与分配有效性验证4.1 输出结果设计光有数字还不行要有对照仿真跑完之后最忌讳上来就给一张分配表。利益分配的结果要能说服人必须有一组对照关系独立运行总成本、联盟运行总成本、合作带来的成本节约额、每个主体的独立收益、每个主体的Shapley分配值、分配后各主体最终收益。这六项缺一不可因为它们共同回答三个问题合作到底省了多少每个主体分到了多少增量每个主体最终收益是否都不低于独立收益有个研究生跟我讨论项目时困惑“为什么我Shapley值算出来光伏电站分到的增量是负数”。我让他把子联盟收益表打出来一看才发现他在算v(S)时把系统固定运维成本算进了每个子联盟而这个固定成本不因联盟构成变化导致边际贡献被系统性扭曲。记住一句话v(S)里只应该包含“加入这个联盟才产生的可变动收益/成本”固定成本不应参与边际贡献计算。4.2 一个三主体案例的收益分配结果示例为了让分配过程更具体我构造一个简化案例。假设系统里只有三个主体A是光伏电站B是储能运营商C是燃气轮机综合供能商。单位统一为万元/日各主体独立运行收益、联盟运行总收益、增量收益和两种分配方案的结果如下表主体独立运行收益Shapley增量分配分配后总收益核仁增量分配分配后总收益A 光伏2.00.92.90.82.8B 储能1.20.71.91.02.2C 燃气轮机2.51.43.91.23.7合计5.73.08.73.08.7这个例子联盟总收益8.7万元/日对比独立运行总收益5.7万元/日合作增量正好3.0万元/日。Shapley值和核仁分配后每个主体分到的最终收益都高于独立运行收益这就是所谓的满足“个体理性”。用这个表格做验证比你说一百句“我们的分配很公平”都管用。注意一个细节分配后总收益一栏三家都比独立运行时高但增量大小有差异。光伏A在Shapley下分到0.9在核仁下分到0.8储能B反过来。这说明不同分配规则对“谁贡献大”的判定不同。真实项目里业主往往会问“为什么这个数跟那个数不一样”这时候你要能解释背后的博弈论含义而不是丢一句“算法不同”。如果业务上对某个主体有政策倾斜可以把Shapley值再乘一个权重系数做成加权Shapley这也是常见做法。4.3 可视化与收敛性分析别只会画bar图可视化这一环节很多人的做法就是画一个柱状图对比收益能看但信息量有限。我的习惯是至少输出三张图第一张是24小时功率平衡图画出每个时段购电、光伏出力、燃气轮机出力、储能充放电和负荷的堆叠曲线这张图能一眼看出调度策略是否合理光伏大发时段有没有发生严重的弃电储能是否在谷段充电、峰段放电。第二张是各主体独立收益与分配收益对比图用成对柱状图呈现注意把独立收益和分配后收益做成两组并列而不是只画分配后的单柱否则“个体理性”这个关键信息无法从图上读出来。第三张是当主体数量较多、用蒙特卡洛采样逼近Shapley值时画出分配值随采样次数增加的收敛曲线证明你的采样次数已经足够支撑结论稳定这是在答辩和评审时非常加分的一个细节。4.4 验证分配稳定性的实操要点稳定性验证不能只看个体理性。我把项目里常用的验证项列在这里照着一项项做基本就不会被评审问倒个体理性每个主体分配后收益不低于独立收益。联盟理性对大联盟N所有主体分配收益之和等于联盟总收益不多不少。子联盟稳定性核检验没有一个子集SS中成员分到的收益之和小于S单独运行的收益v(S)。如果存在这样的S这部分主体就有动机脱离大联盟Shapley值分配结果应该在核内而核仁的结果天然在这个核内。帕累托改进正向性增量收益全部分配完毕没有剩余。实际操作中子联盟稳定性检验的计算量也不小因为你又要枚举所有S。我的做法是主体少的时候全枚举主体多的时候只针对规模小于等于3的子联盟做抽样检验大部分实际问题里小规模联盟的脱离威胁最大。5. 常见问题与排查技巧实录5.1 求解卡死或MIPGap不收敛怎么办这是出现频率最高的问题。明明模型很简单但Gurobi跑了几千秒还不收敛多半是下面几个原因二进制变量太多且约束松弛太松可行域出现大量对称的次优解Big-M参数取太大让松弛问题病态目标函数里存在巨额惩罚项导致求解器把大量时间花在反复试探不可行域。我的排查顺序是先把所有二进制变量打印出来看数量如果超过200个考虑能不能用时段聚合或者对称破缺约束砍一些然后把Big-M从1000000收紧到实际物理量的10倍左右最后在sdpsettings里加上mipgaptol, 1e-3很多时候1e-4和1e-3在实际工程里差别微乎其微但求解时间能差出一倍以上。不要一上来就追求全局最优先接受一个物理上合理的次优解再逐步收紧容差。5.2 Shapley值计算爆炸的降级方案前面已经说过子联盟数量是2的n次方n15的时候就是32768次优化求解即使每次1秒也要跑9个小时。我在实际项目里给n较大的系统用过几种降级方案按推荐程度排序优先尝试Sobol序列的低差异采样随机抽取一部分子联盟计算边际贡献再用样本均值估计Shapley值。我是对比过的Sobol比纯随机采样的收敛速度高不少同样500个子联盟Sobol的方差能小一个量级。采样同时配合Matlab Parallel Toolbox把不同子联盟的优化任务分发到各个worker上速度能再提若干倍。如果你连并行都不想配那就按业务理解把主体合并把多个同质主体合并成一个博弈方比如三台同型号燃机合成一个主体n降下来计算量直接少一个量级。这个方法兼顾精度和速度我最常推荐。5.3 分配结果不满足个体理性的排查路线如果算完之后发现某个主体分到的钱比独立运行还少先从这三个方向检查第一检查独立运行阶段的约束条件是不是比联盟运行阶段更宽比如独立运行时储能调度可以不考虑备用约束而联盟运行阶段你却给它加了备用容量约束那联盟运行收益自然偏低分配基数就错了。第二检查子联盟v(S)的口径当S不包含某主体时是否把该主体原来承担的负荷或者固定成本也错误地剔除了。第三检查是否存在部分主体之间具有强互补性比如储能和光伏的协同增量很大单独加一个其他边缘主体反而摊薄了平均边际贡献这时候增量分配不均其实是正常的需要通过核仁方案重新分配来平衡。5.4 “No suitable solver”和路径问题不少人在新电脑上配好Yalmip和Gurobi后一运行就报No suitable solver。优先检查三件事Gurobi的Matlab接口是否已经addpath到搜索路径这个最容易被忽略Yalmip更新后是否执行过yalmip(clear)旧缓存会把求解器信息锁死sdpsettings里指定的solver名称是否和你装的版本一致Gurobi新版求解器的调用名仍然是gurobi但部分旧版本Yalmip不识别最好升级Yalmip到最新版。5.5 结果波动太大或场景变了结论就翻如果你在多个典型日场景下反复跑同一套模型发现分配结论忽高忽低甚至排序规律都变了那问题大概率不在博弈分配而在上游场景输入。光伏出力曲线、负荷曲线这些参数本身具有随机性你传进去的“典型日”不够典型优化出来的结果自然不稳定。这时候我建议引入随机优化思路把场景集从1个扩到几十个每个场景带一个概率权重目标函数变成期望成本。如果嫌计算量大可以先做主成分分析聚类出几个典型场景再对每个场景单独做合作博弈分配最后按概率加权汇总。文章开头热门词里的bilstm、transformer这些序列预测方法也完全可以用来提前预测光伏出力和负荷曲线让上游输入更稳下游分配结论才能立得住。再拓展一步如果你的主体之间信息交互是分布式的不希望把数据统一交到一个中央调度中心那可以用ADMM一类的分布式求解方法。和本文的三阶段框架并不冲突只是把第二阶段“联盟优化”替换成分布式迭代流程每个主体保留自己的隐私数据通过交换边界变量迭代收敛到全局最优解。利益分配部分依然用Shapley值或核仁。最后分享一点我自己的体会做这个项目最难的不是把Matlab代码跑通而是想明白“合作收益到底从哪里来”。评审和业主都会盯着这个问题问你需要能清楚地回答收益来自多能互补、削峰填谷、减少弃风弃光、提高设备利用率这些物理来源在优化模型里都有对应的数学刻画。你是不是真的把联盟收益算清新了谈判桌上对方一眼就能看出来。把这套“三阶段建模、Shapley分配、个体理性验证”的流程完整走一遍无论你是写论文还是做项目手里的交付物都会扎实很多。