www.gusucode.com > ident_featured 案例代码 matlab源码程序 > ident_featured/OnlinePolynomialModelEstimationExample.m
%% Online ARMAX Polynomial Model Estimation % % This example shows how to implement an online polynomial model estimator. % You estimate two ARMAX models for a nonlinear chemical reaction process. % These models capture the behavior of the process at two operating % conditions. The model behavior is identified online and used to adjust % the gains of an adaptive PI controller during system operation. % Copyright 2013-2014 The MathWorks, Inc. %% Continuously 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: % % <<../cstr.png>> %% % This is a jacketed diabatic (i.e., nondiabatic) tank reactor described % extensively in Bequette's book "Process Dynamics: Modeling, Analysis and % Simulation", published by Prentice-Hall, 1998. 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 rate. After stirring, the end product streams out % of the vessel at the same rate as reagent A is fed into the tank (the % volume in the reactor tank is constant). Details of the operation of the % CSTR and its 2-state nonlinear model used in this example are explained in % the example <docid:ident_examples.example-ex36617320>. % % The inputs of the CSTR model are: % % $$ \begin{array} {ll} % u_1(t) = C_{Af}(t) \; & \textnormal{Concentration of A in inlet feed % stream} [kgmol/m^3] \\ % u_2(t) = T_f(t) \; & \textnormal{Inlet feed stream temperature} [K] \\ % u_3(t) = T_j(t) \; & \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(t) = x_1(t) = C_{A}(t) \; & \textnormal{Concentration of A in reactor tank} [kgmol/m^3] \\ % y_2(t) = x_2(t) = T(t) \; & \textnormal{Reactor temperature} [K] \\ % \end{array} $$ % % The control objective is to maintain the concentration of reagent A, % $C_{A}(t)$ at the desired level $C_{Aref}(t)$, which changes over time. % The jacket temperature $T_j(t)$ is manipulated by a PI controller in % order to reject disturbances arising from the inlet feed stream % temperature $T_f(t)$. The input of the PI controller is the tracking % error signal, $C_{Aref}(t)-C_A(t)$. The inlet feed stream concentration, % $C_{Af}(t)$, is assumed to be constant. The Simulink model |iddemo_cstr| % implements the CSTR plant as the block |CSTR|. open_system('iddemo_cstr'); %% % The $T_f(t)$ feed temperature input consists of a white noise disturbance % on top of a constant offset. The noise power is 0.0075 [K^2]. This noise % level causes up to 2% deviation from the desired $C_{Aref}(t)$. % % The $C_{Aref}(t)$ signal in this example contains a step change from 1.5 % [kgmol/m^3] to 2 [kgmol/m^3] at time $t=400$. In addition to this % step change, $C_{Aref}(t)$ also contains a white noise perturbation for t % in the [0,200) and [400,600) ranges. The power of this white noise % signal is 0.015. The noise power is adjusted empirically to approximately % give a signal-to-noise ratio of 10. Not having sufficient excitation in % the reference signal in closed-loop identification can lead to not having % sufficient information to identify a unique model. The implementation of % $C_{Aref}(t)$ is in the |iddemo_cstr/CA Reference| block. %% Online Estimation for Adaptive Control % % It is known from the nonlinear model that the CSTR output $C_{A}(t)$ is % more sensitive to the control input $T_j(t)$ at higher $C_{A}(t)$ levels. % The Recursive Polynomial Model Estimator block is used to detect this % change in sensitivity. This information is used to adjust the gains of % the PI controller as $C_{A}(t)$ varies. The aim is to avoid having a a % high gain control loop which may lead to instability. % % You estimate a discrete transfer-function from $T_j(t)$ to $C_{A}(t)$ % online with the Recursive Polynomial Model Estimator block. The adaptive % control algorithm uses the DC gain of this transfer function. The % tracking error $C_{Aref}(t)-C_A(t)$, is divided by the normalized DC % gain of the estimated transfer function. This normalization is done to % have a gain of 1 on the tracking error at the initial operating point, % for which the PI controller is designed. For instance, the error signal % is divided by 2 if the DC gain becomes 2 times its original value. % This corresponds to dividing the PI controller gains by 2. This adpative % controller is implemented in |iddemo_cstr/Adaptive PI Controller|. %% Recursive Polynomial Model Estimator Block Inputs % % The 'Recursive Polynomial Model Estimator' block is found under the % |System Identification Toolbox/Estimators| library in Simulink. You % use this block to estimate linear models with ARMAX structure. ARMAX % models have the form: % % $$ % \begin{array} {l}A(q)\bar{y}(t)=B(q)\bar{u}(t)+C(q)\bar{e}(t) \\[0.1in] % A(q) = 1+a_1z^{-1}+a_2z^{-2}+\cdots+a_{na}z^{-na} \\ % B(q) = (b_01+b_1z^{-1}+b_2z^{-2}+\cdots+a_{nb-1}z^{-nb+1})z^{-nk} \\ % C(q) = 1+c_1z^{-1}+c_2z^{-2}+\cdots+c_{nc}z^{-nc} \\ % \end{array} % $$ % % * The |Inputs| and |Output| inport of the recursive polynomial model % estimator block correspond to $\bar{u}(t)$ and $\bar{y}(t)$ respectively. % For the CSTR model $\bar{y}$ and $\bar{u}$ are deviations from the jacket % temperature and A concentration trim operating points: % $\bar{y}=C_A(t)-\bar{C}_A(t)$, $\bar{u}=T_j(t)-\bar{T}_j(t)$. It is good % to scale $\bar{u}$ and $\bar{y}$ to have a peak amplitude of 1 to improve % the numerical condition of the estimation problem. The trim operating % points, $\bar{C}_A(t)$ and $\bar{T}_j(t)$, are not known exactly before % system operation. They are estimated and extracted from the measured % signals by using a first-order moving average filter. These preprocessing % filters are implemented in the |iddemo_cstr/Preprocess Tj| and % |iddemo_cstr/Preprocess CA| blocks. open_system('iddemo_cstr/Preprocess Tj'); %% % * The optional |Enable| inport of the Recursive Polynomial Model % Estimator block controls the parameter estimation in the block. Parameter % estimation is disabled when the Enable signal is zero. Parameter % estimation is enabled for all other values of the |Enable| signal. In % this example the estimation is disabled for the time intervals % $t\in[200,400)$ and $t\in[600,800)$. During these intervals the measured % input $T_j(t)$ does not contain sufficient excitation for closed-loop % system identification. %% Recursive Polynomial Model Estimator Block Setup % % Configure the block parameters to estimate a second-order ARMAX model. In % the *Model Parameters* tab, specify: % % * *Model Structure:* ARMAX. Choose ARMAX since the current and past % values of the disturbances acting on the system, $T_f(t)$, are expected % to impact the CSTR system output $C_A(t)$. % % * *Initial Estimate:* None. By default, the software uses a value of 0 % for all estimated parameters. % % * *Number of parameters in A(q) (na)*: |2|. The nonlinear model has % 2 states. % % * *Number of parameters in B(q) (nb)*: |2|. % % * *Number of parameters in C(q) (nc)*: |2|. The estimated model % corresponds to a second order model since the maximum of na, nb, and nc % are 2. % % * *Input Delay (nk)*: |1|. Like most physical systems, the CSTR system % does not have direct feedthrough. Also, there are no extra time delays % between its I/Os. % % * *Parameter Covariance Matrix*: |1e4|. Specify a high covariance value % because the initial guess values are highly uncertain. % % * *Sample Time*: |0.1|. The CSTR model is known to have a bandwidth of % about 0.25Hz. $T_s=0.1$ chosen such that 1/0.1 is greater than 20 times % this bandwidth (5Hz). %% % <<../rarmax_dialogParams1.png>> %% % Click *Algorithm and Block Options* to set the estimation options: % % * *Estimation Method*: |Forgetting Factor| % % * *Forgetting Factor*: |1-5e-3|. Since the estimated parameters are % expected to change with the operating point, set the forgetting factor to % a value less than 1. Choose $\lambda = 1-5e-3$ which corresponds to a % memory time constant of $T_0=\frac{Ts}{1-\lambda}=100$ seconds. A 100 % second memory time ensures that a significant amount data used for % identification is coming from the 200 second identification period at % each operating point. % % * Select the *Output estimation error* check box. You use this block % output to validate the estimation. % % * Select the *Add enable port* check box. You only want to adapt the % estimated model parameters when extra noise is injected in the reference % port. The parameter estimation n is disabled through this port when the % extra noise is no longer injected. % % * *External reset*: |None|. %% % <<../rarmax_dialogParams2.png>> %% Recursive Polynomial Model Estimator Block Outputs % % At every time step, the recursive polynomial model estimator provides an % estimate for $A(q)$, $B(q)$, $C(q)$, and the estimation error $\bar{e}$. % The |Error| outport of the polynomial model estimator block contains % $\bar{e}(t)$ and is also known as the one-step-ahead prediction error. % The |Parameters| outport of the block contains the A(q), B(q), and C(q) % polynomial coefficients in a bus signal. Given the chosen polynomial % orders ($na=2$, $nb=2$, $nc=2$, $nk=1$) the |Parameters| bus elements % contain: % % $$ \begin{array} {l} % A(t) = [1 \; a_1(t) \; a_2(t)] \\ % B(t) = [0 \; b_0(t) \; b_1(t)] \\ % C(t) = [1 \; c_1(t) \; c_2(t)] \\ % \end{array} $$ % % The estimated parameters in the A(q), B(q), and C(q) polynomials change % during simulation as follows: sim('iddemo_cstr'); close_system('iddemo_cstr/Preprocess Tj'); open_system('iddemo_cstr/ABC'); %% % The parameter estimates quickly change from their initial values of 0 due % to the high value chosen for the initial parameter covariance matrix. The % parameters in the $A(q)$ and $B(q)$ polynomials approach their values at % $t=200$ rapidly. However, the parameters in the $C(q)$ polynomial show % some fluctuations. One reason behind these fluctuations is that the % disturbance $T_f(t)$ to CSTR output $C_A(t)$ is not fully modelled by the % ARMAX structure. The error model $C(q)$ is not important for the control % problem studied here since the $T_j(t)$ to $C_A(t)$ relationship is % captured by the transfer function $\frac{B(q)}{A(q)}$. Therefore, the % fluctuation in $C(q)$ is not a concern for this identification and % control problem. % % The parameter estimates are held constant for $t\in[200,400)$ since the % estimator block was disabled for this interval (0 signal to the |Enable| % inport). The parameter estimation is enabled at $t=400$ when the CSTR % tank starts switching to its new operating point. The parameters of % $A(q)$ and $B(q)$ converge to their new values by $t=600$, and then held % constant by setting the |Enable| port to 0. The convergence of $A(q)$ and % $B(q)$ is slower at this operating point. This slow convergence is % because of the smaller eigenvalues of the parameter covariance matrix at % t=400 compared to the initial 1e4 values set at t=0. The parameter % covariance, which is a measure of confidence in the estimates, is updated % with each time step. The algorithm quickly changed the parameter % estimates when the confidence in estimates were low at t=0. The improved % parameter estimates capture the system behavior better, resulting in % smaller one-step-ahead prediction errors and smaller eigenvalues in the % parameter covariance matrix (increased confidence). The system behavior % changes in t=400. However, the block is slower to change the parameter % estimates due to the increased confidence in the estimates. You can use % the External Reset option of the Recursive Polynomial Model Estimator % block to provide a new value for parameter covariance at t=400. To see % the value of the parameter covariance, select the *Output parameter % covariance matrix* check box in the Recursive Polynomial Model Estimator % block. %% Validating the Estimated Model % % The |Error| output of the |Recursive Polynomial Model Estimator| block % gives the one-step-ahead error for the estimated model. close_system('iddemo_cstr/ABC'); open_system('iddemo_cstr/1-Step Ahead Prediction Error'); %% % % The one-step-ahead error is higher when there are no extra perturbations % injected in the $T_j(t)$ channel for system identification. These higher % errors may be caused by the lack of sufficient information in the % $T_j(t)$ input channel that the estimator block relies on. However, even % this higher error is low and bounded when compared to the measured % fluctuations in $C_A(t)$. This gives confidence in the estimated % parameter values. % % A more rigorous check of the estimated model is to simulate the estimated % model and compare with the actual model output. The % |iddemo_cstr/Time-Varying ARMAX| block implements the time-varying ARMAX % model estimated by the Online Polynomial Model Estimator block. The error % between the output of the CSTR system and the estimated time-varying % ARMAX model output is: % close_system('iddemo_cstr/1-Step Ahead Prediction Error'); open_system('iddemo_cstr/Simulation Error'); %% % % The simulation error is again bounded and low when compared to the % fluctuations in the $C_A(t)$. This further provides confidence that the % estimated linear models are able to predict the nonlinear CSTR model % behavior. % %% % The identified models can be further analyzed in MATLAB. The model % estimates for the operating points $C_A=1.5 [kgmol/m^3]$ and $C_A=2 % [kgmol/m^3]$ can be obtained by looking at the estimated A(q), B(q), and % C(q) polynomials at $t = 200$ and $t = 600$ respectively. Bode plots of % these models are: Ts = 0.1; tidx = find(t>=200,1); P200 = idpoly(AHat(:,:,tidx),BHat(:,:,tidx),CHat(:,:,tidx),1,1,[],Ts); tidx = find(t>=600,1); P600 = idpoly(AHat(:,:,tidx),BHat(:,:,tidx),CHat(:,:,tidx),1,1,[],Ts); bodemag(P200,'b',P600,'r--',{10^-1,20}); legend('Estimated Model at C_A=1.5 [kgmol/m^3]', ... 'Estimated Model at C_A=2.0 [kgmol/m^3]', ... 'Location', 'Best'); %% % The estimated model has a higher gain at higher concentration levels. % This is in agreement with prior knowledge about the nonlinear CSTR % plant. The transfer function at $C_A(t)=2 [kgmol/m^3]$ has a $6dB$ higher % gain (double the amplitude) at low frequencies. %% Summary % % You estimated two ARMAX models to capture the behavior of the nonlinear % CSTR plant at two operating conditions. The estimation was done during % closed-loop operation with an adaptive controller. You looked at two % signals to validate the estimation results: One step ahead prediction % errors and the errors between the CSTR plant output and the simulation of % the estimation model. Both of these errors signals were bounded and small % compared to the CSTR plant output. This provided confidence in the % estimated ARMAX model parameters. %% % bdclose('iddemo_cstr');