www.gusucode.com > ident_featured 案例代码 matlab源码程序 > ident_featured/idnlgreydemo2.m

    %% Two Tank System: C MEX-File Modeling of Time-Continuous SISO System
% This example shows how to perform IDNLGREY modeling based on C MEX model
% files. It uses a simple system where nonlinear state space modeling
% really pays off.

%   Copyright 2005-2014 The MathWorks, Inc.

%% A Two Tank System
% The objective is to model the liquid level of the lower tank of a
% laboratory scale two tank system, schematically shown in Figure 1.
%
% <<../twotank.png>>
%
% *Figure 1:* Schematic view of a two tank system.

%% Input-Output Data
% We start the modeling job by loading the available input-output data,
% which was simulated using the below IDNLGREY model structure, with noise
% added to the output. The twotankdata.mat file contains one data set with
% 3000 input-output samples, generated using a sampling rate (Ts) of 0.2
% seconds. The input u(t) is the voltage [V] applied to a pump, which
% generates an inflow to the upper tank. A rather small hole at the bottom
% of this upper tank yields an outflow that goes into the lower tank, and
% the output y(t) of the two tank system is then the liquid level [m] of
% the lower tank. We create an IDDATA object z to hold the tank data. For
% bookkeeping and documentation purposes we also specify channel names and
% units. This step is optional.
load(fullfile(matlabroot, 'toolbox', 'ident', 'iddemos', 'data', 'twotankdata'));
z = iddata(y, u, 0.2, 'Name', 'Two tanks');
set(z, 'InputName', 'Pump voltage', 'InputUnit', 'V',             ...
       'OutputName', 'Lower tank water level', 'OutputUnit', 'm', ...
       'Tstart', 0, 'TimeUnit', 's');

%%
% The input-output data that will be used for estimation are shown in a
% plot window.
figure('Name', [z.Name ': input-output data']);
plot(z);

%%
% *Figure 2:* Input-output data from a two tank system.

%% Modeling the Two Tank System
% The next step is to specify a model structure describing the two tank
% system. To do this, let x1(t) and x2(t) denote the water level in the
% upper and the lower tank, respectively. For each tank, fundamental
% physics (mass balance) states that the change of water volume depends on
% the difference between in- and outflow as (i = 1, 2):
%
%    d/dt (Ai*xi(t)) = Qini(t) - Qouti(t)
%
% where Ai [m^2] is the cross-sectional area of tank i and Qini(t) and
% Qouti(t) [m^3/s] are the inflow to and the outflow from tank i at time t.
%
% For the upper tank, the inflow is assumed to be proportional to the
% voltage applied to the pump, i.e., Qin1(t) = k*u(t). Since the outlet
% hole of the upper tank is small, Bernoulli's law can be applied, stating
% that the outflow is proportional to the square root of the water level,
% or more precisely that:
%
%    Qout1(t) = a1*sqrt(2*g*x1(t))
%
% where a1 is the cross-sectional area of the outlet hole and g is the
% gravity constant. For the lower tank, the inflow equals the outflow from
% the upper tank, i.e., Qin2(t) = Qout1(t), and the outflow is given by
% Bernoulli's law:
%
%    Qout2(t) = a2*sqrt(2*g*x2(t))
%
% where a2 is the cross-sectional area of the outlet hole.

%%
% Put altogether these facts lead to the following state-space structure:
%
%    d/dt x1(t) = 1/A1*(k*u(t) - a1*sqrt(2*g*x1(t)))
%    d/dt x2(t) = 1/A2*(a1*sqrt(2*g*x1(t)) - a2*sqrt(2*g*x2(t)))
%          y(t) = x2(t)

%% Two Tank C MEX Model File
% These equations are next put into a C MEX-file with 6 parameters (or
% constants), A1, k, a1, g, A2 and a2. The C MEX-file is normally a bit
% more involved than the corresponding file written using MATLAB language,
% but C MEX modeling generally gives a distinct advantage in terms of
% execution speed, especially for more complex models. A template C
% MEX-file is provided (see below) to help the user to structure the code.
% For most applications, it suffices to define the number of outputs and to
% enter the code lines that describe dx and y into this template. An
% IDNLGREY C MEX-file should always be structured to return two outputs:
%
%    dx: the right-hand side(s) of the state-space equation(s)
%     y: the right-hand side(s) of the output equation(s)
%
% and it should take 3+Npo(+1) input arguments specified as follows:
%
%     t: the current time
%     x: the state vector at time t ([] for static models)
%     u: the input vector at time t ([] for time-series models)
%     p1, p2, ..., pNpo: the individual parameters (which can be real
%        scalars, column vectors or 2-dimensional matrices); Npo is here
%        the number of parameter objects, which for models with scalar
%        parameters coincide with the number of parameters Np
%     FileArgument: optional inputs to the model file
%
% In our two tank system there are 6 scalar parameters and hence the number
% of input arguments to the C MEX modeling file should be 3+Npo = 3+6 = 9.
% The trailing 10:th argument can here be omitted as no optional
% FileArgument is employed in this application.

%%
% Writing a C MEX modeling file is normally done in four steps:
%
%    1. Inclusion of C-libraries and definitions of the number of outputs.
%    2. Writing the function computing the right-hand side(s) of the state
%       equation(s), compute_dx.
%    3. Writing the function computing the right-hand side(s) of the output
%       equation(s), compute_y.
%    4. Writing the main interface function, which includes basic error
%       checking functionality, code for creating and handling input and
%       output arguments, and calls to compute_dx and compute_y.
%
%  Let us view the C MEX source file (except for some comments) for the two
%  tank system and based on this discuss these four items in some more
%  detail.
%
% <<../twotanksmexcode.png>>
%
% *Figure 3:* C MEX source code for the two tank system.

%%
% 1. Two C-libraries mex.h and math.h are normally included to provide
% access to a number of MEX-related as well as mathematical functions. The
% number of outputs is also declared per modeling file using a standard
% C-define:
%
%    /* Include libraries. */
%    #include "mex.h"
%    #include "math.h"
%
%    /* Specify the number of outputs here. */
%    #define NY 1

%%
% 2-3. Next in the file we find the functions for updating the states,
% compute_dx, and the output, compute_y. Both these functions hold argument
% lists, with the output to be computed (dx or y) at position 1, after
% which follows all variables and parameters required to compute the
% right-hand side(s) of the state and the output equations, respectively.
%
% The first step in these functions is to unpack the model parameters that
% will be used in the subsequent equations. Any valid variable name (except
% for those used in the input argument list) can be used to provide
% physically meaningful names of the individual parameters.
%
% As is the case in C, the first element of an array is stored at position
% 0. Hence, dx[0] in C corresponds to dx(1) in MATLAB(R) (or just dx in case
% it is a scalar), the input u[0] corresponds to u (or u(1)), the parameter
% A1[0] corresponds to A1, and so on.
%
% The two tank model file involves square root computations. This is
% enabled through the inclusion of the mathematical C library math.h. The
% math library realizes the most common trigonometric functions (sin, cos,
% tan, asin, acos, atan, etc.), exponential (exp) and logarithms (log,
% log10), square root (sqrt) and power of functions (pow), and absolute
% value computations (fabs). The math.h library must be included whenever
% any math.h function is used; otherwise it can be omitted. See "Tutorials
% on Nonlinear Grey Box Model Identification: Creating IDNLGREY Model
% Files" for further details about the C math library.

%%
% 4. The main interface function should almost always have the same
% content and for most applications no modification whatsoever is needed.
% In principle, the only part that might be considered for changes is where
% the calls to compute_dx and compute_y are made. For static systems, one
% can leave out the call to compute_dx. In other situations, it might be
% desired to only pass the variables and parameters referred in the state
% and output equations. For example, in the output equation of the two tank
% system, where only one state is used, one could very well shorten the
% input argument list to:
%
%    void compute_y(double *y, double *x)
%
% and call compute_y as:
%
%    compute_y(y, x);
%
% The input argument lists of compute_dx and compute_y might also be
% extended to include further variables inferred in the interface function,
% like the number of states and the number of parameters.
%
% Once the model source file has been completed it must be compiled, which
% can be done from the MATLAB command prompt using the mex command; see
% "help mex". (This step is omitted here.)

%%
% When developing model specific C MEX-files it is often useful to start
% the work by copying the IDNLGREY C MEX template file. This template
% contains skeleton source code as well as detailed instructions on how to
% customize the code for a particular application. The location of the
% template file is displayed by typing the following at the MATLAB command
% prompt.
%
%    fullfile(matlabroot, 'toolbox', 'ident', 'nlident', 'IDNLGREY_MODEL_TEMPLATE.c')
%
% Also see "Creatin IDNLGREY Model Files" example for more details on
% IDNLGREY C MEX model files.

%% Creating a Two Tank IDNLGREY Model Object
% The next step is to create an IDNLGREY object describing the two tank
% system. For convenience we also set some bookkeeping information about
% the inputs and outputs (name and units).
FileName      = 'twotanks_c';          % File describing the model structure.
Order         = [1 1 2];               % Model orders [ny nu nx].
Parameters    = {0.5; 0.0035; 0.019; ...
                 9.81; 0.25; 0.016};   % Initial parameters.
InitialStates = [0; 0.1];              % Initial initial states.
Ts            = 0;                     % Time-continuous system.
nlgr = idnlgrey(FileName, Order, Parameters, InitialStates, Ts, ...
                'Name', 'Two tanks');
set(nlgr, 'InputName', 'Pump voltage', 'InputUnit', 'V',             ...
          'OutputName', 'Lower tank water level', 'OutputUnit', 'm', ...                          ...
          'TimeUnit', 's');

%%
% We continue to add information about the names and the units of the
% states and the model parameters via the commands SETINIT and SETPAR.
% Furthermore, both states x1(t) and x2(t) are tank levels that cannot be
% negative, and thus we also specify that x1(0) and x2(0) >= 0 via the
% 'Minimum' property. In fact, we also know that all model parameters ought
% to be strictly positive. We therefore set the 'Minimum' property of all
% parameters to some small positive value (eps(0)). These settings implies
% that constraint estimation will be carried out in the upcoming estimation
% step (i.e., the estimated model will be a model such that all entered
% constraints are honored).
nlgr = setinit(nlgr, 'Name', {'Upper tank water level' 'Lower tank water level'});
nlgr = setinit(nlgr, 'Unit', {'m' 'm'});
nlgr = setinit(nlgr, 'Minimum', {0 0});   % Positive levels!
nlgr = setpar(nlgr, 'Name', {'Upper tank area'        ...
                      'Pump constant'          ...
                      'Upper tank outlet area' ...
                      'Gravity constant'       ...
                      'Lower tank area'        ...
                      'Lower tank outlet area'});
nlgr = setpar(nlgr, 'Unit', {'m^2' 'm^3/(s*V)' 'm^2'  'm/(s^2)' 'm^2' 'm^2'});
nlgr = setpar(nlgr, 'Minimum', num2cell(eps(0)*ones(6,1)));   % All parameters > 0!

%%
% The cross-sectional areas (A1 and A2) of the two tanks can rather
% accurately be determined. We therefore treat these and g as constants and
% verify that the 'Fixed' field is properly set for all 6 parameters
% through the command GETPAR. All in all, this means that 3 of the model
% parameters will be estimated.
nlgr.Parameters(1).Fixed = true;
nlgr.Parameters(4).Fixed = true;
nlgr.Parameters(5).Fixed = true;
getpar(nlgr, 'Fixed')

%% Performance of the Initial Two Tank Model
% Before estimating the free parameters k, a1 and a2 we simulate the
% system using the initial parameter values. We use the default
% differential equation solver (a Runge-Kutta 45 solver with adaptive step
% length adjustment) and set the absolute and relative error tolerances to
% rather small values (1e-6 and 1e-5, respectively). Notice that the
% COMPARE command, when called with two input arguments, as default will
% estimate all initial state(s) regardless of whether any initial state has
% been defined to be 'Fixed'. In order to only estimate the free initial
% state(s), call COMPARE with a third and a fourth input argument as
% follows: compare(z, nlgr, 'init', 'm'); as both initial states of the
% tank model by default are 'Fixed', no initial state estimation will be
% performed by this command. 
nlgr.SimulationOptions.AbsTol = 1e-6;
nlgr.SimulationOptions.RelTol = 1e-5;

compare(z, nlgr);

%%
% *Figure 4:* Comparison between true output and the simulated output
% of the initial two tank model.

%%
% The simulated and true outputs are shown in a plot window, and as can be
% seen the fit is not so impressive.

%% Parameter Estimation
% In order to improve the fit, the 3 free parameters are next estimated
% using NLGREYEST. (Since, by default, the 'Fixed' fields of all initial
% states are false, no estimation of the initial states will be done in
% this call to the estimator.)
nlgr = nlgreyest(z, nlgr, nlgreyestOptions('Display', 'on'));

%% Performance of the Estimated Two Tank Model
% To investigate the performance of the estimated model, a simulation of it
% is performed (the initial states are here reestimated).
compare(z, nlgr);

%%
% *Figure 5:* Comparison between true output and the simulated output of
% the estimated two tank model.

%%
% The agreement between the true and the simulated outputs is quite good. A
% remaining question is, however, if the two tank system can be accurately
% described using a simpler and linear model structure. To answer this, let
% us try to fit the data to some standard linear model structures, and then
% use COMPARE to see how well these models capture the dynamics of the
% tanks.
nk = delayest(z);
arx22 = arx(z, [2 2 nk]);   % Second order linear ARX model.
arx33 = arx(z, [3 3 nk]);   % Third order linear ARX model.
arx44 = arx(z, [4 4 nk]);   % Third order linear ARX model.
oe22  = oe(z,  [2 2 nk]);   % Second order linear OE model.
oe33  = oe(z,  [3 3 nk]);   % Third order linear OE model.
oe44  = oe(z,  [4 4 nk]);   % Fourth order linear OE model.
sslin = ssest(z);           % State-space model (order determined automatically)
compare(z, nlgr, 'b', arx22, 'm-', arx33, 'm:', arx44, 'm--', ...
        oe22, 'g-', oe33, 'g:', oe44, 'g--', sslin, 'r-');

%%
% *Figure 6:* Comparison between true output and the simulated outputs of a
% number of estimated two tank models.

%%
% The comparison plot clearly reveals that the linear models cannot pick up
% all dynamics of the two tank system. The estimated nonlinear IDNLGREY
% model on the other hand shows an excellent fit to the true output.
% In addition, the IDNLGREY model parameters are also well in line with
% those used to generate the true output. In the following display
% computations, we are using the command GETPVEC, which returns a parameter
% vector created from the structure array holding the model parameters of
% an IDNLGREY object.
disp('   True      Estimated parameter vector');
ptrue = [0.5; 0.005; 0.02; 9.81; 0.25; 0.015];
fprintf('   %1.4f    %1.4f\n', [ptrue'; getpvec(nlgr)']);

%%
% The prediction errors obtained using PE are small and look very much like
% random noise.
figure;
pe(z, nlgr);

%%
% *Figure 7:* Prediction errors obtained for the estimated IDNLGREY two
% tank model.

%%
% Let us also investigate what happens if the input voltage is increased
% from 5 to 6, 7, 8, 9 and 10 V in a step-wise manner. We do this by
% calling STEP with different specified step amplitudes starting from a
% fixed offset of 5 Volts. The step response configuration is facilitated
% by a dedicated option-set created by |stepDataOptions|:
figure('Name', [nlgr.Name ': step responses']);
t = (-20:0.1:80)';
Opt = stepDataOptions('InputOffset',5,'StepAmplitude',6);
step(nlgr, t, 'b', Opt);
hold on

Opt.StepAmplitude = 7;  step(nlgr, t, 'g', Opt);
Opt.StepAmplitude = 8;  step(nlgr, t, 'r', Opt);
Opt.StepAmplitude = 9;  step(nlgr, t, 'm', Opt);
Opt.StepAmplitude = 10; step(nlgr, t, 'k', Opt);

grid on;
legend('5 -> 6 V', '5 -> 7 V', '5 -> 8 V', '5 -> 9 V', '5 -> 10 V', ...
       'Location', 'Best');
%%
% *Figure 8:* Step responses obtained for the estimated IDNLGREY two
% tank model.

%%
% By finally using the PRESENT command, we get summary information about
% the estimated model:
present(nlgr);

%% Conclusions
% In this example we have shown:
%
%    1. how to use C MEX-files for IDNLGREY modeling, and
%    2. provided a rather simple example where nonlinear state-space
%       modeling shows good potential

%% Additional Information
% For more information on identification of dynamic systems with System Identification Toolbox(TM) 
% visit the
% <http://www.mathworks.com/products/sysid/ System Identification Toolbox>
% product information page.