form_jac.m

来自「计算机潮流计算 应用的是牛拉法」· M 代码 · 共 55 行

function J=form_jac(bus,Y,n,m)
H=zeros(n-1,n-1);  %H是(n-1)*(n-1)阶矩阵    
for i=1:n-1
    for j=1:n-1
        if(i~=j)
            delta=bus(i,3)-bus(j,3);  %节点ij的角度差
            H(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*sin(delta)-imag(Y(i,j))*cos(delta));
            
        else
            H(i,i)=bus(i,2)^2*imag(Y(i,i));
            for k=1:n
                delta=bus(i,3)-bus(k,3); %节点ik的角度差
                H(i,i)=H(i,i)+bus(i,2)*bus(k,2)*(real(Y(i,k))*sin(delta)-imag(Y(i,k))*cos(delta));
            end
        end
    end
end


N=zeros(n-1,m);   %N是(n-1)*m阶矩阵 
for i=1:n-1
    for j=1:m
        if(i~=j)
            delta=bus(i,3)-bus(j,3);  %节点ij的角度差
            N(i,j)=-bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(delta)+imag(Y(i,j))*sin(delta));
        else
            N(i,i)=-bus(i,2)^2*real(Y(i,i)); 
            for k=1:n
                delta=bus(i,3)-bus(k,3); %节点ik的角度差
                N(i,i)=N(i,i)-bus(i,2)*bus(k,2)*(real(Y(i,k))*cos(delta)+imag(Y(i,k))*sin(delta));
            end
        end
    end
end
    
K=zeros(m,n-1);   %K是m*(n-1)阶矩阵 
for i=1:m
    for j=1:n-1
        if(i~=j)
            delta=bus(i,3)-bus(j,3);  %节点ij的角度差
            K(i,j)=bus(i,2)*bus(j,2)*(real(Y(i,j))*cos(delta)+imag(Y(i,j))*sin(delta));
        else
            K(i,i)=N(i,i)+2*bus(i,2)^2*real(Y(i,i));
        end
    end
end


L=H([1:m],[1:m]);   %L是m*m阶矩阵 
for i=1:m
    L(i,i)=-H(i,i)+2*bus(i,2)^2*imag(Y(i,i));
end


J=[H,N;K,L];