%
%  set the maximum control effort
%
  max_control_effort = 2.0;
%  
%  enter the model of the plant transfer function
%
  num_Gp_model = 100;
  den_Gp_model = [1 2 10];
%
%  enter the true plant transfer function
%
  num_Gp = [100];
  den_Gp = [1 2 10];
%
%  now construct the model transfer function and display it
%
  Gp = tf(num_Gp_model,den_Gp_model);
  disp('Gp is : '); Gp  
%
%=========================================================
%  This is model matching stuff
%
%  get the closed loop transfer function
%
%  zero position error ITAE
%
   wo = 5;
   Go = tf(wo^2,[1 1.4*wo wo^2]);
%  Go = tf(wo^3,[1 1.75*wo 2.15*wo^2 wo^3]);
%  Go = tf(wo^4,[1 2.1*wo 3.4*wo^2 2.7*wo^3 wo^4]);
%
%  zero velocity error ITAE
%
%  Go = tf([3.2*wo wo^2],[1 3.2*wo wo^2]);
%  Go = tf([3.25*wo^2 wo^3],[1 1.75*wo 3.25*wo^2 wo^3]);
%  Go = tf([5.14*wo^3 wo^4],[1 2.41*wo 4.93*wo^2 5.14*wo^3 wo^4]);
%
%  deadbeat systems
%
%   Go = tf(wo^2, [1 1.82*wo wo^2]);
%  Go = tf(wo^3, [1 1.90*wo 2.20*wo^2 wo^3]);
%  Go = tf(wo^4, [1 2.20*wo 3.50*wo^2 2.80*wo^3 wo^4]);
%
%  determine the controller
%
  Gc = Go/(Gp*(1-Go));  
%
%============================================================
%
%  Now get rid of pole/zero cancellations in the controller
%  and put it in the correct form for Simulink
%
%  tol determine how close poles/zers need to be before we call them the
%  same
%
  tol = 1e-3;
  disp('Gc is :');
  Gc = minreal(Gc,tol)
  [num_Gc,den_Gc] = tfdata(Gc,'v') ;
%
%=============================================================
%
%  determine the prefilter, for now it is just 1
%
  num_Gpre = 1;
  den_Gpre = 1;
%
%==============================================================
%
%  set the step amplitude and final time
%
  Amp = 1.0; %cm
  Tf = 2.5;  % seconds
%
%  now simulate the system in Simulink
%
  sim('closedloop');
%
%  determine the estimated closed loop response
%
  t = linspace(m_time(1), m_time(end), 1000);
  y = lsim(Go,Amp*ones(1,length(t)),t);
%
%  determine the open loop response
%
  yo = lsim(Gp,Amp*ones(1,length(t)),t);
%
%  plot the results
%  
  figure;
  subplot(3,1,1); 
  plot(t,yo,'k-','Linewidth', 2);
  grid; xlabel('Time (sec)'); ylabel('y (cm)'); title('Open Loop Plant');
  subplot(3,1,2);
  plot(m_time,m_y,'b-', t,y,'r-','Linewidth',2); 
  grid; xlabel('Time (sec)'); 
  ylabel('y (cm)'); legend('Model Matching','Model');
  title(['Closedloop Model Matching, \omega_0 = ', num2str(wo,3) 'rad/sec']);
  subplot(3,1,3); hold on
  plot(m_time,m_u,'m-','Linewidth',2); 
  plot(m_time,-max_control_effort*ones(1,length(m_time)),'--','Linewidth',2);
  plot(m_time,max_control_effort*ones(1,length(m_time)),'--','Linewidth',2);
  title('Control Effort');
  hold off
  grid; xlabel('Time (sec)'); ylabel('u (cm)');
 % axis([0 0.1 -0.5 0.5]);
   