引例

布丰投针实验

法国数学家布丰（1707-1788）最早设计了投针试验。

这一方法的步骤是：

1） 取一张白纸，在上面画上许多条间距为a的平行线。

2） 取一根长度为l（l≤a） 的针，随机地向画有平行直线的纸上掷n次，观察针与直线相交的次数，记为m。

3）计算针与直线相交的概率．

18世纪，法国数学家布丰提出的“投针问题”，记载于布丰1777年出版的著作中：“在平面上画有一组间距为a的平行线，将一根长度为l（l≤a）的针任意掷在这个平面上，求此针与平行线中任一条相交的概率。”

布丰本人证明了，这个概率是：
\[
p=\frac{2l}{\pi a}
\]

(其中π为圆周率）

由于它与π有关，于是人们想到利用投针试验来估计圆周率的值。

布丰惊奇地发现：有利的扔出与不利的扔出两者次数的比，是一个包含π的表示式．如果针的长度等于a/2，那么有利扔出的概率为1/π．扔的次数越多，由此能求出越为精确的π的值。

下面是利用这个公式，用概率的方法得到圆周率的近似值的一些资料。

证明

代码

%%  蒙特卡罗用于布丰投针实验

%% （1）预备知识
%  rand(m,n)函数产生由在[0,1]之间均匀分布的随机数组成的m行n列的矩阵（或称为数组）。
rand(5,4)
%     0.8300    0.1048    0.2396    0.4398
%     0.5663    0.1196    0.8559    0.5817
%     0.9281    0.2574    0.3013    0.9355
%     0.3910    0.3173    0.2108    0.1676
%     0.3645    0.4372    0.8819    0.9232
rand(3) % 若只给一个输入，则会生成一个方阵
%     0.1709    0.4951    0.0541
%     0.9374    0.8500    0.6155
%     0.2400    0.3156    0.5741
% a + rand(m,n)*(b-a) 可以输出在[a,b]之间均匀分布的随机数组成的m行n列的矩阵。
-2 + rand(3,2) * (2 - (-2))  % 输出在[-2,2]之间均匀分布的随机数组成的3行2列的矩阵。
%    -1.2743    0.6013
%    -1.3084    0.0766
%     1.5075    0.7563
% a + rand(m,n)*(b-a)等价于unifrnd(a,b,m,n)
unifrnd(-2,2,3,2)

%% （2）代码部分
l =  0.520;     % 针的长度（任意给的）
a = 1.314;    % 平行线的宽度(大于针的长度l即可)
n = 1000000;    % 做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就要加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)])

%% （3） 由于一次模拟的结果具有偶然性，因此我们可以重复100次后再来求一个平均的pi
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)])

蒙特卡罗方法概述

定义

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

提出

? 蒙特卡罗?法于20世纪40年代美国在第?次世界?战中研制原?弹的“曼哈顿计划”计划的成员S.M.乌拉姆和J. 冯·诺伊曼?先提出。数学家冯·诺伊曼?驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种?法，为它蒙上了? 层神秘?彩。在这之前，蒙特卡罗?法就已经存在。1777年，法国Buffon提出?投针实验的?法求圆周率，这被认 为是蒙特卡罗?法的起源。

原理

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

讨论

1. 蒙特卡罗是?种算法吗 ?

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

2. 蒙特卡罗与计算机仿真的关系

   计算机仿真（模拟）早期称为蒙特卡罗?法，是??利?随机数实验求解随机问题的?法，其主要应?在复杂问题的数 值模拟上。但随着计算机的性能的提?以及各种新兴产业的发展，?前计算机仿真涉及的内容要?得多，例如过程控制、动 画仿真、图像静态模拟等都属于计算机仿真的应?（例如?计算机模拟原?弹爆炸的过程、使?计算机?成特效??等）。 在数学建模中，我们不?刻意的去区分两者的区别，?家只需要知道如何使?计算机对问题进?模拟即可。

3. 蒙特卡罗与枚举法

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

已经学过的例?

第?期 视频 (AOL) 第七讲 : 多元回归分析

探究内?性对回归结果的影响

第?期 视频 (AOL) 番外篇 : 基于熵权法对印的模型的修正

信息熵和标准差的关系

蒙特卡罗?法 的应?实例

三?问题

? 你参加?档电视节?，节?组提供了ABC三扇?， 主持?告诉你，其中?扇?后边有辆汽?，其它两扇? 后是空的。 假如你选择了B?，这时，主持?打开了C?，让你 看到C?后什么都没有，然后问你要不要改选A?？

有约束的?线性规划问题

代码

%%  蒙特卡罗求解有约束的非线性规划问题
% max f(x) = x1*x2*x3
% s.t.
% (1) -x1+2*x2+2*x3>=0
% (2) x1+2*x2+2*x3<=72
% (3) x2<=20 & x2>=10
% (4) x1-x2 == 10

%% （1）预备知识
%  (1) format long g  可以将Matlab的计算结果显示为一般的长数字格式（默认会保留四位小数，或使用科学计数法）
5/7
5895*514100
format long g
5/7
5895*514100
%  (2)unifrnd(a,b,m,n)可以输出在[a,b]之间均匀分布的随机数组成的m行n列的矩阵。(等价于 a + rand(m,n)*(b-a))
unifrnd(0,5,4,3)
%           4.07361843196589          3.16179623112705          4.78753417717149
%            4.5289596853781         0.487702024997048          4.82444267599638
%           0.63493408146753          1.39249109433524         0.788065408387741
%            4.5668792806951          2.73440759602492          4.85296390880308

%% （2）代码部分
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: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 %计算tic和toc中间部分的代码的运行时间

%% （3）缩小范围重新模拟得到更加精确的取值
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: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 %计算tic和toc中间部分的代码的运行时间

书店 买书 问题 ( 0 - 1 规划)(无脑)

%% 书店买书问题的蒙特卡罗的模拟
%% （1）预备知识
% (1)unique函数: 剔除一个矩阵或者向量的重复值，并将结果按照从小到大的顺序排列
% adj.  唯一的; 独一无二的   [ju'ni:k]
unique([1 2 5; 6 8 9;2 4 6])
unique([5 6 8 8 4 1 6 2 2 4 8 4 5 6])

% (2)randi([a,b],m,n)函数可在指定区间[a,b]内随机取出大小为m*n的整数矩阵
randi([-5,5],2,6)

%% （2）代码求解
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  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 = k([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   % 18+39+48+17+47+20
min_result

导弹追踪问题

%%  蒙特卡罗用于模拟导弹追击问题
clear;clc
%% （1）预备知识
% mod(m,n)表示求m/n的余数
mod(8,3)
mod(1000,50)

% 设置横纵坐标的范围并标上字符
x = 1:0.01:3;
y = x .^ 2;
plot(x,y)  % 画出x和y的图形
axis([0 3 0 10])  % 设置横坐标范围为[0, 3] 纵坐标范围为[0, 10]
pause(3)  % 暂停3秒后再继续接下来的命令
text(2,4,'清风')  % 在坐标为(2,4)的点上标上字符串：清风
close % 关闭图形窗口

%% (2) 代码求解
% 1. 不画追击的示意图
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 = 1
    x(1)=x(1)+3*v*dt*cos_alpha;   y(1)=y(1)+3*v*dt*sin_alpha; % 计算新的导弹的位置
    if d>50  % 导弹的有效射程为50个单位
        disp('导弹没有击中B船');
        break;  % 退出循环
    end
    if d<=50 & dd<0.001   % 导弹飞行的距离小于50个单位且导弹和B船的距离小于0.001（表示击中）
        disp(['导弹飞行',num2str(d),'单位后击中B船'])
        disp(['导弹飞行的时间为',num2str(t*60),'分钟'])
    end
end

% 2. 画追击的示意图
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:2
    plot(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:2
            plot(x(i),y(i),'.k','MarkerSize',1);
            hold on; % 不关闭图形，继续画图
        end
        pause(0.001);  % 暂停0.001s后再继续下面的操作
    end
    if d>50  % 导弹的有效射程为50个单位
        disp('导弹没有击中B船');
        break;  % 退出循环
    end
    if d<=50 & dd<0.001   % 导弹飞行的距离小于50个单位且导弹和B船的距离小于0.001（表示击中）
        disp(['导弹飞行',num2str(d),'个单位后击中B船'])
        disp(['导弹飞行的时间为',num2str(t*60),'分钟'])
    end
end

作业

作业1

%%  蒙特卡罗求解非线性规划问题
% min f(x) =2*(x1^2)+x2^2-x1*x2-8*x1-3*x2
% s.t.
% (1) 3*x1+x2>9
% (2) x1+2*x2<16
% (3) x1>0 & x2>0

%% （1）初次寻找最小值的代码
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: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;  % 并且将此时的x1 x2 保存到相应的变量中
        end
    end
end
disp(strcat('蒙特卡罗模拟得到的最小值为',num2str(fmin)))
disp('最小值处x1 x2的取值为：')
disp(X)
toc %计算tic和toc中间部分的代码的运行时间

%% （2）缩小范围重新模拟得到更加精确的取值
clc,clear;
tic %计算tic和toc中间部分的代码的运行时间
n=10000000; %生成的随机数组数
x1=unifrnd(2,3,n,1);  % 生成在[2,3]之间均匀分布的随机数组成的n行1列的向量构成x1
x2=unifrnd(2,3,n,1);  % 生成在[2,3]之间均匀分布的随机数组成的n行1列的向量构成x2
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;  % 并且将此时的x1 x2 保存到相应的变量中
        end
    end
end
disp(strcat('蒙特卡罗模拟得到的最小值为',num2str(fmin)))
disp('最小值处x1 x2的取值为：')
disp(X)
toc %计算tic和toc中间部分的代码的运行时间

作业2

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

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

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

% find函数的用法
% 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中最小元素的位置

%% （2）代码部分
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; % 更新所有电子管的寿命
    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

原文地址：https://www.cnblogs.com/pxlsdz/p/12385321.html