文章目录

二：实例代码
- rand,randi,unifrnd,randperm
- normrnd,exprnd
- plot
- 数值显示
- 随机数种子
- 2.1 圆周率随机模拟计算
- 2.2 三门问题
- 2.3 排队问题
- 2.4 非线性规划
- 2.5 TSP
Reference

二：实例代码

rand,randi,unifrnd,randperm

（1）rand(m,n)函数产生由在[0,1]之间均匀分布的随机数组成的m行n列的矩阵（或称为数组）。
（2）a + rand(m,n)(b-a) 可以输出在[a,b]之间均匀分布的随机数组成的m行n列的矩阵。
eg:-2 + rand(3,2) * (2 - (-2)) -->输出在[-2,2]之间均匀分布的随机数组成的3行2列的矩阵。
（3）a + rand(m,n)(b-a)等价于unifrnd(a,b,m,n)
（4）randi([a,b],m,n)函数可在指定区间[a,b]内随机取出大小为m*n的整数矩阵
（5）randi([a,b])可以生成[a,b]区间内的一个随机整数
（6）randperm(n)：生成1-n的随机序列

normrnd,exprnd

（1）normrnd(MU,SIGMA):生成一个服从正态分布(MU参数代表均值,SIGMA参数代表标准差）
（2）exprnd(M)表示生成一个均值为M的指数分布随机数(其对应的参数为1/M)

plot

（1）axis([xl,xr,yl,yr]) box on：画出二维坐标轴
（2）plot(x,y,’***’)：
①’r.’：画出红色小点
②’-o’：画出直线
（3）hold on：表示在相同图上作图
（4）text(x,y,’ ')：在坐标上标注文本

数值显示

(1) format long g 可以将Matlab的计算结果显示为一般的长数字格式

随机数种子

rng(seed) 使用非负整数 seed 为随机数生成函数提供种子，以使 rand、randi 和 randn 生成可预测的数字序列。

rng(‘shuffle’) 根据当前时间为随机数生成函数提供种子。rand、randi 和 randn会在每次调用 rng 时生成不同的数字序列。

2.1 圆周率随机模拟计算

问题：平面上画有等距离为a的一组平行线，向该平面投以长为l(l<a)的针，试求针与平行线相交的概率。
分析：设x为针与最近的一平行线的距离，σ表示针与平行线的夹角，0<x<L/2,0<σ<pi。若使针与平行线相交，一定有x<=L/2sin(σ)，因此：p=2L/pi*a=m/n，并通过模拟x和σ得到pi的近似解。

①使用代码：由于一次模拟的结构具有偶然性，因此可以重复多次求平均值。

clc,clear;
result = zeros(100,1);
l = 0.51;  a = 1.231;
n = 10000;
for num=1:100m = 0;  %相交个数%产生n个x phi的随机数进行模拟x = rand(1,n) * a/2;phi = rand(1,n) * pi;for i=1:nif x(i) <= 1/2 * sin(phi(i))m = m+1; endendp = m/n;mypi1 = (2*l)/(a*p);result(num) = mypi1;
end
mymeanpi = mean(result);
disp(['蒙特卡罗方法得到pi为：', num2str(mymeanpi)])

②随机运行一次的画图代码

l =  0.520;     a = 1.314;
n = 10000;
m = 0;
x = rand(1, n) * a / 2 ;
phi = rand(1, n) * pi;
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);
disp(['蒙特卡罗方法得到pi为：', num2str(mypi)])

2.2 三门问题

问题：你参加⼀档电视节⽬，节⽬组提供了ABC三扇⻔，主持⼈告诉你，其中⼀扇⻔后边有辆汽⻋，其它两扇⻔后是空的。假如你选择了B⻔，这时，主持⼈打开了C⻔，让你看到C⻔后什么都没有，然后问你要不要改选A⻔？
分析：验证改变主意获奖的概率是不改变主意获奖概率的两倍。
①三门问题可能有三种情况：不获奖、不改变选择获奖、改变选择获奖

clc,clear;
n = 100000;
a = 0; b = 0; c = 0; %分别是不改变获奖、改变获奖、不获奖的次数
for i=1:nx = randi([1,3]);  %奖品所在门编号y = randi([1,3]);  %选择的门编号change = randi([0,1]);  %是否改变if x == y   %选择的门有奖品if change == 0a = a+1;else c = c+1;endelse if change == 0c = c+1;else b = b+1;endend
end
disp(['蒙特卡罗方法得到的不改变主意时的获奖概率为：', num2str(a/n)]);
disp(['蒙特卡罗方法得到的改变主意时的获奖概率为：', num2str(b/n)]);
disp(['蒙特卡罗方法得到的没有获奖的概率为：', num2str(c/n)]);

结果：

蒙特卡罗方法得到的不改变主意时的获奖概率为：0.16643
蒙特卡罗方法得到的改变主意时的获奖概率为：0.33205
蒙特卡罗方法得到的没有获奖的概率为：0.50152

2.3 排队问题

问题：假设某银⾏⼯作时间只有⼀个服务窗⼝，⼯作⼈员只能逐个的接待顾客。当来的顾客较多时，⼀部分顾客就需要排队等待。
假设：
1：顾客到来的间隔时间服从参数为0.1的指数分布
2：每个顾客的服务时间服从均值为10，⽅差为4的正态分布(单位为分钟，若服务时间⼩于1分钟，则按1分钟计算) 3) 排队按先到先服务的规则，且不限制队伍的⻓度，每天⼯作时⻓为8⼩时。
试回答下⾯的问题：
1：模拟⼀个⼯作⽇，在这个⼯作⽇共接待了多
少客户，客户平均等待的时间为多少?
2：模拟100个⼯作⽇，计算出平均每⽇接待客户的个数以及每⽇客户的平均等待时⻓。
符号：
ci：第i个客户到达的时间
bi：第i个客户开始服务的时间
ei：第i个客户结束服务的时间
xi：第i个客户和第i-1个客户到达的间隔时间
yi：第i个客户的服务时间
waiti：第i个客户的等待时间
分析：
1)c(i) = c(i-1)+xi (c0=0)
2)b(i) = max(e(i-1),c(i))
3)e(i) = b(i)+y(i) (e0=0)
4)wait(i) = b(i) - c(i)
①随机模拟一天的排队情况

clc,clear;
tic
i = 1;
w = 0; %总等待时间
c0 = 0;    e0 = 0;
x(1) = exprnd(10); %M=0.1
c(1) = c0+x(1);
b(1) = max(e0,c(1));
while b(i) <= 480   %一天工作8小时y(i) = normrnd(10,2);  %正太随机服务时间if y(i) < 1y(i) = 1;ende(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 = i-1;
t = w/n;   %客户的平均等待时间
disp(['银行一天8小时一共服务的客户人数为: ',num2str(n)])
disp(['客户的平均等待时间为: ',num2str(t)])
toc

②随机模拟day天计算平均等待时长与平均接待人数

clc,clear;
day = 10000;
n = zeros(100,1);
t = zeros(100,1);
for k=1:dayi = 1;w = 0;  %总等待时间c0 = 0;  e0 = 0;x(1) = exprnd(10); %M=0.1c(1) = c0+x(1);b(1) = max(e0,c(1));while b(i) <= 480      %一天工作8小时y(i) = normrnd(10,2);  %正太随机服务时间if y(i) < 1y(i) = 1;ende(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));endn(k) = i-1;t(k) = w/n(k);
end
disp([num2str(day),'个工作日中，银行每日平均服务的客户人数为: ',num2str(mean(n))])
disp([num2str(day),'个工作日中，银行每日客户的平均等待时间为: ',num2str(mean(t))])

结果：

银行一天8小时一共服务的客户人数为: 49
客户的平均等待时间为: 35.177
时间已过 0.009927 秒。
100个工作日中，银行每日平均服务的客户人数为: 43.49
100个工作日中，银行每日客户的平均等待时间为: 30.2145
时间已过 0.052765 秒。

2.4 非线性规划

问题：利用蒙特卡洛模拟算法解非线性规划的近似解，并用作非线性规划matlab方法的初值得到精确解。
思路：
①从约束条件中求出每个变量的大致范围
②在上述范围中随机生成若干试验点，并验证是否满足所有约束条件。若满足则分到可行集合中，并在其中找到最大值或最小值。
例题1：
max f = x1x2x3
s.t. -x1+2x2+2x3 >= 0
x1+2x2+2x3 <= 72
10 <= x2 <=20
x1-x2 = 10
步骤一：等式化简求变量的大致范围
x2 ∈ [ 10 , 20 ]
x1 = x2+10 -->x1 ∈ [ 20 , 30 ]
x3 >= x1-2x2/2 --> x3 >= -10
x3 <= 72-x1-2x2/2 --> x3 <= 16
可得：x3 ∈ [ -10 , 16 ]
步骤二：在范围内随机生成数据，验证约束条件，并在满足条件的集合中找到最优值
unifrnd(a,b,n,m)

clc,clear;
tic
n = 1000000;
&生成范围内随机数
x1 = unifrnd(20,30,n,1);
x2 = x1-10;
x3 = unifrnd(-10,16,n,1);
fmax = -inf;
for i=1:nx = [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 > fmaxfmax = result;X = x;endend
end
disp(strcat('蒙特卡罗模拟得到的最大值为',num2str(fmax)))
disp('最大值处x1 x2 x3的取值为：')
disp(X)
toc

结果：

蒙特卡罗模拟得到的最大值为3445.6022
最大值处x1 x2 x3的取值为：
22.5897459433039 12.5897459433039 12.1153738601129

通过蒙特卡洛模拟得到的近似解可用于matlab方法的初值设定，并可以通过缩小范围重新获得更加精确的取值。

例题2：
min f = 2x1^2 + x2^2 - x1x2 - 8x1 - 3x2
s.t. 3x1+x2 > 9
x1+2x2 < 16
x1>0 x2>0
步骤一：等式化简求变量的大致范围
x1+2x2 < 16 --> x1 < 16-2x2 --> x1 < 16
x2 < 16-x2/2 --> x2 < 8
3x1+x2 > 9 --> x2 > 9-3*x1 --> x2 > -39
x1 > (9-x2)/3 --> x1 > 1/3
由此可得： x1 ∈ [0,16],x2 ∈ [0,8]
步骤二：在范围内随机生成数据，验证约束条件，并在满足条件的集合中找到最优值

clc,clear;
format long g
tic
n = 1000000;
%生成范围内随机数
x1 = unifrnd(0,16,n,1);
x2 = unifrnd(0,8,n,1);
fmin = inf;
for i=1:nx = [x1(i),x2(i)];if (3*x(1)+x(2) > 9) & (x(1)+2*x(2) < 16)result = 2*x(1)*x(1)+x(2)^2-x(1)*x(2)-8*x(1)-3*x(2);if result < fminfmin = result;X = x;endend
end
disp(strcat('蒙特卡罗模拟得到的最小值为',num2str(fmin)))
disp('最小值处x1 x2的取值为：')
disp(X)
toc

2.5 TSP

问题：一个售货员访问n个城市，折n个城市是一个完全图，售货员需要恰好经过所有城市一次，并且回到最初的城市，且希望旅行的路程最短。
分析：TSP是个NP难问题，若采用枚举法，则会有O((n-1)!)的时间复杂度。因此蒙特卡洛模拟法不一定可以得到精确解，且对于较大n的运行时间较长而无法适用。
思路：随机生成1-n的序列即售货员走过城市的顺序，计算旅行的总路程。若小于当前最短路程则替换。

clc,clear;
%%1.记录节点坐标，并在作图标注
coord = [0.6683 0.2536;0.6195 0.2634;0.4 0.4439;0.2439 0.1463;0.1707 0.2293;0.2293 0.7810;0.5171 0.9414;0.8732 0.6536;0.6878 0.5219;0.8488 0.3609];
n = size(coord,1);
figure(1)
plot(coord(:,1),coord(:,2),'o');
for i=1:ntext(coord(i,1)+0.01,coord(i,2)+0.01,num2str(i));
end
hold on
%%2.计算结点间距离
d = zeros(n,n);
for i=1:nfor j=i+1:ncoord_i = coord(i,:); coord_j = coord(j,:);x_i = coord_i(1); y_i = coord_i(2);x_j = coord_j(1); y_j = coord_j(2);d(i,j) = sqrt((x_i-x_j)^2+(y_i-y_j)^2);end
end
d = d+d';
%%3.生成随机序列，模拟计算最短路径
min_result = inf;
min_path = [1:n];
N = 10000000;
for i=1:Nresult = 0;path = randperm(n);for j=1:n-1result = result+d(path(j),path(j+1));endresult = result+d(path(n),path(1));if result < min_resultmin_result = resultmin_path = path;end
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

结果：

min_path
5 4 2 1 10 9 8 7 6 3

Reference

清风数学建模：
https://www.bilibili.com/video/BV1DW411s7wi

国赛培训——随机化算法——蒙特卡洛模拟相关推荐

计算机检测维修与数据恢复国赛培训线上直播课火热进行中
红警(中国)维修连锁,已于2020年2月份正式上线"网上教学"模块,利用"互联网+教育"进行直播教学. 学生在线学习线上直播课一.培训目标通过理论与实战 ...
最新2022年高职大数据国赛任务书详解与模拟练习
2022高职大数据竞赛模拟练习-模拟数据说明 2022高职大数据竞赛模拟练习-离线数据处理任务一:数据抽取
国赛培训——评价模型——TOSIS法
文章目录例题 (1)原始矩阵正向化 (2)标准化处理 (3)计算得分并归一化代码 Reference 例题题目:评价下表中20条河流的水质情况注:含氧量越高越好:PH值越接近7越好:细菌总数越 ...
2022数学建模国赛备赛阶段性记录（1-1）
数学建模国赛培训记录,主要使用软件为MATLAB,主要内容为在数学建模竞赛中常用的操作.数学与模型以及部分练习题的解析. 一.常规操作 1.基本运算 MATLAB内四则运算相当于计算机的加减乘除,对应 ...
数学建模竞赛常用算法介绍及对应国赛获奖论文分类整理分享
数学建模竞赛中应当掌握的算法: 数学建模国赛每年的题型都类似,除非是个人专业性很强,否则作者不太建议选择华为出的题,剩余的题型每年都类似,是有迹可循的,毕竟站在巨人的肩膀上看的更远.下面就介绍一些数模 ...
备战数学建模国赛，快速搞定算法模型！
全世界只有3.14 % 的人关注了青少年数学之旅说到数学建模,大家的第一反应就是国赛.美赛等数学建模比赛,但这只是冰山一角,不过这个反应却也很正常,因为很多小伙伴接触数学建模的契机,大部分还是因为 ...
2022国赛数学建模思路 - 案例：集成算法AdaBoost
2022 高教社杯(国赛数学建模)思路解析 2022高教社杯ABCD赛题思路解析: https://blog.csdn.net/dc_sinor/article/details/126211983 集 ...
2020年第四届计算机检测维修与数据恢复国赛模拟比赛
第四届计算机检测维修与数据恢复模拟比赛在我司实训车间举行,本次模拟比赛按照 2019-2020年度广东省职业院校技能大赛(中职组)"计算机检测维修与数据恢复"项目样题任务书进行出 ...
Algorithm：数学建模大赛(国赛和美赛)的简介/内容、数学建模做题流程、历年题目类型及思想、常用算法、常用工具之详细攻略
Algorithm:数学建模大赛(国赛和美赛)的简介/内容.数学建模做题流程.历年题目类型及思想.常用算法.常用工具之详细攻略目录国内数学建模大赛简介 1.本科生数学建模大赛 2.研究生数学建模大 ...

国赛培训——随机化算法——蒙特卡洛模拟