From b0de5d18f4ae563395117482a975fb6298a44b9c Mon Sep 17 00:00:00 2001 From: Bernardo Bahia Monteiro Date: Thu, 26 Oct 2023 21:27:09 -0400 Subject: [PATCH 01/13] [docs] Improve formatting of matrices in newmark_ss --- sharpy/linear/src/lingebm.py | 35 +++++++++++++++++++++++++++-------- 1 file changed, 27 insertions(+), 8 deletions(-) diff --git a/sharpy/linear/src/lingebm.py b/sharpy/linear/src/lingebm.py index 83cfefdc4..51ea90c66 100644 --- a/sharpy/linear/src/lingebm.py +++ b/sharpy/linear/src/lingebm.py @@ -1422,14 +1422,33 @@ def newmark_ss(Minv, C, K, dt, num_damp=1e-4): Substituting the former relation onto the later ones, rearranging terms, and writing it in state-space form: .. math:: - \begin{bmatrix} I + M^{-1}K \Delta t^2\beta \quad \Delta t^2\beta M^{-1}C \\ (\gamma \Delta t M^{-1}K) - \quad (I + \gamma \Delta t M^{-1}C) \end{bmatrix} \begin{Bmatrix} \mathbf{\dot q}_{n+1} \\ - \mathbf{\ddot q}_{n+1} \end{Bmatrix} = - \begin{bmatrix} (I - \Delta t^2(1/2-\beta)M^{-1}K \quad (\Delta t - \Delta t^2(1/2-\beta)M^{-1}C \\ - (-(1-\gamma)\Delta t M^{-1}K \quad (I - (1-\gamma)\Delta tM^{-1}C \end{bmatrix} - \begin{Bmatrix} \mathbf{q}_{n} \\ \mathbf{\dot q}_{n} \end{Bmatrix} + - \begin{Bmatrix} (\Delta t^2(1/2-\beta) \\ (1-\gamma)\Delta t \end{Bmatrix} M^{-1}F_n+ - \begin{Bmatrix} (\Delta t^2\beta) \\ (\gamma \Delta t) \end{Bmatrix}M^{-1}F_{n+1} + \begin{bmatrix} + I + M^{-1}K \Delta t^2\beta & \Delta t^2\beta M^{-1}C \\ + (\gamma \Delta t M^{-1}K) & (I + \gamma \Delta t M^{-1}C) + \end{bmatrix} + \begin{Bmatrix} + \mathbf{\dot q}_{n+1} \\ + \mathbf{\ddot q}_{n+1} + \end{Bmatrix} + = + \begin{bmatrix} + (I - \Delta t^2(1/2-\beta)M^{-1}K & (\Delta t - \Delta t^2(1/2-\beta)M^{-1}C \\ + (-(1-\gamma)\Delta t M^{-1}K & (I - (1-\gamma)\Delta tM^{-1}C + \end{bmatrix} + \begin{Bmatrix} + \mathbf{q}_{n} \\ \ + mathbf{\dot q}_{n} + \end{Bmatrix} + + + \begin{Bmatrix} + (\Delta t^2(1/2-\beta) \\ + (1-\gamma)\Delta t + \end{Bmatrix} + M^{-1}F_n+ + \begin{Bmatrix} + (\Delta t^2\beta) \\ + (\gamma \Delta t) + \end{Bmatrix}M^{-1}F_{n+1} To understand SHARPy code, it is convenient to apply the following change of notation: From fbb1dee822aad0a0a3c30331274352ac837699ec Mon Sep 17 00:00:00 2001 From: Bernardo Bahia Monteiro Date: Thu, 26 Oct 2023 21:31:53 -0400 Subject: [PATCH 02/13] [docs] Remove incorrect dots since the system is discrete (newmark_ss) --- sharpy/linear/src/lingebm.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/sharpy/linear/src/lingebm.py b/sharpy/linear/src/lingebm.py index 51ea90c66..c3194d779 100644 --- a/sharpy/linear/src/lingebm.py +++ b/sharpy/linear/src/lingebm.py @@ -1427,8 +1427,8 @@ def newmark_ss(Minv, C, K, dt, num_damp=1e-4): (\gamma \Delta t M^{-1}K) & (I + \gamma \Delta t M^{-1}C) \end{bmatrix} \begin{Bmatrix} - \mathbf{\dot q}_{n+1} \\ - \mathbf{\ddot q}_{n+1} + \mathbf{q}_{n+1} \\ + \mathbf{\dot q}_{n+1} \end{Bmatrix} = \begin{bmatrix} From 93142c579a2ab3bb576db4b36c2a26cf8f079b90 Mon Sep 17 00:00:00 2001 From: Bernardo Bahia Monteiro Date: Thu, 26 Oct 2023 21:51:13 -0400 Subject: [PATCH 03/13] [docs] Fix and improve formatting of equations - Fix equations to match the ones in the code, correct parenthesis, etc - Bring M^{-1} term to the inside of the matrix multiplying the forces, for consistency of dimensions. Replace curly brackets by square brackets since it is a matrix and not a vector --- sharpy/linear/src/lingebm.py | 39 ++++++++++++++++++------------------ 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/sharpy/linear/src/lingebm.py b/sharpy/linear/src/lingebm.py index c3194d779..b569f6725 100644 --- a/sharpy/linear/src/lingebm.py +++ b/sharpy/linear/src/lingebm.py @@ -1399,7 +1399,7 @@ def newmark_ss(Minv, C, K, dt, num_damp=1e-4): The following steps describe how to apply the Newmark-beta scheme to a state-space formulation. The original idea is based on [1]. - The equation of a second order system dynamics reads: + The equation of a second order dynamical system reads: .. math:: M\mathbf{\ddot q} + C\mathbf{\dot q} + K\mathbf{q} = F @@ -1415,7 +1415,7 @@ def newmark_ss(Minv, C, K, dt, num_damp=1e-4): .. math:: \mathbf{q}_{n+1} &= \mathbf{q}_n + \mathbf{\dot q}_n\Delta t + - (\frac{1}{2}-\beta)\mathbf{\ddot q}_n \Delta t^2 + \beta \mathbf{\ddot q}_{n+1} \Delta t^2 + O(\Delta t^3) \\ + (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) @@ -1423,8 +1423,8 @@ def newmark_ss(Minv, C, K, dt, num_damp=1e-4): .. math:: \begin{bmatrix} - I + M^{-1}K \Delta t^2\beta & \Delta t^2\beta M^{-1}C \\ - (\gamma \Delta t M^{-1}K) & (I + \gamma \Delta t M^{-1}C) + I + \beta\Delta t^2 M^{-1}K & \beta\Delta t^2 M^{-1}C \\ + \gamma \Delta t M^{-1}K & I + \gamma \Delta t M^{-1}C \end{bmatrix} \begin{Bmatrix} \mathbf{q}_{n+1} \\ @@ -1432,33 +1432,34 @@ def newmark_ss(Minv, C, K, dt, num_damp=1e-4): \end{Bmatrix} = \begin{bmatrix} - (I - \Delta t^2(1/2-\beta)M^{-1}K & (\Delta t - \Delta t^2(1/2-\beta)M^{-1}C \\ - (-(1-\gamma)\Delta t M^{-1}K & (I - (1-\gamma)\Delta tM^{-1}C + I - \Delta t^2(1/2-\beta)M^{-1}K & \Delta t I - (1/2-\beta)\Delta t^2 M^{-1}C \\ + -(1-\gamma)\Delta t M^{-1}K & I - (1-\gamma)\Delta t M^{-1}C \end{bmatrix} \begin{Bmatrix} \mathbf{q}_{n} \\ \ mathbf{\dot q}_{n} \end{Bmatrix} + - \begin{Bmatrix} - (\Delta t^2(1/2-\beta) \\ - (1-\gamma)\Delta t - \end{Bmatrix} - M^{-1}F_n+ - \begin{Bmatrix} - (\Delta t^2\beta) \\ - (\gamma \Delta t) - \end{Bmatrix}M^{-1}F_{n+1} + \begin{bmatrix} + (\Delta t^2(1/2-\beta)\, M^{-1} \\ + (1-\gamma)\Delta t \, M^{-1} + \end{bmatrix} + F_n + + \begin{bmatrix} + (\beta \Delta t^2) \, M^{-1}\\ + (\gamma \Delta t) \, M^{-1} + \end{bmatrix} + F_{n+1} To understand SHARPy code, it is convenient to apply the following change of notation: .. math:: \textrm{th1} = \gamma \\ \textrm{th2} = \beta \\ - \textrm{a0} = \Delta t^2 (1/2 -\beta) \\ - \textrm{b0} = \Delta t (1 -\gamma) \\ - \textrm{a1} = \Delta t^2 \beta \\ - \textrm{b1} = \Delta t \gamma \\ + \textrm{a0} = (1/2 -\beta) \Delta t^2 \\ + \textrm{b0} = (1 -\gamma) \Delta t \\ + \textrm{a1} = \beta \Delta t^2 \\ + \textrm{b1} = \gamma \Delta t \\ Finally: From 8c77819d25e812db2f1f78b23c4e67184313f68a Mon Sep 17 00:00:00 2001 From: Bernardo Bahia Monteiro Date: Thu, 26 Oct 2023 21:53:59 -0400 Subject: [PATCH 04/13] [docs] organize text in more logical order SHARPy specific notation goes last --- sharpy/linear/src/lingebm.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/sharpy/linear/src/lingebm.py b/sharpy/linear/src/lingebm.py index b569f6725..bebfcfa44 100644 --- a/sharpy/linear/src/lingebm.py +++ b/sharpy/linear/src/lingebm.py @@ -1451,16 +1451,6 @@ def newmark_ss(Minv, C, K, dt, num_damp=1e-4): \end{bmatrix} F_{n+1} - To understand SHARPy code, it is convenient to apply the following change of notation: - - .. math:: - \textrm{th1} = \gamma \\ - \textrm{th2} = \beta \\ - \textrm{a0} = (1/2 -\beta) \Delta t^2 \\ - \textrm{b0} = (1 -\gamma) \Delta t \\ - \textrm{a1} = \beta \Delta t^2 \\ - \textrm{b1} = \gamma \Delta t \\ - Finally: .. math:: @@ -1472,6 +1462,16 @@ def newmark_ss(Minv, C, K, dt, num_damp=1e-4): To finally isolate the vector at :math:`n+1`, instead of inverting the :math:`A_{ss1}` matrix, several systems are solved. Moreover, the output equation is simply :math:`y=x`. + To understand SHARPy code, it is convenient to apply the following change of notation: + + .. math:: + \textrm{th1} = \gamma \\ + \textrm{th2} = \beta \\ + \textrm{a0} = (1/2 -\beta) \Delta t^2 \\ + \textrm{b0} = (1 -\gamma) \Delta t \\ + \textrm{a1} = \beta \Delta t^2 \\ + \textrm{b1} = \gamma \Delta t \\ + Args: Minv (np.array): Inverse mass matrix :math:`\mathbf{M^{-1}}` C (np.array): Damping matrix :math:`\mathbf{C}` From d8359a0ffa2fb5b184ae1be81b0ec4f8551bd79d Mon Sep 17 00:00:00 2001 From: Bernardo Bahia Monteiro Date: Thu, 26 Oct 2023 22:06:11 -0400 Subject: [PATCH 05/13] Avoid factorizing Ass1 matrix multiple times --- sharpy/linear/src/lingebm.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/sharpy/linear/src/lingebm.py b/sharpy/linear/src/lingebm.py index bebfcfa44..84479950b 100644 --- a/sharpy/linear/src/lingebm.py +++ b/sharpy/linear/src/lingebm.py @@ -1508,10 +1508,12 @@ def newmark_ss(Minv, C, K, dt, num_damp=1e-4): [-b0 * MinvK, Imat - b0 * MinvC]]) Ass1 = np.block([[Imat + a1 * MinvK, a1 * MinvC], [b1 * MinvK, Imat + b1 * MinvC]]) - Ass = np.linalg.solve(Ass1, Ass0) - Bss0 = np.linalg.solve(Ass1, np.block([[a0 * Minv], [b0 * Minv]])) - Bss1 = np.linalg.solve(Ass1, np.block([[a1 * Minv], [b1 * Minv]])) + # Factorize Ass1 once, solve multiple times + Ass1_factors = sc.linalg.lu_factor(Ass1) + Ass = sc.linalg.lu_solve(Ass1_factors, Ass0) + Bss0 = sc.linalg.lu_solve(Ass1_factors, np.block([[a0 * Minv], [b0 * Minv]])) + Bss1 = sc.linalg.lu_solve(Ass1_factors, np.block([[a1 * Minv], [b1 * Minv]])) # eliminate predictior term Bss1 return libss.SSconv(Ass, Bss0, Bss1, C=np.eye(2 * N), D=np.zeros((2 * N, N))) From 3530a23a024168827408992cc7545891f0ad495d Mon Sep 17 00:00:00 2001 From: Bernardo Bahia Monteiro Date: Thu, 26 Oct 2023 22:24:06 -0400 Subject: [PATCH 06/13] Factor M inside newmark_ss to improve numerical stability and performance --- sharpy/linear/src/lingebm.py | 30 +++++++++++++++-------- tests/linear/rom/test_springmasssystem.py | 2 +- 2 files changed, 21 insertions(+), 11 deletions(-) diff --git a/sharpy/linear/src/lingebm.py b/sharpy/linear/src/lingebm.py index 84479950b..ffb180c68 100644 --- a/sharpy/linear/src/lingebm.py +++ b/sharpy/linear/src/lingebm.py @@ -152,7 +152,7 @@ def __init__(self, tsinfo, structure=None, custom_settings=dict()): # Store structure at linearisation and linearisation conditions self.structure = structure self.tsstruct0 = tsinfo - self.Minv = None + self.Minv = None # Not used anymore since M is factorized inside newmark_ss self.scaled_reference_matrices = dict() # keep reference values prior to time scaling @@ -932,7 +932,7 @@ def assemble(self, Nmodes=None): Ccut = np.dot(Phi.T, np.dot(self.Cstr, Phi)) Ass, Bss, Css, Dss = newmark_ss( - np.linalg.inv(np.dot(self.U[:, :Nmodes].T, np.dot(self.Mstr, self.U[:, :Nmodes]))), + np.dot(self.U[:, :Nmodes].T, np.dot(self.Mstr, self.U[:, :Nmodes])), Ccut, np.dot(self.U[:, :Nmodes].T, np.dot(self.Kstr, self.U[:, :Nmodes])), self.dt, @@ -983,10 +983,8 @@ def assemble(self, Nmodes=None): else: # Full system - self.Minv = np.linalg.inv(self.Mstr) - Ass, Bss, Css, Dss = newmark_ss( - self.Minv, self.Cstr, self.Kstr, + self.Mstr, self.Cstr, self.Kstr, self.dt, self.newmark_damp) self.Kin = None self.Kout = None @@ -1354,7 +1352,7 @@ def cont2disc(self, dt=None): self.dlti = True -def newmark_ss(Minv, C, K, dt, num_damp=1e-4): +def newmark_ss(M, C, K, dt, num_damp=1e-4, M_is_SPD=False): r""" Produces a discrete-time state-space model of the structural equations @@ -1473,11 +1471,12 @@ def newmark_ss(Minv, C, K, dt, num_damp=1e-4): \textrm{b1} = \gamma \Delta t \\ Args: - Minv (np.array): Inverse mass matrix :math:`\mathbf{M^{-1}}` + M (np.array): Mass matrix :math:`\mathbf{M}` C (np.array): Damping matrix :math:`\mathbf{C}` K (np.array): Stiffness matrix :math:`\mathbf{K}` dt (float): Timestep increment num_damp (float): Numerical damping. Default ``1e-4`` + M_is_SPD (bool): whether to factorized M using Cholesky (only works for SPD matrices) or LU decomposition. Default: False Returns: tuple: the A, B, C, D matrices of the state space packed in a tuple with the predictor and delay term removed. @@ -1497,11 +1496,22 @@ def newmark_ss(Minv, C, K, dt, num_damp=1e-4): b1 = th1 * dt b0 = dt - b1 - # relevant matrices + # Identity matrix N = K.shape[0] Imat = np.eye(N) - MinvK = np.dot(Minv, K) - MinvC = np.dot(Minv, C) + + # Factorize M and obtain relevant matrices + # Even though the inverse needs to be explicitly calculated, using the matrix + # factorization is more numerically stable and faster than multiplying by the + # inverse + M_factors = sc.linalg.cho_factor(M) if M_is_SPD else sc.linalg.lu_factor(M) + + def M_solve(b): + return sc.linalg.cho_solve(M_factors, b) if M_is_SPD else sc.linalg.lu_solve(M_factors, b) + + Minv = M_solve(Imat) + MinvK = M_solve(K) + MinvC = M_solve(C) # build StateSpace Ass0 = np.block([[Imat - a0 * MinvK, dt * Imat - a0 * MinvC], diff --git a/tests/linear/rom/test_springmasssystem.py b/tests/linear/rom/test_springmasssystem.py index 77adcd7bd..31d40b379 100644 --- a/tests/linear/rom/test_springmasssystem.py +++ b/tests/linear/rom/test_springmasssystem.py @@ -134,7 +134,7 @@ def build_system(self, system_inputs, system_time): else: # Discrete time system dt = 1e-2 - Adt, Bdt, Cdt, Ddt = lingebm.newmark_ss(Minv, C, k, dt=dt, num_damp=0) + Adt, Bdt, Cdt, Ddt = lingebm.newmark_ss(m, C, k, dt=dt, num_damp=0) system = libss.StateSpace(Adt, Bdt, Cdt, Ddt, dt=dt) From 0f85982cebeb6dfa0bc327b543234e0a6968a83f Mon Sep 17 00:00:00 2001 From: Bernardo Bahia Monteiro Date: Thu, 26 Oct 2023 22:48:38 -0400 Subject: [PATCH 07/13] Remove dead code and improve readability by using @ instead of np.dot This is a numpy best practice when multiplying matrices --- sharpy/linear/src/lingebm.py | 12 +++--------- 1 file changed, 3 insertions(+), 9 deletions(-) diff --git a/sharpy/linear/src/lingebm.py b/sharpy/linear/src/lingebm.py index ffb180c68..f801fa8fc 100644 --- a/sharpy/linear/src/lingebm.py +++ b/sharpy/linear/src/lingebm.py @@ -925,16 +925,10 @@ def assemble(self, Nmodes=None): if self.proj_modes == 'undamped': Phi = self.U[:, :Nmodes] - if self.Ccut is None: - # Ccut = np.zeros((Nmodes, Nmodes)) - Ccut = np.dot(Phi.T, np.dot(self.Cstr, Phi)) - else: - Ccut = np.dot(Phi.T, np.dot(self.Cstr, Phi)) - Ass, Bss, Css, Dss = newmark_ss( - np.dot(self.U[:, :Nmodes].T, np.dot(self.Mstr, self.U[:, :Nmodes])), - Ccut, - np.dot(self.U[:, :Nmodes].T, np.dot(self.Kstr, self.U[:, :Nmodes])), + Phi.T @ self.Mstr @ Phi, + Phi.T @ self.Cstr @ Phi, + Phi.T @ self.Kstr @ Phi, self.dt, self.newmark_damp) From 775f978e83111736d33e995045ec57a0caccd995 Mon Sep 17 00:00:00 2001 From: Bernardo Bahia Monteiro Date: Fri, 27 Oct 2023 15:17:42 -0400 Subject: [PATCH 08/13] [docs] Improve docs for newmark_ss --- sharpy/linear/src/lingebm.py | 78 ++++++++++++++++++++++-------------- 1 file changed, 48 insertions(+), 30 deletions(-) diff --git a/sharpy/linear/src/lingebm.py b/sharpy/linear/src/lingebm.py index f801fa8fc..fab4a2055 100644 --- a/sharpy/linear/src/lingebm.py +++ b/sharpy/linear/src/lingebm.py @@ -1412,57 +1412,75 @@ def newmark_ss(M, C, K, dt, num_damp=1e-4, M_is_SPD=False): \gamma \mathbf{\ddot q}_{n+1} \Delta t + O(\Delta t^3) Substituting the former relation onto the later ones, rearranging terms, and writing it in state-space form: - .. math:: + A_{ss1} \begin{Bmatrix} \mathbf{\dot q}_{n+1} \\ \mathbf{\ddot q}_{n+1} \end{Bmatrix} = + A_{ss0} \begin{Bmatrix} \mathbf{\dot q}_{n} \\ \mathbf{\ddot q}_{n} \end{Bmatrix} + + B_{ss0} F_n + + B_{ss1} F_{n+1} + where + .. math:: + A_{ss1} &= \begin{bmatrix} I + \beta\Delta t^2 M^{-1}K & \beta\Delta t^2 M^{-1}C \\ \gamma \Delta t M^{-1}K & I + \gamma \Delta t M^{-1}C \end{bmatrix} - \begin{Bmatrix} - \mathbf{q}_{n+1} \\ - \mathbf{\dot q}_{n+1} - \end{Bmatrix} - = + \\ + A_{ss0} &= \begin{bmatrix} I - \Delta t^2(1/2-\beta)M^{-1}K & \Delta t I - (1/2-\beta)\Delta t^2 M^{-1}C \\ -(1-\gamma)\Delta t M^{-1}K & I - (1-\gamma)\Delta t M^{-1}C \end{bmatrix} - \begin{Bmatrix} - \mathbf{q}_{n} \\ \ - mathbf{\dot q}_{n} - \end{Bmatrix} - + + \\ + B_{ss0} &= \begin{bmatrix} (\Delta t^2(1/2-\beta)\, M^{-1} \\ (1-\gamma)\Delta t \, M^{-1} \end{bmatrix} - F_n + + \\ + B_{ss1} &= \begin{bmatrix} (\beta \Delta t^2) \, M^{-1}\\ (\gamma \Delta t) \, M^{-1} \end{bmatrix} - F_{n+1} - Finally: - .. math:: - A_{ss1} \begin{Bmatrix} \mathbf{\dot q}_{n+1} \\ \mathbf{\ddot q}_{n+1} \end{Bmatrix} = - A_{ss0} \begin{Bmatrix} \mathbf{\dot q}_{n} \\ \mathbf{\ddot q}_{n} \end{Bmatrix} + - \begin{Bmatrix} (\Delta t^2(1/2-\beta) \\ (1-\gamma)\Delta t \end{Bmatrix} M^{-1}F_n+ - \begin{Bmatrix} (\Delta t^2\beta) \\ (\gamma \Delta t) \end{Bmatrix}M^{-1}F_{n+1} - - To finally isolate the vector at :math:`n+1`, instead of inverting the :math:`A_{ss1}` matrix, several systems are - solved. Moreover, the output equation is simply :math:`y=x`. + This is not in standard space-state form because the state update equation depends of the input at the + ..math:`n+1` timestep. This term can be eliminated by defining the state + ..math:: + X_n = + \begin{Bmatrix} + \mathbf{\dot q}_{n} \\ + \mathbf{\ddot q}_{n} + \end{Bmatrix} + B_{ss1} F_n - To understand SHARPy code, it is convenient to apply the following change of notation: - + Then + ..math:: + X_{n+1} &= A_{ss1}^{-1}(A_{ss0} X_n + (A_{ss0}B_{ss1} + B_{ss0}) F_n) \\ + \begin{Bmatrix} + \mathbf{\dot q}_{n} \\ + \mathbf{\ddot q}_{n} + \end{Bmatrix} + &= X_n - B_{ss1} F_n + + see also libss.conv for more details on the elimination of the term + multiplying ..math:`F_{n+1}` in the state equation. + + This function retuns a tuple with the discrete state-space matrices ..math:` (A,B,C,D)` where + ..math:: + A &= A_{ss1}^{-1}A_{ss0} \\ + B &= A_{ss1}^{-1}(B_{ss0} + A_{ss0}B{ss1}) \\ + C &= I \\ + D &= -B_{ss1} + + + The following notation is used in the code: .. math:: - \textrm{th1} = \gamma \\ - \textrm{th2} = \beta \\ - \textrm{a0} = (1/2 -\beta) \Delta t^2 \\ - \textrm{b0} = (1 -\gamma) \Delta t \\ - \textrm{a1} = \beta \Delta t^2 \\ - \textrm{b1} = \gamma \Delta t \\ + \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 \\ Args: M (np.array): Mass matrix :math:`\mathbf{M}` From ed20f3f86bc0fd3466be2638b36a86c9df83bfad Mon Sep 17 00:00:00 2001 From: Bernardo Bahia Monteiro Date: Fri, 27 Oct 2023 15:33:57 -0400 Subject: [PATCH 09/13] [docs] Improve wording and notation in newmark_ss --- sharpy/linear/src/lingebm.py | 25 ++++++++++++------------- 1 file changed, 12 insertions(+), 13 deletions(-) diff --git a/sharpy/linear/src/lingebm.py b/sharpy/linear/src/lingebm.py index fab4a2055..a899a4ce4 100644 --- a/sharpy/linear/src/lingebm.py +++ b/sharpy/linear/src/lingebm.py @@ -1348,28 +1348,22 @@ def cont2disc(self, dt=None): def newmark_ss(M, C, K, dt, num_damp=1e-4, M_is_SPD=False): r""" - Produces a discrete-time state-space model of the structural equations - + Produces a discrete-time state-space model of the 2nd order ordinary differential equation (ODE) given by .. math:: + \mathbf{M}\mathbf{\ddot{x}} + \mathbf{C}\mathbf{\dot{x}} + \mathbf{K}\mathbf{x} = \mathbf{F} - \mathbf{\ddot{x}} &= \mathbf{M}^{-1}( -\mathbf{C}\,\mathbf{\dot{x}}-\mathbf{K}\,\mathbf{x}+\mathbf{F} ) \\ - \mathbf{y} &= \mathbf{x} - - - based on the Newmark-:math:`\beta` integration scheme. The output state-space model - has form: + This ODE is discretized based on the Newmark-:math:`\beta` integration scheme. + The output state-space model has the form: .. math:: - \mathbf{X}_{n+1} &= \mathbf{A}\,\mathbf{X}_n + \mathbf{B}\,\mathbf{F}_n \\ \mathbf{Y} &= \mathbf{C}\,\mathbf{X} + \mathbf{D}\,\mathbf{F} - - with :math:`\mathbf{X} = [\mathbf{x}, \mathbf{\dot{x}}]^T` + with :math:`\mathbf{Y} = [\mathbf{x}, \mathbf{\dot{x}}]^T` Note that as the state-space representation only requires the input force - :math:`\mathbf{F}` to be evaluated at time-step :math:`n`,the :math:`\mathbf{C}` and :math:`\mathbf{D}` matrices - are, in general, fully populated. + :math:`\mathbf{F}` to be evaluated at time-step :math:`n`, thus the pass-through matrix + :math:`\mathbf{D}` is not zero. The Newmark-:math:`\beta` integration scheme is carried out following the modifications presented by Geradin [1] that render it unconditionally stable. The displacement and velocities are estimated as: @@ -1465,6 +1459,11 @@ def newmark_ss(M, C, K, dt, num_damp=1e-4, M_is_SPD=False): see also libss.conv for more details on the elimination of the term multiplying ..math:`F_{n+1}` in the state equation. + This system is identified with a standard discrete state space model + ..math:: + X_{n+1} = A X_n + B F_n + Y_n = C X_n + D F_n + This function retuns a tuple with the discrete state-space matrices ..math:` (A,B,C,D)` where ..math:: A &= A_{ss1}^{-1}A_{ss0} \\ From 5871ce22d0e78b08adfee2fcc1d83a34c1615ac8 Mon Sep 17 00:00:00 2001 From: Bernardo Bahia Monteiro Date: Fri, 27 Oct 2023 21:11:39 +0000 Subject: [PATCH 10/13] [docs] Use autofunction instead of automodule to print function signature --- docs/source/includes/linear/src/lingebm/newmark_ss.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/includes/linear/src/lingebm/newmark_ss.rst b/docs/source/includes/linear/src/lingebm/newmark_ss.rst index 3543f29c7..295cff4e6 100644 --- a/docs/source/includes/linear/src/lingebm/newmark_ss.rst +++ b/docs/source/includes/linear/src/lingebm/newmark_ss.rst @@ -1,4 +1,4 @@ newmark_ss ---------- -.. automodule:: sharpy.linear.src.lingebm.newmark_ss \ No newline at end of file +.. autofunction:: sharpy.linear.src.lingebm.newmark_ss From d68af9079ec1736c6818f5eff7d192cc3641ab81 Mon Sep 17 00:00:00 2001 From: Bernardo Bahia Monteiro Date: Fri, 27 Oct 2023 22:24:03 +0000 Subject: [PATCH 11/13] [docs] Uniformize notation and correct mistakes in docstring of newmark_ss --- sharpy/linear/src/lingebm.py | 186 ++++++++++++++++++++--------------- 1 file changed, 108 insertions(+), 78 deletions(-) diff --git a/sharpy/linear/src/lingebm.py b/sharpy/linear/src/lingebm.py index a899a4ce4..16ea81b06 100644 --- a/sharpy/linear/src/lingebm.py +++ b/sharpy/linear/src/lingebm.py @@ -1348,56 +1348,52 @@ def cont2disc(self, dt=None): def newmark_ss(M, C, K, dt, num_damp=1e-4, M_is_SPD=False): r""" - Produces a discrete-time state-space model of the 2nd order ordinary differential equation (ODE) given by + Produces a discrete-time state-space model of the 2nd order ordinary differential equation (ODE) given by: + .. math:: - \mathbf{M}\mathbf{\ddot{x}} + \mathbf{C}\mathbf{\dot{x}} + \mathbf{K}\mathbf{x} = \mathbf{F} + \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-:math:`\beta` integration scheme. The output state-space model has the form: + .. math:: - \mathbf{X}_{n+1} &= \mathbf{A}\,\mathbf{X}_n + \mathbf{B}\,\mathbf{F}_n \\ - \mathbf{Y} &= \mathbf{C}\,\mathbf{X} + \mathbf{D}\,\mathbf{F} + \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 - with :math:`\mathbf{Y} = [\mathbf{x}, \mathbf{\dot{x}}]^T` + where :math:`\mathbf{y} = \begin{Bmatrix} \mathbf{q} \\ \mathbf{\dot{q}} \end{Bmatrix}` Note that as the state-space representation only requires the input force - :math:`\mathbf{F}` to be evaluated at time-step :math:`n`, thus the pass-through matrix - :math:`\mathbf{D}` is not zero. - - The Newmark-:math:`\beta` integration scheme is carried out following the modifications presented by - Geradin [1] that render it unconditionally stable. The displacement and velocities are estimated as: - - .. math:: - x_{n+1} &= x_n + \Delta t \dot{x}_n + \left(\frac{1}{2}-\theta_2\right)\Delta t^2 \ddot{x}_n + \theta_2\Delta t - \ddot{x}_{n+1} \\ - \dot{x}_{n+1} &= \dot{x}_n + (1-\theta_1)\Delta t \ddot{x}_n + \theta_1\Delta t \ddot{x}_{n+1} - - The stencil is unconditionally stable if the tuning parameters :math:`\theta_1` and :math:`\theta_2` are chosen as: + :math:`\mathbf{f}` to be evaluated at time-step :math:`n`, thus the pass-through matrix + :math:`\mathbf{D_{ss}}` is not zero. - .. math:: - \theta_1 &= \frac{1}{2} + \alpha \\ - \theta_2 &= \frac{1}{4} \left(\theta_1 + \frac{1}{2}\right)^2 \\ - \theta_2 &= \frac{5}{80} + \frac{1}{4} (\theta_1 + \theta_1^2) \text{TBC SOURCE} + This function retuns a tuple with the discrete state-space matrices :math:`(\mathbf{A_{ss},B_{ss},C_{ss},D_{ss}})`. - where :math:`\alpha>0` accounts for small positive algorithmic damping. + Theory + ------ - The following steps describe how to apply the Newmark-beta scheme to a state-space formulation. The original idea - is based on [1]. + The following steps describe how to apply the Newmark-:math:`\beta` scheme + to the ODE in order to generate the discrete time-state space-model. It + folows the development of [1]. - The equation of a second order dynamical system reads: + .. admonition:: Notation - .. math:: - M\mathbf{\ddot q} + C\mathbf{\dot q} + K\mathbf{q} = F + 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. - Applying that equation to the time steps :math:`n` and :math:`n+1`, rearranging terms and multiplying by - :math:`M^{-1}`: + Evaluating the ODE to the time steps :math:`t_n` and :math:`t_{n+1}` and + isolating the acceleration term: .. math:: - \mathbf{\ddot q}_{n} = - M^{-1}C\mathbf{\dot q}_{n} - M^{-1}K\mathbf{q}_{n} + M^{-1}F_{n} \\ - \mathbf{\ddot q}_{n+1} = - M^{-1}C\mathbf{\dot q}_{n+1} - M^{-1}K\mathbf{q}_{n+1} + M^{-1}F_{n+1} + \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} \\ - The relations of the Newmark-beta scheme are: + The update equations of the Newmark-beta scheme are [1]: .. math:: \mathbf{q}_{n+1} &= \mathbf{q}_n + \mathbf{\dot q}_n\Delta t + @@ -1405,74 +1401,108 @@ def newmark_ss(M, C, K, dt, num_damp=1e-4, M_is_SPD=False): \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) - Substituting the former relation onto the later ones, rearranging terms, and writing it in state-space form: + where :math:`\Delta t = t_{n+1} - t_n`. + + The stencil is unconditionally stable if the tuning parameters + :math:`\gamma` and :math:`\beta` are chosen as: + + .. math:: + \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) + + where :math:`\alpha>0` accounts for small positive algorithmic damping (:math:`\alpha` is ``num_damp`` in the code). + + Substituting the former relations onto the later ones, rearranging terms, and writing it in state-space form: + .. math:: - A_{ss1} \begin{Bmatrix} \mathbf{\dot q}_{n+1} \\ \mathbf{\ddot q}_{n+1} \end{Bmatrix} = - A_{ss0} \begin{Bmatrix} \mathbf{\dot q}_{n} \\ \mathbf{\ddot q}_{n} \end{Bmatrix} + - B_{ss0} F_n + - B_{ss1} F_{n+1} + \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} + where + .. math:: - A_{ss1} &= + \mathbf{A_{ss1}} &= \begin{bmatrix} - I + \beta\Delta t^2 M^{-1}K & \beta\Delta t^2 M^{-1}C \\ - \gamma \Delta t M^{-1}K & I + \gamma \Delta t M^{-1}C + \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} \\ - A_{ss0} &= + \mathbf{A_{ss0}} &= \begin{bmatrix} - I - \Delta t^2(1/2-\beta)M^{-1}K & \Delta t I - (1/2-\beta)\Delta t^2 M^{-1}C \\ - -(1-\gamma)\Delta t M^{-1}K & I - (1-\gamma)\Delta t M^{-1}C + \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} \\ - B_{ss0} &= + \mathbf{B_{ss0}} &= \begin{bmatrix} - (\Delta t^2(1/2-\beta)\, M^{-1} \\ - (1-\gamma)\Delta t \, M^{-1} + (\Delta t^2(1/2-\beta) \mathbf{M}^{-1} \\ + (1-\gamma)\Delta t \mathbf{M}^{-1} \end{bmatrix} \\ - B_{ss1} &= + \mathbf{B_{ss1}} &= \begin{bmatrix} - (\beta \Delta t^2) \, M^{-1}\\ - (\gamma \Delta t) \, M^{-1} + (\beta \Delta t^2) \mathbf{M}^{-1}\\ + (\gamma \Delta t) \mathbf{M}^{-1} \end{bmatrix} - This is not in standard space-state form because the state update equation depends of the input at the - ..math:`n+1` timestep. This term can be eliminated by defining the state - ..math:: - X_n = - \begin{Bmatrix} - \mathbf{\dot q}_{n} \\ - \mathbf{\ddot q}_{n} - \end{Bmatrix} + B_{ss1} F_n + This is not in standard space-state form because the state update equation + depends of the input at :math:`t_{n+1}`. This term can be eliminated by + defining the state + + .. math:: \mathbf{x}_n = + \begin{Bmatrix} + \mathbf{q}_{n} \\ + \mathbf{\dot q}_{n} + \end{Bmatrix} + - \mathbf{A_{ss1}}^{-1}\mathbf{B_{ss1}} \mathbf{f}_n Then - ..math:: - X_{n+1} &= A_{ss1}^{-1}(A_{ss0} X_n + (A_{ss0}B_{ss1} + B_{ss0}) F_n) \\ + + .. math:: + \mathbf{x}_{n+1} &= \mathbf{A_{ss1}}^{-1}[ + \mathbf{A_{ss0}} \mathbf{x}_n + ( + \mathbf{A_{ss0}}\mathbf{A_{ss1}}^{-1}\mathbf{B_{ss1}} + + \mathbf{B_{ss0}} + ) + \mathbf{f}_n + ] \\ \begin{Bmatrix} \mathbf{\dot q}_{n} \\ \mathbf{\ddot q}_{n} \end{Bmatrix} - &= X_n - B_{ss1} F_n - - see also libss.conv for more details on the elimination of the term - multiplying ..math:`F_{n+1}` in the state equation. - - This system is identified with a standard discrete state space model - ..math:: - X_{n+1} = A X_n + B F_n - Y_n = C X_n + D F_n - - This function retuns a tuple with the discrete state-space matrices ..math:` (A,B,C,D)` where - ..math:: - A &= A_{ss1}^{-1}A_{ss0} \\ - B &= A_{ss1}^{-1}(B_{ss0} + A_{ss0}B{ss1}) \\ - C &= I \\ - D &= -B_{ss1} + &= \mathbf{x}_n + \mathbf{B_{ss1}} \mathbf{f}_n + + See also :func:`sharpy.linear.src.libss.SSconv` for more details on the elimination of the term + multiplying :math:`\mathbf{f}_{n+1}` in the state equation. + + This system is identified with a standard discrete-time state-space model + + .. math:: + \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 + + where + + .. math:: + \mathbf{A_{ss}} &= \mathbf{A_{ss1}}^{-1}\mathbf{A_{ss0}} \\ + \mathbf{B_{ss}} &= \mathbf{A_{ss1}}^{-1}(\mathbf{B_{ss0}} + + \mathbf{A_{ss0}}\mathbf{A_{ss1}}^{-1}\mathbf{B_{ss1}}) \\ + \mathbf{C_{ss}} &= \mathbf{I} \\ + \mathbf{D_{ss}} &= \mathbf{B_{ss1}} - The following notation is used in the code: + .. admonition:: Notation is used in the code + .. math:: \texttt{th1} &= \gamma \\ \texttt{th2} &= \beta \\ @@ -1487,10 +1517,10 @@ def newmark_ss(M, C, K, dt, num_damp=1e-4, M_is_SPD=False): K (np.array): Stiffness matrix :math:`\mathbf{K}` dt (float): Timestep increment num_damp (float): Numerical damping. Default ``1e-4`` - M_is_SPD (bool): whether to factorized M using Cholesky (only works for SPD matrices) or LU decomposition. Default: False + M_is_SPD (bool): whether to factorized M using Cholesky (only works for SPD matrices) or LU decomposition. Default: ``False`` Returns: - tuple: the A, B, C, D matrices of the state space packed in a tuple with the predictor and delay term removed. + tuple: with the :math:`(\mathbf{A_{ss},B_{ss},C_{ss},D_{ss}})` matrices of the discrete-time state-space representation References: [1] - Geradin M., Rixen D. - Mechanical Vibrations: Theory and application to structural dynamics From d6f0ffe26ad58b3bfa0c7f0e28fb373e10e7f2b9 Mon Sep 17 00:00:00 2001 From: Bernardo Bahia Monteiro Date: Fri, 27 Oct 2023 22:48:18 +0000 Subject: [PATCH 12/13] [docs] Suggest >16GB memory for Docker installation --- docs/source/content/installation.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/source/content/installation.md b/docs/source/content/installation.md index 72912462c..90415ab54 100644 --- a/docs/source/content/installation.md +++ b/docs/source/content/installation.md @@ -207,6 +207,10 @@ to your taste. ## Using SHARPy from a Docker container +> **Tip** To install the Python environment, miniconda needs approximatelly 16GB of +> memory. It is recommended to have at least that amount available for the +> container in which you are building SHARPy. + Docker containers are similar to lightweight virtual machines. The SHARPy container distributed through [Docker Hub](https://hub.docker.com/) is a CentOS 8 machine with the libraries compiled with `gfortran` and `g++` and an From db53dcfd85389b6cb8491ad63b7c2a5b4cd9b7aa Mon Sep 17 00:00:00 2001 From: Bernardo Bahia Monteiro Date: Fri, 27 Oct 2023 22:48:56 +0000 Subject: [PATCH 13/13] [docs] Create section for testing in Docker so it is easier to find --- docs/source/content/installation.md | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/source/content/installation.md b/docs/source/content/installation.md index 90415ab54..412c9f736 100644 --- a/docs/source/content/installation.md +++ b/docs/source/content/installation.md @@ -272,14 +272,16 @@ docker cp sharpy:/my_file.txt my_file.txt # copy from container to host ``` The `sharpy:` part is the `--name` argument you wrote in the `docker run` command. +**Enjoy!** + +## Testing with Docker + You can run the test suite once inside the container as: ``` cd sharpy_dir python -m unittest ``` -**Enjoy!** - ## Running SHARPy