蒙特卡罗经典例子

陌上花开 ~

已于 2022-07-27 23:12:36 修改

阅读量2.8k

点赞数 2

分类专栏：数学建模文章标签：概率论深度学习机器学习

于 2021-08-29 20:38:31 首次发布

本文链接：https://blog.csdn.net/qq_45550375/article/details/119827373

版权

蒙特卡罗模拟信息技术概率问题优化问题数值计算

关键词由CSDN通过智能技术生成

数学建模专栏收录该内容

14 篇文章 66 订阅

订阅专栏

布丰投针实验求π

代码：

l =  0.520;   % 针长
a = 1.314;    % 平行线的宽度(大于针长即可)
n = 10000;    % 做n次投针试验
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;  
        plot(phi(i), x(i), 'r.')   %横坐标为phi，纵坐标为x , 用红色的小点进行标记
        hold on  % 在原来的图形上继续绘制
    end
end
p = m / n;    % 针和平行线相交出现的频率
mypi = (2 * l) / (a * p);  % 根据公式计算得到的pi
disp(['蒙特卡罗方法得到pi为：', num2str(mypi)])

结果：

蒙特卡罗方法得到pi为：3.1358

在这里插入图片描述
将上述代码重复100次后求平均值：

result = zeros(100,1);  % 初始化保存100次结果的矩阵
l =  0.520;     a = 1.314;
n = 1000000;    
for num = 1:100
    m = 0;  
    x = rand(1, n) * a / 2 ;
    phi = rand(1, n) * pi;
    for i=1:n
        if x(i) <= l / 2 * sin(phi (i))
            m = m + 1;
        end
    end
    p = m / n;
    mypi = (2 * l) / (a * p);
    result(num) = mypi;  % 把求出来的myphi保存到结果矩阵中
end
mymeanpi = mean(result);  % 计算result矩阵中保存的100次结果的均值
disp(['蒙特卡罗方法得到pi为：', num2str(mymeanpi)])

结果：

蒙特卡罗方法得到pi为：3.1412

错排问题求e

在这里插入图片描述

预备知识：

randperm(5)  % 生成1-5组成的一个随机序列
%      3     5     1     2     4
%      1     4     5     3     2

% 假设a是一个向量，那么find(a)可以用来返回这个向量中非零元素的下标，如果a中所有元素都为0，则返回空值
find([1,5,6,0,8,0,-5])  %      1     2     3     5     7
find([0,0,0,0,0])  %   空的 1×0 double 行矢量

% (3) 矩阵（或向量）和常量的比较运算可返回逻辑矩阵（或向量）(元素全为0和1)
[1,5,6,0,8,0,-5] > 0      %    1   1   1   0   1   0   0
[1,5,6,0,8,0,-5] == 0    %    0   0   0   1   0   1   0

% (4) isempty(A)函数可以用来判断A是否为空, 如果A为空， isempty(A) 返回逻辑值1(true)，否则返回逻辑值0(false)。
isempty(find([0,0,0,0,0]))   %    1
isempty(find([0,1,0,0,0]))   %    0
isempty([0,0,0,0,0])  %它不是空矩阵（空矩阵是指里面没有元素）

求解：

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

结果：

自然常数e的蒙特卡罗模拟值为：2.7167
历时 5.235055 秒

三门问题

在这里插入图片描述
蒙特卡罗求在成功的条件下的概率：

n = 100000;  % 蒙特卡罗模拟重复次数
a = 0;  % 不改变主意时能赢得汽车的次数
b = 0;  % 改变主意时能赢得汽车的次数
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.33452
蒙特卡罗方法得到的改变主意时的获奖概率为：0.66548

蒙特卡罗考虑失败情况的概率(无条件概率)：

n = 100000;  % 蒙特卡罗模拟重复次数
a = 0;  % 不改变主意时能赢得汽车的次数
b = 0;  % 改变主意时能赢得汽车的次数
c = 0;  % 没有获奖的次数
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~=y
    if x == y   % 如果x和y相同，那么我们只有不改变主意时才能赢
        if change == 0  % 不改变主意
        	a = a + 1; 
        else  % 改变了主意
            c= c+1;
        end
    else  % x ~= y ，如果x和y不同，那么我们只有改变主意时才能赢
         if change == 0  % 不改变主意
        	c = c + 1; 
         else  % 改变了主意
            b= b + 1;
         end
    end
end
disp(['蒙特卡罗方法得到的不改变主意时的获奖概率为：', num2str(a/n)]);
disp(['蒙特卡罗方法得到的改变主意时的获奖概率为：', num2str(b/n)]);
disp(['蒙特卡罗方法得到的没有获奖的概率为：', num2str(c/n)]);

结果：

蒙特卡罗方法得到的不改变主意时的获奖概率为：0.16723
蒙特卡罗方法得到的改变主意时的获奖概率为：0.33271
蒙特卡罗方法得到的没有获奖的概率为：0.50006

装备价值预测

现在有一把神器，初始为 1 级，可免费领取（即价值为 0），可花费金币对其升级，每次 10000 金币，最多升到 5 级，问 5 级神器价值多少金币
在这里插入图片描述
预备知识：

% 以一定的概率产生随机数  randsrc(m,n,[alphabet; prob])
% m和n表示生成的随机数矩阵的行数和列数
% alphabet表示需要产生的随机数的数字，用一个行向量表示
% prob表示这些数字出现的概率大小，用一个行向量表示，向量长度和alphabet向量要完全相同, 且这些概率的和要为1

% 例子：要产生1、4、 6这三个数。它们分别出现的概率为 0.1、0.2、0.7
% 要设计程序使得按照这个概率产生10个随机数

alphabet = [1 4 6]; prob = [0.1 0.2 0.7];
randsrc(10,1,[alphabet; prob])

求解：

clear;clc
tic  
% success矩阵储存升级的成功率
success = [0.65 0.2  0.1  0.05  0;
           0.25 0.4  0.2  0.1   0.05;
           0.1  0.2  0.4  0.2   0.1;
           0    0.1  0.3  0.4   0.2] ;
n = 10000; 
MONEY = zeros(n,1);  % 初始化用来存储每次蒙特卡罗计算出来的表示强化费用的向量
for i = 1:n
    rank = 1; 
    money = 0; 
    alphabet = [1 2 3 4 5];   % 用来表示五个等级
    while rank ~= 5  % 只要等级不是5级， 就一直循环下去
        prob =success(rank,:);    % 令生成随机数的概率为第rank行
        rank = randsrc(1,1,[alphabet; prob]);   % 生成一个在1-5中的随机数，表示强化后的等级
        money = money + 10000;  % 更新强化的费用
    end
    MONEY(i) = money;  % 将这次蒙特卡罗的结果保存到MONEY向量中
end
disp(['将武器升级到5级的平均花费为：',num2str(mean(MONEY))])
toc  %计算tic和toc中间部分的代码的运行时间

结果：

将武器升级到5级的平均花费为：160335
历时 2.784793 秒。

模拟排队问题

在这里插入图片描述
问题1求解：

tic  %计算tic和toc中间部分的代码的运行时间
i = 1;  
w = 0;  % w用来表示所有客户等待的总时间，初始化为0
e0 = 0;  c0 = 0;   
x(1) = exprnd(10);  
c(1) = c0 + x(1);  
b(1) = c(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; 
t = w/n; % 客户的平均等待时间
disp(['银行一天8小时一共服务的客户人数为: ',num2str(n)])
disp(['客户的平均等待时间为: ',num2str(t)])
toc  %计算代码运行时间

问题2求解：

tic  
day = 100;  % 假设模拟100天
n = zeros(day,1); % 每日接待客户数
t = zeros(day,1); % 每日客户平均等待时长
for k = 1:day
    i = 1;  
    w = 0;  
    e0 = 0;  c0 = 0;   
    x(1) = exprnd(10);  
    c(1) = c0 + x(1);  
    b(1) = c(1); 
    while b(i) <= 480  
        y(i) = normrnd(10,2); 
        if y(i) < 1  
            y(i) = 1;
        end
        e(i) = b(i) + y(i); 
        wait(i) = b(i) - c(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  %计算代码运行时间

解有约束的非线性规划问题

例一

求满足以下约束条件下函数 $f(x)=x_{1}x_{2}x_{3}$ 的最大值
$t.\begin{cases}-x_{1}+2x_{2}+2x_{3}>=0\\ x_{1}+x_{2}+2x_{3}<=72\\ 10<=x_{2}<=20\\ x_{1}-x_{2}=10\end{cases}$

首先得到一个比较精准的解

代码：

clc,clear;
tic 
n=10000000; %生成的随机数组数
x1=unifrnd(20,30,n,1);  % 生成在[20,30]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=x1 - 10;
x3=unifrnd(-10,16,n,1); 
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

结果：

蒙特卡罗模拟得到的最大值为3445.574
最大值处x1 x2 x3的取值为：
          22.6113760117399          12.6113760117399          12.0829260000158

历时 1.979151 秒。

在比较精准的解的基础上，再次蒙特卡罗：

求解：

clc,clear;
tic 
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:n
    x = [x1(i), x2(i), x3(i)]; 
    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; 
        end
    end
end
disp(strcat('蒙特卡罗模拟得到的最大值为',num2str(fmax)))
disp('最大值处x1 x2 x3的取值为：')
disp(X)
toc

结果：

蒙特卡罗模拟得到的最大值为3445.6014
最大值处x1 x2 x3的取值为：
          22.5887179336818          12.5887179336818          12.1169119252507

历时 0.831883 秒。

例二

在这里插入图片描述
初次寻找最小值：

clc,clear;
format long g   %可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
tic 
n=10000000; 
x1=unifrnd(0,16,n,1);  % 生成在[0,16]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=unifrnd(0,8,n,1); 
fmin=+inf; % 初始化函数f的最小值为正无穷
for i=1:n
    x = [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;  
        end
    end
end
disp(strcat('蒙特卡罗模拟得到的最小值为',num2str(fmin)))
disp('最小值处x1 x2的取值为：')
disp(X)
toc %计算tic和toc中间部分的代码的运行时间

结果：

蒙特卡罗模拟得到的最小值为-15.1429
最小值处x1 x2的取值为：
          2.71485111081151          2.85830659834113

历时 2.435691 秒。

缩小范围重新模拟：

clc,clear;
tic %计算tic和toc中间部分的代码的运行时间
n=10000000; 
x1=unifrnd(2,3,n,1);  % 生成在[2,3]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=unifrnd(2,3,n,1);  
fmin=+inf; % 初始化函数f的最小值为正无穷
for i=1:n
    x = [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;  
        end
    end
end
disp(strcat('蒙特卡罗模拟得到的最小值为',num2str(fmin)))
disp('最小值处x1 x2的取值为：')
disp(X)
toc %计算tic和toc中间部分的代码的运行时间

结果：

蒙特卡罗模拟得到的最小值为-15.1429
最小值处x1 x2的取值为：
          2.71414872138368          2.85692847449316

历时 0.954697 秒。

买书问题（01规划）

在这里插入图片描述
求解：

min_money = +Inf;  % 初始化最小的花费为无穷大
min_result = randi([1, 6],1,5);  % 初始化五本书都在哪一家书店购买，后续我们不断对其更新
n = 100000;  % 蒙特卡罗模拟次数
M = [18	 39	29	48	59
        24	45	23	54	44
        22	45	23	53	53
        28	47	17	57	47
        24	42	24	47	59
        27	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);  
    end
    if money < min_money  
        min_money = money;  % 我们更新最小的花费
        min_result = result; % 用这组数据更新最小花费的结果
    end
end
min_money 
min_result

结果：

min_money =

   189


min_result =

     1     1     4     1     4

导弹追踪问题

位于坐标原点的A船向位于其正东方20个单位的B船发射导弹，导弹始终对准B船，B船以时速V单位（常数）沿东北方向逃逸。若导弹的速度为3V，导弹的射程是50个单位，画出导弹运行的曲线，导弹是否能在射程内击中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; 
d=0; 
m=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船的位置
    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 = 1
    x(1)=x(1)+3*v*dt*cos_alpha;   y(1)=y(1)+3*v*dt*sin_alpha; % 计算新的导弹的位置
    if d>50
        disp('导弹没有击中B船');
        break;  
    end
    if d<=50 && dd<0.001   
        disp(['导弹飞行',num2str(d),'单位后击中B船'])
        disp(['导弹飞行的时间为',num2str(t*60),'分钟'])
    end
end

结果：

导弹飞行27.8019单位后击中B船
导弹飞行的时间为2.7802分钟

画上追击的示意图：

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; 
d=0; 
m=sqrt(2)/2;   
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 导弹与B船的距离
for i=1:2
    plot(x(i),y(i),'.k','MarkerSize',1);  % 画出导弹和B船所在的坐标，点的大小为1，颜色为黑色(k)，点的大小为1
    grid on;  % 打开网格线
    hold on;  % 不关闭图形，继续画图
end
axis([0 30 0 10])  % 固定x轴的范围为0-30  固定y轴的范围为0-10
k = 0;  % 引入一个变量，为了控制画图的速度
while(dd>=0.001)  %这里的临界值要结合dt来取，否则可能导致错过交界处的情况
    t=t+dt; 
    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:2
            plot(x(i),y(i),'.k','MarkerSize',1);
            hold on; % 不关闭图形，继续画图
        end
%         pause(0.0001);  % 暂停0.0001s后再继续下面的操作
    end
    if d>50  
        disp('导弹没有击中B船');
        break;  
    end
    if d<=50 && dd<0.001   
        disp(['导弹飞行',num2str(d),'个单位后击中B船'])
        disp(['导弹飞行的时间为',num2str(t*60),'分钟'])
    end
end

结果：
在这里插入图片描述

导弹飞行27.8019个单位后击中B船
导弹飞行的时间为2.7802分钟

TSP(旅行商问题)

在这里插入图片描述
例题：

预备知识：

plot([1,2],[5,10],'-o') % 画出一条线段，x范围是[1, 2] ，y范围是[5,10]
text(1.5,7.5,'csp') % 在坐标(1.5,7.5)处标上文本：csp
close %关闭一个图形

% randperm函数的用法
randperm(5)  % 生成1-5组成的一个随机序列
%      3     5     1     2     4
%      1     4     5     3     2

求解：

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列
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);   % 初始化两个城市的距离矩阵
% 计算城市i和j的距离
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);   
    end
end
d = d+d'   % 生成距离矩阵的对称的一面

min_result = +inf;  % 初始化最短距离
min_path = [1:n];   % 初始化最短的路径就是1-2-3-...-n

N = 10000000;  
for k = 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_result
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)  
    hold on
end

结果：
在这里插入图片描述

d =

  列 1 至 6

         0    0.0498    0.3289    0.4378    0.4982    0.6710
    0.0498         0    0.2842    0.3934    0.4501    0.6323
    0.3289    0.2842         0    0.3361    0.3141    0.3601
    0.4378    0.3934    0.3361         0    0.1107    0.6149
    0.4982    0.4501    0.3141    0.1107         0    0.5349
    0.6710    0.6323    0.3601    0.6149    0.5349         0
    0.7042    0.6857    0.5111    0.8407    0.7919    0.3397
    0.4494    0.4654    0.5176    0.8083    0.8207    0.6528
    0.2690    0.2674    0.2982    0.5815    0.5941    0.5171
    0.2100    0.2492    0.4564    0.6418    0.6908    0.7375

  列 7 至 10

    0.7042    0.4494    0.2690    0.2100
    0.6857    0.4654    0.2674    0.2492
    0.5111    0.5176    0.2982    0.4564
    0.8407    0.8083    0.5815    0.6418
    0.7919    0.8207    0.5941    0.6908
    0.3397    0.6528    0.5171    0.7375
         0    0.4579    0.4529    0.6686
    0.4579         0    0.2274    0.2937
    0.4529    0.2274         0    0.2277
    0.6686    0.2937    0.2277         0


min_result =

    2.6907


min_path =

     8     9    10     1     2     4     5     3     6     7     8

最优维修方案

某设备上安装有四只型号规格完全相同的电子管，已知电子管寿命为1000–2000小时之间的均匀分布(假定为整数)，当电子管损坏时有两种维修方案，一是每次更换损坏的那一只；二是当其中一只损坏时四只同时更换，已知更换时间为换一只时需1小时，4只同时换为2小时，更换时机器因停止运转每小时的损失为20元，又每只电子管价格10元，试用模拟方法决定哪一个方案经济合理？

预备知识：

% randi([a,b],m,n)  随机生成m*n的矩阵，矩阵中的每个元素都是[a,b]中的随机整数
randi([1, 5],3,2)
randi([1, 5])  % 不写m*n代表只生成1个随机数

% find函数的用法

a = [2 3 5 1 7 5];
find(a)  % 找到a中所有非0元素的位置
find(a == 5)  % 找到a中等于5的元素的位置
find(a == 5,1)  % 找到a中第一个等于5的元素的位置
find(a == min(a))   % 找到a中最小元素的位置

求解：

clear;clc
T = 100000000;   % 模拟的总时间（单位为小时）
t = 0;   % 初始化当前时刻
c1 = 0; c2 = 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;   % 初始化
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

结果：

所以第二种方案更划算

邮轮定价问题（2015 年电工杯 B 题）

在2021国赛集训期间，我们用蒙特卡罗实现了这道题的第四问，具体如下：

针对问题四，建立航行最大收益预期模型并且计算第 8 次航行的预期收益。由于实际预定人数是由意愿预定人数转化而来，其转化率与票价有关，当票价超过了顾客的预期价格，意愿预定则不能转化为实际预定。假设预期价格服从均匀分布，通过均匀分布的分布函数找出价格影响实际预定人数的关系，进而确定实际预定人数的函数，从而构建最大收益的非线性规划模型。最后通过蒙特卡洛模拟计算出第 8 次航行的预期售票收益为 1526686

代码实现：

%蒙特卡罗模拟第 8 次航行
%每种舱位每周预定价格在价格区间[vmax,vmin]内服从均匀分布
% vmax,vmin,Mt,H
clc;
[n,m]=size(vmax);
x=0;y=0;
res=0;z2=0;res2=0;res3=0;cnt=1;
z=zeros(n,1);
rr=zeros(n,1);
for k=1:100000
    sum1=0;sum2=0;
    for i=1:n
        z(i)=randi([vmin(i),vmax(i)]); %产生在[vmax,vmin]内的随机数
        while(i>1&&abs(z(i)-z(i-1))/z(i-1)>0.2)% 20%的约束
            z(i)=randi([vmin(i),vmax(i)]);
        end
        rr(i)=randi([max(0,Mt(i)-40),Mt(i)]);
        p=ceil((vmax(i)-z(i))/(vmax(i)-vmin(i))*rr(i));
        
        sum1=sum1+p;
        sum2=sum2+p*z(i);
        if(sum1>H) %需求和 不能超过 邮轮特定舱位的上限值
            break;
        end
    end
    if(sum1>H) %需求和 不能超过 邮轮特定舱位的上限值
            continue;
    end
    if(res<sum2)
        res=sum2;
        z2=z;
        res2=sum1;
        res3=rr;
        x(cnt)=cnt;
        y(cnt)=res;
        cnt=cnt+1;
    end
end
plot(x, y, ':r*');
legend('更新值');%生成图例
xlabel('成功更新最优解次数');ylabel('收益值');
% cnt
res2
res
z2
rr

变量截图：
在这里插入图片描述
运行结果：

逻辑解释：

蒙特卡罗过程中，每种舱位每周预定价格在价格区间[vmax,vmin]内服从均匀分布
for k=1:100000即蒙特卡罗循环100000次
randi([vmin(i),vmax(i)])产生在[vmax,vmin]内的随机数while(i>1&&abs(z(i)-z(i-1))/z(i-1)>0.2)实现题干中20%的约束的要求
if(sum1>H)保证需求和不能超过邮轮特定舱位的上限值
ceil((vmax(i)-z(i))/(vmax(i)-vmin(i))*rr(i))来更新最优解

感想：

总之可以通过枚举我们想要的条件，在时间复杂度允许的范围内用蒙特卡罗进行模拟，实现对问题的求解