【2023年第十三届MathorCup高校数学建模挑战赛】思路总结分析

兜兜里有好多糖

已于 2023-11-21 20:31:29 修改

阅读量994

点赞数 10

分类专栏：数学建模笔记文章标签：算法大数据数据分析模拟退火算法人工智能

于 2023-04-17 17:51:22 首次发布

本文链接：https://blog.csdn.net/Gorege__Hu/article/details/130205272

版权

数学建模笔记专栏收录该内容

10 篇文章 32 订阅

订阅专栏

【写在前面的话】
我们选择A题，分析A题题目可以得知属于一种组合优化模型，类似于旅行商问题，0-1背包问题等等。该类问题通常采用遗传算法，粒子群算法，模拟退火算法等算法进行求解。由于本题需要我们建立出数学模型之后通过转换为QUBO模型，从而建立量子退火模型，从而可以实现在量子计算机中求解。所以我们直接选择采用模拟退火模型。

第一问太简单，我TM直接穷举🤣，代码给你们看看（别太荒谬哈哈）：

clc;clear;
load bank.mat
max=0;
i_max=0;
j_max=0;
for j=1:2:200
    for i=1:10
        banifit=1000000*bank(i,j)*0.08*(1-bank(i,j+1))-1000000*bank(i,j)*bank(i,j+1);
        if(banifit>max)
            max=banifit;
            i_max=i;
            j_max=j;
        end
    end
end

开玩笑的，正儿八经还得用模拟退火模型求解嘿嘿
第一问的代码（这么详细的注释，我哭死，你怎么还不哭😡）：

clc;
clear;
tic
load bank.mat
x1=bank(:,1:2:end);
x1=x1(:);
x1=x1';
x2=bank(:,2:2:end);
x2=x2(:);
x2=x2';
n=length(x1);
%% 参数初始化
T0 = 1000;   % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 500;  % 最大迭代次数
Lk = 20;  % 每个温度下的迭代次数
alfa = 0.95;  % 温度衰减系数

%%  生成初始解
% 初始解中每个元素都为0，表示一个阈值都不选
way0 = zeros(1,n); 
load_matrix=way0;
% 计算相应的银行收益率
profit0=sum(1000000*0.08*way0.*x1.*(1-x2)-1000000*way0.*x1.*x2);

%% 记录一些中间过程的参数，方便最后画图
max_profit = profit0;     % 初始化找到的最佳的解对应的利润为profit0
PROFIT = zeros(maxgen,1); % 记录每一次外层循环结束后找到的max_profit (方便画图）

%% 模拟退火过程
for iter = 1 : maxgen  % 外循环, 我这里采用的是指定最大迭代次数
    for i = 1 : Lk  %  内循环，在每个温度下开始迭代       
        [way1,load_matrix] = gen_new_way(way0,load_matrix,n);  % 调用我们自己写的gen_new_way函数生成新的阈值方案
        profit1 =sum(1000000*0.08*way1.*x1.*(1-x2)-1000000*way1.*x1.*x2); % 计算新阈值方案的利润
        if profit1 > profit0 
            way0 = way1; % 更新当前阈值方案为新的阈值方案
            profit0 = profit1;
        else
            p = exp(-(profit0 - profit1)/T); % 根据Metropolis准则计算一个概率
            if rand(1) < p   % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
                way0 = way1;  % 更新当前阈值方案为新阈值方案
                profit0 = profit1;
            end
        end
        % 判断是否要更新找到的最佳的解
        if profit0 > max_profit  % 如果当前解更好，则对其进行更新
            max_profit = profit0;  % 更新最大的利润
            best_way = way0;  % 更新找到的最好的阈值方案
        end
    end
    PROFIT(iter) = max_profit; % 保存本轮外循环结束后找到的最大的利润
    T = alfa*T;   % 温度下降     
end

%%结果
disp('最佳的阈值方案是：'); disp(best_way)
disp('此时最优值是：'); disp(max_profit)
plot(1:maxgen,PROFIT,'r-');
xlabel('迭代次数');
ylabel('银行最大收益值');
hold on
toc

主函数中用到的gen_new_way函数：

function [way1, delta_weight] = gen_new_way(way0, n, weight)
% way0：原来的装载方案；  n是物品个数  ; weight是物品的重量
% way1：新的装载方案； delta_weight：重量的变化量
       i = randi([1,n],1);  % 随机选取一个物品i
       if way0(i) == 1  % 如果物品i在原来的背包中
           way0(i) = 0;  % 从背包里取出物品i
           empty_ind = find(way0==0); % 看看背包中哪几个位置是空的
           length_empth_ind = length(empty_ind);  % 计算空位置的个数
           ind = randi([1, length_empth_ind], 1); % 在这几个空位置中取一个随机下标 
           j = empty_ind(ind);  % j就是随机选择出来的背包中的某个空位置
           way0(j) = 1;  % 装入物品j
           delta_weight = weight(j) - weight(i);  % 计算重量的变化量
       else  % 如果物品i原来就不在背包中
           way0(i) = 1; % 在背包中放入物品i
           delta_weight = weight(i);  % 计算重量的变化量
           if rand(1) < 0.5 % 判断是否需要再取出一个物品（这里的0.5是需要再取出物品的概率，可以自己调整）
               fill_ind = find(way0==1); % 看看背包中哪几个位置有物品
               length_fill_ind = length(fill_ind);  % 计算有物品的位置的个数
               ind = randi([1, length_fill_ind], 1); % 在这几个有物品的位置中取一个随机下标 
               j = fill_ind(ind);  % j就是随机选择出来的背包中的某个有物品的位置
               way0(j) = 0;  % 取出物品j
               delta_weight = weight(i) - weight(j);  % 计算重量的变化量
           end
       end
       way1 = way0;  % 将调整后的way0赋值给way1
end

第二问就是第一问的一个拓展，主函数差不多，要变得是收益函数的形式和新解的产生方法（正经起来了不是？）：

clc;
clear;
tic
load bank.mat
x1=bank(:,1:2:6);
x1=x1(:);
x1=x1';
x2=bank(:,2:2:6);
x2=x2(:);
x2=x2';
n=length(x1);
%% 参数初始化
T0 = 1000;   % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 50;  % 最大迭代次数
Lk = 20;  % 每个温度下的迭代次数
alfa = 0.95;  % 温度衰减系数

%%  生成初始解
% 初始解中每个元素都为0，表示一个阈值都不选
way0 = zeros(1,n);
% 随机选取三个阈值
i_1= randi([1,10],1); 
i_2= randi([11,20],1); 
i_3= randi([21,30],1); 
way0(i_1)=1;
way0(i_2)=1;
way0(i_3)=1;
% 计算三个阈值的组合通过率和组合坏账率
T_tresh_ind=find((way0)~=0);
T_tresh=x1(T_tresh_ind(1))*x1(T_tresh_ind(2))*x1(T_tresh_ind(3));
H_tresh_ind=find((way0)~=0);
H_tresh=1/3*(x2(H_tresh_ind(1))+x2(H_tresh_ind(2))+x2(H_tresh_ind(3)));
% 计算相应的银行收益率
profit0=sum(1000000*0.08*T_tresh*(1-H_tresh)-1000000*T_tresh*H_tresh);

%% 记录一些中间过程的参数，方便最后画图
max_profit = profit0;     % 初始化找到的最佳的解对应的利润为profit0
PROFIT = zeros(maxgen,1); % 记录每一次外层循环结束后找到的max_profit 

%% 模拟退火过程
for iter = 1 : maxgen  % 外循环, 我这里采用的是指定最大迭代次数
    for i = 1 : Lk  %  内循环，在每个温度下开始迭代       
        [way1] = gen_new_way_2(way0,n);  % 调用我们自己写的gen_new_way函数生成新的装载方案
        T_tresh_ind=find((way1)~=0);
        T_tresh=x1(T_tresh_ind(1))*x1(T_tresh_ind(2))*x1(T_tresh_ind(3));
        H_tresh_ind=find((way1)~=0);
        H_tresh=1/3*(x2(H_tresh_ind(1))+x2(H_tresh_ind(2))+x2(H_tresh_ind(3)));
        profit1 =sum(1000000*0.08*T_tresh*(1-H_tresh)-1000000*T_tresh*H_tresh);
        if profit1 > profit0  
            way0 = way1; % 更新当前阈值方案为新的阈值方案
            profit0 = profit1;
        else
            p = exp(-(profit0 - profit1)/T); % 根据Metropolis准则计算一个概率
            if rand(1) < p   % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
                way0 = way1;  % 更新当前阈值方案为新的阈值方案
                profit0 = profit1;
            end
        end
        % 判断是否要更新找到的最佳的解
        if profit0 > max_profit  % 如果当前解更好，则对其进行更新
            max_profit = profit0;  % 更新最大的利润
            best_way = way0;  % 更新找到的最好的阈值方案
        end
    end
    PROFIT(iter) = max_profit; % 保存本轮外循环结束后找到的最大的利润
    T = alfa*T;   % 温度下降     
end

%% 结果
disp('最佳的阈值方案是：'); disp(best_way)
disp('此时最优值是：'); disp(max_profit)
plot(1:maxgen,PROFIT,'r-');
xlabel('迭代次数');
ylabel('银行最大收益值');
hold on
toc

主函数中用到的gen_new_way_2函数：

function [way1] = gen_new_way_2(way0,n)
% way0：原来的阈值方案；  
% way1：新的阈值方案；
       way0_split=reshape(way0,10,[]);
       way0_1=way0_split(:,1);
       way0_2=way0_split(:,2);
       way0_3=way0_split(:,3);
       i=randi([1,n],1);
       if i<=10
           if way0(i) == 1  % 如果阈值i在原来的选择范围中
               empty_ind = find(way0_1==0); % 查看哪几个阈值是没有被选中的
               length_empth_ind = length(empty_ind);  % 计算空位置的个数
               ind = randi([1, length_empth_ind], 1); % 在这几个空位置中取一个随机下标 
               j = empty_ind(ind);  % j就是随机选择出来的阈值中的某个空位置
               way0(j) = 1;  % 将该阈值纳入考虑范畴
               way0(i) = 0;  % 从已选择的阈值中删除第i个阈值
           else
               empty_ind = find(way0_1==1);
               j = empty_ind(1);
               way0(j)=0;
               way0(i)=1;
           end
       elseif i<=20
           if way0(i) == 1 
               empty_ind = find(way0_2==0); 
               length_empth_ind = length(empty_ind); 
               ind = randi([1, length_empth_ind], 1); 
               j = empty_ind(ind); 
               way0(j+10) = 1; 
               way0(i) = 0; 
           else
               empty_ind = find(way0_2==1);
               j = empty_ind(1);
               way0(j+10)=0;
               way0(i)=1;
           end
       else
           if way0(i) == 1  
               empty_ind = find(way0_3==0);
               length_empth_ind = length(empty_ind); 
               ind = randi([1, length_empth_ind], 1);
               j = empty_ind(ind); 
               way0(j+20) = 1; 
               way0(i) = 0; 
           else
               empty_ind = find(way0_3==1);
               j = empty_ind(1);
               way0(j+20)=0;
               way0(i)=1;
           end
       end
       way1 = way0;  % 将调整后的阈值选择赋值给way1
end

其实我感觉我第二问的产生新解的函数写的有点复杂了，有点冗余，第三问虽然问题更复杂了，涉及的数据更多了，但是代码又简化了很多嘿嘿：

clc;
clear;
tic
load bank.mat
x1=bank(:,1:2:end);
x1=x1(:);
x1=x1';
x2=bank(:,2:2:end);
x2=x2(:);
x2=x2';
%% 参数初始化
n=length(x1);
T0 = 1000;   % 初始温度
T = T0; % 迭代中温度发生改变，第一次迭代时温度T0
maxgen = 500;  % 最大迭代次数
Lk = 20;  % 每个温度下的迭代次数
alfa = 0.96;  % 温度衰减系数

%%  生成初始解
% 初始解中每个元素都为0，表示一个阈值都不选取
way0 = zeros(1,n); 
% 随机选取三个阈值 
i_1= randi([1,10],1);  
i_2= randi([11,20],1); 
i_3= randi([21,30],1); 
way0(i_1)=1;
way0(i_2)=1;
way0(i_3)=1;
% 计算三个阈值的组合通过率和组合坏账率
T_tresh_ind=find((way0)~=0);
T_tresh=x1(T_tresh_ind(1))*x1(T_tresh_ind(2))*x1(T_tresh_ind(3));
H_tresh_ind=find((way0)~=0);
H_tresh=1/3*(x2(H_tresh_ind(1))+x2(H_tresh_ind(2))+x2(H_tresh_ind(3)));
% 计算相应的银行收益率
profit0=sum(1000000*0.08*T_tresh*(1-H_tresh)-1000000*T_tresh*H_tresh);

%% 记录一些中间过程的参数，方便最后画图
max_profit = profit0;     % 初始化找到的最佳的解对应的利润为profit0
PROFIT = zeros(maxgen,1); % 记录每一次外层循环结束后找到的max_profit 

%% 模拟退火过程
for iter = 1 : maxgen  % 外循环, 我这里采用的是指定最大迭代次数
    for i = 1 : Lk  %  内循环，在每个温度下开始迭代       
        [way1] = gen_new_way_3(way0,n);  % 调用我们自己写的gen_new_way函数生成新的阈值方案
        % 计算三个阈值的组合通过率和组合坏账率和相应的银行收益率
        T_tresh_ind=find((way1)~=0);
        T_tresh=x1(T_tresh_ind(1))*x1(T_tresh_ind(2))*x1(T_tresh_ind(3));
        H_tresh_ind=find((way1)~=0);
        H_tresh=1/3*(x2(H_tresh_ind(1))+x2(H_tresh_ind(2))+x2(H_tresh_ind(3)));
        profit1 =sum(1000000*0.08*T_tresh*(1-H_tresh)-1000000*T_tresh*H_tresh);
        %权衡是否采用新的阈值方案
        if profit1 > profit0  % 如果没有超重,而且新装载方案的利润更高
            way0 = way1; % 更新当前装载方案为新装载方案
            profit0 = profit1;
        else
            p = exp(-(profit0 - profit1)/T); % 根据Metropolis准则计算一个概率
            if rand(1) < p   % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
                way0 = way1;  % 更新当前装载方案为新装载方案
                profit0 = profit1;
            end
        end
        % 判断是否要更新找到的最佳解
        if profit0 > max_profit  % 如果当前解更好，则对其进行更新
            max_profit = profit0;  % 更新最大的利润
            best_way = way0;  % 更新找到的最好的装载方案
        end
    end
    PROFIT(iter) = max_profit; % 保存本轮外循环结束后找到的最大的利润
    T = alfa*T;   % 温度下降     
end

%% 结果
disp('最佳的阈值组合是：'); disp(best_way)
disp('此时最优值是：'); disp(max_profit)
plot(1:maxgen,PROFIT,'r-');
xlabel('迭代次数');
ylabel('银行最大收益值');
hold on
toc

主函数中用到的gen_new_way_3函数：

function [way1] = gen_new_way_3(way0,n)
       empty_ind = find(way0==1); % 查找目前的阈值选择的是哪几个，生成其下标索引
       length_empth_ind = length(empty_ind);  
       ind = randi([1, length_empth_ind], 1); 
       j = empty_ind(ind);%随机选取目前下标中的一个
       way0(j) = 0;%将该阈值剔除出我们的选择范围
       i=randi([1,n],1);%随机选择一个新的阈值
       while(panduan(i,empty_ind)==1)%判断该新阈值是否与原阈值属于同一个银行卡
            i=randi([1,n],1);%如果已有，则重新生成，保证一个银行卡只选择一个阈值
       end
       way0(i) = 1; %确保该新阈值是有效阈值之后，将其纳入选择范围
       way1 = way0;  % 将调整后的阈值选择赋值给way1
end

其中用到的panduan函数：

function [result] = panduan(i,index)
%输入随机生成的阈值下标和参数为1的索引向量
%判断新生成的阈值下标是否处于原阈值下标属于同一个银行卡
%如果属于同一个银行卡，则返回1，如果不属于，则返回0
    if i>=index(1)-mod(index(1),10)+1 && i<=index(1)-mod(index(1),10)+10
        result=1;
    elseif i>=index(2)-mod(index(2),10)+1 && i<=index(2)-mod(index(2),10)+10
        result=1;
    elseif i>=index(3)-mod(index(3),10)+1 && i<=index(3)-mod(index(3),10)+10
        result=1;
    else
        result=0;
    end
end

第二三问运行结果每次可能都会有些许不同，因为并不是只有一个最优解，所以他会在几个最优解之间跳动，这是正常现象。
还不点赞？！怎么回事？！

兜兜里有好多糖

关注

10
点赞
踩
24

收藏

觉得还不错? 一键收藏
打赏
4
评论
【2023年第十三届MathorCup高校数学建模挑战赛】思路总结分析

我们选择A题，分析A题题目可以得知属于一种组合优化模型，类似于旅行商问题，0-1背包问题等等。该类问题通常采用遗传算法，粒子群算法，模拟退火算法等算法进行求解。由于本题需要我们建立出数学模型之后通过转换为QUBO模型，从而建立量子退火模型，从而可以实现在量子计算机中求解。第二三问运行结果每次可能都会有些许不同，因为并不是只有一个最优解，所以他会在几个最优解之间跳动，这是正常现象。第二问就是第一问的一个拓展，主函数差不多，要变得是收益函数的形式和新解的产生方法（正经起来了不是？
复制链接

扫一扫