Matlab-潮流计算程序N节点.doc

%=%=%=%潮流计算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 t=t+1; for jj=1:14 new_data1(t,jj)=data(ii,jj); end; a(1,t)=ii; s=s+1; %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; for jj=1:14 new_data1(t,jj)=data(ii,jj); end; a(1,t)=ii; r=r+1; %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 new_data1(t,jj)=data(ii,jj); end; a(1,t)=ii; w=w+1; end;end;%creat a new_data2x,y=size(data2)for ii=1:x for jj=1:2 for mm=1:number if data2(ii,jj)=a(1,mm) new_data2(ii,jj)=mm; end; end; end;end;for ii=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=new_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_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); 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_data2(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)./(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(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; 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(ii,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,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 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); % 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(ii,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).*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); 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 every 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,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=);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 ii=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); PE(ii,jj)=-G(ii,jj).*e(ii,1)-B(ii,jj).*f(ii,1); else ss=0; qq=0; for num=1:(n+m+1) 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; 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);%TEST+1 end; end;end;来源：(/s/blog_4c4e78800100cmlq.html) - Matlab潮流计算程序（改进）_小周_新浪博客 copy=3.14159;disp(=copy=)for ii=1:n for jj=1:m+n if (ii=jj) QE(ii,jj)=B(ii,jj).*e(ii,1)-G(ii,jj).*f(ii,1);%TEST+1 QF(ii,jj)=G(ii,jj).*e(ii,1)+B(ii,jj).*f(ii,1);%TEST+1 else ss=0; qq=0; for num=1:(n+m+1) 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; 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 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) if (ii=jj) UE(ii,jj)=0; UF(ii,jj)=0; else ss=0; qq=0; for num=1:(n+m+1) 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; UF(ii,jj)=-2.*f(ii,1); 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;%disp(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_edelta_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(substract_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_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) Old_Angle(ii,1)=Angle(jj,1); Old_Angle(ii,2)=ii; end;end;end;disp(=Angle=BUS-NO.);disp(Angle=);disp(Old_Angle);disp(=Try-times=)disp(Try-times=)disp(try_time);%disp(p=);%disp(p);% for jj=1:number% if a(1,jj)=118% man=jj % end;%end;%disp(man=);%disp(man)sum4=0;for jj=1:number Y_conj(number,jj)=conj(Y(number,jj); sum4=sum4+Y_conj(number,jj).*(e(jj,1)-f(jj,1)*i);end;%sum4=sum4*e(number,1);disp(=Balance P Q=BUS-NO);%disp(sum4);Blance_Q(1,1)=imag(sum4)*100;Blance_Q(1,2)=a(