☰
基于纳什博弈的多微网电热双层共享策略Matlab复现全解析
2026/9/24 21:47:56 网站建设 项目流程

如果你论文库里放着一篇题目很长、公式很多的SCI,要求你用Matlab把里面的多微网调度结果一模一样地复现出来,你会从哪一行代码开始?我先说我的答案:我不会先打开Matlab,而是先拿一张A4纸,把“纳什博弈”“多微网”“电热双层共享”这几个关键词拆开,想清楚每个概念对应代码里的哪个变量、哪个函数、哪个循环。

我最近就按这个思路复现过一个基于纳什博弈的多微网电热双层共享策略,这个过程比我想象的费时间,但也把这类论文的套路摸清楚了。很多所谓SCI复现的难点,并不在代码本身,而在于论文的数学模型和Matlab表达之间那条缝。这篇博客就是我的复现笔记,把从拆公式到写代码、调参、验证纳什均衡的完整链路记录下来,给后面做博弈论多主体调度、电热耦合系统或分布式优化的同学一个能直接上手的方法。

1. 先弄清论文在解决什么问题:共享策略的博弈视角

1.1 多微网为什么要“共享”电和热

这种文章的背景通常是:一个区域里分布着好几个微网,每个微网有光伏、风机、热电联产机组、电锅炉、储能或者蓄热罐。单独运行时,有的微网白天电多热少,有的晚上热多电缺。如果各微网之间可以共享电和热,就能减少从外部电网买电、减少燃气锅炉耗气,整体经济运行指标会变好。

“共享”听起来简单,落在数学模型上就是一组交易变量。以电为例,共享意味着微网A在一个时段既可以向微网B卖电,也可以从微网C买电,电力潮流通过公共联络线在微网之间流动。热力共享类似,只不过共享的是热功率,通常通过区域热网或者站房间的热力管道连接。电热双层共享,就是在同一个优化框架里同时考虑电力交易和热力交易,并且两者不是独立的,因为CHP机组可以把燃气同时变成电和热,电锅炉又能把电变成热。

复现的第一步就是要把这些设备关系理清楚:哪些设备是电源,哪些是热源,哪些是电热耦合设备。我用一张简单的表记录每个微网内部设备的输入输出关系,这是后面写约束条件的基础。

设备输入输出耦合关系
CHP机组天然气电、热电热之比由热电比决定
光伏/风电可再生能源电出力上限来自预测曲线
电锅炉电热用电功率按效率转化为热功率
电储能电电充电功率、放电功率和SOC递推
蓄热罐热热储热功率、放热功率和蓄热状态递推

这张表一旦列出来,论文里那些公式的物理含义就通了一大半。

1.2 “双层”到底指什么

“电热双层共享”里的“双层”需要先判断是指两层能量载体,还是指双层优化结构。我复现的这篇论文,主要指电和热两个品类的共享层,也就是电力共享层和热力共享层。电力共享包括微网之间通过联络线的功率交换和电价形成机制,热力共享包括通过热网管道传递的热功率交换和热价形成机制。

这两个“层”虽然可以放在同一个模型里写,但处理逻辑很不一样。电力共享对时间尺度的敏感度很高,功率平衡要求瞬时满足,线路容量也会限制交换功率;热力共享则有明显的热惯性,管道里的热水不是瞬时到达的,温度梯度、热损失、传输延时都可能被建模进去。如果原论文用的是热功率平衡而没有细化热网传热,复现起来相对简单;如果包含节点温度混合和管道动态,就需要额外引入网络变量。

另外还要注意,不要把“双层”硬套成“上层优化、下层优化”的Stackelberg模型。虽然确实有论文是三层结构,但常见的多微网电热共享是单层博弈决策:每个微网同时决定自身设备出力和交易功率,共享价格在迭代中形成。拿到论文先找“upper level”“lower level”或者“layer”出现的位置,再决定你的代码主循环应该是嵌套优化还是简单迭代。

1.3 为什么用纳什博弈而不是集中优化

这是做复现之前必须想明白的问题。集中优化的做法是把所有微网的变量扔进同一个优化问题,以整个系统总成本最小化为目标,求解一次得到全局最优调度。这个方法从数学上最干净,代码也最好写,但它有一个致命前提:所有微网愿意把自己的负荷、设备参数、成本信息完全交给一个中心调度者。

实际场景中,不同微网可能是不同运营商,彼此之间既有合作又有竞争,数据隐私和利益分配都很敏感。纳什博弈的思路是反过来:每个微网是一个独立的主体,只关心自己的运行成本和收益,在给定其他微网策略的前提下,选择让自己收益最大的策略。当所有微网都没有动机单独改变自己的策略时,就达到纳什均衡。

复现时,这意味着不能把所有目标函数相加求最小,而是要把“每个微网独立最小化自身目标”写成 N 个子问题,再用共享价格或者共享功率把这些子问题连接起来。这种结构上的差别,决定了你的代码核心是循环嵌套优化,而不是一次矩阵求解。

2. 复现前的模型落地:从数学公式到计算对象

2.1 微网内部设备与能量平衡

开始写Matlab之前,我建议先手写出每个微网的变量、约束和目标函数,格式越接近矩阵就越能减少后面的调试时间。

一个典型微网 i 在时段 t 的主要变量有:

  • CHP输出电功率 P_chp(i,t),输出热功率 H_chp(i,t)
  • 电锅炉消耗电功率 P_eb(i,t),输出热功率 H_eb(i,t)
  • 电储能充电功率 P_ch(i,t)、放电功率 P_dis(i,t)
  • 蓄热罐储/放热功率 H_st_ch(i,t)、H_st_dis(i,t)
  • 与外部共享的购电功率 P_buy(i,t)、售电功率 P_sell(i,t)
  • 与热网交换的购热/售热功率 H_buy(i,t)、H_sell(i,t)

在这些变量基础上,每个时段都要满足两个平衡约束。电力平衡的物理含义是:光伏出力、CHP电出力、电储能放电和购电之和,等于本地电负荷、电锅炉用电、电储能充电和售电之和。热力平衡的物理含义是:CHP热出力、电锅炉热出力、蓄热罐放热和购热之和,等于本地热负荷和售热之和。

在Matlab里,这类约束就是等式约束矩阵 Aeq 的某一行。建议把每个微网的变量排成一个统一顺序的向量,比如 x = [P_chp; H_chp; P_eb; P_dis; P_ch; H_st_dis; H_st_ch; P_buy; P_sell; H_buy; H_sell],然后所有约束都按这个顺序写。这样写子问题代码时不会乱。

2.2 电力共享层与热力共享层的耦合关系

电热双层共享最有意思的地方是两个共享层不是简单并列存在,它们通过CHP机组和电锅炉耦合在一起。CHP机组的热电比决定了一个微网在产电的同时必然产热,如果本地热负荷小,多余的热可以卖给其他微网或者由蓄热罐吸收。电锅炉正好反过来,如果某个微网电多了但热不够,可以用电锅炉把富余的电转化成热,从而减少从热网购买热量。

因此,在写目标函数时不能把电和热分开优化两次,必须放到同一个子问题里同时决策。这也是为什么用分解算法求解时要特别小心:如果把电热耦合变量强行掰开,迭代可能不收敛。

耦合关系还体现在共享价格上。电力共享价格会影响CHP机组的发电倾向,发电倾向又会影响CHP产热,CHP产热接着影响该微网在热力共享层的买卖热策略。我复现的时候明显感觉到,电热价格必须放在同一个迭代框架里更新,分开更新很容易出现“电收敛了,热还在振荡”的情况。

2.3 目标函数、可行域和边界条件的矩阵化

典型的多微网经济调度目标是一个二次函数加线性项,常写为:

min sum_t [ a_i * P_chp(i,t)^2 + b_i * P_chp(i,t) + c_i + C_om + 购能费用 - 售能收益 ]

其中购能费用包括购电费用 lambda_e * P_buy 和购热费用 lambda_h * H_buy,售能收益类似。在Matlab中,如果所有变量连续且目标二次,就直接写成 x'Qx + f'x,用 quadprog 求解。如果有机组启停的0-1变量,就要用 intlinprog,或者用YALMIP调用Gurobi/CPLEX。

约束条件里最容易遗漏的是CHP可行域。很多复现代码只写了出力上下限,忽略了CHP电出力和热出力之间的可行域约束,比如热出力不能超过当前电出力对应的最大热输出,也不能低于最小热电比对应的热输出。这类约束写成线性不等式就是 C*x <= d 的一行。

还要注意联络线容量约束,也就是 P_buy + P_sell 的上下限。有些共享策略要求某个微网在同一时段不能同时买卖电力,这个约束在连续模型里其实很难严格表达,常见做法是用一个很小的时间尺度或者引入0-1变量。复现时如果论文没有明确要求,直接用联络线容量限制即可,不必额外加二进制变量,否则求解难度会明显上升。

3. 纳什博弈求解框架与Matlab迭代实现

3.1 把每个微网看成独立优化子问题

纳什博弈的数值实现,本质上就是反复求解一组互相耦合的独立优化问题。每个微网 i 的优化子问题可以写成:

给定其他所有微网的共享变量后,微网 i 最小化自己的目标函数,同时满足自己的设备约束和共享变量约束。

这就有两个细节需要注意。第一,共享变量既可以是共享功率,也可以是共享价格。如果是功率共享,微网 i 决策卖出多少功率,会影响其他微网能买到的功率,这种问题通常属于广义纳什均衡问题,需要额外处理共享约束。如果是价格共享,每个微网面对同一组价格信号做独立决策,再把所有微网的决策汇总起来更新价格,这种实现更简单,也是复现时最常见的做法。

我建议你复现时先看论文里的“价格”是不是真的作为变量存在。如果论文有拉格朗日乘子的影子,大概率属于价格驱动型共享,实现起来压力小很多。

3.2 Gauss-Seidel风格的迭代流程与价格更新

我复现时用的主循环是经典的三步走:求解子问题、汇总共享不平衡、更新共享价格。

开始时先给共享价格一个初值,比如所有时段统一取当地电网购电价格和天然气折算热价的某个比例。然后进入迭代:

  1. 对每个微网 i,在当前共享电价 lambda_e 和共享热价 lambda_h 下,调用子问题优化函数,得到该微网的购电/售电、购热/售热和设备出力。
  2. 汇总所有微网的共享电功率和共享热功率,计算系统层面的售购偏差。
  3. 根据偏差修正共享价格,偏差过大时抬高价格或压低价格。
  4. 判断价格和策略是否同时收敛,不满足收敛条件就回到第2步。

这个框架很像拍卖市场里的“价格叫价”过程。用Matlab写成伪代码就是:

for k = 1:max_iter lambda_e_old = lambda_e; lambda_h_old = lambda_h; for i = 1:N [x{i}, obj_i(i)] = solve_microgrid(params{i}, lambda_e, lambda_h); Psell(i,:) = x{i}.Psell; Pbuy(i,:) = x{i}.Pbuy; Hsell(i,:) = x{i}.Hsell; Hbuy(i,:) = x{i}.Hbuy; end imbalance_e = sum(Psell,1) - sum(Pbuy,1); imbalance_h = sum(Hsell,1) - sum(Hbuy,1); lambda_e = lambda_e + alpha_e * imbalance_e; lambda_h = lambda_h + alpha_h * imbalance_h; if norm(lambda_e - lambda_e_old, inf) < tol && ... norm(lambda_h - lambda_h_old, inf) < tol break; end end

这里如果所有微网都基于上一轮的价格同时求解,就是Jacobi式并行更新;如果微网按顺序求解,并且后面微网能看到前面微网刚更新的结果,就是Gauss-Seidel式异步更新。后者通常收敛更快,但要求你严格定义“共享信息”的更新顺序,稍不注意结果会有差异。

3.3 收敛判据与纳什均衡验证

很多初学复现的同学只看前后两轮的价格差,价差小于阈值就算收敛。这样其实不够。价格收敛不代表每个微网都已经达到最优反应,尤其当目标函数不是严格凸的时候,价格可能到了一个稳定影子但策略还在小幅波动。

我的做法是同时检查三个条件:

  • 共享价格变化量小于阈值
  • 共享功率偏差绝对值之和小于阈值
  • 每个微网的目标函数变化量小于阈值

前两个条件判断市场出清情况,第三个条件判断最优反应是否稳定。更严格的做法是在迭代结束后做一次“单边偏离检验”:把其他微网的共享策略固定住,重新优化某一个微网,看它的目标函数是否还能下降。如果下降量明显超过设定误差,说明当前结果不是纳什均衡,需要继续迭代或者调整初值和步长重新求解。

单边偏离检验是论文里不会细写、但审稿人经常问的东西。复现代码里加上这一段,在写报告或者回答问题时会非常加分。

4. Matlab代码结构模块拆解与关键函数实现

4.1 参数初始化与输入数据结构

复现这种论文,最忌讳把所有东西都塞进一个脚本。我的代码组织方式是一个主脚本加几个功能函数,数据全部用结构体传入引出。

参数初始化的函数大概这样:

function params = init_params(T, N) params.T = T; % 调度时段数 params.N = N; % 微网数量 for i = 1:N params.mg(i).Pchp_max = 2.0; % MW params.mg(i).Pchp_min = 0.3; params.mg(i).Hchp_max = 1.6; params.mg(i).Peb_max = 0.8; params.mg(i).Peb_eff = 0.95; params.mg(i).load_e = load_data(i, T); params.mg(i).load_h = heat_data(i, T); params.mg(i).pv = pv_data(i, T); params.mg(i).R_line = 0.05; end end

用结构体的好处是,后面写子问题优化函数时,只需传 params.mg(i) 进去,函数内部不会污染外部变量,也方便批量修改微网数量。所有输入曲线统一用24维列向量,时段从1到T。

4.2 子问题优化函数:用quadprog还是YALMIP

针对单个微网的子问题,我提供一个用quadprog求解的简化版本。这个版本忽略储能和0-1变量,只保留核心的电热耦合设备,目的是展示变量排列和约束矩阵写法。

function [x, obj] = solve_microgrid(par, lambda_e, lambda_h) % 变量顺序: % 1 Pchp, 2 Hchp, 3 Peb, 4 Hsell, 5 Hbuy, 6 Psell, 7 Pbuy H = zeros(7); H(1,1) = 2 * par.a_chp; % 二次成本系数 f = [par.b_chp; par.c_h; 0; -lambda_h; lambda_h; -lambda_e; lambda_e]; Aeq = [1 0 -par.Peb_eff 0 0 -1 1; % 电平衡 0 1 par.Peb_eff -1 1 0 0]; % 热平衡 beq = [par.load_e - par.pv; par.load_h]; A = [1 -par.ratio_chp_max 0 0 0 0 0; -1 par.ratio_chp_min 0 0 0 0 0]; b = [0; 0]; lb = [par.Pchp_min; 0; 0; 0; 0; 0; 0]; ub = [par.Pchp_max; par.Hchp_max; par.Peb_max; ... par.Hsell_max; par.Hbuy_max; par.Psell_max; par.Pbuy_max]; options = optimoptions('quadprog','Display','off'); [x, obj] = quadprog(H, f, A, b, Aeq, beq, lb, ub, [], options); end

如果论文模型里含0-1变量,比如储能的充放状态、机组启停、买卖互斥状态,quadprog就不适用了。我建议直接用YALMIP建模,求解器选Gurobi或者CPLEX。YALMIP的表达方式更接近论文公式,能大幅降低建模出错率,缺点是需要额外安装工具箱。Matlab自带的intlinprog也可以做,但变量一多,性能差距就出来了。

4.3 主迭代脚本、数据流与结果导出

主脚本的内容就是前面写的迭代循环。每轮迭代把每个微网的优化结果存到一个 cell 数组里,既有按微网索引的一层,也有按迭代次数的另一层。这样调试时可以很方便地画出“不同迭代轮次下各微网买卖电量”的变化过程,对判断收敛非常有帮助。

数据流上,我会在每个子问题函数内部只返回结构体,不返回一大堆散变量。比如:

sol.Pchp = x(1); sol.Hchp = x(2); sol.Peb = x(3); sol.Hsell = x(4); sol.Hbuy = x(5); sol.Psell = x(6); sol.Pbuy = x(7); sol.obj = obj;

主脚本汇总后,再把结果用 writetable 或者 writematrix 导出到Excel,方便和论文里的数据表逐项对照。

4.4 对标论文图表的可视化代码思路

SCI复现的最后一关永远是图表。很多同学的代码能算出结果,但画出来的图跟论文对不上,原因往往是横纵轴物理量不对、坐标范围不一致、或者某一组曲线用的是“净功率”而不是“购售功率”。

我画图时习惯把每类变量分成一张图。以电功率共享为例,用堆叠面积图表示各微网的购售电曲线,用柱状图表示某个典型时段的共享功率分配,用折线图展示共享电价收敛过程。热力共享类似,但要注意热功率单位是kW还是MW,不统一的话图线形状完全不一样。

Matlab里常用到的就是 plot、stairs、area、bar 这几个函数。关键在于设置线宽、字体、坐标轴范围时尽量向论文靠拢,同时保留一个能一键输出 .fig 和 .png 的导出函数。复现论文时不需要完全复刻配色,但趋势和数值必须对得上。

5. 复现过程中最容易翻车的五个细节

5.1 初值设置:纳什博弈对初始点非常敏感

这是我踩得最深的一个坑。集中优化问题通常只要可行域是凸的,初始点不影响全局最优解;纳什博弈不一样,不同共享价格初值有可能收敛到不同均衡,甚至不收敛。

我试过把共享电价初值设成0,结果迭代前几步各微网都拼命卖电,价格直接冲上天,然后开始振荡。后来我改用当地电网分时电价作为初值,并在电价更新方程里加入阻尼系数,收敛情况明显改善。

遇到不收敛时,先不要怀疑算法框架写错,先检查价格初值是否落在接近市场出清价的区间。可以先用集中优化求一次最优调度,把最优共享功率对应的拉格朗日乘子拿出来当初始价格,纳什博弈的收敛速度会快很多。

5.2 共享变量更新顺序不同,结果完全不同

Gauss-Seidel式更新如果顺序乱了,结果可能对不上原论文。有些作者默认每个微网按编号顺序决策,后面的微网在做优化时能看到前面微网的最新共享功率;有些作者则用Jacobi式并行更新,所有微网都基于上一轮的信息。两种更新方式在最优点一致的情况下也可能产生不同的收敛路径。

复现时最好在论文方法部分找到对更新过程的描述,哪怕没有明确写,也可以通过对比实验结果反推。如果论文用了“parallel”或者“simultaneously”,就用Jacobi;如果用了“sequential”或者“order”,就用Gauss-Seidel。

5.3 价格变量发散:阻尼更新是常规解

价格更新不是越大步长越快收敛。步长太大,价格会在均衡点附近来回震荡;步长太小,收敛需要几百轮迭代。我使用的步长通常取0.01到0.05之间,具体大小跟目标函数系数和功率量纲有关。

更稳妥的做法是引入阻尼系数,即每轮更新后的新价格等于上一轮价格和本轮计算价格的加权平均:

lambda_new = beta * lambda_calculated + (1 - beta) * lambda_old

beta 取0.5到0.8。这样能显著抑制振荡。我一般把step和阻尼系数都设为可调参数,调试时用二分思想找到适合当前数据规模的值。

5.4 中文注释乱码与Matlab版本兼容问题

这是一个纯工程问题,但能卡住半天。新版Matlab默认按UTF-8读取脚本,但很多旧脚本是用GBK保存的,直接打开后中文注释全变乱码。我的处理办法是用记事本打开脚本,另存为UTF-8编码,再回到Matlab里运行。

另外,老代码里很多优化算法在不同版本中会被移除。比如quadprog的active-set算法在部分新版本中已经不再支持,如果论文代码是几年前的,运行时可能会报错。遇到这种情况,优先看函数文档确认当前版本支持的算法,再把options里的Algorithm改成'interior-point-convex'。

5.5 “复现”不等于“照搬公式”:参数对标最花时间

最容易被忽视的是基础数据。论文里每个微网的负荷曲线、光伏曲线、设备容量、天然气价格、上网电价,这些参数如果没有在正文里给全,复现结果一定对不上。

我的建议是先做一个“参数核验表”,把论文正文、附录、图表注释里能挖到的所有参数列出来,填到Matlab的params里。缺的参数根据常见取值范围补齐,然后在代码里做全局替换,方便后期调整。对标论文曲线时,先对比总体趋势,再比较典型时段的数值,如果趋势一致且数量级相同,基本就算复现成功。

6. 复现之外:这套策略的扩展方向与实用建议

6.1 从静态场景扩展到多时段动态调度

我一开始复现时只做的单一典型日24小时优化,后来发现论文可能还考虑了不同季节或者多日连续运行。多时段动态模型会让储能约束变成跨时段耦合约束,比如蓄热罐的容量递推方程不能简单拆到每个时刻独立求解,必须把整个调度周期的SOC作为子问题的一部分。

这个扩展在Matlab里并不难,只是决策变量的维度变大。只要把之前的变量向量复制 T 份,并把储能递推约束写进 A 矩阵,求解器仍然能处理。关键是共享价格也变成 T 维向量,迭代更新时要按每个时段分别计算共享偏差。

6.2 换用ADMM或KKT条件转换求解

如果论文里的纳什均衡最后是用KKT条件统一求解的,那就不是简单的迭代子问题,而是把所有微网的KKT条件拼成一个大型互补问题或者带均衡约束的优化问题。这类问题的Matlab实现门槛更高,需要专门处理互补条件和最优性条件。

另一种常见替代方案是ADMM。ADMM和纳什博弈虽然理论出发点不同,但在多主体分布式优化里经常能收到类似的收敛效果。如果你迭代算法一直不收敛,可以尝试把共享功率约束上的拉格朗日乘子用ADMM更新,公式更稳定,代码实现也不复杂。不过要小心,ADMM得到的是全局最优的分布式实现,而纳什均衡只保证单边最优性,两者结果不完全等价。

6.3 写代码前先做的一件小事

最后分享一个非常容易出现、但能省好几个小时的细节:在你开始写任何代码前,先给每个变量和参数建立命名规范,最好在代码文件头部注释里写清楚每个变量对应的物理单位。

比如Pchp,单位是MW;Hchp,单位是MWth;价格lambda_e,单位是元/MWh;目标函数里的成本系数a_chp,单位是元/MWh^2。很多对不上图的情况,不是模型错,而是单位差了一个1000倍系数。把单位写进变量名或者注释里,检查数质量会清楚很多。

我每次复现完一个复杂模型,最后都会把单个微网的子问题单独跑一遍,让它退化成独立运行模式,确认结果和常识一致后,再打开共享策略和博弈迭代。这个“从小到大”的调试顺序帮我避开了大量无效排错,也让我对论文模型的理解更深。希望这份复现思路也能让你的Matlab代码少走一点弯路。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询