数学建模 - 自动化车床管理中的故障检查与刀具更换模型

自动化车床管理中的故障检查与刀具更换模型

摘要

本文针对自动化车床管理中的故障检查与排除问题，以效益最好（即使总损失费用最少）为目标，建立不同条件下的检查间隔模型和刀具更换间隔模型，并提出改进的检查方式。

为了估计刀具的寿命分布，本文对附件中100次刀具故障记录的数据进行K-S法正态检验，得出有99.4%的概率可以认为刀具寿命服从正态分布。再运用极大似然法估计其正态分布参数，得其服从分布。

对于最优刀具更换间隔的求解，以效益最好为目标，将零件平均损失期望作为目标函数，利用刀具寿命的分布，解得效益最好时的刀具更换间隔为426个零件。

对于检查间隔的求解，需要在以下两种条件下分别讨论：

在工序正常时只产生合格品、工序故障时只产生不合格品的条件下，认为发生故障后的第一次检查即能检查出故障，依此构造出以检查间隔为自变量的零件平均损失期望函数。若以等间隔的方式进行检查，解得效益最好时的检查间隔为15个零件，且零件平均损失费期望最小值为4.51元/个。

在工序正常时生产2%的不合格品、工序故障时生产60%不合格品的条件下，因故障与否不一定能被准确查出，所以根据工序运行状态与检查结果分四种情况讨论：工序故障且检查结果为有故障、工序故障但检查结果为无故障、工序无故障且检查结果为无故障、工序无故障但检查结果为有故障。在此基础上推导出关于检查间隔的零件平均损失期望函数，求解得出使得效益最好的检查间隔为20个零件，且零件平均损失费期望最小值为6.87元/个。

最后，为了获得更高的效益，将检查方案由一次检查一个零件调整为一次检查多个零件再判断是否故障。保持求得的最优检查间隔（20个零件）不变，求得最佳的检查方案为：工作人员每次检查3个零件，若大于等于2个零件不合格即可判断车床有故障。另外，提出非等间隔的检查方式作为模型改进方向。

一、 问题重述

随着我国工业生产的飞速发展，制造业的生产技术已经进入自动化生产时代。但是，自动车床一旦发生故障而又未能及时检查出来，将会产生大量的不合格品，不仅给企业带来严重的经济损失，而且造成资源的严重浪费。因此，对自动化车床进行高效的管理显得尤为必要。

现已知某自动化车床连续加工零件，由于各种因素出现故障导致工厂遭受损失，效益变低。根据故障类型、发生的频率，将自动车床发生的故障分为刀具损坏故障，其它故障。分别占故障总数的95%和5%。为了及时检查出机器是否发生故障，工作人员通过检查零件是否出现故障。根据积累的100次刀具故障记录和生产工序的费用参数，解决以下三个问题：

（1）假定工序故障时产出的零件均为不合格品，正常时产出的零件均为合格品，试制定检查间隔与刀具更换策略使得工厂遭受的损失最少；

（2）假设工序正常时产出的零件不合格品概率为2%，工序故障时产出的零件有40%为合格品。工序正常而误认为有故障造成车床停机会造成1500元/次的损失费。试改进问题一的检查间隔和刀具更换策略；

（3）在（2）的假设条件下，考虑改进的检查方式以获得更高的效益。

二、 问题分析

题目中工序出现故障是完全随机，且生产任一零件时出现故障的机会均相同。由此，在理论上我们可以将问题转化为概率模型。

首先，根据刀具故障占95%这一条件，可以认为故障的发生与刀具的寿命密切相关。为了估计刀具的寿命，通过对附件中100次刀具故障记录的数据进行分析与处理，将100次刀具故障记录的数据绘制出频率分布图，可以发现其图像较符合正态分布。为了验证该猜想，对其做进一步的正态分布检验，发现其正态分布的检验效果很好。据此，可以进一步用极大似然估计法确定其分布中的参数。

由于发生故障后调节使恢复正常时需花费3000元/次，大于更换刀具费用1000元/次。则认为其具体做法为：在刀具更换时先对工序进行检修，若有故障则对其进行修理，花费3000元，若无故障则更换刀具，花费1000元。

由于需要求解的参数有刀具更换间隔和检查间隔两个参数，如果同时找其最优解会十分复杂，因此我们采取分级优化的方式，即先求得一个参数的最优解，再以该参数为基础，求得另一个参数的最优解。由于刀具的更换主要考虑的因素为刀具的寿命，且更换刀具可以避免刀具故障的发生，因此，刀具更换间隔对检查间隔的依赖很小，且对效益的影响程度很高。所以，我们采用先求解刀具更换间隔、再以此为基础求解检查间隔的方法。

因为更换刀具可以避免刀具故障，故对故障的发生点分两种情况讨论：故障发生在刀具更换间隔内和故障发生在刀具更换间隔外。利用问题中不同的假设条件，对这两种情况作详细讨论，即可求解出效益最好的检查间隔和刀具更换间隔。

对于检查间隔的求解，需要在以下两种条件下分别进行讨论：

在工序正常时只产生合格品、工序故障时只产生不合格品的条件下，检查出不合格零件即可认为车床发生故障。根据刀具寿命的分布，可以得出工厂损失期望的表达式，利用MATLAB即可解出最优的检查间隔。

当考虑到工序正常也会生产不合格品，工序故障也会生产合格品时，检查时可能会产生误判。工序正常时认为工序故障会造成停机，工序故障时认为工序正

常会造成不合格品比例增大。故可根据工序正常与否、检查结果正确与否分四种情况进行讨论，即：工序故障且检查结果为有故障、工序故障但检查结果为无故障、工序无故障且检查结果为无故障、工序无故障但检查结果为有故障四种情况。再利用刀具寿命的分布，可以列出总损失期望的表达式，同样利用MATLAB即可解出最优的检查间隔。

问题三假设与问题二相同，但由于前面两个问题在检查时均通过检查一个零件来判断车床是否发生故障，因此这种检查方式很容易导致判断错误。故考虑改进检查方案，可以通过增加检查零件的个数，来减小判断错误出现的概率。建模时确定检查零件中不合格品的比例，当不合格品的比例超过该值后即可判断车床出现故障。基于问题二的检查间隔和刀具更换策略，可以得出损失期望的表达式，并求解得出最优的比例系数。

三、 问题假设

（1）假设自动化车床故障种类中其它故障任何措施都无法避免； （2）假设工人每次检查能准确检测合格品与不合格品； （3）假设当故障发生后更换刀具，按维修情况来处理；

（4）假设只有工序处于故障状态下生产的不合格品才需要考虑零件损失费； （5）假设两次检查期间故障发生时刻服从均匀分布； （6）假设生产加工时每个零件合格与否相互独立；

四、 符号与变量说明

符号 检查间隔 刀具更换周期 刀具寿命 刀具寿命的分布函数 刀具寿命的概率密度分布 自动车床连续加工零件个数 刀具运行一次总费用的期望 故障时产出的零件损失费用 发现故障调节使恢复正常的费用 说明 单位 个 个 个 无 无 个 元 元 元 未发现故障时更换一把新刀具的费用 刀具更换间隔内刀具连续工作的平均长度 工序正常而误认有故障仃机产生的损失费用 元/次 个 元

五、 模型建立与求解

建立车床管理模型，使得检查间隔和刀具更换周期达到最好的效益。重点在于确切的估计出刀具寿命的分布。根据刀具寿命的概率密度函数，估计出连续加工零件时每一个零件出现故障的概率。结合生产工序的费用参数，确定损失的期望。优化检查间隔和刀具更换间隔，使得损失的期望最小。 5.1建模准备

5.1.1确定刀具寿命分布 1、样本数据分析

为了估计刀具寿命的分布，根据题目附件提供的100次刀具故障记录的数据，运用spss软件对其进行分析和处理，结果如下：

表1:100次刀具故障记录数据统计表 N 100 均值 600.00 标准差 196.629 M2统计值 38276 M2统计值^1/2 195.64

图1 ：故障件完成零件个数直方图

根据直方图发现，刀具故障完成零件数的频数分布与正态分布函数曲线有很好的拟合度。 2、 样本数据正态检验：

（1）对故障间完成的零件数运用spss进行正态P-P分析，P-P图拟合结果为：

图2 ：故障间完成的零件的正态P-P图

P-P图拟合结果为一条近似的直线，拟合结果较好。故采用正态分布描述刀具寿命的分布。

(2)对刀具故障数据运用spss进行K-S检验，检验结果为：

表2： 单样本 Kolmogorov-Smirnov 检验 刀具寿命 样本容量N 100 均值 600.0000 正态参数 标准差 196.62917 绝对值 0.042 最极端差别 正 0.042

负 Kolmogorov-Smirnov Z 渐近显著性(双侧) -0.041 0.421 0.994

3、正态分布参数估计

假设刀具寿命服从正态分布，且假设题目所给样本来自于正态总体，用极大似然估计法估计正态参数以确定其分布：

设是样本观测值，由于正态总体的概率密度为

构造似然函数为：

对似然函数取对数得：

由似然方程组

解得

代入数据计算得:

故刀具寿命服从参数为的正态分布，概率密度为：

5.2刀具更换间隔模型

刀具更换策略通过预测故障发生时间，在连续工作生产一定数目的零件后，定期更换刀具以减少损失。根据刀具寿命的分布，故障可能发生在刀具更换之前，也可能发生在刀具更换之后。

以两次刀具更换间隔为一个周期，该模型的目标函数应为刀具运行周期内的零件平均损失。

假设为零件数，刀具寿命X的分布函数，刀具的平均寿命为，刀具发生故障后带来的损失为，未发生故障时采取预防性刀具更换策略的花费为，则总费用期望：

零件的平均损失费用：

一个周期内刀具连续工作的平均长度为：

则目标函数为：

其中，。

使得上式（11）取得最小值的即为最优的刀具更换间隔。 5.3检查间隔模型

无论工序正常时产出不合格零件，还是工序故障产出合格零件，均会影响检查结果。工作人员通过检查一个零件的方式判断车床是否出现故障，需考虑误判的可能性。故在不同假设条件下建立检查间隔模型。 5.3.1问题一的检查间隔模型

根据问题一的条件，工序故障时产出的零件均为不合格品，工序正常产出的均为合格品。

若刀具更换间隔内检查出故障，费用分为三部分：检查费、故障时产出的零件损失费、发现故障进行调节使恢复正常的费用；

若刀具更换间隔内未检查出故障，费用分为两部分：检查费、未发现故障时更换一把新刀具的费用。

假设为车床连续加工零件个数，刀具更换间隔为，检查间隔为，刀具寿命的概率密度函数为。则：

1、刀具更换的间隔内检查出故障时损失费用为

根据假设（5）故障在检查的间隔内均匀分布。故当检查出故障时，平均故障运行时间为，故障情况下加工零件的损失费用为：代入零件损失费元/次，则。

设为刀具的工作寿命，检查次数为，检查费用为代入检查费用为元/次，则。发现故障进行调节使回复正常的平均费用元/次。

刀具发生故障后损失为：

其中，为加工零件的数量，为刀具的工作寿命。 2、刀具更换间隔内未发生故障损失费用

检查费用为代入元/次，则，未发生故障更换一把新刀具的费用为。 刀具未发生故障损失为：

其中，为加工零件的数量，为刀具的工作寿命。 3、目标函数

刀具运行一个周期内损失期望即为目标函数：

其中：

为加工零件的数量，为刀具的工作寿命。

在车床实际加工零件过程中，刀具工作时间越长越容易发生故障。为了及时发现故障减少损失，检查间隔应随着加工零件数量的增多而减小。为了简化模型，方便工作人员实际操作，故选择固定检查间隔的方式。

因此，使得上式（16）取得最小值的即为最优的检查间隔。

5.3.2问题二的检查间隔模型

根据问题二的条件，工序正常时产出的零件不全是合格品，有2%为不合格品；而工序故障时产出的零件有40%为合格品，60%为不合格品。

工作人员采取固定检查间隔对车床检查，当检查出一个不合格零件即视为机器故障。由于检查结果可能会导致工作人员对车床运行状态的误判。故将工序的运作状态和检查结果分为四种情况：

1、刀具更换间隔内工序发生故障，并在该段时间内检测到工序发生故障

由于检查时存在误判的可能性，在故障发生后工作人员可能无法第一时间检查到故障发生。故检测到故障的时刻与故障发生的时刻之间的间隔会影响损失期望值。

图3 ：故障点与检查点位置示意图

假设在第i-1次检查与第i次检查之间，车床发生故障。由假设(5)知，故障在此区间内发生的概率是均等的，则该段区间内平均故障时间为。

假设在第i次检查之后的每次检查中，零件不合格的概率为60%。且每次检查相互独立，检查到零件不合格就停止检查。

因此，从发生故障到检查出故障的检查次数服从参数为0.6的几何分布，其期望。故障后平均的故障运行时间为。

假设为刀具寿命，检查次数为。若，则认为该故障无法被工作人员检查出来，但根据刀具更换策略可以通过系统检修排除故障。则在工序发生故障，并正确检测到工序故障情况下，工厂损失费用为：

工序故障情况生产零件的损失费用为：，将带入得出。检查费用，发现故障后排除故障的费用元/次。

损失期望表达式为：

2、刀具更换间隔内工序发生故障，但在该段时间内未检测到故障

当工序发生故障，而被检查的零件均为合格品时，则故障在此刀具更换周期内未被检查出来。但考虑到实际情况，本文制定的刀具更换策略是通过先检查设备是否故障，确保设备无故障后再更换刀具。

因此，若在一个更换刀具的周期内故障未被检查出来，在更换刀具时该故障能被解决。

支出费用包括检查费，故障情况下产出零件损失费，在加工零件数为时，发现故障进行调节使恢复正常的平均费用。

假设刀具寿命为，故障发生后生产不合格零件的概率保持不变，则故障运行时间为。故障情况下零件损失费用为：

将参数值带入式（17）得综上可得情况二的支出费用为：

将各项费用的表达式代入上式（21）的支出费用:

3、刀具更换间隔内工序未发生故障，并在该段时间内检测到工序无故障

工序无故障时，检查的零件均为合格品，则认为工序未发生故障。但考虑到该工序正常时产出的零件不全是合格品，有2%为不合格品。则损失费用为：

检查费为，未发现故障时更换一把新刀具的费用。故 将各项费用的表达式代入上式（20）：

4、刀具更换间隔内工序未发生故障，但在该段时间内检测到工序发生故障

由于工序正常时产出的零件不全是合格品，有2%为不合格品，故工序正常工作时，错误的检测到工序发生故障的概率，损失费用主要包括检查费、工序正常而误认有故障停机产生的损失费，未发现故障时更换一把新刀具的费用。

结合误判的概率可得情况四的支出费用为：

将各项费用的表达式代入上式（）的支出费用：

综上可得，总支出费用为：

5.4改进检查方式后的检查间隔模型

检查出一个不合格零件即判断为车床出现故障，易造成误判。故可采用每次检查个零件的检查方式。若其中不合格零件数大于个，则认为工序发生故障，否则认为工序正常工作。

假设生产过程中每个零件合格与否是相互独立，每次检查一个零件花费元/次。

工序在故障情况下生产时，零件不合格的概率为60%。。检查个零件，其中不合格零件数服从参数为的二项分布，其中不合格零件大于等于个的概率为。从发生故障到检查出故障的检查次数服从参数为的几何分布，其期望。

工序在正常情况下生产时，零件不合格的概率为2%，则其中不合格零件数服从参数为的二项分布，故其中不合格零件数大于等于的概率为，从发生故障到检查出故障的检查次数服从参数为的几何分布，其期望。

将检查方式调整后，问题二中四种工作情况所对应的误判概率会发生相应的变化：

1、刀具更换间隔内工序发生故障，并在该段时间内检测到工序发生故障 损失支出调整费用为：

2、刀具更换间隔内工序发生故障，但在该段时间内未检测到故障 损失支出费用调整为：

3、刀具更换间隔内工序未发生故障，并在该段时间内检测到工序无故障 损失支出费用调整为：

4、刀具更换间隔内工序未发生故障，但在该段时间内检测到工序发生故障 损失支出费用调整为：

综上可得调整后的总支出费用为：

六、 模型求解

6.1求解刀具更换间隔

对式（11）求导可得：

其中，。

使得最小的应满足上式（33）。用MATLAB编程求解上式（33）得：。 6.2求解检查间隔模型

6.2.1问题一的检查间隔求解

为了简化模型，现假设检查间隔是等间隔检查，为定值。上式（16） 转化为函数极值问题，令可得：

故检查间隔为：

确定刀具更换间隔N=426的条件下运用MATLAB求得：当s=15时损失费最少，为1921元。

表3 ：固定刀具更换间隔条件下检查间隔的优化结果 （个） （个） （元） （元/个） 426 15 1921 4.51

6.2.2问题二的检查间隔求解

刀具更换间隔N=426时，采用定期检查方式对上式（27）进行优化，运用MATLAB解得：当s=20时损失费最少， 元。

表4 ： 固定刀具更换间隔条件下检查间隔优化结果 （个） （个） （元） （元/个） 426 20 2927 6.87 与表3结果进行对比可以发现零件的平均损失费用增加，这是因为工序在故障状态下生产零件，也会生产40%的合格品。而且检查方式为检查一个零件，故

障发生后不能及时被工作人员发现，因此增加了工厂的损失。由此可见，模型结果与实际相符。

6.3改进检查方式后的检查间隔求解

采用每次检查个零件的检查方式。若其中不合格零件数大于个，则认为工序发生故障，否则认为工序正常工作。

刀具更换间隔N=426 、s=20时，运用MATLAB求得最优的与， 元。

表5 ：固定刀具更换间隔与检查间隔条件下检查个数 （个） （个） （个） （个） （元） （元/个） 426 20 3 2 2809 6.59 七、 模型评价与改进方向

7.1模型评价

优点：检查间隔模型较符合实际实用性好，具有实际指导意义，论文中数学推导严谨，理论性强；

缺点：模型是基于先优化刀具更换间隔再优化检查间隔，未考虑两者同时对零件平均损失期望的影响。 7.2模型改进方向

由于刀具的寿命服从参数为的正态分布，概率密度函数为

则刀具的寿命的分布函数为 令设备寿命大于n的概率为。

令设备运行到时刻时仍在正常工作，则设备在内失效的概率为：

定义为设备的失效率，很小的时候，描述了在时刻仍在正常工作的条件下，在内失效的概率。

将随机变量的分布函数和概率密度代入失效率的表达式可得：

失效率的变化趋势如图（）所示，

图4 ：失效率与加工零件个数的函数图像

在刀具预防性更换间隔内，刀具的失效率单调增大。随着加工零件的数量增多，刀具损坏的概率越大，因此理想的检查方式应该是随着加工零件的数量的增多而调整检查的间隔，在刀具故障概率较小的时候检查的间隔可以较大，而在刀具故障较大的时候为了尽量快的检查到故障，检查间隔应该较小。

因此合理的检查间隔应该随着检查零件数变化，会有更好的检查效果。 假设用一次函数来描述检查间隔的变化，，由于每次检查零件的数量调整为3个，故每次检查的费用由10元/次调整为30元/次。

在工序发生故障的情况下，此时每次检查发现故障的概率，从发生故障后到发现故障时，检查次数的期望为，

在工序未故障的情况下，此时每次检查发现故障的概率，从发生故障后到发现故障时，检查次数的期望为，则四种情况损失费期望函数调整如下：

1、刀具更换间隔内工序发生故障，并在该段时间内检测到工序发生故障；

调整后的损失费用为： （40）

2、刀具更换间隔内工序发生故障，但在该段时间内未检测到故障；

调整后的损失费用为:

（41）

3、刀具更换间隔内工序未发生故障，并在该段时间内检测到工序无故障

调整后的支出费用为：

（42）

4、刀具更换间隔内工序未发生故障，但在该段时间内检测到工序发生故障

调整后的损失费用为：

（43）

综上可得总损失费用为：

（44）

八、 参考文献

[1] 徐全智、杨晋浩，数学模型，北京：高等教育出版社，2010

[2] 姜启源、谢金星、叶俊，数学模型，北京：高等教育出版社，2012 [3] 周凯、宋军全、邬学军，数学建模竞赛入门与提高，浙江：浙江大学出版社，2013

九、 附件清单

附件1：概率密度函数、分布函数、失效率函数图像绘制MATLAB程序（draw.m） 附件2：样本二阶矩求解MATLAB程序（M2.m）

附件3：刀具更换间隔求解MATLAB程序（find_N.m）

附件4：问题一检查间隔求解MATLAB程序（find_1_s.m） 附件5：问题二检查间隔求解MATLAB程序（find_2_s.m）

附件6：问题三参数、求解MATLAB程序（find_3_h_g.m）

十、 附件

附件1：概率密度函数、分布函数、失效率函数图像绘制MATLAB程序（draw.m） syms x;

f=1/(195.64*sqrt(2*pi))*exp(-(x-600)^2/(2*195.64^2)); %概率密度 F=int(f,x,0,x); %分布函数 p=f/(1-F); figure(1);

ezplot(f,[0,1200]);ylabel('f(x)');title('正态分布概率密度函数'); figure(2);

ezplot(F,[0,1200]);ylabel('F(x)');title('正态分布分布函数'); figure(3);

ezplot(p,[0,1200]);ylabel('f/(1-F)');title('失效率');

附件2：样本二阶矩求解MATLAB程序（M2.m） a=[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];

m2=sum((a-600).^2)/100 m2_root=m2.^0.5

附件3：刀具更换间隔求解MATLAB程序（find_N.m） %正态分布 syms x;

f=1/(195.64*sqrt(2*pi))*exp(-(x-600)^2/(2*195.64^2)); %概率密度 F=int(f,x,0,x); %分布函数 c1=3000;c2=1000; p=int(F,x,0,x);

y=f/(1-F)*(x-p)-F-c2/(c1-c2); yy=inline(y); solve(y) %fzero(yy,0)

附件4：问题一检查间隔求解MATLAB程序（find_1_s.m）

%正态分布，求最小效益时的s syms x; t=426;

f=1/(195.64*sqrt(2*pi))*exp(-(x-600)^2/(2*195.64^2)); %概率密度 y=inline(f);

F=quad(y,0,t); %分布函数 F1=1-F; ?(x)=1-F(x) f2=x*f; %f(x)=x*f(x) y2=inline(f2); F2=quad(y2,0,t); min=9999999; for s=1:426 x=t;

fee=(100*s+3000)*F+10/s*F2+(10*t/s+1000)*F1; if fee

附件5：问题二检查间隔求解MATLAB程序（find_2_s.m） %第二问，正态分布，求最小效益时的s syms x; t=426;

f=1/(195.64*sqrt(2*pi))*exp(-(x-600)^2/(2*195.64^2)); %概率密度 y=inline(f);

f2=x*f; %f(x)=x*f(x) y2=inline(f2);

F=quad(y,0,t); %分布函数 F1=1-F; ?(x)=1-F(x) min=9999999; for s=1:426

F2=quad(y,0,t-5/3*s); F3=quad(y2,0,t-5/3*s); F4=quad(y,t-5/3*s,t); F5=quad(y2,t-5/3*s,t);

fee1=(260*s+50/3+3000)*F2+10/s*F3; fee2=(120*t+10*t/s+3000)*F4-120*F5; fee3=(10*t/s+1000+0.02*1500*t/s)*F1; fee=fee1+fee2+fee3; if fee

mins=s; end end min mins

附件6：问题三参数、求解MATLAB程序（find_3_h_g.m） %第三问，定参数h、g t=426;s=20; %求组合

comb=zeros(s,s); for n=1:s

comb(n,n)=1; comb(n,1)=n; for r=2:n-1

comb(n,r)=comb(n-1,r)+comb(n-1,r-1); end end

%求二项分布 p1=zeros(s,s); p2=zeros(s,s); for h=1:s

for g=1:h i=[g:h];

p1(h,g)=sum(comb(h,g:h).*(0.6.^i).*(0.4.^(h-i))); p2(h,g)=sum(comb(h,g:h).*(0.02.^i).*(0.98.^(h-i))); end end

%分布函数 syms x;

f=1/(195.64*sqrt(2*pi))*exp(-(x-600)^2/(2*195.64^2)); %概率密度 y=inline(f);

f2=x*f; %f(x)=x*f(x) y2=inline(f2);

F=quad(y,0,t); %分布函数 F1=1-F; ?(x)=1-F(x) %求最优解 min=9999999; for h=1:s

for g=1:h

if (p1(h,g)>0.6)&&(p2(h,g)<0.02) Y=1/p1(h,g); F2=quad(y,0,t-Y*s); F3=quad(y2,0,t-Y*s);

F4=quad(y,t-Y*s,t); F5=quad(y2,t-Y*s,t);

fee1=(120*(0.5+Y)*s+10*h*Y+3000)*F2+10*h/s*F3; fee2=(120*t+10*h*t/s+3000)*F4-120*F5; fee3=(10*h*t/s+1000+p2(h,g)*1500*t/s)*F1; fee=fee1+fee2+fee3; if fee

本文来源：https://www.bwwdw.com/article/ye7p.html