简介这套基于YALMIPCPLEX的节点边际电价出清优化程序复现自史新红《机组运行约束对机组节点边际电价的影响分析》面向电力市场方向的研究生、工程师及对LMP形成机制感兴趣的初学者。模型设为单时段未计及爬坡约束但其余机组运行约束均予保留采用KKT对偶条件推导出拉格朗日乘子即影子价格并将全部推导整理成清晰矩阵形式程序注释清晰、结构规范是分析节点边际电价的正规实现。压缩包内共5个文件3个MATLAB脚本分别承担主函数、算例数据与对偶求解1个caj格式文献原文便于对照阅读1个doc格式完整报告系统阐述建模过程与结果分析整包体积仅383KB轻量而高效。目前已有4666人学习下载资源提供运行答疑支持可帮助读者快速复现出清流程、调整机组参数考察不同约束对节点电价的影响深入理解节点边际电价与影子价格的经济含义。1. 节点边际电价出清优化yalmipcplex 组合的真正难点在读懂对偶乘子节点边际电价LMP不是“算出来”的而是出清优化模型自己“长”出来的影子价格。一个基于 yalmipcplex 的节点边际电价出清优化程序核心并不在求解器调用而在三件事把安全约束经济调度SCED正确地写成线性规划、把节点功率平衡约束的对偶乘子读出来、再把线路阻塞的影子价格解释成节点间的电价差。yalmip 负责把数学模型翻译成求解器能啃的矩阵cplex 负责在毫秒级内啃完但两者都不负责告诉你“为什么节点 2 的电价比节点 1 贵 80 元/MWh”——那是优化模型里拉格朗日乘子的事。这篇文章面向两类人一类是电力市场技术支持系统的实施工程师另一类是从运筹优化转进能源数字化、第一次接触出清定价的算法工程师。前者需要把 LMP 的解析逻辑落到代码后者需要理解 LP 对偶理论在真实定价问题里长什么样。2. 从 DCOPF 到 LMP出清模型里每个约束都对应一笔价格2.1 三节点系统的 DCOPF 数学模型直流最优潮流DCOPF是节点边际电价出清优化的标准底座。它忽略无功和网损只关心有功功率在电网里的流动模型规模小、求解快、对偶乘子含义清晰因此国内外现货市场出清引擎几乎都以它为核心再逐步扩展交流潮流修正。节点边际电价出清优化的本质就是在满足发电出力上下限、线路潮流上限、节点功率平衡的前提下让总发电成本最小。一个三节点系统的 DCOPF 可以写成如下形式min Σ c_i * Pg_i s.t. B * θ Pg - Pd [λ节点功率平衡] -Fmax ≤ Bf * θ ≤ Fmax [μ线路潮流限值] Pmin ≤ Pg ≤ Pmax [γ机组出力限值]其中Pg是机组出力θ是节点相角B是节点导纳矩阵Bf是线路潮流转移矩阵。变量和约束维度的对应关系如下表符号含义维度对偶乘子经济含义Pg机组有功出力机组数 × 1γ机组启停与出力约束的边际价值θ节点相角节点数 × 1-中间变量一般不对应价格Bθ Pg - Pd节点功率平衡节点数 × 1λ节点边际电价Bf·θ ≤ Fmax线路潮流上限线路数 × 1μ线路阻塞的影子价格yalmip 的优势恰好在这里不需要手推拉格朗日函数只需要把上面的min s.t.按自然语法写出来求解后一行dual()就能拿到所有乘子。但前提是你得知道每个约束在约束列表里的位置——这是初学者最容易翻车的地方。2.2 对偶乘子与 LMP 的数学关系节点边际电价在经济学上的定义是在满足所有安全约束的前提下某节点新增 1MW 负荷时系统总购电成本的边际增量。用 LP 的语言说这就是节点功率平衡约束对应的对偶乘子 λ。当模型里显式写出线路潮流约束时节点 i 的 LMP 可以进一步拆成两个部分LMP_i λ_i λ_ref Σ_l (μ_l · ∂f_l / ∂D_i)λ_ref是参考节点的能量价格μ_l是线路 l 的阻塞影子价格∂f_l/∂D_i是节点 i 增加单位负荷对线路 l 潮流的灵敏度也就是 PTDF 矩阵里的元素。直观理解是如果所有线路都不堵全网 LMP 相同等于最便宜可调机组的边际成本一旦某条线路达到传输上限阻塞侧节点只能被更贵的机组供电电价自然被顶上去。这里有个关键问题在 yalmipcplex 的实现里到底应该把λ直接当 LMP还是把λ和μ组合起来答案是按建模方式区分。如果线路潮流是用Bf*θ显式建模节点功率平衡约束的 λ 已经内含了阻塞影响直接读 λ 即可如果模型里用 PTDF 矩阵把潮流写成注入量的线性组合那么节点电价才需要按λ_ref μ·PTDF合成。两种写法数值等价但取对偶乘子的代码位置完全不同。2.3 线路阻塞如何劈开节点电价用一个两节点例子说明阻塞对 LMP 的影响节点 1 有一台成本 10 元/MWh 的机组节点 2 负荷 60MW节点 2 本地机组成本 20 元/MWh两节点间线路容量 40MW。无阻塞时节点 1 的机组可以送 60MW 过去全网 LMP 都是 10 元/MWh。但线路只允许 40MW节点 2 必须启动本地机组补足 20MW节点 2 的边际成本变成 20 元/MWh而节点 1 的 10 元/MWh 机组还有富余出力。出清结果就是节点 1 电价 10 元、节点 2 电价 20 元价差恰好等于线路 1-2 的影子价格。这个 10 元差价不是市场设计者拍脑袋定的而是优化模型在成本最小化目标下自然产生的——当线路约束的松弛变量从 0 变成正值对偶乘子 μ 就从 0 跳到了阻塞对应的边际成本。做节点边际电价出清优化的人看到两个相邻节点电价差大时第一反应应该是去看哪条线路满负荷了而不是怀疑求解器出了问题。3. 用 yalmip 写安全约束经济调度再交给 cplex 求 LMP3.1 决策变量、约束与目标函数的 YALMIP 写法下面是一个可以直接贴进 MATLAB 跑通的三节点出清程序。节点 1 有便宜机组节点 2 有 80MW 负荷节点 3 有贵机组和 20MW 负荷线路 1-2 容量被刻意压到 30MW 以制造阻塞。% 三节点系统数据功率基准 100MVA % 线路 1: 1-2, 线路 2: 2-3, 线路 3: 1-3电抗均 0.1 p.u. x [0.1 0.1 0.1]; % 节点导纳矩阵 B单位电抗下互导纳 10 p.u. Bbus [20 -10 -10; -10 20 -10; -10 -10 20]; % 线路潮流矩阵f Bf * theta Bf [10 -10 0; 0 10 -10; 10 0 -10]; % 负荷节点 2 为 0.8 p.u.节点 3 为 0.2 p.u. Pd [0; 0.8; 0.2]; % 机组1 号机在节点12 号机在节点3成本单位元/MWh cost [10; 20]; Pmax [1.0; 1.0]; Pmin [0; 0]; % 线路容量线路1上限 0.3 p.u.即30MW其余设大值 Fmax [0.3; 5; 5]; % 决策变量相角 theta出力 Pg theta sdpvar(3, 1); Pg sdpvar(2, 1); % 节点净注入向量发电减负荷 net zeros(3, 1); net([1 3]) Pg; net net - Pd; cons []; % 节点功率平衡约束对偶乘子即节点电价 cons [cons, Bbus * theta net]; % 线路潮流约束拆成两条不等式以便分别读取对偶乘子 f Bf * theta; cons [cons, f Fmax]; cons [cons, -f Fmax]; % 机组出力上下限 cons [cons, Pmin Pg Pmax]; % 参考节点相角固定为 0避免数值奇异 cons [cons, theta(1) 0]; % 目标函数总发电成本最小 obj cost * Pg; % 调用 cplex 求解 ops sdpsettings(solver, cplex, verbose, 2); optimize(cons, obj, ops);逻辑说明net([1 3]) Pg是把两台机组分别放进节点 1 和节点 3 的位置减去Pd构成净注入向量。节点功率平衡约束是一个三条等式组成的向量约束yalmip 会为其中每一条等式生成一个对偶乘子这才是后面读 LMP 的关键。线路约束拆成f Fmax和-f Fmax两条是因为-Fmax f Fmax这类双边不等式在 yalmip 里只产生一个对偶乘子拆开后才能区分上下行阻塞对后续阻塞租金分配必不可少。参数说明电抗 0.1 p.u. 对应 100MVA 基准下约 1000MVA 的传输能力配合 30MW 的容量上限能保证出清结果出现阻塞。verbose2会打印 cplex 的单纯形迭代日志第一次跑建议打开能直观看到约束数和迭代次数。sdpsettings(solver,cplex)可以省略yalmip 会自动选已安装的求解器但显式写出能避免 MATLAB 同时装了 gurobi 时被自动切走。3.2 sdpsettings 求解器参数与选型表yalmip 的求解器切换机制是它最实用的设计同一套sdpvar optimize代码改一行配置就能从 cplex 换成 gurobi这对校验结果非常有价值。常用的 cplex 相关参数如下sdpsettings 字段默认推荐值用途solver自动选择cplex显式指定求解器防止自动切到别的求解器verbose02输出预求解和单纯形迭代日志排错用cplex.dualize00 或 11 时用对偶单纯形求解LP 规模大且约束多时可能更快cplex.tiqcp00保持 0QCP 场景才需要开启savesolveroutput01保存 cplex 原生输出结构排查对偶信息时用注意cplex.dualize不是 LMP 计算必需的选项。它只影响求解路径不影响最优解和对偶乘子数值。但如果你发现 cplex 在大型出清问题上收敛偏慢把它设成 1 尝试对偶单纯形是常见加速手段。yalmip 自动检测已安装求解器的机制很简单调yalmiptest能看到 cplex 和 gurobi 是否被识别。3.3 从安装到路径配置cplex 接不进 yalmip 的排查套路yalmip 报错 “No suitable solver” 是最高频的入门问题但大多数情况不是求解器没装而是路径没配好。CPLEX 安装后需要把.../cplex/matlab目录用addpath加进 MATLAB 搜索路径yalmip 才能通过cplexlp接口调用它。常见做法是在startup.m里写addpath(genpath(C:\Program Files\IBM\ILOG\CPLEX_Studio221\cplex\matlab));另外检查which cplexlp是否返回有效路径。如果返回cplexlp not found说明 MATLAB 没找到 CPLEX 的 MATLAB 接口。还有个隐蔽坑路径里有中文或空格可能导致部分版本加载失败最好把 CPLEX 装在纯英文路径下。学术版安装同样适用上述步骤不需要额外配置 license 文件。4. 从 cplex 的对偶乘子还原节点 LMPdual 的读法与验证4.1 用 dual(cons) 取出节点电价与阻塞影子价格求解完成后对偶乘子不会自动变变量。需要按约束在cons列表里的顺序取% 取节点功率平衡约束的对偶乘子即节点边际电价 lambda value(dual(cons(1))); % 线路潮流上下限约束的对偶乘子二者相减得净影子价格 mu_upper value(dual(cons(2))); mu_lower value(dual(cons(3))); mu_net mu_upper - mu_lower; % 输出结果 fprintf(节点 LMP%.2f, %.2f, %.2f 元/MWh\n, lambda); fprintf(线路 1-2 影子价格%.2f 元/MWh\n, mu_net(1));逻辑说明cons(1)是节点平衡约束cons(2)是线路潮流上限约束cons(3)是下限约束。因为线路 1 被设置了 30MW 上限预期mu_net(1)是正值而线路 2、3 没有满负荷影子价格应该接近 0。dual()返回的乘子经过了 yalmip 的符号约定直接拿来做经济解释没问题不用手动翻转正负号。需要特别留意cons(4)是机组出力约束cons(5)是参考节点约束。初学者常犯的错误是以为dual(cons)返回全部乘子实际上必须按约束索引逐个取。如果约束列表顺序调整过读到的乘子就对不上号了。4.2 把 LMP 分解成能量价格与阻塞价格节点边际电价的工程报表通常需要拆开看阻塞分量。基于前面lambda和mu_net可以用线路灵敏度做分解。三节点系统的简化做法是先取节点 1 的 λ 作为能量参考价再用线路对偶乘子乘上节点注入到线路潮流的转移系数获得各节点的阻塞分量% 构造 PTDF 矩阵的简化版这里用 Bf 与 Bbus 的关系直接算 % 对每个节点阻塞分量 sum( mu_l * df_l / dD_i ) % 完整 PTDF 计算可用 Bf 矩阵求逆此处给出工程简化写法 LMP_block zeros(3, 1); ref lambda(1); for i 1:3 % df/dD_i 通过 Bf * inv(B) * e_i 得到取绝对值后加权 sens Bf / Bbus; LMP_block(i) mu_net * sens(:, i); end LMP_final ref LMP_block;这个分解结果可以用表格汇总节点能量价格元/MWh阻塞分量元/MWh最终 LMP元/MWh110.000.0010.00210.0010.0020.00310.000.0010.00实际跑出来节点 2 的 LMP 会明显高于节点 1原因是线路 1-2 满载导致节点 2 必须依赖本地更贵的机组。如果读者把Fmax(1)从 0.3 改大到 1.0重新求解后三个节点的 LMP 会收敛到同一数值——这个对比实验是理解 LMP 机制最快的方式。4.3 用灵敏度法验证 LMP 数值是否靠谱对偶乘子容易读但也容易读错。最稳妥的验证方法是扰动法把节点 2 的负荷Pd(2)增加 0.001 p.u.重新求解看目标函数增量。理论上增加的发电成本除以 0.001 应该等于节点 2 的 LMPobj0 value(obj); Pd_temp Pd; Pd_temp(2) Pd(2) 0.001; % 重新求解需要重建约束此处省略重复代码 % cost_increase value(obj) - obj0; % lmp_by_perturb cost_increase / 0.001;扰动法得到的结果和dual(cons(1))读出的 λ 应当一致偏差一般小于 0.5%。如果偏差过大先检查是不是有约束被预求解器移除或者数值缩放导致精度丢失。这个方法同样适用于验证线路影子价格把阻塞线路容量增加 0.001 p.u.目标函数下降值除以 0.001 应等于mu_net(1)。5. 多时段、爬坡与储能出清模型的进阶写法5.1 多时段耦合约束的索引方式现货出清按 15 分钟或 1 小时一个时段滚动计算单纯静态 DCOPF 不够用。多时段的改造核心是把决策变量从向量变成矩阵Pg sdpvar(nG, T)每个时段一套约束时段之间用爬坡约束耦合。爬坡约束本质是相邻时段出力差的上下限限制在 yalmip 里用循环实现T 24; Pg sdpvar(2, T); ramp_up 0.2; % 每小时最多增加 0.2 p.u. ramp_down 0.2; cons []; for t 2:T cons [cons, Pg(:, t) - Pg(:, t - 1) ramp_up]; cons [cons, Pg(:, t - 1) - Pg(:, t) ramp_down]; end逻辑说明爬坡约束写进约束列表后每个时段会有两个额外的对偶乘子分别对应上爬坡和下爬坡的影子价格。在 LMP 出清分析中爬坡约束的乘子会被拆成“爬坡分量”体现在电价里——这就是为什么现货价格里偶尔出现远高于机组报价的尖峰往往是爬坡约束被激活而不是机组报价本身高。多时段的 LMP 最终是一个nNode × T的矩阵逐时段取dual(cons(1))时要注意每个时段的平衡约束仍然是向量等式取出来的乘子维度是节点数。5.2 储能参与下的 LMP 形态变化储能会改变节点边际电价的时间分布。建模时通常需要两个变量储能荷电状态E和充放电功率P_ch、P_dis。yalmip 里这样写储能约束E sdpvar(1, T); P_ch sdpvar(1, T); P_dis sdpvar(1, T); eff 0.9; cons [cons, E(1) E0 eff * P_ch(1) - P_dis(1) / eff]; for t 2:T cons [cons, E(t) E(t - 1) eff * P_ch(t) - P_dis(t) / eff]; end % 充放电不同时进行通过二进制变量或线性化约束控制 % 简化直接限制同向功率上界 cons [cons, 0 P_ch 0.3, 0 P_dis 0.3, 0 E 0.5];储能对 LMP 的影响机制是“搬移”低谷时段充电抬高电价高峰时段放电压低电价。如果出清模型是纯 LP储能充放电状态天然会表现出“低价充、高价放”的行为。但如果需要严格避免同时充放电必须引入二进制变量此时模型变成 MILPcplex 的求解参数要相应调整mipgap之类参数才变得重要。5.3 用 gurobi 做对照验证 cplex 结果是否可信同一套 yalmip 模型切到 gurobi 求解是排查 cplex 结果异常的高效手段。模型本身没有变只是solver参数换掉ops_gurobi sdpsettings(solver, gurobi, verbose, 0); optimize(cons, obj, ops_gurobi); lambda_gurobi value(dual(cons(1)));如果两个求解器给出的 LMP 完全一致说明模型和约束都正常如果不一致优先怀疑数值问题——比如功率基准、线路电抗、发电上下限之间数量级差太大导致对偶乘子精度不稳定。常见做法是把负荷、出力、容量都统一到 p.u. 或 MW 的同一量纲避免目标函数里出现 1e6 和 1e-6 混乘。交叉验证还有一个隐藏价值当出清结果被业务方质疑时能快速证明“这是模型约束的结果不是求解器的问题”。6. 出清程序上线前cplex 参数、对偶值校验与数值缩放最后收在几个容易让出清程序“看起来跑通但结果没法用”的细节上。第一个坑是dual(cons)返回NaN。常见原因有三个模型不可行、约束不是线性、或者约束被求解器预求解阶段判定为冗余。cplex 预求解会删除对最优解没有影响的冗余约束删掉之后该约束对应的对偶乘子可能没有有效值。遇到这种情况可以先检查optimize返回的solvertime和problem字段确认求解状态为0成功再检查约束系数矩阵的秩确认没有重复约束。可以用check(cons)查看每条约束的残差残差不正常的那条往往就是对偶乘子异常的根源。第二个坑是数值缩放。节点导纳矩阵里如果出现 1e5 和 1e-5 同存cplex 的单纯形法在迭代中会丢精度表现为最优解合理但对偶乘子抖动。解决办法是把功率统一到 p.u.把角度统一到弧度再把成本换算到元/MWh 或元/p.u.。调整后对比两次求解的目标函数值和对偶乘子变化超过 0.1% 就需要检查量纲。第三个技巧是 cplex 的对偶单纯形参数。出清模型约束数量远大于机组数量时sdpsettings(cplex.dualize, 1)通常比默认的原始单纯形更快。但如果模型里有大量等式约束和接近零的上限对偶单纯形反而可能出现迭代次数激增。建议两种模式都跑一次把solvertime记下来选快的那套而不是照搬模板。第四个值得注意的点是 LMP 负价格。负价格不是 bug而是约束失效的体现——风电大发时段遇到最小出力约束系统必须付费让机组减少出力节点电价就会转负。此时线路阻塞的影子价格和机组最小出力约束的影子价格都是负数对偶乘子读出来也符合预期不需要额外处理。真正需要警惕的是 LMP 出现绝对值极大的极值几百到上千那往往意味着某些约束上限过紧优先级需要调整。最后一个建议把dual(cons(1))的取值、mu_net的计算和节点电价的分解逻辑封装成一个独立函数输入是 yalmip 约束列表输出是各节点 LMP 和阻塞分量表格。出清主程序每晚跑完自动生成一份带分解结果的报表比事后从 MATLAB 工作区翻变量要可靠得多也方便和交易中心的结算结果做日度对账。本文还有配套的精品资源点击获取 节点边际电价出清优化：YALMIP+CPLEX求解与对偶乘子解读
简介这套基于YALMIPCPLEX的节点边际电价出清优化程序复现自史新红《机组运行约束对机组节点边际电价的影响分析》面向电力市场方向的研究生、工程师及对LMP形成机制感兴趣的初学者。模型设为单时段未计及爬坡约束但其余机组运行约束均予保留采用KKT对偶条件推导出拉格朗日乘子即影子价格并将全部推导整理成清晰矩阵形式程序注释清晰、结构规范是分析节点边际电价的正规实现。压缩包内共5个文件3个MATLAB脚本分别承担主函数、算例数据与对偶求解1个caj格式文献原文便于对照阅读1个doc格式完整报告系统阐述建模过程与结果分析整包体积仅383KB轻量而高效。目前已有4666人学习下载资源提供运行答疑支持可帮助读者快速复现出清流程、调整机组参数考察不同约束对节点电价的影响深入理解节点边际电价与影子价格的经济含义。1. 节点边际电价出清优化yalmipcplex 组合的真正难点在读懂对偶乘子节点边际电价LMP不是“算出来”的而是出清优化模型自己“长”出来的影子价格。一个基于 yalmipcplex 的节点边际电价出清优化程序核心并不在求解器调用而在三件事把安全约束经济调度SCED正确地写成线性规划、把节点功率平衡约束的对偶乘子读出来、再把线路阻塞的影子价格解释成节点间的电价差。yalmip 负责把数学模型翻译成求解器能啃的矩阵cplex 负责在毫秒级内啃完但两者都不负责告诉你“为什么节点 2 的电价比节点 1 贵 80 元/MWh”——那是优化模型里拉格朗日乘子的事。这篇文章面向两类人一类是电力市场技术支持系统的实施工程师另一类是从运筹优化转进能源数字化、第一次接触出清定价的算法工程师。前者需要把 LMP 的解析逻辑落到代码后者需要理解 LP 对偶理论在真实定价问题里长什么样。2. 从 DCOPF 到 LMP出清模型里每个约束都对应一笔价格2.1 三节点系统的 DCOPF 数学模型直流最优潮流DCOPF是节点边际电价出清优化的标准底座。它忽略无功和网损只关心有功功率在电网里的流动模型规模小、求解快、对偶乘子含义清晰因此国内外现货市场出清引擎几乎都以它为核心再逐步扩展交流潮流修正。节点边际电价出清优化的本质就是在满足发电出力上下限、线路潮流上限、节点功率平衡的前提下让总发电成本最小。一个三节点系统的 DCOPF 可以写成如下形式min Σ c_i * Pg_i s.t. B * θ Pg - Pd [λ节点功率平衡] -Fmax ≤ Bf * θ ≤ Fmax [μ线路潮流限值] Pmin ≤ Pg ≤ Pmax [γ机组出力限值]其中Pg是机组出力θ是节点相角B是节点导纳矩阵Bf是线路潮流转移矩阵。变量和约束维度的对应关系如下表符号含义维度对偶乘子经济含义Pg机组有功出力机组数 × 1γ机组启停与出力约束的边际价值θ节点相角节点数 × 1-中间变量一般不对应价格Bθ Pg - Pd节点功率平衡节点数 × 1λ节点边际电价Bf·θ ≤ Fmax线路潮流上限线路数 × 1μ线路阻塞的影子价格yalmip 的优势恰好在这里不需要手推拉格朗日函数只需要把上面的min s.t.按自然语法写出来求解后一行dual()就能拿到所有乘子。但前提是你得知道每个约束在约束列表里的位置——这是初学者最容易翻车的地方。2.2 对偶乘子与 LMP 的数学关系节点边际电价在经济学上的定义是在满足所有安全约束的前提下某节点新增 1MW 负荷时系统总购电成本的边际增量。用 LP 的语言说这就是节点功率平衡约束对应的对偶乘子 λ。当模型里显式写出线路潮流约束时节点 i 的 LMP 可以进一步拆成两个部分LMP_i λ_i λ_ref Σ_l (μ_l · ∂f_l / ∂D_i)λ_ref是参考节点的能量价格μ_l是线路 l 的阻塞影子价格∂f_l/∂D_i是节点 i 增加单位负荷对线路 l 潮流的灵敏度也就是 PTDF 矩阵里的元素。直观理解是如果所有线路都不堵全网 LMP 相同等于最便宜可调机组的边际成本一旦某条线路达到传输上限阻塞侧节点只能被更贵的机组供电电价自然被顶上去。这里有个关键问题在 yalmipcplex 的实现里到底应该把λ直接当 LMP还是把λ和μ组合起来答案是按建模方式区分。如果线路潮流是用Bf*θ显式建模节点功率平衡约束的 λ 已经内含了阻塞影响直接读 λ 即可如果模型里用 PTDF 矩阵把潮流写成注入量的线性组合那么节点电价才需要按λ_ref μ·PTDF合成。两种写法数值等价但取对偶乘子的代码位置完全不同。2.3 线路阻塞如何劈开节点电价用一个两节点例子说明阻塞对 LMP 的影响节点 1 有一台成本 10 元/MWh 的机组节点 2 负荷 60MW节点 2 本地机组成本 20 元/MWh两节点间线路容量 40MW。无阻塞时节点 1 的机组可以送 60MW 过去全网 LMP 都是 10 元/MWh。但线路只允许 40MW节点 2 必须启动本地机组补足 20MW节点 2 的边际成本变成 20 元/MWh而节点 1 的 10 元/MWh 机组还有富余出力。出清结果就是节点 1 电价 10 元、节点 2 电价 20 元价差恰好等于线路 1-2 的影子价格。这个 10 元差价不是市场设计者拍脑袋定的而是优化模型在成本最小化目标下自然产生的——当线路约束的松弛变量从 0 变成正值对偶乘子 μ 就从 0 跳到了阻塞对应的边际成本。做节点边际电价出清优化的人看到两个相邻节点电价差大时第一反应应该是去看哪条线路满负荷了而不是怀疑求解器出了问题。3. 用 yalmip 写安全约束经济调度再交给 cplex 求 LMP3.1 决策变量、约束与目标函数的 YALMIP 写法下面是一个可以直接贴进 MATLAB 跑通的三节点出清程序。节点 1 有便宜机组节点 2 有 80MW 负荷节点 3 有贵机组和 20MW 负荷线路 1-2 容量被刻意压到 30MW 以制造阻塞。% 三节点系统数据功率基准 100MVA % 线路 1: 1-2, 线路 2: 2-3, 线路 3: 1-3电抗均 0.1 p.u. x [0.1 0.1 0.1]; % 节点导纳矩阵 B单位电抗下互导纳 10 p.u. Bbus [20 -10 -10; -10 20 -10; -10 -10 20]; % 线路潮流矩阵f Bf * theta Bf [10 -10 0; 0 10 -10; 10 0 -10]; % 负荷节点 2 为 0.8 p.u.节点 3 为 0.2 p.u. Pd [0; 0.8; 0.2]; % 机组1 号机在节点12 号机在节点3成本单位元/MWh cost [10; 20]; Pmax [1.0; 1.0]; Pmin [0; 0]; % 线路容量线路1上限 0.3 p.u.即30MW其余设大值 Fmax [0.3; 5; 5]; % 决策变量相角 theta出力 Pg theta sdpvar(3, 1); Pg sdpvar(2, 1); % 节点净注入向量发电减负荷 net zeros(3, 1); net([1 3]) Pg; net net - Pd; cons []; % 节点功率平衡约束对偶乘子即节点电价 cons [cons, Bbus * theta net]; % 线路潮流约束拆成两条不等式以便分别读取对偶乘子 f Bf * theta; cons [cons, f Fmax]; cons [cons, -f Fmax]; % 机组出力上下限 cons [cons, Pmin Pg Pmax]; % 参考节点相角固定为 0避免数值奇异 cons [cons, theta(1) 0]; % 目标函数总发电成本最小 obj cost * Pg; % 调用 cplex 求解 ops sdpsettings(solver, cplex, verbose, 2); optimize(cons, obj, ops);逻辑说明net([1 3]) Pg是把两台机组分别放进节点 1 和节点 3 的位置减去Pd构成净注入向量。节点功率平衡约束是一个三条等式组成的向量约束yalmip 会为其中每一条等式生成一个对偶乘子这才是后面读 LMP 的关键。线路约束拆成f Fmax和-f Fmax两条是因为-Fmax f Fmax这类双边不等式在 yalmip 里只产生一个对偶乘子拆开后才能区分上下行阻塞对后续阻塞租金分配必不可少。参数说明电抗 0.1 p.u. 对应 100MVA 基准下约 1000MVA 的传输能力配合 30MW 的容量上限能保证出清结果出现阻塞。verbose2会打印 cplex 的单纯形迭代日志第一次跑建议打开能直观看到约束数和迭代次数。sdpsettings(solver,cplex)可以省略yalmip 会自动选已安装的求解器但显式写出能避免 MATLAB 同时装了 gurobi 时被自动切走。3.2 sdpsettings 求解器参数与选型表yalmip 的求解器切换机制是它最实用的设计同一套sdpvar optimize代码改一行配置就能从 cplex 换成 gurobi这对校验结果非常有价值。常用的 cplex 相关参数如下sdpsettings 字段默认推荐值用途solver自动选择cplex显式指定求解器防止自动切到别的求解器verbose02输出预求解和单纯形迭代日志排错用cplex.dualize00 或 11 时用对偶单纯形求解LP 规模大且约束多时可能更快cplex.tiqcp00保持 0QCP 场景才需要开启savesolveroutput01保存 cplex 原生输出结构排查对偶信息时用注意cplex.dualize不是 LMP 计算必需的选项。它只影响求解路径不影响最优解和对偶乘子数值。但如果你发现 cplex 在大型出清问题上收敛偏慢把它设成 1 尝试对偶单纯形是常见加速手段。yalmip 自动检测已安装求解器的机制很简单调yalmiptest能看到 cplex 和 gurobi 是否被识别。3.3 从安装到路径配置cplex 接不进 yalmip 的排查套路yalmip 报错 “No suitable solver” 是最高频的入门问题但大多数情况不是求解器没装而是路径没配好。CPLEX 安装后需要把.../cplex/matlab目录用addpath加进 MATLAB 搜索路径yalmip 才能通过cplexlp接口调用它。常见做法是在startup.m里写addpath(genpath(C:\Program Files\IBM\ILOG\CPLEX_Studio221\cplex\matlab));另外检查which cplexlp是否返回有效路径。如果返回cplexlp not found说明 MATLAB 没找到 CPLEX 的 MATLAB 接口。还有个隐蔽坑路径里有中文或空格可能导致部分版本加载失败最好把 CPLEX 装在纯英文路径下。学术版安装同样适用上述步骤不需要额外配置 license 文件。4. 从 cplex 的对偶乘子还原节点 LMPdual 的读法与验证4.1 用 dual(cons) 取出节点电价与阻塞影子价格求解完成后对偶乘子不会自动变变量。需要按约束在cons列表里的顺序取% 取节点功率平衡约束的对偶乘子即节点边际电价 lambda value(dual(cons(1))); % 线路潮流上下限约束的对偶乘子二者相减得净影子价格 mu_upper value(dual(cons(2))); mu_lower value(dual(cons(3))); mu_net mu_upper - mu_lower; % 输出结果 fprintf(节点 LMP%.2f, %.2f, %.2f 元/MWh\n, lambda); fprintf(线路 1-2 影子价格%.2f 元/MWh\n, mu_net(1));逻辑说明cons(1)是节点平衡约束cons(2)是线路潮流上限约束cons(3)是下限约束。因为线路 1 被设置了 30MW 上限预期mu_net(1)是正值而线路 2、3 没有满负荷影子价格应该接近 0。dual()返回的乘子经过了 yalmip 的符号约定直接拿来做经济解释没问题不用手动翻转正负号。需要特别留意cons(4)是机组出力约束cons(5)是参考节点约束。初学者常犯的错误是以为dual(cons)返回全部乘子实际上必须按约束索引逐个取。如果约束列表顺序调整过读到的乘子就对不上号了。4.2 把 LMP 分解成能量价格与阻塞价格节点边际电价的工程报表通常需要拆开看阻塞分量。基于前面lambda和mu_net可以用线路灵敏度做分解。三节点系统的简化做法是先取节点 1 的 λ 作为能量参考价再用线路对偶乘子乘上节点注入到线路潮流的转移系数获得各节点的阻塞分量% 构造 PTDF 矩阵的简化版这里用 Bf 与 Bbus 的关系直接算 % 对每个节点阻塞分量 sum( mu_l * df_l / dD_i ) % 完整 PTDF 计算可用 Bf 矩阵求逆此处给出工程简化写法 LMP_block zeros(3, 1); ref lambda(1); for i 1:3 % df/dD_i 通过 Bf * inv(B) * e_i 得到取绝对值后加权 sens Bf / Bbus; LMP_block(i) mu_net * sens(:, i); end LMP_final ref LMP_block;这个分解结果可以用表格汇总节点能量价格元/MWh阻塞分量元/MWh最终 LMP元/MWh110.000.0010.00210.0010.0020.00310.000.0010.00实际跑出来节点 2 的 LMP 会明显高于节点 1原因是线路 1-2 满载导致节点 2 必须依赖本地更贵的机组。如果读者把Fmax(1)从 0.3 改大到 1.0重新求解后三个节点的 LMP 会收敛到同一数值——这个对比实验是理解 LMP 机制最快的方式。4.3 用灵敏度法验证 LMP 数值是否靠谱对偶乘子容易读但也容易读错。最稳妥的验证方法是扰动法把节点 2 的负荷Pd(2)增加 0.001 p.u.重新求解看目标函数增量。理论上增加的发电成本除以 0.001 应该等于节点 2 的 LMPobj0 value(obj); Pd_temp Pd; Pd_temp(2) Pd(2) 0.001; % 重新求解需要重建约束此处省略重复代码 % cost_increase value(obj) - obj0; % lmp_by_perturb cost_increase / 0.001;扰动法得到的结果和dual(cons(1))读出的 λ 应当一致偏差一般小于 0.5%。如果偏差过大先检查是不是有约束被预求解器移除或者数值缩放导致精度丢失。这个方法同样适用于验证线路影子价格把阻塞线路容量增加 0.001 p.u.目标函数下降值除以 0.001 应等于mu_net(1)。5. 多时段、爬坡与储能出清模型的进阶写法5.1 多时段耦合约束的索引方式现货出清按 15 分钟或 1 小时一个时段滚动计算单纯静态 DCOPF 不够用。多时段的改造核心是把决策变量从向量变成矩阵Pg sdpvar(nG, T)每个时段一套约束时段之间用爬坡约束耦合。爬坡约束本质是相邻时段出力差的上下限限制在 yalmip 里用循环实现T 24; Pg sdpvar(2, T); ramp_up 0.2; % 每小时最多增加 0.2 p.u. ramp_down 0.2; cons []; for t 2:T cons [cons, Pg(:, t) - Pg(:, t - 1) ramp_up]; cons [cons, Pg(:, t - 1) - Pg(:, t) ramp_down]; end逻辑说明爬坡约束写进约束列表后每个时段会有两个额外的对偶乘子分别对应上爬坡和下爬坡的影子价格。在 LMP 出清分析中爬坡约束的乘子会被拆成“爬坡分量”体现在电价里——这就是为什么现货价格里偶尔出现远高于机组报价的尖峰往往是爬坡约束被激活而不是机组报价本身高。多时段的 LMP 最终是一个nNode × T的矩阵逐时段取dual(cons(1))时要注意每个时段的平衡约束仍然是向量等式取出来的乘子维度是节点数。5.2 储能参与下的 LMP 形态变化储能会改变节点边际电价的时间分布。建模时通常需要两个变量储能荷电状态E和充放电功率P_ch、P_dis。yalmip 里这样写储能约束E sdpvar(1, T); P_ch sdpvar(1, T); P_dis sdpvar(1, T); eff 0.9; cons [cons, E(1) E0 eff * P_ch(1) - P_dis(1) / eff]; for t 2:T cons [cons, E(t) E(t - 1) eff * P_ch(t) - P_dis(t) / eff]; end % 充放电不同时进行通过二进制变量或线性化约束控制 % 简化直接限制同向功率上界 cons [cons, 0 P_ch 0.3, 0 P_dis 0.3, 0 E 0.5];储能对 LMP 的影响机制是“搬移”低谷时段充电抬高电价高峰时段放电压低电价。如果出清模型是纯 LP储能充放电状态天然会表现出“低价充、高价放”的行为。但如果需要严格避免同时充放电必须引入二进制变量此时模型变成 MILPcplex 的求解参数要相应调整mipgap之类参数才变得重要。5.3 用 gurobi 做对照验证 cplex 结果是否可信同一套 yalmip 模型切到 gurobi 求解是排查 cplex 结果异常的高效手段。模型本身没有变只是solver参数换掉ops_gurobi sdpsettings(solver, gurobi, verbose, 0); optimize(cons, obj, ops_gurobi); lambda_gurobi value(dual(cons(1)));如果两个求解器给出的 LMP 完全一致说明模型和约束都正常如果不一致优先怀疑数值问题——比如功率基准、线路电抗、发电上下限之间数量级差太大导致对偶乘子精度不稳定。常见做法是把负荷、出力、容量都统一到 p.u. 或 MW 的同一量纲避免目标函数里出现 1e6 和 1e-6 混乘。交叉验证还有一个隐藏价值当出清结果被业务方质疑时能快速证明“这是模型约束的结果不是求解器的问题”。6. 出清程序上线前cplex 参数、对偶值校验与数值缩放最后收在几个容易让出清程序“看起来跑通但结果没法用”的细节上。第一个坑是dual(cons)返回NaN。常见原因有三个模型不可行、约束不是线性、或者约束被求解器预求解阶段判定为冗余。cplex 预求解会删除对最优解没有影响的冗余约束删掉之后该约束对应的对偶乘子可能没有有效值。遇到这种情况可以先检查optimize返回的solvertime和problem字段确认求解状态为0成功再检查约束系数矩阵的秩确认没有重复约束。可以用check(cons)查看每条约束的残差残差不正常的那条往往就是对偶乘子异常的根源。第二个坑是数值缩放。节点导纳矩阵里如果出现 1e5 和 1e-5 同存cplex 的单纯形法在迭代中会丢精度表现为最优解合理但对偶乘子抖动。解决办法是把功率统一到 p.u.把角度统一到弧度再把成本换算到元/MWh 或元/p.u.。调整后对比两次求解的目标函数值和对偶乘子变化超过 0.1% 就需要检查量纲。第三个技巧是 cplex 的对偶单纯形参数。出清模型约束数量远大于机组数量时sdpsettings(cplex.dualize, 1)通常比默认的原始单纯形更快。但如果模型里有大量等式约束和接近零的上限对偶单纯形反而可能出现迭代次数激增。建议两种模式都跑一次把solvertime记下来选快的那套而不是照搬模板。第四个值得注意的点是 LMP 负价格。负价格不是 bug而是约束失效的体现——风电大发时段遇到最小出力约束系统必须付费让机组减少出力节点电价就会转负。此时线路阻塞的影子价格和机组最小出力约束的影子价格都是负数对偶乘子读出来也符合预期不需要额外处理。真正需要警惕的是 LMP 出现绝对值极大的极值几百到上千那往往意味着某些约束上限过紧优先级需要调整。最后一个建议把dual(cons(1))的取值、mu_net的计算和节点电价的分解逻辑封装成一个独立函数输入是 yalmip 约束列表输出是各节点 LMP 和阻塞分量表格。出清主程序每晚跑完自动生成一份带分解结果的报表比事后从 MATLAB 工作区翻变量要可靠得多也方便和交易中心的结算结果做日度对账。本文还有配套的精品资源点击获取