求MATLAB簡單潮流計算程式,求MATLAB簡單潮流計算程式。。。。

時間 2021-09-02 08:26:16

1樓:匿名使用者

function lianxuchaoliu

clear;

clc;

n=9;%節點數;

nl=9;%支路數;

isb=1;%平衡節點號;

pr=0.00001;%誤差精度;

b1=[1 4 0.0576i 0 1.05 1; 4 5 0.

017+0.092i 0.158i 1 0; 5 6 0.

039+0.17i 0.358i 1 0; 3 6 0.

0586i 0 1.05 1; 6 7 0.0119+0.

1008i 0.209i 1 0; 7 8 0.0085+0.

072i 0.149i 1 0; 8 2 0.0625i 0 1.

05 1; 8 9 0.032+0.161i 0.

306i 1 0; 9 4 0.01+0.085i 0.

176i 1 0];

%依次是支路首端;末端,支路阻抗;對地電納;支路變比;折算到哪一側標誌(高壓側為1;低壓側為0);

b2=[0 0 1.05 1.05 0 1; 1.

63 0 1.05 1.05 0 3; 0.

85 0 1.05 1.05 0 3; 0 0 1 0 0 2; 0 0.

9+0.3i 1 0 0 2; 0 0 1 0 0 2; 0 1+0.35i 1 0 0 2; 0 0 1 0 0 2; 0 1.

25+0.5i 1 0 0 2];

%依次是節點的發電機功率;負荷功率;節點電壓初值;pv節點電壓v給定值;節點無功補償裝置容量;節點分類標號(平衡1;pq2;pv3);

y=zeros(n);%求導納陣;

for i=1:nl

if b1(i,6)==0

p=b1(i,1);q=b1(i,2);

else

p=b1(i,2);q=b1(i,1);

endy(p,q)=y(p,q)-1./(b1(i,3)*b1(i,5));

y(q,p)=y(p,q);

y(q,q)=y(q,q)+1./(b1(i,3)*b1(i,5)^2)+b1(i,4)./2;

y(p,p)=y(p,p)+1./b1(i,3)+b1(i,4)./2;

end%disp('系統的導納陣為:');

%disp(y);

g=real(y);b=imag(y);

for i=1:n

e(i)=real(b2(i,3));

f(i)=imag(b2(i,3));

v(i)=b2(i,4);

endfor i=1:n

s(i)=b2(i,1)-b2(i,2);

b(i,i)=b(i,i)+b2(i,5);

endp=real(s);q=imag(s);

w=zeros(2*n,1);jac=zeros(2*n);

t=0;

while t==0

for i=1:n

if b2(i,6)~=isb

c=0;d=0;

for j=1:n

c=c+g(i,j)*e(j)-b(i,j)*f(j);

d=d+g(i,j)*f(j)+b(i,j)*e(j);

endif b2(i,6)==2%p,q節點;

w(2*i)=p(i)-e(i)*c-f(i)*d;

w(2*i-1)=q(i)-f(i)*c+e(i)*d;

else if b2(i,6)==3%p,v節點;

w(2*i)=p(i)-e(i)*c-f(i)*d;

w(2*i-1)=v(i)^2-(e(i)^2+f(i)^2);

endendelse

w(2*i-1)=0;

w(2*i)=0;

endend%disp(w);

w1=w(3:2*n);

for i=1:n

for j=1:n

if b2(i,6)~=isb

if b2(i,6)==2%p,q節點;

if j~=i

jac(2*i,2*j-1)=-1*(g(i,j)*e(i)+b(i,j)*f(i));

jac(2*i-1,2*j)=(g(i,j)*e(i)+b(i,j)*f(i));

jac(2*i,2*j)=b(i,j)*e(i)-g(i,j)*f(i);

jac(2*i-1,2*j-1)=b(i,j)*e(i)-g(i,j)*f(i);

else if j==i

m=0;h=0;

for r=1:n

m=m+g(i,r)*e(r)-b(i,r)*f(r);

h=h+g(i,r)*f(r)+b(i,r)*e(r);

endjac(2*i,2*j-1)=-1*m-g(i,i)*e(i)-b(i,i)*f(i);

jac(2*i-1,2*j)=-1*m+g(i,i)*e(i)+b(i,i)*f(i);

jac(2*i,2*j)=-1*h+b(i,i)*e(i)-g(i,i)*f(i);

jac(2*i-1,2*j-1)=h+b(i,i)*e(i)-g(i,i)*f(i);

endendelse if b2(i,6)==3%p,v節點;

if j~=i

jac(2*i,2*j-1)=-1*(g(i,j)*e(i)+b(i,j)*f(i));

jac(2*i-1,2*j)=0;

jac(2*i,2*j)=b(i,j)*e(i)-g(i,j)*f(i);

jac(2*i-1,2*j-1)=0;

else if j==i

m=0;h=0;

for r=1:n

m=m+g(i,r)*e(r)-b(i,r)*f(r);

h=h+g(i,r)*f(r)+b(i,r)*e(r);

endjac(2*i,2*j-1)=-1*m-g(i,i)*e(i)-b(i,i)*f(i);

jac(2*i-1,2*j)=-2*f(i);

jac(2*i,2*j)=-1*h+b(i,i)*e(i)-g(i,i)*f(i);

jac(2*i-1,2*j-1)=-2*e(i);

endendendendelse

jac(2*i-1,2*j-1)=0;

jac(2*i,2*j)=0;

jac(2*i-1,2*j)=0;

jac(2*i,2*j-1)=0;

endendend%disp(jac);

jac1=jac(3:2*n,3:2*n);

for k=1:2*n-2

m=0;

for i=k+1:2*n-2

m=jac1(i,k)./jac1(k,k);

w1(i)=w1(i)-m*w1(k);

for j=k+1:2*n-2

jac1(i,j)=jac1(i,j)-m*jac1(k,j);

endendende(2*n-2)=w1(2*n-2)./jac1(2*n-2,2*n-2);

for i=2*n-3:-1:1

c=0;

for k=i+1:2*n-2

c=c+jac1(i,k)*e(k);

e(i)=(w1(i)-c)./jac1(i,i);

endend%disp(e);

for i=1:n-1

e(i+1)=e(i+1)-e(2*i-1);

f(i+1)=f(i+1)-e(2*i);

end%disp(e);

%disp(f);

for i=1:n-1

b(2*i-1)=abs(e(2*i-1));

b(2*i)=abs(e(2*i));

endkb=max(b);

%disp(kb);

if kb1%以切向量中分量最大(絕對值最大)的狀態變數選定為連續引數;

for i=1:2*n-2

jd(i)=abs(d(i));

endfor i=1:2*n-3

if jd(i)>=jd(i+1)

zd=jd(i);ek=i;

else if jd(i)<=jd(i+1)

zd=jd(i+1);ek=i+1;

endendendfor j=1:2*n-1

jj1(2*n-1,j)=0;

endjj1(2*n-1,ek)=1;

w2=zeros(2*n-1,1);

for i=1:2*n-1

w2(i,1)=0;

endif d(ek)>0

w2(ek,1)=1;

else if d(ek)<0

w2(ek,1)=-1;

endendendfor k=1:2*n-1

m=0;

for i=k+1:2*n-1

m=jj1(i,k)./jj1(k,k);

w2(i)=w2(i)-m*w2(k);

for j=k+1:2*n-1

jj1(i,j)=jj1(i,j)-m*jj1(k,j);

endendendd(2*n-1)=w2(2*n-1)./jj1(2*n-1,2*n-1);

for i=2*n-2:-1:1

c=0;

for k=i+1:2*n-1

c=c+jj1(i,k)*d(k);

d(i)=(w2(i)-c)./jj1(i,i);

endend%disp(d);

for i=1:n-1

e(i+1)=e(i+1)+l*d(2*i-1);

f(i+1)=f(i+1)+l*d(2*i);

endt=t+l*d(2*n-1);

disp(d(2*n-1));

%disp(e);

%disp(f);

%disp(t);

%對預估的近似解進行校正;

tt=0;

while tt==0

for i=1:n

if b2(i,6)~=isb

c=0;d=0;

for j=1:n

c=c+g(i,j)*e(j)-b(i,j)*f(j);

d=d+g(i,j)*f(j)+b(i,j)*e(j);

endif b2(i,6)==2% p,q節點;

if b2(i,2)~=0

wj(2*i)=p(i)-t*real(b2(i,2))-e(i)*c-f(i)*d;

wj(2*i-1)=q(i)-t*imag(b2(i,2))-f(i)*c+e(i)*d;

else

wj(2*i)=p(i)-e(i)*c-f(i)*d;

wj(2*i-1)=q(i)-f(i)*c+e(i)*d;

endelse if b2(i,6)==3%p,v節點;

wj(2*i)=p(i)+t*1.5*1.05-e(i)*c-f(i)*d;

wj(2*i-1)=0;

endendelse

wj(2*i)=0;

wj(2*i-1)=0;

endendwj1=wj(3:2*n);

for i=1:n

for j=1:n

if b2(i,6)~=isb

if b2(i,6)==2%p,q節點;

if j~=i

jjac(2*i,2*j-1)=-1*(g(i,j)*e(i)+b(i,j)*f(i));

jjac(2*i-1,2*j)=(g(i,j)*e(i)+b(i,j)*f(i));

jjac(2*i,2*j)=b(i,j)*e(i)-g(i,j)*f(i);

jjac(2*i-1,2*j-1)=b(i,j)*e(i)-g(i,j)*f(i);

else if j==i

m=0;h=0;

for r=1:n

m=m+g(i,r)*e(r)-b(i,r)*f(r);

h=h+g(i,r)*f(r)+b(i,r)*e(r);

endjjac(2*i,2*j-1)=-1*m-g(i,i)*e(i)-b(i,i)*f(i);

jjac(2*i-1,2*j)=-1*m+g(i,i)*e(i)+b(i,i)*f(i);

jjac(2*i,2*j)=-1*h+b(i,i)*e(i)-g(i,i)*f(i);

jjac(2*i-1,2*j-1)=h+b(i,i)*e(i)-g(i,i)*f(i);

endendelse if b2(i,6)==3%p,v節點;

if j~=i

jjac(2*i,2*j-1)=-1*(g(i,j)*e(i)+b(i,j)*f(i));

jjac(2*i-1,2*j)=0;

jjac(2*i,2*j)=b(i,j)*e(i)-g(i,j)*f(i);

jjac(2*i-1,2*j-1)=0;

else if j==i

m=0;h=0;

for r=1:n

m=m+g(i,r)*e(r)-b(i,r)*f(r);

h=h+g(i,r)*f(r)+b(i,r)*e(r);

endjjac(2*i,2*j-1)=-1*m-g(i,i)*e(i)-b(i,i)*f(i);

jjac(2*i-1,2*j)=-2*f(i);

jjac(2*i,2*j)=-1*h+b(i,i)*e(i)-g(i,i)*f(i);

jjac(2*i-1,2*j-1)=-2*e(i);

endendendendelse

jjac(2*i-1,2*j-1)=0;

jjac(2*i,2*j)=0;

jjac(2*i-1,2*j)=0;

jjac(2*i,2*j-1)=0;

endendendjjac1=jjac(3:2*n,3:2*n);

for k=1:2*n-2

m=0;

for i=k+1:2*n-2

m=jjac1(i,k)./jjac1(k,k);

wj1(i)=wj1(i)-m*wj1(k);

for j=k+1:2*n-2

jjac1(i,j)=jjac1(i,j)-m*jjac1(k,j);

endendende1(2*n-2)=wj1(2*n-2)./jjac1(2*n-2,2*n-2);

for i=2*n-3:-1:1

c=0;

for k=i+1:2*n-2

c=c+jjac1(i,k)*e1(k);

e1(i)=(wj1(i)-c)./jjac1(i,i);

endend%disp(e1);

for i=1:n-1

e(i+1)=e(i+1)-e1(2*i-1);

f(i+1)=f(i+1)-e1(2*i);

endfor i=1:n-1

bx(2*i-1)=abs(e1(2*i-1));

bx(2*i)=abs(e1(2*i));

endkb1=max(bx);

if kb10

ttt=0;ccc=ccc+1;

else

ttt=1;

endend

disp(t);

disp(ccc);

%disp(e);

%disp(f);

for i=1:n

fz1(i)=sqrt(e(i)^2+f(i)^2);

enddisp(fz1);

求matlab程式,給定4組資料x1 2,3,7,6,5 1,5,8,

建立模型y f x1,x2,x3,x4 因為y是關於x1,x2,x3,x4的線性函式,所以有f x1,x2,x3,x4 a0 a1 x1 a2 x2 a3 x3 a4 x4,x1 2,3,7,6,5 x2 1,5,8,6,7 x3 4,7,9,12,10 x4 5,8,10,13,11 y 23,2...

設計matlab程式計算圓域上的二重積分

這個可以用matlab的符號積分或者數值積分解決,下面提供4種方法 1 直角座標系符號積分 syms x y int y int sin pi x 2 y 2 y,sqrt 1 x 2 sqrt 1 x 2 先對y積分 i vpa int int y,x,1,1 i 2.0 2 極座標系符號積分 s...

求簡單的金工實習數控車床簡單程式

funac數控車程式設計如下 o9001 n10 g50 x100 z10 設立座標系,定義對刀點的位置 n20 g00 x16 z2 m03 移到倒角延長線,z 軸2mm 處 n30 g01 u10 w 5 g98 f120 倒3 45 角 n40 z 48 加工 26 外圓 n50 u34 w ...