function [rsys,k,theta,tau1,tau2] = reduce_plant_model(sys, order)
% REDUCE_PLANT_MODEL  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 RHP-zeros (i.e.,
% negative numerator time constants).  Individual model parameters are also
% returned.
%
% [RSYS,K,THETA,TAU1,TAU2] = REDUCE_PLANT_MODEL(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 REMOVE_POSITIVE_NUMERATOR_TIME_CONSTANTS, PROCESS_MODEL.

% Author(s): Bora Eryilmaz, The MathWorks, Inc.

% Check the plant model for SIMC model reduction compatibility.
check_model(sys);

% Check for positive, real numerator time constants.
[Z,P,K] = zpkdata(sys, 'v');
if any(real(Z)<0 & ~imag(Z))
  error('Positive numerator time constants should be removed from the model.')
end

% Inverse response (negative numerator) time constants in descreasing order.
Zrealpos = Z(real(Z)>0 & ~imag(Z));
Tj0inv = sort(+1./Zrealpos, 'descend');

% Lag time constants in descreasing order.
Prealneg = P(real(P)<0 & ~imag(P));
taui0 = sort(-1./Prealneg, 'descend');

% Default second-order plus time delay model parameters.
k = K * prod(-Z(Z~=0)) / prod(-P(P~=0));
theta = totaldelay(sys) + sum(Tj0inv);
tau1 = 0;
tau2 = 0;

% Augment lag times with up to 2 integrators for integrating systems.
nInt = min(order,length(find(P==0)));
taui0 = [Inf(nInt,1); taui0];

% Fold removed lag times into tau1, tau2, and theta.
if numel(taui0) >= 1
  tau1 = taui0(1);
end
if order == 1
  if numel(taui0) >= 2
    tau1 = tau1 + taui0(2)/2;
    theta = theta + taui0(2)/2 + sum(taui0(3:end));
  end
elseif order == 2
  if numel(taui0) >= 2
    tau2 = taui0(2);
    if numel(taui0) >= 3
      tau2 = tau2 + taui0(3)/2;
      theta = theta + taui0(3)/2 + sum(taui0(4:end));
    end
  end
else
  error('Desired model order should be 1 or 2.')
end

% Reorder time constants so that tau1 >= tau2.
if tau2 > tau1
  temp = tau1; tau1 = tau2; tau2 = temp;
end

% Construct process model.
rsys = process_model(k,theta,tau1,tau2);
end