您的当前位置：首页正文

Matlab 潮流计算程序N节点

来源：六九路网

%====================================================================

%潮流计算MATLAB 粗略程序

%====================================================================

%creat a new_data

t=0;

s=0;

r=0;

w=0;

number=input('How many node are there=');

% Convert Pq to a new array

for 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; end;

%record the number of the PQ node

end;

%Convert pv to a new array

for 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; end;

end;

%Convert set_v to a new array

%record the number of the PV node

for 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;

%creat a new_data2

[x,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;

for ii=1:x

for jj=3:14

new_data2(ii,jj)=data2(ii,jj);

end;

%creat a Y

Y=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;

%creat B, G

for ii=1:number

for jj=1:number

G(ii,jj)= real(Y(ii,jj));

B(ii,jj)= imag(Y(ii,jj));

end;

%creat Initial_P Initial_Q Initial_V

for 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 correcting

end;

Initial_p_q_v=[set_P;set_Q;set_V];

disp(Initial_p_q_v);

%creat Initial_e,Initial_f

for ii=1:number-1

e(ii,1)=1;

f(ii,1)=0.0;%change f to test used to be 1.0

end;

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 CALULATION

for try_time=1:25

%Creat every node consume P Q and U

n=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 Jacobian

disp(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;

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;

%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;

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;

%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 result

substract_result=Initial_p_q_v-sum;

%disp('substract_result');

%disp(substract_result);

%calculate delta_f_e

delta_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;

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 angle

PI=3.141592

for 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;

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(1,number);

Blance_P(1,1)=real(sum4)*100;

Blance_P(1,2)=a(1,number);

disp('Q Of the Balance node= ');

disp(Blance_Q);

disp('P Of the Balance node= ')

disp(Blance_P);

disp('=================================BUS-NO');

%calculate the Q of the P-V node

Q_PV_node=zeros(number,2);

Y_conj=conj(Y);

for ii=(s+1):(s+r)

for jj=1:number

Q_PV_node(ii,1)=Q_PV_node(ii,1)+(e(ii,1)+f(ii,1)*i).*(Y_conj(ii,jj).*(e(jj,1)-f(jj,1)*i));

end;

for ii=(s+1):(s+r);

Q_PV_node(ii,1)=Q_PV_node(ii,1).*100+new_data1(ii,8)*i;

end;

disp('This program is from blog.sina.com.cn/breadwinner') ;

for ii=1:number

old_number=a(1,ii);

Q_PV_node_old(old_number,1)=Q_PV_node(ii,1);

end;

for ii=1:number

Q_PV_node_old(ii,1)=imag(Q_PV_node_old(ii,1));

end;

for ii=1:number

Q_PV_node_old(ii,2)=ii;

end;

disp('Q gen=');

disp(Q_PV_node_old);

因篇幅问题不能全部显示，请点此查看更多更全内容

查看全文