数模国赛2020A复现
文章目录
1. 前言
优秀论文[1]和[2]都建立了偏微分方程模型，以原件的厚度 x x x作为一个变量，以时间 t t t作为另一个变量，最后用有限元法（有限差分法？）求解偏微分方程。其中[1]的求解效果较好，代码较全，[2]稍微差一些，又改用了其他模型。但是无论哪篇论文，运算量都很大，而且参数求法都是先粗略搜索，再精细搜索，程序运行时间是个谜，少则一两个小时，多则几天也有可能， 比赛时心里没底的话，很悬，能准确求出结果也需要运气。
总的来说，研究这个问题有助于掌握偏微分建模法以及有限元求解的方法。
[1]中的参数求出来都特别小，数量级是 1 0 − 10 10^{-10} 10−10，我猜测是因为代码中把长度单位统一为 m m m的缘故，这么小的数量级不好计算，可见单位的选择也要有所考量。
2. 模型建立
2.1 偏微分方程模型
如上图所示，中间的长方形代表元件，横向平移表示在时间尺度上的变化，元件上下表面在回焊炉内被加热，以元件厚度中心为
x
x
x轴0点，假设
T
(
x
,
t
)
T(x,t)
T(x,t)表示
t
t
t时刻
x
x
x处的温度，建立以下偏微分方程模型：
{
∂
T
(
x
,
t
)
∂
t
−
α
∂
2
T
(
x
,
t
)
∂
x
2
=
0
∂
T
(
x
,
t
)
∂
x
∣
x
=
−
d
2
=
σ
(
T
(
x
,
t
)
−
T
a
i
r
(
t
)
)
∂
T
(
x
,
t
)
∂
x
∣
x
=
d
2
=
−
σ
(
T
(
x
,
t
)
−
T
a
i
r
(
t
)
)
\begin{cases} \dfrac{\partial T(x,t)}{\partial t} - \alpha\dfrac{\partial^2T(x,t)}{\partial x^2} = 0 \\[3ex] \left.\dfrac{\partial T(x,t)}{\partial x}\right|_{x=-\frac{d}{2}}=\sigma (T(x,t) - T_{air}(t)) \\[3ex] \left.\dfrac{\partial T(x,t)}{\partial x}\right|_{x=\frac{d}{2}}=-\sigma (T(x,t) - T_{air}(t)) \end{cases}
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧∂t∂T(x,t)−α∂x2∂2T(x,t)=0∂x∂T(x,t)∣∣∣∣x=−2d=σ(T(x,t)−Tair(t))∂x∂T(x,t)∣∣∣∣x=2d=−σ(T(x,t)−Tair(t))
第一个式子是温度传播规律，假设元件厚度为
d
d
d，则第2个式子和第3个式子是边界温度传播的条件。
σ
\sigma
σ前面符号相反即可，负号加到第2个式子前面也可。
T a i r ( t ) T_{air}(t) Tair(t)表示回焊炉内温度的变化情况，只和时间有关。
2.2离散化处理（有限差分法）
取时间步长
Δ
t
\Delta t
Δt，空间步长
Δ
x
\Delta x
Δx，将连续的平面区域离散化，建立二维平面网格，网格点坐标为[1]：
{
x
i
=
−
d
2
+
i
Δ
x
,
(
i
=
0
,
1
,
2
,
⋯
,
⌊
d
Δ
x
⌋
)
t
j
=
j
Δ
t
,
(
j
=
0
,
1
,
2
,
⋯
,
⌊
t
m
Δ
t
⌋
)
\begin{cases} x_i = -\frac{d}{2} + i\Delta x, & (i = 0,1,2,\cdots,\lfloor \frac{d}{\Delta x} \rfloor) \\[3ex] t_j = j \Delta t, & (j=0,1,2,\cdots,\lfloor \frac{t_m}{\Delta t} \rfloor) \end{cases}
⎩⎪⎨⎪⎧xi=−2d+iΔx,tj=jΔt,(i=0,1,2,⋯,⌊Δxd⌋)(j=0,1,2,⋯,⌊Δttm⌋)
其中
t
m
t_m
tm表示总时间，则
T
[
i
]
[
j
]
=
T
(
x
i
,
t
j
)
T[i][j]=T(x_i,t_j)
T[i][j]=T(xi,tj)。
将上述第1个公式离散化：
T
[
i
]
[
j
]
−
T
[
i
]
[
j
−
1
]
Δ
t
=
α
T
[
i
+
1
]
[
j
]
−
2
T
[
i
]
[
j
]
+
T
[
i
−
1
]
[
j
]
Δ
x
2
\frac{T[i][j] - T_[i][j-1]}{\Delta t}=\alpha \frac{T[i+1][j]-2T[i][j]+T[i-1][j]}{\Delta x^2}
ΔtT[i][j]−T[i][j−1]=αΔx2T[i+1][j]−2T[i][j]+T[i−1][j]
令
A
=
α
Δ
t
Δ
x
2
A = \alpha\frac{\Delta t}{\Delta x^2}
A=αΔx2Δt
得到：
T
[
i
]
[
j
−
1
]
=
−
A
T
[
i
−
1
]
[
j
]
+
(
2
A
+
1
)
T
[
i
]
[
j
]
−
A
T
[
i
+
1
]
[
j
]
(1)
T[i][j-1]=-AT[i-1][j]+(2A+1)T[i][j]-AT[i+1][j] \tag{1}
T[i][j−1]=−AT[i−1][j]+(2A+1)T[i][j]−AT[i+1][j](1)
将第2个式子离散化，得到：
T
[
1
]
[
j
]
−
T
[
0
]
[
j
]
Δ
x
=
σ
(
T
[
0
]
[
j
]
−
T
a
i
r
(
j
Δ
t
)
)
\frac{T[1][j] - T[0][j]}{\Delta x} = \sigma(T[0][j] - T_{air}(j\Delta t))
ΔxT[1][j]−T[0][j]=σ(T[0][j]−Tair(jΔt))
化简得到：
σ
Δ
x
T
a
i
r
(
j
Δ
t
)
=
(
σ
Δ
x
+
1
)
T
[
0
]
[
j
]
−
T
[
1
]
[
j
]
(2)
\sigma\Delta xT_{air}(j\Delta t) = (\sigma\Delta x + 1)T[0][j] - T[1][j] \tag{2}
σΔxTair(jΔt)=(σΔx+1)T[0][j]−T[1][j](2)
将第3个式子离散化，得到：
T
[
m
]
[
j
]
−
T
[
m
−
1
]
[
j
]
Δ
x
=
−
σ
(
T
[
m
]
[
j
]
−
T
a
i
r
(
j
Δ
t
)
)
\frac{T[m][j] - T[m-1][j]}{\Delta x} = -\sigma(T[m][j] - T_{air}(j\Delta t))
ΔxT[m][j]−T[m−1][j]=−σ(T[m][j]−Tair(jΔt))
其中
m
=
⌊
d
Δ
x
⌋
m=\lfloor \frac{d}{\Delta x} \rfloor
m=⌊Δxd⌋。
化简得：
σ
Δ
x
T
a
i
r
(
j
Δ
t
)
=
(
σ
Δ
x
+
1
)
T
[
m
]
[
j
]
−
T
[
m
−
1
]
[
j
]
(3)
\sigma\Delta xT_{air}(j\Delta t) = (\sigma\Delta x + 1)T[m][j] - T[m-1][j] \tag{3}
σΔxTair(jΔt)=(σΔx+1)T[m][j]−T[m−1][j](3)
综合(1)(2)(3)，将方程写为矩阵形式：
记为：
Ω
T
j
=
B
j
−
1
\Omega T_j=B_{j-1}
ΩTj=Bj−1
初始得到
T
0
T_0
T0，构建
Ω
,
B
j
−
1
\Omega,B_{j-1}
Ω,Bj−1，由此计算出
T
j
T_j
Tj；再计算
Ω
,
B
j
\Omega,B_j
Ω,Bj，从而计算
T
j
+
1
T_{j+1}
Tj+1，迭代这个过程。
2.3 参数讨论
模型中有两个参数需要确定，分别是 σ \sigma σ和 α \alpha α。事实上，一共有10个小温区分属为5个大温区，第 i i i个大温区选择参数 σ i , α i \sigma_i,\alpha_i σi,αi，即每个大温区的参数不一样，这样通过模型去拟合所给的附件中的数据，效果更好，但是参数量变多，求解更加困难。
2.4 回焊炉内的温度分布
温区内保持均温，温区之间呈现线性关系，图像如下，其中纵坐标为温度 T T T，横坐标为距离 s s s：
画图代码：
import numpy as np
import matplotlib.pyplot as plt
#绘制T_air随距离变化的图像，分段函数
T1 = 175 #小温区1~5
T2 = 195 #6
T3 = 235 #7
T4 = 255 #8~9
def Tair_s(s): #当前距离和空气温度的函数
if s < 25:
T = (T1 - 25) / (25 - 0) * (s - 0) + 25
elif s < 197.5:
T = T1
elif s < 202.5:
T = (T2 - T1) / (202.5 - 197.5) * (s - 197.5) + T1
elif s < 233:
T = T2
elif s < 238:
T = (T3 - T2) / (238 - 233) * (s - 233) + T2
elif s < 268.5:
T = T3
elif s < 273.5:
T = (T4 - T3) / (273.5 - 268.5) * (s - 268.5) + T3
elif s < 339.5:
T = T4
elif s < 344.5:
T = (25 - T4) / (344.5 - 339.5) * (s - 339.5) + T4
else:
T = 25
return T
x = np.linspace(0, 450, 10000)
y = np.array([Tair_s(s) for s in x]) #根据x计算各点函数值
plt.figure()
plt.plot(x, y)
#plt.ylim(-0.2, 1.2) #限制y的范围
#plt.xlim(0, 2) #限制x的范围
plt.show()
2.5 模型参数拟合
代码参考[1]，先大概估计参数的范围（主要是根据所建立的模型的物理量关系，这里 α \alpha α的值比较小，数量级大致为1e-11），使用python求解。获得整体温度分布图如下图，可以看到是两边对称的：
其中心温度，即x=0的温度随时间t的变化如下图所示：
和附件中给的数据对比如下：
模型值和所给的附件中的数据差别还是比较大的，可见随意选择的参数并不能很好的符合。
代码如下：
#有限元差分求解代码
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei']
##固定参数
delta_x = 1e-6 #单位m
delta_t = 0.5 #单位s
thickness = 0.15e-3 #元件厚度，m
t_m = 373 #附件中的总时间
v = 70 / 60 #传送带前进速度
T1 = 175 #小温区1~5
T2 = 195 #6
T3 = 235 #7
T4 = 255 #8~9
size_x = round(thickness / delta_x + 1)
size_t = round(t_m / delta_t + 1)
##待求参数，先根据[2]求参的结果，画图检验一下
sigma = 18657485
alpha = [4.437e-11, 5.621e-11, 7.449e-11, 4.997e-11, 2.401e-11] #第1个温区的alpha参数，以此类推
t_interval = [1, 348, 409, 470, 592, 747] #5个温区的截止时间，以此为第1个温区的……
##待求值
T = np.zeros((size_x, size_t))
T[:,0] = 25
def Tair_s(s): #当前距离和空气温度的函数
if s < 25:
T = (T1 - 25) / (25 - 0) * (s - 0) + 25
elif s < 197.5:
T = T1
elif s < 202.5:
T = (T2 - T1) / (202.5 - 197.5) * (s - 197.5) + T1
elif s < 233:
T = T2
elif s < 238:
T = (T3 - T2) / (238 - 233) * (s - 233) + T2
elif s < 268.5:
T = T3
elif s < 273.5:
T = (T4 - T3) / (273.5 - 268.5) * (s - 268.5) + T3
elif s < 339.5:
T = T4
elif s < 344.5:
T = (25 - T4) / (344.5 - 339.5) * (s - 339.5) + T4
else:
T = 25
return T
for i in range(0, len(alpha)):
A = alpha[i] * delta_t / 2 / (delta_x * delta_x)
Omega = np.eye(size_x, size_x) #构建Omega矩阵
Omega[0, 0] = 1 + sigma * delta_x
Omega[0, 1] = -1
Omega[size_x-1, size_x-2] = -1
Omega[size_x-1, size_x-1] = 1 + sigma * delta_x
for k in range(1, size_x - 1):
Omega[k, k-1] = -A
Omega[k, k] = 2 * A + 1
Omega[k, k+1] = -A
B = np.zeros((size_x, 1)) #构建B矩阵
for j in range(t_interval[i], t_interval[i+1]):
s = j * delta_t * v #当前前进距离
T_air = Tair_s(s)
B[0] = sigma * delta_x * T_air
B[size_x - 1] = sigma * delta_x * T_air
for j1 in range(1, size_x - 1):
B[j1] = T[j1, j-1]
T[:, [j]] = np.linalg.inv(Omega).dot(B) #注意矩阵乘法使用dot,以及[j]，后者影响维度
##画三维图
x = np.linspace(-thickness/2, thickness/2, size_x) #厚度
y = np.linspace(0, t_m, size_t) #时间
X, Y = np.meshgrid(x, y)
fig1 = plt.figure()
ax = plt.axes(projection='3d') #ax是一个变量，像是str， ndarray等
ax.plot_surface(X, Y, T.T, cmap ='jet', edgecolor ='green') #T需要转置一下，否则维度对不上
ax.set_xlabel('x', fontsize = 20)
ax.set_ylabel('t', fontsize = 20)
ax.set_zlabel('T', fontsize = 20);
###画二维图
fig2 = plt.figure()
file_path = '附件.xls' #读取附件
data_1 = pd.read_excel(file_path, sheet_name = "Sheet1")
data_2 = data_1.values #把dataframe(data_1)转化为ndarray(data_2)
fj_x = data_2[:, 0]
fj_y = data_2[:, 1]
plot1 = plt.plot(y[38:], T[75, 38:], color = 'green', label='模型值')
plot2 = plt.plot(fj_x, fj_y, color = 'red', label = '附件值')
plt.title('工件中心温度随时间的变化关系', fontsize = 20)
plt.legend(loc = 1, fontsize = 20)
plt.xlabel('t', fontsize = 20)
plt.ylabel('T', fontsize = 20)
plt.show()
以上代码中的参数取值分别为：
σ
=
8657485
\sigma =8657485
σ=8657485
α = [ 4.437 e − 11 , 5.621 e − 11 , 7.449 e − 11 , 4.997 e − 11 , 2.401 e − 11 ] {\alpha}=[4.437e-11, 5.621e-11, 7.449e-11, 4.997e-11, 2.401e-11] α=[4.437e−11,5.621e−11,7.449e−11,4.997e−11,2.401e−11]
更改 σ \sigma σ的值发现影响不大对于温度的曲线，那么关键参数只有 α \alpha α了。
考虑到模型的复杂度，使得向“梯度下降”一类的求参方法比较困难（不知这类方法是否可行），只能通过搜索来了，所谓搜索，无非就是5重for循环暴力枚举，复杂度不用说了，能不能在规定时间内跑出来看脸。（-_-!）
不过还是有些技巧的，我们可以先通过粗略搜索，框定一个大概的范围，然后再这个范围内不断缩小搜索精度。这是一门技术活，需要非凡的意志和耐心。
我们把以下函数作为损失函数：
L
=
1
2
∑
i
=
1
n
(
T
^
[
i
]
−
T
[
i
]
)
2
L=\frac{1}{2}\sum_{i=1}^n(\hat{T}[i]-T[i])^2
L=21i=1∑n(T^[i]−T[i])2
其中
T
^
\hat{T}
T^为附件给出的值，
T
T
T为模型预测值，
n
n
n为附件中的数据总数，注意
i
=
1
i=1
i=1对应的时间为19s。
我们先以
1
×
1
0
−
11
1\times10^{-11}
1×10−11的精度开始搜索，建立以下优化模型：
arg
min
α
L
=
1
2
∑
i
=
1
n
(
T
^
[
i
]
−
T
[
i
]
)
2
\arg \mathop{\min}_{\alpha}L=\frac{1}{2}\sum_{i=1}^n(\hat{T}[i]-T[i])^2
argminαL=21i=1∑n(T^[i]−T[i])2
s . t . α i = j e − 11 , ( j = 1 , 2 , ⋯ , 9 i = 1 , 2 , 3 , 4 , 5 ) s.t. \quad \alpha_i=je-11,(j=1,2,\cdots,9\quad i=1,2,3,4,5) s.t.αi=je−11,(j=1,2,⋯,9i=1,2,3,4,5)
结果发现5重for循环总共计算1e5次，时间复杂度太高。转而想到这5个参数可以分别拟合，比如先拟合 α 1 \alpha_1 α1，得到其最优值后再拟合 α 2 \alpha_2 α2，以此类推如下图：
代码如下：
#5重for循环求参
#有限元差分求解代码
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei']
##固定参数
delta_x = 1e-6 #单位m
delta_t = 0.5 #单位s
thickness = 0.15e-3 #元件厚度，m
t_m = 373 #附件中的总时间
v = 70 / 60 #传送带前进速度
T1 = 175 #小温区1~5
T2 = 195 #6
T3 = 235 #7
T4 = 255 #8~9
size_x = round(thickness / delta_x + 1)
size_t = round(t_m / delta_t + 1)
##待求参数，先根据[2]求参的结果，画图检验一下
sigma = 8657485
#alpha = [4.437e-11, 5.621e-11, 7.449e-11, 4.997e-11, 2.401e-11] #第1个温区的alpha参数，以此类推
t_interval = [1, 348, 409, 470, 592, 747] #5个温区的截止时间，以此为第1个温区的……
##待求值
T = np.zeros((size_x, size_t))
T[:,0] = 25
def Tair_s(s): #当前距离和空气温度的函数
if s < 25:
T = (T1 - 25) / (25 - 0) * (s - 0) + 25
elif s < 197.5:
T = T1
elif s < 202.5:
T = (T2 - T1) / (202.5 - 197.5) * (s - 197.5) + T1
elif s < 233:
T = T2
elif s < 238:
T = (T3 - T2) / (238 - 233) * (s - 233) + T2
elif s < 268.5:
T = T3
elif s < 273.5:
T = (T4 - T3) / (273.5 - 268.5) * (s - 268.5) + T3
elif s < 339.5:
T = T4
elif s < 344.5:
T = (25 - T4) / (344.5 - 339.5) * (s - 339.5) + T4
else:
T = 25
return T
def compu_loss(): #计算损失函数的值
delta_T = T[75, 38:347] - fj_y[0:309]
L = 1/2 * np.sum(delta_T * delta_T)
return L
def compu_T(): #计算当前的温度T
for i in range(0, len(alpha)):
A = alpha[i] * delta_t / 2 / (delta_x * delta_x)
Omega = np.eye(size_x, size_x) # 构建Omega矩阵
Omega[0, 0] = 1 + sigma * delta_x
Omega[0, 1] = -1
Omega[size_x - 1, size_x - 2] = -1
Omega[size_x - 1, size_x - 1] = 1 + sigma * delta_x
for k in range(1, size_x - 1):
Omega[k, k - 1] = -A
Omega[k, k] = 2 * A + 1
Omega[k, k + 1] = -A
B = np.zeros((size_x, 1)) # 构建B矩阵
for j in range(t_interval[i], t_interval[i + 1]):
s = j * delta_t * v # 当前前进距离
T_air = Tair_s(s)
B[0] = sigma * delta_x * T_air
B[size_x - 1] = sigma * delta_x * T_air
for j1 in range(1, size_x - 1):
B[j1] = T[j1, j - 1]
T[:, [j]] = np.linalg.inv(Omega).dot(B) # 注意矩阵乘法使用dot,以及[j]，后者影响维度
file_path = '附件.xls' #读取附件
data_1 = pd.read_excel(file_path, sheet_name = "Sheet1")
data_2 = data_1.values #把dataframe(data_1)转化为ndarray(data_2)
fj_x = data_2[:, 0]
fj_y = data_2[:, 1]
accuracy = 1e-11
alpha_op = [] #alpha的最优值
L_min = 1e99 #损失函数的最小值
alpha = [0,0,0,0,0]
flag = 1
for alpha1 in range(1, 16):
alpha1 = alpha1 * accuracy
alpha[0] = alpha1
for alpha2 in range(1, 2):
alpha2 = alpha2 * accuracy
alpha[1] = alpha2
for alpha3 in range(1, 2):
alpha3 = alpha3 * accuracy
alpha[2] = alpha3
for alpha4 in range(1, 2):
alpha4 = alpha4 * accuracy
alpha[3] = alpha4
for alpha5 in range(1, 2):
alpha5 = alpha5 * accuracy
alpha[4] = alpha5
print(flag)
flag = flag+1
compu_T()
L = compu_loss()
if L < L_min:
L_min = L
#alpha_op = alpha #不能直接等于，否则alpha_op和alpha共用一块内存,一个变化，另一个也变
alpha_op = alpha.copy()
print(alpha_op)
alpha = alpha_op
compu_T()
# ##画三维图
# x = np.linspace(-thickness/2, thickness/2, size_x) #厚度
y = np.linspace(0, t_m, size_t) #时间
#
# X, Y = np.meshgrid(x, y)
#
# fig1 = plt.figure()
# ax = plt.axes(projection='3d') #ax是一个变量，像是str， ndarray等
# ax.plot_surface(X, Y, T.T, cmap ='jet', edgecolor ='green') #T需要转置一下，否则维度对不上
# ax.set_xlabel('x', fontsize = 20)
# ax.set_ylabel('t', fontsize = 20)
# ax.set_zlabel('T', fontsize = 20);
###画二维图
fig2 = plt.figure()
plot1 = plt.plot(y[38:], T[75, 38:], color = 'green', label='模型值')
plot2 = plt.plot(fj_x, fj_y, color = 'red', label = '附件值')
plt.title('工件中心温度随时间的变化关系', fontsize = 20)
plt.legend(loc = 1, fontsize = 20)
plt.xlabel('t', fontsize = 20)
plt.ylabel('T', fontsize = 20)
plt.show()
得到最优值后
α
1
\alpha_1
α1取最优值，遍历
α
2
\alpha_2
α2，固定其他变量，以此类推最终得到：
α
=
[
9
,
11
,
15
,
10
,
5
]
×
1
0
−
11
\alpha=[9,11,15,10,5]\times10^{-11}
α=[9,11,15,10,5]×10−11
拟合效果如下，可见效果非常不错：
全阶段代码：
#5重for循环求参
#有限元差分求解代码
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei']
##固定参数
delta_x = 1e-6 #单位m
delta_t = 0.5 #单位s
thickness = 0.15e-3 #元件厚度，m
t_m = 373 #附件中的总时间
v = 70 / 60 #传送带前进速度
T1 = 175 #小温区1~5
T2 = 195 #6
T3 = 235 #7
T4 = 255 #8~9
size_x = round(thickness / delta_x + 1)
size_t = round(t_m / delta_t + 1)
##待求参数，先根据[2]求参的结果，画图检验一下
sigma = 8657485
#alpha = [4.437e-11, 5.621e-11, 7.449e-11, 4.997e-11, 2.401e-11] #第1个温区的alpha参数，以此类推
t_interval = [1, 348, 409, 470, 592, 747] #5个温区的截止时间，以此为第1个温区的……
##待求值
T = np.zeros((size_x, size_t))
T[:,0] = 25
def Tair_s(s): #当前距离和空气温度的函数
if s < 25:
T = (T1 - 25) / (25 - 0) * (s - 0) + 25
elif s < 197.5:
T = T1
elif s < 202.5:
T = (T2 - T1) / (202.5 - 197.5) * (s - 197.5) + T1
elif s < 233:
T = T2
elif s < 238:
T = (T3 - T2) / (238 - 233) * (s - 233) + T2
elif s < 268.5:
T = T3
elif s < 273.5:
T = (T4 - T3) / (273.5 - 268.5) * (s - 268.5) + T3
elif s < 339.5:
T = T4
elif s < 344.5:
T = (25 - T4) / (344.5 - 339.5) * (s - 339.5) + T4
else:
T = 25
return T
def compu_loss(): #计算损失函数的值
delta_T = T[75, 38:747] - fj_y[0:709] #根据不同的alpha选择不同区间
L = 1/2 * np.sum(delta_T * delta_T)
return L
def compu_T(): #计算当前的温度T
for i in range(0, len(alpha)):
A = alpha[i] * delta_t / 2 / (delta_x * delta_x)
Omega = np.eye(size_x, size_x) # 构建Omega矩阵
Omega[0, 0] = 1 + sigma * delta_x
Omega[0, 1] = -1
Omega[size_x - 1, size_x - 2] = -1
Omega[size_x - 1, size_x - 1] = 1 + sigma * delta_x
for k in range(1, size_x - 1):
Omega[k, k - 1] = -A
Omega[k, k] = 2 * A + 1
Omega[k, k + 1] = -A
B = np.zeros((size_x, 1)) # 构建B矩阵
for j in range(t_interval[i], t_interval[i + 1]):
s = j * delta_t * v # 当前前进距离
T_air = Tair_s(s)
B[0] = sigma * delta_x * T_air
B[size_x - 1] = sigma * delta_x * T_air
for j1 in range(1, size_x - 1):
B[j1] = T[j1, j - 1]
T[:, [j]] = np.linalg.inv(Omega).dot(B) # 注意矩阵乘法使用dot,以及[j]，后者影响维度
file_path = '附件.xls' #读取附件
data_1 = pd.read_excel(file_path, sheet_name = "Sheet1")
data_2 = data_1.values #把dataframe(data_1)转化为ndarray(data_2)
fj_x = data_2[:, 0]
fj_y = data_2[:, 1]
accuracy = 1e-11
alpha_op = [] #alpha的最优值
L_min = 1e99 #损失函数的最小值
alpha = [0,0,0,0,0]
flag = 1
for alpha1 in range(9, 10):
alpha1 = alpha1 * accuracy
alpha[0] = alpha1
for alpha2 in range(11, 12):
alpha2 = alpha2 * accuracy
alpha[1] = alpha2
for alpha3 in range(15, 16):
alpha3 = alpha3 * accuracy
alpha[2] = alpha3
for alpha4 in range(10, 11):
alpha4 = alpha4 * accuracy
alpha[3] = alpha4
for alpha5 in range(1, 20):
alpha5 = alpha5 * accuracy
alpha[4] = alpha5
print(flag)
flag = flag+1
compu_T()
L = compu_loss()
if L < L_min:
L_min = L
#alpha_op = alpha #不能直接等于，否则alpha_op和alpha共用一块内存,一个变化，另一个也变
alpha_op = alpha.copy()
print(alpha_op)
alpha = alpha_op
compu_T()
# ##画三维图
# x = np.linspace(-thickness/2, thickness/2, size_x) #厚度
y = np.linspace(0, t_m, size_t) #时间
#
# X, Y = np.meshgrid(x, y)
#
# fig1 = plt.figure()
# ax = plt.axes(projection='3d') #ax是一个变量，像是str， ndarray等
# ax.plot_surface(X, Y, T.T, cmap ='jet', edgecolor ='green') #T需要转置一下，否则维度对不上
# ax.set_xlabel('x', fontsize = 20)
# ax.set_ylabel('t', fontsize = 20)
# ax.set_zlabel('T', fontsize = 20);
###画二维图
fig2 = plt.figure()
plot1 = plt.plot(y[38:], T[75, 38:], color = 'green', label='模型值')
plot2 = plt.plot(fj_x, fj_y, color = 'red', label = '附件值')
plt.title('工件中心温度随时间的变化关系', fontsize = 20)
plt.legend(loc = 1, fontsize = 20)
plt.xlabel('t', fontsize = 20)
plt.ylabel('T', fontsize = 20)
plt.show()
3. 模型求解
3.1 问题一
把给出的温区温度和求出的参数代入模型中。
要求的几个点的温度距离起点分别为：
[
111.25
,
217.75
,
253.25
,
304
]
c
m
[111.25,217.75,253.25,304]cm
[111.25,217.75,253.25,304]cm
除以速度(1.3cm/s)得到从起点出发的时间为：
[
85.57692308
,
167.5
,
194.80769231
,
233.84615385
]
s
[ 85.57692308, 167.5 , 194.80769231, 233.84615385]s
[85.57692308,167.5,194.80769231,233.84615385]s
乘以2取整得到
T
T
T对应下标的数据即为答案：
(
T
[
171
]
,
T
[
335
]
,
T
[
390
]
,
T
[
468
]
)
(T[171], T[335] , T[390], T[468])
(T[171],T[335],T[390],T[468])
为：
[
130.3
,
166.2
,
184.2
,
228
]
[130.3,166.2,184.2,228]
[130.3,166.2,184.2,228]
温度曲线如下所示：
3.2 问题二
思路如下：
把第二题各个温区的温度代入模型中，从小到大遍历速度，每对应一个速度，得到 T [ 76 , : ] T[76, :] T[76,:]，检验其是否符合制程界限的条件，如果符合，增加速度，这样直到找到满足制程界限的最大速度，注意 T [ 76 ] [ i ] T[76][i] T[76][i]和 T [ 76 ] [ i + 1 ] T[76][i+1] T[76][i+1]之间的间隔是0.5s。
最终得到结果在72cm/min附近，四篇有论文中有三篇在76，77附近，结果还是有点差异的。
5个约束对应的四张图如下：
求解代码如下：
#问题2，求满足制程界限的最大允许速度
#有限元差分求解代码
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
# from pylab import *
#
# mpl.rcParams['font.sans-serif'] = ['SimHei']
##固定参数
delta_x = 1e-6 #单位m
delta_t = 0.5 #单位s
thickness = 0.15e-3 #元件厚度，m
t_m = 373 #附件中的总时间
#v = 78 / 60 #传送带前进速度
T1 = 182 #小温区1~5
T2 = 203 #6
T3 = 237 #7
T4 = 254 #8~9
size_x = round(thickness / delta_x + 1)
size_t = round(t_m / delta_t + 1)
##待求参数，先根据[2]求参的结果，画图检验一下
sigma = 8657485
alpha = [9e-11,11e-11,15e-11,10e-11,5e-11] #第1个温区的alpha参数，以此类推
t_interval = [1, 348, 409, 470, 592, 747] #5个温区的截止时间，以此为第1个温区的……
##待求值
T = np.zeros((size_x, size_t))
T[:,0] = 25
def Tair_s(s): #当前距离和空气温度的函数
if s < 25:
T = (T1 - 25) / (25 - 0) * (s - 0) + 25
elif s < 197.5:
T = T1
elif s < 202.5:
T = (T2 - T1) / (202.5 - 197.5) * (s - 197.5) + T1
elif s < 233:
T = T2
elif s < 238:
T = (T3 - T2) / (238 - 233) * (s - 233) + T2
elif s < 268.5:
T = T3
elif s < 273.5:
T = (T4 - T3) / (273.5 - 268.5) * (s - 268.5) + T3
elif s < 339.5:
T = T4
elif s < 344.5:
T = (25 - T4) / (344.5 - 339.5) * (s - 339.5) + T4
else:
T = 25
return T
def if_satisfy(): #是否满足制程界限
max_rising_slope = 0
max_descending_slope = 0
peak_value = 0
i150begin = 0
i190end = 0
i217begin = 0
i217end = 0
for i in range(1, len(T[76,:])):
if T[76, i] - T[76, i-1] > max_rising_slope:
max_rising_slope = T[76, i] - T[76, i-1]
elif T[76, i] - T[76, i-1] < max_descending_slope:
max_descending_slope = T[76, i] - T[76, i-1]
if T[76, i] > peak_value:
peak_value = T[76, i]
if T[76, i-1] < 150 and T[76, i] >= 150:
i150begin = i
if T[76, i-1] < 190 and T[76, i] >= 190:
i190end = i
if T[76, i - 1] < 217 and T[76, i] >= 217:
i217begin = i
if T[76, i - 1] > 217 and T[76, i] <= 217:
i217end = i
m_ris_s_l.append(max_rising_slope*2)
m_des_s_l.append(max_descending_slope*2)
t_b_150_190.append((i190end - i150begin)*0.5)
t_u_217.append((i217end - i217begin)*0.5)
peak_l.append(peak_value)
if peak_value <240 or peak_value > 250:
return False
if max_rising_slope > 1.5:
return False
if max_descending_slope < -1.5:
return False
if (i190end - i150begin)*0.5 > 120 or (i190end - i150begin)*0.5 < 60:
return False
if (i217end - i217begin)*0.5 > 90 or (i217end - i217begin)*0.5 < 40:
return False
return True
V_list = []
m_ris_s_l = []
m_des_s_l = []
t_b_150_190 = []
t_u_217 = []
peak_l = []
for V in range(70, 101):
V_list.append(V)
v = V/60
for i in range(0, len(alpha)):
A = alpha[i] * delta_t / 2 / (delta_x * delta_x)
Omega = np.eye(size_x, size_x) #构建Omega矩阵
Omega[0, 0] = 1 + sigma * delta_x
Omega[0, 1] = -1
Omega[size_x-1, size_x-2] = -1
Omega[size_x-1, size_x-1] = 1 + sigma * delta_x
for k in range(1, size_x - 1):
Omega[k, k-1] = -A
Omega[k, k] = 2 * A + 1
Omega[k, k+1] = -A
B = np.zeros((size_x, 1)) #构建B矩阵
for j in range(t_interval[i], t_interval[i+1]):
s = j * delta_t * v #当前前进距离
T_air = Tair_s(s)
B[0] = sigma * delta_x * T_air
B[size_x - 1] = sigma * delta_x * T_air
for j1 in range(1, size_x - 1):
B[j1] = T[j1, j-1]
T[:, [j]] = np.linalg.inv(Omega).dot(B) #注意矩阵乘法使用dot,以及[j]，后者影响维度
if if_satisfy() == False:
print(V)
fig1 = plt.figure()
plt.plot(V_list, m_ris_s_l, color = "green", label = "最大温度上升斜率")
plt.plot(V_list, m_des_s_l, color = "red", label = "最小温度下降斜率")
plt.legend(loc = 1, fontsize = 20)
plt.title('最大上升斜率和最小下降斜率和速度的关系', fontsize = 20)
plt.xlabel('V(cm/min)', fontsize = 20)
plt.ylabel('delta_T', fontsize = 20)
fig2 = plt.figure()
plt.plot(V_list, t_b_150_190, color= 'red')
plt.title('温度在150-190的时间和速度的关系', fontsize = 20)
plt.xlabel('V', fontsize = 20)
plt.ylabel('t', fontsize = 20)
fig3 = plt.figure()
plt.plot(V_list, t_u_217, color= 'red')
plt.title('温度大于217的时间和速度的关系', fontsize = 20)
plt.xlabel('V', fontsize = 20)
plt.ylabel('t', fontsize = 20)
fig4 = plt.figure()
plt.plot(V_list, peak_l, color = 'red')
plt.title('温度峰值和速度的关系', fontsize = 20)
plt.xlabel('V', fontsize = 20)
plt.ylabel('T', fontsize = 20)
3.3 问题三
将求解代码封装为类，改写如下，可以看到，代码清楚多了。
#问题三，模拟退火求最小面积
#有限元差分求解代码
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
###封装成类，面向对象编程
class LuWenQuXian: #拼音“炉温曲线” (^_^)
def __init__(self, T1, T2, T3, T4, V):
self.v = V / 60 #单位cm/s
self.T1 = T1
self.T2 = T2
self.T3 = T3
self.T4 = T4
self.delta_x = 1e-6
self.delta_t = 0.5
self.thickness = 0.15e-3
self.t_m = 373
self.size_x = round(self.thickness / self.delta_x + 1)
self.size_t = round(self.t_m / self.delta_t + 1)
self.sigma = 8657485
self.alpha = [9e-11,11e-11,15e-11,10e-11,5e-11]
self.t_interval = [1, 348, 409, 470, 592, 747]
self.T = np.zeros((self.size_x, self.size_t))
self.T[:,0] = 25
def Tair_s(self, s):
if s < 25:
T = (self.T1 - 25) / (25 - 0) * (s - 0) + 25
elif s < 197.5:
T = self.T1
elif s < 202.5:
T = (self.T2 - self.T1) / (202.5 - 197.5) * (s - 197.5) + self.T1
elif s < 233:
T = self.T2
elif s < 238:
T = (self.T3 - self.T2) / (238 - 233) * (s - 233) + self.T2
elif s < 268.5:
T = self.T3
elif s < 273.5:
T = (self.T4 - self.T3) / (273.5 - 268.5) * (s - 268.5) + self.T3
elif s < 339.5:
T = self.T4
elif s < 344.5:
T = (25 - self.T4) / (344.5 - 339.5) * (s - 339.5) + self.T4
else:
T = 25
return T
def if_satisfy(self): # 是否满足制程界限
max_rising_slope = 0
max_descending_slope = 0
peak_value = 0
i150begin = 0
i190end = 0
i217begin = 0
i217end = 0
for i in range(1, len(self.T[76, :])):
if self.T[76, i] - self.T[76, i - 1] > max_rising_slope:
max_rising_slope = self.T[76, i] - self.T[76, i - 1]
elif self.T[76, i] - self.T[76, i - 1] < max_descending_slope:
max_descending_slope = self.T[76, i] - self.T[76, i - 1]
if self.T[76, i] > peak_value:
peak_value = self.T[76, i]
if self.T[76, i - 1] < 150 and self.T[76, i] >= 150:
i150begin = i
if self.T[76, i - 1] < 190 and self.T[76, i] >= 190:
i190end = i
if self.T[76, i - 1] < 217 and self.T[76, i] >= 217:
i217begin = i
if self.T[76, i - 1] > 217 and self.T[76, i] <= 217:
i217end = i
if peak_value < 240 or peak_value > 250:
return False
if max_rising_slope > 1.5:
return False
if max_descending_slope < -1.5:
return False
if (i190end - i150begin) * 0.5 > 120 or (i190end - i150begin) * 0.5 < 60:
return False
if (i217end - i217begin) * 0.5 > 90 or (i217end - i217begin) * 0.5 < 40:
return False
return True
def compu_T(self): # 计算当前的温度T
for i in range(0, len(self.alpha)):
A = self.alpha[i] * self.delta_t / 2 / (self.delta_x * self.delta_x)
Omega = np.eye(self.size_x, self.size_x) # 构建Omega矩阵
Omega[0, 0] = 1 + self.sigma * self.delta_x
Omega[0, 1] = -1
Omega[self.size_x - 1, self.size_x - 2] = -1
Omega[self.size_x - 1, self.size_x - 1] = 1 + self.sigma * self.delta_x
for k in range(1, self.size_x - 1):
Omega[k, k - 1] = -A
Omega[k, k] = 2 * A + 1
Omega[k, k + 1] = -A
B = np.zeros((self.size_x, 1)) # 构建B矩阵
for j in range(self.t_interval[i], self.t_interval[i + 1]):
s = j * self.delta_t * self.v # 当前前进距离
T_air = self.Tair_s(s)
B[0] = self.sigma * self.delta_x * T_air
B[self.size_x - 1] = self.sigma * self.delta_x * T_air
for j1 in range(1, self.size_x - 1):
B[j1] = self.T[j1, j - 1]
self.T[:, [j]] = np.linalg.inv(Omega).dot(B) # 注意矩阵乘法使用dot,以及[j]，后者影响维度
def plot_T_xANDt(self): #绘制温度和时间已经厚度关系的三维图
##画三维图
x = np.linspace(-self.thickness / 2, self.thickness / 2, self.size_x) # 厚度
y = np.linspace(0, self.t_m, self.size_t) # 时间
X, Y = np.meshgrid(x, y) #生成网格
plt.figure()
ax = plt.axes(projection='3d') # ax是一个变量，像是str， ndarray等
ax.plot_surface(X, Y, self.T.T, cmap='jet', edgecolor='green') # T需要转置一下，否则维度对不上
ax.set_xlabel('x', fontsize=20)
ax.set_ylabel('t', fontsize=20)
ax.set_zlabel('T', fontsize=20);
def plot_T_t(self): #绘制元件中心温度和时间关系的二维图
###画二维图
y = np.linspace(0, self.t_m, self.size_t) # 时间
plt.figure()
plt.plot(y[38:], self.T[75, 38:], color='green', label='模型值')
plt.title('工件中心温度随时间的变化关系', fontsize=20)
plt.legend(loc=1, fontsize=20)
plt.xlabel('t', fontsize=20)
plt.ylabel('T', fontsize=20)
plt.show()
if __name__ == '__main__':
luwenquxian1 = LuWenQuXian(182, 203, 237, 254, 78)
luwenquxian1.compu_T()
luwenquxian1.plot_T_xANDt()
luwenquxian1.plot_T_t()
接下来编写“模拟退火类”，使用退火算法求解问题，其中5个参数应该满足的条件是：
{
165
≤
T
1
≤
185
185
≤
T
2
≤
205
225
≤
T
1
≤
245
245
≤
T
1
≤
265
65
≤
v
≤
100
\begin{cases} 165\leq T_1 \leq 185\\[1ex] 185\leq T_2 \leq 205\\[1ex] 225\leq T_1 \leq 245\\[1ex] 245\leq T_1 \leq 265\\[1ex] 65\leq v \leq 100\\[1ex] \end{cases}
⎩⎪⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎪⎪⎧165≤T1≤185185≤T2≤205225≤T1≤245245≤T1≤26565≤v≤100
除此之外，还应该满足制程界限的约束。
完整代码如下：
#问题三，模拟退火求最小面积
#有限元差分求解代码
import math
import numpy as np
import matplotlib.pyplot as plt
from random import random
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
###封装成类，面向对象编程
class LuWenQuXian: #拼音“炉温曲线” (^_^)
def __init__(self, T1, T2, T3, T4, V):
self.v = V / 60 #单位cm/s
self.T1 = T1
self.T2 = T2
self.T3 = T3
self.T4 = T4
self.delta_x = 1e-6
self.delta_t = 0.5
self.thickness = 0.15e-3
self.t_m = 373
self.size_x = round(self.thickness / self.delta_x + 1)
self.size_t = round(self.t_m / self.delta_t + 1)
self.sigma = 8657485
self.alpha = [9e-11,11e-11,15e-11,10e-11,5e-11]
self.t_interval = [1, 348, 409, 470, 592, 747]
self.T = np.zeros((self.size_x, self.size_t))
self.T[:,0] = 25
def Tair_s(self, s):
if s < 25:
T = (self.T1 - 25) / (25 - 0) * (s - 0) + 25
elif s < 197.5:
T = self.T1
elif s < 202.5:
T = (self.T2 - self.T1) / (202.5 - 197.5) * (s - 197.5) + self.T1
elif s < 233:
T = self.T2
elif s < 238:
T = (self.T3 - self.T2) / (238 - 233) * (s - 233) + self.T2
elif s < 268.5:
T = self.T3
elif s < 273.5:
T = (self.T4 - self.T3) / (273.5 - 268.5) * (s - 268.5) + self.T3
elif s < 339.5:
T = self.T4
elif s < 344.5:
T = (25 - self.T4) / (344.5 - 339.5) * (s - 339.5) + self.T4
else:
T = 25
return T
def if_satisfy(self): # 是否满足制程界限
max_rising_slope = 0
max_descending_slope = 0
peak_value = 0
i150begin = 0
i190end = 0
i217begin = 0
i217end = 0
for i in range(1, len(self.T[76, :])):
if self.T[76, i] - self.T[76, i - 1] > max_rising_slope:
max_rising_slope = self.T[76, i] - self.T[76, i - 1]
elif self.T[76, i] - self.T[76, i - 1] < max_descending_slope:
max_descending_slope = self.T[76, i] - self.T[76, i - 1]
if self.T[76, i] > peak_value:
peak_value = self.T[76, i]
if self.T[76, i - 1] < 150 and self.T[76, i] >= 150:
i150begin = i
if self.T[76, i - 1] < 190 and self.T[76, i] >= 190:
i190end = i
if self.T[76, i - 1] < 217 and self.T[76, i] >= 217:
i217begin = i
if self.T[76, i - 1] > 217 and self.T[76, i] <= 217:
i217end = i
if peak_value < 240 or peak_value > 250:
return False
if max_rising_slope > 1.5:
return False
if max_descending_slope < -1.5:
return False
if (i190end - i150begin) * 0.5 > 120 or (i190end - i150begin) * 0.5 < 60:
return False
if (i217end - i217begin) * 0.5 > 90 or (i217end - i217begin) * 0.5 < 40:
return False
return True
def if_param_in_range(self): #判断5个参数是否在范围内
if self.T1 < 165 or self.T1 > 185:
return False
if self.T2 < 185 or self.T2 > 205:
return False
if self.T3 < 225 or self.T1 > 245:
return False
if self.T2 < 245 or self.T2 > 265:
return False
if self.v < 65 or self.v > 100:
return False
return True
def compu_T(self): # 计算当前的温度T
for i in range(0, len(self.alpha)):
A = self.alpha[i] * self.delta_t / 2 / (self.delta_x * self.delta_x)
Omega = np.eye(self.size_x, self.size_x) # 构建Omega矩阵
Omega[0, 0] = 1 + self.sigma * self.delta_x
Omega[0, 1] = -1
Omega[self.size_x - 1, self.size_x - 2] = -1
Omega[self.size_x - 1, self.size_x - 1] = 1 + self.sigma * self.delta_x
for k in range(1, self.size_x - 1):
Omega[k, k - 1] = -A
Omega[k, k] = 2 * A + 1
Omega[k, k + 1] = -A
B = np.zeros((self.size_x, 1)) # 构建B矩阵
for j in range(self.t_interval[i], self.t_interval[i + 1]):
s = j * self.delta_t * self.v # 当前前进距离
T_air = self.Tair_s(s)
B[0] = self.sigma * self.delta_x * T_air
B[self.size_x - 1] = self.sigma * self.delta_x * T_air
for j1 in range(1, self.size_x - 1):
B[j1] = self.T[j1, j - 1]
self.T[:, [j]] = np.linalg.inv(Omega).dot(B) # 注意矩阵乘法使用dot,以及[j]，后者影响维度
def compu_S(self): #计算第三问中要求的面积
begini = 0
endi = 0
peak_value = 0
for i in range(1, len(self.T[76, :])):
if self.T[76, i-1] < 217 and self.T[76, i] >= 217:
begini = i
if self.T[76, i] > peak_value:
peak_value = self.T[76, i]
else:
endi = i
break
if peak_value <= 217:
return 0
S = 0
for i in range(begini, endi):
S += (self.T[76, i] - 217 + self.T[76, i+1] - 217) / 2 * self.delta_t
return S
def plot_T_xANDt(self): #绘制温度和时间已经厚度关系的三维图
##画三维图
x = np.linspace(-self.thickness / 2, self.thickness / 2, self.size_x) # 厚度
y = np.linspace(0, self.t_m, self.size_t) # 时间
X, Y = np.meshgrid(x, y) #生成网格
plt.figure()
ax = plt.axes(projection='3d') # ax是一个变量，像是str， ndarray等
ax.plot_surface(X, Y, self.T.T, cmap='jet', edgecolor='green') # T需要转置一下，否则维度对不上
ax.set_xlabel('x', fontsize=20)
ax.set_ylabel('t', fontsize=20)
ax.set_zlabel('T', fontsize=20);
def plot_T_t(self): #绘制元件中心温度和时间关系的二维图
###画二维图
y = np.linspace(0, self.t_m, self.size_t) # 时间
plt.figure()
plt.plot(y[38:], self.T[75, 38:], color='green', label='模型值')
plt.title('工件中心温度随时间的变化关系', fontsize=20)
plt.legend(loc=1, fontsize=20)
plt.xlabel('t', fontsize=20)
plt.ylabel('T', fontsize=20)
plt.show()
class SA:
def __init__(self, iter = 5, T0 = 5, Tf = 0.01, alpha = 0.9):
self.iter = iter
self.alpha = alpha
self.T0 = T0
self.Tf = Tf
self.T = T0
self.param = np.zeros((self.iter, 5))
for i in range(self.iter):
T1 = random() * 20 + 165 #165 <= T1 <= 185
T2 = random() * 20 + 185
T3 = random() * 20 + 225
T4 = random() * 20 + 245
v = random() * 35 + 65
self.param[i, 0] = T1
self.param[i, 1] = T2
self.param[i, 2] = T3
self.param[i, 3] = T4
self.param[i, 4] = v
self.history = {'S':[], 'T': []}
def if_param_new_in_ran(self, param_new): #判断generate_new生成的新参数是否在范围内
if param_new[0, 0] < 165 or param_new[0, 0] > 185:
return False
if param_new[0, 1] < 185 or param_new[0, 1] > 205:
return False
if param_new[0, 2] < 225 or param_new[0, 2] > 245:
return False
if param_new[0, 3] < 245 or param_new[0, 3] > 265:
return False
if param_new[0, 4] < 65 or param_new[0, 4] > 100:
return False
return True
def generate_new(self, param_old):
param_new = np.zeros(param_old.shape)
while self.if_param_new_in_ran(param_new) is not True:
param_new[0, 0] = param_old[0, 0] + self.T * (random() - random())
param_new[0, 1] = param_old[0, 1] + self.T * (random() - random())
param_new[0, 2] = param_old[0, 2] + self.T * (random() - random())
param_new[0, 3] = param_old[0, 3] + self.T * (random() - random())
param_new[0, 4] = param_old[0, 4] + self.T * (random() - random())
return param_new
def Metropolis(self, f, f_new):
if f_new <= f:
return True
else:
p = math.exp((f - f_new) / self.T)
if random() < p:
return True
else:
return False
def best(self):
f_list = []
for i in range(self.iter):
luwenquxian = LuWenQuXian(self.param[i, 0], self.param[i, 1], self.param[i, 2], self.param[i, 3], self.param[i, 4])
luwenquxian.compu_T()
f = luwenquxian.compu_S()
f_list.append(f)
f_best = min(f_list)
idx = f_list.index(f_best)
return f_best, idx
def run(self):
count = 0
while self.T > self. Tf:
print(count)
for i in range(self.iter):
print(i,)
luwenquxian = LuWenQuXian(self.param[i, 0], self.param[i, 1], self.param[i, 2], self.param[i, 3],
self.param[i, 4])
luwenquxian.compu_T()
f = luwenquxian.compu_S()
print('flag1')
param_new = self.generate_new(self.param[[i], :]) #没有输出‘flag2’，说明卡在这里了
print('flag2')
luwenquxian1 = LuWenQuXian(param_new[0, 0], param_new[0, 1], param_new[0, 2], param_new[0, 3], param_new[0, 4])
luwenquxian1.compu_T()
f_new = luwenquxian1.compu_S()
if self.Metropolis(f, f_new):
self.param[[i], :] = param_new
ft, _ = self.best()
self.history['S'].append(ft)
self.history['T'].append(self.T)
self.T = self.T * self.alpha
count += 1
f_best, idx = self.best()
print(f"S = {f_best}, T1 = {self.param[idx, 0]}, T2 = {self.param[idx, 1]}, T3 = {self.param[idx, 2]}, T4 = {self.param[idx, 3]}, v = {self.param[idx, 4]}")
if __name__ == '__main__':
sa = SA()
sa.run()
# luwenquxian1 = LuWenQuXian(182, 203, 237, 254, 78)
# luwenquxian1.compu_T()
# luwenquxian1.plot_T_xANDt()
# luwenquxian1.plot_T_t()
运行效果还是有些出入的，相差比较大，目前没有发现问题，欢迎读者交流。
3.4 问题四
重新建立目标函数，仍然使用模拟退火算法，修改以上代码即可。
4. 总结与反思
这篇文章断断续续写了总共有三四天的样子，从看懂别人的模型到自己编程实现，从头到尾走了一遍整个过程，工作量很大，但是收获也是很大，同时也感觉自己能力上的不足，希望能对读者有帮助。
参考资料
[1] 优秀论文A195
[1] 优秀论文A195
[2] 优秀论文A070
[3] 通过Python绘制分段函数
[4] Three-dimensional Plotting in Python using Matplotlib
[5] Three-Dimensional Plotting in Matplotlib
[6] 优秀论文A212
4万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
