# Regression revisited

# Purpose
Does the equation need to be decoupled?

# Setup


```python
# %load imports.py
## Local packages:

%matplotlib inline
%load_ext autoreload
%autoreload 2
%config Completer.use_jedi = False  ## (To fix autocomplete)

## 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
#if os.name == 'nt':
#    plt.style.use('presentation.mplstyle')  # Windows

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

from src.symbols import *
from src.parameters import *
import src.symbols as symbols
from src import prime_system
from src.models import regression
from src.visualization.plot import track_plot
from src.equation import Equation

## Load models:
# (Uncomment these for faster loading):
#import src.models.vmm_simple  as vmm
import src.models.vmm_martin  as vmm 
#import src.models.vmm_abkowitz as vmm 

```


```python
u1d,v1d,r1d = sp.symbols('u1d, v1d, r1d')
u_,v_,r_ = sp.symbols('u, v, r')

X_qs_, Y_qs_, N_qs_ = sp.symbols('X_qs, Y_qs, N_qs')


X_eq = vmm.X_eom.subs([
    (X_force,sp.solve(vmm.fx_eq,X_force)[0]),
])

Y_eq = vmm.Y_eom.subs([
    (Y_force,sp.solve(vmm.fy_eq,Y_force)[0]),
])


N_eq = vmm.N_eom.subs([
    (N_force,sp.solve(vmm.mz_eq,N_force)[0]),
])


subs = [
    #(X_qs,X_qs_), 
    #(Y_qs,Y_qs_), 
    #(N_qs,N_qs_),
   

    (u1d,u1d),
    (v1d,v1d),
    (r1d,r1d),

    (u,u_),
    (v,v_),
    (r,r_),
    
]
eq_X_ = X_eq.subs(subs)
eq_Y_ = Y_eq.subs(subs)
eq_N_ = N_eq.subs(subs)

A,b = sp.linear_eq_to_matrix([eq_X_,eq_Y_,eq_N_],[u1d,v1d,r1d])
```


```python
A
```


```python
b
```

Simulation model can be expressed as:


```python
acceleration = sp.matrices.MutableDenseMatrix([u1d,v1d,r1d])
sp.Eq(sp.UnevaluatedExpr(A)*sp.UnevaluatedExpr(acceleration),sp.UnevaluatedExpr(b))
```

The acceleartion can be calculated by calculating the inverse of the mass matrix:


```python
sp.Eq(sp.UnevaluatedExpr(acceleration),sp.UnevaluatedExpr(sp.simplify(A.inv()))*sp.UnevaluatedExpr(b))
```

This expression can be formulated in a more clean way by introducing the following term:


```python
A_inv = A.inv()
S = sp.symbols('S')
eq_S=sp.Eq(S,-sp.fraction(A_inv[1,1])[1])
eq_S
```


```python
A_inv_S = A_inv.subs(eq_S.rhs,S)
sp.Eq(sp.UnevaluatedExpr(acceleration),sp.UnevaluatedExpr(A_inv_S)*sp.UnevaluatedExpr(b))
```

Time integration of this expression will simulate the system.

This expression can also be formulated as a regression problem to identify the hydrodynamic coefficients within $X_{qs}$, $Y_{qs}$, $N_{qs}$ which are the quasi static force model. This is usually the part that differs between various simulation models, such as the Abkowitz model, the SeaMan model etc.

The regression problem is expressed as:


```python
coeff_matrix = sp.matrices.MutableDenseMatrix([A_coeff*X_X,B_coeff*X_Y,C_coeff*X_N])
eq_regression = sp.Eq(acceleration, coeff_matrix)
eq_regression
```

where $A_{coeff}$, $B_{coeff}$ and $C_{coeff}$ are coefficient vectors that should be regressed based on the feature vectors: $X_X$, $X_Y$ and $X_N$. These coefficient vectors cannot be used directly in the simulation model above. In order to derive how these coefficient vectors should be transformed, the acceleration predicted with this regression can be inserted into the equation of motion:


```python
Math(vlatex(sp.Eq(sp.UnevaluatedExpr(A)*sp.UnevaluatedExpr(coeff_matrix),sp.UnevaluatedExpr(b))))
```

The right hand side of this equation is instead expressed with matrices:


```python
eq_X_qs = sp.Eq(X_qs_*X_X, b[0])
Math(vlatex(eq_X_qs))
```


```python
eq_Y_qs = sp.Eq(X_qs_*X_Y, b[1])
Math(vlatex(eq_Y_qs))
```


```python
eq_N_qs = sp.Eq(N_qs_*X_N, b[2])
Math(vlatex(eq_N_qs))
```


```python
subs=[
    (eq_X_qs.rhs,eq_X_qs.lhs),
    (eq_Y_qs.rhs,eq_Y_qs.lhs),
    (eq_N_qs.rhs,eq_N_qs.lhs),
    
]
Math(vlatex(sp.Eq(sp.UnevaluatedExpr(A)*sp.UnevaluatedExpr(coeff_matrix),sp.UnevaluatedExpr(b.subs(subs)))))
```

So that the parameter vectors $X_{qs}$, $Y_{qs}$, $N_{qs}$ can be calculated from $A_{coeff}$, $B_{coeff}$ and $C_{coeff}$ in the following way.


```python
Math(vlatex(sp.Eq(sp.UnevaluatedExpr(A)*sp.UnevaluatedExpr(sp.matrices.MutableDenseMatrix([A_coeff,B_coeff,C_coeff])),
                  sp.UnevaluatedExpr(sp.matrices.MutableDenseMatrix([X_qs_,Y_qs_,N_qs_])))))
```

Please not that the transformed coefficient vectors: $X_{qs}$, $Y_{qs}$, $N_{qs}$ now also contain inertial forces such as as centrigal force etc. (I think). 


```python
subs = {value:key for key,value in p.items()}
A_ = A*sp.matrices.MutableDenseMatrix([A_coeff,B_coeff,C_coeff])
A_lambda=lambdify(A_.subs(subs))
```

## Load test


```python
#id=22773
#id=22616
id=22774
#id=22770


df, units, meta_data = mdl.load(id=id, dir_path='../data/processed/kalman')
df.index = df.index.total_seconds()
df = df.iloc[0:-100].copy()
df.index-=df.index[0]
df['t'] = df.index
df.sort_index(inplace=True)
df['-delta'] = -df['delta']
df['V'] = np.sqrt(df['u']**2 + df['v']**2)

df['thrust'] = df['Prop/PS/Thrust'] + df['Prop/SB/Thrust']
df['U'] = df['V']
df['beta'] = -np.arctan2(df['v'],df['u'])
```


```python
meta_data['rho']=1000
meta_data['mass'] = meta_data['Volume']*meta_data['rho']
```

# Ship parameters


```python
T_ = (meta_data.TA + meta_data.TF)/2
L_ = meta_data.lpp
m_ = meta_data.mass
rho_ = meta_data.rho
B_ = meta_data.beam
CB_ = m_/(T_*B_*L_*rho_)
I_z_ = m_*meta_data.KZZ**2
#I_z_ = 900


ship_parameters = {
        'T' : T_,
        'L' : L_,
        'CB' :CB_,
        'B' : B_,
        'rho' : rho_,
        'x_G' : 0,  # motions are expressed at CG
        'm' : m_,
        'I_z': I_z_, 
        'volume':meta_data.Volume,
    }

ps = prime_system.PrimeSystem(**ship_parameters)  # model

scale_factor = meta_data.scale_factor
ps_ship = prime_system.PrimeSystem(L=ship_parameters['L']*scale_factor, rho=meta_data['rho'])  # ship
```

# Brix parameters


```python
def calculate_prime(row, ship_parameters):
    return run(function=row['brix_lambda'], inputs=ship_parameters)


mask = df_parameters['brix_lambda'].notnull()
df_parameters.loc[mask,'brix_prime'] = df_parameters.loc[mask].apply(calculate_prime, ship_parameters=ship_parameters, axis=1)

df_parameters['brix_prime'].fillna(0, inplace=True)
#df_parameters['brix_SI'].fillna(0, inplace=True)


```


```python
df_parameters['prime'] = df_parameters['brix_prime']

def to_SI(row):
    
    prime = row['prime']
    denominator = row['denominator']
    lamda = lambdify(denominator)
    try:
        denominator_value = run(lamda, inputs=ship_parameters)
        SI = prime*denominator_value
    except:
        return None
    else:
        return SI
    
df_parameters['brix_SI'] = df_parameters.apply(to_SI, axis=1)
    
```

# Regression


```python
from statsmodels.sandbox.regression.predstd import wls_prediction_std
def show_pred(X,y,results, label):
    
    display(results.summary())
    
    X_ = X
    y_ = y
    y_pred = results.predict(X_)
    
    prstd, iv_l, iv_u = wls_prediction_std(results, exog=X_, alpha=0.05)
    #iv_l*=-1 
    #iv_u*=-1
    
    fig,ax=plt.subplots()
    ax.plot(X_.index,y_, label='Numerical gradient from model test')
    ax.plot(X_.index,y_pred, '--', label='OLS')
    ax.set_ylabel(label)
    
    ax.fill_between(X_.index, y1=iv_l, y2=iv_u, zorder=-10, color='grey', alpha=0.5, label=r'5\% confidence')
    ax.legend();
```

## N


```python
N_ = sp.symbols('N_')

diff_eq_N = regression.DiffEqToMatrix(ode=vmm.N_qs_eq.subs(N_qs,N_), 
                                      label=N_, base_features=[delta,u,v,r])
```


```python
Math(vlatex(diff_eq_N.acceleration_equation))
```


```python
X = diff_eq_N.calculate_features(data=df)
y = diff_eq_N.calculate_label(y=df['r1d'])

model_N = sm.OLS(y,X)
results_N = model_N.fit()

show_pred(X=X,y=y,results=results_N, label=r'$\dot{r}$')
```

## Y


```python
Y_ = sp.symbols('Y_')
diff_eq_Y = regression.DiffEqToMatrix(ode=vmm.Y_qs_eq.subs(Y_qs,Y_), 
                                      label=Y_, base_features=[delta,u,v,r])
```


```python
Math(vlatex(diff_eq_Y.acceleration_equation))
```


```python
X = diff_eq_Y.calculate_features(data=df)
y = diff_eq_Y.calculate_label(y=df['v1d'])


model_Y = sm.OLS(y,X)
results_Y = model_Y.fit()

show_pred(X=X,y=y,results=results_Y, label=r'$\dot{v}$')
```

## X


```python
X_ = sp.symbols('X_')
diff_eq_X = regression.DiffEqToMatrix(ode=vmm.X_qs_eq.subs(X_qs,X_), 
                                      label=X_, base_features=[delta,u,v,r,thrust])
```


```python
diff_eq_X.acceleration_equation
```


```python
X = diff_eq_X.calculate_features(data=df)
y = diff_eq_X.calculate_label(y=df['u1d'])

model_X = sm.OLS(y,X)
results_X = model_X.fit()

show_pred(X=X,y=y,results=results_X, label=r'$\dot{u}}$')
```


```python
results_summary_X = regression.results_summary_to_dataframe(results_X)
results_summary_Y = regression.results_summary_to_dataframe(results_Y)
results_summary_N = regression.results_summary_to_dataframe(results_N)
```


```python
A_coeff_ = results_summary_X['coeff']
B_coeff_ = results_summary_Y['coeff']
C_coeff_ = results_summary_N['coeff']

coeffs = run(A_lambda,A_coeff=A_coeff_.values, B_coeff=B_coeff_.values, C_coeff=C_coeff_.values, 
    **df_parameters['brix_SI'], **ship_parameters)
```


```python
results_summary_X['decoupled'] = coeffs[0][0]
results_summary_Y['decoupled'] = coeffs[1][0]
results_summary_N['decoupled'] = coeffs[2][0]

```


```python
A_coeff_*(m_-df_parameters.loc['Xudot','brix_SI'])
```


```python
A_coeff_
```

## Add the regressed parameters
Hydrodynamic derivatives that depend on acceleration cannot be obtained from the VCT regression. They are however essential if a time simulation should be conducted. These values have then been taken from Brix semi empirical formulas for the simulations below.


```python
df_parameters_all = df_parameters.copy()
for other in [results_summary_X, results_summary_Y, results_summary_N]:
    df_parameters_all = df_parameters_all.combine_first(other)

df_parameters_all.rename(columns={'decoupled':'regressed'}, inplace=True)
df_parameters_all.drop(columns=['brix_lambda'], inplace=True)

df_parameters_all['SI'] = df_parameters_all['regressed'].combine_first(df_parameters_all['brix_SI'])  # prefer regressed
df_parameters_all['SI'].fillna(0,inplace=True)
```

# Simulate forces


```python
df['X_qs_pred'] = run(function=vmm.simulator.X_qs_lambda, inputs=df, **df_parameters_all['SI'])
df['Y_qs_pred'] = run(function=vmm.simulator.Y_qs_lambda, inputs=df, **df_parameters_all['SI'])
df['N_qs_pred'] = run(function=vmm.simulator.N_qs_lambda, inputs=df, **df_parameters_all['SI'])

```


```python
X = diff_eq_X.calculate_features(data=df)
X_forces = X*results_summary_X['coeff']
X_forces.index = df.index

X = diff_eq_Y.calculate_features(data=df)
Y_forces = X*results_summary_Y['coeff']
Y_forces.index = df.index

X = diff_eq_N.calculate_features(data=df)
N_forces = X*results_summary_N['coeff']
N_forces.index = df.index

```


```python
display(px.line(X_forces, y=X_forces.columns, width=1000, height=400))
display(px.line(Y_forces, y=Y_forces.columns, width=1000, height=400))
display(px.line(N_forces, y=N_forces.columns, width=1000, height=400))
```

# Simulation

The way that the regression is formulated, inertial forces, such as centrifugal force will be included into the derivatives (I think) which means that centrifugal force : $-m \cdot r \cdot u$ will be included into $Y_{ur}$ coefficient. This coefficient is therefore not pure hydrodynamic, and can potentially be counted twice..?
The coefficients are recalculated below to avooid this:


```python
parameters=df_parameters_all['SI'].copy()
x_G_ = ship_parameters['x_G']
parameters['Xrr']+=(-m_*x_G_)
parameters['Xvr']+=(-m_)
parameters['Yur']+=m_
parameters['Nur']+=m_*x_G_


```


```python
result_vmm = vmm.simulator.simulate(df_=df, parameters=parameters, ship_parameters=ship_parameters)
```


```python
result_vmm.plot_compare()
```


```python

```

270e54a9c834ff12ae1a26de546ee81c59a5daee
ipynb
Jupyter Notebook
notebooks/11.01_regression_revisited.ipynb
martinlarsalbert/wPCC
16e0d4cc850d503247916c9f5bd9f0ddb07f8930
MIT
notebooks/11.01_regression_revisited.ipynb
martinlarsalbert/wPCC
16e0d4cc850d503247916c9f5bd9f0ddb07f8930
MIT
notebooks/11.01_regression_revisited.ipynb
martinlarsalbert/wPCC
16e0d4cc850d503247916c9f5bd9f0ddb07f8930
MIT
Qwen/Qwen-72B
1. YES
2. YES