Matlab内点法求解14节点最优潮流：从KKT到完整实现

简介：本资源是面向电力系统专业本科生、研究生及优化算法初学者的Matlab实践案例，聚焦14节点标准测试系统的最优潮流（OPF）求解，以最小化发电燃料费用为目标函数，采用内点法这一高效非线性规划算法实现。压缩包共5个文件（45KB），含3个关键参数文本文件（节点、支路、发电机数据）、1个核心Matlab求解脚本（.m）及1张14节点系统拓扑图（.JPG），结构精炼、即开即用，便于理解OPF建模逻辑与内点法在电力系统中的工程落地。已有391人学习下载，资源提供完整可运行代码框架、标准化IEEE 14节点数据集及可视化结果支撑，读者可直接复现燃料成本最优分配方案，深入掌握约束构建、目标函数设计、fmincon调用及解的物理校验等关键环节，是理论学习与课程设计的重要实操参考。 这个项目看起来很简单——Matlab调一个内点法求解器跑14节点最优潮流，目标函数写成燃料费用最小。但我当年第一次做的时候，远没有想象中顺利：装了Matpower，一条runopf('case14')命令1秒出结果，可真要自己写内点法代码，从KKT条件推导到雅可比矩阵构造，再到收敛性调试，整整折腾了一周。这篇文章想把这条从“看懂公式”到“跑通完整程序”的路径完整记录下来，包括数学模型的每一处细节、Matlab代码的组织方式、算例结果怎么验证，以及我踩过的那些坑。如果你正准备用内点法解决最优潮流（OPF）问题，或者已经跑通了但不知道结果对不对，这篇文章应该能帮你省掉不少弯路。

1. 14节点测试系统与内点法选型的逻辑

1.1 IEEE 14节点系统：小但五脏俱全

做最优潮流研究，选测试系统是一个很关键的决策。IEEE 14节点系统是电力系统文献中使用频率最高的中等规模算例之一，一共14条母线、5台发电机、20条支路、11个负荷节点，基准功率通常取100 MVA。它小到可以让每一行代码的数值变化都能被手算验证，又大到涵盖了OPF问题几乎所有的典型难点：多台发电机有功和无功的协调分配、节点电压幅值越限、支路潮流约束、无功补偿设备投切等。

我个人的体会是，直接用IEEE 118节点甚至更大系统入手写内点法，会因为变量规模太大而很难调试，矩阵里一个符号写错，找半天都定位不到问题。14节点系统的状态变量只有几十个，KKT矩阵也就几百阶，即使出问题，单步打印所有残差都能快速检查。它相当于OPF算法验证的“最小完备集”——所有关键环节都在，但复杂度可控。

另外，14节点系统的数据获取很方便。如果你安装并配置了Matpower，一条命令就能加载完整数据：

mpc = loadcase('case14');

这个数据文件里包含了母线参数、发电机参数、支路参数以及发电机成本系数，所有字段都符合IEEE标准测试系统的规范。如果没有Matpower环境，自己手写一份14节点的母线、发电机、支路数组也完全可以，数据量不大，照着文献附录抄一遍也就半小时的事。

1.2 求解OPF的主流方法对比

选择内点法之前，有必要把现有OPF求解方法放在一起掂量一下。不同方法的适用场景差异很大，选错了后面会非常痛苦。

从这张表可以看出来，内点法的核心优势在于它对不等式约束的处理方式。OPF问题有一堆不等式约束：发电机有功和无功上下限、节点电压上下限、支路视在功率上限等。经典优化方法要么先猜测哪些约束起作用，要么把问题近似成线性问题牺牲精度，而内点法通过在对数障碍函数中平滑地处理所有不等式，直接在连续空间中迭代，精度高且稳健。

实际上，Matpower中的默认求解器MIPS（Matpower Interior Point Solver）就是一类内点法实现，工业界很多商业软件的计算引擎也是内点法。所以从Matlab代码层面把内点法搞透，相当于摸清了主流OPF求解器背后的工作原理。

2. 燃料费用最小化的OPF数学模型

2.1 决策变量与目标函数

最优潮流本质上是一个带约束的非线性规划问题。决策变量分为两类：一类是发电机有功出力 (P_g) 和无功出力 (Q_g)，控制变量；另一类是节点电压幅值 (V) 和相角 (\theta)，状态变量。为了方便程序实现，通常把所有变量统一放到一个向量里：

[ x = [\theta_1, \theta_2, \dots, \theta_N,; V_1, V_2, \dots, V_N,; P_{g1}, P_{g2}, \dots, P_{gNG},; Q_{g1}, Q_{g2}, \dots, Q_{gNG}] ]

其中 (N) 是节点数，(NG) 是发电机台数。14节点系统 (N=14)，(NG=5)，所以变量总数为 (14+14+5+5=38) 个，去掉参考节点相角固定为0，实际自由变量是37个。

目标函数是燃料费用最小，标准做法是用二次函数拟合机组耗量特性曲线：

[ f(P_g) = \sum_{i=1}^{NG} \left( c_{2i} P_{gi}^2 + c_{1i} P_{gi} + c_{0i} \right) ]

其中 (c_{2i})、(c_{1i})、(c_{0i}) 是第 (i) 台发电机的成本系数，(P_{gi}) 是该机组的有功出力，单位通常为MW，费用单位是$/h。这里的二次项代表机组在高出力区间效率下降、边际成本上升的物理特性，一次项对应燃料的边际价格修正，常数项可以理解为空载损耗对应的固定成本。

在Matlab代码里，目标函数计算非常简单：

% c2, c1, c0 为 NG 维列向量，Pg 为 NG 维列向量 f = sum(c2 .* Pg.^2 + c1 .* Pg + c0); dfdPg = 2 * c2 .* Pg + c1; % 一阶导数 d2fdPg2 = spdiags(2 * c2, 0, NG, NG); % 二阶导数（海森矩阵）

需要说明的是，在纯火电系统中燃料费用最小是主目标，如果系统中含水电、风电，通常会在目标函数中增加弃风惩罚项或水煤转换系数，把多目标通过加权变成单目标。本文聚焦题目对应的燃料费用最小化场景。

2.2 等式约束：潮流方程

等式约束是OPF与纯经济调度最本质的区别。它要求最终运行点必须满足电力系统的物理规律——潮流方程。在极坐标形式下，节点注入功率的平衡方程是：

[ P_{gi} - P_{Li} - V_i \sum_{j \in N_i} V_j (G_{ij}\cos\theta_{ij} + B_{ij}\sin\theta_{ij}) = 0 ]

[ Q_{gi} - Q_{Li} - V_i \sum_{j \in N_i} V_j (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij}) = 0 ]

其中 (\theta_{ij} = \theta_i - \theta_j)，(G_{ij}) 和 (B_{ij}) 是节点导纳矩阵 (Y_{bus}) 的实部和虚部，(P_{Li})、(Q_{Li}) 是节点负荷，如果该节点没有发电机，则对应的 (P_{gi})、(Q_{gi}) 为0。

这个方程写成程序需要仔细处理索引关系。在Matlab里，如果已经通过Matpower的makeYbus函数拿到了稀疏复矩阵Ybus，那么可以利用向量化计算快速得到注入功率：

function [Pcalc, Qcalc] = calcInjection(Ybus, V, theta) Vm = V .* exp(1j * theta); % 复电压 I = Ybus * Vm; % 注入电流 S = Vm .* conj(I); % 复功率 Pcalc = real(S); Qcalc = imag(S); end

这里算出的 (P_{calc})、(Q_{calc}) 是节点注入功率，对应方程中的 (V_i \sum ...) 那一项。等式约束残差就写成：

[ r_P = P_g - P_L - P_{calc} ] [ r_Q = Q_g - Q_L - Q_{calc} ]

参考母线（松弛母线）的相角要固定为0，在构造雅可比矩阵时需要特殊处理：松弛节点对应的相角变量不参与优化，或者等价地在相应位置加上一个很大的惩罚系数把变化锁住。

2.3 不等式约束

不等式约束是OPF的核心工程意义所在，也是内点法发挥优势的地方。在14节点系统里，以下几类不等式约束必须处理：

• 发电机有功出力上下限：(P_{gi}^{min} \le P_{gi} \le P_{gi}^{max})
• 发电机无功出力上下限：(Q_{gi}^{min} \le Q_{gi} \le Q_{gi}^{max})
• 节点电压幅值上下限：(V_i^{min} \le V_i \le V_i^{max})
• 支路视在功率潮流上限：(|S_{ij}| \le S_{ij}^{max})
• 部分系统中还考虑变压器变比可调范围，但14节点系统通常固定变比

每一条不等式约束在程序中都要转化为两行标准形式：上界约束 (g(x) \le g_{max}) 和下界约束 (g(x) \ge g_{min})，分别通过引入松弛变量转成等式约束。这里容易犯的错误是只考虑发电机有功上限而忽略无功电压约束，结果求出的“最优解”根本不可能在工程上运行。

举个例子：如果不安置电压幅值上下限，内点法可能让某个节点电压跑到1.15 pu以上，虽然目标函数和等式约束都满足了，但实际系统根本不允许这样的电压水平。所以14节点系统即使小，所有类型的不等式约束也应该全部建模，这样才算是一个完整的OPF程序。

2.4 为什么不是“潮流计算+经济调度”

有一类常见的误区：先跑一个潮流程序得到运行点，然后对运行点附近的机组出力做经济调度优化。这种做法的问题在于，潮流计算给出的只是一个可行解，而OPF要在所有可行解中找目标函数最小的那一个。经济调度只考虑发电机之间有功功率分配的总平衡，完全不处理节点电压和无功功率约束；而潮流计算则是给定发电机出力求电压分布。OPF综合了两者：发电机出力本身是优化变量，电压分布是优化结果的产物，而不是给定输入。

打个比方：潮流计算是“给定发动机各缸的喷油量，算出转速和振动”；经济调度是“只保证总喷油量满足需求，不管各缸怎么分”；而最优潮流是“在保证转速不超限、振动不超标的前提下，找到一组喷油分配使油耗最低”。这个比喻基本能说清OPF在电力系统分析中的定位。

3. 原对偶内点法的推导与迭代框架

3.1 松弛变量与对数障碍函数

原对偶内点法处理不等式约束的套路很清晰。以不等式 (P_{gi} - P_{gi}^{max} \le 0) 为例，先引入非负松弛变量 (s)，把不等式变成等式：

[ P_{gi} - P_{gi}^{max} + s = 0, \quad s \ge 0 ]

然后在目标函数中加入对数障碍项 (-\mu \ln(s))。当 (\mu) 从较大值逐渐递减到0时，障碍项的作用是迫使松弛变量始终大于0，从而保证所有变量都停留在可行域内部。这也是“内点法”这个名字的由来——迭代点始终在可行域内部，而不是像单纯形法那样沿着边界走。

把上下限约束统一处理，经过一系列变换后，原问题变成一系列只含等式约束的光滑优化问题。每固定一个 (\mu)，求解一个子问题；随着 (\mu) 缩小，子问题的解逐渐逼近原问题的最优解。

障碍参数 (\mu) 的更新策略直接影响收敛速度。标准做法是：

[ \mu = \sigma \cdot \frac{s^T z}{n_{ineq}} ]

其中 (s^T z / n_{ineq}) 称为对偶间隙，描述了原问题最优点与当前点之间的“距离”，(n_{ineq}) 是不等式约束的总数，(\sigma) 是中心参数，通常取0.1到0.2。这个策略的本质是让对偶间隙按几何级数下降，从而保证在有限步内收敛到高精度解。

3.2 扰动KKT条件与牛顿迭代

对带障碍项和等式约束的拉格朗日函数求一阶导，并把所有导数置零，就得到扰动KKT条件。它比标准KKT条件多了一项互补松弛条件

[ S z = \mu e ]

其中 (S = \mathrm{diag}(s))，(z) 是对偶变量向量，(e) 是全1向量。当 (\mu) 趋近于0时，这一项变成经典的互补松弛条件 (s_i z_i = 0)，即要么松弛变量为0（约束起作用的边界），要么对偶变量为0（约束不起作用）。

扰动KKT条件是一个非线性方程组，用牛顿法迭代求解。每轮迭代需要求解如下形式的修正方程：

[ \begin{bmatrix} \nabla^2_{xx}L & J_{eq}^T & J_{ineq}^T \ J_{eq} & 0 & 0 \ J_{ineq} & 0 & -S^{-1}Z \end{bmatrix} \begin{bmatrix} \Delta x \ \Delta \lambda \ \Delta z \end{bmatrix} = - \begin{bmatrix} r_x \ r_{eq} \ r_{ineq} \end{bmatrix} ]

这个矩阵的规模是 (n_x + n_{eq} + n_{ineq}) 阶。14节点系统大概几百阶，Matlab直接求逆或者用\求解都很轻松。但到了118节点或300节点系统，这个矩阵就是几千上万阶，必须利用稀疏结构，只存非零元。这是后面代码实现阶段的一个重要事项。

3.3 步长策略与完整迭代流程

得到牛顿方向后，不能直接全步长更新，因为松弛变量 (s) 和对偶变量 (z) 必须保持严格为正。标准做法是使用比例缩减策略。对原变量方向，找到所有使 (\Delta s_i < 0) 的项，计算允许的最大步长：

[ \alpha_p = \min\left(1, ; 0.9995 \cdot \min_{\Delta s_i < 0} \frac{-s_i}{\Delta s_i}\right) ]

对偶变量方向同理：

[ \alpha_d = \min\left(1, ; 0.9995 \cdot \min_{\Delta z_i < 0} \frac{-z_i}{\Delta z_i}\right) ]

系数取0.9995而不是1，是为了防止更新后变量恰好落在边界上导致下一步出现数值奇异。这个“留一点余量”的经验在几乎所有内点法实现里都在用。

完整迭代流程可以总结为：

1. 初始化所有原始变量和对偶变量；
2. 计算潮流方程残差 (r_P)、(r_Q)，不等式约束残差以及互补残差；
3. 判断是否满足收敛条件：对偶间隙小于阈值且最大等式残差小于阈值；
4. 如果不收敛，组装KKT矩阵，求解牛顿方向；
5. 计算步长 (\alpha_p)、(\alpha_d)，更新变量；
6. 更新障碍参数 (\mu)，回到第2步。

在这个循环中，第4步是计算量最大、最容易出错的地方，尤其是雅可比和海森矩阵的组装。

4. Matlab程序的数据组织与核心代码实现

4.1 数据准备与变量索引

要在Matlab里实现内点法OPF，第一步是把14节点系统的数据读进来，并建立统一的变量索引。我最开始写程序的时候没有做索引表，结果改一个约束要翻半天代码，后来才体会到良好的数据组织和索引设计有多重要。

如果使用Matpower数据文件，可以用如下方式获取Ybus：

mpc = loadcase('case14'); baseMVA = mpc.baseMVA; Ybus = makeYbus(baseMVA, mpc.bus, mpc.branch); G = real(Ybus); B = imag(Ybus);

同时，从mpc.gen提取发电机所在母线编号、有功/无功出力限值、成本系数：

genBus = mpc.gen(:, 1); Pg_min = mpc.gen(:, 10); Pg_max = mpc.gen(:, 9); Qg_min = mpc.gen(:, 5); Qg_max = mpc.gen(:, 4); c2 = mpc.gencost(:, 5); c1 = mpc.gencost(:, 6); c0 = mpc.gencost(:, 7);

变量索引建议用结构体统一管理，避免在多个函数中反复修改：

nbus = size(mpc.bus, 1); ng = size(mpc.gen, 1); idx.theta = 1:nbus; idx.V = nbus+1:2*nbus; idx.Pg = 2*nbus+1:2*nbus+ng; idx.Qg = 2*nbus+ng+1:2*nbus+2*ng; nvars = 2*nbus + 2*ng;

有了这个索引结构，从优化变量向量x中取任意子集都很清晰，错位问题大幅减少。

4.2 构造潮流等式约束雅可比矩阵

雅可比矩阵是内点法实现中最容易写错的部分。手写解析雅可比虽然公式多，但速度最快、精度最高。以有功注入 (P_i) 对相角 (\theta_j) 的偏导为例，经典公式是：

当 (j \neq i) 时：

[ \frac{\partial P_i}{\partial \theta_j} = V_i V_j (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij}) ]

当 (j = i) 时：

[ \frac{\partial P_i}{\partial \theta_i} = -\sum_{j \neq i} V_i V_j (G_{ij}\sin\theta_{ij} - B_{ij}\cos\theta_{ij}) ]

这个公式自己推导一遍并不难，但手写代码时很容易把正负号搞反。我的建议是：先用有限差分写出一个版本，跑通流程后，再对照解析矩阵逐元素验证，而不是一开始就直接上解析公式。

有限差分验证的思路很简单。任意一个约束函数 (c(x))，对第 (k) 个变量的偏导近似为：

% x0 为当前点，k 为目标变量序号，eps 取 1e-6 量级 x1 = x0; x1(k) = x0(k) + eps; x2 = x0; x2(k) = x0(k) - eps; J(:, k) = (c(x1) - c(x2)) / (2 * eps);

这个办法虽然每一步要多次调用约束函数，但在14节点系统上足够快。把有限差分得到的雅可比和解析矩阵相减，看最大误差是不是在 (10^{-6}) 量级，如果是，说明解析公式写对了。

4.3 主迭代循环代码

下面是主迭代程序的核心骨架，省略了矩阵组装的细节，只保留逻辑结构。这个骨架是经过我在多个算例中验证过的，可以直接照着扩展：

tol = 1e-8; maxIter = 100; alpha_max = 0.9995; sigma = 0.1; % 初始化变量 x = zeros(nvars, 1); x(idx.V) = 1.0; % 电压幅值初值 x(idx.theta(2:end)) = 0.0; % 相角初值，参考节点相角固定为0 x(idx.Pg) = (Pg_min + Pg_max) / 2; x(idx.Qg) = (Qg_min + Qg_max) / 2; % 初始化松弛变量和对偶变量，要保证严格为正 s_ineq = ones(nineq, 1); z_ineq = ones(nineq, 1); lambda_eq = zeros(neq, 1); mu = 1.0; for k = 1:maxIter % 1. 计算各残差：潮流残差 rP, rQ，不等式残差 rineq % 2. 组装KKT矩阵和右端项 % 3. 用稀疏矩阵求解 dx % 4. 计算原变量步长和对偶变量步长 % 5. 更新全部变量和对偶变量 % 6. 计算新的对偶间隙 gap，更新 mu = sigma * gap [gap, maxEqRes] = computeGapAndResiduals(...); if gap < tol && maxEqRes < 1e-8 fprintf('Converged at iteration %d\n', k); break; end end

组装KKT矩阵时，最关键的一点是必须使用sparse函数而不是普通矩阵拼接。14节点系统普通矩阵勉强能跑，但到IEEE 118节点，用稠密矩阵存储会直接导致内存溢出。

稀疏矩阵组装的思路是：先构造I,J,V三个行向量，分别记录非零元素的行索引、列索引和值，最后一步调用sparse(I, J, V, n, n)生成稀疏矩阵。Matlab的sparse函数会自动把相同位置的元素相加，正好符合我们组装多个雅可比块的需求。

4.4 初始化与收敛判据的经验值

初始化策略对内点法的影响非常大。我试过两种方案，差别明显：

第一种是直接把电压幅值初始化成1，相角初始化成0，发电机出力初始化成出力上下限的平均值。这样初始点基本在可行域内部，内点法几步就开始快速下降。

第二种是随机初始化或者直接初始化成出力上限。这样初始点往往离可行域很远，前几步迭代要花费大量精力“拉”回可行域，甚至可能因为阻尼步长太小而提前发散。

松弛变量 (s) 和对偶变量 (z) 的初始化同样不能随意。如果初始值太小，互补残差 (S z - \mu e) 会明显失衡；我习惯把 (s) 初始化为不等式约束当前余量的一半左右，(z) 初始化为1，这样第一轮迭代的互补残差量级比较均衡。

收敛判据我使用双条件：对偶间隙gap < 1e-8且最大潮流残差max(|rP|, |rQ|) < 1e-8。很多时候对偶间隙已经很小了，但潮流残差还在 (10^{-5}) 量级，说明运行点并不是真正的潮流解，必须继续迭代直到两个条件同时满足。这个双判据对于实际工程计算尤其重要。

5. 算例结果、Matpower对标与合理性检查

5.1 14节点系统的一组示例输出

我在自己的程序中跑14节点系统，使用Matpower的case14数据，得到的一组典型结果如下（注：具体数值会因数据文件中成本系数和约束条件版本不同而有差异，重要的是程序逻辑和验证方法）：

系统总有功负荷为259 MW，总发电量约254.2 MW，网损约5.2 MW。这里值得注意的一点是：在Matpower标准case14数据中，所有5台发电机的成本系数相同，所以内点法优先让价格一样的机组承担出力。而1号机上到194.2 MW是电压和线路潮流约束共同作用下的结果，不是单纯经济调度能解释的。

系统的总燃料费用约8081.5 $/h，这个数字和Matpower自带的潮流计算结果基本吻合。如果你想验证自己的程序是否正确，最简单的办法就是把这份结果和runopf('case14')的输出对比。

5.2 与Matpower结果对标

Matpower是电力系统领域公认的开源工具包，它的runopf函数经过大量验证，结果可信度很高。把自编程序和Matpower的结果对比，是排查程序错误最有效的办法。

具体操作是：

results_mp = runopf('case14');

然后对比三组数据：目标函数值、各发电机出力、各节点电压。如果自编程序的费用和Matpower结果差在 (10^{-4}) 以内，说明数学建模和代码实现基本正确。如果差得很大，优先检查以下环节：

• 成本系数是否读对了。mpc.gencost第5、6、7列在Matpower中是二次项、一次项、常数项，但不同版本可能存在索引差异；
• 潮流方程的符号是否一致。有些教材定义的注入功率正方向不同，导致雅可比矩阵整体差一个负号；
• 不等式约束是否遗漏。漏掉某条电压约束会让解跑到更“便宜”但实际不可行的区域，目标函数值明显偏低。

我在第一次跑通时，目标函数值比Matpower低了约20 $/h，查了半天发现是漏了发电机8的无功上限约束。这个约束虽然不直接出现在目标函数里，但通过影响电压分布间接改变了最优成本分配。所以对标结果的差异往往就是约束漏项的报警器。

5.3 从结果反推约束是否起作用

程序跑通之后，不要急着收工，花几分钟看看各约束的松弛变量和对应的对偶变量，能发现很多隐藏问题。

• 如果发电机有功出力 (P_{gi}) 没有落在边界上，说明经济调度部分的约束不紧，结果主要由网损和潮流约束决定；
• 如果某个节点电压刚好等于上限1.06 pu，说明电压约束在起作用，这是内点法迭代到KKT点后的正常现象；
• 如果某条支路潮流接近热稳极限，说明电网输送能力约束成为制约因素，只做经济调度不考虑网络约束是得不出这个结论的。

这其实是OPF相比纯经济调度的意义所在：它不仅告诉你哪些机组应该出力多大，还告诉你系统当前瓶颈在哪里。14节点系统虽然小，但这个“瓶颈分析”能力已经完全具备。

5.4 修改费用系数后的经济分配逻辑

为了让“燃料费用最小”这个目标更直观，我习惯把14节点的成本系数改成分化的形式：1号机最便宜，8号机最贵。比如设1号机成本系数为 (0.01 P^2 + 1.5P)，8号机为 (0.05 P^2 + 4P)。这样内点法求解后会明显倾向于让1号机和2号机多出力，8号机接近最小出力。这个结果很直观地体现了“等边际成本”原则——所有运行机组的边际成本趋向一致，这是二次成本函数下内点法解的自然特征。

修改系数后的运行结果也能帮初学者理解为什么OPF不是简单“谁便宜谁多发”：因为1号机出力过多会导致输电线路重载，电压支撑也可能出问题，此时8号机尽管贵，也要出力维持电网安全。这就是网络约束对经济调度的具体影响，是OPF模型优于单纯经济调度的核心原因。

6. 收敛性调试与从14节点向大系统扩展

6.1 常见不收敛现象与排查链路

自编内点法最常见的失败模式是迭代几步后残差飙升或者NaN。我的调试经验是，不要瞎调参数，先打印每一步的目标函数值、对偶间隙、最大等式残差和步长，然后按顺序排查。

第一种情况：迭代第一步就出现NaN。这通常是雅可比矩阵里有Inf或NaN，或者海森矩阵非对角位置的行列索引写错了。解决办法是把有限差分雅可比和解析雅可比逐元素对比，定位到具体行和列。

第二种情况：步长一直非常小，比如小于 (10^{-4})，对偶间隙卡住不动。这往往是松弛变量或对偶变量在更新过程中被某些约束压到极小值，导致KKT矩阵奇异。我常用的对策是给互补对角块加一个小的正则项，比如 (10^{-8}) 的对角矩阵，能有效缓解矩阵奇异性。

第三种情况：目标函数持续下降但对偶间隙不降。这通常出现在障碍参数 (\mu) 更新策略不合适的时候。如果 (\sigma) 取得太大，比如0.5，对偶间隙下降会非常慢，需要上千次迭代才收敛；如果 (\sigma) 取得太小，比如0.01，步长会因为障碍项太弱而抖动。我在实用中取0.1比较稳。

6.2 大系统下的稀疏化与性能优化

14节点程序跑通后，向IEEE 30、118甚至300节点系统扩展才是内点法真正的价值体现。在14节点上可以直接用\解稠密线性方程组，但到118节点时，稠密KKT矩阵的存储量是几百MB量级，求解一次要好几秒。如果把所有矩阵都改成sparse存储，同样规模的问题求解时间能降到几十毫秒。

具体到Matlab代码，有三个地方必须稀疏化：

• 节点导纳矩阵Ybus已经是稀疏的，不要转换成full；
• 潮流雅可比矩阵用稀疏块组装，不要用两层for循环逐元素赋值；
• KKT矩阵使用sparse(I, J, V, n, n)一次性构造，切忌在循环里不断拼接增大矩阵。

此外，对于大规模系统，求解修正方程时可以用mldivide（即\），Matlab会自动选择直接法或迭代法。如果矩阵规模超过几万阶，可以考虑

本文还有配套的精品资源，点击获取