clear global
clear
clc
global gen_n h PT Tj f Pe E Xdp bus_n wn Ys ftx Yss Yg
global vxnum vynum dnum wnum gxnum gynum
load psat145ieee.mat
gen_num=Syn.con(:,1);
[A,ix]=sortrows(gen_num);
Syn.con=Syn.con(ix,:);
Syn.pm=Syn.pm(ix,:);
gen_num=Syn.con(:,1)';
load_num=PQ.con(:,1)';
bus_n=Bus.n;
gen_n=Syn.n;
load_n=PQ.n;
f=60; wn=2*pi*f;
Xdp=Syn.con(:,9); %列向量
Tj=Syn.con(:,18);
PT=Syn.pm; %列向量
default_n=[7];
default_bra_num=[7 6];
default_bra_b=Line.con(16,10);
Vm=DAE.V;
alf=DAE.a;
V=Vm.*(cos(alf)+j*sin(alf));
num=gen_num;
S=Bus.Pg(num)+j*Bus.Qg(num);
Ep=V(num)+j*Xdp.*conj(S./V(num)); %求相量E'
E=abs(Ep);
d0=angle(Ep); %列向量
w0=ones(gen_n,1); %列向量
num=load_num;
Y1=Line.Y;
Y1=Y1+sparse(num,num,conj(PQ.P0+j*PQ.Q0)./(Vm(num).^2),bus_n,bus_n); %包含负荷节点导纳Y
Y1=Y1+sparse(gen_num,gen_num,-j*1./Xdp,bus_n,bus_n);
Y2=Y1;
Y2(default_n,default_n)=Y2(default_n,default_n)+1e20; %三相接地故障
Y3=Y1;
m=default_bra_num(1); n=default_bra_num(2);
Y3(m,m)=Y3(m,m)+Y1(m,n)-j*default_bra_b/2; Y3(m,n)=Y3(m,n)-Y1(m,n);
Y3(n,m)=Y3(m,n); Y3(n,n)=Y3(n,n)+Y1(m,n)-j*default_bra_b/2;%切除故障
for i=1:bus_n
for j=1:bus_n
Ys1{i,j}=[real(Y1(i,j)) -imag(Y1(i,j));imag(Y1(i,j)) real(Y1(i,j))];
Ys2{i,j}=[real(Y2(i,j)) -imag(Y2(i,j));imag(Y2(i,j)) real(Y2(i,j))];
Ys3{i,j}=[real(Y3(i,j)) -imag(Y3(i,j));imag(Y3(i,j)) real(Y3(i,j))];
end
end
Ys1=sparse(cell2mat(Ys1)); Ys2=sparse(cell2mat(Ys2)); Ys3=sparse(cell2mat(Ys3));
vxnum=1:2:2*bus_n-1; vxnum=vxnum';
vynum=2:2:2*bus_n; vynum=vynum';
dnum=1:2:2*gen_n-1; dnum=dnum';
wnum=2:2:2*gen_n; wnum=wnum';
gxnum=2*gen_num-1; gxnum=gxnum'; %发电机节点
gynum=2*gen_num; gynum=gynum';
Yg=sparse(gxnum,gynum,1./Xdp,2*bus_n,2*bus_n)+sparse(gynum,gxnum,-1./Xdp,2*bus_n,2*bus_n);
x(dnum,1)=d0; x(wnum,1)=w0; y(vxnum,1)=real(V); y(vynum,1)=imag(V);
ftx=1;
if ftx==1
x(wnum,1)=w0-1;
end
h1=0.18;
t_back1=0;
t_back2=0.1;
t_end=1.5;
t(1)=0;
h=h1;
Ys=Ys1;
xresult(:,1)=x(dnum);
i=1;
t0=cputime;
while t(i)<t_end
t(i+1)=t(i)+h1;
if t(i+1)>t_back1 & t(i)<t_back1
h=t_back1-t(i);
t(i+1)=t_back1;
end
if t(i+1)==t_back1+h1
Ys=Ys2; h=h1; Yss=Ys2; y=scpf(x);
end
if t(i+1)>t_back2 & t(i)<t_back2
h=t_back2-t(i);
t(i+1)=t_back2;
end
if t(i+1)==t_back2+h1
Ys=Ys3; h=h1; Yss=Ys3; y=scpf(x);
end
if t(i+1)>t_end&& t(i)<t_end
h=t_end-t(i);
t(i+1)=t_end;
end
hhh(i)=h;
[x,y]=dw(Tj,Ys,x,y,gen_n,h,ftx); %解微分方程
xresult(:,i+1)=x(dnum);
i=i+1;
end
et = cputime-t0
clear x y
x(dnum,1)=d0; x(wnum,1)=w0; y(vxnum,1)=real(V); y(vynum,1)=imag(V);
ftx=2;
t_back1=0;
t_back2=0.1;
t_end=1.5;
t(1)=0;
h=h1;
Ys=Ys1;
xxresult(:,1)=x(dnum);
i=1;
while t(i)<t_end
t(i+1)=t(i)+h1;
if t(i+1)>t_back1 & t(i)<t_back1
h=t_back1-t(i);
t(i+1)=t_back1;
end
if t(i+1)==t_back1+h1
Ys=Ys2; h=h1; Yss=Ys2; y=scpf(x);
end
if t(i+1)>t_back2 & t(i)<t_back2
h=t_back2-t(i);
t(i+1)=t_back2;
end
if t(i+1)==t_back2+h1
Ys=Ys3; h=h1; Yss=Ys3; y=scpf(x);
end
if t(i+1)>t_end& t(i)<t_end
h=t_end-t(i);
t(i+1)=t_end;
end
hhh(i)=h;
[x,y]=dw(Tj,Ys,x,y,gen_n,h,ftx); %解微分方程
xxresult(:,i+1)=x(dnum);
i=i+1;
end
% load 3dresult.mat
% xx(3).t=t;
% xx(3).xresult=xresult;
% xx(3).xxresult=xxresult;
% xx(3).h=h1*ones(length(t),1);
load result.mat
ind=int16(t*1000+1);
%subplot(324)
plot(t',abs(xresult(20,:)-xresult(50,:)-result1(20,ind)+result1(50,ind))*180/pi,'r',...
t',abs(xxresult(20,:)-xxresult(50,:)-result1(20,ind)+result1(50,ind))*180/pi,'b:')
legend('explicit RKN method','classic RK method',2)
legend('boxoff')
set(gca,'FontSize',9);
xlabel('时间/s','fontsize',9);
ylabel(['误差/( ','\circ',')'],'fontsize',9);
% plot(t',(result(20,ind)-result(50,ind))*180/pi,'-k',t',(xresult(20,:)-xresult(50,:))*180/pi,':or',...
% t',(xxresult(20,:)-xxresult(50,:))*180/pi,':xb')
% legend('精确值','3级4阶辛RKN法','传统4阶RK法')
% xlabel('时间(s)');
% ylabel('功角(゜)');
% ylim([-10 150])
% tx=xresult;
% subplot(311)
% plot(t',xresult(20,:)'-xresult(50,:)','r',t',xresult(26,:)'-xresult(50,:)','g')
% xlim([0 t_end])
% subplot(312)
% plot(t',xresult(20,:)'-xresult(26,:)')
% xlim([0 t_end])
% load result3.mat
% ind=int16(t*1000+1);
% subplot(324)
% plot(t',abs(xresult(20,:)-xresult(50,:)-result(20,ind)+result(50,ind))*180/pi,'r')
% xlabel('时间(s)');
% ylabel('误差(゜)');
% if ftx==1
% title('3级4阶辛rkn法');
% end
% if ftx==2
% title('传统4阶rk法');
% end
% a=[0.03 1.125 1.422
% 0.04 0.828 1.047
% 0.05 0.672 0.859
% 0.06 0.578 0.734
% 0.07 0.515 0.641
% 0.08 0.468 0.594
% 0.09 0.422 0.532
% 0.10 0.343 0.438
% 0.11 0.313 0.391
% 0.12 0.296 0.375
% 0.13 0.271 0.344
% 0.14 0.250 0.312
% 0.15 0.250 0.312
% 0.16 0.234 0.297
% 0.17 0.234 0.297
% 0.18 0.218 0.273];
%
% subplot(324)
% plot(a(:,1),a(:,2),'r-',a(:,1),a(:,3),'b--','LineWidth',2.0)
% xlabel('步长(s)')
% ylabel('CPU时间(s)')
% legend('3级4阶辛RKN法','传统4阶RK法')
% xlim([0.03 0.18])
% subplot(325)
% plot(a(:,2),a(:,3))
% load 3dresult.mat
% for i=1:13
% ind=int16(xx(i).t*1000+1);
% plot3(xx(i).h,xx(i).t,abs(xx(i).xresult(20,:)-xx(i).xresult(50,:)-result(20,ind)+result(50,ind))*180/pi,'r',...
% xx(i).h,xx(i).t,abs(xx(i).xxresult(20,:)-xx(i).xxresult(50,:)-result(20,ind)+result(50,ind))*180/pi,'b')
% hold on
% end
% xlabel('步长(s)')
% ylabel('时间(s)')
% zlabel('误差(゜)')
% clear
% clc
% c1=0.757434; c2=0.5;
% la=0.5*c1^2;
% a=diag([la,la,la]);
% a(2,1)=(c2^2-c1^2)/2;
% a(3,1)=(c1+c2)*(1-2*c1)/2;
% a(3,2)=(1-2*c1)*(1-c1-c2)/2;
% inva=inv(a)
% c1=0.138824; c2=0.5;
% la=0.5*c1^2;
% b=diag([la,la,la]);
% b(2,1)=(c2^2-c1^2)/2;
% b(3,1)=(c1+c2)*(1-2*c1)/2;
% b(3,2)=(1-2*c1)*(1-c1-c2)/2;
% invb=inv(b)
评论0