sysselfljwd.m

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

%该程序仿真两区域电力系统本身的联结稳定性。按推导出来的条件先测试系统是否联结稳定，然后看鲁棒程度的大小。
%系统数据
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=0.3
A11=[A1 zeros(8,1) at1;d1' 0 1;m12' 0 0];
A111=[zeros(9,10);m12' 0 0];
A12=[zeros(9,10);-m21' 0 0];
A21=[zeros(9,10);-m12' 0 0];
A22=[A2 zeros(8,1) at2;d2' 0 1;m21' 0 0];
A222=[zeros(9,10);m21' 0 0];
A1111=[zeros(9,10);(alfa11-1)*m12' 0 0];
A110=A11+A1111;
A2222=[zeros(9,10);(alfa12-1)*m21' 0 0];
A220=A22+A2222;
A120=[zeros(9,10);-alfa11*m21' 0 0];
A210=[zeros(9,10);-alfa12*m12' 0 0];
A=[A110 A120;A210 A220];
A0=[A110 zeros(10,10);zeros(10,10) A220];
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 4];
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;
% A11=[A110(1:9,1:10);zeros(1,10)];
% A111=[zeros(9,10);A110(10,1:10)];
% A22=[A220(1:9,1:10);zeros(1,10)];
% A222=[zeros(9,10);A220(10,1:10)];
A110=A11+A1111;
A220=A22+A2222;
A=[A110 A120;A210 A220];
A0=[A110 zeros(10,10);zeros(10,10) A220];
% A110=[A11(1:9,1:10);zeros(1,8) 0 0];
% A220=[A22(1:9,1:10);zeros(1,8) 0 0];
% A0=[A110,zeros(10,10);zeros(10,10),A220];
% %判断联结稳定性
e110=1;e120=1;e210=1;e220=1;%基本互联矩阵
% I=eye(10);
% mu1=10;
% mu2=10;
% 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;
% E1=eig(A111'*P1+P1*A111);
% E1max=max(E1)
% if E1max<0
%     E1max=0;
% end
% E2=eig(A222'*P2+P2*A222);
% E2max=max(E2)
% if E2max<0
%     E2max=0;
% end
% w11=-(mu1-e110*(alfa11-1)*E1max-e120*alfa11*norm(P1*A12))/Q1
% w12=(e120*alfa11*norm(P1*A12))/Q2
% w21=(e210*alfa12*norm(P2*A21))/Q1
% w22=-(mu2-e220*(alfa12-1)*E2max-e210*alfa12*norm(P2*A21))/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=100*ones(1,i);w2=10*ones(1,i);u11=[w1;w2];
A00=[A110,e120*A120;e210*A210,A220];
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('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(221),plot(Ts11,y11(:,3),'-b',Ts00,y00(:,3),'-r'),title('Deviation of frequency f1'),xlabel('t/s');
subplot(222),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} output Response curve'),xlabel('Time(sec)'),ylabel('\Delta{\itf1/Hz}');
% subplot(222),plot(Ts11,y11(:,5),'-b',Ts00,y00(:,5),'-r'),title('(b) \Delta{\itPe1} output Response curve'),xlabel('Time(sec)'),ylabel('\Delta{\itPe1}/p.u.');
% subplot(221),plot(Ts11,y11(:,8),'-b',Ts00,y00(:,8),'-r'),title('(c) \Delta{\itf2} output Response curve'),xlabel('Time(sec)'),ylabel('\Delta{\itf2/Hz}');
% subplot(222),plot(Ts11,y11(:,10),'-b',Ts00,y00(:,10),'-r'),title('(d) \Delta{\itPe2} output Response curve'),xlabel('Time(sec)'),ylabel('\Delta{\itPe2}/p.u.');

%求不确定参数eij范数界
B11b=[A111,zeros(10,10);zeros(10,10),zeros(10,10)];
B12b=[zeros(10,10),A12;zeros(10,10),zeros(10,10)];
B21b=[zeros(10,10),zeros(10,10);A21,zeros(10,10)];
B22b=[zeros(10,10),zeros(10,10);zeros(10,10),A222];
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)/mu;
% P21b=(B21b'*Pb+Pb*B21b)/mu;
% Pr=[P12b,P21b];
% Pr1=[zeros(20,20),P12b,P21b,zeros(20,20)];
% Qb0=norm(Pr)
% Qb1=norm(Pr1)
% Qb=Qb0^2
% eij2=1/Qb
% %第二个判断条件，看是否小于1
% QQ1=norm(P12b);
% QQ2=norm(P21b);
% QQ=e120*QQ1+e210*QQ2
P11b=B11b'*Pb+Pb*B11b;
P22b=B22b'*Pb+Pb*B22b;
P12b=B12b'*Pb+Pb*B12b;
P21b=B21b'*Pb+Pb*B21b;
Pr=[P11b,P12b,P21b,P22b];
Qb=norm(Pr);
Qb0=Qb^2;
QQ11=norm(P11b);
QQ22=norm(P22b);
QQ12=norm(P12b);
QQ21=norm(P21b);
L=e110*alfa11*QQ11+e120*alfa11*QQ12+e210*alfa12*QQ21+e220*alfa12*QQ22;
aa=eig(e110*QQ11+e220*QQ22);
aamin=min(aa);
R=mu+aamin;
if L<R
    disp('The system is connectively stable.')
else
    disp('The system is not connectively stable.')
end
aa1=eig(QQ11+QQ22);
aa1min=min(aa1);
alfamax=(mu+aa1min)/(QQ12+QQ21+QQ11+QQ22)
alfa2=(mu+aamin)^2/Qb0