简介这是一套面向电力市场与电力系统优化方向的程序包复现了《机组运行约束对机组节点边际电价的影响分析》中的核心算例适合研究生、高年级本科生及科研人员快速上手节点边际电价出清建模。程序基于Yalmip与Cplex求解单时段出清优化借助KKT对偶条件提取机组运行约束对应的拉格朗日乘子也就是影子价格并整理为矩阵形式同时清晰区分了节点边际电价与系统边际电价的差异。需要说明的是计算基于单时段模型未计及爬坡约束。包内共5个文件包括3个脚本程序、1篇CAJ格式原论文和1份Word版报告压缩包仅383KB轻量且结构简明。目前已有逾4600人学习下载。随附报告对建模思路、约束处理与结果讨论进行了完整说明程序注释详细便于对照复现作者还提供运行答疑适合作为课程设计或论文复现的基础工具。1. 节点边际电价出清yalmipcplex这套组合到底在算什么如果你正在做电力市场结算、现货出清仿真或者节点边际电价Locational Marginal PriceLMP测算大概率绕不开“基于yalmipcplex的电力市场-节点边际电价出清优化”这套东西。它解决的是一件事给定电网拓扑、机组报价、负荷和线路容量怎么在最小化总购电成本的同时算清楚每个节点应该收多少钱。这里有个反直觉的结论LMP并不是发电机组报价的简单平均值而是出清优化这个数学问题在功率平衡约束上拉格朗日乘子的“影子价格”。所以你把优化模型写对、把求解器调顺LMP自然就从对偶变量里读出来了。适合谁电力市场方向的研究生、交易中心或售电公司的结算测算人员以及想从linprog这种通用求解器往工业级CPLEX迁移的工程师。节点边际电价出清这件事数学建模本身不难难的是让yalmipcplex在你自己的数据上不翻车地把结果跑出来。2. 从出清模型到LMP先讲清节点边际电价在优化里怎么产生2.1 出清优化模型目标函数、约束和LMP的“影子价格”来源节点边际电价出清本质上是一个经济调度问题加网络约束。把时间维度先拿掉单时段的直流潮流DC-OPF模型长这样目标是最小化所有机组的发电成本约束有三类——功率平衡、发电机出力上下限、线路潮流限值。用数学语言写目标函数是 min sum_i (a_i * P_i^2 b_i * P_i)其中P_i是第i台机组出力a_i、b_i是成本系数线性项或二次项都可以。功率平衡约束是 sum_i P_i sum_d D_d也就是总出力必须等于总负荷。线路潮流约束用直流潮流近似通常写成 F_l sum_i PTDF_{l,i} * (P_i - D_i)要求 -F_max F_l F_max。这里PTDF是功率传输分布因子Power Transfer Distribution Factor它把每节点的净注入功率映射到每条线路的潮流上。LMP是怎么来的根据拉格朗日对偶理论功率平衡等式约束对应的拉格朗日乘子λ就是系统能量价格每条线路潮流约束对应的乘子μ_l就是阻塞价格分量。节点k的LMP λ sum_l (PTDF_{l,k} * μ_l)注意符号约定。当没有线路阻塞时所有节点LMP相等有阻塞时LMP会按PTDF加权拆分。这就是为什么你要用优化求解器而不是手工核算——因为LMP本质上是这个二次规划问题的对偶信息。2.2 为什么用YALMIP建模以及CPLEX在MATLAB里的角色我一般不会手写KKT条件或拉格朗日函数的显式表达式去求对偶变量太容易出错。YALMIP是MATLAB里的建模语言它最大的价值是让你用sdpvar声明变量直接用、运算符写约束然后它自动把模型翻译成求解器需要的稀疏矩阵格式。相比之下直接调CPLEX的MATLAB接口需要手动拼A矩阵、b向量、lb和ub一旦约束数量超过几十条拼错一个索引就要排半天错。CPLEX则是求解器引擎。它处理LP、QP、MIP的速度和数值稳定性在工业界是公认的。在同一个模型里YALMIP负责“把问题说清楚”CPLEX负责“真的把它算出来”。你可以在sdpsettings里把solver指定为cplexYALMIP就会自动生成CPLEX能吃的模型并且把CPLEX返回的对偶变量保存在约束句柄上。这也是这套方案比纯手写求逆或单纯用MATLAB自带linprog更有吸引力的地方你不需要自己推导对偶只需要知道YALMIP的dual(约束)这个函数怎么用。3. 把模型写成代码YALMIPCPLEX跑通LMP出清的最小实现3.1 数据准备母线、线路、机组参数怎么组织最小实现不需要多大的电网我建议先做一个3母线2机组的例子3条母线、2台发电机、2条线路、1个负荷节点。这样你能手工验算结果确认LMP逻辑是对的。母线数据写成一个结构体包含母线编号和负荷线路数据包含起始母线、终止母线、电抗、容量限值机组数据包含所在母线、成本系数、出力上下限。常见做法是把数据硬编码在脚本里但只要你换数据集就会意识到硬编码的痛苦。我一般会把数据准备和模型求解分开先写出一个函数load_case()返回上述三个结构体后续要换IEEE 30节点或118节点数据时只需要换这个函数。PTDF矩阵怎么来最可靠的路径是用节点导纳矩阵B去掉参考母线后的降阶版本取逆得到B_inv然后线路潮流 (1/x_l) * (theta_from - theta_to)而theta B_inv * (P_net)。写成矩阵就是PTDF。在MATLAB里这一步用几个矩阵运算就能完成注意参考母线slack bus那一行要删掉。3.2 最小代码块目标函数、功率平衡、线路潮流约束这里我给出一个完整的、可以直接抄走改数据的MATLAB脚本框架。假设你已经在MATLAB路径里配好了YALMIP和CPLEX。% 基于yalmipcplex的节点边际电价出清最小示例3母线2机组 % 步骤1定义数据 % 线路: 1-2, 2-3, 电抗均为0.1 pu, 容量限制为100 MW % 机组: g1在母线1, g2在母线2, 负荷在母线3 bus_load [0; 0; 150]; % 各母线负荷, MW gen_bus [1; 2]; % 发电机所在母线编号 gen_cost [20 0.05; 30 0.01]; % 每行: b线性成本$/MWh, a二次成本$/MW^2h gen_lim [0 100; 0 80]; % 出力上下限, MW line_from [1; 2]; line_to [2; 3]; line_x [0.1; 0.1]; line_cap [100; 100]; % 线路潮流上限, MW nbus 3; nline 2; ngen 2; % 步骤2计算PTDF矩阵直流潮流 % 构建节点导纳矩阵B不含参考母线这里设母线1为参考 B zeros(nbus); for l 1:nline f line_from(l); t line_to(l); B(f,f) B(f,f) 1/line_x(l); B(t,t) B(t,t) 1/line_x(l); B(f,t) B(f,t) - 1/line_x(l); B(t,f) B(t,f) - 1/line_x(l); end B_red B(2:end, 2:end); % 去参考母线 B_inv inv(B_red); % PTDF: 线路l对节点j的注入功率灵敏度 PTDF zeros(nline, nbus); for l 1:nline f line_from(l); t line_to(l); % 参考母线对应列置0非参考母线用B_inv的对应行 PTDF(l,1) 0; for k 2:nbus e_f zeros(nbus-1,1); e_t zeros(nbus-1,1); if f 1, e_f(f-1) 1; end if t 1, e_t(t-1) 1; end PTDF(l,k) (e_f - e_t) * B_inv / line_x(l); end end % 步骤3定义变量和约束 P sdpvar(ngen, 1); % 发电机出力变量 F sdpvar(nline, 1); % 线路潮流变量 net_inj zeros(nbus, 1); for i 1:ngen net_inj(gen_bus(i)) net_inj(gen_bus(i)) P(i); end net_inj net_inj - bus_load; % 净注入功率 Constraints []; for l 1:nline Constraints [Constraints, F(l) PTDF(l,:) * net_inj]; % 直流潮流等式 end for i 1:ngen Constraints [Constraints, gen_lim(i,1) P(i) gen_lim(i,2)]; end Constraints [Constraints, sum(P) sum(bus_load)]; % 功率平衡 Constraints [Constraints, -line_cap F line_cap]; % 线路容量 % 步骤4目标函数总发电成本最小 Objective 0; for i 1:ngen Objective Objective gen_cost(i,1) * P(i) gen_cost(i,2) * P(i)^2; end % 步骤5求解 options sdpsettings(solver,cplex, verbose, 1); optimize(Constraints, Objective, options); % 步骤6提取LMP功率平衡约束的对偶变量 lambda dual(Constraints(ismember(Constraints, [sum(P) sum(bus_load)]))); % 注意yalmip的dual取的是等式约束的拉格朗日乘子对不等式约束取的是影子价格 lmp zeros(nbus,1); for k 1:nbus lmp(k) lambda sum(PTDF(:,k) .* dual(Constraints(ismember(Constraints, ... [F(1) -line_cap(1); F(1) line_cap(1); F(2) -line_cap(2); F(2) line_cap(2)])))); end这段代码有几点要说明。第一PTDF的构造用了逐线路循环对每个节点算灵敏度这里为了可读性牺牲了效率真正跑IEEE 300节点时你会想把它向量化。第二Constraint句柄索引用了ismember方式这在YALMIP里可行但更稳妥的做法是把功率平衡约束拆出来单独命名Balance [sum(P) sum(bus_load)]然后加进Constraint列表时也保留变量名。第三目标函数里的二次项是严格凸的CPLEX会走QP求解路径LMP对应的是功率平衡约束的对偶变量二次成本不影响对偶提取方式。3.3 求解与LMP提取CPLEX返回的对偶变量怎么读上面的代码块里提取LMP的方式比较绕我实际项目里会写成更清晰的方式。先定义单独的约束句柄Balance_con (sum(P) sum(bus_load)); % 等式约束句柄 Line_con []; for l 1:nline Line_con [Line_con; -line_cap(l) F(l) line_cap(l)]; end Constraints [P_lim, Balance_con, Line_con]; optimize(Constraints, Objective, options); lambda dual(Balance_con); mu_low dual(Line_con(1:2:end)); % 线路下限约束的对偶 mu_up dual(Line_con(2:2:end)); % 线路上限约束的对偶这里的关键是YALMIP的dual函数返回的是对应约束的拉格朗日乘子。对等式约束它是能量价格的来源对不等式约束它只在约束活跃紧约束时才非零。LMP的节点价格要按PTDF把线路阻塞的对偶值折算回每个节点。手动算一遍你就会知道当线路22-3阻塞时节点3的LMP一定会高于节点1和2因为它是负荷端需要更贵的边际机组多出力。4. CPLEX配置与求解器参数让出清优化不翻车的几个必调项4.1 cplex配置到matlab路径、许可证和常见的安装坑“cplex配置到matlab”这个关键词是检索量很高的词也是很多新手第一次卡住的地方。安装完IBM ILOG CPLEX后你要把它提供的MATLAB接口路径加进MATLAB search path。通常路径长这样C:\Program Files\IBM\ILOG\CPLEX_Studio221\cplex\matlab\x64_win64。在MATLAB里执行addpath(genpath(...))然后保存path。常见翻车点有三个一是版本位数不匹配CPLEX Studio自带x64和win64两个目录MATLAB必须也是64位二是环境变量没有配可能在命令行工具里能跑但MATLAB调用不到三是license服务没启动运行cplexlicm或检查cplex.getVersion会直接报错。验证配置是否成功的标准动作我一般是跑yalmiptest或者直接在MATLAB里输入which cplexlp。如果返回了cplexlp的路径说明CPLEX工具箱已经认得如果返回错误先检查path、再检查许可证。还有一个隐藏坑YALMIP版本太老不一定认识新版本CPLEX。我习惯把YALMIP和CPLEX都更新到比较新的版本然后先跑一个optimize(sdpvar(1)0, sdpvar(1), sdpsettings(solver,cplex))能出结果再往下走。4.2 求解器参数MIP gap、输出精简、数值容差怎么调一旦model能求解下一步就是调参数。节点边际电价出清一般是LP或凸QPCPLEX默认参数通常已经不错但在大规模算例里你还是要主动设这几项。参数名作用我的经验值备注cplex.lpmethod选择线性规划求解算法0自动/ 2dual simplex大规模约束多dual simplex比primal更稳cplex.simplex.tolerances.optimality最优性容差1e-9精度差会造成LMP数值抖动cplex.display.func控制终端输出1看每次迭代的对偶进展cplex.mip.tolerances.mipgap整数最优gap1e-4如果模型有启停整数变量这个值决定停在哪verboseYALMIP输出级别12会输出太多0则出了问题无迹可寻设置方式是在sdpsettings里用点号比如options sdpsettings(solver,cplex, ... cplex.lpmethod, 2, ... cplex.simplex.tolerances.optimality, 1e-9, ... verbose, 1);注意参数名里的点号在YALMIP里是合法的你不需要转成下划线。如果你发现CPLEX求解时间异常先看是不是把每一条机组组合状态的0-1变量都塞进去了那是正常的NP-Hard如果纯LP检查约束里有没有大M法引入的冗余。前面这个3母线例子用默认参数是毫秒级完成的。5. 避坑电力市场出清优化中5个高频踩坑记录5.1 现象LMP为负或者节点价格差异离谱我最早在6母线系统上算LMP结果节点3的LMP是-150节点1却飙到300。查到最后发现是线路容量约束写反了方向。我在YALMIP里用了F -line_cap而不是F -line_cap导致约束根本没有起作用对偶变量取值完全错乱。解决方式很简单每次跑完要检查dual(Line_con)如果线路实际潮流远低于容量却有非零对偶那一定是约束方向写反了。另一个常见原因是PTDF矩阵的符号约定不一致直流潮流的正方向要和线路from/to定义对齐不然阻塞判断正好反了。5.2 现象同样的数据linprog能解但CPLEX卡死或报数值问题有一次我在一个200节点、300条线路的算例上CPLEX报了Numerical difficulties结果也不收敛。后来发现是线路电抗值直接用真实欧姆数数值跨度从0.001到0.1但PTDF直接拿去算了。解决方法是把整个系统标幺化电抗改成标幺值容量也按基准容量归一。电网数据最好在进YALMIP之前就做预处理以100 MVA为基准把欧姆、MW都转成标幺值。还有约束里如果存在成百上千条相似的上下限不等式YALMIP会生成稀疏矩阵注意不要用全零行去凑约束那会加重CPLEX预处理负担。5.3 现象dual函数读出来的对偶变量全是NaN或0这是YALMIP使用中非常经典的坑。dual(约束句柄)必须在optimize成功返回后立即调用而且必须引用同一个句柄对象。如果你把同一个约束写成Constraints [Constraints, P 0]后再重新赋值给新变量句柄索引可能就丢了。还有YALMIP里形如lb x ub的双边不等式对偶变量是分开存放在两个约束里的我踩坑就是直接对双边约束调用dual得到的是空的或零。正确做法是拆成x lb和x ub两个单独约束或者分别dual(Constraints(1))、dual(Constraints(2))。另外如果模型求解器返回的是infeasibledual自然没有意义先检查可行性。5.4 现象目标函数有二次项但LMP对偶值明显不对二次成本a*P^2存在时目标函数是凸的CPLEX用QP求解功率平衡约束的对偶变量依然是能量价格没问题。但是有的老版本YALMIP或CPLEX接口如果你把二次项写成sdpvar*Q*sdpvar而Q不是严格半正定CPLEX会直接拒绝或改走MIP近似。我遇到的是目标函数写成sum(a*P^2 b*P)没有任何问题。如果你为了让成本曲线分段线性而引入了大量整数变量那模型的凸性被破坏LMP的对偶意义就会模糊——这不是代码bug是模型设计问题。遇到这种情况老老实实用连续凸的二次成本。5.5 现象CPLEX报错“No solver available”或“cplex not found”这个常常发生在配置阶段表现为optimize报错说找不到cplex。原因通常是YALMIP缓存了旧solver列表或者MATLAB path里没有CPLEX接口目录。解决步骤clear classes然后重新addpath再yalmiptest确认。如果还不行看看环境变量PATH里有没有CPLEX的bin目录Windows上经常需要手动加。我自己曾经因为装了CPLEX 12.10后又装了别的软件改了环境变量结果CPLEX许可证找不到了报错信息却误导我以为是solver路径问题。最后用cplex.setup检查解决。6. 验证与进阶用参考LMP反推模型对错再往多时段扩展6.1 用2母线手工验算LMP逻辑在跑任何大算例之前先建一个2母线1线路的小系统机组A在母线1成本20 $/MWh机组B在母线2成本35 $/MWh母线2负荷150 MW线路容量100 MWh。无阻塞时机组A出力150两台LMP都等于20。阻塞后线路只能送100 MW机组B必须在母线2出力50系统边际成本升高到35此时节点1的LMP20 阻塞影子价格节点2的LMP35节点2的本地机组成为边际机组。你手工算出的LMP差异和代码里dual读出来的值应该对上。这个验证我每次换算法或换数据集都会跑一遍确保PTDF符号和dual提取逻辑没变。6.2 从单时段扩展到多时段出清单时段DC-OPF是电力市场的基本积木。下一步是多时段出清常见做法是给每个时段变量加下标t功率平衡和线路容量每个时段都成立如果还考虑机组爬坡约束需要在相邻时段间加不等式|P_i(t)-P_i(t-1)| ramp_i。YALMIP里可以定义sdpvar(ngen, T)的矩阵变量循环添加约束。CPLEX对这类带稀疏结构的大模型依然高效但要注意不要在每个时段都复制一整份PTDF公共数据提到循环外。我一般会把所有PTDF和线路容量做成稀疏矩阵再用repmat构造块结构这样内存占用小很多。6.3 最后的一个习惯先看对偶维度和约束编号我现在每次跑完出清都会先打印length(dual(Constraints))和约束列表确认对偶变量个数和约束个数一致。然后再画一张LMP随负荷变化的曲线看趋势是否合理。前两年有一次我在提取阻塞价格时用错了PTDF列顺序导致某个节点LMP和负荷变化方向相反还以为是市场力排了三天才发现是索引错位。这个习惯救了我不下五次。如果你也从这套yalmipcplex的方案起步记住节点边际电价出清的程序里最危险的不是公式而是索引和约束句柄。希望这些踩坑记录能帮你少走这段弯路。本文还有配套的精品资源点击获取