% Copyright c 1998-2002 The Board of Trustees of the University of Illinois
% 		  All rights reserved.
% Developed by:	Large Scale Systems Research Laboratory
%               Professor Richard Braatz, Director
%               Department of Chemical Engineering
%		University of Illinois
%		http://brahms.scs.uiuc.edu
% 
% Permission hereby granted, free of charge, to any person obtaining a copy
% of this software and associated documentation files (the "Software"), to
% deal with the Software without restriction, including without limitation
% the rights to use, copy, modify, merge, publish, distribute, sublicense,
% and/or sell copies of the Software, and to permit persons to whom the 
% Software is furnished to do so, subject to the following conditions:
% 		1. Redistributions of source code must retain the above copyright
%		   notice, this list of conditions and the following disclaimers.
%		2. Redistributions in binary form must reproduce the above 
%		   copyright notice, this list of conditions and the following 
%		   disclaimers in the documentation and/or other materials 
%		   provided with the distribution.
%		3. Neither the names of Large Scale Research Systems Laboratory,
%		   University of Illinois, nor the names of its contributors may
%		   be used to endorse or promote products derived from this 
%		   Software without specific prior written permission.
%
% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
% OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL 
% THE CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR
% OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, 
% ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
% DEALINGS IN THE SOFTWARE.
%
%-------------------------------------------------------------------------------
%  Program:      Computer Simulation of an Adhesive Coating Operation 
%  Authors:      Michael Simkus, Jeffrey G. Virmond, and Richard D. Braatz 
%  Last modified:	 January 2, 2002
%
%  Description:	  This program simulates the operation of an adhesive coater.
%                 It enables the user to simulate either open loop or closed
%                 loop operations.  Both Eulerian and Lagrangian frames of
%                 reference are available for the open loop simulation.
%                 The objective of the control algorithm is to minimize the 
%                 effect of disturbances on the changes in coating thicknesses.
%
%  Note:	  Figure and page references are to Avery Final Report 
%		  by R. D. Braatz, M. L. Tyler, and M. Morari.  The process
%		  including industrial implementation are described in detail
%		  in the journal paper:
% 
%	R.D. Braatz, M.L. Tyler, M. Morari, F.R. Pranckh, and L. Sartor.
%	Identification and Cross-directional Control of Coating Processes,
%	AIChE Journal, 38:1329-1339, 1992.
%
%
%	Users of this code should cite the above paper.
%------------------------------------------------------------------------------- 
clear
clf                                  
format compact
time_in_past = 9;						%Length of the time interval plotted
number_of_time_steps = 15;			%Simulation time
n = 12;                       	%Number of actuators
r = 0 ;                          %setpoint
k_hat = 1; disp(['The process model gain k_hat = ' num2str(k_hat)])             
P_hat = k_hat*(0.25*eye(12)-.02*ones(12,12));	%Open loop process model
d = [0.187 -0.284 -0.381 -0.013 0.182 0.016 ...
     -0.027 0.072 0.041 -0.043 0.056 0.194]' ;
variance = 0.000046; 	         %Variance of white measurement noise
noise = randn(n,number_of_time_steps) * sqrt(variance);  
x = NaN*ones(n,time_in_past);    %Allocates memory to facilitate calculations
                                 
choice_one = menu('Which simulation would you like to run ?','Open Loop','Closed Loop');

if (choice_one ~= 1) & (choice_one ~= 2), 
  disp('running open loop simulation'),
  choice_one = 1; 
end

if (choice_one == 1),				%Open loop simulation
  choice_two = menu(['Which frame of reference for the open loop ' ...
		     				'simulation would you like to view?'], ...
		     				'Lagrangian (along with the flow)', ...
                     'Eulerian (from outside the flow)');%these are defined in any fluids text
  if (choice_two ~= 1) & (choice_two ~= 2), 
    disp('will use a Lagrangian viewpoint'), 
    choice_two = 1; 
  end
  if (choice_two == 2),				%Eulerian frame of reference
    M = moviein(number_of_time_steps);
    for counter = 1:number_of_time_steps,               
      v(:,counter) = noise(:,counter) + d;
      u(:,counter) = [0 0 0 0 0 0 0 0 0 0 0 0 ]' ;
      for counter2 = 1:time_in_past-1,
        x(:,counter2) = x(:,counter2+1);
      end								%Repositions coating thicknesses
      x(:,time_in_past) = P_hat * u(:,counter) + v(:,counter);   
      
      if counter == 1,				%Defines rectangular grid for 3D plot 
        [dummy_position,dummy_lane] = meshgrid(0:time_in_past-1,1:n);
      end
      surf(0:time_in_past-1,1:n,x)	%Generates 3D plot of coating thicknesses
      ylabel('sensor lane #'),xlabel('position'),zlabel('deviations in coating thickness')
      M(:,counter) = getframe;  	%Sets each column vector to a movie frame
    end
    decide = menu('Would you like to see the movie now?','Yes','No');
    if decide == 1,
      movie(M,1,0.4)             %Runs movie  
    end
  else   								%Lagrangian viewpoint  
    for counter = 1:time_in_past,
      v(:,counter) = noise(:,counter) + d;
      if counter == 1,
        x(:,counter) = v(:,counter);
      else       
        x(:,counter) = P_hat * u(:,counter-1) + v(:,counter);
      end   
        u(:,counter) = [0 0 0 0 0 0 0 0 0 0 0 0 ]' ;
     end
      figure(1)
      [dummy_time,dummy_lane] = meshgrid(0:time_in_past-1,1:n); %Defines rectangular grid for 3D plot
      surf(0:time_in_past-1,1:n,x)	%Generates 3D plot of coating thicknesses
      view(-20,20), title('Open loop simulation')
      xlabel('time'),ylabel('sensor lane #'),zlabel('deviations in coating thickness')
      print -dpsc figure1.ps
  end
else										%Run closed loop simulation
  gamma = input(['Please enter a value for gamma ' ...
                 '(the filter parameter, default=0.2835): ']);
  if isempty(gamma), gamma = 0.2835; end
  k = input(['Please enter the true value for the ' ...
             'process gain k (default is 1): ']);
  if isempty(k), k = 1; end 
  P = (k/k_hat)*P_hat;
  for counter = 1:time_in_past,
    if counter == 1,
      y_m(:,counter) = noise(:,counter) + d;
    else     
      y_m(:,counter) = P * u(:,counter-1) + noise(:,counter) + d;
    end
    if counter == 1,
      y(:,counter) = y_m(:,counter);
      delta_u(:,counter) = inv(P_hat'*P_hat)*P_hat'*(r - y(:,counter));
      u(:,counter) = delta_u(:,counter);
    else
      y(:,counter) = (1-gamma)*(y(:,counter-1) + ...
			    P_hat*delta_u(:,counter-1)) + gamma*y_m(:,counter);
      delta_u(:,counter) = inv(P_hat'*P_hat)*P_hat'*(r-y(:,counter));
      u(:,counter) = delta_u(:,counter) + u(:,counter-1);
    end
  end
  [z,x] = meshgrid(0:time_in_past-1,1:n); %Defines rectangular grid for 3D plot
  figure(2)
  surf(0:time_in_past-1,1:n,y_m)	%plot closed loop time response
  xlabel('time'),ylabel('sensor lane #'),zlabel('deviations in coating thickness')
  print -dpsc figure2.ps
  figure(3)
  surf(0:time_in_past-1,1:n,u)	%plot closed loop time response
  xlabel('time'),ylabel('sensor lane #'),zlabel('manipulated variable changes') 
  print -dpsc figure3.ps
end