七、随机数据输入模型构建范例

基本过程：

1）  生成数据的直方图

2）  根据直方图与各种pdf曲线对比，选择分布簇

3）  进行参数估计

4）  进行假设检验，如果检验不通过重复步骤2）-4）

5）  获得该数据的随机分布函数

 

7.1   指数分布数据处理范例

在对银行营业部客户服务水平的仿真过程中，需要分析顾客到达时间间隔服从怎样的随机分布函数。通过秒表测时收集了连续的100位顾客到达营业厅的时刻点数据，处理得出到达时间间隔数据如下：

8.04

4.77

12.42

9.22

3.91

8.50

14.74

9.20

5.82

28.56

10.33

4.62

7.66

7.24

39.22

2.14

15.46

16.16

16.45

1.18

22.38

3.12

0.59

3.76

3.01

4.00

4.59

7.07

1.55

29.58

10.38

9.59

12.41

20.05

0.78

9.34

8.59

2.54

1.34

11.50

1.34

3.82

9.47

15.83

14.00

25.77

3.59

3.62

0.08

2.87

21.83

18.41

21.70

21.95

1.91

5.97

2.53

9.21

2.20

4.86

0.96

18.79

6.33

7.16

5.80

14.65

0.86

4.33

4.31

0.07

0.04

1.33

17.15

22.81

6.59

5.94

3.51

4.22

5.49

13.27

26.06

0.29

1.04

16.29

1.66

5.81

3.78

37.63

0.44

1.77

2.24

19.26

9.32

16.39

11.56

2.67

8.16

2.40

3.53

5.14

分析顾客到达时间间隔服从什么分布？开发并检验一个合适的模型。

求解过程：

1）  数据的直方图

使用Minitab生成数据直方图如下：

2）  选择分布簇

对照pdf曲线图和上面的直方图，可以发现，该直方图同指数分布曲线比较吻合，因此，假设顾客到达时间间隔服从指数分布，及其概率密度函数pdf为：

其中为该分布的均值。

3）参数估计，

在Minitab中获得直方图的拟合曲线和参数的估计值，见下图。

从图上可以看出，为9.018，即选择的指数分布密度函数pdf如下：

　

   （2）

概率分布函数cdf如下：

4）检验

  运用检验如下假设。

  H0：x的概率密度函数服从式（2）的指数分布

  检验计算表

No

Min

Max

fi

pi

npi

npi-fi

(npi-fi)2/npi

Fvalue

1

0.000

3.922

36

0.353

35.270

-0.730

0.015

0.353

2

3.922

7.845

21

0.228

22.830

1.830

0.147

0.581

3

7.845

11.767

15

0.148

14.778

-0.222

0.003

0.729

4

11.767

15.690

7

0.096

9.566

2.566

0.688

0.824

5

15.690

19.612

9

0.062

6.192

-2.808

1.273

0.886

6

19.612

23.535

6

0.040

4.008

-0.636

0.036

0.926

7

23.535

27.457

2

0.026

2.594

0.952

8

27.457

31.380

2

0.017

1.679

0.969

9

31.380

35.302

0

0.011

1.087

0.980

10

35.302

39.225

2

0.020

1.995

1

 =SUM(ABOVE) 2.162

 

  因为>2.162，故在水平0.05下接受H0，认为x服从均值为9.018的指数分布。

 

所以，该银行客户到达时间间隔服从均值为9.018Min的指数分布。

 

7.2   三角分布数据处理范式

  记录机床操作工作业时间数据如下，单位为秒，分析该作业时间服从什么分布？开发并检验一个合适的模型。

31.1

34.8

21.1

26.6

32.5

36.6

38.6

33.4

27.8

42.3

23.2

29.3

37.9

26.9

30.4

42.1

31.7

31.8

24.3

29.3

28.7

30.8

35.2

25.6

43.1

46.4

35.3

45.9

31.0

34.5

33.3

32.4

32.4

30.2

42.3

49.0

32.6

41.8

37.7

33.1

29.7

30.1

31.4

24.9

30.6

29.5

29.8

26.5

29.1

23.2

33.0

29.3

29.6

42.6

25.8

26.3

45.2

23.0

37.8

27.6

33.5

23.0

47.2

33.2

43.4

27.3

35.2

43.0

31.4

29.3

30.1

33.2

31.5

28.1

30.5

41.5

33.1

31.1

35.6

32.7

41.4

26.6

27.8

37.8

43.2

35.7

32.9

30.7

31.8

38.4

42.4

41.6

38.1

30.2

34.1

26.3

30.8

38.7

35.4

38.4

 

求解过程：

3）  数据的直方图

使用Minitab生成数据直方图如下：

4）  选择分布簇

对照pdf曲线图和上面的直方图，可以发现，该直方图同三角分布曲线比较吻合，因此，假设作业时间服从三角分布，及其概率密度函数pdf为：

其中a，b，c为该分布的参数，a为最小取值，b为众数，c为最大取值。

3）参数估计，

该分布的参数a和c的估计值直接通过观察数据即可获得，其中a=21.2，c=49.0。而b可由下式计算获得。

已有数据的均值E(x)为33.34，则b=29.92，即选择的三角分布密度函数pdf如下：

   （2）

概率分布函数cdf如下：

4）检验

  运用检验如下假设。

  H0：x的概率密度函数服从式（2）的三角分布

检验计算表

No

Min

Max

fi

pi

npi

npi-fi

(npi-fi)^2/npi

Fvalue

1

21.1

23.9

6

0.032

3.162

-2.350

0.437

0.032

2

23.9

26.7

9

0.095

9.487

0.126

3

26.7

29.5

12

0.158

15.812

3.812

0.919

0.285

4

29.5

32.3

13

0.189

18.890

5.890

1.837

0.474

5

32.3

35.0

17

0.161

16.087

-0.913

0.052

0.634

6

35.0

37.8

8

0.132

13.162

5.162

2.024

0.766

7

37.8

40.6

8

0.102

10.237

2.237

0.489

0.868

8

40.6

43.4

12

0.073

7.312

-4.688

3.005

0.942

9

43.4

46.2

3

0.044

4.387

-0.150

0.004

0.985

10

46.2

49.0

3

0.015

1.462

1.000

 =SUM(ABOVE) 8.767

 

 

  因为>8.767，故在水平0.05下接受H0，认为x服从参数为[21.2，29.9，49]的三角分布。

 

所以，该工序作业时间服从参数为[21.2，29.9，49]的三角分布。

 

 

 

 

八、示例1：设备故障间隔时间模型的识别

发动机曲轴生产线有5台设备，设为A/B/C/D/E，如图一所示。该产线生产两种产品P1和P2，同一机器对产品的加工时间相等，当机器加工不同产品进行换型时，各台设备的换型时间差异较大，如表一所示。设备间以滚道输送链连接并进行曲轴的输送，设备间滚道和生产线上下料滚道各可以存放6件，滚道移位时间为1Sec。其中设备D具有随机故障的发生，故障发生时刻和维修时间数据如表二所示，完整的历史数据见：输入模型之设备故障间隔时间模型识别数据.xls。为了进行仿真建模，需要识别故障间隔时间服从的分布以及维修时间服从的分布。对数据表进行分析时，假设缸体生产线每天工作时间为8小时，即从：8：00-16：00。

试建立故障间隔和维修时间的随机输入模型，设计该系统的生产排程规则，并用Witness建立该系统的仿真模型。

图一   产线布局图

表一   加工和换型时间数据表

设备

加工时间Sec

换型时间Min

A

46.8

2

B

47

14

C

47.1

2

D

47.6

80

E

45

44.3

　

表二   数据处理局部表格

8.1 设备故障间隔时间输入模型的构建

（1）间隔时间数据的获得

上表中，A-K列数据为现场记录数据，L和M列是为了计算两次故障间隔发生时间而进行的计算数据。其中L列系统时间为K列发生时刻数据转换获得的，M列故障间隔时间为C列日期和K列发生时刻数据处理获得的，具体计算函数如下。

以第2行的数据为例，L列系统时间L2计算表达式为：=ToTime(K2)，即将K列数据转换为系统时间值。

ToTime函数为自定义函数，函数体如下：

Public Function ToTime(sStr As String) As Date

   ToTime = CDate(sStr)

End Function

因为K列发生时刻数据为字符串，为了计算为字符串在每天中的时间值（一天24小时的时间值为0-1之间的小数），使用了Cdate函数进行转换，具体参看xls。

以第3行的数据为例，M列间隔时间M3计算表达式为：=CalInterval(C2,L2,C3,L3)，即根据前后两次故障发生的日期和时间，计算处两次故障发生的时间间隔。

CalInterval函数为自定义函数，函数体如下：

Public Function CalInterval(nDate1 As Integer, rTime1 As Date, nDate2 As Integer, rTime2 As Date) As Double

   If nDate2 < nDate1 Then

        CalInterval = nDate2 * 480 + (rTime2 - rTime1) * 1440

   Else

        CalInterval = (nDate2 - nDate1) * 480 + (rTime2 - rTime1) * 1440

   End If

End Function

nDate1和rTime1分别为前一次发生故障的日期和时间，nDate21和rTime21分别为当次发生故障的日期和时间。

如果前一次和当次故障都在同一个月中发生，故障间隔计算表达式为：

CalInterval = (nDate2 - nDate1) * 480 + (rTime2 - rTime1) * 1440

上述表达式表示两种情况：如果两次故障均发生在同一天，故障间隔等于故障发生时刻的系统时间之差乘以1440分钟；如果两次故障不是发生在同一天，则故障间隔等于故障发生日期之差乘以480（因为每天正常班次时间为8小时=480分钟）再加上两次故障发生时刻系统时间之差。

如果前一次和当次故障不在同一个月中发生，故障间隔计算表达式为：

CalInterval = nDate2 * 480 + (rTime2 - rTime1) * 1440

不再同一月发生的条件由表达式nDate2 < nDate1控制，因为从表中数据分析可以看出，如果当次故障同前一次故障不是在同一月发生，则当次故障发生的日期肯定要小于前次故障发生的日期。

（2）间隔时间直方图的生成

将M列的数据复制并粘贴为行的形式，然后将行数据复制到Matlab窗口，令Y=[ 98     405  207  101  142  506  2031       1208       1573       643  139  472  363  424  1166       109  124  126  104  207  198  247  403  1548       42    213  783  657  67       162  112  65    145  884  561  3437       638  728  138  190  423  117  767  210  81    58    122  186  629  2088       285  1515       1679       1689       1528       413  142  226  676  452  288  4     317  103  688  204  1979       1100       553  1197       607       1279       92]

在Matlab窗口输入命令：hist(y,9)，即可获得下列直方图。

从直方图中可以看出，该随机分布近似服从指数分布。

（3）选择分布簇

对照pdf曲线图和上面的直方图，可以发现，该直方图同指数分布曲线比较吻合，因此，假设顾客到达时间间隔服从指数分布，及其概率密度函数pdf为：

其中为该分布的均值。

（4）参数估计，

     可以计算出数据的均值为589.9，即为589.9，选择的指数分布密度函数pdf如下：

   （2）

概率分布函数cdf如下：

（5）卡方检验

卡方检验计算如下表所示。

 

检验统计量为3.871，查表>3.871，所以接受设备故障间隔时间服从均值为589.9的指数分布。

注：或者在Matlab中使用kstest(z,[z expcdf(z,589)])进行检验，其中z为故障间隔时间的向量，获得的结果也为0，即不能拒绝该分布服从均值为589的指数分布。

8.2 故障维修时间的输入模型构建

（1）数据分析和直方图的生成

从数据表二中的维修时间中可以看出，对维修时间的记录很粗糙，基本上维修时间记录为5的整数倍，因此可以将维修时间模型开发为为离散分布模型。首先识别出重大故障发生的频率，即维修时间为130Min的故障比率为1/74。然后将130的数据从数据列表中移除，这样剩下的维修时间如下：

15   20    20    20                                      

25   35    35    15    50    18    15    20    75    30

40   40    28    30    20    45    80    30    80    15

20   40    45    15    25    10    25    20    15    10

15   25    45           20    15    15    30    45    25

20   35    15    10    50    35    20    20    30    55

60   21    25    53    20    18    20    13    20    20

55   15    15    25    25    45    25    55    12    20

对其中的非5整数倍数据进行整数倍处理，处理规则类似四舍五入法，例如12将转化为10，13将转化为15。得到新的数列如下：

15   20    20    20                                      

25   35    35    15    50    20    15    20    75    30

40   40    30    30    20    45    80    30    80    15

20   40    45    15    25    10    25    20    15    10

15   25    45           20    15    15    30    45    25

20   35    15    10    50    35    20    20    30    55

60   20    25    55    20    20    20    15    20    20

55   15    15    25    25    45    25    55    10    20

在Matlab命令窗口使用Hist(Z,15)生成数列直方图如下图。

（2）经验分布表

根据直方图，很容易获得每种维修时间出现的频次，获得经验分布表如下。

维修时间

10

15

20

25

30

35

40

45

50

55

60

65

70

75

80

出现频次

4

13

19

9

6

4

3

5

2

5

1

0

0

1

2

 

8.3 建立仿真模型

 

（1）模型中的元素列表