SISO for a Ground Motion Event

import mdof
from mdof import modal, transform
import sdof
import numpy as np
from numpy import linspace, sqrt, pi
from mdof.utilities.config import Config
from mdof.utilities.printing import print_modes, plot_transfer, plot_io, plot_pred
import matplotlib.pyplot as plt
import scienceplots
plt.style.use(["science"])

Unknown system with one input and one output

parameter	value
\(\omega_n\)	natural period
\(\zeta\)	damping ratio

Configure

Inputs

inputs  = np.loadtxt("uploads/elcentro.txt")
dt = 0.02

# ta = np.arange(0, inputs.size*dt, dt)
t = linspace(0, (inputs.size-1)*dt, len(inputs))

Outputs

outputs = None

# Example SDOF system
mass = 1       # mass
k = 30         # stiffness
zeta = 0.01    # damping ratio

Generate output if one was not given:

if outputs is None:
    omega_n = sqrt(k/mass)  # natural frequency (rad/s)
    Tn = 2*pi/omega_n  # natural periods (s)
    c = 2*zeta*mass*omega_n    # damping coefficient
    print(f"natural period: {Tn:<3.5}s")
    print(f"damping ratio:  {zeta}")
    displ, veloc, outputs = sdof.integrate(inputs,dt,k,c,mass)

natural period: 1.1471s
damping ratio:  0.01

Analysis with System Identification

Transfer Function Methods

# Set parameters
conf = Config()
conf.damping = zeta
conf.period_band = (0.1,3) # Period band (s)
conf.pseudo = True # use pseudo accelerations (Sa)

# A place to store models and their predictions
transfer_models = {}

# Generate a transfer function representation of the system
transfer_models["Fourier Transform"] = transform.fourier_transfer(inputs, outputs, dt, **conf)
transfer_models["Response Spectrum"] = transform.response_transfer(inputs, outputs, dt, **conf)

# Determing the fundamental frequency
fourier_periods, fourier_amplitudes = modal.spectrum_modes(*transfer_models["Fourier Transform"])
response_periods, response_amplitudes = modal.spectrum_modes(*transfer_models["Response Spectrum"], height=10)

State Space Methods

# Generate a state space realization of the system
conf.order = 2
conf.horizon = 300
realization = mdof.system(method="srim", inputs=inputs, outputs=outputs, **conf)
# Obtain natural period and damping ratio from the state space model
ss_modes = modal.system_modes(realization, dt, **conf)
print_modes(ss_modes, Tn=Tn, zeta=zeta)

Spectral quantities:
       T(s)        ζ        EMACO      MPC       EMACO*MPC     T % error    ζ % error
      1.148      0.00998    0.0        1.0        0.0          0.0999       -0.1995
Mean Period(s): 1.1482934446028537
Standard Dev(s): 0.0

Visualize Transfer Functions vs. State Space Methodds

plot_transfer(transfer_models, title="Transfer Functions")
color_cycle = plt.rcParams['axes.prop_cycle'].by_key()['color']
plt.gcf().axes[0].vlines([fourier_periods[0], response_periods[0]],ymin=0,ymax=1,linestyles='--',colors=color_cycle[:2])
plt.gcf().axes[0].vlines(1/ss_modes[next(iter(ss_modes.keys()))]["freq"],ymin=0,ymax=1,linestyles='--',colors='r',label="State Space")
plt.gcf().axes[0].legend()
plt.gcf().axes[0].text(fourier_periods[0]-1.05,0.95,r"$T_{fourier}$ = "+str(np.round(fourier_periods[0],3)),fontsize=15)
plt.gcf().axes[0].text(response_periods[0]-1.1,0.8,r"$T_{response}$ = "+str(np.round(response_periods[0],3)),fontsize=15)
plt.gcf().axes[0].set_ylim((0,1.1));

assert np.isclose(1/Tn, ss_modes[next(iter(ss_modes.keys()))]["freq"], atol=1e-2), (1/Tn, 1/ss_modes[next(iter(ss_modes.keys()))]["freq"])
assert np.isclose(Tn, fourier_periods[0], atol=1e-2)

plot_io(inputs=inputs, outputs=outputs, t=t, title="SDOF Response")

# Reproduce the response with the state space model
from control import ss, forced_response
out_mdof = forced_response(ss(*realization,dt), U=inputs, squeeze=False, return_x=False).outputs
plot_pred(ytrue=outputs, models=out_mdof, t=t, title="State Space Model Displacement Response")