2023高教社杯数学建模竞赛a题代码
# 此文件用于存储常量
import numpy as np
import math

def cosToSin(x):
    return np.sqrt(1-x**2)

lat = np.deg2rad(39.4) # 维度
lot = np.deg2rad(98.5) # 经度
H = 3 # 海拔高度，单位km
G_0 = 1.366 # 太阳常数
ST = [9, 10.5, 12, 13.5, 15] # 当地时间
omega = [] # 太阳时角
D = [306, 337, 0, 31, 61, 92, 122, 153, 184, 214, 245, 275] # 从春分第0天起算的天数，默认一年365天，春分为3月21日
sin_delta = [] # 太阳赤纬角的sin值
cos_delta = [] # 太阳赤纬角的cos值
sin_alpha = np.zeros((12, 5)) # 太阳高度角的sin值
cos_alpha = np.zeros((12, 5)) # 太阳高度角的cos值，如何确定cos正负？
DNI = np.zeros((12, 5)) # 法向直接辐射辐照度
solar_angle = [] # 太阳入射角
cos_theta = np.zeros((12, 5)) # 太阳入射角
eta_ref = 0.92 # 镜面反射率

for _ in range(5):
    omega.append(np.pi/12*(ST[_]-12))

for _ in range(12):
    sin_delta.append(np.sin(2*np.pi*D[_]/365)*np.sin(2*np.pi*23.45/360))
    cos_delta.append(cosToSin(sin_delta[_]))

for month in range(12):
    for time in range(5):
        sin_alpha[month][time] = cos_delta[month]*math.cos(lat)*math.cos(omega[time]) + sin_delta[month]*math.sin(lat)
        cos_alpha[month][time] = cosToSin(sin_alpha[month][time])

cos_gamma = np.zeros((12, 5))
for month in range(12):
    for time in range(5):
        cos_gamma[month][time] = (sin_delta[month]-sin_alpha[month][time]*math.sin(lat))/(cos_alpha[month][time]*math.cos(lat))

a = 0.4237 - 0.00821*(6-H)**2
b = 0.5055 + 0.00595*(6.5-H)**2
c = 0.2711 + 0.01858*(2.5-H)**2

for month in range(12):
    for time in range(5):
        DNI[month][time] = G_0*(a+b*math.exp(-c/sin_alpha[month][time]))


for month in range(12):
        for time in range(5):
            cos_theta[month][time] = cos_delta[month]*cos_alpha[month][time] # 太阳入射角

max_circle = 25
constance.py
Azreal7/2023mcm_a
import numpy as np
import math
import constance as c

def distance_3d(a:tuple, b:tuple):
    return np.sqrt((a[0]-b[0])**2+(a[1]-b[1])**2+(a[2]-b[2])**2)

def cosToSin(x):
    return np.sqrt(1-x**2)

def distance_2d(a, b):
    return np.sqrt((a[0]-b[0])**2+(a[1]-b[1])**2)

# 适用于month*time*i的三维array求每月平均值
def get_average_mni(input, len):
    average_mni = []
    for month in range(12):
        losses = 0
        for time in range(5):
            for i in range(len):
                losses += input[month][time][i]
        average_mni.append(losses/(5*len))
    return average_mni

# 计算阴影遮挡损失, 输入为：所有定日镜的坐标,cos太阳方位角,cos太阳高度角,镜面高度，吸收塔高度,镜面长度,镜子间距
def get_shadow_effi(x, y, cos_delta, cos_alpha, h1, h2, len1, average_distance):
    al = []
    t = np.zeros((12, 5, len(x)))
    loss = np.zeros((12, 5, len(x)))
    for i in range(len(x)):
        al.append(distance_2d((0, 0), (x[i], y[i])))
        al[i] = np.arctan(h2-h1/al[i]) # 塔与定日镜tan, a1
    for month in range(12):
        for time in range(5):
            for i in range(len(x)):
                t[month][time][i] = (np.pi/2-np.arccos(c.cos_theta[month][time])-al[i])/2 # a2
                x1 = np.pi/2-(t[month][time][i]+al[i]) # a3
                a4 = np.pi/2+np.arccos(c.cos_theta[month][time])-x1
                loss[month][time][i] = average_distance*np.sin(np.pi/2-np.arccos(c.cos_theta[month][time]))/(len1*np.sin(a4))
                if loss[month][time][i] > 1 or int(distance_2d((0, 0), (x[i], y[i]))) == 107:
                    loss[month][time][i] = 1
    return loss, t

# 计算集热器截断效率，输入为集热器的高,集热器的直径,定日镜长,定日镜宽,数据组长度,集热器中心与定日镜中心距离
def get_trunc_effi(h_hotter, L_hotter, h_m, w_m, len, d_HR):
    effi_trunc = np.zeros(len) # 集热器截断效率
    for i in range(len):
        effi_trunc[i] = h_hotter/(h_m+2*d_HR[i]*np.tan(4.65e-3))*L_hotter/(w_m+2*d_HR[i]*np.tan(4.65e-3))
        if effi_trunc[i] > 1:
           effi_trunc[i] = 1
    return effi_trunc

# 计算集热器中心与定日镜中心距离和大气透射率，输入为坐标,集热器重心高度，定日镜高度，数据长度
def get_d_HR(x, y, h1, h2, len):
    eta_at = []
    d_HR = []
    for i in range(len):
        d_HR.append(distance_3d((0, 0, h1), (x[i], y[i], h2)))
        eta_at.append(0.99321-0.0001176*d_HR[i]+1.97e-8*d_HR[i]**2)
    return d_HR, eta_at

# 计算光学效率，输入为阴影遮挡效率，余弦效率，大气透射率，集热器截断效率,数组长度和镜面反射率
def get_light_effi(loss, effi_cos, eta_at, effi_trunc, len, eta_ref=0.92):
    light_effi = np.zeros((12, 5, len)) # 光学效率
    for month in range(12):
        for time in range(5):
            for i in range(len):
                light_effi[month][time][i] = (loss[month][time][i])*effi_cos[month][time][i]*eta_at[i]*effi_trunc[i]*eta_ref
    return light_effi

# 计算余弦效率，输入为t和长度
def get_effi_cos(t, len):
    effi_cos = np.zeros((12, 5, len))
    for month in range(12):
        for time in range(5):
            for i in range(len):
                effi_cos[month][time][i] = math.cos(t[month][time][i])
    return effi_cos

def get_heat_W(light_effi, S, DNI, len, shadow_loss):
    E_field = np.zeros((12, 5))
    for month in range(12):
        for time in range(5):
            for i in range(len):
                E_field[month][time] += S*(shadow_loss[month][time][i])*light_effi[month][time][i]
            E_field[month][time] *= DNI[month][time]
    return E_field

# 根据数量和间距得到每个半径分布数量
def get_nums(w, number, distance):
    min_dis = 5+w # 镜面之间最短距离
    r = 100 # 圆的半径
    nums = []
    while number > 0:
        C_circle = 2*np.pi*r
        num1 = int(C_circle/min_dis)
        if number >= num1:
            nums.append(num1)
            number -= num1
            r += distance
        else:
            nums.append(number)
            number = 0
    return nums

# 根据数量获取坐标系
def get_xy(nums, distance):
    r = 100
    nums_len = len(nums)
    x = []
    y = []
    for i in range(nums_len):
        deg = 0
        deg_plus = 2*math.pi/nums[i]
        for _ in range(nums[i]):
            x.append(np.cos(deg)*r)
            y.append(np.sin(deg)*r)
            deg += deg_plus
        r += distance
    return x, y

# 以下函数为第二题优化版
def new_get_shadow_loss(distance, h_mirror, length, len):
    r = 100
    al = []
    t = np.zeros((12, 5, len))
    loss = np.zeros((12, 5, len))
    for i in range(len):
        al.append(r)
        al[i] = np.arctan(84-h_mirror/al[i]) # alpha
        r += distance
    for month in range(12):
        for time in range(5):
            for i in range(len):
                t[month][time][i] = (np.pi/2-np.arccos(c.cos_theta[month][time])-al[i])/2
                if i == 0:
                    loss[month][time][i] = 1
                    continue
                h1 = distance*np.tan(al[i])
                c0 = h1
                b = length*np.sin(al[i])+length*np.cos(al[i])*np.tan(np.arccos(c.cos_theta[month][time]))
                x_a = (b+c0)/(np.tan(al[i])+np.tan(np.arccos(c.cos_theta[month][time])))*-1
                y_a = np.tan(np.arccos(c.cos_theta[month][time]))*x_a+b
                x_b = -1*length*np.cos(al[i])-distance+length*np.cos(al[i])
                y_b = -1*np.tan(al[i])*length-c0
                D_ab = np.sqrt((x_a-x_b)**2+(y_a-y_b)**2)
                loss[month][time][i] = D_ab/length
                if loss[month][time][i] > 1:
                    loss[month][time][i] = 1
    return loss, t

def new_get_d_HR(h2, len, average_distance):
    r = 100
    eta_at = []
    d_HR = []
    for i in range(len):
        d_HR.append(distance_3d((0, 0, 84), (r, 0, h2)))
        eta_at.append(0.99321-0.0001176*d_HR[i]+1.97e-8*d_HR[i]**2)
        r += average_distance
    return d_HR, eta_at

def new_get_heat_W(light_effi, S, DNI, len, shadow_loss, nums):
    E_field = np.zeros((12, 5))
    for month in range(12):
        for time in range(5):
            for i in range(len):
                E_field[month][time] += S*(shadow_loss[month][time][i])*light_effi[month][time][i]*nums[i]
            E_field[month][time] *= DNI[month][time]
    return E_field

# len:圈数,a_d:圈间距,h_mirror:镜子高度,length:镜子长度,width:镜子宽度, height_diff:两镜高度差
def q3_get_shadow_loss(distance, h_mirror, length, width, height, len=c.max_circle):
    r = 100
    al = []
    t = np.zeros((12, 5, len))
    loss = np.zeros((12, 5, len))
    for i in range(len):
        al.append(r)
        al[i] = np.arctan(84-h_mirror[i]/al[i]) # alpha
        r += distance[i]
    for month in range(12):
        for time in range(5):
            for i in range(len):
                t[month][time][i] = (np.pi/2-np.arccos(c.cos_theta[month][time])-al[i])/2
                if i == 0:
                    loss[month][time][i] = 1
                    continue
                h1 = distance[i-1]*np.tan(al[i])
                c0 = h1 - (height[i]-height[i-1])
                b = length[i]*np.sin(al[i])+length[i]*np.cos(al[i])*np.tan(np.arccos(c.cos_theta[month][time]))
                x_a = (b+c0)/(np.tan(al[i])+np.tan(np.arccos(c.cos_theta[month][time])))*-1
                y_a = np.tan(np.arccos(c.cos_theta[month][time]))*x_a+b
                x_b = -1*length[i]*np.cos(al[i])-distance[i-1]+length[i-1]*np.cos(al[i])
                y_b = -1*np.tan(al[i])*length[i-1]-c0
                D_ab = np.sqrt((x_a-x_b)**2+(y_a-y_b)**2)
                if width[i] > width[i-1]:
                    loss[month][time][i] = D_ab*width[i]/(length[i-1]*width[i-1])
                else:
                    loss[month][time][i] = D_ab/length[i-1]
                if loss[month][time][i] > 1:
                    loss[month][time][i] = 1
    return loss, t

# distance:间距, number:数量, width:宽度
def q3_get_new_nums(distances, number, width):
    r = 100
    nums = np.zeros(c.max_circle)
    for i in range(c.max_circle):
        min_dis = width[i] + 5
        C_circle = 2*np.pi*r
        if number == 0:
            break
        nums[i] = int(C_circle/min_dis)
        if nums[i] > number:
            nums[i] = number
            number = 0
            break
        number -= nums[i]
        r += distances[i]
    return nums

def q3_get_d_HR(height, distances, len=c.max_circle):
    r = 100
    eta_at = []
    d_HR = []
    for i in range(len):
        d_HR.append(distance_3d((0, 0, 84), (r, 0, height[i])))
        eta_at.append(0.99321-0.0001176*d_HR[i]+1.97e-8*d_HR[i]**2)
        r += distances[i]
    return d_HR, eta_at

def q3_get_xy(nums, distances):
    r = 100
    nums_len = len(nums)
    x = []
    y = []
    for i in range(nums_len):
        deg = 0
        if nums[i] == 0:
            break
        deg_plus = 2*math.pi/nums[i]
        for _ in range(int(nums[i])):
            x.append(np.cos(deg)*r)
            y.append(np.sin(deg)*r)
            deg += deg_plus
        r += distances[i]
    return x, y

def q3_get_trunc_effi(h_hotter, L_hotter, h_m, w_m, d_HR, len=c.max_circle):
    effi_trunc = np.zeros(len) # 集热器截断效率
    for i in range(len):
        effi_trunc[i] = h_hotter/(h_m[i]+2*d_HR[i]*np.tan(4.65e-3))*L_hotter/(w_m[i]+2*d_HR[i]*np.tan(4.65e-3))
        if effi_trunc[i] > 1:
           effi_trunc[i] = 1
    return effi_trunc

def q3_get_heat_W(light_effi, S, DNI, len, shadow_loss, nums):
    E_field = np.zeros((12, 5))
    for month in range(12):
        for time in range(5):
            for i in range(len):
                E_field[month][time] += S[i]*(shadow_loss[month][time][i])*light_effi[month][time][i]*nums[i]
            E_field[month][time] *= DNI[month][time]
    return E_field




function.py
Azreal7/2023mcm_a
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import function as f
import constance as c

df = pd.read_excel(r"附件.xlsx")

x = df['x坐标 (m)'].tolist()
y = df['y坐标 (m)'].tolist()

# plt.scatter(x, y)
# plt.show()

d_HR, eta_at = f.get_d_HR(x, y, 84, 4, len(x))

# 此部分用于区分x,y中不同的圆，以计算圆之间的平均距离
distance = []
for i in range(len(x)):
    distance.append(int(f.distance_2d((0, 0), (x[i], y[i]))))
distance = sorted(list(set(distance)))
average_distance = 0 # 圆环之间的平均距离
for i in range(1, len(distance)):
    average_distance += distance[i] - distance[i-1]
average_distance /= len(distance) - 1

# 计算阴影遮挡效率
loss, t = f.get_shadow_effi(x, y, c.cos_delta, c.cos_alpha, 4, 80, 6, average_distance)
# print("loss:" + loss)

# 计算平均阴影遮挡效率
average_effi = f.get_average_mni(loss, len(x))
print( average_effi)

# 计算集热器截断效率
effi_trunc = f.get_trunc_effi(8, 7, 6, 6, len(x), d_HR)
print(np.mean(effi_trunc))

# 计算余弦效率
effi_cos = f.get_effi_cos(t, len(x))
average_cos = f.get_average_mni(effi_cos, len(x))
print(average_cos)

# 求光学效率，镜面反射率默认0.92
light_effi = f.get_light_effi(loss, effi_cos, eta_at, effi_trunc, len(x))
average_light = f.get_average_mni(light_effi, len(x))
print(average_light)

# print(average_light)
E_field = f.get_heat_W(light_effi, 6*6, c.DNI, len(x), loss)
density = []
for i in range(12):
    E = 0
    for time in range(5):
        E += E_field[i][time]
    density.append(E/(36*len(x)*5))
print(density)
q1.py
Azreal7/2023mcm_a
import random
import time
import math
import function as f
import numpy as np
import constance as c
from matplotlib import pyplot as plt
import pandas as pd


class sa_tsp:
    # 接受的输入参数：起始温度start_temp, 结束温度end_temp, 内循环次数num
    def __init__(self, s_t, e_t, num):
        self.start_temp = s_t
        self.end_temp = e_t
        self.num = num

    # 获取当前状态的功率密度
    def get_power_density(self, state):
        power, _ = self.get_power(state)
        return power / (state[0] * state[1] * state[3])

    # 获取当前状态的功率值和摆放的圈数
    def get_power(self, state):
        nums = f.get_nums(state[0], state[3], state[4])
        S = state[0]*state[1]
        x, y = f.get_xy(nums, state[4])
        loss, t = f.new_get_shadow_loss(state[4], 4, 6, len(nums))
        d_HR, eta_at = f.new_get_d_HR(state[2], len(nums), state[4])
        effi_trunc = f.get_trunc_effi(8, 7, state[0], state[1], len(nums), d_HR)
        effi_cos = f.get_effi_cos(t, len(nums))
        light_effi = f.get_light_effi(loss, effi_cos, eta_at, effi_trunc, len(nums))
        E_field = f.new_get_heat_W(light_effi, S, c.DNI, len(nums), loss, nums)
        return np.mean(E_field), len(nums)

    def get_power_and_print(self, state):
        nums = f.get_nums(state[0], state[3], state[4])
        S = state[0]*state[1]
        x, y = f.get_xy(nums, state[4])
        loss, t = f.new_get_shadow_loss(state[4], 4, 6, len(nums))
        f.get_average_mni(loss, len(nums))
        d_HR, eta_at = f.new_get_d_HR(state[2], len(nums), state[4])
        effi_trunc = f.get_trunc_effi(8, 7, state[0], state[1], len(nums), d_HR)
        print(np.mean(effi_trunc))
        effi_cos = f.get_effi_cos(t, len(nums))
        f.get_average_mni(effi_cos, len(nums))
        light_effi = f.get_light_effi(loss, effi_cos, eta_at, effi_trunc, len(nums))
        f.get_average_mni(light_effi, len(nums))
        density = []
        E_field = f.new_get_heat_W(light_effi, S, c.DNI, len(nums), loss, nums)
        for i in range(12):
            E = 0
            for time in range(5):
                E += E_field[i][time]
            density.append(E/(36*len(x)*5))
        print(density)
        return np.mean(E_field), len(nums)
    # 获取下一个状态
    # 如果符合约束条件，则返回新状态，否则放回原状态
    def get_next_state(self, state):
        next_state = state.copy()
        select = random.randint(0, 4)  # 生成0到4之间的随机整数
        flag = random.choice([-1, 1])
        vary_list = [4, 4, 2, 100, 10]
        vary_num = random.random()
        next_state[select] = next_state[select] + flag * vary_list[select] * vary_num
        next_state[3] = int(next_state[3])      # 数目只能是整数
        power, num = self.get_power(next_state)
        if self.is_valid(next_state, power, num):
            return next_state
        else:
            return state

    # 判断是否满足约束条件
    def is_valid(self, state, power, num):
        return 2 <= state[0] <= 8 and 2 <= state[1] <= 8 and 2 <= state[2] <= 6 and state[1] <= 2 * state[2] \
               and state[0] >= state[1] and state[4] >= state[0] + 5 and power >= 60000 and 100 + (num-1) * state[4] <= 350

    # 降温函数
    def update_temp(self, temp, num):
        # return self.start_temp / math.log(1+num)
        return 0.9 * temp

    # 模拟退火算法
    def simulated_annealing(self, initial_state):
        state = initial_state
        best_power_density, best_state = self.get_power_density(initial_state), initial_state
        cur_temp = self.start_temp
        change_num = 0

        # 外层循环，这里的终止条件采取的是设置温度的阈值
        while cur_temp > self.end_temp:
            # 内层循环，这里的终止条件采取的是设置内循环数目
            loop_num = self.num
            for i in range(loop_num):
                new_state = self.get_next_state(state)
                if state == new_state:
                    continue
                initial_power_density, new_power_density = self.get_power_density(state), self.get_power_density(new_state)
                delta = initial_power_density - new_power_density
                rec_p = math.exp(-(delta / cur_temp))
                if rec_p > random.uniform(0, 1):
                    state = new_state

                # print("state:", state, "pd:", new_power_density)

                # 判断是否更优
                if new_power_density > best_power_density:
                    best_state = new_state
                    best_power_density = new_power_density
                    # print("temp:", cur_temp, ", best_power_density:", best_power_density, "state:", best_state)
            # print("temp:", cur_temp, ", best_power_density:", best_power_density, "state:", best_state)
            # 更新温度

            cur_temp = self.update_temp(cur_temp, change_num)
        return best_power_density, best_state


if __name__ == '__main__':

    start_time = time.time()

    initial_state = [5.5, 5.5, 4, 2900, 10.5]

    # 模拟退火
    test = sa_tsp(10, 0.1, 100)
    best_power_density, best_state = test.simulated_annealing(initial_state)

    end_time = time.time()

    # print(best_power_density)
    # print(best_state)
    print(test.get_power_and_print(best_state))
    print("time:", end_time-start_time, "s")




q2.py
Azreal7/2023mcm_a
import function as f
import constance as c
import numpy as np
import math
import time
import random
from matplotlib import pyplot as plt
import pandas as pd


class sa_tsp:
    # 接受的输入参数：起始温度start_temp, 结束温度end_temp, 内循环次数num
    def __init__(self, s_t, e_t, num):
        self.start_temp = s_t
        self.end_temp = e_t
        self.num = num

    # 获取当前状态的功率密度
    def get_power_density(self, state):
        power, num = self.get_power(state)
        area = 0
        for i in range(len(num)):
            area += state[1][i] * state[2][i] * num[i]
        return power / area

    # 获取当前状态的功率值和摆放的圈数数组
    # state:数量 长度 宽度 高度 距离
    def get_power(self, state):
        nums = f.q3_get_new_nums(state[4], state[0], state[2])
        S = np.zeros(c.max_circle)
        for i in range(c.max_circle):
            S[i] = state[1][i] * state[2][i]
        loss, t = f.q3_get_shadow_loss(state[4], state[3], state[1], state[2], state[3])
        d_HR, eta_at = f.q3_get_d_HR(state[3], state[4])
        effi_trunc = f.q3_get_trunc_effi(8, 7, state[1], state[2], d_HR)
        effi_cos = f.get_effi_cos(t, c.max_circle)
        light_effi = f.get_light_effi(loss, effi_cos, eta_at, effi_trunc, c.max_circle)
        E_field = f.q3_get_heat_W(light_effi, S, c.DNI, c.max_circle, loss, nums)
        return np.mean(E_field), nums


    # 获取下一个状态
    # 如果符合约束条件，则返回新状态，否则放回原状态
    def get_next_state(self, state):
        next_state = [state[0]]
        for i in range(1, 5):
            next_state.append(state[i].copy())
        vary_list = [100, 4, 4, 2, 10]
        flag = random.choice([-1, 1])
        vary_num = random.random()
        select_num = random.random()
        if select_num >= 0.5:
            next_state[0] = next_state[0] + flag * vary_list[i] * vary_num
        row_select = random.randint(0, c.max_circle - 1)
        flag = random.choice([-1, 1])
        vary_num = random.random()
        select_num = random.random()
        c_select = random.randint(1, 4)
        next_state[c_select][row_select] = next_state[c_select][row_select] + flag * vary_list[c_select] * vary_num
        next_state[0] = int(next_state[0])      # 数目只能是整数
        power, num = self.get_power(next_state)
        if self.is_valid(next_state, power, num):
            return next_state
        else:
            return state

    # 判断是否满足约束条件
    def is_valid(self, state, power, num):
        valid = True
        for i in range(len(num)):
            if num[i] == 0:
                break
            valid = valid and 2 <= state[1][i] <= 8 and 2 <= state[2][i] <= 8 and 2 <= state[3][i] <= 6 \
                    and state[1][i] <= state[3][i] * 2 and state[1][i] <= state[2][i] \
                    and state[4][i] >= state[2][i] + 5
        r = 100
        for i in range(len(num)-1):
            r += state[4][i]
        return valid and r <= 350 and power >= 60000

    # 降温函数
    def update_temp(self, temp, num):
        return 0.9 * temp

    # 模拟退火算法
    def simulated_annealing(self, initial_state):
        state = initial_state
        best_power_density, best_state = self.get_power_density(initial_state), initial_state
        cur_temp = self.start_temp
        change_num = 0

        # 外层循环，这里的终止条件采取的是设置温度的阈值
        while cur_temp > self.end_temp:
            # 内层循环，这里的终止条件采取的是设置内循环数目
            loop_num = self.num
            for i in range(loop_num):
                new_state = self.get_next_state(state)
                if state == new_state:
                    continue
                initial_power_density, new_power_density = self.get_power_density(state), self.get_power_density(new_state)
                delta = initial_power_density - new_power_density
                rec_p = math.exp(-(delta / cur_temp))
                if rec_p > random.uniform(0, 1):
                    state = new_state

                # 判断是否更优
                if new_power_density > best_power_density:
                    best_state = new_state
                    best_power_density = new_power_density
                    # print("temp:", cur_temp, ", best_power_density:", best_power_density, "state:", best_state)

            print("temp:", cur_temp, ", best_power_density:", best_power_density)
            # 更新温度
            cur_temp = self.update_temp(cur_temp, change_num)
        return best_power_density, best_state


if __name__ == '__main__':

    start_time = time.time()

    initial_state = [2900, [5.4] * c.max_circle, [5.4] * c.max_circle, [4] * c.max_circle, [10.4] * c.max_circle]

    # 模拟退火
    test = sa_tsp(100, 0.1, 100)
    best_power_density, best_state = test.simulated_annealing(initial_state)

    end_time = time.time()

    print(best_power_density)
    print(best_state)
    print("time:", end_time-start_time, "s")
q3.py
Azreal7/2023mcm_a