newmark_ss

sharpy.linear.src.lingebm.newmark_ss(M, C, K, dt, num_damp=0.0001, M_is_SPD=False)

Produces a discrete-time state-space model of the 2nd order ordinary differential equation (ODE) given by:

\[\mathbf{M}\mathbf{\ddot{q}} + \mathbf{C}\mathbf{\dot{q}} + \mathbf{K}\mathbf{q} = \mathbf{f(t)}\]

This ODE is discretized based on the Newmark-\(\beta\) integration scheme.

The output state-space model has the form:

\[\begin{split}\mathbf{x}_{n+1} &= \mathbf{A_{ss}}\mathbf{x}_n + \mathbf{B_{ss}}\mathbf{f}_n \\ \mathbf{y}_n &= \mathbf{C_{ss}}\mathbf{x}_n + \mathbf{D_{ss}}\mathbf{f}_n\end{split}\]

where \(\mathbf{y} = \begin{Bmatrix} \mathbf{q} \\ \mathbf{\dot{q}} \end{Bmatrix}\)

Note that as the state-space representation only requires the input force \(\mathbf{f}\) to be evaluated at time-step \(n\), thus the pass-through matrix \(\mathbf{D_{ss}}\) is not zero.

This function retuns a tuple with the discrete state-space matrices \((\mathbf{A_{ss},B_{ss},C_{ss},D_{ss}})\).

Theory

The following steps describe how to apply the Newmark-\(\beta\) scheme to the ODE in order to generate the discrete time-state space-model. It folows the development of [1].

Notation

Bold upper case letters represent matrices, bold lower case letters represent vectors. Non-bold symbols are scalars. Curly brackets indicate (block) vectors and square brackets indicate (block) matrices.

Evaluating the ODE to the time steps \(t_n\) and \(t_{n+1}\) and isolating the acceleration term:

\[\begin{split}\mathbf{\ddot q}_{n} &= - \mathbf{M}^{-1}\mathbf{C}\mathbf{\dot{q}}_{n} - \mathbf{M}^{-1}\mathbf{K}\mathbf{q}_{n} + \mathbf{M}^{-1}\mathbf{f}_{n} \\ \mathbf{\ddot q}_{n+1} &= - \mathbf{M}^{-1}\mathbf{C}\mathbf{\dot{q}}_{n+1} - \mathbf{M}^{-1}\mathbf{K}\mathbf{q}_{n+1} + \mathbf{M}^{-1}\mathbf{f}_{n+1} \\\end{split}\]

The update equations of the Newmark-beta scheme are [1]:

\[\begin{split}\mathbf{q}_{n+1} &= \mathbf{q}_n + \mathbf{\dot q}_n\Delta t + (1/2-\beta)\mathbf{\ddot q}_n \Delta t^2 + \beta \mathbf{\ddot q}_{n+1} \Delta t^2 + O(\Delta t^3) \\ \mathbf{\dot q}_{n+1} &= \mathbf{\dot q}_n + (1-\gamma)\mathbf{\ddot q}_n \Delta t + \gamma \mathbf{\ddot q}_{n+1} \Delta t + O(\Delta t^3)\end{split}\]

where \(\Delta t = t_{n+1} - t_n\).

The stencil is unconditionally stable if the tuning parameters \(\gamma\) and \(\beta\) are chosen as:

\[\begin{split}\gamma &= \frac{1}{2} + \alpha \\ \beta &= \frac{(1+\alpha)^2}{4} = \frac{(1/2+\gamma)^2}{4} = \frac{1}{16} + \frac{1}{4}(\gamma + \gamma^2)\end{split}\]

where \(\alpha>0\) accounts for small positive algorithmic damping (\(\alpha\) is num_damp in the code).

Substituting the former relations onto the later ones, rearranging terms, and writing it in state-space form:

\[\begin{split}\mathbf{A_{ss1}} \begin{Bmatrix} \mathbf{q}_{n+1} \\ \mathbf{\dot q}_{n+1} \end{Bmatrix} = \mathbf{A_{ss0}} \begin{Bmatrix} \mathbf{q}_{n} \\ \mathbf{\dot q}_{n} \end{Bmatrix} + \mathbf{B_{ss0}} \mathbf{f}_n + \mathbf{B_{ss1}} \mathbf{f}_{n+1}\end{split}\]

where

\[\begin{split}\mathbf{A_{ss1}} &= \begin{bmatrix} \mathbf{I} + \beta\Delta t^2 \mathbf{M}^{-1}\mathbf{K} & \beta\Delta t^2 \mathbf{M}^{-1}\mathbf{C} \\ \gamma \Delta t \mathbf{M}^{-1}\mathbf{K} & \mathbf{I}+ \gamma \Delta t \mathbf{M}^{-1}\mathbf{C} \end{bmatrix} \\ \mathbf{A_{ss0}} &= \begin{bmatrix} \mathbf{I} - \Delta t^2(1/2-\beta)\mathbf{M}^{-1}\mathbf{K} & \Delta t \mathbf{I} - (1/2-\beta)\Delta t^2 \mathbf{M}^{-1}\mathbf{C} \\ -(1-\gamma)\Delta t \mathbf{M}^{-1}\mathbf{K} & \mathbf{I} - (1-\gamma)\Delta t \mathbf{M}^{-1}\mathbf{C} \end{bmatrix} \\ \mathbf{B_{ss0}} &= \begin{bmatrix} (\Delta t^2(1/2-\beta) \mathbf{M}^{-1} \\ (1-\gamma)\Delta t \mathbf{M}^{-1} \end{bmatrix} \\ \mathbf{B_{ss1}} &= \begin{bmatrix} (\beta \Delta t^2) \mathbf{M}^{-1}\\ (\gamma \Delta t) \mathbf{M}^{-1} \end{bmatrix}\end{split}\]

This is not in standard space-state form because the state update equation depends of the input at \(t_{n+1}\). This term can be eliminated by defining the state

\[\begin{split}\mathbf{x}_n = \begin{Bmatrix} \mathbf{q}_{n} \\ \mathbf{\dot q}_{n} \end{Bmatrix} - \mathbf{A}_{\mathbf{ss1}}^{-1}\mathbf{B_{ss1}} \mathbf{f}_n\end{split}\]

Then

\[\begin{split}\mathbf{x}_{n+1} &= \mathbf{A}_{\mathbf{ss1}}^{-1}[ \mathbf{A_{ss0}} \mathbf{x}_n + ( \mathbf{A_{ss0}}\mathbf{A}_{\mathbf{ss1}}^{-1}\mathbf{B_{ss1}} + \mathbf{B_{ss0}} ) \mathbf{f}_n ] \\ \begin{Bmatrix} \mathbf{q}_{n} \\ \mathbf{\dot q}_{n} \end{Bmatrix} &= \mathbf{x}_n + \mathbf{A}_{\mathbf{ss1}}^{-1}\mathbf{B_{ss1}} \mathbf{f}_n\end{split}\]

See also sharpy.linear.src.libss.SSconv() for more details on the elimination of the term multiplying \(\mathbf{f}_{n+1}\) in the state equation.

This system is identified with a standard discrete-time state-space model

\[\begin{split}\mathbf{x}_{n+1} &= \mathbf{A_{ss}} \mathbf{x}_n + \mathbf{B_{ss}} \mathbf{f}_n \\ \mathbf{y_n} &= \mathbf{C_{ss}} \mathbf{x}_n + \mathbf{D_{ss}} \mathbf{f}_n\end{split}\]

where

\[\begin{split}\mathbf{A_{ss}} &= \mathbf{A}_{\mathbf{ss1}}^{-1}\mathbf{A_{ss0}} \\ \mathbf{B_{ss}} &= \mathbf{A_{\mathbf{ss1}}^{-1}(\mathbf{B_{ss0}} + \mathbf{A_{ss0}}\mathbf{A_{\mathbf{ss1}}^{-1}\mathbf{B_{ss1}}) \\ \mathbf{C_{ss}} &= \mathbf{I} \\ \mathbf{D_{ss}} &= \mathbf{A_{\mathbf{ss1}}^{-1}\mathbf{B_{ss1}}\end{split}\]

Notation is used in the code

\[\begin{split}\texttt{th1} &= \gamma \\ \texttt{th2} &= \beta \\ \texttt{a0} &= (1/2 -\beta) \Delta t^2 \\ \texttt{b0} &= (1 -\gamma) \Delta t \\ \texttt{a1} &= \beta \Delta t^2 \\ \texttt{b1} &= \gamma \Delta t \\\end{split}\]

param M:

Mass matrix \(\mathbf{M}\)

type M:

np.array

param C:

Damping matrix \(\mathbf{C}\)

type C:

np.array

param K:

Stiffness matrix \(\mathbf{K}\)

type K:

np.array

param dt:

Timestep increment

type dt:

float

param num_damp:

Numerical damping. Default 1e-4

type num_damp:

float

param M_is_SPD:

whether to factorized M using Cholesky (only works for SPD matrices) or LU decomposition. Default: False

type M_is_SPD:

bool

returns:

with the \((\mathbf{A_{ss},B_{ss},C_{ss},D_{ss}})\) matrices of the discrete-time state-space representation

rtype:

tuple

References

[1] - Geradin M., Rixen D. - Mechanical Vibrations: Theory and application to structural dynamics