The Behavioral Toolbox

The data is a vector time series containing the observed variables of the system. Since these variables are not partitioned into inputs and outputs, which variables are inputs is not known. Also, how many variables are inputs is not known. It is possible however to fining the number of inputs from the given data. The function of the toolbox that does this is w2c. In addition to finding the number of inputs \(m\), it finds also the order \(n\) and the lag \(\ell\) of the underlying data-generating system.

Here is an example:

• specify the simulation parameters—number of variables q, number of inputs m, lag ell, order n, and number of data samples Td

q = 3; m = 1; ell = 2; n = 3; Td = 20; 

• generate a random linear time-invariant system with the specified parameters

Rp = randB(q, [m, ell, n]); 

• generate a random trajectory of the system with the specified number of samples

wd = R2w(Rp, q, Td); 

• call the w2c function

[~, mh, ellh, nh] = w2c(wd);

• check if the computed result is correct

[m == mh, ell == ellh, n == nh] % -> [true true true]  

But what about the assumption? A “sufficiently long” random trajectory almost surely fully specifies the system. Quantifying what “sufficiently long” means, however, requires knowledge of \(m\) and \(n\). There is no way out of the dilemma, either we have to assume that the data fully specifies the system, or we have to know \(m\) and \(n\), in which case we can check the condition.

Exercise: Change the simulation parameters and find out empirically when the data is “sufficiently long”.

We’ve seen that we can find the number of inputs \(m\) from the data. Can we find also the input variables \(u\)? No, in general, we can not do this because different sets of variables can act as inputs. The function of the toolbox that finds all possibilities of partitioning the variables into inputs and outputs is BT2IO.

BT2IO (as well as other functions of the toolbox) takes as an input argument a matrix \(B_T\) that can be obtained from the data using w2BT. \(B_T\) defines a non-parametric representation of the finite-horizon behavior \(\mathcal{B}|_T\). The horizon \(T\) is a problem dependent. For BT2IO the required horizon is \(T = 2\ell+1\).

Continuing the example above,

• first, we obtain the matrix \(B_T\):

T = 2 * ell + 1; BT = w2BT(wd, T, [m ell n]);

• then, we call BT2IO:

IO = BT2IO(BT, q)

Its output

IO =

     1     2     3
     1     3     2
     2     1     3
     2     3     1
     3     1     2
     3     2     1

is a matrix, the rows of which indicate the possible input/output partitionings, e.g.,

ud = wd(:, IO(1, 1:m)); yd = wd(:, IO(:, m+1:end));

is the first input/output partitioning.

Exercise: In the example above, all partitionings of the variables are valid input/output partitionings. Create an example where this is not the case.

In the examples above, we’ve used the function randB for obtaining a random data-generating system B and the function B2w for obtaining a random trajectory of that system.

• randB constructs a kernel representation R corresponding to the specified structure (number of variables, number of inputs, lag, and order) with nonzero elements generated as random uniformly distributed in the interval \([0,1]\). The kernel representation R is then converted to an input/state/output representation B.
• B2w uses the function lsim from the Control System Toolbox of Matlab in order to simulate a trajectory of B under random initial conditions and a random input.

The system produced by randB may be BIBO unstable:

B = R2ss(Rp, q); any(abs(eig(B)) > 1) % -> true

When it is unstable, the trajectory obtained from random input and initial conditions is almost surely exponentially diverging

T = 300; 
u = rand(T, m); y = lsim(B, u, [], rand(n, 1)); 
w = [u y]; max(w(:)) % -> large

Generating a “sufficiently long well behaved” random trajectory from unstable system requires special care (e.g., design of a stabilizing controller).

An alternative way to generate a random trajectory of the system is using the nonparameteric representation \(B_T\) of the finite-horizon behavior. For a random vector \(g\in\mathbb{R}^{mT+n}\), we obtain a random trajectory \(w = B_T g\). When \(||g||\) is bounded (e.g., its elements are uniformly distributed in the interval \([0,1]\)), since by construction the matrix \(B_T\) is orthonormal, the resulting trajectory \(w\) is also bounded.

BT = R2BT(Rp, q, T); w = BT * rand(size(BT, 2), 1); max(w) % -> small

Thus, using the nonparameteric representation \(B_T\) of the finite-horizon behavior, standard linear algebra (computation of an orthonormal basis) allows us to obtain “sufficiently long well behaved” trajectories of unstable systems.

The distance of a signal \(w\) to a system \(\mathcal{B}\) has the familiar definition and geometrical interpretation: \(\min_{\hat w \in \mathcal{B}} \| w - \hat w \|\). It’s computation is implemented in the function dist. As a specification of \(\mathcal{B}\), dist accepts an lti object as well as the finite-horizon representation \(B_T\) with horizon \(T\) the length of \(w\).

Finding the distance of a signal to a system is useful in particular for checking if the signal is a valid trajectory of the system:

dist(wd, BT) % -> 0

The function dist returns also the projection \(\hat w\) of \(w\) in \(\mathcal{B}\). The signal \(\hat w\) is useful because it the best approximation of \(w\) by a trajectory of the system. It is called the smoothed signal. Here is an example.

• We add a perturbation on the trajectory

  wn = randn(size(wd)); wd_noisy = wd + 0.1 * norm(wd) * wn / norm(wd); 
  

• and compute the projection wdh of the perturbed signal wd_noisy on the system

  [d, wdh] = dist(wd_noisy, BT);
  

• The distance of the best approximation to the original trajectory is smaller than the distance of the perturbed trajectory to the original trajectory:

  100 * [norm(wd - wd_noisy) norm(wd - wdh)] / norm(wd)
  

How to check if two systems are equal? By construction the system

Bp = ss2ss(B, rand(n));

is equal to B. However, the Control Systems Toolbox of Matlab does not allow us to infer it by

try 
    B == Bp % -> Operator '==' is not supported for operands of type 'ss'.
catch me
    me
end    

Instead, we can check the norm of difference:

norm(Bp - B) % -> Inf

however, it is restricted to stable systems.

The function Bdist computes a distance between systems defined as the angle between the corresponding behaviors:

Bdist(Bp, B) % -> 0

This distance measure is applicable also to unstable systems.