Robust Control Toolbox

Gain-Scheduled Control of a Chemical Reactor

This example shows how to design and tune a gain-scheduled controller for a chemical reactor transitioning from low to high conversion rate. For background, see Seborg, D.E. et al., "Process Dynamics and Control", 2nd Ed., 2004, Wiley, pp. 34-36.

Continuous Stirred Tank Reactor

The process considered here is a continuous stirred tank reactor (CSTR) during transition from low to high conversion rate (high to low residual concentration). Because the chemical reaction is exothermic (produces heat), the reactor temperature must be controlled to prevent a thermal runaway. The control task is complicated by the fact that the process dynamics are nonlinear and transition from stable to unstable and back to stable as the conversion rate increases. The reactor dynamics are modeled in Simulink. The controlled variables are the residual concentration Cr and the reactor temperature Tr, and the manipulated variable is the temperature Tc of the coolant circulating in the reactor's cooling jacket.


We want to transition from a residual concentration of 8.57 kmol/m^3 initially down to 2 kmol/m^3. To understand how the process dynamics evolve with the residual concentration Cr, find the equilibrium conditions for five values of Cr between 8.57 and 2 and linearize the process dynamics around each equilibrium. Log the reactor and coolant temperatures at each equilibrium point.

clear G
CrEQ = linspace(8.57,2,5)';  % concentrations
TrEQ = zeros(5,1);           % reactor temperatures
TcEQ = zeros(5,1);           % coolant temperatures
for k=1:5
   opspec = operspec('rct_CSTR_OL');
   % Set desired residual concentration
   opspec.Outputs(1).y = CrEQ(k);
   opspec.Outputs(1).Known = true;
   % Compute equilibrium condition
   [op,report] = findop('rct_CSTR_OL',opspec,...
   % Log temperatures
   TrEQ(k) = report.Outputs(2).y;
   TcEQ(k) = op.Inputs.u;
   % Linearize process dynamics
   G(:,:,k) = linearize('rct_CSTR_OL', 'rct_CSTR_OL/CSTR', op);
G.InputName = {'Cf','Tf','Tc'};
G.OutputName = {'Cr','Tr'};
G.SamplingGrid = struct('Cr',CrEQ);
G.TimeUnit = 'min';

Plot the reactor and coolant temperatures at equilibrium as a function of concentration.

subplot(311), plot(CrEQ,'b-*'), grid, title('Residual concentration'), ylabel('CrEQ')
subplot(312), plot(TrEQ,'b-*'), grid, title('Reactor temperature'), ylabel('TrEQ')
subplot(313), plot(TcEQ,'b-*'), grid, title('Coolant temperature'), ylabel('TcEQ')

An open-loop control strategy consists of following the coolant temperature profile above to smoothly transition between the Cr=8.57 and Cr=2 equilibria. However, this strategy is doomed by the fact that the reaction is unstable in the mid range and must be properly cooled to avoid thermal runaway. This is confirmed by inspecting the poles of the linearized models for the five equilibrium points considered above (three out of the five models are unstable).

ans(:,:,1) =

  -0.5225 + 0.0000i
  -0.8952 + 0.0000i

ans(:,:,2) =

   0.1733 + 0.0000i
  -0.8866 + 0.0000i

ans(:,:,3) =

   0.5114 + 0.0000i
  -0.8229 + 0.0000i

ans(:,:,4) =

   0.0453 + 0.0000i
  -0.4991 + 0.0000i

ans(:,:,5) =

  -1.1077 + 1.0901i
  -1.1077 - 1.0901i

Feedback Control Strategy

To prevent thermal runaway while ramping down the residual concentration, use feedback control to adjust the coolant temperature Tc based on measurements of the residual concentration Cr and reactor temperature Tr. For this application, we use a cascade control architecture where the inner loop regulates the reactor temperature and the outer loop tracks the concentration setpoint. Both feedback loops are digital with a sampling period of 0.5 minutes.


The target concentration Cref ramps down from 8.57 kmol/m^3 at t=10 to 2 kmol/m^3 at t=36 (the transition lasts 26 minutes). The corresponding profile Tref for the reactor temperature is obtained by interpolating the equilibrium values TrEQ from trim analysis. The controller computes the coolant temperature adjustment dTc relative to the initial equilibrium value TcEQ(1)=297.98 for Cr=8.57. Note that the model is set up so that initially, the output TrSP of the "Concentration controller" block matches the reactor temperature, the adjustment dTc is zero, and the coolant temperature Tc is at its equilibrium value TcEQ(1).

t = [0 10:36 45];
C = interp1([0 10 36 45],[8.57 8.57 2 2],t);
subplot(211), plot(t,C), grid, set(gca,'ylim',[0 10])
title('Target residual concentration'), ylabel('Cref')
subplot(212), plot(t,interp1(CrEQ,TrEQ,C))
title('Corresponding reactor temperature at equilibrium'), ylabel('Tref'), grid

Control Objectives

Use TuningGoal objects to capture the design requirements. First, Cr should follow setpoints Cref with a response time of about 5 minutes.

R1 = TuningGoal.Tracking('Cref','Cr',5);

The inner loop (temperature) should stabilize the reaction dynamics with sufficient damping and fast enough decay.

MinDecay = 0.2;
MinDamping = 0.5;
% Constrain closed-loop poles of inner loop with the outer loop open
R2 = TuningGoal.Poles('Tc',MinDecay,MinDamping);
R2.Openings = 'TrSP';

The Rate Limit block at the controller output specifies that the coolant temperature Tc cannot vary faster than 10 degrees per minute. This is a severe limitation on the controller authority which, when ignored, can lead to poor performance or instability. To take this rate limit into account, observe that Cref varies at a rate of 0.25 kmol/m^3/min. To ensure that Tc does not vary faster than 10 degrees/min, the gain from Cref to Tc should be less than 10/0.25=40.

R3 = TuningGoal.Gain('Cref','Tc',40);

Finally, require at least 7 dB of gain margin and 45 degrees of phase margin at the plant input Tc.

R4 = TuningGoal.Margins('Tc',7,45);

To achieve these requirements, we use a PI controller in the outer loop (see Figure 1) and a lead compensator in the inner loop (see Figure 2). Note that due to the slow sampling rate, a pure gain in the inner loop is not enough to adequately stabilize the chemical reaction at the mid-range concentration Cr = 5.28 kmol/m^3/min.

Figure 1: PI controller for concentration loop.

Figure 2: Lead compensator for temperature loop.

Plot the Bode response of the linearized models obtained for 5 concentration values between 8.57 and 2 kmol/m^3. The reaction dynamics vary substantially with concentration, suggesting that the PI controller and lead compensator gains should be scheduled as a function of the residual concentration Cr.

clf, bode(G(:,3),{0.01,10})

Gain-Scheduled Controller

Because the target bandwidth is within a decade of the Nyquist frequency, it is easier to tune the controller directly in the discrete domain. Discretize the linearized process dynamics with sample time of 0.5 minutes. Use the ZOH method to reflect how the digital controller interacts with the continuous-time plant.

Ts = 0.5;  % in minutes
Gd = c2d(G,Ts);

Use low-order polynomials to model the dependence of the controller gains Kp,Ki,Kt,a,b on the scheduling variable Cr. First-order polynomials are a natural starting point, but here better results are obtained with quadratic polynomials in Cr, for example,

$$ K_p(C_r) = K_{p0} + K_{p1} C_r + K_{p2} C_r^2 . $$

Using systune, you can automatically tune the coefficients $K_{p0}, K_{p1}, K_{p2}, K_{i0}, \ldots$ to meet the requirements R1-R4 at the five equilibrium points computed above. This amounts to tuning the gain-scheduled controller at five design points along the Cref trajectory. Use gainsurf to model the quadratic gain schedules $K_p(C_r),K_i(C_r),\ldots$ at the five concentration values CrEQ.

Kp = gainsurf('Kp',0,CrEQ,CrEQ.^2);
Ki = gainsurf('Ki',-1,CrEQ,CrEQ.^2);
Kt = gainsurf('Kt',0,CrEQ,CrEQ.^2);
a = gainsurf('a',0,CrEQ,CrEQ.^2);
b = gainsurf('b',0,CrEQ,CrEQ.^2);

Use these gains to build the PI and lead controllers.

PI = pid(Kp,Ki,'Ts',Ts,'TimeUnit','min');
PI.u = 'ECr';   PI.y = 'TrSP';

LEAD = Kt * tf([1 -a],[1 -b],Ts,'TimeUnit','min');
LEAD.u = 'ETr';   LEAD.y = 'Tc';

Use connect to build a closed-loop model of the overall control system at the five design points. Mark the controller outputs TrSP and Tc as "analysis points" so that loops can be opened and stability margins evaluated at these locations. The closed-loop model T0 is a 5-by-1 array of linear models depending on the tunable coefficients $K_{p0}, K_{p1}, K_{p2}, K_{i0}, \ldots$. Each model is discrete and sampled every half minute.

S1 = sumblk('ECr = Cref - Cr');
S2 = sumblk('ETr = TrSP - Tr');
T0 = connect(Gd(:,'Tc'),LEAD,PI,S1,S2,'Cref','Cr',{'TrSP','Tc'});

You can now use systune to tune the controller coefficients against the requirements R1-R4. Make the stability margin requirement a hard constraints and optimize the remaining requirements.

T = systune(T0,[R1 R2 R3],R4);
Final: Soft = 1.21, Hard = 0.99996, Iterations = 289

The resulting design satisfies the hard constraint (Hard<1) and nearly satisfies the remaining requirements (Soft close to 1). To validate this design, simulate the responses to a ramp in concentration with the same slope as Cref. Each plot shows the linear responses at the five design points CrEQ.

t = 0:Ts:20;
uC = interp1([0 2 5 20],(-0.25)*[0 0 3 3],t);
subplot(211), lsim(getIOTransfer(T,'Cref','Cr'),uC)
grid, set(gca,'ylim',[-1.5 0.5]), title('Residual concentration')
subplot(212), lsim(getIOTransfer(T,'Cref','Tc'),uC)
grid, title('Coolant temperature variation')

Note that rate of change of the coolant temperature remains within the physical limits (10 degrees per minute or 5 degrees per sample period).

Controller Validation

Each controller gain in the Simulink model is set up as a quadratic function of the scheduling variable Cr. For example, Figure 3 shows how the Kp block implements the formula

$$ K_p(C_r) = K_{p0} + K_{p1} C_r + K_{p2} C_r^2 . $$

The controller gains are updated at each sampling period based on the measured value of Cr, which ensures that the control signal varies smoothly as a function of Cr.

Figure 3: Scheduling of the proportional gain Kp as a function of Cr.

To validate the gain-scheduled controller in Simulink, extract the tuned coefficients using gainsurfdata. Use the same variable names as in the Simulink model, for example, Kp_0,Kp_1,Kp_2 for the proportional gain $K_p$.

[Kp_0,Kp_1,Kp_2] = gainsurfdata(setBlockValue(Kp,T));
[Ki_0,Ki_1,Ki_2] = gainsurfdata(setBlockValue(Ki,T));
[Kt_0,Kt_1,Kt_2] = gainsurfdata(setBlockValue(Kt,T));
[a_0,a_1,a_2] = gainsurfdata(setBlockValue(a,T));
[b_0,b_1,b_2] = gainsurfdata(setBlockValue(b,T));

The simulation results appear in Figure 4. The gain-scheduled controller successfully drives the reaction through the transition with adequate response time and no saturation of the rate limits (controller output du matches effective temperature variation dTc). The reactor temperature stays close to its equilibrium value Tref, indicating that the controller keeps the reaction near equilibrium while preventing thermal runaway.

Figure 4: Transition with gain-scheduled cascade controller.

Controller Implementation

Inspect how each gain varies with Cr during the transition.

SG = G.SamplingGrid;
subplot(321), Kp.SamplingGrid = SG; view(setBlockValue(Kp,T)), ylabel('Kp')
subplot(322), Ki.SamplingGrid = SG; view(setBlockValue(Ki,T)), ylabel('Ki')
subplot(323), Kt.SamplingGrid = SG; view(setBlockValue(Kt,T)), ylabel('Kt')
subplot(324), a.SamplingGrid = SG; view(setBlockValue(a,T)), ylabel('a')
subplot(325), b.SamplingGrid = SG; view(setBlockValue(b,T)), ylabel('b')

You can generate code for this digital gain-scheduled controller "as is". If preferred, you can also convert the quadratic gain formulas into classic gain-scheduling lookup tables.