Stochastic Chebyshev Goal Programming Mixed Integer Linear Model for Sustainable Global Production Planning

Production planning is a necessary process that directly affects the efficiency of production systems in most industries. The complexity of the current production planning problem depends on increased options in production, uncertainties in demand and production resources. In this study, a stochastic multi-objective mixed-integer optimization model is developed to ensure production efficiency in uncertainty conditions and satisfy the requirements of sustainable development. The efficiency of the production system is ensured through objective functions that optimize backorder quantity, machine uptime and customer satisfaction. The other three objective functions of the proposed model are related to optimization of profits, emissions, and employment changing. The objective functions respectively represent the three elements of sustainable development: economy, environment, and sociality. The proposed model also assures the production manager’s discretion over whether or not to adopt production options such as backorder, overtime, and employment of temporary workers. At the same time, the resource limits of the above options can also be adjusted according to the situation of each production facility via the model’s parameters. The solutions that compromise the above objective functions are determined with the Chebyshev goal programming approach together with the weights of the goals. The model is applied to the multinational production system of a Southeast Asian supplier in the textile industry. The goal programming solution of the model shows an improvement in many aspects compared to this supplier’s manufacturing practices under the same production conditions. Last but not least, the study develops different scenarios based on different random distributions of uncertainty demand and different weights between the objective functions. The analysis and evaluation of these scenarios provide a reference basis for managers to adjust the production system in different situations. Analysis of uncertain demand with more complex random distributions as well as making predictions about the effectiveness of scenarios through the advantages of machine learning can be considered in future studies.


Introduction
The textile industry plays a vital role in many Asian countries' economies especially in emerging nations. It is a diverse sector, which has its own chain that includes converting chemical or natural fibers into consumer goods such as garments or industrial textiles. It will require force to transfer from the existing Industry 3.0 to Industry 4.0 and empower smart production. Sustainable fabrication is significant for the effective use of waste reduction, resources allocation efficiency and related costs. Textiles are severely interlinked with environmental, governance and social issues. In the past, efforts of the manufacturer have been primarily concentrated on enhancing the social attributes of textiles, such as creating suitable working conditions, setting minimum wages, maintaining work safety, Production plans can be developed with forecasted demand using DEA-based approaches, which were proposed by Du et al. [6]. In that study, the production system performance is measured by the CCR efficiency score, which was introduced by Charnes, Cooper and Rhodes [7]. Another heuristic-based approach is recommended to determine the effective inventory level to save inventory-related costs, such as holding cost and backorder penalty [8]. The results of the synthesis paper published in 2019 by Bagshaw stated that different quantitative models are developed for different production conditions and problem characteristics [9]. In other words, a universal model or method for all operational production problems does not exist. This study also demonstrates that linear programming planning being an effective quantitative analysis approach for the improvement of production systems [9,10]. An optimal control formulation, introduced by Devizón et. al., integrates two dynamical systems with typical inventory process parameters such as production rate, inventory level, and costs [11]. Another linear planning model developed to minimize costs in the mining industry is suggested by Ozsan et al. The solution of this model is under different production conditions due to limitations on resources [12]. In research published in 2018, Campo developed an optimized mathematical model for the whole production planning problem with the objective function and constraints that represent the nature of the garment industry. In which, the model includes variables that affect specific production systems of the garment industry such as physical variations of fabrics, skills of new workers, and so on [3]. For the furniture industry, Gremani's mathematical model solves the problem of the trade-off between inventory and production of semi-finished or finished products. The results of the study show that while inventory of semi-finished or chopped materials is more cost-effective, it is advantageous to determine the quantity of finished product inventory [13]. The special requirements for their highly skilled workforce and long-time flexibility are characteristic elements of the automotive industry, according to Sillekens et al. Therefore, the workforce planning has been integrated into the mathematical model that incorporates a heuristic algorithm to make useful decisions for automotive production planning [14]. Another application of linear programming model in production planning is in the food industry. The mathematical model of Munhoz and Morabito is concerned with the acidity parameters of the product that are typical of the fruit juice production process [15].
The review paper by Ali Cheraghalikhani, Farid Khoshalhan, and Hadi Mokhtari has categorized production planning models into 6 groups, as showen in Table 1 [16]. Accordingly, the group 1 and group 2 models solve the problem of production planning under certain conditions. Meanwhile, the studies of the remaining groups have considered the uncertainties when developing models or optimization algorithms. While group 3 and group 4 approach uncertainty solutions using fuzzy theory, the studies in group 5 and group 6 apply the theories and probability distributions. This study also pointed out important issues in the production planning problem where the number of related studies decreased over time. Machine utilization, labor characteristics, and multiple facilities are the issues on the list. In the conclusion of this review, the authors propose potential research directions for the production planning problem. Among them, there are multi-objective stochastic models with more than three important issues considered. Table 1. Production planning model classification.

Single Objective Multiple Objective
Deterministic models Group 1 Group 2 Uncertainty models Fuzzy models Group 3 Group 4 Stochastic models Group 5 Group 6 According to research by Mula et al., uncertainty factors in the production system are listed as demand, environment, system resource, lead time, and yield. The study also synthesizes common approaches to uncertainty including stochastics model, dynamic programming, fuzzy theory, and simulation-based approaches [17].
Uncertainty in production plans often appears as demand uncertainty or fluctuations in the production process, such as production times or material loss. Bakir and Byrne analyzed the difference of the solutions provided by deterministic model and stochastic model, in which the uncertainty factor is market demand. The analytical results show that this difference is dependent on the variance of uncertain demand [18].
The uncertainty of parameters in the linear planning model is also approached in terms of fuzzy numbers. Accordingly, some of the parameters describing the production system limit are recognized by the authors as uncertain and described as fuzzy parameters [19][20][21]. In addition, parameters intended to determine the value of the target functions such as purchasing materials cost or production cost, can also be uncertain. Therefore, the linear programming model in Vasant's study in 2006 also introduced a fuzzy objective function [22]. In a recently published study, the authors used fuzzy numbers to describe the uncertainty of unit production time, inventory time, and packing time in the automotive industry. The results of the linear planning model are assessed to have improved the production system in terms of shortening inventory times as well as increasing economic efficiency [23]. Zimerman's method has been applied in the fuzzy mathematical linear programming model, which minimizes all kinds of costs by Kalaf et al. [24]. Baykasoglu and Gocken converted the fuzzy numbers to a non-linear crisp equivalent before using a Tabu Search algorithm to find an optimized solution [25]. Meanwhile, another choice is a genetic algorithm for solving the crisp model. Moreover, this study also uses a combination of fuzzy AHP and Topis to determine the solution that is reconciled [26,27].
In addition to the fuzzy number approaches, production planning uncertainty is also examined through stochastic models with random distributions and the theory of probability. Orcun presented an integrated non-linear programming model with uncertain demand. The author then used a heuristics approach to evaluate the trade-off between the production system's time and safety stock level [28]. A hybrid programming model is developed through a combination of stochastic programming and linear programming for the refinery's production plan, in which the weighting factor is used to control the uncertainty of the hybrid model [29]. Kazemi Zanjani used a scenario tree to control demand uncertainty in a multi-product production planning study. The scenarios in this stochastic model are defined based on the probability distribution of demand during production planning [30]. With similar approaches, Zhang created scenarios for demand uncertainty based on optimistic, pessimistic, and neutral forecast values and their probabilities [31]. For the probabilities, Kazemi Zanjani uses the mean values of the uncertainty parameters as a method for processing the proposed stochastic model [32].
As mentioned above, production plans not only target financial and economic benefits but also consider objectives such as employment, production time, etc. Consequently, the presence of studies that place more than one objective function in a linear programming model is growing. Komsiyah's research has proposed a multi-objective model with the fuzzy goal programming approach in the furniture industry. However, the objective functions of this model revolve around different cost categories that can be combined into a total cost objective function [33]. Meanwhile, in addition to costs, Tirkolaee's study has pointed out another objective function, which is to maximize customer satisfaction. In that study, the authors applied weighted goal programming to solve the problem of multi-objective production planning with assumptions that did not consider overtime and outsourcing [34]. A study applied in Hong Kong has developed an optimization model with three target functions: Maximizing profit, minimizing costs due to defective products, and maximizing equipment utilization at a certain production facility [35]. A similar study has been developed for multinational production systems in which ensuring stability in employment and optimizing the efficiency of the import/export quota becomes the objective functions of the mathematical model [36].
In recent years, the impact of climate change on social activities has become more and more apparent. In particular, industry and manufacturing are often considered the source that generates more emissions compared to other sectors. Therefore, environmental factors also appear in the optimization models to help managers make sustainable production decisions. These environmental impacts in the manufacturing sector are often considered by researchers in terms of production energy consumption. In 2018, Hahn and Brandenburg proposed an optimization model in which the environmental impact factors are converted into a form of costs [37]. Taking a similar approach, Tsai's model describes the environmental impacts such as carbon taxes and the cost savings from reused energy [38]. The optimization model of Modarres and Izadpanahi published in 2016 has the relationship between the uncertainty level of the parameters and the value of the cost of electrical energy. Goal programming is again chosen as the approach to solve the multiobjective model of production planning [39]. In order to increase the efficiency of planning processes, optimization models and algorithms are increasingly integrated both quantitatively and extensively with project management tools [40]. Research by Valenko and Klanšek presented an integrated approach to optimization and Microsoft Project software that increased the cost-effectiveness of the planning process [41].
In summary, the effectiveness of decisions in the production planning process is supported by the results of numerous studies with different approaches. The approach that uses optimization mathematical models has been showing useful applications in the global production planning problems of many industries. Regarding production options, a few studies consider overtime, backorder or part-time worker production in addition to regular production. For the model's objective functions, most studies focus on economic factors such as profits, costs, and penalties. Meanwhile, studies in recent years consider the environmental impact of production systems as an objective function. In addition, indicators showing the system performance such as maximizing utilization of equipment and customer satisfaction are also considered in the models. Moreover, the uncertainty of the parameters in the model has also shown to have a significant effect on the optimization results by researchers. As the results summarize in Table A1, this study develops a model that integrates the above factors simultaneously.

Sets, Parameters and Decision Variables
In this optimization models, parameters and decision variables are formulated based on the set of product families, manufacturing facilities, the time period and the objective function as Table 2.

Set
Indices Description The set of the product families J j The set of the factories T t The set of the planning periods K K The set of objective functions As shown in Table 3, the model parameters can be classified as five groups. The first group is the parameters related to the selling price of the product along with the costs in production, labor, inventory, backorder penalties, recruitment and layoffs. The second parameter group describes the upper and lower boundaries of inventory, backorder, workforce, and machine capacity. The third group includes parameters related to demand, the standard deviation of demand and its probability under uncertainty. The parameters that provide information about the processing time of product families by the machine and the operator are classified in the fourth group. The fifth group includes the parameters of total worker time, overtime limits and allowable rates of workforce change in production facilities. In addition, such parameters as batch size, production specific electricity use of product families; emission factor of manufacturing facilities; goal value, and weight of objective functions fall into this group. The unit inventory cost to hold a product at the end of each period K VND/unit BOC ij The unit backorder cost for a product at the end of each period K VND/unit HC j The hiring cost for one full-time worker K VND/man The fraction of workforce allowable for variation in each period % α j The fraction of overtime hours in factory j used in each period % BZ i The production batch size of product i Units/batch EF j CO2 Emission factor of factory j tCO2/mwh SE Specific electricity use for production mwh/unit z The inverse distribution function of a standard normal distribution with cumulative probability β The proportion of demand met from stock (service level Type II) % G k The aspiration level of objective function k WE k The weight of objective function k The decision variables of the model assist managers in the optimal quantity of products that can be sold, S it , based on information about demand under uncertain conditions. In Table 4, the above production quantities supplied by sources including regular production, temporary worker production, overtime production, inventory, and backorder are denoted x ijt , y ijt , v ijt , I ijt , B ijt respectively. In addition, the model also makes decisions about the workforce level (W jt ) as well as the amount of recruitment (H jt ) and layoffs (F jt ) in manufacturing facilities in each period. Furthermore, the total overtime (OT jt ), temporary labor use time (TPT jt ), and machine uptime (MPT jt ) are also calculated by the model. Finally, deviation variables (o k and u k ) are used to determine the difference between the goal value and the achieved value of the objective functions obtained by the solution. The maximal deviation, ω, is used for Chebyshev Goal Programming to ensure balance between objective functions. The quantity of product i sold at period t Units In this paper, the authors aim for six objective functions that reflect the efficiency of the production system in both economic aspects and factors of sustainable development such as environment and society. First, this model maximizes the profitability of the production system (VND). In Equation (1), the first term represents the estimated total revenue while the rest shows the costs. Those costs are production costs, labor costs, inventory costs, backorder penalties, and cost of hiring/firing workers are presented in the order. The second object function of this model is to minimize emissions (tCO2) that are estimated based on parameters such as emission factor and specific electricity use for production as Equation (2).
Mathematics 2021, 9, 483 8 of 22 Besides environmental impacts, social influences are also a part of ensuring sustainable development. In Equation (3), this study aims to minimize changes in the workforce level (man). This objective function both ensures production stability through ensuring the skill of workers, while also minimizing employment disturbances in the area where the manufacturing facilities are located.
Backorder is one of the production planning options that this model considers. Although this option can bring some advantages to the production system such as saving inventory space, a large backorder could lead to the risk of customers turning to alternative suppliers. Therefore, as Equation (4), the next objective function is to minimize the amount of system production backorder (unit).
In order to improve the efficiency of investment by increasing the utilization rate of production machines, the model maximizes the total operating time of production machinery (hour) in all manufacturing facilities through Equation (5).
In Equation (6), the model's sixth objective function is to maximize customer satisfaction level (%), which is determined by the ratio between the amount of product the system can supply and the expected value of demand. This higher ratio can improve the popularity of the product in the market, but it can also increase costs, leading to a decrease in net profit.

Constraints
The constraints presented in this section describe the resource limitation and policy of the production system. In Equations (7) and (8), the quantity of products expected to be sold does not exceed the expected value of demand and is above the minimum known demand. The expected values of demand are estimated through known minimum, mostlikely, and forecasted maximum demand, combined with their probabilities. This expected sales quantity is balanced with regular production, overtime production, production by temporary workers, inventory, and backorder, as shown in Equation (9). The binary decision variable ϕ ijt in Equation (10) is used to determine the minimum inventory level in Equation (11). The right-hand side of Equation (11) describes the safety stock of a product if it is regularly produced at a given facility over a certain period of time with service level type II (β). In which, z β is the inverse distribution function of a standard normal distribution with cumulative probability β [42]. Equations (12) and (13) ensure that inventory level and backorder do not exceed the maximum allowed limit.
Mathematics 2021, 9, 483 9 of 22 Equation (14) defines change in the workforce by recruiting new workers and layoffs. In particular, the number of workers at a time cannot be outside the upper and lower borders according to the policies of each manufacturing facility as Equation (15). Moreover, Equations (16) and (17) ensure that the number of newly recruited or fired workers cannot exceed a permitted percentage compared to the workforce level in the preceding period. Equation (18) is a non-linear variation of the equation H jt .F jt = 0, which ensures that manufacturing facilities cannot hire and fire workers in the same period.
The constraints (19) ensure that the total time required for regular production does not exceed the total time of the workforce. Similarly, constraints (20)-(23) represent the limits of the production time by temporary workers and overtime. In Equations (24) and (25), the total time required for the aforementioned production types is also limited by the design capacity of the machinery system at the production facility. Finally, all decision variables in this model are considered as non-negative variables.

Chebyshev Goal Programming (CGP) Application
In order to find compromise solutions for all six objective functions outlined above, this study uses the Chebyshev goal programming approach, also known as Minmax Goal Programming. The difference of this variant compared to other variants of goal programming such as weighted goal programming or lexicographic goal programming is to assist in finding balanced solutions between the objective functions instead of extreme solutions. In the first step of this approach, the goal values of the objective functions G k are set either through solving the optimization problem for each objective function individually or consulting decision makers. Next, the deviation variables in both negative u k and positive o k directions from the goal value of each objective function are added. The result of the above step is to form new constraints of the model that are presented in Equations (26)- (31).
The goal values are then used as a normalized constant according to the percentage normalization method to convert the deviation variables to the same unit. Additionally, each deviation variable of each objective function has a weighting parameter WE k that represents its importance relative to the rest. These deviation variables are controlled by the maximal deviation variable ω as constraints (32)- (37). Finally, the only objective function of the model is to minimize the maximum deviation variable with Equation (38).

Case Study of Textile Industry
In this study, a multiple objective optimization model was developed, which applies for the multinational textile supplier's production planning process. Their production system consists of two manufacturing facilities in Vietnam (VN-1, VN-2) and another in Cambodia (KH). These facilities employ an average of 500 full-time employees to produce five product families including large-sized handbags (LH), handbags (HB), backpacks (BP), wallets (W), and key bags (KB) with a total production in 2019 of approximately one million products as presented in Figure 1. According to production managers, the average required production time of the above product families for temporary and full-time workers is different as described in Figure 2.
normalization method to convert the deviation variables to the same unit. Additionally, each deviation variable of each objective function has a weighting parameter that represents its importance relative to the rest. These deviation variables are controlled by the maximal deviation variable as constraints (32)- (37). Finally, the only objective function of the model is to minimize the maximum deviation variable with Equation (38). (38)

Case Study of Textile Industry
In this study, a multiple objective optimization model was developed, which applies for the multinational textile supplier's production planning process. Their production system consists of two manufacturing facilities in Vietnam (VN-1, VN-2) and another in Cambodia (KH). These facilities employ an average of 500 full-time employees to produce five product families including large-sized handbags (LH), handbags (HB), backpacks (BP), wallets (W), and key bags (KB) with a total production in 2019 of approximately one million products as presented in Figure 1. According to production managers, the average required production time of the above product families for temporary and full-time workers is different as described in Figure 2. In Table 5, manufacturing facilities have the same total number of periodic working hours, but different labor policies such as allowable rate of the workforce changing, upper/lower limits of workforce level, maximum inventory level, overtime limit, and machine capacity.  In Table 5, manufacturing facilities have the same total number of periodic working hours, but different labor policies such as allowable rate of the workforce changing, upper/lower limits of workforce level, maximum inventory level, overtime limit, and machine capacity.
To estimate the production system emissions, this study is based on the textile industry's average electrical energy consumption [43] and the emissions factor of the country in which the products are manufactured [44]. Under uncertain conditions, the expected value of demand for product families is estimated using the three-point estimation technique with the program evaluation and review technique (PERT) distribution. Accordingly, the probabilities for the known minimum, most-likely and forecasted maximum demand are 1/6, 4/6, and 1/6 respectively. The optimization problem is solved through IBM CPLEX software version 12.8 on hardware CPU i7 intel and 16 GB RAM. Table 6 describes the dimensions of the optimization problem. By solving the model with each objective function individually, the extreme solutions and their CPU times form the payoff matrix as Table 7.

Goal Programming Solution
Based on the extreme solutions found above, this study defines the goal values as the first step in applying CGP to the model. Meanwhile, the weighting parameters of the objective functions use the non-negative power of ten to represent the dominant nature or the priority between the objective functions. Table 8 shows that 1% deviation of the goal value of profit is ten times more important than the goal value related to emissions. Similarly, every 1% deviation of the emission-related goal value is more important than a thousand times the goal value relative to machine operation time. In other words, the objective function related to machine operation time is dominated by the emission-related objective function.
As shown in Table 9, the solution obtained from the CGP model has some differences in the value of the objective functions compared to the supplier's current production plan in 2019 which is named CP solution in this article. According to these results, CGP's solution resulted in improvements in profitability, emissions, workforce changing and backorder of 21.02%, 37.16%, 30.77% and 23.44% respectively compared to the CP solution. Although the expected value of profit increases, the CGP model also sacrifices the objective functions of machine operation time and customer satisfaction level.  In other words, the CGP solution focuses on effectiveness rather than efficiency. In Figure 3, the production plan of the HB, W, and KB product families is determined by the model to be lower than the CP production plan and uncertain demand. This decision leads to the consequences of low inventory levels as well as reduced production with temporary workers or overtime thereby reducing total costs. Meanwhile, the solution of the model concentrates resources to maintain the production level of more effectiveness product families such as LHB and PB. The CGP model's production plan of these two product families to chase their uncertain demand. In other words, the CGP solution focuses on effectiveness rather than efficiency. In Figure 3, the production plan of the HB, W, and KB product families is determined by the model to be lower than the CP production plan and uncertain demand. This decision leads to the consequences of low inventory levels as well as reduced production with temporary workers or overtime thereby reducing total costs. Meanwhile, the solution of the model concentrates resources to maintain the production level of more effectiveness product families such as LHB and PB. The CGP model's production plan of these two product families to chase their uncertain demand. As one of the factors affecting production levels, the workforce level policy also has the difference between the current operation and the solution of the model as shown in Figure 4. For CP solutions, the workforce is replenished in the first two quarters to ensure production before being maintained in all manufacturing facilities. Meanwhile, the CGP solution shows that the workforce has different changes in manufacturing facilities. The labor force level increased slightly at the VN-1 facility in the first quarter and increased by 50% at the KH facility in the first two quarters. Especially, at the VN-2 facility, the number of full-time workers is downwardly adjusted in the last two quarters of the year to maintain economic effectiveness. In general, the CGP solution suggests a lower workforce level and more variation over time than the manager's current policy.   As one of the factors affecting production levels, the workforce level policy also has the difference between the current operation and the solution of the model as shown in Figure 4. For CP solutions, the workforce is replenished in the first two quarters to ensure production before being maintained in all manufacturing facilities. Meanwhile, the CGP solution shows that the workforce has different changes in manufacturing facilities. The labor force level increased slightly at the VN-1 facility in the first quarter and increased by 50% at the KH facility in the first two quarters. Especially, at the VN-2 facility, the number of full-time workers is downwardly adjusted in the last two quarters of the year to maintain economic effectiveness. In general, the CGP solution suggests a lower workforce level and more variation over time than the manager's current policy.  According to detailed production plan, shown in Figure 5, CP solution's regular production plan shows variability in both volume and variety of product families at each manufacturing facility. In other words, each manufacturing facility was in charge of three to four product families with volumes that fluctuate significantly over time. For CGP solution, besides lower total production, the role of production facilities is also more clearly defined to provide stability in production. In which, the product families which involves handbag as LHB and HB are manufactured at facilities in Vietnam. Meanwhile, the remaining production lines are set up at a manufacturing facility in Cambodia. This arrangement helps the production system save time and the cost of re-establishing the production lines.  Figures 6 and 7 shows the CP solution's production plan by temporary workers and overtime. It can be observed that these two forms of production occur commonly in all manufacturing facilities. In contrast, because of the selection of an effectiveness level of  According to detailed production plan, shown in Figure 5, CP solution's regular production plan shows variability in both volume and variety of product families at each manufacturing facility. In other words, each manufacturing facility was in charge of three to four product families with volumes that fluctuate significantly over time. For CGP solution, besides lower total production, the role of production facilities is also more clearly defined to provide stability in production. In which, the product families which involves handbag as LHB and HB are manufactured at facilities in Vietnam. Meanwhile, the remaining production lines are set up at a manufacturing facility in Cambodia. This arrangement helps the production system save time and the cost of re-establishing the production lines.  According to detailed production plan, shown in Figure 5, CP solution's regular production plan shows variability in both volume and variety of product families at each manufacturing facility. In other words, each manufacturing facility was in charge of three to four product families with volumes that fluctuate significantly over time. For CGP solution, besides lower total production, the role of production facilities is also more clearly defined to provide stability in production. In which, the product families which involves handbag as LHB and HB are manufactured at facilities in Vietnam. Meanwhile, the remaining production lines are set up at a manufacturing facility in Cambodia. This arrangement helps the production system save time and the cost of re-establishing the production lines.  Figures 6 and 7 shows the CP solution's production plan by temporary workers and overtime. It can be observed that these two forms of production occur commonly in all manufacturing facilities. In contrast, because of the selection of an effectiveness level of    Figures 6 and 7 shows the CP solution's production plan by temporary workers and overtime. It can be observed that these two forms of production occur commonly in all manufacturing facilities. In contrast, because of the selection of an effectiveness level of production, CGP profits do not use the temporary workforce in production and only approximately 25,000 units are produced by overtime. This overtime production is equivalent to 11.63% compared to CP solution. production, CGP profits do not use the temporary workforce in production and only approximately 25,000 units are produced by overtime. This overtime production is equivalent to 11.63% compared to CP solution.  Since production levels are determined more efficiently, the inventory levels of the CGP solution are also lower with the same safety stock policy as shown in Figure 8. In addition, the variety of product families of CP solution also has more disadvantages in warehouse management at manufacturing facilities.  production, CGP profits do not use the temporary workforce in production and only approximately 25,000 units are produced by overtime. This overtime production is equivalent to 11.63% compared to CP solution.  Since production levels are determined more efficiently, the inventory levels of the CGP solution are also lower with the same safety stock policy as shown in Figure 8. In addition, the variety of product families of CP solution also has more disadvantages in warehouse management at manufacturing facilities. Since production levels are determined more efficiently, the inventory levels of the CGP solution are also lower with the same safety stock policy as shown in Figure 8. In addition, the variety of product families of CP solution also has more disadvantages in warehouse management at manufacturing facilities. For backorder, both solutions show the backorder that appeared at facility KH where there was the largest planned output as shown in Figure 9. Accordingly, the backorder of the CGP solution is concentrated in the last quarters of the year and is approximately 76.76% equivalent to the CP solution.   For backorder, both solutions show the backorder that appeared at facility KH where there was the largest planned output as shown in Figure 9. Accordingly, the backorder of the CGP solution is concentrated in the last quarters of the year and is approximately 76.76% equivalent to the CP solution.
For backorder, both solutions show the backorder that appeared at facility KH where there was the largest planned output as shown in Figure 9. Accordingly, the backorder of the CGP solution is concentrated in the last quarters of the year and is approximately 76.76% equivalent to the CP solution.
To summarize, the solution that the optimization model provides with the CGP multi-objective approach shows advantages in terms of economic efficiency, lower backorder, and environmental and social sustainability factors. On the other hands, this solution sacrifices the utilization of production machines as well as the risk of the absence of goods when market demand increases, because of a low and more secured production level. However, the optimization model in this study still ensures decision-makers' right to make decisions about satisfaction and sacrifice objectives through weighting parameters.  To summarize, the solution that the optimization model provides with the CGP multiobjective approach shows advantages in terms of economic efficiency, lower backorder, and environmental and social sustainability factors. On the other hands, this solution sacrifices the utilization of production machines as well as the risk of the absence of goods when market demand increases, because of a low and more secured production level. However, the optimization model in this study still ensures decision-makers' right to make decisions about satisfaction and sacrifice objectives through weighting parameters.

Experiments
In this part, the authors examine the changes of the CGP solution according to two factors: Weights of the objective function and probability distributions of demand. For weighting, according to combinatorics, the number of weight combinations increases with a factorial function of the number of objective functions. Thus, there are seven hundred and twenty weight combinations for the six objective functions. However, examining all combinations is ineffective because combinations with similar first, second, and third priority of objective function do not make significant differences in results. Therefore, this study selected seven weighting combinations that have significant differences in priority among the objective functions. For demand's probability distributions, the expected value of demand is estimated using the three-point estimation method in which PERT distribution and triangular distribution are widely used. The probabilities for the known minimum, most-likely and forecasted maximum demand, which are provide by those distribution, are presented in Table 10. Accordingly, the authors develop fourteen scenarios based on the two mentioned factors as described and denoted in Table 11. In it, the scenario where its solution has been discussed in the previous section is named S-Base.  Table 12 shows the value of the objective functions for scenarios along with the supplier's current operations. It can be seen that there are no significant differences in the results between pairs of scenarios with the same weight combinations such as (S-1, S-8), (S-4, S-10), etc. This implies that the difference of the PERT and Triangular distributions does not clearly affect the optimization results. On the other hand, the optimization results of the model clearly show the influence of the weight combinations, which represent the will of the decision-makers on the production system. These scenarios all show their own strengths and weaknesses compared to the rest. However, the solutions of the scenarios S-4, S-6, S-11, and S-13 have low practical application because of low economic efficiency compared to CP solution. The results of scenarios S-5 and S-12 have the same level of customer satisfaction, change in workforce, and machine operation time as CP solution, but greater profits and emissions. The tradeoffs for profit and environmental impact that only creates an improvement in backorder can be difficult to convince production managers. The S-2 and S-9 scenarios are almost identical to the supplier's current policies, but generate more profits and lower emissions. However, these two scenarios have much higher backorder than the rest. Therefore, when applying them, suppliers need to pay certain attention to the loyalty of customers. Breakthrough scenarios from current operating policy are S-Base, S-1, S-3, S-7, S-8, and S-10 with marked improvement in profitability and sustainable development factors. Among them, the S-3 and S-10 scenarios sacrifice a small profits and environmental impact to achieve less variability in workforce, less backorder, and greater customer satisfaction level.
Finally, the authors normalize the units between the objective functions of the scenarios to scores in the interval [0, 1] in order to analyze the trade-offs between them. The scores presented in Figure 10 show that the S-3 and S-10 scenarios have the greatest number of objective functions whose scores are greater than 0.7. Next is the group of S-Base, S-1, S-7, and S-8 scenarios that have three objective functions with scores between 0.4 and 0.7, besides one objective function that achieves upper extreme scores. The remaining scenarios also have maximum scores at a certain objective function, but the rest are below 0.5.  In conclusion, the optimization model in this study provides managers with solutions that can improve the production planning process. These solutions can be determined under the influence of both subjective and objective factors. Subjective factors represent the will of production managers or decision-makers in choosing the priority between the objective functions. At the same time, market surveys to estimate the objective variability of demand along with its probabilities also increase the model's effectiveness in practice.

Conclusions
In this study, the authors develop a stochastic multi-objective mixed-integer programming model for global sustainable multi-product production planning. The variables of this model assist decision making regarding the quantity of production options including regular production, temporary worker production, overtime production, backorder for each period, and production facilities in many countries. Production resource policies can be adjusted by the managers through model parameters such as facilities' maximum inventory level, inventory service level, workforce level limitation, overtime limitation, backorder limitation, machine capacity, and so on. The most optimistic, pessimistic, and most-likely values of the demand, along with their PERT probabilities are the basis for the model to calculate the expected value of the demand and the safety stock. The objective functions that the proposed model considers are, respectively, maximizing total profit, minimizing production emissions, minimizing change in workforce, minimizing backorder, maximizing machine operation time, and maximizing customer satisfaction. The multi-objective solution in this study is determined by the Chebyshev variant of the In conclusion, the optimization model in this study provides managers with solutions that can improve the production planning process. These solutions can be determined under the influence of both subjective and objective factors. Subjective factors represent the will of production managers or decision-makers in choosing the priority between the objective functions. At the same time, market surveys to estimate the objective variability of demand along with its probabilities also increase the model's effectiveness in practice.

Conclusions
In this study, the authors develop a stochastic multi-objective mixed-integer programming model for global sustainable multi-product production planning. The variables of this model assist decision making regarding the quantity of production options including regular production, temporary worker production, overtime production, backorder for each period, and production facilities in many countries. Production resource policies can be adjusted by the managers through model parameters such as facilities' maximum inventory level, inventory service level, workforce level limitation, overtime limitation, backorder limitation, machine capacity, and so on. The most optimistic, pessimistic, and most-likely values of the demand, along with their PERT probabilities are the basis for the model to calculate the expected value of the demand and the safety stock. The objective functions that the proposed model considers are, respectively, maximizing total profit, minimizing production emissions, minimizing change in workforce, minimizing backorder, maximizing machine operation time, and maximizing customer satisfaction. The multi-objective solution in this study is determined by the Chebyshev variant of the goal programming approach. This variation is able to help goal programming solutions to avoid extreme values of high priority objective functions. Next, scenarios are developed based on the change in the weights of the objective functions and the probability distribution of demand. Analysis of these scenarios helps decision-makers to understand the overall picture of the production system when these two factors fluctuate.
The main contribution of this study was to develop an integrated multiple objectives model as a theoretical reinforcement for the production planning problem. The model considers typical current production system characteristics such as multi-product, multifacility, multi-period, and uncertainty factors. The proposed model also provides assistance for managers in deciding whether or not the production options are applicable. The options mentioned include regular production, temporary worker production, overtime production, and backorder. The applicability of these production options can be adjusted through the parameters to which the model relates to resource constraints. In addition, the model's objective functions are defined to ensure both efficient and sustainable development of the production system. Profit maximization, emission minimization, and employment volatility minimization are the objective functions for which the model ensures a balance between the three pillars of sustainable development: Economic, environmental, and social. At the same time, the objective functions that minimize backorder, maximize machine uptime, and maximize customer satisfaction towards optimizing production system efficiency. In order to achieve solutions that are weighted and balanced at the same time for these objective functions, this study proposes the Chebyshev goal programming approach. Another contribution of this study is the development and evaluation of scenarios that depend on the random distributions of uncertain demand and the weights of the objective functions from the point of view of the production managers. This assessment is a useful reference for managers to adjust their production plans for uncertainties and production changes according to the global situation.
Future studies may consider the complexity of other random distributions in terms of the uncertainties of demand as well as production resources. In addition, the advantages of machine learning are also suitable approaches to predict the effectiveness of scenarios as their number may increase exponentially.