This notebook loads data from a ship roll decay test simulation and tries to predict the parameters in a matematical model describing the dynamics of this system. If you are not familiar with the ship roll decay test, you can check out this 
[Video](https://youtu.be/gN123b8gM-o).


```python
# This Python 3 environment comes with many helpful analytics libraries installed
# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python
# For example, here's several helpful packages to load

import numpy as np # linear algebra
import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)
import matplotlib.pyplot as plt
import sympy as sp
import sympy.physics.mechanics as me
from sympy.physics.vector.printing import vlatex
from IPython.display import Math
import statsmodels.api as sm
from sklearn.metrics import r2_score
from scipy.integrate import solve_ivp

# Input data files are available in the read-only "../input/" directory
# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory

import os
for dirname, _, filenames in os.walk('/kaggle/input'):
    for filename in filenames:
        print(os.path.join(dirname, filename))

# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using "Save & Run All" 
# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session
```

# Loading data from the simulation


```python
data = pd.read_csv('/kaggle/input/kvlcc2-roll-decay/fnpf_kvlcc2_rolldecay_0kn.csv', index_col=0)
```

This file contains time series for roll angle: $\phi$ (**phi**) and its first and second time derivatives $\dot{\phi}$ (**phi1d**) and $\ddot{\phi}$ (**phi2d**).


```python
data.head()
```

The ship is forces into an initial roll angle at the start of the test. When the ship is released it starts an oscillation motion that is decaying over time, due to the hydrodynamic roll damping of the ship.


```python
fig,ax=plt.subplots()
fig.set_size_inches(20,5)
data.plot(y='phi', ax=ax);
ax.set_ylabel(r'Roll angle $\phi$ $[rad]$')
ax.grid(True)
```

The ship is now behaving very much like the classical spring-mass-damper system in dynamics:



This system can be described by the following differential equation:


```python
phi = me.dynamicsymbols('phi')  # Roll angle
phi_dot = phi.diff()
phi_dot_dot = phi_dot.diff()

A_44 = sp.symbols('A_44')  # description='Total mass moment of inertia', unit='kg*m**2'
B_44 = sp.Function('B_{44}')(phi_dot)
C_44 = sp.Function('C_{44}')(phi)

lhs = A_44*phi_dot_dot + B_44 + C_44
rhs = 0
roll_equation_himeno = sp.Eq(lhs=lhs, rhs=rhs)
Math(vlatex(roll_equation_himeno))
```

where $A_{44}$ is the total roll intertia (representing the mass $m$), $B_{44}(\dot{\phi})$ is a function for the roll damping ($c$ in the figure above) and $C_{44}(\phi)$ is the stiffness (represented by the spring force $k$ in the figure). The $A_{44}=1$ is assumed to give this equation one unique solution. The equation is now ritten as:


```python
roll_equation_normalized = roll_equation_himeno.subs(A_44,1)
Math(vlatex(roll_equation_normalized))
```

## Stiffness model
The following polynomial will be used for the stiffness model:


```python
C_1, C_3, C_5 = sp.symbols('C_1 C_3 C_5')
restoring_equation_cubic = sp.Eq(C_44, C_1 * phi + C_3 * phi ** 3 + C_5 * phi ** 5)
restoring_equation_cubic
```

## Damping model
The following polynomial will be used for the damping model:


```python
B_1, B_2, B_3 = sp.symbols('B_1 B_2 B_3')
b44_cubic_equation = sp.Eq(B_44, B_1 * phi_dot + B_2 * phi_dot * sp.Abs(phi_dot) + B_3 * phi_dot ** 3)
Math(vlatex(b44_cubic_equation))
```

The model can now be rewritten with these polynomials:


```python
roll_decay_equation_cubic = roll_equation_normalized.subs([
    (C_44, restoring_equation_cubic.rhs),
    (B_44, b44_cubic_equation.rhs),
    
])
Math(vlatex(roll_decay_equation_cubic))
```


```python
def calculate_features(data):
    phi=data['phi']
    phi1d=data['phi1d']
    data_features = -np.array([phi1d, np.abs(phi1d)*phi1d, phi1d**3, phi, phi**3, phi**5]).T
    return data_features

def generate_features_and_label(data):
    
    data_features = calculate_features(data=data)
    
    columns = ['B_1','B_2','B_3','C_1','C_3','C_5']
    X = pd.DataFrame(data=data_features, index=data.index, columns=columns)
    
    y = data['phi2d']
    
    return X,y
  
X,y = generate_features_and_label(data=data)
```

# Regression
Regression to estimate the parameters $B_1$, $B_2$,..., from the time series.

## Linear regression


```python
model = sm.OLS(y,X)
results = model.fit()
```


```python
results.summary()
```


```python
y_pred = results.predict(X)
```

## Simpler model


```python
def generate_features_and_label_simple(data):
    
    phi=data['phi']
    phi1d=data['phi1d']
    phi2d=data['phi2d']
    
    X = -np.array([phi1d, phi, phi**3, phi**5]).T
    y = phi2d.values
    
    return X,y
  
X2,y2 = generate_features_and_label_simple(data=data)
```


```python
model2 = sm.OLS(y2,X2)
results2 = model2.fit()
```


```python
results2.summary()
```


```python
y_pred2 = results2.predict(X2)
```


```python
fig,ax=plt.subplots()
fig.set_size_inches(20,5)
data.plot(y='phi2d', ax=ax, label=r'$\ddot{\phi}$');

ax.plot(data.index, y_pred, '--', label=r'$\ddot{\phi}_{pred}$')
ax.plot(data.index, y_pred2, ':', label=r'$\ddot{\phi}_{pred2}$')

ax.set_ylabel(r'Roll angle acceleration $\ddot{\phi}$ $[rad/s^2]$')

ax.grid(True)
ax.legend();
```

## Simulate


```python
def updater(t, states):
    phi = states[0]
    phi1d = states[1]
        
    s = {
        'phi':phi,
        'phi1d':phi1d
    }
    x = calculate_features(data=s)
    phi2d = results.predict(x)[0]
        
    d_states_dt = np.array([phi1d, phi2d])
    return d_states_dt

y0 = [data.iloc[0]['phi'],
      data.iloc[0]['phi1d'],
     ]
t = data.index
t_span = [t[0], t[-1]]
result = solve_ivp(fun=updater, t_span=t_span, y0=y0, t_eval=t)
df_result = pd.DataFrame(index=result.t, data=result.y.T, columns = ['phi','phi1d'])
```


```python
fig,ax=plt.subplots()
fig.set_size_inches(20,5)
data.plot(y='phi', ax=ax, label=r'$\phi$');

df_result.plot(y='phi', style='--', label=r'$\phi_{sim}$', ax=ax)

ax.set_ylabel(r'Roll angle $\phi$ $[rad]$')

ax.grid(True)
ax.legend();


```


```python
r2_score(y_true=data['phi'], y_pred=df_result['phi'])
```


```python
error = data['phi']-df_result['phi']
error.hist(bins=100);
```

51f7b0a0a7e57edd48126ba6ac581bd839612580
ipynb
Jupyter Notebook
kernels/starter-notebook-kvlcc2-roll-decay/starter-notebook-kvlcc2-roll-decay.ipynb
martinlarsalbert/kaggle
5f75b0b7bf6adf1f5c9c20c2c3d4e1f6670716ac
MIT
kernels/starter-notebook-kvlcc2-roll-decay/starter-notebook-kvlcc2-roll-decay.ipynb
martinlarsalbert/kaggle
5f75b0b7bf6adf1f5c9c20c2c3d4e1f6670716ac
MIT
kernels/starter-notebook-kvlcc2-roll-decay/starter-notebook-kvlcc2-roll-decay.ipynb
martinlarsalbert/kaggle
5f75b0b7bf6adf1f5c9c20c2c3d4e1f6670716ac
MIT
Qwen/Qwen-72B
1. YES
2. YES