问题3

# TODO import

import re
import os
import sys
import hmz
import pathlib
import mitosheet
import numpy as np
import pandas as pd
import matlab.engine

import scipy
from scipy.integrate import odeint
from scipy.optimize import minimize

import time
from time import time, sleep

import copy
import random

import sympy
from sympy import limit
from sympy import diff
from sympy import integrals

import sklearn
import graphviz
from sklearn import tree
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import train_test_split
from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error as MSE
from sklearn.metrics import mean_absolute_error as MAE
from sklearn.metrics import classification_report, roc_auc_score

import sko
from sko.GA import GA

import numba
from numba import jit

import plotly
import plotly.express as px
import plotly.graph_objects as go
import plotly.figure_factory as ff
plotly.offline.init_notebook_mode()

import cufflinks as cf
cf.set_config_file(
    offline=True, 
    world_readable=True,
    theme='pearl',  # cf.getThemes()
)

import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei']  # KaiTi
plt.rcParams['axes.unicode_minus'] = False

from IPython.core.interactiveshell import InteractiveShell 
InteractiveShell.ast_node_interactivity = 'all'
# InteractiveShell.ast_node_interactivity = 'last'

import cv2 as cv

# import torch
# import torchvision
# import torch.nn as nn
# import torch.nn.functional as F
# import torch.utils.data as Data
# from torch.utils.data import DataLoader
# from torch.utils.data.dataset import Dataset

import pylatex
import latexify

import warnings
warnings.filterwarnings("ignore")

# TODO 日志、计时
# from colorama.Fore import RED, RESET
from colorama import Fore
import logging
fmt = '%(asctime)s - %(levelname)8s - %(message)s'
formatter = logging.Formatter(fmt)
handler_control = logging.StreamHandler()  # stdout to console
handler_control.setLevel('INFO')  # 设置 INFO 级别
handler_control.setFormatter(formatter)

logger = logging.getLogger()
# logger.setLevel('INFO')
logger.addHandler(handler_control)

def timeit(text):
    def func_deco(func):
        """ 用来计时的装饰器函数 """
        def func_wrapper(*args, **kwargs):
            from time import time
            t0 = time()
#             logging.info(text + "开始计时")
            print(Fore.RED, text + "开始计时：", Fore.RESET)
            res = func(*args, **kwargs)
            t1 = time()
#             logging.info(text + "用时：" + str(t1 - t0) + "s")
            print(Fore.RED, text + "结束计时，用时：", str(t1 - t0), "s", Fore.RESET)
            return res
        return func_wrapper
    return func_deco

# TODO DIR
ROOTDIR = pathlib.Path(os.path.abspath('.'))
IMG_HTML = ROOTDIR / 'img-html'
IMG_SVG = ROOTDIR / 'img-svg'
DATA_RAW = ROOTDIR / 'data-raw'
DATA_COOKED = ROOTDIR / 'data-processed'

# TODO 附件4参数
浮子质量 = 4866  # kg
浮子底半径 = 1  # m
浮子圆柱部分高度 = 3  # m
浮子圆锥部分高度 = 0.8  # m
振子质量 = 2433  # kg
振子半径 = 0.5  # m
振子高度 = 0.5  # m
海水的密度 = 1025  # kg/m^3
重力加速度 = 9.8  # m/s^2
弹簧刚度 = 80000  # N/m
弹簧原长 = 0.5  # m
扭转弹簧刚度 = 250000  # N·m
静水恢复力矩系数 = 8890.7  # N·m

# TODO 附件3参数

class question1234:
    """设置具体问题几的参数"""
    def __init__(self, question):
        global 入射波浪频率
        global 垂荡附加质量
        global 纵摇附加转动惯量
        global 垂荡兴波阻尼系数
        global 纵摇兴波阻尼系数
        global 垂荡激励力振幅
        global 纵摇激励力矩振幅
        global 波浪频率
        global 波浪周期
        
        if question == 1:
            # 问题1：参数
            # 纵摇附加转动惯量 = 6779.315  # kg·m^2
            # 纵摇兴波阻尼系数 = 151.4388  # N·m·s
            # 纵摇激励力矩振幅 = 1230  # N·m
            纵摇附加转动惯量 = None  # 问题1未使用，为避免使用/误用，初始化为 None
            纵摇兴波阻尼系数 = None  # 问题1未使用，为避免使用/误用，初始化为 None
            纵摇激励力矩振幅 = None  # 问题1未使用，为避免使用/误用，初始化为 None
            入射波浪频率 = 1.4005  # s^{-1}
            垂荡附加质量 = 1335.535  # kg
            垂荡兴波阻尼系数 = 656.3616  # N·s/m
            垂荡激励力振幅 = 6250  # N
            波浪频率 = 入射波浪频率
            波浪周期 = 1 / 波浪频率
        elif question == 2:
            # 问题2：参数
            # 纵摇附加转动惯量 = 7131.29
            # 纵摇兴波阻尼系数 = 2992.724
            # 纵摇激励力矩振幅 = 2560
            纵摇附加转动惯量 = None  # 问题2未使用，为避免使用/误用，初始化为 None
            纵摇兴波阻尼系数 = None  # 问题2未使用，为避免使用/误用，初始化为 None
            纵摇激励力矩振幅 = None  # 问题2未使用，为避免使用/误用，初始化为 None
            入射波浪频率 = 2.2143
            垂荡附加质量 = 1165.992
            垂荡兴波阻尼系数 = 167.8395
            垂荡激励力振幅 = 4890
            波浪频率 = 入射波浪频率
            波浪周期 = 1 / 波浪频率
        elif question == 3:
            # 问题3：参数
            入射波浪频率 = 1.7152
            垂荡附加质量 = 1028.876
            纵摇附加转动惯量 = 7001.914
            垂荡兴波阻尼系数 = 683.4558
            纵摇兴波阻尼系数 = 654.3383
            垂荡激励力振幅 = 3640
            纵摇激励力矩振幅 = 1690
            波浪频率 = 入射波浪频率
            波浪周期 = 1 / 波浪频率
        elif question == 4:
            # 问题4：参数
            入射波浪频率 = 1.9806
            垂荡附加质量 = 1091.099
            纵摇附加转动惯量 = 7142.493
            垂荡兴波阻尼系数 = 528.5018
            纵摇兴波阻尼系数 = 1655.909
            垂荡激励力振幅 = 1760
            纵摇激励力矩振幅 = 2140
            波浪频率 = 入射波浪频率
            波浪周期 = 1 / 波浪频率
        
        return None

class trange:
    """设置时间区间和间隔的参数"""
    def __init__(self, left, right, step):
        global t_left
        global t_right
        global t_step
        
        t_left = left
        t_right = right
        t_step = step
        
        return None

_ = question1234(3)
_ = trange(0, 200, 0.2)

# TODO 跟变量有关的参数函数（这个后面有用，但是用处不大）

def S浮子底面积_func(r=浮子底半径):
    return np.pi * r**2
S浮子底面积 = S浮子底面积_func()

def V排_func(h, pprint=True):
    """
    :param h: 圆柱壳体的入水深度
    :param pprint: 是否打印状态
    :return: V排 (m^3)
    """
    if h >= 0:
        print("圆锥壳体完全浸没")
        V排 = (1/3 * S浮子底面积 * 浮子圆锥部分高度) + (S浮子底面积 * h)
    else:
        print("圆锥壳体漂浮")
        depth = 浮子圆锥部分高度 + h
        r = 浮子底半径 * depth / 浮子圆锥部分高度
        V排 = 1/3 * S浮子底面积_func(r) * depth
    return V排
# print("浮子入水体积：", V排_func(3))
# print("浮子入水体积：", V排_func(2.4147))
# print("浮子入水体积：", V排_func(0))
# print("浮子入水体积：", V排_func(-0.001))
# print("浮子入水体积：", V排_func(-0.8))

def F静水恢复力_func(h, pprint=False):
    """ 类似(就是)浮力 方向向上
    :param h: 圆柱壳体的入水深度
    :param pprint: 是否打印状态
    :return: F静水恢复力 (N)
    """
    F静水恢复力 = 海水的密度 * 重力加速度 * V排_func(h, pprint)
    return F静水恢复力
# F静水恢复力_func(2.4147)

def F兴波阻尼力_func(v, k=垂荡兴波阻尼系数):
    """ 方向同速度方向
    :param v: 速度
    :return:
    """
    F兴波阻尼力 = k * v
    return F兴波阻尼力

def F波浪激励力_func(t, omega=入射波浪频率, f=垂荡激励力振幅):
    """ 方向向上
    :param t: 时间
    :return: F波浪激励力 (N)
    """
    F波浪激励力 = f * np.cos(omega * t)
    return F波浪激励力
# F波浪激励力_func(0)

def F附加惯性力_func(m=垂荡附加质量, g=重力加速度):
    """ 方向向下 """
    F附加惯性力 = m * g
    return F附加惯性力
F附加惯性力 = F附加惯性力_func()

def F重力_func(m=浮子质量+振子质量, g=重力加速度):
    """ 方向向下 """
    F重力 = m * g
    return F重力
F重力 = F重力_func()

def c直线阻尼器的阻尼系数_func1():
    c直线阻尼器的阻尼系数 = 10000  # N·s/m
    return c直线阻尼器的阻尼系数

def c直线阻尼器的阻尼系数_func2(v浮子, v振子, k=10000, a=0.5):
    c直线阻尼器的阻尼系数 = k * abs(v浮子 - v振子)**a  # N·s/m
    return c直线阻尼器的阻尼系数

浮子和振子的垂荡和纵摇的微分方程组

$$\begin{align} \nonumber m_2 \frac{d^2 X_2(t)}{dt^2} cos(\theta_2(t)) + c \frac{d X_2(t)}{dt} cos(\theta_2(t)) + k X_2(t) cos(\theta_2(t)) & = c \frac{d X_1(t)}{dt} + k X_1(t)\\ \nonumber m_1' \frac{d^2 X_1(t)}{dt^2} + m_2 \frac{d^2 X_2(t)}{dt^2} cos(\theta_2(t)) & = F(t)\\ \nonumber j_2 \frac{d^2 \theta_2(t)}{dt^2} + C \frac{d \theta_2(t)}{dt} + K \theta_2(t) & = C \frac{d \theta_1(t)}{dt} + K \theta_1(t)\\ \nonumber j_1' \frac{d^2 \theta_1(t)}{dt^2} + j_2 \frac{d^2 \theta_2(t)}{dt^2} & = M(t)\\ \end{align}$$

$m$ 质量

$j$ 转动惯量

$c$ 力的阻尼系数

$C$ 力矩的阻尼系数

$k$ 刚度系数

$K$ 扭转弹簧的刚度系数

$F(t) = F_{波浪激励力} - F_{兴波阻尼力} - F_{静水恢复力}$

$= f cos(\omega t) + k_{垂荡兴波阻尼系数} v + \rho g S h$

$= f cos(\omega t) + k_{垂荡兴波阻尼系数} \frac{d X_1(t)}{dt} + k_{系数} X_1(t)$

$= f cos(\omega t) + k_1 y_2 + k_2 y_1$

$M(t) = M_{波浪激励力矩} - M_{兴波阻尼力矩} - M_{静水恢复力矩}$

$ = L cos(\omega t) + K_{纵摇兴波阻尼力矩系数} \frac{d \theta_1(t)}{dt} + K_{静水恢复力矩系数} \theta_1(t)$

$ = L cos(\omega t) + K_1 y_6 + K_2 y_5$

$y' = \left[ \begin{matrix} y_1' \\ y_2' \\ y_3' \\ y_4' \\y_5' \\ y_6' \\ y_7' \\ y_8' \end{matrix} \right] = \left[ \begin{matrix} \dot{X_1} \\ \ddot{X_1} \\ \dot{X_2} \\ \ddot{X_2} \\ \dot{\theta_1} \\ \ddot{\theta_1} \\ \dot{\theta_2} \\ \ddot{\theta_2} \end{matrix} \right] = \left[ \begin{matrix} \dot{X_1} \\ (f cos(\omega t) + k_1 \dot{X_1} + k_2 X_1 - m_2 \ddot{X_2} cos(\theta_2)) / m_1' \\ \dot{X_2} \\ (c(\dot{X_1} - \dot{X_2} cos(\theta_2)) + k(X_1 - X_2 cos(\theta_2))) / m_2 cos(\theta_2) \\ \dot{\theta_1} \\ (L cos(\omega t) + K_1 \dot{\theta_1} + K_2 \theta_1 - j_2 \ddot{\theta_2}) / j_1' \\ \dot{\theta_2} \\ (C(\dot{\theta_1} - \dot{\theta_2}) + K(\theta_1 - \theta_2)) / j_2 \\ \end{matrix} \right]$

$= \left[ \begin{matrix} y_2 \\ (f cos(\omega t) + k_1 y_2 + k_2 y_1 - m_2 \ddot{X_2} cos(y_7)) / m_1' \\ y_4 \\ (c(y_2 - y_4 cos(y_7)) + k(y_1 - y_3 cos(y_7))) / m_2 cos(y_7) \\ y_6 \\ (L cos(\omega t) + K_1 y_6 + K_2 y_5 - j_2 \ddot{\theta_2}) / j_1' \\ y_8 \\ (C(y_6 - y_8) + K(y_5 - y_7)) / j_2 \\ \end{matrix} \right]$

$j_1 = 33000$

$j_2 = \left(\left(0.469588 + 0.5 + X_1 cos(\theta_2) - X_2 \left(t\right)\right)^3 - \left(0.469558 + X_1 cos(\theta_2) - X_2 \left(t\right)\right)^3\right) m_2 / 1.5$

$j_2 = ((0.969588 + y_3)^3 - (0.469558 + y_3)^3) m_2 / 1.5$

浮子与振子的垂荡位移与速度和纵摇角位移与角速度

# TODO 问题3：参数
m1, m2, me = 浮子质量, 振子质量, 垂荡附加质量
j2_func = lambda y3: ((0.969588 + y3)**3 - (0.469558 + y3)**3) * m2 / 1.5  # y3
j1, je = 33000, 纵摇附加转动惯量
omega = 入射波浪频率

c = 10000
f = 垂荡激励力振幅
k = 弹簧刚度
k1 = -垂荡兴波阻尼系数
k2 = -海水的密度 * 重力加速度 * S浮子底面积

C = 1000
L = 纵摇激励力矩振幅
K = 扭转弹簧刚度
K1 = -纵摇兴波阻尼系数
K2 = -静水恢复力矩系数

def ode_func(y, t):
    dy1 = y[2-1]
    dy3 = y[4-1]
    dy4 = c * (y[2-1] - y[4-1] * np.cos(y[7-1])) + k * (y[1-1] - y[3-1] * np.cos(y[7-1]))
    dy2 = f * np.cos(omega * t) + k1 * y[2-1] + k2 * y[1-1] - dy4
    dy2 /= m1 + me
    dy4 /= m2 * np.cos(y[7-1])
    
    j2 = j2_func(y[3-1])
    
    dy5 = y[6-1]
    dy7 = y[8-1]
    dy8 = C * (y[6-1] - y[8-1]) + K * (y[5-1] - y[7-1])
    dy6 = L * np.cos(omega * t) + K1 * y[6-1] + K2 * y[5-1] - dy8
    dy6 /= j1 + je
    dy8 /= j2
    return [dy1, dy2, dy3, dy4, dy5, dy6, dy7, dy8]

t = np.linspace(t_left, t_right, num=int(t_right / t_step) + 1)
y0 = [0 for _ in range(8)]
res3 = odeint(ode_func, y0, t)

保存结果

# TODO 保存结果

def get_result3_df(t=t, result3=res3):
    columns = ['时间 (s)', '浮子垂荡位移 (m)', '浮子垂荡速度 (m/s)', '振子垂荡位移 (m)', '振子垂荡速度 (m/s)', 
               '浮子纵摇角位移 (rad)', '浮子纵摇角速度 (rad/s)', '振子纵摇角位移 (rad)', '振子纵摇角速度 (rad/s)']
    shijian = pd.DataFrame(t, columns=columns[:1])
    result3 = pd.DataFrame(result3, columns=columns[1:])
    result3_df = pd.concat([shijian, result3], axis=1)
    # 时间，浮子位移，浮子速度，振子位移，振子速度，浮子角位移，浮子角速度，振子角位移，振子角速度
    # -> 时间，浮子位移，浮子速度，浮子角位移，浮子角速度，振子位移，振子速度，振子角位移，振子角速度
    result3_df = result3_df.iloc[:, [0, 1, 2, 5, 6, 3, 4, 7, 8]]
    return result3_df

def save_result3_df(result3_df=get_result3_df(), save=True):
    # 前 40 周期，间隔 0.2s
    result3 = result3_df.iloc[:int(40 * 波浪周期 / t_step) + 1:int(0.2 / t_step), :]
    if save:
        result3.to_csv('result3.csv', encoding='utf_8_sig')
    
    #  10 s、20 s、40 s、60 s、100 s
    idx = list(map(lambda x: int(x / t_step), [10, 20, 40, 60, 100]))
    result3_paper = result3_df.iloc[idx, :]
    if save:
        result3_paper.to_csv('result3-paper.csv', encoding='utf_8_sig')
    
    return result3_df

result3_df = save_result3_df()

绘图

result3_df.iplot(x='时间 (s)')

# 位移 速度
def plot_plotly_xv(title, t=t, y=res3, svg_name=None):
    trace1 = go.Scatter(x=t, y=y[:, 0], name="$浮子垂荡位移 ~ X_1$", yaxis='y1')
    trace2 = go.Scatter(x=t, y=y[:, 1], name="$浮子垂荡速度 ~ V_1$", yaxis='y2')
    trace3 = go.Scatter(x=t, y=y[:, 2], name="$振子垂荡位移 ~ X_2$", yaxis='y1')
    trace4 = go.Scatter(x=t, y=y[:, 3], name="$振子垂荡速度 ~ V_2$", yaxis='y2')
    fig = go.Figure(data=[trace1, trace2, trace3, trace4])
    fig.update_layout(
        width=1000, 
        height=600,
        xaxis=dict(title='$时间 (s)$'),
        yaxis=dict(title='$垂荡位移 (m)$'),
        yaxis2=dict(title='$垂荡速度 (m/s)$', anchor='x', overlaying='y', side='right'),
        legend=dict(y=1.22, yanchor="top", x=1, xanchor="right"),
        title=title,
        template='plotly_white',
    )
    if svg_name is not None:
        fig.write_image(IMG_SVG / svg_name)
    fig.show()
    return None
def plot_plotly_rw(title, t=t, y=res3, svg_name=None):
    trace1 = go.Scatter(x=t, y=y[:, 4], name="$浮子纵摇角位移~\\theta_1$", yaxis='y1', line=dict(color='rgb(232,137,189)'))
    trace2 = go.Scatter(x=t, y=y[:, 5], name="$浮子纵摇角速度~ \omega_1$", yaxis='y2', line=dict(color='rgb(103,194,163)'))
    trace3 = go.Scatter(x=t, y=y[:, 6], name="$振子纵摇角位移~\\theta_2$", yaxis='y1', line=dict(color='rgb(252,140,99)'))
    trace4 = go.Scatter(x=t, y=y[:, 7], name="$振子纵摇角速度~ \omega_2$", yaxis='y2', line=dict(color='rgb(142,160,201)'))
    fig = go.Figure(data=[trace1, trace2, trace3, trace4])
    fig.update_layout(
        width=1000, 
        height=600,
        xaxis=dict(title='$时间 (s)$'),
        yaxis=dict(title='$纵摇角位移 (rad)$'),
        yaxis2=dict(title='$纵摇角速度 (rad/s)$', anchor='x', overlaying='y', side='right'),
        legend=dict(y=1.22, yanchor="top", x=1, xanchor="right"),
        title=title,
        template='plotly_white',
    )
    if svg_name is not None:
        fig.write_image(IMG_SVG / svg_name)
    fig.show()
    return None

plot_plotly_xv('浮子和振子的垂荡位移、垂荡速度', svg_name='问题3-浮子和振子的垂荡位移、垂荡速度.svg')
plot_plotly_rw('浮子和振子的纵摇角位移、纵摇角速度', svg_name='问题3-浮子和振子的纵摇角位移、纵摇角速度.svg')

# 相对位移 相对速度
def plot_plotly_diff_xvrw(t=t, y=res3, save=False, show=False):
    trace1 = go.Scatter(x=t, y=y[:, 0]-y[:, 2], name="$浮子和振子的相对垂荡位移 ~ X$", line=dict(color='#8ECFC9'))
    trace2 = go.Scatter(x=t, y=y[:, 1]-y[:, 3], name="$浮子和振子的相对垂荡速度 ~ X$", line=dict(color='#FA7F6F'))
    trace3 = go.Scatter(x=t, y=y[:, 4]-y[:, 6], name="$浮子和振子的相对纵摇角位移~\\theta$", line=dict(color='#82B0D2'))
    trace4 = go.Scatter(x=t, y=y[:, 5]-y[:, 7], name="$浮子和振子的相对纵摇角速度~ \omega$", line=dict(color='#BEB8DC'))
    fig1 = go.Figure(data=[trace1])
    fig2 = go.Figure(data=[trace2])
    fig3 = go.Figure(data=[trace3])
    fig4 = go.Figure(data=[trace4])
    fig1.update_layout(
        width=1000, height=600,
        xaxis=dict(title='$时间 (s)$'),
        yaxis=dict(title='相对垂荡位移 (m)'),
        title='$浮子和振子的相对垂荡位移$',
    )
    fig2.update_layout(
        width=1000, height=600,
        xaxis=dict(title='$时间 (s)$'),
        yaxis=dict(title='相对垂荡速度 (m/s)'),
        title='$浮子和振子的相对垂荡速度$',
    )
    fig3.update_layout(
        width=1000, height=600,
        xaxis=dict(title='$时间 (s)$'),
        yaxis=dict(title='相对纵摇角位移 (rad)'),
        title='$浮子和振子的相对纵摇角位移$',
    )
    fig4.update_layout(
        width=1000, height=600,
        xaxis=dict(title='$时间 (s)$'),
        yaxis=dict(title='相对纵摇角速度 (rad/s)'),
        title='$浮子和振子的相对纵摇角速度$',
#         legend=dict(y=1.22, yanchor="top", x=1, xanchor="right"),
#         template='plotly_white',
    )
    
    if save:
        fig1.write_image(IMG_SVG / svg_name)
        fig2.write_image(IMG_SVG / svg_name)
        fig3.write_image(IMG_SVG / svg_name)
        fig4.write_image(IMG_SVG / svg_name)
    if show:
        fig1.show()
        fig2.show()
        fig3.show()
        fig4.show()
    return None

plot_plotly_diff_xvrw(save=False, show=True)