"""Clean linear-time-invariant state-space representation.
A :class:`StateSpace` holds the validated ``A``, ``B``, ``C``, ``D`` matrices of a
continuous- or discrete-time system together with the state, input, and output
names. The averaged matrices of an MBC transform are exported directly into this
form, and a zero-order-hold discretization is provided for use in a digital twin.
"""
from __future__ import annotations
from dataclasses import dataclass, field
from typing import TYPE_CHECKING
import numpy as np
from scipy.linalg import expm
if TYPE_CHECKING:
import numpy.typing as npt
from vane.io.mbc_transform import MBCResult
__all__ = ["StateSpace", "state_space_from_mbc"]
[docs]
@dataclass
class StateSpace:
"""A linear time-invariant state-space system.
Parameters
----------
a : npt.NDArray[np.float64]
State matrix, shape ``(n_states, n_states)``.
b, c, d : npt.NDArray[np.float64] or None
Input, output, and feedthrough matrices, shapes ``(n_states, n_inputs)``,
``(n_outputs, n_states)``, ``(n_outputs, n_inputs)``; any may be ``None``.
state_names, input_names, output_names : list[str]
Names of the states, inputs, and outputs; each is either empty or matches
the corresponding dimension.
dt : float or None
Sample time; ``None`` for a continuous-time system, otherwise the discrete
sample time in seconds.
Raises
------
ValueError
If the matrix shapes or name lengths are inconsistent.
"""
a: npt.NDArray[np.float64]
b: npt.NDArray[np.float64] | None = None
c: npt.NDArray[np.float64] | None = None
d: npt.NDArray[np.float64] | None = None
state_names: list[str] = field(default_factory=list)
input_names: list[str] = field(default_factory=list)
output_names: list[str] = field(default_factory=list)
dt: float | None = None
def __post_init__(self) -> None:
"""Validate matrix shapes, the sample time, and name lengths."""
if self.a.ndim != 2 or self.a.shape[0] != self.a.shape[1]:
msg = f"A must be square, got shape {self.a.shape}"
raise ValueError(msg)
n = self.a.shape[0]
if self.b is not None:
self._require_matrix(self.b, "B")
if self.b.shape[0] != n:
msg = f"B must have {n} rows, got {self.b.shape}"
raise ValueError(msg)
if self.c is not None:
self._require_matrix(self.c, "C")
if self.c.shape[1] != n:
msg = f"C must have {n} columns, got {self.c.shape}"
raise ValueError(msg)
self._validate_feedthrough()
if self.dt is not None and not self.dt > 0.0:
msg = f"Discrete sample time must be positive, got {self.dt}"
raise ValueError(msg)
self._check_names(self.state_names, n, "state_names")
self._check_names(self.input_names, self.n_inputs, "input_names")
self._check_names(self.output_names, self.n_outputs, "output_names")
def _validate_feedthrough(self) -> None:
"""Validate D against whichever input/output dimensions are known."""
if self.d is None:
return
self._require_matrix(self.d, "D")
if self.b is not None and self.d.shape[1] != self.b.shape[1]:
msg = f"D must have {self.b.shape[1]} columns (inputs), got {self.d.shape}"
raise ValueError(msg)
if self.c is not None and self.d.shape[0] != self.c.shape[0]:
msg = f"D must have {self.c.shape[0]} rows (outputs), got {self.d.shape}"
raise ValueError(msg)
@staticmethod
def _require_matrix(matrix: npt.NDArray[np.float64], label: str) -> None:
"""Raise if ``matrix`` is not two-dimensional."""
if matrix.ndim != 2:
msg = f"{label} must be a 2-D matrix, got a {matrix.ndim}-D array"
raise ValueError(msg)
@staticmethod
def _check_names(names: list[str], expected: int, label: str) -> None:
"""Raise if a non-empty name list does not match the expected length."""
if names and len(names) != expected:
msg = f"{label} has {len(names)} entries, expected {expected}"
raise ValueError(msg)
@property
def n_states(self) -> int:
"""Return the number of states."""
return int(self.a.shape[0])
@property
def n_inputs(self) -> int:
"""Return the number of inputs (0 if B is absent)."""
return 0 if self.b is None else int(self.b.shape[1])
@property
def n_outputs(self) -> int:
"""Return the number of outputs (0 if C is absent)."""
return 0 if self.c is None else int(self.c.shape[0])
@property
def is_discrete(self) -> bool:
"""Return whether the system is discrete-time."""
return self.dt is not None
[docs]
def is_stable(self) -> bool:
"""Return whether the system is asymptotically stable.
Returns
-------
bool
For a continuous system, whether every eigenvalue of ``A`` has
negative real part; for a discrete system, whether every eigenvalue
lies strictly inside the unit circle.
"""
eigenvalues = np.linalg.eigvals(self.a)
if self.is_discrete:
return bool(np.all(np.abs(eigenvalues) < 1.0))
return bool(np.all(eigenvalues.real < 0.0))
[docs]
def discretized(self, dt: float) -> StateSpace:
"""Return a zero-order-hold discretization at sample time ``dt``.
Parameters
----------
dt : float
Sample time in seconds; must be positive.
Returns
-------
StateSpace
The discrete-time system; ``C`` and ``D`` are unchanged.
Raises
------
ValueError
If ``dt`` is not positive or the system is already discrete.
"""
if not dt > 0.0:
msg = f"Sample time must be positive, got {dt}"
raise ValueError(msg)
if self.is_discrete:
msg = "System is already discrete"
raise ValueError(msg)
ad, bd = _zero_order_hold(self.a, self.b, dt)
return StateSpace(
a=ad,
b=bd,
c=self.c,
d=self.d,
state_names=list(self.state_names),
input_names=list(self.input_names),
output_names=list(self.output_names),
dt=dt,
)
[docs]
def state_space_from_mbc(result: MBCResult) -> StateSpace:
"""Export the averaged MBC matrices as a continuous-time state-space system.
Parameters
----------
result : MBCResult
An MBC transform result whose ``avg_a`` is present.
Returns
-------
StateSpace
The continuous-time system with state names from the result.
Raises
------
ValueError
If ``result.avg_a`` is ``None``.
"""
if result.avg_a is None:
msg = "MBCResult has no averaged A matrix to export"
raise ValueError(msg)
return StateSpace(
a=result.avg_a,
b=result.avg_b,
c=result.avg_c,
d=result.avg_d,
state_names=list(result.state_descriptions),
)
def _zero_order_hold(
a: npt.NDArray[np.float64], b: npt.NDArray[np.float64] | None, dt: float
) -> tuple[npt.NDArray[np.float64], npt.NDArray[np.float64] | None]:
"""Zero-order-hold discretization of ``(A, B)`` via the augmented exponential."""
n = a.shape[0]
if b is None:
ad: npt.NDArray[np.float64] = np.asarray(expm(a * dt), dtype=np.float64)
return ad, None
m = b.shape[1]
augmented = np.zeros((n + m, n + m), dtype=np.float64)
augmented[:n, :n] = a
augmented[:n, n:] = b
exponential = np.asarray(expm(augmented * dt), dtype=np.float64)
return exponential[:n, :n], exponential[:n, n:]