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

    %% Linear Approximation of Complex Systems by Identification
%
% This example shows how to obtain linear approximations of a complex,
% nonlinear system by means of linear model identification. The approach is
% based on selection of an input signal that excites the system. A linear
% approximation is obtained by fitting a linear model to the simulated
% response of the nonlinear model for the chosen input signal.
%
% This example uses Simulink(R), Control System Toolbox(TM) and Simulink
% Control Design(TM). 

% Copyright 2011-2015 The MathWorks, Inc.

%% Introduction
% In many situations, a linear model is obtained by simplification of a
% more complex nonlinear system under certain local conditions. For
% example, a high fidelity model of the aircraft dynamics may be described
% by a detailed Simulink model. A common approach taken to speed up the
% simulation of such systems, study their local behavior about an operating
% point or design compensators is to linearize them. If we perform the
% linearization of the original model analytically about an operating
% point, this, in general, would yield a model of an order as high as, or
% close to, the number of states in the original model. This order may be
% unnecessarily high for the class of inputs it is supposed to be used with
% for analysis or control system design. Hence we can consider an
% alternative approach centered around collecting input-output data from
% the simulation of the system and using it to derive a linear model of
% just the right order.

%% Analytical Linearization of F14 Model
% Consider the F14 model. This is a already a linear model but contains
% derivative blocks and sources of disturbance that may affect the nature
% of its output. We can "linearize" it between its one input port and the
% two output ports as follows:
open_system('idF14Model')
IO = getlinio('idF14Model')
syslin = linearize('idF14Model',IO)

%% 
% syslin is a model with 2 outputs, one input and 9 states. This is because
% the linearization path from the "Stick Input" input to the two outputs
% has 9 states in the original system. We can verify this by using
% |operpoint|:
operpoint('idF14Model')

%%
% Could this order be reduced while still maintaining fidelity to the
% responses for the chosen square-wave ("Stick input") input?

%% Preparing Identification Data
% Simulate the model and log the input square wave (u) and the outputs
% "Angle of attack" (y1) and "Pilot G force" (y2) for 0:30 second time
% span. This data, after interpolation to a uniformly spaced time vector
% (sample time of 0.0444 seconds), is stored in the "idF14SimData.mat"
% file.
load idF14SimData
Z = iddata([y1, y2],u,Ts,'Tstart',0);
Z.InputName = 'Stick input'; Z.InputUnit = 'rad/s';
Z.OutputName = {'Angle of attack', 'Pilot G force'};
Z.OutputUnit = {'rad', 'g'};
t = Z.SamplingInstants;

subplot(311)
plot(t,Z.y(:,1)), ylabel('Angle of attack (rad)') 
title('Logged Input-Output Data')

subplot(312)
plot(t,Z.y(:,2)), ylabel('Pilot G force (g)')

subplot(313)
plot(t,Z.u), ylabel('Stick input (rad/s)')
axis([0 30 -1.2 1.2])
xlabel('Time (seconds)')

%% Estimation of State Space Models 
% Estimate state-space models of orders 2 to 4 using the |ssest| command.
% We configure the estimation to use "simulation" focus and choose to not
% estimate the disturbance component of the model. 
opt = ssestOptions('Focus','simulation');
syslin2 = ssest(Z, 2, 'DisturbanceModel', 'none', opt);
syslin3 = ssest(Z, 3, 'DisturbanceModel', 'none', opt);
syslin4 = ssest(Z, 4, 'DisturbanceModel', 'none', opt);

%%
% Compare the performance of the linearized model |syslin| and the three
% identified models to the data. Note that syslin is an SS model while
% |syslin2|, |syslin3| and |syslin4| are IDSS models. 
syslin.InputName = Z.InputName;
syslin.OutputName = Z.OutputName; % reconcile names to facilitate comparison
clf
compare(Z, syslin, syslin2, syslin3, syslin4)

%%
% The plot shows that the 3rd order model (|syslin3|) works pretty well as
% a linear approximation of aircraft dynamic about its default (t=0)
% operating conditions. It fits the data slightly better than the one
% obtained by analytical linearization (|syslin|). If the original
% idF14Model is linear, why doesn't its linearization result, |syslin|,
% produce a 100% fit to the data? There are two reasons for this:
% 
% # The measured outputs are affected by the wind gusts which means that
%   the logged outputs are not simply a function of the Stick input. There
%   are disturbances affecting it.
% # The "Pilot G force" block uses derivative blocks whose linearization
%   depends upon the value of the time constant "c". "c" should be small
%   (we use 1e-5) but not zero. A non-zero value for c introduces an
%   approximation error in the linearization.
% 
% Let us view the parameters of the model |syslin3| which seems to capture
% the responses pretty well: 
syslin3

%% Model Simplification by Reduction and Estimation
% We can also take the approach of reducing the order of the linearized
% model |syslin| and refining the parameters of the reduced model to best
% fit the data Z. To figure out a good value for the reduced order, we use
% |hsvd|:
[S, BalData] = hsvd(syslin);
clf; bar(S)

%%
% The bar chart shows that the singular values are quite small for states 4
% and above. Hence 3rd order might be optimal for reduction.
sysr = balred(syslin,3,BalData)
opt2 = bodeoptions; opt2.PhaseMatching = 'on';
clf; bodeplot(sysr,syslin,opt2)

%%
% The bode plot shows good fidelity up to 10 rad/s. |sysr| is able to
% emulate the responses as good as the original 9-state model as shown by
% the |compare| plot:
compare(Z, sysr, syslin)

%%
% Let us refine the parameters of |sysr| to improve its fit to data. For
% this estimation, we choose the "Levenberg-Marquardt" search method and
% change the maximum number of allowable iterations to 10. These choices
% were made based on some trial and error. We also turn on the display of
% estimation progress.
opt.Display = 'on'; 
opt.SearchMethod = 'lm';
opt.SearchOption.MaxIter = 10;
sysr2 = ssest(Z, sysr, opt)
compare(Z, sysr2)

%%
% The refined model |sysr2| fits the response of the F14 model quite well
% (about 99% for the first output, and about 97% for the second).

%% Conclusions
% We showed an alternative approach to analytical linearization for
% obtaining linear approximations of complex systems. The results are
% derived for a particular input signal and, strictly speaking, are
% applicable to that input only. To improve the applicability of results to
% various input profiles, we could perform several simulations using the
% various types of inputs. We could then merge the resulting datasets into
% one multi-experiment dataset (see |iddata/merge|) and use it for
% estimations. In this example, we used a complex, but linear system for
% convenience. The real benefit of this approach would be seen for
% nonlinear systems.

%%
% We also showed an approach for reducing the order of a linear
% system while keeping the reduced model faithful to the simulated response
% of the original Simulink model. The role of the reduced model |sysr| was to
% provide an initial guess for the estimated model |sysr2|. The approach
% also highlights the fact that any linear system, including those of a
% different class, can be used as the initial model for estimations.

bdclose('idF14Model')

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