牛顿拉夫逊法 matlab 程序,牛顿拉夫逊法程序设计

该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

function gsj

clear all;

n=0;

fidin=fopen('start3.txt');

while ~feof(fidin) % 判断是否为文件末尾

tline=fgetl(fidin); % 从文件读行

if double(tline(1))>=48&&double(tline(1))<=57||double(tline(1))==43||double(tline(1))==45

% 判断首字符是否是数值 +号对应43，-号对应45

n=n+1;

if n==1

PQ=tline;

PQ=sscanf(PQ,'%f') %PQ数组中的偶数元素为Q，奇数元素为P

n1=size(PQ)/2;

numofPQ=n1(1) %统计PQ节点个数

end

if n==2

PV=tline; %PV节点为节点2，总之PV节点的编号在PQ节点编号前面

PV=sscanf(PV,'%f') %PV数组中偶数元素为电压模值U，奇数元素为P

n2=size(PV)/2;

numofPV=n2(1) %统计PV节点个数

end

if n==3

QVmaxmin=tline; %PV节点为节点2，总之PV节点的编号在PQ节点编号前面

QVmaxmin=sscanf(QVmaxmin,'%f') %PV数组中偶数元素为电压模值U，奇数元素为P

end

if n==4

S=tline; %平衡节点为节点1

S=sscanf(S,'%f')

end

if n==5

Yb0=tline; %自导纳

Yb0=sscanf(Yb0,'%f')

end

if n>=6

Yb=tline; %

Yb=sscanf(Yb,'%f');

num=numofPQ+numofPV+1 ; %num为系统节点数

for j=1:(numofPQ+numofPV+1) %1对应的平衡节点

YBre(n-5,j)=Yb(2*j-1); %节点导纳阵实部

YBim(n-5,j)=Yb(2*j); %虚部

if (n-5)==j

YBre(n-5,j)=YBre(n-5,j)+Yb0(2*j-1);

YBim(n-5,j)=YBim(n-5,j)+Yb0(2*j);

end

end

end

end

end

YBB=YBre+i*YBim %YBB为节点导纳矩阵

fclose(fidin);

%%%%%%%%%%%%%%%%上面部分为数据读入及导纳矩阵形成，下面为定义部分%%%%%%%%%%%%%%%%%%%%%%%%%

e=ones(num-1,1); %对除平衡节点外的节点电压e和f赋初值，从标号1到标号num-1的数组内放置的

f=zeros(num-1,1); %是从第2节点到第num节点的元素

P=zeros(num-1,1); %注释同上

Q=zeros(num-1,1);

Jacobi=zeros(2*(num-1));

I=zeros(1,num-1); %存放注入电流

detaPQU=zeros(2*(num-1),1); %存放P,Q,V的不平衡量

detaFE= zeros(2*(num-1),1); %存放f,e的不平衡量,注意：f在前，e在后

number=1; %存放迭代次数

flagPVPQ=zeros(num-1,1); %标记PV是否转换为PQ

%%%%%%%%%%%%%%%%%%%%%%%%%%％％％％％％下面为潮流计算部分％％％％％％%%%%%%%%%%%%%%%%%%%%%%%%%

e1=S(1);f1=S(2); %1号节点为平衡节点

flowComputer=fopen('flowComputer.txt','wt'); % 创建flowComputer.txt文件

while number<=100

number=number %显示迭代次数

fprintf(flowComputer,'迭代次数：%g\n',number);

fprintf(flowComputer,'P,Q,U 不平衡量：\n');

for j=2:num

sumP=e(j-1)*(YBre(j,1)*e1-YBim(j,1)*f1)+f(j-1)*(YBre(j,1)*f1+YBim(j,1)*e1);

%if j>numofPV+1 %判断是否为PQ节点

sumQ=f(j-1)*(YBre(j,1)*e1-YBim(j,1)*f1)-e(j-1)*(YBre(j,1)*f1+YBim(j,1)*e1);

%end

n=2;

while n<=num %带入电压迭代值求P、Q的迭代值

sumP=sumP+ e(j-1)*(YBre(j,n)*e(n-1)-YBim(j,n)*f(n-1))+ f(j-1)*(YBre(j,n)*f(n-1)+YBim(j,n)*e(n-1));

% if j>numofPV+1 %判断是否为PQ节点

sumQ=sumQ+ f(j-1)*(YBre(j,n)*e(n-1)-YBim(j,n)*f(n-1))- e(j-1)*(YBre(j,n)*f(n-1)+YBim(j,n)*e(n-1));

%无论PQ节点还是PV节点，其有功和无功在后面的节点注入电流计算中均会用到

% end

n=n+1;

end

if j<=(numofPV+1)

P(j-1)=sumP;

Q(j-1)=sumQ;

detaPQU(2*(j-1)-1)=PV(2*(j-1)-1)-sumP; %PV节点的P和U的不平衡量

if sumQ>=QVmaxmin(2*(j-1))&&sumQ<=QVmaxmin(2*(j-1)-1) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

detaPQU(2*(j-1)) =PV(2*(j-1))^2-e(j-1)^2-f(j-1)^2;

elseif sumQ

detaPQU(2*(j-1))=QVmaxmin(2*(j-1))-sumQ; %QVmaxmin(2*(j-1)-1)对应的是QVmin

flagPVPQ(1,j-1)=1;

else

detaPQU(2*(j-1))=QVmaxmin(2*(j-1)-1)-sumQ; %QVmaxmin(2*(j-1)-1)对应的是QVmax

flagPVPQ(1,j-1)=1;

end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fprintf(flowComputer,'%g\n',detaPQU(2*(j-1)-1));

fprintf(flowComputer,'%g\n',detaPQU(2*(j-1)));

else

P(j-1)=sumP;

Q(j-1)=sumQ;

detaPQU(2*(j-1)-1)=PQ(2*(j-1-numofPV)-1)-sumP; %PQ节点的P和Q的不平衡量

detaPQU(2*(j-1)) =PQ(2*(j-1-numofPV))-sumQ;

fprintf(flowComputer,'%g\n',detaPQU(2*(j-1)-1));

fprintf(flowComputer,'%g\n',detaPQU(2*(j-1)));

end

end

detaPQU=detaPQU

displayPQ=P+i*Q

%for j=2:num

%fprintf(flowComputer,'%g\n',detaPQU(2*(j-1)-1)+i*detaPQU(2*(j-1)));

%end

fprintf(flowComputer,'\n');

%%%%%%%%%%%以上为赋初值及求不平衡量部分,以下为Jacobi矩阵元素的求解%%%%%%%%%%%%%%

for j=2:num

I(j-1)=(P(j-1)-i*Q(j-1))/(e(j-1)-i*f(j-1));%除平衡节点外的节点都有注入电流,I(i)=a(i)+j*b(i)

end

for j=2:num

for k=2:num

if k==j

Jacobi(2*(j-1)-1,2*(j-1)-1)=-YBim(j,j)*e(j-1)+YBre(j,j)*f(j-1)+imag(I(j-1));%H(i,i)

Jacobi(2*(j-1)-1,2*(j-1)) = YBre(j,j)*e(j-1)+YBim(j,j)*f(j-1)+real(I(j-1));%N(i,i)

if j<=1+numofPV&&flagPVPQ(j-1,1)==0 %判断PV节点是否要转换为PQ节点

Jacobi(2*(j-1),2*(j-1)-1) = 2*f(j-1); %R(i,i)

Jacobi(2*(j-1),2*(j-1)) = 2*e(j-1); %S(i,i)

else %PQ

Jacobi(2*(j-1),2*(j-1)-1) =-YBre(j,j)*e(j-1)-YBim(j,j)*f(j-1)+real(I(j-1));%J(i,i)

Jacobi(2*(j-1),2*(j-1)) =-YBim(j,j)*e(j-1)+YBre(j,j)*f(j-1)-imag(I(j-1));%H(i,i)

end

else

if YBre(j,k)==0&&YBim(j,k)==0 %判断i，j两个节点之间是否有直接的联系

continue;

else

Jacobi(2*(j-1)-1,2*k-3) =-YBim(j,k)*e(j-1)+YBre(j,k)*f(j-1); %H(i,j)

Jacobi(2*(j-1)-1,2*k-2) = YBre(j,k)*e(j-1)+YBim(j,k)*f(j-1); %N(i,j)

if j<=1+numofPV %PV

Jacobi(2*(j-1),2*k-3) = 0; %R(i,j)

Jacobi(2*(j-1),2*k-2) = 0; %S(i,j)

else %PQ

Jacobi(2*(j-1),2*k-3) =-Jacobi(2*(j-1)-1,2*k-2); %J(i,j)

Jacobi(2*(j-1),2*k-2) = Jacobi(2*(j-1)-1,2*k-3); %H(i,j)

end

end

end

end

end

%Jacobi=Jacobi

invofJ=inv(Jacobi)

detaFE=invofJ*detaPQU

fprintf(flowComputer,'电压不平衡量(f在前，e在后)：\n');

for j=1:num-1

fprintf(flowComputer,'%g\n', detaFE(2*j-1));

fprintf(flowComputer,'%g\n', detaFE(2*j));

end

fprintf(flowComputer,'\n');

%absdetaFE=abs(detaFE)

%maximume=max(abs(detaFE))

if max(abs(detaFE))<=1e-3 %如果e和f的不平衡量达到指定精度，可停止迭代

display('This is an end');

break;

else

number=number+1;

fprintf(flowComputer,'经过不平衡量修正后的电压：\n');

for j=2:num

f(j-1)=f(j-1)+detaFE(2*(j-1)-1); %注意：f在前，e在后

e(j-1)=e(j-1)+detaFE(2*(j-1));

fprintf(flowComputer,'%g',e(j-1));

fprintf(flowComputer,'%gi\n',f(j-1));

end

fprintf(flowComputer,'\n');

u=[e+i*f ,abs(e+i*f)]

end

end %对应while循环

fclose(flowComputer);

%%%%%%%%%%%%%%%%%%%平衡节点功率及线路功率的计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

stable=(YBre(1,1)-i*YBim(1,1))*(S(1)-i*S(2));

for j=2:num

if YBre(1,j)==0&&YBim(1,j)==0

continue;

else

stable=stable+(YBre(1,j)-i*YBim(1,j))*(e(j-1)-i*f(j-1));

end

end

Ss=stable*(S(1)+i*S(2)) %平衡节点功率

powerOfline=zeros(num); %存放线路功率

for j=1:num

for k=1:num

if j~=k

if j==1

powerOfline(j,k)=(e1^2+f1^2)*(YBre(j,k)-i*YBim(j,k))-(e1+i*f1)*(e(k-1)-i*f(k-1))*(YBre(j,k)-i*YBim(j,k));

elseif k==1

powerOfline(j,k)=(e(j-1)^2+f(j-1)^2)*(YBre(j,k)-i*YBim(j,k))-(e(j-1)+i*f(j-1))*(e1-i*f1)*(YBre(j,k)-i*YBim(j,k));

else

powerOfline(j,k)=(e(j-1)^2+f(j-1)^2)*(YBre(j,k)-i*YBim(j,k))-(e(j-1)+i*f(j-1))*(e(k-1)-i*f(k-1))*(YBre(j,k)-i*YBim(j,k));

end

end

end

end

powerOfline= powerOfline

powerloss=0; %存放功率损耗

for j=1:num

for k=1:num

if j~=k

powerloss(j,k)=powerOfline(j,k)+powerOfline(k,j);

end

end

end

powerloss=powerloss