Matlab 潮流计算程序N节点

1、%=%=%=%潮流计算MATLAB 粗略程序              %=%=%=%creat a new_datat=0;s=0;r=0;w=0;number=input('How many node are there=');% Convert Pq to a new arrayfor ii=1:number    if data(ii,4)=1     

2、   t=t+1;            for jj=1:14                new_data1(t,jj)=data(ii,jj);            end;  &

3、#160;         a(1,t)=ii;            s=s+1;                           &

4、#160;            %record the number of the PQ node    end;end;%Convert pv to a new arrayfor ii=1:number    if data(ii,4)=2        t=t+1;       &#

5、160;    for jj=1:14                new_data1(t,jj)=data(ii,jj);            end;            a(1,t)=ii

6、;            r=r+1;                                     

7、0;  %record the number of the PV node    end;end;%Convert set_v to a new arrayfor ii=1:number    if data(ii,4)=3        t=t+1;            for jj=1:14    &#

8、160;           new_data1(t,jj)=data(ii,jj);            end;            a(1,t)=ii;          

9、;  w=w+1;    end;end;%creat a new_data2x,y=size(data2)             for ii=1:x    for jj=1:2        for mm=1:number        

10、60;   if data2(ii,jj)=a(1,mm)                new_data2(ii,jj)=mm;            end;        end;    end;end;for i

11、i=1:x    for jj=3:14        new_data2(ii,jj)=data2(ii,jj);    end;end;%creat a YY=zeros(number,number);YY=zeros(number,number);yy=zeros(number,number);for ii=1:x   % for jj=1:14        iii=n

12、ew_data2(ii,1);        jjj=new_data2(ii,2);        if new_data2(ii,5)=2            sub=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)-new_

13、data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*i;            Y(iii,jjj)=-sub./new_data2(ii,14);            YY(iii,jjj)=sub./new_data2(ii,14);  

14、0;         Y(jjj,iii)=-sub/new_data2(ii,14);             YY(jjj,iii)=sub./new_data2(ii,14);            yy(iii,jjj)=(1.-new_data2(ii,14)./(new_dat

15、a2(ii,14).*new_data2(ii,14).*sub;            yy(jjj,iii)=(new_data2(ii,14)-1)./(new_data2(ii,14).*sub;        else            Y(iii,jjj)=-new_data2(ii,6)./

16、(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)+new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*i;             YY(iii,jjj)=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(i

17、i,6).*new_data2(ii,6)-new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*i;            Y(jjj,iii)=-new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)+new_data2(ii,7)./(new_data2(i

18、i,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)*i;            YY(jjj,iii)=new_data2(ii,6)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii,6)-new_data2(ii,7)./(new_data2(ii,7).*new_data2(ii,7)+new_data2(ii,6).*new_data2(ii

19、,6)*i;            yy(iii,jjj)=new_data2(ii,8)./2.*i;            yy(jjj,iii)=new_data2(ii,8)./2.*i;        end;    %end;end;for iii=1:number&

20、#160;   Y(iii,iii)=0;end;%for ii=1:x  %  for jj=1:14    for     iii=1:number        for jj=1:number   % if iii=jj    Y(iii,iii)=Y(iii,iii)+YY(iii,jj)+yy(iii,jj);    

21、  % end;   end;end;%creat B, Gfor ii=1:number        for jj=1:number            G(ii,jj)= real(Y(ii,jj);             B(ii,jj)= imag(Y(i

22、i,jj);    end;end;%creat Initial_P Initial_Q Initial_Vfor ii=1:(s+r)    set_P(ii,1)=(new_data1(ii,9)-new_data1(ii,7)./100;end;for ii=1:s;    set_Q(ii,1)=(new_data1(ii,10)-new_data1(ii,8)./100;end;for ii=1:r    set_V(ii,1)=new_data1(ii+s,12)

23、.*new_data1(ii+s,12);%try to modify for sike of correctingend;Initial_p_q_v=set_P;set_Q;set_V;disp(Initial_p_q_v);%creat Initial_e,Initial_ffor ii=1:number-1    e(ii,1)=1;    f(ii,1)=0.0;%change f to test used to be 1.0end;  e(number,1)=new_data1(number,12); 

24、f(number,1)=0; % e(64,1)=0.88;%test 118ieee % f(64,1)=0.39395826829394;   % f(14,1)=0;   % e(10,1)=1.045;   %e(11,1)=1.01;    %e(12,1)=1.07;    %e(13,1)=1.09;%/%/%/%/% Start NEWTOWN CALULATIONfor try_time=1:25  %Creat ev

25、ery node consume P Q and Un=s;m=r;for ii=1:(n+m)    sum1=0;    for jj=1:(n+m+1)        sum1=sum1+e(ii,1).*(G(ii,jj).*e(jj,1)-B(ii,jj).*f(jj,1)+f(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*e(jj,1);    end;    p(ii

26、,1)=sum1;end;for ii=1:n    sum2=0;    for jj=1:(n+m+1)        sum2=sum2+f(ii,1).*(G(ii,jj).*e(jj,1)-B(ii,jj).*f(jj,1)-e(ii,1).*(G(ii,jj).*f(jj,1)+B(ii,jj).*e(jj,1);    end;    q(ii,1)=sum2;end;disp('q=

27、');disp(q);u=zeros(n+m),1);for ii=(n+1):(n+m)    u(ii,1)=e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);end;for ii=n+1:(n+m)     extra_u(ii-n),1)=u(ii,1);end;disp('extra_u=');disp(extra_u);sum=p;q;extra_u;disp(sum)disp(s);disp(p);%creat Jacobiandisp(n);disp(m);for i

28、i=1:(n+m)    for jj=1:(n+m)        if (ii=jj)            PF(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);               &#

29、160;      PE(ii,jj)=-G(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1);                             else       

30、60;    ss=0;            qq=0;                       for num=1:(n+m+1)        

31、60;       ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);                qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);            end; 

32、                      PF(ii,jj)=-ss+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1            PE(ii,jj)=-qq-G(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1)

33、;%TEST+1                   end;    end;end;来源：( - Matlab 潮流计算程序（改进）_小周_新浪博客 copy=3.14159;disp('=copy=')for ii=1:n    for jj=1:m+n      &#

34、160; if (ii=jj)                     QE(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1                  

35、60;    QF(ii,jj)=G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1                   else            ss=0;      

36、0;     qq=0;                      for num=1:(n+m+1)                ss=ss+G(ii,num).*f(num,1)+B(ii,nu

37、m).*e(num,1);                qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);            end;             &

38、#160;                    QF(ii,jj)=-qq+G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1            QE(ii,jj)=ss+B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1

39、                  end;    end;end;%disp('QF');%disp(QF);%disp('QE');%disp(QE);UE=zeros(n+m),(n+m);UF=zeros(n+m),(n+m);for ii=n+1:n+m    for jj=1:(n+m)   &

40、#160;    if (ii=jj)                      UE(ii,jj)=0;            UF(ii,jj)=0;        else&

41、#160;           ss=0;            qq=0;                       for num=1:(n+m+1) &

42、#160;              ss=ss+G(ii,num).*f(num,1)+B(ii,num).*e(num,1);                qq=qq+G(ii,num).*e(num,1)-B(ii,num).*f(num,1);     

43、60;      end;                                  UF(ii,jj)=-2.*f(ii,1);       

44、;     UE(ii,jj)=-2.*e(ii,1);        end;    end;end;for ii=(n+1):(n+m)    for jj=1:(n+m)    extra_UE(ii-n),jj)=UE(ii,jj);    extra_UF(ii-n),jj)=UF(ii,jj);    end;end;%dis

45、p('extra_UE');%disp(extra_UE);%disp('extra_Uf');%disp(extra_UF);Jacobian=PF,PE;QF,QE;extra_UF,extra_UE;%disp('Jacobian=');%disp(Jacobian);%creat substract resultsubstract_result=Initial_p_q_v-sum;%disp('substract_result');%disp(substract_result);%calculate delta_f_ede

46、lta_f_e=-inv(Jacobian)*substract_result;%disp(delta_f_e);for ii=1:number-1;    f(ii,1)=f(ii,1)+delta_f_e(ii,1);    e(ii,1)=e(ii,1)+delta_f_e(ii+number-1,1);end;if max(substract_result)<1e-4    break;end ;end;%disp('substract_result');%disp(subs

47、tract_result);%disp('e=');%disp(e);%disp('f=');%disp(f);for ii=1:number   uuu(ii,1)= e(ii,1).*e(ii,1)+f(ii,1).*f(ii,1);U_RESULT(ii,1)=sqrt(uuu(ii,1);end;for ii=1:number  for  jj=1:number    if ii=a(1,jj)        Old

48、_Uresult(ii,1)=U_RESULT(jj,1)    end;end;end;for ii=1:number           Old_Uresult(ii,2)=ii;end;%disp('U_result');%disp(U_RESULT);disp('=');disp('The last result is :')disp('=U=BUS-NO.');disp('U=')disp(Old_Uresult);%calculate the anglePI=3.141592for ii=1:number       Angle(ii,1)=atan(f(ii,1)./e(ii,1)./PI*180;   end;for ii=1:number  for  jj=1:number    if ii=a(1,jj)