利用matlab基础语法实现牛顿法潮流计算

问题描述

现代电力系统分析课程留了用计算机进行潮流计算的作业，可以使用各种软件实现，matlab就包含在其中，按老师的意思应该是使用matlab中的某个app。
可是我头铁，app是什么，当然要自己写啊！
以此题为例

整体思路

根据Y阵，PQ节点的P、Q值，PV节点的P、V值，平衡节点的δ、V值，利用极坐标的牛顿法，求出每个节点的δ、V值和P、Q值，和ΔP、ΔQ值，根据ΔP、ΔQ值判断是否收敛，同时得出下一个P、Q值，反复迭代求解。

具体步骤

1. 输入Y阵Ymat，PQ、PV节点参数Ps、Qs、Ps、V，平衡节点参数，此为目标值。
2. 给出初始V、δ值，其中V节点、平衡节点为1中的值，在迭代过程中不再变化。
3. 根据δ、U值求其对应的P、Q值。
4. 根据P、Q、δ、V、Ymat求Jacobi矩阵J。
5. 将复数矩阵Ymat进行分离，实部组成G阵，虚部组成B阵。根据如下公式计算J阵。
   $$\begin{gathered} H_{ii}=\{\begin{array}{cc}-V_i{V}_j(G_{ij}sin{\delta }_{ij}-B_{ij}cos{\delta }_{ij})& i̸=j\\ V_i^2{B}_{ii}+Q_i& i=j\end{array} \\ {N}_{ij}=\{\begin{array}{cc}-V_i{V}_j(G_{ij}cos{\delta }_{ij}+B_{ij}sin{\delta }_{ij})& i̸=j\\ -V_i^2{G}_{ii}-P_i& i=j\end{array} \\ {K}_{ij}=\{\begin{array}{cc}{V}_i{V}_j(G_{ij}cos{\delta }_{ij}+B_{ij}sin{\delta }_{ij})& i̸=j\\ V_i^2{G}_{ii}-P_i& i=j\end{array} \\ {L}_{ij}=\{\begin{array}{cc}-V_i{V}_j(G_{ij}sin{\delta }_{ij}-B_{ij}cos{\delta }_{ij})& i̸=j\\ V_i^2{G}_{ii}-Q_i& i=j\end{array} \\ J=\left[\begin{array}{cc}H& N\\ K& L\end{array}\right], \end{gathered}$$
6. 根据Ps与P，Qs与Q的差求出修正量DP、DQ。
7. 根据J阵与DP、DQ，利用牛顿法公式$$\left[\begin{array}{c}\Delta P\\ \Delta Q\end{array}\right]=-J\left[\begin{array}{c}\Delta \delta \\ V_{D_2}^{-1}\Delta V\end{array}\right],其中:V_{D_2}=\left[\begin{array}{cccc}{V}_1& & & \\ & V_2& & \\ & & ...& \\ & & & V_m\end{array}\right]$$求Δδ、ΔV。
8. 根据Δδ、ΔV、δ、V求得下一个δ、V。
9. 重复3~8。

程序代码

%chaoliujtz.m

%季天泽
%华中科技大学电气与电子工程学院
%强电磁与新技术国家重点实验室
clear
clc

%   节点必须按照PQ、PV、平衡节点来排序
%   PQ节点有m个，共有n+1个节点，即PV节点有n-m个
%   每个PQ节点有1个关于Q的方程式，每个PV或PQ节点有1个关于P的方程式，平衡节点没有方程式
%   所以P有n行，Q有m行，U与δ有n+1行

Ymat=[
    1.042093-i*8.242876      -0.588235+i*2.352941     i*3.666667       -0.453858+i*1.891074
    -0.588235+i*2.352941     1.069005-i*4.727377      0               -0.480769+i*2.403846
    i*3.666667               0                       -i*3.333333      0
    -0.453858+i*1.891074     -0.480769+i*2.403846     0               0.934627-i*4.261590
    ];
z=max(size(Ymat));
%Ps、Qs代表实际的P、Q
Ps(1)=-0.3;
Ps(2)=-0.55;
Ps(3)=0.5;
Qs(1)=-0.18;
Qs(2)=-0.13;
n=size(Ps);
m=size(Qs);
n=n(2);
m=m(2);
%delta=[ones(1,10)*-1,0]
delta(z)=0;
%输入U时，注意PV节点和平衡节点直接输入的是目标值的标幺值，而PQ节点输入的仅是迭代初值
U=[1,1,1.1,1.05];
output=[];
error=[];
G=real(Ymat);
B=imag(Ymat);
for ii=1:500
    p1=(U*(G.*cos(delta-delta')+B.*sin(delta-delta')));
    q1=(U*(G.*sin(delta-delta')-B.*cos(delta-delta')));
    P=U(1:n).*p1(1:n);
    Q=U(1:m).*q1(1:m);
    clear p1 q1;
    DP=Ps-P;
    DQ=Qs-Q;
    J = Jacobi(P,Q,Ymat,U,delta);
    invJ=inv(J);
    %修正方程
    M=-invJ*[DP';DQ'];
    Ddelta=M(1:n)';
    DU=(diag(U(1:m))*M(n+1:m+n))';
    Ddelta(n+1)=0;
    DU(n+1)=0;
    delta=delta+Ddelta;
    U=U+DU;
    error=[error;[DP,DQ,DU]];
    output=[output;U,delta];
end
plot(1:500,output)

%Jacobi.m


function J = Jacobi(P,Q,Y,U,delta)
%   U,Y,S都是复矩阵，获取之后进行拆分,m为PQ节点的个数，n为节点总数，n-1-m为PV节点数
%   此处显示详细说明
%   节点按照PQ、PV、平衡节点来排序
%   PQ节点有m个，共有n+1个节点，即PV节点有n-m个
%   每个PQ节点有1个关于Q的方程式，每个PV或PQ节点有1个关于P的方程式，平衡节点没有方程式
%   所以P有n行，Q有m行
n=size(P);
m=size(Q);
n=n(2);
m=m(2);
G=real(Y);
B=imag(Y);
Q(n+1)=0;
P(n+1)=0;
for i=1:n
    for j=1:n
        if i~=j
            H(i,j)=-U(i)*U(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
        elseif i==j
            H(i,j)=U(i)*U(j)*B(i,j)+Q(i);
        end
    end
end
for i=1:n
    for j=1:m
        if i~=j 
            N(i,j)=-U(i)*U(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
        elseif i==j
            N(i,j)=-U(i)*U(j)*G(i,j)-P(i);
        end
    end
end
for i=1:m
    for j=1:n
        if i~=j 
            K(i,j)=U(i)*U(j)*(G(i,j)*cos(delta(i)-delta(j))+B(i,j)*sin(delta(i)-delta(j)));
        elseif i==j
            K(i,j)=U(i)*U(j)*G(i,j)-P(i);
        end
    end
end
for i=1:m
    for j=1:m
        if i~=j 
            L(i,j)=-U(i)*U(j)*(G(i,j)*sin(delta(i)-delta(j))-B(i,j)*cos(delta(i)-delta(j)));
        elseif i==j
            L(i,j)=U(i)*U(j)*B(i,j)-Q(i);
        end
    end
end 
J=[
    H, N;
    K, L
   ];