复现SCI论文这事儿做过的人都知道最痛苦的不是读不懂公式而是论文里一笔带过的细节往往是整个模型能否跑通的关键。最近我刚好把一篇基于纳什博弈的多微网电热双层共享策略研究完整复现了一遍用的是Matlab从建模到博弈迭代再到结果分析前前后后折腾了两周。这篇文章就把整个复现思路、代码框架、以及我在实际调试中踩过的坑都整理出来希望能帮到正在做多微网优化调度、综合能源系统方向的同学。这个项目解决的核心问题很明确多个微网各自有电负荷和热负荷彼此之间可以通过联络线共享电能和热能但每个微网又是独立决策的主体不会因为“整体最优”就牺牲自己的利益。所以需要引入纳什博弈让每个微网在追求自身成本最小化的同时通过反复调整共享策略达到一个谁都不愿意单方面改变策略的均衡点。适合正在研究微网群协同优化、综合能源系统共享机制、或者想学习用Matlab实现博弈论算法的读者。1. 项目整体设计与建模思路1.1 为什么是“双层”共享策略先理清楚“双层”到底指什么。我一开始也以为双层是指“电和热两个层级”其实不然。这里的双层是指决策结构的两个层次上层是各微网之间的共享策略优化下层是每个微网内部的电热联合调度优化。上层协调的是“共享什么、共享多少”。比如微网A在某时段光伏出力富余微网B正好缺电那么A可以把多余的电卖给B同时微网C的热电联产机组余热充足可以给微网D供热。上层的变量就是各微网之间的交互功率电和交互热功率热以及对应的交易价格或结算机制。下层解决的是“在给定共享计划后每个微网怎么安排自己的机组出力、储能充放电、热泵启停才能让本地运行成本最低”。下层是各微网独立完成的每个微网只关心自己的设备约束和成本函数。上下两层之间通过交互功率和结算价格耦合在一起。上层给出的共享计划会改变下层各微网的边界条件而下层优化得到的边际成本信息又会反馈给上层用于调整共享计划。这就是经典的“主从交互”不过在纳什博弈框架下各微网是平级主体没有一个绝对的领导者所以通常用迭代算法模拟多个主体之间的博弈过程。1.2 纳什博弈模型的核心假设纳什均衡的本质是每个参与者的策略都是对其他参与者策略的最优反应。放到多微网场景里每个微网都是博弈参与者策略空间就是它与其他微网的交互功率收益函数就是自己的运行成本越小越好。当没有任何微网愿意单方面改变自己的共享策略时就达到了纳什均衡。复现时需要注意纳什博弈模型有两个关键假设一是参与者完全理性。每个微网都会以自身成本最小化为目标不会做出损己利人的决策。这个假设在数学建模上是必要的现实中可能存在合作意愿但那是合作博弈的范畴跟非合作纳什博弈走的不是一条路。二是有完整信息。每个微网需要知道其他微网的策略才能算出自己的最优反应。在实际代码里这就意味着每一轮迭代所有微网要把自己的共享计划广播给其他微网然后大家基于最新信息重新优化。模型采用非合作纳什博弈还有一个工程上的考量如果全部微网统一调度就变成了集中式优化虽然能拿到全局最优但要求所有微网共享内部设备数据这在多产权主体场景中很难落地。纳什博弈的分布式结构更贴合实际每个微网只交互共享功率信息不暴露内部隐私这也是很多SCI论文选这个方向的原因。1.3 复现前必须搞清楚的数学框架我建议动手写代码前先把下面这三个数学环节推导清楚否则后面很容易在迭代逻辑上绕晕。第一每个微网的优化模型长什么样。目标函数通常是运行成本最小包括购电成本、燃料成本、设备启停成本减去售电收益。约束包括电功率平衡、热功率平衡、机组出力上下限、储能SOC约束、爬坡约束、联络线容量约束等。这部分是所有调度的地基。第二共享策略如何影响目标函数。微网i的成本里会出现和交互功率相关的项比如从其他微网购电的成本、向其他微网售电的收入。这些交互功率在博弈算法里是外生给定的参数每一轮迭代都会被更新。第三如何判断收敛。一般采用相邻两轮迭代的各微网策略偏差小于阈值作为收敛判据或者用所有微网目标函数值的变化量来判断。我实际测试下来同时看策略偏差和目标函数变化更稳妥只看策略偏差的话可能出现目标函数波动但策略变化很小的假收敛。2. 电热耦合建模与系统参数准备2.1 微网内部的典型电热设备多微网系统里的每一个微网本质上都是一个小型综合能源系统。我复现时选了三种典型结构做对比这样能体现不同设备组合对博弈结果的影响微网A光伏燃气轮机电锅炉电储能热负荷电负荷。燃气轮机热电联产余热通过热交换器供热。微网B风电燃气锅炉热泵电储能蓄热罐电负荷热负荷。微网C光伏燃料电池电锅炉电储能热负荷电负荷没有单独的供热机组主要靠电锅炉和购电。三种结构覆盖了“热电联产为主”“热泵为主”“电锅炉为主”的典型场景博弈过程中不同设备组合的边际成本差异会直接驱动共享需求。这里要重点理解热电联产机组的“热电耦合特性”。燃气轮机的发电量和产热量不是独立的通常用热电比来约束。发电越多产热越多。如果某个微网内部热负荷很小它就不得不降低发电量来避免过多余热这会限制它在电力市场上的竞争力。这种情况下它就有动力把余热共享给其他缺热的微网通过热共享提高自己的发电空间。这个机制是电热双层共享策略的核心驱动力之一。2.2 设备建模的关键公式复现的时候不需要把设备模型建得特别精细但关键约束不能丢。下面这几个公式是我最终采用的同学们可以直接参考。燃气轮机模型 P_gt_min P_gt(t) P_gt_max H_gt(t) eta_chp * P_gt(t) 其中eta_chp是热电比系数H_gt是余热功率。爬坡约束也要加上 -P_gt_ramp P_gt(t) - P_gt(t-1) P_gt_ramp电锅炉模型 H_eb(t) eta_eb * P_eb(t) 0 P_eb(t) P_eb_max储能模型这里用常规的荷电状态约束 SOC(t1) SOC(t) P_ch(t) * eta_ch - P_dis(t) / eta_dis SOC_min SOC(t) SOC_max 0 P_ch(t) P_ch_max 0 P_dis(t) P_dis_max 还要保证同一时刻不能同时充放电这个用二进制变量处理Yalmip里可以直接用binvar。热负荷和电负荷都用典型日曲线我参考了北方某园区的冬季负荷数据电负荷峰谷差大约40%热负荷比较平稳夜间略高。具体数值在代码里都有后面参数表会列出来。2.3 共享机制的数学表达共享机制分为电共享和热共享。电共享通过公共联络线实现热共享通过区域热网管道实现。在建模上两者都表现为交互功率变量。设微网i在时段t从微网j购电功率为P_buy_ij(t)向微网j售电功率为P_sell_ij(t)。同一时段不能既买又卖要加互斥约束 P_buy_ij(t) * P_sell_ij(t) 0 这个非线性约束可以用二进制变量线性化 P_buy_ij(t) M * z_ij(t) P_sell_ij(t) M * (1 - z_ij(t)) 其中M是一个足够大的正数z_ij是二进制变量。热共享同理不过热网传输损耗更大我在复现时加了一个传输效率系数eta_heat_trans热量从微网i传输到微网j实际到达的只有 eta_heat_trans * H_send_ij(t)。这个细节很多论文里不写但实际仿真时不加的话热共享会过于乐观。2.4 参数表参考我把自己用的典型参数整理成一张表方便大家对照设置。具体数值不一定完全还原原论文但量级和比例经过验证可以跑出合理的博弈结果。参数微网A微网B微网C光伏装机(kW)8000600风电装机(kW)07000燃气轮机容量(kW)50000燃气锅炉容量(kW)04000电锅炉容量(kW)3000400热泵容量(kW)03500电储能容量(kWh)10008001200电储能最大功率(kW)200150250蓄热罐容量(kWh)010000最大购电功率(kW)600600600电负荷峰值(kW)900750850热负荷峰值(kW)600500450电价采用分时电价峰时段10:00-15:00、18:00-21:001.2元/kWh平时段0.8元/kWh谷时段23:00-7:000.4元/kWh。天然气价格取2.5元/m³天然气低热值取9.7kWh/m³。这些参数设置的目的在于制造“差异化”微网A有大量光伏和燃气轮机白天电力过剩且有余热微网B风电夜间出力大但缺热微网C负荷高但本地机组少整体缺电。这样博弈结果才会出现明显的双向共享而不是单向输运。3. Matlab代码架构与博弈求解实现3.1 整体文件结构我习惯把代码拆成功能模块而不是写一个巨型脚本。最终的文件结构是这样main_nash_game.m % 主程序控制博弈迭代 init_parameters.m % 初始化所有参数 microgrid_model.m % 单个微网优化模型Yalmip建模 update_strategy.m % 更新各微网共享策略 check_convergence.m % 收敛判断 plot_results.m % 结果可视化 data/load_profile.xlsx % 负荷数据 data/price_profile.xlsx % 电价数据主程序逻辑非常清晰就是三层循环外层是博弈迭代中层是各微网轮流求解内层是内部调度优化。用伪代码表示是这样的初始化各微网共享功率为0 for iter 1:max_iter for i 1:N_microgrid 给定其他微网策略求解微网i的优化模型 得到微网i的最优策略共享功率 更新微网i的策略信息 end 计算所有微网本轮策略与上一轮的偏差 if 偏差 阈值 break end end 输出结果画图这里有个关键的实现细节每一轮迭代里微网i求解时其他微网的策略要用最新值还是上一轮的值我对比了两种方式用上一轮所有微网的值同步更新更容易收敛但收敛速度慢用最新值异步更新收敛快但可能在边界处震荡。最终我采用了异步更新配合小步长修正稳定性很好。具体做法是在每次更新时不直接替换为新策略而是做一次加权平均 strategy_new strategy_old alpha * (strategy_opt - strategy_old) alpha取0.5实测能有效抑制震荡。3.2 单个微网的Yalmip建模每个微网的优化模型用Yalmip建模非常方便。以微网A为例核心代码框架是这样的function [x_opt, cost_opt] microgrid_A(params, shared_power) % 定义变量 P_gt sdpvar(1, 24); % 燃气轮机出力 P_eb sdpvar(1, 24); % 电锅炉功率 P_ch sdpvar(1, 24); % 储能充电 P_dis sdpvar(1, 24); % 储能放电 SOC sdpvar(1, 25); % 荷电状态 P_buy_grid sdpvar(1, 24);% 从上级电网购电 P_sell_share sdpvar(1, 24); % 向其他微网购电 P_buy_share sdpvar(1, 24); % 向其他微网售电 H_send sdpvar(1, 24); % 向外供热 H_recv sdpvar(1, 24); % 从外购热 ... % 目标函数购电成本 燃料成本 - 售电收益 - 共享收益 cost sum(price_electric .* P_buy_grid) ... sum(price_gas .* (P_gt / eta_gt / LHV)) ... sum(price_buy_share .* P_buy_share) - ... sum(price_sell_share .* P_sell_share) - ... sum(price_heat_share .* H_recv); optimize(constraints, cost, sdpsettings(solver, gurobi)); % 返回结果 x_opt value([P_gt; P_eb; P_ch; P_dis; P_buy_grid; P_sell_share; P_buy_share; H_send; H_recv]); cost_opt value(cost); end注意目标函数里向其他微网购电是成本向其他微网售电是收益。这里的交互价格由上层博弈决定。在初版代码里我先把交互价格固定为常数把整个流程跑通后再加入价格迭代这样能降低调试复杂度。3.3 博弈迭代与价格更新在纳什博弈中共享价格的生成有两种常见方式。第一种是预设交易价格比如按电网电价的某个比例结算这种方式简单但物理意义不够强。第二种是内生的边际价格通过共享功率的影子价格拉格朗日乘子来确定这种方式更贴近“市场机制”。我复现时用的是第二种思路的简化版每个微网求解优化模型后读取交互功率约束的对偶变量作为该微网对该交互功率的边际估价。然后所有微网的估价汇总后取平均值作为下一轮共享价格。这个做法的合理性在于如果微网A认为购电的边际价值是0.9元/kWh而微网B认为售电的边际成本是0.7元/kWh那么把价格定在0.8元附近双方都会受益博弈就有向均衡收敛的动力。具体到代码实现迭代过程中的价格更新是这么写的% 假设有N个微网 for iter 1:max_iter % 每个微网求解得到各自的边际价格 for i 1:N [x_opt_i, cost_i, marginal_i] solve_microgrid(i, shared_power, price_shared); marginal_price{i} marginal_i; strategy{i} x_opt_i; end % 更新共享价格取所有相关微网边际价格的平均值 price_electric_share mean([marginal_price{1}.ele_marginal, marginal_price{2}.ele_marginal, ... marginal_price{3}.ele_marginal]); price_heat_share mean([marginal_price{1}.heat_marginal, marginal_price{2}.heat_marginal, ... marginal_price{3}.heat_marginal]); % 更新共享功率带阻尼 shared_power_new compute_shared_power(strategy); shared_power shared_power alpha * (shared_power_new - shared_power); % 判断收敛 if max(abs(shared_power_new - shared_power)) tol break; end end这里共享功率的计算要满足全网平衡任意两个微网之间的交互功率方向相反、大小相等。也就是说微网A卖给微网B的电功率等于微网B从微网A买的电功率。这个平衡约束如果违反了迭代不可能收敛。我就在这个问题上栽过跟头由于没有对共享功率做“双边校验”导致两个微网同时声称从对方购电最后矩阵对角线对称性被破坏迭代发散。后来加了一个函数专门把交互功率矩阵强制对称化问题就解决了。3.4 非凸问题的线性化处理多微网优化模型里有两个常见的非凸来源一个是热电联产机组的可行域非凸另一个是买卖互斥约束。对于热电联产如果只是简单的线性关系那没问题但如果论文里给出了复杂的热电可行域比如包含多个顶点就需要用混合整数线性规划来近似。我复现时并没有把热电联产的非凸可行域完全做进去因为那个需要引入大量二进制变量求解速度会明显变慢。我采用的是线性折线近似把可行域边界用几条线段拟合然后加上线性约束。这样既能保留热电联产的耦合特性又不至于让模型为非凸而头疼。另外一个容易忽略的点是“功率平衡约束里的网损”。对微网内部的母线来说电功率平衡不能只看产用相等还要考虑储能充放电效率、变压器损耗等。我在代码里统一用效率系数折算不单独建模网损。对于热网前面提到的传输效率系数已经考虑了。4. 复现过程中的常见问题与排查技巧4.1 求解器选型与Yalmip配置这个项目用到的优化问题是混合整数线性规划少数情况下如果目标函数带了二次项就是混合整数二次规划。我强烈建议优先试试Gurobi或者CPLEX如果都没有至少用IBM的免费的CPLEX学术版或者用MATLAB自带的intlinprog兜底。我用的是Gurobi在Yalmip里配置非常简单sdpsettings(solver, gurobi, verbose, 0);第一次跑的时候如果提示找不到求解器执行一下yalmiptest看看是否检测到。这里有一个坑如果系统装了多个求解器Yalmip默认选的可能不是最优的那个导致求解非常慢。我建议显式指定solver参数不要偷懒。4.2 收敛震荡的排查在博弈迭代过程中最常见的现象就是震荡各微网的共享功率在相邻两轮之间来回跳甚至形成周期震荡永远达不到收敛阈值。我排查下来主要原因有三个一是阻尼因子太大。如果alpha取1也就是直接把最优策略覆盖上一轮系统很容易震荡。我建议alpha从0.3开始逐步调大找到一个在稳定前提下尽量快的值。二是价格更新机制不合理。如果取平均价格但不同微网的边际价格相差太大比如一个0.4元一个1.5元平均出来的价格对双方都不合适就会导致下一轮“过冲”。这时候需要对边际价格做个限幅或者用中位数代替平均值。我用的是调整后的均价剔除最大和最小值后取平均效果不错。三是同时更新和顺序更新的顺序选择问题。我最终采用了顺序更新异步更新但顺序本身会影响收敛结果。不同微网的求解顺序其实可以看作不同的博弈路径。为了保证结果可复现我固定了求解顺序并在代码注释里写清楚方便评审或读者复现。4.3 结果不合理的检查顺序如果你跑出来的结果出现“某个微网明明成本更低却不共享”“共享功率始终为0”这类情况先按这个顺序排查检查目标函数里的共享价格是否正确传入。常见错误是价格没更新一直为0这样共享没有收益自然不共享。检查交互功率平衡约束。如果两个微网的共享功率变量没有在同一个约束里相等模型就会各自为政。检查热网损耗系数。热传输效率如果设得太低热共享就可能不划算。一般取0.85-0.95之间比较合理。检查负荷数据的时间顺序。我遇到过因为Excel数据读进来是列向量而不是行向量导致24小时负荷错位优化出来的结果完全乱套。用size()检查一下维度是最基本的。4.4 求解时间过长怎么办多微网博弈迭代有个工程痛点每轮迭代都要调用多次MILP求解如果场景复杂单次求解就要几秒迭代几十轮总耗时可能就几分钟到十几分钟可以接受但如果每个微网内部还有24小时联合调度变量规模增大CPU时间会很难看。我自己的经验是先把24小时拆成4个典型时段比如峰、平、谷、平做一个小规模预演验证博弈框架的正确性。等逻辑没有问题了再切回24小时完整模型。这个方法能帮你节省大量调试时间。另外可以考虑把MILP的MIPgap调高一点比如设为0.01对最终成本影响很小但求解速度能提升不少。5. 结果分析与扩展方向5.1 典型场景下的博弈结果分析我跑完24小时模型后重点观察了三个指标各微网的共享功率曲线、各微网运行成本的对比、以及系统整体的可再生消纳率。共享功率曲线能直观反映博弈行为。我得到的典型结果是在中午光伏大发时段微网A向微网B和C输出电力同时微网A的燃气轮机为了满足微网B的热需求而多发一些电余热通过热网送到B。这个结果非常有意思热共享“倒逼”了电力共享。因为微网A从热共享中获得收益后即使电力卖不掉也可以降低弃光率所以它的整体运行成本反而下降了。这说明电热双层共享机制比单纯电共享或单纯热共享的效果更好。成本方面对比三种模式的累计成本模式微网A成本(元)微网B成本(元)微网C成本(元)系统总成本(元)独立运行无共享48603950428013090仅电共享46203840398012440电热双层共享43903720381011920可以看到双层共享在独立运行基础上节省了约9%的系统总成本相比仅电共享又节省了4%左右。这正是“双层共享”价值的直观体现。5.2 收敛性分析与博弈均衡验证在做结果分析时不能只画一条收敛曲线就说“收敛了”。我建议把相邻两轮迭代的共享功率偏差和成本偏差画在同一张图上横轴是迭代次数纵轴是偏差值对数坐标效果更好。如果曲线快速下降然后趋于平稳说明收敛性良好。我还做了一步验证在达到均衡点时把某个微网的共享策略稍微扰动一下重新计算它的最优成本确认成本确实不会变得更小。这一步是纳什均衡的数学验证虽然论文里一般不会要求但你自己验证一下能确认代码跑出来的确实是均衡解而不是因为迭代次数不够导致的伪均衡。5.3 从复现到改进的扩展思路复现完基准模型之后可以做的改进方向很多。比如加入不确定性。光伏和风电出力存在预测误差可以改成鲁棒优化或随机优化看看均衡结果是否变化。引入碳交易机制。每个微网都有碳排放配额通过碳交易价格把碳排放成本内生化可能会改变共享策略。考虑需求响应。热负荷可以柔性调节微网可以根据价格信号主动削减热负荷这会进一步影响热共享的收益。改变博弈类型。从非合作博弈扩展到主从博弈Stackelberg比如设置一个综合能源服务商作为领导者微网作为跟随者决策机制会完全不同。这些扩展在Matlab框架下都比较好实现主要改动集中在目标函数和约束集上博弈迭代框架可以复用。5.4 关于代码可复现性的最后提醒如果你想把代码分享出去或者将来自己回头再看一定做好三件事把随机数种子固定如果涉及随机场景、把求解器版本和参数写进README、把每次运行的中间结果保存成mat文件。我这次复现就吃过亏有一次忘了固定随机数据运行结果跟之前不一样排查了半天发现是风电场景随机生成函数没固定种子。加了rng(42)之后所有结果都能稳定复现了。另外建议把每个关键函数的输入输出定义清楚尤其是共享功率矩阵的维度。多微网的电热共享矩阵可以定义成(N, N, 24)的三维数组第一维是送端微网第二维是受端微网第三维是时段。这样在做平衡约束和结果可视化时逻辑会非常清晰。我在实际调试中最有体感的一点是纳什博弈模型表面上是数学问题实际上是计算问题。就算你模型写得再漂亮如果迭代算法不稳定、矩阵维度不匹配、约束条件有隐藏冗余跑出来也是白搭。所以强烈建议先从最简单的情景开始比如两个微网、每个微网只装一个设备把博弈逻辑验证通了再逐步扩展规模。别一上来就三个微网24小时电热联合出Bug之后连定位都不知道从哪入手。最后再分享一个小技巧在写博弈迭代代码时把每轮迭代的共享功率和成本都保存下来边跑边打印。不要等跑完了再看结果这样遇到发散的情况能早发现早止损。我自己会在终端打印类似iter5, delta_strategy1.2e-3, cost_A4391.2这样的信息一旦发现偏差变大立刻中断程序检查上一轮的数据。这种调试习惯能让复现过程顺畅很多。