谁知道风电场并网的潮流计算MATLAB程序

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; %record the number of the PQ node
end;
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; %record the number of the PV node
end;
end;
%Convert set_v to a new array
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;
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;
end;
end;
end;for ii=1:x
for jj=3:14
new_data2(ii,jj)=data2(ii,jj);
end;
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)/2i;
yy(jjj,iii)=new_data2(ii,8)/2i;
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_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)=00;%change f to test used to be 10
end; e(number,1)=new_data1(number,12);
f(number,1)=0;
% e(64,1)=088;%test 118ieee
% f(64,1)=039395826829394;
% f(14,1)=0;
% e(10,1)=1045;
%e(11,1)=101;
%e(12,1)=107;
%e(13,1)=109;
%//////////////////////////////////////////////////////////////////////////
%/////////////////////////////////////////////////////////////////////////
%//////////////////////////////////////////////////////////////////////////
%//////////////////////////////////////////////////////////////////////////
% 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;
end;
end;copy=314159;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)=-2f(ii,1);
UE(ii,jj)=-2e(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 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;
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 angle
PI=3141592
for ii=1:number
Angle(ii,1)=atan(f(ii,1)/e(ii,1))/PI180;
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=sum4e(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 nodeQ_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;
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 blogsinacomcn/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); 1、矩阵的稀疏性没有体现，delta_f_e=-inv(Jacobian)substract_result;最好用一个函数来求稀疏线性方程组的解，在该函数中可以考虑Jacobian阵的结构特点（用all函数），然后进行合适的节点编号，使得注入元达到最小，而后进行前代和回代过程；
2、可以考虑用极坐标来形成Jacobian矩阵，这样的话，矩阵在结构上是位子对称的，很容易用矩阵分块技术来形成结构对称的分块对角矩阵，而后用递归或者是并行求解该Jacobian方程，这样的方法比较适用于大规模的电力系统仿真计算中。

以上就是关于谁知道风电场并网的潮流计算MATLAB程序全部的内容，包括:谁知道风电场并网的潮流计算MATLAB程序、、等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！