Merge pull request #70 from ShipMMG/dev

v0.0.8

Author: taiga4112

Project.toml

@@ -1,22 +1,29 @@
 name = "ShipMMG"
 uuid = "37f2b0bf-0c13-4883-8808-e75eb56597e7"
 authors = ["Taiga MITSUYUKI <[email protected]>"]
-version = "0.0.7"
+version = "0.0.8"
 
 [deps]
 AdvancedPS = "576499cb-2369-40b2-a588-c64705576edc"
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Dierckx = "39dd38d3-220a-591b-8e3c-4c3a8c710a94"
 DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 KernelDensity = "5ab0869b-81aa-558d-bb23-cbf5423bbe9b"
 Libtask = "6f1fad26-d15e-5dc8-ae53-837a1d7b8c9f"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ParameterizedFunctions = "65888b18-ceab-5e60-b2b9-181511a3b968"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
+ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
+Revise = "295af30f-e4ad-537b-8983-00126c2a3abe"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
+StatsPlots = "f3b207a7-027a-5e70-b257-86293d7955fd"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 
 [compat]
-AdvancedPS = "^0.4"
+AdvancedPS = "^0.6"
 Dierckx = "0.5"
 DifferentialEquations = "^7"
 Distributions = "^0.25"
@@ -25,7 +32,7 @@ KernelDensity = "^0.6"
 Libtask = "^0.8"
 ParameterizedFunctions = "^5.12"
 Parameters = "0.12"
-Turing = "^0.24"
+Turing = "^0.36"
 julia = "^1.6"
 
 [extras]

src/ShipMMG.jl

@@ -7,6 +7,7 @@ using Parameters
 using Distributions
 using Turing
 using ForwardDiff
+using LinearAlgebra
 
 include("data.jl")
 export calc_position,

src/kt.jl

@@ -90,7 +90,7 @@ function kt_simulate(
     spl_δ = Spline1D(time_list, δ_list)
 
     X0 = [u0; v0; r0; x0; y0; Ψ0; δ_list[1]]
-    p = [K, T, spl_δ]
+    p = (K, T, spl_δ)
     prob = ODEProblem(kt_model!, X0, (time_list[1], time_list[end]), p)
     sol = solve(prob, algorithm, reltol=reltol, abstol=abstol)
     sol_timelist = sol(time_list)
@@ -333,7 +333,7 @@ function create_model_for_mcmc_sample_kt(
     u0 = [r_obs[1]; δ_obs[1]]
     K_start = 0.10
     T_start = 60.0
-    p = [K_start, T_start]
+    p = (K_start, T_start)
     prob1 = ODEProblem(KT!, u0, (time_obs[1], time_obs[end]), p)
 
     # create probabilistic model
@@ -342,7 +342,7 @@ function create_model_for_mcmc_sample_kt(
         K ~ K_prior_dist
         T ~ T_prior_dist
 
-        p = [K, T]
+        p = (K, T)
         prob = remake(prob1, p=p)
         sol = solve(prob, Tsit5())
         predicted = sol(time_obs)

src/mmg.jl

@@ -1103,7 +1103,7 @@ function create_model_for_mcmc_sample_mmg(
     N_vrr_dash_start = 0.055
     N_rrr_dash_start = -0.013
 
-    p = [
+    p = (
         R_0_dash_start,
         X_vv_dash_start,
         X_vr_dash_start,
@@ -1121,7 +1121,7 @@ function create_model_for_mcmc_sample_mmg(
         N_vvr_dash_start,
         N_vrr_dash_start,
         N_rrr_dash_start,
-    ]
+    )
 
     u0 = 2.29 * 0.512
     v0 = 0.0
@@ -1152,7 +1152,7 @@ function create_model_for_mcmc_sample_mmg(
         N_vrr_dash ~ N_vrr_dash_prior_dist
         N_rrr_dash ~ N_rrr_dash_prior_dist
 
-        p = [
+        p = (
             R_0_dash,
             X_vv_dash,
             X_vr_dash,
@@ -1170,7 +1170,259 @@ function create_model_for_mcmc_sample_mmg(
             N_vvr_dash,
             N_vrr_dash,
             N_rrr_dash,
-        ]
+        )
+        prob = remake(prob1, p=p)
+        sol = solve(prob, solver, abstol=abstol, reltol=reltol)
+        predicted = sol(time_obs)
+        for i in eachindex(predicted)
+            obs[1][i] ~ Normal(predicted[i][1], σ_u) # u
+            obs[2][i] ~ Normal(predicted[i][2], σ_v) # v
+            obs[3][i] ~ Normal(predicted[i][3], σ_r) # r
+        end
+    end
+
+    return fitMMG(time_obs, [u_obs, v_obs, r_obs], prob1)
+end
+
+function create_model_for_mcmc_sample_mmg(
+    data::ShipData,
+    basic_params::Mmg3DofBasicParams,
+    k_0,
+    k_1,
+    k_2,
+    multivariate_prior_dist;
+    ρ=1025.0,
+    σ_u_prior_dist=Uniform(0.00, 0.20),
+    σ_v_prior_dist=Uniform(0.00, 0.20),
+    σ_r_prior_dist=Uniform(0.00, 0.20),
+    solver=Tsit5(),
+    abstol=1e-6,
+    reltol=1e-3
+)
+    time_obs = data.time
+    u_obs = data.u
+    v_obs = data.v
+    r_obs = data.r
+    δ_obs = data.δ
+    n_p_obs = data.n_p
+
+    @unpack L_pp,
+    B,
+    d,
+    x_G,
+    D_p,
+    m,
+    I_zG,
+    A_R,
+    η,
+    m_x,
+    m_y,
+    J_z,
+    f_α,
+    ϵ,
+    t_R,
+    x_R,
+    a_H,
+    x_H,
+    γ_R_minus,
+    γ_R_plus,
+    l_R,
+    κ,
+    t_P,
+    w_P0,
+    x_P = basic_params
+
+    # create sytem model
+    spl_δ = Spline1D(time_obs, δ_obs)
+    spl_n_p = Spline1D(time_obs, n_p_obs)
+    function MMG!(dX, X, p, t)
+        u, v, r, δ, n_p = X
+        R_0_dash,
+        X_vv_dash,
+        X_vr_dash,
+        X_rr_dash,
+        X_vvvv_dash,
+        Y_v_dash,
+        Y_r_dash,
+        Y_vvv_dash,
+        Y_vvr_dash,
+        Y_vrr_dash,
+        Y_rrr_dash,
+        N_v_dash,
+        N_r_dash,
+        N_vvv_dash,
+        N_vvr_dash,
+        N_vrr_dash,
+        N_rrr_dash = p
+
+        U = sqrt(u^2 + (v - r * x_G)^2)
+
+        β = 0.0
+        if U == 0.0
+            β = 0.0
+        else
+            β = asin(-(v - r * x_G) / U)
+        end
+
+        v_dash = 0.0
+        if U == 0.0
+            v_dash = 0.0
+        else
+            v_dash = v / U
+        end
+
+        r_dash = 0.0
+        if U == 0.0
+            r_dash = 0.0
+        else
+            r_dash = r * L_pp / U
+        end
+
+        w_P = w_P0 * exp(-4.0 * (β - x_P * r_dash)^2)
+
+        J = 0.0
+        if n_p == 0.0
+            J = 0.0
+        else
+            J = (1 - w_P) * u / (n_p * D_p)
+        end
+        K_T = k_0 + k_1 * J + k_2 * J^2
+        β_R = β - l_R * r_dash
+        γ_R = γ_R_minus
+
+        if β_R < 0.0
+            γ_R = γ_R_minus
+        else
+            γ_R = γ_R_plus
+        end
+
+        v_R = U * γ_R * β_R
+
+        u_R = 0.0
+        if J == 0.0
+            u_R = sqrt(η * (κ * ϵ * 8.0 * k_0 * n_p^2 * D_p^4 / pi)^2)
+        else
+            u_R =
+                u *
+                (1.0 - w_P) *
+                ϵ *
+                sqrt(η * (1.0 + κ * (sqrt(1.0 + 8.0 * K_T / (pi * J^2)) - 1))^2 + (1 - η))
+        end
+
+        U_R = sqrt(u_R^2 + v_R^2)
+        α_R = δ - atan(v_R, u_R)
+        F_N = 0.5 * A_R * ρ * f_α * (U_R^2) * sin(α_R)
+
+        X_H = (
+            0.5 *
+            ρ *
+            L_pp *
+            d *
+            (U^2) *
+            (
+                -R_0_dash +
+                X_vv_dash * (v_dash^2) +
+                X_vr_dash * v_dash * r_dash +
+                X_rr_dash * (r_dash^2) +
+                X_vvvv_dash * (v_dash^4)
+            )
+        )
+        X_R = -(1.0 - t_R) * F_N * sin(δ)
+        X_P = (1.0 - t_P) * ρ * K_T * n_p^2 * D_p^4
+        Y_H = (
+            0.5 *
+            ρ *
+            L_pp *
+            d *
+            (U^2) *
+            (
+                Y_v_dash * v_dash +
+                Y_r_dash * r_dash +
+                Y_vvv_dash * (v_dash^3) +
+                Y_vvr_dash * (v_dash^2) * r_dash +
+                Y_vrr_dash * v_dash * (r_dash^2) +
+                Y_rrr_dash * (r_dash^3)
+            )
+        )
+        Y_R = -(1 + a_H) * F_N * cos(δ)
+        N_H = (
+            0.5 *
+            ρ *
+            (L_pp^2) *
+            d *
+            (U^2) *
+            (
+                N_v_dash * v_dash +
+                N_r_dash * r_dash +
+                N_vvv_dash * (v_dash^3) +
+                N_vvr_dash * (v_dash^2) * r_dash +
+                N_vrr_dash * v_dash * (r_dash^2) +
+                N_rrr_dash * (r_dash^3)
+            )
+        )
+        N_R = -(x_R + a_H * x_H) * F_N * cos(δ)
+        dX[1] = du = ((X_H + X_R + X_P) + (m + m_y) * v * r + x_G * m * (r^2)) / (m + m_x)
+        dX[2] =
+            dv =
+                (
+                    (x_G^2) * (m^2) * u * r - (N_H + N_R) * x_G * m +
+                    ((Y_H + Y_R) - (m + m_x) * u * r) * (I_zG + J_z + (x_G^2) * m)
+                ) / ((I_zG + J_z + (x_G^2) * m) * (m + m_y) - (x_G^2) * (m^2))
+        dX[3] = dr = (N_H + N_R - x_G * m * (dv + u * r)) / (I_zG + J_z + (x_G^2) * m)
+        dX[4] = dδ = derivative(spl_δ, t)
+        dX[5] = dn_p = derivative(spl_n_p, t)
+    end
+
+    R_0_dash_start = 0.022
+    X_vv_dash_start = -0.040
+    X_vr_dash_start = 0.002
+    X_rr_dash_start = 0.011
+    X_vvvv_dash_start = 0.771
+    Y_v_dash_start = -0.315
+    Y_r_dash_start = 0.083
+    Y_vvv_dash_start = -1.607
+    Y_vvr_dash_start = 0.379
+    Y_vrr_dash_start = -0.391
+    Y_rrr_dash_start = 0.008
+    N_v_dash_start = -0.137
+    N_r_dash_start = -0.049
+    N_vvv_dash_start = -0.030
+    N_vvr_dash_start = -0.294
+    N_vrr_dash_start = 0.055
+    N_rrr_dash_start = -0.013
+
+    p = (
+        R_0_dash_start,
+        X_vv_dash_start,
+        X_vr_dash_start,
+        X_rr_dash_start,
+        X_vvvv_dash_start,
+        Y_v_dash_start,
+        Y_r_dash_start,
+        Y_vvv_dash_start,
+        Y_vvr_dash_start,
+        Y_vrr_dash_start,
+        Y_rrr_dash_start,
+        N_v_dash_start,
+        N_r_dash_start,
+        N_vvv_dash_start,
+        N_vvr_dash_start,
+        N_vrr_dash_start,
+        N_rrr_dash_start,
+    )
+
+    u0 = 2.29 * 0.512
+    v0 = 0.0
+    r0 = 0.0
+    X0 = [u_obs[1]; v_obs[1]; r_obs[1]; δ_obs[1]; n_p_obs[1]]
+    prob1 = ODEProblem(MMG!, X0, (time_obs[1], time_obs[end]), p)
+
+    # create probabilistic model
+    @model function fitMMG(time_obs, obs, prob1)
+        σ_u ~ σ_u_prior_dist
+        σ_v ~ σ_v_prior_dist
+        σ_r ~ σ_r_prior_dist
+        p ~ multivariate_prior_dist;
         prob = remake(prob1, p=p)
         sol = solve(prob, solver, abstol=abstol, reltol=reltol)
         predicted = sol(time_obs)

src/mmg_wind.jl

@@ -1142,7 +1142,7 @@ function create_model_for_mcmc_sample_mmg(
     N_vrr_dash_start = 0.055
     N_rrr_dash_start = -0.013
 
-    p = [
+    p = (
         R_0_dash_start,
         X_vv_dash_start,
         X_vr_dash_start,
@@ -1160,7 +1160,7 @@ function create_model_for_mcmc_sample_mmg(
         N_vvr_dash_start,
         N_vrr_dash_start,
         N_rrr_dash_start,
-    ]
+    )
 
     u0 = 2.29 * 0.512
     v0 = 0.0
@@ -1191,7 +1191,7 @@ function create_model_for_mcmc_sample_mmg(
         N_vrr_dash ~ N_vrr_dash_prior_dist
         N_rrr_dash ~ N_rrr_dash_prior_dist
 
-        p = [
+        p = (
             R_0_dash,
             X_vv_dash,
             X_vr_dash,
@@ -1209,7 +1209,7 @@ function create_model_for_mcmc_sample_mmg(
             N_vvr_dash,
             N_vrr_dash,
             N_rrr_dash,
-        ]
+        )
         prob = remake(prob1, p=p)
         sol = solve(prob, solver, abstol=abstol, reltol=reltol)
         predicted = sol(time_obs)