简介这套基于YALMIPCPLEX的电力市场出清优化包面向电力市场方向的研究生与科研人员完整复现《机组运行约束对机组节点边际电价的影响分析_史新红》的核心方法。程序采用KKT对偶条件求解拉格朗日乘子影子价格以矩阵形式实现节点边际电价出清模型单时段设定、未考虑爬坡约束边界条件清晰适合作为教学与科研基线。压缩包共5个文件含3个MATLAB程序主程序与辅助函数、1篇CAJ原始文献和1份DOC完整报告合计约383KB目录结构紧凑。目前已有4668人学习/浏览。通过对照报告与代码注释可系统掌握节点边际电价的正规建模流程、机组约束作用机制与影子价格求解细节程序清晰美观、注释充分运行中遇到问题还可答疑对开展电力市场出清分析具有直接参考价值。1. 一个 zip 装下的电力市场出清节点边际电价与 yalmipcplex 到底在解决什么拿到「基于yalmipcplex的电力市场-节点边际电价出清优化.zip」这类代码包很多人第一件事是装好 YALMIP 和 CPLEX然后双击主脚本直接 run。实际结果往往三选一求解器报许可证错误、模型直接不可行、或者出清价格算出来是负的。问题大多不在求解器而在对出清模型本身的理解——节点边际电价LMP不是一个全局统一定价它是网络约束下的影子价格线路一堵同一个市场里不同节点就会出不同的价。YALMIP 负责把优化问题用接近数学表达式的写法描述出来CPLEX 负责把大规模线性或混合整数优化快速解掉这两者组合几乎是电力市场出清方向最主流的技术栈。本文用一套三节点系统从模型拆解、完整可运行代码、CPLEX 参数设置到必须规避的五个坑把这条链路一次讲透。适合手里有 MATLAB、想真正跑通并看懂 LMP 出清结果的研究生和一线从业者。2. 出清优化模型拆解LMP 的三分量与 DC OPF 约束族2.1 为什么节点电价不一样能量、阻塞与损耗节点边际电价在现货市场里的物理含义是在满足所有运行约束的前提下某节点新增一兆瓦负荷时全系统发电成本的最小增量。这个增量不是一个全局常数因为输电网络会让不同节点的“边际成本”产生空间差异。LMP 在工程上通常被拆成三个分量能量分量全网统一的系统能量边际成本由最贵的在线边际机组决定相当于无阻塞时的统一出清价。阻塞分量由于线路潮流越限而抬升或压低的电价偏移阻塞越严重不同节点的价格差越大。损耗分量与网络损耗有关直流潮流模型DC OPF里通常忽略交流模型或带网损修正的模型才会单独算。最常见的出清实现是安全约束经济调度SCED目标是最小化总购电成本约束覆盖功率平衡、机组上下限、爬坡和线路潮流。节点边际电价就是这套优化问题的拉格朗日乘子组合而不是某个机组报价的简单搬运。这也是为什么很多新手直接把最低报价当成出清价、发现对不上的原因。2.2 DC OPF 出清模型目标函数与四条约束DC 最优潮流DC OPF是出清优化的标准内核。它把交流潮流简化成线性模型忽略无功和网损只保留有功功率和相角关系从而让整个问题变成线性规划既能保证求解速度和数值稳定性又方便从对偶乘子中提取价格信号。真实电力市场做日前出清普遍采用 DC OPF 加安全校核的方式AC 潮流在事后校核和辅助服务市场里用。标准模型写成如下形式目标函数minimize ∑ C_i · P_i其中 C_i 是机组 i 的报价$/MWhP_i 是出清出力MW。约束族分四类功率平衡全网发电等于负荷加网损DC 模型里通常忽略网损写成 ∑P ∑D。机组上下限P_i_min ≤ P_i ≤ P_i_max对应机组物理出力范围。爬坡约束多时段模型中相邻时段出力差不能超过爬坡速率|P_i,t - P_i,t-1| ≤ R_i。线路潮流约束直流模型里线路潮流 F_l (θ_f - θ_t)/x_l且 |F_l| ≤ F_l_maxθ 是节点相角x_l 是线路电抗。出清这个词的实质就是把这四条约束和一个目标函数一起丢给求解器求出的最优解就是机组出清出力配套的对偶乘子就是节点电价。模型里把线路容量改大阻塞消失LMP 就会落回同一个值把某条线路容量调小价格立刻分化。这个“调参数看价格变化”的过程是理解 LMP 最有效的实验。2.3 为什么是 YALMIP CPLEX而不是直接写 linprogMATLAB 自带的 linprog 能解线性规划但需要把问题手抄成 A、b、Aeq、beq 矩阵出清模型的约束一旦上了爬坡和线路容量矩阵维度很容易错位改一个参数要动一大片代码。YALMIP 的价值在于用符号变量描述优化问题代码和教科书公式几乎一一对应新增一类约束只是加一行约束表达式改模型成本非常低。CPLEX 的价值在于求解能力和对偶信息。出清优化看起来是普通 LP实际项目里会掺入机组启停 0/1 变量变成 MILPCPLEX 在 MIP 求解上的成熟度和数值稳定性明显优于 MATLAB 自带求解器。更实际的一点是CPLEX 返回的对偶乘子质量高YALMIP 的 dual 函数能直接提取不需要自己从原始解反推。如果哪天换了 Gurobi 或 MOSEK同样的 YALMIP 模型只改一个 sdpsettings 参数就能切换这就是建模工具箱带来的最大红利。3. 三节点出清最小实现数据、代码与结果逐个对上3.1 测试系统设计让阻塞主动出现在 1-2 线路教学用的三节点系统不需要复杂但要设计成“无阻塞一个价、有阻塞价差明显”这样跑出来的结果才有对比价值。我常用的参数如下项目参数数值线路 1-2电抗 / 容量0.1 p.u. / 80 MW线路 1-3电抗 / 容量0.1 p.u. / 150 MW线路 2-3电抗 / 容量0.1 p.u. / 150 MW机组 G1所在节点 / 报价 / 上限节点 1 / 14 $/MWh / 150 MW机组 G2所在节点 / 报价 / 上限节点 3 / 20 $/MWh / 150 MW负荷所在节点 / 功率节点 2 / 150 MW三条线路电抗一致功率会优先走 1-2 这条直接路径。当 G1 满发满足全部负荷时1-2 线路上的潮流会接近甚至超过 80 MW阻塞自然出现此时必须让节点 3 的 G2 多发把便宜电挤出一部分节点 2 的 LMP 就会高于节点 1。这个系统虽然只有两个机组但已经具备 LMP 分化的完整要素。注意线路容量用的是 MW电抗用标幺值 p.u.基准容量取 100 MVA。代码里建模时要把这个基准关系理清楚具体写在 3.2 节的变量说明里。3.2 完整建模代码YALMIP 变量、约束与 CPLEX 求解以下是完整脚本直接从 MATLAB 新建脚本复制即可运行。代码分成数据、建模、求解、输出四段每段都留了注释。%% 三节点 DC OPF 出清YALMIP CPLEX clear; clc; %% 1. 系统数据 S_base 100; % 基准容量 MVA潮流约束里要用 % 线路: [首端 末端 电抗(p.u.) 容量(MW)] lines [ 1 2 0.1 80 1 3 0.1 150 2 3 0.1 150 ]; % 机组: [节点 报价($/MWh) Pmin(MW) Pmax(MW)] gen [ 1 14 0 150 3 20 0 150 ]; % 负荷: [节点 功率(MW)] loads [2 150]; n_bus 3; %% 2. 节点电纳矩阵 B标幺值 B zeros(n_bus, n_bus); for k 1:size(lines, 1) f lines(k,1); t lines(k,2); x lines(k,3); B(f,f) B(f,f) 1/x; B(t,t) B(t,t) 1/x; B(f,t) B(f,t) - 1/x; B(t,f) B(t,f) - 1/x; end %% 3. 决策变量 theta sdpvar(n_bus, 1); % 节点相角弧度 Pg sdpvar(size(gen,1), 1); % 机组出力 MW %% 4. 节点功率注入 Dvec zeros(n_bus, 1); for k 1:size(loads,1) Dvec(loads(k,1)) loads(k,2); end Inj zeros(n_bus, 1); for k 1:size(gen,1) Inj(gen(k,1)) Inj(gen(k,1)) Pg(k); end %% 5. 约束 Cons []; % 功率平衡S_base*B*theta 是 MW两边量纲统一 Bal (S_base * B * theta Inj - Dvec); Cons [Cons, Bal]; % 参考节点相角 Cons [Cons, theta(1) 0]; % 机组上下限 for k 1:size(gen,1) Cons [Cons, gen(k,3) Pg(k) gen(k,4)]; end % 线路潮流约束单向不等式分开写方便后续提取对偶乘子 Flow cell(size(lines,1), 1); for k 1:size(lines,1) f lines(k,1); t lines(k,2); x lines(k,3); Flow{k} (theta(f) - theta(t)) / x * S_base; % MW Cons [Cons, Flow{k} lines(k,4)]; Cons [Cons, -lines(k,4) Flow{k}]; end %% 6. 目标函数总购电成本最小 Objective gen(:,2) * Pg; %% 7. 求解 ops sdpsettings(solver, cplex, verbose, 1); sol optimize(Cons, Objective, ops); %% 8. 输出 if sol.problem ~ 0 disp([求解失败: sol.info]); else fprintf(最优总成本: %.2f $/h\n, value(Objective)); for k 1:size(gen,1) fprintf(G%d bus %d: %.2f MW, 报价 %.0f $/MWh\n, ... k, gen(k,1), value(Pg(k)), gen(k,2)); end for k 1:size(lines,1) fprintf(线路 %d-%d 潮流: %.2f MW (上限 %.0f MW)\n, ... lines(k,1), lines(k,2), value(Flow{k}), lines(k,4)); end end代码里的几个关键点值得展开说明。功率平衡约束写成S_base * B * theta Inj - Dvec而不是直接B * theta ...。原因是电抗、导纳都是标幺值B 乘以 theta 得到的是标幺功率不乘回基准容量就对不上 MW 单位的注入向量。这一行是新手最容易抄错然后得到不可行解或量纲错乱结果的地方。线路潮流表达式(theta(f) - theta(t)) / x * S_base同样处理了标幺到有名值的转换。注意线路潮流方向定义为首端流向末端容量约束写成双向不等式其中每条线路的上限约束单独命名这样后面用 dual 提取阻塞对应的对偶乘子时乘子序号和方向不会混乱。ops sdpsettings(solver, cplex, verbose, 1)里 verbose 控制求解日志输出。0 是静默1 只输出关键步骤2 会打印 CPLEX 每一轮迭代定位问题时用 2正常跑用 0 或 1 即可。3.3 结果对照机组出力、线路潮流与能量价格以这套数据跑出来的结果正确方向应该是G1 出力被 1-2 线路压到 90 MW 附近G2 出力约 60 MW 补足负荷缺口1-2 线路潮流卡在 80 MW 上限总购电成本明显高于无阻塞时的 2100 $/h150 MW × 14 $/MWh。我做这个示例时最常用的验证方式是先把lines里 1-2 的容量改成 150 MW 跑一遍确认无阻塞场景下 G1 出 150 MW、全网 LMP 回到 14 $/MWh再改回 80 MW 跑第二遍看价格是否分化。两个场景的对比足以验证模型逻辑正确。如果需要从模型里直接读能量价格可以对功率平衡约束取对偶lambda_energy dual(Bal);Bal是向量约束dual 返回一个列向量每个元素对应一个节点的平衡约束乘子。需要注意的是三个节点的 KCL 方程线性相关乘子本身不唯一CPLEX 会返回一组可行对偶解但取哪个作为系统能量价格有约定。我在第 6 章会给出更可靠的数值校验方法这条先留意。4. CPLEX 求解器接入与调参从能跑到跑稳4.1 让 YALMIP 找到 CPLEX路径、测试与许可证CPLEX 本身是独立求解器要让它被 MATLAB 里的 YALMIP 调用需要把 CPLEX 安装目录下的 MATLAB 接口文件夹加入 MATLAB 搜索路径。常见做法是在 MATLAB 里执行addpath(genpath(CPLEX安装目录/cplex/matlab))然后保存路径让下次启动仍然生效。装好之后先跑一个命令确认链路通畅yalmiptest输出列表里如果 cplex 一栏显示 found说明 YALMIP 已经能识别。如果显示 not found优先检查两件事路径是否真的加进去了CPLEX 版本和 MATLAB 版本是否匹配。新版 CPLEX 对 MATLAB 版本非常敏感版本差距大时接口目录可能无法加载这时换个与 MATLAB 版本对应的 CPLEX 版本是最快的解决方式。许可证是另一个独立问题。CPLEX 通常需要许可证文件学术版本或商业版本安装后要确保环境变量ILOG_LICENSE_FILE指向正确的许可证路径。如果跑yalmiptest时 cplex 存在但实际optimize一调用就报许可错误大多是许可证没被 MATLAB 进程读到重新设置环境变量并重启 MATLAB 即可。4.2 sdpsettings 参数表MIP gap、LP 算法、线程与时间出清模型从 LP 到 MILPCPLEX 的可调参数里有几个几乎每次都要碰。以下是常用参数速查表参数名默认值作用建议verbose0控制求解日志详细程度调问题时设 2正常跑设 0cplex.lpmethod0LP 算法0 自动1 单纯形2 对偶单纯形4 Barrier大规模出清常设 4对偶信息也稳定cplex.mip.tolerances.mipgap0.0001MIP 相对间隙决定求解终止精度工程场景可放宽到 0.001 提速cplex.threads0并行线程数0 为自动内存够就设物理核数 8~16cplex.timelimit1e75求解时间上限秒线上出清必须设置如 600savesolveroutput0保存 CPLEX 原始输出到 sol.solveroutput定位异常时设 1实际使用时我会在脚本里写一个统一的 ops 配置ops sdpsettings(solver, cplex, ... verbose, 1, ... cplex.lpmethod, 4, ... cplex.threads, 8, ... cplex.timelimit, 600);Barrier 算法lpmethod4对大型 SCED 的稀疏结构特别友好内存占用和迭代次数通常优于单纯形。但要注意Barrier 求解完成后 CPLEX 默认会交叉到单纯形以获得顶点解这对提取对偶乘子更友好不用刻意关闭交叉选项。MIP gap 在纯 LP 模型里不生效一旦加了机组启停变量它直接决定求解时间。0.001 的间隙对绝大多数出清分析已经够用硬要往 1e-6 抠会让 MILP 多跑几个小时而价格几乎没变化。4.3 规模化建模从三节点到上千节点的习惯把三节点换成 IEEE 标准节点或实际电网数据时最容易翻车的不是求解器而是建模时的编码习惯。我见过大量代码在 3 节点上跑得飞快换到几百节点后建模阶段就卡住几十分钟问题出在 for 循环里反复拼接约束矩阵。YALMIP 里约束拼接用[Cons, Cons, newConstraint]在节点少时没感觉节点多了会产生大量临时对象。更好的做法是先把约束放入 cell 数组最后一次性合并ConsList cell(n_lines * 2 n_gen n_bus 1, 1); idx 1; % ... 填充 ConsList{idx}idx 自增 Cons [ConsList{:}];多时段 SCED 的变量用二维 sdpvarPg sdpvar(n_gen, n_t)爬坡约束写成矩阵形式而不是逐时段循环RampCons [-R_max diff(Pg, 1, 2) R_max];diff一次生成所有相邻时段差值代码简洁且约束数量可控。数据规模上CPLEX 处理几千个节点、上万条约束的 LP 完全是常规负载真正的瓶颈往往在建模层而不是求解层。5. 出清模型必踩的五个坑现象、原因与解决办法5.1 CPLEX 报许可证错误模型一个都跑不了现象optimize返回sol.problem -1命令行提示 CPLEX Error比如CPLEX Error 1016: License not available。原因CPLEX 安装了但 MATLAB 进程拿不到许可证。常见于许可证服务未启动、环境变量ILOG_LICENSE_FILE没设置、或学术许可证绑定的机器信息与当前环境不一致。解决先在系统命令行执行 CPLEX 自带的交互程序确认许可证本身可用再检查环境变量是否设置到了正确路径。最后在 MATLAB 里重启会话让环境变量生效。如果许可证没问题检查 CPLEX 插件目录是否与 MATLAB 版本匹配版本错位也会报出一些看起来像许可证错误的异常。5.2 一加线路容量约束就不可行现象不带线路约束时模型正常求解加上Flow limit后 CPLEX 直接报 infeasible。原因最常见是功率平衡与容量约束在数学上存在冲突。比如负荷 150 MW 全部压在单条线路上而这条线路容量只有 80 MW无论机组怎么出都没有可行解另一个常见原因是线路潮流表达式的单位没统一潮流算出来比实际大 100 倍导致容量约束把所有通道全掐死。解决把线路容量逐一放大观察哪一条从 infeasible 变成 feasible定位到具体瓶颈同时用value(Flow{k})检查潮流量级是否正常。如果单位对不上回到 3.2 节检查S_base有没有乘到潮流表达式里。5.3 对偶乘子符号反了LMP 算成负数现象节点边际电价低于最低机组报价或者负荷节点的价格比发电节点还低出现明显的负阻塞分量。原因YALMIP 的dual函数返回的对偶乘子符号与公式里的符号约定相关而不同求解器的返回符号并不完全一致。很多人套用论文里的 LMP 公式把乘子符号直接代进去结果恰好颠倒。解决先用一个无阻塞的两节点系统手算验证。两节点单线路、单负荷、单机组节点 1 发电节点 2 负荷LMP 应当等于机组报价且两个节点相同。跑通这个对照后再检查带阻塞场景的乘子符号。我自己的习惯是永远用数值差分法做最终校验方法见第 6 章不依赖求解器返回来猜符号。5.4 相角没设参考节点雅可比直接奇异现象求解结果里 theta 数值异常大或者 CPLEX 报 numerical issues有时甚至直接不可行。原因直流潮流的 B 矩阵秩为 n-1必须有参考节点相角为零约束才能满秩。没加theta(1) 0时任意相角整体偏移都不改变潮流求解器在数值上无法收敛。解决在约束里显式加入Cons [Cons, theta(1) 0]。注意参考节点的选择会影响部分对偶乘子的数值但对 LMP 的分量分解没有物理影响。5.5 多时段模型爬坡约束写错下标结果看起来合理实际漏约束现象爬坡速率设得很小但相邻时段出力跳变明显超过速率结果却不报错。原因循环里写成了Pg(:,k) - Pg(:,k1)但下标范围多走了一位或者用了diff之后忘记把首时段排除。更隐蔽的原因是爬坡约束里用了绝对值函数abs(...)YALMIP 会产生额外变量改变问题结构求解器给出的是等价 MIP 的解但数值容差下偶尔会留出小毛刺。解决输出结果后用一段简单脚本自动检查相邻时段差值超过速率就打印出来比肉眼看图可靠得多。爬坡约束尽量写成-R P_next - P_prev R避免使用abs让模型保持线性结构。5.6 用 sol 结构与调试命令把不可行问题切成两半当模型不可行时不要盯着代码发愣按这套流程走能省半天时间。先看sol.problem与sol.infoCPLEX 返回的消息通常直接指出冲突的约束类别然后去掉线路约束看是否恢复可行如果恢复问题在网络约束与发电能力的配合上如果还不可行再去掉负荷的绝对值改为逐步加载找到第一个不可行的窗口。YALMIP 提供assign与check可以手动指定变量值后逐条检验约束残差对定位下界冲突特别有用assign(Pg, [100; 50]); assign(theta, zeros(3,1)); check(Cons);check输出每条约束的残差和是否满足残差最大的几条就是冲突源。配合 5.2 的单位检查绝大多数不可行问题都可以在十分钟内定位。6. 进阶多时段 SCED、LMP 数值校验与结算口径6.1 多时段 SCED 的变量与爬坡约束写法从单时段扩展到 24 时段核心是把机组出力变成矩阵变量并加入爬坡和启停。变量定义改成Pg sdpvar(n_gen, 24)目标函数变成sum(sum(gen_cost .* Pg))爬坡约束用 diff 一次性生成代码量增加不到十行。如果是机组组合还要引入on binvar(n_gen, 24)启停变量并把出力上下限改写成on .* Pmin Pg on .* Pmax。这一步会让问题从 LP 变成 MILP求解时间从秒级跳到分钟级属于正常现象。6.2 用数值差分给 LMP 做体检符号约定太容易翻车所以我现在每个出清模型跑通后都会加一段独立的数值校验脚本。原理很简单把某节点负荷增加 1 MW重新求解出清优化目标函数的增量就是该节点的 LMP 数值近似。LMP_num zeros(n_bus, 1); for i 1:n_bus Dvec_test Dvec; Dvec_test(i) Dvec_test(i) 1; % 该节点负荷 1 MW % 重新求解得到新的目标值 % LMP_num(i) new_objective - base_objective end对三节点系统这个循环跑三次即可完成校验。把数值差分结果和模型里提取的 LMP 做对比两者偏差应在一个很小的范围内。如果差分法是正的而模型提取是负的不用怀疑提取公式的符号一定有问题以差分法为准去修正。这个习惯是我交了学费才养成的。最早做某跨平台系统的出清模块时对偶乘子提取公式少了一个转置算出的阻塞分量方向和实际完全相反整整排查了一个下午最后就是靠这十几行差分脚本定位的。从那以后差分校验成了每个出清模型的出厂测试项宁可多跑几次 LP也不靠眼睛判断价格方向。希望帮到你。本文还有配套的精品资源点击获取
