There are two files below, save the first one as y2.m and the other as y2p.m and run y2.m in matlab. Please let me know if you have any questions. You might also want to look at the following related papers %[1] Q. Zou and S. Devasia "Preview-based Stable-Inversion for Output Tracking," % ASME J. of Dynamic Systems, Measurement and Control, % Vol. 121 (4), pp. 625-630, December 1999. %[2] S. Devasia "Approximated Stable Inversion for Nonlinear Systems with % Nonhyperbolic Internal Dynamics," IEEE Trans. on Automatic Control, % Vol. 44 (7), pp. 1419-1425, July 99. %[3] J.S. Dewey, K. K. Leang and S. Devasia "Experimental and Theoretical Results % in Output-Trajectory Redesign for Flexible Structures," ASME Journal of % Dynamic Systems, Measurement, and Control, Vol. 120 (4), pp. 456-461, Dec. 1998. %[4] S. Devasia and B. Paden "Stable Inversion for Nonlinear Nonminimum-Phase % Time-Varying Systems," IEEE Trans. on Automatic Control, % Vol. 43 (2), pp. 283-288, Feb. 1998. %[5] S. Devasia, B. Paden and C. Rossi "Exact-Output Tracking Theory for Systems % with Parameter Jumps," International Journal of Control, % Vol. 67 (1), pp. 117-131, May 1997. %[6] S. Devasia, D. Chen and B. Paden "Nonlinear Inversion-Based Output Tracking," % IEEE Transactions on Automatic Control, Vol. 41 (7), pp. 930-942, July 1996. %[7] D. Croft and S. Devasia "Vibration Compensation for High Speed Scanning % Tunneling Microscopy," Review of Scientific Instruments published by % the American Institute of Physics, Vol. 70 (12), pp. 4600-4605, December 1999. %[8] D. Croft, G. Shedd and S. Devasia “Creep, Hysteresis, and Vibration Compensation for % Piezoactuators: Atomic Force Microscopy Application,” ASME Journal of Dynamic Systems, % Measurement and Control, Vol. 123 (35), pp. 35-43, March 2001. and a recent article on should inversion be used in the presence of uncertainty % [9] S. Devasia “Should Model-based Inverse Inputs be used as Feedforward under Plant Uncertainty?” IEEE Trans. on Automatic Control, Vol. 47(11), pp. 1865-1871, Nov 2002 Regards Santosh Devasia %file 1 save as y2.m ******************* % This should be saved as y2.m %%%%% clear clf %% kgk=.1; nsteps=5; % del=0.005; % integration step size t1=-4:del:-del; tm=0:del:2*pi; t2=((max(tm)):del:10); T_time=[t1 tm t2]; ym=-cos(tm); ym=kgk*(ones(size(ym))+ym); ym=[zeros(size(t1)) ym zeros(size(t2))]; % desired output profile T_time = -4:del:max(T_time)+del; %figure(1); clf; plot(T_time,ym) % % Initializing g_x1=[]; g_x2=[]; % the integration stuff % [mm,nn]=size(ym); x1=zeros(size(ym)); x2=zeros(size(ym)); f1=zeros(size(ym)); f2=zeros(size(ym)); % ***************************************** % Starting the loop for the zero dynamics for jk=1:nsteps, 1+nsteps-jk for jj=1:nn, f1(jj)=-ym(nn+1-jj); f2(jj)=x1(jj)*x1(jj); end [x1,yx1]=lsim(-1,1,1,0,f1',(T_time+4)); [x2,yx2]=lsim(-1,1,1,0,f2',(T_time+4)); for j=1:nn, x1(j)=yx1(nn+1-j); end g_x1=[g_x1 ; x1' ]; g_x2=[g_x2 ; x2' ]; %figure(2); clf ; %subplot(221), plot(T_time,g_x1') %subplot(222), plot(T_time,g_x2') end z1=x1; z2=x2; yd=+kgk*sin(tm); yd=[zeros(size(t1)) yd zeros(size(t2))]; ydd=+kgk*cos(tm); ydd=[zeros(size(t1)) ydd zeros(size(t2))]; u1=zeros(size(yd)); x1=zeros(size(yd)); x2=zeros(size(yd)); x3=zeros(size(yd)); x4=zeros(size(yd)); % % for j=1:nn, x1(1,j)=ym(j)+3*z1(j); x2(1,j)=4*ym(j)+yd(j)+6*z1(j); x3(1,j)=z1(j); x4(1,j)=z2(j); al=10*x1(j) -7*x2(j) -12*x3(j) +x1(j)*x1(j)*x1(j); be=2 + sin(x4(j))*sin(x4(j)); u1(j)=(1/be)*(ydd(j)-al); end % figure(1); clf plot(T_time,u1); xlabel('time') ylabel('Inverse Input') % %%% the integration routine without truncation save ss1 T_time u1 YY0=[0 ; 0; 0; 0]; [tt,yy]=ode23('y2p',[-3.9 9.99], YY0); figure(2); clf subplot(211), plot(tt,yy) xlabel('time') ylabel('States') yd=yy(:,1) - 3*yy(:,3); subplot(212), plot(tt,yd,'-.',T_time,ym,'-'); xlabel('time') ylabel('Output without truncation') %%% the integration routine with truncation %save dsdata_4 % results of truncation yy0=[0 ; 0; 0; 0]; [tt,yy]=ode23('y2p',[0, 9.99],yy0); figure(3); clf subplot(211), plot(tt,yy) xlabel('time') ylabel('States') yd=yy(:,1) - 3*yy(:,3); subplot(212), plot(tt,yd,'-.',T_time,ym,'-'); xlabel('time') ylabel('Output with truncation') %save dsdata_0 %%% %%%%%%%%%%%%%%%%%%*********************************** % save this as y2p.m ************************ %%%%%%%%%%%%%%%%%%*********************************** % save this as y2p.m function yp = fun(t,y); load ss1 yp = zeros(4,1); u=interp1(T_time',u1',t); yp(1)=-y(1) + y(2); yp(2)=-3*y(2) + y(1)*y(1)*y(1) +(2 +(sin(y(4)))^2)*u; yp(3)= y(1) - 2*y(3); yp(4)= -y(4) +y(3)*y(3);