蒙特卡洛方法是一种很常用的方法,笔者在学习该方法时寻找相关资料做了一些学习笔记,方便之后查阅。笔记通过一些蒙特卡洛方法的实例,用MATLAB编程解决这些实际问题来学习蒙特卡洛方法。
蒙特卡洛模型
一、蒙特卡洛模拟
1.蒙特卡洛方法概述
定义:蒙特卡洛方法也称统计模拟法,是一种随机模拟方法,以概率统计理论为基础,使用(伪)随机数解决计算问题。将所求问题与一定的概率模型相联系,用计算机实现统计模拟或抽样,以获得问题的近似解。
名称由来:20世纪40年代J.冯·诺伊曼、S.M.乌拉姆等人在“曼哈顿计划”中研究核武器工作中涉及多次实验结果的分析,统计分析的概率统计特征与赌博的统计特性之间有密切联系,因此借用摩纳哥的赌城蒙特卡洛命名。1777年蒲丰提出投针求圆周率的方法则被认为是蒙特卡洛方法的起源。
原理:大数定律,当重复次数足够多时,频率近似等于概率。
蒙特卡洛不是一个算法,而是一个方法,没有通用的代码,需要结合问题抽象出概率模型进行应用。蒙特卡洛方法可以看作枚举法的一个变异,当问题的解有无数多个或非常多时,即使借助现有的计算机也无法枚举出所有的结果,这时就需要用蒙特卡洛方法将一个不确定的问题转化成很多确定的问题并得到一个近似的解。
2. 试验与实例
2.1 蒲丰投针试验
18世纪,法国数学家蒲丰 (buffon,另译“布丰”) 最早设计了投针试验,记载于其1777年出版的著作中:“在平面上画有一组间距为
a
a
a的平行线,将一根长为
l
(
l
≤
a
)
l(l\leq a)
l(l≤a)的针任意掷在这个平面上,求此针与平行线中任意一条相交的概率。”
“蒲丰投针问题”是概率论中一个经典的几何概型问题,这同时是一种计算圆周率的方法——随机投针法。
将随机投针法计算圆周率
π
\pi
π的步骤概括如下:
(1) 在一张白纸上画很多条等距且间距为
a
a
a的平行线
(2) 将长度为 l ( l ≤ a ) l(l\leq a) l(l≤a)的针随机投掷在这张白纸上,设共投掷了n次,其中观察到有m次针与直线相交
(3) 计算针与直线相交的概率
蒲丰自己已经证明,这个概率是
p
=
2
l
π
a
p=\frac{2l}{\pi a}
p=πa2l.根据概率的定义,可以用频率估计概率。针与直线相交的频率为
m
n
\frac{m}{n}
nm,当试验的次数足够大时,频率可以近似等于概率,从而求出
π
\pi
π.
历史上一些科学家用该方法计算的圆周率:
理论上试验的次数越大精度越高,但实际上有一些偏差,能看出这种求
π
\pi
π的方法具有很大的随机性。
相交概率的证明
用几何概型来解决该问题,可以从面积、线段两个角度入手,设针与直线相交的概率为
p
p
p
方法一:面积法
将针与直线的夹角设为
θ
\theta
θ,中点
M
M
M与最近一条直线之间的距离为
x
x
x,当
x
<
l
2
sin
θ
x<\frac{l}{2}\sin\theta
x<2lsinθ 时,针与直线相交。
以
θ
\theta
θ为横坐标,
x
x
x为纵坐标,
θ
\theta
θ的范围是
0
≤
θ
≤
π
0\leq\theta\leq\pi
0≤θ≤π,
x
x
x的范围是
0
≤
x
≤
a
2
0\leq x\leq \frac{a}{2}
0≤x≤2a.针与直线相交的概率是阴影部分的面积除以整个矩形外框的面积:
p
=
∫
0
π
1
2
l
sin
θ
d
θ
1
2
π
a
=
2
l
π
a
p= \frac{\int_0^\pi \frac{1}{2}l\sin\theta d\theta}{\frac{1}{2}\pi a}=\frac{2l}{\pi a}
p=21πa∫0π21lsinθdθ=πa2l
方法二:线段法
当针与直线的夹角为
θ
\theta
θ时,将针AB投影在垂直于直线的
y
y
y方向,想象针的上端点A可以在线段PQ上移动,则只有在红色区域内移动时才能与直线相交。夹角为
θ
\theta
θ的针与直线相交的概率为
p
(
θ
)
=
l
sin
θ
a
p(\theta)= \frac{l\sin\theta}{a}
p(θ)=alsinθ
θ
\theta
θ的范围是
0
≤
θ
≤
π
0\leq\theta\leq\pi
0≤θ≤π,所有针的夹角在这个范围内等可能分布,总体的概率等价于对
θ
\theta
θ进行加权平均
p
=
∫
0
π
1
a
l
sin
θ
d
θ
π
−
0
=
2
l
π
a
p= \frac{\int_0^\pi \frac{1}{a}l\sin\theta d\theta}{\pi -0}=\frac{2l}{\pi a}
p=π−0∫0πa1lsinθdθ=πa2l
两种方法得到的结果是相同的,针与直线相交的概率为:
p
=
2
l
π
a
p= \frac{2l}{\pi a}
p=πa2l
蒲丰投针求 π \pi π编程实现(MATLAB)
编程需要用到的MATLAB语法
输出分布在[a,b]之间均匀分布的随机数组成的m行n列矩阵:
//两种方法,替换a,b,m,n为数字即可
a + rand(m,n) * (b - a)
unifrnd(a,b,m,n)
蒲丰投针——蒙特卡洛模拟代码
单次求 π \pi π以及绘图的代码:
figure
l = 1; % 针的长度(任意给)
a = 1.5; % 平行线间距离(大于l即可)
n = 1000000; % 做n次投针试验,n越大求出来的pi越准确
m = 0; % 记录针与平行线相交次数
x = rand(1,n) * a / 2 ; % 在[0,a/2]内服从均匀分布随机产生n个数,表示针的中点和最近的一条平行线的距离
% 或可以替换成 x = unifrnd(0,a/2,1,n)
theta = rand(1,n) * pi; % 在[0,pi]内服从均匀分布随机产生n个数,表示针和平行线的夹角
% 或可以替换成 x = unifrnd(0,pi,1,n)
for i=1:n % 开始循环,判断每根针是否和直线相交
if x(i) <= l / 2 * sin(theta(i)) % 如果相交
m = m + 1; % 相交次数m加一
% %画图需要将n调小在1000左右,否则很慢
% plot(theta(i), x(i), 'r.') %用红色小点标记
% hold on % 保持原来的图窗
end
end
% fplot(@(theta)l/2*sin(theta),[0 pi],"Color",'r',"LineWidth",1.5) %画出轮廓线
% xlabel('θ') %标记坐标轴
% ylabel('x')
% xticks(0:pi:pi) %刻度范围和间距
% xticklabels({'0','\pi'}) %坐标刻度标记
% yticks(0:a/2:a/2);
% yticklabels({'0','0.5a'})
% set(gca,'FontSize',15,'Fontname', 'Times New Roman'); %调整坐标区字体
% axis([0,pi, 0,a/2]); % x轴在0-pi,y轴在0-a/2
% box on; % 画坐标轴框
p = m / n; % 针和平行线相交的频率
pi_needle = (2 * l) / (a * p); % 根据投针概率公式计算得到的pi
disp(['投针-蒙特卡罗方法得到pi为:', num2str(pi_needle)])
n为5000时的绘图结果如下,与理论面积法的图形比较符合。
当n取1000000时,注释画图代码,某次的运行结果:
投针-蒙特卡罗方法得到pi为:3.1444
为了使得求出的 π \pi π更加精确,对该模型进一步优化,采取求100个 π \pi π平均值的策略,优化模型代码如下:
%进一步优化:重复100次后求平均的pi
result = zeros(100,1); % 提前定义保存100次运算结果的矩阵
l = 1; a = 1.5;
n = 1000000;
for num = 1:100 %100次循环
m = 0; %每次的相交次数
x = rand(1,n) * a / 2;
theta = rand(1, n) * pi;
for i=1:n %100次中每次的试验次数
if x(i) <= l / 2 * sin(theta (i))
m = m + 1;
end
end
p = m / n; %频率
pi_needle = (2 * l) / (a * p);
result(num) = pi_needle; % 将每次的结果保存到结果矩阵result中
end
pi_mean = mean(result); % 计算result矩阵中保存的100次结果的均值
disp(['优化投针-蒙特卡罗方法得到pi为:', num2str(pi_mean)])
某次运行结果:
优化投针-蒙特卡罗方法得到pi为:3.1415
这种优化方式的合理性在于,单次仿真依赖于随机数生成器,其结果证明不了什么。多次模拟求平均值可以使得结果更接近真实值。
2.2 三门问题
简述:参加一档电视节目,节目组提供了ABC三扇门,其中一扇门后面有奖品(一辆汽车),其他两扇门后面是空的,主持人知道哪一扇门后面是奖品。假如你选择了B门,这时,主持人会打开另外两扇门其中一扇没有奖品的C门,问你要不要改变主意选择A门。由于电视节目中主持人的名字叫蒙提霍尔,因此“三门问题”也称“蒙提霍尔问题”。
该问题的答案是,换门的话,赢得奖品的概率是 2 / 3 2/3 2/3,不换门的话,赢得奖品的概率是 1 / 3. 1/3. 1/3.
三门问题的概率论解释
一种直观解释
可以采用穷举的思想考虑这个问题。有三扇门A,B,C,奖品在某扇门背后的概率是一样的,为
1
/
3
1/3
1/3,我们在这里假设奖品在A门后,先提前想好是否要改变主意,即换或不换。
选择某扇门的概率是相等的,均为
1
/
3.
1/3.
1/3.当选择了有奖品的门A时,主持人会在B,C中随机选择打开一扇门,则打开B,C的概率相同,均为
1
/
2.
1/2.
1/2.若改变主意,赢的概率为
2
/
3
2/3
2/3;不改变主意,赢的概率为
1
/
3.
1/3.
1/3.我们发现,在这里改变主意使得最终赢的概率
1
/
3
1/3
1/3变成 了刚开始时选错的概率
2
/
3.
2/3.
2/3.
条件概率解释
我们事先不知道奖品在哪扇门后,假设我们选择了门A,主持人打开了门B.设奖品在A,B,C门后分别为事件A,B,C,则在主持人没有开门前,有
P
(
A
)
=
P
(
B
)
=
P
(
C
)
=
1
/
3
P(A)=P(B)=P(C)=1/3
P(A)=P(B)=P(C)=1/3,这是先验概率。设主持人打开B门为事件BB,则奖品在A后的条件下,主持人打开B门的概率为
P
(
B
B
∣
A
)
=
1
/
2
P(BB|A)=1/2
P(BB∣A)=1/2;奖品在C后的条件下主持人打开B的概率为
P
(
B
B
∣
C
)
=
1
P(BB|C)=1
P(BB∣C)=1;奖品在B后的概率为
P
(
B
)
=
0
P(B)=0
P(B)=0.此时可以用全概率公式求出主持人打开B门的概率
P
(
B
B
)
=
P
(
B
B
∣
A
)
P
(
A
)
+
P
(
B
B
∣
B
)
P
(
B
)
+
P
(
B
B
∣
C
)
P
(
C
)
=
1
/
2
×
1
/
3
+
0
+
1
×
1
/
3
=
1
/
2
P(BB)=P(BB|A)P(A)+P(BB|B)P(B)+P(BB|C)P(C)=1/2\times 1/3+0+1\times 1/3=1/2
P(BB)=P(BB∣A)P(A)+P(BB∣B)P(B)+P(BB∣C)P(C)=1/2×1/3+0+1×1/3=1/2
再用贝叶斯公式求出主持人打开B门的条件下A的概率,也就是由于主持人的动作我们对先验概率做出调整,得到A后验概率
P
(
A
∣
B
B
)
=
P
(
A
)
P
(
B
B
∣
A
)
P
(
B
B
)
=
1
/
3
×
1
/
2
1
/
2
=
1
/
3
P(A|BB)=\frac{P(A)P(BB|A)}{P(BB)}=\frac{1/3\times1/2}{1/2}=1/3
P(A∣BB)=P(BB)P(A)P(BB∣A)=1/21/3×1/2=1/3
同理可以求出主持人打开B门的条件下,奖品在C后的概率为
P
(
C
∣
B
B
)
=
P
(
C
)
P
(
B
B
∣
C
)
P
(
B
B
)
=
1
/
3
×
1
1
/
2
=
2
/
3
P(C|BB)=\frac{P(C)P(BB|C)}{P(BB)}=\frac{1/3\times1}{1/2}=2/3
P(C∣BB)=P(BB)P(C)P(BB∣C)=1/21/3×1=2/3
那么我们改变主意,换选C门赢得奖品的概率要大一些。
蒙特卡洛模拟三门问题
编程需要用到的MATLAB语法
输出分布在[a,b]之间均匀分布的整数组成的m行n列矩阵:
randi([a,b],m,n) %没有,m,n默认随机生成一个[a,b]间的整数
连接两个字符型数组
//两种方法
['string1','string2']
strcat('string1','string2')
三门问题——蒙特卡洛模拟代码
代码1.不考虑失败的情况,仅对(不)改变主意能赢的次数做统计
%% 不考虑失败情况
n = 100000; % 试验的重复次数
a = 0; % 不改变主意时能赢得奖品的次数
b = 0; % 改变主意时能赢得奖品的次数
k = 3; % 门的个数(假设有3扇门)
for i= 1 : n % 开始模拟n次
x = randi([1,k]); % 随机生成一个1-k间的整数x表示奖品在第x扇门后
y = randi([1,k]); % 随机生成一个1-k间的整数y表示自己选的门
% 下面分为两种情况讨论:x=y和x~=y
if x == y % 如果x和y相同,只有不改变主意时才能赢
a = a + 1;
else % x ~= y ,如果x和y不同,只有改变主意时才能赢
b = b + 1;
end
end
disp(['蒙特卡罗方法得到不改变主意时的获奖概率:', num2str(a/n)]);
disp(['蒙特卡罗方法得到改变主意时的获奖概率:', num2str(b/n)]);
某次运行结果:
蒙特卡罗方法得到的不改变主意时的获奖概率为:0.33252
蒙特卡罗方法得到的改变主意时的获奖概率为:0.66748
代码2.考虑失败的情况,假设在每次试验中是否改变主意是随机的,得到失败的概率、改变主意赢的概率、不改变主意赢的概率
%% 考虑失败情况
n = 100000; % n代表蒙特卡罗模拟重复次数
a = 0; % 不改变主意时能赢得奖品的次数
b = 0; % 改变主意时能赢得奖品的次数
c = 0; % c表示没有获奖的次数
k = 3; % 门的个数(假设有3扇门)
for i = 1 : n % 开始模拟n次
x = randi([1,k]); % 随机生成1-k之间的整数x表示奖品在第x扇门后
y = randi([1,k]); % 随机生成1-k之间的整数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.16708
蒙特卡罗方法得到的改变主意时的获奖概率为:0.33312
蒙特卡罗方法得到的没有获奖的概率为:0.4998
2.3 模拟排队问题
假设某银行工作时间只有一个服务窗口,工作人员只能逐个接待顾客。当来的顾客较多时,一部分顾客就需要排队等待。
假设:
(1)顾客到来的间隔时间服从参数为0.1的指数分布
(2)每个顾客的服务时间服从均值为10,方差为4的正态分布(单位为分钟,若服务时间小于1分钟,则按1分钟计算)
(3)排队按先到先服务的规则,且不限制队伍的长度,每天工作时长为8小时。
试回答下面的问题:
(1)模拟一个工作日,在这个工作日共接待了多少客户,客户平均等待的时间为多少?
(2)模拟100个工作日,计算出平均每日接待客户的个数以及每日客户的平均等待时长。
排队问题理论模型
指数分布的概率密度函数为
f
(
x
)
=
λ
e
−
λ
x
f(x)=\lambda e^{-\lambda x}
f(x)=λe−λx
其中
λ
\lambda
λ是指数分布的参数,该指数分布的期望值为
1
λ
\frac{1}{\lambda}
λ1.
首先找出各个变量的迭代关系。设每个客户到达的时刻点为
a
a
a(arrive),开始服务的时刻点为
b
b
b(begin),结束服务离开的时刻点为
l
l
l(leave).
观察可以得出以下规律:
(1)客户开始服务的时刻点取决于上一个客户结束服务的时刻点以及该客户到达的时刻点,取二者中的较大值:
b
i
=
m
a
x
(
l
i
−
1
,
a
i
)
b_i=max(l_{i-1},a_i)
bi=max(li−1,ai)
(2)将客户到达的事件间隔设为
x
i
x_i
xi,该随机变量服从指数分布,则客户到达时刻可以表示为:
a
i
=
a
i
−
1
+
x
i
a_i=a_{i-1}+x_i
ai=ai−1+xi
(3)设每个客户服务的服务时间为
y
i
y_i
yi,该 随机变量服从均值为10,方差为4的正态分布,则客户服务结束离开的时刻为:
l
i
=
b
i
+
y
i
l_i=b_i+y_i
li=bi+yi
蒙特卡洛模拟排队问题
编程需要用到的MATLAB语法
生成一个服从均值为mu,标准差为sigma的正态分布随机数:
normrnd(MU,SIGMA)
生成一个均值为mu的指数分布随机数:
exprnd(mu)% 其参数为lambda,则mu=1/lambda
求一个向量x中各元素的均值:
mean(x)
排队问题——蒙特卡洛模拟代码
问题一、模拟一个工作日,假设银行每天的工作时间为 8 h = 480 m i n 8h=480min 8h=480min
i = 1; % 表示第i个客户
x(i) = exprnd(10); % 从银行开门到第一个顾客到达的时间间隔
a(i) = x(1); % 第一个顾客到达的时间
b(i) = a(1); wait(i) = 0; % 第一个顾客不用等待
l(i) = normrnd(10,2); % 第一个顾客结束服务离开的时间
while b(i) <= 480 % 第i个顾客开始服务的时间大于营业时间时终止
y(i) = normrnd(10,2); % 第i个顾客的服务时间服从均值为10,方差为4(标准差为2)的正态分布
if y(i) < 1 % 根据题目假设,服务时间小于1min记作1min
y(i) = 1;
end
l(i) = b(i) + y(i); % 第i个顾客离开的时间
wait(i) = l(i) - a(i); % 第i个顾客等待的时间
x(i + 1) = exprnd(10); % 第i+1个顾客到达时间与上一个顾客到达时间的间隔
a(i + 1) = a(i) + x(i+1); % 第i+1位顾客到达的时间,时间间隔服从均值为10的指数分布
b(i + 1) = max([a(i + 1),l(i)]); % 第i+1位顾客开始的时间是其到达时间和上一位顾客开始时间的较大值
i = i + 1; % 更新到下一位顾客
end
n = i - 1; % 服务的总人数
wait(end) = []; % 删掉最后一个超出营业时间的顾客的等待时间
w_mean = mean(wait); % 顾客等待时间的平均值
disp(['一个工作日服务的顾客总数:',num2str(n)])
disp(['顾客等待的平均时间:',num2str(w_mean)])
某次的运行结果:
一个工作日服务的顾客总数:45
顾客等待的平均时间:43.8722
问题二、模拟100个工作日
% 100个工作日取平均
for day = 1:100
i = 1; % 表示第i个客户
x(i) = exprnd(10); % 从银行开门到第一个顾客到达的时间间隔
a(i) = x(1); % 第一个顾客到达的时间
b(i) = a(1); wait(i) = 0; % 第一个顾客不用等待
l(i) = normrnd(10,2); % 第一个顾客结束服务离开的时间
while b(i) <= 480 % 第i个顾客开始服务的时间大于营业时间时终止
y(i) = normrnd(10,2); % 第i个顾客的服务时间服从均值为10,方差为4(标准差为2)的正态分布
if y(i) < 1 % 根据题目假设,服务时间小于1min记作1min
y(i) = 1;
end
l(i) = b(i) + y(i); % 第i个顾客离开的时间
wait(i) = l(i) - a(i); % 第i个顾客等待的时间
x(i + 1) = exprnd(10); % 第i+1个顾客到达时间与上一个顾客到达时间的间隔
a(i + 1) = a(i) + x(i+1); % 第i+1位顾客到达的时间,时间间隔服从均值为10的指数分布
b(i + 1) = max([a(i + 1),l(i)]); % 第i+1位顾客开始的时间是其到达时间和上一位顾客开始时间的较大值
i = i + 1; % 更新到下一位顾客
end
n = i - 1; % 服务的总人数
wait(end) = []; % 删掉最后一个超出营业时间的顾客的等待时间
w_mean = mean(wait); % 顾客等待时间的平均值
N(day) = n; % 储存每天的人数
W(day) = w_mean; % 储存每个工作日顾客等待的平均值
end
disp(['100个工作日每天服务顾客数的平均值:',num2str(mean(N))])
disp(['100个工作日每天顾客等待平均时间的平均值:',num2str(mean(W))])
某次的运行结果:
100个工作日每天服务顾客数的平均值:42.93
100个工作日每天顾客等待平均时间的平均值:40.3217
2.4 有约束的非线性规划问题
数学规划问题
{
min
f
(
x
)
目标函数
s
.
t
.
g
(
x
)
≤
0
不等式约束
h
(
x
)
=
0
等式约束
\begin{cases} \min{f(x)} \text{ 目标函数}\\ s.t. \quad g(x) \leq 0\text{ 不等式约束}\\ \quad \quad h(x)=0\text{ 等式约束} \end{cases}
⎩⎪⎨⎪⎧minf(x) 目标函数s.t.g(x)≤0 不等式约束h(x)=0 等式约束
线性规划(LP):
f
(
x
)
,
g
i
(
x
)
,
h
i
(
x
)
f(x),g_i(x),h_i(x)
f(x),gi(x),hi(x)为线性函数
非线性规划(NLP):
f
(
x
)
,
g
i
(
x
)
,
h
i
(
x
)
f(x),g_i(x),h_i(x)
f(x),gi(x),hi(x)至少一个为非线性。若不存在
g
i
(
x
)
,
h
i
(
x
)
g_i(x),h_i(x)
gi(x),hi(x)即
X
=
R
n
X=R^n
X=Rn为无约束非线性规划(无约束最优化问题)
蒙特卡洛模拟有约束的非线性规划问题
求目标函数
f
(
x
)
=
x
1
x
2
x
3
f(x)=x_1x_2x_3
f(x)=x1x2x3的最大值,约束条件如下:
s
.
t
.
{
−
x
1
+
2
x
2
+
2
x
3
≥
0
x
1
+
2
x
2
+
2
x
3
≤
72
10
≤
x
2
≤
20
x
1
−
x
2
=
10
s.t. \begin{cases} -x_1+2x_2+2x_3\geq 0\\ x_1+2x_2+2x_3\leq 72\\ 10\leq x_2\leq 20\\ x_1-x_2=10 \end{cases}
s.t.⎩⎪⎪⎪⎨⎪⎪⎪⎧−x1+2x2+2x3≥0x1+2x2+2x3≤7210≤x2≤20x1−x2=10
首先分别求出每一个自变量的所有可能取值范围,
10
≤
x
2
≤
20
10\leq x_2\leq 20
10≤x2≤20,
x
1
−
x
2
=
10
x_1-x_2=10
x1−x2=10,则
20
≤
x
1
≤
30
20\leq x_1\leq 30
20≤x1≤30,
x
3
≥
1
2
(
x
1
−
2
x
2
)
=
1
2
(
10
−
x
2
)
≥
−
5
x_3\geq \frac{1}{2}(x_1-2x_2)=\frac{1}{2}(10-x_2)\geq-5
x3≥21(x1−2x2)=21(10−x2)≥−5,
x
3
≤
1
2
(
72
−
x
1
−
2
x
2
)
=
1
2
(
62
−
3
x
2
)
≤
16
x_3\leq \frac{1}{2}(72-x_1-2x_2)=\frac{1}{2}(62-3x_2)\leq 16
x3≤21(72−x1−2x2)=21(62−3x2)≤16
即
{
10
≤
x
2
≤
20
20
≤
x
1
≤
30
−
5
≤
x
3
≤
16
\begin{cases} 10\leq x_2\leq 20\\ 20\leq x_1\leq 30\\ -5\leq x_3\leq 16 \end{cases}
⎩⎪⎨⎪⎧10≤x2≤2020≤x1≤30−5≤x3≤16
先进行初次尝试
n = 10000000; % 生成的随机数的组数
x2 = unifrnd(10,20,n,1); % [10,20]之间均匀分布的n个随机数
x1 = x2 + 10; % 约束
x3 = unifrnd(-5,16,n,1);% 约束,[-10,16]之间均匀分布的随机数
fmax = -Inf; %(求最大值)设初始目标函数最大值为负无穷
for i = 1:n
if ((-x1(i) + 2*x2(i) + 2*x3(i))>=0) && ((x1(i) + 2*x2(i) +2*x3(i))<=72) % 判断生成的解是否满足条件
x = [x1(i),x2(i),x3(i)]; % 一个解向量
f_x = x1(i)*x2(i)*x3(i); % 目标函数值
if f_x >= fmax % 判断新目标函数值是否最大
fmax = f_x; % 更新最大目标函数值
X = x; % 更新最优解
end
end
end
format long % 显示更精确的值
disp("用蒙特卡洛方法求出的非线性规划问题的最大目标函数值为:")
fmax
disp("用蒙特卡洛方法求出的非线性规划问题的最优解为:")
X
某次的运行结果如下:
用蒙特卡洛方法求出的非线性规划问题的最大目标函数值为:
fmax = 3.445484294987000e+03
用蒙特卡洛方法求出的非线性规划问题的最优解为:
X = 22.632824076330230 12.632824076330232 12.050666809610465
接下来运用得到的结果将变量的取值范围进一步缩小,使得结果更加精确:
% 根据之前的结果缩小解的范围
n = 10000000; % 生成的随机数的组数
x2 = unifrnd(12,13,n,1); % [10,20]之间均匀分布的n个随机数
x1 = x2 + 10; % 约束
x3 = unifrnd(11,13,n,1);% 约束,[-10,16]之间均匀分布的随机数
fmax = -Inf; %(求最大值)设初始目标函数最大值为负无穷
for i = 1:n
if ((-x1(i) + 2*x2(i) + 2*x3(i))>=0) && ((x1(i) + 2*x2(i) +2*x3(i))<=72) % 判断生成的解是否满足条件
x = [x1(i),x2(i),x3(i)]; % 一个解向量
f_x = x1(i)*x2(i)*x3(i); % 目标函数值
if f_x >= fmax % 判断新目标函数值是否最大
fmax = f_x; % 更新最大目标函数值
X = x; % 更新最优解
end
end
end
format long % 显示更精确的值
disp("用蒙特卡洛方法求出的非线性规划问题的最大目标函数值为:")
fmax
disp("用蒙特卡洛方法求出的非线性规划问题的最优解为:")
X
某次的运行结果如下:
用蒙特卡洛方法求出的非线性规划问题的最大目标函数值为:
fmax = 3.445603424713144e+03
用蒙特卡洛方法求出的非线性规划问题的最优解为:
X = 22.583486684921731 12.583486684921731 12.124764237429215
2.5 0-1规划问题(书店买书)
以书店买书问题为例:某同学要从六家线上商城选购五本书籍B1,B2,B3,B4,B5,每本书籍在不同商家的售价以及每个商家的单次运费如下表所示,请给该同学制定最省钱的选购方案。(注:在同一个店买多本书也只会收取一次运费)
编程需要用到的MATLAB语法
删去一个矩阵或者向量的重复值,并将结果按照从小到大的顺序排列,结果以一个列向量的形式输出
unique(A) % 将矩阵或向量A中的元素去重并从小到大排列后以列向量的形式输出
蒙特卡洛模拟书店买书问题
这个问题可能的结果有
6
5
=
7776
6^5=7776
65=7776种,因此将蒙特卡洛模拟的次数设为
n
=
100000
n=100000
n=100000很可能找到最佳方案,如果n较小,则结果会有一定的波动。
cost_min = Inf % 花费初始化为无穷大
result_min = randi([1, 6],1,5); % 初始化买书方案,五本书都在哪一家书店购买,有5本书6家书店,用于储存最优结果
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 i = 1:n
result = randi([1,6],1,5); % 随机生成的方案
index = unique(result); % 在哪些商店购买了商品,用于计算运费
sum_f = sum(freight(index)); % 计算买书花费的运费
sum_b = 0; % 初始化买书的花费为0
for k = 1:5
sum_b = sum_b + M(result(k),k); % 买书花的钱
end
cost = sum_f + sum_b; %总价格
if cost < cost_min
cost_min = cost; % 更新最低价格
result_min = result; % 更新买书方案
end
end
cost_min
result_min