集训中用到的代码记录，方便后续再次调用

模拟退火（不透明物体染色）

MATLAB代码如下,这是最基础的模拟算法，可以根据实际修改的地方有更新解的方式、循环停止的方式、添加约束条件

%% SA 模拟退火: 
% clear; clc
%% 参数初始化
narvs = 3; % 变量个数
T0 = 100;   % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 50 ;  % 每个温度下的迭代次数
alfa = 0.95;  % 温度衰减系数
T_min = 0.1;
x_lb = 0.00;    %所求解的解的下界
x_ub = 0.5;      %解的上界

% data1=xlsread("C:/Users/asus/Desktop/集训模型5/附件5.xlsx",1,'B53:D62');  %读入数据表格 

%  随机生成一个初始解
 x0 = zeros(1,narvs);
 for i = 1: narvs
     x0(i) = x_lb + (x_ub-x_lb)*rand(1);    
 end

y0 = Obj_fun2(x0,u,data1,data2,data3); % 计算当前解的函数值，自己设定函数值
min_y = y0;     % 初始化找到的最佳的解对应的函数值为y0
best_x = x0;

%% 模拟退火过程
while(T>T_min)
for iter = 1 : maxgen  % 内循环, 我这里采用的是指定最大迭代次数

%%%%%%%%%%产生新解
        y = randn(1,narvs);  % 生成1行narvs列的N(0,1)随机数
        z = y / sqrt(sum(y.^2)); % 根据新解的产生规则计算z
        x_new = x0 + z*T; % 根据新解的产生规则计算x_new的值
        % 如果这个新解的位置超出了定义域，就对其进行调整
        for j = 1: narvs
            if x_new(j) < x_lb
                r = rand(1);
                x_new(j) = r*x_lb+(1-r)*x0(j);
            elseif x_new(j) > x_ub
                r = rand(1);
                x_new(j) = r*x_ub+(1-r)*x0(j);
            end
        end
        x1 = x_new;    % 将调整后的x_new赋值给新解x1
        y1 = Obj_fun2(x1,u,data1,data2,data3);  % 计算新解的函数值
        w1 =60*2*x1(1)+65*x1(3)+63*2*x1(3);
%                    disp("w1=");disp(w1)

        if y1 < y0    % 如果新解函数值小于当前解的函数值
            x0 = x1; % 更新当前解为新解
%             disp('y1<y0更新当前解为新解')
            y0 = y1;
 
        else
            p = exp(-( y1-y0)/T); % 根据Metropolis准则计算一个概率
            rrr = rand(1);
            if  rrr< p   % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
                x0 = x1;  % 更新当前解为新解
%                 disp('y1>y0更新当前解为新解')
                y0 = y1;
            end
        end
        % 判断是否要更新找到的最佳的解
        if  y0<min_y  % 如果当前解更好，则对其进行更新
            best_x = x0;  % 更新找到的最好的x
            min_y = y0;
%             disp('更新找到的最佳的解!!!!!')
        end
      end
    T = alfa*T;   % 温度下降
    sst = strcat('当前温度为',num2str(T));
    disp(sst)
    ss = strcat('最小的成本为',num2str(cost_min),'其色差为',num2str(min_y));
    disp(ss)  
end
disp('最佳的位置是：'); disp(best_x)
disp('此时对应色差是：'); disp(min_y)

遗传算法（火车作为分配）

python代码如下

import random

#数据预处理
import numpy as np
import pandas as pd
###读取xlsx文件
file_path = "C:/Users/asus/Desktop/集训模型2/数据.xlsx"     #要读取的文件的名称
file_dict = pd.read_excel(file_path,sheet_name=None)   #读取文件
def yuchuli(file_dict):
    Old = []
    V = []
    All = []
    pianhao = []
    feipianhao = []
    V_num = []
    for sheet_name, df in file_dict.items():      #对文件内容进行处理
        if sheet_name =="问题一数据":      #对某一页的数据进行处理
            df.columns = range(df.shape[1])       #去掉表头
            cols = df.iloc[:, :]    #读取数据列表的某一行，某一列
            cols = cols.values
            for i in range(0,len(cols)):
                bianhao = i + 1
                if cols[i][1] == '无偏好':
                    feipianhao.append(i + 1)
                    bed_type = "无偏好"
                    position1 = '无偏好'
                    position2 = '无偏好'
                else:

                    bed_type = cols[i][1]
                    position1 = cols[i][2]
                    position2 = cols[i][3]
                is_elderly = cols[i][4]
                if is_elderly == '是':
                    Old.append((bianhao,bed_type, position1, position2, is_elderly))
                else:
                    V_num.append(i+1)
                    if bed_type != '无偏好':
                        pianhao.append(i + 1)
                    V.append((bianhao,bed_type, position1, position2, is_elderly))
                All.append((bianhao,bed_type, position1, position2,is_elderly))
    return Old, V, All,pianhao,feipianhao,V_num

def is_number_in_2d_collection(collection, number):
    for row in collection:
        for item in row:
            if int(item) == number:
                return True
    return False

#创建初始种群
def create_answer(Old,V,n,g,baojian):
    print("开始创建初始种群！")
    result = []
    result_labels = []
    is_n = 0
    # for i in range(n):
    while(is_n<n):
        #print("开始产生第{}个解".format(is_n))
        B = np.ones((baojian,6))*240
        old_labels = np.ones((baojian,6))*240
        #优先安排老人住下铺
        for o in Old:
            flag = 1
            while(flag!=0):
                k = random.randint(1,baojian)    #随机选一个包间k
                #print("k0=",k)
                if k < g+1:     #如果分配到软卧
                    l = random.sample([0,2],1)   #在下铺随机选一个
                    if B[k-1][l] == 240:
                        B[k-1][l] = o[0]
                        old_labels[k-1][l] = 1
                        flag = 0
                        #print("第{}个老人已经分配铺位".format(o[0]))
                else:     #如果被分配到硬卧
                    l = random.sample([0, 3], 1)  # 在下铺随机选一个
                    if B[k-1][l] == 240:
                        B[k-1][l] = o[0]
                        old_labels[k-1][l] = 1
                        flag = 0
                        #print("第{}个老人已经分配铺位".format(o[0]))
                if o[0] == 234 and is_number_in_2d_collection(B, o[0]):
                    break

        #再安排非老年人
        #print("再安排非老年人")
        for v in V:
            #print("非老年人开始安排铺位！")
            flag = 1
            while (flag!=0):
                k = random.randint(1, baojian)  # 随机选一个包间k
                #print("k1= ",k)
                if k < g+1:  # 如果分配到软卧
                    l = random.sample([0,1,2,3], 1)  # 在随机选一个铺位
                    if B[k-1][l] == 240:
                        #print("v k l ",v[0],k,l)
                        B[k-1][l] = v[0]
                        old_labels[k-1][l] = 0
                        flag = 0
                else:     #如果被分配到硬卧
                    l = random.sample([0,1,2,3,4,5], 1)  # 随机选一个铺位
                    if B[k-1][l] == 240:
                        B[k-1][l] = v[0]
                        old_labels[k-1][l] = 0
                        flag = 0
                if v[0] == 233 and is_number_in_2d_collection(B, v[0]):
                    break

        if is_n != 0 and np.any(np.all(result == B, axis=1)):
            continue
        else:
            result_labels.append(old_labels)
            result.append(B)
            is_n += 1
        print("is_n=",is_n)
    #返回n种解以及解得标签
    return result, result_labels

######计算适应度函数
soft_label = ['下','上','下','上']
hard_label = ['下','中','上','下','中','上']
def error_level(new_answer,All,d):
    error = []
    #print("000ALL=",All)
    ##计算每一种方案的反感度
    for answer in new_answer:    ##遍历解空间的每一种方案
        dislike = 0
        for j in range(len(answer)):          ##遍历每一个包间
            for i in range(6):    ##遍历每一个位置
                s = int(answer[j][i])
                #print("提取在包间的{}位乘客的信息".format(s))
                if s !=240:
                    #print("筛选后提取在包间的{}位乘客的信息".format(s))
                    #print("All=",All)
                    people = All[s-1]    ##提取在包间的i位乘客的信息
                    if people[1] != '无偏好':   ##分析有偏好的人的反感度
                        if j < d:
                            if people[1] != '软卧':
                                dislike += 2
                            if people[2] != soft_label[i] and people[3]!=soft_label[i]:
                                dislike += 1
                        else:
                            if people[1] != '硬卧':
                                dislike +=2
                            if people[2] != hard_label[i] and people[3] != hard_label[i]:
                                dislike += 1
        error.append(dislike)
    return error

#######设计一个根据适应度来选择优胜者的函数
def choice_selected(old_answer,All,n):

    result = []
    error_c = error_level(old_answer, All)   ##所有解的反感度的列表
    error_cc = error_c
    error_c.sort()
    for i in range(n):
        min_index = error_cc.index(error_c[i])
        result.append(old_answer[min_index])
    return result
#####重新写一个基因重组的算法
######基因重组
def gene_recombination(old_answer,n,pianhao,peo):
    # 将数组转化为一维的
    for i in range(len(old_answer)):
        array_2d = np.array(old_answer[i])
        old_answer[i] = np.reshape(array_2d, -1)
        #print("第{}个解为{}".format(i, old_answer[i]))
    #print("基因重组前解的个数=",len(old_answer))
    for i in range(n//2):
        rand = random.uniform(0,1)
        new_answer1 = old_answer[i]
        new_answer2 = old_answer[i + n // 2]
        old_answer.append(new_answer1)
        old_answer.append(new_answer2)
        if rand < peo:
            x = random.sample(pianhao,1)
            index_x1 = list(old_answer[i]).index(x[0])
            index_x2 = list(old_answer[i+n//2]).index(x[0])
            old_answer[i][index_x1] = old_answer[i][index_x2]
            old_answer[i][index_x2] = x[0]
            old_answer[i+n//2][index_x2] = old_answer[i+n//2][index_x1]
            old_answer[i + n // 2][index_x1] = x[0]

    # 将数组恢复为二维的
    #print("基因重组后解的个数=", len(old_answer))
    for i in range(len(old_answer)):
        #print("基因重组后解的个数=",len(old_answer))
        #print("第{}个解为{}".format(i, old_answer[i]))
        flat_array = np.array(old_answer[i])
        old_answer[i] = np.reshape(flat_array, (42, 6))
    #print("进行基因重组后的种群=",old_answer)
    return old_answer
#######基因变异
def variation(old_answer,pianhao,feipianhao,pro,V_num):
    #将数组转化为一维的
    #print("基因变异前解的空间数量=",len(old_answer))
    for i in range(len(old_answer)):
        array_2d = np.array(old_answer[i])
        old_answer[i] = np.reshape(array_2d,-1)
        #print("第{}个解为{}".format(i,old_answer[i]))
    for i in range(len(old_answer)):
        new_answer = old_answer[i]
        old_answer.append(new_answer)
        #print("种群规模：",len(old_answer))
        rand = random.uniform(0, 1)
        if rand < pro:    #如果概率小于变异概率则发生变异
            bianyi = random.sample(V_num,48)
            # rand1 = random.uniform(0,1)
            # if rand1 < 0.95:
            #     rand1 = 0.95
            # bianyi1 = random.sample(pianhao, round(48*rand1))
            # bianyi2 = random.sample(feipianhao,48-round(48*rand1))
            # bianyi = bianyi1+bianyi2
            indexjj = {}
            for j in bianyi:
                #print("i={},j={}".format(i,j))
                #print("old_answer=",old_answer[i])
                indexj = list(old_answer[i]).index(j)
                #print('old_answer[i]=',len(old_answer[i]))
                indexjj[j] = indexj
            #print("第{}个解需进行变异，变异前乘客号为{}".format(i,indexjj))
            values = list(indexjj.values())
            random.shuffle(values)
            indexjj = {k:v for k,v in zip(indexjj.keys(),values)}
            #new_dict = {k: v for k, v in zip(my_dict.keys(), values)}
            #print("第{}个解需进行变异，变异后的乘客号为{}".format(i, indexjj))
            for v in values:    ####把打乱后的数重新插入解中
                for key, value in indexjj.items():
                    if value == v:
                        #print("键:", key)
                        #key = key
                        old_answer[i][v] = key
                        break
    #将数组恢复为二维的
    for i in range(len(old_answer)):
        #print("基因变异后解的空间数量=", len(old_answer))
        flat_array = np.array(old_answer[i])
        old_answer[i] = np.reshape(flat_array,(42,6))
    return old_answer
#
Old, V, All_pasenger ,pianhao,feipianhao,V_num= yuchuli(file_dict)

print("老人合集=",len(Old))
print("非老人合集=",V)
print("所有乘客合集=",All_pasenger)
n = 20 #设置种群规模
middle_answer, middle_answer_labels = create_answer(Old,V,n,9,42)
error9 = error_level(middle_answer,All_pasenger,9)
error9 = [x - 120 for x in error9]
print("软卧车厢为9时反感度为：",error9)
middle_answer_21, middle_answer_labels_21 = create_answer(Old,V,n,21,46)
error21 = error_level(middle_answer_21,All_pasenger,21)
error21 = [x - 160 for x in error21]
print("软卧车厢为21时反感度为：",error21)

import matplotlib.pyplot as plt

# 创建数据
x = [i for i in range(1,21)]
# 绘制曲线1
plt.plot(x, error9, label='soft=36')

# 绘制曲线2
plt.plot(x, error21, label='soft=84')

# 设置图例
plt.legend()

# 设置图标标题和轴标签
plt.title("dislike")
plt.xlabel("ways")
plt.ylabel("dislike")
# 显示图形
plt.show()

一些常用的常规操作

MATLAB求导

syms lada beta
% % 定义表达式
CP=0.5176.*((1./(lada+0.08.*beta)-0.035./(beta.^3+1)).*116-0.4.*beta-5)...
.*exp(-1.*(1./(lada+0.08.*beta)-0.035./(beta.^3+1)).*21)+0.068;
CP=(0.44-0.0167*beta)*sin((pi*(lada-3))/(15-0.3*beta))-0.00184*(lada-3)*beta;
% 求 x 的偏导数
dy_dlada = diff(CP, lada);
% 求 z 的偏导数
dy_dbeta = diff(CP, beta);
% 输出结果
disp(['偏导数 dy/dlada = ' char(dy_dlada)]);
disp(['偏导数 dy/dbeta = ' char(dy_dbeta)]);

MATLAB求微分方程

dsolve('DT=a*T-b*T*ln(T)','t')
dsolve('DT=a*T*ln(K/T)-(m*v+2*n*v^2)*T','t')

 MATLAB求积分

syms x a
y=int(a/(sqrt(1+x^2)),x)