Linear Exact Modelling

A Tool to Explore Interrelations Between Data and Dynamic Systems

Martin Rickli
Automatic Control Lab, ETH Zurich, Switzerland
rickli@aut.ee.ethz.ch

Introduction

Given vector-valued data
  w(t) = (p0 + p1 t + p2 t^2 + ... + pk t^k ) exp(lambda t)
find the operator Theta(d/dt), such that
  Theta(d/dt) w(t) = 0

lambda may be a complex number but Theta has only real valued coefficients.

The algorithm is recursive in processing the number of data and proceeds in several steps

 Theta_new = Psort x Delta x Gamma x P x Theta_old
 
 It keeps track of controllable rows (rows of Theta that have full row rank for all lambda).

Content

Startup

Depending on your installation the method to load a package may be different.

<<ExactModeling.m	(* loads all necessary functions *)

(* late braking news from the author ... *)

ExactModelingInfo[]

Symbols (Functions) defined in ExactModeling

?ExactModeling`*

?ResetExactModeling

?ProcessNewMeasurement

Options

Options[ExactModeling]

?Gauss

?Orthogonal

?MaxDifference

?NoiseLevel

?Normalize

?Balance

?TriangularizeModel

Specialities

?SetMode

?StoreResults

?GetUpdateMatrices

?SavePrefs

Examples

The following examples should give you an impression what this package may be used for. We only cover basic features. Stay tuned for a more demanding application.

Pure Exponential Data (step by step)

p1 = {1,1};
p2 = {1,1/2};
p3 = {1,1/3};

(* At the beginning all rows are controllable:
   Length[ControlSet] = Dimension[p] *)

ResetExactModeling[s,{1,2},PrintLevel -> 9,
  Method -> Gauss, WorkingPrecision -> Infinity,
  NoiseLevel -> 0]

th = IdentityMatrix[2];

(* first data w1 = p1 Exp[0 t] *)
error = CalculateError[th,p1,0]

{P,G,Ps,th1} = 
  ProcessNewMeasurement[th,p1,0,s];

(* check result *)

error = CalculateError[th1,p1,0]

(* second data w2 = p2 Exp[1 t] *)
{P,G,Ps,th2} = 
  ProcessNewMeasurement[th1,p2,1,s];

(* third data w3 = p3 Exp[2 t] *)
{P,G,Ps,th3} = 
  ProcessNewMeasurement[th2,p3,2,s];

TableForm[th3]

CalculateTransferMatrix[th3]

Polynomial time exponential data

Let us compute an example where the data is a true polynomial-time exponential time series, i.e. p(t) is a true polynomial.
The update scheme can only annihilate pure exponential errors vectors. We therefore have to compute "related time series" w'.

p[1,0] = {0,0,1};
p[1,1] = {0,1,2};
p[1,2] = {1,2,3};
p1 = p[1,0] + p[1,1] t + (p[1,2]/2) t^2;

lam = 2;
th0 = IdentityMatrix[3];
ResetExactModeling[s,{1,2,3}, PrintLevel -> 3,
WorkingPrecision -> Infinity, NoiseLevel -> 0]

CheckMeasurement[th0,p1,lam]

(* -> build "new" (related) data *)

p11 = RelatedTimeSeries[p1,lam,t]
CheckMeasurement[th0,p11,lam]

p12 = RelatedTimeSeries[p11,lam,t]
CheckMeasurement[th0,p12,lam]

{pm,gam,ps,th01} =
  ProcessNewMeasurement[th0,p12,lam,s, Method -> Gauss,
    PrintLevel -> 3];

CheckMeasurement[th01,p1,lam]
CheckMeasurement[th01,p11,lam]

{pm,gam,ps,th02} =
  ProcessNewMeasurement[th01,p11,lam,s];

{pm,gam,ps,th03} =
  ProcessNewMeasurement[th02,p1,lam,s];

gs3 = CalculateTransferMatrix[th03]

Polynomial time exponential data (automatic)

Manually computing related time series is too tedious so we've built it in.

p[1,0] = {0,0,1};
p[1,1] = {0,1,2};
p[1,2] = {1,2,3};
p1 = p[1,0] + p[1,1] t + (p[1,2]/2) t^2;
lam = 2;
th0 = IdentityMatrix[3];
ResetExactModeling[s,{1,2,3}, PrintLevel->2,
  WorkingPrecision -> Infinity, NoiseLevel -> 0]

{pm,gam,ps,thpoly} =
  ProcessNewMeasurement[th0,p1,lam,s]

gsp = CalculateTransferMatrix[thpoly]

Transfer Matrix Data  Z(lambda)

The problem of rational interpolation, i.e. given frequency response data find a suitable dynamical system that reproduces the data, can be solved as well. The data vector w(t) has the form
w(t) = [1  Z(lam)] ' exp(lam t)
The function RatInt solves this problem. You don't have to build the data vectors w yourself. Input parameters are the list of frequencies (lam) and corresponding measurements (zlam).

(* Example: Z(s) = 1/(s+1) , lambda = 0, 1, 2 *)

?RatInt

theta = RatInt[{1, 1/2, 1/3},{0,1,2},s,
            PrintLevel->2, WorkingPrecision -> Infinity,
            NoiseLevel -> 0, Method -> Gauss]

TableForm[theta]

(* Transfer matrix *)

gs = CalculateTransferMatrix[theta]

Transfer Matrix To Theta Model (exact)

Instead of having frequency response data it is also possible to specify a transfer matrix and ask for a corresponding behavioral representation. This may be seen as an academic problem, but it's very useful for testing purposes. Note the use of complex data. Silently we assume complex conjugate data occurs as well giving rise to a description with real coefficients only.

fun = 1/(s-1)-1/(s+2)+1/(s+10);
lam = {1/10 + 3*I, -1/10 + I, 1/10, 1/10 + 2*I};

th = TF2Theta[fun,lam,s,PrintLevel->2, Method -> Gauss,
      Sorting->None, WorkingPrecision->Infinity,
      NoiseLevel -> 0]

gs = CTM[th]
Factor@Together@fun

Transfer Matrix To Theta Model (finite precision)

Working with finite precision numbers may lead to unexpected results especially in the context of _exact_ modeling!

fun = 1/(s-1)-1/(s+2)+1/(s+10);
lam = {1/10,1/10+I,1/10+2*I,1/10+3*I};

thn = TF2Theta[fun,lam,s,PrintLevel->0, Sorting->None,
   WorkingPrecision->16, NoiseLevel -> 10.^-11]

gsn = CTM[thn]
Simplify@Together@fun

Transfer Matrix To Theta Model (data driven sorting)

One challenge in algorithm design is handling very small and very big numbers. In our framework, the order of a polynom is very crucial and depends on the decision which numbers are taken to be zero.
Given a set of data which are to be processed recursively the order of calculation may highly influence the result. We suggest to apply a test on each error vector and proceed with the data wich gives the biggest error. This test may be activated with option Sorting -> MaxError.

fun = 1/(s-1)-1/(s+2)+1/(s+10);
lam = {.1,.1+I,.1+2*I,.1+3*I};
lam = {1/10,1/10+I,1/10+2*I,1/10+3*I};

thn2 = TF2Theta[fun,lam,s,PrintLevel->0, Sorting->MaxError,
   WorkingPrecision->16, NoiseLevel -> 10.^-7]

gsn2 = CTM[thn2]
Simplify@Together@fun

MIMO Impulse Response Data (automatic)

The realization problem may be reformulated as a rational interpolation problem where impulse response data corresponds to frequency response data (and derivations thereof) at frequency zero.
The function MarkovModel was implemented to tackle this kind of problem.

(* Markov parameters of a 2x3 system *)

A4 = {{8,11,6},{16,21,7}};
A3 = {{4,6,3},{8,10,2}};
A2 = {{2,2,1},{4,4,0}};
A1 = {{1,1,1},{2,2,0}};
A0 = 0 A1;

th3=MarkovModel[A0,A1,A2,A3,q,PrintLevel->3,
     Method->Gauss, WorkingPrecision->Infinity,
     NoiseLevel -> 0, Normalize->NormHRC]

g3 = CalculateTransferMatrix[th3]