Getting Started

This page takes you from a fresh installation to a model you have built, evaluated, simulated and linearized yourself — in about twenty lines of MATLAB. It assumes no prior contact with multilinear algebra.

Read Requirements and Installation first if the toolbox is not on your path yet. The explicit-model walkthrough below needs nothing beyond the toolbox itself, and none of the gated capabilities on the requirements page. The one exception is flagged where it arises: simulating the implicit model in the last section needs the Optimization Toolbox, which the requirements page lists as required rather than optional.

Once installed, confirm your environment before anything else:

environmentChecker

What a multilinear model is

A multilinear function is linear in each of its variables on its own, but not in all of them at once. The product x1u1x_{1}u_{1} is the smallest example: fix u1u_{1} and it is linear in x1x_{1}, fix x1x_{1} and it is linear in u1u_{1}, yet the function itself is not linear. Multilinear functions are a subclass of the polynomials and a superclass of both the linear and the Boolean functions, so a single representation covers continuous dynamics, switching logic and the mixture of the two.

Used as the right-hand sides of a state-space model, they give a multilinear time-invariant (MTI) model. The analogy to carry with you is this. A linear time-invariant model is a set of matrices, and linear algebra is the theory behind it. An MTI model is a set of tensors, and multilinear algebra is the theory behind it [1].

Stored in full, those tensors have one dimension of length two per variable, so their size doubles with every variable added. The toolbox never stores them that way. It keeps them decomposed, as sums of outer products — a canonical polyadic (CP) decomposition, normalized to the CPN form [3] used throughout — which is what makes models of realistic size tractable. Two decompositions are available, CPNTensor and TTTensor; you rarely handle either directly.

Two model classes sit on top of them:

The reason both exist is worth knowing early: the explicit class is not closed under composition. Connect two explicit MTI models and the result may be polynomial rather than multilinear. Implicit models are closed, so anything assembled from components — a network, a plant with its controller — is built as an mdss. Implicit modelling also means you need not decide in advance which variable each equation solves for, and conservation laws can be written as constraints rather than substituted away.

Your first model

Take the discrete-time bilinear model

x1(k+1)=0.9x1(k)+0.1x1(k)u1(k),y1(k)=2x1(k). x_{1}(k+1) = 0.9\,x_{1}(k) + 0.1\,x_{1}(k)\,u_{1}(k), \qquad y_{1}(k) = 2\,x_{1}(k) .

It is as small as a genuinely multilinear model gets: one state, one input, one output, and one term — x1u1x_{1}u_{1} — that no linear model can represent.

1. Build it

A model is written down as two matrices. The structure matrix says which variables appear in each term, one row per variable and one column per term; the parameter matrix says what each term is worth in each equation, one row per equation. There are three terms in total — x1x_{1}, x1u1x_{1}u_{1}, and x1x_{1} again for the output — but the output can reuse the first column, so two columns suffice:

%                   x1   x1*u1
structureMatrix = [  1     1 ;   % x1
                     0     1 ];  % u1

parameterMatrix = [ 0.9   0.1 ;  % x1(k+1)
                    2     0   ]; % y1

sys = mss(structureMatrix, parameterMatrix, ...
          1, ...   % stateIndex           (1st row of the structure matrix)
          2, ...   % inputIndex           (2nd row of the structure matrix)
          1, ...   % timeStepSize         (>0 = discrete-time)
          1, ...   % stateEquationIndex   (1st row of the parameter matrix)
          2);      % outputEquationIndex  (2nd row of the parameter matrix)

A structure entry of 0 means the variable is absent from that term, not that it is multiplied by zero. Column 1 is therefore x1x_{1} alone and column 2 is x1u1x_{1}u_{1}. The two index arguments in the middle map rows of the structure matrix onto states and inputs, and the last two map rows of the parameter matrix onto the state and output equations — the model does not infer them from the matrix shapes, because a model may be laid out in any order.

2. See what the constructor worked out

The output equation is nonzero only in column 1, so the constructor knows the two equations share that column and records which columns each one uses:

sys.nState                % 1
sys.nInput                % 1
sys.nOutput               % 1
sys.columnIndexStateEq    % [1; 2]
sys.columnIndexOutputEq   % 1        -- shared with the state equation

The two matrices you typed are not kept as such. They are folded into a CPNTensor, which splits the structure and the parameters into separate channels over disjoint supports; sys.structureMatrix and sys.parameterMatrix reconstruct a view of them when you ask.

3. Evaluate it at a point

The model is a function, and can be called as one — no simulation involved. functionValue evaluates the state equation and outputValue the output equation:

sys.functionValue(2, 3)   % 2.4   = 0.9*2 + 0.1*2*3
sys.outputValue(2, 3)     % 4     = 2*2

4. Simulate it

msim is the lsim of this toolbox. The model’s timeStepSize decides the scheme on its own — positive, as here, marches the difference equation sample by sample; zero would integrate with a MATLAB ODE solver. Give it an input trajectory with one column per input, the matching time vector, and an initial state:

k = (0:10)';
u = double(k >= 3);              % unit step at k = 3

[y, t, x] = msim(sys, u, k, 1);  % x0 = 1
plot(t, x, '.-'), grid on, xlabel('k'), ylabel('x_1')

The trajectories come back with one row per time point, in the order lsim uses: outputs first, then time, then states.

For the first three samples the input is zero and the state decays by a factor 0.9 per step, from 1 to 0.729. From k = 3 on it stops moving — it sits at 0.729 for the rest of the run.

5. See why it is not linear

That flat tail is the whole point of the example, and it is not a solver artefact. Factor the state equation:

x1(k+1)=(0.9+0.1u1(k))x1(k). x_{1}(k+1) = \bigl(0.9 + 0.1\,u_{1}(k)\bigr)\,x_{1}(k) .

The input does not merely drive this system, it changes the system’s own decay factor. At u1=0u_{1} = 0 the factor is 0.9 and the state decays; at u1=1u_{1} = 1 it is exactly 1 and the state is frozen; at u1=2u_{1} = 2 it is 1.1 and the state grows without bound. Simulate the same model with twice the step and you do not get twice the response — you get a different qualitative behaviour:

[~, ~, x2] = msim(sys, 2*u, k, 1);
x(end)     % 0.7290   -- frozen
x2(end)    % 1.4206   -- growing

linearize shows the same fact as numbers. It returns a standard MATLAB state-space object — a sparss, at the model’s own sample time — or the four matrices directly:

[A0, B0, C0, D0] = linearize(sys, 0.729, 0);   % A0 = 0.9   B0 = 0.0729
[A1, B1, C1, D1] = linearize(sys, 0.729, 1);   % A1 = 1.0   B1 = 0.0729

Two linearizations of one model, two different A matrices. That is the gap the toolbox exists to close: an MTI model carries the dependence itself, rather than being replaced by one linear approximation per operating point. The operating point need not be an equilibrium, and nothing checks that it is one.

The same model, written as equations

Typing structure matrices is how the representation works, not how you are expected to work. An implicit model can be parsed straight from its equations written as text, in residual form, with stringToMdss:

dsys = stringParser.stringToMdss( ...
    {"0 = xp1 - 0.9*x1 - 0.1*x1*u1"; ...
     "0 = y1  - 2*x1"}, ...
    1);                              % timeStepSize

The variable names carry the roles: xp1 is the state change x1(k+1)x_{1}(k+1), x1 a state, u1 an input, and y1 an algebraic variable — which is what an output becomes once causality is no longer assumed. All five role prefixes are configurable; z is the fifth, for binary signals.

mdss simulates through the same msim name, with a different signature, because an implicit model needs starting guesses for its algebraic and binary variables ([] asks the solver to find them):

Note: This step needs the Optimization Toolbox

Simulating an implicit model solves its algebraic sub-problems with lsqnonlin, so unlike everything above it needs the Optimization Toolbox — listed as required, not optional, on the Requirements page. If environmentChecker reported it missing, the rest of this page still works; this one command does not.

[yd, td, xd] = msim(dsys, u, k, 1);
max(abs(xd - x))    % 0

Zero to the last bit: the two models are the same model, reached two ways. Also available are symbolicToMdss for Symbolic Math Toolbox equations and polynomialToMdss.

What else the toolbox does

Most commands share, or nearly share, their Control System Toolbox syntax; some differ only by an m prefix, and several are overloaded methods with identical signatures. Beyond what you have used above:

Where MTI models are used

Application of MTI models began with HVAC systems [2], where projects have demonstrated real-time applicability. More recent work covers power networks, which opens the way to large-scale multi-energy systems. More information is on the project website, mti.systems.

One confusion is worth heading off, because the words are almost identical. Multilinear models are not multi-linear models. The latter are collections of linear models — one LTI model per operating point. They can be represented as tensors too, and within this toolbox, but they suffer the curse of dimensionality the decomposed representation is there to avoid.

Where to go next

References

  1. G. Lichtenberg (2012): Hybrid Tensor Systems, Habilitation, TU Hamburg.
  2. G. Pangalos, A. Eichler, G. Lichtenberg (2015): Hybrid Multilinear Modeling and Applications. https://doi.org/10.1007/978-3-319-11457-6_5
  3. K. Kruppa, G. Pangalos, G. Lichtenberg (2014): Multilinear approximation of nonlinear state space models. https://doi.org/10.3182/20140824-6-za-1003.00455
  4. L. Schnelle, G. Lichtenberg, C. Warnecke (2022): Using Low-rank Multilinear Parameter Identification for Anomaly Detection of Building Systems. https://doi.org/10.1016/j.ifacol.2022.07.173
  5. G. Lichtenberg, G. Pangalos, C. Cateriano Yáñez, A. Luxa, N. Jöres, L. Schnelle, C. Kaufmann (2022): Implicit multilinear modeling: An introduction with application to energy systems. https://doi.org/10.1515/auto-2021-0133

MyToolbox Documentation | Generated automatically by CI/CD pipeline