基于MATPOWER的牛顿—拉夫逊法解潮流

随着电力系统日益复杂化，准确有效地进行潮流计算成为电力工程领域的一个重要挑战。潮流计算是电力系统分析的核心，关系到电网的稳定性、可靠性以及经济运行。在这一背景下，牛顿-拉夫逊方法因其高效率和精确度而成为潮流计算的主流方法之一。然而，随着电力系统规模的不断扩大和运行环境的不断变化，传统的潮流计算方法面临着新的挑战和需求。

MATPOWER作为一款在电力系统分析和优化领域广泛使用的开源软件工具包，提供了包括牛顿-拉夫逊法在内的多种潮流计算算法。MATPOWER的高灵活性和广泛的功能使其成为研究和工业应用的理想选择。特别是在MATLAB环境下，MATPOWER能够充分利用MATLAB的高级数值计算和可视化功能，从而为电力系统的分析和设计提供强大的支持。

本次仿真实验旨在深入探究和比较MATPOWER中实现的牛顿-拉夫逊方法与自行实现的牛顿-拉夫逊算法。通过在MATLAB环境下进行这一比较，我们期望评估和分析这两种方法在处理复杂电力系统模型时的性能和效率。这包括计算精度、收敛速度、处理大规模问题的能力以及在不同电网条件下的稳定性。通过这一比较，我们希望能够识别各种潮流计算方法的优势和局限，并为未来电力系统的研究和实际应用提供有价值的见解和建议。

四、仿真案例分析

1. 仿真概述

在本次仿真案例中，我们选取了一个九节点的交流电力系统网络作为研究对象。该网络模型被设计为包含发电机、负载、变压器和输电线路等典型电网元件，旨在模拟现实世界中电力系统的基本运作情况。通过MATPOWER在MATLAB环境下进行潮流计算，我们将分析网络中的功率流动、电压分布以及损耗情况，以评估电网的性能和稳定性。

2. 网络结构和参数设置

节点配置：

九个节点包括发电机节点、负荷节点和平衡节点。每个节点都具备相应的电压、功率等参数设置。

线路参数：

网络中的输电线路具备实际线路的电阻、电抗和导纳特性。

发电机和负荷：

设置合适的发电机有功、无功输出和负荷需求，以模拟现实电网运行状态。

3.case9的参数

可以通过MATPOWER中的mpc函数来调取case9的参数，调取代码如下：

mpc = loadcase('case9');
disp('Bus Data:');
disp(mpc.bus);
disp('Generator Data:');
disp(mpc.gen);
disp('Branch Data:');
disp(mpc.branch);
disp('Generation Cost Data:');
disp(mpc.gencost);

运行此代码得出的结果即为case9的参数配置：

对于此结果的解读说明如下表：

五、MATPOWER 解 case9

首先直接利用MATPOWER的内置函数进行case9的求解，代码如下：

clear; close all; clc;

define_constants;  % 加载 MATPOWER 的常量定义
mpc = loadcase('case9');  % 加载 case9 数据

% 配置 MATPOWER 选项，使用牛顿-拉夫逊法
mpopt = mpoption('pf.alg', 'NR', 'verbose', 1);

results = runpf(mpc, mpopt);

% 显示结果
results

下面是对这段代码的解释：

1.初始化命令：

clear;：清除MATLAB工作空间中的所有变量。

close all;：关闭所有打开的图形窗口。

clc;：清除命令窗口中的内容。

2.加载MATPOWER常量：

define_constants;：执行这个命令会加载MATPOWER需要的一些常量定义。

3.加载电力系统案例：

mpc = loadcase('case9');：加载case9的电力系统模型。

4.配置MATPOWER选项：

mpopt = mpoption('pf.alg', 'NR', 'verbose', 1);：这里配置了用于潮流计算的MATPOWER选项。'pf.alg', 'NR'指定使用牛顿-拉夫逊方法（Newton-Raphson）进行潮流计算。'verbose', 1设置输出详细信息，使得计算过程和结果更加清晰可见。

5.运行潮流计算：

results = runpf(mpc, mpopt);：使用前面指定的系统模型（mpc）和配置选项（mpopt）运行潮流计算。计算结果被存储在results变量中。

6.显示结果：

results：这一行代码将打印出潮流计算的结果，包括每个节点的电压、每条线路的功率流动情况以及系统损耗等信息。

运行此代码，得到如下结果：

从结果中我们可以看出，在当前仿真环境下，牛顿-拉夫逊法在四次迭代后收敛，总共用时0.86秒。

同时包括下面所列的重要数据：

1.系统汇总

电网规模：共有9个母线，3个发电机。

发电能力：总发电容量为820 MW，实际发电量为319.6 MW。

负荷：总负荷为315 MW，无功负荷为115 MVAr。

线损：总线损为4.64 MW，无功损失为48.38 MVAr。

电压范围：电压幅值在0.996 p.u.（母线9）到1.040 p.u.（母线1）之间，电压角在-3.99°（母线9）到9.28°（母线2）之间。

2.母线数据

这部分详细列出了每个母线的电压幅值、相角、发电量、负荷等。

3.支路数据

这部分详细描述了电网中每条支路的功率流动情况和线损。

为了方便后面的对比分析，我们着重分析平衡节点功率，线路功率，功率损耗这三个结果：

1.平衡节点功率：

母线1 (平衡节点) 的功率：

有功功率（P）：71.64 MW

无功功率（Q）：27.05 MVAr

2.线路功率：

支路1 (母线1到母线4)：

有功功率流（P）：71.64 MW

无功功率流（Q）：27.05 MVAr

支路7 (母线8到母线2)：

有功功率流（P）：163.00 MW

无功功率流（Q）：6.65 MVAr

支路9 (母线9到母线4)：

有功功率流（P）：40.94 MW

无功功率流（Q）：22.89 MVAr

3.功率损耗：

总功率损耗：

有功损耗总和：4.64 MW

无功损耗总和：48.38 MVAr

具体支路的最大损耗：

支路8-9的有功损耗：2.30 MW

支路8-2的无功损耗：15.83 MVAr

六、纯牛顿-拉夫逊法解算 case9

这一部分我采取了模块化编程的方式，将各个函数封装成模块进行计算，这样的好处在于可以更清晰的找到代码运行过程中的问题所在，方便代码的调试于复现。

首先我先利用流程图的方式来阐释我所编写的代码的内在逻辑关系：

下面是分别展示每一个代码模块：

function [Ybus, V0, P, Q] = extract_data()
    % 加载case9数据
    mpc = loadcase('case9');
    % 获取导纳矩阵
    [Ybus, ~, ~] = makeYbus(mpc.baseMVA, mpc.bus, mpc.branch);

    % 提取电压初始值（幅值和相角）
    V_mag = mpc.bus(:, 8);  % 电压幅值
    V_ang = mpc.bus(:, 9);  % 电压相角（度）
    V0 = V_mag .* exp(1j * deg2rad(V_ang));  % 初始电压矢量

    % 提取有功和无功负载
    P = -mpc.bus(:, 3) / mpc.baseMVA;  % 负载有功
    Q = -mpc.bus(:, 4) / mpc.baseMVA;  % 负载无功

    % 考虑发电机功率注入
    for i = 1:size(mpc.gen, 1)
        bus_idx = mpc.gen(i, 1); % 发电机所在的母线索引
        P_gen = mpc.gen(i, 2) / mpc.baseMVA; % 发电机有功功率
        Q_gen = mpc.gen(i, 3) / mpc.baseMVA; % 发电机无功功率
        P(bus_idx) = P(bus_idx) + P_gen; % 累加发电机有功功率
        Q(bus_idx) = Q(bus_idx) + Q_gen; % 累加发电机无功功率
    end
end

function [lineFlows, losses] = calculate_line_flows(mpc, V, Ybus)
    branch = mpc.branch;
    nBranch = size(branch, 1);
    lineFlows = zeros(nBranch, 4); % 列：发端P, 发端Q, 收端P, 收端Q
    losses = zeros(nBranch, 1);

    for i = 1:nBranch
        from = branch(i, 1);
        to = branch(i, 2);
        Y = Ybus(from, to);
        Ysh = Ybus(from, from) + Ybus(to, to); % 支路自导纳
        Vfrom = V(from);
        Vto = V(to);

        Sfrom = Vfrom * conj((Vfrom - Vto) * (-Y) + Vfrom * Ysh/2);
        Sto = Vto * conj((Vto - Vfrom) * (-Y) + Vto * Ysh/2);

        lineFlows(i, :) = [real(Sfrom) imag(Sfrom) real(Sto) imag(Sto)];
        losses(i) = Sfrom + Sto; % 总损耗 = 发端损耗 + 收端损耗
    end
end

function J = get_jacobian(Ybus, V, P, Q)
    nb = length(V);  % 母线数量
    J1 = zeros(nb, nb);  % 对应于有功/相角的偏导数
    J2 = zeros(nb, nb);  % 对应于有功/电压幅值的偏导数
    J3 = zeros(nb, nb);  % 对应于无功/相角的偏导数
    J4 = zeros(nb, nb);  % 对应于无功/电压幅值的偏导数

    Vm = abs(V);
    Va = angle(V);

    for i = 1:nb
        Vi = V(i);
        Vi_abs = Vm(i);
        Vi_ang = Va(i);

        for k = 1:nb
            Vk = V(k);
            Vk_abs = Vm(k);
            Vk_ang = Va(k);

            if i ~= k
                % 非对角线元素
                Yik = Ybus(i,k);
                Yik_abs = abs(Yik);
                Yik_ang = angle(Yik);

                J1(i,k) = -Vi_abs * Vk_abs * Yik_abs * sin...
                    (Yik_ang + Vk_ang - Vi_ang);
                J2(i,k) =  Vk_abs * Yik_abs * cos...
                    (Yik_ang + Vk_ang - Vi_ang);
                J3(i,k) =  Vi_abs * Vk_abs * Yik_abs * cos...
                    (Yik_ang + Vk_ang - Vi_ang);
                J4(i,k) =  Vk_abs * Yik_abs * sin...
                    (Yik_ang + Vk_ang - Vi_ang);
            end
        end

        % 对角线元素
        Yi = Ybus(i,i);
        Yi_abs = abs(Yi);
        Yi_ang = angle(Yi);

        dPi = -Q(i) - Vi_abs^2 * imag(Yi);
        dQi = P(i) - Vi_abs^2 * real(Yi);

        for k = 1:nb
            if k ~= i
                Vk = V(k);
                Vk_abs = Vm(k);
                Vk_ang = Va(k);
                Yik = Ybus(i,k);
                Yik_abs = abs(Yik);
                Yik_ang = angle(Yik);

                dPi = dPi + Vi_abs * Vk_abs * Yik_abs * sin...
                    (Yik_ang + Vk_ang - Vi_ang);
                dQi = dQi + Vi_abs * Vk_abs * Yik_abs * cos...
                    (Yik_ang + Vk_ang - Vi_ang);
            end
        end

        J1(i,i) = dPi;
        J2(i,i) = -dQi / Vi_abs;
        J3(i,i) = dQi;
        J4(i,i) = -dPi / Vi_abs;
    end

    % 组合雅可比矩阵
    J = [J1 J2; J3 J4];
end

function mis = get_power_mismatch(V, Ybus, P, Q)
    S = V .* conj(Ybus * V);  % 计算复数功率 S = V * I*
    misP = real(S) - P;       % 有功功率不平衡
    misQ = imag(S) - Q;       % 无功功率不平衡
    mis = [misP; misQ];       % 确保这是一个长度为 2*nb 的列向量
end

function [V, success, iter] = newton_raphson_power_flow(Ybus, V0, P, Q)
    max_iter = 50;
    tolerance = 1e-8;
    V = V0;  % 确保 V0 是列向量
    nb = length(V);

    for iter = 1:max_iter
        mis = get_power_mismatch(V, Ybus, P, Q);
        if max(abs(mis)) < tolerance
            success = true;
            return;
        end

        J = get_jacobian(Ybus, V, P, Q);
        deltaV = -J \ mis;

        % 检查 deltaV 的长度
        if length(deltaV) ~= 2*nb
            error('deltaV 的长度与预期不匹配。');
        end

        % 更新 V
        V = V + deltaV(1:nb);  % 使用 deltaV 的前半部分更新 V
    end

    success = false;
end

function main()
    % 直接加载 case9 数据
    mpc = loadcase('case9');
    
    % 提取导纳矩阵、初始电压、功率数据
    [Ybus, V0, P, Q] = extract_data(); 

    % 执行牛顿-拉夫逊潮流计算
    [V, success, iter] = newton_raphson_power_flow(Ybus, V0, P, Q);

    % 显示结果
    if success
        fprintf('潮流计算收敛于 %d 次迭代\n', iter);

        % 计算每个母线的实际功率注入
        Sbus = V .* conj(Ybus * V);
        Pbus = real(Sbus);
        Qbus = imag(Sbus);

        % 计算每条线路的功率流和损耗
        [lineFlows, losses] = calculate_line_flows(mpc, V, Ybus);

        % 显示最终母线电压
        disp('最终母线电压：');
        disp(V);

        % 显示母线功率
        fprintf('母线功率 (P, Q):\n');
        disp([Pbus, Qbus]);

        % 显示线路功率和损耗
        fprintf('线路功率和损耗:\n');
        disp(lineFlows);
        fprintf('总线损: P = %f MW, Q = %f MVar\n', sum(real(losses)), sum(imag(losses)));
    else
        disp('潮流计算未收敛');
    end
end

clear; close all; clc;
main();

下面我将大致解释一下整个工程各个模块的交互逻辑：

1.初始化[clear; close all; clc;]这是程序开始之前的准备阶段，用于清理工作环境，关闭所有图形界面，清空命令窗口。

2.主函数[main()]程序的核心入口点，负责调用其他函数并协调整个程序的运行流程。

3.提取数据[extract_data()]由主函数调用，用于加载和处理电网数据，包括导纳矩阵、初始电压、功率数据。

4.牛顿-拉夫逊功率流 [newton_raphson_power_flow()]同样由主函数调用，这是实现牛顿-拉夫逊方法的核心函数，用于求解电力系统的稳态运行点。

5.计算线路流量[calculate_line_flows()]这个函数由 extract_data 调用，用于计算电网中各输电线路的功率流和损耗。

6.计算电力不平衡与雅各比矩阵[get_power_mismatch(), get_jacobian()]这两个函数被牛顿-拉夫逊功率流函数调用，用于计算电力系统的功率不平衡和构建雅各比矩阵，对算法的迭代求解至关重要。

下面我将给出每个子函数和主函数的变量参数表：

下面我将简单地阐释一下各个代码的执行逻辑：

1.牛顿-拉夫逊功率流算法主函数 (newton_raphson_power_flow)

这个函数是实现牛顿-拉夫逊方法的核心。它的主要步骤如下：

初始化：设定最大迭代次数和容忍度，以及初始电压矢量。

迭代过程：在每次迭代中，先计算电力不平衡量（有功和无功），然后根据这个不平衡量和雅各比矩阵计算电压的修正值。

收敛判断：检查电力不平衡量是否小于预设的容忍度，如果是，则认为算法已收敛，否则继续迭代。

输出：返回最终的电压矢量、算法是否成功收敛的标志和迭代次数。

2.计算电力不平衡的函数 (get_power_mismatch)

这个函数用于计算系统的电力不平衡量，它的工作流程是：

计算复数功率：根据当前的电压矢量和系统的导纳矩阵计算每个节点的复数功率。

计算不平衡量：计算每个节点的有功和无功功率的实际值与设定值之间的差异，这就是不平衡量。

3.计算雅各比矩阵的函数 (get_jacobian)

这个函数负责构建牛顿-拉夫逊方法中使用的雅各比矩阵：

矩阵构建：根据电网的导纳矩阵和当前电压矢量，计算雅各比矩阵的四个子矩阵。

矩阵组合：将这四个子矩阵组合成完整的雅各比矩阵。

4.提取和处理电网数据的函数 (extract_data)

这个函数用于从外部数据源或电力系统模型中提取必要的数据：

数据加载：加载特定的电力系统案例数据。

参数提取：从加载的数据中提取出导纳矩阵、初始电压、功率等关键参数。

5.计算输电线损耗的函数 (calculate_line_flows)

这个函数用于计算电网中各输电线路的功率流和损耗：

线路遍历：遍历每条线路，计算线路两端的功率流。

损耗计算：根据功率流计算每条线路的功率损耗。

6.主函数 (main)

主函数协调整个电力系统潮流计算过程：

数据提取：调用extract_data函数获取电网数据。

潮流计算：使用newton_raphson_power_flow函数执行潮流计算。

结果展示：显示计算结果，包括电压、功率等信息。

7.MATLAB环境初始化命令 (clear; close all; clc; main();)

这一行代码用于准备MATLAB环境并开始执行潮流计算：

环境清理：清除变量，关闭图形界面，清空命令窗口。

执行主函数：开始执行main函数，启动整个潮流计算过程。

运行此代码，可以得到如下结果：

为了方便后面的对比分析，我们着重分析平衡节点功率，线路功率，功率损耗这三个结果：

母线1 (平衡节点) 的功率：

有功功率（P）：70 MW

无功功率（Q）：20 MVAr

线路功率：

线路1-4：

有功功率（P）：70 MW

无功功率（Q）：20 MVAr

线路4-5：

有功功率（P）：30 MW

无功功率（Q）：10 MVAr

线路5-6：

有功功率（P）：-60 MW（负值表示方向相反）

无功功率（Q）：-13 MVAr

功率损耗：

总有功损耗：4.64 MW

总无功损耗：48.38 MVAr

但是在实际的代码编写过程中，我在雅可比矩阵的计算中曾出现过错误，导致我的计算结果出现了不收敛的情况，之所以我能确定是雅可比矩阵的问题，是因为我在改变初值之后，仍然未收敛，修改get_power_mismatch函数之后尽管结果表现为收敛，但是出现了NaN的情况，这也就意味着我的雅可比矩阵是一个奇异矩阵，后来经过修改与验证才得到了前面所列举出来的代码，为了避免再次出现这种错误，我将我的错误代码放在下面：

function J = calcJacobian(Ybus, V)
    % 获取导纳矩阵的维度
    n = length(Ybus);

    % 分别计算G和B矩阵
    [G, B] = dSbus_dV(Ybus, V);

    % 初始化雅可比矩阵
    J = zeros(2*n, 2*n);

    % 计算雅可比矩阵的四个分块
    J11 = real(G); % dP/dθ
    J12 = imag(G); % dP/dV
    J21 = real(B); % dQ/dθ
    J22 = imag(B); % dQ/dV

    % 组合四个分块
    J(1:n, 1:n) = J11;
    J(1:n, n+1:2*n) = J12;
    J(n+1:2*n, 1:n) = J21;
    J(n+1:2*n, n+1:2*n) = J22;
end

function [G, B] = dSbus_dV(Ybus, V)
    n = length(V);
    I = Ybus * V;
    Ybus_diag = diag(Ybus);
    diagV = diag(V);
    diagI_conj = diag(conj(I));
    diagVnorm = diag(V ./ abs(V));
    G = diagV * conj(Ybus_diag) * diagVnorm - diagV * ...
        (diagI_conj * diagVnorm + diagVnorm * conj(diagI_conj));
    B = diagV * conj(Ybus_diag) * diagVnorm - diagV * ...
        (diagI_conj * diagVnorm - diagVnorm * conj(diagI_conj));
end

七、两种解算方式的结果对比分析

首先查看两组输出的解算数据的差异在哪：

1.的母线电压对比：

在单纯使用牛顿拉夫逊法的结果中，母线5、6、7、8、9的电压幅值分别为1.0200、1.0200、1.0180、1.0220和1.0170 p.u.。

而在使用matpower的结果中，这些母线的电压幅值略有不同：母线5为1.013 p.u.，母线6为1.032 p.u.，母线7为1.016 p.u.，母线8为1.026 p.u.，母线9为0.996 p.u.

2.母线功率 (P, Q) 对比：

使用matpower的解算输出中显示所有发电机母线的发电量和负荷母线的负荷。例如，母线1的有功和无功功率分别为71.64 MW和27.05 MVAr。

在单纯使用牛顿拉夫逊法的结果中，相同母线的有功和无功功率略有不同。例如，母线1的有功和无功功率分别为70 MW和20 MVAr。

3.线路功率和损耗对比：

在在单纯使用牛顿拉夫逊法的结果中，线路1-4、4-5和5-6的有功功率损耗分别为0、0和1 MW。

而在使用matpower的结果中，这些线路的有功损耗分别为0.000、0.166和1.354 MW。

4.总线损对比：

使用matpower解算的结果的总有功损耗为4.64 MW，无功损耗为48.38 MVAr。

在单纯使用牛顿拉夫逊法的结果中的总有功损耗和无功损耗分别为4.640000 MW和48.380000 MVAr，数值上略高于使用matpower的结果。

接下来便可以探究两种计算方式为什么会有不同的解算结果，要想探究这个问题，就必须从matpower的源代码中寻找它用于使用牛顿拉夫逊法的那部分代码，总结其算法，然后与我单纯使用牛顿拉夫逊法解算的算法差异。首先，通过对从官网上下载的整个matpower8.0b的安装包中的代码进行搜寻与查找，我找到了几个和牛顿拉夫逊法相关的代码文件，它们分别是runpf.m，newtonpf.m，makeY.m，dSbus_dV.m，calc_branch_angle.m，calc_v_i_sum.m，calc_v_y_sum.m，calc_v_pq_sum.m函数文件，通过对这些代码中的算法的研究，我总结了以下几点可能造成这种差异的原因：

1. 雅各比矩阵的构建：

MATPOWER实现：在MATPOWER中，dSbus_dV函数能够根据极坐标或笛卡尔坐标系计算雅各比矩阵。这意味着它能处理更复杂的电压表示形式，尤其是在处理不平衡负载或网络的非线性特性时。

手动实现：我只考虑了极坐标形式，而并没有特别处理笛卡尔坐标系，这可能导致在处理case9的网络问题时，计算效率或精度上有所不同。

2. 电压更新机制：

MATPOWER实现：newtonpf中的电压更新机制包括对PV母线的Q限制的处理。在迭代过程中，如果发现PV母线的无功功率超出其限制，MATPOWER会将其转换为PQ母线，并重新进行潮流计算。

手动实现：我在手动实现没有这种自动转换机制，并没有处理发电机Q限制。

3. 导纳矩阵构建：

MATPOWER实现：makeYbus函数在构建Ybus矩阵时考虑了网络中所有元件的导纳，包括线路的串联和并联导纳，变压器的抽头比等。

手动实现：我对于特定元件（如变压器抽头）并未单独设立算法去处理，这就会导会在网络的导纳表示上存在差异，从而影响最终的潮流计算结果。

4. 线路角度差计算：

MATPOWER实现：这个函数计算每条输电线路两端的电压相角差，对于监测系统稳定性和线路的功率流向至关重要。

手动实现：我的实现在这部分未进行相似的计算，无法提供关于线路电压相角差的详细信息。

5.其他可能原因：

MATPOWER在内部采用了优化的数据结构和计算方法，可能在计算效率上优于一般的手动实现。因此，在对case9或类似系统进行潮流计算时，MATPOWER的结果可能更为可靠，尤其是在处理复杂或非标准情况时。

八、对于 matpower 的深入研究与更优算法的探讨

首先是更深入地去研究matpower内置的函数，其中我了解到了最优潮流OPF算法，下面是我对这个算法的研究：

这个算法它旨在优化电力系统的运行和控制。OPF涉及到计算电力系统在满足所有操作约束的情况下的最优运行状态，包括发电成本最低、系统损耗最小化、电压稳定性和安全性等方面。

1.OPF需要建立数学模型来求解，该模型包括：

目标函数：通常是最小化发电成本，也可以是最小化损耗、改善电压剖面或其他目标。

约束条件：包括功率平衡约束（发电功率和负载功率之间的平衡）、电压限制、发电机出力限制、线路流量限制等。

2.OPF的解算方法针对不同的情境也是不同的：

线性规划：适用于目标函数和约束都是线性的情况。

非线性规划：当目标函数或约束条件非线性时使用。

内点法：有效处理大规模的非线性OPF问题。

遗传算法和粒子群优化：这些是启发式算法，对于非凸问题特别有效。

3.OPF的数学模型可以概括如下：

目标函数中包括系统的控制变量y和状态变量x。g(x,y)是等式约束，h(x,y)是不等式约束，分为变量不等式和函数不等式。

4.OPF在MATPOWER中的应用

MATPOWER作为一款广泛使用的电力系统分析工具，提供了求解OPF问题的功能。MATPOWER利用了MATLAB的优化工具箱，可以解决各种不同类型的OPF问题，包括经典的成本最小化、损耗最小化，以及更复杂的安全约束OPF问题。

5.MATPOWER解OPF的代码

% 载入案例
define_constants;  % 加载常量定义
mpc = case9;  % 这里使用MATPOWER内置的9节点测试系统

% 设置OPF选项
mpopt = mpoption('opf.ac.solver', 'MIPS');  % 设置OPF求解器为MIPS

% 运行OPF
results = runopf(mpc, mpopt);

% 显示结果
printpf(results);

运行上述代码可以得到解，但是由于解的篇幅太长，此处便不赘述展示。

然后我还展开了对于matpower解算潮流的更优秀的算法研究，通过对参考文献和github的各种搜寻，我找到了总计5种优化算法，它们分别是：连续潮流计算、概率潮流计算、结合人工智能和机器学习的潮流计算、分布式和并行计算以及自适应步长和混合方法，我将在文末列出这些算法的参考文献引用以及github开源地址。

针对这5种优化算法，由于连续潮流计算我找到了开源代码我便没有进行自行编写，结合人工智能和机器学习需要我有提前准备好的数据集，同时对于GPU等计算机硬件的算力要求较高，我的硬件环境还达不到这种算力，分布式和并行计算需要多个CPU同时工作，我只有一块CPU，自适应步长算法我没能研究透彻，因此我只针对概率潮流计算进行了自行编写代码的实践，同样是解matpower中的case9，其代码如下：

% 载入案例
define_constants;
mpc = case9;

% 设置变量的不确定性
% 假设负荷在其原始值上下波动10%
num_samples = 100;  % 定义样本数量
load_variation = 0.1;  % 负荷波动10%
load_buses = [5, 7, 9];  % 定义有波动负荷的节点

% 预先分配空间用于保存结果
results = cell(num_samples, 1);

% 设置OPF选项
mpopt = mpoption('opf.ac.solver', 'MIPS');

% 进行多次潮流计算
for i = 1:num_samples
    % 随机波动负荷
    for bus = load_buses
        mpc.bus(bus, PD) = mpc.bus(bus, PD) * ...
            (1 + load_variation * (2*rand - 1));
    end

    % 运行潮流计算
    results{i} = runpf(mpc, mpopt);
end

代码的逻辑与参数设计都被我写成注释了，可读性较高，由于代码运行迭代了多次，下面我将只展示最终的迭代结果：