Control System Toolbox 

This example shows the limitations of PI control for processes with long dead time and illustrates the benefits of a control strategy called "Smith Predictor."
The example is inspired by:
A. Ingimundarson and T. Hagglund, "Robust Tuning Procedures of DeadTime Compensating Controllers," Control Engineering Practice, 9, 2001, pp. 11951208.
On this page… 

The process openloop response is modeled as a firstorder plus dead time with a 40.2 second time constant and 93.9 second time delay:
s = tf('s'); P = exp(93.9*s) * 5.6/(40.2*s+1); P.InputName = 'u'; P.OutputName = 'y'; P
P = From input "u" to output "y": 5.6 exp(93.9*s) *  40.2 s + 1 Continuoustime transfer function.
Note that the delay is more than twice the time constant. This model is representative of many chemical processes. Its step response is shown below.
step(P), grid on
ProportionalIntegral (PI) control is a commonly used technique in Process Control. The corresponding control architecture is shown below.
Compensator C is a PI controller in standard form with two tuning parameters: proportional gain Kp and an integral time Ti. We use the PIDTUNE command to design a PI controller with the open loop bandwidth at 0.006 rad/s:
Cpi = pidtune(P,pidstd(1,1),0.006); Cpi
Cpi = 1 1 Kp * (1 +  * ) Ti s with Kp = 0.0501, Ti = 47.3 Continuoustime PI controller in standard form
To evaluate the performance of the PI controller, close the feedback loop and simulate the responses to step changes in the reference signal ysp and output disturbance signal d. Because of the delay in the feedback path, it is necessary to convert P or Cpi to the statespace representation using the SS command:
Tpi = feedback([P*Cpi,1],1,1,1); % closedloop model [ysp;d]>y Tpi.InputName = {'ysp' 'd'}; step(Tpi), grid on
The closedloop response has acceptable overshoot but is somewhat sluggish (it settles in about 600 seconds). Increasing the proportional gain Kp speeds up the response but also significantly increases overshoot and quickly leads to instability:
Kp3 = [0.06;0.08;0.1]; % try three increasing values of Kp Ti3 = repmat(Cpi.Ti,3,1); % Ti remains the same C3 = pidstd(Kp3,Ti3); % corresponding three PI controllers T3 = feedback(P*C3,1); T3.InputName = 'ysp'; step(T3) title('Loss of stability when increasing Kp')
The performance of the PI controller is severely limited by the long dead time. This is because the PI controller has no knowledge of the dead time and reacts too "impatiently" when the actual output y does not match the desired setpoint ysp. Everyone has experienced a similar phenomenon in showers where the water temperature takes a long time to adjust. There, impatience typically leads to alternate scolding by burning hot and freezing cold water. A better strategy consists of waiting for a change in temperature setting to take effect before making further adjustments. And once we have learned what knob setting delivers our favorite temperature, we can get the right temperature in just the time it takes the shower to react. This "optimal" control strategy is the basic idea behind the Smith Predictor scheme.
The Smith Predictor control structure is sketched below.
The Smith Predictor uses an internal model Gp to predict the delayfree response yp of the process (e.g., what water temperature a given knob setting will deliver). It then compares this prediction yp with the desired setpoint ysp to decide what adjustments are needed (control u). To prevent drifting and reject external disturbances, the Smith predictor also compares the actual process output with a prediction y1 that takes the dead time into account. The gap dy=yy1 is fed back through a filter F and contributes to the overall error signal e. Note that dy amounts to the perceived temperature mismatch after waiting long enough for the shower to react.
Deploying the Smith Predictor scheme requires
A model Gp of the process dynamics and an estimate tau of the process dead time
Adequate settings for the compensator and filter dynamics (C and F)
Based on the process model, we use:
For F, use a firstorder filter with a 20 second time constant to capture lowfrequency disturbances.
F = 1/(20*s+1); F.InputName = 'dy'; F.OutputName = 'dp';
For C, we redesign the PI controller with the overall plant seen by the PI controller, which includes dynamics from P, Gp, F and dead time. With the help of the Smith Predictor control structure we are able to increase the openloop bandwidth to achieve faster response and increase the phase margin to reduce the overshoot.
% Process P = exp(93.9*s) * 5.6/(40.2*s+1); P.InputName = 'u'; P.OutputName = 'y0'; % Prediction model Gp = 5.6/(40.2*s+1); Gp.InputName = 'u'; Gp.OutputName = 'yp'; Dp = exp(93.9*s); Dp.InputName = 'yp'; Dp.OutputName = 'y1'; % Overall plant S1 = sumblk('ym = yp + dp'); S2 = sumblk('dy = y0  y1'); Plant = connect(P,Gp,Dp,F,S1,S2,'u','ym'); % Design PI controller with 0.08 rad/s bandwidth and 90 degrees phase margin Options = pidtuneOptions('PhaseMargin',90); C = pidtune(Plant,pidstd(1,1),0.08,Options); C.InputName = 'e'; C.OutputName = 'u'; C
C = 1 1 Kp * (1 +  * ) Ti s with Kp = 0.574, Ti = 40.2 Continuoustime PI controller in standard form
Comparison of PI Controller vs. Smith Predictor
To compare the performance of the two designs, first derive the closedloop transfer function from ysp,d to y for the Smith Predictor architecture. To facilitate the task of connecting all the blocks involved, name all their input and output channels and let CONNECT do the wiring:
% Assemble closedloop model from [y_sp,d] to y Sum1 = sumblk('e = ysp  yp  dp'); Sum2 = sumblk('y = y0 + d'); Sum3 = sumblk('dy = y  y1'); T = connect(P,Gp,Dp,C,F,Sum1,Sum2,Sum3,{'ysp','d'},'y');
Use STEP to compare the Smith Predictor (blue) with the PI controller (red):
step(T,'b',Tpi,'r') grid on legend('Smith Predictor','PI Controller')
The Smith Predictor provides much faster response with no overshoot. The difference is also visible in the frequency domain by plotting the closedloop Bode response from ysp to y. Note the higher bandwidth for the Smith Predictor.
bode(T(1,1),'b',Tpi(1,1),'r',{1e3,1}) grid on legend('Smith Predictor','PI Controller')
In the previous analysis, the internal model
matched the process model P exactly. In practical situations, the internal model is only an approximation of the true process dynamics, so it is important to understand how robust the Smith Predictor is to uncertainty on the process dynamics and dead time.
Consider two perturbed plant models representative of the range of uncertainty on the process parameters:
P1 = exp(90*s) * 5/(38*s+1); P2 = exp(100*s) * 6/(42*s+1); bode(P,P1,P2), grid on title('Nominal and Perturbed Process Models')
To analyze robustness, collect the nominal and perturbed models into an array of process models, rebuild the closedloop transfer functions for the PI and Smith Predictor designs, and simulate the closedloop responses:
Plants = stack(1,P,P1,P2); % array of process models T1 = connect(Plants,Gp,Dp,C,F,Sum1,Sum2,Sum3,{'ysp','d'},'y'); % Smith Tpi = feedback([Plants*Cpi,1],1,1,1); % PI step(T1,'b',Tpi,'r') grid on legend('Smith Predictor 1','PI Controller')
Both designs are sensitive to model mismatch, as confirmed by the closedloop Bode plots:
bode(T1(1,1),'b',Tpi(1,1),'r') grid on legend('Smith Predictor 1','PI Controller')
To reduce the Smith Predictor's sensitivity to modeling errors, check the stability margins for the inner and outer loops. The inner loop C has openloop transfer C*Gp so the stability margin are obtained by
margin(C * Gp)
title('Stability Margins for the Inner Loop (C)')
The inner loop has comfortable gain and phase margins so focus on the outer loop next. Use CONNECT to derive the openloop transfer function L from ysp to dp with the inner loop closed:
Sum1o = sumblk('e = ysp  yp'); % open the loop at dp L = connect(P,Gp,Dp,C,F,Sum1o,Sum2,Sum3,{'ysp','d'},'dp'); bodemag(L(1,1))
This transfer function is essentially zero, which is to be expected when the process and prediction models match exactly. To get insight into the stability margins for the outer loop, we need to work with one of the perturbed process models, e.g., P1:
H = connect(Plants(:,:,2),Gp,Dp,C,Sum1o,Sum2,Sum3,{'ysp','d'},'dy'); H = H(1,1); % openloop transfer ysp > dy L = F * H; margin(L) title('Stability Margins for the Outer Loop (F)') grid on, set(gca,'xlim',[1e2 1])
This gain curve has a hump near 0.04 rad/s that lowers the gain margin and increases the hump in the closedloop step response. To fix this issue, pick a filter F that rolls off earlier and more quickly:
F = (1+10*s)/(1+100*s); F.InputName = 'dy'; F.OutputName = 'dp';
Verify that the gain margin has improved near the 0.04 rad/s phase crossing:
L = F * H; margin(L) title('Stability Margins for the Outer Loop with Modified F') grid on, set(gca,'xlim',[1e2 1])
Finally, simulate the closedloop responses with the modified filter:
T2 = connect(Plants,Gp,Dp,C,F,Sum1,Sum2,Sum3,{'ysp','d'},'y'); step(T2,'b',Tpi,'r') grid on legend('Smith Predictor 2','PI Controller')
The modified design provides more consistent performance at the expense of a slightly slower nominal response.
Improving Disturbance Rejection
Formulas for the closedloop transfer function from d to y show that the optimal choice for F is
where tau is the internal model's dead time. This choice achieves perfect disturbance rejection regardless of the mismatch between P and Gp. Unfortunately, such "negative delay" is not causal and cannot be implemented. In the paper:
Huang, H.P., et al., "A Modified Smith Predictor with an Approximate Inverse of Dead Time," AiChE Journal, 36 (1990), pp. 10251031
the authors suggest using the phase lead approximation:
where B is a lowpass filter with the same time constant as the internal model Gp. You can test this scheme as follows:
Define B(s) and F(s)
B = 0.05/(40*s+1); tau = totaldelay(Dp); F = (1+B)/(1+B*exp(tau*s)); F.InputName = 'dy'; F.OutputName = 'dp';
Redesign PI controller with reduced bandwidth
Plant = connect(P,Gp,Dp,F,S1,S2,'u','ym'); C = pidtune(Plant,pidstd(1,1),0.02,pidtuneOptions('PhaseMargin',90)); C.InputName = 'e'; C.OutputName = 'u'; C
C = 1 1 Kp * (1 +  * ) Ti s with Kp = 0.144, Ti = 40.1 Continuoustime PI controller in standard form
Computed closedloop model T3
T3 = connect(Plants,Gp,Dp,C,F,Sum1,Sum2,Sum3,{'ysp','d'},'y');
Compare T3 with T2 and Tpi
step(T2,'b',T3,'g',Tpi,'r') grid on legend('Smith Predictor 2','Smith Predictor 3','PI Controller')
This comparison shows that our last design speeds up disturbance rejection at the expense of slower setpoint tracking.