数模国赛2020A题

题解 Python, 数模 评论

第一问

忽略电路板厚度，假定温区内部恒温

通过热传导方程得到不同温区之间的空隙中，温度在空间上线性分布

通过牛顿冷却定律和比热容公式解微分方程可知，物体在恒温环境下温度随时间的函数

然后由附件数据分段拟合参数

下面是拟合的代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
import numpy as np
import pylab as plt
import pandas as pd
from scipy.optimize import curve_fit
exc = pd.read_excel("input.xlsx", usecols=range(0, 2))
data = exc.values


pos = np.array([[0, 0], [25, 55.5], [60.5, 91], [96, 126.5], [131.5, 162], [167, 197.5], [202.5, 233], [238, 268.5], [273.5, 304], [309, 339.5], [344.5, 375], [380, 410.5], [435.5, 460.5]]
               )

ranx = np.array([[0, pos[2][0]], [pos[2][0], pos[5][1]], [pos[5][1], pos[6][0]],
                 [pos[6][0], pos[6][1]], [pos[6][1], pos[7][0]], [
                     pos[7][0], pos[7][1]],
                 [pos[7][1], pos[8][0]], [pos[8][0], pos[9][1]], [
                     pos[9][1], pos[10][1]],
                 [pos[10][1], pos[12][0]]])


def getab(i, tem):  # 计算第i个小温区后的温度曲线T=ax+b
    x1 = pos[i][1]
    x2 = pos[i+1][0]
    t1 = tem[i]
    t2 = tem[i+1]
    a = (t2-t1)/(x2-x1)
    b = t1-a*x1
    return (a, b)


def getfx(x, tem):  # 获取x位置的环境温度曲线f(x)=ax+b
    for i in range(0, 12):
        if (pos[i][0] <= x and x <= pos[i][1]):
            return (0, tem[i])
        elif (pos[i][1] <= x and x <= pos[i+1][0]):
            return getab(i, tem=tem)
    return (-1, -1)


def f(tt, lmd):  # 用于拟合的函数
    def gettem(t):  # 获取x时刻的零件温度，已知初值T(t0)=T0
        # a, b = (0, 175)
        a, b = getfx(t[0]*v, tem=tem)
        c = (T0-a*v*t0+a*v/lmd-b)*np.exp(lmd*t0)
        return a*v*t-a*v/lmd+b+c*np.exp(-lmd*t)
    return gettem(tt)


if __name__ == '__main__':
    tem = np.array([25, 173, 173, 173, 173, 173,
                    198, 230, 257, 257, 25, 25, 25])
    v = 70/60  # 传送带速度，单位cm/s

    tT = data
    t = tT.T[0]
    T = tT.T[1]
    lmds = []
    TT = []
    x = t*v
    lj, j = 0, 0
    for ranxi in ranx:
        while (j < len(x) and x[j] <= ranxi[1]):
            j += 1
        t0 = t[lj]
        T0 = T[lj]
        popt, pcov = curve_fit(f, t[lj:j], T[lj:j])
        lmds.extend(popt)
        TT.extend(f(t[lj:j], lmd=popt))
        lj = j
        if (j >= len(x)):
            break

    plt.scatter(t, T, s=1)
    plt.plot(t, TT, color='r')
    plt.savefig('getlmd.png', dpi=500)
    plt.show()
    print(lmds)

第一问代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
import numpy as np
import pylab as plt
import pandas as pd
from getlmd import ranx, pos, getfx


# lmd = 0.01442337549312217
lmds = [0.008643609716781697, 0.02103327687676206, 0.010922351732360641, 0.014168584576603811, 0.01604373036324419,
        0.0269012209873633, 0.023523032208845447, 0.0190525180588788, 0.00079443137454158, 0.009429030064290005]


def gettem(t, t0, T0, lmd, v, tem):  # 获取t时段的零件温度，已知初值T(t0)=T0
    a, b = getfx(t[0]*v, tem=tem)
    c = (T0-a*v*t0+a*v/lmd-b)*np.exp(lmd*t0)
    return a*v*t-a*v/lmd+b+c*np.exp(-lmd*t)


def gettT(v, tem, dt=0.5, fil=lambda x: x[1] >= 30):  # 计算温度曲线
    t = np.arange(dt, pos[-1][0]/v, dt)
    T = []
    x = t*v
    lj, j = 0, 0
    t0, T0 = 0, 25
    for (i, ranxi) in enumerate(ranx):
        while (j < len(x) and x[j] <= ranxi[1]):
            j += 1
        Ti = gettem(t[lj:j], t0, T0, lmds[i], v=v, tem=tem)
        T.extend(Ti)
        lj = j
        t0, T0 = t[j-1], T[j-1]
        if (j >= len(x)):
            break

    tT = np.array([t, T]).T
    if (fil != None):
        tT_f = filter(fil, tT)
        tT = np.array(list(tT_f))
        t, T = tT.T
    return (t, T, tT)


if __name__ == '__main__':
    tem = np.array([25, 173, 173, 173, 173, 173,
                    198, 230, 257, 257, 25, 25, 25])
    v = 78/60  # 传送带速度，单位cm/s

    t, T, tT = gettT(v=v, tem=tem)
    plt.plot(t, T)
    plt.savefig('sub1.png', dpi=500)
    plt.show()
    res = pd.DataFrame(tT, columns=["时间", "温度"])
    res.to_csv('result.csv', index=False, header=False)
    printtim = [(pos[3][0]+pos[3][1])/2/v, (pos[6][0]+pos[6][1]) /
                2/v, (pos[7][0]+pos[7][1])/2/v, pos[8][1]/v]
    lmdid = [1, 3, 5, 7]
    j = 1
    for (i, timei) in enumerate(printtim):
        while (j < len(t) and t[j] < timei):
            j += 1
        if j >= len(t):
            break
        print(gettem(np.array([timei]), t[j-1], T[j-1],
                     lmd=lmds[lmdid[i]], v=v, tem=tem))

得到的各点温度：

1
2
3
4
[135.71310782]
[169.75603785]
[191.00134962]
[223.82106867]

炉温曲线：

第二问

大力枚举速度v，先大范围粗糙枚举，再小范围精确枚举

需要复用许多第一问的代码，所以函数化编程很重要，同时也要封装好来

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
import numpy as np
import pylab as plt
from sub1 import gettT


def checkv(t, T, tT, dt=0.5):
    maxT = np.max(T)
    if not(240 <= maxT and maxT <= 250):
        return False
    for i in range(1, len(T)):
        if np.abs((T[i]-T[i-1])/dt) > 3:
            print((i, T[i], T[i-1], t[i], t[i-1]))
            return False
        if (T[i-1] < 150 and T[i] >= 150):
            t1 = t[i]
        if (T[i-1] < 190 and T[i] >= 190):
            t2 = t[i]
        if (T[i-1] < 217 and T[i] >= 217):
            t3 = t[i]
        if (T[i-1] >= 217 and T[i] < 217):
            t4 = t[i]
    if not (60 <= t2-t1 and t2-t1 <= 120):
        return False
    if not (40 <= t4-t3 and t4-t3 <= 90):
        return False
    return True


def printlim(t, T, tT, dt=0.5):  # 输出各制程的极限值
    maxT = np.max(T)
    print("maxT=", maxT)
    maxk = 0
    mink = 0
    for i in range(1, len(T)):
        maxk = max((T[i]-T[i-1])/dt, maxk)
        mink = min((T[i]-T[i-1])/dt, mink)
        if (T[i-1] < 150 and T[i] >= 150):
            t1 = t[i]
        if (T[i-1] < 190 and T[i] >= 190):
            t2 = t[i]
        if (T[i-1] < 217 and T[i] >= 217):
            t3 = t[i]
        if (T[i-1] >= 217 and T[i] < 217):
            t4 = t[i]
    print("maxk=", maxk)
    print("mink=", mink)
    print("t2-t1=", t2-t1)
    print("t4-t3=", t4-t3)


if __name__ == '__main__':
    tem = np.array([25, 182, 182, 182, 182, 182,
                    203, 237, 254, 254, 25, 25, 25])
    maxv = 0
    dt = 0.1
    for v in np.linspace(65/60, 100/60, num=35*2+1):  # 间隔0.5
        t, T, tT = gettT(v=v, dt=dt, tem=tem)
        # print("v={}".format(v))
        if (checkv(t, T, tT, dt)):
            maxv = v
    print(maxv*60)
    for v in np.linspace(70/60, 75/60, num=5*100+1):  # 间隔0.01
        t, T, tT = gettT(v=v, dt=dt, tem=tem)
        # print("v={}".format(v))
        if (checkv(t, T, tT, dt)):
            maxv = v
    print(maxv*60)
    t, T, tT = gettT(v=maxv, dt=dt, tem=tem)
    printlim(t, T, tT, dt=dt)
    plt.plot(t, T)
    plt.savefig('sub2.png', dpi=500)
    plt.show()

得到的各参数：

1
2
3
4
5
6
7
v1=72.0
v2=72.33
maxT= 240.00007138887304
maxk= 2.658684542907963
mink= -1.8523553812633509
t2-t1= 103.69999999999999
t4-t3= 82.90000000000006

在最大速度下的炉温曲线：

第三问

把阴影面积作为估值函数，做模拟退火

调了scikit-opt的库

然后还画了最小面积随迭代次数变化的曲线

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
from sko.SA import SA, SABoltzmann, SAFast
from sub1 import gettT
from sub2 import checkv
import numpy as np
import pandas as pd
import pylab as plt
exc = pd.read_excel("input.xlsx", usecols=range(0, 2))
data = exc.values
data = data.T


def getarea(dt, T):
    res = 0
    maxT = max(T)
    id = T.tolist().index(maxT)
    for (i, Ti) in enumerate(T):
        if (Ti >= 217 and i <= id):
            res += dt*(Ti-217)
    return res


def X2v(X):
    return (X[0]*35/20+65)/60


def X2tem(X):
    return [25]+[175+X[1]-10]*5+[195+X[2]-10] + \
           [235+X[3]-10]+[255+X[4]-10]*2+[25]*3


def cost(X):  # 计算模拟退火的代价函数即阴影面积
    if (X.min() < 0 or X.max() > 20):  # 判断输入向量是否合法
        return 1e10+np.abs(X.min()-0)+np.abs(X.max()-20)
    v = X2v(X)
    tem = X2tem(X)
    t, T, tT = gettT(v=v, tem=tem, dt=0.1)
    if (not checkv(t, T, tT, dt=0.1)):  # 判断是否满足制程界限
        return 1e5
    return getarea(dt=0.1, T=T)


if __name__ == '__main__':
    sa = SA(func=cost, x0=[0]*5, lb=0, ub=20, max_stay_counter=5000, L=300)
    best_x, best_y = sa.run()
    print("best_x={}\nminArea={}".format(best_x, best_y))
    v = X2v(best_x)
    tem = X2tem(best_x)
    print("v={}\ntem={}\n".format(v*60, tem))
    t, T, tT = gettT(v=v, tem=tem, dt=0.1)
    plt.plot(t, T)
    plt.show()
    by_f = filter(lambda x: x < 1e3, sa.best_y_history)
    by = list(by_f)
    plt.plot(pd.DataFrame(by).cummin(axis=0))  # 画出最小面积随迭代次数变化的曲线
    plt.show()
    plt.savefig('sub3_2.png', dpi=500)

1
2
3
4
best_x=[1.06929985e+01 1.97873914e+01 5.68152357e+00 8.90629441e-03 1.99983037e+01]
minArea=446.4150721594619
v=83.71274730906958
tem=[25, 184.78739135418277, 184.78739135418277, 184.78739135418277, 184.78739135418277, 184.78739135418277, 190.68152356719315, 225.00890629441466, 264.99830371980664, 264.99830371980664, 25, 25, 25]

第四问

把炉温曲线中相同温度的点的时间做平均值，与峰值时间的平方和作为对称性的衡量标准

然后本问是多目标优化，即阴影面积和对称度

这里采用的是线性组合作为估值函数，经过多次尝试，发现面积变化始终不大，而对称性比较敏感，最终选择了1:10的比例赋权

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
from sko.SA import SA
from sub1 import gettT
from sub2 import checkv
from sub3 import getarea, X2tem, X2v
import numpy as np
import pandas as pd
import pylab as plt
from bisect import bisect
exc = pd.read_excel("input.xlsx", usecols=range(0, 2))
data = exc.values
data = data.T


def getsy(tT):  # 对称度估值
    tT_f = filter(lambda x: x[1] >= 127, tT)
    tT = np.array(list(tT_f))
    t, T = tT.T.tolist()
    maxT = max(T)
    id = T.index(maxT)
    mt = t[id]
    res = 0
    for i in range(id+1, len(T)):
        j = bisect(T, T[i], 0, id)
        res += ((t[i]+t[j])/2-mt)**2
    return res


def cost(X):
    if (X.min() < 0 or X.max() > 20):
        return 1e100+np.abs(X.min()-0)+np.abs(X.max()-20)
    v = X2v(X)
    tem = X2tem(X)
    t, T, tT = gettT(v=v, tem=tem, dt=dt)
    if (not checkv(t, T, tT, dt=dt)):
        return 1e50
    return getarea(dt=dt, T=T)+getsy(tT)*1e1


if __name__ == '__main__':
    dt = 0.1
    sa = SA(func=cost, x0=[0]*5, lb=0, ub=20, max_stay_counter=5000, L=300)
    best_x, best_y = sa.run()
    print("best_x={}\nbest_y={}".format(best_x, best_y))
    v = X2v(best_x)
    tem = X2tem(best_x)
    print("v={}\ntem={}\n".format(v*60, tem))
    t, T, tT = gettT(v=v, tem=tem, dt=0.1)
    print("Area=", getarea(dt, T))
    print("Sy=", getsy(tT))
    plt.plot(t, T)
    plt.plot([0, t.max()], [217, 217], '--')
    plt.savefig('sub4.png', dpi=500)
    plt.show()
    by_f = filter(lambda x: x < 1e30, sa.best_y_history)
    by = list(by_f)
    plt.plot(pd.DataFrame(by).cummin(axis=0))  # 画出最小面积随迭代次数变化的曲线
    plt.show()
    plt.savefig('sub4_2.png', dpi=500)

1
2
3
4
5
6
best_x=[12.12066253  2.29087199  0.94809455 19.92827282 19.99803738]
best_y=167663.96859958651
v=86.21115943067932
tem=[25, 167.2908719854776, 167.2908719854776, 167.2908719854776, 167.2908719854776, 167.2908719854776, 185.94809454727962, 244.92827282463267, 264.998037383657, 264.998037383657, 25, 25, 25]
Area= 446.41859958757425
Sy= 16721.754999999896

linkfqyA link to the past.

咸鱼一条