机组运行约束的影子价格如何塑造节点边际电价？——5节点系统LMP分解实战

简介：一套面向电力市场与毕业设计场景的完整源程序包，围绕《机组运行约束对机组节点边际电价的影响分析》论文展开，适合电气工程、电力经济方向的高年级本科生与研究生使用。内容基于节点电价出清全时段优化模型，推导了系统备用约束下机组节点边际电价与机组运行约束影子价格、备用约束影子价格的关系式，系统分析机组功率上下限、爬坡滑坡约束对节点电价以及电能价格、阻塞价格分量的影响，并采用5节点算例验证所提方法与传统节点电价计算结果的一致性。压缩包共6个文件，约333KB，含m源程序、lp优化模型、xlsx算例数据、docx分析报告及2张结果对比图，文件结构清晰，可直接运行复现算例结果，便于对照论文学习。已有163人学习，适合想深入理解节点电价形成机制、评估机组物理参数对收益影响的研究者。

1. 相同的报量、不同的出清价：机组运行约束的影子价格才是 LMP 的幕后推手

在电力市场出清里，同样的 10 MW 申报，放在不同时段、不同机组上，最终节点边际电价（LMP）可能差出几十元/MWh。很多人只盯负荷和阻塞，却忽略机组自身的功率上下限、爬坡滑坡约束——这些不等式约束一旦激活，会通过对偶变量直接写入电价。这份源码对应知网论文《机组运行约束对机组节点边际电价的影响分析》，把备用约束下的 LMP 关系式完整推导出来，并用 5 节点系统做了全时段出清验证。对做电力市场结算、写毕业设计或者研究阻塞费用分摊的人来说，这是一套可以直接复算的基准程序，尤其适合用来验证自己对影子价格方向的理解是否正确。

2. 机组运行约束的数学刻画：从功率上下限到爬坡滑坡的影子价格

2.1 出清模型的对偶形式：LMP 是如何从约束里“挤”出来的

全时段安全约束经济调度（SCED）写成线性规划后，节点电价并非直接来自目标函数，而是来自等式约束的拉格朗日乘子。离散化到若干个时段 T，目标函数通常取发电成本最小：

min Σ_t Σ_i c_i · p_{i,t} + 启动/空载常数项

约束包括：节点功率平衡、机组出力上下限、爬坡滑坡约束、线路潮流限值、系统备用约束。把不等式约束写成标准形 A·x ≤ b 后，对偶变量便是资源的边际价值。比如机组 i 在时段 t 被顶上最大出力，其上限约束对应的影子价格为 μ_{i,t}^max。这个影子价格的含义是：如果该机组再多给 1 MW 上限，系统总购电成本能下降多少。节点边际电价不是某个单一乘子，而是这些影子价格通过机组对节点注入的灵敏度折算回节点的结果。

我在源码包中看到的模型把各时段叠加成一个大的稀疏线性规划，而不是逐时段独立出清。这个处理很关键，因为爬坡约束天然跨时段。如果只在单个时段内做经济调度，爬坡滑坡约束根本无法进入约束矩阵，影子价格恒为零。论文中推导的关系式可以写成如下形式：

λ_{n,t} = 系统电能价格分量 + 阻塞价格分量 + 机组运行约束影子价格折算项

2.2 备用约束下的 LMP 与影子价格关系式

考虑系统备用后，模型会多出“备用容量之和 ≥ 负荷 + 备用需求”的约束。备用约束的影子价格是全局性的，它会通过拉格朗日乘子抬升所有节点的电能价格分量，而不是只影响某个节点。机组功率上限约束和下限约束的影子价格则带有明显机组属性，需要通过注入转移因子映射到节点价格。

设机组 g 对节点 n 的注入转移因子为 GSF_{g,n}，线路 l 对节点 n 的潮流灵敏度为 PTDF_{l,n}，则节点 n 在时段 t 的 LMP 可以分解为：

λ_{n,t} = λ_t^E + Σ_l μ_{l,t}^line · PTDF_{l,n} + Σ_g ( μ_{g,t}^up - μ_{g,t}^down ) · GSF_{g,n}

其中 λ_t^E 是系统功率平衡约束的乘子，μ_{l,t}^line 是线路阻塞影子价格，μ_{g,t}^up 和 μ_{g,t}^down 分别是机组功率上、下限约束的影子价格。系统备用约束的影子价格则并入 λ_t^E 中参与全局价格形成。这个式子看起来不复杂，但实际使用时最容易踩坑的是符号方向：linprog 默认不等式约束为 A·x ≤ b，其拉格朗日乘子非负；如果模型里把功率下限写成 -p ≤ -p_min，乘子符号就要跟着翻转。论文附录里的推导对这部分做了详细展开，程序包里的 lmp.m 也正是按这个方向实现的。

2.3 爬坡约束的全时段建模要小心“相邻时段”的对偶成对出现

爬坡约束写成差分形式：

-RD_g ≤ p_{g,t} - p_{g,t-1} ≤ RU_g

在矩阵里，每个时段对每一台机构造两行不等式行。注意变量顺序必须统一：如果决策变量 x 按“时段主序、机组次序”展开，那么差分矩阵的行要先遍历机组再遍历时段，或者相反，这决定了对偶变量如何回读。源码中的 dlsc5.lp 文件把约束命名为 ramp_up_g1_t2 这类清晰的名字，读 LP 文件时可以直接按名字定位。

爬坡影子价格成对出现：通常只有“爬坡上限”或“滑坡上限”其中一个乘子非零，两个同时非零意味着目标函数出现退化。实践中我发现很多自编程序在提取影子价格时只取了不等式行的前半段，后半段滑坡约束的乘子被漏掉，导致 LMP 分解结果不闭合。验证方法是把所有约束的乘子带回到 KKT 条件，看目标函数对每个决策变量的梯度是否为零；如果残差大于 1e-6，基本就是乘子映射错位。

3. 5节点算例数据准备与 LMP 计算代码实现

3.1 读懂 dlsc5.lp 与 5节点数据.xlsx 的结构

程序包里的 5节点数据.xlsx 是原始输入，dlsc5.lp 是写好的线性规划文本，两者约束顺序必须一一对应。拿到手我先看 Excel 里的表名和列名，因为它直接决定后续矩阵拼接的索引。典型结构如下：

全时段出清时，负荷曲线会按 T 个时段展开，每台机组在每个时段都对应一个出力变量。Excel 里的固定参数如 pmin、pmax 一般不随时间变化，但在模型里要扩展成 T 段的列向量。我一般会先写一个load_case5.m把 Excel 读进来，再构造稀疏矩阵，避免在循环里频繁调用 xlsread，否则 24 时段、50 台机组的数据拼装能跑出 5 分钟以上的性能问题。

3.2 用 MATLAB linprog 做全时段出清并提取对偶变量

我处理这套算例时用的是 MATLAB 的 linprog，求解器选项选 dual-simplex。之所以不用默认的 interior-point，是因为对偶单纯形法在回读影子价格时更稳定，尤其是面对退化约束时不会随意换基导致影子价格跳变。下面是可运行的最小化代码骨架：

% lmp_full_dispatch.m % T 个时段, G 台机组, N 个节点 % 决策变量 x 按 时段主序、机组次序 展开 ngen = size(gen, 1); T = 24; x0 = zeros(T * ngen, 1); % 目标函数：最小化能量成本 % gen.rho 是每台机组的边际成本列向量, 单位元/MWh f = repmat(gen.rho, T, 1); % 节点功率平衡: Aeq * x = beq % Aeq 的行数是 N * T, 每个时段一个节点功率平衡行 % 机组 g 在节点 bus(g), 则 Aeq(bus(g) + (t-1)*N, (t-1)*ngen + g) = 1 [Aeq, beq] = build_balance(gen, load, T, N); % 爬坡约束: -RD <= x(t,i) - x(t-1,i) <= RU % 每个时段每个机组生成两行不等式 [Ain, bin] = build_ramp(gen, T); % 出力上下限 lb = repmat(gen.pmin, T, 1); ub = repmat(gen.pmax, T, 1); % 线性规划求解 opts = optimoptions('linprog', 'Algorithm', 'dual-simplex', 'Display', 'iter'); [x, fval, exitflag, output] = linprog(f, Ain, bin, Aeq, beq, lb, ub, opts); % 提取对偶变量 dual_eqlin = output.lambda.eqlin; % 节点功率平衡乘子 dual_ineqlin = output.lambda.ineqlin; % 爬坡与备用约束乘子 dual_lower = output.lambda.lower; % 下限约束影子价格 dual_upper = output.lambda.upper; % 上限约束影子价格

dual_eqlin的长度应该是 N × T，前 N 个值对应第 1 时段各节点的功率平衡乘子。如果不计网损，这些值就是该时段各节点的 LMP 雏形。dual_lower和dual_upper与 lb、ub 对应，它们正是 2.2 节公式中的 μ_{g,t}^down 和 μ_{g,t}^up。爬坡约束的乘子则需要从dual_ineqlin中按build_ramp的行顺序切分；我通常会在build_ramp里保存一个constraint_names元胞数组，把每个约束对应的机组、时段和类型记录清楚，省去后面逐个猜。

3.3 用 LP 文件核对约束顺序

dual_ineqlin的长度只告诉你不等式约束的总数，不会告诉你每个数对应哪个约束。最可靠的方式是直接读 dlsc5.lp 文件，按约束名和顺序对齐。用命令行最快：

grep -n "^c" dlsc5.lp | head -40

LP 文件里每个约束行以c开头，约束名一般写在行首。比如ramp_up_g2_t5表示第 2 台机组第 5 时段的爬坡上限约束。这个名称顺序与求解器输出的对偶变量顺序一致。将grep结果保存成文本，再在 MATLAB 里读入，映射到dual_ineqlin对应位置：

% map_dual_ineq.m fid = fopen('constraint_order.txt', 'r'); names = textscan(fid, '%s', 'Delimiter', '\n'); names = names{1}; fclose(fid); for k = 1:length(names) if contains(names{k}, 'ramp_up') fprintf('索引 %d: %s, 影子价格 = %.4f\n', ... k, names{k}, dual_ineqlin(k)); end end

这里需要严格保证constraint_order.txt的行顺序与求解器内部约束顺序一致。建议在生成Ain时同步把每个不等式的名称写进日志文件，而不是事后从 LP 文件猜。很多“LMP 分解对不上”的问题，最后查下来都是这一行顺序错位。

4. 全时段出清结果分解：电能价格与阻塞价格分量怎么算

4.1 从对偶变量重构 LMP 三分量

拿到dual_eqlin之后，节点电价还只是功率平衡乘子，没有包含机组上下限约束和阻塞约束的贡献。需要把 2.2 节的公式落成代码。先把线路潮流约束的乘子提取出来，再计算节点的 PTDF 矩阵。对于 5 节点系统，PTDF 可以基于电抗矩阵手算，也可以直接用 MATLAB 的自带函数：

% decompose_lmp.m % base_index 为参考节点编号 PTDF = compute_ptdf(branch, base_index); % 节点数 × 线路数 % 线路约束影子价格: 剔除备用约束后的 remaining 行 mu_line = dual_ineqlin(line_ineq_index); % 阻塞价格分量 congestion_component = PTDF * mu_line; % 机组上下限约束影子价格映射到节点 % gen_matrix 维度 N × G, gen_matrix(n,g)=1 表示机组 g 在节点 n umap = gen_matrix * (dual_upper - dual_lower) / gen_matrix_sample_time; % 系统电能价格分量取参考节点平衡乘子 energy_component = dual_eqlin(base_index); % 最终 LMP lmp = energy_component + congestion_component + umap;

代码里的umap是把每台机组的上下限影子价格按其所在节点做聚合，再换算成节点增量。这里要注意单位：如果模型里出力单位是 p.u.，影子价格单位是元/p.u.，需要乘以 baseMVA 回到元/MWh。我在第一次复算时就是忘了这个换算，导致所有机组约束的影子价格都比论文结果大了 100 倍，排查了很久才定位到。

4.2 多案例对比：哪些约束在哪些时段被激活

论文里做了多案例分析，我按同样的思路把约束分了三类：仅功率上下限、上下限加爬坡、上下限加爬坡加备用。下面摘取 t=10 时段的一组典型影子价格，用于说明数量级和方向：

从表中可以看到，功率上限的影子价格出现在 G1 这类低成本机组上，因为低成本机组被顶到上限后，系统不得不调用更贵的机组，影子价格为正。下限约束影子价格为负，含义是“如果允许 G3 进一步压低出力，总购电成本能下降”。爬坡约束影子价格集中在负荷快速上升的黄昏时段，这与实际系统晚高峰爬坡紧张的情况一致。备用约束影子价格对所有节点一视同仁，它不改变节点间相对价格，只会平移整个价格平面。

4.3 与传统单时段出清方法的一致性验证

论文中的计算方法要和传统计算方法做交叉验证。传统方法不考虑爬坡约束，逐时段独立出清。全时段模型的正确性验证可以这样做：把爬坡上下限放大到 1e6，此时爬坡约束退化为恒满足，对偶乘子应为零，LMP 应与单时段结果一致：

% verify_against_single_period.m % 将爬坡限值放大, 重新出清 gen.RU = 1e6; gen.RD = 1e6; [x_free, ~, ~, output_free] = solve_full_dispatch(gen, load, branch); % 与单时段无爬坡结果做差 lmp_single = solve_single_period(gen, load, branch); lmp_full_free = reshape(output_free.lambda.eqlin, N, T); max_diff = max(abs(lmp_full_free(:) - lmp_single(:))); fprintf('最大偏差: %.6e\n', max_diff);

如果max_diff在 1e-6 到 1e-8 之间，说明全时段模型在无爬坡时能退回到单时段结果，矩阵拼接没有错误。如果偏差很大，优先检查爬坡约束的差分方向：是t时段减t-1时段还是反过来。方向写反会导致价格严重错位，偏差可能达到几十元/MWh。这套验证步骤在拿到任意算例后都值得先跑一遍，再去看价格分解。

5. 进阶用法：用影子价格定位“卡脖子”机组并评估收益改进空间

5.1 从对偶变量中快速找出激活约束

全时段模型的对偶向量通常很长，但大量元素接近零。直接用容差过滤，再映射回约束名：

% find_active_constraints.m tol = 1e-6; active_upper = find(dual_upper > tol); active_lower = find(dual_lower < -tol); for i = 1:length(active_upper) idx = active_upper(i); t = floor((idx - 1) / ngen) + 1; g = idx - (t - 1) * ngen; fprintf('时段 %d, 机组 %d: 出力上限激活, 影子价格 %.4f\n', ... t, g, dual_upper(idx)); end

这里我用floor把一维索引拆回时段和机组编号，前提是决策变量严格按“时段主序、机组次序”展开。如果源码里的变量顺序倒过来，这段映射就要改成mod方式。筛选active_lower时用< -tol而不是> tol，因为下限约束在标准形里通常被改写成-p ≤ -p_min，乘子符号为负。建议在程序里把正负约定写成注释，否则隔几个月回来看很容易混淆。

5.2 用影子价格做灵活性改造的投资估算

机组爬坡约束影子价格可以直接用来估算爬坡能力提升的收益。例如某个时段的爬坡约束影子价格为 8.3 元/MWh，意味着爬坡能力每增加 1 MWh/时段，该时段系统成本下降约 8.3 元。我把这个值乘以改造后的爬坡增量，再乘以高频使用时段数，可以快速判断投资价值：

# sensitivity.py shadow_price = 8.3 # 元/MWh, 爬坡约束影子价格 delta_ramp = 10 # MWh/h, 计划增加的爬坡能力 hours_per_month = 2 * 30 # 晚高峰每天 2 小时 benefit = shadow_price * delta_ramp * hours_per_month print(f"月收益增量: {benefit:.0f} 元")

这种估算的上限很明确：一旦爬坡能力真正增加，约束不再激活，影子价格会降到零，实际收益远小于初始值。所以更严谨的做法是把改造后的爬坡参数重新代入模型求解，用新的 LMP 与旧 LMP 的差值做收益分析，而不是直接用原影子价格外推。我在做市场成员侧评估时，通常会把这两种方法的结果都放出来：用影子价格快速筛选值得改造的机组，再对候选机组做全时段重算，避免高估灵活性改造价值。

本文还有配套的精品资源，点击获取