2019数学建模国赛A题 Python代码

import numpy as np import pandas as pd import matplotlib.pyplot as plt import math data = '附件3-弹性模量与压力.xlsx' data = pd.read_excel('附件3-弹性模量与压力.xlsx') # 导入附件3数据 plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签 plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号 x_pandas_list = data[u'燃油密度(mm3/ms)'] # 附件3 燃油密度列数据 y_pandas_list = data[u'压力(MPa)'] # 附件3 压力列数据 coefficient1_numpy_list = np.polyfit(x_pandas_list, y_pandas_list, 2) # [ 10784.65349198 -15693.86029255 5648.09019803]y=压强 x=密度拟合系数 def ρ_P(x): # 密度转压力函数 return coefficient1_numpy_list[0] * x ** 2 + coefficient1_numpy_list[1] * x + coefficient1_numpy_list[2] coefficient2_numpy_list = np.polyfit(y_pandas_list, x_pandas_list, 2) # [-6.59843062e-07 5.23409641e-04 8.04251284e-01]y=密度 x=压强拟合系数 def P_ρ(x): # 压力转密度函数 return coefficient2_numpy_list[0] * x ** 2 + coefficient2_numpy_list[1] * x + coefficient2_numpy_list[2] plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签 plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号 data = '附件1-凸轮边缘曲线.xlsx' data = pd.read_excel(data) x_pandas_list1 = data[u'极角（rad）'] y_pandas_list1 = data[u'极径（mm）'] cof = np.polyfit(x_pandas_list1, y_pandas_list1, 6) # 凸轮 x-极角 y-极径拟合系数 def y_(x): return cof[0] * x ** 6 + cof[1] * x ** 5 + cof[2] * x ** 4 + cof[3] * x ** 3 + cof[4] * x ** 2 + cof[5] * x + cof[6] ω = 0.06 # 弧度制下的角速度rad/ms T = 2 * math.pi / ω # 凸轮运动周期时间 单位ms t = 0 # 记录时刻 t_ = 0.1 # 每隔0.01ms时间段刷新一次状态 α0 = 0 # 初始时刻凸轮旋转角度0 α = α0 # 记录凸轮旋转角度时刻 hMax = 0 # 记录每一时刻高度的最大值 P_camj = 0.5 # 初始时刻压力 ρ_camj = P_ρ(0.5) # 初始时刻密度 d_camj = 5 # 柱塞腔直径5mm r_camj = d_camj / 2 # 柱塞腔半径2.5mm s_camj = math.pi * r_camj ** 2 v_part = 20 # 向上残余容积为20mm3 h_part = v_part / s_camj # 向上残余高度 m0_camj = 92.4217168684373 # 初始时刻柱塞腔油气质量 m_camj = m0_camj # 每一时刻柱塞腔油气质量 C = 0.85 # 流量系数 d_a = 1.4 # 小孔处的直径1.4mm r_a = d_a / 2 # 小孔处的半径1.4mm A = S_a = math.pi * r_a ** 2 # 小孔处的面积1.5393804002589984mm2 l_primary = 500 # 内腔长度500mm d_primary = 10 # 内腔直径10mm r_primary = d_primary / 2 # 内腔半径5mm S_primary = math.pi * r_primary ** 2 # 内腔横截面积78.53981633974483mm2 V_primary = S_primary * l_primary # 内腔体积39269.90816987242mm3 ρ0 = 0.85 # 初始100MPa压力下的密度 ρ160 = P_ρ(160) # 160MPa压力下的密度 0.8711048440368855mg/mm3 m0 = ρ0 * V_primary # 初始高压油管油气33379.42194439156mg m = m0 # 记录每一时刻高压油管内的质量 初始为m0 P = 100 # 记录每一时刻高压油管内的压强 ρ = 0.85 # 录每一时刻高压油管内的密度 P_list = [] # 压强时刻表 def cam_h(t): # x为时刻 返回值为活塞距离柱塞腔底部的高度 global α, hMax # 全局变量可做修改 hmax = 0 # 记录这一时刻下活塞向上运动的最大值 for i in np.arange(0, 2 * math.pi, 0.1): # 弧度制遍历 h_point = y_(i) * math.cos(ω * t + i) # 每一时间点下不同弧度坐标位置 h_point y坐标 if h_point > hmax: # 把每一个极角遍历出的最大值传递给hmax hmax = h_point hMax = 7.247588972988529 - hmax + h_part return hMax def P_cam(t): # 柱塞腔返回t时刻下的压力 return ρ_P(ρ_cam(t)) def ρ_cam(t): # 柱塞腔返回t时刻下燃油的密度 global m_camj, ρ_camj, P_camj if t % T <= T / 2: m_camj = m0_camj P_camj = 0.5 else: if P_camj <= 100: ρ_camj = m_camj / (s_camj * cam_h(t)) # 更新密度 P_camj = ρ_P(ρ_camj) # 更新压强 else: det_m = C * A * math.sqrt(2 * (P_camj - 100) / ρ_camj) * ρ_camj * t_ m_camj = m_camj - det_m ρ_camj = m_camj / (s_camj * cam_h(t)) # 更新密度 P_camj = ρ_P(ρ_camj) # 更新压强 P_list = [] # 压强时刻表 # 进油函数 参数t为时刻 def enter(t): global m, P, ρ # 全局变量，可做修改 print(t) if P_camj > 100.1769: Q = C * A * math.sqrt(2 * (P_camj - 100.1769) / P_ρ(P_camj)) # 单位时间进油量 det_m = Q * P_ρ(P_camj) * t_ # 0.01为步长 质量改变量 m = m + det_m # 更新质量 rou = m / V_primary # 更新密度 P = ρ_P(rou) # 更新压强 else: pass # 针阀拟合 data = '附件2-针阀运动曲线.xlsx' data = pd.read_excel('附件2-针阀运动曲线.xlsx') # 导入附件3数据 plt.rcParams['font.sans-serif'] = ['SimHei'] # 用来正常显示中文标签 plt.rcParams['axes.unicode_minus'] = False # 用来正常显示负号 x_pandas_list = data[u'时间(ms)'] # 附件2 时间数据 y_pandas_list = data[u'距离(mm)'] # 附件2 距离数据 y1_pandas_list = y_pandas_list[0:46] x1_pandas_list = x_pandas_list[0:46] cof1 = np.polyfit(x1_pandas_list, y1_pandas_list, 6) # 第一段分段函数系数[15.39735971 -1.93086567 0.026043 ] y2_pandas_list = y_pandas_list[200:246] x2_pandas_list = x_pandas_list[200:246] cof2 = np.polyfit(x2_pandas_list, y2_pandas_list, 6) # 第三段分段函数系数[ 15.44372154 -73.71211028 87.92142494] A_s = 0 def H(t): if 0 < t % 200 <= 0.45: return cof1[0] * (t % 200) ** 6 + cof1[1] * (t % 200) ** 5 + cof1[2] * (t % 200) ** 4 + cof1[3] \ * (t % 200) ** 3 + cof1[4] * (t % 200) ** 2 + cof1[5] * (t % 200) + cof1[6] elif 0.45 < t % 200 < 2: return 2 elif 2 <= t % 200 < 2.45: return cof2[0] * (t % 200) ** 6 + cof2[1] * (t % 200) ** 5 + cof2[2] * (t % 200) ** 4 + cof2[3] \ * (t % 200) ** 3 + cof2[4] * (t % 200) ** 2 + cof2[5] * (t % 200) + cof2[6] else: return 0 def Aj_1(t): # 阻塞模型下的流出面积 global A_s if t % 200 < 0.33: A_s = math.pi * (1.25 + H(t) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi elif 0.33 <= t % 200 < math.sqrt(6) - 0.33: A_s = 0.7 ** 2 * math.pi elif math.sqrt(6) - 0.33 <= t % 200 < math.sqrt(6): A_s = math.pi * (1.25 + (H(math.sqrt(6) - t)) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi else: A_s = 0 def out_1(t): global m, P, ρ_camj # 全局变量，可做修改 Q = C * A_s * math.sqrt(2 * P / ρ_camj) # 单位时间进油量 det_m = Q * ρ_camj * t_ # 0.01为步长 质量改变量 m = m - det_m # 更新质量 ρ_camj = m / V_primary # 更新密度 P = ρ_P(ρ_camj) # 更新压强 def Aj_2(t): # 阻塞模型下的流出面积 global A_s if t % 200 < 0.33 + 50: A_s = math.pi * (1.25 + H(t) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi elif 0.33 + 50 <= t % 200 < math.sqrt(6) - 0.33 + 50: A_s = 0.7 ** 2 * math.pi elif math.sqrt(6) - 0.33 + 50 <= t % 200 < math.sqrt(6) + 50: A_s = math.pi * (1.25 + (H(math.sqrt(6) - t)) * math.tan(math.pi / 20)) ** 2 - 1.25 ** 2 * math.pi else: A_s = 0 def out_2(t): global m, P, ρ_camj # 全局变量，可做修改 Q = C * A_s * math.sqrt(2 * P / ρ_camj) # 单位时间进油量 det_m = Q * ρ_camj * t_ # 0.01为步长 质量改变量 m = m - det_m # 更新质量 ρ_camj = m / V_primary # 更新密�