简介电力市场节点边际电价出清优化方案围绕史新红《机组运行约束对机组节点边际电价的影响分析》进行复现以单时段出清模型为核心在忽略爬坡约束的前提下完整考虑机组运行约束适合电力市场、电力经济相关专业学生及研究者用于学习节点边际电价LMP与系统边际电价的形成机理。采用KKT条件将对偶问题转化为矩阵形式求解拉格朗日乘子影子价格程序结构清晰、注释完整便于对照报告理解影子价格在出清优化中的经济含义。压缩包共5个文件包含3个Matlab程序m、1篇参考文献caj和1份完整报告doc总大小仅383KB轻量易下载。目前已有4668人学习资源内附出清过程详细报告并针对运行问题提供答疑支持可帮助读者快速理解边际电价出清的规范做法并复现算例。1. 打开这个zip之前先弄懂节点边际电价为什么是一堆影子价格我做电力市场出清仿真这几年最常被问的一句话是节点边际电价LMP到底是怎么算出来的很多人第一反应是“用报价加权平均”或者“按负荷比例分摊成本”这两个直觉在单节点系统里能凑合一旦把电网拓扑和线路阻塞放进来误差能到几十美元。真实情况是LMP不是市场规则里拍脑袋定的而是安全约束经济调度SCED这个线性规划问题里功率平衡约束和线路阻塞约束的拉格朗日乘子组合出来的结果。换句话说价格不是算出来的是从优化问题的对偶问题里“问”出来的。这个名为“基于yalmipcplex的电力市场-节点边际电价出清优化”的zip本质上是把这条链路完整打通用YALMIP在MATLAB里把出清模型写成约束和目标交给CPLEX求解再从求解结果里抽出每个节点的边际电价。它适合三类人正在做电力市场仿真的研究生、售电或发电侧做策略的算法工程师、想从通用优化转行电力调度的开发者。读完你不仅能跑通一个三节点系统的出清还能自己拆LMP、判断结果对不对、避开一堆让人挠头的坑。2. SCED与LMP直流潮流下的安全约束经济调度怎么定价2.1 没有线路阻塞时电网其实是一个“大池子”先看最简单的情况假设整个电网是一个点所有发电机和负荷接在同一个母线上不考虑线路容量。那么出清问题就是在满足总负荷的前提下让发电总成本最小。这是一个再简单不过的线性规划目标是最小化所有机组的发电成本之和约束只有一个——总出力等于总负荷。这个问题的对偶乘子拉格朗日乘子就是统一出清价所有节点完全一样。为什么因为在这个假设下不管电从哪个机组来、送到哪个节点去成本只跟机组自身的报价有关跟“位置”无关。这就是电力市场里常说的“全网一口价”。但现实电网不是这样。输电线路有容量上限电能必须遵循基尔霍夫定律流动。于是同样1MW负荷加在送端节点和加在受端节点引起的发电成本增量完全不同。这个“不同”就是节点边际电价存在的原因。我一般会给新人打一个比方LMP是什么是“在这个节点多接1MW负荷整个系统最优发电成本会增加多少钱”。这个增量包含了机组边际报价也包含了线路阻塞带来的机会成本。2.2 线路越过容量价格就开始“撕裂”用一个最小例子说明。两个节点一条线路连着。送端节点1有一台便宜机组报价10美元/MWh受端节点2有一台贵机组报价20美元/MWh负荷全在节点2。线路容量很大时便宜机组直接送电全网价格10美元。现在把线路容量调小比如只能送60MW而负荷是100MW。便宜机组最多往节点2送60MW剩下的40MW必须让节点2的贵机组发。这时候全网不再是统一价。节点1的LMP是10美元——因为在节点1多1MW负荷本地便宜机组直接多发1MW就行。节点2的LMP变成20美元——因为在节点2多1MW负荷线路已经满了只能让贵机组加出力。两个节点的价差正好等于阻塞线路的“稀缺租金”。这个租金不是任何人主观定的是线路潮流约束的对偶乘子。安全约束经济调度SCED的核心就是把这些约束全部放进优化问题里让价格自己“长出来”。2.3 为什么偏偏是YALMIP CPLEX理论上你自己用MATLAB写一个单纯形法或者调linprog也能解这个LP。但出清模型有个特殊要求你得知道每个约束对应的影子价格而且要把影子价格准确地映射到“功率平衡乘子”“线路阻塞乘子”这些业务概念上。自己写求解器变量和约束多了以后对偶信息的跟踪极其容易出错。YALMIP的价值在于建模层和求解层分离。你只需要声明变量是sdpvar然后像写数学公式一样写约束YALMIP负责把模型转换成标准LP/QP形式传给底层求解器。CPLEX则在求解层提供了成熟的单纯形和对偶单纯形算法数值稳定大规模稀疏问题也能扛。更关键的是YALMIP的dual函数能直接取回约束对应的对偶乘子这比从CPLEX原生接口里翻对偶解要友好得多。我见过不少直接从CPLEX Java/C接口取乘子的工程代码里全是矩阵索引和符号转换维护成本极高。另外YALMIP的好处是求解器可替换。同一个人样模型把sdpsettings里的solver改成其他商业求解器或开源求解器也能跑只是性能和对偶信息获取的方便程度不同。我一般在教学和验证阶段固定用CPLEX因为它的LP对偶解默认返回完整且预求解不会轻易“吃掉”你关注的约束乘子。3. YALMIP写三节点出清建模、求解和第一个LMP3.1 直流潮流的近似逻辑与适用边界市场出清仿真里最常用的潮流模型是直流潮流DCOPF不是交流潮流。交流潮流要同时算有功、无功、电压幅值和相角问题是非凸的求解慢、可能不收敛而且LMP的定义在交流模型下要复杂得多——无功价格、电压约束的影子价格都会掺进来。直流潮流做的事情很粗暴忽略无功、忽略网损、假设电压幅值都是1只保留有功和相角的关系。这个关系写出来很简单线路有功潮流P_ij约等于两端相角差除以电抗也就是P_ij(θ_i-θ_j)/x_ij。x是线路电抗单位用标幺值。这个线性关系让整个出清问题变成一个纯线性规划凸、可收敛、对偶乘子有清晰的价格含义。代价是精度——直流潮流不适用于电压问题、无功优化场景但在节点边际电价仿真、市场力分析、阻塞管理这些场景里它是绝对的主流选择。你去看各地的现货市场出清系统核心引擎基本都是某种形式的DCOPF或加入了网损近似的扩展版。需要注意一个边界x_ij不能为0。电抗为零意味线路导纳无穷大B矩阵奇异潮流方程无解。我见过有人在测试数据里为了省事把所有线路电抗写成0结果YALMIP一跑就报“NaN”或者不可行最后查了半天发现是电抗单位写错了。3.2 三节点系统的参数与完整求解代码为了把原理落到代码上我构造一个最简单但能体现阻塞定价的三节点链式系统。节点1是送端节点2是中转节点3是负荷中心。节点1有一台便宜机组节点2有一台贵机组负荷100MW全部在节点3。线路1-2容量设为60MW线路2-3容量设为150MW。这样一定会出现阻塞——便宜机组的电最多送60MW到节点2节点3的缺口只能让节点2的贵机组顶上。系统参数如下对象参数数值线路1-2电抗(p.u.) / 容量(MW)0.1 / 60线路2-3电抗(p.u.) / 容量(MW)0.1 / 150机组G1所在节点 / 报价($/MWh) / 出力上下限(MW)1 / 10 / 0~100机组G2所在节点 / 报价($/MWh) / 出力上下限(MW)2 / 20 / 0~100负荷D3所在节点 / 大小(MW)3 / 100下面是完整的YALMIP CPLEX求解代码。这段代码我建议你原样跑通一次再逐步改参数观察价格变化。%% 三节点直流潮流出清YALMIP CPLEX % 节点1送端2中转3受端(负荷) % 线路L12连接1-2L23连接2-3 clear; clc; %% 系统数据 mpc.nb 3; % 节点数 mpc.nl 2; % 线路数 % 线路: [起点, 终点, 电抗(p.u.), 容量(MW)] mpc.line [1 2 0.1 60; 2 3 0.1 150]; % 机组: [节点, 报价($/MWh), 出力下限(MW), 出力上限(MW)] mpc.gen [1 10 0 100; 2 20 0 100]; mpc.load [0 0 100]; % 节点3负荷100MW %% 构建直流潮流B矩阵和线路-节点关联矩阵Bf B zeros(mpc.nb, mpc.nb); for k 1:mpc.nl f mpc.line(k,1); t mpc.line(k,2); x mpc.line(k,3); B(f,f) B(f,f) 1/x; B(f,t) B(f,t) - 1/x; B(t,f) B(t,f) - 1/x; B(t,t) B(t,t) 1/x; end Bf zeros(mpc.nl, mpc.nb); for k 1:mpc.nl f mpc.line(k,1); t mpc.line(k,2); x mpc.line(k,3); Bf(k,f) 1/x; Bf(k,t) -1/x; end ref 3; % 参考节点选3相角置0 %% 定义优化变量 P sdpvar(2, 1); % 两台机组出力 theta sdpvar(3, 1); % 节点相角 %% 约束 cons_plim (mpc.gen(:,3) P mpc.gen(:,4)); % 机组出力上下限 cons_ref (theta(ref) 0); % 参考节点相角为0 F Bf * theta; % 线路有功潮流表达式 cons_fline (F mpc.line(:,4)); % 线路潮流上限 cons_fdown (-mpc.line(:,4) F); % 线路潮流下限 cons_balance (sum(P) sum(mpc.load)); % 功率平衡 % 合并所有约束传给求解器 all_cons [cons_plim, cons_ref, cons_fline, cons_fdown, cons_balance]; %% 目标函数与求解 obj mpc.gen(:,2) * P; % 总发电成本线性成本 options sdpsettings(solver, cplex, verbose, 1); sol optimize(all_cons, obj, options); %% 读取结果 P_opt value(P); F_opt value(F); theta_opt value(theta); fprintf(G1出力: %.2f MW, G2出力: %.2f MW\n, P_opt(1), P_opt(2)); fprintf(线路L12潮流: %.2f MW, 线路L23潮流: %.2f MW\n, F_opt(1), F_opt(2)); fprintf(总发电成本: %.2f 美元\n, value(obj));逻辑上说得很直白先构建B矩阵和Bf矩阵B矩阵用来算节点相角和功率注入的关系Bf矩阵用来从相角直接得到每条线路的潮流。然后声明两台机组的出力P和三个节点的相角theta为优化变量。约束分四组机组出力上下限、参考节点相角、线路潮流上下限、功率平衡。目标函数是线性成本之和。参数说明要强调几个地方。第一参考节点必须固定否则所有相角同时加一个常数潮流不变但矩阵求解不唯一YALMIP不会报错但结果看着像随机数。第二F Bf * theta不是约束是表达式它把相角和潮流联系起来后续的线路容量约束都作用在这个表达式的值上。第三sdpsettings里solver指定用cplexverbose1让求解器打印迭代日志如果CPLEX没有正确安装或授权YALMIP会报错这时可以把solver改成其他支持的求解器先验证模型逻辑。跑完这段代码你会得到G1出力60MW、G2出力40MW、总成本1400美元。这个结果本身就是阻塞定价的起点——便宜机组只发到线路容量允许的上限剩下的全靠贵机组。3.3 第一次看结果哪些约束在“卡脖子”求解完成后不要只看出力数字要看哪条约束是积极约束active constraint也就是dual值不等于0的约束。积极约束就是“卡脖子”的位置线路L12的潮流上限约束和功率平衡约束都会有非零对偶乘子。机组G1的出力上限60MW不会触发它的上限是100MW触发的是线路容量60MW。我习惯把结果分成两种场景做对比。场景A线路L12容量60MW价格分层场景B把线路L12容量改成150MW也就是不阻塞。两张场景的LMP对比如下节点场景AL12容量60MW场景BL12容量150MW节点110 美元/MWh10 美元/MWh节点220 美元/MWh10 美元/MWh节点320 美元/MWh10 美元/MWh机组出力G160, G240G1100, G20场景B里全网一口价10美元跟2.1节的理论预测完全一致。场景A里节点1和节点2/3差了10美元差的这10美元就是线路L12阻塞的影子价格。这个“一堵就分层”的现象是所有节点边际电价仿真的核心。后面第四章我要讲的LMP分解就是把这10美元从数学上拆出来。4. LMP的提取与分解从对偶乘子到节点电价4.1 取对偶乘子YALMIP的dual()与符号约定YALMIP里取对偶乘子非常简单直接调用dual(约束句柄)。但要命的是符号约定对于形如“表达式 0”的约束YALMIP返回的乘子是非负的对于等式约束乘子可以是任意实数符号由等式方向决定。我初学时在这里栽过跟头——把功率平衡约束的dual取出来直接当成LMP结果算出来的价格全反了。标准做法是把功率平衡等式约束写成sum(P) sum(mpc.load)然后dual(cons_balance)就是能量价格基础分量λ。线路潮流约束F Fmax返回的是非负乘子μ_up-Fmax F返回的也是非负乘子μ_dn。不需要手动取负号只需要在汇总公式里按KKT条件的方向组合。4.2 拆LMP能量分量、阻塞分量与参考节点有了对偶乘子还需要一个矩阵节点转移因子PTDF。PTDF的含义是某个节点注入1MW功率参考节点吸收每条线路上会增加多少潮流。直流潮流下PTDF矩阵可以由Bf和B矩阵直接算出来。下面这段代码接着上一节的系统数据计算PTDF并把三个节点的LMP拆成能量分量和阻塞分量。%% 计算PTDF矩阵 nonref [1 2]; % 非参考节点参考节点3的列保持0 Bred B(nonref, nonref); % 去掉参考节点行列的B矩阵 PTDF zeros(mpc.nl, mpc.nb); PTDF(:, nonref) Bf(:, nonref) / Bred; % 线路注入灵敏度参考节点列为0 %% 取对偶乘子 lambda dual(cons_balance); % 功率平衡乘子能量价格基础 mu_up dual(cons_fline); % 线路FFmax的乘子非负 mu_dn dual(cons_fdown); % 线路-FmaxF的乘子非负 %% 组合LMP LMP zeros(mpc.nb, 1); for i 1:mpc.nb LMP(i) lambda - sum((mu_up - mu_dn) .* PTDF(:, i)); end %% 打印结果 fprintf(能量分量(全网统一): %.2f $/MWh\n, lambda); for i 1:mpc.nb blocking lambda - LMP(i); % 阻塞分量含符号 fprintf(节点%d LMP%.2f $/MWh, 阻塞分量%.2f $/MWh\n, ... i, LMP(i), blocking); end这里PTDF(:,i)是一个列向量表示节点i的注入对各条线路潮流的灵敏度。用mu_up减去mu_dn再和PTDF逐元素相乘后求和得到该节点阻塞分量的修正值。公式为什么是“减去”而不是“加上”回到第二章的双节点例子送端节点1的LMP10全网能量分量为20阻塞修正值是-10方向必须能让送端更便宜、受端更贵。如果你拿到的数据算出来方向相反优先检查PTDF的参考节点是否选对了以及线路约束写成“F Fmax”时YALMIP返回的dual是否真的非负。在三节点例子里这个代码会输出lambda20美元节点1的LMP10美元阻塞分量-10节点2和节点3的LMP20美元阻塞分量0。这跟3.3节表格完全一致。阻塞分量只作用于阻塞线路的“送端”这是直流潮流的典型特征。注意阻塞盈余的概念也由此而来节点3的用户按20美元买电节点1的机组按10美元卖电中间每度电的10美元差额就是阻塞盈余在真实市场里会通过阻塞费用结算机制分配。关于参考节点再补充一句LMP的数值总量不依赖参考节点的选择三个节点算出来都一致。但分解出来的能量分量和阻塞分量会随参考节点变化——参考节点的阻塞分量恒为0能量分量会被“吸收”一部分。所以做分解时不要在报告里只说“能量价是多少”一定要注明参考节点是哪个。4.3 扰动校核用电价的定义反推你的结果这是我最推荐的验证习惯。LMP的严格定义是“某节点增加单位负荷引起的最优成本增量”。与其信公式不如直接让计算机“扰动”一下。做法是把求解过程封装成一个函数然后对每个节点做一次1MW负荷扰动重新求解比较成本差。function [P, F, obj_val] solve_sced(mpc) % 把第三章的建模和求解包成函数返回出力、潮流和目标值 % mpc结构里有line/gen/load/nb/nl等字段 % 具体建模代码同第三章略 end %% 扰动校核对每个节点加1MW负荷求成本增量 c0 0; % 先求基础场景的总成本 [~, ~, c0] solve_sced(mpc); for i 1:mpc.nb mpc2 mpc; mpc2.load(i) mpc2.load(i) 1; % 只在节点i加1MW [~, ~, c1] solve_sced(mpc2); fprintf(节点%d扰动法LMP%.2f $/MWh\n, i, c1 - c0); end在三节点例子里这个扰动校核的结果会和4.2节的LMP完全吻合节点1加1MW负荷成本增加10美元节点2和节点3加1MW负荷成本增加20美元。逻辑也很好验证——节点1的负荷增加不会加重线路L12的负担本地负荷吸收本地出力便宜机组直接加发节点2和节点3的负荷增加则必须靠贵机组。扰动法看着笨但它是唯一不考虑对偶理论也能解释得通的价格验证手段。我每次写新模型都会保留这个校核脚本作为LMP结果的“金标准”。5. 出清仿真的五类翻车现场与排查路径5.1 INFEASIBLE电抗为零、参考缺失与负荷不平衡现象optimize返回求解状态为不可行或者YALMIP直接报NaN。原因有三大类按我实际遇到频次排序第一线路电抗x写成0B矩阵奇异第二参考节点相角约束没有加模型自由度过大导致数值病态第三总发电出力上限小于总负荷数学上无解。还有一类隐蔽的负荷数据写在某个节点但该节点没有发电机也没有外来线路形成孤岛。解决先查B矩阵的行列式行列式接近0立刻定位电抗异常检查所有节点是否至少有发电机或线路连通到参考节点。快速定位不可行约束的办法是给功率平衡加一个人工松弛变量比如把平衡约束改成sum(P)slacksum(load)让slack无界求解后看slack落在哪个节点就知道哪里功率缺了。5.2 dual()全为零对偶信息被预求解消掉了现象模型求解成功出力、潮流都对但dual(cons_balance)返回0。LMP算出来全是0怎么看都不对。原因CPLEX预求解阶段可能判定某条约束冗余并直接删掉对偶乘子随之丢失。尤其当功率平衡约束的表达式写法不够“标准”时容易触发。另外YALMIP里多个约束合并成复合约束后dual返回的向量顺序和约束顺序如果不匹配也容易让人误读成零。解决把需要取乘子的约束单独定义成独立的句柄不要和其他约束揉在一起。求解前在sdpsettings里关闭部分预求解options.cplex.preprocessing.reduce 0。如果还是取不到用扰动法4.3节交叉验证——不需要信dual也能得到价格。这是最保险的兜底手段。5.3 全网一口价是不是哪里没设容量现象明明加了线路容量约束三节点系统还是全网统一价跟第二章“大池子”一样。原因线路容量设得太大了。容量上限在200MW线路实际潮流才100MW约束虽然存在但从未被激活乘子当然为0价格不会分层。很多人第一次跑通代码后发现价格都一样以为自己优化错了其实是数据没卡到阻塞点。解决逐步缩小关键线路容量观察LMP分叉的临界点。这个临界值其实有物理含义线路容量紧到一定程度市场上的便宜机组送不过去受端价格开始爬升。我在3.3节里把L12容量从150降到60就是故意制造这个分叉。调参时每次都对比“潮流值vs容量上限”两个数字相等的那条线就是当前最拥堵的线。5.4 LMP是负的或高于机组报价稀缺定价不是bug现象某节点LMP出现负值或者高于全网最高机组报价20美元达到35美元。原因这是线性规划里约束影子价格的正常表现不是求解器抽风。负价格通常出现在机组有最小技术出力比如最小出力50MW但负荷只有30MW或者必须开机维持电压支撑的场景价格高于最高报价则说明某条受约束的线路把稀缺性放大受端愿意付出的边际代价远超本地最贵机组的报价。解决不要“修”结果先理解是哪条约束在定价。把该节点的扰动法LMP手动算一遍你会发现它和影子价格完全一致。如果业务上不能接受负价格那是市场规则需要加价格下限不是优化模型的问题。仿真时保留这个现象即可。5.5 多解与退化同一份数据不同的价格现象代码没变、数据没变不同机器上跑出来G1/G2出力的分配有差异甚至LMP在某些节点上有微小波动。原因LP存在多个最优解尤其是机组报价都是线性没有二次项时成本曲线是一条直线多个出力组合可以达到相同的总成本。单纯形法会沿着顶点走起点和预求解差异会让它落在不同的最优顶点上。价格本身不一定变但个别节点的阻塞分量可能因为参考约束的取舍而不同。解决给目标函数加一个极小的二次项比如0.0001×P^2让问题变成凸二次规划最优解唯一。这在实际现货市场出清里也是常用手段——报价段内部按线性处理但出清程序里会加微小的平滑项避免解漂移。代价是价格含义从纯线性变成带微弱二次修正但量级完全可以忽略。6. 从单时段到多时段SCED扩展成时序价格的进阶技巧6.1 给模型加时间轴爬坡约束与24小时价格单时段模型解决好“某一时刻的电价”但现货市场价格是时序的调度员关心的是机组能不能在15分钟或1小时内从50MW爬升到80MW。给模型加时间轴的套路是固定的所有变量加一个时间维度P变成2×T的矩阵theta变成3×T的矩阵功率平衡、线路容量约束对每个时段单独写目标函数变成所有时段的成本之和。额外要加的约束是爬坡约束T 24; % 24个时段 P sdpvar(2, T); % 每台机组每个时段的出力 R_up 30; % 每时段最大爬坡上限MW R_dn 30; % 每时段最大爬坡下限MW cons_ramp []; for t 2:T cons_ramp [cons_ramp, P(:,t) - P(:,t-1) R_up]; cons_ramp [cons_ramp, P(:,t-1) - P(:,t) R_dn]; end这里P(:,t)-P(:,t-1)用矩阵形式一次性约束两台机组比写两个标量约束干净。爬坡约束会对价格产生实质影响低谷时段便宜机组被强制压低出力高峰时段又爬不上去受端节点在爬坡受限的时段会看到价格尖峰。每个时段的功率平衡约束dual取出来就是该时段的能量价格序列把它们连成曲线就是典型的24小时出清价。注意多时段模型里LMP的取值要区分“该时段的lambda”和“跨时段爬坡约束的影子价格”后者会影响多个时段的价格形成分析时不要混在一起。6.2 MIP与LP的两阶段出清机组启停后怎么取LMP这是最容易被忽视的进阶问题。一旦引入机组启停变量变成0/1整数模型从LP变成MIP。从数学上严格说MIP的对偶乘子不再有经典LP影子价格的解释直接从混合整数解里取dual得到的“价格”是失真且没有经济含义的。行业里标准的做法是两阶段出清第一阶段用MIP求解机组启停和出力计划目标函数包含启停成本第二阶段把第一步解出来的0/1状态变量固定成常数重新求解一次纯LP从第二次LP里取对偶乘子作为LMP。YALMIP里做这个操作也很方便——binvar声明启停变量第一阶段optimize之后把启停变量的值用value提取出来替换成约束值再二次optimize。第二阶段模型不包含整数变量CPLEX自动走LP路径dual()就能正常取回乘子。我去年在某模拟项目X里第一次做带启停的出清时直接对MIP结果取dual算出来的价格低得离谱后来查文献才发现这个坑。从此固定了“先MIP定状态、再LP定价”的习惯。6.3 校验脚本与工程习惯最后说一个长期受益的习惯每次建新模型第一件事不是跑结果而是写一个自动校核脚本包含三类检查。一是潮流检查把解出来的P和theta代回Bf*theta确认潮流方程成立。二是平衡检查每时段的发电机总出力减去负荷误差不超过1e-6。三是价格检查用4.3节的扰动法随机抽三个节点验证LMP。这三关过了再谈报告和曲线。这套习惯帮我挡掉了无数个“看起来很美但实际算错”的版本。说到底出清优化最玄学的地方不在YALMIP语法也不在CPLEX配置而在“你有没有真正理解什么是影子价格”。我踩过的坑从符号搞反到预求解吞乘子没有一个能靠调参解决全是回到定义本身才找到答案。希望这篇笔记能帮你少走这些弯路把时间花在分析和业务上。本文还有配套的精品资源点击获取