System Manipulation, Reduction, and Correction (system)¶
System manipulation, reduction, and corrections functions.
@author: Joseph C. Slater
-
system.c2d(A, B, C, D, dt)[source]¶ Convert continuous state system to discrete time.
Converts a set of continuous state space system matrices to their discrete counterpart.
Parameters: A, B, C, D : float arrays
State space system matrices
dt : float
Time step of discrete system
Returns: Ad, Bd, C, D : float arrays
Discrete state space system matrices
Notes
Note
Zero-order hold solution
See also
Examples
>>> import numpy as np >>> np.set_printoptions(precision=4, suppress=True) >>> import vibrationtesting as vt >>> A1 = np.array([[ 0., 0. , 1. , 0. ]]) >>> A2 = np.array([[ 0., 0. , 0. , 1. ]]) >>> A3 = np.array([[-1.4, 1.2, -0.0058, 0.0014]]) >>> A4 = np.array([[ 0.8, -1.4, 0.0016, -0.0038]]) >>> A = np.array([[ 0., 0. , 1. , 0. ], ... [ 0., 0. , 0. , 1. ], ... [-1.4, 1.2, -0.0058, 0.0014], ... [ 0.8, -1.4, 0.0016, -0.0038]]) >>> B = np.array([[ 0.], ... [ 0.], ... [ 0.], ... [ 1.]]) >>> C = np.array([[-1.4, 1.2, -0.0058, 0.0014]]) >>> D = np.array([[-0.2]]) >>> Ad, Bd, *_ = vt.c2d(A, B, C, D, 0.01) >>> print(Ad) [[ 0.9999 0.0001 0.01 0. ] [ 0. 0.9999 0. 0.01 ] [-0.014 0.012 0.9999 0.0001] [ 0.008 -0.014 0.0001 0.9999]] >>> print(Bd) [[ 0. ] [ 0. ] [ 0. ] [ 0.01]]
-
system.d2c(Ad, Bd, C, D, dt)[source]¶ Return continuous A, B, C, D from discrete form.
Converts a set of digital state space system matrices to their continuous counterpart.
Parameters: Ad, Bd, C, D : float arrays
Discrete state space system matrices
dt : float
Time step of discrete system
Returns: A, B, C, D : float arrays
State space system matrices
Notes
Note
Zero-order hold solution
See also
Examples
>>> import vibrationtesting as vt >>> Ad = np.array([[ 0.9999,0.0001,0.01,0.], ... [ 0.,0.9999,0.,0.01], ... [-0.014,0.012,0.9999,0.0001], ... [ 0.008, -0.014,0.0001,0.9999]]) >>> Bd = np.array([[ 0. ], ... [ 0. ], ... [ 0. ], ... [ 0.01]]) >>> C = np.array([[-1.4, 1.2, -0.0058, 0.0014]]) >>> D = np.array([[-0.2]]) >>> A, B, *_ = vt.d2c(Ad, Bd, C, D, 0.01) >>> print(A) [[-0.003 0.004 1.0001 -0.0001] [-0.004 -0.003 -0. 1.0001] [-1.4001 1.2001 -0.003 0.004 ] [ 0.8001 -1.4001 0.006 -0.003 ]]
-
system.guyan(M, K, master=None, fraction=None)[source]¶ Guyan reduced model.
Applies Guyan Reductions to second order system of equations of the form
\[M \ddot{x} + K x = 0\]which are reduced to the form
\[M_r \ddot{x}_m + K_r x_m = 0\]where
\[ \begin{align}\begin{aligned}x = T x_m\\M_r= T^T M T\\K_r= T^T K T\end{aligned}\end{align} \]Parameters: M, K : float arrays
Mass and Stiffness matrices
master : float array or list, optional
List of retained degrees of freedom (0 indexing)
fraction : float, optional
Fraction of degrees of freedom (0< fraction <1.0) to retain in model. If both master and fraction are neglected, fraction is set to 0.25.
Returns: Mred, Kred, T : float arrays
Reduced Mass matric, Reduced Stiffness matrix, Transformation matrix
master_dofs : int list
List of master degrees of freedom (0 indexing)
truncated_dofs : int list
List of truncated degrees of freedom (0 indexing)
Notes
Reduced coordinate system forces can be obtained by Fr = T.T @ F.
Reduced damping matrix can be obtained using Cr = T.T @ C @ T.
If mode shapes are obtained for the reduced system, full system mode shapes are phi = T @ phi_r.
Code is not as efficient as possible. Using block submatrices would be more efficient.
Examples
>>> import vibrationtesting as vt >>> M = np.array([[4, 0, 0], ... [0, 4, 0], ... [0, 0, 4]]) >>> K = np.array([[8, -4, 0], ... [-4, 8, -4], ... [0, -4, 4]]) >>> Mred, Kred, T, master, truncated_dofs = vt.guyan(M, K, fraction = 0.5)
-
system.mode_expansion_from_model(Psi, omega, M, K, measured)[source]¶ Deflection extrapolation to full FEM model coordinates, matrix method.
Provided an equation of the form:
\(\begin{pmatrix}-\begin{bmatrix}M_{mm}&M_{mu}\\ M_{um}&M_{uu} \end{bmatrix} \omega_i^2 +\begin{bmatrix}K_{mm}&K_{mu}\\ K_{um}&K_{uu}\end{bmatrix}\end{pmatrix}\) \(\begin{bmatrix}\Psi_{i_m}\\ \Psi_{i_u}\end{bmatrix}= 0\)
Where:
- \(M\) and \(K\) are the mass and stiffness matrices, likely from a finite element model
- \(\Psi_i\) and \(\omega_i\) represent a mode/frequency pair
- subscripts \(m\) and \(u\) represent measure and unmeasured of the mode
Determines the unknown portion of the mode (or operational deflection) shape, \(\Psi_{i_u}\) by direct algebraic solution, aka
\(\Psi_{i_u} = - (K_{uu}- M_{ss} \omega_i^2) ^{-1} (K_{um}-M_{um}\omega_i^2)\Psi_{i_m}\)
Parameters: Psi : float array
mode shape or shapes, 2-D array columns of which are mode shapes
omega : float or 1-D float array
natural (or driving) frequencies
M, K : float arrays
Mass and Stiffness matrices
measured : float or integer array or list
List of measured degrees of freedom (0 indexed)
Returns: Psi_full : float array
Complete mode shape
Notes
See also
incomplete multi-mode update. Would require each at a different frequency.
Examples
>>> import vibrationtesting as vt >>> M = np.array([[4, 0, 0], ... [0, 4, 0], ... [0, 0, 4]]) >>> K = np.array([[8, -4, 0], ... [-4, 8, -4], ... [0, -4, 4]]) >>> measured = np.array([[0, 2]]) >>> omega, zeta, Psi = vt.sos_modal(M, K) >>> Psi_measured = np.array([[-0.15], [-0.37]]) >>> Psi_full = vt.mode_expansion_from_model(Psi_measured, omega[0], M, K, ... measured) >>> print(np.hstack((Psi[:,0].reshape(-1,1), Psi_full))) [[-0.164 -0.15 ] [-0.2955 0.2886] [-0.3685 -0.37 ]]
-
system.model_correction_direct(Psi, omega, M, K, method='Baruch')[source]¶ Direct model updating using model data.
Parameters: Psi : float array
Expanded mode shapes from experimental measurements. Must be real-valued. See mode_expansion_from_model.
omega : float array
Natural frequencies identified from modal analysis (diagonal matrix or vector).
M, K : float arrays
Analytical mass and stiffness matrices
method : string, optional
Returns: Mc, Kc : float arrays
Corrected mass and stiffness matrices.
Notes
[R12] Baruch, M. and Bar-Itzhack, I.Y., “Optimal Weighted Orthogonalization of Measured Modes,” AIAA Journal, 16(4), 1978, pp. 346-351. [R22] Berman, A. and Nagy, E.J., 1983, “Improvements of a Large Analytical Model using Test Data,” AIAA Journal, 21(8), 1983, pp. 1168-1173.
-
system.real_modes(Psi, autorotate=True)[source]¶ Real modes from complex modes.
Assuming a transformation
\[Psi_{real} = Psi_{complex} T\]exists, where \(T\) is a complex transformation matrix, find \(Psi_{real}\) using linear algebra.
Parameters: Psi : complex float array
Complex mode shapes (displacement)
autorotate : Boolean, optional
Attempt to rotate to near-real first
Returns: Psi : float array
Real modes
Notes
Note
Rotation of modes should be performed to get them as close to real as possible first.
Warning
Current autorotate bases the rotation on de-rotating the first element of each vector. User can use their own pre-process by doing to and setting autorotate to False.
Examples
>>> import vibrationtesting as vt >>> M = np.array([[4, 0, 0], ... [0, 4, 0], ... [0, 0, 4]]) >>> Cso = np.array([[.1,0,0], ... [0,0,0], ... [0,0,0]]) >>> K = np.array([[8, -4, 0], ... [-4, 8, -4], ... [0, -4, 4]]) >>> Bt = np.array([[1],[0],[0]]) >>> Ca = np.array([[1,0,0]]) >>> Cd = Cv = np.zeros_like(Ca) >>> A, B, C, D = vt.so2ss(M, Cso, K, Bt, Cd, Cv, Ca) >>> Am, Bm, Cm, Dm, eigenvalues, modes = vt.ss_modal(A, B, C, D) >>> Psi = vt.real_modes(modes[:,0::2])
-
system.rsolve(B, C, **kwargs)[source]¶ Solve right Gauss elimination equation.
Given \(A B = C\) return \(A = C B^{-1}\)
This uses scipy.linalg.solve with a little matrix manipulation first. All keyword arguments of scipy.linalg.solve may be used.
Parameters: B, C : float arrays Returns: A : float array Notes
See also
scipy.linalg.solve
Examples
>>> import numpy as np >>> import vibrationtesting as vt >>> B = np.array([[ 8, -4, 0], ... [-4, 8, -4], ... [ 0, -4, 4]]) >>> C = np.array([[ 32, -16, 0], ... [-16, 36, -20], ... [ 4, -24, 20]]) >>> A = vt.rsolve(B, C) >>> print(np.round(rsolve(B, C))) [[ 4. 0. 0.] [-0. 4. -1.] [ 0. -1. 4.]]
-
system.serep(M, K, master)[source]¶ System Equivalent Reduction Expansion Process reduced model.
Reduce size of second order system of equations by SEREP processs while returning expansion matrix
Equation of the form: \(M \ddot{x} + K x = 0\) is reduced to the form \(M_r \ddot{x}_m + Kr x_m = 0\) where \(x = T x_m\), \(M_r= T^T M T\), \(K_r= T^T K T\)
Parameters: M, K : float arrays
Mass and Stiffness matrices
master : float array or list
List of retained degrees of freedom
Returns: Mred, Kred, T : float arrays
Reduced Mass matric, reduced stiffness matrix, Transformation matrix
truncated_dofs : int list
List of truncated degrees of freedom, zero indexed
Notes
Reduced coordinate system forces can be obtained by Fr = T.T @ F
Reduced damping matrix can be obtained using Cr = T.T*@ C @ T.
If mode shapes are obtained for the reduced system, full system mode shapes are phi = T @ phi_r
See also
Examples
>>> import vibrationtesting as vt >>> M = np.array([[4, 0, 0], ... [0, 4, 0], ... [0, 0, 4]]) >>> K = np.array([[8, -4, 0], ... [-4, 8, -4], ... [0, -4, 4]]) >>> retained = np.array([[1, 2]]) >>> Mred, Kred, T, truncated_dofs = vt.serep(M, K, retained)
-
system.slice(Matrix, a, b)[source]¶ Slice a matrix properly- like Octave.
Addresses the confounding inconsistency that M[a,b] acts differently if a and b are the same length or different lengths.
Parameters: Matrix : float array
Arbitrary array
a, b : int lists or arrays
list of rows and columns to be selected from Matrix
Returns: Matrix : float array
Properly sliced matrix- no casting allowed.
-
system.so2ss(M, C, K, Bt, Cd, Cv, Ca)[source]¶ Convert second order system to state space.
Given second order linear matrix equation of the form \(M\\ddot{x} + C \\dot{x} + K x= \\tilde{B} u\) and \(y = C_d x + C_v \\dot{x} + C_a\\ddot{x}\) returns the state space form equations \(\\dot{z} = A z + B u\), \(y = C z + D u\)
Parameters: M, C, K, Bt, Cd, Cv, Cd : float arrays
Mass , damping, stiffness, input, displacement sensor, velocimeter, and accelerometer matrices
Returns: A, B, C, D : float arrays
State matrices
Examples
>>> import vibrationtesting as vt >>> M = np.array([[2, 1], ... [1, 3]]) >>> K = np.array([[2, -1], ... [-1, 3]]) >>> C = np.array([[0.01, 0.001], ... [0.001, 0.01]]) >>> Bt = np.array([[0], [1]]) >>> Cd = Cv = np.zeros((1,2)) >>> Ca = np.array([[1, 0]]) >>> A, B, Css, D = vt.so2ss(M, C, K, Bt, Cd, Cv, Ca) >>> print('A: {}'.format(A)) A: [[ 0. 0. 1. 0. ] [ 0. 0. 0. 1. ] [-1.4 1.2 -0.0058 0.0014] [ 0.8 -1.4 0.0016 -0.0038]] >>> print('B: ', B) B: [[ 0. ] [ 0. ] [-0.2] [ 0.4]] >>> print('C: ', Css) C: [[-1.4 1.2 -0.0058 0.0014]] >>> print('D: ', D) D: [[-0.2]]
-
system.sos_frf(M, C, K, Bt, Cd, Cv, Ca, omega_low, omega_high, in_index, out_index, num_freqs=1000)[source]¶ Return FRF of second order system.
Given second order linear matrix equation of the form \(M\\ddot{x} + C \\dot{x} + K x= \\tilde{B} u\) and \(y = C_d x + C_v \\dot{x} + C_a\\ddot{x}\) converts to state space form and returns the requested frequency response function
Parameters: M, C, K, Bt, Cd, Cv, Cd : float arrays
Mass, damping, stiffness, input, displacement sensor, velocimeter, and accelerometer matrices
omega_low, omega_high : floats
low and high frequencies for evaluation
in_index, out_index : ints
input and output numbers (starting at 1)
num_freqs : int
number of frequencies at which to return FRF
Returns: omega : float array
frequency vector
H : float array
frequency response function
Examples not working for second order system
Need to make one for second order expansion
Examples
>>> import vibrationtesting as vt >>> A = np.array([[ 0., 0. , 1. , 0. ], ... [ 0., 0. , 0. , 1. ], ... [-1.4, 1.2, -0.0058, 0.0014], ... [ 0.8, -1.4, 0.0016, -0.0038]]) >>> B = np.array([[ 0.], ... [ 0.], ... [ 0.], ... [ 1.]]) >>> C = np.array([[-1.4, 1.2, -0.0058, 0.0014]]) >>> D = np.array([[-0.2]]) >>> omega, H = vt.ssfrf(A, B, C, D, 0, 3.5, 1, 1) >>> vt.frfplot(omega, H)
-
system.sos_modal(M, K, C=False, damp_diag=0.03, shift=1)[source]¶ Eigen analysis of proportionally damped system.
Optimally find mass normalized mode shapes and natural frequencies of a system modelled by \(M\ddot{x}+Kx=0\).
If provided, obtain damping ratios presuming \(C\) can be decoupled.
Provides a warning if diagonalization of damping matrix fails worse than relative error of damp_diag.
Parameters: M, K : float arrays
Mass and stiffness matrices
C : float array, optional
Damping matrix
damp_diag : float, optional
Maximum amount of off-diagonal error allowed in assuming C can be diagonalized
shift : float, optional
Shift used in eigensolution. Should be approximately equal to the first non-zero eigenvalue.
Returns: omega : float array (1xN)
Vector of natural frequencies (rad/sec)
zeta : float array (1xN)
Vector of damping ratios
Psi : float array (NxN)
Matrix of mass normalized mode shapes by column
Examples
>>> import numpy as np >>> import vibrationtesting as vt >>> M = np.array([[4, 0, 0], ... [0, 4, 0], ... [0, 0, 4]]) >>> K = np.array([[8, -4, 0], ... [-4, 8, -4], ... [0, -4, 4]]) >>> omega, zeta, Psi = vt.sos_modal(M, K, K/10) >>> print(omega) [ 0.445 1.247 1.8019] >>> print(Psi.T@K@Psi) [[ 0.1981 0. -0. ] [ 0. 1.555 -0. ] [-0. -0. 3.247 ]]
Check that it works for rigid body modes.
>>> K2 = K-np.eye(K.shape[0])@M*(Psi.T@K@Psi)[0,0] >>> omega, zeta, Psi = vt.sos_modal(M, K2) >>> print(omega) [ 0. 1.1649 1.7461] >>> print(Psi) [[-0.164 0.3685 -0.2955] [-0.2955 0.164 0.3685] [-0.3685 -0.2955 -0.164 ]] >>> print(np.diag(Psi.T@K2@Psi)) [ 0. 1.3569 3.0489]
How about non-proportional damping
>>> C = K/10 >>> C[0,0] = 2 * C[0,0] >>> omega, zeta, Psi = vt.sos_modal(M, K2, C) Damping matrix cannot be completely diagonalized. Off diagonal error of 22%. >>> print(omega) [ 0. 1.1649 1.7461] >>> print(zeta) [ 0. 0.1134 0.113 ] >>> print(Psi.T@C@Psi) [[ 0.0413 -0.0483 0.0388] [-0.0483 0.2641 -0.0871] [ 0.0388 -0.0871 0.3946]]
-
system.ss_modal(A, B=None, C=None, D=None)[source]¶ State space modes, frequencies, damping ratios, and modal matrices.
Parameters: A, B, C, D : float arrays
State space system matrices
Returns: Am, Bm, Cm, Dm : float arrays
Modal state space system matrices
Examples
>>> import vibrationtesting as vt >>> M = np.array([[4, 0, 0], ... [0, 4, 0], ... [0, 0, 4]]) >>> Cso = np.array([[.1,0,0], ... [0,0,0], ... [0,0,0]]) >>> K = np.array([[8, -4, 0], ... [-4, 8, -4], ... [0, -4, 4]]) >>> Bt = np.array([[1],[0],[0]]) >>> Ca = np.array([[1,0,0]]) >>> Cd = Cv = np.zeros_like(Ca) >>> A, B, C, D = vt.so2ss(M, Cso, K, Bt, Cd, Cv, Ca) >>> Am, Bm, Cm, Dm, eigenvalues, modes = vt.ss_modal(A, B, C, D) >>> print(Am+(0.00001+0.00001j)) [[-0.0013+0.4451j 0.0000+0.j 0.0000+0.j 0.0000+0.j 0.0000+0.j 0.0000+0.j ] [ 0.0000+0.j -0.0013-0.445j 0.0000+0.j 0.0000+0.j 0.0000+0.j 0.0000+0.j ] [ 0.0000+0.j 0.0000+0.j -0.0068+1.247j 0.0000+0.j 0.0000+0.j 0.0000+0.j ] [ 0.0000+0.j 0.0000+0.j 0.0000+0.j -0.0068-1.247j 0.0000+0.j 0.0000+0.j ] [ 0.0000+0.j 0.0000+0.j 0.0000+0.j 0.0000+0.j -0.0044+1.8019j 0.0000+0.j ] [ 0.0000+0.j 0.0000+0.j 0.0000+0.j 0.0000+0.j 0.0000+0.j -0.0044-1.8019j]] >>> print(Cm) [[ 0.0594-0.0001j 0.0594+0.0001j 0.0039-0.717j 0.0039+0.717j 0.0241-0.9307j 0.0241+0.9307j]]
-
system.ssfrf(A, B, C, D, omega_low, omega_high, in_index, out_index, num_freqs=1000)[source]¶ Return FRF of state space system.
Obtains the computed FRF of a state space system between selected input and output over frequency range of interest.
Parameters: A, B, C, D : float arrays
state system matrices
omega_low, omega_high : floats
low and high frequencies for evaluation
in_index, out_index : ints
input and output numbers (starting at 1)
num_freqs : int
number of frequencies at which to return FRF
Returns: omega : float array
frequency vector
H : float array
frequency response function
Examples
>>> import vibrationtesting as vt >>> A = np.array([[ 0., 0. , 1. , 0. ], ... [ 0., 0. , 0. , 1. ], ... [-1.4, 1.2, -0.0058, 0.0014], ... [ 0.8, -1.4, 0.0016, -0.0038]]) >>> B = np.array([[ 0.], ... [ 0.], ... [ 0.], ... [ 1.]]) >>> C = np.array([[-1.4, 1.2, -0.0058, 0.0014]]) >>> D = np.array([[-0.2]]) >>> omega, H = vt.ssfrf(A, B, C, D, 0, 3.5, 1, 1) >>> vt.frfplot(omega, H)