sysselfljwd0.m

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

%该程序仿真两区域电力系统本身的联结稳定性。按推导出来的条件先测试系统是否联结稳定，然后看鲁棒程度的大小。
%系统数据(考虑alfa1i固定的情况）
A1=[-.2 0 0 0 0 0 0 -4;4.75 -5 0 0 0 0 0 0;0 0.1667 -0.1667 0 0 0 0 0;0 0 2 -2 0 0 0 0;...
   0 -0.08 -0.0747 -0.112 -3.994 10 -0.928 -9.1011;0 0 0 0 0.2 -0.5 0 0;0 0 0 0 1.3194 0 -1.3889 -0.2778;...
   0 0.01 0.0093 0.014 -0.0632 0 0.116 -0.1124];
A2=[-.2 0 0 0 0 0 0 -4;4.75 -5 0 0 0 0 0  0;0 0.1667 -0.1667 0 0 0 0 0;0 0 2 -2 0 0 0 0;...
   0 -0.1 -0.0933 -0.14 -4.096 10 -0.7442 -9.1079;0 0 0 0 0.2 -0.5 0 0;0 0 0 0 1.3194 0 -1.3889 -0.2778;...
   0 0.0125 0.0117 0.0175 -0.0506 0 0.0928 -0.1115];
at1=[0 0 0 0 0.6667 0 0 -0.0833]';at2=at1;
d1=[0 0 0 0 0 0 0 10]';d2=d1;
f1=at1;f2=at2;
k21=22.2144;k12=k21;
m12=[0 0 0 0 0 0 0 k21]';m21=m12;
b1=[1.6 0 0 0 6 0 0 0]';b2=[2 0 0 0 5 0 0 0]';
c1=[0 0.3 0.28 0.42 0 0 0 0;0 0 0 0 -1.52 0 2.78 0.217;zeros(1,7) 1];c2=c1;
alfa11=1
alfa12=1
A11=[A1 zeros(8,1) at1;d1' 0 1;m12' 0 0];
A12=[zeros(9,10);-m21' 0 0];
A21=[zeros(9,10);-alfa12*m12' 0 0];
A22=[A2 zeros(8,1) at2;d2' 0 1;alfa12*m21' 0 0];
B=[b1 zeros(8,1);zeros(2,2);zeros(8,1) b2;zeros(2,2)];
B11=[b1;0;0];
B22=[b2;0;0];
F=[f1 zeros(8,1);zeros(2,2);zeros(8,1) f2;zeros(2,2)];
C=[c1 zeros(3,12);zeros(1,8) 1 zeros(1,11);zeros(1,9) 1 zeros(1,10);zeros(3,10) c2 zeros(3,2);zeros(1,18) 1 0;zeros(1,19) 1];
D=zeros(10,2);

%由于子系统不稳定，致使Lyapunov方程无解，故给系统加LQ控制使其稳定
v=[0 0 0 0 0 0 0 4 4 1];
Q11=diag(v);
R11=1;
[K1,S1,E1]=lqr(A11,B11,Q11,R11);
Q22=Q11;
R22=R11;
[K2,S2,E2]=lqr(A22,B22,Q22,R22);
A11=A11-B11*K1;
A22=A22-B22*K2
A=[A11 A12;A21 A22];
A0=[A11 zeros(10,10);zeros(10,10) A22];

% %判断联结稳定性
%基本互联矩阵
e120=0.001
e210=0.001
I=eye(10);
mu1=1;
mu2=1;
II1=mu1*I;
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);
n1=10;n2=10;n=n1+n2;
w11=-mu1/(2*Q1)
w12=(2*(n-n1)*e120^2*norm(P1*A12)^2)/(Q2*n2*mu1)
w21=(2*(n-n2)*e210^2*norm(P2*A21)^2)/(Q1*n1*mu2)
w22=-mu2/(2*Q2)
% w11=-(mu1-e110*abs(alfa11-1)*norm(A111'*P1+P1*A111))/(2*Q1)
% w12=((n-n1)*e120^2*alfa11^2*norm(A12'*P1+P1*A12)^2)/(2*n2*Q2*(mu1-e110*abs(alfa11-1)*norm(A111'*P1+P1*A111)))
% w21=((n-n2)*e210^2*alfa12^2*norm(A21'*P2+P2*A21)^2)/(2*n1*Q1*(mu2-e220*abs(alfa12-1)*norm(A222'*P2+P2*A222)))
% w22=-(mu2-e220*abs(alfa12-1)*norm(A222'*P2+P2*A222))/(2*Q2)
W=[w11 w12;w21 w22]
Weig=eig(W)
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

%仿真系统在阶跃负载扰动下的响应
t=0:0.01:60;i=length(t);
w1=10*ones(1,i);w2=1*ones(1,i);u11=[w1;w2];
A00=[A11,e120*A12;e210*A21,A22];
sys11=ss(A00,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('Deviation of Pt1'),xlabel('t/s');
% subplot(222),plot(Ts11,y11(:,2),'-b',Ts00,y00(:,2),'-r'),title('Deviation of PH1'),xlabel('t/s');
% subplot(223),plot(Ts11,y11(:,3),'-b',Ts00,y00(:,3),'-r'),title('Deviation of frequency f1'),xlabel('t/s');
% subplot(224),plot(Ts11,y11(:,5),'-b',Ts00,y00(:,5),'-r'),title('Deviation of Pe1'),xlabel('t/s');
% subplot(221),plot(Ts11,y11(:,6),'-b',Ts00,y00(:,6),'-r'),title('Deviation of Pt2'),xlabel('t/s');
% subplot(222),plot(Ts11,y11(:,7),'-b',Ts00,y00(:,7),'-r'),title('Deviation of PH2'),xlabel('t/s');
% subplot(223),plot(Ts11,y11(:,8),'-b',Ts00,y00(:,8),'-r'),title('Deviation of frequency f2'),xlabel('t/s');
% subplot(224),plot(Ts11,y11(:,10),'-b',Ts00,y00(:,10),'-r'),title('Deviation of Pe2'),xlabel('t/s');


subplot(221),plot(Ts11,y11(:,3),'-b',Ts00,y00(:,3),'-r'),title('(a) \Delta{\itf1}响应曲线'),xlabel('Time(sec)'),ylabel('\Delta{\itf1/Hz}');
subplot(222),plot(Ts11,y11(:,5),'-b',Ts00,y00(:,5),'-r'),title('(b) \Delta{\itPe1}响应曲线'),xlabel('Time(sec)'),ylabel('\Delta{\itPe1}/p.u.');
subplot(223),plot(Ts11,y11(:,8),'-b',Ts00,y00(:,8),'-r'),title('(c) \Delta{\itf2}响应曲线'),xlabel('Time(sec)'),ylabel('\Delta{\itf2/Hz}');
subplot(224),plot(Ts11,y11(:,10),'-b',Ts00,y00(:,10),'-r'),title('(d) \Delta{\itPe2}响应曲线'),xlabel('Time(sec)'),ylabel('\Delta{\itPe2}/p.u.');

%求不确定参数eij范数界
B12b=[zeros(10,10),A12;zeros(10,10),zeros(10,10)];
B21b=[zeros(10,10),zeros(10,10);A21,zeros(10,10)];
Bb=[A11,zeros(10,10);zeros(10,10),A22];
Ib=eye(20);
mu=1;
IIb=mu*Ib;
Pb=lyap(Bb',IIb);
P12b=B12b'*Pb+Pb*B12b;
P21b=B21b'*Pb+Pb*B21b;
% Pr=[P12b,P21b];
Pr1=[zeros(20,20),P12b,P21b,zeros(20,20)];
% Qb0=norm(Pr)
Qb1=norm(Pr1);
Qb=Qb1^2;
eij2=mu^2/Qb
%第二个判断条件，看是否小于1
QQ1=norm(P12b);
QQ2=norm(P21b);
QQ=e120*QQ1+e210*QQ2;
if QQ<mu
    disp('The system is connectively stable.')
else
    disp('The system is not connectively stable.')
end
emax=mu/(QQ1+QQ2)