4.1 Steps for implementing the PSO
PSO is a study of the process of searching for food in a flock of birds in nature. The algorithm simulates the behaviour of a flock of birds searching for food and the cooperation between birds to achieve the optimal purpose of the flock. The PSO is a very simple algorithm that performs effective optimization for various functions. As an optimization algorithm based on the theory of population intelligence, it operates on the principle that the population intelligence generated by the cooperative and competitive relationship between the particles in the population guides the optimization search[30]. In PSO, each particle represents a bird of the natural bird community, and the range of the bird in search of food is the search space of the particle. All particles have a fitness value, which is determined by the objective function. Velocity and position are the two basic parameters of a particle. The velocity represents the speed of the particle's flight and the position represents the direction of the flight. Each particle searches for its optimal solution in the search space and records it as \(pbest\) to share with the other particles in the population. After that, the optimal solutions found by each particle are compared and the best result is recorded as the current population optimal solution \(gbest\). The velocities and positions of all particles in the population are adjusted by \(pbest\) and \(gbest\). In this paper, algorithms are developed using MATLAB to minimize the preventive maintenance cost of machine tools based on the above principles.
The first step is to choose the size of the particle population, which is initialized with N particles. The position of the i-th particle is denoted as \({X}_{i},i=\text{1,2},\dots ,N\), the velocity of the i-th particle is denoted as \({V}_{i},i=\text{1,2},\dots ,N\), the optimal position passed by the i-th particle is \({pbest}_{i},i=\text{1,2},\dots ,N\), and the optimal position passed by the whole particle population is \(gbest\). The learning factors \({C}_{1}\) and \({C}_{2}\) are used to adjust the maximum step of particle motion during each iteration. \({C}_{1}\) is mainly used to adjust the maximum step of particle motion so that the particles can move toward the optimal position of an individual. \({C}_{2}\) moves the particles toward the optimal position of the whole group. The convergence speed is adjusted using inertia weights \(W\) [31]. After the initialization and parameters are determined, the reliability constraint, i.e., the upper limit of the preventive maintenance interval, needs to be defined in the algorithm. All the particles will be iteratively searched for the best in this interval, define this interval is:
$$\begin{array}{c}{Xlimit}_{min}\le {Xlimit}_{i}\le {Xlimit}_{max}\left(1\right)\end{array}$$
Where:
\({Xlimit}_{min}\) the minimum allowable position;
\({Xlimit}_{max}\) the maximum allowable position.
Similarly, \({V}_{i}\) is to be constrained by defining the velocity interval is:
$$\begin{array}{c}-\left({V}_{{i}_{min}}\right)\le {V}_{i}\le {V}_{{i}_{max}}\left(2\right)\end{array}$$
The maximum velocity of the particle is calculated by formula (3):
$$\begin{array}{c}{V}_{{i}_{max}}=\frac{Xlimitmax-Xlimitmin}{t}\left(3\right)\end{array}$$
Where:
\({V}_{i}\) velocity of the i-th particle;
\({V}_{{i}_{min}}\) the minimum allowable velocity;
\({V}_{{i}_{max}}\) the maximum allowable velocity;
\(t\) number of iterations.
\({V}_{i}\) and \({X}_{i}\) are calculated by formula (4) and formula (5):
$$\begin{array}{c}{V}_{i}\left(t+1\right)={V}_{i}\left(t\right)\times W+{C}_{1}\times {R}_{1}\times \left({pbest}_{i}\left(t\right)-{X}_{i}\left(t\right)\right)+{C}_{2}\times {R}_{2}\times \left(gbest\left(t\right)-{X}_{i}\left(t\right)\right)\left(4\right)\end{array}$$
$$\begin{array}{c}{X}_{i}\left(t+1\right)={X}_{i}\left(t\right)+{V}_{i}\left(t+1\right)\left(5\right)\end{array}$$
Where:
\({V}_{i}\left(t+1\right)\) updated velocity;
\({V}_{i}\left(t\right)\) velocity of the i-th particle;
\(W\) inertia weights;
\({C}_{1}\) cognitive constant;
\({R}_{1}\) random number between 0 and 1;
\({C}_{2}\) social constant;
\({R}_{2}\) random number between 0 and 1;
\({X}_{i}\left(t+1\right)\) updated position;
\({X}_{i}\left(t\right)\) position of the i-th particle.
\({R}_{1}\) and \({R}_{2}\) are independent of each other, the diversity of the population is determined by these two parameters; The inertia weight \(W\) affects the flight speed of the particle and the search range of the solution space. Typically, larger inertia weights have an enhanced effect on the global search capability of the algorithm, while smaller inertia weights have an enhanced effect on the local search capability of the algorithm[32]. The inertia weight \(W\) is calculated by formula (6):
$$\begin{array}{c}W={W}_{max}-\frac{\left({W}_{max}-{W}_{min}\right)\times t}{{t}_{max}}\left(6\right)\end{array}$$
Where:
\({W}_{min}\) initial inertial weight;
\({W}_{max}\) final inertial weight;
\({t}_{max}\) Maximum number of iterations.
In this paper, the minimum value of inertia weights is 0.4 and the maximum value is 0.9, and these values have good convergence accuracy. The values of \({C}_{1}\) and \({C}_{2}\) also help the particles converge faster as they move toward the optimal solution. The learning factor values are chosen by trial runs on the algorithm between 0 and 4, which is the typical range used for general particle swarm optimization problems. It can be observed that the convergence of the PSO is faster when the values of \({C}_{1}\) and \({C}_{2}\) are taken as 2. The steps for the implementing of the particle swarm optimization algorithm are shown in Fig. 5 using the following flowchart.
4.2 Machine tool maintenance cost optimization model
The probability of the machine completing the machining task at the specified time and conditions, is denoted as \(R\left(t\right)\). \(R\left(t\right)\), also known as the reliability function, is a probabilistic description of the time-dependent behaviour of the product[32]. Assuming that the operating time of the machine tool before failure is T and the specified operating time is t, the expression for the reliability of the machine tool at moment \(t\) is:
$$\begin{array}{c}R\left(t\right)=P\left(T>t\right) t>0\left(7\right)\end{array}$$
Reliability \(R\left(t\right)\) is a function of time, and the greater the reliability, the greater the reliability of the product to perform the specified function, and 0 ≤ \(R\left(t\right)\) ≤ 1. When the product is put into use, the reliability is 1, i.e. \(R\left(0\right)\) = 1; With the use of time reliability decreases with increasing use time. The possibility that the machine will not complete the machining task given the specified time and conditions. That can be described by the distribution function \(F\left(t\right)\) as follows:
$$\begin{array}{c}F\left(t\right)=P\left(T\le t\right) t>0\left(8\right)\end{array}$$
From the above equation, it can be seen that \(F\left(t\right)\) is the possibility that describes when the product fails at a time less than or equal to \(t\). In reliability theory, \(F\left(t\right)\) is called the failure distribution function. According to the definition of the failure distribution function and reliability, the relationship between \(F\left(t\right)\)and \(R\left(t\right)\) is derived as:
$$\begin{array}{c}F\left(t\right)=1-R\left(t\right) t>0\left(9\right)\end{array}$$
The time between failures of the machine is related to the level of reliability of the machine itself, and the reliability parameters of the machine are necessary to make maintenance interval decisions. By collecting the reliability data of each component during machine operation, parameter estimation, distribution fitting and hypothesis testing are performed to obtain the distribution of time between failures, then the reliability parameters of each component are obtained. There are many common equipment failure distribution functions, including exponential distribution, Weibull distribution, logarithmic distribution and normal distribution, etc. Equipment fault distribution of the machine tool selected by Weibull distribution. The probability of failure density function of the Weibull distribution[33] is calculated by formula (10):
$$\begin{array}{c}f\left(t\right)=\frac{\beta }{\theta }{\left(\frac{t-\gamma }{\theta }\right)}^{\beta -1}exp\left[{-\left(\frac{t-\gamma }{\theta }\right)}^{\beta }\right]\left(10\right)\end{array}$$
Where:
\(\beta\) shape factor;
\(\theta\) scale factor;
\(\gamma\) location factor.
Further, the reliability function R(t) is calculated by formula (11):
$$\begin{array}{c}R\left(t\right)= exp\left[{-\left(\frac{t-\gamma }{\theta }\right)}^{\beta }\right]\left(11\right)\end{array}$$
The cumulative failure probability density function F(t) is calculated by formula (12):
$$\begin{array}{c}F\left(t\right)=1-exp\left[{-\left(\frac{t-\gamma }{\theta }\right)}^{\beta }\right]\left(12\right)\end{array}$$
The maintenance costs of rotary ultrasonic vibration assisted EDM machines is calculated by formula (13):
$$\begin{array}{c}C\left(T\right)={C}_{p}+{C}_{c}\left(13\right)\end{array}$$
Where, \({C}_{p}\) is preventive maintenance(PM) costs, \({C}_{c}\) is the cost of repair after failure, and \(C\left(T\right)\) is the total maintenance cost. Assuming that the average cost of PM cost for the m-th component of the machine tool is \({C}_{pm}\), the average cost of repairing the m-th component after failure is \({C}_{cm}\), the reliability level of the accessory after each repair is the same, the probability of performing preventive maintenance during the preventive maintenance interval\(T\)is \(R\left(T\right)\) and the probability of performing post-failure restorative repair is\(F\left(T\right)\), so the total maintenance cost during the preventive maintenance interval is calculated by formula (14):
$$\begin{array}{c}C\left(T\right)=\sum _{m=1}^{M}{C}_{pm}R\left(T\right)+\sum _{m=1}^{M}{C}_{cm}F\left(T\right)\left(14\right)\end{array}$$
Substitute formula (11) and formula (12) into formula (14), C(T) is calculated by formula (15):
$$\begin{array}{c}C\left(T\right)=\sum _{m=1}^{M}{C}_{pm}{exp}\left[-{\left(\frac{T-{\gamma }_{m}}{{\theta }_{m}}\right)}^{{\beta }_{m}}\right]+\sum _{m=1}^{M}{C}_{cm}[1-{exp}\left[-{\left(\frac{T-{\gamma }_{m}}{{\theta }_{m}}\right)}^{{\beta }_{m}}\right]]\left(15\right)\end{array}$$
Also, in order to ensure that T is within the permissible pre-repair interval for each component, the upper limit of the preventive maintenance interval for the m-th component is calculated by formula (16):
$$\begin{array}{c}T<{\theta }_{m}{\left({ln}\left[\frac{1}{1-F\left(T\right)}\right]\right)}^{\frac{1}{{\beta }_{m}}}\left(16\right)\end{array}$$
The maximum preventive maintenance interval for each system of the machine tool is calculated based on the above equation, ensuring that the optimized results maintain the minimum level of reliability required. Considering the economic impact, choose the right maintenance interval to make the maintenance cost as cheap as possible, thus the total maintenance cost during the preventive maintenance interval can be used as the objective function and the following mathematical model can be constructed:
$$Min：C\left(T\right)=\frac{\sum _{m=1}^{M}{C}_{pm}{exp}\left[-{\left(\frac{T-{\gamma }_{m}}{{\theta }_{m}}\right)}^{{\beta }_{m}}\right]+\sum _{m=1}^{M}{C}_{cm}[1-{exp}\left[-{\left(\frac{T-{\gamma }_{m}}{{\theta }_{m}}\right)}^{{\beta }_{m}}\right]]}{T}$$
$$R\left(T\right)= exp\left[{-\left(\frac{T-{\gamma }_{m}}{\theta }\right)}^{\beta }\right]$$
$$Subject to:T<{\theta }_{m}{\left({ln}\left[\frac{1}{1-F\left(T\right)}\right]\right)}^{\frac{1}{{\beta }_{m}}}$$
$$m=\text{1,2},\dots ,M$$
Where:
\({C}_{pm}\) average PM cost of the m-th component;
\({C}_{cm}\) average cost of repairing the mth component after failure;
\({\beta }_{m}\) shape factor;
\({\theta }_{m}\) scale factor;
\({\gamma }_{m}\) location factor;
\(T\) preventive maintenance intervals;
\(F\left(T\right)\) maximum allowable failure probability function;
\(R\left(T\right)\) reliability function;
\(\left[R\right]\) lower limit of machine reliability.
When \(C\left(T\right)\) is the smallest, \(T\) is the optimal preventive maintenance interval, and \(T\) must be less than the upper limit of the preventive maintenance interval of each component. The model takes the total maintenance cost during the PM interval of the machine tools as the optimization objective and the maintenance interval as the decision variable. The maintenance interval optimization based on this model can reduce the maintenance cost while ensuring the reliability of machine tools.
4.3 Minimization of machine tool maintenance costs under reliability constraints
The expected service life of the machine tool is 87600 hours, the lower limit of machine reliability is 0.875. The Weibull distribution parameters, reliability, failure probability, and maintenance cost for each system are shown in Table 1.
Table 1
The machine tool systems main data
| |
Servo control
|
X-Y axis
|
Z-axis
|
Spindle rotation
|
Ultrasonic vibration
|
|
\({\theta }_{m}\)
|
1400
|
1475
|
1300
|
1400
|
1350
|
|
\({\beta }_{m}\)
|
3
|
4
|
3
|
2
|
3
|
|
\({\gamma }_{m}\)
|
0
|
0
|
0
|
0
|
0
|
|
\(R\left(T\right)\)
|
0.87
|
0.88
|
0.86
|
0.87
|
0.87
|
|
\(F\left(T\right)\)
|
0.13
|
0.12
|
0.14
|
0.13
|
0.13
|
|
\({C}_{pm}\)
|
1000
|
1100
|
1000
|
1150
|
1200
|
|
\({C}_{cm}\)
|
2000
|
2000
|
1500
|
2500
|
2500
|
Through the expression of the reliability function \(R\left(T\right)\) and the data set given in Table 1, the variation curve of the reliability with the period of the premaintenance interval is obtained by substituting the data into formula (11), and the graph is shown in Fig. 6. The graph shows that the reliability decreases as the premaintenance interval rises, and drops to a minimum of 0.6 after the interval reaches 1000 hours. According to the reliability constraint of 0.875 for the machine tool can be known from the image that the preventive maintenance interval cycle is at [0,534] hours. According to the algorithm proposed in Section 4.1, it is known that the range of the maximum optimality-seeking interval for the particles is [0,534], and the values are substituted into formula (1) to add reliability constraints to the algorithm.
According to the solution idea of this paper, the particle swarm algorithm mentioned in Section 4.1 was developed by applying MATLAB software. To reduce the maintenance cost, the number of particles is selected as 50. The maximum number of iterations is set to 600. The rest of the initialization parameters are shown in Table 2. The position and velocity of each particle in the population are updated according to formula (4)-(6), and \({pbest}_{i}\) and \(gbest\) are stored in the particle swarm algorithm memory and are not displayed at each iteration.
Table 2
PSO initialization parameters
|
no.
|
Parameters
|
Value
|
|
1
|
\({W}_{min}\)
|
0.4
|
|
2
|
\({W}_{max}\)
|
0.9
|
|
3
|
\({C}_{1}\)
|
2
|
|
4
|
\({C}_{2}\)
|
2
|
|
5
|
\({R}_{1}\)
|
Between 0 and 1
|
|
6
|
\({R}_{2}\)
|
Between 0 and 1
|
After that, the iterative optimization search is repeated until the optimization is completed after the specified iteration stopping criterion is reached, the image of the maintenance cost variation with the maintenance cycle is generated. The extreme value point with the lowest maintenance cost under the reliability constraint is represented in the figure. The upper limit of the preventive maintenance interval for each system is shown in Table 3, which ensures that the optimized preventive maintenance interval is less than the minimum value for each system. If the final optimized period is greater than the minimum value of each system interval period, the optimization result is invalid.
Table 3
Preventive maintenance interval for each system
| |
Servo control
|
X-Y axis
|
Z-axis
|
Spindle rotation
|
Ultrasonic vibration
|
|
PM interval(h)
|
624
|
652
|
565
|
534
|
441
|
The setting of iteration termination criteria (iteration number \(t\)) also affects the optimization results. The optimal set of solutions obtained by solving MATLAB with 300–600 iterations respectively is shown in Fig. 7.
The specific optimization values corresponding to different numbers of iterations in the above figure are shown in Table 4. From Fig. 7 and Table 4, it can be seen that a different number of iterations affects the optimization results. As the numbers of iterations continues to rise, the maintenance cost of the machine tool decreases, and the corresponding preventive maintenance interval increases. When the numbers of iterations are 700, the iteration stopping criterion has been exceeded, the maintenance cost is ¥20761 and the preventive maintenance interval is 486 hours. Although the maintenance cost is reduced compared with the cost of 600 iterations, it can be seen that the preventive maintenance interval is higher than the minimum value of the ultrasonic vibration system maintenance interval and the optimization result is invalid.
Table 4
Cost and interval corresponding to different number of iterations
|
Number of iterations
|
Total Maintenance Costs(¥)
|
PM interval(h)
|
|
300
|
47218
|
212
|
|
400
|
35744
|
279
|
|
500
|
28731
|
349
|
|
600
|
24088
|
417
|
|
700
|
20761
|
486
|
After 600 iterations of particle updates, the program arrived at the final optimized resulting optimal preventive maintenance interval for the machine tool system, i.e., \(gbest\), of 417 hours, which is less than the maintenance cycle of each system, and the corresponding minimum preventive maintenance cost of ¥24088, as shown in Fig. 7. This means it will cost ¥24088 to perform preventive maintenance tasks every 417 hours during the planning cycle.
In this section, the particle swarm algorithm mentioned in Section 4.1 is adopted to solve model, and a program is developed using MATLAB to analyze the relationship between the PM interval and the corresponding maintenance cost for rotating ultrasonic vibration-assisted EDM machines with different numbers of iterations. The method does not involve complex mathematical steps. It takes into account reliability constraints by incorporating them into the algorithm, which not only helps to select the best maintenance schedule, but also helps to simplify the maintenance cost optimization steps and reduce the downtime loss due to machine downtime caused by failures.
4.1 Steps for implementing the PSO
PSO is a study of the process of searching for food in a flock of birds in nature. The algorithm simulates the behaviour of a flock of birds searching for food and the cooperation between birds to achieve the optimal purpose of the flock. The PSO is a very simple algorithm that performs effective optimization for various functions. As an optimization algorithm based on the theory of population intelligence, it operates on the principle that the population intelligence generated by the cooperative and competitive relationship between the particles in the population guides the optimization search[30]. In PSO, each particle represents a bird of the natural bird community, and the range of the bird in search of food is the search space of the particle. All particles have a fitness value, which is determined by the objective function. Velocity and position are the two basic parameters of a particle. The velocity represents the speed of the particle's flight and the position represents the direction of the flight. Each particle searches for its optimal solution in the search space and records it as \(pbest\) to share with the other particles in the population. After that, the optimal solutions found by each particle are compared and the best result is recorded as the current population optimal solution \(gbest\). The velocities and positions of all particles in the population are adjusted by \(pbest\) and \(gbest\). In this paper, algorithms are developed using MATLAB to minimize the preventive maintenance cost of machine tools based on the above principles.
The first step is to choose the size of the particle population, which is initialized with N particles. The position of the i-th particle is denoted as \({X}_{i},i=\text{1,2},\dots ,N\), the velocity of the i-th particle is denoted as \({V}_{i},i=\text{1,2},\dots ,N\), the optimal position passed by the i-th particle is \({pbest}_{i},i=\text{1,2},\dots ,N\), and the optimal position passed by the whole particle population is \(gbest\). The learning factors \({C}_{1}\) and \({C}_{2}\) are used to adjust the maximum step of particle motion during each iteration. \({C}_{1}\) is mainly used to adjust the maximum step of particle motion so that the particles can move toward the optimal position of an individual. \({C}_{2}\) moves the particles toward the optimal position of the whole group. The convergence speed is adjusted using inertia weights \(W\) [31]. After the initialization and parameters are determined, the reliability constraint, i.e., the upper limit of the preventive maintenance interval, needs to be defined in the algorithm. All the particles will be iteratively searched for the best in this interval, define this interval is:
$$\begin{array}{c}{Xlimit}_{min}\le {Xlimit}_{i}\le {Xlimit}_{max}\left(1\right)\end{array}$$
Where:
\({Xlimit}_{min}\) the minimum allowable position;
\({Xlimit}_{max}\) the maximum allowable position.
Similarly, \({V}_{i}\) is to be constrained by defining the velocity interval is:
$$\begin{array}{c}-\left({V}_{{i}_{min}}\right)\le {V}_{i}\le {V}_{{i}_{max}}\left(2\right)\end{array}$$
The maximum velocity of the particle is calculated by formula (3):
$$\begin{array}{c}{V}_{{i}_{max}}=\frac{Xlimitmax-Xlimitmin}{t}\left(3\right)\end{array}$$
Where:
\({V}_{i}\) velocity of the i-th particle;
\({V}_{{i}_{min}}\) the minimum allowable velocity;
\({V}_{{i}_{max}}\) the maximum allowable velocity;
\(t\) number of iterations.
\({V}_{i}\) and \({X}_{i}\) are calculated by formula (4) and formula (5):
$$\begin{array}{c}{V}_{i}\left(t+1\right)={V}_{i}\left(t\right)\times W+{C}_{1}\times {R}_{1}\times \left({pbest}_{i}\left(t\right)-{X}_{i}\left(t\right)\right)+{C}_{2}\times {R}_{2}\times \left(gbest\left(t\right)-{X}_{i}\left(t\right)\right)\left(4\right)\end{array}$$
$$\begin{array}{c}{X}_{i}\left(t+1\right)={X}_{i}\left(t\right)+{V}_{i}\left(t+1\right)\left(5\right)\end{array}$$
Where:
\({V}_{i}\left(t+1\right)\) updated velocity;
\({V}_{i}\left(t\right)\) velocity of the i-th particle;
\(W\) inertia weights;
\({C}_{1}\) cognitive constant;
\({R}_{1}\) random number between 0 and 1;
\({C}_{2}\) social constant;
\({R}_{2}\) random number between 0 and 1;
\({X}_{i}\left(t+1\right)\) updated position;
\({X}_{i}\left(t\right)\) position of the i-th particle.
\({R}_{1}\) and \({R}_{2}\) are independent of each other, the diversity of the population is determined by these two parameters; The inertia weight \(W\) affects the flight speed of the particle and the search range of the solution space. Typically, larger inertia weights have an enhanced effect on the global search capability of the algorithm, while smaller inertia weights have an enhanced effect on the local search capability of the algorithm[32]. The inertia weight \(W\) is calculated by formula (6):
$$\begin{array}{c}W={W}_{max}-\frac{\left({W}_{max}-{W}_{min}\right)\times t}{{t}_{max}}\left(6\right)\end{array}$$
Where:
\({W}_{min}\) initial inertial weight;
\({W}_{max}\) final inertial weight;
\({t}_{max}\) Maximum number of iterations.
In this paper, the minimum value of inertia weights is 0.4 and the maximum value is 0.9, and these values have good convergence accuracy. The values of \({C}_{1}\) and \({C}_{2}\) also help the particles converge faster as they move toward the optimal solution. The learning factor values are chosen by trial runs on the algorithm between 0 and 4, which is the typical range used for general particle swarm optimization problems. It can be observed that the convergence of the PSO is faster when the values of \({C}_{1}\) and \({C}_{2}\) are taken as 2. The steps for the implementing of the particle swarm optimization algorithm are shown in Fig. 5 using the following flowchart.
4.2 Machine tool maintenance cost optimization model
The probability of the machine completing the machining task at the specified time and conditions, is denoted as \(R\left(t\right)\). \(R\left(t\right)\), also known as the reliability function, is a probabilistic description of the time-dependent behaviour of the product[32]. Assuming that the operating time of the machine tool before failure is T and the specified operating time is t, the expression for the reliability of the machine tool at moment \(t\) is:
$$\begin{array}{c}R\left(t\right)=P\left(T>t\right) t>0\left(7\right)\end{array}$$
Reliability \(R\left(t\right)\) is a function of time, and the greater the reliability, the greater the reliability of the product to perform the specified function, and 0 ≤ \(R\left(t\right)\) ≤ 1. When the product is put into use, the reliability is 1, i.e. \(R\left(0\right)\) = 1; With the use of time reliability decreases with increasing use time. The possibility that the machine will not complete the machining task given the specified time and conditions. That can be described by the distribution function \(F\left(t\right)\) as follows:
$$\begin{array}{c}F\left(t\right)=P\left(T\le t\right) t>0\left(8\right)\end{array}$$
From the above equation, it can be seen that \(F\left(t\right)\) is the possibility that describes when the product fails at a time less than or equal to \(t\). In reliability theory, \(F\left(t\right)\) is called the failure distribution function. According to the definition of the failure distribution function and reliability, the relationship between \(F\left(t\right)\)and \(R\left(t\right)\) is derived as:
$$\begin{array}{c}F\left(t\right)=1-R\left(t\right) t>0\left(9\right)\end{array}$$
The time between failures of the machine is related to the level of reliability of the machine itself, and the reliability parameters of the machine are necessary to make maintenance interval decisions. By collecting the reliability data of each component during machine operation, parameter estimation, distribution fitting and hypothesis testing are performed to obtain the distribution of time between failures, then the reliability parameters of each component are obtained. There are many common equipment failure distribution functions, including exponential distribution, Weibull distribution, logarithmic distribution and normal distribution, etc. Equipment fault distribution of the machine tool selected by Weibull distribution. The probability of failure density function of the Weibull distribution[33] is calculated by formula (10):
$$\begin{array}{c}f\left(t\right)=\frac{\beta }{\theta }{\left(\frac{t-\gamma }{\theta }\right)}^{\beta -1}exp\left[{-\left(\frac{t-\gamma }{\theta }\right)}^{\beta }\right]\left(10\right)\end{array}$$
Where:
\(\beta\) shape factor;
\(\theta\) scale factor;
\(\gamma\) location factor.
Further, the reliability function R(t) is calculated by formula (11):
$$\begin{array}{c}R\left(t\right)= exp\left[{-\left(\frac{t-\gamma }{\theta }\right)}^{\beta }\right]\left(11\right)\end{array}$$
The cumulative failure probability density function F(t) is calculated by formula (12):
$$\begin{array}{c}F\left(t\right)=1-exp\left[{-\left(\frac{t-\gamma }{\theta }\right)}^{\beta }\right]\left(12\right)\end{array}$$
The maintenance costs of rotary ultrasonic vibration assisted EDM machines is calculated by formula (13):
$$\begin{array}{c}C\left(T\right)={C}_{p}+{C}_{c}\left(13\right)\end{array}$$
Where, \({C}_{p}\) is preventive maintenance(PM) costs, \({C}_{c}\) is the cost of repair after failure, and \(C\left(T\right)\) is the total maintenance cost. Assuming that the average cost of PM cost for the m-th component of the machine tool is \({C}_{pm}\), the average cost of repairing the m-th component after failure is \({C}_{cm}\), the reliability level of the accessory after each repair is the same, the probability of performing preventive maintenance during the preventive maintenance interval\(T\)is \(R\left(T\right)\) and the probability of performing post-failure restorative repair is\(F\left(T\right)\), so the total maintenance cost during the preventive maintenance interval is calculated by formula (14):
$$\begin{array}{c}C\left(T\right)=\sum _{m=1}^{M}{C}_{pm}R\left(T\right)+\sum _{m=1}^{M}{C}_{cm}F\left(T\right)\left(14\right)\end{array}$$
Substitute formula (11) and formula (12) into formula (14), C(T) is calculated by formula (15):
$$\begin{array}{c}C\left(T\right)=\sum _{m=1}^{M}{C}_{pm}{exp}\left[-{\left(\frac{T-{\gamma }_{m}}{{\theta }_{m}}\right)}^{{\beta }_{m}}\right]+\sum _{m=1}^{M}{C}_{cm}[1-{exp}\left[-{\left(\frac{T-{\gamma }_{m}}{{\theta }_{m}}\right)}^{{\beta }_{m}}\right]]\left(15\right)\end{array}$$
Also, in order to ensure that T is within the permissible pre-repair interval for each component, the upper limit of the preventive maintenance interval for the m-th component is calculated by formula (16):
$$\begin{array}{c}T<{\theta }_{m}{\left({ln}\left[\frac{1}{1-F\left(T\right)}\right]\right)}^{\frac{1}{{\beta }_{m}}}\left(16\right)\end{array}$$
The maximum preventive maintenance interval for each system of the machine tool is calculated based on the above equation, ensuring that the optimized results maintain the minimum level of reliability required. Considering the economic impact, choose the right maintenance interval to make the maintenance cost as cheap as possible, thus the total maintenance cost during the preventive maintenance interval can be used as the objective function and the following mathematical model can be constructed:
$$Min：C\left(T\right)=\frac{\sum _{m=1}^{M}{C}_{pm}{exp}\left[-{\left(\frac{T-{\gamma }_{m}}{{\theta }_{m}}\right)}^{{\beta }_{m}}\right]+\sum _{m=1}^{M}{C}_{cm}[1-{exp}\left[-{\left(\frac{T-{\gamma }_{m}}{{\theta }_{m}}\right)}^{{\beta }_{m}}\right]]}{T}$$
$$R\left(T\right)= exp\left[{-\left(\frac{T-{\gamma }_{m}}{\theta }\right)}^{\beta }\right]$$
$$Subject to:T<{\theta }_{m}{\left({ln}\left[\frac{1}{1-F\left(T\right)}\right]\right)}^{\frac{1}{{\beta }_{m}}}$$
$$m=\text{1,2},\dots ,M$$
Where:
\({C}_{pm}\) average PM cost of the m-th component;
\({C}_{cm}\) average cost of repairing the mth component after failure;
\({\beta }_{m}\) shape factor;
\({\theta }_{m}\) scale factor;
\({\gamma }_{m}\) location factor;
\(T\) preventive maintenance intervals;
\(F\left(T\right)\) maximum allowable failure probability function;
\(R\left(T\right)\) reliability function;
\(\left[R\right]\) lower limit of machine reliability.
When \(C\left(T\right)\) is the smallest, \(T\) is the optimal preventive maintenance interval, and \(T\) must be less than the upper limit of the preventive maintenance interval of each component. The model takes the total maintenance cost during the PM interval of the machine tools as the optimization objective and the maintenance interval as the decision variable. The maintenance interval optimization based on this model can reduce the maintenance cost while ensuring the reliability of machine tools.
4.3 Minimization of machine tool maintenance costs under reliability constraints
The expected service life of the machine tool is 87600 hours, the lower limit of machine reliability is 0.875. The Weibull distribution parameters, reliability, failure probability, and maintenance cost for each system are shown in Table 1.
Table 1
The machine tool systems main data
| |
Servo control
|
X-Y axis
|
Z-axis
|
Spindle rotation
|
Ultrasonic vibration
|
|
\({\theta }_{m}\)
|
1400
|
1475
|
1300
|
1400
|
1350
|
|
\({\beta }_{m}\)
|
3
|
4
|
3
|
2
|
3
|
|
\({\gamma }_{m}\)
|
0
|
0
|
0
|
0
|
0
|
|
\(R\left(T\right)\)
|
0.87
|
0.88
|
0.86
|
0.87
|
0.87
|
|
\(F\left(T\right)\)
|
0.13
|
0.12
|
0.14
|
0.13
|
0.13
|
|
\({C}_{pm}\)
|
1000
|
1100
|
1000
|
1150
|
1200
|
|
\({C}_{cm}\)
|
2000
|
2000
|
1500
|
2500
|
2500
|
Through the expression of the reliability function \(R\left(T\right)\) and the data set given in Table 1, the variation curve of the reliability with the period of the premaintenance interval is obtained by substituting the data into formula (11), and the graph is shown in Fig. 6. The graph shows that the reliability decreases as the premaintenance interval rises, and drops to a minimum of 0.6 after the interval reaches 1000 hours. According to the reliability constraint of 0.875 for the machine tool can be known from the image that the preventive maintenance interval cycle is at [0,534] hours. According to the algorithm proposed in Section 4.1, it is known that the range of the maximum optimality-seeking interval for the particles is [0,534], and the values are substituted into formula (1) to add reliability constraints to the algorithm.
According to the solution idea of this paper, the particle swarm algorithm mentioned in Section 4.1 was developed by applying MATLAB software. To reduce the maintenance cost, the number of particles is selected as 50. The maximum number of iterations is set to 600. The rest of the initialization parameters are shown in Table 2. The position and velocity of each particle in the population are updated according to formula (4)-(6), and \({pbest}_{i}\) and \(gbest\) are stored in the particle swarm algorithm memory and are not displayed at each iteration.
Table 2
PSO initialization parameters
|
no.
|
Parameters
|
Value
|
|
1
|
\({W}_{min}\)
|
0.4
|
|
2
|
\({W}_{max}\)
|
0.9
|
|
3
|
\({C}_{1}\)
|
2
|
|
4
|
\({C}_{2}\)
|
2
|
|
5
|
\({R}_{1}\)
|
Between 0 and 1
|
|
6
|
\({R}_{2}\)
|
Between 0 and 1
|
After that, the iterative optimization search is repeated until the optimization is completed after the specified iteration stopping criterion is reached, the image of the maintenance cost variation with the maintenance cycle is generated. The extreme value point with the lowest maintenance cost under the reliability constraint is represented in the figure. The upper limit of the preventive maintenance interval for each system is shown in Table 3, which ensures that the optimized preventive maintenance interval is less than the minimum value for each system. If the final optimized period is greater than the minimum value of each system interval period, the optimization result is invalid.
Table 3
Preventive maintenance interval for each system
| |
Servo control
|
X-Y axis
|
Z-axis
|
Spindle rotation
|
Ultrasonic vibration
|
|
PM interval(h)
|
624
|
652
|
565
|
534
|
441
|
The setting of iteration termination criteria (iteration number \(t\)) also affects the optimization results. The optimal set of solutions obtained by solving MATLAB with 300–600 iterations respectively is shown in Fig. 7.
The specific optimization values corresponding to different numbers of iterations in the above figure are shown in Table 4. From Fig. 7 and Table 4, it can be seen that a different number of iterations affects the optimization results. As the numbers of iterations continues to rise, the maintenance cost of the machine tool decreases, and the corresponding preventive maintenance interval increases. When the numbers of iterations are 700, the iteration stopping criterion has been exceeded, the maintenance cost is ¥20761 and the preventive maintenance interval is 486 hours. Although the maintenance cost is reduced compared with the cost of 600 iterations, it can be seen that the preventive maintenance interval is higher than the minimum value of the ultrasonic vibration system maintenance interval and the optimization result is invalid.
Table 4
Cost and interval corresponding to different number of iterations
|
Number of iterations
|
Total Maintenance Costs(¥)
|
PM interval(h)
|
|
300
|
47218
|
212
|
|
400
|
35744
|
279
|
|
500
|
28731
|
349
|
|
600
|
24088
|
417
|
|
700
|
20761
|
486
|
After 600 iterations of particle updates, the program arrived at the final optimized resulting optimal preventive maintenance interval for the machine tool system, i.e., \(gbest\), of 417 hours, which is less than the maintenance cycle of each system, and the corresponding minimum preventive maintenance cost of ¥24088, as shown in Fig. 7. This means it will cost ¥24088 to perform preventive maintenance tasks every 417 hours during the planning cycle.
In this section, the particle swarm algorithm mentioned in Section 4.1 is adopted to solve model, and a program is developed using MATLAB to analyze the relationship between the PM interval and the corresponding maintenance cost for rotating ultrasonic vibration-assisted EDM machines with different numbers of iterations. The method does not involve complex mathematical steps. It takes into account reliability constraints by incorporating them into the algorithm, which not only helps to select the best maintenance schedule, but also helps to simplify the maintenance cost optimization steps and reduce the downtime loss due to machine downtime caused by failures.