# SPDX-License-Identifier: GPL-3.0-or-later
# Copyright (C) 2015-2025 Simon Marwitz, Volkmar Zabel, Andrei Udrea et al.
"""Data-driven SSI (SSIData, SSIDataMC) for operational modal analysis."""
import numpy as np
import scipy.linalg
import os
from .Helpers import lq_decomp, validate_array, simplePbar, ConfigFile
from .PreProcessingTools import PreProcessSignals
from .ModalBase import ModalBase
import logging
logger = logging.getLogger(__name__)
logger.setLevel(level=logging.INFO)
[docs]
class SSIDataMC(ModalBase):
"""Reference-based Data-driven Stochastic Subspace Identification (SSI-Data/MC).
Identifies modal parameters from raw time-series data by constructing a
block-Hankel matrix, computing its LQ decomposition, and extracting state-
space models via SVD. The standard workflow is:
1. :meth:`build_block_hankel` — build the projection matrix from raw data.
2. :meth:`compute_modal_params` — run the multi-order modal identification.
3. Pass the result to :class:`~pyOMA.core.StabilDiagram.StabilCalc` for
stabilisation-diagram analysis.
Parameters
----------
prep_signals : PreProcessSignals
Pre-processed signal object providing raw time-series and channel
metadata.
.. TODO::
* define unit tests to check functionality after changes
* add switch to keep synthesized time-histories
"""
def __init__(self, *args, **kwargs):
"""
Parameters
----------
*args, **kwargs
Passed to :class:`~pyOMA.core.ModalBase.ModalBase`.
"""
super().__init__(*args, **kwargs)
self.state = [False, False, False, False]
self.num_block_rows = None
self.num_blocks = 1
self.P_i_ref = None
self.P_i_1 = None
self.Y_i_i = None
self.S = None
self.U = None
self.V_T = None
self.max_model_order = None
self.modal_contributions = None
[docs]
@classmethod
def init_from_config(cls, conf_file, prep_signals):
cfg = ConfigFile(conf_file)
num_block_rows = cfg.int('Number of Block-Columns')
max_model_order = cfg.int('Maximum Model Order')
ssi_object = cls(prep_signals)
ssi_object.build_block_hankel(num_block_rows)
ssi_object.compute_modal_params(max_model_order)
return ssi_object
[docs]
def write_config(self, conf_file):
ConfigFile.write(conf_file, {
'Number of Block-Columns': self.num_block_rows,
'Maximum Model Order': self.max_model_order,
})
[docs]
def build_block_hankel(self, num_block_rows=None, reduced_projection=True):
'''
Builds a Block-Hankel Matrix of the measured time series with varying
time lags and estimates the subspace matrix from its LQ decomposition.
::
<- num_time samples - num_block_rows-> _
[ y_0 y_1 ... y_(j-1) ]^
[ y_1 y_2 ... y_j ]num_block_rows (=i)*n_l
[ ... ... ... ... ]v
[ y_(2i-1) y_(2i) ... y_(2i+j-2) ]_
The notation mostly follows Peeters 1999.
Parameters
-------
num_block_rows: integer, required
The number of block rows of the Subspace matrix
'''
if num_block_rows is None:
num_block_rows = self.num_block_rows
if not isinstance(num_block_rows, int):
raise TypeError(f"num_block_rows must be an int, got {type(num_block_rows)}")
self.num_block_rows = num_block_rows
signals = self.prep_signals.signals
total_time_steps = self.prep_signals.total_time_steps
ref_channels = sorted(self.prep_signals.ref_channels)
n_l = self.prep_signals.num_analised_channels
n_r = self.prep_signals.num_ref_channels
q = num_block_rows
p = num_block_rows
N = int(total_time_steps - 2 * p)
logger.info(f'Building Block-Hankel matrix with {p} block-columns and {q} block rows')
Hankel_matrix = self._assemble_hankel_matrix(signals, ref_channels, n_l, n_r, p, q, N)
logger.debug(Hankel_matrix.shape)
l, q = lq_decomp(Hankel_matrix, mode='full')
logger.info('Estimating subspace matrix...')
L_red, Q_red, P_i_ref, P_i_1, Y_i_i, U, S, V_T = \
self._extract_subspace_matrices(l, q, Hankel_matrix, n_l, n_r, p, reduced_projection)
self.L_red = L_red
self.Q_red = Q_red
self.P_i_1 = P_i_1
self.P_i_ref = P_i_ref
self.Y_i_i = Y_i_i
self.S = S
self.U = U
self.V_T = V_T
self.state[0] = True
@staticmethod
def _assemble_hankel_matrix(signals, ref_channels, n_l, n_r, p, q, N):
"""Fill and return the normalised block-Hankel data matrix.
Parameters
----------
signals : numpy.ndarray, shape (total_time_steps, n_l)
ref_channels : list of int
n_l, n_r : int
p, q : int
Number of block rows for future/past halves (both equal to num_block_rows).
N : int
Number of columns (time windows).
Returns
-------
Hankel_matrix : numpy.ndarray, shape ((q*n_r + (p+1)*n_l), N)
"""
Y_minus = np.zeros((q * n_r, N))
Y_plus = np.zeros(((p + 1) * n_l, N))
for ii in range(q):
Y_minus[(q - ii - 1) * n_r:(q - ii) * n_r, :] = signals[(ii):(ii + N), ref_channels].T
for ii in range(p + 1):
Y_plus[ii * n_l:(ii + 1) * n_l, :] = signals[(q + ii):(q + ii + N)].T
Hankel_matrix = np.vstack((Y_minus, Y_plus))
Hankel_matrix /= np.sqrt(N)
return Hankel_matrix
@staticmethod
def _extract_subspace_matrices(l, q, Hankel_matrix, n_l, n_r, p, reduced_projection):
"""Extract the projection / subspace matrices from the LQ factors.
Returns
-------
L_red, Q_red, P_i_ref, P_i_1, Y_i_i : numpy.ndarray
U, S, V_T : numpy.ndarray
SVD factors of either the reduced L block (reduced_projection=True)
or P_i_ref (reduced_projection=False). V_T is None when
reduced_projection=True.
"""
a = n_r * p
b = n_r
c = n_l - n_r
d = n_l * p
P_i_ref = l[a:a + b + c + d, :a] @ q[:a, :]
if reduced_projection:
U, S, V_T = np.linalg.svd(l[a:a + b + c + d, :a], full_matrices=False)
V_T = None
else:
U, S, V_T = np.linalg.svd(P_i_ref, full_matrices=False)
P_i_1 = l[a + b + c:a + b + c + d, :a + b] @ q[:a + b, :]
# Y_i_i = l[a : a + b + c, : a + b + c] @ q[ : a + b + c, : ]
Y_i_i = Hankel_matrix[a:a + b + c, :]
L_red = l[a:, :a + b]
Q_red = q[:a + b, :]
return L_red, Q_red, P_i_ref, P_i_1, Y_i_i, U, S, V_T
[docs]
def compute_modal_params(self, max_model_order=None,
j=None, validation_blocks=None,
synth_sig=True):
'''
Perform a multi-order computation of modal parameters. Successively
calls
* estimate_state(order,)
* modal_analysis(A,C)
* synthesize_signals(A, C, Q, R, S, j)
at ascending model orders, up to max_model_order.
See the explanations in the the respective methods, for a detailed
explanation of parameters.
Parameters
----------
max_model_order: integer, optional
Maximum model order, where to interrupt the algorithm. If not given,
it is determined from the previously computed subspace matrix.
'''
if not self.state[0]:
raise RuntimeError("Call build_block_hankel() first.")
if max_model_order is not None:
if not isinstance(max_model_order, int):
raise TypeError(f"Expected int for 'max_model_order', got {type(max_model_order).__name__!r}.")
else:
max_model_order = self.S.shape[0]
if max_model_order > self.S.shape[0]:
raise ValueError(f"max_model_order must be <= {self.S.shape[0]}, got {max_model_order}.")
# num_block_rows = self.num_block_rows
num_analised_channels = self.prep_signals.num_analised_channels
# num_ref_channels = self.prep_signals.num_ref_channels
# sampling_rate = self.prep_signals.sampling_rate
if j is None and validation_blocks is None:
j = self.prep_signals.total_time_steps
elif validation_blocks is None:
if j > self.prep_signals.signals.shape[0]:
raise ValueError(f"j must be <= {self.prep_signals.signals.shape[0]}, got {j}.")
logger.info('Computing modal parameters...')
eigenvalues = np.zeros(
(max_model_order, max_model_order), dtype=np.complex128)
modal_frequencies = np.zeros((max_model_order, max_model_order))
modal_damping = np.zeros((max_model_order, max_model_order))
mode_shapes = np.zeros(
(num_analised_channels,
max_model_order,
max_model_order),
dtype=complex)
if synth_sig:
modal_contributions = np.zeros((max_model_order, max_model_order))
else:
modal_contributions = None
pbar = simplePbar(max_model_order - 1)
for order in range(1, max_model_order):
next(pbar)
A, C, Q, R, S = self.estimate_state(order)
f, d, phi, lamda, = self.modal_analysis(A, C)
modal_frequencies[order,:order] = f
modal_damping[order,:order] = d
mode_shapes[:phi.shape[0],:order, order] = phi
eigenvalues[order,:order] = lamda
if synth_sig:
_, delta = self.synthesize_signals(A, C, Q, R, S, j=j, validation_blocks=validation_blocks)
modal_contributions[order,:order] = delta
self.max_model_order = max_model_order
self.modal_contributions = modal_contributions
self.modal_frequencies = modal_frequencies
self.modal_damping = modal_damping
self.mode_shapes = mode_shapes
self.eigenvalues = eigenvalues
self.state[2] = True
[docs]
def estimate_state(self, order,):
'''
Estimate the state matrices A, C and noise covariances Q, R and S from
the subspace / projection matrix. Several methods exist, e.g.
* Peeters 1999 Reference Based Stochastic Subspace Identification for OMA
* DeCock 2007 Subspace Identification Methods
* the algorithm used in BRSSICovRef.
Here, the first algorithm, a residual-based computation of Q, R and S,
is implemented.
Parameters
----------
order: integer, required
The model order, at which to truncate the singular values of the
projection Matrix P_i_ref
Returns
-------
A: numpy.ndarray
State matrix: Array of shape (order, order)
C: numpy.ndarray
Output matrix: Array of shape (num_analised_channels, order)
Q: numpy.ndarray
state noise covariance matrix: Symmetric array of shape (order, order)
R: numpy.ndarray
signal noise covariance matrix: Array of shape (num_analised_channels, num_analised_channels)
S: numpy.ndarray
system noise - signal noise covariance matrix: Array of shape (order, num_analised_channels)
'''
num_block_rows = self.num_block_rows
num_analised_channels = self.prep_signals.num_analised_channels
# num_ref_channels = self.prep_signals.num_ref_channels
U = self.U[:,:order]
S = self.S[:order]
P_i_1 = self.P_i_1 # n_l * p x n_r * (p + 1) x N
Y_i_i = self.Y_i_i # n_l x n_r *p + n_l x N
# compute state-space model
S_2 = np.power(S, 0.5)
O = U * S_2[np.newaxis,:]
O_i_1 = O[:num_analised_channels * num_block_rows,:order]
O_i = O[:,:order]
if self.V_T is not None:
V_T = self.V_T[:order,:]
X_i = S_2[:, np.newaxis] * V_T
else:
P_i_ref = self.P_i_ref # (p + 1) * n_l x n_r * p x N
X_i = np.linalg.pinv(O_i) @ P_i_ref
X_i_1 = np.linalg.pinv(O_i_1) @ P_i_1
X_i_1_Y_i = np.vstack((X_i_1, Y_i_i))
AC = X_i_1_Y_i @ np.linalg.pinv(X_i)
A = AC[:order,:]
C = AC[order:,:]
roh_w_v = X_i_1_Y_i - AC @ X_i
QSR = roh_w_v @ roh_w_v.T
Q = QSR[:order,:order]
S = QSR[:order, order:order + num_analised_channels]
R = QSR[order:order + num_analised_channels,
order:order + num_analised_channels]
return A, C, Q, R, S
[docs]
def modal_analysis(self, A, C, rescale_fun=None):
'''
Computes the modal parameters from a given state space model as described
by Peeters 1999 and Döhler 2012. Mode shapes are scaled to unit modal
displacements. Complex conjugate and real modes are removed prior to
further processing. Typically, order // 2 modes are in the returned arrays.
Parameters
----------
A: numpy.ndarray
State matrix: Array of shape (order, order)
C: numpy.ndarray
Output matrix: Array of shape (num_analised_channels, order)
Returns
-------
modal_frequencies: numpy.ndarray, shape (order,)
Array holding the modal frequencies for each mode
modal_damping: numpy.ndarray, shape (order,)
Array holding the modal damping ratios (0,100) for each mode
mode_shapes: numpy.ndarray, shape (n_l, order,)
Complex array holding the mode shapes
eigenvalues: numpy.ndarray, shape (order,)
Complex array holding the eigenvalues for each mode
'''
# collect variables
accel_channels = self.prep_signals.accel_channels
velo_channels = self.prep_signals.velo_channels
sampling_rate = self.prep_signals.sampling_rate
n_l = self.num_analised_channels
order = A.shape[0]
if order != A.shape[1]:
raise RuntimeError(f"Internal error: A must be square, got shape {A.shape}.")
# allocate output arrays
modal_frequencies = np.full((order), np.nan)
modal_damping = np.full((order), np.nan)
mode_shapes = np.full((n_l, order), np.nan, dtype=complex)
eigenvalues = np.full((order), np.nan, dtype=complex)
# compute modal model
eigvals, eigvecs_r = np.linalg.eig(A)
Phi = C.dot(eigvecs_r)
conj_indices = self.remove_conjugates(eigvals, eigvecs_r, inds_only=True)
for i, ind in enumerate(conj_indices):
lambda_i = eigvals[ind]
mode_shape_i = Phi[:, ind]
a_i = np.abs(np.arctan2(np.imag(lambda_i), np.real(lambda_i)))
b_i = np.log(np.abs(lambda_i))
freq_i = np.sqrt(a_i ** 2 + b_i ** 2) * sampling_rate / 2 / np.pi
damping_i = 100 * np.abs(b_i) / np.sqrt(a_i ** 2 + b_i ** 2)
if rescale_fun is not None:
mode_shape_i = rescale_fun(mode_shape_i)
# scale modeshapes to modal displacements
mode_shape_i = self.integrate_quantities(
mode_shape_i, accel_channels, velo_channels, freq_i * 2 * np.pi)
# rotate mode shape in complex plane
mode_shape_i = self.rescale_mode_shape(mode_shape_i)
modal_frequencies[i] = freq_i
modal_damping[i] = damping_i
mode_shapes[:mode_shape_i.shape[0], i] = mode_shape_i
eigenvalues[i] = lambda_i
argsort = np.argsort(modal_frequencies)
return modal_frequencies[argsort], modal_damping[argsort], mode_shapes[:, argsort], eigenvalues[argsort],
@staticmethod
def _solve_kalman_gain(CPCR, APCS):
"""Solve for the steady-state Kalman gain from the DARE solution.
Raises a more informative LinAlgError than the bare numpy one when
the innovation covariance CPCR is singular: this typically happens
with noise-free/synthetic data, where the estimated noise
covariances (Q, R, S) are degenerate; real measurement data with
genuine ambient/process noise usually avoids it.
"""
try:
return np.linalg.solve(CPCR.T, APCS.T).T
except np.linalg.LinAlgError as exc:
raise np.linalg.LinAlgError(
f"{exc}: singular innovation covariance while solving for "
"the Kalman gain. This typically happens with noise-free/"
"synthetic data, where the estimated noise covariances "
"(Q, R, S) are degenerate; real measurement data with "
"genuine ambient/process noise usually avoids it."
) from exc
[docs]
def synthesize_signals(self, A, C, Q, R, S, j=None, **kwargs):
'''
Computes the modal response signals and the contribution of each mode.
The algorithm follows Peeters 1999 and the Lyapunov equation is solved
as a discrete-time algebraic Riccati equation (DARE). For long signals,
the computation may become time-consuming, thus only time steps up to j
may be used to synthesize the signal.
Parameters
----------
A: numpy.ndarray
State matrix: Array of shape (order, order)
C: numpy.ndarray
Output matrix: Array of shape (num_analised_channels, order)
Q: numpy.ndarray
state noise covariance matrix: Symmetric array of shape (order, order)
R: numpy.ndarray
signal noise covariance matrix: Array of shape (num_analised_channels, num_analised_channels)
S: numpy.ndarray
system noise - signal noise covariance matrix: Array of shape (order, num_analised_channels)
j: integer, optional
length of signal to synthesize (number of timesteps)
Returns
-------
sig_synth: numpy.ndarray, shape (num_analised_channels, j, order // 2)
Array holding the modally decomposed input signals for
each channel n_l and all modes
modal_contributions: numpy.ndarray, shape (order, )
Array holding the contributions of each mode to the input
signals.
'''
order = A.shape[0]
if order != A.shape[1]:
raise RuntimeError(f"Internal error: A must be square, got shape {A.shape}.")
if j is None:
j = self.prep_signals.total_time_steps
signals = self.prep_signals.signals[:j,:]
n_l = self.prep_signals.num_analised_channels
modal_contributions = np.zeros((order))
sig_synth = np.zeros((n_l, j, order // 2))
try:
P = scipy.linalg.solve_discrete_are(
a=A.T, b=C.T, q=Q, r=R, s=S, balanced=True)
except Exception:
logger.warning('Correlations of residuals are not symmetric. Skiping Modal Contributions')
return sig_synth, modal_contributions
APCS = A @ P @ C.T + S
CPCR = C @ P @ C.T + R
K = self._solve_kalman_gain(CPCR, APCS)
eigvals, eigvecs_r = np.linalg.eig(A)
conj_indices = self.remove_conjugates(eigvals, eigvecs_r, inds_only=True)
A_0 = np.diag(eigvals)
C_0 = C @ eigvecs_r
K_0 = np.linalg.solve(eigvecs_r, K)
states = np.zeros((order, j), dtype=complex)
AKC = A_0 - K_0 @ C_0
K_0m = K_0 @ signals.T
for k in range(j - 1):
states[:, k + 1] = K_0m[:, k] + AKC @ states[:, k]
Y = signals.T
# norm = 1 / np.einsum('ji,ji->j', Y, Y)
Sigma_data = np.einsum('ji,ji->j', Y, Y)
Sigma_data_synth = np.zeros((n_l, order))
for i, ind in enumerate(conj_indices):
lambda_i = eigvals[ind]
ident = eigvals == lambda_i.conj()
ident[ind] = 1
C_0I = C_0[:, ident]
this_sig_synth = C_0I @ states[ident,:]
if not np.all(np.isclose(this_sig_synth.imag, 0)):
logger.warning(f'Synthetized signals are complex at mode index {order}:{ind}.')
sig_synth[:,:, i] = this_sig_synth.real
mYT = np.einsum('ji,ji->j', sig_synth[:,:, i], Y)
Sigma_data_synth[:, i] = mYT
# modal_contributions[i] = np.mean(norm * mYT)
total_sig_synth = np.sum(sig_synth, axis=-1) # shape n_l, N
Sigma_synth = np.einsum('ji,ji->j', total_sig_synth, total_sig_synth)
modal_contributions[:] = np.mean(Sigma_data_synth / np.sqrt(Sigma_data * Sigma_synth)[:, np.newaxis], axis=0)
self._sig_synth = sig_synth
self._modal_conributions = modal_contributions
return sig_synth, modal_contributions
[docs]
def save_state(self, fname):
logger.info('Saving results to {}...'.format(fname))
dirname, _ = os.path.split(fname)
if not os.path.isdir(dirname):
os.makedirs(dirname)
out_dict = {'self.state': self.state}
out_dict['self.setup_name'] = self.setup_name
out_dict['self.start_time'] = self.start_time
if self.state[0]: # subspace matrix
out_dict['self.num_block_rows'] = self.num_block_rows
out_dict['self.num_blocks'] = self.num_blocks
# out_dict['self.P_i_1'] = self.P_i_1
# out_dict['self.P_i_ref'] = self.P_i_ref
out_dict['self.L_red'] = self.L_red
out_dict['self.Q_red'] = self.Q_red
out_dict['self.Y_i_i'] = self.Y_i_i
out_dict['self.S'] = self.S
out_dict['self.U'] = self.U
out_dict['self.V_T'] = self.V_T
if self.state[2]: # modal params
out_dict['self.modal_frequencies'] = self.modal_frequencies
out_dict['self.modal_damping'] = self.modal_damping
out_dict['self.eigenvalues'] = self.eigenvalues
out_dict['self.mode_shapes'] = self.mode_shapes
out_dict['self.modal_contributions'] = self.modal_contributions
out_dict['self.max_model_order'] = self.max_model_order
np.savez_compressed(fname, **out_dict)
[docs]
@classmethod
def load_state(cls, fname, prep_signals):
logger.info('Loading results from {}'.format(fname))
in_dict = np.load(fname, allow_pickle=True)
if 'self.state' in in_dict:
# bool(...): entries loaded straight out of the .npz archive are
# numpy.bool_, not plain Python bool.
state = [bool(s) for s in in_dict['self.state']]
else:
return
if not isinstance(prep_signals, PreProcessSignals):
raise TypeError(f"Expected PreProcessSignals for 'prep_signals', got {type(prep_signals).__name__!r}.")
setup_name = str(in_dict['self.setup_name'].item())
if setup_name != prep_signals.setup_name:
raise ValueError(f"'setup_name' must be one of {[prep_signals.setup_name]}, got {setup_name!r}.")
start_time = prep_signals.start_time
if start_time != prep_signals.start_time:
raise ValueError(f"'start_time' must match prep_signals.start_time {prep_signals.start_time!r}, got {start_time!r}.")
# prep_signals = in_dict['self.prep_signals'].item()
ssi_object = cls(prep_signals)
if state[0]: # subspace matrix
ssi_object.num_block_rows = validate_array(in_dict['self.num_block_rows'])
ssi_object.num_blocks = validate_array(in_dict['self.num_blocks'])
ssi_object.L_red = validate_array(in_dict['self.L_red'])
ssi_object.Q_red = validate_array(in_dict['self.Q_red'])
# rebuild projections from reduced L and Q matrices
# saves ~ 95 % storage (Xe2 MB) compared to storing projections directly
a = prep_signals.num_ref_channels * ssi_object.num_block_rows
b = prep_signals.num_ref_channels
c = prep_signals.num_analised_channels - prep_signals.num_ref_channels
d = prep_signals.num_analised_channels * ssi_object.num_block_rows
ssi_object.P_i_ref = ssi_object.L_red[:b + c + d,: a] @ ssi_object.Q_red[:a,:]
ssi_object.P_i_1 = ssi_object.L_red[b + c:b + c + d,: ] @ ssi_object.Q_red[:,: ]
ssi_object.Y_i_i = validate_array(in_dict['self.Y_i_i'])
ssi_object.S = validate_array(in_dict['self.S'])
ssi_object.U = validate_array(in_dict['self.U'])
ssi_object.V_T = validate_array(in_dict['self.V_T'])
ssi_object.S = validate_array(in_dict['self.S'])
ssi_object.U = validate_array(in_dict['self.U'])
ssi_object.V_T = validate_array(in_dict['self.V_T'])
if state[2]: # modal params
ssi_object.modal_frequencies = validate_array(in_dict['self.modal_frequencies'])
ssi_object.modal_damping = validate_array(in_dict['self.modal_damping'])
ssi_object.eigenvalues = validate_array(in_dict['self.eigenvalues'])
ssi_object.mode_shapes = validate_array(in_dict['self.mode_shapes'])
ssi_object.modal_contributions = validate_array(in_dict['self.modal_contributions'])
ssi_object.max_model_order = validate_array(in_dict['self.max_model_order'])
ssi_object.state = state
return ssi_object
[docs]
class SSIData(SSIDataMC):
"""Data-driven SSI without correlation synthesis (non-Monte-Carlo variant).
Identical workflow to :class:`SSIDataMC` but skips the synthesis step
(``synth_sig=False``), making it faster when variance estimation is not
required.
"""
[docs]
def compute_modal_params(self, max_model_order): # pylint: disable=arguments-differ
'''
Perform a multi-order computation of modal parameters. Successively
calls
* estimate_state(order,)
* modal_analysis(A,C)
at ascending model orders, up to max_model_order.
See the explanations in the the respective methods, for a detailed
explanation of parameters.
Parameters
----------
max_model_order: integer, optional
Maximum model order, where to interrupt the algorithm. If not given,
it is determined from the previously computed subspace matrix.
'''
super().compute_modal_params(max_model_order, synth_sig=False)
[docs]
def estimate_state(self, order, max_modes=None, algo='svd'):
'''
Compute the state matrix A and output matrix C from the singular values
and vectors of the projection matrix, truncated at the requested order. Estimation of the
state matrix can be performed by QR decomposition or Singular Value decomposition
of the shifted observability matrix. If max_modes is specified, the singular
value decomposition is truncated additionally, also known as Crystal Clear SSI.
Parameters
----------
order: integer, required
Model order, at which the state matrices should be estimated
max_modes: integer, optional
Maximum number of modes, that are known to be present in the signal,
to suppress noise modes
algo: str, optional
Algorithm to use for estimation of A. Either 'svd' or 'qr'.
Returns
-------
A: numpy.ndarray
State matrix: Array of shape (order, order)
C: numpy.ndarray
Output matrix: Array of shape (num_analised_channels, order)
'''
if order > self.S.shape[0]:
raise RuntimeError(f'Order cannot be higher than {self.S.shape[0]}. Consider using more block_rows/block_columns.')
if algo not in ['svd', 'qr']:
raise ValueError(f"'algo' must be one of {['svd', 'qr']}, got {algo!r}.")
n_l = self.num_analised_channels
num_block_rows = self.num_block_rows
U = self.U[:,:order]
S = self.S[:order]
# compute state-space model
S_2 = np.power(S, 0.5)
O = U * S_2[np.newaxis,:]
On_up = O[:n_l * num_block_rows,:order]
On_down = O[n_l:n_l * (num_block_rows + 1),:order]
if algo == 'svd':
if max_modes is not None:
[u, s, v_t] = np.linalg.svd(On_up, 0)
s = 1. / s[:max_modes]
# On_up_i = np.dot(np.transpose(v_t[:max_modes, :]), np.multiply(
# s[:, np.newaxis], np.transpose(u[:, :max_modes])))
On_up_i = v_t[:max_modes,:].T @ (s[:, np.newaxis] * u[:,:max_modes].T)
else:
On_up_i = np.linalg.pinv(On_up) # , rcond=1e-12)
A = On_up_i @ On_down
elif algo == 'qr':
Q, R = np.linalg.qr(On_up)
S = Q.T.dot(On_down)
A = np.linalg.solve(R, S)
C = O[:n_l,:order] # output matrix
return A, C, None, None, None
[docs]
class SSIDataCV(SSIDataMC):
@staticmethod
def _coerce_blocks_array(blocks, num_blocks, name):
"""Validate and coerce a blocks argument to a numpy array."""
if blocks is None:
return np.arange(num_blocks)
if isinstance(blocks, (list, tuple)):
blocks = np.array(blocks)
elif not isinstance(blocks, np.ndarray):
raise RuntimeError(f"Argument {name!r} must be an iterable but is type {type(blocks)}")
if blocks.max() >= num_blocks:
raise ValueError(f"{name}.max() must be < {num_blocks}, got {blocks.max()}.")
return blocks
def _compute_lq_factors(self, hankel_matrices, training_blocks, K, K2, N_b, N, n_r, n_l, q, p, pbar):
"""Compute per-block LQ decompositions and assemble the combined unique Q matrix."""
n_training_blocks = training_blocks.shape[0]
R_matrices = np.empty((n_r * q + n_l * (p + 1), K * n_training_blocks))
Q_matrices = []
Q_unique_matrices = np.empty((K2, N))
for i in range(n_training_blocks):
i_block = training_blocks[i]
next(pbar)
L, Q = lq_decomp(hankel_matrices[i_block], mode='reduced', unique=True)
R_matrices[:, i * K:(i + 1) * K] = L
Q_matrices.append(Q)
logger.debug(f'R shapes: actual: {R_matrices.shape} expected: {(n_r * p + n_l * (p + 1), K * n_training_blocks)}')
R_full_breve, Q_full_breve = lq_decomp(R_matrices, mode='reduced', unique=True)
_ = [next(pbar) for _ in range(n_training_blocks)]
del R_matrices
logger.debug(f'Q_breve shapes: actual: {Q_full_breve.shape} expected: ,{(K2, K * n_training_blocks)}')
Q_breve_matrices = np.hsplit(Q_full_breve, np.arange(K, n_training_blocks * K, K))
logger.debug(f'Q_breve_j shapes: actual: {Q_breve_matrices[0].shape}, expected: {(K2, K)}')
del Q_full_breve
for i in range(n_training_blocks):
next(pbar)
Q_unique_matrices[:, i * N_b:(i + 1) * N_b] = Q_breve_matrices[i] @ Q_matrices[i]
return R_full_breve, Q_unique_matrices
[docs]
def build_block_hankel(self, num_block_rows=None, num_blocks=1, training_blocks=None, reduced_projection=True):
'''
Builds serveral Block-Hankel Matrices of the measured time series with varying
time lags and estimates the subspace matrices from their LQ decompositions.
Uniqueness of the subspace estimates is ensured by an intermediate LQ
decomposition, where the diagonals of the L matrices are constrained to
positive values. A subspace matrix estimate is computed by the mean over
the training blocks leaving any remainig blocks for validation.
Note: Blocks are not completely i.i.d. as we borrow p+q timesteps from the
previous block for the projection of a full block (assembly of Hankel matrix)
.. TODO::
* investigate correct scaling of the subspace matrices
[sqrt(N_b), sqrt(N_b * num_blocks), sqrt(N_b*n_training_blocks)] ?
* use sparse SVD (scipy.sparse.svds) to save memory and cpu time
Parameters
-------
num_block_rows: integer, required
The number of block rows of the Subspace matrix
num_blocks: integer, optional
The number of blocks, used for cross-validation
training_blocks: list, optional
The selected blocks to use for system identification (=training)
'''
if num_block_rows is None:
num_block_rows = self.num_block_rows
if not isinstance(num_block_rows, int):
raise TypeError(f"Expected int for 'num_block_rows', got {type(num_block_rows).__name__!r}.")
if not isinstance(num_blocks, int):
raise TypeError(f"Expected int for 'num_blocks', got {type(num_blocks).__name__!r}.")
training_blocks = self._coerce_blocks_array(training_blocks, num_blocks, 'training_blocks')
n_training_blocks = training_blocks.shape[0]
self.num_block_rows = num_block_rows
self.num_blocks = num_blocks
signals = self.prep_signals.signals
total_time_steps = self.prep_signals.total_time_steps
ref_channels = sorted(self.prep_signals.ref_channels)
n_l = self.prep_signals.num_analised_channels
n_r = self.prep_signals.num_ref_channels
q = num_block_rows
p = num_block_rows
logger.info(f'Building Block-Hankel matrix from {n_training_blocks} out of {num_blocks} signal blocks with {p} block-columns and {q} block rows.')
N_b = int(np.floor((total_time_steps - q - p) / num_blocks))
if N_b < n_r * q:
raise RuntimeError(f'Block-length ({N_b}) must not be smaller than the number of reference channels * number of block rows (={n_r * q}).')
# might omit some timesteps in favor of equally sized blocks
N = N_b * num_blocks
# shorten signals by omitted samples to have them available for Kalman-Filter startup later
N_offset = total_time_steps - q - p - N
signals = signals[N_offset:,:]
K = min((q * n_r) + (p + 1) * n_l, N_b)
K2 = min(((q * n_r) + (p + 1) * n_l, K * n_training_blocks))
Y_minus = np.zeros((q * n_r, N))
Y_plus = np.zeros(((p + 1) * n_l, N))
for ii in range(q):
Y_minus[(q - ii - 1) * n_r:(q - ii) * n_r,:] = signals[ii:ii + N, ref_channels].T
for ii in range(p + 1):
Y_plus[ii * n_l:(ii + 1) * n_l,:] = signals[q + ii:q + ii + N,:].T
Hankel_matrix = np.vstack((Y_minus, Y_plus))
Hankel_matrix /= np.sqrt(N)
logger.debug(Hankel_matrix.shape)
hankel_matrices = np.hsplit(Hankel_matrix, np.arange(N_b, N_b * num_blocks, N_b))
np.hstack([hankel_matrices[i_block] for i_block in training_blocks])
pbar = simplePbar(n_training_blocks * 3)
L, Q = self._compute_lq_factors(hankel_matrices, training_blocks, K, K2, N_b, N, n_r, n_l, q, p, pbar)
logger.info('Estimating subspace matrix...')
a = n_r * q
b = n_r
c = n_l - n_r
d = n_l * p
P_i_ref = L[a:a + b + c + d,: a] @ Q[:a,:] # (p + 1) * n_l x n_r * q x N
if reduced_projection:
[U, S, V_T] = np.linalg.svd(L[a:a + b + c + d,: a], full_matrices=False)
V_T = None
else:
[U, S, V_T] = np.linalg.svd(P_i_ref, full_matrices=False)
P_i_1 = L[a + b + c:a + b + c + d,: a + b] @ Q[: a + b,: ] # n_l * p x n_r * (q + 1) x N
Y_i_i = L[a: a + b + c,: a + b + c] @ Q[: a + b + c,: ] # n_l x n_r *q + n_l x N
self.L_red = L[a:,:a + b] # (p + 1) n_l x n_r * (q + 1)
self.Q_red = Q[:a + b,:] # n_r * (q + 1) x N
self.P_i_1 = P_i_1
self.P_i_ref = P_i_ref
self.Y_i_i = Y_i_i
self.S = S
self.U = U
self.V_T = V_T
# self.max_model_order = self.S.shape[0]
self.state[0] = True
def _run_kalman_block(self, i, validation_blocks, block_starts, signals, N_b, N_0_offset,
N_offset, order, n_l, AKC, K_0, start_states):
"""Run the Kalman filter for a single validation block and return states and measurements."""
i_block = validation_blocks[i]
start_state = start_states[i]
if i_block == 0:
_N_offset = N_0_offset
elif start_state is not None:
_N_offset = 1
else:
_N_offset = N_offset
states = np.zeros((order, N_b + _N_offset), dtype=complex)
block_start = block_starts[i] - _N_offset
block_end = block_starts[i] + N_b
signals_block = signals[block_start:block_end,:]
K_0m = K_0 @ signals_block.T
if start_state is not None:
states[:, 0] = start_state
for k in range(N_b + _N_offset - 1):
states[:, k + 1] = K_0m[:, k] + AKC @ states[:, k]
start_states[i + 1] = states[:, -1]
states = states[:, _N_offset:]
Y = signals_block[_N_offset:,:].T
return states, Y
def _decompose_modes(self, states, Y, sig_synth, eigvals, C_0, conj_indices, order, n_l):
"""Compute per-mode synthesized signals and return Sigma_data_synth."""
Sigma_data = np.einsum('ji,ji->j', Y, Y)
Sigma_data_synth = np.zeros((n_l, order))
for i, ind in enumerate(conj_indices):
lambda_i = eigvals[ind]
ident = eigvals == lambda_i.conj()
ident[ind] = 1
C_0I = C_0[:, ident]
this_sig_synth = C_0I @ states[ident,:]
if not np.all(np.isclose(this_sig_synth.imag, 0)):
logger.warning(f'Synthetized signals are complex at mode index {order}:{ind}.')
sig_synth[:,:, i] = this_sig_synth.real
Sigma_data_synth[:, i] = np.einsum('ji,ji->j', sig_synth[:,:, i], Y)
total_sig_synth = np.sum(sig_synth, axis=-1)
Sigma_synth = np.einsum('ji,ji->j', total_sig_synth, total_sig_synth)
contrib = np.mean(Sigma_data_synth / np.sqrt(Sigma_data * Sigma_synth)[:, np.newaxis], axis=0)
return sig_synth, contrib
[docs]
def synthesize_signals(self, A, C, Q, R, S, validation_blocks=None, N_offset=None, **kwargs):
'''
Computes the modal response signals and the contribution of each mode.
The algorithm follows Peeters 1999 and the Lyapunov equation is solved
as a discrete-time algebraic Riccati equation (DARE). For long signals,
the computation may become time-consuming, thus only time steps up to j
may be used to synthesize the signal.
Parameters
----------
A: numpy.ndarray
State matrix: Array of shape (order, order)
C: numpy.ndarray
Output matrix: Array of shape (num_analised_channels, order)
Q: numpy.ndarray
state noise covariance matrix: Symmetric array of shape (order, order)
R: numpy.ndarray
signal noise covariance matrix: Array of shape (num_analised_channels, num_analised_channels)
S: numpy.ndarray
system noise - signal noise covariance matrix: Array of shape (order, num_analised_channels)
validation_blocks: list, optional
The selected blocks to be synthethized and used for system validation.
N_offset: integer, optional
The number of samples to be used from any previous block for
Kalman-Filter startup.
Returns
-------
sig_synth: numpy.ndarray, shape (num_analised_channels, j, order // 2)
Array holding the modally decomposed input signals for
each channel n_l and all modes
modal_contributions: numpy.ndarray, shape (order, )
Array holding the contributions of each mode to the input
signals.
'''
order = A.shape[0]
if order != A.shape[1]:
raise RuntimeError(f"Internal error: A must be square, got shape {A.shape}.")
num_blocks = self.num_blocks
validation_blocks = self._coerce_blocks_array(validation_blocks, num_blocks, 'validation_blocks')
n_validation_blocks = validation_blocks.shape[0]
n_l = self.prep_signals.num_analised_channels
signals = self.prep_signals.signals
total_time_steps = self.prep_signals.total_time_steps
q = self.num_block_rows
p = self.num_block_rows
N_b = int(np.floor((total_time_steps - q - p) / num_blocks))
# might omit some timesteps in favor of equally sized blocks
N = N_b * num_blocks
# blocks start at N_0_offset + p + q
# (in training we virtually borrow p + q timesteps from the previous block)
N_0_offset = total_time_steps - N
if N_offset is None:
N_offset = N_b // 5
if 0 in validation_blocks and N_0_offset < N_offset:
logger.info(f"Block '0' is in the validation dataset, but only has {N_0_offset} startup-samples (recommended/chosen: {N_offset}) from any previous block for the Kalman Filter. Expect a degraded performance.")
modal_contributions = np.zeros((order))
all_sig_synth = []
try:
P = scipy.linalg.solve_discrete_are(
a=A.T, b=C.T, q=Q, r=R, s=S, balanced=True)
except Exception:
logger.warning('Correlations of residuals are not symmetric. Skiping Modal Contributions')
return all_sig_synth, modal_contributions
APCS = A @ P @ C.T + S
CPCR = C @ P @ C.T + R
K = self._solve_kalman_gain(CPCR, APCS)
eigvals, eigvecs_r = np.linalg.eig(A)
conj_indices = self.remove_conjugates(eigvals, eigvecs_r, inds_only=True)
C_0 = C @ eigvecs_r
K_0 = np.linalg.solve(eigvecs_r, K)
AKC = np.diag(eigvals) - K_0 @ C_0
block_starts = validation_blocks * N_b + N_0_offset
start_states = [None for _ in range(num_blocks + 1)]
for i in np.argsort(validation_blocks):
sig_synth = np.zeros((n_l, N_b, order // 2))
states, Y = self._run_kalman_block(
i, validation_blocks, block_starts, signals, N_b, N_0_offset,
N_offset, order, n_l, AKC, K_0, start_states)
sig_synth, contrib = self._decompose_modes(
states, Y, sig_synth, eigvals, C_0, conj_indices, order, n_l)
modal_contributions += contrib
all_sig_synth.append(sig_synth)
modal_contributions /= n_validation_blocks
self._sig_synth = all_sig_synth
self._modal_contributions = modal_contributions
return all_sig_synth, modal_contributions
def plot_sig_synth(modal_data, modelist=None, channel_inds=None, ref_channel_inds=None, axes=None, i_block=None):
import matplotlib.pyplot as plt
prep_signals = modal_data.prep_signals
sig_synth = modal_data._sig_synth
if isinstance(sig_synth, list): # multi-block for cross-validation SSIDataCV
sig_synth = sig_synth[i_block]
num_blocks = modal_data.num_blocks
total_time_steps = prep_signals.total_time_steps
q = modal_data.num_block_rows
p = modal_data.num_block_rows
N_b = int(np.floor((total_time_steps - q - p) / num_blocks))
N = N_b * num_blocks
N_0_offset = total_time_steps - q - p - N
# N_offset = N_b // 15
block_starts = i_block * N_b + N_0_offset + p + q
t = prep_signals.t[block_starts:block_starts + N_b]
signals = prep_signals.signals.T[:, block_starts:block_starts + N_b]
else:
t = prep_signals.t
signals = prep_signals.signals.T
# ref_channels = modal_data.prep_signals.ref_channels
sampling_rate = modal_data.prep_signals.sampling_rate
channel_headers = modal_data.prep_signals.channel_headers
modal_contributions = modal_data._modal_contributions
if channel_inds is None:
channel_inds = np.arange(modal_data.prep_signals.num_analised_channels)
num_channels = len(channel_inds)
if ref_channel_inds is None:
ref_channel_inds = np.arange(modal_data.prep_signals.num_ref_channels)
# num_ref_channels = len(ref_channel_inds)
# Plot signals for each mode and each channel
num_modes = sig_synth.shape[-1]
if modelist is None:
modelist = list(range(num_modes))
num_plots = len(modelist) + 2
fig1, axes = plt.subplots(num_plots, 1, sharex='col', sharey='col', squeeze=False)
for ip, i in enumerate(modelist):
rho = modal_contributions[i]
this_signals = sig_synth[:,:, i]
for j in range(num_channels):
i_l = channel_inds[j]
color = str(np.linspace(0, 1, len(channel_inds) + 2)[j + 1])
ls = 'solid'
axes[ip, 0].plot(t, this_signals[i_l,:], color=color, ls=ls)
axes[ip, 0].set_ylabel(f'$\\delta_{{{i + 1}}}$={rho:1.2f}',
rotation=0, labelpad=40, va='center', ha='left')
for j in range(num_channels):
i_l = channel_inds[j]
color = str(np.linspace(0, 1, len(channel_inds) + 2)[j + 1])
ls = 'solid'
this_signals = signals[i_l,:]
this_signals_synth = np.sum(sig_synth, axis=2)[i_l,:]
axes[-1, 0].plot(t, this_signals, color=color, ls=ls,
label=f'{channel_headers[i_l]}')
axes[-2, 0].plot(t, this_signals_synth, color=color, ls=ls,)
axes[-1, 0].set_ylabel('Measured', rotation=0, labelpad=50, va='center', ha='left')
axes[-2, 0].set_ylabel(f'$\\sum\\delta$={np.sum(modal_contributions):1.2f}',
rotation=0, labelpad=50, va='center', ha='left')
axes[-1, 0].set_xlabel('$t$ [\\si{\\second}]')
for ax in axes.flat:
ax.set_yticks([])
fig1.legend(title='Channels')
# Plot power spectral density functions for each channel and all modes
fig2, axes = plt.subplots(num_channels, 1, sharex='col', sharey='col', squeeze=False)
ft_freq = np.fft.rfftfreq(len(t), d=(1 / sampling_rate))
for j in range(num_channels):
i_l = channel_inds[j]
this_signals = signals[i_l,:]
this_signals_synth = sig_synth[i_l,:,:]
ft_meas = np.fft.rfft(this_signals * np.hanning(len(t)))
if j == 0: label = f'Inp.'
else: label = None
axes[j, 0].plot(ft_freq, 10 * np.log10(np.abs(ft_meas)), ls='solid', color='k', label=label)
for ip, i in enumerate(modelist):
ft_synth = np.fft.rfft(this_signals_synth[:, i] * np.hanning(len(t)))
color = str(np.linspace(0, 1, len(modelist) + 2)[ip + 1])
ls = ['-', '--', ':', '-.'][i % 4]
if j == 0: label = f'm={i+1}'
else: label = None
axes[j, 0].plot(ft_freq, 10 * np.log10(np.abs(ft_synth)), color=color, ls=ls, label=label)
axes[j, 0].set_ylabel(f'{channel_headers[i_l]}',
rotation=0, labelpad=20, va='center', ha='center')
axes[-1, 0].set_xlabel('$f$ [\\si{\\hertz}]')
for ax in axes.flat:
ax.set_yticks([])
ax.set_xlim(0, 1 / 2 * sampling_rate)
fig2.legend(title='Mode')
return fig1, fig2
def main():
pass
if __name__ == '__main__':
main()
