Tutorial 2.4: Deriving the EOM of a 4DOF floating wind turbine

Tutorial 2.4: Deriving the EOM of a 4DOF floating wind turbine#

In this tutorial you will learn to derive the Equations of Motion of a floating offshore wind turbine. The tutorial is quite elaborate, but contains large portions of example code to fully understand so that it can be applied in other offshore application projects

We will schematize the floating wind turbine as follows:

A rigid body, assumed rectangular, for all floaters together
- \(L\) = 50 m, \(B\) = 50m, \(H\) = 10 m, average density \(\rho\) = 500 kg/m3
- With horizontal anchor stiffness of 10 MN/m, and no vertical/rotational anchor stiffness. The anchor is attatched at the centre of the floater.
- With 3 DOF: surge, heave, roll
A turbine point mass connected via a rigid beam (assumed of negligible mass) and a rotational spring
- \(L_{tower}\) = 200 m, \(m_{turbine}\) = 12.000 kg
- A rotational spring of stiffness 15 GNm/rad
- Located at 3/4 from the left of the substructure.
Loaded by
- Wind
- Waves

Assume seawater of \(\rho_{water}\) = 1025 kg/m3. The full system is without damping. We also assume small displacements. Additionally this means that the water forces are not dependent on the ship position in time.

Part 1: Kinematic equations#

We first start by defining the variables.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
from sympy import *

var("t x_s z_s phi_s phi_t")      
# x_s substructure surge, z_s heave, phi_s roll, phi_t tower angle (to vertical)
x_s = Function("x_s")(t)
z_s = Function("z_s")(t)
phi_s = Function("phi_s")(t)
phi_t = Function("phi_t")(t)

# Define the kinematic relations here

# contants needed
var("x_st z_st L_t")
# dx_st hotizontal distance substructure centre to tower base
# dz_st vertica distance, L_t tower length

# positions of nacelle
x_t = x_s + x_st*cos(phi_s) - z_st*sin(phi_s) - L_t*sin(phi_t)
z_t = z_s + x_st*sin(phi_s) + z_st*cos(phi_s) + L_t*cos(phi_t)

# for energy relations later
phi_st = phi_s - phi_t

The velocities can then be obtained using:

# Compute/define the velocities here
x_s_dot = diff(x_s, t)
z_s_dot = diff(z_s, t)
phi_s_dot = diff(phi_s, t)
x_t_dot = diff(x_t, t)
z_t_dot = diff(z_t, t)
phi_t_dot = diff(phi_t, t) # will not be sued for point masses
phi_st_dot = diff(phi_st, t) # velocity will not be used

Part 2: Energy equations#

Kinetic energy:#

There are 2 masses in this system. The substructure also has some rotational inertia.

var("rho_s B_s H_s W_s m_t")
# substructure density, breadth, height, width; turbine mass

# Define the kinetic energy here (T)
m_s = rho_s*B_s*H_s*W_s
J_s = 1/12*m_s*(B_s**2 + H_s**2)
T_s = 1/2*m_s*(x_s_dot**2 + z_s_dot**2) + 1/2*J_s*(phi_s_dot**2)
T_t = 1/2*m_t*(x_t_dot**2 + z_t_dot**2) + 1/2*0*(phi_t_dot**2) # point mass so J_t = 0
T = T_s+T_t

Potential energy:#

var("rho_w g kr_st k_h")
# water density, rotational spring stiffness
# Define the potential energy here (V)

k_s = B_s*W_s*rho_w

draft_s = (B_s*W_s*H_s*rho_s)/(B_s*W_s*rho_w)
KB_s = draft_s/2 # COB, due to constant shape
nabla_s = B_s*W_s*draft_s # Submerged volume, taken in neutral position
J_sub = 1/12*W_s*B_s**3
BM_s = J_sub/nabla_s
KG_s = H_s/2 # COG, due to uniform weight
GM_s = KB_s + BM_s - KG_s
kr_s = rho_w*g*nabla_s*GM_s # Nm/rad

V_s = m_s*g*z_s + 1/2*k_h*x_s**2 + 1/2*k_s*z_s**2 + 1/2*kr_s*phi_s**2
V_t = m_t*g*z_t + 1/2*kr_st*phi_st**2 # need relative angle for this spring
V = V_s + V_t

Work by external force#

The work done by an external force working in the horizontal directio. The forces act on block 1 and 2 respectively

F_wave = Function("F_wave")(t)
M_wave = Function("M_wave")(t)
F_wind = Function("F_wind")(t)

# Define your external work here (W)
W = F_wave*z_s + M_wave*phi_s + F_wind*x_t

Step 3: Construct the Lagrangian#

# Define your Lagrangian here (L)
L = T - V - W

Step 4: Obtaining the EoM#

In order to obtain the EoMs we have to take derivatives w.r.t. \(x_1\) and its velocity, as well as \(x_2\).

# Compute the EOM here
EOM_x_s = diff( diff(L, x_s_dot), t) - diff(L, x_s)
EOM_z_s = diff( diff(L, z_s_dot), t) - diff(L, z_s)
EOM_phi_s = diff( diff(L, phi_s_dot), t) - diff(L, phi_s)
EOM_phi_t = diff( diff(L, phi_t_dot), t) - diff(L, phi_t)

The equations are nonlinear since they depend on the cosines and sindes of the angles. We want to obtain a linear function, so let’s linearize. We consider a perturbation around the point we want to linearize.

On a technical note: Using SymPy for function substitution can be tricky, hence we will use temporary symbols called tmp in order to do the substitution.

# dictionaries for substitution
var("x_s_0 x_s_epsilon")
tmp1_x_s = symbols("tmp1_x_s")
psi_x_s = Function("psi_x_s")(t) # perturbation function
tmp2_x_s = symbols("tmp2_x_s")

var("z_s_0 z_s_epsilon")
psi_z_s = Function("psi_z_s")(t) # perturbation function
tmp1_z_s = symbols("tmp1_z_s")
tmp2_z_s = symbols("tmp2_z_s")

var("phi_s_0 phi_s_epsilon")
psi_phi_s = Function("psi_phi_s")(t) # perturbation function
tmp1_phi_s = symbols("tmp1_phi_s")
tmp2_phi_s = symbols("tmp2_phi_s")

var("phi_t_0 phi_t_epsilon")
psi_phi_t = Function("psi_phi_t")(t) # perturbation function
tmp1_phi_t = symbols("tmp1_phi_t")
tmp2_phi_t = symbols("tmp2_phi_t")

subs1_dict = {x_s: tmp1_x_s, z_s: tmp1_z_s, phi_s: tmp1_phi_s, phi_t: tmp1_phi_t}
subs2_dict = {tmp1_x_s: x_s_0 + x_s_epsilon*psi_x_s,
              tmp1_z_s: z_s_0 + z_s_epsilon*psi_z_s,
              tmp1_phi_s: phi_s_0 + phi_s_epsilon*psi_phi_s,
              tmp1_phi_t: phi_t_0 + phi_t_epsilon*psi_phi_t}
subs3_dict = {diff(x_s, (t, 2)): tmp2_x_s, x_s: tmp1_x_s,
              diff(z_s, (t, 2)): tmp2_z_s, z_s: tmp1_z_s,
              diff(phi_s, (t, 2)): tmp2_phi_s, phi_s: tmp1_phi_s,
              diff(phi_t, (t, 2)): tmp2_phi_t, phi_t: tmp1_phi_t}
subs4_dict = {tmp2_x_s: diff(x_s_0 + x_s_epsilon*psi_x_s, (t, 2)), tmp1_x_s: x_s_0 + x_s_epsilon*psi_x_s,
              tmp2_z_s: diff(z_s_0 + z_s_epsilon*psi_z_s, (t, 2)), tmp1_z_s: z_s_0 + z_s_epsilon*psi_z_s,
              tmp2_phi_s: diff(phi_s_0 + phi_s_epsilon*psi_phi_s, (t, 2)), tmp1_phi_s: phi_s_0 + phi_s_epsilon*psi_phi_s,
              tmp2_phi_t: diff(phi_t_0 + phi_t_epsilon*psi_phi_t, (t, 2)), tmp1_phi_t: phi_t_0 + phi_t_epsilon*psi_phi_t}
epsilons_dict = {x_s_epsilon: 1, z_s_epsilon: 1, phi_s_epsilon: 1, phi_t_epsilon: 1}


startpos_dict = {x_s_0: 0, z_s_0: 0, phi_s_0: 0, phi_t_0: 0}

# x_s
EOM_psi_x_s = EOM_x_s.evalf(subs=subs1_dict)
EOM_psi_x_s = EOM_psi_x_s.evalf(subs=subs2_dict)
EOM_lin_x_s = series(EOM_psi_x_s, x_s_epsilon, n=2)

EOM_psi2_x_s = EOM_x_s.evalf(subs=subs3_dict)
EOM_psi2_x_s = EOM_psi2_x_s.evalf(subs=subs4_dict)
EOM_lin_x_s = series(EOM_psi2_x_s, x_s_epsilon, n=2)

EOM_lin_x_s = EOM_lin_x_s.removeO().evalf(subs=epsilons_dict)
EOM_lin_x_s_simplified = EOM_lin_x_s.evalf(subs=startpos_dict) # makes symbolic calcs easier
EOM_lin_x_s_iso = solve(EOM_lin_x_s_simplified, diff(psi_x_s, (t, 2)))

x_s_dotdot = EOM_lin_x_s_iso[0].evalf(subs=startpos_dict)

# z_s
EOM_psi_z_s = EOM_z_s.evalf(subs=subs1_dict)
EOM_psi_z_s = EOM_psi_z_s.evalf(subs=subs2_dict)
EOM_lin_z_s = series(EOM_psi_z_s, z_s_epsilon, n=2)

EOM_psi2_z_s = EOM_z_s.evalf(subs=subs3_dict)
EOM_psi2_z_s = EOM_psi2_z_s.evalf(subs=subs4_dict)
EOM_lin_z_s = series(EOM_psi2_z_s, z_s_epsilon, n=2)

EOM_lin_z_s = EOM_lin_z_s.removeO().evalf(subs=epsilons_dict)
EOM_lin_z_s_simplified = EOM_lin_z_s.evalf(subs=startpos_dict) # makes symbolic calcs easier
EOM_lin_z_s_iso = solve(EOM_lin_z_s_simplified, diff(psi_z_s, (t, 2)))

z_s_dotdot = EOM_lin_z_s_iso[0].evalf(subs=startpos_dict)

# phi_s

EOM_psi_phi_s = EOM_phi_s.evalf(subs=subs1_dict)
EOM_psi_phi_s = EOM_psi_phi_s.evalf(subs=subs2_dict)
EOM_lin_phi_s = series(EOM_psi_phi_s, phi_s_epsilon, n=2)

EOM_psi2_phi_s = EOM_phi_s.evalf(subs=subs3_dict)
EOM_psi2_phi_s = EOM_psi2_phi_s.evalf(subs=subs4_dict)
EOM_lin_phi_s = series(EOM_psi2_phi_s, phi_s_epsilon, n=2)

EOM_lin_phi_s = EOM_lin_phi_s.removeO().evalf(subs=epsilons_dict)
EOM_lin_phi_s_simplified = EOM_lin_phi_s.evalf(subs=startpos_dict) # makes symbolic calcs easier
EOM_lin_phi_s_iso = solve(EOM_lin_phi_s_simplified, diff(psi_phi_s, (t, 2)))

phi_s_dotdot = EOM_lin_phi_s_iso[0].evalf(subs=startpos_dict)

# phi_t

EOM_psi_phi_t = EOM_phi_t.evalf(subs=subs1_dict)
EOM_psi_phi_t = EOM_psi_phi_t.evalf(subs=subs2_dict)
EOM_lin_phi_t = series(EOM_psi_phi_t, phi_t_epsilon, n=2)

EOM_psi2_phi_t = EOM_phi_t.evalf(subs=subs3_dict)
EOM_psi2_phi_t = EOM_psi2_phi_t.evalf(subs=subs4_dict)
EOM_lin_phi_t = series(EOM_psi2_phi_t, phi_t_epsilon, n=2)

EOM_lin_phi_t = EOM_lin_phi_t.removeO().evalf(subs=epsilons_dict)
EOM_lin_phi_t_simplified = EOM_lin_phi_t.evalf(subs=startpos_dict) # makes symbolic calcs easier
EOM_lin_phi_t_iso = solve(EOM_lin_phi_t_simplified, diff(psi_phi_t, (t, 2)))

phi_t_dotdot = EOM_lin_phi_t_iso[0].evalf(subs=startpos_dict)

Now we isolate it for the acceleration. Then we convert the EOM to matrices for easier interpretation.

var("acc1 acc2 acc3 acc4 vel1 vel2 vel3 vel4")

dict_values = { Derivative(psi_x_s, (t,2)): acc1,
                Derivative(psi_z_s, (t,2)): acc2,
                Derivative(psi_phi_s, (t,2)): acc3,
                Derivative(psi_phi_t, (t,2)): acc4,
                Derivative(psi_x_s, t): vel1,
                Derivative(psi_z_s, t): vel2,
                Derivative(psi_phi_s, t): vel3,
                Derivative(psi_phi_t, t): vel4}

EOM_1 = EOM_lin_x_s_simplified.evalf(subs=dict_values)
EOM_2 = EOM_lin_z_s_simplified.evalf(subs=dict_values)
EOM_3 = EOM_lin_phi_s_simplified.evalf(subs=dict_values)
EOM_4 = EOM_lin_phi_t_simplified.evalf(subs=dict_values)

MTRX = linear_eq_to_matrix([EOM_1, EOM_2, EOM_3, EOM_4],
                           [acc1, acc2, acc3, acc4])
# Note: The results per line are the same af from for example EOM_lin_phi_t_iso

M = MTRX[0]
F = MTRX[1]

\[\begin{split}\displaystyle \left[\begin{matrix}1.0 B_{s} H_{s} W_{s} \rho_{s} + 1.0 m_{t} & 0 & - 1.0 m_{t} x_{st} \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} - 1.0 m_{t} z_{st} \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} & - 1.0 L_{t} m_{t} \cos{\left(\psi_{\phi t}{\left(t \right)} \right)}\\0 & 1.0 B_{s} H_{s} W_{s} \rho_{s} + 1.0 m_{t} & 1.0 m_{t} x_{st} \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} - 1.0 m_{t} z_{st} \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} & - 1.0 L_{t} m_{t} \sin{\left(\psi_{\phi t}{\left(t \right)} \right)}\\- 1.0 m_{t} x_{st} \psi_{\phi s}{\left(t \right)} - 1.0 m_{t} z_{st} & 1.0 m_{t} x_{st} - 1.0 m_{t} z_{st} \psi_{\phi s}{\left(t \right)} & 0.0833333333333333 B_{s}^{3} H_{s} W_{s} \rho_{s} + 0.0833333333333333 B_{s} H_{s}^{3} W_{s} \rho_{s} + 1.0 m_{t} x_{st}^{2} + 1.0 m_{t} z_{st}^{2} & 1.0 L_{t} m_{t} x_{st} \psi_{\phi s}{\left(t \right)} \cos{\left(\psi_{\phi t}{\left(t \right)} \right)} - 1.0 L_{t} m_{t} x_{st} \sin{\left(\psi_{\phi t}{\left(t \right)} \right)} + 1.0 L_{t} m_{t} z_{st} \psi_{\phi s}{\left(t \right)} \sin{\left(\psi_{\phi t}{\left(t \right)} \right)} + 1.0 L_{t} m_{t} z_{st} \cos{\left(\psi_{\phi t}{\left(t \right)} \right)}\\- 1.0 L_{t} m_{t} & - 1.0 L_{t} m_{t} \psi_{\phi t}{\left(t \right)} & - 1.0 L_{t} m_{t} x_{st} \psi_{\phi t}{\left(t \right)} \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} + 1.0 L_{t} m_{t} x_{st} \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} + 1.0 L_{t} m_{t} z_{st} \psi_{\phi t}{\left(t \right)} \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} + 1.0 L_{t} m_{t} z_{st} \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} & 1.0 L_{t}^{2} m_{t}\end{matrix}\right]\end{split}\]

\[\begin{split}\displaystyle \left[\begin{matrix}- 1.0 L_{t} m_{t} vel_{4}^{2} \sin{\left(\psi_{\phi t}{\left(t \right)} \right)} - 1.0 k_{h} \psi_{x s}{\left(t \right)} + 1.0 m_{t} vel_{3}^{2} x_{st} \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} - 1.0 m_{t} vel_{3}^{2} z_{st} \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} - F_{wind}{\left(t \right)}\\- B_{s} H_{s} W_{s} g \rho_{s} - 1.0 B_{s} W_{s} \rho_{w} \psi_{z s}{\left(t \right)} + 1.0 L_{t} m_{t} vel_{4}^{2} \cos{\left(\psi_{\phi t}{\left(t \right)} \right)} - g m_{t} + 1.0 m_{t} vel_{3}^{2} x_{st} \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} + 1.0 m_{t} vel_{3}^{2} z_{st} \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} - F_{wave}{\left(t \right)}\\- 0.0833333333333333 B_{s}^{3} W_{s} g \rho_{w} \psi_{\phi s}{\left(t \right)} - \frac{0.5 B_{s} H_{s}^{2} W_{s} g \rho_{s}^{2} \psi_{\phi s}{\left(t \right)}}{\rho_{w}} + 0.5 B_{s} H_{s}^{2} W_{s} g \rho_{s} \psi_{\phi s}{\left(t \right)} + 1.0 L_{t} m_{t} vel_{4}^{2} x_{st} \psi_{\phi s}{\left(t \right)} \sin{\left(\psi_{\phi t}{\left(t \right)} \right)} + 1.0 L_{t} m_{t} vel_{4}^{2} x_{st} \cos{\left(\psi_{\phi t}{\left(t \right)} \right)} - 1.0 L_{t} m_{t} vel_{4}^{2} z_{st} \psi_{\phi s}{\left(t \right)} \cos{\left(\psi_{\phi t}{\left(t \right)} \right)} + 1.0 L_{t} m_{t} vel_{4}^{2} z_{st} \sin{\left(\psi_{\phi t}{\left(t \right)} \right)} - g m_{t} x_{st} + g m_{t} z_{st} \psi_{\phi s}{\left(t \right)} - 1.0 kr_{st} \psi_{\phi s}{\left(t \right)} + 1.0 kr_{st} \psi_{\phi t}{\left(t \right)} + x_{st} F_{wind}{\left(t \right)} \psi_{\phi s}{\left(t \right)} + z_{st} F_{wind}{\left(t \right)} - M_{wave}{\left(t \right)}\\L_{t} g m_{t} \psi_{\phi t}{\left(t \right)} - 1.0 L_{t} m_{t} vel_{3}^{2} x_{st} \psi_{\phi t}{\left(t \right)} \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} - 1.0 L_{t} m_{t} vel_{3}^{2} x_{st} \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} - 1.0 L_{t} m_{t} vel_{3}^{2} z_{st} \psi_{\phi t}{\left(t \right)} \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} + 1.0 L_{t} m_{t} vel_{3}^{2} z_{st} \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} + L_{t} F_{wind}{\left(t \right)} + 1.0 kr_{st} \psi_{\phi s}{\left(t \right)} - 1.0 kr_{st} \psi_{\phi t}{\left(t \right)}\end{matrix}\right]\end{split}\]

As one can see, this system is already linearized. We can therefore immediately start solving the system. In case of a non-linear example, the linearization must be performed for each EOM separately.

Step 4: Solve the equation#

For students looking for a challenge

The student is invited to try and plot a 100s trace of the horizontal and vertical nacelle movement. All relevant info can be find in the previous tutorial on how to use solve_ivp for multiple DOF, and on how to plot the relevant output. The floater is subjected to the following force (in MN):

\[ F_h(t) = 1.25 \cos(0.1 t) + 0.5 \cos(0.5 t - 0.75) + 0.1 \cos(1.3 t + 2.28) \]

Hint: Think about how the nacelle moves. Which DOF contribute to its movement?
Hint: Looking at the force in 1 direction (for the horizontal DOF), which movements do you excpect? Are the DOF coupled? How could one see that in the DOF-matrices?

In practice an irregular wave consits of waves generated from a spectrum, each with their respective amplitude, frequency and/or phase difference. More information is given here. Ask yourself then:

Which forces act on the structure?
Assuming small deformations, how does a water wave result in a force and moment in time on the structure Think about wave potentials, their pressures, and integration over the area on which the water acts

# Filling in values

L = 50
B = 50
H = 10
rho_ship = 400
L_tower = 200
m_turbine = 12e3
k_r = 15e9
x_r = L/4 # Distance from centre of floater
z_r = H/2
rho_water = 1025
g_val = 9.81
k_anchor = 10e6

def F_h(t):
    return (1.25 * cos(0.1*t) + 0.5*cos(0.5*t - 0.75) + 0.1*cos(1.3*t + 2.28))*10e3

dict_fillin = {B_s: L, H_s:H, W_s:B, m_t: m_turbine, z_st: z_r, x_st: x_r,
               L_t: L_tower, rho_s: rho_ship, g: g_val, rho_w: rho_water,
               F_wind: 0, kr_st: k_r, M_wave: 0, k_h: k_anchor}

M_filled = M.evalf(subs=dict_fillin)
M_filled

\[\begin{split}\displaystyle \left[\begin{matrix}10012000.0 & 0 & - 150000.0 \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} - 60000.0 \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} & - 2400000.0 \cos{\left(\psi_{\phi t}{\left(t \right)} \right)}\\0 & 10012000.0 & - 60000.0 \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} + 150000.0 \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} & - 2400000.0 \sin{\left(\psi_{\phi t}{\left(t \right)} \right)}\\- 150000.0 \psi_{\phi s}{\left(t \right)} - 60000.0 & 150000.0 - 60000.0 \psi_{\phi s}{\left(t \right)} & 2168841666.66667 & 12000000.0 \psi_{\phi s}{\left(t \right)} \sin{\left(\psi_{\phi t}{\left(t \right)} \right)} + 30000000.0 \psi_{\phi s}{\left(t \right)} \cos{\left(\psi_{\phi t}{\left(t \right)} \right)} - 30000000.0 \sin{\left(\psi_{\phi t}{\left(t \right)} \right)} + 12000000.0 \cos{\left(\psi_{\phi t}{\left(t \right)} \right)}\\-2400000.0 & - 2400000.0 \psi_{\phi t}{\left(t \right)} & 12000000.0 \psi_{\phi t}{\left(t \right)} \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} - 30000000.0 \psi_{\phi t}{\left(t \right)} \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} + 30000000.0 \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} + 12000000.0 \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} & 480000000.0\end{matrix}\right]\end{split}\]

F_filled = F.evalf(subs=dict_fillin)
F_filled

\[\begin{split}\displaystyle \left[\begin{matrix}- 60000.0 vel_{3}^{2} \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} + 150000.0 vel_{3}^{2} \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} - 2400000.0 vel_{4}^{2} \sin{\left(\psi_{\phi t}{\left(t \right)} \right)} - 10000000.0 \psi_{x s}{\left(t \right)}\\150000.0 vel_{3}^{2} \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} + 60000.0 vel_{3}^{2} \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} + 2400000.0 vel_{4}^{2} \cos{\left(\psi_{\phi t}{\left(t \right)} \right)} - F_{wave}{\left(t \right)} - 2562500.0 \psi_{z s}{\left(t \right)} - 98217720.0\\30000000.0 vel_{4}^{2} \psi_{\phi s}{\left(t \right)} \sin{\left(\psi_{\phi t}{\left(t \right)} \right)} - 12000000.0 vel_{4}^{2} \psi_{\phi s}{\left(t \right)} \cos{\left(\psi_{\phi t}{\left(t \right)} \right)} + 12000000.0 vel_{4}^{2} \sin{\left(\psi_{\phi t}{\left(t \right)} \right)} + 30000000.0 vel_{4}^{2} \cos{\left(\psi_{\phi t}{\left(t \right)} \right)} - 19937435409.1463 \psi_{\phi s}{\left(t \right)} + 15000000000.0 \psi_{\phi t}{\left(t \right)} - 1471500.0\\- 30000000.0 vel_{3}^{2} \psi_{\phi t}{\left(t \right)} \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} - 12000000.0 vel_{3}^{2} \psi_{\phi t}{\left(t \right)} \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} + 12000000.0 vel_{3}^{2} \sin{\left(\psi_{\phi s}{\left(t \right)} \right)} - 30000000.0 vel_{3}^{2} \cos{\left(\psi_{\phi s}{\left(t \right)} \right)} + 15000000000.0 \psi_{\phi s}{\left(t \right)} - 14976456000.0 \psi_{\phi t}{\left(t \right)}\end{matrix}\right]\end{split}\]

# Initial conditions and solve_ivp
q0 = np.zeros(8)

def qdot(t,q):

    print(t)
    
    vt = q[int(len(q)/2):]

    F_x_s = F_h(t)

    dict_values = { psi_x_s: q[0], psi_z_s: q[1],psi_phi_s: q[2], psi_phi_t: q[3], 
                    vel1: q[4], vel2: q[4], vel3: q[6], vel4: q[7],
                    F_wave: F_x_s}  
    Mass_matr =  M_filled.evalf(subs=dict_values)
    Force_vect = F_filled.evalf(subs=dict_values)
    at = Mass_matr.inv()*Force_vect

    return list(vt) + list(np.transpose(at)[0])

sol = solve_ivp(fun=qdot,t_span=[0,100],y0=q0) # Takes around 6 minutes to run 100 seconds of simulation

0
1e-06
9999999999999998e-05
9999999999999997e-05

999999999999999e-05
888888888888888e-05
999999999999999e-05

9.999999999999999e-05
0.00030000000000000003

00039999999999999996
0009
0009888888888888888

0.0011
0.0011

0031000000000000003
0041
0091

0.00998888888888889
0.0111

0111
031100000000000003
0411

09110000000000001
09998888888888889
1111

0.1111
0.271793766300492

352140649450738
753875065201968
8252945168910755

9145688315024599
9145688315024599
072664136876536

1517117895635738
5469500529987636
6172146331650197

7050453583728395
7050453583728395
02285817335736

07700284428481
3477261989220604
3958547953020162

4560155407769606
4560155407769606
5643048826318608

6184495535593109
8891729081965614
9373015045765167

1.9974622500514614
1.9974622500514614

1058234275283247
1600040162667566
430906959958915

4790674832819652
539268137435778
539268137435778

0784767992070643
118984073784866
3215204466738735

2.3575269129652523
2.4025349958294764

4025349958294764
4835495449850793
524056819562881

7265931924518885
7625996587432673
8076077416074914

8076077416074914
8824048990159237
91980347772014

3.1067963712412205
3.140039552311635

3.181593528649653
3.181593528649653

2500888302205486
2843364810059965
4555747349332355

486017091186967
5240700365041313
5240700365041313

592783821142421
627140713461566
7989251750572905

3.82946463489653
3.8676389596955802

8676389596955802
964310880731274
012646841249121

254326643838356
297291942076441
3509985648740495

4.3509985648740495
3.93272719801245

965271317170885
12799191296306
156920018881668

193080151279929
193080151279929
258168389596799

4.290712508755234
4.4534331045474085

4.482361210466017
4.518521342864278

518521342864278
584440748214873
61740045089017

782198964266656
811496477755809
848118369617251

848118369617251
944457820777153
992627546357103

5.233476174256857
5.276293708105703

329815625416759
329815625416759
91024619257417

4.941310104052629
5.0966296614449265

124242027203557
158757484401845
158757484401845

5.220885307358764
5.251949218837224

407268776229521
434881141988152
46939659918644

46939659918644
5341784649324754
566569397805493

5.728524062170582
5.7573160025021535

7933059279166175
7933059279166175
874202309218608

5.9146504998696035
6.1168914531245795

152845400369909
19778783442657
19778783442657

5.854376915879278
5.8849124098606085

03758987976726
064732541083998
0986608677299206

0986608677299206
1597318556925815
190267349673912

6.342944819580563
6.370087480897301

6.404015807543224
6.404015807543224

473634240775971
508443457392345
682489540474215

713431066355436
752107973706963
752107973706963

81899212603397
852434202197473
019644583014989

049370872938103
086528735341996
086528735341996

146784211243605
17691194919441
327550638948434

3543308504602605
387806114850044
387806114850044

4575784125584175
492464561412604
666895305683539

697905215776149
736667603391913
736667603391913

7.8006476011986665
7.832637600102044

992587594618928
021023149199706
056567592425681

056567592425681
11556143600056
145058357788

292542966725199
318762452758477
351536810300077

351536810300077
42546599975627
462430594484367

8.647253568124848
8.680110985660933

72118275758104
72118275758104
41193545740932

442134780963944
593131398737054
619975241896718

8.653530045846297
8.653530045846297

713928692955541
744128016510164
895124634283274

921968477442938
955523281392518
955523281392518

017080905716087
047859717877873
201753778686799

229112722830607
263311403010368
263311403010368

33701759883865
373870696752789
558136186323493

9.590894495580507
9.631842382151774

631842382151774
322306210037075
35180361355043

4992906311172
525510545351292
558285438143907

558285438143907
617280245170614
646777648683969

9.794264666250738
9.82048458048483

853259473277445
853259473277445
92606715796806

962471000313368
144490212039907
176849183013513

217297896730521
217297896730521
913392412915202

943458882734083
093791231828478
120516982778591

153924171466235
153924171466235
214057111103992

244123580922873
394455930017267
421181680967381

454588869655025
454588869655025
516810591400748

10.547921452273611
10.703475756637923

731129855191577
765697478383647
765697478383647

837187914414494
872933132429917
051659222507034

083432749631855
123149658537882
123149658537882

822777330801431
851317257010324
994016888054785

11.019385711351578
11.05109674047257

05109674047257
108176592890354
136716519099247

11.279416150143708
11.304784973440501

336496002561493
336496002561493
408304515120227

11.444208771399595
11.623730052796429

655644947266977
695538565355163
695538565355163

11.39460381820015
11.423657726019476

568927265116116
59475296095552
627035080754773

627035080754773
68514289639343
714196804212756

859466343309396
8852920391488
917574158948053

917574158948053
978324455095503
00869960316923

160575343537857
187575475158946
221325639685308

221325639685308
284481003675658
316058685670832

12.473947095646707
12.502016146309085

537102459637058
537102459637058
277147379143871

305058248873154
444612597519564
469422259501147

500434336978127
500434336978127
556256076436691

584166946165974
723721294812384
748530956793967

779543034270947
779543034270947
860515922674871

901002366876833
103434587886644
139422538288388

184407476290568
184407476290568
835433015553525

12.863378006194814
13.003102959401257

027942951082403
058992940683835
058992940683835

13.114882921966412
13.1428279126077

13.282552865814145
13.30739285749529

338442847096722
338442847096722
407005508069044

441286838555206
61269349098601
643165784751487

681256151958332
681256151958332
395276293985958

423693017430574
565776634653663
591035944382213

13.622610081542899

13.622610081542899
13.679443528432135

70786025187675
849943869099839
87520317882839

906777315989075
906777315989075
968278245905262

999028710863355
152781035653822
180114782283239

21428196557001
21428196557001
276370707446139

307415078384205
46263693307453
49023192946392

524725674950659
524725674950659
269673745436634

297369635369945
435849085036505
460467653866116

49124086490313
49124086490313
546632644769755

574328534703065
712807984369626
737426553199237

14.76819976423625
14.76819976423625

847588080798392
887282239079463
085753030484815

121036726734657
165141347046957
165141347046957

823329567931891
850894469779712
988718979018815

15.013221113994655
15.043848782714456

043848782714456
098978586410096
126543488257918

15.26436799749702
15.28887013247286

31949780119266
31949780119266
391011126152375

42676778863223
605551101031512
637334801013608

677064425991226
677064425991226
374586269566182

402130503752943
539851674686744
564335438408309

594940143060265
594940143060265
650028611433786

15.677572845620546
15.815294016554349

839777780275913
87038248492787
87038248492787

15.933794107847547
15.965499919307385

124028976606578
152211920126433
187440599526255

187440599526255
245706397140477
274839295947586

16.42050378998314
16.446399700033908

16.478769587597363
16.478769587597363

534656065808882
562599304914645
702315500443447

72715393520412
758201978654967
758201978654967

16.83331774106212
16.870875622265697

058665028283578
09204981157565
133780790690732

133780790690732
81330015512416
840849243358754

16.978594684531735
17.00308276296249

17.03369286100093
17.03369286100093

088791037470124
116340125704717
2540855668777

27857364530845
309183743346892
309183743346892

385133586692696
4231085083656
612983116730113

646738602661582
688932960075917
688932960075917

36312902596323
390101667271395
52496487381223

17.548940554975047
17.578910156428567

578910156428567
632855439044903
65982808035307

794691286893904
81866696805672
84863656951024

17.84863656951024
17.913478185275988

94589899315886
108003032573222
13682152846911

17284464833897
17284464833897
90456713647376

932532419955518
072358837364312
097216867125876

12828940432783
12828940432783
18421997129135

18.212185254773107
18.352011672181902

376869701943466
40794223914542
40794223914542

466156116463335
495263055122294
64079774841708

666670582780597
699011625734993
699011625734993

757652294057618
786972628218933
933574299025498

959636818279996
992214967348122
992214967348122

19.04626632114089
19.073291998037273

208420382519193
232443206427092
262471736311962

262471736311962
340743854400113
37987991344419

575560208664566
610347816703744
653832326752717

653832326752717
31595263834016
342693089354263

476395344424763
50016463421507
52987624645296

52987624645296
58335714848116
610097599495262

74379985456576
76756914435607
79728075659396

79728075659396
87438713568573
912940325231617

105706272961044
139975774779607
182812652052814

182812652052814
850372153688106
87691785223518

009646344970548
033242521456835
062737742064694

20.062737742064694
20.11582913915884

142374837705916
275103330441283
29869950692757

20.32819472753543
20.32819472753543

396120380338044
43008320673935
59989733874589

20.630086517769275
20.667822991548505

667822991548505
382716665025974
409977633771245

546282477497606
570514449715628
60080441498815

60080441498815
655326352478696
682587321223966

20.818892164950327
20.84312413716835

873414102440872
873414102440872
93486387252667

96558875756957
119213182784062
146524191711084

18066295286986
18066295286986
928568683267365

95614597368061
094032425746835
11854557278083

21.149187006573328
21.149187006573328

20434158739982
231918877813065
36980532987929

394318476913288
424959910705784
424959910705784

21.482528548733544

511312867747424
65523446281682
680820524162492

71280310084458
71280310084458
772791289920356

802785384458243
952755857147682
97941727451469

012744046223457
012744046223457
765603506310697

21.792003709043758
21.92400472270905

947471569582877
976805128175165
976805128175165

029525508460807
055885698603632
18768664931774

21111792944469
240407029603382
240407029603382

30891237297415
343165044659532
514428403086452

544875222362347
58293374645722
58293374645722

22.29339251341244
22.31988525531697

22.45234896483961
22.47589806875475

505334448648668
505334448648668
558319932457724

22.584812674362254
22.717276383884897

740825487800034
770261867693954
770261867693954

845231812929487
882716785547252
070141648636085

103461624296322
14511159387162
14511159387162

821726765456493
84745921433776
976121458744107

99899474663857
027586356506646
027586356506646

079016872321436
104732130228832
233308419765812

256166426794607
284738935580602
284738935580602

3591353037843
396333487886146
582324408395387

23.615389460930363
23.656720776599084

656720776599084
33711434102027
363302043740106

494240557339282
517518515312467
54661596277895

54661596277895
59899136821862
625179070938454

75611758453763
779395542510816
8084929899773

23.8084929899773

23.87828904855791
23.913187077848214

087677224299735
118697694780007
157473282880346

157473282880346
859826289433745
885492939161967

013826187803076
036640987561498
06515948725952

06515948725952
116492786715966
14215943644419

270492685085298
29330748484372
321825984541743

321825984541743
387475973850307
420300968504588

584425941775994
61360371480202
650075931084558

650075931084558
374484075839778
400813121488795

53245834973388
55586194586634
585116441031914

585116441031914
63777453232995
664103577978967

79574880622405
81915240235651
848406897522086

848406897522086
910614318042988
94171802830344

25.0972365796057
25.124884322059433

1594440001266
1594440001266
90038497034717

926374006759712
05631918882242
079420554522457

25.108297261647504
25.108297261647504

160275334472587
18626437088513
31620955294784

339310918647875
368187625772922
368187625772922

25.426981219809193
25.456378016827326

603362001918
629492488156345
66215559595427

25.66215559595427
25.715287596485524

741853596751152
87468359807929
898297820537625

25.927815598610543
25.927815598610543

983499566887343
01134155102574
150551471717733

175299902062978
206235439994533
206235439994533

26.264311208957345
26.29334909343875

438538515845778
46434996871814
49661428480859

26.49661428480859
26.54902536250917

575230901359458
7062585956109
72955240792227

75866967331148
75866967331148
827815470547332

86238836916526
035252862254893
065984327693048

104398659490744
104398659490744
8111302960124

83736060736286
968512164115158
99182799642668

020972786816078
020972786816078
073433409516998

27.099663720867458
27.230815277619755

254131109931276
283275900320675
283275900320675

358063105468233
395456708042012
58242472091091

61566347875427
65721192605847
65721192605847

334129594741203
359556441951465
48669067800278

50929231996746
53754437242331
53754437242331

588398066843837
6138249140541
740959150105414

763560792070095
791812844525943
791812844525943

86566620003107
902592877783636
087226266546455

28.12004998010429
28.161079622051584

161079622051584
84354396105112
869409519313713

27.99873731062666
28.02172891797118

050468427151838
050468427151838
102199543677017

28.128065101939605
28.257392893252554

280384500597076
309124009777733
309124009777733

38090566970449
416796499667868
59625064948476

628153609452205
668032309411515
668032309411515

359053016255977
3840175194951
508840035690707

531030705236592
55876904216895
55876904216895

608698048647195
633662551886317
758485068081924

28.78067573762781

28.80841407456017
28.80841407456017

28.876896060400924
28.9111370533213

082342017923192
11277845607464
150824003763947

150824003763947
85949513875369
88503567085045

012738331334255
035441026531377
063819395527776

063819395527776
114900459721298
14044099181806

268143652301863
290846347498984
319224716495384

319224716495384
387219726264146
42121723114853

591204755570438
62142475991211
6591997653392

6591997653392
369006473866367
39389735255186

518351745979317
54047697147753
5681335033503

29.5681335033503
29.617915260721283

642806139406773
767260532834232
789385758332447

817042290205215
817042290205215
883606690885035

916888891224943
083299892924487
11288407100441

149864293604306
149864293604306
867800242829194

893179219141185
020074100701134
042633190756234

30.070832053325113
30.070832053325113

12159000594909
146968982261082
27386386382103

29642295387613
32462181644501
32462181644501

391267680248472
424590612150205
591205271658858

620825655571508
65785113546232
65785113546232

30.37436711904096
30.39923977033894

30.523603026828823
30.545712050204802

573348329424775
573348329424775
623093632020726

647966283318706
772329539808588
794438563184567

30.82207484240454
30.82207484240454

30.88813097388317
30.92115903962249

086299368319075
115657648976246
152355499797707

31.152355499797707
30.87266190261236

897955432716273
024423083235824
046906221105967

075010143443645
075010143443645
125597203651466

150890733755375
27735838427493
299841522145073

32794544448275
32794544448275
394145625121958

42724571544156
592746167039575
62216846954589

658946347678782
658946347678782
377724878326354

31.40261459524816
31.527063179857173

549187372676553
576842613700776
576842613700776

62662204754438
651511764466186
7759603490752

31.79808454189458
31.825739782918802

825739782918802
89112364089405
92381556988168

087275214819805
116334707253245
15265907279505

15265907279505
87570529574259
900688052154486

02560183421396
04780872880231
07556734703775

07556734703775
125532859861536
15051561627343

32.2754293983329
32.297636292921254

32.32539491115669
32.32539491115669

39240109699268
425904189910675
59341965450064

32.623200181538856
32.66042584033663

66042584033663
37587650700439
40111730492825

527321294547505
5497575593687
577802890395205

577802890395205
628284486242904
65352528416676

77972927378602
80216553860721
83021086963372

83021086963372
89731512682575
930867255421774

33.09862789840186
33.12845201270943

16573215559389
16573215559389
88008161913398

90501699388411
02969386763477
05185864519044

079564617135034
079564617135034
12940845603358

15433037548285
278939972729205
30109279001745

32878381162775
32878381162775
396204981745065

42991556680372
59846849209701
6284334565936

665889662214326
665889662214326
38001403303256

405629143734956
53370469724697
556473684538

33.584934918651776
33.584934918651776

63616514005658
66178025075898
78985580427099

33.81262479156202
33.8410860256758

33.8410860256758
33.907698098385424

94100413474024
1075343165143
13713968216302

17414638922392
17414638922392
891120397033156

33.91613758271183
34.041223511105215

06346100948626
09125788246257
09125788246257

14129225381993
166309439498605
29139536789199

31363286627303
34142973924934
34142973924934

406411635559195
438902583714125
60135732448877

63023816729314
66633922079862
66633922079862

39270671358409
41834520075146
546537636588326

569327402959324
59781461092307
59781461092307

34.64909158525782
34.67473007242519

80292250826206
825712274633055
854199482596805

854199482596805
91766417882602
949396526940625

10805826751366
136264799171094
17152296374288

17152296374288
90338743754653
92798141502139

0509513023957
07281261570669
100139257345425

100139257345425
14932721229515
17392118977001

29689107714432
31875239045531
346079032094046

346079032094046
40826843451686
43936313572827

594836641785314
62247637619546
65702604420813

65702604420813
3965425118311
42177425169963

54793295104226
570361164258735
59839643077932

59839643077932
64885991051637
674091650384895

800250349727534
822678562944006
85071382946459

85071382946459
91293431107913
94404455188639

09959575592273
12724930330697
161816237537266

161816237537266
899985115425046
924620758405275

04779897330641
069697322622176
09707025926687

09707025926687
14634154522733
17097718820756

2941554031087
31605375242446
343426689069155

343426689069155
405583478665804
436661873464125

592053847455745
61967908727648
654210637052394

654210637052394
39295475981306
41771879518501

541538972044776
563551447930955
59106704278868

59106704278868
640595113532584
665359148904535

7891793257643
81119180165048
838707396508205

838707396508205
902183607674374
93392171325746

09261224117288
12082389058006
15608845233905

15608845233905
88793301362874
912545822189

37.035609864990334
37.05748791704391

08483548211087
08483548211087
1340610992314

37.158673907791666

281737950593
30361600264657
33096356771353

33096356771353
39578853226846
428201014545934

590263425933266
61907452129101
655088390488196

655088390488196
3799951119476
40451088406463

37.52708974464979
37.54888154208716

57612128888386
57612128888386
62515283311793

37.64966860523496
37.77224746582012

---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/core/evalf.py:1481, in evalf(x, prec, options)
   1480 try:
-> 1481     rf = evalf_table[type(x)]
   1482     r = rf(x, prec, options)

KeyError: psi_phi_t

During handling of the above exception, another exception occurred:

KeyboardInterrupt                         Traceback (most recent call last)
Cell In[20], line 21
     17     at = Mass_matr.inv()*Force_vect
     19     return list(vt) + list(np.transpose(at)[0])
---> 21 sol = solve_ivp(fun=qdot,t_span=[0,100],y0=q0) # Takes around 6 minutes to run 100 seconds of simulation

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/scipy/integrate/_ivp/ivp.py:591, in solve_ivp(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, args, **options)
    589 status = None
    590 while status is None:
--> 591     message = solver.step()
    593     if solver.status == 'finished':
    594         status = 0

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/scipy/integrate/_ivp/base.py:181, in OdeSolver.step(self)
    179 else:
    180     t = self.t
--> 181     success, message = self._step_impl()
    183     if not success:
    184         self.status = 'failed'

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/scipy/integrate/_ivp/rk.py:144, in RungeKutta._step_impl(self)
    141 h = t_new - t
    142 h_abs = np.abs(h)
--> 144 y_new, f_new = rk_step(self.fun, t, y, self.f, h, self.A,
    145                        self.B, self.C, self.K)
    146 scale = atol + np.maximum(np.abs(y), np.abs(y_new)) * rtol
    147 error_norm = self._estimate_error_norm(self.K, h, scale)

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/scipy/integrate/_ivp/rk.py:64, in rk_step(fun, t, y, f, h, A, B, C, K)
     62 for s, (a, c) in enumerate(zip(A[1:], C[1:]), start=1):
     63     dy = np.dot(K[:s].T, a[:s]) * h
---> 64     K[s] = fun(t + c * h, y + dy)
     66 y_new = y + h * np.dot(K[:-1].T, B)
     67 f_new = fun(t + h, y_new)

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/scipy/integrate/_ivp/base.py:138, in OdeSolver.__init__.<locals>.fun(t, y)
    136 def fun(t, y):
    137     self.nfev += 1
--> 138     return self.fun_single(t, y)

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/scipy/integrate/_ivp/base.py:20, in check_arguments.<locals>.fun_wrapped(t, y)
     19 def fun_wrapped(t, y):
---> 20     return np.asarray(fun(t, y), dtype=dtype)

Cell In[20], line 16, in qdot(t, q)
     12 dict_values = { psi_x_s: q[0], psi_z_s: q[1],psi_phi_s: q[2], psi_phi_t: q[3], 
     13                 vel1: q[4], vel2: q[4], vel3: q[6], vel4: q[7],
     14                 F_wave: F_x_s}  
     15 Mass_matr =  M_filled.evalf(subs=dict_values)
---> 16 Force_vect = F_filled.evalf(subs=dict_values)
     17 at = Mass_matr.inv()*Force_vect
     19 return list(vt) + list(np.transpose(at)[0])

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/matrices/common.py:2031, in MatrixOperations.evalf(self, n, subs, maxn, chop, strict, quad, verbose)
   2028 """Apply evalf() to each element of self."""
   2029 options = {'subs':subs, 'maxn':maxn, 'chop':chop, 'strict':strict,
   2030         'quad':quad, 'verbose':verbose}
-> 2031 return self.applyfunc(lambda i: i.evalf(n, **options))

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/matrices/common.py:1989, in MatrixOperations.applyfunc(self, f)
   1986 if not callable(f):
   1987     raise TypeError("`f` must be callable.")
-> 1989 return self._eval_applyfunc(f)

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/matrices/common.py:1931, in MatrixOperations._eval_applyfunc(self, f)
   1930 def _eval_applyfunc(self, f):
-> 1931     out = self._new(self.rows, self.cols, [f(x) for x in self])
   1932     return out

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/matrices/common.py:1931, in <listcomp>(.0)
   1930 def _eval_applyfunc(self, f):
-> 1931     out = self._new(self.rows, self.cols, [f(x) for x in self])
   1932     return out

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/matrices/common.py:2031, in MatrixOperations.evalf.<locals>.<lambda>(i)
   2028 """Apply evalf() to each element of self."""
   2029 options = {'subs':subs, 'maxn':maxn, 'chop':chop, 'strict':strict,
   2030         'quad':quad, 'verbose':verbose}
-> 2031 return self.applyfunc(lambda i: i.evalf(n, **options))

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/core/evalf.py:1647, in EvalfMixin.evalf(self, n, subs, maxn, chop, strict, quad, verbose)
   1645     options['quad'] = quad
   1646 try:
-> 1647     result = evalf(self, prec + 4, options)
   1648 except NotImplementedError:
   1649     # Fall back to the ordinary evalf
   1650     if hasattr(self, 'subs') and subs is not None:  # issue 20291

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/core/evalf.py:1482, in evalf(x, prec, options)
   1480 try:
   1481     rf = evalf_table[type(x)]
-> 1482     r = rf(x, prec, options)
   1483 except KeyError:
   1484     # Fall back to ordinary evalf if possible
   1485     if 'subs' in options:

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/core/evalf.py:600, in evalf_add(v, prec, options)
    597 while 1:
    598     options['maxprec'] = min(oldmaxprec, 2*prec)
--> 600     terms = [evalf(arg, prec + 10, options) for arg in v.args]
    601     n = terms.count(S.ComplexInfinity)
    602     if n >= 2:

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/core/evalf.py:600, in <listcomp>(.0)
    597 while 1:
    598     options['maxprec'] = min(oldmaxprec, 2*prec)
--> 600     terms = [evalf(arg, prec + 10, options) for arg in v.args]
    601     n = terms.count(S.ComplexInfinity)
    602     if n >= 2:

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/core/evalf.py:1482, in evalf(x, prec, options)
   1480 try:
   1481     rf = evalf_table[type(x)]
-> 1482     r = rf(x, prec, options)
   1483 except KeyError:
   1484     # Fall back to ordinary evalf if possible
   1485     if 'subs' in options:

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/core/evalf.py:648, in evalf_mul(v, prec, options)
    646 from .numbers import Float
    647 for arg in args:
--> 648     result = evalf(arg, prec, options)
    649     if result is S.ComplexInfinity:
    650         special.append(result)

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/core/evalf.py:1486, in evalf(x, prec, options)
   1483 except KeyError:
   1484     # Fall back to ordinary evalf if possible
   1485     if 'subs' in options:
-> 1486         x = x.subs(evalf_subs(prec, options['subs']))
   1487     xe = x._eval_evalf(prec)
   1488     if xe is None:

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/core/evalf.py:1023, in evalf_subs(prec, subs)
   1021 newsubs = {}
   1022 for a, b in subs.items():
-> 1023     b = S(b)
   1024     if b.is_Float:
   1025         b = b._eval_evalf(prec)

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/core/sympify.py:383, in sympify(a, locals, convert_xor, strict, rational, evaluate)
    381         conv = _sympy_converter.get(superclass)
    382     if conv is not None:
--> 383         return conv(a)
    385 if cls is type(None):
    386     if strict:

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/core/numbers.py:1069, in Float.__new__(cls, num, dps, precision)
   1067 elif num is S.NaN:
   1068     return num
-> 1069 elif _is_numpy_instance(num):  # support for numpy datatypes
   1070     num = _convert_numpy_types(num)
   1071 elif isinstance(num, mpmath.mpf):

File ~/.virtualenvs/.venv/lib/python3.9/site-packages/sympy/core/sympify.py:73, in _is_numpy_instance(a)
     68 """
     69 Checks if an object is an instance of a type from the numpy module.
     70 """
     71 # This check avoids unnecessarily importing NumPy.  We check the whole
     72 # __mro__ in case any base type is a numpy type.
---> 73 return any(type_.__module__ == 'numpy'
     74            for type_ in type(a).__mro__)

KeyboardInterrupt: 

# Plot
# From earlier: x_t = x_s + x_st*cos(phi_s) - z_st*sin(phi_s) - L_t*sin(phi_t)
# To see changes: do -x_r
x_s = sol.y[0]
phi_s = sol.y[2]
phi_t = sol.y[3]
move_nacelle = x_s + x_r*np.cos(phi_s) - z_r*np.sin(phi_s) - L_tower*np.sin(phi_t) - x_r

plt.plot(sol.t,move_nacelle, label="Nacelle movement at 3 frequencies")
plt.xlabel("Time [s]")
plt.ylabel("Excursion from initial position [m]")
plt.title("Nacelle from Lagrangian equations")
plt.legend();