Mode shapes using ERA
Closed this issue · 5 comments
Dear @twmacro
Thanks for this brilliant package.
Is it possible to obtain the final mode shapes using Eigensystem Realization Algorithm?
Best regards
Ayubirad
Hi @Ayubirad,
Thank you for asking! The answer is yes, but I'm realizing now that I should polish it up a bit. I'm working on the ERA algorithm now, so I'll add this to the update. I'm guessing that will be done within a couple weeks or so.
The "C_modal" attribute is the final complex state-space mode-shape matrix. The rows correspond to the measurements. To get the final real mode-shapes currently, you'll need to do a small manipulation of C_modal:
phi_era = era.ERA(...).C_modal[:, ::2].real # assumes underdamped modes
Here is a complete example:
import numpy as np
import scipy.linalg as la
from pyyeti import era, ode
M = np.identity(3)
K = np.array(
[
[4185.1498, 576.6947, 3646.8923],
[576.6947, 2104.9252, -28.0450],
[3646.8923, -28.0450, 3451.5583],
]
)
D = np.array(
[
[4.96765646, 0.97182432, 4.0162425],
[0.97182432, 6.71403672, -0.86138453],
[4.0162425, -0.86138453, 4.28850828],
]
)
(w2, phi) = la.eigh(K, M)
omega = np.sqrt(w2)
freq_hz = omega / (2 * np.pi)
modal_damping = phi.T @ D @ phi
zeta = np.diag(modal_damping) / (2 * omega)
dt = 0.01
t = np.arange(0, 1, dt)
F = np.zeros((3, len(t)))
ts = ode.SolveUnc(M, D, K, dt)
sol = ts.tsolve(force=F, v0=[1, 1, 1])
era_fit = era.ERA(
sol.v,
sr=1 / dt,
auto=True,
input_labels=["x", "y", "z"],
)
phi_era = era_fit.C_modal[:, ::2].real
print("\nMode shape ratio era/known =\n", phi_era / phi)
Since there is no noise added, the mode shapes "phi_era" are exactly correct, but with arbitrary scale. The output of the above code is:
Current fit includes all modes:
Mode Freq. (Hz) Zeta MSV EMAC MPC MAC
--------------------------------------------------------------------
* 1 1.33603 0.02000 *0.817 *1.000 *1.000 *1.000
* 2 7.39083 0.07500 *0.839 *1.000 *1.000 *1.000
* 3 13.79671 0.05000 *1.000 *1.000 *1.000 *1.000
Auto-selected modes fit:
Mode Freq. (Hz) Zeta MSV EMAC MPC MAC
--------------------------------------------------------------------
1 1.33603 0.02000 0.817 1.000 1.000 1.000
2 7.39083 0.07500 0.839 1.000 1.000 1.000
3 13.79671 0.05000 1.000 1.000 1.000 1.000
Mode shape ratio era/known =
[[-0.05585 0.64635 0.84434]
[-0.05585 0.64635 0.84434]
[-0.05585 0.64635 0.84434]]
I hope this helps! :-)
Regards,
Tim
Hi @Ayubirad,
After thinking about the above a bit more, I believe the above calculation for "phi_era" is true only for cases where the damping is diagonalized by the real (undamped) mode shapes (phi_era). But, for lightly damped structures, the equation is probably a decent approximation. In any case, "C_modal" are the complex mode shapes as computed by ERA.
Regards,
Tim
Hi @Ayubirad,
Try #3! :-) I realized my approach above isn't robust ... if the real part happens to be small, we could be looking at small noisy numbers. Here is a better approach and the method I plan to incorporate in the next update. It's very simple: scale each complex mode so the abs-max value is 1, then take the real part:
C_modal = era_fit.C_modal[:, ::2] # assumes underdamped modes
rows = np.abs(C_modal).argmax(axis=0)
cols = np.arange(C_modal.shape[1])
phi_era = (C_modal / C_modal[rows, cols]).real
We can even calculate a "realness" parameter for each mode shape (shown below). For the example above, the modes are real, so the realness parameter would be 1 for all modes. But, if the damping D is modified so it is not diagonalized by the real modes, the realness parameter drops. Here is an updated example with a modified D (notice also this is using a partially updated ERA ... it has "EMACO" and "CMIO"):
import numpy as np
import scipy.linalg as la
from pyyeti import era, ode
M = np.identity(3)
K = np.array(
[
[4185.1498, 576.6947, 3646.8923],
[576.6947, 2104.9252, -28.0450],
[3646.8923, -28.0450, 3451.5583],
]
)
# D = np.array(
# [
# [4.96765646, 0.97182432, 4.0162425],
# [0.97182432, 6.71403672, -0.86138453],
# [4.0162425, -0.86138453, 4.28850828],
# ]
# )
# modified D:
D = np.array(
[
[4.96765646, 1.2, 4.0162425],
[1.2, 6.71403672, -0.86138453],
[4.0162425, -0.86138453, 4.28850828],
]
)
(w2, phi) = la.eigh(K, M)
omega = np.sqrt(w2)
freq_hz = omega / (2 * np.pi)
modal_damping = phi.T @ D @ phi
zeta = np.diag(modal_damping) / (2 * omega)
dt = 0.01
t = np.arange(0, 1, dt)
F = np.zeros((3, len(t)))
ts = ode.SolveUnc(M, D, K, dt)
sol = ts.tsolve(force=F, v0=[1, 1, 1])
era_fit = era.ERA(
sol.v,
sr=1 / dt,
auto=True,
input_labels=["x", "y", "z"],
)
C_modal = era_fit.C_modal[:, ::2]
rows = np.abs(C_modal).argmax(axis=0)
cols = np.arange(C_modal.shape[1])
phi_era_scaled = C_modal / C_modal[rows, cols]
phi_era = phi_era_scaled.real
realness = 1 - abs(phi_era_scaled.imag / phi_era).max(axis=0)
rat = phi_era / phi
print("\nMode shape ratio era/known =\n", rat / abs(rat).max(axis=0))
print(f"\n{realness = }")
# Note that since there was no noise added, ERA will recover the exact
# frequencies, damping, and complex modes. This can be confirmed by
# comparing these values to the output of `ode.eigss`:
print("\n")
print("Confirm exact answers from ERA:")
print("-------------------------------")
s = ode.eigss(ode.make_A(M, D, K), delcc=False)
print("Complex mode shape ratio =\n", s.ur[:3, ::2] / era_fit.C_modal[:, ::2])
print(f"Natural frequency ratio = {era_fit.freqs / s.wn}")
print(f"Damping ratio = {era_fit.zeta / s.zeta}")
And here is the output:
Current fit includes all modes:
Mode Freq (Hz) Zeta MSV EMAC EMACO MPC CMIO MAC
---- ---------- ---------- ------ ------ ------ ------ ------ ------
* 1 1.33604 0.01641 *0.824 *1.000 *1.000 *1.000 *1.000 *1.000
* 2 7.39081 0.07537 *0.838 *1.000 *1.000 *1.000 *1.000 *1.000
* 3 13.79667 0.05015 *1.000 *1.000 *1.000 *1.000 *1.000 *1.000
Auto-selected modes fit:
Mode Freq (Hz) Zeta MSV EMAC EMACO MPC CMIO MAC
---- ---------- ---------- ------ ------ ------ ------ ------ ------
1 1.33604 0.01641 0.824 1.000 1.000 1.000 1.000 1.000
2 7.39081 0.07537 0.838 1.000 1.000 1.000 1.000 1.000
3 13.79667 0.05015 1.000 1.000 1.000 1.000 1.000 1.000
Mode shape ratio era/known =
[[ 1. -1. -0.9999]
[ 1. -0.9971 -0.9992]
[ 1. -0.9977 -1. ]]
realness = array([ 0.997 , 0.9827, 0.9645])
Confirm exact answers from ERA:
-------------------------------
Complex mode shape ratio =
[[-0.4698-2.7875j -1.5419+0.0903j 1.1551-0.1847j]
[-0.4698-2.7875j -1.5419+0.0903j 1.1551-0.1847j]
[-0.4698-2.7875j -1.5419+0.0903j 1.1551-0.1847j]]
Natural frequency ratio = [ 1. 1. 1.]
Damping ratio = [ 1. 1. 1.]
Dear @twmacro
I wanted to express my sincerest gratitude for your exceptional work in clarifying the ERA procedure. Furthermore, I'm eagerly anticipating the updated package. I was wondering if there might be a possibility of obtaining early access to it, prior to its official release! 👀
Furthermore, I have been studying a paper titled "Innovative Stabilization Diagram for Automated Structural Modal Identification based on ERA and Hierarchical Cluster Analysis" [1]. With the aid of the anticipated updated package, it would undoubtedly facilitate the provision of this diagram effortlessly 👍
Thank you again for your outstanding work, and please let me know if there's anything I can do to assist you in your endeavors.
Sincerely,
Ayubirad
[1] https://link.springer.com/article/10.1007/s13349-021-00514-8
Hi @Ayubirad,
As of version 1.2.7, I added the "phi" attribute following the code shown above. I also included the calc for "realness" shown above, but I don't think it's that useful as is. The normal MPC indicator is enough for this type of check I think. Still, I may update the calculation of this factor for learning purposes.
I plan to check out the reference you indicate as well.
Best regards,
Tim