【学习笔记】MATLAB与数学建模——蒙特卡罗模拟&仿真

目录

引例：投针试验

试验描述：

试验分析：

代码实现：

蒙特卡罗模拟&仿真

预备知识

蒙特卡罗模拟应用实例

实例一、三门问题

问题描述

问题分析与代码实现

实例二、排队问题1-港口卸货

问题描述

问题分析与代码实现

实例三、排队问题2-银行排队

问题描述

问题分析与代码实现

实例四、有约束的非线性规划

问题描述

问题分析与代码实现

实例五、书店选择(0-1规划)

问题描述

问题分析与代码实现

实例六 、导弹追踪

问题描述

问题分析与代码实现

实例七、旅行商问题

问题描述

问题分析与代码实现

实例八、加油站存储策略

问题描述

问题分析与代码实现

---

 

引例：投针试验

 

试验描述：

(1). 取一张白纸，在上面画出许多条间距为a的平行线；

(2). 取一根长度为L(L<a)的针，随机地向纸上丢掷n次，记录针与直线相交的次数，记为m；

(3). 计算得到针与直线相交的概率(大数定律)

试验分析：

  如图所示，Q为针的中点，x为中点Q到最近一条平行线的距离，为针与这条平行线的夹角，则显然有如下关系：

而如果针与直线相交，则满足：

图示为：

 

 根据以上关系，经过证明可得针与直线相交的概率即为蓝色区域面积与矩形面积的比值，经过积分计算得到: （本文不作计算）,即，由此我们可以用投针试验来估计的值：

代码实现：

L=1; %设置针的长度
a=2; %设置平行线之间的距离
n=10000; %设置单次模拟的试验次数，次数越多，求解的PI越准确
m=0; %用来记录针与平行线相交的次数
N=1000; %设置模拟的次数，避免偶然性
x=rand(1,n)*a/2; %在0到a/2上服从均匀分布随机产生n个数，表示针的中点与最近平行线的距离
Angle=rand(1,n)*pi %在0到PI上服从均匀分布随机产生n个数，表示针与最近平行线的夹角
%开始进行模拟
result=zeros(1,N); %结果矩阵储存每次模拟结果
for j=1:N %设置多次模拟，避免偶然性
    m=0;
    p=0; %每次循环要初始化
    for i=1:n %开始投掷
        if x(i) <= L/2 * sin(Angle(i)) %如果针与平行线相交
            m = m + 1 %就计数
        end
    end
    p=m/n; %针与平行线相交的频率
    Single_Cal_PI=(2*L)/(a*p); %该次模拟得出的PI值
    result(j)=Single_Cal_PI; %将模拟结果储存到结果矩阵中
end
Cal_PI=mean(result); %将结果矩阵中记录的结果求均值
disp(['经过模拟得到的PI值为 ',num2str(Cal_PI)]);
 %输出结果，num2str函数将数值转化为字符串

蒙特卡罗模拟&仿真

蒙特卡罗模拟是一种随机模拟方法，如果我们所求解的问题与概率模型存在一定的关系，我们便可以借助计算机多次模拟事件的发生，以获得问题的近似解。

蒙特卡罗模拟不能说是一种方法，而是一种思想，因此针对不同的问题我们要设计不同的代码。

预备知识

预备知识包括一些可能不熟悉，但会用到的函数

%%预备知识
randi([low,up],m,n);
%生成一个m*n阶矩阵，每个元素是取值范围为[low up]的整数

randi([low,up])
%生成取值范围为[low up]的整数

['字符串1','字符串2'];
strcat('字符串1','字符串2')；
%两种字符串连接方式

num2str();
%将数值转化为字符串

disp('字符串');
%输出函数

normrnd(MIU,SIGMA);
%生成一个服从正态分布(MIU参数代表均值,SIGMA参数代表标准差)的随机数

exprnd(M);
%生成一个均值为λ=M的指数分布随机数(其对应的参数为θ=1/M)

mean();
%用来求解均值的函数

tic函数和toc函数;
%用来返回代码运行的时间
%例如我们要计算一段代码的运行时间，就可以在这段代码前加上tic，在这段代码后加上toc 

format long g;
%设置计算结果为长数字型

unifrnd(low,up,m,n);
%生成一个m*n阶矩阵，矩阵元素取值满足[low up]的均匀分布

unique()；
%剔除一个矩阵或者向量的重复值，并将结果按照从小到大的顺序排列

mod(a,b);
%求a/b的余数

randperm(k);  
%生成1-k组成的一个随机序列

蒙特卡罗模拟应用实例

实例一、三门问题

问题描述

你参加了一个电视节目，面前有三道门，其中一扇门后有一个大奖，此时你选择了B门，若此时主持人打开C门，发现为空，问你是否要更换为A门？

问题分析与代码实现

%实例一、三门问题
N=1000000; %模拟次数
Ch=0; %记录改变选项赢得奖品次数
NoCh=0; %记录不改变选项赢得奖品次数

for i=1:N
    x=randi([1 3]); %随机生成1到3的整数，表示奖品所在位置
    y=randi([1 3]); %表示初始选择
    if x==y
        NoCh=NoCh+1;
    else
        Ch=Ch+1;
    end
end

barh(Ch,'red');
hold on;
barh(NoCh,'blue');
title('两种抉择后赢的次数');
legend('改变选择后赢的次数','不改变选择后赢的次数')

结果可视化如下：

 可以明显地看到，改变选择后赢的次数远远大于不改变的，实际上，改变后赢的可能性为不改变的2倍。

实例二、排队问题1-港口卸货

问题描述

只有一个港口可供卸货，每两艘船到达港口的时间间隔为15分钟到145分钟，轮船卸货时间取决于货物量，时间长为45分钟到90分钟。

通过模拟仿真，回答如下问题：

1、每艘船在港口呆的平均时间；

2、每艘船等待的平均时间；

3、港口的空闲时间比例；

4、排队时最长的队伍。

问题分析与代码实现

假设以1分钟为最小计时单位。

假设到达服从一致分布(可以根据不同情况调整到达时间的分布)

代码实现：

%实例二、排队问题1-港口卸货
clear;clc;
%随机生成数据(一致分布)
N=5; %设置总船只数，这里我设置为5，可更改
Between=randi([15 145],1,N); %生成均匀一致分布的到达时间间隔
Unload=randi([45 90],1,N); %生成随机均匀一致分布的卸货时间长度
Arrive=cumsum(Between); %计算船只到达时间

Start=zeros(1,N); %初始化开始卸货时间点
Idle=zeros(1,N); %初始化港口空闲时间长度
Wait=zeros(1,N); %初始化等待时间长度
Finish=zeros(1,N); %初始化卸货完毕时间点
Total=zeros(1,N); %初始化轮船在港口总用时长度

%第一艘船建模计算
Idle(1)=Arrive(1); %第一艘船未到达时间为港口空闲时间
Wait(1)=0; %等待时间
Start(1)=Arrive(1)+Wait(1); %第一艘开始卸货时间
Finish(1)=Start(1)+Unload(1); %第一艘结束卸货时间
Total(1)=Finish(1)-Arrive(1); %第一艘在港口的总时间

%对剩下的船只进行模拟仿真
for ind=2:N
    Idle(ind)=max(0,Arrive(ind)-Finish(ind-1)); %第ind艘船到达前的空闲时间，要大于等于0
    Wait(ind)=max(0,Finish(ind-1)-Arrive(ind)); %第ind艘船到达后的等待时间，要大于等于0
    Start(ind)=Arrive(ind)+Wait(ind); %第ind艘船开始卸货的时间
    Finish(ind)=Start(ind)+Unload(ind); %第ind艘船完成卸货的时间
    Total(ind)=Finish(ind)-Arrive(ind); %第ind艘船在港口的总时间
    %可视化
    for t=1:10:Finish(ind)-Arrive(ind)
        bar(Finish,'green');
        hold on;
        bar(Start,'red');
        bar(Arrive,'White');
        legend('卸货时间','等待时间','到达时间','Location','Northeastoutside');
        axis([0,N+1,Arrive(1),Finish(ind-1)+t]);
        grid on;
        xlabel('船只编号');
        ylabel('时间经历');
        hold off;
        drawnow;
    end
end

%计算港口排队系统的性能指标
Total_mean=mean(Total); %平均在港口时间
Wait_mean=mean(Wait); %平均等待时间
Idle_percent=sum(Idle)./Finish(end) %空闲率

结果可视化如下：

 在实际运行过程中，也可以看到动态的变化过程。

 

实例三、排队问题2-银行排队

问题描述

银行只有一个服务窗口，当顾客较多，一部分顾客就需要排队等待。

假设：顾客到来的时间服从参数为0.1的指数分布；

每个顾客的服务时间服从均值为10，方差为4的正态分布；

排队按照先到先服务原则，每天工作八小时。

建立模型，求解客户平均等待时长为多少？

问题分析与代码实现

该问题与实例二——港口卸货非常类似，区别在于时间分布不同，代码实现如下：

%实例三、排队问题2-银行排队问题
clear;clc;
N=100; %取N天为样本，即模拟100次
Total=zeros(1,N); %所有顾客等待时间和

for j=1:1:N
    i=1; %表示第i个客户，从1开始
    Arrive(1)=exprnd(10); %到达时间间隔，服从参数为0.1的指数分布
    Wait(1)=0; %等待时间
    Start(1)=Arrive(1)+Wait(1); %开始办理业务时间
    Process(1)=normrnd(10,2) %办理业务时长，服从均值为10，方差为4(标准差为2)的正态分布
    End(1)=Start(1)+Process(1); %结束时间等于开始时间加上办理业务时间
    
    while Start(i)<=480
        i=i+1;
        Arrive(i)=Arrive(i-1)+exprnd(10); %第i个客户到达时间
        Start(i)=max(Arrive(i),End(i-1)); %第i个客户开始时间
        Wait(i)=max(0,End(i-1)-Arrive(i)); %第i个客户等待的时间
        Process(i)=normrnd(10,2); %第i个客户办理业务的时间
        End(i)=Start(i)+Process(i); %第i个客户结束的时间
        Total(j)=Total(j)+Wait(i); %等待的总时长
        Average_Wait(j)=(Total(j))/i; %平均每个客户等待时长
    end

        Final_Average=mean(Average_Wait); %N次模拟得出的总平均等待时长
end

%可视化
x=1:1:N;
plot(x,Average_Wait,'-.r.');
xlabel('天数');
ylabel('顾客平均排队时长(分钟)');
title('每天顾客平均等待时长');

结果可视化如下:

 

实例四、有约束的非线性规划

问题描述

求解如下一个有约束的非线性规划问题：

问题分析与代码实现

对于求解该类型问题，我们先通过约束条件确定每个自变量的大概取值范围，在这些范围中随机生成若组试验点，若小组整体都满足约束条件，则划分到可行组，并从中找到函数的最值。

%实例四、有约束的非线性规划问题
% max f = x1*x2*x3
% s.t.
% -x1+2*x2+2*x3>=0
% x1+2*x2+2*x3<=72
% x2<=20 & x2>=10
% x1-x2 == 10
%有约束条件可以得到x1,x2,x3的大概取值范围
%其中x3满足: 5-x2/2<=x3<=31-3*x2/2 ，因此取区间[-5 16]

clear;clc;
tic
format long g; %设置输出数值类型
N=1000000; %N组数据点
x2=unifrnd(10,20,1,N); %均匀分布
x1=x2+10;
x3=unifrnd(-5,16,1,N);
f=-inf;

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>f
            f=result;
            final_X=x;
        end
    end
end

disp(['经过蒙特卡罗模拟后，得出的最大目标函数值为：',num2str(f)]);
disp(['最终自变量取值为：x1=',num2str(final_X(1)),' x2=',num2str(final_X(2)),' x3=',num2str(final_X3))]);
toc

输出结果为：

实例五、书店选择(0-1规划)

问题描述

某同学要从六家线上商城选购五本书籍B1,B2,B3,B4,B5，每本书在不同商家的售价以及每个商家的单次运费如下表所示，请给该同学制定最省钱的选购方案。

	B1	B2	B3	B4	B5	运费
A商城	18	39	29	48	59	10
B商城	24	45	23	54	44	15
C商城	22	45	23	53	53	15
D商城	28	47	17	57	47	10
E商城	24	42	24	47	59	10
F商城	27	48	20	55	53	15

问题分析与代码实现

%实例五、书店选择(0-1规划)
clear;clc;
Final_Shop=randi([1 6],1,5); %初始化每本书对应选择的书店
Final_Cost=inf; %初始化花费
N=10000; %模拟的次数
Price=[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]; %每个店铺各个书的售价
Trans=[10 15 15 10 10 15]; %每个店铺的运费

for k = 1:N  % 开始循环模拟
    Shop = randi([1, 6],1,5); %重新随机每本书对应选择的书店
    index = unique(Shop);  %去除重复元素并排序，确定有哪些书店
    Cost = sum(Trans(index)); %计算买书花费的运费
    %计算总花费=运费+五本书的售价
    for i=1:5   
        Cost = Cost + Price(Shop(i),i); %每本书对应的花费 
    end
    
    if Cost < Final_Cost  % 判断刚刚随机生成的这组数据的花费是否小于最小花费，如果小于的话
        Final_Cost = Cost;  % 我们更新最小的花费
        Final_Shop = Shop;% 用这组数据更新最小花费的结果
    end
end

disp('经过蒙特卡洛模拟：')
disp(['这五本书分别在这些书店购买：书店',num2str(Final_Shop(1)),' 书店',num2str(Final_Shop(2)),' 书店',num2str(Final_Shop(3)),' 书店',num2str(Final_Shop(4)),' 书店',num2str(Final_Shop(5))]);
disp(['总共花费',num2str(Final_Cost),'元']);

输出结果为：

实例六 、导弹追踪

问题描述

位于坐标原点的A船向位于其正东方20个单位的B船发射导弹，导弹始终对准B船，B船以时速v单位沿东北方向逃逸。若导弹的速度为3v，导弹的射程为50，画出导弹运行的曲线，导弹是否能在射程内集中B船？

问题分析与代码实现

本问题在使用仿真时，实际上是微元法的体现，将连续的时间均等分割为大量的离散时间，每一段时间都尽可能短。

时间 t 时的情形如下(P时导弹的坐标)：

代码实现为：

%实例六、导弹追踪
clear;clc
tic
v=200; %任意给定v的值
dt=1e-8; %定义时间间隔
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定义为一个常量，更加简洁
Distance=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 导弹与B船的距离

plot(x(1),y(1),'r.','MarkerSize',0.5);
grid on;
hold on;
plot(x(2),y(2),'b.','MarkerSize',1);
hold on;

axis([0 30 0 10])  %固定x轴的范围为0-30  固定y轴的范围为0-10
k=0;  %引入一个变量  为了控制画图的速度（因为Matlab中画图的速度超级慢）

while(Distance>=1e-5)  %设置击中的临界距离
    t=t+dt; %更新导弹击落B船的时间
    d=d+3*v*dt; %更新导弹飞行的距离
    x(2)=20+t*v*m;  y(2)=t*v*m; %计算B船的新位置
    Distance=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); %更新导弹与B船的距离
    tan_alpha=(y(2)-y(1))/(x(2)-x(1)); %根据导弹与B的位置计算导弹的方向，即tan(α)
    cos_alpha=(x(2)-x(1))/Distance;   
    sin_alpha=(y(2)-y(1))/Distance;  
    x(1)=x(1)+3*v*dt*cos_alpha;   y(1)=y(1)+3*v*dt*sin_alpha; %计算导弹的新位置
    
    k = k +1 ;  
    
    if mod(k,2000) == 0   % 每刷新1000次时间就画出下一个导弹和B船所在的坐标  mod(m,n)表示求m/n的余数
        plot(x(1),y(1),'r.','MarkerSize',0.5);
        hold on;
        plot(x(2),y(2),'b.','MarkerSize',1);
        hold on;
        pause(0.001);  % 暂停0.001s后再继续下面的操作
    end
    
    if d>50  % 导弹的有效射程为50个单位
        disp('导弹没有击中B船');
        break;  % 退出循环
    end
    
    if d<=50 && Distance<1e-5  % 导弹飞行的距离小于50个单位且导弹和B船的距离小于临界距离
        disp(['导弹飞行',num2str(d),'单位后击中B船'])
        disp(['导弹飞行的时间为',num2str(t*60),'分钟'])
    end
end

legend('导弹运行轨迹','B运行轨迹','Location','NortheastOutside');
toc

实际可以看到动态图，最终结果可视化如下：

 输出结果为：

 

实例七、旅行商问题

问题描述

一个售货员必须访问10个城市，这10个城市可用完全图表示，售货员需要恰好访问所有城市一次，并且回到最终的城市。售货员希望旅行距离之和最少。

完全图：一种简单的无向图，其中每对不同的顶点之间都恰有一条边相连

问题分析与代码实现

%实例七、旅行商问题
clear;clc;
%我们这里随机生成10个城市的坐标
Position=unifrnd(1,25,10,2);
N = size(Position,1);  % 城市的数目
plot(Position(:,1),Position(:,2),'.','Markersize',12); %画出城市的分布散点图

for i = 1:N
    text(Position(i,1)+0.3,Position(i,2),num2str(i)) %在图上标上城市的编号
end
hold on;

Distance = zeros(N); %初始化两个城市的距离矩阵全为0
for i = 2:N  
    for j = 1:i  
        Position_i = Position(i,:); %取出第i个城市的坐标
        x_i = Position_i(1);     y_i = Position_i(2); %城市i的横坐标为x_i，纵坐标为y_i
        Position_j = Position(j,:); %取出第j个城市的坐标
        x_j = Position_j(1);     y_j = Position_j(2); %城市j的横坐标为x_j，纵坐标为y_j
        Distance(i,j) = sqrt((x_i-x_j)^2 + (y_i-y_j)^2);   %计算从城市i到城市j的距离
    end
end
%此时Distance是一个下三角矩阵，表示的是城市i到城市j的距离
%又因为城市i到j的距离 等于 城市j到i的距离
Distance = Distance+Distance'; %生成距离矩阵的对称的一面

Min_Cost = inf; %初始化最短距离
Min_Path = [1:N]; %初始化最短的路径：1-2-3-...-n

Num = 10000000;  %蒙特卡罗模拟的次数
for i=1:Num  
    Cost = 0; %初始化走过的路程
    Path = randperm(N); %将1到n的序列重新随机
    
    for j=1:N-1  
        Cost = Distance(Path(j),Path(j+1)) + Cost; %按照序列求得距离之和
    end
    Cost = Distance(Path(1),Path(N)) + Cost; %加上从最后城市返回到最初城市的距离
    
    if Cost < Min_Cost  %判断是否是更优解，如果是，则更新
        Min_Path = Path;
        Min_Cost = Cost;
    end
end

Min_Path = [Min_Path,Min_Path(1)]; %在最短路线后补充上最初的城市
N = N+1; %补充了最初的城市

for i = 1:N-1 
     j = i+1;
    %按顺序取出城市路线
    Position_i=Position(Min_Path(i),:);
    x_i = Position_i(1);     y_i = Position_i(2); 
    Position_j=Position(Min_Path(j),:);
    x_j = Position_j(1);     y_j = Position_j(2);
   
    plot([x_i,x_j],[y_i,y_j],'--') %按照顺序依次连线，最终显示出路径
    pause(0.5)  % 暂停0.5s再画下一条线段
    hold on
end

结果可视化如下：

输出结果为：

 

实例八、加油站存储策略

问题描述

你作为加油站的顾问，需要为加油站设计补货策略：补充频率和补充量。已知每次补货的运输费用：d(花费固定，与油量无关)，还需要考虑油的存储费用(费用与油量相关)。如何设计方案，使得加油站成本最小，同时能够满足顾客的需求？

每日油量需求的历史数据：

需要的油量	频次	频率
1000——1099	10	0.01
1100——1199	20	0.02
1200——1299	50	0.02
1300——1399	120	0.12
1400——1499	200	0.2
1500——1599	270	0.27
1600——1699	180	0.18
1700——1799	80	0.08
1800——1899	40	0.04
1900——1999	30	0.03
总天数	1000	1

问题分析与代码实现

根据经济批量补货公式：

最佳补货周期：

最佳补货量：

r：每天的需求量

d：每次补货的运输费用

s：每天单位油量的花费

首先，根据题目中油量需求的历史数据，用频率估计概率，并结合均匀分布与插值算法，模拟估计出每日需求量(本文省略计算过程，只体现最后结果)：

%%模拟每日需求量
function y=gasoline_demand_function(x) %自定义函数
%x=rand(100，1); %提醒：其中可以放向量

y=(x<0.01).*(x+0.2)*5000....
+(x>=0.01&x<0.03).*(x+0.2)*5000....
+(x>=0.03&x<0.08).*(x+0.545)*2000....
+(x>=0.08&x<0.20).*(x+1.42)*833.33....
+(x>=0.20&x<0.40).*(x+2.5)*500....
+(x>=0.40&x<0.67).*(x+3.515)*370.37....
+(x>=0.67&x<0.85).*(x+2.12)*555.55....
+(x>=0.85&x<0.93).*(x+0.472)*1250....
+(x>=0.93&x<0.97).*(x-0.23)*2500....
+(x>=0.97).*(x-0.6)*5000
%当代码过长，一行写不下时，在每行结束加上 ' .... '
end

模拟过程：

%%补货模型
format short; %设置输出格式 5字长定点数
N=30; %模拟天数
I=zeros(1,N+1); %每日结束后的库存量
D=zeros(1,N); %每天油的运送补充量
d=1000; %设置每次运费
s=0.1; %设置每天单位油量的存储费用
C=zeros(1,N); %每天结束后总的成本
q=zeros(1,N); %每天油的需求量

%模拟主体
q=gasoline_demand_function(rand(N,1)); %随机每日需求量
r=mean(q); %平均需求量
T=round(sqrt(2*d/(s*r))); %最优补货时间 round()函数 表示四舍五入
Q=r*T; %最优补货量
D(1:T:end)=Q; %即每T天，补充量为Q，其余为0
C(1:T:end)=d; %即每T天，都会有一次运费
for ind=1:N
    I(ind+1)=max(I(ind)+D(ind)-q(ind),0); %计算每天结束后的库存量，满足大于等于0
    C(ind)=C(ind)+I(ind+1)*s; %加上每天的存储成本
end

%可视化
plot(1:N,q','--',1:N,I(2:end),'o-',1:N,C,'-','LineWidth',2);
%画出每日需求量、每日库存量、总成本
legend('每日需求量','每日库存量','总成本') %图例

结果可视化如下：

 问题解决啦！