diff --git a/docs/_static/images/gallery/Example5.png b/docs/_static/images/gallery/Example5.png new file mode 100755 index 000000000..42e861209 Binary files /dev/null and b/docs/_static/images/gallery/Example5.png differ diff --git a/docs/_static/images/gallery/Example6.png b/docs/_static/images/gallery/Example6.png new file mode 100755 index 000000000..700c96125 Binary files /dev/null and b/docs/_static/images/gallery/Example6.png differ diff --git a/docs/_static/images/gallery/safeway.png b/docs/_static/images/gallery/safeway.png new file mode 100755 index 000000000..74c5c19dd Binary files /dev/null and b/docs/_static/images/gallery/safeway.png differ diff --git a/docs/_templates/home.html b/docs/_templates/home.html index fcf429411..ca1e18068 100644 --- a/docs/_templates/home.html +++ b/docs/_templates/home.html @@ -281,7 +281,7 @@

Gallery

{% for example in examples %} + href="{{example.link}}">
{{ example.title }}
{{ example.description }}
diff --git a/docs/classes/Kalman_LDS.md b/docs/classes/Kalman_LDS.md deleted file mode 100644 index 6d0f2ea43..000000000 --- a/docs/classes/Kalman_LDS.md +++ /dev/null @@ -1,224 +0,0 @@ ---- -title: Kalman's Mathematical Description of Linear Dynamical Systems -author: Chrystal Chern -date: Friday, December 15, 2023 -... - - - -## 0. Question and Scope - -```{=tex} -\begin{centering} -``` -*How much of the physical world can be determined from a given amount of experimental data?* -```{=tex} -\end{centering} -``` - -- Real, finite dimensional, continuous-time, linear time-invariant (LTI) dynamical systems (LTV before Section 5). -- Impulse response. - -## 1. Intro and Summary - -- **1&2&3**: Intro; define a dynamical system -- **4**: Present the problem of system realization from impulse response -- **5**: Present the *canonical structure theorem*: if a system realization is controllable and observable, $\implies$ it is an irreducible realization of the impulse response. -- **6&7**: Define controllability and observability and compute the canonical structure of an LTI system. -- **8**: Define MIMO state variables for LTI systems from transfer functions. -- **9**: Present common errors made by equating transfer functions and state space system realizations. - - -## 2. Axioms of Dynamical Systems - -A dynamical system consists of the following (developed for linear time-varying (LTV) systems by Kalman, but only presented for LTI systems here): - -### *state*, $\begin{bmatrix} \bm{x}(t) \\ \dot{\bm{x}}(t) \end{bmatrix} \in \Sigma = \mathbb{R}^{n}$ - -- Physical interpretation: instantaneous position and momentum. -- Definition: minimal information needed about the system's history which suffices to predict the effect of the past upon the future. -- Properties of state: $\begin{bmatrix} \bm{x}(t>\tau) \\ \dot{\bm{x}}(t>\tau) \end{bmatrix}$ is fully defined by $\begin{bmatrix} \bm{x}(\tau) \\ \dot{\bm{x}}(\tau) \end{bmatrix}$ and $\bm{u}(t\geq\tau)$. - -### *time*, $t \in \Theta = \mathbb{R}$ - -### *input*, $u \in \Omega =$ piecewise continuous functions - -- Physical interpretation: forcing function. - -### *transition function*, $\varphi: \Omega \times \Theta \times \Theta \times \Sigma \to \Sigma$ - -$$ -\bm{x}(t) = \varphi(\bm{u}, t; t_{o}, \bm{x}_{o}) -$$ - -- Linear on $\Omega \times \Theta$. -- Continuous with respect to $\Sigma, \Theta, \Omega$ and their induced products. - -### *output function*, $\psi: \Theta \times \Sigma \to \mathbb{R}$ - -- Physical interpretation: variables that are directly observed. -- Linear on $\Sigma$. -- Continuous with respect to $\Sigma, \Theta, \Omega$ and their induced products. - -Every linear time-invariant ($\bm{A,B,C}$ linear in $\Sigma$ and independent of time) dynamical system is thus governed by the equations[^1] - -$$ -\boxed{ - \begin{aligned} - \frac{d}{dt}\bm{x}(t) &= \bm{Ax}(t) + \bm{Bu}(t) \\ - \bm{y}(t) &= \bm{Cx}(t) - \end{aligned} -} -$${@eq:goveqs} - -The general solution to the differential equation is - -$$ -\begin{aligned} -\bm{x}(t) &= \Phi_{\bm{A}}(t,t_{o})\bm{x}_{o} + \int_{t_{o}}^{t}{\Phi_{\bm{A}}(t,\tau)\bm{Bu}(\tau)d\tau} -\\ -\Phi_{\bm{A}}(t,\sigma) &= e^{\bm{A}(t-\sigma)} -\end{aligned} -$${#eq:gensol} - -and it can be verified that $\Phi_{\bm{A}}(t,\sigma) = \Phi_{\bm{A}}(t,\tau)\Phi_{\bm{A}}(\tau,\sigma) ~~\forall{t,\tau,\sigma} \in \mathbb{R}$. - -[^1]: Sometimes you'll see a system with a *feed-through* term, $\bm{D}$ in the output function, $\bm{y}(t) = \bm{Cx}(t) + \bm{Du}(t)$. This term can always be added or removed because inputs $\bm{u}(t)$ and outputs $\bm{y}(t)$ are deterministic quantities. - -## 3. Equivalent Dynamical Systems - -The state vector $\bm{x}(t)$ is an abstract quantity in $\mathbb{R}^{n}$. Therefore it can be expressed in any $\mathbb{R}^{n}$ coordinates. The system with the state vector $\bar{\bm{x}}(t)$ is equivalent to the system with the state vector $\bm{x}(t)$ if there exists a nonsingular matrix $\bm{T} \in \mathbb{R}^{n\times n}$ such that - -$$ -\bar{\bm{x}}(t) = \bm{Tx}(t)~. -$$ - -The equivalent governing equations for the system are as follows: - -$$ -\begin{aligned} -\frac{d}{dt}\bar{\bm{x}}(t) &= \bm{TAT^{-1}\bar{x}}(t) + \bm{TBu}(t) \\ -\bm{y}(t) &= \bm{CT^{-1}\bar{x}}(t) -\end{aligned} -$$ - -## 4. System Realization from Impulse Response - -The *impulse response matrix* of a system is a time-dependent array of the output at each $\bm{y}$ coordinate in response to an pulse input, $\delta(t-t_{o})$, at each $\bm{u}$ coordinate. - -That is, given the input $u_{ij}(t) = \delta_{ij}\delta(t-t_{o})$, the output is $y_{ij} = S_{ij}(t,t_{o})$, where (by plugging in (@eq:gensol),) - -$$ -\bm{S}(t,\tau) = \bm{C} \Phi_{\bm{A}}(t,\tau) \bm{B} ~. -$$ - -Given $\bm{S}(t,t_{o})$ for a dynamical system, the output can be found for any input by the convolution integral: - -$$ -\bm{y}(t) = \int_{t_{o}}^{t}{\bm{S}(t,\tau)\bm{u}(\tau)d\tau} ~. -$$ - ->The central question of this paper is, -> ->```{=tex} ->\begin{centering} ->``` ->*When and how does the impulse-response matrix determine the dynamical equations of the system?* ->```{=tex} ->\end{centering} ->``` - -$$ -\begin{aligned} -\\ -\\ -\end{aligned} -$$ - -## 5. Kalman Canonical Staircase State Space Realization (K-CSSSR) - -Kalman presents his canonical form for state space system realization of a dynamical system. He proves that any system realization can be transformed into one which isolates the state vector into four parts: - -- controllable and unobservable ($\bm{x}_{cuo}$) -- controllable and observable ($\bm{x}_{co}$) -- uncontrollable and unobservable ($\bm{x}_{ucuo}$) -- ucontrollable and observable ($\bm{x}_{uco}$) - -The governing equations are as follows: - -$$ -\boxed{ - \begin{aligned} - \frac{d}{dt}\begin{bmatrix} - \bm{x}_{cuo}(t) \\ - \bm{x}_{co}(t) \\ - \bm{x}_{ucuo}(t) \\ - \bm{x}_{uco}(t) \\ - \end{bmatrix} &= - \begin{bmatrix} - \bm{A}_{cuo} & \bm{A}_{12} & \bm{A}_{13} & \bm{A}_{14} \\ - 0 & \bm{A}_{co} & 0 & \bm{A}_{24}\\ - 0 & 0 & \bm{A}_{ucuo} & \bm{A}_{34} \\ - 0 & 0 & 0 & \bm{A}_{uco} \\ - \end{bmatrix} - \begin{bmatrix} - \bm{x}_{cuo}(t) \\ - \bm{x}_{co}(t) \\ - \bm{x}_{ucuo}(t) \\ - \bm{x}_{uco}(t) \\ - \end{bmatrix} + - \begin{bmatrix} - \bm{B}_{cuo}(t) \\ - \bm{B}_{co}(t) \\ - 0 \\ - 0 \\ - \end{bmatrix} - \bm{u}(t) \\ - \bm{y}(t) &= - \begin{bmatrix} - 0 & \bm{C}_{co} & 0 & \bm{C}_{uco} - \end{bmatrix} - \bm{x}(t) - \end{aligned} -} -$${#eq:canonical} - -The impulse response matrix of a linear dynamical system depends only on the dynamics of the controllable and observable modes $\bm{x}_{co}(t)$. - -$$ -\bm{S}(t,\tau) = \bm{C}_{co}\Phi_{\bm{A}_{co}}(t,\tau)\bm{B}_{co} -$$ - -Equivalently, any controllable and observable realization of an impulse response matrix is *irreducible* and equivalent to any other controllable and observable realization. - -If a dynamical system is controllable and observable, it has a unique impulse response matrix. - ->The answer to the central question of this paper is, ->```{=tex} ->\begin{centering} ->``` ->*The impulse-response matrix fully determines the dynamical equations of a controllable and observable system.* ->```{=tex} ->\end{centering} ->``` - - -## 6. Definition of Controllability and Observability - -The following statements are equivalent: - -- The system is controllable. -- The pair, $\{\bm{A},\bm{B}\}$, is controllable. -- The matrix $\bm{P} = \begin{bmatrix} \bm{B} & \bm{AB} & \cdots & \bm{A^{n-1}B} \end{bmatrix}$ has rank $n$. - -The following statements are equivalent: - -- The system is observable. -- The pair, $\{\bm{A},\bm{C}\}$, is observable. -- The matrix $\bm{Q} = \begin{bmatrix} \bm{C} \\ \bm{CA} \\ \cdots \\ \bm{CA^{n-1}} \end{bmatrix}$ has rank $n$. - -## 7. Computing the Canonical Form (@eq:canonical) - -The transformation matrix that separates us into controllable and uncontrollable modes is $\bm{M} = \begin{bmatrix} \bm{M}_{c} & \bm{M}_{uc} \end{bmatrix}$, where the columns of $\bm{M}_{c}$ are the linearly independent columns of $\bm{P}$ and the columns of $\bm{M}_{uc}$ are chosen such that $\bm{M}$ is full rank. - -The transformation matrix that separates us into observable and unobservable modes is $\bm{N} = \begin{bmatrix} \bm{N}_{o} \\ \bm{N}_{uo} \end{bmatrix}$, where the rows of $\bm{N}_{o}$ are the linearly independent rows of $\bm{Q}$ and the rows of $\bm{N}_{uo}$ are chosen such that $\bm{N}$ is full rank. \ No newline at end of file diff --git a/docs/classes/era.rst b/docs/classes/era.rst deleted file mode 100644 index 314a38e75..000000000 --- a/docs/classes/era.rst +++ /dev/null @@ -1,141 +0,0 @@ - -Eigensystem Realization Algorithm (ERA) ---------------------------------------- - - -As shown in the previous section, a structural system’s dynamic behavior -can be represented by the four coefficients -(:math:`\mathbf{A},\mathbf{B},\mathbf{C},\mathbf{D}`) of its discrete -LTI state-space representation. The `Ho-Kalman -Algorithm `__, or -`Eigensystem Realization Algorithm `__, -produces a *reduced order model* for these four coefficients, -(:math:`\mathbf{\tilde{A}},\mathbf{\tilde{B}},\mathbf{\tilde{C}},\mathbf{\tilde{D}}`), -based on an impulse input and its corresponding response output. Then, -modal properties can be extracted from :math:`\mathbf{\tilde{A}}` and -:math:`\mathbf{\tilde{C}}`. - -With the discrete LTI model, a unit impulse input with zero initial -conditions produces an output of constants -(:math:`\mathbf{D,CB,CAB,...,CA}^{k-1}\mathbf{B}`). These constants are -called *Markov parameters* because they must be unique for a given -system – there is only one possible output for a unit impulse input. - -.. math:: - - - \begin{aligned} - \mathbf{x}_{k+1} &= \mathbf{Ax}_{k} + \mathbf{Bu}_{k} \\ - \mathbf{y}_{k} &= \mathbf{Cx}_{k} + \mathbf{Du}_{k} \\ - \end{aligned} - -.. math:: - - - \begin{aligned} - \mathbf{u}_{0},\mathbf{u}_{1},\mathbf{u}_{2},...,\mathbf{u}_{k} &= \mathbf{I,0,0,...,0} \\ - \mathbf{x}_{0},\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{k} &= \mathbf{0,B,AB,...,A}^{k-1}\mathbf{B} \\ - \mathbf{y}_{0},\mathbf{y}_{1},\mathbf{y}_{2},...,\mathbf{y}_{k} &= \mathbf{D,CB,CAB,...,CA}^{k-1}\mathbf{B} \\ - \end{aligned} - -Knowing that the impulse response output data directly give the Markov -parameters, the data can then be stacked into the generalized blockwise -Hankel matrix :math:`\mathbf{H}`: - -.. math:: - - - \mathbf{H} - = - \begin{bmatrix} - \mathbf{y}_{1} & \mathbf{y}_{2} & \cdots & \mathbf{y}_{m_{c}} \\ - \mathbf{y}_{2} & \mathbf{y}_{3} & \cdots & \mathbf{y}_{m_{c}+1} \\ - \vdots & \vdots & \ddots & \vdots \\ - \mathbf{y}_{m_{o}} & \mathbf{y}_{m_{o}+1} & \cdots & \mathbf{y}_{m_{o}+m_{c}-1} \\ - \end{bmatrix} - = - \begin{bmatrix} - \mathbf{CB} & \mathbf{CAB} & \cdots & \mathbf{CA}^{m_{c}-1}\mathbf{B} \\ - \mathbf{CAB} & \mathbf{CA}^{2}\mathbf{B} & \cdots & \mathbf{CA}^{m_{c}}\mathbf{B} \\ - \vdots & \vdots & \ddots & \vdots \\ - \mathbf{CA}^{m_{o}-1}\mathbf{B} & \mathbf{CA}^{m_{o}}\mathbf{B} & \cdots & \mathbf{CA}^{m_{c}+m_{o}-2}\mathbf{B} \\ - \end{bmatrix} - = - \mathbf{\mathcal{OC}} - -where :math:`\mathbf{\mathcal{O}}` and :math:`\mathbf{\mathcal{C}}` are -the observability and controllability matrices of the system: - -.. math:: - - - \mathbf{\mathcal{O}} = \begin{bmatrix} - \mathbf{C} \\ \mathbf{CA} \\ \mathbf{CA}^{2} \\ \vdots \\ \mathbf{CA}^{m_{o}-1} - \end{bmatrix}, \hspace{1cm} - \mathbf{\mathcal{C}} = \begin{bmatrix} \mathbf{B} & \mathbf{AB} & \mathbf{A}^{2}\mathbf{B} & \cdots & \mathbf{A}^{m_{c}-1}\mathbf{B} \end{bmatrix} - -The shifted Hankel matrix, :math:`\mathbf{H'}` (one time step ahead of -:math:`\mathbf{H}`), is shown below: - -.. math:: - - - \mathbf{H'} - = - \begin{bmatrix} - \mathbf{y}_{2} & \mathbf{y}_{3} & \cdots & \mathbf{y}_{m_{c}+1} \\ - \mathbf{y}_{3} & \mathbf{y}_{4} & \cdots & \mathbf{y}_{m_{c}+2} \\ - \vdots & \vdots & \ddots & \vdots \\ - \mathbf{y}_{m_{o}+1} & \mathbf{y}_{m_{o}+2} & \cdots & \mathbf{y}_{m_{o}+m_{c}} \\ - \end{bmatrix} - = - \begin{bmatrix} - \mathbf{CAB} & \mathbf{CA}^{2}\mathbf{B} & \cdots & \mathbf{CA}^{m_{c}}\mathbf{B} \\ - \mathbf{CA}^{2}\mathbf{B} & \mathbf{CA}^{3}\mathbf{B} & \cdots & \mathbf{CA}^{m_{c}+1}\mathbf{B} \\ - \vdots & \vdots & \ddots & \vdots \\ - \mathbf{CA}^{m_{o}}\mathbf{B} & \mathbf{CA}^{m_{o}+1}\mathbf{B} & \cdots & \mathbf{CA}^{m_{c}+m_{o}-1}\mathbf{B}\\ - \end{bmatrix} - = - \mathbf{\mathcal{O}A\mathcal{C}} - -By taking the dominant terms of the singular value decomposition -of :math:`\mathbf{H}`, and transforming the relationship between -:math:`\mathbf{H} = \mathbf{\mathcal{OC}}` and -:math:`\mathbf{H'} = \mathbf{\mathcal{O}A\mathcal{C}}`, a reduced-order -model is constructed as follows: - -.. math:: - - - \mathbf{H} = \mathbf{U}\Sigma\mathbf{V}^{H} = - \begin{bmatrix} \mathbf{\tilde{U}} & \mathbf{U}_{t} \end{bmatrix} - \begin{bmatrix} \tilde{\Sigma} & \mathbf{0} \\ \mathbf{0} & \Sigma_{t} \end{bmatrix} - \begin{bmatrix} \mathbf{\tilde{V}}^{H} \\ \mathbf{V}_{t}^{H} \end{bmatrix} - \approx \mathbf{\tilde{U}}\tilde{\Sigma}\mathbf{\tilde{V}}^{H} - -where the superscript :math:`H` denotes conjugate transpose and the -subscript :math:`t` indicates elements to be truncated such that only -the first :math:`r` dominant singular values in :math:`\tilde{\Sigma}` -are retained, - -.. math:: - - - \begin{aligned} - \mathbf{\tilde{A}} &= \tilde{\Sigma}^{-1/2}\mathbf{\tilde{U}}^{H}\mathbf{H'\tilde{V}}\tilde{\Sigma}^{-1/2} \\ - \mathbf{\tilde{B}} &= \tilde{\Sigma}^{1/2}\mathbf{\tilde{V}}^{H} - \begin{bmatrix} \mathbf{I}_{q} \\ \mathbf{0} \end{bmatrix} \\ - \mathbf{\tilde{C}} &= \begin{bmatrix} \mathbf{I}_{p} & \mathbf{0} \end{bmatrix} \mathbf{\tilde{U}}\tilde{\Sigma}^{1/2} \\ - \mathbf{\tilde{D}} &= \mathbf{y}_{0} - \end{aligned} - -where :math:`p` indicates the number of outputs and :math:`q` the number -of inputs, and - -.. math:: - - - \begin{aligned} - \mathbf{\tilde{x}}_{k+1} &= \mathbf{\tilde{A}\tilde{x}}_{k} + \mathbf{\tilde{B}u}_{k} \\ - \mathbf{y}_{k} &= \mathbf{\tilde{C}\tilde{x}}_{k} + \mathbf{\tilde{D}u}_{k}. \\ - \end{aligned} diff --git a/docs/classes/fourier.md b/docs/classes/fourier.md deleted file mode 100644 index 6abd53031..000000000 --- a/docs/classes/fourier.md +++ /dev/null @@ -1,384 +0,0 @@ - - - -## Fourier series expansion - -Certain **periodic** functions can be expressed as a Fourier series expansion. With $P$ as the period of the function, - -The sine-cosine form represents $f(x)$ in the orthogonal basis $\left\{\cos{\frac{2\pi n x}{P}}, \sin{\frac{2\pi n x}{P}} ~|~ n \in \mathbb{Z}_{+} \right\}$. -$$\begin{aligned} - f(x) &= \frac{a_{0}}{2} - + \sum_{n=1}^{\infty}\left({ - a_{n}\cos{\frac{2{\pi}nx}{P}} - + b_{n}\sin{\frac{2{\pi}nx}{P}} - }\right) \\ - a_{0} &= \frac{2}{P} \int_{-\frac{P}{2}}^{\frac{P}{2}} {f(t)dt}~, \\ - a_{n} &= \frac{2}{P} \int_{-\frac{P}{2}}^{\frac{P}{2}} - {f(t) \cos{\frac{2{\pi}nt}{P}} dt}~, \text{ and} \\ - b_{n} &= \frac{2}{P} \int_{-\frac{P}{2}}^{\frac{P}{2}} - {f(t) \sin{\frac{2{\pi}nt}{P}} dt}~. -\end{aligned}$$ - -The exponential form represents $f(x)$ in the orthogonal basis $\left\{\exp{i \frac{2\pi n x}{P}} ~|~ n \in \mathbb{Z} \right\}$. -$$\begin{aligned} -f(x) &= \sum_{n=-\infty}^{n=\infty} {c_{n}e^{i\frac{2{\pi}n}{P}x}} \\ -c_{n} &= \frac{1}{P} \int_{-\frac{P}{2}}^{\frac{P}{2}} {f(t) e^{-i\frac{2{\pi}n}{P}t} dt} -\end{aligned}$$ - - -## Fourier transform - -Plug $c_{n}$ into the exponential form: -$$\begin{aligned} - f(x) &= \sum_{n=-\infty}^{n=\infty} {\frac{1}{P}} - \int_{-\frac{P}{2}}^{\frac{P}{2}} - f(t) e^{-i\frac{2{\pi}n}{P}t} dt~ - e^{i\frac{2{\pi}n}{P}x} \\ - &= \sum_{n=-\infty}^{n=\infty} {\frac{1}{P}} - \int_{-\frac{P}{2}}^{\frac{P}{2}} - f(t) e^{i\frac{2{\pi}n}{P}(x-t)} dt \\ -\end{aligned}$$ - -To generalize to certain **non-periodic** functions, take the limit as $P \rightarrow \infty$, and change from the variable of summation $n$ to the variable of integration $\omega = \frac{2\pi{}n}{P}$ . - -$$\begin{aligned} - f(x) &= \underset{P \rightarrow \infty}\lim{} - \sum_{n=-\infty}^{n=\infty} \left( {\frac{1}{P}} - \int_{-\frac{P}{2}}^{\frac{P}{2}} - f(t) e^{i\frac{2{\pi}n}{P}(x-t)} dt \right) \\ - &= \underset{P \rightarrow \infty}\lim{} - \int_{n=-\infty}^{n=\infty} \left( - \int_{-\frac{P}{2}}^{\frac{P}{2}} - f(t) e^{i\frac{2{\pi}n}{P}(x-t)} dt \right) {\frac{dn}{P}} \\ - &= \int_{\omega=-\infty}^{\omega=\infty} \left( - \int_{-\infty}^{\infty} - f(t) e^{i\omega(x-t)} dt \right) \frac{d\omega}{2\pi} \\ - &= \frac{1}{2\pi}\int_{-\infty}^{\infty} e^{i\omega x} \left( - \int_{-\infty}^{\infty} - f(t) e^{-i\omega t} dt \right) d\omega \\ -\end{aligned}$$ - -The Fourier transform $\mathcal{F}(f)$ is expressed as - -$$\begin{aligned} -\mathcal{F}(f) &= \phi(\omega) \\ -\phi(\omega) &= \int_{-\infty}^{\infty}{f(t)e^{-i\omega t} dt} -\end{aligned}$$ - -and exists for all $f$ that are piecewise continuous and absolutely integrable. - -The inverse Fourier transform $\mathcal{F}^{-1}(\phi)$ is expressed as - -$$\begin{aligned} -\mathcal{F}^{-1}(\phi) &= f(t) \\ -f(t) &= \frac{1}{2\pi}\int_{-\infty}^{\infty}{\phi(\omega)e^{i\omega x} d\omega} -\end{aligned}$$ - - -## Discrete Fourier Transform (DFT) - -The Discrete Fourier Transform fits a summation approximation over $t \in \{0,1,2,\dots,n-1\}$ of the integral $\phi(\omega) = \int_{-\infty}^{\infty}{f(t)e^{-i\omega t} dt}$ over the discrete domain $\omega_{n} = e^{i\frac{2\pi}{n}}, n \in \{0,1,2,\dots,n-1\}$. - -$$ -\phi(n) = -\sum_{t=0}^{t=n-1}f(t)e^{i\frac{2\pi}{n}t}~, -$$ - -This corresponds to frequencies $\omega \in {0, \frac{2\pi}{n}, \frac{4\pi}{n}}$ - -$$ -\phi(\omega_{n}) = -\sum_{t=0}^{t=n-1}f(t)\omega_{n}^{t}~, -$$ - -This transformation can be expressed as a Hermitian matrix multplication. - - -$$ - \begin{bmatrix} - c_{0} \\ - c_{1} \\ - c_{2} \\ - \vdots \\ - c_{n-1} \\ - \end{bmatrix} - = - \begin{bmatrix} - 1 & 1 & 1 & \cdots & 1 \\ - 1 & \omega_{n} & \omega_{n}^{2} & \cdots & \omega_{n}^{n-1} \\ - 1 & \omega_{n}^{2} & \omega_{n}^{4} & \cdots & \omega_{n}^{2(n-1)} \\ - \vdots & \vdots & \vdots & \ddots & \vdots \\ - 1 & \omega_{n}^{n-1} & \omega_{n}^{2(n-1)} & \cdots & \omega_{n}^{(n-1)^{2}} \\ - \end{bmatrix} - \begin{bmatrix} - y_{0} \\ - y_{1} \\ - y_{2} \\ - \vdots \\ - y_{n-1} \\ - \end{bmatrix} -$$ - - - -## Fast Fourier Transform (FFT) - -The bare-bones (radix-2) Cooley-Tukey algorithm requires the series length to be a power of 2. - -Base case: $n=1$ - -$$ - c_{0} - = F_{1}y_{0} - = \begin{bmatrix} - 1 - \end{bmatrix} - y_{0} -$$ - -The next case: $n=2$ -$$ -\begin{aligned} - \begin{bmatrix} - c_{0} \\ - c_{1} - \end{bmatrix} - &= F_{2} - \begin{bmatrix} - y_{0} \\ - y_{1} - \end{bmatrix} - = \begin{bmatrix} - 1 & 1 \\ - 1 & \omega_{2} - \end{bmatrix} - \begin{bmatrix} - y_{0} \\ - y_{1} - \end{bmatrix} - = \begin{bmatrix} - 1 & 1 \\ - 1 & -1 - \end{bmatrix} - \begin{bmatrix} - y_{0} \\ - y_{1} - \end{bmatrix} \\ - &= \begin{bmatrix} - I_{1} & D_{1} \\ - I_{1} & -D_{1} - \end{bmatrix} - \begin{bmatrix} - F_{1} & 0 \\ - 0 & F_{1} - \end{bmatrix} - \begin{bmatrix} - y_{0} \\ - y_{1} - \end{bmatrix} - = \begin{bmatrix} - 1 & 1 \\ - 1 & -1 - \end{bmatrix} - \begin{bmatrix} - 1 & 0 \\ - 0 & 1 - \end{bmatrix} - \begin{bmatrix} - y_{0} \\ - y_{1} - \end{bmatrix} -\end{aligned} -$$ - -For $n = 4, 8, 16, 32, ...$ - -$$ - \begin{bmatrix} - c_{0} \\ - c_{1} \\ - \vdots \\ - c_{n-1} - \end{bmatrix} - = F_{n} - \begin{bmatrix} - y_{0} \\ - y_{1} \\ - \vdots \\ - y_{n-1} - \end{bmatrix} - = \begin{bmatrix} - I_{n/2} & D_{n/2} \\ - I_{n/2} & -D_{n/2} - \end{bmatrix} - \begin{bmatrix} - F_{n/2} & 0 \\ - 0 & F_{n/2} - \end{bmatrix} - \begin{bmatrix} - y_{even} \\ - y_{odd} - \end{bmatrix} -$$ - - - - - - - diff --git a/docs/classes/index.rst b/docs/classes/index.rst deleted file mode 100644 index d32bb7479..000000000 --- a/docs/classes/index.rst +++ /dev/null @@ -1,13 +0,0 @@ -Classes -======= - -.. toctree:: - :maxdepth: 2 - - statespace.rst - modes.rst - era.rst - okid.rst - srim.rst - - diff --git a/docs/classes/modes.md b/docs/classes/modes.md deleted file mode 100644 index 1e1d94aa6..000000000 --- a/docs/classes/modes.md +++ /dev/null @@ -1,88 +0,0 @@ -# Modal Properties from State Space Realization - -Structural system dynamics are in continuous time. With the time-domain system identification methods in the `mdof` package such as OKID-ERA and SRIM, a structural system's discrete-time state space realization is obtained from measured data. The following process recovers the structure's modal properties (i.e., natural frequencies, modal damping ratios, and mode shapes) from the discrete state space realization coefficients, $\mathbf{A}$ and $\mathbf{C}$. - -## Eigendecompositions of $\mathbf{A}_{c}$ and $\mathbf{A}$ - -The relationship between the eigendecompositions of $\mathbf{A}_{c}$, the continuous-time state transition matrix, and $\mathbf{A}$, the discrete-time state transition matrix, is shown below. - -$$\begin{aligned} -\mathbf{A}_{c} &= \Phi\Lambda\Phi^{-1} \\ -\mathbf{A} &= \Psi\Gamma\Psi^{-1} \\ -\end{aligned}$$ - -$$\begin{aligned} -\mathbf{A} = e^{\mathbf{A}_{c}\Delta t} = e^{\Phi\Lambda\Phi^{-1}\Delta t} = \Phi e^{\Lambda\Delta t}\Phi^{-1} -\end{aligned}$$ - -$$\begin{aligned} -\Psi = \Phi, \quad \Gamma = e^{\Lambda\Delta t} \\ -\end{aligned}$$ - -where - -$$\Psi = -\begin{bmatrix} -\psi_{1} & \psi_{2} & \cdots & \psi_{r} -\end{bmatrix}, \quad{} -\Phi = -\begin{bmatrix} -\phi_{1} & \phi_{2} & \cdots & \phi_{r} -\end{bmatrix} -$$ - -$$ -\Gamma = -\begin{bmatrix} -\gamma_{1} & 0 & \cdots & 0 \\ -0 & \gamma_{2} & \cdots & 0 \\ -\vdots & \vdots & \ddots & \vdots \\ -0 & 0 & \cdots & \gamma_{r} -\end{bmatrix}, \quad -\Lambda = \begin{bmatrix} -\lambda_{1} & 0 & \cdots & 0 \\ -0 & \lambda_{2} & \cdots & 0 \\ -\vdots & \vdots & \ddots & \vdots \\ -0 & 0 & \cdots & \lambda_{r} -\end{bmatrix}$$ - -- $\gamma_{j}$ and $\psi_{j}$ are the eigenvalues and eigenvectors of $\mathbf{A}$, - -- $\lambda_{j}$ and $\phi_{j}$ are the eigenvalues and eigenvectors of $\mathbf{A}_{c}$, - -- $j \in [1,2,\dots,r]$, and $r$ is the model order. - - -## Natural Frequencies and Modal Damping Ratios - -From the eigendecompositions of $\mathbf{A}_{c}$ and $\mathbf{A}$, we have - -$$\begin{aligned} -\Gamma = e^{\Lambda\Delta t} -\implies -\lambda_{j} = (\ln{\gamma_{j}})/\Delta t~. -\end{aligned}$$ - -Assuming modal damping, these eigenvalues $\lambda_{j}$ contain the damped natural frequencies and modal damping factors: - -$$\begin{aligned} -\lambda_{j} &= -\zeta_{j}\omega_{j} \pm i\left(\omega_{j}\sqrt{1-\zeta_{j}^{2}}\right), \hspace{0.5cm} i=\sqrt{-1} \\ \\ -\lambda_{j}\bar{\lambda}_{j} &= \zeta_{j}^{2}\omega_{j}^{2} + \omega_{j}^{2}\left(1-\zeta_{j}^{2}\right) \\ - &= \omega_{j}^{2}\left(\zeta_{j}^{2} + 1 - \zeta_{j}^{2}\right) \\ - &= \omega_{j}^{2} \\ \\ -\omega_{j} &= \sqrt{\lambda_{j}\bar{\lambda}_{j}} = | \lambda_{j} |, \\ -\zeta_{j} &= -\frac{\text{Re}(\lambda_{j})}{\omega_{j}},\\ -\end{aligned}$$ - -where the overline symbol $~\overline{\cdot}~$ indicates the complex conjugate, $\text{Re}(\cdot)$ indicates the real part of the complex number, $\omega_{j}$ are the modal frequencies, and $\zeta_{j}$ are the modal damping ratios. - -## Mode Shapes - -The eigenvectors of the continuous state transition matrix $\mathbf{A}_{c}$ transform the modal coordinates into the state coordinates. We define the mode shapes as the transformation from modal coordinates to the output coordinates. - -$$\begin{aligned} -\mathbf{x}(t) &= \mathbf{\Phi q}(t) \\ -\mathbf{y}(t) &= \mathbf{Cx}(t) + \mathbf{Du}(t) \\ -\mathbf{\Phi}_{\text{modal}} &= \mathbf{C\Phi} -\end{aligned}$$ - diff --git a/docs/classes/modes.rst b/docs/classes/modes.rst deleted file mode 100644 index 01891b35b..000000000 --- a/docs/classes/modes.rst +++ /dev/null @@ -1,124 +0,0 @@ -Modal Properties from State Space Realization -============================================= - -Structural system dynamics are in continuous time. With the time-domain -system identification methods in the ``mdof`` package such as OKID-ERA -and SRIM, a structural system’s discrete-time state space realization is -obtained from measured data. The following process recovers the -structure’s modal properties (i.e., natural frequencies, modal damping -ratios, and mode shapes) from the discrete state space realization -coefficients, :math:`\mathbf{A}` and :math:`\mathbf{C}`. - -Eigendecompositions of :math:`\mathbf{A}_{c}` and :math:`\mathbf{A}` --------------------------------------------------------------------- - -The relationship between the eigendecompositions of -:math:`\mathbf{A}_{c}`, the continuous-time state transition matrix, and -:math:`\mathbf{A}`, the discrete-time state transition matrix, is shown -below. - -.. math:: - - \begin{aligned} - \mathbf{A}_{c} &= \Phi\Lambda\Phi^{-1} \\ - \mathbf{A} &= \Psi\Gamma\Psi^{-1} \\ - \end{aligned} - -.. math:: - - \begin{aligned} - \mathbf{A} = e^{\mathbf{A}_{c}\Delta t} = e^{\Phi\Lambda\Phi^{-1}\Delta t} = \Phi e^{\Lambda\Delta t}\Phi^{-1} - \end{aligned} - -.. math:: - - \begin{aligned} - \Psi = \Phi, \quad \Gamma = e^{\Lambda\Delta t} \\ - \end{aligned} - -where - -.. math:: - - \Psi = - \begin{bmatrix} - \psi_{1} & \psi_{2} & \cdots & \psi_{r} - \end{bmatrix}, \quad{} - \Phi = - \begin{bmatrix} - \phi_{1} & \phi_{2} & \cdots & \phi_{r} - \end{bmatrix} - -.. math:: - - - \Gamma = - \begin{bmatrix} - \gamma_{1} & 0 & \cdots & 0 \\ - 0 & \gamma_{2} & \cdots & 0 \\ - \vdots & \vdots & \ddots & \vdots \\ - 0 & 0 & \cdots & \gamma_{r} - \end{bmatrix}, \quad - \Lambda = \begin{bmatrix} - \lambda_{1} & 0 & \cdots & 0 \\ - 0 & \lambda_{2} & \cdots & 0 \\ - \vdots & \vdots & \ddots & \vdots \\ - 0 & 0 & \cdots & \lambda_{r} - \end{bmatrix} - -- :math:`\gamma_{j}` and :math:`\psi_{j}` are the eigenvalues and - eigenvectors of :math:`\mathbf{A}`, - -- :math:`\lambda_{j}` and :math:`\phi_{j}` are the eigenvalues and - eigenvectors of :math:`\mathbf{A}_{c}`, - -- :math:`j \in [1,2,\dots,r]`, and :math:`r` is the model order. - -Natural Frequencies and Modal Damping Ratios --------------------------------------------- - -From the eigendecompositions of :math:`\mathbf{A}_{c}` and -:math:`\mathbf{A}`, we have - -.. math:: - - \begin{aligned} - \Gamma = e^{\Lambda\Delta t} - \implies - \lambda_{j} = (\ln{\gamma_{j}})/\Delta t~. - \end{aligned} - -Assuming modal damping, these eigenvalues :math:`\lambda_{j}` contain -the damped natural frequencies and modal damping factors: - -.. math:: - - \begin{aligned} - \lambda_{j} &= -\zeta_{j}\omega_{j} \pm i\left(\omega_{j}\sqrt{1-\zeta_{j}^{2}}\right), \hspace{0.5cm} i=\sqrt{-1} \\ \\ - \lambda_{j}\bar{\lambda}_{j} &= \zeta_{j}^{2}\omega_{j}^{2} + \omega_{j}^{2}\left(1-\zeta_{j}^{2}\right) \\ - &= \omega_{j}^{2}\left(\zeta_{j}^{2} + 1 - \zeta_{j}^{2}\right) \\ - &= \omega_{j}^{2} \\ \\ - \omega_{j} &= \sqrt{\lambda_{j}\bar{\lambda}_{j}} = | \lambda_{j} |, \\ - \zeta_{j} &= -\frac{\text{Re}(\lambda_{j})}{\omega_{j}},\\ - \end{aligned} - -where the overline symbol :math:`~\overline{\cdot}~` indicates the -complex conjugate, :math:`\text{Re}(\cdot)` indicates the real part of -the complex number, :math:`\omega_{j}` are the modal frequencies, and -:math:`\zeta_{j}` are the modal damping ratios. - -Mode Shapes ------------ - -The eigenvectors of the continuous state transition matrix -:math:`\mathbf{A}_{c}` transform the modal coordinates into the state -coordinates. We define the mode shapes as the transformation from modal -coordinates to the output coordinates. - -.. math:: - - \begin{aligned} - \mathbf{x}(t) &= \mathbf{\Phi q}(t) \\ - \mathbf{y}(t) &= \mathbf{Cx}(t) + \mathbf{Du}(t) \\ - \mathbf{\Phi}_{\text{modal}} &= \mathbf{C\Phi} - \end{aligned} diff --git a/docs/classes/okid.rst b/docs/classes/okid.rst deleted file mode 100644 index 00d35ffc8..000000000 --- a/docs/classes/okid.rst +++ /dev/null @@ -1,111 +0,0 @@ - -Observer Kalman Filter Identification (OKID) ---------------------------------------------- - -Structural dynamics are noisy, hard to measure, and lightly damped, and -ERA is intended only to characterize impulse responses rather than time -histories. However, available data from ambient or small excitations -during structure service can be de-noised and used to estimate impulse -response data. Then, ERA can be used to obtain a reduced order model -even if the available data are not a clean impulse response. This -process is called `Observer Kalman -Identification `__, or OKID-ERA when -combined with ERA. - -When noise is incorporated into the discrete LTI state-space -representation of a structural system, it becomes a *linear Gaussian -model* of a *hidden Markov process*. - -Because the data are assumed to follow a linear Gaussian model, Kalman -filtering can estimate an impulse response that is most consistent with -the input-output data. The estimated model after filtering is the same -as that of ERA: - -.. math:: - - - \begin{aligned} - \mathbf{x}_{k+1} &= \mathbf{Ax}_{k} + \mathbf{Bu}_{k} \\ - \mathbf{y}_{k} &= \mathbf{Cx}_{k} + \mathbf{Du}_{k} \\ - \end{aligned} - -Since the input is no longer an impulse, the state-space evolution -includes more terms than ERA. - -.. math:: - - - \begin{aligned} - \mathbf{u}_{0},\mathbf{u}_{1},\mathbf{u}_{2},...,\mathbf{u}_{k} :=& \text{given input} \\ - \mathbf{x}_{0},\mathbf{x}_{1},\mathbf{x}_{2},...,\mathbf{x}_{k} =& \mathbf{0},(\mathbf{Bu}_{0}),(\mathbf{ABu}_{0}+\mathbf{Bu}_{1}),...,(\mathbf{A}^{k-1}\mathbf{Bu}_{0}+\mathbf{A}^{k-2}\mathbf{Bu}_{1}+...+\mathbf{Bu}_{k-1}) \\ - \mathbf{y}_{0},\mathbf{y}_{1},\mathbf{y}_{2},...,\mathbf{y}_{k} =& \mathbf{Du}_0,(\mathbf{CBu}_{0}+\mathbf{Du}_{1}),(\mathbf{CABu}_{0}+\mathbf{CBu}_{1}+\mathbf{Du}_{2}),..., \\ - & (\mathbf{CA}^{k-1}\mathbf{Bu}_{0}+\mathbf{CA}^{k-2}\mathbf{Bu}_{1}+...+\mathbf{Du}_{k}). - \end{aligned} - -The output data can be expressed in terms of the Markov parameters and -an upper triangular *data matrix* :math:`\mathscr{B}` built from the -input data; however, inverting :math:`\mathscr{B}` is often -computationally expensive or ill-conditioned. - -.. math:: - - - \underbrace{\begin{bmatrix} \mathbf{y}_{0} & \mathbf{y}_{1} & \mathbf{y}_{2} & \cdots & \mathbf{y}_{m} \end{bmatrix}}_{\mathbf{S}} - = - \underbrace{\begin{bmatrix} \mathbf{y}_{0} & \mathbf{y}_{1} & \mathbf{y}_{2} & \cdots & \mathbf{y}_{m} \end{bmatrix}_{\delta}}_{\mathbf{S}_{\delta}} - \underbrace{\begin{bmatrix} - \mathbf{u}_{0} & \mathbf{u}_{1} & \cdots & \mathbf{u}_{m} \\ - \mathbf{0} & \mathbf{u}_{0} & \cdots & \mathbf{u}_{m-1} \\ - \vdots & \vdots & \ddots & \vdots \\ - \mathbf{0} & \mathbf{0} & \cdots & \mathbf{u}_{0} \\ - \end{bmatrix}}_{\mathscr{B}} - -where the subscript :math:`\delta` indicates that the response comes -from an impulse input. - -The Kalman filter is applied by augmenting the system with the outputs -:math:`\mathbf{y}_{i}` to form the *augmented data matrix* -:math:`\mathscr{V}`: - -.. math:: - - - \mathscr{V} - = - \begin{bmatrix} - \mathbf{u}_{0} & \mathbf{u}_{1} & \cdots & \mathbf{u}_{l} & \cdots & \mathbf{u}_{m} \\ - \mathbf{0} & \mathbf{v}_{0} & \cdots & \mathbf{v}_{l-1} & \cdots & \mathbf{v}_{m-1} \\ - \vdots & \vdots & \ddots & \vdots & \ddots & \vdots \\ - \mathbf{0} & \mathbf{0} & \cdots & \mathbf{v}_{0} & \cdots & \mathbf{v}_{m-l} \\ - \end{bmatrix}, \hspace{1cm} - \mathbf{v}_{i} = \begin{bmatrix} \mathbf{u}_{i} \\ \mathbf{y}_{i} \end{bmatrix}. - -Then, the Markov parameters (i.e., the impulse response) can be -estimated as a function of the input and output data as follows: - -.. math:: \hat{\mathbf{S}}_\delta = \mathbf{S}\mathscr{V}^{\dagger} - -where the superscript :math:`\dagger` indicates pseudo-inverse, - -Extract the estimated, or *observer*, Markov parameters from the block -columns of :math:`\hat{\mathbf{S}}_\delta`: - -.. math:: - - - \begin{aligned} - \hat{\mathbf{S}}_{\delta 0} &\in \mathbb{R}^{p\times q} \\ - \hat{\mathbf{S}}_{\delta k} &= - \begin{bmatrix} \hat{\mathbf{S}}_{\delta k}^{(1)} & \hat{\mathbf{S}}_{\delta k}^{(2)} \end{bmatrix} , \hspace{0.5cm} k\in[1,2,...] \\ - \hat{\mathbf{S}}_{\delta k}^{(1)} &\in\mathbb{R}^{p\times q}, \hspace{0.5cm} - \hat{\mathbf{S}}_{\delta k}^{(2)} \in\mathbb{R}^{p\times p} - \end{aligned} - -Reconstruct the system Markov parameters: - -.. math:: - - - \mathbf{y}_{\delta 0} = \hat{\mathbf{S}}_{\delta 0} = \mathbf{D}, \hspace{0.5cm} - \mathbf{y}_{\delta k} = \hat{\mathbf{S}}_{\delta k}^{(1)} - + \sum_{i=1}^{k}{\hat{\mathbf{S}}_{\delta k}^{(2)}}\mathbf{y}_{\delta (k-i)}. diff --git a/docs/classes/psd.md b/docs/classes/psd.md deleted file mode 100644 index 663473938..000000000 --- a/docs/classes/psd.md +++ /dev/null @@ -1,100 +0,0 @@ - - -# Power Spectral Density - -The power spectral density is the norm of the Fourier transform. In a way, it measures energy content at each frequency of response. - -It is also the averaged(?) Fourier transform of the autocorrelation. - -## Norm of Fourier Transform - -The power spectral density, $P(f)$, of a signal $y(t)$, is - -$$\begin{aligned} -P(f) &= \frac{1}{T} -\int_{-T}^{T}{ - \left| Y_{T}(f) \right|^{2}df -} -\end{aligned}$$ - -Where the Fourier transform $Y_{T}(f)$ of $y_{T}(t)$ is defined as: -$$\begin{aligned} -Y_{T}(f) &= \mathcal{F}\left\{y_{T}(t)\right\} = \int_{-T}^{T}{ - e^{-i2\pi ft}y_{T}(t)dt -} -\end{aligned}$$ - -where $f \in \mathbb{R}_{+}$ is frequency in Hertz, and $t \in \mathbb{R}_{+}$ is time in seconds. - -**question.** -$$\begin{aligned} -y_{T}(t) = \begin{bmatrix} -y(t) \\ -y(t+\Delta t) \\ -y(t+2\Delta t) \\ -\vdots \\ -y(t+(T-1)\Delta t) \\ -\end{bmatrix} -\in \mathbb{R}^{T} ~~\textbf{??} -\end{aligned}$$ - -## Fourier Transform of Autocorrelation - -Autocorrelation, discrete: - -$$\begin{aligned} -R_{yy} &\approx \frac{1}{N}\mathbf{Y}_{p}(k)\mathbf{Y}_{p}^{*}(k), ~~N \gg 0 -\end{aligned}$$ - -$$ -\begin{aligned} -\mathbf{Y}_{p}(k) &= \begin{bmatrix} \mathbf{y}_{p}(k) & \mathbf{y}_{p}(k+1) & \cdots & \mathbf{y}_{p}(k+N-1) \end{bmatrix} \\ -&= \begin{bmatrix} -\mathbf{y}(k) & \mathbf{y}(k+1) & \cdots & \mathbf{y}(k+N-1)\\ -\mathbf{y}(k+1) & \mathbf{y}(k+2) & \cdots & \mathbf{y}(k+N) \\ -\vdots & \vdots & \ddots & \vdots \\ -\mathbf{y}(k+p-1) & \mathbf{y}(k+p) & \cdots & \mathbf{y}(k+N+p-2) -\end{bmatrix} -\end{aligned} -$$ - -```{=tex} -\pagebreak -``` - -Autocorrelation, continuous: - -$$\begin{aligned} -R_{yy}(\tau) -&= \mathbb{E}\left\{\mathbf{y}(t+\tau)\mathbf{y}^{*}(t)\right\} \\ -&= \int_{-\infty}^{\infty}{\mathbf{y}(t+\tau)\mathbf{y}^{*}(t)f_{y}(t)dt} -\\ -&= \int_{-\infty}^{\infty}{\mathbf{y}(t+\tau)\mathbf{y}^{*}(t)dt} ~~\textbf{??} -\\ -&= \underset{T \rightarrow \infty}{\lim}\frac{1}{T}\int_{-\infty}^{\infty}{\mathbf{y}_{T}^{*}(t-\tau)\mathbf{y}_{T}(t)dt} ~~\textbf{??} -\end{aligned}$$ - -Fourier transform of this autocorrelation: -$$\begin{aligned} -\mathcal{F}\left\{R_{yy}(\tau)\right\} &= \mathbf{R}_{yy}(f) \\ -&= \int_{-\infty}^{\infty}e^{-i2\pi ft}R_{yy}(\tau)d\tau \\ -&= \int_{-\infty}^{\infty}e^{-i2\pi ft}\left[ - \underset{T \rightarrow \infty}{\lim}\frac{1}{T}\int_{-\infty}^{\infty}{\mathbf{y}_{T}^{*}(t-\tau)\mathbf{y}_{T}(t)dt} -\right](\tau)d\tau \\ -&= \underset{N \rightarrow \infty}{\lim}\frac{\left( \Delta t \right)^{2}}{T}\left| \sum_{n=-N}^{N} y_{n} e^{-i2\pi fn\Delta t} \right|^{2} \\ -&\approx \underset{T \rightarrow \infty}{\lim}{\frac{1}{T}\left| Y_{T}(f) \right|^{2}} -\end{aligned}$$ - -Averaged...**??** Fourier transform of this autocorrelation: -$$\begin{aligned} -P &= \int_{-\infty}^{\infty}{\mathbf{R}_{yy}(f) df} -&= \int_{-\infty}^{\infty}{ - \underset{T \rightarrow \infty}{\lim} - \frac{1}{T} - \left| Y_{T}(f) \right|^{2}df - } -\end{aligned}$$ \ No newline at end of file diff --git a/docs/classes/srim.md b/docs/classes/srim.md deleted file mode 100644 index bc553e058..000000000 --- a/docs/classes/srim.md +++ /dev/null @@ -1,299 +0,0 @@ -# System Realization by Information Matrix (SRIM) - -For discrete-time systems, the correlation between inputs, outputs, and state yield information about the system's state evolution and response. In fact, the state equations can be estimated by manipulating correlation matrices through the method, [System Realization by Information Matrix](https://doi.org/10.2514/2.4068) (SRIM). - -## Discrete-Time System Matrix Equation -We begin with discrete-time state equations that correspond to the structure's dynamics (see [Discrete LTI State-Space Representation](https://chrysatlchern.github.io/mdof/theory/statespace.html#discrete-lti-state-space-representation)). - -$$ -\begin{aligned} -\mathbf{x}(k+1) &= \mathbf{A}\mathbf{x}(k) + \mathbf{B}\mathbf{u}(k) \\ -\mathbf{y}(k) &= \mathbf{C}\mathbf{x}(k) + \mathbf{D}\mathbf{u}(k) -\end{aligned} -$$ - -By noting the state evolution - -$$ -\begin{aligned} -\mathbf{x}(k+1) &= \mathbf{Ax}(k)+\mathbf{B}\mathbf{u}(k)\\ -\mathbf{x}(k+2) &= \mathbf{A}^2\mathbf{x}(k) + \mathbf{AB}\mathbf{u}(k) + \mathbf{B}\mathbf{u}(k+1)\\ -\mathbf{x}(k+3) &= \mathbf{A}^{3}\mathbf{x}(k) + \mathbf{A}^{2}\mathbf{Bu}(k) + \mathbf{A}\mathbf{B}\mathbf{u}(k+1) + \mathbf{B}\mathbf{u}(k+2), -\end{aligned} -$$ - -we can generalize the response for the timepoint $k+p-1$: - -$$ -\begin{aligned} -\mathbf{x}(k+p) &= \mathbf{A}^{p}\mathbf{x}(k) + \sum_{i=1}^{p}\mathbf{A}^{p-i}\mathbf{Bu}(k+i-1) -\\ -\mathbf{x}(k+p-1) &= \mathbf{A}^{p-1}\mathbf{x}(k) + \sum_{i=1}^{p-1}\mathbf{A}^{p-i-1}\mathbf{Bu}(k+i-1) -\\ -\mathbf{y}(k+p-1) &= \mathbf{CA}^{p-1}\mathbf{x}(k) + \sum_{i=1}^{p-1}\mathbf{CA}^{p-i-1}\mathbf{Bu}(k+i-1)+\mathbf{Du}(k+p-1)~. -\end{aligned} -$$ - -Then, we can vertically stack $p$ successive time-points into a column vector and express this vector as $\mathbf{y}_{p}(k)$: - -$$ -\begin{aligned} -\mathbf{y}_{p}(k) &= \mathcal{O}_{p}\mathbf{x}(k) + \mathcal{T}_{p}\mathbf{u}_{p}(k) \\ -\begin{bmatrix} -\mathbf{y}(k) \\ -\mathbf{y}(k+1) \\ -\vdots \\ -\mathbf{y}(k+p-1) -\end{bmatrix} -=& -\begin{bmatrix} -\mathbf{C} \\ -\mathbf{CA} \\ -\mathbf{CA}^{2} \\ -\vdots \\ -\mathbf{CA}^{p-1} -\end{bmatrix} -\mathbf{x}(k) -~+ \\ -& -\begin{bmatrix} -\mathbf{D} \\ -\mathbf{CB} & \mathbf{D} \\ -\mathbf{CAB} & \mathbf{CB} & \mathbf{D} \\ -\vdots & \vdots & \vdots & \ddots \\ -\mathbf{CA}^{p-2}\mathbf{B} & \mathbf{CA}^{p-3}\mathbf{B} & \mathbf{CA}^{p-4}\mathbf{B} & \cdots & \mathbf{D} -\end{bmatrix} -\begin{bmatrix} -\mathbf{u}(k) \\ -\mathbf{u}(k+1) \\ -\vdots \\ -\mathbf{u}(k+p-1) -\end{bmatrix}~. -\end{aligned} -$$ - - - -Finally, we horizontally stack $N$ successive timepoints of these column vectors in a matrix, to get the matrix equation - -$$ -\boxed{\mathbf{Y}_{p}(k) = \mathcal{O}_{p}\mathbf{X}(k) + \mathcal{T}_{p}\mathbf{U}_{p}(k)} ~, -$$ - -where -$$ -\begin{aligned} -\mathbf{Y}_{p}(k) &= \begin{bmatrix} \mathbf{y}_{p}(k) & \mathbf{y}_{p}(k+1) & \cdots & \mathbf{y}_{p}(k+N-1) \end{bmatrix} \\ -&= \begin{bmatrix} -\mathbf{y}(k) & \mathbf{y}(k+1) & \cdots & \mathbf{y}(k+N-1)\\ -\mathbf{y}(k+1) & \mathbf{y}(k+2) & \cdots & \mathbf{y}(k+N) \\ -\vdots & \vdots & \ddots & \vdots \\ -\mathbf{y}(k+p-1) & \mathbf{y}(k+p) & \cdots & \mathbf{y}(k+N+p-2) -\end{bmatrix} -\end{aligned} -$$ - -$$ -\mathbf{X}(k) = \begin{bmatrix} \mathbf{x}(k) & \mathbf{x}(k+1) & \cdots & \mathbf{x}(k+N-1) \end{bmatrix} -$$ - -$$ -\begin{aligned} -\mathbf{U}_{p}(k) &= \begin{bmatrix} \mathbf{u}_{p}(k) & \mathbf{u}_{p}(k+1) & \cdots & \mathbf{u}_{p}(k+N-1) \end{bmatrix} \\ -&= \begin{bmatrix} -\mathbf{u}(k) & \mathbf{u}(k+1) & \cdots & \mathbf{u}(k+N-1)\\ -\mathbf{u}(k+1) & \mathbf{u}(k+2) & \cdots & \mathbf{u}(k+N) \\ -\vdots & \vdots & \ddots & \vdots \\ -\mathbf{u}(k+p-1) & \mathbf{u}(k+p) & \cdots & \mathbf{u}(k+N+p-2) -\end{bmatrix}~. -\end{aligned} -$$ - - -## Observability Matrix from Information Matrix -By post-multiplying the matrix equation by $\frac{1}{N}\mathbf{U}_{p}^{T}(k)$, $\frac{1}{N}\mathbf{Y}_{p}^{T}(k)$ or $\frac{1}{N}\mathbf{X}_{p}^{T}(k)$, we obtain relationships between correlation matrices $\mathbf{R}_{yy}$, $\mathbf{R}_{yu}$, $\mathbf{R}_{uu}$, and $\mathbf{R}_{xx}$ (See [Appendix](#appendix-manipulation-of-discrete-time-system-matrix-equation-into-correlation-matrix-relationships)). - -$$ -\mathbf{R}_{yy} - \mathbf{R}_{yu}\mathbf{R}_{uu}^{-1}\mathbf{R}_{yu}^{T} = \mathcal{O}_{p}\mathbf{R}_{xx}\mathcal{O}_{p}^{T} ~, -$$ - -where - -$$\begin{aligned} -\mathbf{R}_{yy} &= \frac{1}{N}\mathbf{Y}_{p}(k)\mathbf{Y}_{p}^{T}(k), \quad{} -\mathbf{R}_{yu} = \frac{1}{N}\mathbf{Y}_{p}(k)\mathbf{U}_{p}^{T}(k) \\ -\mathbf{R}_{uu} &= \frac{1}{N}\mathbf{U}_{p}(k)\mathbf{U}_{p}^{T}(k) , \quad{} -\mathbf{R}_{xx} = \frac{1}{N}\mathbf{X}(k)\mathbf{X}^{T}(k) ~. -\end{aligned}$$ - -The left side of the equation is found from input and output measurements, and is called the *information matrix* of the data. Its singular value decomposition is computed to yield the *observability matrix* $\mathcal{O}_{p}$. - -$$ -\mathbf{R}_{yy} - \mathbf{R}_{yu}\mathbf{R}_{uu}^{-1}\mathbf{R}_{yu}^{T} = \mathbf{U} \Sigma \mathbf{U}^{T} = \mathcal{O}_{p}\mathbf{R}_{xx}\mathcal{O}_{p}^{T} ~. -$$ - -## State Equation Matrices from Observability Matrix - -Now, the state equation matrices $\mathbf{A}$ and $\mathbf{C}$ can be obtained from the observability matrix $\mathcal{O}_p$. - -$$ -\begin{aligned} -\mathcal{O}_{p} -= -\begin{bmatrix} -\mathbf{C} \\ -\mathbf{CA} \\ -\mathbf{CA}^{2} \\ -\vdots \\ -\mathbf{CA}^{p-1} -\end{bmatrix} -, \quad{} -\mathcal{O}_{p}(:-1) -= -\begin{bmatrix} -\mathbf{C} \\ -\mathbf{CA} \\ -\mathbf{CA}^{2} \\ -\vdots \\ -\mathbf{CA}^{p-2} -\end{bmatrix} -, \quad{} -\mathcal{O}_{p}(1:) -= -\begin{bmatrix} -\mathbf{CA} \\ -\mathbf{CA}^{2} \\ -\mathbf{CA}^{3} \\ -\vdots \\ -\mathbf{CA}^{p-1} -\end{bmatrix} -\end{aligned} -$$ - -$$ -\mathbf{A} = \mathcal{O}_{p}(:-1)^{+}\mathcal{O}_{p}(1:) -$$ - -$$ -\mathbf{C} = \mathcal{O}_{p}(0) -$$ - -## Output Error Minimization - -$$ -\Phi = \begin{bmatrix} -\mathbf{C} & \mathcal{U}_{p}(0) & \mathbf{0}_{p\times r} \\ -\mathbf{CA} & \mathcal{U}_{p}(1) & \mathbf{C}\mathcal{U}_{r}(0) \\ -\mathbf{CA^{2}} & \mathcal{U}_{p}(2) & \mathbf{CA}\mathcal{U}_{r}(0) + \mathbf{C}\mathcal{U}_{r}(1) \\ -\vdots & \vdots & \vdots \\ -\mathbf{CA}^{ns-1} & \mathcal{U}_{p}(ns-1) & \sum_{k=0}^{ns-2}\mathbf{CA}^{ns-k-2}\mathcal{U}_{r}(k) -\end{bmatrix} \in \mathbb{R}^{(ns*p) \times (pr+pq+rq)} -$$ - - - -## Appendix: Manipulation of discrete-time system matrix equation into correlation matrix relationships - -In ([Juang 1997](https://doi.org/10.2514/2.4068)), the discrete-time system matrix equation is manipulated into a form describing the relationship between correlation matrices $\mathbf{R}_{yy}$, $\mathbf{R}_{yu}$, $\mathbf{R}_{uu}$, and $\mathbf{R}_{xx}$. - -Post-multiplying the [discrete-time system matrix equation](#discrete-time-system-matrix-equation) by $\frac{1}{N}\mathbf{U}_{p}^{T}(k)$: - -$$\begin{aligned} -\frac{1}{N}\mathbf{Y}_{p}(k)\mathbf{U}_{p}^{T}(k) &= \mathcal{O}_{p}\frac{1}{N}\mathbf{X}(k)\mathbf{U}_{p}^{T}(k) + \mathcal{T}_{p}\frac{1}{N}\mathbf{U}_{p}(k)\mathbf{U}_{p}^{T}(k) -\\ -\mathbf{R}_{yu} &= \mathcal{O}_{p}\mathbf{R}_{xu} + \mathcal{T}_{p}\mathbf{R}_{uu} -\\ -\mathcal{T}_{p} &= \left( \mathbf{R}_{yu} - \mathcal{O}_{p}\mathbf{R}_{xu} \right)\mathbf{R}_{uu}^{-1} -\end{aligned}$$ - -Post-multiplying by $\frac{1}{N}\mathbf{Y}_{p}^{T}(k)$: - -$$\begin{aligned} -\frac{1}{N}\mathbf{Y}_{p}(k)\mathbf{Y}_{p}^{T}(k) &= \mathcal{O}_{p}\frac{1}{N}\mathbf{X}(k)\mathbf{Y}_{p}^{T}(k) + \mathcal{T}_{p}\frac{1}{N}\mathbf{U}_{p}(k)\mathbf{Y}_{p}^{T}(k) -\\ -\mathbf{R}_{yy} &= \mathcal{O}_{p}\mathbf{R}_{yx}^{T} + \mathcal{T}_{p}\mathbf{R}_{yu}^{T} -\\ -\mathbf{R}_{yy} &= \mathcal{O}_{p}\mathbf{R}_{yx}^{T} + \left( \mathbf{R}_{yu} - \mathcal{O}_{p}\mathbf{R}_{xu} \right)\mathbf{R}_{uu}^{-1}\mathbf{R}_{yu}^{T} -\end{aligned}$$ - -Post-multiplying by $\frac{1}{N}\mathbf{X}_{p}^{T}(k)$: - -$$\begin{aligned} -\frac{1}{N}\mathbf{Y}_{p}(k)\mathbf{X}_{p}^{T}(k) &= \mathcal{O}_{p}\frac{1}{N}\mathbf{X}(k)\mathbf{X}_{p}^{T}(k) + \mathcal{T}_{p}\frac{1}{N}\mathbf{U}_{p}(k)\mathbf{X}_{p}^{T}(k) -\\ -\mathbf{R}_{yx} &= \mathcal{O}_{p}\mathbf{R}_{xx} + \mathcal{T}_{p}\mathbf{R}_{xu}^{T} -\\ -\mathbf{R}_{yx} &= \mathcal{O}_{p}\mathbf{R}_{xx} + \left( \mathbf{R}_{yu} - \mathcal{O}_{p}\mathbf{R}_{xu} \right)\mathbf{R}_{uu}^{-1}\mathbf{R}_{xu}^{T} -\end{aligned}$$ - -Substituting the equation for $\mathbf{R}_{yx}$ into the equation for $\mathbf{R}_{yy}$: - -$$\begin{aligned} -\mathbf{R}_{yy} =& ~\mathcal{O}_{p} -\left(\mathcal{O}_{p}\mathbf{R}_{xx} + \left( \mathbf{R}_{yu} - \mathcal{O}_{p}\mathbf{R}_{xu} \right)\mathbf{R}_{uu}^{-1}\mathbf{R}_{xu}^{T}\right)^{T} -\\ -&+ -\left( \mathbf{R}_{yu} - \mathcal{O}_{p}\mathbf{R}_{xu} \right)\mathbf{R}_{uu}^{-1}\mathbf{R}_{yu}^{T} -\\ -=& ~\mathcal{O}_{p}\mathbf{R}_{xx}\mathcal{O}_{p}^{T} - + \mathcal{O}_{p}\mathbf{R}_{xu}\mathbf{R}_{uu}^{-1} \left( \mathbf{R}_{yu}^{T} - \mathbf{R}_{xu}^{T}\mathcal{O}_{p}^{T} \right) -\\ -&+ -\left( \mathbf{R}_{yu} - \mathcal{O}_{p}\mathbf{R}_{xu} \right)\mathbf{R}_{uu}^{-1}\mathbf{R}_{yu}^{T} -\\ -=& ~\mathcal{O}_{p}\mathbf{R}_{xx}\mathcal{O}_{p}^{T} - + \mathcal{O}_{p}\mathbf{R}_{xu}\mathbf{R}_{uu}^{-1} \mathbf{R}_{yu}^{T} - \mathcal{O}_{p}\mathbf{R}_{xu}\mathbf{R}_{uu}^{-1} \mathbf{R}_{xu}^{T}\mathcal{O}_{p}^{T} -\\ -&+ - \mathbf{R}_{yu}\mathbf{R}_{uu}^{-1}\mathbf{R}_{yu}^{T} - \mathcal{O}_{p}\mathbf{R}_{xu} \mathbf{R}_{uu}^{-1}\mathbf{R}_{yu}^{T} -\\ -=& ~\mathcal{O}_{p}\mathbf{R}_{xx}\mathcal{O}_{p}^{T} - - \mathcal{O}_{p}\mathbf{R}_{xu}\mathbf{R}_{uu}^{-1} \mathbf{R}_{xu}^{T}\mathcal{O}_{p}^{T} + - \mathbf{R}_{yu}\mathbf{R}_{uu}^{-1}\mathbf{R}_{yu}^{T} -\end{aligned}$$ - -Moving all of the terms that can be composed with measured data to the left side: - -$$\begin{aligned} -\mathbf{R}_{yy} - \mathbf{R}_{yu}\mathbf{R}_{uu}^{-1}\mathbf{R}_{yu}^{T} -&= \mathcal{O}_{p}\mathbf{R}_{xx}\mathcal{O}_{p}^{T} - \mathcal{O}_{p}\mathbf{R}_{xu}\mathbf{R}_{uu}^{-1} \mathbf{R}_{xu}^{T}\mathcal{O}_{p}^{T} \\ -&= \mathcal{O}_{p}\left( \mathbf{R}_{xx} - \mathbf{R}_{xu}\mathbf{R}_{uu}^{-1} \mathbf{R}_{xu}^{T} \right) \mathcal{O}_{p}^{T} -\end{aligned}$$ - -We make the assumption that all current and future input data is uncorrelated with the current state, which means that the average of the products $\mathbf{x}(k)\mathbf{u}(k+i), ~~ i \in [0,1,2,\dots]$ is zero. This gives: - -$$\begin{aligned} -\mathbf{R}_{xu} &= -\frac{1}{N} -\begin{bmatrix} -\sum_{j=0}^{N-1}\mathbf{x}(k+j)\mathbf{u}(k+j) \\ -\sum_{j=0}^{N-1}\mathbf{x}(k+j)\mathbf{u}(k+j+1) \\ -\sum_{j=0}^{N-1}\mathbf{x}(k+j)\mathbf{u}(k+j+2) \\ -\vdots \\ -\sum_{j=0}^{N-1}\mathbf{x}(k+j)\mathbf{u}(k+j+p-1) -\end{bmatrix}^{T} \\ -&= -\mathbf{0} -\end{aligned}$$ - -in order to yield: - -$$ -\mathbf{R}_{yy} - \mathbf{R}_{yu}\mathbf{R}_{uu}^{-1}\mathbf{R}_{yu}^{T} = \mathcal{O}_{p}\mathbf{R}_{xx}\mathcal{O}_{p}^{T}~. -$$ \ No newline at end of file diff --git a/docs/classes/srim.rst b/docs/classes/srim.rst deleted file mode 100644 index c097777c9..000000000 --- a/docs/classes/srim.rst +++ /dev/null @@ -1,354 +0,0 @@ -System Realization by Information Matrix (SRIM) -=============================================== - -For discrete-time systems, the correlation between inputs, outputs, and -state yield information about the system's state evolution and response. -In fact, the state equations can be estimated by manipulating -correlation matrices through the method, `System Realization by -Information Matrix `__ (SRIM). - -Discrete-Time System Matrix Equation ------------------------------------- - -We begin with discrete-time state equations that correspond to the -structure's dynamics (see `Discrete LTI State-Space -Representation `__). - -.. math:: - - - \begin{aligned} - \bm{x}(k+1) &= \bm{Ax}(k) + \bm{Bu}(k) \\ - \bm{y}(k) &= \bm{Cx}(k) + \bm{Du}(k) - \end{aligned} - -By noting the state evolution - -.. math:: - - - \begin{aligned} - \bm{x}(k+1) &= \bm{Ax}(k)+\bm{B}\bm{U}_{p}(k)\\ - \bm{x}(k+2) &= \bm{A}^2\bm{X}(k) + \bm{AB}\bm{U}_{p}(k) + \bm{Bu}(k+1)\\ - \bm{x}(k+3) &= \bm{A}^{3}\bm{X}(k) + \bm{A}^{2}\bm{Bu}(k) + \bm{ABu}(k+1) + \bm{Bu}(k+2), - \end{aligned} - -we can generalize the response for the timepoint :math:`k+p-1`: - -.. math:: - - - \begin{aligned} - \bm{x}(k+p) &= \bm{A}^{p}\bm{x}(k) + \sum_{i=1}^{p}\bm{A}^{p-i}\bm{Bu}(k+i-1) - \\ - \bm{x}(k+p-1) &= \bm{A}^{p-1}\bm{x}(k) + \sum_{i=1}^{p-1}\bm{A}^{p-i-1}\bm{Bu}(k+i-1) - \\ - \bm{y}(k+p-1) &= \bm{CA}^{p-1}\bm{x}(k) + \sum_{i=1}^{p-1}\bm{CA}^{p-i-1}\bm{Bu}(k+i-1)+\bm{Du}(k+p-1)~. - \end{aligned} - -Then, we can vertically stack :math:`p` successive time-points into a -column vector and express this vector as :math:`\bm{y}_{p}(k)`: - -.. math:: - - - \begin{aligned} - \bm{y}_{p}(k) &= \mathcal{O}_{p}\bm{x}(k) + \mathcal{T}_{p}\bm{u}_{p}(k) \\ - \begin{bmatrix} - \bm{y}(k) \\ - \bm{y}(k+1) \\ - \vdots \\ - \bm{y}(k+p-1) - \end{bmatrix} - =& - \begin{bmatrix} - \bm{C} \\ - \bm{CA} \\ - \bm{CA}^{2} \\ - \vdots \\ - \bm{CA}^{p-1} - \end{bmatrix} - \bm{x}(k) - ~+ \\ - & - \begin{bmatrix} - \bm{D} \\ - \bm{CB} & \bm{D} \\ - \bm{CAB} & \bm{CB} & \bm{D} \\ - \vdots & \vdots & \vdots & \ddots \\ - \bm{CA}^{p-2}\bm{B} & \bm{CA}^{p-3}\bm{B} & \bm{CA}^{p-4}\bm{B} & \cdots & \bm{D} - \end{bmatrix} - \begin{bmatrix} - \bm{u}(k) \\ - \bm{u}(k+1) \\ - \vdots \\ - \bm{u}(k+p-1) - \end{bmatrix}~. - \end{aligned} - -.. raw:: html - - - -Finally, we horizontally stack :math:`N` successive timepoints of these -column vectors in a matrix, to get the matrix equation - -.. math:: - - - \boxed{\bm{Y}_{p}(k) = \mathcal{O}_{p}\bm{X}(k) + \mathcal{T}_{p}\bm{U}_{p}(k)} ~, - -where - -.. math:: - - - \begin{aligned} - \bm{Y}_{p}(k) &= \begin{bmatrix} \bm{y}_{p}(k) & \bm{y}_{p}(k+1) & \cdots & \bm{y}_{p}(k+N-1) \end{bmatrix} \\ - &= \begin{bmatrix} - \bm{y}(k) & \bm{y}(k+1) & \cdots & \bm{y}(k+N-1)\\ - \bm{y}(k+1) & \bm{y}(k+2) & \cdots & \bm{y}(k+N) \\ - \vdots & \vdots & \ddots & \vdots \\ - \bm{y}(k+p-1) & \bm{y}(k+p) & \cdots & \bm{y}(k+N+p-2) - \end{bmatrix} - \end{aligned} - -.. math:: - - - \bm{X}(k) = \begin{bmatrix} \bm{x}(k) & \bm{x}(k+1) & \cdots & \bm{x}(k+N-1) \end{bmatrix} - -.. math:: - - - \begin{aligned} - \bm{U}_{p}(k) &= \begin{bmatrix} \bm{u}_{p}(k) & \bm{u}_{p}(k+1) & \cdots & \bm{u}_{p}(k+N-1) \end{bmatrix} \\ - &= \begin{bmatrix} - \bm{u}(k) & \bm{u}(k+1) & \cdots & \bm{u}(k+N-1)\\ - \bm{u}(k+1) & \bm{u}(k+2) & \cdots & \bm{u}(k+N) \\ - \vdots & \vdots & \ddots & \vdots \\ - \bm{u}(k+p-1) & \bm{u}(k+p) & \cdots & \bm{u}(k+N+p-2) - \end{bmatrix}~. - \end{aligned} - -Observability Matrix from Information Matrix --------------------------------------------- - -By post-multiplying the matrix equation by -:math:`\frac{1}{N}\bm{U}_{p}^{T}(k)`, -:math:`\frac{1}{N}\bm{Y}_{p}^{T}(k)` or -:math:`\frac{1}{N}\bm{X}_{p}^{T}(k)`, we obtain relationships -between correlation matrices :math:`\bm{R}_{yy}`, -:math:`\bm{R}_{yu}`, :math:`\bm{R}_{uu}`, and -:math:`\bm{R}_{xx}` (See -`Appendix <#appendix-manipulation-of-discrete-time-system-matrix-equation-into-correlation-matrix-relationships>`__). - -.. math:: - - - \bm{R}_{yy} - \bm{R}_{yu}\bm{R}_{uu}^{-1}\bm{R}_{yu}^{T} = \mathcal{O}_{p}\bm{R}_{xx}\mathcal{O}_{p}^{T} ~, - -where - -.. math:: - - \begin{aligned} - \bm{R}_{yy} &= \frac{1}{N}\bm{Y}_{p}(k)\bm{Y}_{p}^{T}(k), \quad{} - \bm{R}_{yu} = \frac{1}{N}\bm{Y}_{p}(k)\bm{U}_{p}^{T}(k) \\ - \bm{R}_{uu} &= \frac{1}{N}\bm{U}_{p}(k)\bm{U}_{p}^{T}(k) , \quad{} - \bm{R}_{xx} = \frac{1}{N}\bm{X}(k)\bm{X}^{T}(k) ~. - \end{aligned} - -The left side of the equation is found from input and output -measurements, and is called the *information matrix* of the data. Its -singular value decomposition is computed to yield the *observability -matrix* :math:`\mathcal{O}_{p}`. - -.. math:: - - - \bm{R}_{yy} - \bm{R}_{yu}\bm{R}_{uu}^{-1}\bm{R}_{yu}^{T} = \bm{U} \Sigma \bm{U}^{T} = \mathcal{O}_{p}\bm{R}_{xx}\mathcal{O}_{p}^{T} ~. - -State Equation Matrices from Observability Matrix -------------------------------------------------- - -Now, the state equation matrices :math:`\bm{A}` and -:math:`\bm{C}` can be obtained from the observability matrix -:math:`\mathcal{O}_p`. - -.. math:: - - - \begin{aligned} - \mathcal{O}_{p} - = - \begin{bmatrix} - \bm{C} \\ - \bm{CA} \\ - \bm{CA}^{2} \\ - \vdots \\ - \bm{CA}^{p-1} - \end{bmatrix} - , \quad{} - \mathcal{O}_{p}(:-1) - = - \begin{bmatrix} - \bm{C} \\ - \bm{CA} \\ - \bm{CA}^{2} \\ - \vdots \\ - \bm{CA}^{p-2} - \end{bmatrix} - , \quad{} - \mathcal{O}_{p}(1:) - = - \begin{bmatrix} - \bm{CA} \\ - \bm{CA}^{2} \\ - \bm{CA}^{3} \\ - \vdots \\ - \bm{CA}^{p-1} - \end{bmatrix} - \end{aligned} - -.. math:: - - - \bm{A} = \mathcal{O}_{p}(:-1)^{+}\mathcal{O}_{p}(1:) - -.. math:: - - - \bm{C} = \mathcal{O}_{p}(0) - -Appendix: Manipulation of discrete-time system matrix equation into correlation matrix relationships ----------------------------------------------------------------------------------------------------- - -In (`Juang 1997 `__), the discrete-time -system matrix equation is manipulated into a form describing the -relationship between correlation matrices :math:`\bm{R}_{yy}`, -:math:`\bm{R}_{yu}`, :math:`\bm{R}_{uu}`, and -:math:`\bm{R}_{xx}`. - -Post-multiplying the `discrete-time system matrix -equation <#discrete-time-system-matrix-equation>`__ by -:math:`\frac{1}{N}\bm{U}_{p}^{T}(k)`: - -.. math:: - - \begin{aligned} - \frac{1}{N}\bm{Y}_{p}(k)\bm{U}_{p}^{T}(k) &= \mathcal{O}_{p}\frac{1}{N}\bm{X}(k)\bm{U}_{p}^{T}(k) + \mathcal{T}_{p}\frac{1}{N}\bm{U}_{p}(k)\bm{U}_{p}^{T}(k) - \\ - \bm{R}_{yu} &= \mathcal{O}_{p}\bm{R}_{xu} + \mathcal{T}_{p}\bm{R}_{uu} - \\ - \mathcal{T}_{p} &= \left( \bm{R}_{yu} - \mathcal{O}_{p}\bm{R}_{xu} \right)\bm{R}_{uu}^{-1} - \end{aligned} - -Post-multiplying by :math:`\frac{1}{N}\bm{Y}_{p}^{T}(k)`: - -.. math:: - - \begin{aligned} - \frac{1}{N}\bm{Y}_{p}(k)\bm{Y}_{p}^{T}(k) &= \mathcal{O}_{p}\frac{1}{N}\bm{X}(k)\bm{Y}_{p}^{T}(k) + \mathcal{T}_{p}\frac{1}{N}\bm{U}_{p}(k)\bm{Y}_{p}^{T}(k) - \\ - \bm{R}_{yy} &= \mathcal{O}_{p}\bm{R}_{yx}^{T} + \mathcal{T}_{p}\bm{R}_{yu}^{T} - \\ - \bm{R}_{yy} &= \mathcal{O}_{p}\bm{R}_{yx}^{T} + \left( \bm{R}_{yu} - \mathcal{O}_{p}\bm{R}_{xu} \right)\bm{R}_{uu}^{-1}\bm{R}_{yu}^{T} - \end{aligned} - -Post-multiplying by :math:`\frac{1}{N}\bm{X}_{p}^{T}(k)`: - -.. math:: - - \begin{aligned} - \frac{1}{N}\bm{Y}_{p}(k)\bm{X}_{p}^{T}(k) &= \mathcal{O}_{p}\frac{1}{N}\bm{X}(k)\bm{X}_{p}^{T}(k) + \mathcal{T}_{p}\frac{1}{N}\bm{U}_{p}(k)\bm{X}_{p}^{T}(k) - \\ - \bm{R}_{yx} &= \mathcal{O}_{p}\bm{R}_{xx} + \mathcal{T}_{p}\bm{R}_{xu}^{T} - \\ - \bm{R}_{yx} &= \mathcal{O}_{p}\bm{R}_{xx} + \left( \bm{R}_{yu} - \mathcal{O}_{p}\bm{R}_{xu} \right)\bm{R}_{uu}^{-1}\bm{R}_{xu}^{T} - \end{aligned} - -Substituting the equation for :math:`\bm{R}_{yx}` into the equation -for :math:`\bm{R}_{yy}`: - -.. math:: - - \begin{aligned} - \bm{R}_{yy} =& ~\mathcal{O}_{p} - \left(\mathcal{O}_{p}\bm{R}_{xx} + \left( \bm{R}_{yu} - \mathcal{O}_{p}\bm{R}_{xu} \right)\bm{R}_{uu}^{-1}\bm{R}_{xu}^{T}\right)^{T} - \\ - &+ - \left( \bm{R}_{yu} - \mathcal{O}_{p}\bm{R}_{xu} \right)\bm{R}_{uu}^{-1}\bm{R}_{yu}^{T} - \\ - =& ~\mathcal{O}_{p}\bm{R}_{xx}\mathcal{O}_{p}^{T} - + \mathcal{O}_{p}\bm{R}_{xu}\bm{R}_{uu}^{-1} \left( \bm{R}_{yu}^{T} - \bm{R}_{xu}^{T}\mathcal{O}_{p}^{T} \right) - \\ - &+ - \left( \bm{R}_{yu} - \mathcal{O}_{p}\bm{R}_{xu} \right)\bm{R}_{uu}^{-1}\bm{R}_{yu}^{T} - \\ - =& ~\mathcal{O}_{p}\bm{R}_{xx}\mathcal{O}_{p}^{T} - + \mathcal{O}_{p}\bm{R}_{xu}\bm{R}_{uu}^{-1} \bm{R}_{yu}^{T} - \mathcal{O}_{p}\bm{R}_{xu}\bm{R}_{uu}^{-1} \bm{R}_{xu}^{T}\mathcal{O}_{p}^{T} - \\ - &+ - \bm{R}_{yu}\bm{R}_{uu}^{-1}\bm{R}_{yu}^{T} - \mathcal{O}_{p}\bm{R}_{xu} \bm{R}_{uu}^{-1}\bm{R}_{yu}^{T} - \\ - =& ~\mathcal{O}_{p}\bm{R}_{xx}\mathcal{O}_{p}^{T} - - \mathcal{O}_{p}\bm{R}_{xu}\bm{R}_{uu}^{-1} \bm{R}_{xu}^{T}\mathcal{O}_{p}^{T} + - \bm{R}_{yu}\bm{R}_{uu}^{-1}\bm{R}_{yu}^{T} - \end{aligned} - -Moving all of the terms that can be composed with measured data to the -left side: - -.. math:: - - \begin{aligned} - \bm{R}_{yy} - \bm{R}_{yu}\bm{R}_{uu}^{-1}\bm{R}_{yu}^{T} - &= \mathcal{O}_{p}\bm{R}_{xx}\mathcal{O}_{p}^{T} - \mathcal{O}_{p}\bm{R}_{xu}\bm{R}_{uu}^{-1} \bm{R}_{xu}^{T}\mathcal{O}_{p}^{T} \\ - &= \mathcal{O}_{p}\left( \bm{R}_{xx} - \bm{R}_{xu}\bm{R}_{uu}^{-1} \bm{R}_{xu}^{T} \right) \mathcal{O}_{p}^{T} - \end{aligned} - -We make the assumption that all current and future input data is -uncorrelated with the current state, which means that the average of the -products :math:`\bm{x}(k)\bm{u}(k+i), ~~ i \in [0,1,2,\dots]` is -zero. This gives: - -.. math:: - - \begin{aligned} - \bm{R}_{xu} &= - \frac{1}{N} - \begin{bmatrix} - \sum_{j=0}^{N-1}\bm{x}(k+j)\bm{u}(k+j) \\ - \sum_{j=0}^{N-1}\bm{x}(k+j)\bm{u}(k+j+1) \\ - \sum_{j=0}^{N-1}\bm{x}(k+j)\bm{u}(k+j+2) \\ - \vdots \\ - \sum_{j=0}^{N-1}\bm{x}(k+j)\bm{u}(k+j+p-1) - \end{bmatrix}^{T} \\ - &= - \bm{0} - \end{aligned} - -in order to yield: - -.. math:: - - - \bm{R}_{yy} - \bm{R}_{yu}\bm{R}_{uu}^{-1}\bm{R}_{yu}^{T} = \mathcal{O}_{p}\bm{R}_{xx}\mathcal{O}_{p}^{T}~. - diff --git a/docs/classes/statespace.rst b/docs/classes/statespace.rst deleted file mode 100644 index dea21d2a4..000000000 --- a/docs/classes/statespace.rst +++ /dev/null @@ -1,149 +0,0 @@ -State Space Model of Structural Dynamics ----------------------------------------- - -.. figure:: figures/si_msmdof.png - :alt: MDOF Structure - - MDOF Structure - -When an multiple degree-of-freedom (MDOF) system is subject to multiple -support excitation, such as in the figure above, the displacement vector -is extended to include the support DOF. An `equation of -motion <#equation-of-motion>`__ is derived as follows. - -Begin by forming a partitioned equation of dynamic equilibrium for all -the DOF: - -Partitioned Equation of Dynamic Equilibrium -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. math:: - - - \begin{bmatrix} \mathbf{m} & \mathbf{m}_{g} \\ \mathbf{m}^T_{g} & \mathbf{m}_{gg} \end{bmatrix} - \begin{bmatrix} \mathbf{\ddot{u}}^{t}_{f} \\ \mathbf{\ddot{u}}_{g} \end{bmatrix} - + - \begin{bmatrix} \mathbf{c} & \mathbf{c}_{g} \\ \mathbf{c}^T_{g} & \mathbf{c}_{gg} \end{bmatrix} - \begin{bmatrix} \mathbf{\dot{u}}^{t}_{f} \\ \mathbf{\dot{u}}_{g} \end{bmatrix} - + - \begin{bmatrix} \mathbf{k} & \mathbf{k}_{g} \\ \mathbf{k}^T_{g} & \mathbf{k}_{gg} \end{bmatrix} - \begin{bmatrix} \mathbf{u}^{t}_{f} \\ \mathbf{u}_{g} \end{bmatrix} - = - \begin{bmatrix} \mathbf{0} \\ \mathbf{p}_{g} \end{bmatrix} - -where the subscript :math:`g` indicates support DOF, the subscript -:math:`f` indicates structural DOF, and the superscript :math:`t` -indicates the total of quasi-static (:math:`\mathbf{u}^{s}_{f}`, due to -static application of support displacements) and dynamic -(:math:`\mathbf{u}_{f}`, evaluated by dynamic analysis) structural -displacements. - -Taking the first half of the partitioned equilibrium, separating the -structural displacements -(:math:`\mathbf{u}^{t}_{f}=\mathbf{u}^{s}_{f}+\mathbf{u}_{f}`), and -moving all :math:`\mathbf{u}_{g}` and :math:`\mathbf{u}^{s}_{f}` terms -to the right side, - -.. math:: - - - \mathbf{m}\mathbf{\ddot{u}}_{f} + \mathbf{c}\mathbf{\dot{u}}_{f} + \mathbf{k}\mathbf{u}_{f} - = -(\mathbf{m}\mathbf{\ddot{u}}^{s}_{f}+\mathbf{m}_{g}\mathbf{\ddot{u}}_{g}) - -(\mathbf{c}\mathbf{\dot{u}}^{s}_{f}+\mathbf{c}_{g}\mathbf{\dot{u}}_{g}) - -(\mathbf{k}\mathbf{u}^{s}_{f}+\mathbf{k}_{g}\mathbf{u}_{g}) - -The term -:math:`(\mathbf{k}\mathbf{u}^{s}_{f}+\mathbf{k}_{g}\mathbf{u}_{g})=\mathbf{0}` -due to static equilibrium, allowing the term to be dropped and giving -:math:`\mathbf{u}^{s}_{f} = \mathbf{-k}^{-1}\mathbf{k}_{g}\mathbf{u}_{g} = \mathbf{\iota u}_{g}`; -the term -:math:`(\mathbf{c}\mathbf{\dot{u}}^{s}_{f}+\mathbf{c}_{g}\mathbf{\dot{u}}_{g})` -is dropped because it is usually small relative to the inertia term; and -the term :math:`\mathbf{m}_{g}\mathbf{\ddot{u}}_{g}` is dropped because -mass is usually neglected at supports. - -The equilibrium equation thus simplifies. - -Equation of Motion -^^^^^^^^^^^^^^^^^^ - -.. math:: - - - \begin{aligned} - \mathbf{M\ddot{u}}_{f}(t) + \mathbf{Z\dot{u}}_{f}(t) + \mathbf{Ku}_{f}(t) &= -\mathbf{M\iota}\mathbf{\ddot{u}}_{g}(t) \\ - \mathbf{m}\mathbf{\ddot{u}}_{f} + \mathbf{c}\mathbf{\dot{u}}_{f} + \mathbf{k}\mathbf{u}_{f} - &= -\mathbf{m}\mathbf{\iota}\mathbf{\ddot{u}}_{g} - \end{aligned} - -Hence, the following equation presents the continuous linear -time-invariant (LTI) state-space representation of a structural system. - -Continuous LTI State-Space Representation -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. math:: - - - \begin{aligned} - \mathbf{\dot{x}} &= \mathbf{A}_{c}\mathbf{x} + \mathbf{B}_{c}\mathbf{u} \\ - \begin{bmatrix} \mathbf{\dot{u}}_{f}(t) \\ \mathbf{\ddot{u}}_{f}(t) \end{bmatrix} - &= - \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{Z} \end{bmatrix} - \begin{bmatrix} \mathbf{u}_{f}(t) \\ \mathbf{\dot{u}}_{f}(t) \end{bmatrix} - + - \begin{bmatrix} \mathbf{0} \\ -\mathbf{\iota} \end{bmatrix} - \mathbf{\ddot{u}}_{g}(t) \\ \\ - \mathbf{y} &= \mathbf{Cx} + \mathbf{Du} \\ - \mathbf{\ddot{u}}_{f}(t) &= - \begin{bmatrix} -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{Z} \end{bmatrix} - \begin{bmatrix} \mathbf{u}_{f}(t) \\ \mathbf{\dot{u}}_{f}(t) \end{bmatrix} - + - \begin{bmatrix} -\mathbf{\iota} \end{bmatrix} - \mathbf{\ddot{u}}_{g}(t) - \end{aligned} - -In order to move from the continuous to the discrete case, the -coefficients :math:`\mathbf{A}_{c}` and :math:`\mathbf{B}_{c}` are -transformed by solving the first-order differential equation with the -signal’s value held constant between time steps (“zero-order hold -method”). The coefficients :math:`\mathbf{C}` and :math:`\mathbf{D}` are -unchanged. The results are shown in the following equation. - -.. figure:: figures/si_discretize.png - :alt: Signal Discretization - - Signal Discretization - -Discrete LTI State-Space Representation -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -.. math:: - - - \begin{aligned} - \mathbf{x}_{k+1} &= \mathbf{Ax}_{k} + \mathbf{Bu}_{k} \\ - \mathbf{y}_{k} &= \mathbf{Cx}_{k} + \mathbf{Du}_{k} \\ - \end{aligned} - -.. math:: - - - \mathbf{x}_{k} = \mathbf{x}(k\Delta t), \hspace{1cm} \mathbf{u}_{k} = \mathbf{u}(k\Delta t), \hspace{1cm} \mathbf{y}_{k} = \mathbf{y}(k\Delta t) - -.. math:: - - - \mathbf{A} = e^{\mathbf{A}_{c}\Delta t}, \hspace{1cm} \mathbf{B} = \int_{0}^{\Delta t}{e^{\mathbf{A}_{c}\tau}}\mathbf{B}_{c}d\tau - -where: - -.. math:: - - - \begin{aligned} - \mathbf{A} & \text{: discrete state transition matrix} \\ - \mathbf{B} & \text{: discrete input influence matrix} \\ - \mathbf{C} & \text{: output influence matrix} \\ - \mathbf{D} & \text{: direct transmission or feed-through matrix} - \end{aligned} diff --git a/docs/conf.py b/docs/conf.py index a4127e3bb..7ecf6ee65 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -47,12 +47,13 @@ 'css/theme-css/'+str(file.name) for file in (Path(__file__).parents[0]/"_static/css/theme-css/").glob("*.css") ] html_additional_pages = {'index': 'home.html'} +g = "https://gallery.stairlab.io" html_context = { 'description': description, 'examples': [ - {"title": "Basics", "link": "examples/00_Overview", "image": "../_static/images/gallery/overview.svg"}, - {"title": "Displacements", "link": "examples/01_SISO_Intro", "image": "../_static/images/gallery/sdof_full.svg"}, - {"title": "Motions", "link": "examples/04_MIMO_Intro", "image": "../_static/images/gallery/2dof_full.svg"}, + {"title": "Basics", "link": f"{g}/examples/example7/", "image": "../_static/images/gallery/safeway.png"}, + {"title": "Displacements", "link": f"{g}/examples/example5/", "image": "../_static/images/gallery/Example5.png"}, + {"title": "Motions", "link": f"{g}/examples/example6/", "image": "../_static/images/gallery/Example6.png"}, ], **globals() } diff --git a/docs/design/index.rst b/docs/design/index.rst new file mode 100644 index 000000000..3596f5684 --- /dev/null +++ b/docs/design/index.rst @@ -0,0 +1,8 @@ +Design +====== + +``veux`` is designed around a set of core abstract classes: + +- A ``Canvas`` abstracts away details about the rendering backend. +- An ``Artist`` abstracts away the drawing process. It owns a canvas. + diff --git a/docs/index.rst b/docs/index.rst index 7bac2943c..3b874465e 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -8,5 +8,5 @@ veux documentation :caption: Contents: library/index - classes/index + design/index diff --git a/docs/library/index.rst b/docs/library/index.rst index a0c4fdb6f..a5e0f3661 100644 --- a/docs/library/index.rst +++ b/docs/library/index.rst @@ -1,5 +1,5 @@ -veux -==== +Library +======= .. currentmodule:: veux @@ -26,8 +26,8 @@ frequency domain and time domain analysis of vibration signals: .. toctree:: :maxdepth: 1 - veux.server veux.canvas + veux.server .. top-level functions .. -------------------