該樓層疑似違規(guī)已被系統(tǒng)折疊?隱藏此樓查看此樓

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

% 判斷首字符是否是數(shù)值 +號(hào)對(duì)應(yīng)43，-號(hào)對(duì)應(yīng)45

n=n+1;

if n==1

PQ=tline;

PQ=sscanf(PQ,'%f') %PQ數(shù)組中的偶數(shù)元素為Q，奇數(shù)元素為P

n1=size(PQ)/2;

numofPQ=n1(1) %統(tǒng)計(jì)PQ節(jié)點(diǎn)個(gè)數(shù)

end

if n==2

PV=tline; %PV節(jié)點(diǎn)為節(jié)點(diǎn)2，總之PV節(jié)點(diǎn)的編號(hào)在PQ節(jié)點(diǎn)編號(hào)前面

PV=sscanf(PV,'%f') %PV數(shù)組中偶數(shù)元素為電壓模值U，奇數(shù)元素為P

n2=size(PV)/2;

numofPV=n2(1) %統(tǒng)計(jì)PV節(jié)點(diǎn)個(gè)數(shù)

end

if n==3

QVmaxmin=tline; %PV節(jié)點(diǎn)為節(jié)點(diǎn)2，總之PV節(jié)點(diǎn)的編號(hào)在PQ節(jié)點(diǎn)編號(hào)前面

QVmaxmin=sscanf(QVmaxmin,'%f') %PV數(shù)組中偶數(shù)元素為電壓模值U，奇數(shù)元素為P

end

if n==4

S=tline; %平衡節(jié)點(diǎn)為節(jié)點(diǎn)1

S=sscanf(S,'%f')

end

if n==5

Yb0=tline; %自導(dǎo)納

Yb0=sscanf(Yb0,'%f')

end

if n>=6

Yb=tline; %

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

num=numofPQ+numofPV+1 ; %num為系統(tǒng)節(jié)點(diǎn)數(shù)

for j=1:(numofPQ+numofPV+1) %1對(duì)應(yīng)的平衡節(jié)點(diǎn)

YBre(n-5,j)=Yb(2*j-1); %節(jié)點(diǎn)導(dǎo)納陣實(shí)部

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為節(jié)點(diǎn)導(dǎo)納矩陣

fclose(fidin);

%%%%%%%%%%%%%%%%上面部分為數(shù)據(jù)讀入及導(dǎo)納矩陣形成，下面為定義部分%%%%%%%%%%%%%%%%%%%%%%%%%

e=ones(num-1,1); %對(duì)除平衡節(jié)點(diǎn)外的節(jié)點(diǎn)電壓e和f賦初值，從標(biāo)號(hào)1到標(biāo)號(hào)num-1的數(shù)組內(nèi)放置的

f=zeros(num-1,1); %是從第2節(jié)點(diǎn)到第num節(jié)點(diǎn)的元素

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; %存放迭代次數(shù)

flagPVPQ=zeros(num-1,1); %標(biāo)記PV是否轉(zhuǎn)換為PQ

%%%%%%%%%%%%%%%%%%%%%%%%%%％％％％％％下面為潮流計(jì)算部分％％％％％％%%%%%%%%%%%%%%%%%%%%%%%%%

e1=S(1);f1=S(2); %1號(hào)節(jié)點(diǎn)為平衡節(jié)點(diǎn)

flowComputer=fopen('flowComputer.txt','wt'); % 創(chuàng)建flowComputer.txt文件

while number<=100

number=number %顯示迭代次數(shù)

fprintf(flowComputer,'迭代次數(shù)：%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節(jié)點(diǎn)

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節(jié)點(diǎn)

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

%無論P(yáng)Q節(jié)點(diǎn)還是PV節(jié)點(diǎn)，其有功和無功在后面的節(jié)點(diǎn)注入電流計(jì)算中均會(huì)用到

% 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節(jié)點(diǎn)的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)對(duì)應(yīng)的是QVmin

flagPVPQ(1,j-1)=1;

else

detaPQU(2*(j-1))=QVmaxmin(2*(j-1)-1)-sumQ; %QVmaxmin(2*(j-1)-1)對(duì)應(yīng)的是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節(jié)點(diǎn)的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));%除平衡節(jié)點(diǎn)外的節(jié)點(diǎn)都有注入電流,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節(jié)點(diǎn)是否要轉(zhuǎn)換為PQ節(jié)點(diǎn)

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兩個(gè)節(jié)點(diǎn)之間是否有直接的聯(lián)系

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的不平衡量達(dá)到指定精度，可停止迭代

display('This is an end');

break;

else

number=number+1;

fprintf(flowComputer,'經(jīng)過不平衡量修正后的電壓：\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 %對(duì)應(yīng)while循環(huán)

fclose(flowComputer);

%%%%%%%%%%%%%%%%%%%平衡節(jié)點(diǎn)功率及線路功率的計(jì)算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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)) %平衡節(jié)點(diǎn)功率

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

總結(jié)

以上是生活随笔為你收集整理的matlab牛顿拉夫逊算法,牛顿拉夫逊法程序设计的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。