备战数学建模41-蒙特卡罗模拟(攻坚战5)

蒙特卡罗⽅法于20世纪40年代美国在第⼆次世界⼤战中研制原⼦弹的“曼哈顿计划”计划的成员S.M.乌拉姆和J.冯·诺伊曼⾸先提出。数学家冯·诺伊曼⽤驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种⽅法，为它蒙上了⼀层神秘⾊彩。在这之前，蒙特卡罗⽅法就已经存在。1777年，法国Buffon提出⽤投针实验的⽅法求圆周率，这被认为是蒙特卡罗⽅法的起源蒙特卡罗⽅法⼜称统计模拟法，是⼀种随机模拟⽅法，以概率和统计理论⽅法为基础的⼀种计算⽅法，是使⽤随机数（或更常⻅的伪随机数）来解决很多计算问题的⽅法。将所求解的问题同⼀定的概率模型相联系，⽤电⼦计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这⼀⽅法的概率统计特征，故借⽤赌城蒙特卡罗命名。

一、了解蒲(布)丰投针实验

二、蒙特卡洛模拟概述

2.1、蒙特卡洛定义

2.2、蒙特卡洛方法的提出及基本原理

2.3、蒙特卡洛方法的讨论

三、蒙特卡洛模拟的应用实例

3.1、蒙特卡洛模拟三门问题

3.2、蒙特卡洛模拟排队论问题

3.3、蒙特卡洛模拟有约束的非线性规划问题

3.4、蒙特卡洛模拟书店买书问题(0-1规划)

3.5、蒙特卡洛模拟导弹追踪问题

3.6、蒙特卡洛模拟旅行商问题(Travling saleman problem,TSP)

四、使用蒙特卡洛模拟法解决问题

4.1、蒙特卡洛模拟求解自然常数e

4.2、蒙特卡洛模拟求解非线性规划问题

4.3、蒙特卡洛模拟求解方案经济性选择问题

一、了解蒲(布)丰投针实验

我们看一下布丰投针实验，就是画间距为a得平行线，在上面投针，针得长度为l，最后计算针与平行线相交的概率，这个概率除了和间距a以及针长l有关外，还和圆周率Π有关。

我们看一下布丰实验得证明过程，类似于投针在如下的x<=1/2sin&的范围内的概率，我们用蒙特卡罗方法求解Π，只需要模拟投针过程，求出p，然后通过(2*l)/(a*p)即可计算出Π的值。

我们使用matlab编程，实现蒙特卡洛模拟布丰投针实验，模拟投针10000次，求出落入指定区域的概率，然后通过公式计算出Π值，具体得代码如下：

l =  0.520;     % 针的长度（任意给的）
a = 1.314;    % 平行线的宽度(大于针的长度l即可)
n = 10000;    % 做n次投针试验，n越大求出来的pi越准确
m = 0;    % 记录针与平行线相交的次数
x = rand(1, n) * a / 2 ;   % 在[0, a/2]内服从均匀分布随机产生n个数， x中每一个元素表示针的中点和最近的一条平行线的距离
phi = rand(1, n) * pi;    % 在[0, pi]内服从均匀分布随机产生n个数，phi中的每一个元素表示针和最近的一条平行线的夹角
axis([0,pi, 0,a/2]);   box on;  % 画一个坐标轴的框架，x轴位于0-pi，y轴位于0-a/2， 并打开图形的边框
for i=1:n  % 开始循环，依次看每根针是否和直线相交if x(i) <= l / 2 * sin(phi (i))     % 如果针和平行线相交m = m + 1;    % 那么m就要加1plot(phi(i), x(i), 'r.')   % 模仿书上的那个图，横坐标为phi，纵坐标为x , 用红色的小点进行标记hold on  % 在原来的图形上继续绘制end
end
p = m / n;    % 针和平行线相交出现的频率
mypi = (2 * l) / (a * p);  % 我们根据公式计算得到的pi
disp(['蒙特卡罗方法得到pi为：', num2str(mypi)])

模拟的效果如下：

二、蒙特卡洛模拟概述

2.1、蒙特卡洛定义

蒙特卡罗⽅法⼜称统计模拟法，是⼀种随机模拟⽅法，以概率和统计理论⽅法为基础的⼀种计算⽅法，是使⽤随机数（或更常⻅的伪随机数）来解决很多计算问题的⽅法。将所求解的问题同⼀定的概率模型相联系，⽤电⼦计算机实现统计模拟或抽样，以获得问题的近似解。为象征性地表明这⼀⽅法的概率统计特征，故借⽤赌城蒙特卡罗命名。

2.2、蒙特卡洛方法的提出及基本原理

由⼤数定理可知，当样本容量⾜够⼤时，事件的发⽣频率即为其概率。

2.3、蒙特卡洛方法的讨论

算法（Algorithm）是指解题⽅案的准确⽽完整的描述，是⼀系列解决问题的清晰指令。蒙特卡罗准确的来说只是⼀种思想，或者是是⼀种⽅法。如果我们所求解的问题与概率模型有⼀定的关联，那么我们就可以使⽤计算机多次模拟事件发⽣，以获得问题的近似解。从数学建模⻆度来看，⼤家千万别认为蒙特卡罗有⼀个通⽤的代码。每个问题对应的代码都是不同的，我们分析清楚题⽬后，就要⾃⼰进⾏编写适⽤于这个题⽬的代码。

枚举法是我们中学就接触的算法，就是把所有可能发⽣情况都考虑进去，最终计算出来⼀个确定结果。这就与蒙特卡罗⽅法的想法很类似，蒙特卡罗法模拟的次数越多，计算的就越准确。由于⽣活中有许多事件发⽣的结果都有⽆限种可能（例如⼀个连续分布的取值），因此我们不可能枚举出所有的结果，这时候就只能通过蒙特卡罗模拟，将⼀个不确定性的问题转化成很多个确定性问题，并得到⼀个近似解，因此蒙特卡罗算法也可以看成是枚举法的⼀种变异。

三、蒙特卡洛模拟的应用实例

3.1、蒙特卡洛模拟三门问题

我们可以看一下三门问题，就是三个门，你选择其中一扇门，主持人给你打开了一个空门，问你要不要改选其它门，这个问题是个概率问题，我们可以通过蒙特卡洛方法进行模拟，然后观察是改选获奖的概率大，还是不改选获奖的概率大。

1）我们考虑两种情况，第一种是默认已经获奖，认为是一个条件概率，即计算改选获奖和不改选获奖的概率。

%在成功的条件下的概率
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~=yif 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)]);

经过上述的10万次模拟开门过程，可以发现应该改变，改变的获奖概率是不改变的两倍。

2)我们考虑第二种情况，就是考虑不获奖的情况，就需要另外用一个变量去记录不获奖的次数，这样根据获奖和不获奖的次数，就可以计算出概率，matlab代码如下：

%考虑失败情况的代码(无条件概率)
n = 100000;  % n代表蒙特卡罗模拟重复次数
a = 0;  % a表示不改变主意时能赢得汽车的次数
b = 0;  % b表示改变主意时能赢得汽车的次数
c = 0;  % c表示没有获奖的次数
for i= 1 : n  % 开始模拟n次x = randi([1,3]);  % 随机生成一个1-3之间的整数x表示汽车出现在第x扇门后y = randi([1,3]);  % 随机生成一个1-3之间的整数y表示自己选的门change = randi([0, 1]); % change =0  不改变主意，change = 1 改变主意% 下面分为两种情况讨论：x=y和x~=yif x == y   % 如果x和y相同，那么我们只有不改变主意时才能赢if change == 0  % 不改变主意a = a + 1; else  % 改变了主意c= c+1;endelse  % x ~= y ，如果x和y不同，那么我们只有改变主意时才能赢if change == 0  % 不改变主意c = c + 1; else  % 改变了主意b= b + 1;endend
end
disp(['蒙特卡罗方法得到的不改变主意时的获奖概率为：', num2str(a/n)]);
disp(['蒙特卡罗方法得到的改变主意时的获奖概率为：', num2str(b/n)]);
disp(['蒙特卡罗方法得到的没有获奖的概率为：', num2str(c/n)]);

通过运行结果我们可以发现，获奖和不获奖各占50%,但是改变主意的获奖率仍然是不改变主意获奖的概率的2倍。

3.2、蒙特卡洛模拟排队论问题

我们先看一下题目，排队论问题就是先到先服务原则，一个先到先服务的串行过程，每个顾客能否服务取决于上一个顾客是否服务结束，我们通过模拟用户到来的时间间隔和每个顾客服务的时间间隔可以求出客户的平均等待时间。

我们在模拟之前需要分析一下题目，主要引入了Ci，bi和ei三个变量，通过排队论题目可以得出第i个客户的到达时间=第i-1个客户的到达时间+时间间隔xi，第i个客户的服务结束时间ei=开始时间+服务持续时间，第i个客户的开始服务时间=max(第i个客户的达到时间，第i-1个客户的服务结束时间)。由这些分析，我们可以使用蒙特卡洛方法进行模拟。

1）我们使用蒙特卡洛方法模拟1个工作日，即480分钟，小于1分钟的，就算作一分钟，客户到达时间间隔假设服从均值为10的指数分布，每个顾客的服务时间通过随机生成的均值为10方差为4的正态分布，最后计算接待客户的总人数和客户的平均等待时间。

%问题1的代码
clc
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个客户到达的时间间隔(均值为10的指数分布)
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;ende(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中间部分的代码的运行时间

运行结果如下，由于每次都是随机模拟的，所以生成的结果大同小异，可以发现这种单个窗口串行的结构使得每位用户平均等待20分钟左右。

我们再来看一下第2问，就是模拟100个工作日，然后计算每天的服务人数和平均等待时间，通过大量的模拟，可以使得模拟结果更加准确，由大数定律可知，当样本容量足够大时，频率就可以近似等于概率。就是外层加个循环，记录100天的，然后求均值即可。

%问题2的代码
clc
clear
tic  %计算tic和toc中间部分的代码的运行时间
day = 100;  % 假设模拟100天
n = zeros(day,1); % 初始化用来保存每日接待客户数结果的矩阵
t = zeros(day,1); % 初始化用来保存每日客户平均等待时长的矩阵
for k = 1:dayi = 1;  % i表示第i个客户，最开始取i=1w = 0;  % w用来表示所有客户等待的总时间，初始化为0e0 = 0;  c0 = 0;   % 初始化e0和c0为0x(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;ende(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)); % 这个新客户开始服务的时间取决于其到达时间和上一个客户结束服务的时间endn(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中间部分的代码的运行时间

模拟的结果如下，每个客户的等待时间达到了30分钟，这个可以说相当可怕，提个鸡肋的建议，多加几个窗口吧，太不容易了。

3.3、蒙特卡洛模拟有约束的非线性规划问题

一般的规划类问题，包括目标函数，决策变量和约束条件，对于规划类问题，用蒙特卡洛方法进行模拟，主要思路如下：需要给出决策变量的大致范围，在这个范围内生成随机数，验证满足条件的决策变量，将这些代入目标函数，找到最大值或则最小值。

对于上面的例题，我们可以先进行如下的推导，可以得到x1,x2,x3三个决策变量的范围，通过在范围内随机生成决策变量，筛选满足条件的决策变量，代入目标函数，求出目标函数最值。

对于上述的例题，使用了1千万组随机数进行模拟，对于满足约束条件的数据，代入目标函数，找到最大值，具体的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:nx = [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保存到一个变量中endend
end
disp(strcat('蒙特卡罗模拟得到的最大值为',num2str(fmax)))
disp('最大值处x1 x2 x3的取值为：')
disp(X)
toc %计算tic和toc中间部分的代码的运行时间

我们可以看一下具体的运行结果，我们通过这个得到的结果，可以对决策变量的范围进行缩小，这样可以模拟出更加准确的结果。

下面根据上述计算出的决策变量的值，对设定的决策变量的范围值进行缩小，这样模拟出来的值会更接近准确值，具体如下:

clc,clear;
tic %计算tic和toc中间部分的代码的运行时间
n=10000000; %生成的随机数组数
x1=unifrnd(22,23,n,1);  % 生成在[22,23]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=x1 - 10;
x3=unifrnd(11,13,n,1);  % 生成在[11,13]之间均匀分布的随机数组成的n行1列的向量构成x3
fmax=-inf; % 初始化函数f的最大值为负无穷（后续只要找到一个比它大的我们就对其更新）
for i=1:nx = [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保存到一个变量中endend
end
disp(strcat('蒙特卡罗模拟得到的最大值为',num2str(fmax)))
disp('最大值处x1 x2 x3的取值为：')
disp(X)
toc %计算tic和toc中间部分的代码的运行时间

运行结果如下：

3.4、蒙特卡洛模拟书店买书问题(0-1规划)

我们看一下下面的买书问题，就是从书店买书，一共需要买5本书，每本书买一次即可，在一家店买多本书也之首一次运费，现在让你设计一个选购方案，使得最省钱。

我们看一下上述规划问题的解题思路，变量i和j分别表示6个商城和5本书，xij表示第i个同学是否在第j家买书，买了为1，不买为0，同时为了约束每本书都只买一次，需要加个约束，另外对于目标函数，主要考虑书的价格和运费，求出总的费用最小。

下面使用蒙特卡洛方法进行模拟整个过程，计算出总的费用=书费+运费，10万次模拟，使得最终的最小值近似等于我们要求得结果。

clc
clear
min_money = +Inf;  % 初始化最小的花费为无穷大，后续只要找到比它小的就更新
min_result = randi([1, 6],1,5);  % 初始化五本书都在哪一家书店购买，后续我们不断对其更新
%若min_result = [5 3 6 2 3]，则解释为：第1本书在第5家店买，第2本书在第3家店买，第3本书在第6家店买，第4本书在第2家店买，第5本书在第3家店买
n = 100000;  % 蒙特卡罗模拟的次数
M = [18     39 29  48  5924    45  23  54  4422    45  23  53  5328    47  17  57  4724    42  24  47  5927    48  20  55  53];  % m_ij  第j本书在第i家店的售价
freight = [10 15 15 10 10 15];  % 第i家店的运费
for k = 1:n  % 开始循环result = randi([1, 6],1,5); % 在1-6这些整数中随机抽取一个1*5的向量，表示这五本书分别在哪家书店购买index = unique(result);  % 在哪些商店购买了商品，因为我们等下要计算运费money = sum(freight(index)); % 计算买书花费的运费% 计算总花费：刚刚计算出来的运费 + 五本书的售价for i = 1:5   money = money + M(result(i),i);  endif money < min_money  % 判断刚刚随机生成的这组数据的花费是否小于最小花费，如果小于的话min_money = money;  % 我们更新最小的花费min_result = result; % 用这组数据更新最小花费的结果end
end
disp(min_money)   % 18+39+48+17+47+20
disp(min_result)

我们看一下，最后总的最小花费为189，买书方案5本书分别在商城1，1，4，1，4购买，最后得花费最小。

3.5、蒙特卡洛模拟导弹追踪问题

我们来看一下这个导弹追踪问题，B船沿着东北方向逃逸，A船始终瞄准B船，向B船发射导弹，计算导弹能否击B船？

我们仔细分析一下这个题目，因为A船得导弹始终对准B船，那么A船设为原点，则B船的坐标很容易得到，导弹的飞行是一个曲线，那么这个切线就是导弹速度的方向，速度方向可以分解为水平和竖直两个方向，这样就可以写出速度公式。

有了上面的公式，我们就可以考虑建立近似的模型，然后使用蒙特卡洛方法进行模拟，我们奖时间间隔划分的很小，就可以模拟一个连续的时间了。首先可以更新B船的位置，然后根据B船的位置可以计算出斜率tana，然后可以推出sina和cosa，这样就可以更新导弹的位置，由此不停地迭代，直到导弹和船的距离小于一个给定值，则认为导弹击中了船。

代码如下：

clear;clc
v=200; % 任意给定B船的速度（后期我们可以再改的）
dt=0.0000001; % 定义时间间隔
x=[0,20]; % 定义导弹和B船的横坐标分别为x(1)和x(2)
y=[0,0]; % 定义导弹和B船的纵坐标分别为y(1)和y(2)
t=0; % 初始化导弹击落B船的时间
d=0; % 初始化导弹飞行的距离
m=sqrt(2)/2;   % 将sqrt(2)/2定义为一个常量，使后面看起来很简洁
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 导弹与B船的距离
while(dd>=0.001)  % 只要两者的距离足够大，就一直循环下去。（两者距离足够小时表示导弹击中，这里的临界值要结合dt来取，否则可能导致错过交界处的情况）t=t+dt; % 更新导弹击落B船的时间d=d+3*v*dt; % 更新导弹飞行的距离x(2)=20+t*v*m;  y(2)=t*v*m;   % 计算新的B船的位置 （注：m=sqrt(2)/2）dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);  % 更新导弹与B船的距离tan_alpha=(y(2)-y(1))/(x(2)-x(1));   % 计算斜率，即tan(α)cos_alpha=sqrt(1/(1+tan_alpha^2));   % sec(α)^2 = (1+tan(α)^2)sin_alpha=sqrt(1-cos_alpha^2);  % sin(α)^2 +cos(α)^2 = 1x(1)=x(1)+3*v*dt*cos_alpha;   y(1)=y(1)+3*v*dt*sin_alpha; % 计算新的导弹的位置if d>50  % 导弹的有效射程为50个单位disp('导弹没有击中B船');break;  % 退出循环endif d<=50 & dd<0.001   % 导弹飞行的距离小于50个单位且导弹和B船的距离小于0.001（表示击中）disp(['导弹飞行',num2str(d),'单位后击中B船'])disp(['导弹飞行的时间为',num2str(t*60),'分钟'])end
end

运行结果如下：

下面是绘制导弹追踪B船的整个过程，代码如下：

clear;clc
v=200; % 任意给定B船的速度（后期我们可以再改的）
dt=0.0000001; % 定义时间间隔
x=[0,20]; % 定义导弹和B船的横坐标分别为x(1)和x(2)
y=[0,0]; % 定义导弹和B船的纵坐标分别为y(1)和y(2)
t=0; % 初始化导弹击落B船的时间
d=0; % 初始化导弹飞行的距离
m=sqrt(2)/2;   % 将sqrt(2)/2定义为一个常量，使后面看起来很简洁
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 导弹与B船的距离
for i=1:2plot(x(i),y(i),'.k','MarkerSize',1);  % 画出导弹和B船所在的坐标，点的大小为1，颜色为黑色(k)，用小点表示grid on;  % 打开网格线hold on;  % 不关闭图形，继续画图
end
axis([0 30 0 10])  % 固定x轴的范围为0-30  固定y轴的范围为0-10
k = 0;  % 引入一个变量  为了控制画图的速度（因为Matlab中画图的速度超级慢）
while(dd>=0.001)  % 只要两者的距离足够大，就一直循环下去。（两者距离足够小时表示导弹击中，这里的临界值要结合dt来取，否则可能导致错过交界处的情况）t=t+dt; % 更新导弹击落B船的时间d=d+3*v*dt; % 更新导弹飞行的距离x(2)=20+t*v*m;  y(2)=t*v*m;   % 计算新的B船的位置 （注：m=sqrt(2)/2）dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);  % 更新导弹与B船的距离tan_alpha=(y(2)-y(1))/(x(2)-x(1));   % 计算斜率，即tan(α)cos_alpha=sqrt(1/(1+tan_alpha^2));   % 利用公式：sec(α)^2 = (1+tan(α)^2)  计算出cos(α)sin_alpha=sqrt(1-cos_alpha^2);  % 利用公式： sin(α)^2 +cos(α)^2 = 1  计算出sin(α)x(1)=x(1)+3*v*dt*cos_alpha;   y(1)=y(1)+3*v*dt*sin_alpha;   % 计算新的导弹的位置k = k +1 ;  if mod(k,500) == 0   % 每刷新500次时间就画出下一个导弹和B船所在的坐标  mod(m,n)表示求m/n的余数for i=1:2plot(x(i),y(i),'.k','MarkerSize',1);hold on; % 不关闭图形，继续画图endpause(0.001);  % 暂停0.001s后再继续下面的操作endif d>50  % 导弹的有效射程为50个单位disp('导弹没有击中B船');break;  % 退出循环endif d<=50 & dd<0.001   % 导弹飞行的距离小于50个单位且导弹和B船的距离小于0.001（表示击中）disp(['导弹飞行',num2str(d),'个单位后击中B船'])disp(['导弹飞行的时间为',num2str(t*60),'分钟'])end
end

3.6、蒙特卡洛模拟旅行商问题(Travling saleman problem,TSP)

旅行商问题也是一个比较热门的问题，就是从一个城市开始走，访问所有城市，所有城市有且只走一次，最后回到原点，找出一种走法，使得费用最低，即边的总权重最小。

模拟走的过程，累加权重，找出最小的，然后绘图，我们此次之使用了10个城市进行模拟，城市数量太多，模拟效果并不好，如下：

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:ntext(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_icoord_j = coord(j,:);   x_j = coord_j(1);     y_j = coord_j(2);  % 城市j的横坐标为x_j，纵坐标为y_jd(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 k = 1:N  % 开始循环result = 0;  % 初始化走过的路程为0path = randperm(n);  % 生成一个1-n的随机打乱的序列for i = 1:n-1  result = d(path(i),path(i+1)) + result;  % 按照这个序列不断的更新走过的路程这个值endresult = d(path(1),path(n)) + result;  % 别忘了加上从最后一个城市返回到最开始那个城市的距离if result < min_result  % 判断这次模拟走过的距离是否小于最短的距离，如果小于就更新最短距离和最短的路径min_path = path;min_result = resultend
end
min_path
min_path = [min_path,min_path(1)];   % 在最短路径的最后面加上一个元素，即第一个点（我们要生成一个封闭的图形）
n = n+1;  % 城市的个数加一个（紧随着上一步）
for i = 1:n-1 j = i+1;coord_i = coord(min_path(i),:);   x_i = coord_i(1);     y_i = coord_i(2); coord_j = coord(min_path(j),:);   x_j = coord_j(1);     y_j = coord_j(2);plot([x_i,x_j],[y_i,y_j],'-')    % 每两个点就作出一条线段，直到所有的城市都走完pause(0.5)  % 暂停0.5s再画下一条线段hold on
end

四、使用蒙特卡洛模拟法解决问题

4.1、蒙特卡洛模拟求解自然常数e

我们使⽤蒙特卡罗的⽅法对这个问题进⾏模拟，并估计出⾃然常数e的值，这个和模拟Π很像。

我们可以用随机生成的数据模拟自己的卡片和打乱顺序后的卡片，最终每个人拿到都不是自己的卡片的次数除以总次数，然后取倒数，就可以得到我们求解的e。

clear;clc
tic  %计算tic和toc中间部分的代码的运行时间
n = 1000000;  % 蒙特卡洛的次数（理论上n取得越大，计算出来的结果越精确）
m = 0;   % 每个人拿到的都不是自己卡片的次数（频数）
people = 100;   % 假设一共有100个人玩这个游戏 （任给的）
for i = 1: n  % 开始循环if isempty(find(randperm(people) - [1:people] == 0))  % 如果每个人拿到的都不是自己的卡片m = m + 1;  % 那么次数就加1end
end
frequency = m / n;  % 每个人拿到的都不是自己卡片的频率（概率）
disp(['自然常数e的蒙特卡罗模拟值为：', num2str(1 / frequency)])  % 注：自然常数真实值约为2.7182
toc  %计算tic和toc中间部分的代码的运行时间

我们用100个人，进行100万次模拟，运行的结果如下所示：

4.2、蒙特卡洛模拟求解非线性规划问题

我们看一下这个非线性规划问题，这个问题看起来不是很复杂，我们使用蒙特卡洛模拟，给出决策变量的大概范围，在范围生成数据进行模拟，满足约束的即为可行解，我们讲可行解代入目标函数，通过大量的模拟，找到一个可行解代入目标函数得到最小值，即为近似可行解。

使用蒙特卡洛进行模拟的matlab代码如下，当然，可以根据模拟的结果，对决策变量的范围进行缩小，然后再次模拟，会得到更加精确的值。

clc,clear;
format long g   %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
tic %计算tic和toc中间部分的代码的运行时间
n=10000000; %生成的随机数组数
x1=unifrnd(0,16,n,1);  % 生成在[0,16]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=unifrnd(0,8,n,1);  % 生成在[0,8]之间均匀分布的随机数组成的n行1列的向量构成x2
fmin=+inf; % 初始化函数f的最小值为正无穷（后续只要找到一个比它小的我们就对其更新）
for i=1:nx = [x1(i), x2(i)];  %构造x向量, 这里千万别写成了：x =[x1, x2]if (3*x(1)+x(2)>9)  &  (x(1)+2*x(2)<16)     % 判断是否满足条件result = 2*(x(1)^2)+x(2)^2-x(1)*x(2)-8*x(1)-3*x(2);  % 如果满足条件就计算函数值if  result  < fmin  % 如果这个函数值小于我们之前计算出来的最小值fmin = result;  % 那么就更新这个函数值为新的最小值X = x;  % 并且将此时的x1 x2 保存到相应的变量中endend
end
disp(strcat('蒙特卡罗模拟得到的最小值为',num2str(fmin)))
disp('最小值处x1 x2的取值为：')
disp(X)
toc %计算tic和toc中间部分的代码的运行时间

4.3、蒙特卡洛模拟求解方案经济性选择问题

我们看一下这个应用题，第一眼看题的时候觉得好搞笑，更换4只成本高而且耗时，而且没坏就换真浪费，直接哪个坏了换哪个不就好了，其实仔细看题会发现，到达寿命后，电子管可能随时坏，如果直接换掉4个，反而可能节约时间。

我们使用蒙特卡洛方法，分别对两种方案进行模拟，对于第一种方案，随机生成四个1000~2000h之间的数字模拟四个电子管寿命，在模拟时间T内，每次找出最短寿命的，更新当前时间，更新方案一的花费，更新经过这些时间后剩余电子管的寿命，同时讲坏的电子管更换为新的寿命。

对于第二种方案，每次都更新四个电子管的寿命，更细时间和方案二的花费。

clear;clc
T = 100000000;   % T表示模拟的总时间（单位为小时）
t = 0;   % 初始化当前时刻为0小时
c1 = 0; c2 = 0;  % 初始化两种方案的总花费都为0%%  方案一
life = randi([1000,2000],1,4);  % 随机生成四个电子管的寿命，假设为整数
while t < T  % 只要现在的时刻没有超过总时刻，就不断循环下去result = min(life);  % 找出寿命最短的那一个电子管的寿命t = t+result+1;  % 现在的时间更改到有电子管损坏的时刻（加上1表示更换电子管需要花费的时间）c1 = c1 + 20 * 1 +10;  % 更新方案一的花费 k = find(life == result,1);   % 找到哪一个电子管是坏的life = life - result -1; % 更新所有电子管的寿命（这里不减去1也是可以的，减少了1也无所谓，对结果的影响很小）    life(k) = randi([1000,2000]);  % 把坏掉的那个电子管的寿命重置
end%%  方案二
t = 0;   % 初始化当前时刻为0小时
while t < T  % 只要现在的时刻没有超过总时刻，就不断循环下去life = randi([1000,2000],1,4); % 随机生成四个电子管的寿命，假设为整数result = min(life); % 找出寿命最小的那一个电子管的寿命t = t+result+2;  % 现在的时间更改到有电子管损坏的时刻（加上2表示更换所有电子管需要花费的时间）c2 =c2 + 20 * 2 +40;  % 更新方案二的花费
end%% 两种方案的花费
c1
c2

通过取较大的时间T进行模拟，可以发现，一次更换四个电子管反而更经济！！！