Generating Pareto Optimal Solutions for Multi-Objective Optimization Problems Using Goal Programming

: Goal programming (GP) is a well-known multi-criteria decision-making tool that is supported by a network of practitioners and researchers who aim to develop its mathematical foundation to cover a wide range of applications. The popularity of GP models stems from their structure, which is based on a satisfying philosophy. This philosophy takes into consideration the preferences of the decision-maker concerning the model parameters. Therefore, the GP model provides the decision-maker with one satisfactory solution that reflects the trade-off between competing objectives. Nevertheless, there is no guarantee regarding the efficiency of this solution. Consequently, this study is designed to improve the quality of decision-making processes by addressing the efficiency issue with the solutions of GP models. The main contribution of this paper is to improve the mathematical framework of the GP model so that it can generate a set of Pareto optimal solutions rather than just one solution. This allows stakeholders to have a complete picture of the feasible space of solutions and select the solution that represents the best compromise according to their preferences. As a result, the proposed methodology is called generational GP. In addition, the study enhances the quality of GP solutions by integrating the notion of the hypervolume subset selection problem with the suggested technique. This, in turn, overcomes the efficiency problem of GP solutions. The performance of the proposed method has been validated through an application to the flow shop scheduling problem. However, our modeling approach is useful for decision-makers in different fields of study. Finally, the merits of the generational GP method are highlighted, with a strong emphasis on potential areas for future research.


Introduction
Most real-world applications take into account multiple competing objectives.There are two primary categories of methods for solving multi-objective optimization problems (MOPs) [1].The first category includes evolutionary algorithms.The second category involves a range of classical techniques.These methods are sorted into three different groups depending on the phase in which decision-makers participate in the decision-making process [2].These groups are known in the literature as a posteriori (generation), a priori, and interactive methods.
The first group involves techniques that generate a set of efficient solutions for the decision-maker.Therefore, the intervention of the stakeholder is limited to the selection process for solutions.The epsilon constraint method is one of the most common techniques in this group.On the contrary, the methods in the second group, i.e., a priori methods, allow decision-makers to express their preferences before starting the solution process.Goal programming (GP) is a popular technique in this group.In GP, model parameters, such as the aspiration levels of goals, are determined based on the preferences of the decision-maker before solving the optimization problem.In interactive methods, stakeholders interact after each step of the solution process to show their preferences [2].Despite the variety of multi-criteria decision-making (MCDM) tools to solve MOPs, this study focuses on developing the mathematical framework of GP to provide the decision-maker with the solution of the best compromise.
The GP technique is one of the most widely used MCDM approaches because of the simplicity of its mathematical framework [3].This framework takes into account the simultaneous optimization of competing objectives.Therefore, the GP model provides the decision-maker with a solution that reflects the trade-off between several conflicting objectives.Moreover, the GP approach is based on the satisficing philosophy, which allows the decision-maker to represent his or her preferences concerning the model parameters.Consequently, the quality of the decision-making processes is enhanced through interacting with the stakeholder before solving the GP model.
The mathematical structure of GP models is based on minimizing a distance function that considers the deviations between the objectives and their target values, which are determined by the decision-maker.The first GP model was initially presented by Charnes et al. [4], and it was considered an extension to linear programming models.Furthermore, several variants of GP models were suggested in the literature to represent the decision-maker's preferences.The most common variants are weighted GP (WGP) and lexicographic GP (LGP) [5].On the one hand, WGP focuses on minimizing the weighted sum of the deviations from all objectives.These weights reflect the relative importance of each objective.On the other hand, the LGP models allow the decision-maker to rank objectives according to their importance.The deviations from the goals at the highest priority level are minimized first.Goals at the next level of importance are optimized while maintaining the minimal values of goals at higher priority levels.Consequently, this variant is based on solving a series of sequential optimization models [6].
Despite the popularity of GP and the variety of its variants to solve different MOPs, there is no guarantee that it provides the decision-maker with an efficient solution.The efficiency issue happens due to the mathematical structure of GP, which is based on the satisfying philosophy [7].The essence behind this philosophy is to provide the decisionmaker with a solution that is as close as possible to the desired target value.This is opposed to the optimization philosophy, which seeks to find an optimal solution.Consequently, several studies were introduced in the literature to handle this problem.The contributions of these studies are classified into two main directions [8].The first one focused on testing only the efficiency of a GP solution, such as the technique introduced by Cabllero et al. [9].The second type of study was designed to provide the decision-maker with an efficient solution in the case of the failure of an efficiency test [8,10,11].
This study aims at developing the mathematical framework of GP models to avoid introducing a dominant solution to the decision-maker.In addition, the study suggests a GP model that generates a set of efficient solutions rather than just one solution.This set covers the feasible space of potential solutions.This, in turn, allows the stakeholder to have a variety of efficient candidate solutions and select the one that ideally reflects his or her preferences.Consequently, the contribution of this paper is twofold.Firstly, the suggested technique is able to produce a set of GP solutions rather than just one solution.Therefore, the framework of GP is adjusted to be similar to the a posteriori approach.This implies that the proposed generational goal programming (GGP) methodology converts the GP technique to a generation method while maintaining the preferences of the decision-maker.The second contribution of this paper is to improve the quality of GP solutions by overcoming the efficiency issues related to them.This is accomplished by combining the proposed technique with the concept of the hypervolume subset selection problem (HYPssp).This allows stakeholders to have a subset of non-dominated solutions that capture the most diverse and well-distributed solution points in the objective space.To assess the reliability of the proposed method, it has been applied to the green permutation flowshop scheduling problem (GPFSP).This application has practical benefits in the fields of manufacturing and industrial engineering, and it is used as an example to show the usefulness of the suggested technique.Nevertheless, the proposed GGP methodology can be compared to the current optimization methods in the literature in addition to applying it to other applications, such as resource allocations among level crossings [12], the vehicle routing problem [13], the ship scheduling problem [14], the complex building design problem [15], and the vessel scheduling problem [16].
The rest of this paper is arranged as follows: Section 2 introduces the proposed GGP method.Section 3 gives an overview of the GPFSP.Moreover, it presents the mathematical programming models used in this study.In Section 4, data settings and computational results are illustrated to test the validity of the proposed technique.The last section concludes the paper.

The proposed GGP methodology
In this study, the suggested GGP technique is introduced in the context of bi-objective minimization problems.However, it can be generalized to any MOP.The purpose of a two-dimensional minimization problem is to find a vector of decision variables in the decision space Χ that minimizes a vector of conflicting objectives in the objective space, i.e., 2 : f Χ →  .Furthermore, this study is interested in points in the objective space that form a subset of 2 .
∈  weakly dominates another point q = (q 1 , q 2 ) 2 1 2 ( ) , q q q = ∈  (also known as p p q  q) if and only if for 1, 2 . Therefore, the desired output from any technique used to solve MOP is to obtain non-dominated solutions.These solutions are incomparable and are used to determine a Pareto frontier [18].
The proposed GGP methodology relies on several steps.The first step is based on calculating the minimum and maximum values for each objective.On the one hand, the minimum value of each objective is considered an ideal value since the suggested technique is illustrated in the context of bi-objective minimization problems.Ideal values are obtained by running two single-objective optimization models.On the other hand, a maximum value is regarded as the worst value that an objective function reaches.Maximum values are estimated in this study by adopting the idea of redundant goals in LGP [19,20].
The issue of goal redundancy occurs in the lexicographic variant of GP when the target values of goals at the highest priority level are equal to their ideal values.This results in redundancy in goals placed at the lowest priority level, which implies that these goals are not achieved [20].This also implies that the goals at the lowest priority level obtain their worst possible values.Consequently, this concept is adopted in this study to estimate the worst value for each objective.On the one hand, the maximum value of the first objective is obtained by placing it at the second priority level in the LGP model.Due to choosing the ideal values of both goals as their target levels, the model achieves the goal at the first priority level, i.e., the second objective, and produces the worst value of the goal at the second rank, i.e., the first objective.On the other hand, the worst value of the second objective is achieved by running the same LGP model after reversing the previous priority order.
After computing the maximum and minimum values for each objective, the next step aims at constructing a group of intervals for each objective to generate a set of aspiration levels for the LGP model.To create these intervals, the range of each objective is determined as the difference between the maximum and minimum values.The length of each interval is the range divided by the number of intervals (n), which is determined by the decision-maker.Furthermore, this number is equivalent to the number of aspiration levels since one aspiration level is randomly generated from each interval.Therefore, the set w of n points of aspiration levels is obtained from these intervals.The x-coordinate of each point represents an aspiration level for the first objective, while the y-coordinate is an aspiration level for the second objective.Moreover, the intervals are constructed such that these points have an increasing order of the first coordinates and a decreasing order of the second coordinates.
After running the LGP model n times to obtain a set of GP solutions, the last step in the proposed methodology is to select a subset of non-dominated solutions.There are several criteria for choosing this subset.This study follows the approach presented by Bringmann et al. [17] to obtain a subset of incomparable solutions.The authors suggested the HYPssp algorithm, which produces a subset of efficient solutions based on the maximum value of the hypervolume performance metric.This indicator, which was introduced by Zitzler and Thiele [21], is designed to measure the volume of the objective space covered by a set of solutions according to a reference value.This paper follows the study of Hughes [22], which estimated the coordinates of the reference point as the maximum (worst) value for each objective.
There are several reasons behind combining the HYPssp algorithm with the proposed methodology.Firstly, GP does not necessarily produce a Pareto-optimal solution [8].As mentioned in the previous section, GP may provide the decision-maker with an inefficient solution since its mathematical structure follows the satisfying philosophy rather than the optimizing one.Secondly, and most importantly, rational decision-makers do not choose a dominant solution [8].Accepting a dominant solution occurs in the traditional variants of GP because other efficient solutions are unavailable to stakeholders due to a lack of knowledge about the feasible region.Consequently, adjusting the structure of GP to generate a set of solutions is considered fundamental for stakeholders who are interested in exploring the feasibility of GP solutions.The HYPssp algorithm ensures that inefficient GP solutions are avoided.Therefore, the decision-maker selects from this subset the most suitable solution for the MOP under consideration.
The proposed GGP methodology can be summarized in the framework of bi-objective minimization problems as follows: 1. Calculate the minimum (ideal) value of each objective using a single-objective optimization model.2. Compute the maximum (worst) value of each objective by considering the redundancy issue of LGP.
3. Set intervals for each objective, and choose an aspiration level at random from each one.4. Construct the set w of points for aspiration levels.5. Use these target values for the goals in the LGP model.6. Run the HYPssp algorithm to obtain a subset of efficient GP solutions.
The advantages of the suggested technique are highlighted by applying it to the GPFSP.This application has practical benefits in the field of industrial engineering.The next section introduces this application and explains in detail the steps of the GGP methodology.

Introduction to the GPFSP
The traditional permutation flowshop scheduling problem (PFSP) is a combinatorial optimization problem that takes objectives related to production time into account.The maximum completion time, i.e., makespan, flow time, and tardiness, are the most common production efficiency-related objectives in the literature.The purpose of the PFSP is to identify the optimal order of processing jobs on machines.Additionally, it makes the assumption that jobs will be processed on machines in the same order [23].Other assumptions include the independence of jobs and the availability of jobs at the start of processing time.Moreover, machines are independent, and their interruption is not allowed during the processing stage.These assumptions are considered in this study.More details about the PFSP and its requirements are available at [24].
The GPFSP is considered an extension of the traditional PFSP.This recent variant considers energy efficiencyrelated objectives, such as energy consumption and carbon emissions.In this study, makespan and total energy consumption (TEC) are considered the competing objectives, and they are assessed using the proposed methodology.Moreover, the speed of machines is the factor that creates the conflict between these objectives.Processing jobs at a higher speed decreases makespan.However, this increases the energy consumed by machines.
Several optimization techniques were introduced to evaluate the trade-off between production and energy efficiency-related objectives in the context of the speed-scaling strategy.Mansouri et al. [25] used the epsilon constraint method and developed constructive heuristics to study the compromise between makespan and TEC for a two-machine sequence-dependent permutation flowshop scheduling problem.Moreover, the compromise between the total flow time and TEC was assessed in the context of the GPFSP [26].The authors employed the augmented epsilon constraint technique to generate Pareto-optimal solutions for small-scale problems.To cope with the complexity of this problem, their study proposed multi-objective iterated greedy algorithms and variable block insertion heuristics for large-scale problems.Furthermore, Saber and Ranjbar [27] studied the conflict between the total tardiness and the total carbon emissions for the GPFSP.The authors suggested a mixed integer mathematical programming model and a multiobjective decomposition-based heuristic algorithm to solve this problem.Recently, in the context of task scheduling for parallel systems, Stewart et al. [28] developed a mixed-integer mathematical programming model to minimize makespan and energy consumption by varying the speed of processors.The authors used the epsilon constraint and the weighted sum scalarization methods to solve this problem.
The next subsection introduces the mathematical programming models used to study the GPFSP.

Mathematical programming models for the GPFSP
As mentioned in Section 2, the proposed GGP methodology relies on calculating the minimum and maximum values for each objective.Model 1 below is a single-objective mixed-integer mathematical programming model that is used twice to compute the minimum (ideal) values for makespan and TEC.These ideal values are used as target levels in the LGP model, i.e., Model 2 below.Before introducing these models, Table 1 below presents the notations, parameters, and decision variables used.It is worth mentioning that Model 1 below is similar to the model introduced by Amiri and Behnamian [29].However, their model is designed to assess the stochastic variant of the GPFSP under scenario analysis.The deterministic version of their model is used in this study.
Model 1. Deterministic mixed-integer mathematical programming model for GPFSP subject to { } 0,1 , , , , The objective functions (1) and (2) aim at minimizing makespan and TEC, respectively.Each objective function is minimized individually to get its ideal value.This implies that Model 1 is used twice.Constraint (3) states that the completion time of the job in the first position on the first machine should equal the processing time on that machine.The inequality is used in the case of machine idle time.Constraint (4) guarantees that a job in sequence position j cannot end its processing on the current machine unless it has already finished its processing on the previous machine.Constraint (5) states that a job in processing on machine m can move to the next position in the sequence after ensuring that the job in the previous position has finished processing on the same machine.In other words, the completion time of any job on any machine is determined by the processing time on that machine in addition to the completion time of its predecessor on the same machine.Constraint (6) states that jobs are processed in the same order, with one speed on each machine.Constraint (7) ensures that each position on each machine has one job with one speed.Furthermore, constraint (8) guarantees that each job is processed at exactly one speed and has one position on each machine.The makespan of the schedule is defined in constraint (9).Idle times on machines are computed using constraint (10).Constraint (11) calculates TEC in kilowatt hours.The non-negativity and binary constraints of the decision variables are defined in constraints ( 12) and (13), respectively.
It is worth noting that the output of Model 1 is twofold.Firstly, it gives the optimal sequence of processing jobs on machines.Secondly, it assigns the optimal processing speed of each job to each machine.Moreover, the LGP version of Model 1 is introduced in Model 2 below.Model 2 is used to estimate the worst value of each objective., , , 0 n p n p ≥ (17 where n 1 and n 2 (p 1 and p 2 ) are the negative (positive) deviational variables of the two goals.These deviational variables are defined to be non-negative, as stated in constraint (17).Constraints ( 15) and ( 16) specify the two goals and their corresponding target values, i.e., C max * and TEC * , respectively.These target values are calculated from Model 1 above, and they represent the ideal values for both goals.The achievement function in equation ( 14) is a vector of positive deviational variables since both objectives have to be minimized.Furthermore, the achievement function places makespan in the first priority order, followed by TEC in the second place.By considering this order, the maximum (worst) value of TEC is computed due to redundancy in the LGP model.The worst value of makespan is estimated by reversing the previous priority order of the two goals.
The proposed GGP methodology is based on generating a set of GP solutions rather than just one solution.These solutions are obtained by varying the aspiration levels for both goals.This is achieved by a random generation from the sets of intervals in constraints ( 18) and (19) below.On the one hand, the first set of intervals guarantees that the target levels for makespan 1  ( ) t g are increasing in order.On the other hand, the construction of the second set of intervals ensures obtaining values in a decreasing order for the aspiration levels of TEC ) [ ) where min C max is the minimum (ideal) value of the makespan, which is calculated from Model 1 above.Max TEC is the maximum (worst) value of the second objective, and it is computed from Model 2.Moreover, the max length C (length TEC ) is the difference between the maximum and minimum values of makespan (TEC) divided by the number of intervals (n).This number is determined by the decision-maker, and it is equal to the number of aspiration levels used for Model 3 below.
The aspiration levels for the first (second) goal are 2 1 ( ) t t g g .In addition, they are computed from the set of intervals in constraints (18) and (19) above.Consequently, constraints (20) and (21) imply that Model 3 is run n times to obtain a set of n GP solutions.Since there is no guarantee about the efficiency of these solutions, the next step involves applying the HYPssp algorithm [17] to select a subset of size l of these solutions.This subset involves efficient solutions that maximize the hypervolume performance indicator.
The next section illustrates the data settings of the GPFSP and the computational results of the proposed GGP technique.

Data settings and computational results
The literature on the flowshop scheduling problem includes several benchmarks to validate the exact and metaheuristic algorithms.The benchmark of Tailard [30] is the most commonly used.Nevertheless, this study adopts the PFSP benchmark instances of Vallada et al. [31] to assess the performance of the proposed GGP methodology.The reason for this choice is based on the variety of small datasets available in this benchmark, which is consistent with the exactness of the suggested technique.The authors generated 240 large instances and the same number for small-to medium-sized instances.A combination of N jobs and M machines is used to construct their proposed set of benchmark instances.In addition, the number of jobs is between 10 and 60, while the number of machines ranges from five to 20, i.e., N = {10, 20, 30, 40, 50, 60}, and M = {5, 10, 15, 20}. 10 instances are created for each combination of the elements of the two sets.The processing times of jobs are calculated using a uniform distribution.This study considers applying the suggested GGP method to the first five instances of the first combination, i.e., the combination of 10 jobs and five machines.
Table 2 below shows the minimum and maximum values calculated for each objective in each instance.On the one hand, the minimum (ideal) value for each objective is computed from Model 1 above.On the other hand, maximum values are obtained by adopting the idea of goal redundancy in LGP, as illustrated in Model 2. Table 3 below illustrates the parameters used for the GPFSP.The regular parameters of the PFSP are adopted from the study of Vallada et al. [31].They include the number of jobs (N), the number of machines (M), and the data for processing time.The study of Mansouri et al. [25] is used as a reference for setting energy parameters.
The uniform and normal distributions are used to generate the aspiration levels for each objective, i.e., 1 t g and 2 ( 1, 2, , ) 1, 2, …, n) for makespan and TEC, respectively.These distributions are used in this study for the purpose of comparison.On the one hand, the parameters of the uniform distribution are the lower and upper bounds of each interval computed from the set of intervals in constraints (18) and (19).On the other hand, the location parameter of the normal distribution is the average of the lower and upper bounds.The standard deviation is computed as the length of the interval divided by six.The computations of these parameters are based on the study of Moghaddam [32].For each of the five instances, Model 3 is used to obtain two sets of GP solutions.This is achieved using normally and uniformly distributed aspiration levels.Due to obtaining inefficient GP solutions under each instance, the HYPssp algorithm [17] is utilized to obtain subsets of non-dominated solutions that give the maximum value of the hypervolume indicator.Columns 3 and 4 of Table 4 below illustrate the normalized hypervolume values calculated from each subset of efficient solutions using both distributions.Observing these values across the five instances concludes that applying the GGP methodology using normally distributed aspiration levels generates subsets of non-dominated solutions that have a relatively larger volume in the objective space.This means that these solutions are more diverse and welldistributed than those generated by a uniform distribution.In addition, the initial size (n) of the set of GP solutions and the size (l) of the subsets of Pareto optimal solutions are provided in columns 1 and 2, respectively.The sizes of the sets of GP solutions are arbitrarily chosen, and the sizes of the corresponding subsets are selected as n / 2.
Tables A1 and A2 in the appendix show detailed computations for the first instance, i.e., the instance 10 × 5-01.Columns 1 and 2 of Table A1 contain the aspiration levels obtained from the uniform distribution and the corresponding GP solutions, respectively.In addition, column 3 involves the aspiration levels generated from the normal distribution, and column 4 presents the corresponding GP solutions.The subsets of efficient solutions are shown in columns 1 and 2 of Table A2.These efficient solutions are produced using the HYPssp algorithm [17].Due to space limitations, the tables for the rest of the instances are omitted.However, they are available upon request.It is worth noting that the mathematical programming models introduced in this study are implemented using the GAMS software version 24.1.1,and CPLEX is used as the solver.The Java code, written by Bringmann et al. [17], is used to run the HYPssp algorithm.The test instances were run on a laptop with a 1.8 GHz Intel Core i5 8th-generation processor and 8 GB of RAM.

Conclusion
This paper presents a new GP approach to improving the quality of decision-making processes.The proposed GGP method aims to overcome the efficiency issue of GP solutions while preserving the preferences of the decisionmaker.Furthermore, this study introduces an enhancement to the mathematical framework of GP by generating a set of solutions rather than just one solution.This set reflects trade-offs between competing objectives and provides the decision-maker with a complete picture of the feasible space of solutions.In addition, the study highlights the importance of dealing with the efficiency problem of GP solutions.This is achieved by integrating the concept of the HYPssp with the proposed technique.This results in obtaining a set of Pareto optimal solutions that are more diverse and representative of the objective space.To validate the performance of the suggested technique, it has been applied to the GPFSP.This application has practical importance in the field of industrial engineering.Concerning future work, the paper recommends the following points for future research: Firstly, apply the GGP technique to different applications in other fields of study to benefit from its advantages.Secondly, conduct the proposed method using other versions of GP, such as the weighted variant.Finally, study the effect of including uncertain parameters on the mathematical structure of the GGP methodology.

Model 2 .
Deterministic  LGP model for the GPFSP Achievement function

Table 1 .
Indexes, parameters, and variables of the mathematical programming models m Power of machine m (KWh)

Table 2 .
Minimum and maximum values of the makespan and total energy consumption

Table 3 .
Summary of the parameters of the GPFSP

Table 4 .
Results of the hypervolume performance metric

Table A2 .
Subsets of efficient GP solutions computed from the previous table for 10 × 5-01 instance