function [C,Kc,taui,taud] = simc_pidtune(sys, tauc)
% SIMC_PIDTUNE  Returns a cascade PID controller based on the SIMC PID tuning
% rules.
%
% SIMC_PIDTUNE(SYS) uses the "tight" rule such that TAUC is equal to the time
% delay in the process model.
%
% [C,KC,TAUI,TAUD] = SIMC_PIDTUNE(SYS,TAUC) returns cascade PID controller
% parameters together with a controller of the form:
%
%          taui*s+1
%   C = Kc -------- (taud*s+1)
%           taui*s
%
%
% [C,KC,TAUI] = SIMC_PIDTUNE(SYS,TAUC) returns cascade PI controller
% parameters together with a controller of the form:
%
%          taui*s+1
%   C = Kc --------
%           taui*s
%
% See also SIMC_REDUCE, PROCESS_MODEL.

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

% Assumes that the system is already in process model format.  Model reduction
% method is called only to extract process parameters, not really for model
% reduction.
try
  % Special handling for pure time delay models, where order(sys) == 0.
  [rsys,k,theta,tau1,tau2] = reduce_plant_model(sys,max(order(sys),1));
catch
 error(lasterror)
end

ni = nargin;
if (ni < 2),
  % Use "tight" rule if closed-loop time constant is not specified.
  tauc = theta;
end

% Closed-loop time constant cannot be zero.
if tauc == 0
  error('The desired closed-loop time constant should be a positive number.')
end

if ~isinf(tau1)
  if tau1~=0
    % Non-integrating process
    Kc = tau1/k/(tauc+theta);
  else
    % Pure time delay
    Kc = 1/k/(tauc+theta);
  end
else
  % Integrating process
  Kc = 1/k/(tauc+theta);
end
taui = min(tau1,4*(tauc+theta));

if ~isinf(tau2)
  % Second-order process
  taud = tau2;
else
  % Double integrating process
  taud = 4*(tauc+theta);
  Kc = Kc/taud;
end

% Construct cascade PID controller model.
if taui~=0
  % PID controller.
  C = zpk( tf(Kc*[taui*taud taui+taud 1], [taui 0]) );
else
  % Pure integral controller.
  C = zpk([], 0, Kc);
end

set(C, 'DisplayFormat', 't');
end
