简介本资源面向电力系统专业本科生、研究生及从事最优潮流研究的工程师提供基于Matlab内点法求解IEEE 14节点系统最优潮流OPF的完整实现方案聚焦燃料费用最小化这一典型经济调度目标。压缩包共5个文件45KB含3个关键参数文本文件节点、支路与发电机数据、1个核心Matlab求解脚本封装fmincon调用与内点算法配置及1张14节点系统拓扑图便于快速建模与结果可视化验证。已有391人学习下载资源结构简洁实用脚本直接调用内点法处理非线性约束优化问题参数文件采用标准格式便于替换扩展拓扑图辅助理解网络结构与潮流分布。读者可即刻运行复现最优发电出力分配、总燃料成本及各节点电压/功率结果掌握OPF建模要点、Matlab优化工具箱实战技巧及电力系统经济调度的核心分析逻辑。 做电力系统优化的同学应该都绕不开最优潮流这块硬骨头。我刚入门的时候天天用Matpower点一下runopf觉得OPF也不过如此后来被要求“不许调库自己写一个内点法求IEEE 14节点系统最优潮流目标函数是燃料费用最小”这才意识到那个“一键出结果”背后全是数学和工程细节。这篇博客就用一篇实战笔记的形式讲清楚怎么从零写一个Matlab内点法程序求解14节点系统的最优潮流。整个过程按“问题建模、算法原理、代码实现、结果验证、踩坑记录”五步展开。适合正在做毕业设计、复现论文方法或者想把OPF底层逻辑吃透的同行参考。1. 先建模14节点系统最优潮流的问题定义与数学形式1.1 为什么要选IEEE 14节点系统在电力系统潮流与优化研究里IEEE 14节点系统是出现频率最高的“入门级标准算例”。它规模不大数据公开又有足够的复杂度来体现算法特性系统里一共14条母线、20条支路、5台发电机总负荷约259MW基准容量100MVA。它包含了普通PQ节点、PV节点、平衡节点还有变压器支路和并联电容器基本覆盖了OPF建模会遇到的大部分设备类型。选它做内点法验证有一个很现实的好处系统规模小状态变量只有几十个雅可比矩阵和海森矩阵都是几十维即使用最笨的稠密矩阵也能在零点几秒内跑完。这让你可以把精力集中在算法逻辑本身而不是矩阵稀疏化和大规模求解上。等你在14节点上把内点法调通了再换IEEE 30节点、118节点核心代码几乎不用改只需要替换数据文件这种“可迁移性”是很值的。另外14节点系统的基准数据在各种工具里都能直接找到尤其是Matpower的case14已经内置了母线、支路、发电机参数。你可以用自己的内点法结果和Matpower的opf结果做交叉验证这是判断程序正确性最快的方法。很多论文里做OPF算法对比也是拿14节点做第一个数值实验因为它能快速暴露算法实现的问题又不会因为系统规模掩盖数值细节。1.2 决策变量、目标函数与机组参数最优潮流本质上是一个带约束的非线性规划问题。目标是最小化系统总燃料费用也就是所有发电机的发电成本之和每台发电机的燃料费用通常用有功出力的二次函数表示min F Σ_i ( a_i * P_gi^2 b_i * P_gi c_i )其中P_gi是第i台发电机的有功出力单位MWa_i、b_i、c_i是费用系数。a_i反映二次成本b_i是一次成本系数c_i是空载费用常数。c_i不影响最优解的取值位置因为常数项的梯度为零但会影响总费用数值代码里一般保留。决策变量分两类一类是节点状态变量包括每个节点的电压幅值V_i和相角θ_i另一类是控制变量包括每台发电机的有功出力P_gi和无功出力Q_gi。在14节点系统里如果不把平衡节点单独处理状态变量总数是14个电压幅值、14个电压相角、5个有功出力、5个无功出力一共38个变量。实际编程时往往固定平衡节点相角为0把它的相角排除在变量之外具体看实现习惯。我这里采用的是一套教学用的典型二次费用系数如下表机组所在母线abcPmax(pu)Pmin(pu)G110.001015.0103.320.10G220.001214.5101.400.10G330.002018.0101.000.05G460.002018.0101.000.05G580.001513.0101.000.05注意这里的a、b对应的P单位是MW也就是把标幺值出力P_pu乘上基准容量100MVA之后代入费用函数。有的文献喜欢在标幺值下写费用函数那a、b的数值会有另一个量级这是正常的。编程前先把单位统一否则后面怎么算都不对。这套系数不是Matpower内置的默认系数Matpower的case14默认用的是线性成本模型需要自己替换成二次多项式。1.3 等式约束与不等式约束怎么列约束条件是OPF和普通经济调度的最大区别。OPF必须满足交流潮流方程也就是每个节点的有功和无功都要平衡。对每个节点i有功等式约束可以写成P_i(V, θ) - P_gi P_di 0其中P_i是注入功率P_di是该节点负荷。无功等式约束同理Q_i(V, θ) - Q_gi Q_di 0这里的V和θ是所有节点的电压向量功率注入方程和牛顿-拉夫逊潮流里的公式完全一样这部分代码可以直接复用潮流计算的mismatch函数。不等式约束则包括四类发电机有功出力上下限、发电机无功出力上下限、节点电压幅值上下限、支路传输功率上限。电压幅值一般限制在0.94到1.06pu之间支路功率约束可以选择视在功率或电流限值看你的数据源怎么定义。每条不等式约束后续都需要引入松弛变量把“小于等于”转成“等式加非负约束”这是内点法能处理它的关键。在14节点系统里不等式约束数量不算多但类型已经足够丰富。建模时我建议把所有不等式约束统一排序统一编号记录每个约束对应的上下限、对应变量索引后面组装KKT系统时就会省很多事。千万不要把约束分散写在不同的代码块里那样调试的时候你会疯掉的。2. 内点法核心逻辑为什么它能高效求解OPF2.1 从“障碍”两个字理解内点法内点法这个名字里最关键的词是“内点”。它不像单纯形法那样沿着可行域的边界找顶点也不像罚函数法那样靠不断试探惩罚系数逼近最优而是从一个严格在可行域内部的点出发靠“障碍项”把迭代点牢牢挡在可行域内部然后逐步逼近边界上的最优点。这个“障碍项”怎么理解呢你可以想象一道物理围墙墙的位置就是不等式约束的边界。优化变量越靠近墙障碍函数的值就越大趋于无穷大这样迭代点就不敢越界。但如果我们始终不许碰墙那最后怎么到达最优解呢方法就是引入一个不断缩小的障碍参数μ。μ大的时候墙离真实边界很远解被压在内部很安全但离最优远μ逐渐变小墙慢慢向真实边界移动解也一步步逼近真正的最优解。当μ趋于0时障碍项消失解就落在边界上了。从数学上看原问题经过松弛和对数变换后变成min f(x) - μ * Σ ln(s_i) s.t. c(x) 0 g(x) s h s 0这里的s是松弛变量。因为ln(s)在s接近0时趋向负无穷等于给目标函数加了一道“软墙”。原对偶内点法的核心就是不断求解这个带μ的系统的KKT条件同时更新μ让解逐渐逼近原问题的最优解。2.2 原对偶内点法的迭代主干要写代码光理解障碍函数还不够还得知道每次迭代具体干什么。原对偶内点法的基本步骤可以用一句话概括把KKT条件写成一组非线性方程组然后用牛顿法迭代求解。具体来说每次迭代做这几件事根据当前变量x、对偶变量λ和μ这里μ指不等式约束对应的对偶变量注意和障碍参数的符号区分计算KKT残差向量。组装牛顿修正方程也就是KKT矩阵这个矩阵是稀疏对称的。求解线性方程组得到变量和对偶变量的修正方向。计算步长保证松弛变量和对偶变量在更新后仍然严格大于0。沿修正方向更新所有变量计算互补间隙gap。如果gap和残差都小于收敛阈值停止迭代否则缩小障碍参数μ继续。这里的“互补间隙”是内点法里最核心的收敛指标。它反映的是原始可行性、对偶可行性以及互补松弛条件的满足程度。每次迭代结束后我们计算gap μ的对偶变量和松弛变量内积之和μ_new σ * gap / n_ineq其中n_ineq是不等式约束个数。当gap小于1e-6时说明已经逼近最优解了。牛顿修正方程的组装是编程里最重的部分。它的矩阵结构大概是[ H Jc A ] [ Δx ] [ r_x ] [ Jc 0 0 ] [ Δλ ] [ r_c ] [ A 0 -Σ ] [ Δμ ] [ r_μ ]H是拉格朗日函数对x的海森矩阵Jc是等式约束雅可比A是不等式约束雅可比Σ是对角阵。14节点系统规模小直接用Matlab的backslash解这个方程就行系统大了再去考虑稀疏分解和预处理那是另一个话题。2.3 三个关键参数μ、σ、步长系数内点法工程实现里最影响成败的是三个参数障碍参数初始值μ0、中心参数σ、步长缩放系数α_max。μ0一般取0.1左右。如果取得太大前几次迭代会花很多时间在“消除障碍”上取得太小初始点可能离可行域太远线性化偏差大容易发散。中心参数σ是每次迭代更新μ时用的比例因子通常在0.1到0.2之间。σ越小μ衰减越快算法收敛快但数值上更容易震荡σ越大收敛越稳但迭代次数变多。我自己的经验是14节点系统取σ0.2稳定和速度都比较均衡。步长计算是内点法一个特别容易出问题的地方。虽然牛顿方向是朝着最优解走的但如果一步跨太大松弛变量可能变成负数违反s0的约束。所以必须算一个“最大安全步长”对修正方向中所有负分量求出允许的最大步长然后乘以一个小于1的安全系数。这个系数我一般取0.9995也就是只走最大安全步长的99.95%留一点余量。如果你设置成1.0迭代点会正好落在边界上导致下一步障碍函数无穷大直接崩掉。3. Matlab代码实现从初始化到收敛的完整流程3.1 数据准备直接基于case14改写代码第一步不是写主函数而是把系统数据准备好。我强烈建议不要自己手敲14条母线的数据直接读取Matpower的case14然后覆盖掉它的发电机费用部分。这样做的好处很明显case14是经过无数人验证的标准数据母线、支路参数不会有错你只需要关注OPF算法本身。具体操作分两步。第一步用Matpower提供的loadcase函数把case14数据读进来然后提取bus、branch、gen三个结构体。第二步把gen里费用相关的列换成我们自己要用的二次系数。Matpower的gen矩阵前6列是出力限值后续列是参与因子和启停状态费用系数通常要另外定义或者直接重写一个gen结构体。我这里给出一个比较清晰的数据定义方式% 系统基准容量 baseMVA 100; % 发电机数据[母线号, Pmax(pu), Pmin(pu), Qmax(pu), Qmin(pu), a, b, c] % 费用函数cost a * P_mw^2 b * P_mw c gen [ 1 3.32 0.10 0.98 -0.20 0.0010 15.0 10; 2 1.40 0.10 0.48 -0.40 0.0012 14.5 10; 3 1.00 0.05 0.40 -0.20 0.0020 18.0 10; 6 1.00 0.05 0.24 -0.06 0.0020 18.0 10; 8 1.00 0.05 0.24 -0.06 0.0015 13.0 10; ]; % 电压限值 Vmin 0.94 * ones(14, 1); Vmax 1.06 * ones(14, 1);母线数据和支路数据直接从case14里取不要自己手抄。mpc loadcase(case14); mpc.gen gen;然后用mpc.bus和mpc.branch作为后续计算的输入。3.2 主循环与修正方程组装主循环是程序的骨架核心逻辑其实不到一百行。下面这段是内点法的主循环框架你可以把它当成模板逐步补全里面的函数% 初始化变量和参数 x initial_point(); lambda rand(n_eq, 1); % 等式约束对偶变量 mu 0.1 * ones(n_ineq, 1); % 不等式约束对偶变量 tau 0.1; sigma 0.2; maxIter 100; tol 1e-6; gap 1; iter 0; while gap tol iter maxIter iter iter 1; % 1. 计算目标函数梯度、等式约束残差、不等式约束残差 [f, grad_f, c, A, J] evaluate_function(x); % 2. 组装KKT残差向量 r_x grad_f J * lambda A * mu; r_c c; r_s ... % 互补松弛残差 % 3. 组装牛顿系统求解修正方向 K assemble_kkt_matrix(J, A, H, slack); rhs [r_x; r_c; r_s]; delta -K \ rhs; % 4. 计算步长保证松弛变量和对偶变量为正 alpha compute_step(delta, slack, mu, 0.9995); % 5. 更新变量 x x alpha * delta.x; lambda lambda alpha * delta.lambda; mu mu alpha * delta.mu; slack slack alpha * delta.slack; % 6. 计算互补间隙更新障碍参数 gap mu * slack; tau sigma * gap / n_ineq; end这里最需要注意的是修正方程的组装。KKT矩阵K是一个分块稀疏矩阵如果你直接暴力塞满整个矩阵14节点没问题但如果你以后想扩展到更大系统最好从一开始就用Matlab的sparse函数构造稀疏矩阵。等式的雅可比J、不等式的雅可比A都是稀疏的海森矩阵H更是只有少数位置非零用稀疏存储能省一大半内存。3.3 写代码时最容易踩的索引坑我在写这个程序时调试时间最多的不是算法逻辑而是变量索引对不上。OPF的变量向量x是一个大列向量里面混着电压幅值、相角、有功出力、无功出力。假如你不把每个子向量的位置定义清楚后面写雅可比和残差时很容易错位而且这种错位通常是静默的程序不会报错结果就是算出来的东西完全不对。我的建议是在文件开头用常量定义索引区间不要写一堆魔法数字。比如nbus 14; ng 5; iV 1 : nbus; itheta nbus1 : 2*nbus; iPg 2*nbus1 : 2*nbusng; iQg 2*nbusng1 : 2*nbus2*ng;后面所有代码都用x(iV)、x(iPg)这样的方式取变量而不是靠记忆第几个位置是什么。这样哪怕系统规模改了只需要改nbus和ng索引区间自动跟着变。再配合把目标函数、约束函数、雅可比矩阵分别写成独立函数每个函数测试完毕再组合调试效率会高很多。4. 结果验证怎么确认算得对、算得好4.1 先看最优解的定性特征程序跑通之后第一步不是急着看总费用而是看最优解符不符合物理直觉。最常用的判断工具是等微增率准则在忽略网络损耗的理想情况下所有没有达到出力上限的机组边际成本应该相等本文还有配套的精品资源点击获取 基于内点法的Matlab最优潮流实现：以IEEE 14节点系统为例
简介本资源面向电力系统专业本科生、研究生及从事最优潮流研究的工程师提供基于Matlab内点法求解IEEE 14节点系统最优潮流OPF的完整实现方案聚焦燃料费用最小化这一典型经济调度目标。压缩包共5个文件45KB含3个关键参数文本文件节点、支路与发电机数据、1个核心Matlab求解脚本封装fmincon调用与内点算法配置及1张14节点系统拓扑图便于快速建模与结果可视化验证。已有391人学习下载资源结构简洁实用脚本直接调用内点法处理非线性约束优化问题参数文件采用标准格式便于替换扩展拓扑图辅助理解网络结构与潮流分布。读者可即刻运行复现最优发电出力分配、总燃料成本及各节点电压/功率结果掌握OPF建模要点、Matlab优化工具箱实战技巧及电力系统经济调度的核心分析逻辑。 做电力系统优化的同学应该都绕不开最优潮流这块硬骨头。我刚入门的时候天天用Matpower点一下runopf觉得OPF也不过如此后来被要求“不许调库自己写一个内点法求IEEE 14节点系统最优潮流目标函数是燃料费用最小”这才意识到那个“一键出结果”背后全是数学和工程细节。这篇博客就用一篇实战笔记的形式讲清楚怎么从零写一个Matlab内点法程序求解14节点系统的最优潮流。整个过程按“问题建模、算法原理、代码实现、结果验证、踩坑记录”五步展开。适合正在做毕业设计、复现论文方法或者想把OPF底层逻辑吃透的同行参考。1. 先建模14节点系统最优潮流的问题定义与数学形式1.1 为什么要选IEEE 14节点系统在电力系统潮流与优化研究里IEEE 14节点系统是出现频率最高的“入门级标准算例”。它规模不大数据公开又有足够的复杂度来体现算法特性系统里一共14条母线、20条支路、5台发电机总负荷约259MW基准容量100MVA。它包含了普通PQ节点、PV节点、平衡节点还有变压器支路和并联电容器基本覆盖了OPF建模会遇到的大部分设备类型。选它做内点法验证有一个很现实的好处系统规模小状态变量只有几十个雅可比矩阵和海森矩阵都是几十维即使用最笨的稠密矩阵也能在零点几秒内跑完。这让你可以把精力集中在算法逻辑本身而不是矩阵稀疏化和大规模求解上。等你在14节点上把内点法调通了再换IEEE 30节点、118节点核心代码几乎不用改只需要替换数据文件这种“可迁移性”是很值的。另外14节点系统的基准数据在各种工具里都能直接找到尤其是Matpower的case14已经内置了母线、支路、发电机参数。你可以用自己的内点法结果和Matpower的opf结果做交叉验证这是判断程序正确性最快的方法。很多论文里做OPF算法对比也是拿14节点做第一个数值实验因为它能快速暴露算法实现的问题又不会因为系统规模掩盖数值细节。1.2 决策变量、目标函数与机组参数最优潮流本质上是一个带约束的非线性规划问题。目标是最小化系统总燃料费用也就是所有发电机的发电成本之和每台发电机的燃料费用通常用有功出力的二次函数表示min F Σ_i ( a_i * P_gi^2 b_i * P_gi c_i )其中P_gi是第i台发电机的有功出力单位MWa_i、b_i、c_i是费用系数。a_i反映二次成本b_i是一次成本系数c_i是空载费用常数。c_i不影响最优解的取值位置因为常数项的梯度为零但会影响总费用数值代码里一般保留。决策变量分两类一类是节点状态变量包括每个节点的电压幅值V_i和相角θ_i另一类是控制变量包括每台发电机的有功出力P_gi和无功出力Q_gi。在14节点系统里如果不把平衡节点单独处理状态变量总数是14个电压幅值、14个电压相角、5个有功出力、5个无功出力一共38个变量。实际编程时往往固定平衡节点相角为0把它的相角排除在变量之外具体看实现习惯。我这里采用的是一套教学用的典型二次费用系数如下表机组所在母线abcPmax(pu)Pmin(pu)G110.001015.0103.320.10G220.001214.5101.400.10G330.002018.0101.000.05G460.002018.0101.000.05G580.001513.0101.000.05注意这里的a、b对应的P单位是MW也就是把标幺值出力P_pu乘上基准容量100MVA之后代入费用函数。有的文献喜欢在标幺值下写费用函数那a、b的数值会有另一个量级这是正常的。编程前先把单位统一否则后面怎么算都不对。这套系数不是Matpower内置的默认系数Matpower的case14默认用的是线性成本模型需要自己替换成二次多项式。1.3 等式约束与不等式约束怎么列约束条件是OPF和普通经济调度的最大区别。OPF必须满足交流潮流方程也就是每个节点的有功和无功都要平衡。对每个节点i有功等式约束可以写成P_i(V, θ) - P_gi P_di 0其中P_i是注入功率P_di是该节点负荷。无功等式约束同理Q_i(V, θ) - Q_gi Q_di 0这里的V和θ是所有节点的电压向量功率注入方程和牛顿-拉夫逊潮流里的公式完全一样这部分代码可以直接复用潮流计算的mismatch函数。不等式约束则包括四类发电机有功出力上下限、发电机无功出力上下限、节点电压幅值上下限、支路传输功率上限。电压幅值一般限制在0.94到1.06pu之间支路功率约束可以选择视在功率或电流限值看你的数据源怎么定义。每条不等式约束后续都需要引入松弛变量把“小于等于”转成“等式加非负约束”这是内点法能处理它的关键。在14节点系统里不等式约束数量不算多但类型已经足够丰富。建模时我建议把所有不等式约束统一排序统一编号记录每个约束对应的上下限、对应变量索引后面组装KKT系统时就会省很多事。千万不要把约束分散写在不同的代码块里那样调试的时候你会疯掉的。2. 内点法核心逻辑为什么它能高效求解OPF2.1 从“障碍”两个字理解内点法内点法这个名字里最关键的词是“内点”。它不像单纯形法那样沿着可行域的边界找顶点也不像罚函数法那样靠不断试探惩罚系数逼近最优而是从一个严格在可行域内部的点出发靠“障碍项”把迭代点牢牢挡在可行域内部然后逐步逼近边界上的最优点。这个“障碍项”怎么理解呢你可以想象一道物理围墙墙的位置就是不等式约束的边界。优化变量越靠近墙障碍函数的值就越大趋于无穷大这样迭代点就不敢越界。但如果我们始终不许碰墙那最后怎么到达最优解呢方法就是引入一个不断缩小的障碍参数μ。μ大的时候墙离真实边界很远解被压在内部很安全但离最优远μ逐渐变小墙慢慢向真实边界移动解也一步步逼近真正的最优解。当μ趋于0时障碍项消失解就落在边界上了。从数学上看原问题经过松弛和对数变换后变成min f(x) - μ * Σ ln(s_i) s.t. c(x) 0 g(x) s h s 0这里的s是松弛变量。因为ln(s)在s接近0时趋向负无穷等于给目标函数加了一道“软墙”。原对偶内点法的核心就是不断求解这个带μ的系统的KKT条件同时更新μ让解逐渐逼近原问题的最优解。2.2 原对偶内点法的迭代主干要写代码光理解障碍函数还不够还得知道每次迭代具体干什么。原对偶内点法的基本步骤可以用一句话概括把KKT条件写成一组非线性方程组然后用牛顿法迭代求解。具体来说每次迭代做这几件事根据当前变量x、对偶变量λ和μ这里μ指不等式约束对应的对偶变量注意和障碍参数的符号区分计算KKT残差向量。组装牛顿修正方程也就是KKT矩阵这个矩阵是稀疏对称的。求解线性方程组得到变量和对偶变量的修正方向。计算步长保证松弛变量和对偶变量在更新后仍然严格大于0。沿修正方向更新所有变量计算互补间隙gap。如果gap和残差都小于收敛阈值停止迭代否则缩小障碍参数μ继续。这里的“互补间隙”是内点法里最核心的收敛指标。它反映的是原始可行性、对偶可行性以及互补松弛条件的满足程度。每次迭代结束后我们计算gap μ的对偶变量和松弛变量内积之和μ_new σ * gap / n_ineq其中n_ineq是不等式约束个数。当gap小于1e-6时说明已经逼近最优解了。牛顿修正方程的组装是编程里最重的部分。它的矩阵结构大概是[ H Jc A ] [ Δx ] [ r_x ] [ Jc 0 0 ] [ Δλ ] [ r_c ] [ A 0 -Σ ] [ Δμ ] [ r_μ ]H是拉格朗日函数对x的海森矩阵Jc是等式约束雅可比A是不等式约束雅可比Σ是对角阵。14节点系统规模小直接用Matlab的backslash解这个方程就行系统大了再去考虑稀疏分解和预处理那是另一个话题。2.3 三个关键参数μ、σ、步长系数内点法工程实现里最影响成败的是三个参数障碍参数初始值μ0、中心参数σ、步长缩放系数α_max。μ0一般取0.1左右。如果取得太大前几次迭代会花很多时间在“消除障碍”上取得太小初始点可能离可行域太远线性化偏差大容易发散。中心参数σ是每次迭代更新μ时用的比例因子通常在0.1到0.2之间。σ越小μ衰减越快算法收敛快但数值上更容易震荡σ越大收敛越稳但迭代次数变多。我自己的经验是14节点系统取σ0.2稳定和速度都比较均衡。步长计算是内点法一个特别容易出问题的地方。虽然牛顿方向是朝着最优解走的但如果一步跨太大松弛变量可能变成负数违反s0的约束。所以必须算一个“最大安全步长”对修正方向中所有负分量求出允许的最大步长然后乘以一个小于1的安全系数。这个系数我一般取0.9995也就是只走最大安全步长的99.95%留一点余量。如果你设置成1.0迭代点会正好落在边界上导致下一步障碍函数无穷大直接崩掉。3. Matlab代码实现从初始化到收敛的完整流程3.1 数据准备直接基于case14改写代码第一步不是写主函数而是把系统数据准备好。我强烈建议不要自己手敲14条母线的数据直接读取Matpower的case14然后覆盖掉它的发电机费用部分。这样做的好处很明显case14是经过无数人验证的标准数据母线、支路参数不会有错你只需要关注OPF算法本身。具体操作分两步。第一步用Matpower提供的loadcase函数把case14数据读进来然后提取bus、branch、gen三个结构体。第二步把gen里费用相关的列换成我们自己要用的二次系数。Matpower的gen矩阵前6列是出力限值后续列是参与因子和启停状态费用系数通常要另外定义或者直接重写一个gen结构体。我这里给出一个比较清晰的数据定义方式% 系统基准容量 baseMVA 100; % 发电机数据[母线号, Pmax(pu), Pmin(pu), Qmax(pu), Qmin(pu), a, b, c] % 费用函数cost a * P_mw^2 b * P_mw c gen [ 1 3.32 0.10 0.98 -0.20 0.0010 15.0 10; 2 1.40 0.10 0.48 -0.40 0.0012 14.5 10; 3 1.00 0.05 0.40 -0.20 0.0020 18.0 10; 6 1.00 0.05 0.24 -0.06 0.0020 18.0 10; 8 1.00 0.05 0.24 -0.06 0.0015 13.0 10; ]; % 电压限值 Vmin 0.94 * ones(14, 1); Vmax 1.06 * ones(14, 1);母线数据和支路数据直接从case14里取不要自己手抄。mpc loadcase(case14); mpc.gen gen;然后用mpc.bus和mpc.branch作为后续计算的输入。3.2 主循环与修正方程组装主循环是程序的骨架核心逻辑其实不到一百行。下面这段是内点法的主循环框架你可以把它当成模板逐步补全里面的函数% 初始化变量和参数 x initial_point(); lambda rand(n_eq, 1); % 等式约束对偶变量 mu 0.1 * ones(n_ineq, 1); % 不等式约束对偶变量 tau 0.1; sigma 0.2; maxIter 100; tol 1e-6; gap 1; iter 0; while gap tol iter maxIter iter iter 1; % 1. 计算目标函数梯度、等式约束残差、不等式约束残差 [f, grad_f, c, A, J] evaluate_function(x); % 2. 组装KKT残差向量 r_x grad_f J * lambda A * mu; r_c c; r_s ... % 互补松弛残差 % 3. 组装牛顿系统求解修正方向 K assemble_kkt_matrix(J, A, H, slack); rhs [r_x; r_c; r_s]; delta -K \ rhs; % 4. 计算步长保证松弛变量和对偶变量为正 alpha compute_step(delta, slack, mu, 0.9995); % 5. 更新变量 x x alpha * delta.x; lambda lambda alpha * delta.lambda; mu mu alpha * delta.mu; slack slack alpha * delta.slack; % 6. 计算互补间隙更新障碍参数 gap mu * slack; tau sigma * gap / n_ineq; end这里最需要注意的是修正方程的组装。KKT矩阵K是一个分块稀疏矩阵如果你直接暴力塞满整个矩阵14节点没问题但如果你以后想扩展到更大系统最好从一开始就用Matlab的sparse函数构造稀疏矩阵。等式的雅可比J、不等式的雅可比A都是稀疏的海森矩阵H更是只有少数位置非零用稀疏存储能省一大半内存。3.3 写代码时最容易踩的索引坑我在写这个程序时调试时间最多的不是算法逻辑而是变量索引对不上。OPF的变量向量x是一个大列向量里面混着电压幅值、相角、有功出力、无功出力。假如你不把每个子向量的位置定义清楚后面写雅可比和残差时很容易错位而且这种错位通常是静默的程序不会报错结果就是算出来的东西完全不对。我的建议是在文件开头用常量定义索引区间不要写一堆魔法数字。比如nbus 14; ng 5; iV 1 : nbus; itheta nbus1 : 2*nbus; iPg 2*nbus1 : 2*nbusng; iQg 2*nbusng1 : 2*nbus2*ng;后面所有代码都用x(iV)、x(iPg)这样的方式取变量而不是靠记忆第几个位置是什么。这样哪怕系统规模改了只需要改nbus和ng索引区间自动跟着变。再配合把目标函数、约束函数、雅可比矩阵分别写成独立函数每个函数测试完毕再组合调试效率会高很多。4. 结果验证怎么确认算得对、算得好4.1 先看最优解的定性特征程序跑通之后第一步不是急着看总费用而是看最优解符不符合物理直觉。最常用的判断工具是等微增率准则在忽略网络损耗的理想情况下所有没有达到出力上限的机组边际成本应该相等本文还有配套的精品资源点击获取