数学建模--蒙特卡洛算法之电子管更换&&刀片寿命问题

目录

1.电子管问题重述

2.电子管问题分析

3.电子管问题求解

4.刀片问题重述

5.刀片问题分析

6.刀片问题求解

---

1.电子管问题重述

某设备上安装有4只型号规格完全相同的电子管，已知电子管寿命服从100～200h之间的均匀分布．

只要有一个电子管损坏，设备就不能工作。

电子管损坏时有两种维修方案：

一、是每次更换损坏的那只；

二、是当其中1只损坏时4只同时更换．

已知更换时间为换1只时需1h，4只同时更换为2h．

更换时机器因停止运转每小时的损失为20元，又每只电子管价格10元，试用模拟方法确定哪一个方案经济合理？

这个里面是出现了两个方案，方案一就是每一次只用去更换掉损坏的那一个；

第二个方案就是当一个损坏的时候，全部更换掉；

这个里面无论是方案一还是方案二，都会产生一个停止运转（更换时间以内）的损失问题，我们要综合考虑这个更换成本和更换损失，计算出这个最经济的方案（就是两个的和最小，即损失加成本的和最小）；

2.电子管问题分析

我们在分析问题的时候就去建立这个对应的数学模型，对于这个出现在模型里面的变量代表的含义，我们一定要搞清楚，这个对于我们的后续代码的编写很重要；

我们的这个思路就是去使用成本加上这个换电子管的损失费用之和相比较，费用低的就是经济的选择，c表示的就是一个电子管的更换成本，f表示每次更换的时候的损失成本，肯呢个有些同学会好奇这个里面的n代表的实际意义是什么，实际上这个n表示的是一个次数，下面我们细讲一下；

我们标题就说了，这个里面的两个问题就是蒙特卡洛算法的问题，而我们的蒙特卡洛算法的实质就是多次模拟，实际上这个关系是这样的：仿真模拟问题里面有很多的案例，这个蒙特卡洛算法就是这个里面的一种仿真模拟的方法，两个的实质很相似，都是多次模拟；

而我们写的这个四个人相互追逐的问题也属于方仿真问题，只不过这个属于过程仿真，需要我们在坐标系里面去解决，他们两个的关系就是下面的图里面的展示结果；

理清了关系，我们介绍思路，这个电子管问题涉及到了时间，我们却要进行多次仿真模拟，这个时候，我们就可以设置一个较长的时间，在这个实际上里面会进行多次的仿真模拟，但是也是会有一些不在这个范围里面，比如一个小写的数据，就1000小时，我们在这个时间里面进行模拟多次，可能有些电子管在这个时间的最后没有完全完成更换，就是这个时间可能不会完全使用完毕，我们只需要计算这个时间里面的完整电子管的相关费用即可，具体到matlab代码里面就是使用length函数求解一个向量的长度，这个长度就是我们这个图片里面的n的数值；

3.电子管问题求解

clc
clear
c=10;  %单个电子管的价格
f=20;  %单位时间的更换损失
ct1=1; %方案1的更换时间
ct2=2; %方案2的更换时间


T=10000;%仿真观测时间

tic %开始计时

for i=1:1000  %仿真次数
    R=100+100*rand(4,100); %模拟产生4乘100只电子管的寿命
    t=cumsum(R,2);          %转换为时间轴的时刻

    %按照方案1计算成本
    tt1=sort(reshape(t,1,[])); %将所有的故障时刻转换到同一时间轴 并排序
    t1=tt1(find(tt1<=T));      %筛选出截止时间T以内的数据
    C1(i)=length(t1)*(c+f*ct1);  %方案1的总成本

    %%按照方案2计算成本
    tt2=cumsum(min(R));     %每次找到寿命最短的电子管更换
    t2=tt2(find(tt2<=T));  %筛选出截止时间T以内的数据
    C2(i)=length(t2)*(4*c+f*ct2); %方案2的成本
end

MC1=mean(C1)
MC2=mean(C2)
toc  %结束计时

这个代码里面有很清楚的注释，下面简单的说明几点：

***length(t1)这些就是我们的分析里面说到的，选出来这个有效时间里面的数据的个数，作为我们的n的具体的数值；

***这个里面实际上是绘制的这个模拟过程的时间轴，就是把每一个灯管做成一个时间轴，每一次进行模拟的时候都会进行累积的，我们最后的mean函数再去求解这个平均值；

***rand函数就是生成的随机数，只不过这个里面是有4灯管，每一个模拟次数是100次，总共的模拟时间就是10000小时；

***R是电子管的寿命的具体的数值，这个我们如何去进行时间轴的转换呢，就是使用的cumsum函数，这个函数的第一个参数R就是寿命的时间长短，2表示行方向上进行求和；

***就是每一个电子管在模拟时间里面多次进行更换，我们cumsum参数2就是行方向求和（默认的话就是列方向求和），这样的话每一个电子管的使用寿命时长就会累加，产生一个类似于时间轴的效果；

4.刀片问题重述

5.刀片问题分析

上面这个刀片问题，很多人无法理解（最开始的时候我自己想了很长时间都想不明白），为了帮助大家理解这个问题，下面通过一个简单的图形帮助大家去理解问题：

就是我们有很多的刀片，我们可能在一个时间去更换他，但是这个刀片的寿命是不一样的，有的长，有的短，我们如果过早地更换，可能这个刀片产生的利润就不是很多，但是时间过长的话，我们的利润多，可能会使得一些刀片损坏而产生损失的费用，比如，我们把这个更换周期设置为T1，这个时候只有刀片11会产生损失，其他的都是产生的利润，更换周期变长之后，这个利润会变多，例如T2，但是损坏的刀片的数量也会变多，这个损失费用也会变大，这个也有优化问题的成分，让我们去求解利润最大的时候的更换时间，相当于是一个最优解；

6.刀片问题求解

我们在这个过程里面有很多的中间变量，为了帮助大家理解，我把这个相关的中间结果也展示在了实时脚本里面；

下面是代码：

clc
clear
N=10000;  %%刀具数目为N
RC=10:10:1000; %%设定状态变化量离散点 刀具使用时间

for j=1:100     %%设定仿真次数

   for i=1:length(RC)  %%不同的更换周期

      R=normrnd(500,100,1,N);  %%产生正态分布随机数
      RR=R;                 %%用于记录R的原始值
      Rc=RC(i);             %%读取第i个刀具更换周期
      r=find(R>Rc);         %%查询寿命大于设定更换周期的刀具序号

      RR(r)=Rc;      %%替换为每把刀具的实际加工时间，RR代表了所有刀具的实际加工时间
      n1=length(r);  %%寿命大于设定时间的刀具
      n2=N-n1;       %%寿命小于设定时间的刀具
      get=sum(RR)*10;   %%计算毛收益
      lost=n2*3000;     %%计算危害损失
      get2(i,j)=get-lost;  %%计算净收益
    end
end

figure
plot(RC,mean(get2'))
xlabel('刀具更换周期/h')
ylabel('净收益')
g=mean(get2');
k=find(g==max(g))
RC(k)

***这个里面的normrnd函数，第一个参数是均值，第二个是方差，1，N表示的就是1行N列的矩阵

***find函数查找的就是寿命长，没有损坏的刀片，把这个更换的时间作为其实际工作时间去更新数据矩阵；

***mean函数求解平均值，默认的话就是对于列求解平均值，其他的注释都很清楚，这里不再赘述，有问题欢迎评论区或者私信留言~~