Data Driven Model Estimation

A system identification package for predicting time-series and to enhance dynamical models from data

Jonas Sjöberg, Chalmers University, Sweden

Håkan Hjalmarsson, Royal Institute of Technology, Sweden

The coming 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.

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

First load packages

\(\($SystemIdentificationDirectory = "c:/users/Home/sjoberg/Mathematica/developments/TechnologyConference";\)\)
\(\(SetDirectory[
        ToFileName[{$SystemIdentificationDirectory, "systemID2006"}]];\)\)
<< ControlSystems`;
<< NumericalControlMethods`
<< NeuralNetworks`;
<< Statistics`ContinuousDistributions`;
<< NumericalControlMethods`Examples`
<< Statistics`NormalDistribution`

    
<< SIfiles.m
<< SIFit.m
<< package.m
<< packageInternal.m
<< InputOutputModels.m
<< blocks.m
<< wiener.m

Functions under development

observabilityForm::usage = "observabilityForm[aa,bb,cc,dd] removes \
unobservable states and returns the system matrices on observability form. \
Works only for single output systems.";

\(observabilityForm[aa_, bb_, cc_, dd_] := 
      Module[{a2, b2, c2, d2, om, tt, ttinv, pars, par, 
          temp}, 

Print["ObservabilityForm is being used"]; 

 (*\ remove\ non - observable\ states\ *) 
{a2, b2, c2, 
            d2} = controlMinimalForm[Transpose[aa], Transpose[cc], 
            Transpose[bb], Transpose[dd]]; 
{a2, c2, b2} = 
          Map[Transpose, {a2, b2, 
              c2}]; 

 (*find\ parameter\ symbols*) 
temp = {}; 
Map[
          If[Head[#] === Symbol, temp = Append[temp, #]] &, {aa, bb, cc, dd}, 
          Infinity]; 
pars = 
          Union[temp]; 
 (*\ 
          dummy\ parameter\ values\ for\ the\ observability\ matrix\ *) 
par = 
          E*N[Range[Length[pars]]]; 
{a3, c3} = {a2, c2} /. 
            Inner[Rule, pars, par, List]; 
om = 
          observabilityMatrix[a3, 
            c3]; 
tt = \((AppendRows @@ 
              NestList[a3\  . \ # &, \(Inverse[om]\)[\([All, {\(-1\)}]\)], 
                Length[a2] - 
                  1])\); 

ttinv = 
          Inverse[tt]; 
 (*\ \(Print[
              Map[MatrixForm, {a2, b2, c2, 
                  dd}]];\)\ *) 
Chop[
            Expand[{ttinv\  . \ a2\ \  . \ tt, ttinv\  . \ b2, c2\  . \ tt, 
                dd}]] /. {\ Times[1.`, a_] -> a, \ 
            Times[\(-1.`\), 
                a_] -> \(-a\)}
]\)

New Functions

sIFit[mod_LinearStateSpace, rest___] := 
      Module[{mod2 = mod, modnl, 
          cov}, 

mod2[[2]] = 
          Fold[DeleteCases[#1, #2] &, 
            mod2[[2]], {(NoiseVariance -> _), (ParameterCovariance -> _)}]; 

modnl = 
          sIFit[linearStateSpace2NonlinearForm[mod2], 
            rest]; 

 (* 
          insert the parameters *) 
mod2 = 
          mod2 /. (Parameters -> _) -> (Parameters -> (Parameters /. 
                      modnl[[2]])); 

cov = 
          ParameterCovariance /. 
            modnl[[2]]; 
mod2[[2]] = 
          Fold[Insert[#1, #2, -1] &, 
            mod2[[2]], {ParameterCovariance -> cov, 
              NoiseVariance -> (NoiseVariance /. 
                    modnl[[2]])}]; 

mod2
]

sIFit[mod_BlockStateSpace, rest___] := 
      Module[{mod2 = mod, modnl, 
          cov}, 

mod2[[2]] = 
          Fold[DeleteCases[#1, #2] &, 
            mod2[[2]], {(NoiseVariance -> _), (ParameterCovariance -> _)}]; 

modnl = 
          sIFit[releaseBlockStructure[mod2], 
            rest]; 

 (* 
          insert the parameters *) 
mod2 = 
          mod2 /. (Parameters -> _) -> (Parameters -> (Parameters /. 
                      modnl[[2]])); 

cov = 
          ParameterCovariance /. 
            modnl[[2]]; 
mod2[[2]] = 
          Fold[Insert[#1, #2, -1] &, 
            mod2[[2]], {ParameterCovariance -> cov, 
              NoiseVariance -> (NoiseVariance /. 
                    modnl[[2]])}]; 

mod2
]

Predict[mod_BlockStateSpace, rest___] := 
      Predict[releaseBlockStructure[mod], rest]

predictionError[mod_BlockStateSpace, rest___] := 
      predictionError[releaseBlockStructure[mod], rest]

comparePlot[mod_BlockStateSpace, rest___] := 
      comparePlot[releaseBlockStructure[mod], rest]

SimulateOutputs[mod_BlockStateSpace, rest___] := 
      SimulateOutputs[releaseBlockStructure[mod], rest]

Predict[mod_LinearStateSpace, rest___] := 
      Predict[linearStateSpace2NonlinearForm[mod], rest]

predictionError[mod_LinearStateSpace, rest___] := 
      predictionError[linearStateSpace2NonlinearForm[mod], rest]

comparePlot[mod_LinearStateSpace, rest___] := 
      comparePlot[linearStateSpace2NonlinearForm[mod], rest]

SimulateOutputs[mod_LinearStateSpace, rest___] := 
      SimulateOutputs[linearStateSpace2NonlinearForm[mod], rest]

Examples

Comparison with subspace algorithms  from the Numerical Control Methods package

From the package, a system model of a steam power system with one input , two outputs and 5 states.

System matrices, A, B, C and D

ss = SteamPowerSystem[];  TraditionalForm[ss]

It has the following Bode diagram

freq = {(10^(-2)), π}; BodePlot[ss, freq]

Generate input signal

We excite the system with a random (white) input signal of length 1000.

Ns = 1000;

u = 
        Table[Random[
            NormalDistribution[0, 1]], {CountInputs[ss]}, {Ns}];

Look at the input signal

SimulationPlot[u, 
        Sampled -> Period[SamplingPeriod[ss]]];

Generate output signal by simulating the system

yu = OutputResponse[ss, u];

Noise-free signal

SimulationPlot[yu, 
        Sampled -> Period[SamplingPeriod[ss]]];

Add noise to the output to imitate a measured signal

The system output is perturbed by white Gaussian noise with variance 1.

σ = 1;
e = σ Table[
            Random[NormalDistribution[0, 1]], {CountOutputs[
                ss]}, {Ns}];
y = yu + e;

SimulationPlot[y, 
        Sampled -> Period[SamplingPeriod[ss]]];

Estimate models using the data

We first estimate a model using the algorithm provided in Control System Professional:

CSPmodel = 
        OutputResponseIdentify[u, y, Sampled -> Period[1], 
          MaxOrder -> 5, Tolerance -> σ];

which warns us that it has problems estimating a good model using this data. Indeed, the model fit is poor: It is a first order model

TraditionalForm[CSPmodel]

Bode plot of CSP estimate

and the fit to the Bode diagram of the true system is poor:

DisplayTogetherGraphicsArray[BodePlot[ss, freq], 
      BodePlot[CSPmodel, freq, 
        PlotStyle -> {{Dashing[{ .02}]}, {Dashing[{ .01,  .025}], 
              Thickness[ .01]}}]]

The new subspace algorithm

z = PackSIData[u, y];
\(\({A\&^, B\&^, K\&^, C\&^, D\&^, Λ\&^, x\&^\_0, G\&^, 
          H\&^, e\&^} = 
        SubspaceIdentification[z, 5, {1}, 50, 50];\)\)

which has a much better model fit

\(\(DisplayTogetherGraphicsArray[BodePlot[ss, freq], 
        BodePlot[G\&^, freq, 
          PlotStyle -> {{Dashing[{ .02}]}, {Dashing[{ .01,  .025}], 
                Thickness[ .01]}}]];\)\)

Non-parametric model estimate of the steam power system

We first set the smoothing parameter to 30 (little smoothing)

NonParametricM = 
        SpectrumAnalysis[z, 30]; 
Bode[NonParametricM, PlotRange -> {{0, Pi}, All}]

which gives a rough estimate.

Smoother estimate

By increasing the smoothing (a lower smoothing parameter) we get a better estimate closer to the true system

NonParametricM = 
        SpectrumAnalysis[z, 10]; 
Bode[NonParametricM, PlotRange -> {{0, Pi}, All}]

Nonlinear identification of an hydrolic actuator

Load data

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

Plot the input signal.

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

Plot the output signal

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

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

Divide the data set into estimation and validation data

ue = u[[All, Range[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}}]

The lmodel can be used for prediction and simulation

Predict[mod1, SIDataPart[zv, {1, 20}], 
      Horizon -> Infinity]

Predict[mod1, SIDataPart[zv, {1, 20}], Horizon -> 1]

Compare model output with measured output

One-step-head prediction on validation data. For a fair validation, validation data is used.

Since the signals are dominated by low frequencies, the difference is hardly visible.

comparePlot[mod1, zv, Horizon -> 1]

Compare simulation with true output

A simulation test is typically tougher. The shortcomings of the lilnear model are clearly visible.

comparePlot[mod1, zv, Horizon -> Infinity]

Go on with nonlinear models

First, the linear ARX model is converted to a nonlinear model format.

mod2 = lBB2nlBB[mod1, CompoundBlock -> False]

The model is easier to crasp with a plot.

plotModel[mod2]

Add nonlinear parts

Since we handle dynamic models, special algorithms are used to initialize the addad nonlinear part

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

plotModel[mod3]

Compare initialized nonlinear model with true output

comparePlot[mod3, zv, Horizon -> Infinity]

Fit all parameters in the model

mod4 = sIFit[mod3, z, 30];

Compare the new model with true output

comparePlot[mod4, zv, Horizon -> Infinity]

An alternative nonlinear model

Take step back, add a nonlinearity to the linear model but include a subset of the regressors.

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

Only the delayed inputs feed into the nonlinearity.

plotModel[mod3]

Fit parameter and check the result

mod4 = sIFit[mod3, z, 30];

This model typically performs better than the previous one, but it depends on the, partly, randomized initialization of the nonlinear part.

comparePlot[mod4, zv, Horizon -> Infinity]

Try a Hammerstein model

Hammerstein models

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[[1]]]

Fit the Hammerstein model and test it on data

mod3 = sIFit[mod2, z, 30];

comparePlot[mod3, zv, Horizon -> Infinity]

Look at the estimated nonlinearity of the Hammerstein model

plotNonlinearity[mod3, z]

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