# PIT on simplified linear Nomoto for many runs

# Purpose
Run PIT on a series of tests the find the best average parameters

# Methodology
* Load time series an concatenate into a large y and X.
* Run OLS regression.

# Setup


```python
# %load imports.py
%matplotlib inline
%load_ext autoreload
%autoreload 2

## External packages:
import pandas as pd
pd.options.display.max_rows = 999
pd.options.display.max_columns = 999
pd.set_option("display.max_columns", None)

import numpy as np
import os
import matplotlib.pyplot as plt

import plotly.express as px 
import plotly.graph_objects as go

import seaborn as sns
import sympy as sp
from sympy.physics.mechanics import (dynamicsymbols, ReferenceFrame,
                                      Particle, Point)
from sympy.physics.vector.printing import vpprint, vlatex
from IPython.display import display, Math, Latex
from src.substitute_dynamic_symbols import run, lambdify

import pyro

import sklearn
import pykalman
from statsmodels.sandbox.regression.predstd import wls_prediction_std
import statsmodels.api as sm

from scipy.integrate import solve_ivp

## Local packages:
from src.data import mdl



```

## Load tests


```python
df_runs = mdl.runs()
```


```python
mask = ~df_runs['sailing']
df_tests = df_runs.loc[mask]
```


```python
df = pd.DataFrame()

for id, row in df_tests.iterrows():

    df_ = mdl.load_run(id=id, dir_path='../data/processed/kalman/')
    df_['time'] = df_.index
    df_['id'] = id
    df = df.append(df_, ignore_index=True, )
    
```


```python
px.line(df, x='y0', y='x0', width=1000, height=1000, line_group='id', color='id')
```

# First order Nomoto model for ship maneuvering dynamics
The first order Nomoto model can be written as:


```python
K, delta, T_1, T_2 = sp.symbols('K delta T_1 T_2')
r = dynamicsymbols('r')

eq_nomoto = sp.Eq(-K*delta,
                 r + T_1*r1d + T_2*r1d.diff())
Math(vlatex(eq_nomoto))
```

where $r$ is yaw rate with its time derivatives and $\delta$ is the rudder angle. $K$, $T_{1}$
 and $T_{1}$ are the coefficients describing the hydrodynamics of the ship.
 
For slow manoeuvres this equation can be further simplified by removing the $\ddot{r}$ term:


```python
eq_nomoto_simple = eq_nomoto.subs(r1d.diff(),0)
Math(vlatex(eq_nomoto_simple))
```


```python
eq_r1 = sp.Eq(r1d,
             sp.solve(eq_nomoto_simple, r1d)[0])
r1d_lambda = lambdify(sp.solve(eq_nomoto_simple, r1d)[0])
Math(vlatex(eq_r1))

```


```python
def nomoto(t,states,df_control,params):
    """
    Simulation model for heading and yaw rate using Nomoto with only K and T_1
    """
    # states:
    psi = states[0]
    r = states[1]
    
    index = df_control.index.get_loc(t, method='nearest')
    delta = float(df_control['delta'].iloc[index])
    r1d_ = r1d_lambda(K=params['K'], T_1=params['T_1'], delta=delta, r=r)
    d_states_dt = np.array([r,r1d_])
    return d_states_dt
```

## Regress nomoto parameters


```python
eq_r1
```


```python
eq_r1.rhs
```


```python
y = df['r1d']
X = df[['delta','r']] 

#y = df_control['r1d_numerical']
#X = df_control[['delta','r_numerical']] 
#X.rename(columns={'r_numerical': 'r'}, inplace=True)

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


```python
df_pred = df.copy()
df_pred['r1d_pred'] = results.predict(X)
prstd, iv_l, iv_u = wls_prediction_std(results, exog=X, alpha=0.10)
df_pred['iv_l'] = iv_l
df_pred['iv_u'] = iv_u


px.line(data_frame=df_pred.resample('1S').mean(), y=['r1d_pred','iv_l','iv_u'], width=1400, height=300,)

```


```python
params={}
params['T_1'] = -1/results.params['r']
params['K'] = -results.params['delta']*params['T_1']
```


```python
def simulate(df, params):
    
    states_0 = np.array([df.iloc[0]['psi'],df.iloc[0]['r']])

    df_sim = df.copy()

    t = df_sim.index

    sol = solve_ivp(fun = nomoto,t_span=[t[0],t[-1]],t_eval=t, y0 = states_0, args=(df_sim,params))
    df_sim['psi'] = sol.y[0,:]
    df_sim['r'] = sol.y[1,:]
    
    return df_sim
    
```


```python
def plot_sim(df, df_sim):
    
    fig,ax=plt.subplots()
    df['psi_deg'] = np.rad2deg(df['psi'])
    df.plot(y='psi_deg', label='model test', ax=ax)
    
    df['-delta_deg'] = -np.rad2deg(df['delta'])
    df.plot(y='-delta_deg', label=r'$-\delta$', ax=ax)
    
    df_sim['psi_deg'] = np.rad2deg(df_sim['psi'])
    df_sim.plot(y='psi_deg', label='simulation', ax=ax)
    
    ax.grid(True)
    ax.set_title('Resimulation with Nomoto model')
    ax.set_ylabel('$\psi$ [deg]')
    ax.set_xlabel('time [s]');
    ax.legend();
    return ax

```


```python
y_lim_min = (-10,10)

for id, test in df.groupby(by='id'):
    
    meta_data = df_tests.loc[id]
    
    test.index = pd.TimedeltaIndex(test['time']).total_seconds()
    
    df_sim = simulate(df=test, params=params)
    ax = plot_sim(df=test, df_sim=df_sim)
    ax.set_title(f'test : {id} {meta_data["test_type"]} {meta_data["comment"]}')
    
    y_lim = ax.get_ylim()
    ax.set_ylim(np.min([y_lim_min[0], y_lim[0]]),
                np.max([y_lim_min[1], y_lim[1]]))
        
```


```python

```

3d7c2a02d0728f02cbce3d7d4a3fbf01d1417081
ipynb
Jupyter Notebook
notebooks/04.02_PIT_nomotos.ipynb
martinlarsalbert/wPCC
16e0d4cc850d503247916c9f5bd9f0ddb07f8930
MIT
notebooks/04.02_PIT_nomotos.ipynb
martinlarsalbert/wPCC
16e0d4cc850d503247916c9f5bd9f0ddb07f8930
MIT
notebooks/04.02_PIT_nomotos.ipynb
martinlarsalbert/wPCC
16e0d4cc850d503247916c9f5bd9f0ddb07f8930
MIT
Qwen/Qwen-72B
1. YES
2. YES