mlgreyest

Identify a multilinear time-invariant (MTI) state-space model from measured data.

Source: src/mlgreyest.m

Overview

mlgreyest(data,r,structuralRestrictions,mtiBase,options) fits a discrete-time explicit multilinear time-invariant (eMTI) model

$$\mathbf{x}(k+1) = \left\langle {\left[{\!\left[{\mathbf{S}_{F},\mathbf{\Phi}_{F}}\right]\!}\right]}\bigg | {\left[ \left(\begin{matrix}1-bv_1\\ v_1\end{matrix}\right), \dots, \left(\begin{matrix}1-bv_{n_v}\\ v_{n_v}\end{matrix}\right) \right] } \right\rangle$$

$$\mathbf{y} = \left\langle{\left[{\!\left[{\mathbf{S}_{G},\mathbf{\Phi}_{G}}\right]\!}\right]}\bigg | {\left[ \left(\begin{matrix}1-bv_1\\ v_1\end{matrix}\right), \dots, \left(\begin{matrix}1-bv_{n_v}\\ v_{n_v}\end{matrix}\right) \right] } \right \rangle$$

to measured data. The model parameters are stored in structure matrices $\mathbf{S}_F\in[-1,1]^{{n_v} \times R_F},\mathbf{S}_G\in[-1,1]^{{n_v}\times R_G}$, and in the parameter matrices $\mathbf{\Phi}_F\in\mathbb{R}^{n_x \times R_F},\mathbf{\Phi}_G\in\mathbb{R}^{n_y \times R_G}$.

The identification result is returned as an mss object where the mtiTensor is a CPNTensor. If the model base $b\in\{0,1\}$ is not passed explicitly to mlgreyest, the default mtiBase of the returned model is the monomial base ($b=0$).

In an mss, the CPNTensor is stored as one structure matrix $\mathbf{S}$ is $\left(n_{x}+n_{u}\right)\times R$ and one parameter matrix $\mathbf{\Phi}$ is $\left(n_{x}+n_{y}\right)\times R$ for the whole model, such that $R\leq R_F+R_G$.

Calling mlgreyest without any Name-Value options (default) will identify values for $\mathbf{\Phi}_F$ and/or $\mathbf{S}_F$ from measured states.

Data modes

State-trajectory mode (default). data.y is the measured state trajectory, so the model order equals the number of measured channels and the initial state is data.y(1,:). Only a state equation is identified. Focus selects whether the freely simulated trajectory ("simulation") or the one-step-ahead prediction from the measured states ("prediction") is fitted.

Output-identification mode (idFromOutputs=true). data.y are output measurements and the state is hidden: it is simulated forward from InitialState, whose length defines the state dimension. A state equation and an output equation are identified. Only the nonlinear methods support this mode, and only Focus="simulation" — without measured states there is nothing to predict from, so a different Focus is overridden with a warning.

Structural restrictions

Two different things are accepted in the third position as structuralRestrictions.

• A numeric or logical matrix is a completely known structure matrix $\mathbf{S}_F$. Only the parameters in $\mathbf{\Phi}_F$ are then computed.

• A CPNTensor restricts the nonlinear, output-identification path per element (only for the identification of a literal-base model). Its six split channels are read as follows:

Syntax

Input Arguments

data

• Type: iddata
• Description: Time-domain measurements. data.u holds the inputs, data.y the measured states (default) or the measured outputs (options.idFromOutputs=true). The ALS requires more than 2*r samples.

r

• Type: scalar, [rStateEq rOutputEq], or []
• Description: CP rank, i.e. the number of product terms.

structuralRestrictions

• Type: [] (default), numeric/logical matrix, or CPNTensor
• Description: See Structural restrictions above. A CPNTensor requires mtiBase = true, options.idFromOutputs = true and a nonlinear options.Method.

mtiBase

• Type: logical
• Default: false
• Description: 0 selects the monomial base, 1 the literal base.

options

Solver

Initialization and stopping

Bounds (nonlinear methods)

Output identification (Nonlinear solvers only)

Output Arguments

sys

• Type: mss
• Description: The identified model. The returned rank can be smaller than r, because trivialReduction is called internally after the identification.

J

• Type: vector (ALS) or scalar (nonlinear)
• Description: The cost. The ALS returns the cost of every accepted update in order, so J(end) is the final one and plot(J) shows the convergence; the nonlinear methods return the solver’s final cost.

x0Est

• Type: column vector or []
• Description: The initial state that was used. With estimateInitialState=true it is the value the solver found; otherwise it echoes InitialState, and in state-trajectory mode the first measured sample.

Examples

The examples below generate their own data with msim, so they can be run as they stand (the iddata constructor needs the System Identification Toolbox).

Identify structure and parameters with the ALS

% a rank-2 monomial model, 2 states and 1 input
S = [0.6 0.2;      % x1
     0.3 0.9;      % x2
     0.4 0.1];     % u1
P = [0.5 0.2;      % x1(k+1)
     0.1 0.4];     % x2(k+1)
sys0 = mss(S,P,[1 2],3,1,[1 2],[],'mtiBase',false);

t  = (0:59)';
u  = 0.5 + 0.4*sin(t/3);
[~,~,x] = sys0.msim(u,t,[0.3 0.6]);
data = iddata(x,u,1);                       % measured STATE trajectory

[sys,J] = mlgreyest(data,2,[],false,'Method',"als",'Display',"off");
fprintf('final cost %.3e\n', J(end));

Add 'Nonnegative',true to keep every identified value non-negative (each subproblem is then solved with lsqnonneg), or 'Focus',"prediction" to fit the one-step-ahead prediction instead of the free simulation.

Fit only the parameters of a known structure

sys = mlgreyest(data,[],S);             % S from the example above

Identify in the literal base

[sys,J] = mlgreyest(data,2,[],true,'Method',"als",'Initialize',"random");

Identify a state and an output equation from outputs only

% state equation on columns 1-2, output equation on columns 3-4
S = [ 1.0   0.0   0.6   1.0 ;   % x1
      0.0   1.0   1.0   0.0 ;   % x2
      0.4  -0.3   0.0   0.2 ];  % u1
P = [ 0.7   0.2   0     0   ;   % x1(k+1)
      0.3   0.5   0     0   ;   % x2(k+1)
      0     0     1.5   0.5 ;   % y1
      0     0     0.8   1.2 ];  % y2
sys0 = mss(S,P,[1 2],3,1,[1 2],[3 4],'mtiBase',true);

t  = (0:59)';
u  = 0.5 + 0.4*sin(t/3);
x0 = [0.3; 0.6];
y = sys0.msim(u,t,x0');
data = iddata(y,u,1);                       % measured OUTPUTS, state hidden

[sys,J,x0Est] = mlgreyest(data,[2 2],[],true, ...
    'Method',"lsqnonlin",'idFromOutputs',true,'InitialState',x0, ...
    'Display',"off",'MaxIterations',80,'MaxFunctionEvaluations',5000);

r = [2 2] requests disjoint column blocks. Use a scalar r to let both equations share the same columns, or name the columns explicitly:

[sys,J] = mlgreyest(data,[],[],true, ...
    'Method',"lsqnonlin",'idFromOutputs',true,'InitialState',x0, ...
    'columnIndexStateEq',[1 2],'columnIndexOutputEq',[2 3 4]);   % may overlap

Restrict the search with a CPNTensor

restriction = sys0.mtiTensor;               % same size as the model to be identified

[sys,J,x0Est] = mlgreyest(data,[],restriction,true, ...
    'Method',"lsqnonlin",'idFromOutputs',true,'InitialState',x0, ...
    'Focus',"simulation",'Display',"off",'MaxFunctionEvaluations',20000);

The tensor is also the warm start, so a restriction tensor that already carries good values converges in very few iterations.

To warm-start from different values on the same free-variable layout, pass a second tensor of the same size as Initialize:

perturbedTensor = copy(sys0.mtiTensor);
perturbedTensor.structureMatrixContinuous(perturbedTensor.structureMatrixContinuous>0)=0.8;
[sys,J] = mlgreyest(data,[],restriction,true, ...
    'Method',"lsqnonlin",'idFromOutputs',true,'InitialState',x0, ...
    'Initialize',perturbedTensor, 'Display', 'iter');

“A CPNTensor Initialize inherits every restriction of the restriction tensor” It supplies only values; which entries are free is decided by the CPNTensor structural restrictions, which it therefore requires.

Estimate the initial state as well

guess = [0.4; 0.4];                          % only a guess now
[sys,J,x0Est] = mlgreyest(data,[],restriction,true, ...
    'Method',"lsqnonlin",'idFromOutputs',true,'InitialState',guess, ...
    'estimateInitialState',true,'MaxFunctionEvaluations',20000);

The initial state becomes an extra block of optimization variables, bounded to $[0,1]$, and is returned as x0Est.

Practical notes

• In output mode x0 is generally not identifiable by itself. With a hidden state and factors that are affine in $x$, the model admits an affine state reparameterization: only the input–output map is pinned down, so a jointly estimated x0Est need not approach the true initial state even when the fit is excellent. Judge the result by the output fit, not by x0Est — unless the restrictions remove the symmetry.
• The returned rank can be smaller than r, because trivialReduction merges duplicate monomial columns and drops unused ones.

See Also

• mss — the identified model class
• msim — simulate the result
• CPNTensor — the tensor representation and the restriction format

MyToolbox Documentation | Generated automatically by CI/CD pipeline