1. INTRODUCTION
Many production units suffer from poor coordination and information exchange between activities, which tend to result in wasteful use of resources (cost, time, materials, etc.) and the emergence of an island culture in their organization. Naturally, all parts of a production unit are interdependent, which is why management decisions need to be made with an integrated approach taking into account all critical factors concerning the decision. This need for integration has long been a subject of interest to researchers and industry managers. According to researchers like Noor al, many other researchers, including Mirzanti et al. (2019) have highlighted the importance of integrating supply, maintenance, and repair activities. Many other researchers, like Boslah et al., have working on integrating production, quality control, and maintenance operations (Berrichi et al., 2009). Therefore, it can be seen that the integrated approach to the management of different domains of production and the integration of information related to different production activities have become popular topics over the past few decades. Having mutual interactions, these domains cannot be treated as islands and it is more reasonable to plan a mechanism for considering all of them together as much as possible. Indeed, many efforts have been made to integrate the plan ning of production, maintenance, and supply procurement activities (Amouzad Mahdiraji and Yousefi Talouki, 2021;Nasiri and Ahmadi Daryakenari, 2021;Sharma and Sharma, 2021).
In this research, we present a new mathematical model for production and maintenance planning and also an efficient metaheuristic algorithm for solving this model. In the remainder of this paper, section 2 reviews the latest studies relevant to the subject, section 3 presents the proposed mathematical model, section 4 explains the fuzzy uncertainty method used in the model, section 5 presents the numerical results obtained with the model, and section 6 concludes the paper.
In a study carried out by Liao (2013), they tried to integrate production planning and maintenance strategies for a deteriorating singlemachine system. In this paper, periodic imperfect preventive maintenance and periodic perfect maintenance were integrated to minimize the expected cost of production and maintenance activities. The objective in this study was to determine the duration of preventive and perfect maintenance cycles such that the expected cost of preventive maintenance is minimized. The results of this study showed that the optimal production plan and the optimal preventive maintenance are interdependent (Aliakbari, 2017).
In a study by Yalaoui et al. (2014), an extended linear programming model was developed for determining the optimum production plan with minimum total cost. The production planning and maintenance model of this study was biobjective, multiproduct, and multiperiod, and had multiple deteriorating production lines. It was also assumed that maintenance operations restore the machines to their As Good As New (AGAN) condition.
Nourelfath et al. (2016) integrated the production planning, maintenance, and quality control operations of an imperfect process in a multiperiod multiproduct system. In this study, the main objective was to minimize all costs while meeting all demands. Ultimately, it was reported raising the level of preventive maintenance lowers the quality control costs. Also, preventive maintenance will not be feasible if its costs are not offset by the reduction in the quality control costs. The results obtained by solving the model also showed a strong relationship between production, preventive maintenance, and quality control operations.
In a study by Ettaye et al. (2017), they integrated maintenance and production planning with a periodic replacement scheme with minimal repairs (preventive maintenance) to prevent unforeseen breakdowns. Because of the complexity of this problem, it takes a lot of time to solve it with exact solution methods, thus these researchers recommended using approximation methods to reach acceptable solutions in reasonably short times.
Ekin (2018) studied the impact of uncertain endogenous factors in the integration of production planning and maintenance programs. For this purpose, a decision model was developed in which uncertainty in the availability of machines serves as a proxy for the aforementioned factors. This model was optimized with a probabilistic simulation method. In the end, it was concluded that there should always be a balance between operating the machines and maintaining them.
Alimian et al. (2019) presented a robust mathematical model for preventive maintenance and production planning. This model takes demand uncertainty into account and utilizes a robust optimization approach to deal with this uncertainty. In this study, maintenance was planned based on a common cause breakdown approach. The results showed that the maintenance strategy can have a significant impact on production costs as well as product return rates.
In a study by Goli et al. (2019), annual production planning was integrated with workforce planning. These researchers formulated a biobjective mathematical model with the objectives of reducing total costs and maximizing customer satisfaction. They also presented two metaheuristic algorithms, namely multiobjective genetic algorithm and multiobjective invasive weed optimization algorithm for solving the model. The results showed the high efficiency of these solution methods.
Glawar et al. (2018) proposed an integrated model of preventive maintenance and production planning. In this model, flexibility and production quality are considered in the planning process. The objective of this model is to reduce all production and maintenance costs. The results of this study showed that combining these two areas of decisionmaking can be very effective in reducing and controlling costs.
In another study by Alimian et al. (2020), they discussed the impact of cell formation on production and maintenance planning and proposed an integrated production and maintenance planning approach, which makes use of dynamic cells and group scheduling. In their mathematical model, these researchers formulated the breakdown of the machines as a rate. After modeling, optimization was performed with GAMS software. The results showed that designing a suitable structure for production cells can greatly affect the breakdown of machines and their maintenance.
Bensmain et al. (2019) introduced a new model for preventive maintenance planning with the impact of maintenance equipment and poor production performance taken into account. In this model, product life costs are fully considered as an objective function. An improved genetic algorithm was used to optimize this mathematical model, the results of which showed the high efficiency of the solution method.
Wang et al. (2020) tried to integrate a conditionbased preventive maintenance method with production planning. In this study, breakdowns and customer demand were considered to be stochastic. A combination of optimization simulation and genetic algorithms was used to solve this problem. The results showed the superiority of the proposed method over similar methods.
Samimi and Sydow (2020) studied human resource management and maintenance in projectbased organizations. In this study, first, previous works in the field of human resource management were reviewed. Then the effect of this management on maintenance in projectbased organizations was studied. The results showed that human resource problems can affect the performance of the organization in various aspects.
After reviewing the previous studies in this field, it was observed that the most important gap in the literature is the integration of production planning, maintenance, and labor scheduling. In reality, these three areas are interrelated and can have a profound impact on the overall performance of a production unit. Thus, the main contribution of this study to the literature is the development of a novel mathematical model for the integration of production planning, maintenance, and labor scheduling. The second innovation of this study is the use of a metaheuristic algorithm called the ant colony algorithm, which has never been used in previous studies in this field.
2. METHODOLOGY
In this study, the labor scheduling aspect of the problem is defined as follows. Consider a factory with S machines (or departments). The goal is to determine a schedule for people to operate these machines during l days of the planning horizon. Each day consists of J shifts with a specified number of hours. The number of workers is K, which is known.
Each machine has a specific breakdown rate and needs to receive maintenance on each shift. The machines will have different production rates before and after the maintenance depending on their breakdown rate and the extent of preventive maintenance performed. Naturally, carrying out preventive maintenance will increase the production rate of the machines.
For a person to be assigned to a machine, he must have a certain level of experience in operating that machine. Also, each shift must have at least a certain number of workers. In addition, labor scheduling is subject to several conditions, which are described below:

•Work shifts are predefined.

•If a person works a night shift, he is not allowed to work the following morning shift

•Each person can take a leave for a maximum of 2 consecutive days. (nonconsecutive leave is allowed)

•If a person works the morning and afternoon shifts, he must not work the night shift.

•Performing tasks repeatedly increases the experience of the person, improving his work rate.

•Over time, fatigue affects a person’s performance, reducing his work rate.
In the proposed mathematical model, the output of each machine depends on three factors: experience of the operator, breakdown rate of the machine, and the machine utilization rate. Since these factors change with conditions, all three of them (breakdown rate, utilization rate, and operator experience) are considered to be fuzzy.
The goal is to provide a work schedule that, while meeting the above conditions, yields the highest production output with a minimum number of workers. Before presenting the proposed mathematical model, we first describe the indices, sets, parameters, and variables used in the model.
Indices and sets

I : Set of days in the planning horizon

i : Index of days (i ∈ I)

J : Set of work shifts

j, j' : Index of work shifts (j, j'∈J)

K : Set of workers

k,k' : Index of workers (k, l ∈ K)

S : Set of machines (departments)

s : Index of machines (s ∈S)
Parameters

M_{ijs} : Minimum number of workers required in shift j from day i to work on machine s

LI_{ks} : Minimum experience of worker k for working on machine s (in terms of output), which is considered to be a fuzzy number

KI_{s} : Maximum efficiency of machine s under ideal conditions (in terms of output)

LE_{s} : Utilization rate of machine s, which is considered to be a fuzzy number.

LF_{s} : Breakdown rate of machine s, which is considered to be a fuzzy number

FG_{js} : Production rate of machine s during shift j before maintenance

Bt_{j} : Duration of maintenance in shift j

RG_{js} : Production rate of machine s during shift j after maintenance

λ : Coefficient of increased breakdown during work

μ : Coefficient of decreased breakdown during rest

δ_{j} : Duration of shift j

U : Maximum working days of a worker in the planning horizon

L : Minimum working days of a worker in the planning horizon

n : Number of days in the planning horizon

n' : Maximum number of night shifts for each worker

t : Maximum number of consecutive working days for each worker
Variables

X_{ijks} : A binary variable showing whether worker k works on machine s on shift j of day i.

Q_{ijks} : Output of worker k working on machine s on shift j of day i
Equation 1 is the objective function of the mathematical model, which minimizes the labor employed divided by the production efficiency. This objective function tries to draw the best possible output from the machines for the given production rate, breakdown rate, and worker experience. Equation 2 computes the production output in terms of the exponential function of shift time and breakdown rate. Equation 3 gives the production output after maintenance, which decreases exponentially and depends on the duration of maintenance. Equation 4 computes the production output in the next shift. This output is affected by the breakdown rate in the previous shift, which increases exponentially. Equation 5 calculates the production output of each worker on each machine in terms of breakdown rate, utilization rate, and worker experience. Equations 2 to 5 have been derived from previous studies on labor scheduling based on fatigue and learning rates.
Equation 6 ensures that each machine on each shift has the required number of workers. Equation 7 states that in each shift, each worker is either working or on leave. Equation 8 states that each worker must work between the minimum and maximum number of days allowed. Equation 9 states that workers who work a night shift (j = 3) are not allowed to work the next morning shift (j = 1). Equation 10 states that a worker can be on leave for a maximum of two consecutive days. Equation 11 states that if a worker works a morning shift and the immediately following afternoon shift, he cannot work the night shift. Equation 12 states that each worker can work a maximum of n’ night shifts.
3. RESULTS AND DISCUSSION
As explained in the model description, the model has three core parameters: utilization rate, breakdown rate, and minimum worker experience. However, there is no specific rule for initializing these parameters and this initialization will be done based on the opinion of experts on human resources. Since people tend to have diverging opinions on such matters, it is necessary to consider the uncertainty involved in these parameters. In this study, the uncertainty in utilization rate, breakdown rate, and minimum worker experience is modeled using the triangular fuzzy uncertainty method.
Suppose $\tilde{a}=({a}^{l},\hspace{0.17em}a,\hspace{0.17em}{a}^{u})\hspace{0.17em}\text{and}\hspace{0.17em}\tilde{b}=({b}^{l},\hspace{0.17em}b,\hspace{0.17em}{b}^{u})$ are two triangular fuzzy numbers. It is said that $\tilde{a}\le \tilde{b}$ if and only if ${a}^{l}\le {b}^{l}\hspace{0.17em}\text{and}\hspace{0.17em}a\le b$.
An rcut of the fuzzy number $\tilde{a}$ is ${a}_{r}=\{x\in \Omega {\mu}_{\tilde{a}}(x)\ge r\}$. Since ${\mu}_{\tilde{a}}$ is semicontinuous, rcut has a bounded and finite value in the form of ${a}_{r}=[{f}_{a}^{1}(r),{g}_{a}^{1}\left(r\right)]$.
The expected value of the fuzzy value $\tilde{a}$, denoted by $EV(\tilde{a})$, is equal to the midpoint of the expected distance, which is calculated as follows:
For the fuzzy numbers $\tilde{a}$ and $\tilde{b}$, the result of arithmetic operation $\tilde{a}\times \tilde{b}$, according to the Lotfizadeh expansion, is a fuzzy set with the following membership function:
Dubois and Prade obtained the following equation for an rcut:
where $\tilde{a}$ and $\tilde{b}$ are fuzzy numbers and λ and γ are real nonnegative numbers.
From equations (13), (14), and (15), it easily follows that:
In this study, the three model parameters that have a degree of uncertainty are considered nondeterministic and defined in the form of triangular fuzzy numbers. Hence, the parameters utilization rate, breakdown rate, and minimum worker experience are formulated as follows.
To modify the mathematical model for fuzzy conditions, we use the fuzzy cut described in Section 41. As explained in that section, a constraint with fuzzy parameters based on Equation 16 is used for this purpose. In the mathematical model, Equation 5 contains fuzzy parameters. Therefore, this equation is rewritten into Equation 20. Other model constraints remain unchanged. It should be noted that in Equation 20, ɑ is a control parameter between 0 and 1, which is usually set to 0.5.
The purpose of validation is to determine whether the model produces correct and reliable outputs. In this study, this evaluation was done by applying the model on a real case in Marun Petrochemical Company with the following specifications.
Ten workers are working in an active production unit. The goal is to schedule the work for a 7 day period with 3 shifts each day. The minimum and maximum number of shifts per worker is 3 and 15 respectively. In each shift, at least 2 people must work in each department. Each shift is 8 hours long. The maximum efficiency of the machine in each department is 80%. The maintenance time in each shift is 30 minutes. The coefficients of increased breakdown and decreased breakdown are 20% and 10% respectively. The parameters worker experience, utilization rate, and breakdown rate for different machines have been given random values Table 2. For each parameter, this random value has been generated using a continuous uniform distribution with the lower bound (L) and upper bound (U) given in Table 1. Also, given the fuzzy nature of the three parameters, three values have been specified for each parameter. These three values are the lower bound, mode, and upper bound of the fuzzy number, respectively.
After entering this information in GAMS and running the optimization model, the outputs shown in Figure 1 were obtained.
As can be seen, GAMS displayed the message “proven optimal solution”, which indicates that the obtained solution is the global optimum. Therefore, the mathematical model and its constraints have created a logical feasibility space, which shows that the optimal solution of the problem can be obtained.
To better illustrate the output, the optimal solution is presented in Tables 2, 3, and 4.
To present the results in these tables, workers were numbered from 1 to 10. The values given in these tables are the number of the worker who must work on a certain machine on a certain shift of a certain day. Table 2, for example, shows that in the first shift of Day 1, workers No. 3, 7, and 10 must work on machine 1 (in department 1).
The numbers given in each cell of Tables 2 to 4 represent the workers who are scheduled to work in the relevant department on the specified shift of the specified day. The presented results show that the output of the mathematical model satisfies all defined constraints, which demonstrates the validity of the formulated model.
This section presents the numerical results related to the evaluation of the performance of the ant colony optimization algorithm (ACO) in solving the formulated model. This performance evaluation is performed in terms of two criteria: solution quality and solution speed. Solution quality was measured in terms of the difference between the objective function value produced by ACO and the one obtained from GAMS. To examine solution speed, we compared the time it takes for these methods to produce the optimal solution. To evaluate the performance of ACO algorithms in solving the problem, this algorithm was coded in MATLAB and 10 problem instances of different sizes were produced. The specifications of these problems are given in Table 5. Other parameters were generated based on the information given in Section 51. More information about these problems is provided in Table 6.
Next, the results of the exact solution of the generated problem instances with GAMS were compared with the results of ACO.
Since it may take a very long time for GAMS to solve largescale problems, we set a time limit of 3600 seconds or 1 hour for reaching a solution. In cases where GAMS could not produce the optimum solution after 1 hour, the solution process was terminated and one of the produced feasible solutions (not necessarily the optimal one) was chosen as the output. Table 6 compares the results obtained from GAMS with those produced by ACO.
As can be seen, while GAMS could not find the optimal solution of the last three problems in less than 1 hour, ACO managed to produce better solutions for these problems in a much shorter time.
The solution speed curves of both methods are drawn in Figure 2.
As the above diagram shows, for the exact solution with GAMS, solution time has increased exponentially with the size of the problem. For ACO, however, solution time has remained very short despite the increasing size of the problem.
The diagram of Figure 3 illustrates the objective function values obtained through the exact solution (GAMS) and metaheuristic solution (ACO) of the model.
In Figure 3, it can be seen that there is not much difference between the qualities of solutions produced by the two methods. Overall, the solutions obtained from ACO are only 0.05% different from the exact solutions of GAMS, which demonstrates the good performance of this algorithm in finding the optimal solution of the formulated problem.
In general, ACO performed very well in terms of solution speed as well as solution accuracy, and can therefore be considered a suitable method for optimizing this problem.
4. CONCLUSION
This study presented a new mathematical model for production planning, maintenance, and labor scheduling. The presented model makes use of a new objective function, which divides the labor force by the production output. Basically, this mathematical model seeks to reduce the labor employed while also increasing the production output. In this model, the parameters related to labor and maintenance are expressed as fuzzy numbers. The fuzzy cut method was used to convert the fuzzy model into its nonfuzzy form. The ant colony optimization algorithm (ACO) was used to solve this mathematical model. The numerical results showed the good performance of this algorithm in terms of solution quality as well as solution speed. The results obtained in the model validation process and when solving numerical problem instances showed that the approach taken in this study, i.e. combining different domains of decision making in production units, can yield good strategies to enhance production planning, maintenance, and labor scheduling all together. To expand on this research, it is recommended to use a robust optimization approach to deal with uncertainty, and also try other metaheuristic methods such as the firefly algorithm or the grey wolf optimization algorithm.