，，，

School of Mechanical Engineering and Automation，Northeastern University，Shenyang 110819，China

（Received 14 June 2019；revised 11 September 2019；accepted 25 September 2019）

Abstract:To improve the accuracy of capacity analysis and prediction for the aircraft assembly stations，an approach for calculating the effective working hour（EWH）of automatic assembly equipment is introduced by using the dynamic mixed Weibull distribution（DMWD）model.Firstly，according to the features of aircraft assembling，a DMWD model considering the dynamic reliability of multiple subsystems and their synthetic effects on the whole equipment is established.A typical automatic drilling&riveting machine is selected as the research object，and the dynamic weights of reliability of three subsystems are modeled and solved.Subsequently the unknown parameters of the DMWD model are estimated based on maximum likelihood estimation（MLE）and Newton-Raphson method.Finally，the EWH of an automatic station is defined and modeled by using the solved dynamic reliability function.Based on the experimental study on a real automatic drilling&riveting machine from a wing panel assembly station，it is shown that the proposed DMWD and EWH models could effectively calculate the equipment reliability with full consideration of its multiple subsystems.The DMWD model is more suitable for improving the solution precision of EWH than the traditional three-parameter Weibull distribution.

Key words：automatic assembly equipment；effective working hours（EWH）；Weibull distribution；reliability

0 Introduction

In recent years，many aircraft manufacturers in China have widely used automatic assembly equipment to promote the efficiency and reduce costs.With the rapid development of aviation manufacturing technology，the precision，stability and complexity of the assembly equipment are becoming higher.However，the practical reliability and availability of automatic assembly equipment are generally low（less than 40%）［1-4］.The main reasons leading to this phenomenon are the absence of reliable methods of capacity estimating and planning for the complex equipment.Although automatic assembly equipment can significantly improve the production efficiency，increasing complexity makes it difficult to analyze the reliability and accurately calculate the effective productivity and output，and the capacity of the equipment depends on the effective reliability model.

Numerous literature studies have shown that the Weibull distribution is the most popular and the most widely used distribution in reliability of mechanical and electrical products［5-6］.The representative studies include：Castet et al.［7］conducted a nonparametric analysis of satellite reliability for 1 584 Earth-orbiting satellites based on mixed Weibull distributions with the maximum likelihood estimation（MLE）method.Djeddi et al.［8］introduced a qualitative and quantitative analysis of operational reliability for an exploited gas turbine based on three-parameter Weibull distribution.Lyu et al.［9］applied the Weibull regression model with the random effects to improve the product reliability with maximizing the mean time to failure（MTTF）.Sürücü et al.［10］con-sidered a three-parameter Weibull distribution to model inter-failure times and used a robust estimation technique to estimate the unknown parameters to monitor the reliability.Okabe et al.［11］proposed a precise method of reliability evaluation of the shortduration AC withstand voltage testby using Weibull distribution function for substation equipment.Apart from Weibull distribution，some research applied other modeling and analyzing methods to estimate the reliability for complex equipment.Beirong et al.［12］established a model for assessing the availability of complex equipment which was subject to random/wear-out failures by using generalized stochastic Petri nets and analytical method.Qiu et al.［13］studied the availability and optimal maintenance policies of a competing-risk/repairable system with multiple failure modes undergoing periodic inspections.G?rkemli et al.［14］proposed a modelling approach for determining the reliability/availability and quantifying the uncertainties of a production system with all the components and their hierarchy by using fuzzy Bayesian method.Kharoufeh et al.［15］proposed two stochastic failure models for the reliability evaluation of manufacturing equipment considering complex operating environmentand time-varying operating conditions.In some studies，fault tree analysis（FTA）technique was presented and applied to analyze and increase equipment reliability［16-17］.

In conclusion，the present study fully recognizes the contributions of aforementioned works in their respective field.However，to the best knowledge of the authors，there is little research available on automatic assembly equipment targeted analysis of capacity and reliability.It is difficult to apply these mentioned methods to this research for the following reasons.The published papers mainly focus on the reliability and availability of the whole equipment.Although some papers explored the availability ofa complex system with multiple failure modes［13］，existing reliability models rarely consider the reliability of multiple subsystems and their synthetic effects on the whole equipment；The mentioned Weibull distribution functions mostly belong to the static model which can hardly reflects the dynamic changes of system/subsystem reliability.

Based on the above analysis，in order to give a quantitative analysis of equipment capacity in this paper，we introduce the effective working hour（EWH）and adopt it as a yardstick to measure the equipment capacity.EWH is the effective output working time of automatic assembly equipment and reflects the duration of certain assembly task.Different from theoretical working hour，EWH is calculated with due consideration of equipment reliability and environmental constraints.To accurately assess EWH，a novel dynamic mixed Weibull distribution（DMWD）model is designed to revise theoretical working hour according to the characteristics of automatic assembly equipment.The main characteristic is that the equipment has multi-types of subsystems with different failure modes.This paper takes automatic drilling&riveting machine in aircraft wing panel assembly station as an example.The validity of the DWMD model is verified by calculating its effective working time.

1 Problem Description

The aim of this section is to give a detailed description of the studied problem，including the assumptions，notations，and the definition of assembly operating element（AOE）.

1.1 Assumptions

Based on the features of the aircraft assembling by automatic assembly equipment，the employment of this method is possible if：

The research subject of this paper is the EWH of automatic assembly equipment，regardless of the effect of manual operation.

Only one type of the aircraft is produced at an automatic assembly station.The difference of working hours in terms of various types is ignored.

Only one automatic assembly equipment is distributed in each station.

Assembly works in each automatic station are the single flow processes.

1.2 Notations

The notations used to define the problem are explained in Table 1.

Table 1 Notations

1.3 Definition of AOE

Generally，a typical aircraft assembly system contains a set of discrete assembly stations，and each station consists of many single flow processes［18］.In an automatic assembly station，each process needs one device to carry out the assigned tasks.In order to improve the accuracy of EWH solution，it is essential to divide processes into basic elements and regard them as statistical units.In elemental form，the operating time of complex equipment can be easily measured and standardized by gathering the beginning and ending points of the operating element.Therefore，we define the AOE as the most basic element of an automatic assembly station.Specifically，each AOE consists of several motions corresponding to a definite duration time.These motions are performed by the qualified automatic assembly equipment at one station.The operating time of an AOE contains two parts：setup time and processing time.For example，a typical AOE（drilling a rivet hole with 5 mm in diameter and 6 mm in depth）on a central wing panel usually needs 18 s to be processed by the automatic assembly equipment.The setup time is 3 s to position the rivet hole，and the processing time is 15 s to drill the rivet hole.According to assembly order（AO），the types and number of AOEs could be obtained.Once all the AOEs are identified and measured，it is convenient to calculate the theoretical working hour by accumulating these AOEs for each automatic station.

2 Dynamic Mixed Weibull Distribution Model

Due to the dynamic characteristics of complex product assembly，the system variables change over time.Thus in this section，a dynamic mixed Weibull distribution（DMWD）model considering the dynamic reliability of multiple subsystems and their synthetic effects on the whole equipment is established.

2.1 Definition of DMWD model

Based on the features and operating environment of automatic assembly，the mathematical definition and description of DMWD model will be given as follows.Firstly we setTDas the inspection interval of automatic assembly equipment，which denotes the maximum achievable number of AOEs within one maintenance cycle.In an inspection interval，the sequence numberxof AOE assembling is set as a dy-namic variable.We assume that when thexth AOE is being assembled，the probability density function（PDF）of equipment operational reliability follows a three-parameter Weibull distribution off(x)～W(m，η，γ；x)，where the shape，scale，and location parameterm，ηandγare independent.Numerically，the classical mathematic formula of PDF and reliability function［19］in terms ofxcan be given by

A mixed Weibull distribution represents a population that consists of several Weibull subpopulations.Therefore we set an automatic assembly equipment containingNtypes of subsystems，in which each subsystem has particular failure mode.Specifically，all the failure functions obey different types of Weibull distributions with dynamic weight coefficientspi(i=1，2，…，N).As for the whole equipment，the probability density of failure functionfi(x)and reliability functionR(x)conform to

In Eq.（2），there are 4N-1 unknown parameters，including 3NWeibull parametersmi，ηi，γiandN-1 dynamic weight coefficientspi.How to estimate these unknown parameters accurately is the basis of DMWD modeling.Thus we will deduce the dynamic weight coefficientspiand estimate the unknown parameters in detail.

2.2 Dynamic weight coefficient solution of subsystems

In this section，how the reliabilities of subsystems affect the whole automatic assembly equipment is analyzed and calculated.As for the aircraft assembly，automatic drilling&riveting machines comprise the majority of automatic assembly equipment.Therefore we select an automatic drilling&riveting machine as the target to analyze its composition and characteristics.In general，a typical automatic drilling&riveting machine contains three subsystems：control system，location&execution system，and measurement&feedback system.The dynamic weights of reliability of three subsystems will be modeled and solved as follows.

2.2.1 Control system

Control system is responsible for the logic control of motion mechanism.Its reliability and stability determines whether the equipment can perform assembly tasks according to the assembly order.In the initial stage of operation，system developing and debugging account for more than 70%of the total work load.Because of higher requirements for the assembly accuracy and efficiency，operators need to invest significant time in simulating，online/off-line programming and testing for control system，which will lead to error accumulation and reduce equipmentreliability.In consequence，the reliability weight of control system should be defined at a high level in the initial stage.When the equipment enters into a steady operation stage，the reliability weight will gradually decline and then get into a steady level.In order to give a quantitative description of this nonlinear decline，we setp1(x)as the dynamic weight of control system’s reliability and hypothesize thatp1(x)belongs to a two-parameter Weibull function，as shown in

Here we set the shape parametermp1＜1to ensure that the weight curve ofp1(x)follows the decreasing interval of Gamma distribution.To solve the initial value ofp1(x)，we should determine the operating time proportion of each subsystem in the initial stage.According to the statistical data of working hours（＜100 h）of three subsystems，we can deduce that the operating time proportion of control system，location&execution system，and measurement&feedback system is about 3∶1∶1.As a result，the initial value of［p1（0），p2（0），p3（0）］is［0.6，0.2，0.2］.Under the condition of pre-determinedp1(0)，withmp1varing from 0.5 to 1，the corresponding scale parameterηp1can be solved.Then the diverse dynamic weight curvep1(x)can be generated.The inspection intervalTDis set as 100 000 based on the technical requirement of aircraft assem-

bly.Fig.1 indicates that as the number ofmp1increases from 0.5 to 1，the weight curve ofp1(x)becomes relatively flatter.Whenmp1=1，p1(x)becomes the exponential distribution.To avoidp1(x)decreasing too rapidly or slowly，we setmp1=0.9，（blue curve in Fig.1）.

Fig.1 Effect of shape parameters on dynamic weight of control system

2.2.2 Location&execution system

Location&execution system constitutes the host structure of automatic assembly equipment，mainly including motion mechanism，locating&clamping mechanism，and end effector.The reliability of location&execution system depends on the quality and stability of hardware structure.After installation and configuration，the host structure is no longer modified.Thus the reliability weight of location&execution system remains relatively stable at a lower value in the earlier stage of a maintenance interval.At the later stage，the mechanical wear and aging will lower the reliability，causing the growth of the reliability weight.Similarly，we setp2(x)as the dynamic weight of location&execution system’s reliability，which follows the translation of two-parameter Weibull function.The initial weightp2(0)is figured out to be 0.2 andp2(x)conforms to

When the shape parametermp2=4，the weight curvep2(x)is close to Gaussian distribution.

The monotonically increasing interval ofp2(x)is subject to the dynamic growth of the reliability of location&execution system.Asmp2is determined，the weight curve stretches along the transverse axis and compresses along the vertical axis with the scale parameterηp2increasing from 5 to 15（see Fig.2）.Asxapproaches to the inspection interval，the upper bound of weight is set as an empirical value of 0.3.Then according to the intervalconstraintp2(x)∈[0.2，0.3]，we can solveηp2=12.931（blue curve in Fig.2）.

Fig.2 Effect of scale parameter on dynamic weight of position and execution system

2.2.3 Measurement&feedback system

Measurement&feedback system is designed to assure the required accuracy and automation by real time measurement from assembly objects.With mechanical wear and aging of equipment，the equipment stability is more dependent on measurement&feedback system.Hence，we deduce that the dynamic weightp3(x)of this system should be growing gradually through a maintenance interval.Based on the solvedp1(x)andp2(x)，p3(x)can be obtained with the constraint of ∑pi=1.Finally the dynamic changes of three dynamic weights can be generated（see Fig.3）.It can be seen that the dynamic weightp3(x)is increasing and then tending to the stability whenxis improving，which accords with the proposed deduction.

Fig.3 Dynamic weights of reliability of three subsystems

3 EWH Solution of Automatic Assembly Equipment

3.1 Parameter estimation of DMWD model

At present，the main parameter estimation methods include graphical procedures（probability plotting）， MLE， andBayesianestimation， etc.Among these methods，MLE is more suitable for multi-parameter estimation of the DMWD model because it is not only more accurate than graphical procedures but also more simplified than Bayesian estimation.Although Bayesian estimation can achieve high accuracy，the multiple parametersofthe DMWD model are difficult to estimate accurately for lack of explicit expressions of integral equation.Thus the MLE method is used to estimate the unknown parameters.Here，we defineθi(mi，ηi，γi)as the parameters of PDF for theithsubsystem.Within one maintenance interval，we assume that a failure would occur onjthAOE of subsystemi.For each subsystem，the failures are independent identically distributed and their observational results of each subsample arexij(j=1，2，…，NSi).Then the failure probability set in the interval[xij，xi(j+1)]isf(xij)dxij，and the maximum likelihood function can be given by

After taking logarithm for Eq.（5），we obtain the likelihood function

To seek the extreme value of Eq.（6），we makeand get the expressions

Because finding exact solution of nonlinear Eq.（7）is difficult by common numerical methods，in this paper we adopt the Newton-Raphson method to get solutions iteratively.Details of the scheme are given as follows.

Step 1The initial parameters to be estimated in Eq.（7）are solved by using graphic method to prevent non-convergence of iteration.

Step 2On the basis of Step 1，the problem is simplified to the parameter estimation of two-parameter Weibull function.Hence the initial shape and scale parametermi0，ηi0can be estimated rapidly.

Step 3mi0，ηi0are taken as the inputs of threeparameter Weibull distribution.Then the graphic method is employed to estimate the initial location parameterγi0.

Step 4In order to calculate the higher-order small quantities Δmi，Δηi，Δγi，we simplify Eq.（7）to the linear equations atθi0(mi0，ηi0，γi0)by means of series expansion.

Step 5The error boundε(m，η，γ)is defined based on the accuracy requirement and compared with Δmi，Δηi，Δγi.If Δmi，Δηi，Δγi＜ε(m，η，γ)，the current results should be the estimated value.Otherwiseθi0+Δθiis substituted to the next iteration instead of Δmi，Δηi，Δγiand iteration is continued until all the parameters)are obtained.

3.2 EWH modeling for automatic stations

At an automatic station，the theoretical working hour can be easily identified and measured.However，the theoretical working hour is static and does not reflect the time extension caused by the decline of reliability under actual assembling conditions.Thus we use the DMWD model to equivalently extend the theoretical working hour so as to calculate the EWH.The detailed derivation of EWH is as follows.Firstly in an automatic station，we setNEas the number of AOE types.NEk(k=1，2，…，NE)is the AOE number of thekth type.HRkis the equipment’s time consumption of switching to thekth type，and the order ofkis accordant with the processing sequence.For thekth type of AOE，the mean measurements is set ashEk.Then the number of theoretical working hourHTkandHTcan be formulated as

According to Eqs.（2）—（7），we can deduce the equipment reliabilityR(x)when assembling thexth AOE.The failure rate isF(x)=1-R(x)，and then we assume that the EWH of assembling the first aircraft at an automatic station meets

Then，we setXas the cumulative number of the finished aircrafts，and then can deduce the EWHHE(X)and the cycle timeTE(X)of assembling theXth aircraft at this station.

where the sequence numberx=l+(X-1)NES，NESis the total AOE’s number of a station andNES= ∑NEk；hDis the number of required working hours per day；「」denotes the top integral function.

4 Experimental Study

In this section，to verify the feasibility of the proposed DMWD and EWH models，a real automatic drilling&riveting machine from the station of integral wing panel assembling is chosen as an example to be tested at the Shenyang Aircraft Company.The parameters used in this paper are collected and set in Section 4.1，and the experimental analysis of the proposed methods is provided in Section 4.2.

4.1 Parameters setting

Firstly，based on the investigation conducted in the assembly workshops，we collected much statistical data to calculate the theoretical working hour of the integral wing panel assembling station.Concretely，the parameters includingNE，NEkandHRkare recorded from the assembly order and process files.The mean time consumptionhEkof each type of AOE is estimated according to the field test data.The number of theoretical working hourHTkis solved by Eq.（8）.Details of these basis data are listed in Table 2.Moreover，the total numberNESof AOE at this station is 3 329，and the total amount of theoretical working hourHT=25.1.

Table 2 Basic data of each type of AOE

Then，in order to get the unknown Weibull parametersof the proposed DMWD model，the sequence numberxof failure AOE from three subsystems is recorded within one maintenance interval（100 thousands of AOE）according to the statistical failure data，as shown in Table 3.As for the automatic drilling&riveting machine，the failure occurs when a subsystem stops performing its required function such as unexpected assembling errors，movement interference，or emergency shutdown under other fault conditions，etc.Based on the data from Table 3，we use MATLAB to estimate the unknown parameters including nine Weibull parameters and their initial values by Newton-Raphson iteration.The results are listed in Table 4.Furthermore，based on the estimated parameters， the reliability functionsR1(x)—R3(x)of three subsystems can be solved（see Fig.4）.Then we combine Eqs.（2）—（4）to calculate the dynamic weight coefficientsp1(x)—p3(x).So the reliability functionR(x)of the DMWD model for the automatic drilling & riveting machine can be solved.Meanwhile we also establish a traditional three-parameter Weibull distribution（TWD）model as the comparison and solve its reliability functionRS(x)with the parameters［0.627 6，4.204 5，0.001 6］through Table 3.The comparison of reliability resultsR(x)andRS(x)are shown in Fig.5.

Table 3 Sequence number of failure AOE from automatic drilling&riveting machine

Table 4 Estimated values of DMWD parameters

Fig.4 Distribution curves of reliability of three subsystems

Fig.5 Comparison of reliability results between DMWD model and traditional Weibull model

Fig.6 Comparison of working hours among DMWD/TWD models and the measured values

4.2 Experimental analysis

Based on the results of reliability function，the EWH models can be solved and verified.According to Eq.（10），the EWHs of assemblingXth aircraft at the integral wing panel station can be calculated.Meanwhile，we collect the measured working time consumption of assembling each corporate jet（X：1—30）at this station from the PDM/ERP systems in 2017.The comparison of working hours among the proposed method （DMWD model），TWD model and the measured values is shown in Fig.6，and the comparison of deviation values between the DMWD and TWD models is shown in Fig.7.Based on Fig.6，the EWH curve is a monotonically increasing convex function with the growth ofX，which is similar to the trend of the measured values.The EWH curve with DMWD model changes more gently than the EWH one with TWD model.Intermsoftheproposed DMWD model and the measured working hours，the similarity coefficient of two data sets reaches up to94.1%，demonstratingthattheEWH with DMWD model can simulate the gradual growing of the actual time consumption.Thus，it is shown that the DMWD model provides significant accuracy in capturing the failure trends of the automatic drilling&riveting machine.

Fig.7 Comparison of deviation values between the proposed method and TWD model

As is shown in Fig.7，the deviation values of DMWD model are obviously lower than those of the TWD one.As for the DMWD model，the absolute value of the maximum deviation is 3.991 h（within 12.4%），less than 4.392 h（15.9%）of the TWD model.Moreover，the mean values of absolute and standard deviation of the proposed method are 1.231 h and 1.541 h，respectively，less than 1.392 h and 1.703 h of the TWD model.In addition，we adopt cubic polynomial fitting to generatethedeviationcurvesoftheDMWD/TWD models.The results show that the fitting curve of DMWD deviation is smoother than that of TWD deviation and closer to zero.Therefore the proposed model is more efficient for calculating and predicting the EWH of automatic assembly stations.

However，F(xiàn)igs.6 and 7 also indicatethat when the cumulative number of the finished aircrafts is relatively small（X＜6），the amount of measured working hour is higher than EWH.On the contrary，most of the measured values are lower than EWH whenX＞6.To the best of the authors’knowledge，there are two main reasons contributing to this difference.Firstly，in the initial stage，operators need to adapt to the new equipment by intensive testing and monitoring.Meanwhile，the new equipment requires a running-in period to reach a steady state.Thus，the amount of measured working hour is higher than EWH at this stage.Secondly，as the number of assembled aircrafts increases，learning effect of operators will work and gradually decrease the switching timeHRk，which partially offsets the growth of working hour by the reliability.Thus，at this stage most of the measured working hours are slightly less than EWH.

5 Conclusions

This paper proposes a novel method for modeling and calculating the EWH for automatic assembly equipment of aircraft assembly stations，by means of reliability analysis with the DMWD model.Based on the experimental study and analysis，it is shown that the proposed method could effectively calculate the equipment reliability with full consideration of its multiple subsystems.Moreover，the DMWD model is more suitable for improving the solution precision of EWH than the traditional threeparameter Weibull distribution.This study is both theoretically and practically significant，since it provides a basis for managers to estimate and predict the capacity of complex equipment or systems.However，it is still difficult to accurately calculate the EWH of some complex operations which need human-machine cooperation.In the further study，we will focus on modeling and solving the nonlinear relationship between the complex equipment and assemblers，and find their synthetic effects on the reliability and EWH.