jmatrix.m

返回电力系统规划中用matpower计算一次潮流后得到的雅可比矩阵.

function J=Jmatrix(branch0,bus0, P0, Q0)

% function [Jaco P Q]=FJaco(Yn,V,a)

%Heis为Jaco'对x的偏导
%H=dP/da;M=dP/dV;N=dQ/da;L=dQ/dV;


branch=branch0;
bus=bus0;
Pin=1000*P0;
Qin=1000*Q0;

% baseMVA=100;
% baseZ=(bus(1,10)^2)/baseMVA;

N=length(bus(:,1));

%%%%%%%%%%%%%   该程序的编号为1和N的节点互换

s1=find(branch(:,1)==1);
branch(s1,1)=N;
s2=find(branch(:,2)==N);
branch(s2,2)=1;

bc=bus(1,:);
bus(1,:)=bus(N,:);
bus(N,:)=bc;

bp=Pin(1);
Pin(1)=Pin(N);
Pin(N)=bp;

bq=Qin(1);
Qin(1)=Qin(N);
Qin(N)=bq;

%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%% 节点导纳矩阵

i=sqrt(-1);
Yn=zeros(N,N);
for k=1:length(branch(:,1))
      Yn( branch(k,1), branch(k,2) )=-1/(branch(k,3) + i*branch(k,4) ); 
      Yn( branch(k,2), branch(k,1) )=Yn( branch(k,1), branch(k,2) );
end

for k=1:N
    Yn(k,k)=-sum( Yn(k,:) );
end

% baseMVA=100;
% baseZ=(bus(1,10)^2)/baseMVA;
% Yn=Yn/baseZ;
%%%%%%%%%%%%%%%%

J11=zeros(N-1, N-1); J12=zeros(N-1, N-1);  
J21=zeros(N-1, N-1);  J22=zeros(N-1, N-1);

bus(:,8)=bus(:,8)*bus(1,10);

for i=1:N-1
    for j=1:N-1
         
        if j==i
            J11(i,i)=imag(Yn(i,i)) *bus(i,8)^2 + Qin(i);
            J12(i,i)=-real(Yn(i,i)) *bus(i,8)^2 - Pin(i);
            J21(i,i)=real(Yn(i,i)) *bus(i,8)^2 - Pin(i);
            J22(i,i)=imag(Yn(i,i)) *bus(i,8)^2 - Qin(i);
        else
            J11(i,j)=-bus(i,8)*bus(j,8)*( real(Yn(i,j))*sin( bus(i,9) - bus(j,9)) - imag(Yn(i,j))*cos( bus(i,9) - bus(j,9))  );
            J12(i,j)=-bus(i,8)*bus(j,8)*( real(Yn(i,j))*cos( bus(i,9) - bus(j,9)) + imag(Yn(i,j))*sin( bus(i,9) - bus(j,9))  );
            J21(i,j)= bus(i,8)*bus(j,8)*( real(Yn(i,j))*cos( bus(i,9) - bus(j,9)) + imag(Yn(i,j))*sin( bus(i,9) - bus(j,9))  );
            J22(i,j)=-bus(i,8)*bus(j,8)*( real(Yn(i,j))*sin( bus(i,9) - bus(j,9)) - imag(Yn(i,j))*cos( bus(i,9) - bus(j,9))  );
        end
        
    end
end
% J11(N,N)=0;  J12(N,N)=0;
% J21(N,N)=0;  J22(N,N)=0; 


J=[J11, J12; J21, J22];     %%%%%    输出N-1阶雅可比矩阵