one.m

产生 S型函数产生 S型函数产生 S型函数产生 S型函数

clear all
v=randn(100,1);%%产生随机数
lam=0.01;T=0.1;
na=2;nb=2;nc=0;nk=1;
%a=[1,-1.5,0.7];b=[0,1,0.5];c=[1];d=[1];f=[1]; 
a=[1,-0.8,-1];b=[0,1,1.5];c=[1];d=[1];f=[1]; 
%%na=2;nb=2;nc=0;nk=1;
%%a=[1,-1.6,0.6];b=[0,1.5,-0.53,0.9];c=[1,0.4,0];d=[1];f=[1];
th0=poly2th(a,b,c,d,f,lam,T);
%u=idinput(600,'prbs');
u=idinput(100,'prbs');
y=idsim([u,v],th0);
z=[y,u];
%idplot(z);hold on;
%%模型参数辨识
adg=0.98;
th1=rarx(z,[na nb nk],'ff',adg);
%plot(th1)
c=th1';
a1=c(1,:); a2=c(2,:); b1=c(3,:); b2=c(4,:);
a11=c(1,100);a12=c(2,100);b11=c(3,100);b12=c(4,100);
y(1)=0; y(2)=0; y(3)=0;du(1)=0; du(2)=0; du(3)=0;Yr=[0 0 0]';Y=[0 0 0]';dU=[0 0 0]';
u(1)=0; u(2)=0; u(3)=0;
[yr,k]=gensig('sqare',50,403,1);%产生方波，周期50
%u=idinput(400,'prbs');
v=randn(400,1);
for k=3:400
    y(k)=-a11*y(k-1)+(-a12)*y(k-2)+b11*u(k-1)+b12*u(k-2)+0.01*v(k);
    C=[a11 a12 b11 b12];    
    A1=1-C(1);A2=C(1)-C(2);A3=C(2);B1=C(3);B2=C(4)-C(3);B3=-C(4);
    %A1=1.7;A2=-0.7;A3=0;
    %z(k)=A1*z(k-1)+A2*z(k-2)+A3*z(k-3)+B1*u(k-1)+B2*u(k-2)+B3*u(k-3);
    E1=[1 0 0];
    g10=A1;g11=A2;g12=A3;
    G1=[g10 g11 g12];
    e21=g10;
    E2=[e21 1 0];
    g20=g11+g10*A1;g21=g12+g10*A2;g22=g10*A3;
    G2=[g20 g21 g22];
    e32=g20;e31=e21;
    E3=[e32 e21 1];                                                                                                                                    
    g30=g21+g20*A1;g31=g22+g20*A2;g32=g20*A3;
    G3=[g30 g31 g32];
    l10=C(3);
    L1=l10;
    h10=C(4);                                                                   
    H1=h10 ;
    l21=e21*C(3)+h10;
    L2=l21;
    h20=e21*C(4);
    H2=h20;
    l32=e32*C(3)+h20;
    L3=l32;                                                                                                                                                                             
    h30=e32*C(4);
    H3=h30;
    L=[L1 L2 L3]';H=[H1 H2 H3]';G=[G1' G2' G3']';
 
    Q=eye(3,3);R=5;
    Lm=inv(L'*Q*L+R);
    Lq=Lm*(L'*Q);
    lq=Lq(1,:);

    Yr=[yr(k+1) yr(k+2) yr(k+3)]';dU=du(k-1) ;Y=[y(k) y(k-1) y(k-2)]';
    u(k)=u(k-1)+lq*[Yr-H*dU-G*Y];
    du(k)=u(k)-u(k-1);
end
figure(2);%第1个图形
plot(y,'-');
hold on
k=1:1:404;
plot(k,yr,':');
%title('y跟踪yr')
figure(3);%第2个图形
plot(u,'-');
%title('u的图像')