# Simplified Linear VMM (Vessel Maneuvering Model)

# Purpose
Implementing according to:
Matusiak, Jerzy. Dynamics of a Rigid Ship. Aalto University, 2017. https://aaltodoc.aalto.fi:443/handle/123456789/24408.

# Methodology
Define the problem using SymPy

# Setup


```python
# %load imports.py
%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

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 import symbols
from src import prime_system
from src.symbols import *
from src.linear_vmm_equations import *

import src.linear_vmm_simplified_equations as eq
import src.models.linear_vmm as model
from src.visualization.plot import track_plot
import src.nonlinear_vmm_equations as nonlinear_vmm_equations

```


```python
symbols.df_parameters
```

## Linearized simplified equation of motion ($x_G=0$)

### X

Nonlinear equation of motion in X-direction:


```python
Math(vlatex(nonlinear_vmm_equations.X_eom))
```

this equation is linearized by assuming : $v=0$ and dropping $r^2$


```python
Math(vlatex(eq.X_eom))
```

### Y

Nonlinear equation of motion in Y-direction:


```python
Math(vlatex(nonlinear_vmm_equations.Y_eom))
```

this equation is linearized by assuming that $u \approx U $ where $U$ is the initial total velocity.


```python
Math(vlatex(eq.Y_eom))
```

### N

Nonlinear equation of motion in N-direction:


```python
Math(vlatex(nonlinear_vmm_equations.N_eom))
```

this equation is linearized by assuming that $u \approx U $ where $U$ is the initial total velocity.


```python
Math(vlatex(eq.N_eom))
```


```python
Math(vlatex(eq.fx_eq))
```


```python
Math(vlatex(eq.fy_eq))
```


```python
Math(vlatex(eq.mz_eq))
```


```python
Math(vlatex(eq.X_eq))
```


```python
Math(vlatex(eq.Y_eq))
```


```python
Math(vlatex(eq.N_eq))
```


```python
r1d_eq = sp.Eq(r1d, sp.solve(N_eq, r1d)[0])
```

## Steady state solution


```python
X_steady = X_eq.subs([
    (u1d,0),
    (v1d,0),
    (r1d,0),
          ])
X_steady
```


```python
Y_steady = Y_eq.subs([
    (u1d,0),
    (v1d,0),
    (r1d,0),
          ])
Y_steady
```


```python
N_steady = N_eq.subs([
    (u1d,0),
    (v1d,0),
    (r1d,0),
          ])
N_steady
```


```python
eqs = [Y_steady,N_steady]
solution = sp.solve(eqs,v,r, dict=True)
solution[0][v]
```


```python
solution[0][r]
```


```python
subs = {value:key for key,value in eq.p.items()}
#u_lambda = lambdify(solution[0][u].subs(subs))
v_lambda = lambdify(solution[0][v].subs(subs))
r_lambda = lambdify(solution[0][r].subs(subs))
```

## Prime system


```python
prime_system.df_prime.loc[['denominator']]
```

## Ship parameters


```python
T_ =10
L_ = 200
CB_ = 0.7
B_ = 30
rho_ = 1025
m_ = T_*B_*L_*CB_*rho_

ship_parameters = {
        'T' : T_,
        'L' : L_,
        'CB' :CB_,
        'B' : B_,
        'rho' : rho_,
        'x_G' : 0,
        'm' : m_,
        'I_z': 0.2*m_*L_**2, 
    }

ps = prime_system.PrimeSystem(**ship_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,'prime'] = df_parameters.loc[mask].apply(calculate_prime, ship_parameters=ship_parameters, axis=1)
df_parameters.loc['Ydelta','prime'] = 0.0005  # Just guessing
df_parameters.loc['Ndelta','prime'] = -df_parameters.loc['Ydelta','prime']/4  # Just guessing
```


```python
df_parameters['prime'].fillna(0, inplace=True)
df_parameters
```

## Steady turning


```python
inputs = dict(df_parameters['prime'])
inputs.update(ps.prime(ship_parameters))

U = 10
delta = np.deg2rad(35)
inputs['delta'] = delta
inputs['u'] = inputs['U'] =  ps._prime(U,'linear_velocity', U=U)


v_steady_prime = run(function=v_lambda, inputs=inputs)
r_steady_prime = run(function=r_lambda, inputs=inputs)

v_steady = ps._unprime(v_steady_prime,'linear_velocity', U=U)
r_steady = ps._unprime(r_steady_prime,'angular_velocity' ,U=U)

```


```python
v_steady
```


```python
r_steady
```


```python
model.u1d_lambda
```


```python
model.v1d_lambda
```


```python
states = {
    'u' : ps._prime(U,'linear_velocity',U=U),
    'v' : ps._prime(v_steady,'linear_velocity',U=U),
    'r' : ps._prime(r_steady,'angular_velocity',U=U),
    #'u1d' : 0,
    #'v1d' : 0,
    #'r1d' : 0,
    
}

control = {
    'delta' : delta,
    'U' : ps._prime(U,'linear_velocity',U=U)
}

inputs = dict(df_parameters['prime'])
inputs.update(ps.prime(ship_parameters))

inputs.update(states)
inputs.update(control)

run(function=model.v1d_lambda, inputs=inputs)
```


```python
run(function=model.u1d_lambda, inputs=inputs)
```


```python
run(function=model.v1d_lambda, inputs=inputs)
```


```python
run(function=model.r1d_lambda, inputs=inputs)
```


```python
t_max = 2*np.pi/r_steady  # Turn one lap
t = np.arange(0,t_max,0.01)


control = {
    'delta' : delta,
}

parameters = dict(df_parameters['prime'])

u_steady = np.sqrt(U**2 - v_steady**2)
y0 = {
    'u' : u_steady, 
    'v' : v_steady,
    'r' : r_steady,
    'x0' : 0,
    'y0' : 0,
    'psi' : 0,
    }

#solution = model.simulate(y0=y0, t=t, df_parameters=df_parameters, df_ship_parameters=df_ship_parameters, control=control, rtol=1e-3, atol=1e-3)
solution = model.simulate(y0=y0, t=t, df_parameters=df_parameters, ship_parameters=ship_parameters, control=control)

columns = ['u','v','r','x0','y0','psi']
df_result_prime = pd.DataFrame(data=solution.y.T, columns=columns)
df_result_prime.index=t[0:len(df_result_prime)]

#df_result_prime['U'] = np.sqrt(df_result_prime['u']**2 + df_result_prime['v']**2)

df_result_prime.plot(y='u')
df_result_prime.plot(y='v')
df_result_prime.plot(y='r')


```


```python
y0
```


```python
solution
```


```python
#units = {key:value[1] for key,value in zip(columns,y0)}
df_result = ps.unprime(values=df_result_prime, U=U)


df_result['beta'] = -np.arctan2(df_result['v'],df_result['u'])
```


```python
df_result.head()
```


```python
track_plot(df=df_result, lpp=ship_parameters['L'], beam=ship_parameters['B'])

df_result.plot(y='u')
df_result.plot(y='v')
df_result.plot(y='r')
df_result.plot(y='psi')
df_result.plot(y='beta')


```

## Steady turnings


```python
inputs = dict(df_parameters['prime'])
inputs.update(ps.prime(ship_parameters))

U = 10
delta = np.deg2rad(35)
inputs['delta'] = delta
inputs['u'] = inputs['U'] =  ps._prime(U,'linear_velocity', U=U)

s = pd.Series(inputs)
N=10
data = np.tile(s.values,(N,1))
df_variation = pd.DataFrame(data, columns=s.index)
df_variation.head()
df_variation['delta'] = np.linspace(0.01,0.3,N)

df_variation['v'] = run(function=v_lambda, inputs=df_variation)
df_variation['r'] = run(function=r_lambda, inputs=df_variation)

```


```python
df_variation.plot(x='delta', y='r')
df_variation.plot(x='delta', y='v')


```


```python
sp.solve(eq.Y_eq,v1d)[0]
```


```python
sp.solve(eq.N_eq,r1d)[0]
```


```python
sp.solve(eq.X_eq,u1d)[0]
```

85060bf0d93cb4e99e16e7642e22b0bab29f0c8c
ipynb
Jupyter Notebook
notebooks/06.01_linear_simplified_VMM.ipynb
martinlarsalbert/wPCC
16e0d4cc850d503247916c9f5bd9f0ddb07f8930
MIT
notebooks/06.01_linear_simplified_VMM.ipynb
martinlarsalbert/wPCC
16e0d4cc850d503247916c9f5bd9f0ddb07f8930
MIT
notebooks/06.01_linear_simplified_VMM.ipynb
martinlarsalbert/wPCC
16e0d4cc850d503247916c9f5bd9f0ddb07f8930
MIT
Qwen/Qwen-72B
1. YES
2. YES