www.gusucode.com > ident_featured 案例代码 matlab源码程序 > ident_featured/iddemo1.m
%% Estimating Simple Models from Real Laboratory Process Data % This example shows how to develop and analyze simple models from a real % laboratory process data. We start with a small description of the % process, learn how to import the data to the toolbox and % preprocess/condition it and then proceed systematically to estimate % parametric and nonparametric models. Once the models have been identified % we compare the estimated models and also validate the model to the actual % output data from the experiment. % Copyright 1986-2015 The MathWorks, Inc. %% System Description % This case study concerns data collected from a laboratory scale % "hairdryer". (Feedback's Process Trainer PT326; See also page 525 in % Ljung, 1999). The process works as follows: Air is fanned through a tube % and heated at the inlet. The air temperature is measured by a % thermocouple at the outlet. The input is the voltage over the heating % device, which is just a mesh of resistor wires. The output is the outlet % air temperature represented by the measured thermocouple voltage. %% Setting up Data for Analysis % First we load the input-output data to the MATLAB(R) Workspace. load dryer2; %% % Vector |y2|, the output, contains 1000 measurements of the themocouple % voltage which is proportional to the temperature in the outlet airstream. % Vector |u2| contains 1000 input data points consisting of the voltage % applied to the heater. The input was generated as a binary random % sequence that switches from one level to the other with probability 0.2. % The sample time is 0.08 seconds. % % The next step is to set up the data as an iddata object dry = iddata(y2,u2,0.08); %% % To get information about the data, just type the name of the |iddata| % object at the MATLAB command window: dry %% % To inspect the properties of the above iddata object, use the |get| % command: get(dry) %% % For better book-keeping, it is good practice to give names to the input % and output channels and Time units. These names would be propagated % throughout the analysis of this iddata object: dry.InputName = 'Heater Voltage'; dry.OutputName = 'Thermocouple Voltage'; dry.TimeUnit = 'seconds'; dry.InputUnit = 'V'; dry.OutputUnit = 'V'; %% % Now that we have the data set ready, we choose the first 300 data points % for model estimation. ze = dry(1:300) %% Preprocessing the Data % Plot the interval from sample 200 to 300: % plot(ze(200:300)); %% % *Figure 1:* A snapshot of the measured hair-dryer data. %% % From the above plot, it can be observed that the data is not zero mean. % So let us remove the constant levels and make the data zero mean. ze = detrend(ze); %% % The same data set after it has been detrended: plot(ze(200:300)) %show samples from 200 to 300 of detrended data %% % *Figure 2:* Detrended estimation data. %% Estimating Nonparametric and Parametric Models % Now that the dataset has been detrended and there are no obvious % outliers, let us first estimate the impulse response of the system % by correlation analysis to get some idea of time constants and the % like: clf mi = impulseest(ze); % non-parametric (FIR) model showConfidence(impulseplot(mi),3); %impulse response with 3 standard %deviations confidence region %% % *Figure 3:* Impulse response of the FIR model estimated using |ze|. %% % The shaded region marks a 99.7% confidence interval. There is a time % delay (dead-time) of 3 samples before the output responds to the input % (significant output outside the confidence interval). %% % The simplest way to get started on a parametric estimation routine is to % build a state-space model where the model-order is automatically % determined, using a prediction error method. Let us estimate a model % using the |ssest| estimation technique: m1 = ssest(ze); %% % |m1| is a continuous-time identified state-space model, represented by an % |idss| object. The estimation algorithm chooses 3 as the optimal order of % the model. To inspect the properties of the estimated model, just enter % the model name at the command window: m1 %% % The display suggests that the model is free-form (all entries of A, B and % C matrices were treated as free parameters) and that the estimated model % fits the data pretty well (over 90% fit). To retrieve the properties of % this model, for example to obtain the |A| matrix of the discrete % state-space object generated above, we can use the dot operator: % % A = m1.a; % % See the "Data and Model Objects in System Identification Toolbox" example % for more information regarding model objects. To find out which % properties of the model object can be retrieved, use |get| command: get(m1) %% % To fetch the values of the state-space matrices and their 1 std % uncertainties, use the |idssdata| command: [A,B,C,D,K,~,dA,dB,dC,dD,dK] = idssdata(m1) %% % The uncertainties are quite large even though the model fit the % estimation data well. This is because the model is over-parameterized, % that is, it has more free parameters than what could be uniquely % identified. The variance of parameters in such cases is not well defined. % However this does not imply that the model is unreliable. We can plot the % time- and frequency-response of this plot and view the variance as % confidence regions as discussed next. %% Analyzing the Estimated Model % % The Bode plot of the generated model can be obtained using the % |bode| function as shown below: h = bodeplot(m1); %% % *Figure 4:* Bode plot of estimated model. %% % Right-click on the plot and pick Characteristics->Confidence Region. % Or, use the |showConfidence| command to view the variance of the % response. showConfidence(h,3) % 3 std (99.7%) confidence region %% % *Figure 5:* Bode plot with 3 std confidence region. %% % Similarly, we can generate the step plot and its associated 3 std % confidence region. We can compare the responses and associated variances % of the parameteric model |m1| with that of the nonparametric model |mi|: showConfidence(stepplot(m1,'b',mi,'r',3),3) %% % *Figure 6:* Step plot of models |m1| and |mi| with confidence regions. %% % We can also consider the Nyquist plot, and mark uncertainty % regions at certain frequencies with ellipses, corresponding to 3 % standard deviations: Opt = nyquistoptions; Opt.ShowFullContour = 'off'; showConfidence(nyquistplot(m1,Opt),3) %% % *Figure 7:* Nyquist plot of estimated model showing the uncertainty regions at certain frequencies. %% % The response plots show that the estimated model |m1| is quite reliable. %% Estimating Models with a Prescribed Structure % System Identification Toolbox can also be used to obtain a model with % a prescribed structure. For example, a difference equation model with 2 % poles, 1 zero and 3 sample delays can be obtained using the |arx| % function as shown below: m2 = arx(ze,[2 2 3]); %% % To look at the model, enter the model name at the command window. m2 %% % A continuous time transfer function with 2 poles, one zero and 0.2 second % transport delay can be estimated using the |tfest| command: m3 = tfest(ze, 2, 1, 0.2) %% Validating the Estimated Model to Experimental Output % How good is an estimated model? One way to find out is to simulate it and % compare the model output with measured output. Select a portion of the % original data that was not used in building the model, say from samples % 800 to 900. Once the validation data has been preprocessed, we use the % |compare| function as shown below to view the quality of prediction: zv = dry(800:900); % select an independent data set for validation zv = detrend(zv); % preprocess the validation data set(gcf,'DefaultLegendLocation','best') compare(zv,m1); % perform comparison of simulated output %% % *Figure 8:* Model's simulated response vs. validation data output. %% % It can be observed here that the agreement is very good. The "Fit" value % shown is calculated as: % % |Fit = 100*(1 - norm(yh - y)/norm(y-mean(y)))| % % where |y| is the measured output (=|zv.y|), and |yh| is the output of the % model |m1|. %% Comparing Estimated Models % To compare the performance of the models that we have estimated, for % example |m1|, |m2| and |m3| with the validation data |zv|, we can again % use the |compare| command: compare(zv,m1,'b',m2,'r',m3,'c'); %% % *Figure 9:* Comparing the responses of models |m1|, |m2|, |m3| on % validation data set |ze|. %% % The pole-zero plots for the models can be obtained using |iopzplot|: h = iopzplot(m1,'b',m2,'r',m3,'c'); %% % *Figure 10:* Poles and zeros of the models |m1|, |m2| and |m3|. %% % The uncertainties in the poles and zeroes can also be obtained. In the % following statement, '3' refers to the number of standard deviations. showConfidence(h,3); %% % *Figure 11:* Pole-zero map with uncertainty regions. %% % The frequency functions above that are obtained from the models can be % compared with one that is obtained using a non-parametric spectral % analysis method (|spa|): gs = spa(ze); %% % The |spa| command produces an IDFRD model. The |bode| function can again be % used for a comparison with the transfer functions of the models obtained. % w = linspace(0.4,pi/m2.Ts,200); opt = bodeoptions; opt.PhaseMatching = 'on'; bodeplot(m1,'b',m2,'r',m3,'c',gs,'g',w,opt); legend('m1','m2','m3','gs') %% % *Figure 12:* Bode responses of |m1|, |m2| and |m3| compared against the % non-parametric spectral estimation model |gs|. %% % The frequency responses from the three models/methods are very close. % This indicates that this response is reliable. % % Also, a Nyquist plot can be analyzed with the uncertainty regions marked % at certain frequencies: showConfidence(nyquistplot(m1,'b',m2,'r',m3,'c',gs,'g'),3) %% % *Figure 13:* Nyquist plots of models |m1|, |m2|, |m3| and |gs|. %% % The non-parametric model |gs| exhibits the most uncertainty in response. %% Additional Information % For more information on identification of dynamic systems with System % Identification Toolbox visit the % <http://www.mathworks.com/products/sysid/ System Identification Toolbox> product % information page.