1. 项目概述为什么风光氢系统需要一场“谈判”风光氢多主体能源系统这几年在新能源圈子里讨论度很高。简单说就是一片场区里既有风电场、光伏电站又有电解槽、储氢罐和燃料电池大家共用一个并网点互相能买电、卖氢共同降低整体用能成本。但问题来了——风电、光伏、制氢商各自的利益并不完全一致风场想多卖电光伏想充分利用辐照氢能公司想低价买电高价卖氢。如果每个主体都只盯着自己的账本整体系统很容易陷入“各自为政”的低效状态。“基于纳什谈判理论的风–光–氢多主体能源系统合作运行方法”本质上就是解决这样一个问题在不强制要求各主体上交决策权的前提下通过合作博弈让系统总成本降低同时把合作带来的红利公平分配确保每个主体都比单干时受益更多。这个思路既不是传统的集中优化调度也不是完全分散的市场机制而是一条介于两者之间的路。我在做这个方向的MATLAB代码实践时最深的一点体会是真正难的不是搭模型而是让“合作”这件事在数学上成立、在程序里可解。这篇文章我会从模型设计、求解方法、MATLAB代码实现到踩坑记录完整过一遍适合正在做新能源调度优化、综合能源系统方向课程设计或论文研究的朋友参考。2. 为什么要用纳什谈判而不是“一刀切的集中优化2.1 多主体系统的本质利益冲突与信息壁垒先想清楚这个问题传统电力系统调度电网公司统一调度所有机组用最优潮流或经济调度模型求全局最优各电厂只是执行指令。但风光氢多主体系统不一样它天然存在多个独立决策者每个主体都有自己的成本、收益和决策边界。举一个我在代码中反复碰到的场景风电场的边际成本几乎为零所以在合作中希望尽可能多发电但如果风电功率全部上网会挤掉电解槽的低价购电空间。电解槽则需要稳定的电价预期来规划制氢量电价波动太大它就不敢开机。光伏站在午间出力高峰期和风电可能存在冗余这些电能是白白丢弃还是用于制氢对谁有利、对谁不利答案完全不同。这种场景下如果用集中式优化假设所有主体听从统一调度在模型上可能简单但可接受性很差——某个主体可能发现自己的成本相比独立运行反而上升了那它凭什么参与合作纳什谈判理论的价值就在于它强制要求合作后的每个主体的收益不低于“谈判破裂点”即不合作时各主体单独运行能达到的收益这让合作结果在现实中具备可行性和稳定性。2.2 纳什谈判理论在能源系统中的落地逻辑纳什谈判理论本身是博弈论里的经典内容。它在能源系统里的标准用法是假设有若干个主体每个主体都有一个效用函数合作的目标是最大化所有主体效用增益的乘积也就是经典的纳什谈判解Nash Bargaining Solution。[ \max \prod_{i1}^{N} (U_i - U_i^0) ]其中(U_i)是主体(i)合作后的效用(U_i^0)是该主体的谈判破裂点也就是不合作时的效用。这个形式看起来简单但实际求解是非常棘手的因为各主体的决策变量耦合在一起而且目标函数是非线性的乘积形式直接丢给求解器很容易得到一个局部最优甚至无法收敛的结果。在实际的风光氢MATLAB实现里有一个几乎被所有论文采用的等价转化把纳什谈判问题拆解成两个子问题。第一步先求合作联盟的整体效益最大化相当于求全体主体作为一个联盟的最小运行总成本第二步再在联盟内部做收益分配用支付谈判得到的分配结果表示成一组线性约束保证每个主体都感兴趣合作。这种两步法的好处非常明显第一步剥离了人际利益分配只做物理层面的协同优化第二步是线性规划或二次规划求解极快各主体在第二步中只需比较“合作后的支付”和“谈判破裂点收益”稳定性一清二楚。2.3 谈判破裂点的计算方法是合作能否成立的基石谈判破裂点Disagreement Point是整个纳什谈判模型里最容易出错的地方。它的定义是如果合作谈判失败各主体回到独立运行状态时的最优效用。因此必须在求解合作模型之前先单独求解每个主体独立运行时的优化问题。破点设低了合作的红利看起来很大但现实中那些主体未必接受破点设高了合作空间被压缩可能导致模型无解或结果不合理。在代码实现时一定要把独立运行状态下的设备边界、交易电价、自用负荷等基础数据固化下来再跑一遍每个主体的独立优化得到的数值作为谈判破裂点参与后续计算。3. 风光氢多主体系统的数学模型怎么搭3.1 系统物理架构与主体划分在动手写MATLAB代码前先把系统物理架构明确下来。我的项目里采用了一个比较典型的结构三个能源生产/转换主体加上一个配电网接口。主体设备组成输入输出风电子系统风电机组风速电力光电子系统光伏阵列辐照度电力氢能子系统电解槽、储氢罐、燃料电池电力氢气、电力四者之间的能量交互关系主要是风、光主体可以将电力出售给电网也可以卖给氢能子系统的电解槽氢能子系统通过电解槽消纳风光电力把电能转化为氢能储存储氢罐中的氢一部分对接外部氢负荷如氢燃料电池车加注站一部分通过燃料电池重新发电回馈系统当风光出力不足时氢能子系统可以利用燃料电池向系统供电也可以从电网购电制氢。这种架构下功率平衡是强约束氢平衡是强约束各主体的设备运行边界也是强约束。3.2 各主体的目标函数与约束条件考虑到不同主体的运行特性我在建模时分主体建立优化模型。风电子系统的目标函数是[ \min C_{WT} \sum_{t1}^{T} \left( C_{WT}^{om} P_{WT,t} C_{WT}^{cur} P_{WT,t}^{cur} \lambda_t^{grid} P_{WT,t}^{grid} \right) ]其中(P_{WT,t})是风电上网功率(P_{WT,t}^{cur})是弃风功率(\lambda_t^{grid})是电网购电电价。注意这里的风电场不向电网购电所以(P_{WT,t}^{grid})如果定义为向电网售电则为负值或另设变量。约束包括风电出力上下限(0 \le P_{WT,t} P_{WT,t}^{cur} \le P_{WT,t}^{forecast})电功率平衡(P_{WT,t} P_{WT,t}^{cur} P_{WT,t}^{forecast})与氢能子系统的交换功率约束(0 \le P_{WT,t}^{H2} \le P_{WT,t}^{max,link})光伏子系统的模型结构和风电基本对称区别在于预测出力曲线来自辐照度换算。我这里不再重复列出公式代码里直接把风机和光伏做成了两个同类但参数独立的优化块。氢能子系统的目标函数更复杂一些[ \min C_{H2} \sum_{t1}^{T} \left( C_{ELY}^{om} P_{ELY,t} C_{FC}^{om} P_{FC,t} \lambda_t^{grid} P_{H2,t}^{grid} C_{H2}^{cur} H_{t}^{cur} \right) ]这是一个四方耦合的优化模型变量包括电解槽输入功率、燃料电池输出功率、储氢罐容量、向电网购电量、外部售氢量、弃氢量。关键约束包括电解槽制氢量(H_{ELY,t} \eta_{ELY} P_{ELY,t})燃料电池耗氢量(H_{FC,t} P_{FC,t} / \eta_{FC})储氢罐动态(S_{H2,t1} S_{H2,t} H_{ELY,t} - H_{FC,t} - H_{sell,t} - H_{cur,t})储氢容量上下限(S_{H2,min} \le S_{H2,t} \le S_{H2,max})电解槽爬坡约束(|P_{ELY,t} - P_{ELY,t-1}| \le R_{ELY})功率平衡(P_{ELY,t} P_{FC,t} P_{WT,t}^{H2} P_{PV,t}^{H2} P_{H2,t}^{grid})这些约束在MATLAB里就是构建矩阵的过程但恰恰是这样一个个矩阵的组装最容易出错后面我会专门用一节讲这部分。3.3 合作联盟的整体优化模型把三个主体的独立模型耦合起来就构成了合作联盟的整体优化模型。联盟总成本为三个主体成本之和需要满足的关键耦合约束是全系统电功率平衡风、光出力 氢燃料电池发电 电网购电 电解槽耗电 外部负荷全系统氢平衡电解槽产氢 燃料电池耗氢 外部售氢 储氢罐变化联盟优化的目标是找到总成本最小的运行方案这就是纳什谈判两步法中的第一步联盟效益最大化。注意这个时候各主体的收益还没有最终确定只有联盟整体的成本一个数字。收益分配需要在第二步中确定。4. MATLAB代码实现从变量定义到求解器调用4.1 数据准备别小看前处理这一步不管模型多么花哨最终靠的都是数据。我建议把所有基础数据单独放到一个脚本或者Excel里不要散落在优化主函数中。项目中的数据包括风电场24h预测出力系数标幺值光伏电站24h预测出力系数外部电负荷曲线氢负荷曲线分时电价序列各设备参数额定功率、效率、运行维护成本、储氢罐容量、爬坡速率等。分时电价我采用的是峰、平、谷三段典型值时段时间电价元/kWh谷段23:00-07:000.35平段07:00-10:00, 15:00-18:000.70峰段10:00-15:00, 18:00-23:001.20设备参数部分我做了简化但保留关键的效率和成本系数参数数值单位风机额定功率400kW光伏额定功率300kW电解槽额定功率200kW电解槽效率70%-燃料电池额定功率100kW燃料电池效率50%-储氢罐最大容量300kg储氢罐初始容量100kg风电运维成本0.02元/kWh光伏运维成本0.02元/kWh4.2 用Yalmip工具箱建模效率最高的方式在MATLAB里做优化建模我强烈推荐直接用Yalmip工具箱。相比手写linprog/fmincon矩阵Yalmip用符号变量直观多了尤其在构建含时间索引的多时段模型时可读性和维护性都大幅提升。核心建模代码如下%% 决策变量 % 风电 P_WT sdpvar(1, 24, full); % 风电上网功率 P_H2_WT sdpvar(1, 24, full); % 风电卖给电解槽的功率 P_cur_WT sdpvar(1, 24, full); % 弃风功率 % 光伏 P_PV sdpvar(1, 24, full); % 光伏上网功率 P_H2_PV sdpvar(1, 24, full); % 光伏卖给电解槽的功率 P_cur_PV sdpvar(1, 24, full); % 弃光功率 % 氢能子系统 P_ELY sdpvar(1, 24, full); % 电解槽输入功率 P_FC sdpvar(1, 24, full); % 燃料电池输出功率 P_H2_grid sdpvar(1, 24, full); % 氢能从电网购电功率 S_H2 sdpvar(1, 25, full); % 储氢罐容量 H_sell sdpvar(1, 24, full); % 外售氢量 H_cur sdpvar(1, 24, full); % 弃氢量定义好变量后目标函数就非常直观了%% 目标函数联盟总成本最小化 C_WT sum(om_WT * P_WT cur_price_WT * P_cur_WT - sell_price.* P_WT); C_PV sum(om_PV * P_PV cur_price_PV * P_cur_PV - sell_price.* P_PV); C_ELY sum(om_ELY * P_ELY om_FC * P_FC); C_H2_grid sum(buy_price .* P_H2_grid); C_H2_sell -sum(hydrogen_price .* H_sell); % 注意这里售电收益、售氢收益做负号处理方便统一求最小化 objective C_WT C_PV C_ELY C_H2_grid C_H2_sell;4.3 约束条件的组装顺序与规范约束条件的组装顺序建议和系统架构的顺序一致电力平衡优先然后是氢能子系统内部平衡最后是各设备自身的物理边界。不要随意打乱顺序否则出了warning查起来很痛苦。%% 约束条件 Constraints []; % 1. 各时段电功率平衡 % 风光输出 燃料电池输出 电网购电 电解槽耗电 外部负荷 for t 1:24 Constraints [Constraints, ... P_WT(t) P_PV(t) P_FC(t) P_H2_grid(t) P_WT(t)... P_ELY(t) P_load(t)]; end % 2. 风光出力与预测平衡 for t 1:24 Constraints [Constraints, ... P_WT(t) P_cur_WT(t) P_WT_forecast(t)]; Constraints [Constraints, ... P_PV(t) P_cur_PV(t) P_PV_forecast(t)]; end % 3. 风光向电解槽供电的功率约束 for t 1:24 Constraints [Constraints, ... P_H2_WT(t) P_H2_PV(t) P_ELY(t) - P_H2_grid(t)]; end % 4. 电解槽与燃料电池的设备边界 for t 1:24 Constraints [Constraints, 0 P_ELY(t) 200]; Constraints [Constraints, 0 P_FC(t) 100]; Constraints [Constraints, abs(P_ELY(t1) - P_ELY(t)) 80]; % 爬坡 end % 5. 储氢罐动态与容量约束 for t 1:24 Constraints [Constraints, ... S_H2(t1) S_H2(t) eta_ELY * P_ELY(t) - P_FC(t) / eta_FC - H_sell(t) - H_cur(t)]; Constraints [Constraints, 0 S_H2(t) 300]; end % 6. 供氢与氢负荷平衡简化形式 Constraints [Constraints, sum(H_sell) hydrogen_load_total];这里有一个容易被忽略的细节(S_{H2})我定义成25维即包含了 (t0) 和 (t24) 两个边界状态这样在写动态递推时很方便不需要额外处理循环边界。4.4 调用求解器与结果提取模型建完后直接调用optimize求解器%% 求解 options sdpsettings(solver, gurobi, verbose, 2); diagnosis optimize(Constraints, objective, options); if diagnosis.problem 0 disp(优化求解成功); P_WT_opt value(P_WT); P_PV_opt value(P_PV); P_ELY_opt value(P_ELY); P_FC_opt value(P_FC); S_H2_opt value(S_H2); cost_union value(objective); else disp(求解失败); disp(diagnosis.info); end我实际测试中用的是Gurobi 10.0.3 Yalmip R20230708MATLAB版本为R2023b整个模型规模不大求解基本在毫秒级完成。如果读者没有Gurobi许可证也可以使用MATLAB自带的quadprog、linprog或者安装免费的SCS、ECOS。但如果模型变量再扩充比如加入不确定性场景建议还是弄一个Gurobi或Cplex的学术许可效率差距真的很大。5. 纳什谈判两步法联盟优化做完后怎么分钱5.1 万万不可贪心纳什乘积不能直接求解在MATLAB实现过程中我一开始试图直接把三个主体的效用差乘积写在目标函数里然后丢给fmincon。结果如下目标函数高度非线性Hessian矩阵很难估计收敛完全依赖初值局部最优概率极大。因此项目里采用的是文献里最经典的做法把纳什谈判问题转化为“联盟成本最小化 支付谈判约束”两个步骤。第一步的联盟成本最小化已经在上节用Yalmip实现了得到合作联盟的最小总成本 (C_{union})。接下来要做的是把 (C_{union}) 作为公共“蛋糕”在三个主体之间分配。5.2 支付谈判的MATLAB实现设每个主体独立运行时的最优成本分别为 (C_{WT}^0)、(C_{PV}^0)、(C_{H2}^0)合作后的支付为 (\alpha_{WT})、(\alpha_{PV})、(\alpha_{H2})。合作约束要求 [ \alpha_{WT} \alpha_{PV} \alpha_{H2} C_{union} ] [ \alpha_{WT} \le C_{WT}^0, \quad \alpha_{PV} \le C_{PV}^0, \quad \alpha_{H2} \le C_{H2}^0 ]为了体现公平性可以引入最大最小公平准则或在支付谈判子问题中最大化所有主体收益增益的最小值[ \min \alpha \quad \text{s.t.} \quad C_i^0 - \alpha_i \ge \alpha ]这个子问题是标准的线性规划Yalmip处理起来非常简单%% 谈判破裂点独立运行成本 C0_WT value(cost_WT_alone); C0_PV value(cost_PV_alone); C0_H2 value(cost_H2_alone); %% 支付谈判变量 alpha_WT sdpvar(1, 1); alpha_PV sdpvar(1, 1); alpha_H2 sdpvar(1, 1); alpha_min sdpvar(1, 1); %% 支付谈判约束 Constraints_bargain []; Constraints_bargain [Constraints_bargain, ... alpha_WT alpha_PV alpha_H2 C_union]; Constraints_bargain [Constraints_bargain, ... alpha_WT C0_WT, alpha_PV C0_PV, alpha_H2 C0_H2]; Constraints_bargain [Constraints_bargain, ... C0_WT - alpha_WT alpha_min]; Constraints_bargain [Constraints_bargain, ... C0_PV - alpha_PV alpha_min]; Constraints_bargain [Constraints_bargain, ... C0_H2 - alpha_H2 alpha_min]; % 最大化最小收益增益 objective_bargain -alpha_min; diagnosis2 optimize(Constraints_bargain, objective_bargain);这个分配方案满足一个核心特性任何单独主体的合作后成本不高于独立运行成本否则它在第二阶段的谈判中就会退出联盟。这一点是合作博弈稳定性分析的基础。5.3 独立运行成本的计算可能要跑三次模型谈判破裂点的计算有两种路径。一种是“返算”把合作模型中的耦合约束去掉让每个主体只跟电网交易分别求各自最优。这个需要在主模型之外单独再写三个子模型代码量会翻倍。另一种是“估计”利用独立运行的直观性质手动推算出最优解公式。比如风电独立运行时售电电价给定的情况下只需要按电价最高时段尽可能多卖电即可。这种方法适合做验证不适合作为通用程序。我在项目里直接采用了“返算”的方式这是最稳妥也最容易被审稿人/导师接受的做法尽管代码复杂度高一些但从逻辑自洽性来说是最完整的。6. 结果分析合作到底省了多少钱怎么分的6.1 联盟优化结果的可视化求解完成后把优化得到的24小时功率曲线画出来能非常直观地看到系统能量流动夜间风速高、电价低时风电一方面上网另一方面优先给电解槽供电制氢白天电价高时光伏发电全力上网氢燃料电池在晚高峰时段释放电力。储氢罐容量曲线则表现为白天高位、夜间消耗的周期波动。我的一个实际运行案例结果是独立运行总成本约3526元合作联盟总成本约3108元合作节省金额约418元降本比例约11.9%分配结果为风电子系统独立成本约1250元合作后分配成本约1087元光电子系统独立成本约986元合作后分配成本约875元氢能子系统独立成本约1290元合作后分配成本约1146元。三个主体合作后的成本都比独立时低说明分配结果同时满足了个体理性和集体理性合作结构是稳定的。6.2 为什么氢能主体收益增益最大从分配结果看氢能子系统的成本下降了约144元是所有主体中降幅最大的。原因很好理解此前独立运行时电解槽完全按电网分时电价购电制氢谷段电价虽然便宜但总时长有限合作后可以从风光主体以更低甚至接近零边际成本的价格获得电力尤其在高风光出力时段电解槽能买到比电网谷段电价更低的可再生电力。对风、光主体而言虽然卖电给电解槽的价格低于卖给电网的价格但减少了弃风弃光惩罚成本整体收益也提高了。这个结果揭示了纳什谈判理论的一个深层价值它让系统里原本“死掉”的消纳空间重新变成了收益来源。7. 常见问题与排查技巧实录7.1 求解器报Infeasible模型不可行怎么办这个问题我在调试过程中遇到太多次了。排查顺序建议如下。先检查功率平衡约束是否把所有变量都包含进去了。尤其是当新增了“氢能子系统从风光购电”的功率变量时很容易出现 (P_{ELY}(t)) 的来源既有风光、又有电网但功率平衡方程漏掉其中某一项的情况。再检查储氢罐的初末状态约束。如果模型中要求 (S_{H2}(0) S_{H2}(24))但氢负荷总量过大导致系统根本无法在24小时内完成一个完整的充放循环这时候模型就不可行。解决办法是放宽末状态约束只要求末容量不低于某一下限。最后检查变量上下限是否有交集。比如电解槽最低技术出力为30kW但某时段系统总可供电功率只有20kW这个时段无论如何都无法满足出力下限。遇到这种问题要么降低技术出力下限要么放宽功率平衡约束。7.2 为什么和文献别人的代码结果对不上一个很常见的原因是目标函数中的正负号不统一。收益最大化问题里有些代码用“负成本”表示收益有些则把售电收益作为负号加到成本函数里。表面上公式一样但拼进总目标函数时如果某一项符号搞反求解结果会完全不同。另一个容易被忽略的问题是时间索引不一致。储氢罐的 (S_{H2}) 到底代表时刻初还是时刻末电解槽的功率是瞬时功率还是小时平均功率如果不统一储能动态方程的时间错位会导致结果出现诡异的振荡。我的建议非常直接把所有变量的单位在代码注释里写得明明白白。功率单位kW电量单位kWh氢气单位kg时间步长1小时。这样排查起来省一半时间。7.3 ADMM/迭代法不收敛的调整经验如果读者后续把模型扩展到含不确定性的分布式求解很可能会用到ADMM。ADMM的核心参数是惩罚系数 (\rho)。太小时原始残差下降慢收敛慢太大时对偶变量振荡结果不光滑。我实测的一个有效调整策略是初始 (\rho) 设为0.1每迭代50次检查原始残差与对偶残差如果原始残差/对偶残差大于10倍则 (\rho) 乘以2如果小于0.1倍则 (\rho) 除以2。这种做法比固定惩罚系数收敛快很多。迭代终止准则建议用相对残差当原始残差和对偶残差同时小于 (10^{-4}) 且迭代次数不超过500次时判定收敛。公式上就是[ |r^k|_2 \le \epsilon^{pri}, \quad |s^k|_2 \le \epsilon^{dual} ]其中 (\epsilon^{pri}) 和 (\epsilon^{dual}) 可以取 (10^{-4}) 数量级具体根据变量量纲调整。7.4 常见问题速查表问题可能原因解决方案求解器报Infeasible功率平衡漏变量逐条核对24时段平衡方程储氢罐容量一直为0动态方程时间索引错位统一时刻初/时刻末定义合作后某主体收益反而降低支付谈判约束未正确加入检查(\alpha_i \le C_i^0)是否成立纳什乘积直接求解不收敛目标函数强非线性改用两步法风光大量弃电但氢能子系统缺电风光与电解槽耦合约束缺失检查(P_{H2,WT}P_{H2,PV}P_{ELY})约束Gurobi报license错误许可证过期重装学术许可证Yalmip报找不到求解器未正确安装工具箱setpath后重启MATLAB8. 项目代码的扩展方向这套基于纳什谈判的光-风-氢合作运行代码在实际使用中其实还可以沿着几个方向直接扩展而且扩展难度不高。考虑碳交易机制时只需在目标函数中加入碳交易成本和碳配额约束本质上是在联盟总成本里增加一项线性成本模型的求解难度不会大幅增加。考虑不确定性的随机优化或鲁棒优化时可以把风速和辐照度场景化形成多场景的纳什谈判模型。这时候第二步支付谈判依然是线性的但第一步联盟效益最大化会变成大场景耦合优化问题对求解器内存有一定要求。考虑动态博弈时可以把单时段周期优化扩展为模型预测控制滚动优化在每一个调度周期内重新求解纳什谈判问题。虽然计算量上去了但系统应对实时波动的能力会明显增强。考虑到多主体信息隐私保护的需求可以改成ADMM分布式求解架构。这也是我目前正在做的扩展方向核心思路是让每个主体只传递边界交换功率和价格信号不暴露各自的成本函数和内部分布式决策信息。9. 实操心得与最后的经验之谈坦白讲纳什谈判理论在风光氢系统中的应用并不是一个数学上特别高深的方向它真正的门槛在于把博弈论的思想转化为一个个可求解、可验证、可复现的工程模型。我在最初做这个项目时最大的弯路就是试图把纳什乘积作为一个整体目标函数直接求解结果绕了好几天才明白两步法的精妙之处。如果你也在做类似的课题我建议你先从三主体的简单场景开始哪怕只有一个典型日的24小时数据先把独立运行、联盟优化、支付谈判三个模块跑通再逐步加设备、加约束、加不确定性。千万别一上来就搞多场景大模型调试起来会非常痛苦。另外把仿真结果做成图非常重要。把24小时功率平衡图、储氢量变化图、成本分配柱状图画出来不仅方便你自己检查结果是否合理在论文写作和项目汇报时也是最有说服力的素材。代码库里我保留了独立运行模型和合作运行模型两套完整求解逻辑方便随时对比算例。如果你在实际复现过程中遇到模型不可行或者分配不稳定的问题欢迎对照这篇文章里的排查顺序逐条核对。