去年冬天我接了一个SCI论文复现项目课题是“基于纳什博弈的多微网主体电热双层共享策略研究”要用Matlab把整套理论从公式变成能跑的代码。说实话接之前我对纳什博弈的理解还停留在课本里“囚徒困境”的水平真正动手之后才发现论文里几行公式背后全是坑——破裂点怎么算、交易价格从哪冒出来、双层结构怎么迭代、log目标函数为什么有的求解器跑不了。这篇文章就是那次复现全过程的经验总结适合正在做综合能源系统、微网优化调度方向的研究生以及所有想把多主体博弈类SCI论文快速跑通的同行。内容上我会从数学模型、Matlab代码实现、踩坑记录三个层面展开。1. 为什么要做“电热双层共享”从三微网算例说起1.1 多微网在什么场景下需要共享电能先解释一个基本问题多个微网为什么要共享能量直接原因是资源分布不均。比如A微网屋顶光伏装得多午后大发但自己用不完B微网有一台热电联产机组冬天电和热都很充裕C微网什么都没有只能从外部电网高价买电。如果三个微网老死不相往来A会弃光B的余热只能排掉C承受高电价——三个主体都吃亏。把大家用一条联络线连起来让电能可以跨微网流动这就是最基础的能源共享。从技术上讲这件事并不复杂复杂的是经济机制A凭什么把便宜的电卖给B卖多少钱合适C加入这个“交易圈”之后大家的收益怎么重新分配这些问题本质上不是物理问题而是利益博弈问题。传统做法是让上级电网做统一调度所有微网听指令就行但现实中的微网往往分属不同业主园区、物业、产线各有各的成本结构和利益诉求不愿意无条件配合集中调度。于是多主体博弈的方法就登场了。1.2 热网加入后问题从“电力交易”升级为“电热耦合调度”电能共享相对直观一张联络线上流过功率就行。但热能共享有两个特殊性第一热能在传输管道里有损耗距离越远损耗越大所以热网往往只在同一园区或相邻微网之间铺设第二最常见的产能设备是热电联产机组CHP它一发就是“电热”捆绑输出热电比在一定范围内可调但不是完全解耦的。所以“共享电”和“共享热”根本没法分开决策。举个例子MG1有一台CHP热负荷重、电负荷轻如果只参与电力交易微网调度方会让CHP少发电避免弃电但这样一来热出力也跟着减少可能满足不了热负荷如果允许热的共享CHP就可以多发电把热卖给隔壁MG2MG2关掉电锅炉两边都划算。这就是“电热双层共享”的核心价值——通过跨微网的热力互济把CHP的潜在产能激活。1.3 “双层”到底指哪两层谁在上面谁在下面我在初读论文时对“双层”这个说法疑惑了很久后来才明白它指的是决策层次不是物理层次。上层是多微网之间的交易层决策变量是每个微网与外界的电能交易量、热能交易量以及对应的交易结算价格下层是单个微网内部的调度层决策变量是CHP出力、锅炉出力、储能充放电功率等。上下层通过交易量和价格耦合在一起上层给的交易量是下层调度的边界条件下层优化出的边际成本又决定了上层愿意接受的交易价格。可以类比成批发市场和零售店的关系上层市场定批发价和进货量下层零售店根据批发价决定上架多少货。SCI论文中常见的做法有两种一种是把下层问题的KKT条件代入上层变成单层数学规划MPEC求解另一种是用交替迭代或分布式算法上下层各算各的通过乘子交换信息。我复现的这篇论文走的是第二种思路因为它更符合多主体隐私保护的需求——每个微网都不需要公开自己的完整成本数据。2. 纳什博弈到优化问题议价模型的两步转化2.1 从非合作博弈到合作博弈为什么共享是谈判问题如果多个微网只是“各报各的价格谁接受谁交易”这就是非合作博弈典型如古诺模型、斯塔克尔伯格模型结果往往是一个纳什均衡——但这个均衡不一定是效率最高的。而“电热共享”想要的是集体合作大家商量好一个统一的调度方案和收益分配方案让每个参与者都比单干时更好这就是合作博弈。合作博弈里有几种解概念比如Shapley值、核、纳什议价解。其中纳什议价解Nash Bargaining Solution在能源系统里用得最多原因是它好算、性质清晰、且满足公平性直觉。它的核心思想可以概括成一句话一群人合伙做生意赚到的总收益怎么分才能让每个人都满意、都愿意留在合作框架内。纳什提出的解法不是平均分而是最大化所有人“效用增量”的乘积。2.2 破裂点不合作时各自能拿到多少纳什议价解有一个关键前提——谈判破裂点disagreement point, d_i也就是某个主体退出合作、选择单干时能获得的效用。在电热共享这个场景里破裂点算起来很简单让每个微网先单独优化一次自己的电热调度得到独立运行的最小成本记为 C_i^0。参与合作后微网 i 的实际运行成本记为 C_i那么合作带来的效用就是U_i C_i^0 - C_iU_i 的物理含义是“参与合作省下的钱”只有 U_i 0微网 i 才有动力留在合作联盟里。这个约束在博弈论里叫个体理性约束对应到代码里就是必须显式保证 C_i 不能超过 C_i^0。2.3 效用乘积最大化的数学转化与“蛋糕”比喻纳什议价解的标准数学形式是最大化所有参与者效用的乘积maximize ∏ U_i ∏ (C_i^0 - C_i)因为对数函数是严格递增的乘积最大化等价于对数之和最大化所以目标函数可以写成maximize Σ ln(C_i^0 - C_i)这个转化有两个好处第一把乘积变成了求和很多数学性质和求解器特性可以直接套用第二对数函数天然惩罚“某一个微网一点好处都没捞到”的情况——如果某个 U_i 趋近于0对数值会趋于负无穷所以算法会自动避开那种让某个参与者吃亏的方案保证分配的公平性。用“蛋糕”打比方总成本下降的量是整块蛋糕纳什议价解既要让整块蛋糕尽量大总效用高又要把蛋糕分得让每个人都觉得划算各 U_i 尽量均匀且不小于零。如果一块蛋糕分下来有人只拿到一口对数乘积就会大幅缩水算法会去调整分配方案让那个人多拿一点——这就是纳什议价解区别于“单纯最小化总成本”的地方。3. 上层竞价与下层调度耦合变量与迭代机制3.1 上层交易层交易变量和结算价格怎么定义在“双层共享”框架里上层的核心变量是各微网之间的电能交易量 P_ij 和热能交易量 H_ij以及交易结算价格 λ。理论上价格应该由市场出清决定但论文里常用一个技巧把交易价格定义为下层调度问题关于交易量的对偶变量拉格朗日乘子。这个乘子的经济含义是“该微网增加一单位交易量时其内部最优成本的变化幅度”本质上就是边际成本。这样做的好处是无需人工定电价价格从每个微网的物理约束和成本函数里自然涌现出来。当某个微网向另一个微网购电时成交价格不会超过买方从外部电网购电的成本也不会低于卖方自己发电的边际成本否则至少有一方不愿交易。这种价格天然落在可行区间内避免了扯皮。3.2 下层调度层每个微网内部的电热联供模型下层是经典的微网经济调度问题。每个微网内部包含若干设备我复现的论文里用到的典型模型有以下几种CHP机组消耗天然气同时输出电功率 P_chp 和热功率 H_chp。简化模型里热电比是常数 r H_chp / P_chp更精细的模型会允许热电比在一定范围内调节。成本是燃料消耗量乘以天然气价格。燃气锅炉只产热成本同样来自天然气是CHP的补充热源。电锅炉消耗电能产热在风电或光伏充足的时段用起来很划算。电储能充放电功率、容量约束、SOC递推方程以及充放电不能同时进行的逻辑约束。热储能类似电储能不过对应的变量是热功率。每个微网的约束条件包括电功率平衡、热功率平衡、设备出力上下限、储能SOC上下限、爬坡约束等。把这些约束组合起来每个微网就是一个混合整数线性规划或凸规划问题——说“混合整数”是因为储能充放电和机组启停需要0-1变量。3.3 ADMM为什么天然适合纳什议价问题双层模型需要求解的是一个耦合问题但纳什议价目标函数有一个非常漂亮的数学性质它是关于各个 U_i 的“和”。求和结构让各微网的子问题可以完全解耦只剩交易变量需要协调。这种结构正是交替方向乘子法ADMM的主场。ADMM的基本流程是每个微网在给定交易价格和交易量的情况下只优化自己的内部调度问题然后上层收集所有微网的交易意愿取平均或做市场化出清更新交易计划再把新的交易计划传下去用拉格朗日乘子更新交易价格。如此循环直到交易计划在微网间达成一致。我复现时的迭代逻辑大概是这样初始化交易变量 z^0、乘子 λ^0、惩罚系数 ρ第 k 轮各微网并行求解给定 λ^k 和 z^k每个微网独立优化自己的成本 C_i 和期望交易量 x_i更新共识交易变量z^{k1} 所有微网 x_i 的平均值或通过一个小的市场出清模型求解更新乘子λ^{k1} λ^k ρ * (x_i^{k1} - z^{k1})计算原残差和对偶残差判断是否收敛。这套流程在论文里展示时只有半页伪代码但落到Matlab里你会发现步长 ρ 怎么取、收敛阈值怎么定、交易变量要不要做可行性投影全是影响成败的细节。4. Matlab代码落地的四个关键模块4.1 模块一数据初始化与算例生成代码的第一步是把论文里的算例数据搬到Matlab里。我建议用结构体组织数据每个微网一个字段设备参数和负荷曲线单独放。比如%% 数据初始化 T 24; % 调度周期小时 price_elec [0.48*ones(1,8), 0.92*ones(1,8), 0.48*ones(1,8)]; % 外部购电价 price_gas 2.3; % 天然气价格元/立方米 % 微网1CHP 光伏 电负荷 热负荷 MG(1).P_chp_max 600; MG(1).P_chp_min 50; MG(1).r_chp 1.2; % 热电比 MG(1).eta_chp 0.45; % 发电效率 MG(1).P_pv [0 0 0 0 20 50 150 250 280 260 220 180 160 180 200 220 180 100 50 20 0 0 0 0]; MG(1).P_load 450 80*sin((1:T)/24*2*pi) randn(T,1)*8; MG(1).H_load 320 60*cos((1:T)/24*2*pi) randn(T,1)*5; % 微网2、微网3类似省略这段代码要注意的点是所有变量单位统一为kW和kWh热负荷与电负荷的曲线要有区分度否则后面结果没有对比价值。论文复现的第一步不是写优化模型而是先把场景调成“不合作时有痛点、合作后有收益”否则模型跑出来总成本压根不降后续博弈就无从谈起。4.2 模块二独立运行基准模型破裂点求解破裂点 C_i^0 的求解是个纯单微网优化问题不需要任何交易变量。我把这一步封装成一个函数输入微网数据和外部能源价格输出独立运行成本 C0 和调度结果。%% 微网独立运行优化 function [C0, result] solve_standalone(MG) P_chp sdpvar(24,1); H_chp sdpvar(24,1); P_dis sdpvar(24,1); % 储能放电 P_ch sdpvar(24,1); SOC sdpvar(24,1); u_ch binvar(24,1); % 储能充电状态 % 约束 Cons []; Cons [Cons, P_chp MG.P_chp_min, P_chp MG.P_chp_max]; Cons [Cons, H_chp MG.r_chp * P_chp]; Cons [Cons, SOC(1) 0.2, SOC(2:end) SOC(1:end-1) 0.9*P_ch/0.95 - P_dis/0.95]; Cons [Cons, SOC 0.1, SOC 0.9]; Cons [Cons, P_ch 200*u_ch, P_dis 200*(1-u_ch)]; Cons [Cons, MG.P_pv P_chp P_dis - P_ch MG.P_load]; Cons [Cons, H_chp MG.H_load]; % 假设热由CHP满足 obj sum(price_gas / MG.eta_chp) * P_chp; % 燃料成本 optimize(Cons, obj, sdpsettings(solver,gurobi)); C0 value(obj); end注意独立运行时如果某些微网只有电锅炉而没有CHP热负荷就得靠电锅炉产热成本会很高。这正是合作后“让MG1的CHP多产热供给MG2”的价值来源。我建议先跑一遍独立运行把每个微网的成本打印出来看看谁最“苦”谁最有富余产能后面解读结果时就一目了然。4.3 模块三纳什议价主模型与求解器选择集中式纳什议价模型是复现代码中最容易出问题的环节。目标函数是max Σ ln(C_i^0 - C_i)在Matlab中我习惯用YALMIP建模但这里有一个非常关键的技术细节ln不是线性函数也不在默认的LP/QP/SOCP框架内它属于指数锥约束。要直接求解这个目标函数求解器必须支持指数锥比如Gurobi9.0以上版本和SCS。如果沿用CPLEXYALMIP要么报错要么会用线性逼近替代——得到的结果完全不靠谱。我用的集中式代码结构如下%% 集中式纳什议价模型用Gurobi求解 C_var sdpvar(N, T); % 每个微网每个时段的运行成本 Cons []; for i 1:N % 各微网的内部约束和独立运行时类似不过多了交易变量P_trade和H_trade Cons [Cons, P_local(i,:) MG(i).P_pv P_chp_i ..., P_trade(i,:)]; Cons [Cons, C_var(i,:) 燃料成本表达式]; end % 交易耦合约束电能和热能在联盟内部守恒 Cons [Cons, sum(P_trade, 1) 0]; Cons [Cons, sum(H_trade, 1) 0]; % 个体理性约束 Cons [Cons, C_var_sum(i) C0(i)]; obj_nash -sum(log(C0 - sum(C_var,2))); % YALMIP默认求最小化取负 optimize(Cons, obj_nash, sdpsettings(solver,gurobi));但这里有个实际问题C0 - sum(C_var,2)在迭代过程中可能因为数值误差变成极小的负数或零导致log出现NaN。我建议先用一个让所有微网合作的“成本下界”来约束 C_var 的范围或者在目标函数里加一个小量 epsilon比如1e-6保证对数输入严格大于零。更稳妥的做法是把目标改成epsilon 1e-3; obj_nash -sum(log(C0 - sum(C_var,2) epsilon));这样能有效避免数值问题又不影响最优解。4.4 模块四ADMM迭代求解与收敛判据集中式求解适合小算例和结果校验但论文复现真正要展示的是分布式求解过程——各微网不需要向中心节点公开内部参数。ADMM在这个环节的价值在于每个微网的子问题独立建模、独立求解中心节点只接收交易变量的值。ADMM的核心循环代码%% ADMM迭代求解 rho 0.5; % 惩罚系数 max_iter 200; tol_pri 1e-4; tol_dual 1e-4; x_prev zeros(2*T, N); % 每个微网的交易变量电能和热能 z zeros(2*T, 1); % 共识交易变量 lambda zeros(2*T, N); % 对偶乘子 for k 1:max_iter % Step 1: 各微网并行求解子问题 for i 1:N x(:,i) solve_ADMM_subproblem(i, z, lambda(:,i), rho, MG, C0); end % Step 2: 更新共识变量交易计划平均值 z_new mean(x, 2); % Step 3: 更新对偶乘子 lambda lambda rho * (x - z_new); % Step 4: 收敛判断 r_k norm(x - z_new, fro); s_k rho * norm(z_new - z); if r_k tol_pri s_k tol_dual break; end z z_new; end子问题的目标函数里除了自身成本还要加上ADMM引入的二次项%% 微网i的ADMM子问题 function [x_i, C_i] solve_ADMM_subproblem(i, z, lambda_i, rho, MG, C0) % 本地变量定义省略... C_i ...; % 本地成本 % 交易变量 x_i 与共识变量 z 的偏差惩罚 obj_local C_i lambda_i * x_i (rho/2) * sum((x_i - z).^2); optimize([Cons_local, x_i 交易量的定义], obj_local, sdpsettings(solver,gurobi)); x_i value(x_i); end这套东西跑起来之后我最大的感受是ADMM的迭代路径很“粘”前期几乎看不到交易量有明显变化但一旦过了某个临界点收敛速度会突然变快。工程上建议每轮迭代都把 r_k 和 s_k 画出来一旦看到残差曲线出现“先平后陡”的形状就说明离收敛不远了。5. 复现踩坑记录五个真正让人崩溃的细节5.1 破裂点算错的三种典型姿势破裂点 C_i^0 是纳什议价模型的基石但代码里非常容易算错。我见过三种典型错误第一种是没关掉交易变量就直接跑联合优化把合作后的成本当成了破裂点导致 U_i 全是零或负数第二种是独立运行模型里漏掉了热负荷约束比如某个微网只有电锅炉代码里忘了加电锅炉产热的平衡方程算出来的 C_i^0 低得离谱第三种是用了外部电网购电价格但忘了计入网络损耗导致独立运行成本偏低。避免方法就一条把每个微网的独立运行模型单独跑一遍人工检查热功率平衡和电功率平衡是否严格成立。C_i^0 算错了后面的对数目标函数全部失去意义而且这个错很隐蔽——程序不会报错只是结果看起来“总是合作无效”。5.2 买电和卖电同时非零的“幽灵交易”在集中式模型里如果不加约束优化器可能会出现微网 i 同时从微网 j 买电、又向微网 k 卖电的“幽灵交易”。物理上这是不可能的同一时段同一节点的净交易量只有一个方向但数学上如果买卖价差存在套利空间优化器就会这么干。解决办法是引入互补约束% P_buy 和 P_sell 不同时非零 P_buy 0; P_buy M * u_buy; P_sell 0; P_sell M * (1 - u_buy);这里 u_buy 是0-1变量M 是足够大的常数。但在ADMM分布式框架里这种双变量建模会让子问题变成混合整数规划收敛性更难保证。一个更干净的替代方案是只用净交易变量 P_trade正为买、负为卖然后在目标函数里同时考虑买入成本和卖出收益用价格差自然抑制“又买又卖”的行为。我最后采用了后者求解效率明显提升。5.3 log目标函数在CPLEX下静默出错这是我踩过最深的坑。当时我偷懒沿用之前的CPLEX求解器YALMIP没有报错模型“正常”跑完但结果显示每个微网的效用都几乎等于零整个合作毫无意义。逐行排查才发现CPLEX根本不支持指数锥约束YALMIP在底层把 log 目标近似成了一段多段线性函数而且近似精度很差。正解是换用Gurobi版本9以上或SCS求解。如果你用的是CVX可以用 exponential cone 形式重写目标cvx_begin variable U(N) maximize(sum(log(U))) subject to U C0 - C_var; U 0; cvx_endCVX会自动选择支持指数锥的求解器如SCS。经验之谈不要为了省事沿用旧求解器数学优化领域“求解器不支持”不是小问题它会静默地毁掉你的结果。5.4 CHP热电比约束写错导致热平衡大面积失衡CHP机组的模型是H_chp r_chp * P_chp看起来简单但很多人忽略了一个问题论文里的热电比有时是基于热功率与电功率的比值有时是基于一次能源分配的比值。如果抄错就会出现“电平衡满足但热平衡完全对不上”的怪现象。还有一个细节CHP有一个运行包络不是任意电出力都能对应任意热出力。低负荷时热电比往往偏高。论文通常会画一个可行域代码里我建议至少加上热出力上下限和电出力上下限两个约束而不是只用一个固定热电比约束% 热电比范围约束变工况简化模型 H_chp r_min * P_chp; H_chp r_max * P_chp;这样既能保留CHP电热耦合的基本特性又不会因为固定热电比把可行域限制得太死。另外CHP的热出力和电出力都应当是连续变量不要图省事把它设成逐时段离散变量否则求解规模会爆炸。5.5 ADMM收敛判据只写一半的翻车现场刚开始写ADMM时我只看了原始残差 r_k没看对偶残差 s_k。结果有一次迭代到了1000轮r_k 已经小于1e-6但 z 和 λ 还在缓慢漂移。原因是对偶残差控制的是乘子更新的幅度如果只控制 r_k乘子可能在最后阶段振荡但幅度很小导致共识变量一直没有完全一致。收敛判据我后来固定用双残差r_k norm(x - z_new, fro); s_k rho * norm(z_new - z); tol_pri sqrt(2*T*N) * 1e-4 1e-4 * max(norm(x,fro), norm(z_new,fro)); tol_dual sqrt(2*T*N) * 1e-4 1e-4 * norm(lambda, fro);此外ADMM的惩罚系数 ρ 也很敏感。ρ 太小原残差收敛慢ρ 太大对偶残差振荡。我建议做一次 ρ 的扫描实验比如取0.1、0.5、1、2画出迭代曲线对比选一个收敛快且稳定的值。工程论文里经常直接用固定 ρ实际复现时自适应调整 ρ 会省下大量时间。6. 三微网24小时算例从结果反推模型是否正确6.1 算例参数与场景描述我用一个三微网系统做测试参数设置尽量贴近实际微网电源设备热源设备储能配置电负荷峰值(kW)热负荷峰值(kW)MG1CHP 600kW 光伏200kWCHP余热无500320MG2风电400kW电锅炉200kW无700250MG3光伏300kW无电储能200kWh 热储能300kWh450150三个微网通过一条公共联络线共享电能铺设了连接MG1与MG2、MG1与MG3的热力管道。天然气价格2.3元/立方米外部购电采用分时电价峰值0.92元/千瓦时谷值0.48元/千瓦时。调度周期24小时步长1小时。这个设置的意图很明确MG1有CHP热产能富余是天然的供能方MG2有电锅炉但自身风电不足热成本很高MG3完全没有热源只能从热网买热或依靠热储能。6.2 核心结果成本下降与效用分配三微网独立运行的总成本是9650元其中MG1约3800元、MG2约4200元、MG3约1650元。合作后总成本降到8580元左右下降约11.1%即省下约1070元。各微网的效用分配如下MG1合作后成本约3300元节省500元效用占比最大。原因是CHP多发电、多发热虽然燃料成本增加了但对外售电、售热的收入远超成本增量。MG2合作后成本约3800元节省400元。它从MG1购买热减少电锅炉的使用节约了大量电费同时风机发电可以卖给MG3。MG3合作后成本约1480元节省170元。它没有热源原本热负荷靠电储能转化折合成本很高合作后直接从MG1的热网买热热储能也得以灵活调配。纳什议价解的一个典型特征是效用最高的参与者不是“最穷的”也不是“最富的”而是边际供给能力最强的那一个。MG1愿意多出力是因为它每多生产一单位热和电自己的边际成本增幅小于市场结算价格MG3虽然配合度高但自身调节能力有限所以效用占比最小。从结果看三个微网成本都下降了个体理性约束全部满足。6.3 结果校验三件套个体理性、帕累托性、与集中式解的一致性复现工作不能跑完出数就算完还得验证结果不是“随机碰巧”的。我常用的校验方法有三个第一个是个体理性校验。检查是否每个 C_i 都严格小于 C_i^0如果出现某个微网成本反而上升说明约束写漏了或者目标函数方向有问题。一般要求最小的 U_i 也要大于总节省金额的5%否则那个微网随时可能退出合作。第二个是帕累托最优校验。用集中式模型跑一个“所有微网作为一个整体最小化总成本”的下界模型得到理论最小总成本 C_sum_lb。合作后总成本应该接近这个下界差距通常小于1%因为纳什议价解本质上是在帕累托前沿上找一个“收益分配公平”的点总成本会略高于纯整体最优但不会差太多。如果合作后总成本离 C_sum_lb 太远说明双层迭代没收敛到位。第三个是与集中式解的一致性校验。我用Gurobi跑集中式纳什议价模型再和ADMM分布式结果对比。如果两者给出的总成本差异超过1%就该怀疑ADMM的收敛判据没有写对或者子问题的二次惩罚项系数不合适。论文里展示的“分布式求解与集中式求解结果几乎一致”就是这么验证出来的。这套代码跑通之后我再回头看那篇论文发现它最难复现的其实不是纳什博弈本身而是“算完结果之后怎么证明你算的对”这一环。很多同学拿到代码跑出一个漂亮的成本下降曲线就觉得大功告成实际上没有经过个体理性和集中式一致性校验的结果大概率只是“看起来能用”而已。最后再分享一个实操中的小习惯我会在Matlab脚本开头加上一句rng(0)把随机数种子固定这样负荷数据、光照数据的随机波动每次都能复现调参时不会因为场景变了而搞不清是代码问题还是数据问题。如果你也正在折腾多微网纳什博弈的复现希望上面的经验和坑能帮你少走几个弯路。
