简介电力市场节点边际电价出清优化的完整复现方案面向电力市场研究人员、高年级本科生及研究生。资源基于史新红论文《机组运行约束对机组节点边际电价的影响分析》在单时段模型下采用YALMIPCPLEX求解器通过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)203025出力下限 (MW)000出力上限 (MW)100100100负荷 (MW)08060支路起点终点电抗 (p.u.)潮流上限 (MW)1120.20502130.30503230.2560定义机组出力向量P_g3×1节点相角向量θ3×1线路潮流向量F3×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)1A(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. YALMIPCPLEX 实现 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.mipgap1e-4混合整数优化的最优间隙越小越精确cplex.timelimit300求解时间上限防止死循环cplex.threads4并行线程数商用机建议设为物理核心数cplex.lpmethod00自动选择LP算法1主单纯形4 barrierverbose2YALMIP输出级别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.problem2线路限值过紧将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);lpmethod1表示主单纯形法对偶单纯形法是2。对于电力出清这种约束矩阵高度稀疏的问题单纯形法速度快并且给出的对偶乘子严格对应基顶点。如果模型是MILPCPLEX会在根节点和每个子节点调用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≤0dual也非负。实际项目中更多使用结算差额来反推影子价格是否正确因为结算差额只依赖LMP和发电成本不依赖约束编排顺序。5.2 扩展多时段和机组组合节点边际电价出清模型可以很容易扩展为多时段经济调度把时间维度加入所有变量将机组爬坡约束加入模型。此时决策变量变为pg(t,n)负荷pd(t,n)线路潮流f(t,k)。目标函数是各时段发电成本之和爬坡约束写成相邻时段的出力差限幅。如果还要处理机组启停就需要引入二进制变量模型变为MILPCPLEX依然可以高效求解。这时价格提取就要小心机组组合的节点边际电价通常从最后一次LP松弛的对偶变量中获得而该LP松弛是在整数变量固定之后求解的。在YALMIP中使用optimize得到MILP解后再对固定整数变量的原LP调用一次optimize提取此时的对偶变量。一个实用技巧是把出清模型写成MATLAB函数输入为负荷向量、线路参数和机组参数输出为LMP和阻塞价格并用单元测试固定几个已知场景。例如三节点系统在某一线路阻塞时节点间LMP差值应该与该线路的阻塞影子价格相关可以手工验证。这样后续换数据、换参数时不会因为某个约束顺序改变而导致对偶变量错位。最后一个小提醒当使用dual()提取对偶时务必在optimize成功之后立即执行不要在其他位置调用。YALMIP不会缓存对偶值重新计算或者修改模型后再调用dual()会得到错误结果或直接报错。我的习惯是将问题写成函数在函数内部完成optimize和dual提取再把结果返回这样能避免很多隐性bug。本文还有配套的精品资源点击获取 基于YALMIP与CPLEX的节点边际电价出清模型实现
简介电力市场节点边际电价出清优化的完整复现方案面向电力市场研究人员、高年级本科生及研究生。资源基于史新红论文《机组运行约束对机组节点边际电价的影响分析》在单时段模型下采用YALMIPCPLEX求解器通过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)203025出力下限 (MW)000出力上限 (MW)100100100负荷 (MW)08060支路起点终点电抗 (p.u.)潮流上限 (MW)1120.20502130.30503230.2560定义机组出力向量P_g3×1节点相角向量θ3×1线路潮流向量F3×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)1A(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. YALMIPCPLEX 实现 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.mipgap1e-4混合整数优化的最优间隙越小越精确cplex.timelimit300求解时间上限防止死循环cplex.threads4并行线程数商用机建议设为物理核心数cplex.lpmethod00自动选择LP算法1主单纯形4 barrierverbose2YALMIP输出级别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.problem2线路限值过紧将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);lpmethod1表示主单纯形法对偶单纯形法是2。对于电力出清这种约束矩阵高度稀疏的问题单纯形法速度快并且给出的对偶乘子严格对应基顶点。如果模型是MILPCPLEX会在根节点和每个子节点调用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≤0dual也非负。实际项目中更多使用结算差额来反推影子价格是否正确因为结算差额只依赖LMP和发电成本不依赖约束编排顺序。5.2 扩展多时段和机组组合节点边际电价出清模型可以很容易扩展为多时段经济调度把时间维度加入所有变量将机组爬坡约束加入模型。此时决策变量变为pg(t,n)负荷pd(t,n)线路潮流f(t,k)。目标函数是各时段发电成本之和爬坡约束写成相邻时段的出力差限幅。如果还要处理机组启停就需要引入二进制变量模型变为MILPCPLEX依然可以高效求解。这时价格提取就要小心机组组合的节点边际电价通常从最后一次LP松弛的对偶变量中获得而该LP松弛是在整数变量固定之后求解的。在YALMIP中使用optimize得到MILP解后再对固定整数变量的原LP调用一次optimize提取此时的对偶变量。一个实用技巧是把出清模型写成MATLAB函数输入为负荷向量、线路参数和机组参数输出为LMP和阻塞价格并用单元测试固定几个已知场景。例如三节点系统在某一线路阻塞时节点间LMP差值应该与该线路的阻塞影子价格相关可以手工验证。这样后续换数据、换参数时不会因为某个约束顺序改变而导致对偶变量错位。最后一个小提醒当使用dual()提取对偶时务必在optimize成功之后立即执行不要在其他位置调用。YALMIP不会缓存对偶值重新计算或者修改模型后再调用dual()会得到错误结果或直接报错。我的习惯是将问题写成函数在函数内部完成optimize和dual提取再把结果返回这样能避免很多隐性bug。本文还有配套的精品资源点击获取