Local Model Structures
local_model_structures.RdParametrization for "local" model classes, in particular, "Data Driven Local Coordinates" as detailed in (McKelvey et al. 2004) and (Ribarits et al. 2005) .
Arguments
- model
stspmod()object, which represents a state space model. Only the case \(m = n > 0\) is implemented, i.e. the output process and the noise process must be of the same dimension.- balance
(character string) For
balance = "lyapunov"orbalance = "minimum phase"the reference model is first balanced by the respective scheme.- sigma_L
(character string) determines the form of the (left) square root of the noise covariance \(\Sigma\). The choice
"chol"gives a lower triangular matrix,"symm"gives a symmetric matrix and"identity"corresponds to the (fixed) identity matrix.
Value
Model template, i.e. a list with slots
h\(((m+s)^2 + m^2)\)-dimensional vector,
H\(((m+s)^2 + m^2, k)\)-dimensional matrix,
class = "stspmod"order = c(m,m,s)and
n.parnumber of free parameters \(=k\).
See also model structures() for more details.
Details
The function tmpl_DDLC and tmpl_GRAM construct model templates which describe
models in a neighborhood of a given reference model.
In a first step the reference state space model is transformed to \(D=I\) and eventually
(depending on the parameter "balance") balanced.
state space models are described by a quadruple \((A,B,C,D=I)\) of matrices which may be embedded into an \((s^2+2ms)\)-dimensional euclidean space. Note that the parameter matrices are not uniqely determined from the ACF or the spectral density of the process, i.e. there is an inherent non identifiablity problem. For minimal models the "equivalence class" of models, which represent the same ACF is given by the set of all models which may be obtained by a state transformation \((A,B,C,D) \rightarrow (TAT^{-1}, TB, CT^{-1}, D)\).
The DDLC parametrization now considers models, \((A,B,C,D=I)\), which are contained in the \(2ms\)-dimensional subspace, which is orthogonal to the \(s^2\)-dimensional tangent space of the set of equivalent models.
The routine tmpl_GRAM considers the \(2ms\)-dimensional subspace, where
models close to the reference models are "approximately" balanced.
Both schemes may fail for "non-generic" models. tmpl_DDLC issues a
warning message and tmpl_GRAM throws an error, in cases where the
\(2ms\)-dimensional subspace is not well defined.
Note that also the parametrization of the left square root L=sigma_L of the
noise covariance is "local", i.e. th = 0 corresponds to the (balanced) reference model.
References
McKelvey T, Helmersson A, Ribarits T (2004). “Data driven local coordinates for multivariable linear systems and their application to system identification.” Automatica, 40(9), 1629 - 1635. doi:10.1016/j.automatica.2004.04.015 .
Ribarits T, Deistler M, Hanzon B (2005). “An analysis of separable least squares data driven local coordinates for maximum likelihood estimation of linear systems.” Automatica, 41(3), 531 - 544. doi:10.1016/j.automatica.2004.11.014 . .
See also
For the computation of Grammians and for balancing of state space models,
see rationalmatrices::balance().
Examples
# create a random state space model with m outputs and s states
m = 3
s = 6
tmpl = tmpl_stsp_full(m, n = m, s, sigma_L = 'symm')
model = r_model(tmpl, bpoles = 1.1, bzeroes = 1.1, sd = 1/s)
model # note that sigma_L is symmetric
#> state space model [3,3] with s = 6 states
#> s[1] s[2] s[3] s[4] s[5] s[6]
#> s[1] -0.393508577 0.12461571 -0.06924699 -0.015220683 -0.07161459 0.04982169
#> s[2] -0.122487510 0.06695231 -0.10280646 0.059401626 -0.04996684 0.01680407
#> s[3] 0.123744942 -0.13837876 -0.24297760 0.153227762 -0.19395870 0.16455487
#> s[4] 0.023668034 -0.06907320 0.10928726 0.308424292 0.01952750 0.32330508
#> s[5] 0.308284473 0.17888572 -0.04409380 -0.076322389 -0.15290357 -0.23584004
#> s[6] -0.163105350 0.13632526 -0.03909642 -0.009370285 0.09134868 0.29838002
#> x[1] 0.071956757 -0.01578565 -0.19892960 0.118843089 -0.09640141 -0.02195228
#> x[2] -0.008326129 0.16998545 0.09406509 -0.150617391 -0.19797936 0.07510532
#> x[3] -0.048403315 -0.10458266 0.06491167 -0.057763783 0.16036459 0.04696967
#> u[1] u[2] u[3]
#> s[1] 0.009792037 -0.18341927 0.28403395
#> s[2] 0.014661062 -0.15528459 -0.07679769
#> s[3] 0.089428255 0.04409296 -0.04988222
#> s[4] 0.176704127 0.04394464 0.09713260
#> s[5] -0.385635619 -0.05439723 -0.05957528
#> s[6] -0.106408031 -0.15011094 -0.05944096
#> x[1] 1.000000000 0.00000000 0.00000000
#> x[2] 0.000000000 1.00000000 0.00000000
#> x[3] 0.000000000 0.00000000 1.00000000
#> Left square root of noise covariance Sigma:
#> u[1] u[2] u[3]
#> u[1] -0.07063615 0.04163931 -0.03794890
#> u[2] 0.04163931 0.01316882 0.33093668
#> u[3] -0.03794890 0.33093668 -0.07177565
model$sigma_L %*% t(model$sigma_L) # noise covariance Sigma
#> [,1] [,2] [,3]
#> [1,] 0.008163417 -0.01495158 0.01918435
#> [2,] -0.014951582 0.11142634 -0.02097532
#> [3,] 0.019184346 -0.02097532 0.11611095
# tmpl_DDLC #############################################
# create a DDLC parametrization of a neighborhood of this model
tmpl = tmpl_DDLC(model, balance = 'lyapunov', sigma_L = 'chol')
# for th = 0, we get the original model (in balanced form)
model = fill_template(numeric(tmpl$n.par), tmpl)
model # note that sigma_L is lower triangular
#> state space model [3,3] with s = 6 states
#> s[1] s[2] s[3] s[4] s[5]
#> s[1] -0.377930599 0.31757105 -0.3241297735 -0.07404308 0.056286812
#> s[2] -0.034259926 -0.01361698 0.0526850557 -0.32231861 -0.009521599
#> s[3] -0.009362975 0.15498744 -0.1722413971 -0.28649675 -0.302613258
#> s[4] 0.334284819 -0.05915110 -0.1919520286 -0.10945703 -0.004483441
#> s[5] 0.023765087 0.01874922 0.3206584043 -0.09936860 0.526304887
#> s[6] 0.004953518 -0.01447720 0.0007736472 0.22019872 0.112936038
#> x[1] -0.199863126 0.09536832 0.0311445852 0.06869977 -0.018603486
#> x[2] -0.130450885 -0.20662250 -0.0659325820 0.02272639 -0.004311740
#> x[3] 0.229541793 0.03544278 -0.0745825899 0.05630364 -0.007114404
#> s[6] u[1] u[2] u[3]
#> s[1] 0.002187688 -0.272201384 0.07570221 -0.149386203
#> s[2] 0.018085937 -0.106723949 0.04504618 0.226561059
#> s[3] 0.044362348 0.060908643 0.13648581 -0.004574977
#> s[4] -0.217370165 -0.043712186 0.01564708 -0.019029603
#> s[5] 0.154031866 -0.003987267 0.01248427 -0.008349216
#> s[6] 0.031307991 -0.001417133 0.01162225 0.002898405
#> x[1] -0.001377618 1.000000000 0.00000000 0.000000000
#> x[2] 0.003036955 0.000000000 1.00000000 0.000000000
#> x[3] 0.007036271 0.000000000 0.00000000 1.000000000
#> Left square root of noise covariance Sigma:
#> u[1] u[2] u[3]
#> u[1] 0.09035163 0.00000000 0.0000000
#> u[2] -0.16548215 0.28989997 0.0000000
#> u[3] 0.21232984 0.04884955 0.2619937
model$sigma_L %*% t(model$sigma_L) # however Sigma is the same as above
#> [,1] [,2] [,3]
#> [1,] 0.008163417 -0.01495158 0.01918435
#> [2,] -0.014951582 0.11142634 -0.02097532
#> [3,] 0.019184346 -0.02097532 0.11611095
#' apply a "small" state transformation T = (diag(s)+eps*X)
eps = sqrt(.Machine$double.eps)
sys = model$sys
d_sys = state_trafo(sys, diag(s) + matrix(rnorm(s^2, sd = eps), nrow = s, ncol = s))
d_pi = (as.vector(unclass(d_sys) - unclass(sys)))/eps
# The vector d_pi is (close to) an element of the tangent space
# of the set of models, which are generated by a state transformation
# of the reference model
# by construction d_pi is (close to) orthogonal to tmpl$H
max(abs(d_pi %*% tmpl$H[1:((m+s)^2), , drop = FALSE]))
#> [1] 1.892004e-08
# the tmpl_DDLC routine may fail in some exceptional cases
m = 1
s = 3
model = stspmod(sys = stsp(A = matrix(0, nrow = s, ncol = s),
B = matrix(rnorm(m*s), nrow = s, ncol = m),
C = matrix(rnorm(m*s), nrow = m, ncol = s),
D = diag(m)),
sigma_L = diag(m))
# For this model "tmpl_DLLC" issues a warning.
junk = tmpl_DDLC(model, sigma_L = 'chol', balance = 'none')
#> Warning: The tangent space of the equivalence class does not have dimension s^2=9 (sv[1]=2.05170976570878, sv[9]=3.2202148341986e-18)
# tmpl_GRAM #############################################
model = fill_template(numeric(tmpl$n.par), tmpl)
tmpl = tmpl_GRAM(model, sigma_L = 'chol')
model = fill_template(numeric(tmpl$n.par), tmpl)
# check grammians
gr = grammians(model$sys, 'lyapunov')
P = gr$P
Q = gr$Q
# P=Q=diag() should hold!
print(round(cbind(P, P-Q), 6))
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
#> [1,] 0.130444 0.000000 0.000000 0.000000 0.000000 0.000000 0 0 0 0
#> [2,] 0.000000 0.066921 0.000000 0.000000 0.000000 0.000000 0 0 0 0
#> [3,] 0.000000 0.000000 0.026716 0.000000 0.000000 0.000000 0 0 0 0
#> [4,] 0.000000 0.000000 0.000000 0.018589 0.000000 0.000000 0 0 0 0
#> [5,] 0.000000 0.000000 0.000000 0.000000 0.004559 0.000000 0 0 0 0
#> [6,] 0.000000 0.000000 0.000000 0.000000 0.000000 0.001123 0 0 0 0
#> [,11] [,12]
#> [1,] 0 0
#> [2,] 0 0
#> [3,] 0 0
#> [4,] 0 0
#> [5,] 0 0
#> [6,] 0 0
# now consider a model close to the reference model
d_th = rnorm(tmpl$n.par, sd = eps)
d_model = fill_template(d_th, tmpl)
d_sys = d_model$sys
gr = grammians(d_sys, 'lyapunov')
d_P = gr$P - P
d_Q = gr$Q - Q
# the "disturbed" system should still be approximately balanced!
print(round(cbind(d_P, d_P - d_Q)/eps, 6) )
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 0.464405 0.000000 0.000000 0.000000 0.000000 0.000000 0 0 0
#> [2,] 0.000000 -0.177384 0.000000 0.000000 0.000000 0.000000 0 0 0
#> [3,] 0.000000 0.000000 -0.080407 0.000000 0.000000 0.000000 0 0 0
#> [4,] 0.000000 0.000000 0.000000 -0.103301 0.000000 0.000000 0 0 0
#> [5,] 0.000000 0.000000 0.000000 0.000000 -0.030998 0.000000 0 0 0
#> [6,] 0.000000 0.000000 0.000000 0.000000 0.000000 0.002177 0 0 0
#> [,10] [,11] [,12]
#> [1,] 0 0 0
#> [2,] 0 0 0
#> [3,] 0 0 0
#> [4,] 0 0 0
#> [5,] 0 0 0
#> [6,] 0 0 0