【数学建模】自动化车床管理(蒙特卡洛模拟法)

自动化车床管理
一道工序用自动化车床连续加工某种零件，由于刀具损坏等原因该工序会出现故障，其中刀
具损坏故障占 95%, 其它故障仅占 5%。工序出现故障是完全随机的, 假定在生产任一零件时出现故障的机会均相同。工作人员通过检查零件来确定工序是否出现故障。
现积累有 100 次刀具故障记录，故障出现时该刀具完成的零件数如附表。现计划在刀具加工
一定件数后定期更换新刀具。
已知生产工序的费用参数如下：
故障时产出的零件损失费用 $f=200$ 元/件；
进行检查的费用 $t=10$ 元/次；
发现故障进行调节使恢复正常的平均费用 $d=3000$ 元/次(包括刀具费)；
未发现故障时更换一把新刀具的费用 $k=1000$ 元/次。

1）假定工序故障时产出的零件均为不合格品，正常时产出的零件均为合格品, 试对该工序
设计效益最好的检查间隔（生产多少零件检查一次）和刀具更换策略。
2）如果该工序正常时产出的零件不全是合格品，有 2%为不合格品；而工序故障时产出的零
件有 40%为合格品，60%为不合格品。工序正常而误认有故障仃机产生的损失费用为 1500 元/次。
对该工序设计效益最好的检查间隔和刀具更换策略。
3）在 2)的情况, 可否改进检查方式获得更高的效益。

附:100 次刀具故障记录(完成的零件数)
459 362 624 542 509 584 433 748 815 505
612 452 434 982 640 742 565 706 593 680
926 653 164 487 734 608 428 1153 593 844
527 552 513 781 474 388 824 538 862 659
775 859 755 649 697 515 628 954 771 609
402 960 885 610 292 837 473 677 358 638
699 634 555 570 84 416 606 1062 484 120
447 654 564 339 280 246 687 539 790 581
621 724 531 512 577 496 468 499 544 645
764 558 378 765 666 763 217 715 310 851

question1

分析问题可知:
如果进行检查要花费$t=10$元/次
检查有两种结果，花费$\{\begin{array}{l}d=3000,发现故障\\ 0,未发现故障\end{array}$
注：说的是 未发现故障时 更换一把新刀具费用是$k=1000$元，前置条件是未发现故障，判断要执行的动作是更换一把新刀具，如果要换才花$k$元，这里假定了假定工序故障时产出的零件均为不合格品，正常时产出的零件均为合格品，因此未发现故障，刀一定没坏，所以不换，花费为$0$元。

我们先对100 次刀具故障记录(完成的零件数) 的数据进行处理：
先用C++跑出 出现的故障记录(完成的零件数) 对应的次数是多少

map<int,int>cnt;
void solve(){
   int x;
   for(int i=0;i<100;i++){
       cin >> x;
       cnt[x]++;
   }
   for(pair<int,int> k:cnt){
       cout << k.first << ' ' << k.second << '\n';
   }
}

改个打印格式再放进MATLAB里面画出柱状图

x =[84 120 164 217 246 280 292 310 339 358 362 378 388 402 416 428 433 434 447 452 459 468 473 474 484 487 496 499 505 509 512 513 515 527 531 538 539 542 544 552 555 558 564 565 570 577 581 584 593 606 608 609 610 612 621 624 628 634 638 640 645 649 653 654 659 666 677 680 687 697 699 706 715 724 734 742 748 755 763 764 765 771 775 781 790 815 824 837 844 851 859 862 885 926 954 960 982 1062 1153];
y =[1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1];
bar(x,y)

可以很明显的看出线条集中在600偏左一点，出现次数为2的周边，查看之前的数据是593
算一下期望

double ex = 0;
    for(pair<int,int> k:cnt){
       ex += (double)k.first * (double)k.second/100;
    }
    cout << ex;

算一下期望得到$E(x)=600$
并且分析上图推测该数据满足正态分布。

C++算标准差/方差：

double ex = 0;
    double ex2 = 0;
    for(pair<int,int> k:cnt){
       ex += (double)k.first * (double)k.second/100;
        ex2 += (double)k.first*(double)k.first * (double)k.second/100;
    }
    cout << ex2-ex*ex << '\n';
    double dx = 0;
    for(pair<int,int> k:cnt){
        dx += ((double)k.first-ex)*((double)k.first-ex) * (double)k.second/100;
    }
    cout << sqrt(dx);

得到方差为38276.4，标准差为195.644
用MATLAB画出正态分布图，如下图

上面这一步应该是拟合

假设检验

略

建立模型

综上所述，最好的检查间隔是600个零件/次，每次检查都换刀具

开始建立模型：
已知:
故障时产出的零件损失费用 $f=200$ 元/件；
进行检查的费用 $t=10$ 元/次；
发现故障进行调节使恢复正常的平均费用 $d=3000$ 元/次(包括刀具费)；
未发现故障时更换一把新刀具的费用 $k=1000$ 元/次。

那么假设完成的零件数在$x$时会故障,如果完成的零件数在$y$时候检查，总花费$g(y)$
$g(y)=\{\begin{array}{l}(y-x)\ast f+d+t,x\le y\\ t,x>y\end{array}$

到目前为止难点就在于未发现故障时要不要更换一把新刀具并且其中刀具损坏故障占 95%, 其它故障仅占 5%。工序出现故障是完全随机的这个条件还没考虑到。
这里先直接假设刀具更换策略和检查策略分开，并且上述难点均分为刀具更换策略。

假设未检查到故障就不换刀,那么检查策略模型:
$\{\begin{array}{l}g(y)=\{\begin{array}{l}(y-x)\ast f+d+t,x\le y\\ t,x>y\end{array}\\ h(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}\\ \min E(g(y))=\int g(y)h(x)dx\\ \mu =600,{\sigma }^2=38276.4\end{array}$
求出检查间隔是$y$个零件/次

刀具更换策略模型:
假设完成的零件数在$x$时会故障,如果完成的零件数在$y$时候换刀
因为刀具损坏故障占 95%所以
$g^{\prime }(y)=\{\begin{array}{l}(y-x)\ast f+d,x\le y\\ 95\%k,x>y\end{array}$

$\{\begin{array}{l}h(x)=\frac{1}{\sqrt{2\pi }\sigma }{e}^{-\frac{(x-\mu {)}^2}{2{\sigma }^2}}\\ \min E(g(y))=\int g^{\prime }(y)h(x)dx\\ \mu =600,{\sigma }^2=38276.4\end{array}$
求出换刀间隔是$y$个零件/次

重新建立模型

1.确立检查间隔$n$(检查间隔相同)
2.确立刀具更换周期$m=sn$,即使不发生故障也要更换
总费用期望最低为目标
已知刀具寿命服从正态分布$N(600,195.644^2)$

1.预防更换费用(刀具寿命超过更换周期$m$)
刀具更换前检查$s-1$次
$c_1=(s-1)t+k$ 概率为$1-F(m)$
第$i$次检查发现故障更换费$(i=1,..,s-1)$
$c_{2i}=it+d+f(n+1)/2$ 概率为$F(in)-F((i-1)n)$
(第$i-1$到第$i$次检查种平均生产$\frac{n+1}{2}$不合格零件)

随机化求解

已知:
故障时产出的零件损失费用 $f=200$ 元/件；
进行检查的费用 $t=10$ 元/次；
发现故障进行调节使恢复正常的平均费用 $d=3000$ 元/次(包括刀具费)；
未发现故障时更换一把新刀具的费用 $k=1000$ 元/次。

定一个$n$和$m$
每次随机一个$x$服从正态分布$N(600,195.644^2)$，相当于换新刀，随机生成刀的寿命
花费$W=\{\begin{array}{l}(\lceil x/n\rceil )t+d+(\lceil x/n\rceil n-x)f,x\le m\\ (m/n)t+k,x>m\end{array}$
生产个数$T=\{\begin{array}{l}x,x\le m\\ m,x>m\end{array}$
计算生成每个零件平均价值$\frac{W}{T}$
随机$k$次$x$得到平均值$\frac{\sum _{i=1}^k\frac{W}{T}}{k}$
取平均值最小的$n$,$m$为解即可

蒙特卡洛模拟
matlab

t = 10; %进行检查的费用
k = 1000; %未发现故障时更换一把新刀具的费用
d = 3000; %发现故障进行调节使恢复正常的平均费用
f = 200; %故障时产出的零件损失费用
ans_n = 0;
ans_m = 0;
ans_wt = 1e9;
x = abs(normrnd(600,195.644,10000));
for n = 1:50
    for s = 1:100
        m=s*n;
        wt = 0;
        for i = 1:length(x);
            w = 0;
            t = 0;
            if x(i) < m
                w = w+ceil(x(i)/n)*t+d+(ceil(x(i)/n)*n-x(i))*f;
                t = t+x(i);
            else
                w = w+s*t+k;
                t = t+m;
            end
            wt = wt + w/t;
            if ans_wt*1000 < wt
                break;
            end
        end
        if ans_wt > wt/1000
            ans_n = n;
            ans_m = m;
            ans_wt = wt/1000;
        end
    end
end
disp(ans_n);
disp(ans_m);