检验与分析¶

In [170]:
# 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")
In [171]:
# 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
In [172]:
# 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'
In [173]:
# 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
In [189]:
# TODO 附件3参数

class parameters:
    """ 随问题参数变化的参数 """
    global m1, m2, me, m1_, m2_, j1, j2, je
    global omega, c, f, k, k1, k2, C, L, K, K1, K2
    
    m1, m2, me = 浮子质量, 振子质量, 垂荡附加质量
    m1_ = m1 + me
    m2_ = 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 = -静水恢复力矩系数

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
            入射波浪频率 = 1.4005  # s^{-1}
            垂荡附加质量 = 1335.535  # kg
            垂荡兴波阻尼系数 = 656.3616  # N·s/m
            垂荡激励力振幅 = 6250  # N
            波浪频率 = 入射波浪频率
            波浪周期 = 1 / 波浪频率
        elif question == 2:
            # 问题2:参数
            纵摇附加转动惯量 = 7131.29
            纵摇兴波阻尼系数 = 2992.724
            纵摇激励力矩振幅 = 2560
            入射波浪频率 = 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 / 波浪频率
        
        _ = parameters()
        
        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(1)
_ = trange(0, 200, 0.2)
In [177]:
# 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直线阻尼器的阻尼系数

检验吻合度¶

$F_{吻合度} = 1 - \cfrac{y1 - y2}{y1}$

问题1参数的不同模型¶

问题1模型和弱化的问题3模型

问题1模型¶

In [190]:
# TODO q1-model
def diff_equation_q1(ys, t, c=c, k=k, k1=k1, k2=k2):
    y1 = ys[2-1]
    y3 = ys[4-1]
    
    y4 = c * (ys[2-1] - ys[4-1]) + k * (ys[1-1] - ys[3-1])
    y2 = f * np.cos(omega * t) + k1 * ys[2-1] + k2 * ys[1-1] - y4
    
    y2 = y2 / m1_
    y4 = y4 / m2
    return [y1, y2, y3, y4]
def get_result1_df(result1, t):
    """
    :param result1: np.ndarray 问题1的结果(有两小问)
    :param t: np.ndarray 时间
    """
    columns = ['时间 (s)', '浮子位移 (m)', '浮子速度 (m/s)', '振子位移 (m)', '振子速度 (m/s)']
    shijian = pd.DataFrame(t, columns=columns[:1])
    result1 = pd.DataFrame(result1, columns=columns[1:])
    result1_df = pd.concat([shijian, result1], axis=1)
    return result1_df
In [191]:
_ = question1234(1)
_ = trange(0, 200, 0.2)

t = np.linspace(t_left, t_right, num=int(t_right / t_step) + 1)
y0 = [0, 0, 0, 0]
result1 = odeint(diff_equation_q1, y0, t)
result1_df = get_result1_df(result1, t)
result1_df.iplot(x='时间 (s)')

问题3模型(弱化)¶

In [192]:
# TODO q3-model
def diff_equation_q3(y, t, c=c, k=k, k1=k1, k2=k2, C=C, K=K, K1=K1, K2=K2):
    """ 弱化的问题3模型 """
    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, 0, 0, 0, 0]
def get_result3_df(t, result3):
    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
In [193]:
_ = question1234(1)
_ = trange(0, 200, 0.2)

t = np.linspace(t_left, t_right, num=int(t_right / t_step) + 1)
y0 = [0, 0, 0, 0, 0, 0, 0, 0]
result3 = odeint(diff_equation_q3, y0, t)
result3_df = get_result3_df(t, result3)
result3_df.iplot(x='时间 (s)')

计算吻合度¶

In [194]:
# TODO 查看二者之差
resault_diff13 = pd.DataFrame(result1_df.values - result3_df.iloc[:, [0, 1, 2, 5, 6]].values)
resault_diff13.iloc[:, 1:].iplot()
In [195]:
def cal_wenhedu(result1, result3):
    wenhe_avg = 1 - abs(result1[1:, :4] - result3[1:, :4]) / abs(result1[1:, :4])
    wenhe_avg = wenhe_avg[:, :]
    return wenhe_avg.mean(0)
cal_wenhedu(result1, result3)
Out[195]:
array([1., 1., 1., 1.])

问题2参数的不同模型¶

问题1模型¶

In [178]:
_ = question1234(2)
_ = trange(0, 200, 0.2)

t = np.linspace(t_left, t_right, num=int(t_right / t_step) + 1)
y0 = [0, 0, 0, 0]
result1 = odeint(diff_equation_q1, y0, t)
result1_df = get_result1_df(result1, t)
result1_df.iplot(x='时间 (s)')

问题3模型(弱化)¶

In [179]:
_ = question1234(2)
_ = trange(0, 200, 0.2)

t = np.linspace(t_left, t_right, num=int(t_right / t_step) + 1)
y0 = [0, 0, 0, 0, 0, 0, 0, 0]
result3 = odeint(diff_equation_q3, y0, t)
result3_df = get_result3_df(t, result3)
result3_df.iplot(x='时间 (s)')

计算吻合度¶

In [180]:
# TODO 查看二者之差
resault_diff13 = pd.DataFrame(result1_df.values - result3_df.iloc[:, [0, 1, 2, 5, 6]].values)
resault_diff13.iloc[:, 1:].iplot()
In [181]:
def cal_wenhedu(result1, result3):
    wenhe_avg = 1 - abs(result1[1:, :4] - result3[1:, :4]) / abs(result1[1:, :4])
    wenhe_avg = wenhe_avg[:, :]
    return wenhe_avg.mean(0)
cal_wenhedu(result1, result3)
Out[181]:
array([1., 1., 1., 1.])

灵敏性分析¶

使用问题3参数进行灵敏度分析!!!

相差¶

In [152]:
# 绘图:位移 速度
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
In [159]:
c = 2000  # 10000
k = 1600  # 80000

C = 3600  # 1000
K = 5000  # 250000
me = 1 * 垂荡附加质量

_ = question1234(1)
_ = trange(0, 100, 0.2)

def ode_func(y, t, k=k, K=K):
    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)
plot_plotly_xv('浮子和振子的垂荡位移', t=t, y=res3, svg_name='灵敏性分析-浮子和振子的垂荡位移.svg')
plot_plotly_rw('浮子和振子的纵摇角位移', t=t, y=res3, svg_name='灵敏性分析-浮子和振子的纵摇角位移.svg')

不同的 c k C K¶

In [160]:
InteractiveShell.ast_node_interactivity = 'last'

def ode_func(y, t, k, K, c=c, C=C):
    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]
In [161]:
# TODO 绘图
InteractiveShell.ast_node_interactivity = 'last'

data_X = []
data_V = []
data_theta = []
data_omega = []

for c, C, k, K in zip([10000 * 0.5, 10000, 10000 * 2], 
                      [1000 * 0.5, 1000, 1000 * 2], 
                      [80000 * 0.5, 80000, 80000 * 2], 
                      [250000 * 0.5, 250000, 250000 * 2]):
    t_left, t_right = 0, 100
    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, args=(k, K))
    
    fuzi_columnes = [f'k={k}-浮子垂荡位移', f'k={k}-浮子垂荡速度', f'k={k}-浮子纵摇角位移', f'k={k}-浮子纵摇角速度']
    zhenzi_columnes = [f'k={k}-振子垂荡位移', f'k={k}-振子垂荡速度', f'k={k}-振子纵摇角位移', f'k={k}-振子纵摇角速度']
    
    data_X.append(go.Scatter(x=t, y=res3[:, 0], name=fuzi_columnes[0]))
    data_V.append(go.Scatter(x=t, y=res3[:, 1], name=fuzi_columnes[1]))
    data_theta.append(go.Scatter(x=t, y=res3[:, 4], name=fuzi_columnes[2]))
    data_omega.append(go.Scatter(x=t, y=res3[:, 5], name=fuzi_columnes[3]))
    
    data_X.append(go.Scatter(x=t, y=res3[:, 2], name=zhenzi_columnes[0]))
    data_V.append(go.Scatter(x=t, y=res3[:, 3], name=zhenzi_columnes[1]))
    data_theta.append(go.Scatter(x=t, y=res3[:, 6], name=zhenzi_columnes[2]))
    data_omega.append(go.Scatter(x=t, y=res3[:, 7], name=zhenzi_columnes[3]))

fig_data_X = go.Figure(data=data_X)
fig_data_X.update_layout(
        width=1000, 
        height=600,
        xaxis=dict(title='$时间 (s)$'),
        yaxis=dict(title='$垂荡位移 (m)$'),
)
fig_data_V = go.Figure(data=data_V)
fig_data_V.update_layout(
        width=1000, 
        height=600,
        xaxis=dict(title='$时间 (s)$'),
        yaxis=dict(title='$垂荡速度 (m)$'),
)
fig_data_theta = go.Figure(data=data_theta)
fig_data_theta.update_layout(
        width=1000, 
        height=600,
        xaxis=dict(title='$时间 (s)$'),
        yaxis=dict(title='$纵摇角位移 (rad)$'),
)
fig_data_omega = go.Figure(data=data_omega)
fig_data_omega.update_layout(
        width=1000, 
        height=600,
        xaxis=dict(title='$时间 (s)$'),
        yaxis=dict(title='$纵摇角速度 (rad)$'),
#         legend=dict(y=1.22, yanchor="top", x=1, xanchor="right"),
#         title=title,
#         template='plotly_white'
)

save = False
if save:
    fig_data_X.write_image(IMG_SVG / "灵敏性分析-垂荡位移.svg")
    fig_data_V.write_image(IMG_SVG / "灵敏性分析-垂荡速度.svg")
    fig_data_theta.write_image(IMG_SVG / "灵敏性分析-纵摇角位移.svg")
    fig_data_omega.write_image(IMG_SVG / "灵敏性分析-纵摇角速度.svg")

fig_data_X.show()
fig_data_V.show()
fig_data_theta.show()
fig_data_omega.show()
In [ ]:
 
In [ ]:
 
In [ ]: