www.gusucode.com > mpc_featured 案例源码 matlab代码程序 > mpc_featured/mpccstr.m
%% Gain-Scheduled MPC Control of Nonlinear Chemical Reactor % This example shows how to use multiple MPC controllers to control a % nonlinear continuous stirred tank reactor (CSTR) as it transitions from % low conversion rate to high conversion rate. % % Multiple MPC Controllers are designed at different operating conditions % and then implemented with the Multiple MPC Controller block in Simulink. % At run time, a scheduling signal is used to switch controller from one to % another. % Copyright 1990-2014 The MathWorks, Inc. %% About the Continuous Stirred Tank Reactor % A Continuously Stirred Tank Reactor (CSTR) is a common chemical system in % the process industry. A schematic of the CSTR system is: % % <<../mpc_cstr.png>> %% % This is a jacketed non-adiabatic tank reactor described extensively in % Seborg's book, "Process Dynamics and Control", published by Wiley, 2004. % The vessel is assumed to be perfectly mixed, and a single first-order % exothermic and irreversible reaction, A --> B, takes place. The inlet % stream of reagent A is fed to the tank at a constant volumetric rate. % The product stream exits continuously at the same volumetric rate and % liquid density is constant. Thus the volume of reacting liquid is % constant. % % The inputs of the CSTR model are: % % $$ \begin{array} {ll} % u_1 = CA_i \; & \textnormal{Concentration of A in inlet feed % stream} [kgmol/m^3] \\ % u_2 = T_i \; & \textnormal{Inlet feed stream temperature} [K] \\ % u_3 = T_c \; & \textnormal{Jacket coolant temperature} [K] \\ % \end{array} $$ % % and the outputs (y(t)), which are also the states of the model (x(t)), are: % % $$ \begin{array} {ll} % y_1 = x_1 = CA \; & \textnormal{Concentration of A in reactor tank} [kgmol/m^3] \\ % y_2 = x_2 = T \; & \textnormal{Reactor temperature} [K] \\ % \end{array} $$ % % The control objective is to maintain the concentration of reagent A, $CA$ % at its desired setpoint, which changes over time when reactor transitions % from low conversion rate to high conversion rate. The coolant % temperature $T_c$ is the manipulated variable used by the MPC controller % to track the reference. The inlet feed stream concentration and % temperature are assumed to be constant. The Simulink model % |mpc_cstr_plant| implements the nonlinear CSTR plant. %% About Gain Scheduled Model Predictive Control % It is well known that the CSTR dynamics are strongly nonlinear with % respect to reactor temperature variations and can be open-loop unstable % during the transition from one operating condition to another. A single % MPC controller designed at a particular operating condition cannot give % satisfactory control performance over a wide operating range. %% % To control the nonlinear CSTR plant with linear MPC control technique, % you have a few options: % % * If a linear plant model cannot be obtained at run time, first you need % to obtain several linear plant models offline at different operating % conditions that cover the typical operating range. Next you can choose % one of the two approaches to implement MPC control strategy: % % (1) Design several MPC controllers offline, one for each plant model. At % run time, use Multiple MPC Controller block that switches MPC controllers % from one to another based on a desired scheduling strategy, as discussed % in this example. Use this approach when the plant models have different % orders or time delays. % % (2) Design one MPC controller offline at a nominal operating point. At % run time, use Adaptive MPC Controller block (updating predictive model at % each control interval) together with Linear Parameter Varying (LPV) % System block (supplying linear plant model with a scheduling strategy). % See <docid:mpc_examples.example-ex18818571> for more details. Use % this approach when all the plant models have the same order and time % delay. % % * If a linear plant model can be obtained at run time, you should use % Adaptive MPC Controller block to achieve nonlinear control. There are % two typical ways to obtain a linear plant model online: % % (1) Use successive linearization. See % <docid:mpc_examples.example-ex07272954> for more details. Use this % approach when a nonlinear plant model is available and can be linearized % at run time. % % (2) Use online estimation to identify a linear model when loop is closed. % See <docid:mpc_examples.example-ex24638454> for more details. Use this % approach when linear plant model cannot be obtained from either an LPV % system or successive linearization. %% Obtain Linear Plant Model at Initial Operating Condition % To run this example, Simulink(R) and Simulink Control Design(R) are % required. if ~mpcchecktoolboxinstalled('simulink') disp('Simulink(R) is required to run this example.') return end if ~mpcchecktoolboxinstalled('slcontrol') disp('Simulink Control Design(R) is required to run this example.') return end %% % First, a linear plant model is obtained at the initial operating % condition, CAi is 10 kgmol/m^3, Ti and Tc are 298.15 K. Functions from % Simulink Control Design such as |operspec|, |findop|, and |linearize|, are % used to generate the linear state-space system from the Simulink model. %% % Create operating point specification. plant_mdl = 'mpc_cstr_plant'; op = operspec(plant_mdl); %% % Feed concentration is known at the initial condition. op.Inputs(1).u = 10; op.Inputs(1).Known = true; %% % Feed temperature is known at the initial condition. op.Inputs(2).u = 298.15; op.Inputs(2).Known = true; %% % Coolant temperature is known at the initial condition. op.Inputs(3).u = 298.15; op.Inputs(3).Known = true; %% % Compute initial condition. [op_point,op_report] = findop(plant_mdl,op); % Obtain nominal values of x, y and u. x0 = [op_report.States(1).x; op_report.States(2).x]; y0 = [op_report.Outputs(1).y; op_report.Outputs(2).y]; u0 = [op_report.Inputs(1).u; op_report.Inputs(2).u; op_report.Inputs(3).u]; %% % Obtain linear model at the initial condition. plant = linearize(plant_mdl,op_point); %% % Verify that the linear model is open-loop stable at this condition. eig(plant) %% Design MPC Controller for Initial Operating Condition % You design an MPC at the initial operating condition. Ts = 0.5; %% % Specify signal types used in MPC. Assume both reactor temperature and % concentration are measurable. plant.InputGroup.UnmeasuredDisturbances = [1 2]; plant.InputGroup.ManipulatedVariables = 3; plant.OutputGroup.Measured = [1 2]; plant.InputName = {'CAi','Ti','Tc'}; plant.OutputName = {'T','CA'}; %% % Create MPC controller with default prediction and control horizons mpcobj = mpc(plant,Ts); %% % Set nominal values in the controller. Note that nominal values for % unmeasured disturbance must be zero. mpcobj.Model.Nominal = struct('X',x0,'U',[0;0;u0(3)],'Y',y0,'DX',[0 0]); %% % Set scale factors because plant input and output signals have different % orders of magnitude. Uscale = [10;30;50]; Yscale = [50;10]; mpcobj.DV(1).ScaleFactor = Uscale(1); mpcobj.DV(2).ScaleFactor = Uscale(2); mpcobj.MV.ScaleFactor = Uscale(3); mpcobj.OV(1).ScaleFactor = Yscale(1); mpcobj.OV(2).ScaleFactor = Yscale(2); %% % The goal will be to track a specified transition in the reactor % concentration. The reactor temperature will be measured and used in % state estimation but the controller will not attempt to regulate it % directly. It will vary as needed to regulate the concentration. Thus, % set its MPC weight to zero. mpcobj.Weights.OV = [0 1]; %% % Plant inputs 1 and 2 are unmeasured disturbances. By default, the % controller assumes integrated white noise with unit magnitude at these % inputs when configuring the state estimator. Try increasing the state % estimator signal-to-noise by a factor of 10 to improve disturbance % rejection performance. Dist = ss(getindist(mpcobj)); Dist.B = eye(2)*10; setindist(mpcobj,'model',Dist); %% % All other MPC parameters are at their default values. %% Test the Controller With a Step Disturbance in Feed Concentration % "mpc_cstr_single" contains a Simulink(R) model with CSTR and MPC % Controller blocks in a feedback configuration. mpc_mdl = 'mpc_cstr_single'; open_system(mpc_mdl) %% % Note that the MPC Controller block is configured to look ahead at (preview) % setpoint changes in the future; that is, anticipating the setpoint % transition. This generally improves setpoint tracking. %% % Define a constant setpoint for the output. CSTR_Setpoints.time = [0; 60]; CSTR_Setpoints.signals.values = [y0 y0]'; %% % Test the response to a 5% increase in feed concentration. set_param([mpc_mdl '/Feed Concentration'],'Value','10.5'); %% % Set plot scales and simulate the response. open_system([mpc_mdl '/Measurements']) open_system([mpc_mdl '/Coolant Temperature']) set_param([mpc_mdl '/Measurements'],'Ymin','305~8','Ymax','320~9') set_param([mpc_mdl '/Coolant Temperature'],'Ymin','295','Ymax','305') sim(mpc_mdl,10); %% % The closed-loop response is satisfactory. %% Simulate Designed MPC Controller Using Full Transition % First, define the desired setpoint transition. After a 10-minute warm-up % period, ramp the concentration setpoint downward at a rate of 0.25 per % minute until it reaches 2.0 kmol/m^3. CSTR_Setpoints.time = [0 10 11:39]'; CSTR_Setpoints.signals.values = [y0(1)*ones(31,1),[y0(2);y0(2);(y0(2):-0.25:2)';2;2]]; %% % Remove the 5% increase in feed concentration used previously. set_param([mpc_mdl '/Feed Concentration'],'Value','10') %% % Set plot scales and simulate the response. set_param([mpc_mdl '/Measurements'],'Ymin','300~0','Ymax','400~10') set_param([mpc_mdl '/Coolant Temperature'],'Ymin','240','Ymax','360') %% % Simulate model. sim(mpc_mdl,60) %% % The closed-loop response is unacceptable. Performance along the full % transition can be improved if other MPC controllers are designed at % different operating conditions along the transition path. In the next % two section, two additional MPC controllers are design at intermediate % and final transition stages respectively. %% Design MPC Controller for Intermediate Operating Condition % Create operating point specification. op = operspec(plant_mdl); %% % Feed concentration is known. op.Inputs(1).u = 10; op.Inputs(1).Known = true; %% % Feed temperature is known. op.Inputs(2).u = 298.15; op.Inputs(2).Known = true; %% % Reactor concentration is known. op.Outputs(2).y = 5.5; op.Outputs(2).Known = true; %% % Find steady state operating condition. [op_point,op_report] = findop(plant_mdl,op); % Obtain nominal values of x, y and u. x0 = [op_report.States(1).x; op_report.States(2).x]; y0 = [op_report.Outputs(1).y; op_report.Outputs(2).y]; u0 = [op_report.Inputs(1).u; op_report.Inputs(2).u; op_report.Inputs(3).u]; %% % Obtain linear model at the initial condition. plant_intermediate = linearize(plant_mdl,op_point); %% % Verify that the linear model is open-loop unstable at this condition. eig(plant_intermediate) %% % Specify signal types used in MPC. Assume both reactor temperature and % concentration are measurable. plant_intermediate.InputGroup.UnmeasuredDisturbances = [1 2]; plant_intermediate.InputGroup.ManipulatedVariables = 3; plant_intermediate.OutputGroup.Measured = [1 2]; plant_intermediate.InputName = {'CAi','Ti','Tc'}; plant_intermediate.OutputName = {'T','CA'}; %% % Create MPC controller with default prediction and control horizons. mpcobj_intermediate = mpc(plant_intermediate,Ts); %% % Set nominal values, scale factors and weights in the controller. mpcobj_intermediate.Model.Nominal = struct('X',x0,'U',[0;0;u0(3)],'Y',y0,'DX',[0 0]); Uscale = [10;30;50]; Yscale = [50;10]; mpcobj_intermediate.DV(1).ScaleFactor = Uscale(1); mpcobj_intermediate.DV(2).ScaleFactor = Uscale(2); mpcobj_intermediate.MV.ScaleFactor = Uscale(3); mpcobj_intermediate.OV(1).ScaleFactor = Yscale(1); mpcobj_intermediate.OV(2).ScaleFactor = Yscale(2); mpcobj_intermediate.Weights.OV = [0 1]; Dist = ss(getindist(mpcobj_intermediate)); Dist.B = eye(2)*10; setindist(mpcobj_intermediate,'model',Dist); %% Design MPC Controller for Final Operating Condition % Create operating point specification. op = operspec(plant_mdl); %% % Feed concentration is known. op.Inputs(1).u = 10; op.Inputs(1).Known = true; %% % Feed temperature is known. op.Inputs(2).u = 298.15; op.Inputs(2).Known = true; %% % Reactor concentration is known. op.Outputs(2).y = 2; op.Outputs(2).Known = true; %% % Find steady-state operating condition. [op_point,op_report] = findop(plant_mdl,op); % Obtain nominal values of x, y and u. x0 = [op_report.States(1).x; op_report.States(2).x]; y0 = [op_report.Outputs(1).y; op_report.Outputs(2).y]; u0 = [op_report.Inputs(1).u; op_report.Inputs(2).u; op_report.Inputs(3).u]; %% % Obtain linear model at the initial condition. plant_final = linearize(plant_mdl,op_point); %% % Verify that the linear model is again open-loop stable at this condition. eig(plant_final) %% % Specify signal types used in MPC. Assume both reactor temperature and % concentration are measurable. plant_final.InputGroup.UnmeasuredDisturbances = [1 2]; plant_final.InputGroup.ManipulatedVariables = 3; plant_final.OutputGroup.Measured = [1 2]; plant_final.InputName = {'CAi','Ti','Tc'}; plant_final.OutputName = {'T','CA'}; %% % Create MPC controller with default prediction and control horizons. mpcobj_final = mpc(plant_final,Ts); %% % Set nominal values, scale factors and weights in the controller. mpcobj_final.Model.Nominal = struct('X',x0,'U',[0;0;u0(3)],'Y',y0,'DX',[0 0]); Uscale = [10;30;50]; Yscale = [50;10]; mpcobj_final.DV(1).ScaleFactor = Uscale(1); mpcobj_final.DV(2).ScaleFactor = Uscale(2); mpcobj_final.MV.ScaleFactor = Uscale(3); mpcobj_final.OV(1).ScaleFactor = Yscale(1); mpcobj_final.OV(2).ScaleFactor = Yscale(2); mpcobj_final.Weights.OV = [0 1]; Dist = ss(getindist(mpcobj_final)); Dist.B = eye(2)*10; setindist(mpcobj_final,'model',Dist); %% Control the CSTR Plant With the Multiple MPC Controllers Block % The following model uses the Multiple MPC Controllers block to implement % three MPC controllers across the operating range. mmpc_mdl = 'mpc_cstr_multiple'; open_system(mmpc_mdl); %% % Note that it has been configured to use the three controllers in a % sequence: mpcobj, mpcobj_intermediate and mpcobj_final. open_system([mmpc_mdl '/Multiple MPC Controllers']); %% % Note also that the two switches specify when to switch from one % controller to another. The rules are: % 1. If CSTR concentration >= 8, use "mpcobj" % 2. If 3 <= CSTR concentration < 8, use "mpcobj_intermediate" % 3. If CSTR concentration < 3, use "mpcobj_final" %% % Simulate with the Multiple MPC Controllers block open_system([mmpc_mdl '/Measurements']); open_system([mmpc_mdl '/MV']); sim(mmpc_mdl) %% % The transition is now well controlled. The major improvement is in the % transition through the open-loop unstable region. The plot of the % switching signal shows when controller transitions occur. The MV % character changes at these times because of the change in dynamic % characteristics introduced by the new prediction model. % %% bdclose(plant_mdl) bdclose(mpc_mdl) bdclose(mmpc_mdl)