www.gusucode.com > control_featured 案例源码程序 matlab代码 > control_featured/GainScheduledProcessExample.m

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

% Copyright 1990-2013 The MathWorks, Inc.

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

CrEQ = linspace(8.57,2,5)';  % concentrations
TrEQ = zeros(5,1);           % reactor temperatures
TcEQ = zeros(5,1);           % coolant temperatures

% Specify trim conditions
opspec = operspec('rct_CSTR_OL',5);
for k=1:5
   % Set desired residual concentration
   opspec(k).Outputs(1).y = CrEQ(k);
   opspec(k).Outputs(1).Known = true;
end

% Compute equilibrium condition and log corresponding temperatures
[op,report] = findop('rct_CSTR_OL',opspec,...
   findopOptions('DisplayReport','off'));
for k=1:5
   TrEQ(k) = report(k).Outputs(2).y;
   TcEQ(k) = op(k).Inputs.u;
end

% Linearize process dynamics at trim conditions
G = linearize('rct_CSTR_OL', 'rct_CSTR_OL/CSTR', op); 
G.InputName = {'Cf','Tf','Tc'};
G.OutputName = {'Cr','Tr'};

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

pole(G)

%%
% The Bode plot further highlights the significant variations in plant dynamics 
% while transitioning from |Cr|=8.57 to |Cr|=2.

clf, bode(G(:,'Tc'),{0.01,10})


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

open_system('rct_CSTR')

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

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('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);

%% Gain-Scheduled Controller
% To achieve these requirements, we use a PI controller in the outer loop
% and a lead compensator in the inner loop. Due to the slow sampling rate, 
% the lead compensator is needed to adequately stabilize the chemical 
% reaction at the mid-range concentration |Cr| = 5.28 kmol/m^3/min.
% Because the reaction dynamics vary substantially with concentration, we
% further schedule the controller gains as a function of concentration.
% This is modeled in Simulink using Lookup Table blocks as shown in 
% Figures 1 and 2.
%
% <<../gsprocess1.png>>
%
% *Figure 1: Gain-scheduled PI controller for concentration loop.*
%
% <<../gsprocess2.png>>
%
% *Figure 2: Gain-scheduled lead compensator for temperature loop.*
%
% Tuning this gain-scheduled controller amounts to tuning the look-up table 
% data over a range of concentration values.
% Rather than tuning individual look-up table entries, parameterize the
% controller gains |Kp,Ki,Kt,a,b| as quadratic polynomials in |Cr|, for example,
%
% $$ K_p(C_r) = K_{p0} + K_{p1} C_r + K_{p2} C_r^2 . $$
%
% Besides reducing the number of variables to tune, this approach ensures
% smooth gain transitions as |Cr| varies. 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 the |tunableSurface| object to parameterize each gain as a quadratic
% function of |Cr|. The "tuning grid" is set to the five concentrations
% |CrEQ| and the basis functions for the quadratic 
% parameterization are $C_r, C_r^2$. Most gains
% are initialized to be identically zero.

TuningGrid = struct('Cr',CrEQ);
ShapeFcn = @(Cr) [Cr , Cr^2];

Kp = tunableSurface('Kp', 0, TuningGrid, ShapeFcn);
Ki = tunableSurface('Ki', -2, TuningGrid, ShapeFcn);
Kt = tunableSurface('Kt', 0, TuningGrid, ShapeFcn);
a = tunableSurface('a', 0, TuningGrid, ShapeFcn);
b = tunableSurface('b', 0, TuningGrid, ShapeFcn);

%% Controller Tuning
% 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;
Gd = c2d(G,Ts);

%%
% Create an |slTuner| interface for tuning the quadratic gain schedules
% introduced above. Use block substitution to replace the nonlinear plant 
% model by the five discretized linear models |Gd| obtained at the design points
% |CrEQ|. Use |setBlockParam| to associate the tunable gain functions 
% |Kp|, |Ki|, |Kt|, |a|, |b| with the Lookup Table blocks of the same name.

BlockSubs = struct('Name','rct_CSTR/CSTR','Value',Gd);
ST0 = slTuner('rct_CSTR',{'Kp','Ki','Kt','a','b'},BlockSubs);
ST0.Ts = Ts;  % sample time for tuning

% Register points of interest 
ST0.addPoint({'Cref','Cr','Tr','TrSP','Tc'})

% Parameterize look-up table blocks
ST0.setBlockParam('Kp',Kp);
ST0.setBlockParam('Ki',Ki);
ST0.setBlockParam('Kt',Kt);
ST0.setBlockParam('a',a);
ST0.setBlockParam('b',b);

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

ST = systune(ST0,[R1 R2 R3],R4);

%%
% 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(ST,'Cref','Cr'),uC)
grid, set(gca,'ylim',[-1.5 0.5]), title('Residual concentration')
subplot(212), lsim(getIOTransfer(ST,'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
% Inspect how each gain varies with |Cr| during the transition.

% Access tuned gain schedules
TGS = getBlockParam(ST);

% Plot gain profiles
clf
subplot(321), viewSurf(TGS.Kp), ylabel('Kp')
subplot(322), viewSurf(TGS.Ki), ylabel('Ki')
subplot(323), viewSurf(TGS.Kt), ylabel('Kt')
subplot(324), viewSurf(TGS.a), ylabel('a')
subplot(325), viewSurf(TGS.b), ylabel('b')

%%
% To validate the gain-scheduled controller in Simulink, first use
% |writeBlockValue| to apply the tuning results to the Simulink model.
% For each Lookup Table block, this evaluates the corresponding quadratic 
% gain formula at the table breakpoints and updates the table data
% accordingly.

writeBlockValue(ST)

%%
% Next push the Play button to simulate the response with the tuned gain
% schedules. The simulation results appear in Figure 3. 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.
%
% <<../gsprocess4.png>>
%
% *Figure 3: Transition with gain-scheduled cascade controller.*

%% Controller Tuning in MATLAB
% Alternatively, you can tune the gain schedules directly in MATLAB without 
% using the |slTuner| interface. First parameterize the gains as quadratic 
% functions of |Cr| as done above.

TuningGrid = struct('Cr',CrEQ);
ShapeFcn = @(Cr) [Cr , Cr^2];

Kp = tunableSurface('Kp', 0, TuningGrid, ShapeFcn);
Ki = tunableSurface('Ki', -2, TuningGrid, ShapeFcn);
Kt = tunableSurface('Kt', 0, TuningGrid, ShapeFcn);
a = tunableSurface('a', 0, TuningGrid, ShapeFcn);
b = tunableSurface('b', 0, TuningGrid, ShapeFcn);

%%
% 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 of |Kp,Ki,Kt,a,b|. Each model 
% is discrete and sampled every half minute.

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

%%
% Finally, use |systune| to tune the gain schedule coefficients.

T = systune(T0,[R1 R2 R3],R4);

%%
% The result is similar to the one obtained above. Confirm by plotting
% the gains as a function of |Cr| using the tuned coefficients in |T|.

clf
subplot(321), viewSurf(setBlockValue(Kp,T)), ylabel('Kp')
subplot(322), viewSurf(setBlockValue(Ki,T)), ylabel('Ki')
subplot(323), viewSurf(setBlockValue(Kt,T)), ylabel('Kt')
subplot(324), viewSurf(setBlockValue(a,T)), ylabel('a')
subplot(325), viewSurf(setBlockValue(b,T)), ylabel('b')

%%
% You can further validate the design by simulating the linear responses 
% at each design point. However, you need to return to Simulink to simulate 
% the nonlinear response of the gain-scheduled controller.