A Tool 4 Linear and Nonlinear System Identification

An AddOn package soon to be finished....

Jonas Sjöberg, Chalmers University, Sweden

Håkan Hjalmarsson, Royal Institute of Technology, Sweden

SIDebuggingMode=False

A System Identification package

This talk describes a system identification package which is under development. Previews of algorithms and supported model structures are given in connection with examples. The package supports algorithms for identification of linear and nonlinear models of dynamical systems which can be used for simulation, prediction, and for control design.

Package description

The symbolic features of Mathematica is used so that a large variety of model structures can be handled. Before the numerical estimation procedure, the symbolic definition of the model structure is simplified and modified before the numerical computations of the parameters. In this way a larger set of model structures can be supported, at the same time as the computational speed is improved.

Extends the possibilities Times Series package and the Neural Network package.

This talk describes ongoing work

Most of the algorithms work, at least on some examples....

GUI is missing, result are not presented in a beutiful way.

The system identification problem

Any application area where you work on data from systems with dynamics.

General model set-up

The problem is to find a model to predict future outputs

(OverHat[y])(t)=g((y^(t-1)),(u^(t-1)))

where (y^(t-1)) and (u^(t-1)) contain all the previous information.

({{(y^(t-1))=[y(t-1),y(t-2),…,y(0)]}, {(u^(t-1))=[u(t-1), u(t-2),…, u(0)]}})

Divide the problem into two parts:

A mapping {(y^(t-1)), (u^(t-1))}⟶x(θ,t)

A second mapping g:x(θ,t)⟶(OverHat[y])(t).

(OverHat[y])=g(θ,x)=(∑_(i=1)^n)(Subscript[c, i])(Subscript[g, i])(x,(Subscript[a, i]))
θ=[(Subscript[c, 1]) (Subscript[a, 1]) (Subscript[c, 2]) (Subscript[a, 2]) …(Subscript[c, n]) (Subscript[a, n])]

The nonlinear parameters (Subscript[a, i]) give the position of the basis functions.

Parameter fitting: minimization

Consider the parameter criterion - which can be defined by the user to a different measure

(Subscript[V, N])(θ)=(1/(2N))(∑_(t=1)^N)|y(t)-g(x(t),θ)(|^2)

and the minimization algorithm

((OverHat[θ])^((i+1)))=((OverHat[θ])^((i)))-(Subscript[μ, i])(Subscript[R, i]^(-1))(Subscript[V, N]^′)((θ^((i))))

((OverHat[θ])^((0))): initial parameter guess
(Subscript[μ, i]): step length to ensure descent updates
(Subscript[R, i]): matrix changing the search direction

Topics: Initial state estimation, stability of model: this is handled by th epackage

General algorithmic issue

Model structures are handled as symbolic expressions describing the equations.

Mathematica is used to generate the necessary derivatives, and to optimize the expressions before numerical computations start. This reduce the computational load.

Who is the tool intended for?

Assume you have obtained a fairly good linear model and that you want to go on and test some nonlinear alternatives: There is special support on extending linear models with nonlinearities.

You want to test many different nonlinear model structures.

You want to fit parameters of a model built on physical insights.

State-of-the-art support for linear system identification, but some new and better algorithms are included.

Tools exist supporting one, or a couple of nonlinear model types, ie, identification of Wiener models. However, in most practical situations where the features of the true system is not known, it more important to be able to test out a variety of model structures so that there is a fair chance to find a good model structures.

Nonlinear Parameterized Basis functions

y={{1}};
x={{1,2}};

fdfrwrd=InitializeFeedForwardNet[x,y,{1},RandomInitialization->True] ;fdfrwrd[[1]]={{{{-1,-1},{1,-1},{0,0}},{{1},{1},{0}}}};

Plot3D[Evaluate[fdfrwrd[{x1,x2}][[1]]],{x1,-8,8},{x2,-8,8},ViewPoint->{3.24,-0.038,0.975},PlotPoints->20]

The sum of two sigmoids in two dimensions, which is a feedforward neural network, can be seen as a local affine function with smooth transitions.

Initialization of basis functions

Incorrect approach (left): The smooth step of the sigmoid is outside the data domain or so sharp that few data points are placed on the transition.

Correct approach (right): The basis functions are placed so that the smooth step takes place within the data domain and is slow enough that many data points are placed on the transition region.

Nonlinear identification of an hydraulic actuator

Load data

<<actuator.dat;
y=Transpose[y];
u=Transpose[u];

Plot the input signal.

ListPlot[Flatten[u],Joined->True]

Plot the output signal

Resonance depends on amplitude - cannot be explained with a linear model

ListPlot[Flatten[y],Joined->True]

Divide the data set into estimation and validation data

ue=u[[All,;;512]];
ye=y[[All,Range[512]]];
uv=u[[All,Range[513,1024]]];
yv=y[[All,Range[513,1024]]];

Pack data into a data object

z=PackSIData[ue,ye];

zv=PackSIData[uv,yv];

Fit a linear ARX model

mod1=ARXFit[z,{{2}},{{2}},{{1}}]

Compare simulation with true output

comparePlot[mod1,zv,Horizon->Infinity]

Go on with nonlinear models

mod2=lBB2nlBB[mod1,CompoundBlock->False]

The model is easier to crasp with a plot.

plotModel[mod2]

Add nonlinear parts

mod3=addBlockParallel[mod2,{1,1},z,FeedForwardNet,{3}];

plotModel[mod3]

The nonlinear model will be at least as good as the nonlinear model it is based on.

Fit all parameters in the model

mod4=sIFit[mod3,z,20];

Compare the new model with true output

comparePlot[mod1,mod4,zv,Horizon->Infinity]

Output Error: an alternative nonlinear model

mod5=OEFit[z,{2},{2},{1}];
mod6=lBB2nlBB[mod5,CompoundBlock->False];
plotModel[mod6]

Here, to add nonlinearities, the linear model is needed to produce simuted model outputs: ys(t). Otherwize the basis function cannot be placed well.

Add nonlinearity, fit the model and evaluate

We can indicate which inputs to include in the nonlinearity.

mod7=addBlockParallel[mod6,{1,1},z,{2,3,4},FeedForwardNet,{3}];

plotModel[mod7]

mod8=sIFit[mod7,z,10];

comparePlot[mod1,mod4,mod6,mod8,zv,Horizon->Infinity]

Special model structures: Hammerstein model

Specify a nonlinearity consisting of a local linear function

mod2=makeHammersteinModel[mod1,z,LocalLinear,8,FreeNots->True];

This illustrates how Hammerstein models are handled in the package. Normally, this plot is not intended for the user.

plotModel[mod2]

Fit the Hammerstein model and test it on data

mod3=sIFit[mod2,z,10];

comparePlot[mod3,zv,Horizon->Infinity]

Look at the estimated nonlinearity

plotNonlinearity[mod3,z]

Linear dynamic part can be evaluated in normal ways: Bode plot, zero-pole plots etc.

More nonlinearities

Look at the model again

plotModel[mod2]

Look inside the compound block

plotModel[mod2,{1,1}]

Add a new nonlinearity

Add a new nonlinearity

mod3=addBlockAfter[mod2,{1,1},z,FeedForwardNet,{3}];

First level looks like before

plotModel[mod3]

The new nonlinear block is placed inside the compound block so that the Hammerstain structure is kept and only its nonlinearity is adapted.

plotModel[mod3,{1,1}]

The new part is initialized so that the new model becomes equal with the previous model. From that all points can be fitted and the new model can be evaluated.

This process of adding and fitting new nonlinearities can be continued as long as you want.

More special model structures:

Wiener model

Wiener-Hammerstein

Hammertstein-Wiener

And much more.

StateSpace models

Most general model format is the state space form:

Linear state space model:

bB2LinearStateSpace[mod5]

State space models can also be defined and manipulated in the block based form:

mod6=linearStateSpace2blockmodel[%,CompoundBlock->False];

plotModel[mod6]

GUI issues

Select linear ARX model with help

ARXSelection[{2,8},{1,8},{1,3},z,zv]

An advanced GUI is needed to help the user to work with all nonlinear possibilities.

Evaluate and validate obtained models

Models can be evaluated in various ways, ie,  linearize the nonlinear model along a trajectory to obtain a set of linear models describing the system+ uncerainty.

Linear models are evaluated with traditional measures.

Algorithmic issues

Estimate of initial state. Methods for linear models, and for nonlinear models as well based on linearization.

Resonable parameter initilization thanks to stepwize procedure but no guarantee against local minima.

Subsets of parameters can be estimated. It is possible to combine blocks containing hard-wired prior knowledge and black box blocks with free parameters.

Symbolic optimization before numerical computations. For example, if there are no parameters influencing the states, they are only simulated at the first iterartion.

Package overview

Non-parametric methods

Spectral estimation, Bode and Nyquist plots

Parameterized models

Linear models

Subspace methods for linear state-space models, very useful for multi-input multi-output systems
Traditional black box model structures: FIR, ARX, OE, ARMAX, BJ and for time-series AR and ARMA models

Non-linear models

The package support of estimation of nonlinear model structures is unprecedented: Nonlinear versions of the linear black box models are supported NLFIR,NLARX, NLOE, NLARMAX,NLBJ and for time-series NLAR and NLARMA models, nonlinear state-space models, Hammerstein and Wiener models.

The user can specify the nonlinearities themselves, but guidance and support is available. The models are estimated by the so called prediction error method.

Conclusions

A new system identification package supporting estimation of dynamic models using data

some standard algorithms for linear systems

some new, improved, algorithms for linear systems

a lot of new possibilities for nonlinear systems

In general, there are problems with local minima. Different initial parameter guesses give different results.

•	Sjöberg, Jonas: On Estimation of Nonlinear Black-Box Models: How to Obtain a Good Initialization. IEEE Workshop on Neural Networks for Signal Processing, Amelia Island Plantation, FL, USA, pp. 72-81. 
•	Sjöberg, Jonas; De Raedt, P.: Nonlinear System Identification: a Software Concept and Examples. Proc. IFAC Symposium on Identification, Fukuoka, Japan.