数学建模暑期集训17：蒙特卡洛法
Posted Z|Star
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了数学建模暑期集训17：蒙特卡洛法相关的知识，希望对你有一定的参考价值。
定义
蒙特卡罗⽅法⼜称统计模拟法，是⼀种随机模拟⽅法，将所求解的问题同⼀定的概率模型相联系，⽤电⼦计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这⼀⽅法的概率统计特征，故借⽤赌城蒙特卡罗命名。
起源
1777年，法国Buffon提出⽤投针实验的⽅法求圆周率，这被认为是蒙特卡罗⽅法的起源。
原理
由⼤数定理可知，当样本容量⾜够⼤时，事件的发⽣频率即为其概率。
蒙特卡洛是一种思想，不是一种具体的算法。
下面将用例题的方式来介绍蒙特卡洛如何在例题中运用。
matlab函数知识补充
本篇博文所用到的matlab的一些函数如下表所示：
| 函数使用示例 | 具体含义 |
|---|---|
| randi([1,5],5,8) | 在区间[1,5]内随机取出大小为5*8的整数矩阵 |
| randi([1,5]) | 在区间[1,5]内随机取出1个整数 |
| normrnd(10,2) | 均值为10 标准差为2（方差为4）的正态分布随机数 |
| exprnd(5) | 均值为5的指数分布随机数（对应的参数为0.2） |
| mean([1,2,3]) | 获得1，2，3的均值 |
| tic …toc | tic函数和toc函数可以用来返回代码运行的时间 |
| format long g | 将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法） |
| unifrnd(0,5,4,3) | 输出在[0,5]之间均匀分布的随机数组成的4行3列的矩阵 |
| plot([1,2],[5,10],’-o’) | 画出一条线段，x范围是[1, 2] ，y范围是[5,10] |
| randperm(5) | 生成1-5组成的一个随机序列(类似于洗牌的操作) |
例题
三门问题
问题背景
你参加⼀档电视节⽬，节⽬组提供了ABC三扇⻔，主持⼈告诉你，其中⼀扇⻔后边有辆汽⻋，其它两扇⻔后是空的。假如你选择了B⻔，这时，主持⼈打开了C⻔，让你看到C⻔后什么都没有，然后问你要不要改选A⻔？
matlab求解
n = 100000; % n代表蒙特卡罗模拟重复次数
a = 0; % a表示不改变主意时能赢得汽车的次数
b = 0; % b表示改变主意时能赢得汽车的次数
for i= 1 : n % 开始模拟n次
x = randi([1,3]); % 随机生成一个1-3之间的整数x表示汽车出现在第x扇门后
y = randi([1,3]); % 随机生成一个1-3之间的整数y表示自己选的门
% 下面分为两种情况讨论：x=y和x~=y
if x == y % 如果x和y相同，那么我们只有不改变主意时才能赢
a = a + 1; b = b + 0;
else % x ~= y ，如果x和y不同，那么我们只有改变主意时才能赢
a = a + 0; b = b +1;
end
end
disp(['蒙特卡罗方法得到的不改变主意时的获奖概率为：', num2str(a/n)]);
disp(['蒙特卡罗方法得到的改变主意时的获奖概率为：', num2str(b/n)]);
求解结果
蒙特卡罗方法得到的不改变主意时的获奖概率为：0.33206
蒙特卡罗方法得到的改变主意时的获奖概率为：0.66794
模拟排队问题
问题背景
假设某银⾏⼯作时间只有⼀个服务窗⼝，⼯作⼈员只能逐个的接待顾客。当来的顾客较多时，⼀部分顾客就需要排队等待。
假设：
- 顾客到来的间隔时间服从参数为0.1的指数分布
- 每个顾客的服务时间服从均值为10，⽅差为4的正态分布(单位为分钟，若服务时间⼩于1分钟，则按1分钟计算)
- 排队按先到先服务的规则，且不限制队伍的⻓度，每天⼯作时⻓为8⼩时。
试回答下⾯的问题：
- 模拟⼀个⼯作⽇，在这个⼯作⽇共接待了多少客户，客户平均等待的时间为多少?
- 模拟100个⼯作⽇，计算出平均每⽇接待客户的个数以及每⽇客户的平均等待时⻓。
模型建立
matlab求解
变量说明：
% x(i)表示第i-1个客户和第i个客户到达的间隔时间，服从参数为0.1的指数分布
% y(i)表示第i个客户的服务持续时间，服从均值为10方差为4(标准差为2)的正态分布 (若小于1则按1计算)
% c(i)表示第i个客户的到达时间，那么c(i) = c(i-1) + x(i)，初始值c0=0
% b(i)表示第i个客户开始服务的时间
% e(i)表示第i个客户结束服务的时间，初始值e0=0
% 第i个客户结束服务的时间 = 第i个客户开始服务的时间 + 第i个客户的服务持续时间
% 即：e(i) = b(i) + y(i）
% 第i个客户开始服务的时间取决于该客户的到达时间和上一个客户结束服务的时间
% 即：b(i) = max(c(i),e(i-1))，初始值b1=c1;
% 第i个客户等待的时间 = 第i个客户开始服务的时间 - 第i个客户到达银行的时间
% 即：wait(i) = b(i) - c(i)
% w表示所有客户等待时间的总和
% 假设一天内银行最终服务了n个顾客，那么客户的平均等待时间t = w/n
问题一：
clear
tic %计算tic和toc中间部分的代码的运行时间
i = 1; % i表示第i个客户，最开始取i=1
w = 0; % w用来表示所有客户等待的总时间，初始化为0
e0 = 0; c0 = 0; % 初始化e0和c0为0
x(1) = exprnd(10); % 第0个客户(假想的)和第1个客户到达的时间间隔
c(1) = c0 + x(1); % 第1个客户到达的时间
b(1) = c(1); % 第1个客户的开始服务的时间
while b(i) <= 480 % 开始设置循环，只要第i个顾客开始服务的时间(时刻)小于480，就可以对其服务（银行每天工作8小时，折换为分钟就是480分钟）
y(i) = normrnd(10,2); % 第i个客户的服务持续时间，服从均值为10方差为4(标准差为2)的正态分布
if y(i) < 1 % 根据题目的意思：若服务持续时间不足一分钟，则按照一分钟计算
y(i) = 1;
end
e(i) = b(i) + y(i); % 第i个客户结束服务的时间 = 第i个客户开始服务的时间 + 第i个客户的服务持续时间
wait(i) = b(i) - c(i); % 第i个客户等待的时间 = 第i个客户开始服务的时间 - 第i个客户到达银行的时间
w = w + wait(i); % 更新所有客户等待的总时间
i = i + 1; % 增加一名新的客户
x(i) = exprnd(10); % 这位新客户和上一个客户到达的时间间隔
c(i) = c(i-1) + x(i); % 这位新客户到达银行的时间 = 上一个客户到达银行的时间 + 这位新客户和上一个客户到达的时间间隔
b(i) = max(c(i),e(i-1)); % 这个新客户开始服务的时间取决于其到达时间和上一个客户结束服务的时间
end
n = i-1; % n表示银行一天8小时一共服务的客户人数
t = w/n; % 客户的平均等待时间
disp(['银行一天8小时一共服务的客户人数为: ',num2str(n)])
disp(['客户的平均等待时间为: ',num2str(t)])
toc %计算tic和toc中间部分的代码的运行时间
问题一求解结果：
银行一天8小时一共服务的客户人数为: 46
客户的平均等待时间为: 35.7262
问题二：
clear
tic %计算tic和toc中间部分的代码的运行时间
day = 100; % 假设模拟100天
n = zeros(day,1); % 初始化用来保存每日接待客户数结果的矩阵
t = zeros(day,1); % 初始化用来保存每日客户平均等待时长的矩阵
for k = 1:day
i = 1; % i表示第i个客户，最开始取i=1
w = 0; % w用来表示所有客户等待的总时间，初始化为0
e0 = 0; c0 = 0; % 初始化e0和c0为0
x(1) = exprnd(10); % 第0个客户(假想的)和第1个客户到达的时间间隔
c(1) = c0 + x(1); % 第1个客户到达的时间
b(1) = c(1); % 第1个客户的开始服务的时间
while b(i) <= 480 % 开始设置循环，只要第i个顾客开始服务的时间(时刻)小于480，就可以对其服务（银行每天工作8小时，折换为分钟就是480分钟）
y(i) = normrnd(10,2); % 第i个客户的服务持续时间，服从均值为10方差为4(标准差为2)的正态分布
if y(i) < 1 % 根据题目的意思：若服务持续时间不足一分钟，则按照一分钟计算
y(i) = 1;
end
e(i) = b(i) + y(i); % 第i个客户结束服务的时间 = 第i个客户开始服务的时间 + 第i个客户的服务持续时间
wait(i) = b(i) - c(i); % 第i个客户等待的时间 = 第i个客户开始服务的时间 - 第i个客户到达银行的时间
w = w + wait(i); % 更新所有客户等待的总时间
i = i + 1; % 增加一名新的客户
x(i) = exprnd(10); % 这位新客户和上一个客户到达的时间间隔
c(i) = c(i-1) + x(i); % 这位新客户到达银行的时间 = 上一个客户到达银行的时间 + 这位新客户和上一个客户到达的时间间隔
b(i) = max(c(i),e(i-1)); % 这个新客户开始服务的时间取决于其到达时间和上一个客户结束服务的时间
end
n(k) = i-1; % n(k)表示银行第k天服务的客户人数
t(k) = w/n(k); % t(k)表示该银行第k天客户的平均等待时间
end
disp([num2str(day),'个工作日中，银行每日平均服务的客户人数为: ',num2str(mean(n))])
disp([num2str(day),'个工作日中，银行每日客户的平均等待时间为: ',num2str(mean(t))])
toc %计算tic和toc中间部分的代码的运行时间
问题二求解结果：
100个工作日中，银行每日平均服务的客户人数为: 43.18
100个工作日中，银行每日客户的平均等待时间为: 31.1212
非线性规划问题
问题背景
matlab求解
clc,clear;
tic %计算tic和toc中间部分的代码的运行时间
n=10000000; %生成的随机数组数
x1=unifrnd(20,30,n,1); % 生成在[20,30]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=x1 - 10;
x3=unifrnd(-10,16,n,1); % 生成在[-10,16]之间均匀分布的随机数组成的n行1列的向量构成x3
fmax=-inf; % 初始化函数f的最大值为负无穷（后续只要找到一个比它大的我们就对其更新）
for i=1:n
x = [x1(i), x2(i), x3(i)]; %构造x向量, 这里千万别写成了：x =[x1, x2, x3]
if (-x(1)+2*x(2)+2*x(3)>=0) & (x(1)+2*x(2)+2*x(3)<=72) % 判断是否满足条件
result = x(1)*x(2)*x(3); % 如果满足条件就计算函数值
if result > fmax % 如果这个函数值大于我们之前计算出来的最大值
fmax = result; % 那么就更新这个函数值为新的最大值
X = x; % 并且将此时的x1 x2 x3保存到一个变量中
end
end
end
disp(strcat('蒙特卡罗模拟得到的最大值为',num2str(fmax)))
disp('最大值处x1 x2 x3的取值为：')
disp(X)
toc %计算tic和toc中间部分的代码的运行时间
求解结果
蒙特卡罗模拟得到的最大值为3445.4538
最大值处x1 x2 x3的取值为：
22.5841 12.5841 12.1233
旅行商问题
问题背景
一个售货员必须访问n个城市，这n个城市是一个完全图，售货员需要恰好访问所有城市的一次，并且回到最终的城市。城市于城市之间有一个旅行费用，售货员希望旅行费用之和最少。
matlab求解
clear;clc
% 只有10个城市的简单情况
coord =[0.6683 0.6195 0.4 0.2439 0.1707 0.2293 0.5171 0.8732 0.6878 0.8488 ;
0.2536 0.2634 0.4439 0.1463 0.2293 0.761 0.9414 0.6536 0.5219 0.3609]' ; % 城市坐标矩阵，n行2列
% 38个城市，TSP数据集网站(http://www.tsp.gatech.edu/world/djtour.html) 上公测的最优结果6656。
% coord = [11003.611100,42102.500000;11108.611100,42373.888900;11133.333300,42885.833300;11155.833300,42712.500000;11183.333300,42933.333300;11297.500000,42853.333300;11310.277800,42929.444400;11416.666700,42983.333300;11423.888900,43000.277800;11438.333300,42057.222200;11461.111100,43252.777800;11485.555600,43187.222200;11503.055600,42855.277800;11511.388900,42106.388900;11522.222200,42841.944400;11569.444400,43136.666700;11583.333300,43150.000000;11595.000000,43148.055600;11600.000000,43150.000000;11690.555600,42686.666700;11715.833300,41836.111100;11751.111100,42814.444400;11770.277800,42651.944400;11785.277800,42884.444400;11822.777800,42673.611100;11846.944400,42660.555600;11963.055600,43290.555600;11973.055600,43026.111100;12058.333300,42195.555600;12149.444400,42477.500000;12286.944400,43355.555600;12300.000000,42433.333300;12355.833300,43156.388900;12363.333300,43189.166700;12372.777800,42711.388900;12386.666700,43334.722200;12421.666700,42895.555600;12645.000000,42973.333300];
n = size(coord,1); % 城市的数目
figure(1) % 新建一个编号为1的图形窗口
plot(coord(:,1),coord(:,2),'o'); % 画出城市的分布散点图
for i = 1:n
text(coord(i,1)+0.01,coord(i,2)+0.01,num2str(i)) % 在图上标上城市的编号（加上0.01表示把文字的标记往右上方偏移一点）
end
hold on % 等一下要接着在这个图形上画图的
d = zeros(n); % 初始化两个城市的距离矩阵全为0
for i = 2:n
for j = 1:i
coord_i = coord(i,:); x_i = coord_i(1); y_i = coord_i(2); % 城市i的横坐标为x_i，纵坐标为y_i
coord_j = coord(j,:); x_j = coord_j(1); y_j = coord_j(2); % 城市j的横坐标为x_j，纵坐标为y_j
d(i,j) = sqrt((x_i-x_j)^2 + (y_i-y_j)^2); % 计算城市i和j的距离
end
end
d = d+d'; % 生成距离矩阵的对称的一面
min_result = +inf; % 假设最短的距离为min_result，初始化为无穷大，后面只要找到比它小的就对其更新
min_path = [1:n]; % 初始化最短的路径就是1-2-3-...-n
N = 10000000; % 蒙特卡罗模拟的次数
for i = 1:N % 开始循环
result = 0; % 初始化走过的路程为0
path = randperm(n); % 生成一个1-n的随机打乱的序列
for i = 1:n-1
result = d(path(i),path(i+1)) + result; % 按照这个序列不断的更新走过的路程这个值
end
result = d(path(1),path(n)) + result; % 别忘了加上从最后一个城市返回到最开始那个城市的距离
if result < min_result % 判断这次模拟走过的距离是否小于最短的距离，如果小于就更新最短距离和最短的路径
min_path = path;
min_result = result
end
end
min_path
min_path = [min_path,min_path(1)]; % 在最短路径的最后面加上一个元素，即以上是关于数学建模暑期集训17：蒙特卡洛法的主要内容，如果未能解决你的问题，请参考以下文章