function [rsys,k,theta,tau1,tau2] = simc_reduce(sys,order,tauc)
% SIMC_REDUCE  Returns a first- or second-order time delay model as
% a low-order approximation of a system model by applying the SIMC model
% reduction rules.  The plant model is assumed to have, if any, real
% LHP-poles (i.e., positive lag time constants) and real zeros.  Individual
% model parameters are also returned.
%
% RSYS = SIMC_REDUCE(SYS, ORDER) internally computes a closed-loop time
% constant TAUC (using the Optimization Toolbox).  Note that the resulting
% time delay THETA can be zero.
%
% RSYS = SIMC_REDUCE(SYS, ORDER, TAUC) generates a low-order model by
% applying the SIMC model reduction rules with a closed-loop time constant
% TAUC.  Note that TAUC has no effect on the reduced model unless the
% system SYS has positive numerator time constants to be removed before
% model reduction.
%
% [RSYS,K,THETA,TAU1,TAU2] = SIMC_REDUCE(SYS,2) returns a second-order time
% delay model, RSYS, as a low-order approximation of the plant model SYS.
% Individual model parameters are also returned.  The approximate model is
% of the form:
%                    k
%   RSYS = --------------------- exp(-theta*s)
%          (tau1*s+1) (tau2*s+1)
%
% [RSYS,K,THETA,TAU1] = REDUCE_PLANT_MODEL(SYS,1) returns a first-order
% time delay model, RSYS, as a low-order approximation of the plant model
% SYS.  Individual model parameters are also returned.  The approximate
% model is of the form:
%
%              k
%   RSYS = ---------- exp(-theta*s)
%          (tau1*s+1)
%
% See also SIMC_PIDTUNE, PROCESS_MODEL.

% Author: Bora Eryilmaz, The MathWorks, Inc.

ni = nargin;
if (ni < 3)
  % Find default closed-loop time constant such that tauc==theta.
  delay = totaldelay(sys);
  if delay ~= 0
    x0 = delay;
  else
    x0 = 1;
  end
  if isempty(ver('optim'))
    error('The Optimization Toolbox is required for estimating TAUC.')
  end
  opt = optimset('fmincon');
  opt = optimset(opt, 'Display', 'none');
  try
    tauc = fmincon( @(x) myfun(sys,order,x), x0, ...
		    [],[],[],[],0,Inf, ...
		    @(x) mycon(sys,order,x), opt );
  catch
 error(lasterror)
end



end

% Nested functions for optimization.
function f = myfun(sys,order,x)
% Cost function
  G = remove_positive_numerator_time_constants(sys, x);
  [Gtmp,k,theta] = reduce_plant_model(G,order);
  f = abs(theta-x);
end
function [c,ceq] = mycon(sys,order,x)
% Constraint function
  G = remove_positive_numerator_time_constants(sys, x);
  [Gtmp,k,theta] = reduce_plant_model(G,order);
  c(1) = -theta;  % theta >= 0;
  c(2) = theta-x; % tauc >= theta
  ceq = [];
end

% Model approximation with closed-time constant, tauc, specified or
% computed using optimization.
try
  G = remove_positive_numerator_time_constants(sys, tauc);
  [rsys,k,theta,tau1,tau2] = reduce_plant_model(G,order);
catch
 error(lasterror)
end

end
