
%===============================================================================
% $Id: ifacconf.tex 19 2011-10-27 09:32:13Z jpuente $  
% Template for IFAC meeting papers
% Copyright (c) 2007-2008 International Federation of Automatic Control
%===============================================================================
\documentclass{ifacconf}

\usepackage{graphicx}      % include this line if your document contains figures
\usepackage{natbib}        % required for bibliography
\usepackage{booktabs}
\usepackage{afterpage}
\usepackage[]{units}
\usepackage{float}
\usepackage{amsmath}
%\usepackage{subfig} 
%\usepackage{hyperref}
\def\beq{\begin{equation}} 
\def\eeq{\end{equation}} 
\def\bb{\begin{itemize}}  % the b stands for bullet
\def\eb{\end{itemize}}
\def\ben{\begin{enumerate}}  %start each number with \item 
\def\een{\end{enumerate}} 
%\def\beqar{\begin{array}} 
%\def\eeqar{\end{array}} 
%\def\beqa{\begin{eqnarray}} 
%\def\eeqa{\end{eqnarray}} 
\def\({\left(} 
\def\){\right)} 
%\def\[{\left[} 
%\def\]{\right]} 
\def\/{\over} 
%\def\p{\partial} 
%\def\pd #1/#2/{\frac{\p #1}{\p #2}} 
%\def\pds #1/#2/#3/{\({\p #1 \over \p #2}\) _{#3}} 
%\def\fr{\frac} 
%\def\*{^{\ast}} 
%\def\beqn{$$} 
%\def\eeqn{$$} 
%\def\eps{\varepsilon} 
\def\q{\quad} 
%\def\qq{\qquad} 
\def\bi{\begin{description}} %remember to write \item[case1] 
\def\ei{\end{description}} 
%\def\bit{\begin{itemize}}    %start each number with \item 
%\def\eit{\end{itemize}} 
%\def\cl{\centerline}
\def\noi{\noindent}
%\def\bc{\begin{center}}
%\def\ec{\end{center}}
%\def\p{\partial}
\def\w{\omega}
\def\x{\par\noi}
\def\xx{\medskip\par\noi}
%\def\sn{\underline\sigma}
\def\eqdef{\buildrel \rm def \over =} 
\newcommand{\twltt}{\tt}

%===============================================================================
\begin{document}
\begin{frontmatter}

\title{Optimal PI-Control and Verification of the SIMC Tuning Rule} 
% Title, preferably not more than 10 words.

%\thanks[footnoteinfo]{Sponsor and financial support acknowledgment
%goes here. Paper titles should be written in uppercase and lowercase
%letters, not all uppercase.}

\author{Chriss Grimholt}  ~~~~\author{Sigurd Skogestad$^*$} 


\address{Department of Chemical Engineering, Norwegian University of Science and Technology (NTNU),    Trondheim, Norway\\
 $^*$e-mail: skoge@ntnu.no}



\begin{abstract}                % Abstract of not more than 250 words.
Optimal PI-settings are derived for first-order with delay processes for specified levels of robustness
($M_s$-value) and compared with the simple SIMC-rule.
Optimality (performance) is defined in terms of the integrated absolute error (IAE) of the output   
for combined step changes in setpoints and input disturbances.
With SIMC, the robustness level is adjusted by changing the tuning parameter $\tau_c$, and 
the SIMC-rule was found to give surprisingly good setting with almost Pareto-optimal
performance. 
The exception is a pure time delay processes where the SIMC-rule gives a pure integral controller with
somewhat sluggish response. 
A simple modification to improve on this, 
is to increase the time constant in the rule by one third of the time delay.
%The SIMC-rules have 

\end{abstract}

\begin{keyword}
Optimality, PI controllers, Pareto-optimal, robustness, performance evaluation, and closed-loop model approximation.
\end{keyword}

\end{frontmatter}
%===============================================================================

\section{Introduction}


 In this paper, we consider a PI-controller 
\beq  c(s) = K_c \cdot \( 1 + { 1 \over \tau_I s}\) 
\label{pidseries} \eeq
where $K_c$ is the controller gain and $\tau_I$ the integral time, applied to a first order with time delay process
\beq  g(s) = { k \over \tau_1 s + 1} \cdot e^{-\theta s} 
\label{modeltrasferfunction} \eeq
where $k$ is the process gain, $\tau_1$ is the process time constant and $\theta$ is the time delay.

\begin{figure}[b]
	\input{blockDiag}
	\caption{Block diagram of feedback control system. We may treat an output disturbance ($d_{ys}$) as a special case of setpoint 			change ($y_s$) }
			%{\footnotesize In this chapter we consider an input (``load'') disturbance ($g_d=g$). }
	\label{1dofblockdiagram}
\end{figure}


The SIMC method for PID controller tuning \citep{skogestad:03} has already found wide industrial usage. The SIMC-rules are analytically derived, and from a first or second order process model we can easily find the PI- and PID-controller setting, respectively. For a first-order process, the SIMC PI-settings are
\beq K_c = {1 \over k} {\tau_1 \over \tau_c + \theta}    \label{kcsispi} \eeq
\beq \tau_I = % \min\{\tau_1, {4 \over k'\ K_c}\}  =
            \min\{\tau_1, 4 (\tau_c + \theta)\}  
\label{tauisispi} \eeq
The SIMC PI-rule has one tuning parameter $\tau_c$ which can be used to trade off between  performance (``tight" control) and robustness (``smooth" control). 
%Although the rule is generally considered to give good trade-off (not only between performance and robustness, but also between set-point and and disturbance response), there have been reported for some lag dominated processes that the rules results sluggish disturbance rejection \citep{haugen:10}. 
The rule was originally derived mainly with simplicity in mind, and although it was mentioned in the original paper that these are "probably the best simple PID tuning in the world", one may want to ask the questions: How good is the SIMC PI-rule? Can it be improved? Is the recommended value $\tau_c=\theta$ a good default value for the tuning parameter $\tau_c$?

To study this, we need to compare the SIMC PI performance with  the ``optimal'' PI-controller for a process on the from \eqref{modeltrasferfunction}.
However, ``optimality" is generally difficult to define as there are many issues to consider, including:
\begin{itemize}
	\item Output performance
	\item Robustness
	\item Input usage
	\item Noise sensitivity
\end{itemize}


This may be considered a multiobjective optimization problem, but fortunately the trade-off space has only one
main dimension, namely high versus low controller gain.
High controller gain favours good output performance, whereas low controller gain favours the three other objectives listed above. 
We can then simplify and say that there are two main objectives:
\begin{enumerate}
	\item Output performance
	\item Robustness, input usage and noise sensitivity
\end{enumerate}
The idea is then to compare the SIMC controller, with $\tau_c$ as a parameter, with the ``Pareto-optimal" PI-controller.
Pareto-optimality applies to multiobjective problems, and means that no further improvement can be made
in objective 1 without sacrificing objective 2. We still need to define more exactly objectives 1 and 2. 

Towards the end the paper, we consider how to find a first-order with delay model, which is required for the SIMC tuning rule.
It can be difficult to obtain open-loop responses for model approximation and \cite{shams:10} presented a new method for using closed-loop data for tuning.
We will modify this method to obtain a first-order with time delay model. 





\section{Evaluation of performance and robustness}

\subsection{Output performance}
Output performance (objective 1) is related to the difference between the measurement $y(t)$ and its setpoint $y_s$ (Figure \ref{1dofblockdiagram}), and may be quantified in many different ways. 
For example, for a setpoint change, we might consider the rise time, overshoot and settling time. However, to quantify the output performance in terms of a single scalar, 
we choose to use the integrated absolute error:
$$ {\rm IAE} = \int_0^{\infty} |y(t)-y_s(t)| dt $$
The IAE-value depends on which disturbance or setpoint change we consider. We choose to consider a weighted average of IAE for a step change in the load disturbance $d$ (IAE$_d$)  and setpoint $y_s$ (IAE$_{y_s}$):
\beq J(c) = 0.5 \left[
{{\rm IAE}_{ys}(c) \over {\rm IAE}^\circ_{ys}} +  
{{\rm IAE}_{d}(c) \over {\rm IAE}^\circ_{d}} \right] \label{eq:jpo} \eeq
Importantly, the weighted cost $J$ is independent of the process gain $k$, the disturbance and setpoint magnitudes, 
and of the unit used for time.

The weighting  factors IAE$^\circ_{ys}$ and IAE$^\circ_{d}$ 
are for reference PI-controllers, which for the given process are IAE-optimal for a setpoint change and a disturbance, respectively. 
To get robust reference controllers, they are required to have $M_s=1.59$. 
Values for  IAE$^\circ_{ys}$ and IAE$^\circ_{d}$,  are given for four different processes in Table~\ref{tabref0}.
We could have used the truly optimal IAE-value as the reference (without the restriction \mbox{$M_s=1.59$}), 
but this would not have changed the results much because the IAE-value is anyway quite close to its minimum at $M_s=1.59$.
Note that two different optimal PI-controllers are used to obtain the two reference values, whereas  a single controller
$c$ is used to find IAE$_{ys}(c)$ and IAE$_d(c)$ when evaluating the weighted IAE-cost $J(c)$.

It may be argued that a two-degree of freedom controller with a setpoint filter may be used to improve
the response for setpoints, but note that a setpoint change is equivalent to an output disturbance 
($d_{ys}$ in Figure~\ref{1dofblockdiagram}) which is not affected by the filter.
Thus, both setpoint changes (output disturbances) and input disturbances should be included when evaluating performance, and
to get a good balance between the two, we weigh them about equally by using the cost in \eqref{eq:jpo}.

%The setpoint performance is often referred to as the ``servo" behavior
%and the disturbance (in this case the input ``load" disturbance) performance is often referred to as
%``regulator" behavior.  $$ {\rm IAE} = \int_0^{\infty} |y(t)-y_s(t)| dt $$


%Fast speed of response and good disturbance rejection  
%(favored by a small value of $\tau_c$). This ``tightness'' can be quantified by the magnitude of the setpoint error,
%$|y(t) - y_s(t)|$, which should be as small as possible.  
%Here, one may consider different ``norms'' of the error, for example, the maximum deviation 
%($\infty$-norm),
%the integrated square deviation (2-norm) and the integrated absolute error (IAE) (1-norm), 
%$$ {\rm IAE} = \int_0^{\infty} |y(t)-y_s(t)| dt $$

\subsection{Robustness, input usage and noise sensitivity}

There are many ways to quantify objective 2 related to the benefits of low gain. 
Robustness my be quantified in terms of sensitivity peak ($M_s$), complementarity  sensitivity peak ($M_t$), gain margin (GM), phase margin (PM), allowable time delay error ($\nicefrac{\Delta \theta}{ \theta}$), etc.
Input usage may also be quantified in many ways, but we have found total variation (TV) to be a good measure:
\beq {\rm TV} = \int_0^{\infty} |\frac{du}{dt}|dt=\sum_0^\infty|u_i-u_{i-1}| \eeq
The noise sensitivity on the outputs ($y$) may be quantified by the complimentary sensitivity peak ($M_t$). However, even more important is usually the noise sensitivity on the inputs ($u$), 
which for a PI-controller is closely related to the magnitude of the controller gain, $K_c$, which should be low to have small noise sensitivity. 
This issue is not considered in this paper, partly because the noise frequency and magnitude shows large variation from case to case. 

In this paper, we choose to quantify all the issues related to objective 2 in terms of the sensitivity peak, $M_s$:
$$M_s = \max_\omega \left| \frac{1}{1+gc(j\omega)} \right|$$ 
which is the peak of the sensitivity function in the frequency domain. In robustness terms, $M_s$ is the inverse of the closest distance between the critical point -1 and the loop transfer function $gc$ in the Nyquist plot. 
For robustness, a small value of $M_s$ is desired, and generally $M_s$ should not exceed 2. A typical ``good'' value is about 1.6, and notice that $M_s <1.6$ guarantees GM $>2.67$ and PM $>36.4^o$  \citep{rivera:86}.
As shown later, the $M_s$-value is also closely related to TV, at least for process and controllers used in this paper.
%The usage of the manipulated variable can be quantified by the total variation (TV),


%
%
%where low input usage indicates smooth and robust control, and high input usage indicates more aggressive and less robust control. We do not use total variation for evaluation of the controller settings in this paper because, as we will show later, it was found to be very correlated with the $M_s$.



%One reason is that we found that the $M_s$-value correlates well with the input usage as given
%by its total variation (TV), which agrees with the findings of \cite{foley:05}.
%Such a correlation is reasonable since a large $M_s$-value corresponds to an oscillatory system with large input
%variations.  

%Good robustness, small input changes and small noise sensitivity 
%(favored by a large value of $\tau_c$). 
%The ``smoothness'' is here quantified by the peak value $M_s\ge 1$ of the frequency-dependent sensitivity function, %$S=1/(1+gc)$.
%In terms of robustness, $1/M_s$ is the closest distance of the loop transfer function $gc$ to the critical $(-1)$-point in the %Nyquist diagram, so 
%$M_s$ should be as small as possible. 
%The value of $M_s$  varies between is 1 (``slow control'' with $\tau_c=\infty$) and infinity (``aggressive
%control'' at the limit to instability), with a typical ``good'' value of about 1.6.
%, which is the peak value as a function of frequency of the sensitivity function. 
%Notice that $M_s <1.7$ guarantees gain margin (GM)$>2.43$ and phase margin (PM)$>34.2^o$  \citep{Rivera:86}.

%\een
%In general, we have a multiobjective optimization problem, so there is no value of $\tau_c$ which is ``optimal''.
%We will consider in more detail the two limiting cases of ``tight'' and ``smooth'' control, and also consider in some detail
%the required input usage.





%, and Table~\ref{tabsimcpo} 
%gives the tunings, costs $J$ and $M_s$-values for the SIMC PI-controller (with $\tau_c=\theta$). 



%$$    {\rm IAE}^o_{ys} = {\rm IAE}_{ys}(c_{ys,opt}, M_s=1.59) $$  
%$$    {\rm IAE}^o_{d} =  {\rm IAE}_{d}(c_{d,opt},   M_s=1.59) $$  
%The reference values IAE$^o_{ys}$ and IAE$^o_{d}$ are given for four processes in Table~\ref{tabref0}.
\input{table_wheight}


\section{Optimal PI Tuning}
%\subsection{Pareto-Optimal Tuning Curve}

For a given first order with time delay process, the Pareto-optimal curve with PI control is generated by solving the the following optimization problem: 

\parbox{\columnwidth}{
For a given set of $M_s$ values $\mathcal{M}_s=\{1.1,...,1.59,...,3\}$, solve:

$$ \min_{c } J(c) = 0.5 \left[
{{\rm IAE}_{ys}(c) \over {\rm IAE}^\circ_{ys}} +  
{{\rm IAE}_{d}(c) \over {\rm IAE}^\circ_{d}} \right]$$
$$ \rm s.t. ~~~~ \it M_s=m$$

for all $m \in \mathcal{M}_s$. }% This was done for four process models

%$$g=\frac{1}{\tau_1 s+1}\exp(-s) $$
%where the process time constant varied between  $\tau_1=0$ (pure time delay), $\tau_1=1$ (delay dominate), $\tau_1=8$ (lag %dominated) to $\tau_1=\infty$ (integrating process). 

%The Pareto-optimal curves for these processes can be seen in Figure \ref{figPo}.

In Table~\ref{tabref0}, the resulting optimal PI-controllers and J-values are given for $M_s=1.59$ for four processes;
\begin{itemize}
	\item pure time delay ($\tau_1/\theta=0$), 
	\item small time constant ($\tau_1/\theta=1$), 
	\item intermediate time constant ($\tau_1/\theta=8$), 
	\item integrating ($\tau_1/\theta=\infty$).
\end{itemize}
We note that $J=1$ for a time delay process, because there is no trade-off between disturbances and setpoints for this process, and because the reference controllers have $M_s=1.59$.
For the integrating process, the optimal value of $J$ is 1.50, mainly because we have to sacrifice setpoint performance.


\subsection{Optimal trade-off between robustness and performance}


\begin{figure*}[p]

	\parbox{1\textwidth}{%
		\begin{center}
	\hspace*{-30mm}
	\includegraphics[width=1.3\textwidth]{fig/stepPo.eps}
		\vspace*{-10mm}
	\caption{Optimal PI control of first-order plus delay processes: 
		Output response (left) and input usage (right) for $M_s=1.20$ (top), $M_s=1.59$ (middle) and $M_s=2.00$ (bottom).
		Setpoint change at $t=0$ and (input) disturbance at $t=20$.}
	\label{stepResponse}
	\end{center}}
%\end{figure*}


%\begin{figure*}[h]
	\vspace*{10mm}
	\parbox{\textwidth}{%
	\begin{center}
		\includegraphics[width=0.81\columnwidth]{fig/poRes0.eps}
		\includegraphics[width=.7954\columnwidth]{fig/poRes1.eps}
		\includegraphics[width=0.81\columnwidth]{fig/poRes8.eps}
		\includegraphics[width=0.81\columnwidth]{fig/poResInf.eps}
	\end{center}
	\caption{Pareto-optimal trade-off between robustness ($M_s$) and performance ($J$) with PI control for four processes} 			 
	 \label{figPo}}
\end{figure*}

Figure \ref{stepResponse} shows the optimal output response (left) and corresponding input usage (right) for three values of $M_s$ (1.2, 1.59, and 2) for the four processes.
We see that as $M_s$ is increased and control gets more aggressive, the output performance is improved whereas the input usage gets worse.

\begin{figure}[htb]
	\begin{center}
	\includegraphics[width=\columnwidth]{fig/sensitivityVsFreq.eps}
	\caption{Optimal PI control of first-order plus delay process: Magnitude of sensitivity function for three $M_s$-values}
	\label{sensitivityVsFreq}
	\end{center}
\end{figure} 

\begin{figure}[!hb]
	\begin{center}
		\includegraphics[width=\columnwidth]{fig/tvYs.eps}
		\includegraphics[width=\columnwidth]{fig/tvD.eps}
	\end{center}
	\caption{Optimal PI control of first-order plus delay process: Input usage (TV) for setpoint change (upper) and input
		 disturbance (lower) as a function of $M_s$.}
	\label{tvVsMs}
\end{figure}






Figure \ref{sensitivityVsFreq} shows the corresponding optimal sensitivity function, $S=1/(1+gc)$, 
for the process \mbox{$g(s)=\frac{1}{8s+1}\exp(-s)$} for the same three $M_s$-values. 
We see how tightening the control (increasing $M_s$) gives better performance (lower $|S|$) at lower frequencies, but 
worse performance at higher frequencies. The worst-case peak for $|S|$ is the $M_s$-value which occurs approximately at the frequency corresponding
to the time delay, $\omega=\nicefrac{1}{\theta}=1$ (the unit depends on the unit used for time).

The main results are given in Figure~\ref{figPo} which gives the Pareto-optimal trade-off between performance ($J$) and robustness ($M_s$) (blue curve) for the four processes.
Note that we have a real trade-off between performance ($J$) and robustness ($M_s$) only when there is a negative slope between these variables, which is in the left region in the figures in Figure~\ref{figPo}. 
Here, we have to compromise  between performance and robustness. 
We never want to be in the region with a zero or positive slope (marked as ``uninteresting''), because if we move to the right in this region, both performance and robustness deteriorate at the same time. 
%In other words,  in this area we can improve both performance ($J$) and robustness ($M_s$) at the same time by moving to the left. 
The minimum point on the curve (dashed black line) represents the maximum $M_s$-value that we would want to use. 
This maximum $M_s$-value varies from about 1.9 for the delay-dominant process to about 2.7 for the integrating process.
$M_s$-values higher than this give oscillating responses which increases the IAE-value.

%with another choice for the tuning parameter (moving to the left). 
%(i.e. using a larger value for $\tau_c$).



\subsection{Correlation between robustness and input usage}

As mentioned earlier, we chose to not consider the input usage when finding the optimal controller. 
One reason is that this is not necessary, because the input usage in terms of total variation (TV) correlates very well with $M_s$. 
This is shown in Figure \ref{tvVsMs} for the four processes. We see that when $M_s$ increases, TV increases for both setpoint change and input disturbance. 
Because the system is subjected to unit step changes in setpoint and input disturbance, the minimum TV-value is 1,
except for the integrating process $\tau_1/\theta=\infty$ where the minimum is 0. 
Also note the TV-value is relatively constant around 1 at low $M_s$-values. It is only when $M_s$ is so large that the system starts oscillating, that TV starts increasing.
This illustrates that $M_s$ is better to use as a robustness measure (objective 2) than TV.
%For the pure time delay process ($\tau_1/\theta=0$) and the small time constant process ($\tau_1/\theta=1$), $M_s$ is better quantifying robustness, as $TV$ is almost constant for $M_s<1.6$.


\subsection{Optimal PI tuning parameters}



\begin{figure}[!b]
	\begin{center}
	\includegraphics[width=\columnwidth]{fig/costFunVsLag.eps}
	\caption{Optimal PI control of first-order plus delay process: Setpoint and disturbance performance as a function of process
			 time constant $\tau_1/\theta$ ($M_s=1.59$).}
	\label{costFunVsLag}
	\end{center}
\end{figure}

\begin{figure*}[!ht]%
       \centering
       \parbox{.47\textwidth}{%
				\includegraphics[width=0.915\columnwidth]{fig/TuningGainVsLagMs.eps}
				\includegraphics[width=0.915\columnwidth]{fig/TuningIntegralVsLagMs.eps}
				\caption{Optimal PI-settings ($K_c, \tau_I$) for first-order with delay processes as a 
						function of $\tau_1/\theta$ for five values of $M_s$. }
				\label{settingsVsLag}} 
				\qquad
     \parbox{.47\textwidth}{%%\begin{minipage}{.5\textwidth}%
     \vspace*{-3mm}
					\includegraphics[width=0.90\columnwidth]{fig/TuningGainVsLagMsZoomFinal.eps}
					\includegraphics[width=0.90\columnwidth]{fig/TuningIntegralVsLagMsZoomFinal.eps}
					\caption{Optimal PI-settings ($K_c, \tau_I$) for first-order with delay processes for small values of $\tau_1/\theta$}
					\label{settingsVsLagZoom}}
%       \end{minipage}%
     \end{figure*}%
     







     


%
%
%\begin{figure}
%\includegraphics[width=0.49\textwidth]{fig/TuningGainVsLagMs.eps}
%\includegraphics[width=0.49\textwidth]{fig/TuningIntegralVsLagMs.eps}
%\caption{Optimal PI-settings ($K_c, \tau_I$) for first-order with delay processes as a function of $\tau_1/\theta$ for five values of $M_s$. }
%\label{settingsVsLag}
%\end{figure}
%
%\begin{figure}
%\includegraphics[width=0.49\textwidth]{fig/TuningGainVsLagMsZoomFinal.eps}
%\includegraphics[width=0.49\textwidth]{fig/TuningIntegralVsLagMsZoomFinal.eps}
%\caption{Optimal PI-settings ($K_c, \tau_I$) for first-order with delay processes for small values of $\tau_1/\theta$}
%\label{settingsVsLagZoom}
%\end{figure}
%\end{figure*}

%
%
%\begin{figure*}
%
%\subfloat[Optimal PI-settings ($K_c, \tau_I$) for first-order with delay processes as a function of $\tau_1/\theta$ for five values of $M_s$. ]{
%\includegraphics[width=0.49\textwidth]{fig/TuningGainVsLagMs.eps}
%\includegraphics[width=0.49\textwidth]{fig/TuningIntegralVsLagMs.eps}
%}{settingsVsLag}
%
%\subfloat[Optimal PI-settings ($K_c, \tau_I$) for first-order with delay processes for small values of $\tau_1/\theta$]{
%\includegraphics[width=0.49\textwidth]{fig/TuningGainVsLagMsZoomFinal.eps}
%\includegraphics[width=0.49\textwidth]{fig/TuningIntegralVsLagMsZoomFinal.eps}
%}{settingsVsLagZoom}
%
%\end{figure*}


The optimal PI controller settings are plotted in Figures \ref{settingsVsLag} and \ref{settingsVsLagZoom} as a function of the scaled process time constant $\tau_1/\theta$ for five values of $M_s$. 
The SIMC-settings for the choice \mbox{$\tau_c=\theta$} are included as a  reference (dashed line). 
The are two main regions: delay dominant ($\tau_1/\theta < 5$) and lag dominant ($\tau_1/\theta > 5)$.  
For lag-dominant processes (Figure \ref{sensitivityVsFreq}), the scaled controller gain $K_ck\theta / \tau_1$ approaches a constant value as we increase the lag $\tau_1/\theta$. 
For example, for $M_s=1.59$, the optimal value is $K_ck\theta / \tau_1=0.4145$ for $\tau_1/\theta = 50$, and 0.4100 for $\tau_1/\theta =\infty$ (integrating process).
The same can be observed for the integral time, which approaches $\tau_I=6.3 \theta$ for an integrating process.
Note that  when $M_s$ increases (less robustness), the controller gain increases and the integral time decrease.
%It might seem strange that there is such a big difference between the SIMC settings and the optimal settings for lag dominated processes (less than 10\%). For processes where $\tau_1/\theta>8$ the $M_s$ value is no longer constant at $M_s=1.59$, but increasing slowly to $M_s=1.7$, so a direct comparison between SIMC and another line can be misleading. 

The delay-dominant region ($\tau_1/\theta<5$) is shown more clearly in Figure \ref{settingsVsLagZoom}. 
Interestingly, the optimal integral time is almost independent of the robustness ($M_s$) in this region. 
We note that  $\tau_I \approx \tau_1$ (dashed line) for small values of $\tau_1/\theta$; more exactly for  $\tau_1/\theta$ between about 1 and 3. This agrees with the well-known IMC rule  \citep{rivera:86}.
However, $\tau_I$ does not approach zero (integrating controller) as $\tau_1/\theta$ becomes small. Rather, for a pure
time delay processes \mbox{($\tau_1/\theta=0$)}, the integral time is approximately $\theta/3$ for all $M_s$-values of interest (which from Figure \ref{figPo} is $M_s<1.9$).
On the other hand, the controller gain ($K_c$) increases with $M_s$ also for small values of $\tau_1/\theta$). 


%This is not the case the lag dominated region (Figure \ref{settingsVsLag}), where integral time is clearly dependant on $M_s$. Also the optimal integral time fits nicely with $\tau_I=\tau_1$.

\subsection{Output disturbance (setpoint) versus input disturbance trade-off}

As already noted from Table~\ref{tabref0}, the optimal controller that minimizes the cost $J(c)$ in (\ref{eq:jpo}), 
puts more emphasis on the (input) disturbance rejection (IAE$_d$) than on the setpoint tracking (IAE$_{ys}$), especially for larger values of $\tau_1/\theta$.
This is further illustrated in Figure~\ref{costFunVsLag}, where we see that for $\tau_1/\theta$ less than about 5, the optimal controller puts about equal weight on setpoints and disturbances, 
whereas for larger values it mostly emphasizes disturbances. For example, for an integrating process, we find
IAE$_d$/IAE$_d^\circ=1.02$ (close to optimal for input disturbance), whereas IAE$_{ys}$/IAE$_{ys}^\circ=1.99$ (twice the optimal). 
% \{(Sigurd: I think 2.1 should be 1.99.
This is a bit surprising since we have equal weights on the two terms in the cost function. 
The reason is that the IAE-optimal controller for setpoint (IAE$_{ys}^\circ$) has no integral action ($\tau_I=\infty$), whereas integral action is required to bring the process back to its setpoint
in reasonable time when there are input disturbances, and this consideration dominates when want to find a single optimal controller.
Even if we increase the emphasis on setpoints, we find a similar pattern. For example, if we use a cost function with only a small weight
(1\%) on input disturbances,
$$ J(c) =  \left[0.99
{{\rm IAE}_{ys}(c) \over {\rm IAE}^\circ_{ys}} +  0.01
{{\rm IAE}_{d}(c) \over {\rm IAE}^\circ_{d}} \right] $$
we find for an integrating process the optimal settings $K_c=0.462$ and $\tau_I = 12.2\theta$, with IAE$_d$/IAE$_d^\circ=1.72$ and IAE$_{ys}$/IAE$_{ys}^\circ=1.91$. % Flere desimaler,   1.8270   1.9889
We notice that there is only a minor improvement in setpoint performance 
(IAE$_{ys}$ decreases about 4\%), %  from 1.99 to 1.83.  % CHRISS: todo, DU M� SJEKKE  (jeg endret fra 7 til 14%)  chriss: ok! should be after should be 1.92, change 1-1.92/1.98=0.303. 
whereas disturbance rejection is much worse (IAE$_d$ increases about 69\%). % changes from 1.02 to 1.83 % CHRISS: todo, Sjekk! chriss: OK!
The conclusion is this that we may put emphasis mainly on input disturbances when tuning PI controllers.

%Notice that there is no trade-off for pure time delay processes (Figure \ref{costFunVsLag}), where the cost function puts equal emphasis on both setpoint change (output disturbance ) and input disturbance. This is because the process only delays the input disturbance, and both input and output disturbances have the same effect on the output $y$. 




\section{Evaluation of the SIMC PI tuning rule}

\begin{figure*}[htb]
	\begin{center}
	\includegraphics[width=.81\columnwidth]{fig/simcRes0.eps}
	\includegraphics[width=.7954\columnwidth]{fig/simcRes1.eps}

	\includegraphics[width=.81\columnwidth]{fig/simcRes8.eps}
	\includegraphics[width=.81\columnwidth]{fig/simcResInf.eps}
	\end{center}
	\caption{Evaluation of optimality of SIMC PI tuning rule for four processes.}
	\label{figSimc}
\end{figure*}     
     



%The idea behind the SIMC tuning rules  \citep{skogestad:03}  was to have a simple analytical derived rule with good closed-loop performance. 

The SIMC PI-settings \citep{skogestad:03} for a first-order with delay model are% (cascade or ideal form)
$$ K_c = {1 \over k} {\tau_1 \over \tau_c + \theta}   \label{kcsispi} $$
$$ \tau_I = % \min\{\tau_1, {4 \over k'\ K_c}\}  =
            \min\{\tau_1, 4 (\tau_c + \theta)\}  $$
where the desired first-order {\em closed-loop} time constant $\tau_c$ is the only tuning parameter. 
For a ``fast and robust'' setting, $\tau_c=\theta$ is recommended.
The corresponding SIMC tuning parameters for $\tau_c=\theta$ 
are given for the four processes in Table~\ref{tabsimcpo}.

\input{table_simc}

\subsection{Effect of the tuning parameter $\tau_c$ on robustness}
By adjusting the tuning parameter $\tau_c$, we can shift the trade-off between performance and robustness. 
In particular, $M_s$ decreases (more robust) as we increase $\tau_c$. 
In Figure~\ref{robVsTc}, we show for four processes, the effect of changing $\tau_c/\theta$ in the region from 0 to 3, on (1) $M_s$-value, (2) gain margin, (3) phase margin and (4) allowed time delay error.
We note that there is a direct correlation between $\tau_c$ and all four robustness measures, and the correlation does not depend strongly on the value of $\tau_1/\theta$.
In particular, this applies to the gain margin (GM), which shows an almost one-to-one (and linear) relationship with $\tau_c$. Thus, for the SIMC-rule, specifying $\tau_c$ is essentially the same
as specifying GM.
%For a first order model with  $\tau_1/\theta=1$, and the recommended tuning parameter $\tau_c=\theta$ gives a $M_s=1.59$. A more aggressive tuning with $\tau_c=0.5\theta$ results in $M_s=1.92$, and a more robust tuning $\tau_c=1.5\theta$ gives a $M_s=1.44$ (Figure \ref{robVsTc}, top).
This linear relationship can be proved mathematically for processes with relatively small values of $\tau_1$, such that the SIMC-rule gives $\tau_I=\tau_1$ (see Appendix) % \ref{derTcGmPmDm}):
\beq GM = \frac{\pi}{2} \(\frac{\tau_c}{\theta}+1\)
\eeq
and Figure \ref{robVsTc} shows that this relationship holds well also for an integrating processes (where the SIMC-rule gives $\tau_I = 4 (\tau_c+\theta)$).

For relatively small values of $\tau_1$, such that $\tau_I=\tau_1$,  a similar linear relationship also exists for the allowable time delay error (see Appendix):
\beq
	\frac{\Delta\theta}{\theta}=\frac{\pi}{2}\( \frac{\tau_c}{\theta}+1\)-1=GM-1
\eeq 
However, Figure \ref{robVsTc} shows that this relationship does not quite hold for an integrating process where the allowable time delay error is somewhat smaller.

%For the process $\tau_1/\theta=8$ in Figure \ref{robVsTc}, there is a shift from $\tau_c<1$ where the robustness is similar to the integrating process, to $\tau_c\geq1$ where the robustness is equal the small time constant processes ($\tau_1/\theta\in\{0,1\}$). This is due to the switching in the integral time, where at small time constants $\tau_I=\tau_1$ and at large time constants \mbox{$\tau_I=4(\tau_c+\theta)$}.


\begin{figure}
	\begin{center}
		\includegraphics[width=.95\columnwidth]{fig/MsVsTc.eps}
		\includegraphics[width=.95\columnwidth]{fig/gmVsTc.eps}
		\includegraphics[width=.95\columnwidth]{fig/pmVsTc.eps}
		\hspace*{-7.5mm}
		\includegraphics[width=1.026\columnwidth]{fig/dmVsTc.eps}
	\end{center}
	\caption{SIMC: Effect of tuning parameter $\tau_c/\theta$ on robustness parameters: (1) $M_s$, (2) GM, (3) PM and (4) 	
			allowable time delay error $\nicefrac{\Delta \theta}{\theta}$}
	\label{robVsTc}
\end{figure}


\subsection{Optimality of the SIMC-rule}



Figure~\ref{figSimc} compares the Pareto-optimal performance (blue curve) with the SIMC PI-controller (green curve) for the four processes. 
The curve for the SIMC controller was generated by varying the tuning parameter $\tau_c$ from a large to a small value. The controllers corresponding to the choices
$\tau_c=1.5\theta$ (smooth tuning),  $\tau_c=\theta$ (default value) and  $\tau_c=0.5 \theta$ (more aggressive tuning) 
are shown by circles. 
Except for the pure time delay process, the differences between the IAE-values ($J$) for SIMC (green curve) and
optimal (blue curve) are within 10\%, which shows that the SIMC PI-rules are close to optimal. In other words, by adjusting $\tau_c$ we can generate the optimal controller for a given
desired robustness. 

%Note that we have a real trade-off between performance ($J$) and robustness ($M_s$) only when there is a negative slope between these
%variables (in the left region in the figures in Figure~\ref{figPO}). 
%We never want to be in the region with a zero or positive slope (to the right in the figures), because here we can improve both
%performance ($J$) and robustness ($M_s$) at the same time with another choice for the tuning parameter 
%(using a larger value for $\tau_c$).
Another important observation from Figure~\ref{figSimc} is 
that the SIMC default recommendation $\tau_c=\theta$ for ``tight'' control (as given by middle of the three circles)  
in all cases is located in the desired trade-off region with a negative slope, well before we reach the minimum.
Also, the recommended choice gives a fairly constant $M_s$-value in the region 1.59 to 1.7.   
From this we conclude that, except for the time delay process,
there is little room to improve on the SIMC PI-rule, at least when performance and
robustness are as defined above ($J$ and $M_s$). 
%This may change if the evaluation criteria are changed, for example, if we put more emphasis on disturbances (IAE$_d$) than on setpoints (IAE_{ys}$.  
%Recall from the derivation of the SIMC rules that the optimal integral time
%for setpoint changes is $\tau_I=\tau_1$ and we introduced the value $\tau_I=4(\tau_c+\theta)$ to get better
%response to input disturbances for integrating processes.   
%Evidently, as seen from the results for the integrating process in Figure~\ref{figPO}, this was a good choice, as it gives
%close-to-optimal performance ($J$) for given values of $M_s$, 

%The IAE-cost $J$ in (\ref{eq:jpo}) is based on equal weighting of
%servo (output disturbance) and regulator (input disturbance) performance. 
%The existence of a trade-off between servo and regulator performance,
%can be quantified by considering how much larger the (Pareto) optimal cost $J_{opt}$ (blue curve) is than 1 at the reference robustness, $M_s=1.59$, see also Table~\ref{tabref0}.
%For a  pure time delay-process, we have that  $J_{opt}=1$ for $M_s=1.59$ and there is no trade-off. 
%The reason is that the setpoint and output disturbance responses are the same.
%On the other hand, for the other extreme of an integrating process, 
%we have a clear trade off since the optimal PI-controller has $J_{opt}=1.51$ (the SIMC PI-controller with $M_s=1.59$ is close to this with $J$ about 1.6).
%The existence of the servo/regulator trade-off for an integrating process, implies that
%one for a given robustness ($M_s$-value)
%can find PI-settings with significantly better regulator (load disturbance) performance or better
%servo (setpoint) performance, but not both at the same time.  

\subsection{Evaluation of SIMC-rule for integral time}

When comparing the optimal PI-settings with the SIMC-rule for an integrating process, it seems that the integral time given by the SIMC-rule is a bit too large (Figure \ref{settingsVsLag}, bottom).
Specifically, for an integrating process and $\tau_c=\theta$  ($M_s=1.70$), the SIMC-rule gives $\tau_I/\theta=8$ whereas the optimal value for $M_s=1.70$ is about 5.5.
This indicates that in terms of the trade-off between setpoint and disturbance performance, the SIMC-rule is putting somewhat too much emphasis on setpoints.
There has also been claims that the SIMC-rule results in 
sluggish disturbance rejection for integrating processes \citep{haugen:10,ruscio:10}. 
To be able to shift the trade-off between setpoint and disturbance performance, 
one may introduce an extra parameter in the tuning rule \citep{alcantara:10}.
%we derived that the integral time should be
%$\tau_I=\tau_1$ for setpoint changes (servo performance), and we introduced the smaller value $\tau_I=4(\tau_c+\theta)$ to get better
%regulator performance for input disturbances for an integrating processes.   
%%Evidently, as seen from the results for the integrating process in Figure~\ref{figPO}, this was a good choice, as it gives
%Thus, with the SIMC approach, 
\cite{haugen:10} suggested to introduce an extra servo/regulator trade-off parameter 
$c$ in the expression for the integral time, 
\beq \tau_I = \min (\tau_1, c(\tau_c+\theta)) \eeq
where $c=4$ gives the original SIMC-rule. 
However, introducing an extra parameter adds complexity and from Figure~\ref{figSimc} the potential benefit does not seem sufficiently large.
Nevertheless, one may consider choosing another (lower) fixed value for $c$, 
and \cite{haugen:10} suggests using $c=2$ to improve the input disturbance performance. 
In Figure \ref{simc4Vs2} we check the optimality of choosing $c=2$ (dashed pink curve) for an integrating process.
We find indeed that performance ($J$) is improved compared to SIMC ($c=4$) if we keep $\tau_c=\theta$.
However, robustness is worse (with $M_s$ close to 2, whereas SIMC gives 1.7). 
More importantly, by decreasing $\tau_c$ for SIMC so that we get the same robustness, the performance with SIMC ($c=4$) 
is even better than with $c=2$.
In fact, we see that the SIMC performance curve (green) is closer to the Pareto-optimal curve in the 
entire optimal performance region  (where $J$ is small).   
In summary, we find that the value $c=4$ in the original SIMC-rule provides a well-balanced servo/regulator trade-off, and to improve disturbance performance for an integrating process we recommended decreasing 
the tuning constant $\tau_c/\theta$, say to around 0.7, rather than changing the value of $c$. 

\begin{figure}[tb]
	\begin{center}
		\hspace*{-9mm}
		\includegraphics[width=1.2\columnwidth]{fig/simc4Vs2.eps}
		\caption{Evaluation of suggested modified SIMC-rule ($c=2$) for integrating process.}
		\label{simc4Vs2}
	\end{center}
\end{figure}

\section{Improved SIMC-rule}

\begin{figure*}[htb]
	\begin{center}
		\includegraphics[width=.81\columnwidth]{fig/simcImpRes0.eps}
		\includegraphics[width=.81\columnwidth]{fig/simcImpRes1.eps}
		\includegraphics[width=.81\columnwidth]{fig/simcImpRes8.eps}
		\includegraphics[width=.81\columnwidth]{fig/simcImpResInf.eps}
	\end{center}
	\caption{Evaluation of optimality of improved SIMC PI tuning rule for four processes.}
	\label{figImpSimc}
\end{figure*}  


%\begin{figure}
%\begin{center}
%\includegraphics[width=\columnwidth]{fig/gmSimcVsSimcImp.eps}
%\includegraphics[width=\columnwidth]{fig/pmSimcVsSimcImp.eps}
%\includegraphics[width=\columnwidth]{fig/dmSimcVsSimcImp.eps}
%\end{center}
%\caption{Robustness comparison between SIMC (red curve) and improved SIMC (green curve)}
%\label{simcVsImpRobust}
%\end{figure}

From the results in Figure~\ref{figSimc}, the main ``problem" with the SIMC-rule is for pure time delay processes, 
where we see that the IAE-value ($J$) for the SIMC controller
is about 40\% higher than the minimum  with the same robustness ($M_s$).
This is further illustrated by the closed-loop simulations in Figure~\ref{fig-delay}.
We see that the response with the SIMC PI-controller 
(denoted SIMC-original in the figure), although being nice and smooth,
is somewhat sluggish initially, because it is actually a pure I-controller 
with $K_c=0, \tau_I=0$ and $K_I=K_c/\tau_I=0.5$. 
%This if for a pure time delay process with $\theta=1$ and $k=1$. 
On the other hand, the IAE-optimal PI-controller (with minimum $J$ for $M_s=1.59$) 
has $K_c$ about 0.2 and $\tau_I$ about 0.32 (and $K_I = 0.62$). 

Fortunately, as already observed, the optimal PI-controller for a pure time delay process (green line in Figure~\ref{figPo}), 
has an almost fixed integral time of approximately $\theta/3$ for all values of $M_s$ between 1.4 and 1.7. (Figure \ref{settingsVsLagZoom}, bottom)
Based on this fact, we propose a simple change to the SIMC-rule, 
namely to replace $\tau_1$ by $\tau_1+\theta/3$ in the rule (PI control), 
which markedly improved the responses for a pure time delay process.
It is important that the change is simple because ``simplicity'' was one of the main objectives 
when originally deriving the SIMC-rule.

A similar change, but with $\theta/2$ rather than $\theta/3$, was originally proposed by \cite{rivera:86} for their ``improved PI" tuning rule, and the effectiveness of this modification is also
clear from the paper of \cite{foley:05}. 
% CUT derivation 05 April 2011
%To derive this modification, we return to the second-order time delay process (\ref{fod2}) and
%the resulting IMC (Smith Predictor) controller in (\ref{sp}). Instead of approximating the delay
%in (\ref{sp}) using $ e^{-\theta s} \approx 1 - \theta s$ (first-order Taylor), which resulted in a
%controller with $\tau_I=0$ for a pure time delay process, we instead apply
% the first-order Pade approximation,
%$$  e^{-\theta s} \approx {1 - (\theta/2)s \over 1 + (\theta/2)s} $$
%and derive from (\ref{sp}) the controller  
%\beq 
%c(s) = {(\tau_1 s + 1)(\tau_2 s+1) \over k} 
%{ (\theta/2)s+1 \over (\tau_c  + \theta) s (\tau_f s + 1)}  
%\label{sp2} \eeq
%where $\tau_f = {\tau_c \theta \over 2(\tau_c+\theta)}$.
%The ``filter'' time constant $\tau_f$ varies between 0 and $\theta/2$, depending on the value of $\tau_c$.
%Thus, the filter term $1/(\tau_f s +1)$ is outside the closed-loop bandwidth and is neglected.
%Next, for a first-order process ($\tau_2=0$), we introduce for the numerator the approximation
%$$ (\tau_1 s + 1) ((\theta/2)s+1) = 
%\tau_1 (\theta/2) s^2 + (\tau_1 + \theta/2)s + 1 \approx
%(\tau_1 + \theta/2)s + 1 $$
%which is exact for the limiting case of a pure time delay process ($\tau_1=0$).
%This results in the same SIMC PI-rules as before, but with $\tau_1$ replaced by
%$\tau_1 + (\theta/2)$ in the rules, which is the ``improved'' PI-rule of \citeasnoun{rivera:86}. 
%%Similarly, for a second-order process ($\tau_2>0$) we introduce the approximation
%%$$ (\tau_2 s + 1) ((\theta/2)s+1) = 
%%\tau_2 (\theta/2) s^2 + (\tau_2 + \theta/2)s + 1 \approx
%%(\tau_2 + \theta/2)s + 1 $$
%%which is exact for $\tau_2=0$ and results in the same SIMC PID-rules as before but with
%%$\tau_2$ replaced by $\tau_2 + \theta/2$ in the rules.
%%$$ K_c = {1 \over k} {\tau_1 + \theta/2 \over \tau_c + \theta},  \q \tau_I = \tau_1 + \theta/2 $$
%For a pure time delay process, this give a PI controller with an integral time $\tau_I=\theta/2$ (independent of the value for the
%tuning parameter $\tau_c$).  
However, as seen in Figure~\ref{fig-delay}, the response with this IMC PI controller also settles rather slowly towards the setpoint,
indicating that the integral time $\theta/2$ is too large.
%The proposed value $\theta/3$ gives a faster settling and has the additional advantrage of being closer
%to the original SIMC-value of 0. 
The conclusion is that we recommend to replace $\tau_1$ by $\tau_1+\theta/3$ in the SIMC-rule to get the improved SIMC-rule:
%\begin{svgraybox}
%{\bf Improved SIMC PI-rule for first-order with delay process.}  
%Replace $\tau_1$ by $\tau_1 + \theta/2$ in the original rules to get:
\beq K_c = {1 \over k} {\tau_1 + {\theta \over 3} \over \tau_c + \theta}  \label{kcsispi-improved} \eeq
\beq \tau_I = % \min\{\tau_1, {4 \over k'\ K_c}\}  =
            \min\{\tau_1+ {\theta \over 3}, 4 (\tau_c + \theta)\}  
\label{tauisispi-improved} \eeq



%\end{svgraybox}
The improvement of this rule for a pure time delay processes is clear from the red curves in 
Figures ~\ref{figImpSimc} (upper left) and ~\ref{fig-delay}; 
for small $M_s$-values the improved SIMC-controller is almost identical to the Pareto-optimal, which confirms
that $\tau_I=\theta/3$ is close to optimal for a pure time delay process. 
One disadvantage with the improved rule, which is not considered in this analysis, is that the larger controller gain ($K_c$) makes the input and output responses more sensitive to measurement noise.

For the process with a small time constant ($\tau_1/\theta=1$), the improved
SIMC-rule (red curve in upper right plot in Figure~\ref{figImpSimc}) is slightly better than the ``original'' SIMC-rule (green curve) for higher $M_s$-values (where we get better performance) 
but slightly worse for lower $M_s$-values. For the two processes with a large time constant ($\tau_1/\theta=8$ and $\tau_1/\theta=\infty$) there is,
as expected, almost no difference between the original and improved SIMC-rules.

\begin{figure}[htb]
	\begin{center}
		\includegraphics[width=\columnwidth]{fig/stepResponse.eps}
	\end{center}
	\caption{Closed-loop setpoint responses for pure time delay process ($\theta=1, k=1, \tau_1=0$) with PI-control.
		%Optimal PI-controller: $K_c=0.197$ and $\tau_I=0.32\ (K_I=K_c/\tau_I=0.62$).
	\footnotesize 
	$\bullet$~IMC~PI:~$K_c=0.29$,  $\tau_I=0.5\ (K_I= K_c/\tau_I = 0.58)$.
	$\bullet$~SIMC~PI~original ($\tau_c=\theta$): $K_c=0, \tau_I=0 (K_I=0.5)$.
	$\bullet$~SIMC~PI~improved ($\tau_c=0.61\theta): K_c=0.207, \tau_I=0.333 (K_I=0.62)$.  
	All three controllers have the same robustness ($M_s=1.59$). 
	For a pure time delay process, the setpoint and disturbance responses are identical, and the input and output are identical.}
	\label{fig-delay}
\end{figure}

%s = tf('s')
%tys=0;td=100;
%thetasim=1, k=1; ys=1; tau1=0;
% g=exp(-thetasim*s)
%theta=1; % used for tuning controller 
%'SIMC:',taui=0.0001; Kc=0.5*taui/(theta*k);taud=0; sim('tunepid');t1=Tid;u1=u;y1=y;e1=y-ys; epsd=0.01; c = Kc*(taui*s+1)*(taud*s+1)/(taui*s*(epsd*taud*s+1)); Ms=ms(g,Kc,taui), IAE=iae(Tid,y-ys) % SIMC
%'ZN:',Kc=0.45; Ki=0.27; taui=Kc/Ki;taud=0; sim('tunepid');t2=Tid;u2=u;y2=y;e2=y-ys; epsd=0.01; c = Kc*(taui*s+1)*(taud*s+1)/(taui*s*(epsd*taud*s+1)); Ms=ms(g,Kc,taui), IAE=iae(Tid,y-ys) % ZN
%'Optimal:',Kc=0.197, taui=0.32,taud=0; sim('tunepid');t3=Tid;u3=u;y3=y;e3=y-ys; epsd=0.01; c = Kc*(taui*s+1)*(taud*s+1)/(taui*s*(epsd*taud*s+1)); Ms=ms(g,Kc,taui), IAE=iae(Tid,y-ys) % Optimal
%'IMC-PI',Kc=0.29, taui=0.5, taud=0; sim('tunepid');t4=Tid;u4=u;y4=y;e4=y-ys; epsd=0.01; c = Kc*(taui*s+1)*(taud*s+1)/(taui*s*(epsd*taud*s+1)); Ms=ms(g,Kc,taui), IAE=iae(Tid,y-ys) % IMC-PI
%'SIMC improved (Ms=1.59):',tauc=0.61*theta;Kc=(1/k)*(tau1+theta/3)/(tauc+theta), taui=tau1+theta/3,taud=0; sim('tunepid');t5=Tid;u5=u;y5=y;e5=y-ys; epsd=0.01; c = Kc*(taui*s+1)*(taud*s+1)/(taui*s*(epsd*taud*s+1)); Ms=ms(g,Kc,taui), IAE=iae(Tid,y-ys) % Improved SIMC (Ms=1.59)
%'SIMC improved (tauc=1):',tauc=1*theta;Kc=(1/k)*(tau1+theta/3)/(tauc+theta), taui=tau1+theta/3,taud=0; sim('tunepid');t6=Tid;u6=u;y6=y;e6=y-ys; epsd=0.01; c = Kc*(taui*s+1)*(taud*s+1)/(taui*s*(epsd*taud*s+1)); Ms=ms(g,Kc,taui), IAE=iae(Tid,y-ys) % Improved SIMC (tauc=1)
%plot(t1,y1,t4,y4,'--',t3,y3,'-.',t5,y5,t6,y6,t2,y2,'.',[0 20],[1 1],':');  axis([0 10 -0.1 1.2])


%{\bf Possible improved SIMC PID-rule for second-order with delay process.} 
%A similar derivation to that following (\ref{sp2}) for the first-order process, 
%gives that we for a second-order process should replace $\tau_2$ by $\tau_2 + \theta/2$ 
%in the original rule (\ref{taudsis}), while keeping $K_c$ and $\tau_I$ unchanged.
%%
%%\beq K_c = {1 \over k} {\tau_1 \over \tau_c + \theta}  \label{kcsis-improved} \eeq
%%\beq \tau_I = % \min\{\tau_1, {4 \over k'\ K_c}\}  =
%%              \min\{\tau_1, 4 (\tau_c + \theta)\}  
%%\label{tauisis-improved} \eeq
%%\beq \tau_D = \tau_2 + \theta/2 \label{taudsis-improved} \eeq
%This is similar to the IMC-PID rule of \citeasnoun{rivera:86}. 
%A more careful analsysis needs to be done to check the possible improvement, and to see if, for example, 
%$\theta/3$ is better as found for the improved SIMC-PI controller. 
%
\section{Obtaining the model from a closed-loop setpoint experiment}

In some cases, open-loop responses may be difficult to obtain, and using closed-loop data may be more effective.
The most famous closed-loop experiment is the Ziegler-Nichols where the system is brought to sustained oscillations
by use of a P-only controller. One disadvantage with the method is that the system is brought to its instability limit. Another
disadvantage is that it does not work for a simple second-order process. Finally, only two pieces of information are used
(the controller gain $K_u$ and the ultimate period $P_u$), so the method cannot possibly work on a wide range of first-order plus delay processes, which we know
are described by three parameters ($k, \tau_1, \theta$).

\cite{yuwana:82}, and more recently \cite{shams:10},  proposed a modification 
to the Ziegler-Nichols closed-loop experiment, which does not suffer from these three disadvantages.
Instead of bringing the system to its limit of stability, one uses a P-controller with a gain that is about half this value, 
such that the resulting overshoot ($D$) to a step change in the setpoint is about 30\% (that is, $D$ is about 0.3). 

We here describe the procedure proposed by \cite{shams:10}.
% which seems to use the most easily available parameters from the closed-loop response.
The system should be at steady-state initially, that is, before the setpoint change is applied.
Then, from the closed-loop setpoint response one obtains the following  parameters
(see Figure~\ref{fig-clpresponse}):

\begin{figure}[tb]
	\includegraphics[width=\columnwidth]{cl-presponse2.pdf}
	\caption{Extracting information from closed-loop setpoint response with P-only controller.}
	\label{fig-clpresponse}
\end{figure}

\bb
	\item Controller gain used in experiment, $K_{c0}$
	\item Setpoint change, $\Delta y_s$
	\item Time from setpoint change to reach first (maximum) peak, $t_p$
	\item Corresponding maximum output change, $\Delta y_p$
	\item Output change at first undershoot, $\Delta y_u$
\eb

This seems to be the information that is most easy (and robust) 
to observe directly, without having to record and analyse all the data before finding
the parameters. Also note that one may stop the experiment already at the first undershoot. 

The undershoot $\Delta y_u$ is used to estimate the steady-state
output change (at infinite time) \citep{shams:10}, 
\beq \Delta y_{\infty} = 0.45(\Delta y_p + \Delta y_u) \eeq
Alternatively, if one has time to wait for the experiment to settle,
one may record $\Delta y_{\infty}$ instead of $\Delta y_u$.

From this information one computes the relative overshoot and the absolute value of the relative steady-state offset,
defined by:
\bb
	\item Overshoot, $D = {\Delta y_p - \Delta y_{\infty} \over \Delta y_{\infty}}$.
	\item Steady-state offset, $B=\left| { \Delta y_s-\Delta y_{\infty} \over \Delta y_{\infty}} \right|$.
\eb

\cite{shams:10} use this information to obtain directly the PI settings. 
Alternatively, we here propose a two-step procedure, where we first from 
$K_{c0}, D, B$ and $t_p$ obtain estimates for the parameters in a first-order plus delay
model. We compute the parameters
$$ A = 1.152D^2 - 1.607D + 1 $$
$$ r = 2A/B $$
and obtain the following first-order plus delay model (see Appendix for derivation):
\beq k = 1/(K_{c0}B)  \label{k-shams} \eeq
\beq \theta = t_p\cdot (0.309 + 0.209e^{-0.61r}) \label{theta-chriss} \eeq
\beq \tau_1 = r \theta \label{tau1-chriss} \eeq
These values may subsequently be used with any tuning method, for example, the SIMC PI-rule.
The closed-loop method may also be used for an unstable process, provided it can be approximated reasonably
well by a stable first-order process. The extension to unstable processes is the reason
for taking the absolute value when obtaining the steady-state offset $B$.

\begin{figure}[hb]
	\includegraphics[width=\columnwidth]{fig/stepOlBareGcl.eps}
	\caption{Closed-loop setpoint identification method: Comparison of resulting open-loop step response with true model.
		%\mbox{$g_0(s) = {(-0.3s + 1)(0.08s + 1 ) \over (2 s+1)(1s+1)(0.4s+1)(0.2s+1)(0.05s+1)^3} $} (solid line),
		%and comparison with various approximations.
		}
	\label{fig-gapprox}
\end{figure}

\xx {\bf Example}. {\it For the process \citep{skogestad:03}
$$ g_0(s) = {(-0.3s + 1)(0.08s + 1 ) \over (2 s+1)(1s+1)(0.4s+1)(0.2s+1)(0.05s+1)^3} $$
the closed-loop setpoint response  with P-only controller with gain $K_{c0}=1.5$ is shown in Figure~\ref{fig-clpresponse}.
The following data is obtained from the closed-loop response
$$ K_{c0}=1.5, \Delta y_s=1, \Delta y_p = 0.79, t_p = 4.4, \Delta y_u = 0.54$$
and we compute
$$ \Delta y_{\infty} = 0.5985, D = 0.32, B = 0.67, A = 0.6038, r = 1.80 $$
which using (\ref{k-shams}) - (\ref{tau1-chriss}) gives the following first-order with delay model approximation,
\beq k = 0.994,  \theta = 1.67,  \tau_1 = 3.00 \label{fod-parameters-cl} \eeq
This gives a good approximation of the open-loop step response, as can seen by comparing the curves for
$g_0$ (blue) and $g_{cl}$ (green) in Figure~\ref{fig-gapprox}.  The approximation is certainly not the best possible, but it should 
be noted that the objective is to use the model for tuning, and the resulting difference in the tuning, and thus
closed-loop response, may be smaller than it appears by comparing the open-loop responses.  
}
%Kc0=1.5
%dys=1
%dyp=0.79
%tp=4.4
%dyu=0.54
%
%dyinf = 0.45*(dyp + dyu)
%D =(dyp -dyinf)/dyinf
%B =abs((dys -dyinf)/dyinf)
%
%A = 1.152*D^2 - 1.607*D + 1.0
%r = 2*A/B
%k = 1/(Kc0*B)
%theta = tp*[0.309 + 0.209*exp(-0.61*r)]
%tau = theta*r




%s = tf('s')
%g = (-0.3*s+1)*(0.08*s+1)/((2*s+1)*(s+1)*(0.4*s+1)*(0.2*s+1)*(0.05*s+1)^3)
%g1 = exp(-1.47*s)/(2.5*s+1)
%g2 = exp(-0.77*s)/((2*s+1)*(1.2*s+1))
%gcl = 0.994*  exp(-1.67*s)/(3.0*s+1)
%step(g,g1,g2,gcl,10)

\section{conclusion}

Comparing the performance of the SIMC-rule with the optimal for a given robustness ($M_s$ value) shows that
the SIMC-rule give settings close to the Pareto-optimal (Figure~\ref{figSimc}).
This means that the room for improving
the SIMC PI-rule is limited, at least for the first-order plus delay processes considered in this paper, and
with a good trade-off between rejecting input and output (setpoint) disturbances.

The tuning parameter $\tau_c$ should be chosen to get the desired trade-off between
fast response (small IAE) on the one side, and smooth input usage and robustness (small $M_s$)
on the other side.
The recommended choice $\tau_c=\theta$ gives robust ($M_s$ about 1.6 to 1.7) 
and somewhat conservative settings when compared with
most other tuning rules. If it is 
desirable to get faster control one may consider reducing $\tau_c$
to about $\theta/2$ (see Figure~\ref{figSimc}). 
More commonly, one may want to have ``smoother'' control with $\tau_c > \theta$ and a smaller controller
gain $K_c$. 

To improve the performance for delay-dominant processes, one may replace
$\tau_1$ by $\tau_1+ {\theta \over 3}$ and
use the ``improved'' SIMC PI-rule in (\ref{kcsispi-improved})--(\ref{tauisispi-improved}).


%The above tuning rules are derived analytically and are thus very suitable for teaching.
%Note that although the same formulas are used to obtain $K_c$ and $\tau_I$ for both PI- and PID-control ,
%the actual values will differ since the effective delay $\theta$ is smaller for
%a second-order model.

%Deleted 13 Mar02
%To get a faster response one may decrease the value of $\tau_c$, 
%and possibly further reduce the integral time. 
%However, one may also want to increase $\tau_c$ to get
%a slower and smoother response with less input usage.
%This results in a smaller controller gain $K_c$, but
%we must require $$K_c \ge |u_d|/y_{max} $$
%(approximately) in order to keep the output deviation less than $y_{max}$ in
%response to a disturbance that requires an input change of magnitude $|u_d|$.




%\begin{ack}
%Place acknowledgements here.
%\end{ack}

%\bibliography{article}             % bib file to produce the bibliography
\begin{thebibliography}{xx}

\harvarditem[Alcantara {\em et al.}]{Alcantara {\em et
al.}}{2010}{alcantara:10}
Alcantara, S., C.~Pedret and R.~Vilanova  (2010). On the model matching
approach to {PID} design: Analytical perspective for robust servo/regulator
tradeoff tuning. {\em J. Process Control} {\bf 20},~596--608.

\harvarditem[Astrom and Hagglund]{Astrom and Hagglund}{1995}{astrom:95}
Astrom, K.J. and T.~Hagglund  (1995). {\em {PID} controllers: {T}heory, design
and tuning}. 2nd ed.. {ISA}-Instrument Society of America.

%\harvarditem[Chien and Fruehauf]{Chien and Fruehauf}{1990}{chien:90}
%Chien, I.L. and P.S. Fruehauf  (1990). Consider {IMC} tuning to improve
%controller performance. {\em Chemical Engineering Progress} pp.~33--41.

\harvarditem[Di Ruscio]{Di Ruscio}{2010}{ruscio:10}
Di Ruscio, D.  (2010). On Tuning PI Controllers for Integrating Plus Time Delay Systems. {\em Modeling, Identification and Control} {\bf 31}(4),~145–-164.

\harvarditem[Foley {\em et al.}]{Foley {\em et al.}}{2005}{foley:05}
Foley, M.W., N.R Ramharack and B.R Copeland  (2005). Comparison of {PI}
controller tuning methods. {\em Ind. Eng. Chem. Res.} {\bf
44}(17),~6741--6750.

\harvarditem[Grimholt]{Grimholt}{2010}{grimholt:10}
Grimholt, C.  (2010). Verification and improvements of the {SIMC} method for
{PI} control. Technical report. 5th year project. Department of Chemical
Engineering. Norwegian University of Science and Technology, Trondheim.
\newblock http://www.nt.ntnu.no/users/skoge
\newblock ~~~~~/diplom/prosjekt10/grimholt.

\harvarditem[Haugen]{Haugen}{2010}{haugen:10}
Haugen, F.  (2010). Comparing {PI} tuning methods in a real benchmark
temperature control system. {\em Modeling, Identifiation and Control} {\bf
31},~79--91.

%\harvarditem[Horn {\em et al.}]{Horn {\em et al.}}{1996}{horn:96}
%Horn, I.G., J.R. Arulandu, J.~Gombas, J.G. VanAntwerp and R.D. Braatz  (1996).
%Improved filter design in {Internal Model Control}. {\em Ind. Eng. Chem.
%Res.} {\bf 35}(10),~3437--3441.

\harvarditem[Rivera {\em et al.}]{Rivera {\em et al.}}{1986}{rivera:86}
Rivera, D.E., M.~Morari and S.~Skogestad  (1986). Internal model control. 4.
{PID} controller design. {\em Ind. Eng. Chem. Res.} {\bf 25}(1),~252--265.

\harvarditem[Seborg {\em et al.}]{Seborg {\em et al.}}{1989}{seborg:89}
Seborg, D.E., T.F. Edgar and D.A. Mellichamp  (1989). {\em Process Dynamics and
Control}. John Wiley \& Sons.

\harvarditem[Shamsuzzoha and Skogestad]{Shamsuzzoha and
Skogestad}{2010}{shams:10}
Shamsuzzoha, M. and S.~Skogestad  (2010). The setpoint overshoot method: A
simple and fast method for closed-loop {PID} tuning. {\em J. Process Control}
{\bf 20},~1220--1234.

\harvarditem[Skogestad]{Skogestad}{2003}{skogestad:03}
Skogestad, S.  (2003). Simple analytic rules for model reduction and {PID}
controller tuning. {\em J. Process Control} {\bf 13},~291--309.

\harvarditem[Skogestad]{Skogestad}{2006}{skogestad:06}
Skogestad, S.  (2006). Tuning for smooth {PID} control with acceptable
disturbance rejection. {\em Ind. Eng. Chem. Res.} {\bf 45},~7817--7822.
%
%\harvarditem[Smith and Corripio]{Smith and Corripio}{1985}{smith:85}
%Smith, C.A. and A.B. Corripio  (1985). {\em Principles and Practice of
%Automatic Process Control}. John Wiley \& Sons.

\harvarditem[Smith]{Smith}{1957}{smith:57}
Smith, O.J.  (1957). Closer control of loops with dead time. {\em Chem. Eng.
Prog.} {\bf 53},~217.

\harvarditem[Tyreus and Luyben]{Tyreus and Luyben}{1992}{tyreus:92}
Tyreus, B.D. and W.L. Luyben  (1992). Tuning {PI} controllers for
integrator/dead time processes. {\em Ind. Eng. Chem. Res.} pp.~2628--2631.

\harvarditem[Yuwana and Seborg]{Yuwana and Seborg}{1982}{yuwana:82}
Yuwana, M. and D.E. Seborg  (1982). A new method for on-line controller tuning.
{\em AIChE Journal} {\bf 28},~434--440.

\harvarditem[Ziegler and Nichols]{Ziegler and Nichols}{1942}{ziegler:42}
Ziegler, J.G. and N.B. Nichols  (1942). Optimum settings for automatic
controllers. {\em Trans. of the A.S.M.E.} {\bf 64},~759--768.




\end{thebibliography}


\appendix

\section{Derivation of the Relationship between $\tau_c$ and $GM$, $PM$ and, $\Delta \theta /\theta$}
\label{derTcGmPmDm}
This derivation applies to processes with a relatively small time constant, or specifically for processes where $\tau_1<4(\tau_c+\theta)$, such that we according the
SIMC-rule have $\tau_I=\tau_1$.
Consider a fist order model,
$$
g=\frac{k}{\tau_1s+1}\exp(-\theta s) 
$$
and a PI controller
$$
c = K_c \(\frac{\tau_Is+1}{\tau_Is}\)
$$
with the SIMC PI-settings
$$
K_c =\frac{1}{k}\frac{\tau_1}{\tau_c+\theta};~\tau_I = \tau_1
$$
This gives the loop transfer function
$$
gc=\frac{1}{(\tau_c+\theta)s}\exp(-\theta s)
$$
The magnitude and phase angles are
$$
|gc(j\omega)|=\frac{1}{(\tau_c+\theta )\omega}
$$
$$
\angle gc(j\omega) = -\frac{\pi}{2}-\omega\theta 
$$
The gain crossover frequency is
$$|gc(j\omega_c)|=\frac{1}{(\tau_c+\theta )\omega_c}=1;~~\omega_c=\frac{1}{\tau_c+\theta}$$
and phase crossover frequency is
$$\angle gc(j\omega_{180}) = -\frac{\pi}{2}-\omega_{180}\theta  = -\pi;~~\omega_{180}=\frac{\pi}{2}\cdot \frac{1}{\theta}$$
The gain margin is then
$$
GM = \frac{1}{|gc(j\omega_{180})|}=\frac{\pi}{2} \(\frac{\tau_c}{\theta}+1\)
$$
The phase margin is
$$
PM = \angle gc(j\omega_c) + \pi= \frac{\pi}{2}-\frac{\theta }{\tau_c+\theta}
$$ 
and the allowable time delay error is
$$
\frac{\Delta\theta}{\theta}=\frac{PM}{\omega_c}\cdot\frac{1}{\theta}=\frac{\pi}{2}\( \frac{\tau_c}{\theta}+1\)-1=GM-1
$$

\section{Estimation of parameters $\tau_1$ and $\theta$ from closed-loop step response.}

\cite{shams:10} discuss at the end of their paper a two-step closed-loop procedure, 
where the first step is to use closed-loop data and some expressions to obtain the parameters $k$, $\tau_1$ and $\theta$.
We use this approach but have modified the expressions.
Our expression for $k$ in (\ref{k-shams}) is given  by their equation (35) by noting that $B = |(1-b)/b|$ where 
$b=\Delta y_{\infty}/\Delta y_s$. 
However, our expressions for $\theta$ and $\tau_1$ in (\ref{theta-chriss})-(\ref{tau1-chriss})
differ somewhat from their equations (36) and (37). 
The reason is that their equations (36) and (37) are not consistent in terms of the time delay estimate, because
the expression for $\tau_1$ in (36) is based on  $\theta = 0.43 t_p$, whereas (37) uses $\theta = 0.305 t_p$.
To correct for this, we first note from (19) in their paper (noting that $\tau_1=\tau_I$ for the delay-dominant case),
that $\tau_1$ and $\theta$ are related by
	     $$ \tau_1 = r \theta$$
where $r = 2A/B$, which is our expression in (\ref{tau1-chriss}). Here, \cite{shams:10} recommend to use 
$\theta=0.44 t_p$ for $\tau_1 < 8 \theta$ and 
$\theta=0.305 t_p$ for $\tau_1 > 8 \theta$.
However, to get better accuracy and a smooth transition, we fitted simulation data for $\theta/t_p$ as a function
of $\tau_1/\theta$ for a wide range of processes with an overshoot of 0.3, and obtained the correlation \citep{grimholt:10}
$$ \theta = t_p\cdot (0.309 + 0.209e^{-0.61(\tau_1/\theta)}) $$
as given in (\ref{theta-chriss}). 
Note here that $(0.309 + 0.209e^{-0.61(\tau_1/\theta)})$ is 0.518 for $r=\tau_1/\theta=0$,  and 0.309 for $r=\infty$.



\end{document}
