sysselfljwd11.m

这是一个研究生课程的matlab课件

%系统参数
A11=[-0.05 6.0 0.0 0.0 -6.0;0.0 -3.33 3.33 0.0 0.0;-5.2083 0.0 -12.5 -12.5 0.0;
    0.423 0.0 0.0 0.0 1.0;0.545 0.0 0.0 0.0 0.0];
A12=[zeros(4,5);-0.545 0 0 0 0];
B11=[0;0;12.5;0;0];
A22=[-0.05 6.0 0.0 0.0 -6.0;0.0 -3.33 3.33 0.0 0.0;-5.2083 0.0 -12.5 -12.5 0.0;
    0.423 0.0 0.0 0.0 1.0;0.545 0.0 0.0 0.0 0.0];
A21=[zeros(4,5);-0.545 0 0 0 0];
B22=[0;0;10.5;0;0];
A=[A11 A12;A21 A22];
B=[B11 zeros(5,1) ; zeros(5,1) B22];
F=[-6.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0;
   0.0 0.0 0.0 0.0 0.0 -6.0 0.0 0.0 0.0 0.0]';
C=eye(10);D=zeros(10,2);
A0=[A11 zeros(5,5);zeros(5,5) A22 ];
%A1=[A11 A12;zeros(5,5) A22 ];

%判断联结稳定性(定理1方法）
%基本互联矩阵
e120=0.2;e210=0.01;
I=eye(5);
mu1=1;
II1=mu1*I;
mu2=1;
II2=mu2*I;
P1=lyap(A11',II1);
P2=lyap(A22',II2);
P10=eig(P1);
P20=eig(P2);
Q1=min(P10);
Q2=min(P20);
Q3=max(P10);
Q4=max(P20);
w11=-mu1/(2*Q1);
n1=5;n2=5;n=n1+n2;
w12=((n-n1)*e120^2*norm(A12'*P1+P1*A12)^2)/(2*n2*Q2*mu1);
w21=((n-n2)*e210^2*norm(A21'*P2+P2*A21)^2)/(2*n1*Q1*mu2);
w22=-mu2/(2*Q2);
W=[w11 w12;w21 w22]
T1=w11;
T2=det(W);
if T1<0 & T2>0
    disp('The system is connectively stable.')
else
    disp('The system is not connectively stable.')
end
%判断联结稳定性(定理2方法）
B12b=[zeros(5,5),A12;zeros(5,5),zeros(5,5)];
B21b=[zeros(5,5),zeros(5,5);A21,zeros(5,5)];
Bb=[A11,zeros(5,5);zeros(5,5),A22];
Ib=eye(10);
mu=100;
IIb=mu*Ib;
Pb=lyap(Bb',IIb);
P12b=(B12b'*Pb+Pb*B12b)/mu;
P21b=(B21b'*Pb+Pb*B21b)/mu;
QQ1=norm(P12b);
QQ2=norm(P21b);
QQ=e120*QQ1+e210*QQ2
%计算鲁棒界
Pr=[zeros(10,10),P12b,P21b,zeros(10,10)];
Qb=norm(Pr)^2;
eij2=1/Qb
eijmax=1/norm(P12b+P21b)
%仿真系统在阶跃负载扰动下的响应
t=0:0.01:60;i=length(t);
w1=0.5*ones(1,i);w2=1*ones(1,i);u11=[w1;w2];
A=[A11 e120*A12;e210*A21 A22];
sys11=ss(A,F,C,D);
[y11,Ts11,x11]=lsim(sys11,u11,t);
t=60:0.01:150;i=length(t);
w1=0*ones(1,i);w2=0*ones(1,i);u00=[w1;w2];
x0=x11(6001,:);
sys00=ss(A0,F,C,D);
[y00,Ts00,x00]=lsim(sys00,u00,t,x0);

subplot(221),plot(Ts11,y11(:,1),'-b',Ts00,y00(:,1),'-r'),title('(a) 区域1的频率偏差曲线'),xlabel('Time(sec)'),ylabel('△f1 / Hz');
% subplot(222),plot(Ts11,y11(:,5),'-b',Ts00,y00(:,5),'-r'),title('(b) Deviation in Pe1'),xlabel('Time(sec)'),ylabel('p.u.');
subplot(222),plot(Ts11,y11(:,6),'-b',Ts00,y00(:,6),'-r'),title('(c) 区域2的频率偏差曲线'),xlabel('Time(sec)'),ylabel('△f2 / Hz');
% subplot(224),plot(Ts11,y11(:,10),'-b',Ts00,y00(:,10),'-r'),title('(d) Deviation in Pe2'),xlabel('Time(sec)'),ylabel('p.u.');

%plot(Ts11,y11(:,1),'-b',Ts00,y00(:,1),'-r'),xlabel('Time(sec)'),ylabel('△f1 / Hz');
%plot(Ts11,y11(:,6),'-b',Ts00,y00(:,6),'-r'),xlabel('Time(sec)'),ylabel('△f2 / Hz');

% t=0:0.01:100;i=length(t);
% w1=1*ones(1,i);
% u11=w1;
% F1=[-6 0 0 0 0]';
% C1=eye(5);D1=zeros(5,1);
% sys11=ss(A11,F1,C1,D1);
% [y11,Ts11,x11]=lsim(sys11,u11,t);
% subplot(221),plot(Ts11,y11(:,1),'-b'),title('(a) Deviation in frequency f1'),xlabel('Time(sec)'),ylabel('Hz');
% subplot(222),plot(Ts11,y11(:,5),'-b'),title('(b) Deviation in Pe1'),xlabel('Time(sec)'),ylabel('p.u.');