简介本资源是面向电力系统专业本科生、研究生及优化算法初学者的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求解方法放在一起掂量一下。不同方法的适用场景差异很大选错了后面会非常痛苦。方法基本思路优点缺点本题适用性经典等微增率法只考虑有功平衡按机组边际成本等值分配简单、适合纯经济调度无法处理电压约束、线路潮流约束不适用线性规划方法将潮流方程和目标函数线性化计算快、理论上能处理大系统线性化误差不可控电压无功问题容易失真一般遗传/粒子群等智能算法随机搜索可行域实现简单、无需梯度信息收敛慢、结果不稳定无法保证KKT最优性不推荐原对偶内点法在可行域内部迭代到KKT点对不等式约束天然友好迭代次数与系统规模弱相关需要高质量的梯度和海森矩阵非常合适从这张表可以看出来内点法的核心优势在于它对不等式约束的处理方式。OPF问题有一堆不等式约束发电机有功和无功上下限、节点电压上下限、支路视在功率上限等。经典优化方法要么先猜测哪些约束起作用要么把问题近似成线性问题牺牲精度而内点法通过在对数障碍函数中平滑地处理所有不等式直接在连续空间中迭代精度高且稳健。实际上Matpower中的默认求解器MIPSMatpower 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节点系统 (N14)(NG5)所以变量总数为 (14145538) 个去掉参考节点相角固定为0实际自由变量是37个。目标函数是燃料费用最小标准做法是用二次函数拟合机组耗量特性曲线[ f(P_g) \sum_{i1}^{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是为了防止更新后变量恰好落在边界上导致下一步出现数值奇异。这个“留一点余量”的经验在几乎所有内点法实现里都在用。完整迭代流程可以总结为初始化所有原始变量和对偶变量计算潮流方程残差 (r_P)、(r_Q)不等式约束残差以及互补残差判断是否满足收敛条件对偶间隙小于阈值且最大等式残差小于阈值如果不收敛组装KKT矩阵求解牛顿方向计算步长 (\alpha_p)、(\alpha_d)更新变量更新障碍参数 (\mu)回到第2步。在这个循环中第4步是计算量最大、最容易出错的地方尤其是雅可比和海森矩阵的组装。4. Matlab程序的数据组织与核心代码实现4.1 数据准备与变量索引要在Matlab里实现内点法OPF第一步是把14节点系统的数据读进来并建立统一的变量索引。我最开始写程序的时候没有做索引表结果改一个约束要翻半天代码后来才体会到良好的数据组织和索引设计有多重要。如果使用Matpower数据文件可以用如下方式获取Ybusmpc 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 nbus1:2*nbus; idx.Pg 2*nbus1:2*nbusng; idx.Qg 2*nbusng1:2*nbus2*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数据得到的一组典型结果如下注具体数值会因数据文件中成本系数和约束条件版本不同而有差异重要的是程序逻辑和验证方法母线(P_g) (MW)(Q_g) (MVAr)(V) (pu)1194.2-16.11.06020.027.41.045360.030.21.01060.018.91.02380.012.51.040系统总有功负荷为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会自动选择直接法或迭代法。如果矩阵规模超过几万阶可以考虑本文还有配套的精品资源点击获取
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求解方法放在一起掂量一下。不同方法的适用场景差异很大选错了后面会非常痛苦。方法基本思路优点缺点本题适用性经典等微增率法只考虑有功平衡按机组边际成本等值分配简单、适合纯经济调度无法处理电压约束、线路潮流约束不适用线性规划方法将潮流方程和目标函数线性化计算快、理论上能处理大系统线性化误差不可控电压无功问题容易失真一般遗传/粒子群等智能算法随机搜索可行域实现简单、无需梯度信息收敛慢、结果不稳定无法保证KKT最优性不推荐原对偶内点法在可行域内部迭代到KKT点对不等式约束天然友好迭代次数与系统规模弱相关需要高质量的梯度和海森矩阵非常合适从这张表可以看出来内点法的核心优势在于它对不等式约束的处理方式。OPF问题有一堆不等式约束发电机有功和无功上下限、节点电压上下限、支路视在功率上限等。经典优化方法要么先猜测哪些约束起作用要么把问题近似成线性问题牺牲精度而内点法通过在对数障碍函数中平滑地处理所有不等式直接在连续空间中迭代精度高且稳健。实际上Matpower中的默认求解器MIPSMatpower 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节点系统 (N14)(NG5)所以变量总数为 (14145538) 个去掉参考节点相角固定为0实际自由变量是37个。目标函数是燃料费用最小标准做法是用二次函数拟合机组耗量特性曲线[ f(P_g) \sum_{i1}^{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是为了防止更新后变量恰好落在边界上导致下一步出现数值奇异。这个“留一点余量”的经验在几乎所有内点法实现里都在用。完整迭代流程可以总结为初始化所有原始变量和对偶变量计算潮流方程残差 (r_P)、(r_Q)不等式约束残差以及互补残差判断是否满足收敛条件对偶间隙小于阈值且最大等式残差小于阈值如果不收敛组装KKT矩阵求解牛顿方向计算步长 (\alpha_p)、(\alpha_d)更新变量更新障碍参数 (\mu)回到第2步。在这个循环中第4步是计算量最大、最容易出错的地方尤其是雅可比和海森矩阵的组装。4. Matlab程序的数据组织与核心代码实现4.1 数据准备与变量索引要在Matlab里实现内点法OPF第一步是把14节点系统的数据读进来并建立统一的变量索引。我最开始写程序的时候没有做索引表结果改一个约束要翻半天代码后来才体会到良好的数据组织和索引设计有多重要。如果使用Matpower数据文件可以用如下方式获取Ybusmpc 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 nbus1:2*nbus; idx.Pg 2*nbus1:2*nbusng; idx.Qg 2*nbusng1:2*nbus2*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数据得到的一组典型结果如下注具体数值会因数据文件中成本系数和约束条件版本不同而有差异重要的是程序逻辑和验证方法母线(P_g) (MW)(Q_g) (MVAr)(V) (pu)1194.2-16.11.06020.027.41.045360.030.21.01060.018.91.02380.012.51.040系统总有功负荷为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会自动选择直接法或迭代法。如果矩阵规模超过几万阶可以考虑本文还有配套的精品资源点击获取
