function [oilfrtf,doilfrtf,doilfrss,oilfrmod,oilfrste]=mkoilfr(Ts,keepout,trunc,perturb,fname);
%MKOILFR Defines Shell heavy oil fractionator model
%
% [oilfrtf,doilfrtf,doilfrss,oilfrmod,oilfrste]=...
%                           mkoilfr(Ts,keepout,trunc,perturb,fname);
%
% This function defines the 'Shell' heavy oil fractionator model.
% This is used as an extended example in the book:
% 'PREDICTIVE CONTROL WITH CONSTRAINTS' by J.M.Maciejowski
%
% Input arguments:
%      Ts: Sampling time (minutes). Default: 4
% keepout: Vector of outputs to retain. Default: [1,2,7].
%          Input 3 (Bottoms reflux duty) is also the last output.
%   trunc: Time at which step response model is truncated. Default: 250 mins.
% perturb: 5x1 vector of perturbations. Scalar => all elements equal.
%          Each element should lie in [-1,1]. Default: zeros(5,1).
%   fname: String defining name of .mat file for storing results. Default:none.
%
% Output arguments (each represents a 5-input, 7-output system):
%   The first 3 are LTI objects. See 'notes' property for some details.
%  oilfrtf: continuous-time, transfer-function form, with time delays.
% doilfrtf: discrete-time, transfer-function form.
% doilfrss: disrete-time, state-space form.
% oilfrmod: discrete-time, state-space, MPC Toolbox MOD format.
% oilfrste: discrete-time, step-response, MPC Toolbox STEP format.
%
% The model was defined in: 
% 'The Shell Process Control Workshop. Process Control Research:
%  Industrial and Academic Perspectives', by
% D.M.Prett and M.Morari (eds), (Butterworth:Stoneham,MA), 1987.
%
% This file uses MPC Toolbx function SS2MOD.

% File written by J.M.Maciejowski, 11 December 1999.
% Copyright (C) J.M.Maciejowski, 1999. All Rights Reserved. 
%
% Modifications: 15.12.99, 31.12.99, 3.1.2000, 6.1.00, 8.1.00.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%% PROCESS INPUT ARGUMENTS:

if nargin==0,
  disp('USAGE:')
  disp(['[oilfrtf,doilfrtf,doilfrss,oilfrmod,oilfrste]=',...
        'mkoilfr(Ts,keepout,trunc,perturb,fname)'])
  return
end
nargchk(1,5,nargin);
if nargin < 2, trunc = 250; end % MINUTES
if nargin < 3, keepout = [1,2,7]; end % top & side end points, bottoms temp.
if nargin < 4, 
  perturb = zeros(5,1);
elseif max(size(perturb))==1,
  perturb = perturb*ones(5,1);
elseif min(size(perturb))>1,
  error('Input argument PERTURB should be scalar or 5x1 vector.')
elseif length(perturb)~=5,
  error('Input argument PERTURB should be scalar or 5x1 vector.')
end
if any(abs(perturb)>1),
  error('Each element of PERTURB must be in the range -1 to +1.')
end
perturb=perturb(:);  % Ensure it's a column vector.
if nargin < 5 | isempty(fname), 
  fname=[];   % Don't save to .mat file
elseif ~ischar(fname),
  error('Input argument FNAME must be a string.')
end

if ~isempty(fname),
  if exist([fname,'.mat'])==2,
    disp(['Warning: File ',fname,'.mat already exists.'])
    yesno = input(['Type ''y'' to continue.',... 
            ' (This will overwrite the existing file.)'],'s');
    if ~strncmp(yesno,'y',1),
      disp('Not creating new oil fractionator model.')
      return
    end
  end
end

%%%%%%%%% CHECK OUTPUT ARGUMENTS:

nargchk(1,5,nargout);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%)

% DEFINE CONTINUOUS-TIME MODEL AS LTI OBJECT:

% Define individual transfer functions first, without time delays,
% and with perturbed gains. Create all outputs first, then delete 
% unwanted ones. (It's inefficient but simple.)

% Output 1:
g11=tf(4.05+2.11*perturb(1),[50,1]);
g12=tf(1.77+0.39*perturb(2),[60,1]);
g13=tf(5.88+0.59*perturb(3),[50,1]);
g14=tf(1.20+0.12*perturb(4),[45,1]);
g15=tf(1.44+0.16*perturb(5),[40,1]);

% Output 2:
g21=tf(5.39+3.29*perturb(1),[50,1]);
g22=tf(5.72+0.57*perturb(2),[60,1]);
g23=tf(6.90+0.89*perturb(3),[40,1]);
g24=tf(1.52+0.13*perturb(4),[25,1]);
g25=tf(1.83+0.13*perturb(5),[20,1]);

% Output 3:
g31=tf(3.66+2.29*perturb(1),[9,1]);
g32=tf(1.65+0.35*perturb(2),[30,1]);
g33=tf(5.53+0.67*perturb(3),[40,1]);
g34=tf(1.16+0.08*perturb(4),[11,1]);
g35=tf(1.27+0.08*perturb(5),[6,1]);

% Output 4:
g41=tf(5.92+2.34*perturb(1),[12,1]);
g42=tf(2.54+0.24*perturb(2),[27,1]);
g43=tf(8.10+0.32*perturb(3),[20,1]);
g44=tf(1.73+0.02*perturb(4),[5,1]);
g45=tf(1.79+0.04*perturb(5),[19,1]);

% Output 5:
g51=tf(4.13+1.71*perturb(1),[8,1]);
g52=tf(2.38+0.93*perturb(2),[19,1]);
g53=tf(6.23+0.30*perturb(3),[10,1]);
g54=tf(1.31+0.03*perturb(4),[2,1]);
g55=tf(1.26+0.02*perturb(5),[22,1]);

% Output 6:
g61=tf(4.06+2.39*perturb(1),[13,1]);
g62=tf(4.18+0.35*perturb(2),[33,1]);
g63=tf(6.53+0.72*perturb(3),[9,1]);
g64=tf(1.19+0.08*perturb(4),[19,1]);
g65=tf(1.17+0.01*perturb(5),[24,1]);

% Output 7:
g71=tf(4.38+3.11*perturb(1),[33,1]);
g72=tf(4.42+0.73*perturb(2),[44,1]);
g73=tf(7.20+1.33*perturb(3),[19,1]);
g74=tf(1.14+0.18*perturb(4),[27,1]);
g75=tf(1.26+0.18*perturb(5),[32,1]);

% Now put them together into a 7x5 multivariable model:
oilfrtf = tf([g11,g12,g13,g14,g15 ; g21,g22,g23,g24,g25;...
            g31,g32,g33,g34,g35 ; g41,g42,g43,g44,g45;...
            g51,g52,g53,g54,g55 ; g61,g62,g63,g64,g65;...
            g71,g72,g73,g74,g75]);

% Now add the time delays:
set(oilfrtf,'ioDelayMatrix',[27,28,27,27,27; 18,14,15,15,15;...
                             2,20,2,0,0; 11,12,2,0,0;...
                             5,7,2,0,0; 8,4,1,0,0; 20,22,0,0,0]);

% Input and output names:
set(oilfrtf,'InputName',{'Top draw';'Side draw';'Bottoms duty';...
			 'Intermediate duty';'Upper duty'});
set(oilfrtf,'OutputName',{'Top end point';'Side end point';...
                          'Top temperature';'Upper reflux temp';...
                          'Side draw temp';'Inter reflux temp';...
                          'Bottoms reflux temp'});

% Create some Notes about the model:
dnum = now;
daystring = datestr(dnum,1);
timestring = datestr(dnum,13);
set(oilfrtf,'Notes',{'Shell heavy oil fractionator model';...
    ['   created by MKOILFR on ',daystring,' at ',timestring,'.'];...
    'with the following input gain perturbations:';...
    num2str(perturb');...
    'Inputs 1,2 and 3 are manipulated variables.';...
    'Input 4 is a measured disturbance.';...
    'Input 5 is an unmeasured disturbance.';...
    'All outputs are measured.';...
    'Time is expressed in minutes.'});

disp([' '])
disp('Continuous-time model created as LTI TF object')
if norm(perturb)==0,
  disp('with nominal gains.')
else
  disp('with the following multiplicative perturbations of the input gains:')
  disp(perturb')
end

% Keep the required outputs only:
oilfrtf = oilfrtf(keepout,:);
disp('The following outputs have been retained:')
disp(keepout);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% NOW DEFINE DISCRETE-TIME MODEL:

if nargout>1,
  doilfrtf = c2d(oilfrtf,Ts);
  doilfrtf = delay2z(doilfrtf); % Put all delays into transfer functions

  % Bring input 3 out as last output (with 1-step delay):
  newtf = tf(1,[1 0],Ts);                 % 1/z
  doilfrtf = [doilfrtf; [0,0,newtf,0,0]]; % Add new output

  % Update notes etc:
  outnames = get(oilfrtf,'OutputName');
  outnames = [outnames;{'Bottoms Duty'}];
  set(doilfrtf,'OutputName',outnames);
  set(doilfrtf,'Notes',get(oilfrtf,'Notes'));

  disp('Discrete-time model created as LTI TF object.')
  disp(['Sample time: ',num2str(Ts),' minutes.'])
  disp('Bottoms Duty (input 3) brought out as additional (last) output')
  disp('     with one-step delay.')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% CREATE MODEL IN MPC TOOLBOX FORMAT:


if nargout>2,
  disp([' '])
  disp('Converting model to state-space form...')
  doilfrss = ss(doilfrtf); % Convert model to state-space form
  doilfrss = minreal(doilfrss); % Remove any nonminimal modes
  a = get(doilfrss,'a'); b = get(doilfrss,'b');
  c = get(doilfrss,'c'); d = get(doilfrss,'d');
  nstates=size(a,1);
  set(doilfrss,'Notes',get(doilfrtf,'Notes'));
  disp(['Number of states in model: ',int2str(nstates)]);
  disp('Discrete-time model created as LTI SS object.')
end

% Create model in 'MOD' format:
% Inputs 1,2 and 3 are manipulated variables.
% Input 4 is a measured disturbance.
% Input 5 is an unmeasured disturbance.
% All the outputs are measured.
if nargout>3,  
  disp([' '])
  disp('Converting model to MPC Toolbox format...')
  infvec = [Ts,nstates,3,1,1,length(keepout)+1,0]; % Information vector for SS2MOD
  oilfrmod = ss2mod(a,b,c,d,infvec);  % MPC Toolbox function.
  disp('Discrete-time model created in MOD format.')
end

% Create model in 'STEP' format: 
if nargout>4,
  oilfrste = ss2step(a,b,c,d,trunc,Ts);
  disp('Discrete-time model created in STEP format,')
  disp(['     step response truncated after ',int2str(trunc),' minutes.'])
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if ~isempty(fname),
  if nargout==1, save(fname,'oilfrtf');
  elseif nargout==2, save(fname,'oilfrtf','doilfrtf');
  elseif nargout==3, save(fname,'oilfrtf','doilfrtf','doilfrss');
  elseif nargout==4, 
    save(fname,'oilfrtf','doilfrtf','doilfrss','oilfrmod'); 
  elseif nargout==5, 
    save(fname,'oilfrtf','doilfrtf','doilfrss','oilfrmod','oilfrste'); 
  end
  disp([' '])
  disp(['File ',fname,'.mat has been written to disk.'])
end

%%%%%%%%%%%%%%%%   End of MKOILFR   %%%%%%%%%%%%%%%%