Credibility based chance constrained programming for parallel machine scheduling under linear deterioration and learning effects with considering setup times dependent on past sequences

,


Introduction
The scheduling of parallel machines is a vital problem that has significant implications not only in the production environment but also in distributed computing and parallel information processing (Alchieri et al., 2018). This scheduling issue represents a critical challenge in production management and operations research due to its combinatorial nature and multiple constraints that must be considered, including machine availability, job processing times, and release dates (Alchieri et al., 2018;Meyr, 2002). In addition, real-world situations may feature additional variables, such as learning and deterioration effects, as well as machine delivery times, that render the problem even more intricate (Kim & Lee, 2021). Developing effective scheduling algorithms that take these factors into account is imperative for ensuring optimal production processes and meeting customer demand (Rostami & Sabripoor, 2022). Uncertainty about influential parameters represents one of the most fundamental issues in real-world problems and environments, where parameters such as processing times and delivery deadlines frequently exhibit a non-deterministic state. The use of fuzzy numbers can help to capture this uncertainty and ambiguity, although it is worth noting that other factors in the workplace can impact problem conditions and basic parameters. In a production environment, raw materials that are queued for processing are often subject to a "deterioration effect of jobs", which may ultimately result in their complete loss (Abedi et al., 2020;Wang & Edwin Cheng, 2007;Wu et al., 2019). Conversely, the abilities and skills of production workers tend to improve over time, thereby accelerating the production process -an effect known as the "learning effect" (Mor et al., 2020;R. Wu et al., 2018). The setup times of machines which are sequence-dependent, constitute another critical parameter that impacts the completion time of tasks in a production system. This type of setup time is typically influenced by the conditions under which it is performed, with one of the more common types being setup time dependent on the past sequence. It is essential to acknowledge that machines within the workplace may exhibit differences in their processing capabilities, stemming from variances in their specifications or brands (Liu et al., 2006;Siemieniuch et al., 2015). In addition to the discourse on workplace conditions, production goals are always a crucial consideration. These objectives can be viewed from the producer's perspective, such as minimizing the maximum completion time of tasks, or from the consumer's perspective, such as minimizing the total delay and early completion. In the literature, the latter two objectives are frequently merged into a single objective termed ET, which corresponds to the JIT concept in the production environment (Pinto et al., 2018). Based on our comprehensive analysis of resources and references, the problem of scheduling non-identical parallel machines that simultaneously considers the effects of task deterioration and learning in a fuzzy environment has yet to be investigated.
The structure of the sections is as follows: In section 2, the necessary basic principles for investigating the problem are presented. Section 3, while defining the problem, a Non-Linear Multi-Objective Mixed Integer Non-Linear Programming (MOMINLP) is first presented. Then, with the help of the credibility index and also chance-constrained programming (CCP), the presented model is made deterministic. Section 4, includes the development of metaheuristic hybrid algorithms for solving the problem under large dimensions. Section 5 shows the computational results and finally, section 6 includes conclusions and future studies.

Literature review
Parallel machine scheduling in deterministic environment without considering the effects of deterioration and learning has been the focus of many researchers, including (Brier & lia dwi jayanti, 2020;Prot et al., 2013;Sariçiçek & Çelik, 2011;Unlu & Mason, 2010;Vallada & Ruiz, 2011). In recent years, one paper used the SA and GA algorithms to study the problem of parallel machine scheduling in a deterministic environment (Wang et al., 2020).
Several researchers have studied this problem considering the effects of task deterioration. (Ji & Cheng, 2008; examined the problem of parallel machine scheduling with the goal of minimizing the total completion time under the condition that task deterioration follows a linear function. (Mazdeh et al., 2010)considered the problem of minimizing the sum of tardiness and the cost of machine deterioration in parallel machine scheduling. (Ruiz-Torres et al., 2013) studied the same problem under the condition that the task deterioration rate is a function of the job sequence, and solved it using SA. Additionally, (Ding et al., 2019;S. J. Yang, 2013;Zhang & Luo, 2013) investigated the problem of parallel machine scheduling under the conditions of task deterioration, and considering factors such as maintenance and rejection. In recent years, (Sekkal & Belkaid, 2020) solved the problem of parallel machine scheduling considering the effects of deterioration and resource constraints using the SA algorithm.
Considering the effect of learning alone in scheduling problems, (Lee et al., 2012) investigated the problem of uniform parallel machine scheduling, while (M. Liu, 2013) considered the joint effect of learning and delivery times of goods. (Przybylski, 2018) took into account the effect of integral-based learning, and (Expósito-Izquierdo et al., 2019) added sequence-dependent setup times. In recent years, the problem of single machine and parallel machine sequencing with learning effect has been solved and evaluated using three algorithms, SA, GSO, and SA-GSO. In another study, (Toksari & Güner, 2009) investigated the problem of scheduling identical parallel machines with the aim of minimizing the total earlinesstardiness and non-linear deterioration coefficient under deterministic environment and common due dates. Another research focuses on the parallel machine scheduling problem, aiming to minimize the total completion time under the influence of position-dependent learning and linear deterioration coefficient, which is solved using a traditional genetic algorithm. Despite numerous studies in the field of parallel machine scheduling under deterministic environments, only a few limited articles have been published in this area under uncertain conditions. (Al-Khamis & M'Hallah, 2011) solved the basic problem of scheduling tasks on parallel machines as a two-level mathematical model without considering the effect of learning and deterioration. (Peng & Liu, 2004) studied the basic problem of parallel machine scheduling under fuzzy processing times and developed three methods for modeling probabilistic problems, namely, EVP, CCP, and DCP, in a fuzzy environment with a new metric called credibility. Finally, they evaluated their proposed methods using simulation.
Recently, two articles have been published in the fuzzy environment, taking into account the effects of learning and decay of tasks and without considering the sequence-dependent setup times. The first article, written by (Rostami et al., 2015), investigates this problem for the first time. They propose an algorithm based on important features for solving optimization problems that consider the effects of learning and decay of tasks. In addition, (Arık & Toksarı, 2018) also examine this problem in the fuzzy environment, considering the effects of decay and learning and solving it using a local search algorithm According to our search and research in this field, so far no research has been conducted on scheduling parallel machines considering the learning effect and the decay effect of tasks simultaneously under conditions where processing times and delivery times are fuzzy and preparation times are dependent on the past sequence in the problem. In this article, the problem of scheduling non-identical parallel machines in a fuzzy environment with the aim of minimizing both the total earliness and tardiness (ET) and the maximum completion time (makespan) is investigated. In this problem, the learning effect and the decay effect simultaneously affect the problem conditions, which can be taken into account in the processing times of tasks. We assume that both the decay effect and the learning effect are functions of task processing priorities.

Mathematical model development
This section provides an overview of the problem and the fundamental principles of fuzzy programming, as well as the index used, and finally, the Augmented ε-constraint method is discussed.

Problem description
There are N tasks that must all be completed at time zero and failure is not allowed. In a production workshop, there are M machines that process tasks at different speeds (Vj). Each machine is operated by workers with a different skill level, represented by a negative coefficient βj. The tasks also deteriorate over time at different rates, represented by another coefficient αi. To start processing each task on each machine, we need to spend some preparation time, which is a sequence of past events (see (Hsu et al., 2011) for more details). Processing times and delivery deadlines are also fuzzy triangular numbers. The objective is to minimize the total earliness and tardiness (ET) and makepan. According to the standard introduced by (Graham et al., 1979), the problem can be formulated as follows: where LE, D, and ST represent the effect of learning, task deterioration, and preparation time, respectively, which depend on past events. Since problem Pm||C max has been proven to be NP-hard by (Liao & Sheen, 2008), adding more difficult conditions to the problem still leaves it NP-hard.

Proposed Mathematical model
In this section, Fuzzy Multi-Objective Mixed Integer Non-linear Programming (FMOMINLP) is presented. Firstly, it is assumed that the sequence for each machine is divided into N hypothetical priorities denoted by r. Then, the task assignment process to these priorities is performed using the model. The model searches the solution space to ultimately reach the optimal solution. Table 1 summarizes the parameters and decision variables of the proposed model.  Earliness of priority r related to machine j xirj 1 if task i is assigned to priority r related to machine j; 0.
MOMIP model: Objective function (2) minimizes the total tardiness and earliness. Objective function (3) minimizes the maximum completion time of jobs. Constraint (4) is designed to calculate the actual processing times of jobs assigned to each priority on each machine. It is clear from this relationship that as the priority number increases, the impact of deterioration on each job increases and processing time increases. Also, as the value of r increases, the effect of learning plays a role and reduces the processing time of that priority. Constraint (5) is to calculate the sequence-dependent setup times on each machine. It is clear that on each machine, the setup time also increases with an increase in r. Constraint (6) is to calculate the completion time of each priority from each machine, which is equal to the completion time of the previous priority on the same machine plus the setup and actual processing times of that priority. Constraint (7) is to convert the Min Max objective function to the simplified Min and calculate the makespan. Constraint (8) calculates the tardiness and earliness of each priority. Constraint (9) requires that each job be assigned to only one priority of one machine. Constraint (10) are designed to calculate the absence of more than one job in each priority of each machine, as well as the absence of idle time on machines. That is, if a job is assigned to a priority of a machine, then no other job can be assigned to that priority from that machine, and there must be a job in its previous priority as well. Finally, constraint (11) defines decision variables. It should be noted that this model is a nonlinear model because it has a series of nonlinear constraints. In this model, constraint (8) is a nonlinear constraint due to the multiplication of two decision variables. So, the solving software does not guarantee reaching the global optimal solution.

Fuzzy Basics
Fuzzy theory was first introduced by (Zadeh, 1965). In this section, the prerequisites for fuzzy numbers are defined, which will be used in subsequent sections (Tanaka, 1992). These definitions and properties are as follows: are two triangular fuzzy numbers, then their sum is also a triangular fuzzy number (Chanas & Kasperski, 2001): Proposition 2: If ξ  is a fuzzy number, r is a crisp number, and Cr represents credibility, Pos represents possibility, and Nec represents necessity, then the relationship between the three measurement tools, i.e., possibility, obligation, and credibility, for an event r ξ ≤  is as follows (B. Liu & Liu, 2002): Proposition 3: According (Zhu & Zhang, 2009)have proven, if 0.5 λ > is a triangular fuzzy number and r is a crisp number, then the product and sum of two fuzzy numbers will be less than or equal to r: where (2) ξ and (3) ξ are the middle value and right side value of the triangular fuzzy number ξ  , respectively. ( , , ) C c c c = whose values i c can be calculated using the following equations (Chen & Wang, 2008): where in Eq. (4):

Credibility-based on Chance Constrained Programming (CCP)
In this section, the model presented above is transformed into a deterministic model by using credibility measurement tools and consulting with CCP. In the model below, the concept of fuzzy chance constrained programming is used to make it deterministic (Peng & Liu, 2004). The difference with the probabilistic case is that instead of using a probability function, a credibility value function for an uncertain event is used. An important point in the model below is that when calculating ET, one should be careful to calculate ( for tasks with a late start and ( for tasks with an early start. Therefore, to establish such conditions, the concept of fuzzy distance between two numbers expressed in proposition 4 has been used. In the model below, two new variables 1 f and 2 f are defined to model the problem with the help of CCP, which are the upper bounds of objective functions 1 and 2, respectively. Also, instead of two variables rj E  and rj T  , a new variable called rj ET   , which is a combination of these two variables, has been used. The two parameters 1 λ and 2 λ are also the confidence levels of the credibility functions of the first and second objective functions.
Credibility-based CCP model:

Constrains
(4), (5), (6), In this model, the objective function (18) minimizes the two objective functions in order. Constraints (19) and (20) are used to define chance constraints in fuzzy problems, considering the credibility value of the two objective functions greater than a predetermined confidence level using credibility tools and seeking the minimum values of i f that satisfy these equations.
These two constraints can be transformed into deterministic constraints with the help of the results in proposition 3. In equation (21), the right-hand side expression represents the distance between two fuzzy numbers, i.e., , r j C  and i d  The distance value for these two numbers can be calculated using the equations in proposition 4, resulting in a triangular fuzzy number. It should be noted that in the above model, it is assumed that the confidence levels are values greater than 0.5, i.e., 0.5 i λ > .

Augmented ε-constraint method
In a multi-objective optimization problem, the Pareto front defines the set of solutions that achieve Pareto optimality. Thus, there is no single optimal solution, but rather a set of optimal solutions representing different trade-offs among the objectives (Jabir et al., 2015;Khoshabi et al., 2020;Petchrompo et al., 2022). In fact, this means that the performance of each objective function cannot be improved without worsening the results of at least one other objective function in the problem (Ghanbari et al., 2022;Spielhofer et al., 2022).
The ε-constraint method is a popular approach suggested and frequently employed for tackling problems with multiple objectives. (Hosseini et al., 2016;Javadi et al., 2020;Laumanns et al., 2006;Z. Yang et al., 2014). This method has some disadvantages, for example, when the number of objectives increases, Sensibility, to the initial conditions, etc (Hemmatesfe et al., 2017;Stanovov et al., 2020). becomes complicated. In order to resolve this problem (Mavrotas, 2009), introduced the augmented ε-constraint method. Numerous studies and articles provide evidence for the efficacy and enhancement resulting from this approach (Esmaili et al., 2011;G & K, 2020;Huy et al., 2023;Torabzadeh et al., 2022;Xing et al., 2019).According to this method, using lexicographic optimization, a payoff table is formed and the two-objective model presented is transformed into a new problem, which is shown in model (24).
subject to It should be noted that the range of objective ( ) is calculated from the payoff table ( ( ) ∈ [ ( ) , ( ) ]), and similarly, the range of my objective is calculated from the payoff table ( = ( ) − ( ) ). In addition, the range of objective ( )should be divided by , which is calculated from equation (25).
Additionally, a small value of ε is considered between 10 and 10 .

Developed Metaheuristic Algorithms
In this section, first the developed NSGA-II algorithm stages for the problem under investigation in this paper are explained, and then the combination algorithm with Variable Neighborhood Search (VNS) will be discussed. The implementation stages of the NSGA-II algorithm can be described as follows:

I. Representation of Solutions (Genotype Space)
The type of solution representation in this paper is considered as two-row matrices with n columns, where n represents the total number of available tasks. The first row of the matrix indicates the machine number that each task (column number) is processed on, and the second row represents the processing priority of each task on that machine. For example, Figure 1 illustrates the real space of the problem and its chromosome representation for 5 tasks and 2 parallel machines.

II. Definition of Fitness Function
To determine the value of each created solution (chromosome), a specific fitness function based on the Pareto rank and crowding distance needs to be developed. However, before that, necessary processing has been carried out to determine the precise value of each objective function due to the fuzzy nature of the objective functions. After calculating the value for each solution i (using the mathematical equations of objective functions provided in section 3.3), the crowding distance of each solution can also be calculated .

III. Implementation of Crossover Operation
In order to perform the crossover operation, a single-point crossover operator has been used. First, one of the columns of the chromosome is randomly selected, and then the columns before the selected column are taken from parent one and the columns after are taken from parent two to form the first offspring.

IV. Implementation of Mutation Operation
To perform mutation on chromosomes, first * mut P n columns are randomly selected (although this value is rounded up to the nearest integer). Then, the values of the selected columns' components are randomly chosen according to a uniform relationship, adopting an integer value between 1 and the maximum value of each row. 0.03 mut P = is the value in this algorithm.

V. Genome Encoding
Given the selected representation structure, almost all problem constraints will always be satisfied. However, constraints related to the order of priority values on each machine may be violated after crossover and mutation operations. To solve this problem, repair strategy is employed where after each execution of crossover and mutation operators, infeasible solutions are repaired. The following steps are taken to repair the solutions: Step 1: The priority values assigned to the tasks of machine j are sorted in ascending order.
Step 2: The rank of each task in the sorted arrangement of Step 1 is considered as its new priority value.
Step 3: The above procedure is applied to all machines to repair all solutions.

VI. Termination Criterion
One of the well-known conditions in literature is the condition of achieving the maximum predetermined number of repetitions, which is also considered in this paper.

VII. Combination with VNS
In this section, we present the features of the VNS algorithm for its integration with the NSGA-II algorithm. This combination can be classified as a Low-level Co-evolutionary Hybridization (LCH) combination. In this combination, after generating solutions using crossover operators, a subset of chromosomes will be influenced by the VNS algorithm with two neighborhood structures (NSS) to extract a better solution (local optimum near them). To evaluate the solutions obtained by the VNS algorithm, a one-sided intensification approach based on the concept of exploiting is used. Exploiting means that 50% of the time, VNS attempts to find solutions with the lowest primary objective function, and 50% of the time, it tries to find solutions with the lower secondary objective function. Therefore, it is expected that the optimal Pareto front will be drawn more tightly, and a better dispersion criterion will be demonstrated. In this algorithm, two neighborhood structures are proposed as follows: First structure: Randomly select two machines, randomly select one job on each machine, swap the two jobs with each other. Second structure: Randomly select two machines, randomly select one job on each machine, swap the two jobs with each other, randomly select another job on each machine, swap the two jobs with each other again. Also, a stopping criterion is defined by specifying a fixed number of iterations. Two options for the number of iterations, 50 and 100, are presented, and 100 is selected as the better option in the parameter settings. Moreover, considering the combination of VNS algorithm with NSGA-II algorithm, other parameters should also be defined, including: VNS algorithm start time (VNSstart), which in this study, will be called after one-third of all NSGA-II iterations, then the VNS algorithm will be invoked. The step length of the VNS algorithm (VNSstep), which means that after starting the VNS algorithm, it will be executed three times. The number of chromosomes (VNSnumber), which means that in each VNS execution, a percentage of the total number of generated offspring will be used for local improvement. This parameter is set to 15% in the parameter settings.

Computational results
To validate and achieve the efficiency of the mathematical model and algorithm presented in the previous section, in this section, randomly generated problems are solved using these methods to discuss and investigate the performance of the presented methods. The reason for generating problems randomly is that the problem under study in this paper has not been investigated in the literature with all the conditions defined simultaneously, so there is no sample problem for it. In randomly generated problems, fuzzy processing times are selected for each operation randomly from uniform distribution function ( ) The decay rate is selected from the uni- In this paper, the confidence level of the first and second objective functions is considered 0.95 and 0.9, respectively. The problems created in this paper are designed in three categories of small, medium, and large sizes. The mathematical programming model is modeled and solved using GAMS software. In the NSGA-II algorithm, the number of iterations is considered 200. The mutation rate, crossover rate, and population size parameters are considered 0.15 m P = ‫و‬ 0.9 C P = ، NP=20 respectively, after adjusting the parameters. The NSGA-II algorithm is coded using MATLAB software. The workflow is as follows: first, a problem is created, then it is solved for each level using the mathematical model, and the Pareto optimal front is obtained. Then, the same problem is solved using the NSGA-II algorithm and the Pareto optimal front is obtained. Finally, the results are compared and the efficiency of the proposed methods is investigated.
The workflow is as follows: first, a problem is created, then a mathematical model is solved for each level and the optimal Pareto front is obtained. Then the same problem is solved using the NSGA-II algorithm. It should be noted that the mathematical model does not guarantee reaching the global optimal solution due to its nonlinearity, so some points on the Pareto front obtained by the mathematical model may be dominated by the solutions obtained from the metaheuristic algorithm. Therefore, two criteria are used here to evaluate the solutions. The first criterion is the error of the Pareto points obtained by algorithm k from the optimal Pareto front obtained by another algorithm (GDk). This index is calculated using equations (23) and (24). It should be noted that in these equations, the value is only applicable to the points that are dominated by other points. In other words, is the number of created Pareto solutions. In these formulas, | | known PF , is the index corresponding to the points created by algorithm k and j is the index corresponding to the points created by another algorithm.
The second criterion is the coverage index. If we denote the Pareto front obtained from the mathematical model by X and the Pareto front obtained from NSGA-II by Y, then the criterion C(X,Y) represents the number of points where the front Y dominates the front X. Table 2 shows the results of comparing the mathematical model and NSGA-II for small-sized problems. The results indicate that NSGA-II has a much higher speed compared to the developed epsilon constraint in generating the Pareto front. However, in this category of problems, the mathematical programming model creates better solutions in terms of both the GD and coverage criteria. A significant point in this table is that for problems with three machines, some of the Pareto points of the metaheuristic algorithm outperform some of the Pareto points of the mathematical model (according to the value of C), indicating that the nonlinear model has failed to create the global Pareto front in this category of problems.
Except for problems with small sizes, all other problems can only reach a suitable Pareto front through NSGA-II. Therefore, to evaluate the obtained solutions using the metaheuristic approach, it is necessary to compare them with solutions that are very close to the Pareto optimal solutions. In this paper, this very good solution is obtained with the help of a large number of NSGA-II algorithm runs and a large population size, and is referred to as the best-known solution (BKS). For obtaining the BKS, the number of iterations is set to 1000 and the population size is set to 100. Table 3 shows the results of these evaluations. In Table 3, the number of population members, solution times, the number of Pareto points created by the algorithm, the spacing metric, the ranges, and average objective functions, GD metric, and finally the average errors of the generated solutions are presented. The average error metric actually calculates the ratio of GD compared to the distance between the average point of two objective functions from the coordinate origin. This metric is used because the values of GD are proportional to the values of objective functions. Therefore, if the values of objective functions are numerically large, then the value of GD also becomes large, which cannot be a suitable metric alone for comparison. The error value is calculated from equation (25): The spacing metric is also calculated using equations (26) and (27) In this table, all metrics have been calculated relative to BKS. The results show that as the size of the problems increases, the solution time increases, the spacing metric decreases, and the error rate remains almost unaffected by the problem size.
Finally, by considering the VNS algorithm and running the algorithm for the case of 16 jobs and 4 machines, Table 4 compares the results with each other. As can be seen from the claim made, the dispersion value of the partial solutions in the one-sided intensification approach is higher than the normal approach.  Furthermore, to analyze the partial solutions created in each of the approaches, Fig. 3 shows the Pareto front resulting from the combination of all partial solutions in three repetitions. As can be seen from the figure, the one-sided intensification approach demonstrates significantly better performance, especially in the corner regions of the Pareto front. On the other hand, to analyze the convergence speed, Fig. 4 is proposed. In this plot, the vertical axis represents the time of birth for each of the final Pareto solutions, and the horizontal axis represents the number of final Pareto solutions. Therefore, the more the plot tends towards a 100% ratio, the slower the convergence and the slower the stabilization of final Pareto solutions. As can be seen from the plot, a significant percentage of partial solutions in each of the approaches have a birth time of more than 80% of the algorithm's execution time, indicating a lack of early convergence in the proposed model.

Conclusions and future research directions
In this paper, we examine the problem of scheduling non-identical parallel machines in a fuzzy environment, considering job degradation, learning effects, and order-dependent setup time effects that simultaneously affect processing time. This article first presents a preliminary definition of the problem and a literature review. A mixed-integer nonlinear mathematical model is then proposed to represent order processing time and delivery dates using triangular fuzzy numbers. Minimizes two objective functions, total early and total late, along with maximum job completion time. The multi-objective fuzzy model is then transformed into a deterministic model using the method of Chance Constrained Programming (CCP) with a measure of reliability that has properties close to the concept of probability. Multi-objective fuzzy models are transformed into deterministic models that can be coded in solver software. Since the proposed mathematical model cannot solve largedimensional problems in a reasonable amount of time, we use the NSGA-II metaheuristic algorithm to solve larger-sized problems. Furthermore, to enhance the power of the NSGA-II algorithm, a hybrid state of this algorithm is presented in the form of a multi-objective genetic algorithm based on the concept of one-way reinforcement.
To evaluate the proposed methods in this article, 20 completely random problems were generated in three sizes, namely small, medium, and large. For small-sized problems that both methods have the ability to solve in a reasonable time, Pareto sets were created by two methods and compared with each other. The results showed that in small-sized problems, the mathematical model creates relatively better solutions than NSGA-II, although in some problems, the mathematical model has weaker performance than NSGA-II. By solving medium and large-sized problems using NSGA-II and comparing the obtained solutions with the mathematical model, it is shown that the proposed hybrid NSGA-II algorithm outperforms the mathematical model in terms of solution quality and computational time.
For further research, considering the effect of non-linear deterioration instead of linear and measuring its difference in the problem use other methods to solve small-scale multi-objective problems. The use of another metaheuristic algorithm and their comparison with the method of solution presented in this paper.