简介：电力市场节点边际电价出清优化的完整复现方案，面向电力市场研究人员、高年级本科生及研究生。资源基于史新红论文《机组运行约束对机组节点边际电价的影响分析》，在单时段模型下采用YALMIP+CPLEX求解器，通过KKT对偶条件解出拉格朗日乘子（即影子价格），并以矩阵形式编程实现，清晰展示了节点边际电价（LMP）出清的正规分析流程。压缩包共5个文件，其中3个MATLAB脚本（主程序、对偶求解与算例）、1个caj格式的论文原文和1份Word版完整报告，整体仅383KB，极为轻便。已有4666人学习下载，深受同行关注。程序中注释详尽、结构清晰，除未考虑爬坡约束外，其余机组约束均覆盖，既适合用作电力市场课程大作业或教学示例，也可供工程人员对照理解机组运行约束对节点电价的影响机理；附带的报告系统解释了不同约束下的电价差异，运行遇到问题还可与作者直接答疑交流。
1. 节点边际电价为什么能出清阻塞
当一条联络线越限时，市场出清价格为什么不再是全网统一价？答案在于节点边际电价（LMP）模型。与传统的系统边际电价（SMP）只依赖系统功率平衡不同，LMP 把电网拓扑和线路潮流约束直接写进优化问题，通过求解直流最优潮流（DC-OPF）得到每个节点的边际价格。这套定价机制已成为美国 PJM、北欧电力市场的主流基础，国内现货市场试点也普遍采用。下面从 DC-OPF 建模入手，用 YALMIP 在 MATLAB 中构建优化模型，调用 CPLEX 求解，最终给出可复现的节点边际电价出清程序，并在最后讨论 LMP 分解与参数调优。适合电力市场方向的研究生、交易员和求解器应用工程师。
2. 从系统边际电价到节点边际电价：DC-OPF的数学建模
2.1 为什么SMP给不出阻塞信号
系统边际电价的计算只要求全系统发电总成本最小，满足总发电等于总负荷，完全不考虑线路潮流约束。在无阻塞时，所有节点共享同一个边际成本，这个价格能引导发电和用电在总量上达到平衡。一旦出现输电阻塞，SMP方案可能要求某台位于负荷中心的机组以较高成本大量出力，而远离负荷的便宜机组却无法满发，这时全网统一价就无法反映不同节点的电能稀缺程度，也无法为阻塞提供正确的经济信号。更严重的是，SMP并不产生阻塞影子价格，市场无法识别需要扩容的线路。
节点边际电价（LMP）则把每条线路的潮流上限作为不等式约束纳入优化，求解出拉格朗日乘子。这些乘子会抬升或者压低不同节点的边际价格。直观上讲，阻塞使电能无法自由流动，被阻塞的区域只能靠本区更贵的电源，因此这些节点的LMP高于其他节点。这种差异化价格正是电网物理约束在价格上的体现。
2.2 DC-OPF的线性化与适用边界
完整的交流最优潮流（AC-OPF）包含电压、无功、有功损耗和非线性三角关系，求解复杂，而且价格形成机制不透明。实际节点边际电价出清中，普遍采用直流潮流近似（DC-OPF），它做了三个假设：
- 支路电阻远小于电抗，忽略有功损耗；
- 节点电压幅值近似为1 p.u.；
- 相角差很小，sinθ≈θ。
于是线路有功潮流可以写为：
P_ij = (θ_i - θ_j) / x_ij其中x_ij是支路电抗。这个公式把潮流表达为相角差的线性函数，所有约束都是线性的，整个出清问题退化为线性规划（LP）。如果机组启动成本、最小开停机时间等整数变量也被纳入，则成为混合整数线性规划（MILP），但出清价格仍然来自LP松弛后的对偶信息。DC-OPF的优点是快速稳定，适合日前、实时市场的大规模迭代计算；缺点是无法处理电压和无功问题，在重载或电压敏感场景需要回到AC-OPF校核。
2.3 出清优化问题的决策变量与约束
以三节点系统为例，假设每个节点均有一台发电机组，节点负荷已知。系统参数如下表所示。
| 参数 | 节点1 | 节点2 | 节点3 |
|---|---|---|---|
| 发电边际成本 ($/MWh) | 20 | 30 | 25 |
| 出力下限 (MW) | 0 | 0 | 0 |
| 出力上限 (MW) | 100 | 100 | 100 |
| 负荷 (MW) | 0 | 80 | 60 |
| 支路 | 起点 | 终点 | 电抗 (p.u.) | 潮流上限 (MW) |
|---|---|---|---|---|
| 1 | 1 | 2 | 0.20 | 50 |
| 2 | 1 | 3 | 0.30 | 50 |
| 3 | 2 | 3 | 0.25 | 60 |
定义机组出力向量P_g（3×1），节点相角向量θ（3×1），线路潮流向量F（3×1）。优化目标为系统总发电成本最小：
min Σ (c_i * P_g,i)约束条件如下：
- 节点功率平衡：对于每个节点，注入功率等于流出功率。用关联矩阵A描述，即 A·F = P_g - P_d；
- 线路潮流定义：F_k = (θ_i - θ_j) / x_k；
- 线路潮流限值：-F_max ≤ F ≤ F_max；
- 机组出力限值：P_min ≤ P_g ≤ P_max；
- 参考节点相角：θ_ref = 0。
关联矩阵A的构造方法很简单：对第k条支路，若起点为i、终点为j，则A(i,k)=1，A(j,k)=-1，其余为0。下面是生成A矩阵和构建节点平衡的MATLAB代码片段，这段代码不依赖任何优化工具箱，只用来准备模型数据。
nb = 3; nl = 3; br = [1 2; 1 3; 2 3]; % 支路起点, 终点 xij = [0.20; 0.30; 0.25]; A = zeros(nb, nl); for k = 1:nl I = br(k,1); J = br(k,2); A(I,k) = 1; A(J,k) = -1; end % A的第k列就对应第k条支路的方向这段代码中，A(i,k)=1表示支路k从节点i流出，A(j,k)=-1表示流入节点j。节点平衡写为A*F == Pg - Pd，即净出力等于线路流出与流入的差值。这里不需要手动列出每个节点的潮流方程，矩阵会自行处理。注意参考节点相角必须固定，否则节点平衡约束线性相关，优化问题会退化出无穷多解。
2.4 节点边际电价来自对偶变量
节点边际电价的本质是节点功率平衡约束的影子价格。在LP问题中，影子价格表示该约束右端项变动一个单位时目标函数的变化量。对节点i来说，如果i节点负荷增加1 MW，系统总成本会上升λ_i，这个λ_i就是该节点的LMP。由于线路约束的存在，不同节点的λ_i往往不同，线路越限越严重，差距越大。后面会看到，使用YALMIP可以很自然地通过dual()函数提取这个值。
注意，系统边际电价SMP可以看作忽略所有线路约束时的特殊LMP：当没有阻塞时，所有节点平衡约束的对偶变量相同，LMP退化为全网统一价。这也是验证模型是否正确的一个间接方法。
3. YALMIP+CPLEX 实现 LMP 出清：核心代码与参数
3.1 YALMIP建模的数据准备
在使用YALMIP之前，需要先安装CPLEX并确保MATLAB能找到其求解器。安装完成后，可以在MATLAB中运行yalmiptest检查是否有CPLEX可用。下面的代码定义3节点系统的完整数据，并构建关联矩阵。为了简洁，假设每个节点一台机组，且发电成本线性，不考虑空载成本。
% 3节点DC-OPF出清，YALMIP + CPLEX clear; clc; nb = 3; ng = nb; c = [20; 30; 25]; % 边际成本 ($/MWh) pmin = zeros(ng,1); pmax = ones(ng,1)*100; % 出力上下限 pd = [0; 80; 60]; % 节点负荷 br = [1 2; 1 3; 2 3]; % 支路起点终点 xij = [0.20; 0.30; 0.25]; % 支路电抗 limit = [50; 50; 60]; % 潮流上限 MW ref = 1; % 参考节点 % 构造关联矩阵 A = zeros(nb, size(br,1)); for k = 1:size(br,1) A(br(k,1), k) = 1; A(br(k,2), k) = -1; end这块代码中，c、pmin、pmax、pd都是列向量，符合YALMIP变量维度的习惯。A矩阵的行是节点，列是支路。注意limit向量对应每条支路的双向限额，即潮流允许在[-limit, limit]之间。
3.2 决策变量与完整模型
YALMIP建模的核心是sdpvar变量。相角theta、机组出力pg、线路潮流f都声明为sdpvar。线路潮流f既可以作为独立变量，也可以通过等式约束捆绑到相角上。下面代码同时给出两种写法，推荐使用变量f加等式约束的方式，这样后续约束表达更清晰。
theta = sdpvar(nb,1); pg = sdpvar(ng,1); f = sdpvar(length(limit),1); % 等式约束：线路潮流与相角关系 line_def = []; for k = 1:size(br,1) I = br(k,1); J = br(k,2); line_def = [line_def, f(k) == (theta(I)-theta(J))/xij(k)]; end % 节点功率平衡 balance = [A*f == pg - pd]; % 线路限值 line_limit = [-limit <= f <= limit]; % 机组限值 gen_limit = [pmin <= pg <= pmax]; % 参考节点相角为0 ref_con = [theta(ref) == 0]; % 目标函数 Objective = sum(c .* pg); % 合并约束 Constraints = [line_def, balance, line_limit, gen_limit, ref_con];在YALMIP中，约束之间的逗号表示逻辑“与”，即所有约束都要满足。line_def使用循环逐个定义每条线路的潮流方程，之所以不用向量化写法，是因为支路起点和终点不是连续的索引，循环可读性更好。balance写成A*f == pg - pd，这其实是nb个等式，YALMIP自动展开成向量约束。ref_con固定参考节点相角，是DC-OPF中必须的条件，否则节点平衡约束矩阵奇异。
3.3 调用CPLEX求解并提取LMP
求解设置使用sdpsettings指定CPLEX作为求解器。下面代码包含常用参数：输出详细性、MIP gap容差、求解器输出保存。对于纯LP问题，mipgap无效，但保留可以让代码在扩展到机组组合时不需要改动。
ops = sdpsettings('solver','cplex', ... 'verbose',2, ... 'cplex.mip.tolerances.mipgap',1e-4, ... 'cplex.timelimit',300); result = optimize(Constraints, Objective, ops); if result.problem == 0 Pg = value(pg); Theta = value(theta); Flow = value(f); LMP = dual(balance); % 节点边际电价，可能符号相反 disp('发电出力: '); disp(Pg); disp('相角: '); disp(Theta); disp('线路潮流: '); disp(Flow); disp('节点LMP: '); disp(LMP); else disp(result.info); endresult.problem == 0表示求解成功。dual(balance)返回每个节点平衡约束的拉格朗日乘子，这就是该节点的边际电价。需要注意，YALMIP的等式约束对偶符号可能与你记忆中的拉格朗日乘子相反。如果出现LMP为负，或明显不合理（比如成本30的节点电价反而低于成本20的节点），可以尝试对dual(balance)取负。实际项目中我一般先打印P_g和LMP一起对比，确认好符号后再封装成函数。
3.4 常用参数速查
CPLEX参数很多，但出清场景下高频使用的是下面几个：
| 参数 | 取值示例 | 作用 |
|---|---|---|
| cplex.mip.tolerances.mipgap | 1e-4 | 混合整数优化的最优间隙，越小越精确 |
| cplex.timelimit | 300 | 求解时间上限，防止死循环 |
| cplex.threads | 4 | 并行线程数，商用机建议设为物理核心数 |
| cplex.lpmethod | 0 | 0自动选择LP算法，1主单纯形，4 barrier |
| verbose | 2 | YALMIP输出级别，2显示求解日志 |
LP问题时，单纯形法对出清价格的对偶变量提取更友好，因为barrier可能返回的dual是内点解，在某些问题上需要跨平台处理。当模型规模达到数千节点时，我一般先用barrier求解，再用单纯形做一次热启动的交叉（crossover），以获得稳定的顶点对偶值。
注意：上述参数通过sdpsettings传入时，参数名必须与CPLEX官方名称一致。YALMIP不负责校验CPLEX参数，写错参数名不会报错，但参数不会被生效。这一点经常被新手踩坑。
4. 求解器调参与常见坑：CPLEX 参数、YALMIP 诊断
4.1 使用YALMIP诊断信息定位问题
YALMIP提供了一套诊断机制，优化结束后首先检查result.problem。==0是成功；==1是求到可行解但可能是次优；==2是不可行；==3是无界；==15是数值问题。自己写脚本时，应该把result.problem的判断写成一个函数，不同错误码给出不同提示。下面是一个常用的诊断代码段：
function check(result) switch result.problem case 0 fprintf('求解成功: %s\n', result.solvertime); case 1 fprintf('求解完成，但解可能不是最优\n'); case 2 fprintf('模型不可行，请检查约束和数据\n'); case 3 fprintf('模型无界，检查目标函数和变量边界\n'); otherwise fprintf('问题代码: %d, 信息: %s\n', result.problem, result.info); end end不可行问题是最常见的。出现不可行时，先检查节点功率平衡是否写错。一个常见错误是pd向量维度与pg不一致，或者负荷方向写反。另一个常见错误是线路潮流上限与电抗值严重不匹配，导致没有任何解能满足所有约束。此时可以先将线路limit改为一个很大的数（比如1e6），如果模型恢复可行，说明是阻塞约束过紧；如果依然不可行，问题出在功率平衡或机组限值。
| 错误现象 | 可能原因 | 排查方法 |
|---|---|---|
| result.problem=2 | 线路限值过紧 | 将limit调大试运行 |
| LMP出现负值 | dual符号取反 | 尝试乘以-1 |
| 求解缓慢 | 缺少参考节点约束 | 检查theta(ref) |
| 不同平台结果不一致 | 数值条件数过大 | 缩放数据 |
4.2 CPLEX参数对出清价格的影响
LP求解器的算法选择会直接影响对偶变量的质量。YALMIP默认让CPLEX自动选择LP算法，但在某些病态矩阵下，自动选择的barrier算法可能给出奇怪的对偶解。我一般固定使用单纯形法：
ops = sdpsettings('solver','cplex', ... 'cplex.lpmethod',1, ... 'cplex.simplex.display',2);lpmethod=1表示主单纯形法，对偶单纯形法是2。对于电力出清这种约束矩阵高度稀疏的问题，单纯形法速度快，并且给出的对偶乘子严格对应基顶点。如果模型是MILP，CPLEX会在根节点和每个子节点调用LP求解，此时lpmethod选项同样作用于LP松弛。在出清计算中，价格来自最后一个LP松弛的对偶，所以LP算法选择很重要。
还要注意CPLEX默认的数值精度。对于节点电价，如果某些线路参数数量级差距太大（比如电抗0.0001和潮流限值10000），会出现数值坏块。建议将标幺值统一在0.01~1范围，潮流限值统一到100的倍数。下面给出一个检查线性约束矩阵规模的技巧：
% 使用YALMIP的export导出约束矩阵 [Model, ~] = export(Constraints, Objective); condest = condest(Model.A); % 估计条件数 fprintf('系数矩阵条件数估计: %.2e\n', condest);如果condest超过1e12，需要缩放数据。缩放方法是把电抗乘以100，或者把潮流限值除以100，让系数矩阵各行量纲接近。
4.3 验证LMP是否正确：无阻塞场景
验证模型最直接的方法是构造一个无阻塞案例。将全部线路潮流上限设为1e6，此时线路约束不会起作用，各节点LMP应当相等，且等于系统边际电价。如果此时dual(balance)不相等，说明对偶符号或约束方向有问题。这里给出一个验证片段：
limit_free = ones(nl,1)*1e6; % 自由流通 ops = sdpsettings('solver','cplex','cplex.lpmethod',1); Constraints2 = [line_def, A*f == pg - pd, ... -limit_free <= f <= limit_free, ... pmin <= pg <= pmax, theta(ref)==0]; optimize(Constraints2, sum(c.*pg), ops); LMP_free = dual(balance);注意这里不能直接复用前面的Constraints，因为前面已经把limit固化进去了。建议把约束构建写成带limit参数的函数，方便在不同场景间切换。我习惯把整个出清模型封装为一个函数lmp_dcopf(pd, limit)，返回发电量、相角、LMP和潮流。这样在测试不同负荷或不同阻塞场景时，只要调用函数即可。
4.4 常见的数值陷阱
电力系统数据中，发电机爬坡率、线路电抗和负荷可能跨越多个数量级。YALMIP对线性约束的默认容差一般为1e-6，但CPLEX内部还有自身的可行容差和最优容差。当节点电价出现不明显的小数漂移或对称节点价格不对称时，可以尝试放松CPLEX的可行容差：
ops = sdpsettings('solver','cplex', ... 'cplex.mip.tolerances.integrality',1e-8, ... 'cplex.simplex.tolerances.feasibility',1e-7);但不要轻易把feasibility调大，否则得到的“可行解”可能实际上不可行，价格对偶也可能失真。另一种做法是提高数据精度，将所有输入都使用double类型，避免用单精度浮点。
5. 从出清到结算：LMP 分解与实用技巧
5.1 把LMP拆成能量价格、阻塞价格和损耗价格
实际电力市场结算时，需要告诉市场参与者当前节点价格为什么高于或低于参考节点。最常见的是将LMP分解为三部分：
LMP_i = 能量价格 + 阻塞价格 + 损耗价格在DC-OPF中忽略损耗，所以损耗价格通常通过能量价格旁边的增量损耗因子近似。更严谨的做法是用AC-OPF或者增加损耗系数。但在DC-OPF框架下，通常只分离能量和阻塞两部分：能量价格取参考节点的LMP，阻塞价格为节点电价减去能量价格。阻塞价格的本质是阻塞线路的影子价格引发的重新调度成本。
YALMIP中可以获取所有约束的对偶变量，包括线路限值约束的对偶。下面代码演示如何从影子价格计算阻塞盈余，并用结算差额验证：
% 分别定义上下界约束，方便提取乘子 line_upper = [f <= limit]; line_lower = [-limit <= f]; ops = sdpsettings('solver','cplex','cplex.lpmethod',1); optimize([line_def, balance, line_upper, line_lower, gen_limit, ref_con], ... sum(c.*pg), ops); mu_upper = dual(line_upper); mu_lower = dual(line_lower); % 方法1：从影子价格计算阻塞盈余（符号需根据目标函数方向验证） surplus_shadow = -mu_upper'*limit + mu_lower'*limit; % 方法2：从结算收入核算，方法2更直观不易出错 surplus_settle = sum(LMP .* pd) - sum(c .* Pg); disp('阻塞盈余(影子价格): '); disp(surplus_shadow); disp('阻塞盈余(结算差额): '); disp(surplus_settle);这里line_upper写成f <= limit，等价于f-limit≤0，其dual非负。line_lower写的-limit <= f，等价于-limit-f≤0，dual也非负。实际项目中，更多使用结算差额来反推影子价格是否正确，因为结算差额只依赖LMP和发电成本，不依赖约束编排顺序。
5.2 扩展多时段和机组组合
节点边际电价出清模型可以很容易扩展为多时段经济调度：把时间维度加入所有变量，将机组爬坡约束加入模型。此时决策变量变为pg(t,n)，负荷pd(t,n)，线路潮流f(t,k)。目标函数是各时段发电成本之和，爬坡约束写成相邻时段的出力差限幅。如果还要处理机组启停，就需要引入二进制变量，模型变为MILP，CPLEX依然可以高效求解。这时价格提取就要小心：机组组合的节点边际电价通常从最后一次LP松弛的对偶变量中获得，而该LP松弛是在整数变量固定之后求解的。在YALMIP中，使用optimize得到MILP解后，再对固定整数变量的原LP调用一次optimize，提取此时的对偶变量。
一个实用技巧是：把出清模型写成MATLAB函数，输入为负荷向量、线路参数和机组参数，输出为LMP和阻塞价格，并用单元测试固定几个已知场景。例如三节点系统在某一线路阻塞时，节点间LMP差值应该与该线路的阻塞影子价格相关，可以手工验证。这样后续换数据、换参数时，不会因为某个约束顺序改变而导致对偶变量错位。
最后一个小提醒：当使用dual()提取对偶时，务必在optimize成功之后立即执行，不要在其他位置调用。YALMIP不会缓存对偶值，重新计算或者修改模型后，再调用dual()会得到错误结果或直接报错。我的习惯是将问题写成函数，在函数内部完成optimize和dual提取，再把结果返回，这样能避免很多隐性bug。
本文还有配套的精品资源，点击获取
