matlab_非线性时滞系统的仿真程序

%----------------------- function for p102main.m ------------------------- function f=p102fun(t,x) global e1 global e2 global u1 global u2 global norm1 global norm2 %--------------------------- initialize ----------------------------- f=zeros(60,1); bb=0.01; lemda1=1.5;lemda2=2; ru1=0.2;ru2=0.2; sigma1=0.1;sigma2=0.2; m1=5;m2=5; m5=0.2;m6=0.2; Mf=4; x1d=[0.5*sin(t)+0.5*sin(0.5*t),0.5*cos(t)+0.25*cos(0.5*t)]'; x2d=[sin(0.5*t)+0.5*sin(1.5*t),0.5*cos(0.5*t)+0.75*cos(1.5*t)]'; xp1=[x(1) x(2)]'; xp2=[x(3) x(4)]'; theta1=x(9:33); % why? theta2=x(34:58); % why? xtal=zeros(10000,1); % why? xtal1=zeros(10000,1); % why? x11=x(1); x12=x(2); % x21 = x(3); % x22 = x(4); %------------------------- error & swith function --------------------- e1=xp1-x1d; e2=xp2-x2d; es1=lemda1*e1(1)+e1(2); % 滤波误差s1; es2=lemda2*e2(1)+e2(2); % 滤波误差s2; s1=es1; s2=es2; %------------------------- NN base function ------------------------ gama1 = lemda1*e1(2)+0.125*cos(0.5*t)+0.5*sin(t); gama2 = lemda2*e2(2)+0.75*1.5*cos(1.5*t)+0.25*cos(0.5*t); z1=[x(1) x(2) x(3) x(4) s1 gama1]'; zh1=ones(25,1); %dh1=ones(5,1); sigma=0.8; % why? sum1=0; for i=1:25 for j=1:6 % z1的维数 zh1(i)=zh1(i)*exp(-(z1(j)-0.5)^2/1.581); % Gaussian函数中的宽度1.581 哪里来的？ end sum1=sum1+zh1(i); end % y1=x(9:33); hz1=zh1'*theta1/sum1; % hz1 hz2 这个量是什么？ %---------------------------- ----------------------------- epsilon1=x(5); epsilon2=x(6); zeta1=x(59); zeta2=x(60); % Nuss1=exp(zeta1^2)*cos((pi/2)*zeta1); p11=0; for i=1:10000 p11=t*(2-sin(x(1))*sin(x(1)))*0.004+p11; % what? end d11=0; for i=1:375 d11=0.1*x(1)^2*0.004+d11; % what? end delay1=4.0*0.05/(2*(1-0.6)*max(s1^2,0.1^2))*d11; p1=p11*0.02+delay1; % ？ % p1=1; Nuss1=exp(zeta1^2)*cos((pi/2)*zeta1); u1=Nuss1*(p1*s1+hz1+epsilon1*tanh(s1/ru1)); %控制律u1 z2=[x(1) x(2) x(3) x(4) s2 gama2 u1]'; % //i rules, j veriables; i=5; j=7;// zh2=ones(25,1); %dh1=ones(5,1); sum2=0; for i=1:25 for j=1:7 % z2的维数 zh2(i)=zh2(i)*exp(-(z2(j)-0.5)^2/1.581); end sum2=sum2+zh2(i); end % y2=x(34:58); hz2=zh2'*theta2/sum2; % hz1 hz2 这个量是什么？ % Nuss2=exp(zeta2^2)*cos((pi/2)*zeta2); p21=0; for i=1:10000 %？ p21=t*(3+sin(x(4)))*0.004/2+p21; % ？ end % p2=p21+0.1*0.001/(2*(1-0.6)*max(s2^2,0.1^2))*0.2*x(4)*sin(x(3))*2; d21=0; for i=1:375 %？ d21=0.2*x(4)*sin(x(3))*0.004+d21; end delay2=4.0*0.05/(2*(1-0.6)*max(s1^2,0.1^2))*d21; p2=p21*0.01+delay2; % ？ % p2=1; Nuss2=exp(zeta2^2)*cos((pi/2)*zeta2); u2=Nuss2*(p2*s2+hz2+epsilon2*tanh(s2/ru2)); %控制律u2 %------------------- dead zone ------------------------- %k1l=0.5;k1r=1.5;k2l=1.5;k2r=2.5;b1l=-0.5;b1r=0.5;b2l=-2.5;b2r=2; k1l=0.5;k1r=0.5; k2l=1.5;k2r=1.5; b1l=-0.5;b1r=0.5; b2l=-2.5;b2r=2; if u1>=b1r % 死区？ w1=k1r*(u1-b1r); % 死区输出？ else if u1<=b1l w1=k1l*(u1-b1l); else w1=0; end end if u2 >= b2r w2=k2r*(u2-b2r); else if u2<=b2l w2=k2l*(u2-b2l); else w2=0; end end %----------------------------object equation ----------------------------- j=1; f(1)=x(2); if t-1-0.5*sin(t)<=0 % 为什么要加这个限制条件？ f(2)=x(3)-0.3*sin(x(3))+w1*(2-(sin(x(1)))^2)+0.1*x11^2+sin(x(3)); f(3)=x(4); f(4)=x(1)^2*w1+w1^2*(x(4)^2+x(1)+0.5*cos(x(3)))+w2*(3+sin(x(2)))+0.2*x(12)*sin(x12)+0.5*sin(x(4)); % 0.2*x(12)*sin(x12) 为什么这么写？ N=(t-1-0.5*sin(t))/0.004; %计算时滞的步数？ for i=1:N xtal(i)=x11; % 是什么？ xtal1(i)=x12; end else f(2)=x(3)-0.3*sin(x(3))+w1*(2-(sin(x(1)))^2)+0.1*xtal(j)^2; f(3)=x(4); f(4)=x(1)^2*w1+w1^2*(x(4)^2+x(1)+0.5*cos(x(3)))+w2*(3+sin(x(2)))+0.2*xtal1(j)*sin(xtal1(j)); % 为什么这么写？ j=j+1; N=(t-1-0.5*sin(t))/0.004; %计算时滞的步数？ for i=375:10000 xtal(i)=x11; xtal1(i)=x12; end end f(5)=m1*(s1*tanh(s1/ru1)-sigma1*x(5)); % ？(30)？ f(6)=m2*(s2*tanh(s2/ru2)-sigma2*x(6)); f(7)=0; f(8)=0; %---------------------------- 范数 ---------------------------- norm1=0; for i=9:33 norm1=norm1+x(i)^2; end norm1=sqrt(norm1); norm2=0; for i=34:58 norm2=norm2+x(i)^2; end norm2=sqrt(norm2); dh1=zeros(25,1); dh2=zeros(25,1); dh1(1:25)=zh1/sum1; dh2(1:25)=zh2/sum2; % ----------------------- adaptive law ----------------------------- % if norm1<Mf | (norm1==Mf & s1*theta1'*zh1<0) f(9:33)= m5*(s1*dh1-sigma1*x(9:33)); else f(9:33)= m5*(s1*dh1-sigma1*x(9:33)-s1*theta1*theta1'*zh1/norm1) end if norm2<Mf | (norm1==Mf & s1*theta2'*zh2<0) f(34:58)= m6*(s2*dh2-sigma1*x(34:58)); else f(34:58)= m6*(s2*dh2-sigma2*x(34:58)-s2*theta2*theta2'*zh2/norm2) end f(59)=p1*s1^2+hz1*s1+epsilon1*s1*tanh(s1/ru1); % 文章(26)式？ f(60)=p2*s2^2+hz2*s2+epsilon2*s2*tanh(s2/ru2); %~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ end ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~