Accelerating the pace of engineering and science

# 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.

```open_system('rct_CSTR_OL')
```

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,...
findopOptions('DisplayReport','off'));
% 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);
end
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')
subplot(312), plot(TrEQ,'b-*'), grid, title('Reactor temperature')
subplot(313), plot(TcEQ,'b-*'), grid, title('Coolant temperature')
```

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).

```pole(G)
```
```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 a feedback controller to adjust the coolant temperature Tc based on measurements of the residual concentration Cr and reactor temperature Tr. This controller is discrete with a sampling period of half a minute.

```open_system('rct_CSTR')
```

The concentration setpoint Cref is a ramp from 8.57 kmol/m^3 to 2 kmol/m^3 in 26 minutes. For the reactor temperature setpoint Tref, use the equilibrium values TrEQ derived above. This means that the controller's main task is to keep the reaction close to equilibrium throughout the transition. The controller output is the coolant temperature adjustment dTc relative to the initial equilibrium value TcEQ(1)=297.98 for Cr=8.57.

```clf
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('Concentration setpoint'), ylabel('Cref')
subplot(212), plot(t,interp1(CrEQ,TrEQ,C))
title('Temperature setpoint'), ylabel('Tref'), grid
```

Plot the Bode response of the five linearized models in G. The reaction dynamics vary substantially with concentration, suggesting that the controller gains should be scheduled as a function of the residual concentration Cr.

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

Control Objectives

Use TuningGoal objects to capture the design requirements. First, Cr should track the setpoint Cref with a response time of about 10 minutes.

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

The closed-loop dynamics should be sufficiently damped and decay reasonably fast.

```R2 = TuningGoal.Poles();
R2.MinDamping = 0.3;
R2.MinDecay = 0.1;
```

The open-loop response at the controller output du should have integral action, a bandwidth of at least 2/10=0.2 rad/min to achieve a 10 min response time, -20 dB/decade roll-off past 2 rad/min, and gain and phase margins of at least 7 dB and 45 degrees.

```R3 = TuningGoal.MinLoopGain('du',0.2,1);
R4 = TuningGoal.MaxLoopGain('du',2,1);
R5 = TuningGoal.Margins('du',7,45);
```

The Rate Limit block at the controller output indicates that the coolant temperature Tc cannot vary faster than plus or minus 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 while Tref varies at a rate of about 2.5 degrees/min. To ensure that du does not vary faster than 10 degrees/min, the gain from Cref to du should be less than 10/0.25=40, and the gain from Tref to du should be less than 10/2.5=4. Use a WeightedGain requirement to capture these two constraints.

```R6 = TuningGoal.WeightedGain({'Cref','Tref'},'du',1,diag([1/40 1/4]));
```

Gain-Scheduled State-Space Controller

Since the target bandwidth is close to the Nyquist frequency, it is more convenient to tune the controller in the discrete domain. Discretize the linearized process dynamics with sample time of half a minute. Use the Tustin method prewarped at 1 rad/min to minimize phase distortion near crossover.

```Ts = 0.5;  % in minutes
Gd = c2d(G,Ts,'tustin',1);
```

At first, it is unclear what controller structure is best suited for this application. To gain insight into the right structure, start with a generic second-order space-space controller whose matrices are functions of Cr. Try a simple affine dependence on Cr:

Using systune, you can automatically tune the matrix coefficients to meet the requirements R1-R6 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 build tunable models of the matrices at the five concentration values CrEQ.

```A = gainsurf('A',zeros(2),CrEQ);
B = gainsurf('B',eye(2),CrEQ);
C = gainsurf('C',[0 0],CrEQ);
D = gainsurf('D',[0 0],CrEQ);
```

Combine A,B,C,D into a tunable model of the gain-scheduled controller. Include a loopswitch block at the controller output du to gain access to the open-loop responses at this location.

```CC1 = loopswitch('du',1) * ss(A,B,C,D,Ts);
CC1.InputName = {'ECr','ETr'};
CC1.OutputName = 'Tc';
CC1.TimeUnit = 'min';
```

Use connect to build a closed-loop model of the overall control system at the five design points. Here T0 is a 5-by-1 array of linear models depending on the tunable coefficients . Each model is discrete sampled every half minute.

```S1 = sumblk('ECr = Cref - Cr');
S2 = sumblk('ETr = Tref - Tr');
T0 = connect(Gd,CC1,S1,S2,{'Cref','Tref','Cf','Tf'},{'Cr','Tr'})
```
```T0 =

5x1 array of generalized discrete-time state-space models.
Each model has 2 outputs, 4 inputs, 4 states, and the following blocks:
A_0: Parametric 2x2 matrix, 1 occurrences.
A_1: Parametric 2x2 matrix, 1 occurrences.
B_0: Parametric 2x2 matrix, 1 occurrences.
B_1: Parametric 2x2 matrix, 1 occurrences.
C_0: Parametric 1x2 matrix, 1 occurrences.
C_1: Parametric 1x2 matrix, 1 occurrences.
D_0: Parametric 1x2 matrix, 1 occurrences.
D_1: Parametric 1x2 matrix, 1 occurrences.
du: Open/closed switch, 1 channels, 1 occurrences.

Type "ss(T0)" to see the current value, "get(T0)" to see all properties, and "T0.Blocks" to interact with the blocks.

```

You can now use systune to tune the gain-scheduled controller coefficients against the requirements R1-R6. Use the stability margin, roll-off, and rate limit requirements as hard constraints and optimize the remaining requirements.

```T1 = systune(T0,[R1 R2 R3],[R4 R5 R6]);
```
```Final: Soft = 1.05, Hard = 0.99957, Iterations = 522
```

The resulting design satisfies the hard constraints (Hard<1) and almost satisfies the remaining requirements (Soft close to 1). To validate this design, simulate the responses to ramps in concentration and temperature with the same slopes as Cref and Tref. 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);
uT = interp1([0 2 5 20],2.38*[0 0 3 3],t);
subplot(211), lsim(getIOTransfer(T1,{'Cref','Tref'},'Cr'),[uC;uT])
grid, set(gca,'ylim',[-1.5 0.5]), title('Residual concentration')
subplot(212), lsim(getIOTransfer(T1,{'Cref','Tref'},'du'),[uC;uT])
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). Verify the open-loop crossover and roll-off at each design point.

```clf, bodemag(getLoopTransfer(T1,'du',-1)), grid
set(gca,'ylim',[-30 60])
```

To simulate the nonlinear response, use Gain and Product blocks to implement the formulas for Note that the controller matrices are updated at each sampling period based on the measured value of Cr. This ensures that the controller varies smoothly as a function of the scheduling variable Cr.

Figure 1: Gain-scheduled state-space controller.

Figure 2: Scheduling of the A matrix as a function of Cr.

Use gainsurfdata to extract the tuned coefficients and simulate the closed-loop response. The simulation results appear below. The gain-scheduled controller successfully drives the reaction through the transition with minimal saturation of the rate limits (controller output du is plotted in yellow, effective temperature variation dTc in magenta).

```[A0,A1] = gainsurfdata(setBlockValue(A,T1));
[B0,B1] = gainsurfdata(setBlockValue(B,T1));
[C0,C1] = gainsurfdata(setBlockValue(C,T1));
[D0,D1] = gainsurfdata(setBlockValue(D,T1));
```

Figure 3: Transition with gain-scheduled state-space controller.

Gain-Scheduled PID Controller

For implementation purposes, it is often desirable to simplify the controller and elucidate what it actually does. To gain insight into the state-space controller inner workings, plot its Bode response at the five design points.

```CC = setBlockValue(CC1,T1);
bode(CC,{1e-3,10}), grid
```

The first channel resembles a PID controller and the second channel is close to a pure gain. This suggests using PID control for the concentration Cr and proportional feedback for the temperature Tr.

Figure 4: Simplified control structure.

As before, use low-order polynomials to model the dependence of the controller gains Kp,Ki,Kd,... on the scheduling variable Cr. First-order polynomials are a natural starting point, but here the best results are obtained for quadratic polynomials in Cr, for example,

Use gainsurf to model these quadratic gain schedules.

```% PID gains for concentration feedback
Kp = gainsurf('Kp',0,CrEQ,CrEQ.^2);
Ki = gainsurf('Ki',-0.1,CrEQ,CrEQ.^2);
Kd = gainsurf('Kd',0,CrEQ,CrEQ.^2);
N = gainsurf('N',10,CrEQ,CrEQ.^2);

% Static gain for temperature feedback
Kt = gainsurf('Kt',0,CrEQ,CrEQ.^2);
```

Combine these gains into the controller of Figure 4 and build a tunable model T0 of the gain-scheduled control system.

```CC2 = loopswitch('du',1) * [pid(Kp,Ki,Kd,1/N,Ts) , Kt];
CC2.InputName = {'ECr','ETr'};
CC2.OutputName = 'Tc';
CC2.TimeUnit = 'min';

S1 = sumblk('ECr = Cref - Cr');
S2 = sumblk('ETr = Tref - Tr');
T0 = connect(Gd,CC2,S1,S2,{'Cref','Tref','Cf','Tf'},{'Cr','Tr'});
```

Use systune to tune the coefficients of

```T2 = systune(T0,[R1 R2 R3],[R4 R5 R6]);
```
```Final: Soft = 1.04, Hard = 0.99966, Iterations = 171
```

The optimization results are similar to those obtained for the state-space controller. Verify the open-loop responses at the five design points.

```clf, bodemag(getLoopTransfer(T2,'du',-1)), grid
set(gca,'ylim',[-30 60])
```

To validate this second design in the Simulink model, extract the tuned coefficients and switch to the PID variant of the controller. This variant implements the quadratic gain formulas as shown in Figure 5. The simulation results appear in Figure 6 and are qualitatively similar to those obtained earlier.

```[Kp_0,Kp_1,Kp_2] = gainsurfdata(setBlockValue(Kp,T2));
[Ki_0,Ki_1,Ki_2] = gainsurfdata(setBlockValue(Ki,T2));
[Kd_0,Kd_1,Kd_2] = gainsurfdata(setBlockValue(Kd,T2));
[N_0,N_1,N_2] = gainsurfdata(setBlockValue(N,T2));
[Kt_0,Kt_1,Kt_2] = gainsurfdata(setBlockValue(Kt,T2));

% Switch to PID controller variant of "Controller" block
MODE = 2;
```

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

Figure 6: Transition with gain-scheduled PID controller.

Inspect how each gain varies with Cr during the transition.

```SG = G.SamplingGrid;
subplot(321), Kp.SamplingGrid = SG; view(setBlockValue(Kp,T2)), ylabel('Kp')
subplot(322), Ki.SamplingGrid = SG; view(setBlockValue(Ki,T2)), ylabel('Ki')
subplot(323), Kd.SamplingGrid = SG; view(setBlockValue(Kd,T2)), ylabel('Kd')
subplot(324), N.SamplingGrid = SG;  view(setBlockValue(N,T2)), ylabel('N')
subplot(325), Kt.SamplingGrid = SG; view(setBlockValue(Kt,T2)), ylabel('Kt')
```

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