内点法求最优潮流（Matlab代码实现）-阿里云开发者社区

💥💥💞💞欢迎来到本博客❤️❤️💥💥

🏆博主优势：🌞🌞🌞博客内容尽量做到思维缜密，逻辑清晰，为了方便读者。

⛳️座右铭：行百里者，半于九十。

📋📋📋本文内容如下：🎁🎁🎁

⛳️赠与读者

👨‍💻做科研，涉及到一个深在的思想系统，需要科研者逻辑缜密，踏实认真，但是不能只是努力，很多时候借力比努力更重要，然后还要有仰望星空的创新点和启发点。建议读者按目录次序逐一浏览，免得骤然跌入幽暗的迷宫找不到来时的路，它不足为你揭示全部问题的答案，但若能解答你胸中升起的一朵朵疑云，也未尝不会酿成晚霞斑斓的别一番景致，万一它给你带来了一场精神世界的苦雨，那就借机洗刷一下原来存放在那儿的“躺平”上的尘埃吧。

或许，雨过云收，神驰的天地更清朗.......🔎🔎🔎

1 内点法

内点法属于约束优化算法。约束优化算法的基本思想是：通过引入效用函数的方法将约束优化问题转换成无约束问题，再利用优化迭代过程不断地更新效用函数，以使得算法收敛。

内点法（罚函数法的一种）的主要思想是：在可行域的边界筑起一道很高的“围墙”，当迭代点靠近边界时，目标函数徒然增大，以示惩罚，阻止迭代点穿越边界，这样就可以将最优解“档”在可行域之内了。

详细讲解：内点法

2 最优潮流概述

2.1 概述

最优潮流和基本潮流的比较潮流计算可以归结为针对一定的扰动变量p(负荷情况)，根据给定的控制变量u(如发电机的有功出力、无功出力或节点电压模值等)，求出相应的状态变量x(如节点电压模值及角度)，这样通过一次潮流计算得到的潮流解决定了电力系统的一个运行状态。这种潮流计算也可以称之为基本潮流(或常规潮流)计算，一次基本潮流计算的结果主要满足了潮流方程式或变量间等式约束条件：

编辑

一次潮流计算所决定的运行状态可能由于某些状态变量或者作为u,x 函数的其它变量在数值上超出了它们所容许的运行限值(即不满足不等式约束条件)，因而在技术上并不是可行的。工程实际上常用的方法是调整某些控制变量的给定值，重新进行前述的基本潮流计算，这样反复进行，直到所有的约束条件都能够得到满足为止。这样便得到了一个技术上可行的潮流解。

由于系统的状态变量及有关函数变量的上下限值间有一定的间距，控制变量也可以在其一定的容许范围内调节，因而对某一种负荷情况，理论上可以同时存在为数众多的、技术上都能满足要求的可行潮流解。每一个可行潮流解对应于系统的某一个特定的运行方式，具有相应总体的经济上或技术上的性能指标(如系统总的燃料消耗量、系统总的网损等)，为了优化系统的运行，就有需要从所有的可行潮流解中挑选出上述性能指标为最佳的一个方案。而这就是本节要讨论的最优潮流所要解决的问题。

因此所谓最优潮流，就是当系统的结构参数及负荷情况给定时，通过控制变量的优选，所找到的能满足所有指定的约束条件，并使系统的某一个性能指标或目标函数达到最优时的潮流分布。

2.2 最优潮流和基本潮流比较

(1)基本潮流计算时控制变量u是事先给定的；而最优潮流中的u则是可变而待优选的变量，为此在最优潮流模型中必然有一个作为u优选准则的目标函数。

(2)最优潮流计算除了满足潮流方程这一等式约束条件之外，还必须满足与运行限制有关的大量不等式约束条件。

(3)进行基本潮流计算是求解非线性代数方程组；而最优潮流计算由于其模型从数学上讲是一个非线性规划问题，因此需要采用最优化方法来求解。

(4)基本潮流计算所完成的仅仅是一种计算功能，即从给定的u求出相应的x；而最优潮流计算则能够根据特定目标函数并在满足相应约束条件的情况下，自动优选控制变量，这便具有指导系统进行优化调整的决策功能。

2.3 最优潮流计算的内点法

内点法则是一种在可行域内部寻优的方法。1984年，Karmarkar提出的内点法具有多项式计算复杂性，在求解大规模线性规划问题时，计算速度比单纯形法快50倍以上。 1986年，Gill将内点法推广到非线性规划领域。

（1）基本思想

内点法最初的基本思想是希望寻优迭代过程始终在可行域内进行，因此，初始点应取在可行域内，并在可行域的边界设置“障碍”使迭代点均为可行域的内点。

（2）存在的困难

但对实际大规模系统，初始可行点的寻找比较困难。

（3）改进方法－跟踪中心轨迹内点法

跟踪中心轨迹内点法对此作了改进，只要求在寻优过程中松弛变量和拉格朗日乘子满足简单的大于零或小于零的条件，可代替原来必须在可行域求解的要求，使计算过程大为简化。

（4）跟踪中心轨迹内点法详介

1）路径跟踪法的基本数学模型

内点法最初的基本思路是希望通过寻优迭代过程始终在可行域内进行，因此，初始点应在可行域内，并在可行域的边界设置“障碍”使迭代点接近边界时其目标函数迅速增大，从而保证迭代点均在可行域的内点。但是对于大规模实际问题而言，寻找初始点往往十分困难。以下介绍的路径跟踪法只要求在寻优过程中松弛变量和拉格朗日乘子满足简单的大于0或者小于0的条件，即可代替原来必须在可行域内求解的要求，使得计算过程大为简化。

一般可以将最优潮流模型简化为如下的非线性优化模型：

编辑

其中扰动因子或者障碍因子u>0。当l或u接近边界时,以上函数将趋于无穷大，因此满足以上障碍目标函数的极小解不可能在边界上找到。这样就通过目标函数的变化把含不等式限制的优化问题A变成只含等式限制优化的问题B了，因此可以直接用拉格朗日乘子法来求解。

编辑

2）路径跟踪法的最优潮流求解思路

路径跟踪法的最优潮流求解过程就是对拉格朗日目标函数求极小值问题

编辑

3 算例及程序实现流程

3.1 算例描述

以系统燃料最省为最优潮流的目标函数，求图2所示简化系统的系统燃料最省的最优潮流计算。线路传输功率边界、发电机有功无功出力上下界和燃料耗费曲线参数分别见表1、表2。所有数据都是以标幺值形式给出，功率基准值为100 MV·A，母线电压上下界分别为1.1和0.9。

3.2 程序具体实现过程

编辑

3.3 案例

编辑

4 Matlab实现

4.1 主函数

% 最优潮流（OPF）的内点法求解
% 说明：本程序包括一个主文件main.m和三个函数文件makeY.m，Coeff.m和dPQ.m，四个文件均应放在MATLAB当前目录下
% 运行main.m，会将计算结果打印到屏幕并保存至文件solution.txt中
clc;
clear;
close all
tic;
%% 算例数据
%节点数据表
%       1  2  3   4        5       6        7        8      9      10     11     12      13      14     15      16    17    18    19
%     节点 区 类 电压     相角    有功     无功     有功    无功   电压   期望    乏值    乏值    并联   并联   远端控 节点  电压  电压
%      号  号 型                  负荷     负荷     出力    出力   基准   电压    上限    下限    电导   电纳  制节点号 号   上限  下限
%  3代表平衡节点，2代表PV节点，0代表PQ节点
N=[ 1  1  0 1.0000   0.00     1.60     0.80     0.00    0.00  230.00 1.032   0.0000  0.0000  0.0000  0.0000    0    1    1.1   0.9
    2  1  0 1.0000   0.00     2.00     1.00     0.00    0.00   18.00 1.025   0.0000  0.0000  0.0000  0.0000    0    2    1.1   0.9
    3  1  0 1.0000   0.00     3.70     1.30     0.00    0.00   13.80 1.025   0.0000  0.0000  0.0000  0.0000    0    3    1.1   0.9
    4  1  2 1.0000   0.00     0.00     0.00     4.50    0.00  230.00 1.026   0.0000  0.0000  0.0000  0.0000    0    4    1.1   0.9
    5  1  3 1.0500   0.00     0.00     0.00     4.50    1.45    0.00 0.996   0.0000  0.0000  0.0000  0.0000    0    5    1.1   0.9];
%支路数据表
  %1  2  3  4  5  6   7      8      9     10    11    12     13     14    15      16      17   18    19   20   21   22
  %首 尾 区 区 电 类 电阻   电抗   电纳  支路额 支路下 支路上  控制   位 变压器最 移相器最 最小  最大  步长 最小 最大 支路
  %点 点 号 号 路 型                     定功率 限功率 限功率 母线号  置 终传输比 终移相角 变比  变比      电压 电压  号
B=[ 1  2  1  1  1  0  0.040  0.2500   0.50      0.    0.    2.      0    0    1.0       0.0    0.0  0.0  0.0  0.0 0.0    1
    1  3  1  1  1  0  0.100  0.3500   0.0       0.    0.    0.65    0    0    1.0       0.0    0.0  0.0  0.0  0.0 0.0    2
    2  3  1  1  1  0  0.080  0.3000   0.50      0.    0.    2.      0    0    1.0       0.0    0.0  0.0  0.0  0.0 0.0    3
    2  4  1  1  1  1  0.0    0.0150   0.0       0.    0.    6.      0    0    1.050     0.0    0.0  0.0  0.0  0.0 0.0    4
    3  5  1  1  1  1  0.0    0.0300   0.0       0.    0.    5.      0    0    1.050     0.0    0.0  0.0  0.0  0.0 0.0    5];
%发电机数据表
%             1    2     3    4    5   6       7         8         9
%            发电 节点 有功 无功 有功 无功   二次      一次       常数
%            机号  号  上界 上界 下界 下界   系数      系数        项
Gen=[ 1    4    8.0  3.0  1.0  -3.0  50.4395  200.4335  1200.6485
      2    5    8.0  5.0  1.0  -2.1  200.55   500.746   1857.201];
[num_node,~] = size(N);    %节点数
[num_branch,~] = size(B);  %支路数
% 不同类型节点个数
num_PQ = 0;
num_PV = 0;
for k = 1:num_node
    if (N(k,3) == 0)
        num_PQ = num_PQ+1;
    end
    if (N(k,3) == 2)
        num_PV = num_PV+1;
    end
end
num_gen = num_node-num_PQ; % 2 发电机个数
num_equa = 2*num_node; %10 等式约束个数m
num_inequa = 2*num_gen+num_node+num_branch; %14 不等式约束个数r
num_state = 2*num_gen+2*num_node; % 14 状态变量个数n
Pload = N(1:num_PQ,6); Qload = N(1:num_PQ,7);
a2 = Gen(:,7); a1 = Gen(:,8); a0 = Gen(:,9); %耗量特性多项式系数
A2 = diag(a2); A1 = diag(a1);              %创建对角矩阵或获取矩阵的对角元素
%% 生成节点导纳矩阵
Y = makeY(N, B, num_node, num_branch);
%% 参数设置
PG = N(num_PQ+1:num_node, 8);  %有功出力
QR = N(num_PQ+1:num_node, 9);  %无功出力
P_ejec = N(:,8)-N(:,6); Q_ejec = N(:,9)-N(:,7);  %有功出力-有功负荷，无功
Xtilde(1:2:num_node*2-1) = (N(:,5))'; Xtilde(2:2:num_node*2) = (N(:,4))';
u = ones(num_inequa,1); l = ones(num_inequa,1);
z = ones(1,num_inequa); w = -0.5*ones(1,num_inequa); y = 1E-10*ones(1,num_equa); y(2:2:num_equa) = -1*y(2:2:num_equa);
epsi = 1E-6;
sigma = 0.1;
num_iteration = 0;
max_iteration = 50;
g_u = [Gen(:,3); Gen(:,4); N(:,18); B(:,12)];
g_l = [Gen(:,5); Gen(:,6); N(:,19); -B(:,12)];
gap_record = zeros(max_iteration,1);
for num_iteration = 1:max_iteration
    gap = l'*z'-u'*w';
    gap_record(num_iteration) = gap;
    if (gap < epsi)
        disp('iterations completed!');
        break;
    end
    miu = gap*sigma/2/num_inequa;
    x = [PG; QR; Xtilde'];  % 状态变量x
    X = [z'; l; w'; u; x; y']; % 待修正的X
    g1 = PG; g2 = QR; % 有功出力g1，无功出力g2
    g3 = (Xtilde(2:2:num_node*2))'; % 节点电压幅值g3
    for k = 1:num_branch % 线路潮流g4
        ii = B(k,1); jj = B(k,2);
        theta = Xtilde(ii*2-1)-Xtilde(jj*2-1);
        g4(k,1) = Xtilde(2*ii)*Xtilde(2*jj)*(real(Y(ii,jj))*cos(theta)+imag(Y(ii,jj))*sin(theta))-Xtilde(2*ii)*Xtilde(2*ii)*real(Y(ii,jj));
    end
    g = [g1; g2; g3; g4];
    
    h = zeros(num_node*2,1);
    for ii = 1:num_node
        for jj = 1:num_node
            theta = Xtilde(ii*2-1)-Xtilde(jj*2-1);
            h(2*ii-1) = h(2*ii-1)-Xtilde(2*ii)*Xtilde(2*jj)*(real(Y(ii,jj))*cos(theta)+imag(Y(ii,jj))*sin(theta));
            h(2*ii) = h(2*ii)-Xtilde(2*ii)*Xtilde(2*jj)*(real(Y(ii,jj))*sin(theta)-imag(Y(ii,jj))*cos(theta));
        end
        h(2*ii-1) = h(2*ii-1)+P_ejec(ii);
        h(2*ii) = h(2*ii)+Q_ejec(ii);
    end
    
    
    L = diag(l); U = diag(u); Z = diag(z); W = diag(w); % Lagrange函数
    %%  Lagrange函数的偏导数
    Ly = h;
    Lz = g-l-g_l;
    Lw = g+u-g_u;
    L_miu_l = L*Z*ones(num_inequa,1)-miu*ones(num_inequa,1);
    L_miu_u = U*W*ones(num_inequa,1)+miu*ones(num_inequa,1);
    
    [A, dh_dx, dg_dx, H_, d2h_dx_y, d2g_dx_c, d2f_dx, temp, L_Z, U_W] = Coeff(num_node, num_PQ, Y, B, X, A2);
    %% 目标函数
    df_dx = [2*A2*PG+a1; zeros(2,1); zeros(length(Xtilde),1)];
    %% 拉格朗日函数
    Lx = df_dx-dh_dx*y'-dg_dx*(z+w)';
    Lx_ = Lx+dg_dx*(L\(L_miu_l+Z*Lz)+U\(L_miu_u-W*Lw));
    
    b = [-L\L_miu_l; Lz; -U\L_miu_u; -Lw; Lx_; -Ly];
    delta_X = dX(H_, dg_dx, dh_dx, L_Z, U_W, A, b, z, l', u', w, x', y);
    
    X = X + delta_X;
    
    %% 更新变量
    z = (X(1:num_inequa))';
    l = X(num_inequa+1:2*num_inequa);
    w = (X(2*num_inequa+1:3*num_inequa))';
    u = X(3*num_inequa+1:4*num_inequa);
    x = X(4*num_inequa+1:5*num_inequa);
    y = (X(5*num_inequa+1:5*num_inequa+num_equa))';
    
    PG = (x(1:num_gen));
    QR = (x(num_gen+1:2*num_gen));
    P_ejec(num_PQ+1:num_node) = PG;
    Q_ejec(num_PQ+1:num_node) = QR;
    Xtilde = (x(2*num_gen+1:2*num_gen+2*num_node))';
end
%% 相角变换到-π到π之间
for k = 1:num_node
    while Xtilde(2*k-1)>pi
        Xtilde(2*k-1) = Xtilde(2*k-1)-2*pi;
    end
    while Xtilde(2*k-1)<-pi
        Xtilde(2*k-1) = Xtilde(2*k-1)+2*pi;
    end
end
%% 以5号平衡节点相角为基准进行折算
for k = 1:num_node
    Xtilde(2*k-1) = Xtilde(2*k-1)-Xtilde(2*num_node-1);
end
if (num_iteration>=50)
    disp('OPF is not convergent!');
end
toc;
%% 保存数据
semilogy(gap_record(1:num_iteration));
grid on;
title('系统最优潮流内点法收敛特性');
xlabel('迭代次数'); ylabel('Gap');
source = [Gen(:,1:2) PG QR A2*(PG.*PG)+A1*PG+a0];
source_sum = sum(source(:,3:5));
voltage = [N(:,1) (Xtilde(2:2:2*num_node))' (Xtilde(1:2:2*num_node-1))'];
branch = [B(:,22) B(:,1:2) g4];
fileID = fopen('solution.txt','w+', 'n', 'UTF-8');
disp( '最优潮流计算结果：');
disp('====================================================');
disp( '|                  有功无功电源出力                |');
disp('====================================================');
disp( ' 发电  节点    有功    无功     燃料');
disp(' 机号   号     出力    出力     费用/$');
disp(['  ',num2str(source(1,1)),'     ',num2str(source(1,2)),'    ',num2str(source(1,3)),'   ',num2str(source(1,4)),'    ',num2str(source(1,5))]);
disp(['  ',num2str(source(2,1)),'     ',num2str(source(2,2)),'    ',num2str(source(2,3)),'   ',num2str(source(2,4)),'   ',num2str(source(2,5))]);
disp(['总计    —   ',num2str(source_sum(1)),'   ',num2str(source_sum(2)),'   ',num2str(source_sum(3))]);
% disp(['断 开 支 路:              ', num2str(o), '                  ',num2str(a)])
disp('====================================================');
disp( '|                     节点电压相量                 |');
disp('====================================================');
disp( '节点    电压      电压');
disp( '号      幅值      相角');
disp([num2str(voltage(1,1)),'       ',num2str(voltage(1,2)),'      ', num2str(voltage(1,3)) ]);
disp([num2str(voltage(2,1)),'       ',num2str(voltage(2,2)),'      ', num2str(voltage(2,3)) ]);
disp([num2str(voltage(3,1)),'       ',num2str(voltage(3,2)),'   ', num2str(voltage(3,3)) ]);
disp([num2str(voltage(4,1)),'       ',num2str(voltage(4,2)),'   ', num2str(voltage(4,3)) ]);
disp([num2str(voltage(5,1)),'       ',num2str(voltage(5,2)),'        ', num2str(voltage(5,3)) ]);
disp('====================================================');
disp('|                    支路有功功率                  |');
disp('====================================================');
disp('支路  从  到   功率');
disp('号 ');
fprintf(fileID, '\r\n%2d    %2d  %2d  %6.4f', branch');
type('solution.txt');

4.2 生成节点导纳矩阵

4.3 拉格朗日函数对所有变量及乘子的偏导数

4.4 算等式/不等式Jacobian矩阵、Hessian矩阵

5 结果

编辑

iterations completed!
时间已过 0.186363 秒。
最优潮流计算结果：
====================================================
|                  有功无功电源出力                |
====================================================
 发电  节点    有功    无功     燃料
 机号   号     出力    出力     费用/$
  1     4    5.5056   1.778    3833.0934
  2     5    2.1568   2.6194   3870.1109
总计    —   7.6624   4.3974   7703.2043
====================================================
|                     节点电压相量                 |
====================================================
节点    电压      电压
号      幅值      相角
1       0.9      -0.0069674
2       1.1      0.40491
3       1.0818   -0.057126
4       1.0697   0.47867
5       1.1        0
====================================================
|                    支路有功功率                  |
====================================================
支路  从  到   功率
号 
 1     1   2  -1.6064
 2     1   3  0.0064
 3     2   3  1.7709
 4     2   4  -5.5056
 5     3   5  -2.1568
>>

6 参考文献

文章中一些内容引自网络，会注明出处或引用为参考文献，难免有未尽之处，如有不妥，请随时联系删除。(文章内容资料获取，更多粉丝福利，MATLAB|Simulink|Python资源获取【请看主页然后私信】

内点法求最优潮流（Matlab代码实现）

⛳️赠与读者

1 内点法