Predictive Maintenance Scheduling with Failure Rate Described by Truncated Normal Distribution

Abstract

The method of risk assessment and planning of technical inspections of machines and optimization of production tasks is the main focus of this study. Any unpredicted failure resulted in the production plans no longer being valid, production processes needing to be rescheduled, costs of unused machine production capacity and losses due to the production of poor-quality products increase, as well as additional costs of human resources, equipment, and materials used during the maintenance. The method reflects the operation of the production system and the nature of the disturbances, allowing for the estimation of unknown parameters related to machine reliability. The machine failure frequency was described with the normal distribution truncated to the positive half of the axis. In production practice, this distribution is commonly used to describe the phenomenon of irregularities. The presented method was an extension of the Six Sigma concept for monitoring and continuous control in order to eliminate and prevent various inconsistencies in processes and resulting products. Reliability characteristics were used to develop predictive schedules. Schedules were assessed using the criteria of solution and quality robustness. Estimation methods of parameters describing disturbances were compared for different job shop scheduling problems. The estimation method based on a maximum likelihood approach allowed for more accurate prediction of scheduling problems. The paper presents a practical example of the application of the proposed method for electric steering gears.

1. Introduction

Six Sigma is a popular technique used to eliminate and prevent inconsistencies in products, processes, and to describe and reduce process failures. For example, Erginel and Hasırcı [1] controlled the failure rate of the screwing backplate to the side panels of the product. They identified failure modes that affected the failure rate of the tightening process. They used the normal distribution and process knowledge to reduce the failure rate.

In line with Six Sigma to stabilize the process, all equipment must be in the best possible technical state, the equipment must be maintained, and the workforce must be trained in operating standards. One of the conditions for the efficient implementation of this methodology is the existence of a real-time data collection system, monitoring key points of the machine or device, where symptoms of impending problems can be detected. In the past, data on the condition of the machine during its operation were collected manually, mainly by the machine operator, using one’s subjective senses and experience, and communicated in the form of oral or written reports [2]. Today, devices can and should be equipped with a set of sensors monitoring both the state of the process (correct material feeding, technological parameters within the assumed range) and the machine itself (noise, vibrations, position of key machine elements, uniformity of movement, etc.). Machine condition checking sensors can be installed by the manufacturer in places provided by the machine designer or added later, in accordance with the selected scope of the machine and process condition monitoring. For machines equipped with an advanced control system, machine health and performance data can be collected directly from this control system, often equipped with standard interfaces. The data collection itself is only the initial stage of the process, during which the data should be filtered and interpreted in order to obtain more synthetic and easier to use information, allowing for operational decisions to be made. The issues of the systemic approach to data acquisition were analyzed, among others, in [3]. In order to stabilize the course of processes and minimize their variability, it was therefore necessary to use the existing sensors through the control system, installation of additional sensors (vibration sensors, accelerometers, etc.), and communication interfaces, as well as the implementation of appropriate data processing and integration procedures, along with convenient sharing of the results in a standardized form. The collected information can be used for online analyses during the operation of the device, as well as for the subsequent offline optimization of process parameters.

Given the data on the condition of the equipment and failure-free operation times, it is possible to describe the condition of the machine by the truncated normal distribution. The causes of possible deviations from the estimated time of failure-free operation are monitored. For the predictions to be reliable, it is necessary to use an appropriate method of estimating the parameters of the truncated normal distribution, which best describe the given failure phenomena [4]. With reliable normal distributions, one can begin the improvement process and the reduction of variation in the system. The presented method is thus an extension of the Six Sigma concept for monitoring and continuous control in order to eliminate and prevent inconsistencies in processes and resulting products. In this paper, three methods of estimating the parameters of the truncated normal distribution to describe the machine failure phenomenon were presented. Reliability of the estimation methods was compared for various job shop scheduling problems.

1.1. Complexity of Scheduling Problems

Production scheduling problems are considered difficult and important issues due to the high computational complexity. Computational complexity increases exponentially with the number of input data typical for a given problem. The complexity of the scheduling problem on parallel machines with sequence-dependent setup times depends on data such as the number of jobs, setup, and production operations or machines [5]. The complexity of the construction project scheduling problem increases with the growing amount of data on various labor productivity, adequacy of equipment, and weather conditions [6]. The problem of planning multi-process production with variable renewable integration depends on the model of energy consumption (energy sources), the number and size of buffers, parameters related to the production line and the production task [7]. The complexity of the berth scheduling problem at marine container terminals depends on the number and variability of input data on services and marine container number and routes [8]. The complexity also grows due to the multidimensionality of a solution space, which depends on a number of criteria adopted to evaluate schedules [5]. The high complexity of the problem causes the need to develop an effective algorithm, both in terms of computing time and power. Computing time and power demand can increase with the complexity of the planning problem, the number of procedures, and tuning parameters used in the algorithms. In addition, when planning production processes, not all input data is known in advance, e.g., machine unavailability times [9], renewable energy sources uncertainty [7], uncertainty of labor productivity, and weather conditions [6]. In the assumed planning horizon, there are disturbances (machine failures) that make the adopted schedule no longer valid. Determining the correct date for a machine inspection to protect the system from variation is still problematic. Too often performed preventive actions are treated as constraints in the production system and lead to a reduction of performance indicators. On the other hand, the machine may suddenly fail despite preventive maintenance, which implies the need to re-optimize the schedule.

1.2. Production and Maintenance Planning Practices

Machine maintenance reduces the frequency of unpredicted production stoppages. The tasks of maintenance personnel include following activities [10]: Corrective maintenance (CM), preventive maintenance (PM), and predictive maintenance (PdM). Corrective actions consist of repairing or replacing damaged machine elements [11]. Preventive actions consist of planning maintenance activities to prevent system failures [12] based on dynamic evaluation of components’ quality, degradation state, and the quality-reliability chain of a production system. Predictive activities may focus on “predicting the date of the next failure based on previous experiences instead of using online prediction” [13]. Better knowledge of the “nature” of disturbance and more reliable predictions of the downtime evolution over time are possible by analyzing historical data on the failure-free operation of components or machines.

Most studies in the field of scheduling production tasks and technical inspections of machines assume that machine inspection time or parameters describing the state of the machine are known at the beginning of the planning period [14,15,16]. The researchers ignore the variability aspect and focuses on the deterministic approach. Researchers also treat a machine inspection as a limitation of availability. However, the process of estimating the maintenance start time is stochastic. Each maintenance operation is scheduled for a specific time which can be described by fuzzy numbers and also depends on the impact on task completion dates [17]. Similarly, Cullum et al. [18] and Bali and Labdelaoui [19] consider the schedule of technical inspections as an optimization problem [18,19]. To select the appropriate inspection date, the cost function that describes the problem (e.g., inspection costs, repair costs, and risk) is minimized. However, it is difficult to quantify the cost of maintenance if the failure type and circumstances are not known. There are also planning methods that use the existing “time windows” to perform the technical inspection of the machine [20]. This approach is in line with the trend of eliminating the “muda”, but is not possible to apply for the bottleneck. Still, others introduce additional times (safety buffers) to the operation duration to minimize the impact of disruptions on a schedule [21]. This approach lowers the performance indicators of the production system. An interesting approach describes the state of a machine using a reliability function. Machine availability deteriorates at a predetermined maximum repair rate. Including the availability component in the objective function places a higher burden on assigning maintenance at intervals that increase the availability of machines. The best allocation of production and maintenance tasks for machines is searched in terms of criteria such as: Makespan or total lateness. After maintenance, the machine is restored to an “as new” [9] or “as bad as old” state [22]. Blokus and Kołowrocki determined a reliability function for an aging series system with the component dependency following the equal load-sharing rule [23]. Wand et al. proposed the reliability equation with three failure modes: Catastrophic (binary state) failure, degradation (continuous processes), and failure due to shocks (impulse processes). Two effects of shocks on performance were considered: A sudden increase in the failure rate and a direct random change in the degradation [24]. Kleiner et al. proposed a strategy for optimal rehabilitation of network pipes in which the water distribution network economics and hydraulic capacity were included. Deterioration of pipes was modeled using aging and stress causes [25]. Kołowrocki and Soszyńska–Budny proposed approaches for multistate systems with ageing components, changes to their structure, and their components’ reliability and safety parameters during the operation processes [26]. Neelacantan et al. modeled the deterioration of the quality of pipes as a result of aging. They noticed the relationship between the pipe diameter and the uptime and included it in the repair cost optimization model [27]. Romaniuk proposed a model for assessing the cost of maintaining the water supply network, where the quality and number of pipelines’ previous failures had a direct impact on reliability [28]. Song et al. analyzed multi-component systems with each component experiencing multiple failure processes due to exposure to degradation and shock loads. Each component of the model could fail and affect all components, potentially causing them to fail more frequently. In the model, the age replacement policy and an inspection-based maintenance policy were applied for the components. The optimal replacement interval or inspection times were searched in order to minimize a cost rate function [29].

Sakib and Wuest [30] noticed the popular practice for engineers and researchers to monitor historical data, model, and simulate using failure probabilities to predict failure-free time and system deterioration over the last decades. Following this spirit, the theory of probability has been proposed to support the data analysis for predictive inspections, machine condition, and remaining usage time or life-cycle management. Historical data on failure-free machine operation times is described by a normal distribution which is commonly used to describe parametric data. Most of the causes of variation are repetitive, and can thus be described by the truncated normal distribution. The normal distribution is symmetrical in relation to its expectation. It is therefore particularly suitable for modeling phenomena in which the probability of occurrence of values above and below the average is similar. For non-negative variables, the truncated normal distribution can be used. It allows for an unequal distribution of the probability mass below and above the expected value, and thus to obtain a distribution with right-hand asymmetry. This distribution, of course, remains a two-parameter one, which gives the possibility of its good fitting to the empirical data. The frequency of the machine failure is described by the normal distribution truncated to the positive half-axis as it describes changing the probability of failure with time.

1.3. Goals and Approaches

Compared to the analyzed literature on joint scheduling of production and maintenance tasks, the distinct nature of the work presented in this paper consists of combining the four: First, the truncated normal distribution was used to model the frequency of machine failures. More precisely, it was assumed that the successive failure-free times were distributed randomly according to the normal distribution truncated to the positive half-axis.

Second, this article is a response to the need to search for methods of estimating the parameters of the truncated normal distribution. Problems with the complementarity and credibility of historical data appeared. Therefore, the estimation method based on the reliability theory used only information about the number of disturbances in historical periods.

Third, the bottleneck is the machine that has the greatest impact on system capacity [31]. The predictive schedule served as an overall plan in the event of the bottleneck failure. The function of predictive planning was to guarantee stable and reliable operation of the production system in the event of disturbances. Thus, schedules were assessed using two criteria: Robustness of the solution and reliability of quality.

Fourth, a practical example of the application of the proposed maintenance planning method for an automotive company was provided.

Based on the review of reference publications, the following research points were identified:

• Methods of achieving the reliability parameters of the truncated normal distribution, even in the case of the absence of complementary and reliable data on historical failure-free times;
• Methods for obtaining the best maintenance and production schedules where the goal is to maximize stability and robustness.

The rest of the paper is organized as follows: The model of failures is presented in Section 2. Reliability characteristics for truncated normal distribution are described in Section 2.4. The maintenance and production scheduling method is described in Section 3. The job shop scheduling problem is described in Section 4. A comparison of the results of experimental tests of three estimation methods of the truncated normal distribution parameters is presented in Section 5. In Section 5.1, the example of estimation of reliability characteristics for the production line of electric steering gears is analyzed. The paper ends with a brief summary of the results and future research objectives (Section 6).

2. A Model of Failures

We considered a production planning model taking into account machine failure frequency. It was assumed that successive failure-fee times truncated normal distributions followed by repair times described by exponential distributions. It was assumed that the parameters of the distributions generally changed with time. Reliability characteristics predictions were built based on historical data on the number of machine failures and uptime in a certain number of historical periods of the same length. Such a system was monitored on r successive time periods, as can be seen in Equation (1):

$$\left[0,\left.T\right),\right[T,\left.2T\right),\dots ,\left[\left(r-1\right)T,rT\right)$$

(1)

of the same durations, for which the information about numbers of detected failures or failure-free times was known. The prediction of system behavior was built for the next period [rT, (r + 1)T]. We assumed that failure-free times X_i,1,…,X_i,Ni in the i-th period [(i − 1)T, iT], i = 1, …, r + 1, had normal distributions with parameters m ∈ R and σ > 0, truncated to the positive half axis. The value N_i represents the number of failures (being, in general, a random variable) detected in the i-th period. It is worth noting that the cumulative distribution function (CDF) F(⋅) of such a distribution has the Equation (2):

$$F\left(t\right)=P\left\{X<t\left|X>0\right.\right\}=\frac{P\left\{0<X<t\right\}}{P\left\{X>0\right\}}=\frac{\mathsf{\Phi }\left(\frac{t-m}{\sigma }\right)-\mathsf{\Phi }\left(-\frac{m}{\sigma }\right)}{1-\mathsf{\Phi }\left(-\frac{m}{\sigma }\right)},t>0,$$

(2)

where $\mathsf{\Phi }\left(\cdot \right)$ stands for the CDF of the standard normal distribution and X is normally distributed with parameters m ∈ R and σ > 0. Differentiating the expression above with respect to variable t, we concluded that the probability density function (PDF) f_i(⋅) of the arbitrary random variable X_i,j, where i = 1, …, N_i, had the following Equation (3):

$$f_i\left(t\right)=\frac{1}{{\sigma }_i\left[1-\mathsf{\Phi }\left(-\frac{m_i}{{\sigma }_i}\right)\right]}\varphi \left(\frac{t-m_i}{{\sigma }_i}\right),t>0,$$

(3)

where $\varphi \left(t\right)=\frac{1}{\sqrt{2\pi }}{e}^{-t^2/2}$ denotes the PDF of the standard normal distribution, and m_i and σ_i are parameters depending on the number of period i.

The mean value E(X_i,k) and the variance Var(X_i,k) of the normal distribution with parameters m_i and σ_i truncated to the positive half axis, with the PDF defined in (3), are given, respectively, by the following Equation (4):

$$E\left(X_{i,k}\right)=m_i+{\sigma }_i\frac{\varphi \left(-\frac{m_i}{{\sigma }_i}\right)}{1-\mathsf{\Phi }\left(-\frac{m_i}{{\sigma }_i}\right)},$$

(4)

$$Var\left(X_{i,k}\right)={\sigma }_i^2\left[1-\frac{\frac{m_i}{{\sigma }_i}\varphi \left(-\frac{m_i}{{\sigma }_i}\right)}{1-\mathsf{\Phi }\left(-\frac{m_i}{{\sigma }_i}\right)}-{\left(\frac{\varphi \left(-\frac{m_i}{{\sigma }_i}\right)}{1-\mathsf{\Phi }\left(-\frac{m_i}{{\sigma }_i}\right)}\right)}^2\right],$$

(5)

where k = 1, …, N_i.

At the end of reliable work period X_i,k, as the failure occurs, a repair time Y_i,k begins immediately, and so on. Repair times Y_i,1, …, Y_i,Ni for i = 1, …, r + 1, are supposed to be exponentially distributed with PDFs g_i(⋅) as seen in Equation (6):

$$g_i\left(t\right)=\left\{\begin{array}{l}{\alpha }_i\exp \left(-{\alpha }_{i}t\right),t>0,\\ 0,t\le 0.\end{array}\right.$$

(6)

As it is well known [29]:

$$E\left(Y_{i,k}\right)=\frac{1}{{\alpha }_i},$$

(7)

$$Var\left(Y_{i,k}\right)=\frac{1}{{\alpha }_i^2},$$

(8)

where σ_i > 0 is known for i = 1, …, r + 1, E(Y_i,k)—the mean value of Y_i,k, and Var(Y_i,k)—the variance of Y_i,k.

We took certain simplifying assumptiosn that each new period of the form [(i − 1)T, iT] started with the beginning of reliable work X_i,1; in other words, we “deleted” the residual repair time Y_i−1,Ni in the i-th period [(i − 1)T, iT]. Thus, we can write:

$${\displaystyle \sum _{k=1}^{N_i}{Z}_{i,k}}={\displaystyle \sum _{k=1}^{N_i}\left(X_{i,k}+Y_{i,k}\right)\approx T,i=1,\dots ,r+1.}$$

(9)

Random variables X_i,k, Y_i,k, for i = 1, …, r + 1, and k = 1, …, N_i are supposed to be totally independent. Thus, the evolution of the system can be observed on successive cycles Z_i,k = X_i,k + Y_i,k, i = 1, …, r + 1, k = 1, …, N_i which are independent random variables with PDFs defined as follows:

$$\begin{gathered} h_i\left(t\right)={\displaystyle \underset{0}{\overset{t}{\int }}{f}_i}\left(t-y\right)g_i\left(y\right)dy={\displaystyle \underset{0}{\overset{t}{\int }}{g}_i}\left(t-y\right)f_i\left(y\right)dy=\frac{{\alpha }_i\exp \left(-{\alpha }_{i}t\right)}{\sqrt{2\pi }{\alpha }_i\left[1-\mathsf{\Phi }\left(-\frac{m_i}{{\sigma }_i}\right)\right]}{\displaystyle \underset{0}{\overset{t}{\int }}\exp \left[-\frac{{\left(y-m_i\right)}^2}{2{\sigma }_i^2}+{\alpha }_{i}y\right]dy} \\ =\frac{{\alpha }_i\exp \left(-{\alpha }_{i}t\right)}{\sqrt{2\pi }{\alpha }_i\left[1-\mathsf{\Phi }\left(-\frac{m_i}{{\sigma }_i}\right)\right]}{\displaystyle \underset{0}{\overset{t}{\int }}\exp \left[-\frac{{\left[y-\left(m_i+{\sigma }_i^2{\alpha }_i\right)\right]}^2}{2{\sigma }_i^2}+\frac{{\sigma }_i^2{\alpha }_i^2}{2}+m_i{\alpha }_i\right]dy} \\ =\frac{{\alpha }_i\exp \left(-{\alpha }_{i}t\right)}{\sqrt{2\pi }{\alpha }_i\left[1-\mathsf{\Phi }\left(-\frac{m_i}{{\sigma }_i}\right)\right]}{\displaystyle \underset{0}{\overset{t}{\int }}\exp \left[-\frac{{\left(y-m_i\right)}^2}{2{\sigma }_i^2}+{\alpha }_{i}y\right]dy} \end{gathered}$$

(10)

and CDFs as seen in Equation (11):

$$H_i\left(t\right)={\displaystyle \underset{0}{\overset{t}{\int }}{h}_i}\left(u\right)du=\frac{{\alpha }_i\exp \left[-{\alpha }_i\left(t-\frac{{\alpha }_i{\sigma }_i^2}{2}-m_i\right)\right]}{{\sigma }_i\left[1-\mathsf{\Phi }\left(-\frac{m_i}{{\sigma }_i}\right)\right]}{\displaystyle \underset{0}{\overset{t}{\int }}du{\displaystyle \underset{0}{\overset{u}{\int }}\varphi \left(\frac{y-\left(m_i+{\alpha }_i{\sigma }_i^2\right)}{{\sigma }_i}\right)}dy}.$$

(11)

2.1. Maximum Likelihood Approach

Suppose that in each of the intervals [0, T], [T, 2T], …, [(r − 1)T, rT], we have given sample values of random variables X_i,k, k = 1, …, N_i for any i = 1, …, r + 1, where n_i is the observed value of N_i. Thus, we have the following observations:

$$\left(x_{1,1},x_{1,2},\dots ,x_{1,n_1}\right),\dots ,\left(x_{r,1},x_{r,2},\dots ,x_{r,n_r}\right)$$

(12)

Let us consider the first historical period $\left[0,T\right)$. The maximum likelihood principle can be used to estimate parameters m_i and σ_i of truncated normal distribution as data x_1,1, x_1,2, …, x_1,n1 which are independent and identically distributed random variables. We defined the likelihood function as follows [32,33]:

$$\begin{array}{l}L\left(m_1,{\sigma }_1\left|x_{1,1},x_{1,2}\dots ,x_{1,n_1}\right.\right)={\displaystyle \prod _{k=1}^{n_1}{f}_1\left(x_{1,k}\right)}={\displaystyle \prod _{k=1}^{n_1}\frac{1}{{\sigma }_1\left[1-\mathsf{\Phi }\left(-\frac{m_1}{{\sigma }_1}\right)\right]}\varphi \left(\frac{x_{1,k}-m_1}{{\sigma }_1}\right)}=\\ =\frac{1}{{\left\{\sqrt{2\pi {\sigma }_1\left[1-\mathsf{\Phi }\left(-\frac{m_1}{{\sigma }_1}\right)\right]}\right\}}^{n_1}}\exp \left[-\frac{1}{2{\sigma }_1^2}{\displaystyle \sum _{k=1}^{n_1}{\left(x_{1,k}-m_1\right)}^2}\right].\end{array}$$

(13)

From Equation (13), we obtain:

$$\ln L\left(m_1,{\sigma }_1\left|x_{1,1},x_{1,2}\dots ,x_{1,n_1}\right.\right)=-n_1\left\{\ln {\sigma }_1+\ln \left[1-\mathsf{\Phi }\left(-\frac{m_1}{{\sigma }_1}\right)\right]+\frac{1}{2}\ln 2\pi \right\}-\frac{1}{2{\sigma }_1^2}{\displaystyle \sum _{k=1}^{n_1}{\left(x_{1,k}-m_1\right)}^2}$$

(14)

Now, differentiating on variables m₁ and σ₁, we obtain the following system of normal equations:

$$-\frac{n_1}{{\sigma }_1\left[1-\mathsf{\Phi }\left(-\frac{m_1}{{\sigma }_1}\right)\right]}\varphi \left(-\frac{m_1}{{\sigma }_1}\right)+\frac{1}{{\sigma }_1^2}{\displaystyle \sum _{k=1}^{n_1}\left(x_{1,k}-m_1\right)}=0,$$

(15)

$$-\frac{n_1}{{\sigma }_1}+\frac{n_1{m}_1}{{\sigma }_1^2\left[1-\mathsf{\Phi }\left(-\frac{m_1}{{\sigma }_1}\right)\right]}\varphi \left(-\frac{m_1}{{\sigma }_1}\right)+\frac{1}{{\sigma }_1^3}{\displaystyle \sum _{k=1}^{n_1}{\left(x_{1,k}-m_1\right)}^2}=0.$$

(16)

Substituting Equation (15) into Equation (16), we obtain:

$$-\frac{n_1}{{\sigma }_1}+\frac{m_1}{{\sigma }_1}\cdot \frac{1}{{\sigma }_1^2}{\displaystyle \sum _{k=1}^{n_1}\left(x_{1,k}-m_1\right)}+\frac{1}{{\sigma }_1^3}{\displaystyle \sum _{k=1}^{n_1}{\left(x_{1,k}-m_1\right)}^2}=0.$$

(17)

and hence:

$${\sigma }_1^2=\frac{m_1{\displaystyle {\sum }_{k=1}^{n_1}\left(x_{1,k}-m_1\right)+}{\displaystyle {\sum }_{k=1}^{n_1}{\left(x_{1,k}-m_1\right)}^2}}{n_1}=\frac{{\displaystyle {\sum }_{k=1}^{n_1}{x}_{1,k}^2-m_1{\displaystyle {\sum }_{k=1}^{n_1}{x}_{1,k}}}}{n_1}.$$

(18)

Now, multiplying Equation (15) by ${\sigma }_1^2\left[1-\mathsf{\Phi }\left(-\frac{m_1}{{\sigma }_1}\right)\right]$ we obtain:

$$-n_1{\sigma }_1\varphi \left(-\frac{m_1}{{\sigma }_1}\right)+\left[1-\mathsf{\Phi }\left(-\frac{m_1}{{\sigma }_1}\right)\right]{\displaystyle {\sum }_{k=1}^{n_1}\left(x_{1,k}-m_1\right)=0}.$$

(19)

Introducing σ₁ given by Equation (18) into Equation (19) we eliminated m₁ numerically, using one of the approximations of the cumulative distribution function and the probability density function of the standard normal distribution. Indeed, using the Maclaurin expansion, we obtain:

$$\varphi \left(x\right)\approx \frac{1}{\sqrt{2\pi }}\left(1-\frac{x^2}{2}+\frac{x^4}{8}\right).$$

(20)

Similarly, in [34] the following logistic approximation for Φ(x) was proposed:

$$\mathsf{\Phi }\left(x\right)\approx \frac{1}{1+e^{-1.702\cdot x}}.$$

(21)

After finding estimators ${\widehat{m}}_1,\dots ,{\widehat{m}}_r$ and ${\widehat{\sigma }}_1,\dots ,{\widehat{\sigma }}_r$ we extrapolated values ${\widehat{m}}_{r+1}$ and ${\widehat{\sigma }}_{r+1}$ for the period [rT, (r + 1)T] for which we had no observation, using the regression method.

2.2. Empirical Moments Approach

Suppose that we had given sample values x_1,1, x_1,2, …, x_1,n1 for period [0, T] as in the previous section. We introduced the sample mean and the variance as follows [32,33]:

$${\overline{x}}_1=\frac{1}{n_1}{\displaystyle {\sum }_{k=1}^{n_1}{x}_{1,k},s^2\left(x_1\right)=}\frac{1}{n_1}{\displaystyle {\sum }_{k=1}^{n_1}{\left(x_{1,k}-{\overline{x}}_1\right)}^2.}$$

(22)

Comparing empirical moments to theoretical ones (see Equations (4) and (5)), we obtained the following system of equations:

$${\overline{x}}_1=m_1+{\sigma }_1\frac{\varphi \left(-\frac{m_1}{{\sigma }_1}\right)}{1-\mathsf{\Phi }\left(-\frac{m_1}{{\sigma }_1}\right)}$$

(23)

and

$$s^2\left(x_1\right)={\sigma }_1^2\left[1-\frac{\frac{m_1}{{\sigma }_1}\varphi \left(-\frac{m_1}{{\sigma }_1}\right)}{1-\mathsf{\Phi }\left(-\frac{m_1}{{\sigma }_1}\right)}-{\left(\frac{\varphi \left(-\frac{m_1}{{\sigma }_1}\right)}{1-\mathsf{\Phi }\left(-\frac{m_1}{{\sigma }_1}\right)}\right)}^2\right]$$

(24)

Using the same approximation (see Equations (20) and (21)), we found estimators for ${\widehat{m}}_1$ and ${\widehat{\sigma }}_1$ from Equations (23) and (24). After finding estimators ${\widehat{m}}_1,\dots ,{\widehat{m}}_r$ and ${\widehat{\sigma }}_1,\dots ,{\widehat{\sigma }}_r$ we built predictions for values ${\widehat{m}}_{r+1}$ and ${\widehat{\sigma }}_{r+1}$ using the regression.

2.3. Renewal Theory Approach

The last proposed method was based on the renewal theory. Recall that if ξ₁, ξ₂,… were nonnegative and independent random variables with the same distribution function B(t); then the following stochastic process:

$$\nu \left(t\right)=\max \left\{n\in N:{\displaystyle {\sum }_{i=1}^n{\xi }_i}\le t\right\},t\ge 0$$

(25)

is called a renewal process generated by random variables ξ₁, ξ₂,… with renewal moments $t_n={\displaystyle {\sum }_{i=1}^n{\xi }_i},n=1,2,\dots$

Following the properties of renewal process, we expressed the first and the second moments of v(t) by means of convolutions of distribution function H₁(t) [34]; thus, we obtained:

$$E\nu \left(t\right)={\displaystyle \sum _{n=1}^{\infty }{H}_1^{n\ast }\left(t\right)}$$

(26)

$$E{\nu }^2\left(t\right)=2{\displaystyle \sum _{n=1}^{\infty }n{H}_1^{n\ast }\left(t\right)-}{\displaystyle \sum _{n=1}^{\infty }{H}_1^{n\ast }\left(t\right),}$$

(27)

where successive convolutions $H_1^{n*}\left(t\right)$ are defined as follows:

$$H_1^{1*}\left(t\right)=H_1\left(t\right),H_1^{n*}\left(t\right)={\displaystyle \underset{0}{\overset{t}{\int }}{H}_1^{\left(n-1\right)*}\left(t-y\right)d{H}_1\left(y\right)},n\ge 2.$$

(28)

We introduced Laplace–Stieltjes transforms of appropriate functions in the following way:

$${\mathsf{\Omega }}_1\left(s\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}{e}^{-st}dE\nu \left(t\right)},{\mathsf{\Omega }}_2\left(s\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}{e}^{-st}dE{\nu }^2\left(t\right)},$$

(29)

$$\widehat{H_1}\left(s\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}{e}^{-st}d{H}_1\left(t\right)},\mathrm{Re}(s)>0$$

(30)

Since the transform of convolution equals the product of transforms, then we obtained from Equation (23):

$${\mathsf{\Omega }}_1\left(s\right)={\displaystyle \underset{0}{\overset{\infty }{\int }}{e}^{-st}}d\left({\displaystyle \sum _{n=1}^{\infty }{H}_1^{n*}\left(t\right)}\right)={\displaystyle \sum _{n=1}^{\infty }{\left({\widehat{H}}_1\left(s\right)\right)}^n=\frac{{\widehat{H}}_1\left(s\right)}{1-{\widehat{H}}_1\left(s\right)}},$$

(31)

since $\left|\widehat{H_1}\left(s\right)\right|<1.$ Similarly we can prove that:

$${\mathsf{\Omega }}_2\left(s\right)=\frac{{\widehat{H}}_1\left(s\right)-{\left[{\widehat{H}}_1\left(s\right)\right]}^2}{{\left[1-{\widehat{H}}_1\left(s\right)\right]}^2}.$$

(32)

Since $\widehat{H_1}\left(s\right)$ can be calculated numerically as a function of unknown parameters m₁ and σ₁, we used a method of the Laplace or Laplace–Stieltjes in order to invert the right sides of Equations (27) and (28) on argument s. The right sides in Equations (27) and (28) are described by R₁(t, m₁, σ₁) and R₂(t, m₁, σ₁), respectively.

We can compare them to the following empirical estimators of Eν(t) and Eν ²(t), respectively:

$$n_1,s^2+n_1^2$$

(33)

where $s^2=\frac{1}{r}{\displaystyle {\sum }_{i=1}^r{\left(n_i-\overline{n}\right)}^2}$ and n₁, …, n_r are numbers of failures physically observed in successive periods [0, T], [T, 2T], …, [(r − 1)T, rT] and $\overline{n}=\frac{1}{r}{\displaystyle {\sum }_{i=1}^r{n}_i.}$ From the system of equations:

$$n_1=R_1\left(t,m_1,{\sigma }_1\right),s^2+n_1^2=R_2\left(t,m_1,{\sigma }_1\right)$$

(34)

we estimate unknown parameters m₁ and σ₁.

For the next period [T, 2T], we used the same algorithm for parameters m₂ and σ₂ with n₂ instead of n₁ as the first estimator in Equation (29). Using ${\widehat{m}}_1,\dots ,{\widehat{m}}_r$ and ${\widehat{\sigma }}_1,\dots ,{\widehat{\sigma }}_r$ we found predictions for values ${\widehat{m}}_{r+1}$ and ${\widehat{\sigma }}_{r+1}$ using the classical regression.

2.4. Predictions of Reliability Characteristics

Consider the interval [(r − 1)T, rT] for which estimators (using one of three methods) ${\widehat{m}}_{r+1}$ and ${\widehat{\sigma }}_{r+1}$ of parameters for truncated normal distribution were achieved. Below are formulas for the reliability characteristics used in the presented model:

The reliability function R(t), [(r − 1)T, rT] gives the probability that, beginning with moment t₀ = rT, the first failure occurs after time t:

$$R\left(t\right)=P\left\{X_{r+1,1}>t\right\}={\displaystyle \underset{t}{\overset{\infty }{\int }}{f}_{r+1}\left(u\right)du=}\frac{1}{{\sigma }_{r+1}\left[1-\mathsf{\Phi }\left(-\frac{m_{r+1}}{{\sigma }_{r+1}}\right)\right]}{\displaystyle \underset{t}{\overset{\infty }{\int }}\varphi \left(\frac{u-m_{r+1}}{{\sigma }_{r+1}}\right)du.}$$

(35)

The mean time to first failure (MTTFF) is as follows:

$$MTTFF={\displaystyle \underset{t}{\overset{\infty }{\int }}R\left(t\right)dt=}\frac{1}{{\sigma }_{r+1}\left[1-\mathsf{\Phi }\left(-\frac{m_{r+1}}{{\sigma }_{r+1}}\right)\right]}{\displaystyle \underset{0}{\overset{\infty }{\int }}dt}{\displaystyle \underset{t}{\overset{\infty }{\int }}\varphi \left(\frac{u-m_{r+1}}{{\sigma }_{r+1}}\right)du.}$$

(36)

The mean time between failures equals mean time to failure plus mean time of repair, or, MTBF = MTTF + MTTR:

$$E\left[X_{r+1,1}+Y_{r+1,1}\right]=m_{r+1}+{\sigma }_{r+1}\frac{\varphi \left(-\frac{m_{r+1}}{{\sigma }_{r+1}}\right)}{1-\mathsf{\Phi }\left(-\frac{m_{r+1}}{{\sigma }_{r+1}}\right)}+\frac{1}{{\alpha }_{r+1}}.$$

(37)

Probability P indicates that in the interval [a, b] ∈ [rT, (r + 1)T] there occurs at least one failure:

$$p=P\left\{a\le X_{r+1,1}\le b\right\}={\displaystyle \underset{a}{\overset{b}{\int }}{f}_{r+1}\left(u\right)du=}\frac{1}{{\sigma }_{r+1}\left[1-\mathsf{\Phi }\left(-\frac{m_{r+1}}{{\sigma }_{r+1}}\right)\right]}{\displaystyle \underset{a}{\overset{b}{\int }}\varphi \left(\frac{u-m_{r+1}}{{\sigma }_{r+1}}\right)du.}$$

(38)

The period of increased probability of failure [a, b + MTTR] (Figure 1), where a is estimated based on the assumption that the probability of the failure-free time of the machine is higher than a equals 30%, and b is estimated based on the assumption that the probability of the failure-free time of the machine is less than b equals 70%:

$$\begin{array}{l}P\left\{X_{r+1,1}>a\right\}=0,7\to 1-F_x\left(a\right)=0,7\to F_x\left(a\right)=0,3\\ \frac{\mathsf{\Phi }\left(\frac{a-m_{r+1}}{{\sigma }_{r+1}}\right)-\mathsf{\Phi }\left(-\frac{m_{r+1}}{{\sigma }_{r+1}}\right)}{1-\mathsf{\Phi }\left(-\frac{m_{r+1}}{{\sigma }_{r+1}}\right)}=0,3.\end{array}$$

(39)

$$\begin{array}{l}P\left\{X_{r+1,1}<b\right\}=0,7\to F_x\left(b\right)=0,7\\ \frac{\mathsf{\Phi }\left(\frac{b-m_{r+1}}{{\sigma }_{r+1}}\right)-\mathsf{\Phi }\left(-\frac{m_{r+1}}{{\sigma }_{r+1}}\right)}{1-\mathsf{\Phi }\left(-\frac{m_{r+1}}{{\sigma }_{r+1}}\right)}=0,7.\end{array}$$

(40)

The obtained reliability characteristics were used in building a robust schedule u, u = 1,…,w. The accuracy of the prediction was computed for the assumption that X_i,1, …, X_i,Ni had normal distributions with parameters m ∈ R and σ > 0, truncated to the positive half-axis and estimated using the described approaches (Section 2.1, Section 2.2 and Section 2.3):

$$\Psi \left(u\right)=\left|\frac{X_{r+1,1}-MTTF}{X_{r+1,1}}\right|$$

(41)

The MTTF is E[X_r+1,1] (37). The prediction of MTTF was effective if ex-post error was less than 0.05.

$$\Psi \left(u\right)\le 0.05$$

The impact of changing real failure times X_r+1,1 on the stability and robustness criteria was investigated for robust schedule u. Different values of reliability predictions achieved by the described approaches were used to generate schedules. Reliability parameters that guarantee the achievement of the most stable and robust schedules in the event of disruptions were adopted.

3. A Predictive Scheduling Model of Production and Maintenance

The process of building stable and robust schedules begins with generating a basic schedule. The scheduling problem relates to the allocation of production jobs to a limited number of machines. Having the basic schedule, the critical machine was identified. The critical machine was the one which was the most heavily loaded in the period described by reliability characteristics [a, b + MTTR]. Basic schedules were obtained using an immune algorithm [35]. Supporters of artificial intelligence methods added the ability to use information to improve the effectiveness of a solution search process. For supporters of the “non-free lunch theorem”, the proposed procedure returned schedules robust to specific types of disturbances, which better took into account the specificity of the processes under consideration. Knowledge gained about reliability characteristics improved the searching process of solutions, which was based not only on the “generate, modify, and test” principle. The best basic schedule was selected for the minimal value of the function:

$$FF\left(u\right)={\varpi }_1\cdot C\left(u\right)+{\varpi }_2\cdot F\left(u\right)+{\varpi }_3\cdot T\left(u\right)+{\varpi }_4\cdot I\left(u\right),$$

(42)

where: ϖ1, …, ϖ4 are weights of criteria, defined by a decision maker, ϖ1, ϖ2, ϖ3, ϖ4 ∈[0,1] and sum of weights is equal to 1.

C(u) is the makespan function:

$$C\left(u\right)=\max \left[t_{z,V_j}\right],$$

(43)

F(u) is the total flow time of jobs:

$$F\left(u\right)={\displaystyle \sum _{j=1}^J\left(t_{z,V_j}-t_{r,1_j}\right)},$$

(44)

T(u) is the total delay of jobs:

$$T\left(u\right)={\displaystyle \sum _{j=1}^J\left[0,D_j\right]},where\left\{\begin{array}{l}0,\hspace{1em}if\hspace{1em}{d}_j-t_{z,V_j}\le 0\\ D_j,\hspace{1em}if\hspace{1em}{d}_j-t_{z,V_j}>0,\end{array}\right.$$

(45)

I (u) is the idle time of machines:

$$I\left(u\right)={\displaystyle \sum _{l=1}^{ł}{I}_l},$$

(46)

where t_z,vj is the end time of operation v_j of job j, v_j = 1, …, V_j, j = 1, …, J, t_r,vj is the start time of operation v_j of task j, d_j is the deadline of job j, D_j is delay in completing job j, and I_l is the idle time of machine l, l = 1,…,L.

The predictive schedule was built by modifying the basic schedule using a heuristic: The minimal impact of disrupted operations on the schedule (MIDOS) [35]. According to the heuristic, a technical inspection of a machine was planned for the period described by [MTTF, MTTF + MTTR]. In the period of high risk of the machine failure [a, b + MTTR], the most flexible operations were assigned. An operation flexibility depended on: (1) The number of changes made in the event of an operation disruption ${\widehat{\nu }}_j$ performed on the critical machine, which equals the number of operations to be rescheduled, $R_{{\widehat{v}}_j}$, (2) the number of parallel machines in which the disrupted operation ${\widehat{\nu }}_j$ can be performed alternatively, $F_{{\widehat{v}}_j}$:

$$MIDO{S}_{{\widehat{v}}_j}=R_{{\widehat{v}}_j}+\left(L-F_{{\widehat{v}}_j}\right)\to \min ,$$

(47)

where L is the number of machines.

The reactive schedule was generated in a situation where the predictive schedule needed to be updated after the disturbance. The newly generated schedule should reproduce as much as possible the previous one by introducing changes caused by the delay of the job disturbed by the machine failure. The stability criterion evaluated the reproducibility of schedule u:

$$SR\left(u*\right)={\displaystyle \sum _{j=1}^J{\displaystyle \sum _{v_j=1}^{V_j}\left|s{t}_{j,v_j}\left(u\right)-s{t}_{j,v_j}\left(u*\right)\right|}}$$

(48)

where $s{t}_{j,v_j}\left(u\right)$ is start time of operation v_j of job j in predictive schedule u; $s{t}_{j,v_j}\left(u*\right)$ is the start time of operation v_j of job j in reactive schedule u*;

After the machine failure, the value of the criterion used to assess the predictive schedule should not deteriorate significantly. The robustness of schedule u is assessed by:

$$QR\left(u^{*}\right)=\left|FF\left(u\right)-FF\left(u^{*}\right)\right|,$$

(49)

where:

$$FF\left(u*\right)={\varpi }_1\cdot C\left(u*\right)+{\varpi }_2\cdot F\left(u*\right)+{\varpi }_3\cdot D\left(u*\right)+{\varpi }_4\cdot I\left(u*\right)$$

(50)

Reactive schedules u* were generated using heuristics: (1) Right shifting (RSh), and (2) re-scheduling operations on the first available parallel machine (RDO). The best schedule was selected using the rule: The minimal impact of rescheduled operations on schedule u (MIROS). Robustness QR(u*) and stability SR(u*) of reactive schedule u* were computed in the MIROS rule:

$$MIROS\left(u*\right)=0.5\left(\raisebox{1ex}{$QR\left(u*\right)$}\!\left/ \!\raisebox{-1ex}{$QR\left(u**\right)$}\right.\right)+0.5\left(\raisebox{1ex}{$SR\left(u*\right)$}\!\left/ \!\raisebox{-1ex}{$SR\left(u**\right)$}\right.\right),$$

(51)

where QR(u**)—the maximum value achieved for the quality robustness criterion, and SR(u**)—the maximum value achieved for the solution robustness criterion.

To sum up, the three stages could be distinguished in the process of developing robust and stable schedules: The construction of basic, predictive, and reactive schedules. The predictive schedules achieved were possible trade-offs between two criteria: Stability (Equation (49)) and robustness (Equation (51)).

4. Mean Time to Failure MTTF Prediction

Computer simulations are used for job shop scheduling problems. Jobs were processed in the exclude-like mode, and operations were non-preemptive. In addition, the job shop system was described by failure-free times of the critical machine for a number of historical periods (horizons) (

Table 1. Failure-free times x_i,k of the critical machine collected for scheduling horizons, i = 1,2,…,35.

The Number of Failures of the Critical Machine (and Failure-Free Times) in Scheduling Horizon i
9(90,90,100,110,100, 110,110,130,130)	7(100,105,105, 120,120,105,105)	9(110,90,105,120, 120,100,110,130,90)	6(105,120,140, 160,180,200)	10(90,130,135,90, 75,80,125,90,80,100)
11(100,90,95,95,90 ,90,85,85,90,90,80)	8(100,140,150, 140,180,100,105,80)	6(50,100,150,150, 280,270)	8(200,100,120, 130,150,90,100,90)	9(90,130,100,90, 130,120,130,90,100)
9(80,90,100,110,120, 125,120,100,110)	7(100,120,125,135, 140,150,160)	10(110,90,95,90,80, 100,90,100,90,90)	11(100,90,80,80, 85,80,90,90,85,70,70)	10(90,95,100,100, 90,150,80,90,80,80)
9(80,90,100,110,120, 125,120,120,130)	7(100,120,125,135, 140,140,160)	8(90,90,100,110,120, 120,140,160)	7(100,110,120, 145,155,150,160)	9(120,110,105,105, 110,100,100,110,120)
9(80,85,100,110,115, 110,120,120,140)	8(80,100,120,125,135, 140,140,150)	9(75,80,85,90,110, 120,130,140,150)	9(70,80,100,110, 115,120,135,130,140)	7(100,110,120, 130,140,140,180)
11(50,60,60,70,80,85, 100,110,115,120,125)	10(55,60,65,65,80,90, 110,120,140,145)	7(90,100,130,140, 150,160,170)	10(70,75,80,85,90, 95,100,110,120,120)	11(55,55,60,80,85,90, 90,100,120,130,130)
11(45,60,60,65,80,90, 100,110,115,135,135)	11(45,55,60,65,70,90, 90,110,120,120,150)	10(40,55,60,90,95, 90,110,125,135,155)	9(80,90,90,95,100, 110,125,135,155)	10(55,60,80,90,90, 90,105,125,135,140)

). The duration of each scheduling horizon i equaled 1000 h. Repair times were described by exponential distributions with given parameters α_i in period i. Suppose the predicted value of α₃₆ was 0.5, the expected distribution value (the mean repair time) was 2 h. After the machine failure, the rescheduled operations could be performed on parallel machines.

Parameters m₁, …, m₃₅ and σ₁, …, σ₃₅ of normal distributions truncated to the positive half-axis described the failure-free times. First, the parameters were estimated using the maximum likelihood method. Parameter σ₁ was numerically estimated using Equation (18) and parameter m₁ was estimated using Equations (19)–(21). After finding successive estimators m₁, …, m₃₅ and σ₁, …, σ₃₅, values of ${\widehat{m}}_{36}$ and ${\widehat{\sigma }}_{36}$ were estimated for the future period [35T, 36T] using the least squares method or Gauss–Newton method (

Table 2. Calculation of parameters ${\widehat{m}}_i$, ${\widehat{\sigma }}_i$ and reliability characteristic mean time to failure (MTTF).

	$\mathbf{The} \mathbf{Prediction} \mathbf{of} {\widehat{\mathit{\sigma }}}_{\mathit{m}+1}$$, {\widehat{\mathit{m}}}_{\mathit{m}+1} \mathbf{and} \mathit{MTTF}\mathbf{Using} \mathbf{the} \mathbf{Gauss-Newton} \mathbf{Method}$
	Maximum Likelihood Approach		Empirical Moments Approach
${\widehat{m}}_{36}$	y = 95.70 + 0.429x − 0.023x2	80.18	y = 99.25 − 0.988x + 0.009x2	76.39
${\widehat{\sigma }}_{36}$	y = 52.27 + 0.09x − 0.023x2	44.38	y = 63.27 − 1.464x + 0.028x2	47.22
The Prediction of ${\widehat{\sigma }}_{m+1}$, ${\widehat{m}}_{m+1}$, and MTTF Using the Least Squares Method
	Maximum Likelihood Approach		Empirical Moments Approach
${\widehat{m}}_{36}$	y = −0.0239x2 + 0.4291x + 95.705	80.17	y = 0.0103x2 – 1.0867x + 99.906	74.13
${\widehat{\sigma }}_{36}$	y = $0.001x^{8\cdot {10}^{-15}}$	44.41	y =0.0281x2 – 1.4997x + 63.531	45.95

). The estimated models were described using the significance test for coefficients, standard error, R², value of coefficient belonging to the interval with probability 95%, and p-value in

Table 3. Evaluation of the estimated models achieved using the Gauss–Newton method with a maximal number of iterations: 50.

$\widehat{\mathit{m}}$ Estimated Using the Maximum Likelihood Approach
	Coefficient	SE	t Statistic	R2	p-Value	95% Confidence Interval (Left End)	95% Confidence Interval (Right End)
a0	95.70475	8.430837	11.35175	0.095	0.000000	78.53170	112.8778
a1	0.42915	1.079850	0.39741		0.693702	−1.77044	2.6287
a2	−0.02390	0.029096	−0.82125		0.417581	−0.08316	0.0354
$\widehat{\sigma }$ Estimated using the maximum likelihood approach
a0	52.27978	6.464445	8.087281	0.040	0.000000	39.11214	65.44743
a1	0.09115	0.827988	0.110085		0.913030	−1.59541	1.77771
a2	−0.00862	0.022310	−0.386394		0.701762	−0.05406	0.03682
$\widehat{m}$ Estimated using the empirical moments approach
a0	99.25615	12.08471	8.213366	0.083	0.000000	74.64040	123.8719
a1	−0.98821	1.54785	−0.638439		0.527731	−4.14108	2.1647
a2	0.00981	0.04171	0.235257		0.815508	−0.07514	0.0948
$\widehat{\sigma }$ Estimated using the empirical moments approach
a0	63.27331	11.33485	5.58219	0.06209691	0.000004	40.18497	86.36164
a1	−1.46463	1.45181	−1.00883		0.320622	−4.42186	1.49260
a2	0.02830	0.03912	0.72350		0.474628	−0.05138	0.10798

. Values of estimators m₁, …, m₃₅ and σ₁, …, σ₃₅ together with the fitted functions are plotted in Figure 2.

The expected value (EX) of normal distributions truncated to the positive half-axis described by parameters achieved using the maximum likelihood method and least squares estimation method was 93.11. The increased probability of the machine failure was indicated at the time period [a, b + MTTR]. Applying Equations (39) and (40), parameters a = 58 and b = 102.28 were achieved.

Next, parameters m₁, …, m₃₅ and σ₁, …, σ₃₅ were computed using the formulas achieved for the empirical moments approach (Equations (23) and (24)). The functions describing m₁, …, m₃₅ and σ₁, …, σ₃₅ are presented in Table 2 and Figure 2. The predicted values of ${\widehat{m}}_{36}$ and ${\widehat{\sigma }}_{36}$ were equal 74.95 and 45.92 for the horizon [35T, 36T]. Mean time to failure was EX = 88.69. The increased probability of the machine failure was indicated at the time window [23, 125 + MTTR] (Equations (39) and (40)).

Experimental tests were carried out for job shop scheduling systems of different sizes in order to assess the accuracy of achieved values of reliability parameters. Each experiment was described by the number of jobs and the number of machines: (11 × 10), (10 × 9), (9 × 8), (8 × 7), and (7 × 6). Processes could have different machine routes, operation times, and deadlines in the experiments. The first machine was critical. Two sets of reliability predictions—MTTF, MTTR, a and b—were calculated using both approaches: Maximum likelihood and sample moments. Experiments were run to answer the question of what the stability and robustness of predictive schedules would be if the actual failure-free time of the critical machine X_36,1 was equal to 80, 82,… or 100 and repair-time equaled 3 h.

5. Computer Simulation Results and Discussion

Three types of schedules were generated to investigate the influence of estimation methods over the ability of achieving accurate predictions of failure-free time: Basic schedules, predictive schedules, and reactive ones. Solutions were searched for the first set of reliability characteristics: I = {a = 58, MTTF = 93, b = 102, MTTR = 3} and the second II = {a = 23, MTTF = 88, b = 125, MTTR = 3} and for changing real failure-free times.

First, the basic schedules were obtained for job shop scheduling problems: (11 × 10), (10 × 9), (9 × 8), (8 × 7), and (7 × 6) using the immune algorithm [21]. The immune algorithm was run six times for each case problem. The best basic schedules were obtained for the task sequences described in

Table 4. The quality of best basic schedules obtained for job shop problems using the MOIA (Multi-objective immune algorithm).

The Size of the Problem	The Job Sequence	Cmax	F	T	I	FFy
7 × 6	{5 2 3 4 6 0 1}	117	273	330	0	1.02
8 × 7	{5 7 2 3 4 0 6 1}	126	335	434	0	1.02
9 × 8	{8 5 2 3 0 4 7 6 1}	137	415	553	0	1.01
10 × 9	{8 5 2 3 4 6 9 0 7 1}	164	530	756	0	0.99
11 × 10	{3 8 2 5 4 7 6 9 0 1 10}	169	623	863	0	1.05

using the objective function of Equation (42).

The MIDOS rule [21] was used to generate predictions for both sets of reliability characteristics (I, II). After a critical machine failure, two RSh and RDO reschedule heuristics were applied. The best reaction schedule was selected using Equation (51). The impact of the machine uptime on two criteria was examined: The solution robustness (SR) and the quality robustness (QR) (Figure 3). The x-axis represents the actual failure-free times of the machine, X_m+1,1 = 70, 72, …, E, where E depends on the size of the scheduling problem and represents the machine end time. The y-axis represents values of criteria SR (48), QR (49), and MIROS (51) of reactive schedules, where the maximum value obtained for the quality robustness and solution robustness criterion equaled to 1, QR(u**) = 1, SR(u**) = 1.

Analyzing the impact of the machine failure time on the SR of predictive schedules generated for the first set of reliability input data for problem (7 × 6), one can notice the following phenomenon: The solution robustness of predictive schedules increased with failure-free operation of the machine until 80. After reaching the peak value, the solution robustness decreased with the failure-free time of the machine (Figure 3a). The same trends were noticeable for the quality robustness and weighted function FF (32) (Figure 3a). All failures happened before the predicted MTTF (equal 93).

Analyzing the instances of (8 × 7) and (9 × 8), the longer the failure-free time of the machine was, the better the solution robustness, quality robustness, and weighted function of the reactive schedule were (Figure 3b,c). Analyzing the instances of (10 × 9) and (10 × 11), when the real failure time approached the predicted value of MTTF, the better the solution robustness, quality robustness, and weighted function of reactive schedule were (Figure 3d,e).

The obtained schedules were also presented for the second set of reliability input data, for all instances of job shop problems (Figure 3). Analyzing the instance of (7 × 6), the longer the failure-free time of the machine was, the better the SR, QR and FF were (Figure 3a). Analyzing the instances of (8 × 7), (9 × 8), (10 × 9), and (10 × 11), when the real failure time approaches the predicted value of MTTF, the better the SR, QR, and FF of reactive schedules were (Figure 3b–e).

In most instances, the phenomenon of achieving lower values of SR, QR, and FF before and after the predicted MTTF was noticeable. This proved that the truncated normal distribution was a useful distribution for failure time modeling. The higher the accuracy of failure-free time prediction was, the lower the values of the SR, QR, and FF were obtained. The best predictions were obtained for the value of accuracy less than or equal to 0.05; that is, for real failure times from 89 to 98 in the event of first dataset and from 84 to 93 in the event of second reliability dataset (Figure 4).

The span of normal distributions describing the failure-free time is very wide in both cases of reliability input data, I (a = 58, b = 102) and II (a = 23 and b = 125). It means that the MIDOS heuristic modifies the basic schedule extensively. In order to investigate the efficiency of the MIDOS rule application, one needs to examine the real failure times before and after the time span described by the safety period [a, b + MTTR], which is not the subject of this paper.

Let us analyze the maximal, minimal, first, and third quintile of the obtained schedules for the solution robustness, quality robustness, and weighted function for the two reliability input data (RID) sets. Analyzing the instance of (7 × 6), one can notice that the reactive schedules obtained for the second RID set absorbed machine disruptions more. The first and third quintiles equaled 0 and 25, respectively, for the weighted function FF. The first and third quintiles equaled 0 and 36, respectively, for the SR and 0 and 14 for QR (Figure 5a).

Analyzing the (8 × 7) instance, the quality of the reactive schedules obtained for first and second RID sets was equal. The first and third quintiles were equal 0 and 49, respectively, for SR for the I set of the RID. The first and third quintiles were equal 0 and 4, respectively, for the QR for the II set of RID (Figure 5b).

Analyzing the instances of (9 × 8), (10 × 9), and (11 × 10), the reactive schedules obtained for the I RID set absorbed machine disruptions more (Figure 5c–e). For example, for the (9 × 8) instance, the first and third quintiles were equal 0 and 31.5, respectively, for the weighted function FF.

Analyzing the results of the simulations, it could be stated that the method based on sample moments allowed for more accurate forecasting in the event of small-size job shop problems. The maximum likelihood approach allowed for more accurate predictions in the event of large-size scheduling problems.

5.1. Estimation of Reliability Characteristics for the Automotive Industry

A production line of electric steering gears was analyzed. Successive failure-free times were supposed to have truncated normal distributions and were followed by exponentially distributed times of repairs. Based on the data about the failure number and failure-free times of the critical machine, the following two reliability characteristics were calculated: Mean time to failure (MTTF), and period of increased probability of failure [a, b + MTTR]. The empirical moments approach was used to estimate the reliability characteristics.

The critical resource measured the average fractions of the cooperating components in the electric steering gear. Data on failure-free operation of the critical resource were collected from March to October 2019. Data on failure modes and maintenance actions taken for the resource of the final steering gear functionality test were also collected. The unknown failure modes resulted from the fact that the production line was in the “infant mortality period” of the life cycle and the maintenance staff was not experienced enough.

The longest mean time between failures was in April. In April, there was only one failure; this was explained by the fact that the planned working time was relatively short, taking into account the entire period studied (Figure 6). In the following months, it was noticeable that the MTBF value stabilized, despite the increase in the planned working time and the monthly number of failures. The values of the MTBF varied from 61.33 to 78.75 h from July to October in 2019. The prediction of the failure-free time of the resource operation for December was needed.

Applying Equations (23) and (24) for the data on failure-free times, successive estimators of truncated normal distributions were obtained: m₁ = 229.23, m₂ = 0.06, m₃ = 395.95, m₄ = 170.79, m₅ = 110.97, m₆ = 199.81, m₇ = 127.54, m₈ = 104.80 h, and σ₁ = 229, σ₁ = 140.75, σ₂ = 0, σ₃ = 394.14, σ₄ = 119.33, σ₅ = 65.32, σ₆ = 213.52, σ₇ = 138.24, σ₈ = 79.66 h. The values ${\widehat{m}}_9$ and ${\widehat{\sigma }}_9$ of the truncated normal distribution were predicted using the regression method. The trend function describing m₁, …, m₈ was f(m) = 4.3469⋅9² + 28.624⋅9 + 149.44. The trend function describing σ₁, …, σ₈ was f(σ) = −6.3399⋅9² + 53.103⋅9 + 66.579. The normal distribution truncated to the positive half-axis was described by parameters ${\widehat{m}}_9$ =54.95 and ${\widehat{\sigma }}_9$ = 30.97, and the expected value equaled to EX = 57.53 for December. The technical inspection of the resource was planned at time 57 h. The increased probability of the machine failure was indicated at the time window [42, 71 + MTTR].

6. Conclusions

A method of risk assessment and planning technical inspection of machines with optimization for planning production tasks based on the probability theory was presented. The aim was to develop a universal method which reflects the operation of the production system and the nature of the disturbances, allowing one to estimate unknown parameters regarding reliability. The frequency of the machine failures was described by the normal distribution truncated to the positive half-axis. Estimation methods of truncated normal distribution parameters were presented. The estimation methods were compared for different sized job shop scheduling problems. Predictive schedules were evaluated using the solution robustness and quality robustness criteria. The estimation method based on sample moments allowed us to obtain more accurate predictions in the case of small-sized job shop scheduling problems. The method based on maximum likelihood allowed for more accurate predictions in the case of large-sized job shop scheduling problems.

Future research will be continued for the production line of electric steering gears with more than one critical resource failing. A method for estimating reliability characteristics will be given for distributions depending on the phases of life cycles of machines. It is also planned to apply other methods of statistical estimation of distribution parameters. Inspired by [36], a model of electric steering gear production system with state transitions of stochastically dependent machines has been proposed. Machine degradation rates depend on operating conditions and deteriorate from a better state to worse one.