Solving Order Planning Problem Using a Heuristic Approach: The Case in a Building Material Distributor

: For building material distributors, order planning is a key process as a result of the increase in construction projects’ scale and complexity. In this paper, the integration of simulation modeling and the response surface methodology (RSM) is presented to solve an order planning problem in the construction supply chain. The interactions of various factors are examined to observe their e ﬀ ects on key system measurements, and a combination of factor levels is determined to achieve the optimal performance. RSM is applied to ﬁnd the possible values of the optimal setting for system responses, which consists of three main steps: central composite design (CCD), Box–Behnken design (BBD), and a comparison of both designs. The model is tested with a realistic case study of a building material distributor in Vietnam to demonstrate its e ﬀ ectiveness. Controllable factors (independent variables), which are the review period (T), order quantity (Q), and safety stock (SS), are found to signiﬁcantly a ﬀ ect system responses, which are the total cost (TC) and customer service level (CSL). The results provide the best settings of factor levels that produce the possible minimum TC and maximum CSL. The developed framework could be applied as a useful reference for decision-makers, purchasing managers, and warehouse managers to obtain the most suitable order policy for a robust order planning process.


Introduction
Buildings are becoming increasingly complex because of the ongoing rise in the size and scope of construction projects. Following that, the rising segmentation of the construction industry goes hand-in-hand with the growth of specialist suppliers or contractors, and the diversity of products, designs, and control activities [1]. Relationships in the construction industry are short-term and normally informal/ad-hoc, which focus on the projects, not the business. Construction demand is inherently unstable as the industry is project-based with defined start and endpoints, and a conventional separation between design and construction. Demand is viewed by temporary coalitions as a series of competitively tendered prototypes. Competencies between construction projects vary from time to time and are distinguished from each other. In the supply chain of the construction industry, order planning is the key process, especially for building material distributors. Adequate order planning is quintessential to a more comprehensive integrated supply chain solution, by which businesses can achieve real-time capabilities to act as an intermediary between manufacturing, sales, and customer service concerns, in order to guarantee on-time, effective and reliable processes of order fulfillment and delivery [2]. For building material distributors, the most common problems The contributions of our paper are three-fold. First, the paper proposed a computational framework to solve the order planning problem by a heuristics approach using RSM, which includes the central composite design (CCD), Box-Behnken design (BBD), and variance dispersion graph (VDG), in order to offer the best possible time-efficient and resource-saving solution in terms of key performances. Most exact methods cannot solve the problem within a satisfactory amount of time, especially for large-scale problems. In contrast, most heuristics allows decision-makers to obtain a feasible and timelier solution. For order planning problems in the construction industry, managing the demand, duration, and progress between projects that are supposed to be implemented simultaneously by a distributor is a complicated task. Thus, to facilitate the practical process and available resources, heuristics can provide decision-makers feasible short-term solutions within a reasonable amount of time. Second, the proposed heuristic model is validated by a realistic case study of a distributor in Vietnam to illustrate the effectiveness of the model's result. Third, the managerial implications of the paper could be a beneficial guide for decision-makers, purchasing managers, and warehouse managers to obtain the most suitable order policy for a robust order planning process. In addition, the practical contribution of this study is the comprehensive insight into order planning problems contributed by the case study of the building material distributor in Vietnam.
The paper is divided into five sections. The introduction and relevant studies are discussed in the first section. In the second section, the research procedure and related methodology are included. This paper proposed a computational framework to solve the order planning problem. Additionally, the case study of a distributor is presented in the third section. In the fourth section, the empirical results are shown, including the results of excel spreadsheet simulation, and experimental design of central composite design, and Box-Behnken design. Moreover, discussions, conclusions, and recommendations are given in the last section.

Literature Review
This paper combined simulation modeling and the response surface methodology (RSM) to solve an order planning problem. A simulation and/or computational experiment are used to model a comprehensive supply chain network that corresponds with a number of supply chain operational elements and management levels [23]. Some studies that proposed simulation-based optimization models to explain supply chain systems are discussed in [24,25]. Supply chain dynamics (response) are typically quantified in terms of order and inventory variance ratios. To experiment with the stochastic supply chain model, a simulation modeling approach was adopted. As a result of modifying the controllable variables, the simulation model was run to generate the various supply chain responses, and thus the functional relationship could be fitted and its response surface defined. RSM is a feasible alternative to modeling, on the other hand, which optimizes stochastic, dynamic, and complex systems such as supply chains. This approach is used to find a functional link between the complex responses of the supply chain and significant controllable variables that influence them [26]. Giddings et al. [27] used RSM for optimality analysis of the cost coefficients in mixed integer linear programming in facility location problems. The study implemented the design of experiments and applied least squares regression to determine cost coefficients that significantly affect the optimal total cost surface within setup coefficient ranges. In a problem of finding the suitable capacity for a factory, Shang et al. [28] aimed at integrating simulation, Taguchi techniques, and RSM. Findings of this hybrid approach is a useful reference for businesses to analyze the dynamic relations among various factors, so as to determine the best combination of their levels that optimizes the impact of demand uncertainty on the performance of the supply chain. Buchholz et al. [29] developed a hybrid model of a process-based simulation ProC/B toolset and RSM for the optimization of the process chain models. In this paper, the experiment of the central composite design was applied to find the optimal values of response variables. The "bullwhip" effect, known also as order variance amplification, is a major cause of supply chain deficiencies. Hassanzadeh et al. [30] studied a three-stage simulation with a single retailer, a single wholesaler, and a singer-producer under both centralized and decentralized chains, in order to analyze the causes of the bullwhip effect from two dimensions of order and inventory variance using RSM. As an extension, a multiobjective minimization problem in which the order and inventory variance ratios are separated into two objective functions was developed by Devika et al. [31]. RSM was used in a hybrid evolutionary approach to analyze interactions between the variance ratios and their corresponding effects, in order to optimize supply chain systems. For measuring design parameters on the bullwhip effect and dynamic responses, Tang et al. [32] analyzed its influences on the supply chain's performance using a hybrid Taguchi and dual RSM. This paper offered practical solutions to the supply chain managers and designers in the trade-off between customer service level and inventory holding cost under an uncertain environment.
To provide an overview of recent contributions, several studies that presented relevant problem characteristics are shown in Table 1. Concerning order planning problems, many studies approached exact methods such as NLP [33], MOLP [34,35], and mathematical models [36][37][38][39]. The widely used method is applying inventory policies such as the economic order quantity (EOQ), economic production quantity (EPQ), (S, T), or (R, Q) model [36,39,40]. The order planning problem is actually the flow shop (FS) problem, which is an NP-hard problem [41]. Most exact methods cannot solve large-scale problems within an acceptable amount of time. Therefore, heuristics [33,42] have widely been applied currently, which produce good feasible solutions timely and effectively. In this paper, to optimize an order planning problem, namely to minimize total relevant costs and maximize customer satisfaction level, RSM was used to determine the best combination of the system parameters. Independent variables, which are the review period (T), order quantity (Q), and safety stock (SS), were treated as factors for response variables, which are the total cost (TC) and customer service level (CSL). Analytical approaches have many limitations to simultaneously solve two or three factors [23]. It is difficult to observe the factor interactions and predict their effects on the overall objective when multiple factors are considered at a time. By using RSM, our research aims to determine the best parameter settings of considered factors to find the possible optimum values of the responses. Findings are expected to offer purchasing managers or decision-makers to obtain the optimal combination of the levels of independent factors, which strongly affect the key system performance measures, i.e., total relevant costs and customer satisfaction level.

Response Surface Methodology (RSM)
Response surface methodology (RSM) was used to design and perform an experimental design via statistical and mathematical approaches. RSM was applied to determine the possible optimal values of factors (i.e., independent variables) to minimize and/or maximize the outputs (response variables). RSM was assumed that there is no correlation among independent variables. There are many stages to perform the process of optimization using RSM. The following parts present the simple and polynomial regression model [43][44][45][46]. The simple regression model is presented in Equation (1) as follows.
where Y denotes response variables, X 1 , X 2 , . . . X k denotes independent variables, k is the total of independent variables, β o is the intercept, β 1 , β 2 , . . . β k is the slope, and ε is the random error. Equation (2) shows the quadratic polynomial regression model.
where Y denotes response variables, X denotes independent variables, β o is the intercept, β 1 is the coefficient for the linear terms, β 2 is the coefficient for the square terms, β 3 is the coefficient for the interaction among terms, and ε is the random error.
In the polynomial regression model, CCD and BBD are two types of experimental designs often referred to by the researchers. CCD is an experimental design with two-level. Minitab provides the rotatability value (α) to ensure the design displays the desirable properties. Meanwhile, BBD is the experiment with a three-level design. BBD does not have the rotatability value (α) and embedded factorial design. Points on the diagrams represent the experimental runs of CCD and BBD are shown in Figure 1 as follows [47]. Appl. Sci. 2020, 10, x FOR PEER REVIEW 6 of 22 (a) CCD experiment (b) BBD experiment

A Computational Framework
In this paper, a computational framework was built (Figure 2), which consists of two parts, described as follows: (1) an Excel spreadsheet was used to develop the order planning model, and calculate the response variables, i.e., total cost and customer service level. Then, (2) Minitab was used to find the feasible values of the optimal setting for those responses using RSM, CCD, and BBD. Besides, a VDG was applied to validate these experimental designs. Three independent variables were being considered to have effects on the response variables in this experiment, consisting of the review period (T, month), order quantity (Q, unit), and safety stock (SS, unit). The combination in the total cost (TC, USD) that includes the holding cost, shortage cost, over-storage cost, transportation cost, and customer service level (CSL, %) were chosen as response variables of the order planning model.

Problem Description
The description of the primary building materials of a distributor in Vietnam is shown in

A Computational Framework
In this paper, a computational framework was built (Figure 2), which consists of two parts, described as follows: (1) an Excel spreadsheet was used to develop the order planning model, and calculate the response variables, i.e., total cost and customer service level. Then, (2) Minitab was used to find the feasible values of the optimal setting for those responses using RSM, CCD, and BBD. Besides, a VDG was applied to validate these experimental designs. Three independent variables were being considered to have effects on the response variables in this experiment, consisting of the review period (T, month), order quantity (Q, unit), and safety stock (SS, unit). The combination in the total cost (TC, USD) that includes the holding cost, shortage cost, over-storage cost, transportation cost, and customer service level (CSL, %) were chosen as response variables of the order planning model.

A Computational Framework
In this paper, a computational framework was built (Figure 2), which consists of two parts, described as follows: (1) an Excel spreadsheet was used to develop the order planning model, and calculate the response variables, i.e., total cost and customer service level. Then, (2) Minitab was used to find the feasible values of the optimal setting for those responses using RSM, CCD, and BBD. Besides, a VDG was applied to validate these experimental designs. Three independent variables were being considered to have effects on the response variables in this experiment, consisting of the review period (T, month), order quantity (Q, unit), and safety stock (SS, unit). The combination in the total cost (TC, USD) that includes the holding cost, shortage cost, over-storage cost, transportation cost, and customer service level (CSL, %) were chosen as response variables of the order planning model.

Problem Description
The description of the primary building materials of a distributor in Vietnam is shown in

Problem Description
The description of the primary building materials of a distributor in Vietnam is shown in

A Computational Framework
In this paper, a computational framework was built (Figure 2), which consists of two parts, described as follows: (1) an Excel spreadsheet was used to develop the order planning model, and calculate the response variables, i.e., total cost and customer service level. Then, (2) Minitab was used to find the feasible values of the optimal setting for those responses using RSM, CCD, and BBD. Besides, a VDG was applied to validate these experimental designs. Three independent variables were being considered to have effects on the response variables in this experiment, consisting of the review period (T, month), order quantity (Q, unit), and safety stock (SS, unit). The combination in the total cost (TC, USD) that includes the holding cost, shortage cost, over-storage cost, transportation cost, and customer service level (CSL, %) were chosen as response variables of the order planning model.

Problem Description
The description of the primary building materials of a distributor in Vietnam is shown in    Currently, the company places orders from a supplier every three or four months. Afterward, the transportation freight charge for 20 ft and 40 ft containers that are used for shipments directly to a warehouse. For inventory management, each unit of products (set, pail, or bag) is charged for holding cost per unit per month. Over-storage cost per unit is charged when ending inventory excesses the warehouse capacity. Shortage cost per unit is charged if satisfied demand is less than the actual demand. Key performances are measure periodically, which are the total cost, calculated by the sum of holding inventory, shortage cost, over-storage cost, and transportation cost, as presented in Equation (3); and the customer service level, which is determined by the accumulation of the total demand satisfied, as shown in Equation (4). Subsequently, the company is facing a problem in which should setup different factors (i.e., referred to independent variables in this paper) that affect the performance of the inventory system. Therefore, designing a computational framework to find the feasible values of the optimal setting of levels for independent variables that strongly affect the minimum total cost of holding cost, shortage cost, over-storage cost, transportation cost, and customer service level is needed.
Total cost (TC) can be calculated from the following Equation (3): where C h is holding cost, C s is shortage cost, C o is over-storage cost, and C t is transportation costs.
Appl. Sci. 2020, 10, 8959 The customer service level (CSL) in this model simplifies by the accumulation of the total demand satisfied, and it can be calculated from the following Equation (4):

Numerical Example
In order to demonstrate the model's effectiveness, the authors chose "Spectite CW100", the key product of the company. Compared to other products, Spectite CW100 is always indispensable for construction projects, yet often running out of stock. Additionally, this product was selected based on ABC analysis (i.e., highest sale volume and high order frequency) [48,49]. The ABC classification divides inventory stock-keeping units (SKUs) into three groups based on their annual monetary purchases. These SKUs are high-value items, i.e., the 15-20% of the total inventory that accounts for 75-80% of the total annual inventory value. The product is distributed to customers and given inventory carrying cost, shortage cost, and over-storage cost of each unit. Based on product behaviors and the current order planning process, the interval ranges of the variables were set by the decision-makers as shown in Table 4 and were chosen as potential values in the factorial design for this problem. Total cost (response) Customer service level (response) CSL CSL (90% and 95%) % Note: calculated by the researchers.

Excel Spreadsheet Simulation
First, the authors pivoted the quantity of historical sale data to count monthly order frequency to generate the probability of demand. Next, the authors used the random function in excel to simulate customer demand and get the values response variables (TC and CSL). The experimental design was created using a central composite design and Box-Behnken design to find the possible independent variables settings for solving the order planning problem. The excel spreadsheet simulation template for the experimental design is described in Figure 4.
In this paper, a case study of a distributor is presented to demonstrate the effectiveness of the developed computational framework to solve the current order planning problem in the company.
The model's parameters setup is described as follows.

Statistical Analysis
In the initial step, the statistical analysis using the full factorial design was performed to test the curvature of the response variables (TC and CSL) and independent variables (T, Q, and SS). Then, the data was analyzed in Minitab using stepwise, and the result shows curvature. The analysis of variance of factorial regression and the Pareto chart of the standardized effects are presented in Figures 5 and 6, respectively. In Figure 5, the analysis of variance of factorial regression, the results suggested that the model had significant curvature because the p-value (p-value = 0.000) was less than the confidence level at 95%. The presence of curvature usually indicates that the factor settings were near a feasible response value. Therefore, two response surface designs, central composite design, and Box-Behnken design were considered afterward to find the possible parameter settings for this case. In the step of full factorial design, for the starting center point, T = 3, Q = 3500, and SS = 1000, were considered. Figure   Figure 4. Excel spreadsheet simulation.

Statistical Analysis
In the initial step, the statistical analysis using the full factorial design was performed to test the curvature of the response variables (TC and CSL) and independent variables (T, Q, and SS). Then, the data was analyzed in Minitab using stepwise, and the result shows curvature. The analysis of variance of factorial regression and the Pareto chart of the standardized effects are presented in Figures 5 and 6, respectively.

. Statistical Analysis
In the initial step, the statistical analysis using the full factorial design was performed to test th rvature of the response variables (TC and CSL) and independent variables (T, Q, and SS). Then e data was analyzed in Minitab using stepwise, and the result shows curvature. The analysis o riance of factorial regression and the Pareto chart of the standardized effects are presented i ures 5 and 6, respectively. In Figure 5, the analysis of variance of factorial regression, the results suggested that the mod d significant curvature because the p-value (p-value = 0.000) was less than the confidence level a %. The presence of curvature usually indicates that the factor settings were near a feasible respons

. Central Composite Design (CCD)
In this section, a central composite design was performed to find the effects of independe riables on the response variables. In this paper, the CCD experiment (one continuous factor an o categorical factors) used the 20-run design with two blocks. Minitab provides the rotatabili lue (α) to ensure the design displays the desirable properties. The rotatability value of the CC periment was chosen based on factorial run and factorial point, α = 1.41 was applied in th perimental design. The CCD experiment running setup and randomized design table are presente Figure 7 and Table 5 below. Figure 7 shows all set-up parameters for the experimental CCD in th search. From the results of Table 5, TC and CSL were calculated based on independent variable hich T, Q, and SS. Note: StdOrder (standard order), RunOrder (run order), and PtType (center poi pe). In Figure 5, the analysis of variance of factorial regression, the results suggested that the model had significant curvature because the p-value (p-value = 0.000) was less than the confidence level at 95%. The presence of curvature usually indicates that the factor settings were near a feasible response value. Therefore, two response surface designs, central composite design, and Box-Behnken design were considered afterward to find the possible parameter settings for this case. In the step of full factorial design, for the starting center point, T = 3, Q = 3500, and SS = 1000, were considered. Figure 6 shows the Pareto chart of the standardized effects from the results and factors. It can be seen that their interaction significantly affected the model's results.

Central Composite Design (CCD)
In this section, a central composite design was performed to find the effects of independent variables on the response variables. In this paper, the CCD experiment (one continuous factor and two categorical factors) used the 20-run design with two blocks. Minitab provides the rotatability value (α) to ensure the design displays the desirable properties. The rotatability value of the CCD experiment was chosen based on factorial run and factorial point, α = 1.41 was applied in this experimental design. The CCD experiment running setup and randomized design table are presented in Figure 7 and Table 5 below. Figure 7 shows all set-up parameters for the experimental CCD in the research. From the results of Table 5, TC and CSL were calculated based on independent variables, which T, Q, and SS. Note: StdOrder (standard order), RunOrder (run order), and PtType (center point type). The initial experimental model was built for the response variables. Figure 8 shows the response surface regression of total cost versus T, Q, and SS. In the first running, the results suggest that the terms SS, the square terms SS*SS, and the interaction among Q*SS, Q*T, and SS*T were not significant (p-value > 0.05), and the paper considered the confidence level of the experimental design at 95%. In the next steps, these terms were eliminated from the CCD model. experimental design. The CCD experiment running setup and randomized design table are presented in Figure 7 and Table 5 below. Figure 7 shows all set-up parameters for the experimental CCD in the research. From the results of Table 5, TC and CSL were calculated based on independent variables, which T, Q, and SS. Note: StdOrder (standard order), RunOrder (run order), and PtType (center point type).    The initial experimental model was built for the response variables. Figure 8 shows the response surface regression of total cost versus T, Q, and SS. In the first running, the results suggest that the terms SS, the square terms SS*SS, and the interaction among Q*SS, Q*T, and SS*T were not significant (p-value > 0.05), and the paper considered the confidence level of the experimental design at 95%. In the next steps, these terms were eliminated from the CCD model. After excluded these terms, the experiment was run once more. The results are shown in Figure  9, the reduced model. From the results, term T and the square term Q*Q were highly significant (pvalue < 0.05), while term Q was barely significant. Thus, term Q was kept for the next analysis to see the effects of these independent variables on the objective functions. After excluded these terms, the experiment was run once more. The results are shown in Figure 9, the reduced model. From the results, term T and the square term Q*Q were highly significant (p-value < 0.05), while term Q was barely significant. Thus, term Q was kept for the next analysis to see the effects of these independent variables on the objective functions. For the second objective function, CSL (%), the CCD experiment was performed consequentially. The response surface regression of the CSL versus T, Q, and SS is presented in Figure 10. In the first running, the results display that the square terms SS*SS, and the interaction among Q*SS and Q*T was not significant (p-value > 0.05). Then, these terms were eliminated from the CCD model. After excluded these terms, the CCD model was run again, the reduced model is shown in Figure  11. The results proposed that only the square term Q*Q was not significant because of the p-value = 0.068 was out of the confidence level at 95%. Meanwhile, the terms Q, SS, T, and the interaction of SS*T dramatically affected the CCD model (p-value < 0.05). After reducing the CCD model into the final form, the residual plots for TC and CSL are displayed in Figure 12, including normal probability plot, histogram, versus fits, and versus order. For the second objective function, CSL (%), the CCD experiment was performed consequentially. The response surface regression of the CSL versus T, Q, and SS is presented in Figure 10. In the first running, the results display that the square terms SS*SS, and the interaction among Q*SS and Q*T was not significant (p-value > 0.05). Then, these terms were eliminated from the CCD model. For the second objective function, CSL (%), the CCD experiment was performed consequentially. The response surface regression of the CSL versus T, Q, and SS is presented in Figure 10. In the first running, the results display that the square terms SS*SS, and the interaction among Q*SS and Q*T was not significant (p-value > 0.05). Then, these terms were eliminated from the CCD model. After excluded these terms, the CCD model was run again, the reduced model is shown in Figure  11. The results proposed that only the square term Q*Q was not significant because of the p-value = 0.068 was out of the confidence level at 95%. Meanwhile, the terms Q, SS, T, and the interaction of SS*T dramatically affected the CCD model (p-value < 0.05). After reducing the CCD model into the final form, the residual plots for TC and CSL are displayed in Figure 12, including normal probability plot, histogram, versus fits, and versus order. After excluded these terms, the CCD model was run again, the reduced model is shown in Figure 11. The results proposed that only the square term Q*Q was not significant because of the p-value = 0.068 was out of the confidence level at 95%. Meanwhile, the terms Q, SS, T, and the interaction of SS*T dramatically affected the CCD model (p-value < 0.05). For the second objective function, CSL (%), the CCD experiment was performed consequentially. The response surface regression of the CSL versus T, Q, and SS is presented in Figure 10. In the first running, the results display that the square terms SS*SS, and the interaction among Q*SS and Q*T was not significant (p-value > 0.05). Then, these terms were eliminated from the CCD model. After excluded these terms, the CCD model was run again, the reduced model is shown in Figure  11. The results proposed that only the square term Q*Q was not significant because of the p-value = 0.068 was out of the confidence level at 95%. Meanwhile, the terms Q, SS, T, and the interaction of SS*T dramatically affected the CCD model (p-value < 0.05). After reducing the CCD model into the final form, the residual plots for TC and CSL are displayed in Figure 12, including normal probability plot, histogram, versus fits, and versus order. After reducing the CCD model into the final form, the residual plots for TC and CSL are displayed in Figure 12, including normal probability plot, histogram, versus fits, and versus order. The results show that the assumption of the model had no concern about the violation of the acceptable outliner. In detail, the normal probability plot nearly followed a straight line, then the residuals could assume from a normal distribution. The versus fits were used to verify that the residuals were scattered randomly about zero. From the plots, the non-constant variance of the residuals was fit. The histogram in the CCD model approximately followed a distribution. The versus order shows the residuals in the order of data collection. As can be seen, the versus order did not show any pattern, hence, there was no time correlation in the residuals.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 22 The results show that the assumption of the model had no concern about the violation of the acceptable outliner. In detail, the normal probability plot nearly followed a straight line, then the residuals could assume from a normal distribution. The versus fits were used to verify that the residuals were scattered randomly about zero. From the plots, the non-constant variance of the residuals was fit. The histogram in the CCD model approximately followed a distribution. The versus order shows the residuals in the order of data collection. As can be seen, the versus order did not show any pattern, hence, there was no time correlation in the residuals.  Figure 13 shown the response optimizer of the CCD experiment to find the feasible response variables with CSL at 90% and 95%. From the result, the CSL set to target that was 90% with the possible minimum total cost 76,330 of the possible settings were T = 3, Q = 3428, and SS = 1336. Meanwhile, when the CSL was increased up to 95%, the possible minimum total cost was 77,090 of the optimal settings that were T = 3, Q = 3520, and SS = 1354. Besides, the surface and contour plots for CSL are displayed in Figure 14. According to the contour plot, the CSL located in an area of dark green color was more satisfying. The trade-off existed in the objective function of the model. The more the customer service level is satisfied, the more of the total cost will be paid by the distributor.  Figure 13 shown the response optimizer of the CCD experiment to find the feasible response variables with CSL at 90% and 95%. From the result, the CSL set to target that was 90% with the possible minimum total cost 76,330 of the possible settings were T = 3, Q = 3428, and SS = 1336. Meanwhile, when the CSL was increased up to 95%, the possible minimum total cost was 77,090 of the optimal settings that were T = 3, Q = 3520, and SS = 1354.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 22 The results show that the assumption of the model had no concern about the violation of the acceptable outliner. In detail, the normal probability plot nearly followed a straight line, then the residuals could assume from a normal distribution. The versus fits were used to verify that the residuals were scattered randomly about zero. From the plots, the non-constant variance of the residuals was fit. The histogram in the CCD model approximately followed a distribution. The versus order shows the residuals in the order of data collection. As can be seen, the versus order did not show any pattern, hence, there was no time correlation in the residuals.  Figure 13 shown the response optimizer of the CCD experiment to find the feasible response variables with CSL at 90% and 95%. From the result, the CSL set to target that was 90% with the possible minimum total cost 76,330 of the possible settings were T = 3, Q = 3428, and SS = 1336. Meanwhile, when the CSL was increased up to 95%, the possible minimum total cost was 77,090 of the optimal settings that were T = 3, Q = 3520, and SS = 1354. Besides, the surface and contour plots for CSL are displayed in Figure 14. According to the contour plot, the CSL located in an area of dark green color was more satisfying. The trade-off existed in the objective function of the model. The more the customer service level is satisfied, the more of the total cost will be paid by the distributor. Besides, the surface and contour plots for CSL are displayed in Figure 14. According to the contour plot, the CSL located in an area of dark green color was more satisfying. The trade-off existed in the objective function of the model. The more the customer service level is satisfied, the more of the total cost will be paid by the distributor.

Box-Behnken Design (BBD)
Following the results from the experiment of the central composite design, Box-Behnken design was performed to compare the effects of independent variables on the response variables. In this experiment, BBD with three continuous factors used the total 15-run design with one block. Compared to the central composite design, the BBD experiment did not have an embedded factorial design and extreme points, and the BBD experiment also did not have the rotatability value (α) in the experimental design. The advantage of the BBD model was that the design usually had fewer design points, hence, they were often less expensive to run than the CCD design with the same number of independent variables. The BBD experiment running setup and randomized design table is presented in Figure 15 and Table 6 as follows. Figure 15 shows all set-up parameters for the experimental BBD. From the results of Table 6, TC and CSL were calculated based on considered independent variables using an excel spreadsheet.

Box-Behnken Design (BBD)
Following the results from the experiment of the central composite design, Box-Behnken design was performed to compare the effects of independent variables on the response variables. In this experiment, BBD with three continuous factors used the total 15-run design with one block. Compared to the central composite design, the BBD experiment did not have an embedded factorial design and extreme points, and the BBD experiment also did not have the rotatability value (α) in the experimental design. The advantage of the BBD model was that the design usually had fewer design points, hence, they were often less expensive to run than the CCD design with the same number of independent variables. The BBD experiment running setup and randomized design table is presented in Figure 15 and Table 6 as follows. Figure 15 shows all set-up parameters for the experimental BBD. From the results of Table 6, TC and CSL were calculated based on considered independent variables using an excel spreadsheet.
The initial model was obtained for the response variables (TC and CSL). Figure 16 shows the response surface regression of total cost versus T, Q, and SS. In the first running, the results suggest that the square terms Q*Q, SS*SS (p-values were 0.976 and 0.915 respectively), and the interaction among T*SS and Q*SS (p-values were 0.721 and 0.902, respectively) were not significant (all p-value > 0.05), while the paper considered the confidence level of the BBD at 95%. Then, these terms were eliminated from the BBD model to obtain the following final model.
After excluded insignificant terms, the experimental BBD was run again, the results of the reduced model are shown in Figure 17. The results display that term SS was statistically significant while the terms T, Q, square term T*T, and the interaction of T*Q significantly affected the model (p-value < 0.05). The authors decided to keep the term SS for further analysis to see the effects of these independent variables on the objective functions.

Box-Behnken Design (BBD)
Following the results from the experiment of the central composite design, Box-Behnken design was performed to compare the effects of independent variables on the response variables. In this experiment, BBD with three continuous factors used the total 15-run design with one block. Compared to the central composite design, the BBD experiment did not have an embedded factorial design and extreme points, and the BBD experiment also did not have the rotatability value (α) in the experimental design. The advantage of the BBD model was that the design usually had fewer design points, hence, they were often less expensive to run than the CCD design with the same number of independent variables. The BBD experiment running setup and randomized design table is presented in Figure 15 and Table 6 as follows. Figure 15 shows all set-up parameters for the experimental BBD. From the results of Table 6, TC and CSL were calculated based on considered independent variables using an excel spreadsheet.    The initial model was obtained for the response variables (TC and CSL). Figure 16 shows the response surface regression of total cost versus T, Q, and SS. In the first running, the results suggest that the square terms Q*Q, SS*SS (p-values were 0.976 and 0.915 respectively), and the interaction among T*SS and Q*SS (p-values were 0.721 and 0.902, respectively) were not significant (all p-value > 0.05), while the paper considered the confidence level of the BBD at 95%. Then, these terms were eliminated from the BBD model to obtain the following final model. After excluded insignificant terms, the experimental BBD was run again, the results of the reduced model are shown in Figure 17. The results display that term SS was statistically significant while the terms T, Q, square term T*T, and the interaction of T*Q significantly affected the model (pvalue < 0.05). The authors decided to keep the term SS for further analysis to see the effects of these independent variables on the objective functions. For the second objective function related to CSL (%), the BBD was performed consequentially. The response surface regression of the CSL versus T, Q, and SS is presented in Figure 18. In the first The initial model was obtained for the response variables (TC and CSL). Figure 16 shows the response surface regression of total cost versus T, Q, and SS. In the first running, the results suggest that the square terms Q*Q, SS*SS (p-values were 0.976 and 0.915 respectively), and the interaction among T*SS and Q*SS (p-values were 0.721 and 0.902, respectively) were not significant (all p-value > 0.05), while the paper considered the confidence level of the BBD at 95%. Then, these terms were eliminated from the BBD model to obtain the following final model. After excluded insignificant terms, the experimental BBD was run again, the results of the reduced model are shown in Figure 17. The results display that term SS was statistically significant while the terms T, Q, square term T*T, and the interaction of T*Q significantly affected the model (pvalue < 0.05). The authors decided to keep the term SS for further analysis to see the effects of these independent variables on the objective functions. For the second objective function related to CSL (%), the BBD was performed consequentially. The response surface regression of the CSL versus T, Q, and SS is presented in Figure 18. In the first running, the results display that the term SS (p-value was 0.585) and the square terms Q*Q and SS*SS (p-values were 0.461 for both) were not significant (p-value > 0.05). In this case, the square terms Q*Q and SS*SS were eliminated from the BBD model in priority, the term SS was kept for further analysis purposes. For the second objective function related to CSL (%), the BBD was performed consequentially. The response surface regression of the CSL versus T, Q, and SS is presented in Figure 18. In the first running, the results display that the term SS (p-value was 0.585) and the square terms Q*Q and SS*SS (p-values were 0.461 for both) were not significant (p-value > 0.05). In this case, the square terms Q*Q and SS*SS were eliminated from the BBD model in priority, the term SS was kept for further analysis purposes. Appl. Sci. 2020, 10, x FOR PEER REVIEW 16 of 22 After excluded these terms, the BBD experiment was run again, the reduced model is shown in Figure 19. The results show that only the term Q, SS, and the square term T*T (p-values were 0.054, 0.564, and 0.067, respectively) were slightly significant because the p-value was slightly out of the confidence level at 95%. Otherwise, only the terms T was highly significant to the BBD experiment (p-value < 0.05). After reducing the BBD model into the final form, the residual plots (normal probability plot, histogram, versus fits, and versus order plot) for the response variables (TC and CSL) are displayed in Figure 20. The results show that all residual plots satisfied the normality and constant variance assumptions, which means the simulated data was random and followed normal distribution.  Figure 21 shows the response optimizer of BBD to determine the feasible response variables, i.e., TC and CSL (CSL, 90%, 95%). From the result, the CSL set to target was 90% with the possible minimum total cost 70,870 of the possible settings were T = 3, Q = 3603, and SS = 1250. Otherwise, when the CSL was at 95%, the possible minimum total cost was 85,920 of the optimal settings were T = 3, Q = 3690, and SS = 1250. After excluded these terms, the BBD experiment was run again, the reduced model is shown in Figure 19. The results show that only the term Q, SS, and the square term T*T (p-values were 0.054, 0.564, and 0.067, respectively) were slightly significant because the p-value was slightly out of the confidence level at 95%. Otherwise, only the terms T was highly significant to the BBD experiment (p-value < 0.05). After excluded these terms, the BBD experiment was run again, the reduced model is shown in Figure 19. The results show that only the term Q, SS, and the square term T*T (p-values were 0.054, 0.564, and 0.067, respectively) were slightly significant because the p-value was slightly out of the confidence level at 95%. Otherwise, only the terms T was highly significant to the BBD experiment (p-value < 0.05). After reducing the BBD model into the final form, the residual plots (normal probability plot, histogram, versus fits, and versus order plot) for the response variables (TC and CSL) are displayed in Figure 20. The results show that all residual plots satisfied the normality and constant variance assumptions, which means the simulated data was random and followed normal distribution.  Figure 21 shows the response optimizer of BBD to determine the feasible response variables, i.e., TC and CSL (CSL, 90%, 95%). From the result, the CSL set to target was 90% with the possible minimum total cost 70,870 of the possible settings were T = 3, Q = 3603, and SS = 1250. Otherwise, when the CSL was at 95%, the possible minimum total cost was 85,920 of the optimal settings were T = 3, Q = 3690, and SS = 1250. After reducing the BBD model into the final form, the residual plots (normal probability plot, histogram, versus fits, and versus order plot) for the response variables (TC and CSL) are displayed in Figure 20. The results show that all residual plots satisfied the normality and constant variance assumptions, which means the simulated data was random and followed normal distribution. After excluded these terms, the BBD experiment was run again, the reduced model is shown in Figure 19. The results show that only the term Q, SS, and the square term T*T (p-values were 0.054, 0.564, and 0.067, respectively) were slightly significant because the p-value was slightly out of the confidence level at 95%. Otherwise, only the terms T was highly significant to the BBD experiment (p-value < 0.05). After reducing the BBD model into the final form, the residual plots (normal probability plot, histogram, versus fits, and versus order plot) for the response variables (TC and CSL) are displayed in Figure 20. The results show that all residual plots satisfied the normality and constant variance assumptions, which means the simulated data was random and followed normal distribution.  Figure 21 shows the response optimizer of BBD to determine the feasible response variables, i.e., TC and CSL (CSL, 90%, 95%). From the result, the CSL set to target was 90% with the possible minimum total cost 70,870 of the possible settings were T = 3, Q = 3603, and SS = 1250. Otherwise, when the CSL was at 95%, the possible minimum total cost was 85,920 of the optimal settings were T = 3, Q = 3690, and SS = 1250.  Figure 21 shows the response optimizer of BBD to determine the feasible response variables, i.e., TC and CSL (CSL, 90%, 95%). From the result, the CSL set to target was 90% with the possible minimum total cost 70,870 of the possible settings were T = 3, Q = 3603, and SS = 1250. Otherwise, when the CSL was at 95%, the possible minimum total cost was 85,920 of the optimal settings were T = 3, Q = 3690, and SS = 1250. Moreover, the surface and contour plots for CSL are shown in Figure 22. The results suggest that the possible minimum total cost was in an area of light green (<74,000), the area of SS less than 1000, and the area of Q was less than 3500 according to the contour plot. The trade-off existed in the objective function of the model, which means, the more the customer service level is satisfied, the more of the total cost will be paid by the distributor.

The Comparision of the CCD and BBD Experiment
This section performs the comparisons of the CCD and BBD experiments. The summary of the response optimizer results of CCD and BBD are displayed in Table 7 below. Next, the overlaid variance dispersion graph for CCD and BBD are drawn as in Figure 23. RSM was applied to build the prediction models, hence, prediction variance was significant. For more than two independent variables in the model, the variance dispersion graph (VDG) was one of the useful tools to validate the surface's results. VDG shows the value of minimum, maximum, and average of the scaled prediction variance (SPV) versus the distance of the design point from the center point [50]. Moreover, the surface and contour plots for CSL are shown in Figure 22. The results suggest that the possible minimum total cost was in an area of light green (<74,000), the area of SS less than 1000, and the area of Q was less than 3500 according to the contour plot. The trade-off existed in the objective function of the model, which means, the more the customer service level is satisfied, the more of the total cost will be paid by the distributor. Moreover, the surface and contour plots for CSL are shown in Figure 22. The results suggest that the possible minimum total cost was in an area of light green (<74,000), the area of SS less than 1000, and the area of Q was less than 3500 according to the contour plot. The trade-off existed in the objective function of the model, which means, the more the customer service level is satisfied, the more of the total cost will be paid by the distributor.

The Comparision of the CCD and BBD Experiment
This section performs the comparisons of the CCD and BBD experiments. The summary of the response optimizer results of CCD and BBD are displayed in Table 7 below. Next, the overlaid variance dispersion graph for CCD and BBD are drawn as in Figure 23. RSM was applied to build the prediction models, hence, prediction variance was significant. For more than two independent variables in the model, the variance dispersion graph (VDG) was one of the useful tools to validate the surface's results. VDG shows the value of minimum, maximum, and average of the scaled prediction variance (SPV) versus the distance of the design point from the center point [50].

The Comparision of the CCD and BBD Experiment
This section performs the comparisons of the CCD and BBD experiments. The summary of the response optimizer results of CCD and BBD are displayed in Table 7 below. Next, the overlaid variance dispersion graph for CCD and BBD are drawn as in Figure 23. RSM was applied to build the prediction models, hence, prediction variance was significant. For more than two independent variables in the model, the variance dispersion graph (VDG) was one of the useful tools to validate the surface's results. VDG shows the value of minimum, maximum, and average of the scaled prediction variance (SPV) versus the distance of the design point from the center point [50].

Discussions and Conclusions
This paper proposed a computational framework for the order planning problem in a distributor. In the proposed model, values of response variables were simulated and calculated, and Minitab was used to determine the possible optimal setting for those responses using the central composite design and Box-Behnken design. In the case study, a combination of the review period (T), order quantity (Q), and safety stock (SS) was selected as independent variables of the model. Additionally, the total cost of the holding cost, shortage cost, over-storage cost, and transportation cost, and customer service level was selected as response variables. The objective of the paper was to minimize the total cost by planning on placing the right order quantities, in the right periods, with the possible safety stock in order to reach customer satisfaction around 90-95 percent using RSM.
For solving an order planning problem of Spectite CW100, the statistical result from the factorial design shows that the curvature had a significant effect on the responses, therefore, two response surface designs, central composite design, and Box-Behnken design, were considered afterward to find the possible parameter settings. Based on the results in the response optimizer, the following values for the response surface design were as follows: According to the terms of the variance in Figure 23, the blue lines show the BBD experiment while the green lines display the CCD experiment. Both of these designs almost had significant differences in the total cost results with 2000 difference, which means, the distance of design point from the center point, and the scaled prediction variance (SPV) shows the CCD design had better SPV (smallest values) than the BBD design, therefore, CCD experiment can be considered as the best model in this paper. On the other hand, based on the results of the response optimizer, BBD design gave the smallest total cost. In other words, this paper provided useful insights for reference decision-makers in the related industry. Based on the objective of the industry, the managers or policymakers can consider the possible experimental designs that give the best performance toward the strategy of sustainable development.

Discussions and Conclusions
This paper proposed a computational framework for the order planning problem in a distributor. In the proposed model, values of response variables were simulated and calculated, and Minitab was used to determine the possible optimal setting for those responses using the central composite design and Box-Behnken design. In the case study, a combination of the review period (T), order quantity (Q), and safety stock (SS) was selected as independent variables of the model. Additionally, the total cost of the holding cost, shortage cost, over-storage cost, and transportation cost, and customer service level was selected as response variables. The objective of the paper was to minimize the total cost by planning on placing the right order quantities, in the right periods, with the possible safety stock in order to reach customer satisfaction around 90-95 percent using RSM.
For solving an order planning problem of Spectite CW100, the statistical result from the factorial design shows that the curvature had a significant effect on the responses, therefore, two response surface designs, central composite design, and Box-Behnken design, were considered afterward to find the possible parameter settings. Based on the results in the response optimizer, the following values for the response surface design were as follows:

1.
Central composite design (CCD), customer service level CSL was at 90%, the minimum total cost was 76,330 with the possible optimal settings were the time period (T) = 3, order quantity (Q) = 3428, safety stock (SS) = 1336; (2) customer service level CSL was at 95%, the minimum total cost was 77,090 with the possible optimal settings were T = 3, Q = 3520, and SS = 1354.

2.
Box-Behnken design (BBD), customer service level CSL was at 90%, the minimum total cost was 70,870 with the possible optimal settings were the time period (T) = 3, order quantity (Q) = 3603, and safety stock (SS) = 1250; (2) customer service level CSL was at 95%, the minimum total cost was 85,920 with the possible optimal settings were T = 3, Q = 3690, and SS = 1250.
According to the scaled prediction variance (SPV), the CCD design had better SPV (smallest values) than the BBD design, therefore, the CCD experiment could be considered as the best model in this paper. On the other hand, based on the results of the response optimizer, BBD design gave the smallest total cost. When the CSL increased from 90% to 95%, the total cost also increased in both designs. The proposed computational framework can be applied to other products of the distributor (Spectite WS, Spectite HP600, and Speccoat PE145), or any ordering process in businesses within the supply chain. In summary, our contributions include (1) the paper proposed a computational framework to solve the order planning problem by a heuristics approach using RSM, (2) the model was tested with a case study of a distributor to show the model's effectiveness, and (3) The insights implication of the model's results was to provide a useful reference for purchasing managers or decision-makers to obtain the most suitable order policy for a robust order planning process, which includes responses that were considered in monetary values. Besides, this could be used as a guideline for other related industries. However, there were some limitations to this paper. To get more robust solutions, future studies should consider considerable data on the products and address this problem under stochastic demand, lead time uncertainty, and unstable markets.