A Multifactory Integrated Production and Distribution Scheduling Problem with Parallel Machines and Immediate Shipments Solved by Improved Whale Optimization Algorithm

This paper discusses the integrated scheduling of production and distribution operations in a multifactory supply chain with maketo-order production system. For the production side of the supply chain, we considered distributed parallel-established factories with identical parallel machines available at each factory. We assumed that the factories could produce all customers’ orders with different production rates and costs. For the distribution side of the supply chain, we considered a limited number of homogeneous vehicles that immediately distribute the finalized orders to the customers. Then, a mixed-integer nonlinear programming model is developed to determine the detailed scheduling of production and distribution that minimizes the total costs of the supply chain including production, distribution, and late delivery costs. To solve the real-world scale problems, we developed a new whale optimization algorithm (WOA). Moreover, we conducted computational experiments by generating several test problems to evaluate the proposed algorithm. Statistical analysis showed that the proposed algorithm has better performance than traditional WOA for different scales of the problem. Moreover, it confirms the capability of the improved whale optimization algorithm (IWOA) to solve the medium-scale instances; however, the results indicate the better performance of genetic algorithm (GA) for the large-scale instances.


Introduction
The increasing competition in different marketplaces and the necessity of appropriate responses to the rapid changes in market requirements have made companies shift from a single production plant to multiple factories adjacent to their customers in different geographical areas.Such strategies reduce the transportation costs of the supply chain and allow companies to utilize economic privileges such as availability and low expense of raw materials, expert labors, energy, and tax allowances.The advent of new joint venture companies in cooperation with global brands in developing countries is an indication of the fact.
Multifactory supply chains are classified in different categories such as ownership structure of factories, interaction mechanism between factories, and shop configuration.In some multifactory supply chains, all plants belong to the same company, while in some others, each factory belongs to different companies.For the first type, all factories work together to fulfill the objectives of the main company while in the second category, factories cooperate to maximize self-benefits.Different mechanisms of interactions in the factories throughout the supply chain may lead to parallel, series, or network structures.Parallel factories can produce various types of products.In series structures, each production plant produces semifinished products for the next factory.The combination of parallel and series structures makes the network one [1].Shop configuration in each factory can be considered the same as basic scheduling problem including single machine, parallel machine, flow shop, job shop, and open shop [2].
Numerous applications of the multifactory supply chains in the last two decades attracted the attention of researchers on multifactory scheduling problems in different industries such as electronics, food and chemical process industries, automotive companies, and electric power generating industries.Furthermore, integrating distribution and production scheduling operations as two costly functions of the supply chains made another subject of interest in applied operations research named integrated production and distribution scheduling problem (IPDSP).Independent optimization of the two following operations, as a traditional solution, would not necessarily lead to global optimality of the whole network.Thus, it is critical to consider interactions between production and distribution operations in the supply chain [3].
Besides, IPDSP in single factories was widely studied in the last decade.The researchers and practitioners extensively consider multifactory scheduling problem with different shop configurations and objective functions as maximizing total tardiness of the orders.Nevertheless, there is a limited amount of literature for the multifactory integrated production and distribution scheduling problem (MIPDSP) to minimize the total costs of the supply chain.
This paper is aimed at integrating the production and immediate distribution scheduling of various orders at parallel factories with identical parallel machines.The best allocation of the facilities to the orders along with the best scheduling of production and distribution pursue minimizing the operational costs including production, transportation, and late delivery.Therefore, the contributions of this paper are as follows: (i) Integrated mathematical formulation of production and immediate distribution scheduling in a parallel multifactory supply chain with a limited number of vehicles (ii) Considering trade-offs between late delivery, production, and transportation costs in make-to-order production systems to make decisions on the detailed scheduling of production and distribution (iii) Considering different orders with different amounts and due dates for any of the customers compatible with their dynamic lot sizing (iv) Considering the limited and equal fleet sizes at any of the factories with homogeneous vehicles (v) Improving newly proposed WOA to solve MIPDSP (vi) Introducing a new procedure used in the metaheuristic algorithms to ensure the feasibility of the solutions.
The remainder of this paper is divided into the following sections.Section (2) consists of a brief literature review.Section (3) defines the problem including assumptions and the mathematical formulation.Section (4) presents the solution approach consisting of a new proposed procedure used with metaheuristic algorithms and improved whale optimization algorithm (IWOA).Section (5) presents computational experiments' results.Conclusions and the trend of further researches are given in Section (6).

Literature Review
The existing papers regarding the scheduling of production and distribution operations can be grouped into three different categories: (i) IPDSP for single-factory supply chains, (ii) production scheduling problem in multifactory supply chains, and (iii) production and distribution scheduling problem in the multifactory environment.In this section, we review the recent papers proposed for each category in brief.
Chen [4] proposed the latest review article on IPDSP.They presented a unified model representation scheme and classified the existing models by the shipping and delivery methods.The IPDSP in the multifactory environment had been declared as a research gap by that time.In a more recent work on deterministic approaches of IPDSP, Guo et al. [5] presented a MINLP model for integrated production and distribution scheduling with parallel machines in a make-toorder supply chain.They developed a harmony search algorithm and a heuristic to solve the problem.Sawik [6] considered two inconsistent objectives including costs and service level in a three-echelon supply chain and addressed a biobjective stochastic MIP model for joint selection of suppliers and scheduling of production and distribution for some numerical examples using the CPLEX solver.Cheng et al. [7] considered integrated scheduling of production and distribution problem to minimize the total cost of production and distribution.In the production side, there are capacitated batch-processing machines.In the distribution side, there is a third party logistic (3PL) provider to deliver the produce batches.They proposed a metaheuristic algorithm based on colony optimization method to solve the production scheduling and a heuristic algorithm for the distribution part to minimize total costs of scheduling for production and distribution operations in a single factory.Behnamian and Ghomi [2] discussed the significance of multifactory supply chains and for the first time provided a comprehensive review of the present papers.
A limited number of papers are available for the multifactory scheduling problem with parallel machines.Cicirello and Smith [8] presented a dynamic multifactory problem with parallel machines and applied wasp-like agents for distributed coordination of agents with parallel machines to decide on the acceptance of an arriving job on a machine queue.Hooker [9] combined MILP and constraint programming to solve the multifactory scheduling problem.They allocated the jobs to factories using MILP and scheduled on facilities using constraint programming (CP) and the two parts of the problem linked using logic-based Benders decomposition.Kerkhove and Vanhoucke [10] studied a real problem for the production 2 Complexity scheduling of knitted fabrics in a multifactory environment with unrelated single machines in each factory.They presented a hybrid metaheuristic, which combines simulated annealing and genetic algorithm to solve the large-scale problems.Behnamian and Ghomi [11] considered the multifactory problem with several parallel identical machines in each factory.The objective function for the problem was to minimize the maximum completion time of the jobs available at the beginning of the time horizon.They presented a MILP model for the problem as well as a heuristic and genetic algorithm for the solution.The presented that the modeling approach has been severely criticized by Yazdani et al. [12].They analyzed and showed that the proposed model and related algorithm are suffering from serious shortcomings.They introduced three new models for minimizing makespan and total completion time in this paper.Furthermore, three efficient metaheuristics are put forward based on an artificial bee colony.A limited literature is available for the integration of production and transportation scheduling with the serial multifactory structure.As the latest research, Karimi and Davoudpour [1] proposed a branch and bound algorithm for the scheduling of production and transportation problem in a serial multifactory production system.They also consider the same problem with batch delivery [13] so that the processed jobs should wait until all the jobs of the same batch are complete.Moreover, they developed their model and considered waiting time to start the process at the receiving factory [14].A time-indexed formulation with an LP relaxation of binary variables is proposed; they considered deadlines for the same problem and suggested a new approach in the imperialistic competitive algorithm (ICA) to solve the large-scale problems.Terrazas-Moreno and Grossmann [15] considered a multiproduct integrated scheduling of production and distribution problem for the petrochemical industry.The partial production of demand is allowed in different complexes.
To the best of our knowledge, the single paper for the integrated scheduling of production and distribution problem in a multifactory environment with a make-to-order production system is proposed by Sun et al. [16].They considered MIPDSP for a particular supply chain with two different transportation modalities: inland and maritime transportation.They also proposed a two-level genetic algorithm to solve the problem, one level for the assignment part of the problem and the other one for scheduling of the assigned jobs.
The MIPDSP studied in this research is strongly different from previous ones in supply chain structure and assumptions.The most appropriate production system for such supply chains is the make-to-order production system with immediate shipments, without any temporary warehouses and intermediate inventory.The real-world instances are project-based procurements in construction projects, automotive industry, knitting and clothing industry, dairy products, etc.

Problem Description
In this paper, we consider a supply chain with a set of parallel factories F, F = 1, … , f in different geographical locations owned by a holding company.The factories are controlled through a central planning system.There are m α , α ∈ F parallel identical machines in each factory.Factories' production speeds are different due to technological differences.There are a set of geographically dispersed customers L, L = 1, … , l .Each customer issues several orders with different sizes.The production cost and time of the orders are different at any of the factories.There are a limited number of capacitated vehicles to deliver the orders immediately after production.The order sizes are smaller than the vehicles' capacity.The orders' delivery cost is proportional to the type of products and traveling distance.An order can wait in the factory until its transportation starts without any holding cost.The vehicles deliver the orders to the customers on or later than their due dates, and the company incurs a timedependent penalty for the late orders.Since all the factories belong to the same company, the most appropriate optimization objective is the minimization of operational costs including production, transportation, and late delivery costs.Figure 1 shows the supply chain configuration.

Assumptions
(i) All parameters and required data are deterministic and available when scheduling is undertaken.
(ii) There are orders from different customers available at time zero.
(iii) The order sizes are smaller than the capacity of the vehicles.
(iv) An equal number of homogeneous vehicles are available at each factory.
(v) All factories can produce all the customers' orders.
(vi) Each order can be processed at most on one machine, and any machine can process at most one order at any time.
(vii) Factories can be idle, but machines cannot be idle when there are unprocessed orders in the plant.3 Complexity (viii) Different orders of a customer can be processed in various plants.
(ix) Preemption is not allowed for production and delivery operations.
(x) The vehicles cannot be idle when there are undelivered orders in a factory.
3.2.Mathematical Formulation.This section presents the mathematical formulation for MIPDSP which is the first step to understanding the characteristic of these problems.Maximum of production completion time and vehicle available time to transport order r of customer k when it is transported immediately after order i of customer j 4 Complexity production and their delivery times are the same.The number of vehicles available at each factory is equal.Consider that f 1 processes order i and f 2 processes order j then changing the allocation of factories will not result in the worst objective function value.Therefore, the factories have the same priority in assigning the orders.Furthermore, if we allocate factory f 1 to process both of the orders and schedule order i before j, in this case, changing the sequence of the orders will not result in a worse objective function value.Hence, it is concluded that the equal number of vehicles has no influence on the sequencing and scheduling of orders to the machines.Considering the lemmas, the optimal sequencing solution of the problem with the limited and equal number of homogeneous vehicles available at each factory without the loss of generality is the same as the optimal sequencing solution with single homogeneous vehicles available at each plant.For the simplification of the proposed model, we assume that there are single vehicles available at each factory and the vehicles at factories are homogeneous.The simplified MINLP model is as follows: The objective function (i) minimizes the total production, transportation, and late delivery costs.The unit production costs for any of the orders are different in factories.For making the problem more realistic, it is assumed that transportation and late delivery costs of the product unit depend on the type of goods.The distribution cost of orders depends on the type of products and the transportation cost per unit of product and distance.The late delivery cost also depends on product type and delivery delay time.The first statement of the objective function represents the total production cost of the orders.The second statement defines total distribution cost, and the last statement shows the total late delivery cost.
Since each order should be processed in one factory and on one machine, Constraints (2) are established to assure it.Constraints (3) ensure that any of the orders is a successor to another order or the first scheduled order on machines, and Constraints (4) ensure each order is the predecessor of another order or the last scheduled order on machines.Nonlinear constraints (5) establish interconnections between machine allocation and orders scheduling subproblems.According to this constraint, if two different machines process two distinct orders, they should not be a successor or a predecessor to each other.Constraints (6) guarantee the selection of at most one order as first order on a machine, and Constraints (7) guarantee the same for the last order.Constraints (8) allocate a machine to the first scheduled jobs on the machines, and constraints (9) do the same for the last scheduled job on machines.Constraints (10) state that an order cannot be both predecessor and successor on a machine at the same time.Constraints (11) control the first scheduled order on a machine not to be proceeded by the last scheduled order, and Constraints (12) control the last order of a production scheduling on a machine not to be proceeded by the first scheduled order on the machine.Constraints (13) determine the production completion of orders, which is equal to the sum of the production start time of the order and its processing duration.Constraints ( 14) and ( 15) set the first order's production to start at the time zero.Constraints ( 16) and (17) ensure that the difference between production start times of two consecutive orders on the same machine is equal to the processing time of the successor.Constraints (18) assure that any order should be a successor or the first order delivered from a factory, and constraints (19) guarantee that any order should be a predecessor or the last order delivered from a factory.Constraints (20) ensure that any two orders transported from a plant not to be a predecessor and a successor to each other at the same time.Constraints (21) control the first scheduled order for transportation in a factory not to be proceeded by the last scheduled order, and Constraints (22) control the last order of shipping schedule in a plant not to be proceeded by the first order.Constraints (23) guarantee that there is at most one order at the beginning of transportation scheduling at each factory and, Constraints (24) guarantee at most one order at the end of transportation scheduling.Constraints (25) and (26) ensure that the orders selected as first and last order for transportation in a factory are in the transportation scheduling of the factory.Constraints ( 27) establish the interconnections between assignment and transportation scheduling subproblems.Constraints (28) calculate the delivery time of the orders that is the sum of transport start and travel times.Constraints (29) determine the traveling time between factory and customer proportional to the distance and vehicles' mean speed.Constraints (30) guarantee that the orders' delivery times are later than their due dates.Constraints (31) and (32) ensure the immediate transportation of produced orders.Constraints (33) determine feasible transportation start time.Transportation of an order should not start earlier than its production completion time and vehicle available time which is already transporting the previous order.Thus, these constraints determine the maximum of production completion time and the vehicle available time for the consecutively transported

Solution Approach
Since the assignment and parallel machine scheduling problems are both NP-hard problems, without loss of generality, the multifactory integrated scheduling of production and distribution problem is NP-hard as well [17,18].As a result, we propose a simulation-based intelligent procedure linked with metaheuristic algorithms to solve the large-scale problems.
It is worth mentioning that every metaheuristic algorithm can be linked to the proposed procedure.In this research, recently proposed WOA is utilized to solve the scheduling problems for the first time.The helix-shape movement inspired from social behavior of humpback whale makes WOA different with existing algorithms.This mechanism provides a superior exploitation capability to the algorithm in solving the complex problems.Furthermore, we improved the whale optimization algorithm which is called IWOA and evaluated its performance within the proposed procedure.We also evaluated the validity of the presented model and procedure by solving several random test problems in different sizes.

Solution Representation.
The proposed problem includes several binary variables like the assignment of orders to factories and assignment of orders to machines.In this research, we use the random key method to represent the solutions.Each solution is represented by a linear chromosome including a certain number of genes with uniform random values in [0,1].The length of the chromosomes is variable based on the number of customers, number of orders, number of factories, and number of machines.The structure of chromosomes is illustrated in Figure 2.

Simulation-Based Intelligent
Procedure.Constraints satisfaction is a critical issue in metaheuristic algorithms for creating and updating the populations.MIPDSP is a complex problem that includes several interrelated subproblems and constraints that need to be feasible during the execution of metaheuristic algorithms.For example, a vehicle should ship an order from a factory that is produced in, an order's production should be scheduled in a plant that is assigned to, and the shipment of an order should be planned in the available time of vehicles.The penalty function is a traditional approach to deal with these kinds of constraints.Using penalty functions for the violating constraints does not necessarily lead to a final feasible solution.Thus, we proposed a four-step intelligent procedure based on random simulation of operational situation and metaheuristic algorithms.In the proposed procedure, simulation of operations guarantees always feasible solutions and metaheuristic algorithms ensure intelligent search within the feasible solutions.Figure 3 illustrates the steps of the procedure.
Step 1.Initial random population The proposed procedure is initiated with a random population in which the solutions are represented in random key alphabet.Each random solution in initial population is denoted by u i, j, α, γ ∈ 0, 1 .
Step 2. Generate priority vectors At this step, the procedure decodes the random variable u i, j, α, γ to a feasible solution x i, j, α, γ ∈ 0, 1 by determining the priority of orders to assign factories and machines.Algorithm 1 represents the pseudocode of decoding mechanism.
Step 3. Operation simulation At this step, all the operational conditions including order assignment to factories, production scheduling on machines, and transportation scheduling are simulated according to priority vectors determined at Step 2. Intermediate variables like production completion times, delivery times, and delays are calculated accordingly.
Step 4. Fitness evaluation The objective function is calculated at this step based on the main and intermediate variables.Moreover, the best solution within the population is determined if the stopping criterion was not satisfied.
Step 5. Evolutionary processes At this step, evolutionary processes are implemented on the population to generate a new population to continue the procedure in Step 2.
To clarify the decoding mechanism, we proposed an illustrative example with two customers and two orders for each.We used Algorithm 1 to generate a random feasible solution to produce the orders at two factories with two parallel machines i, j, α, γ = 1, 2 .Consider the initial random solution which is represented in Figure 4.The algorithm reshapes the chromosome to u i, j, α, γ which is illustrated in Table 2.According to the algorithm, when i = 1 the sum should be 3 and when i = 2 the sum should be 3.6; on the other hand, when j = 1 the sum should be 2.7 and when j = 2 the sum should be 3.6.Thus, the algorithm adjusts the priority vectors I = 2, 1 and J = 2, 1 Then, the algorithm assigns random machines at random factories to the orders with the sequence (1) second order of customer two, (2) second order of the customer one, (3) then the first order of the customer two, and (4) first order of the customer one.We also used the same concept to determine the production and distribution sequencing of the orders.7 Complexity based algorithms [19].This algorithm mimics the social behavior of whales and their bubble-net hunting strategy.This algorithm assumes that the current best position is the target pray or is close to the optimal solution, so the other search agents will try to update their position towards the best search agent.This algorithm includes a new exploitation mechanism known as "bubble-net attacking method" and exploration mechanism called "search for prey."The algorithm starts with a set of random solutions.At each iteration, search agents (whales) update their position towards either the best search agent or a randomly chosen agent.Search agents update their position towards the best position using the following equations: where t indicates the iteration index, A and D are coefficient vectors, X * is best position vector, and A • D is an elementby-element multiplication of two vectors.While updating the position of search agents, the algorithm makes a random selection between two approaches designed for the exploitation phase: (i) shrinking encircling mechanism and (ii) spiral Step 1: Generate priority vectors Step 2: Operation simulation Step 3: Fitness evaluation (i) Fitness evaluation of solutions (ii) Determining the best solution (iii) Evaluating the stopping criterion Step 4: Evolutionary processes Step 0: Initial random population (i) Orders' priority vector (ii) Customers' priority vector Input data i, j, α, γ: indexes and number of orders, customers, factories, machines; npop: number of population; sort: sort descending u n i, j, α, γ : random variable ∈ 0, 1 procedure i, sum j, α, γ = ∑ j ∑ α ∑ γ u n i, j, α, γ ; I → i of sorted i, sum j, α, γ by sum j, α, γ as the priority vector of orders; j, sum i, α, γ = ∑ i ∑ α ∑ γ u n i, j, α, γ ;; j → j of sorted j, sum α, γ by sum α, γ as the priority vector of customers; for i ∈ I for j ∈ J select a random factory α to be assigned to order i of customer j select a random machineγ at factory α α, γ to be assigned to order i of customer j i, j set x n i, j, α, γ = 1; if order i, j is proceeding the order i′, j′ on α, γ then:  8 Complexity updating position.The following equations are designed to perform a random search by updating all search agents towards a random solution in the exploration phase.
While conducting computational experiments using basic WOA, we observed some shortcomings in exploring solutions.Early convergence of WOA resulted in weak solutions compared with the other algorithms like GA and PSO.It seemed that the functions (41) and (42) designed for the random search of the algorithm are not able to ensure a good exploration of the solution region.To resolve the problem, we designed two perturbation mechanisms for both exploration and exploitation phases of the algorithm without disturbing the rationale behind the basic algorithm as follows: (a) On the exploitation phase, while using the shrinking mechanism, we considered adding a random vector to the search agents updated towards the best solution.We inspired this mechanism from the swarm-based firefly algorithm to increase exploration capability of WOA and diversification of solutions.We used (44) instead of (40) while updating the search agents using the shrinking mechanism in the improved algorithm.
(b) On the exploration phase, besides using a random search agent for updating the current population, we considered a partial random movement considering the direction of the best search agent.The proposed improvement mechanism is designed to provide the algorithm with maximum exploration while updating the population.We used the following equations instead of (42) in our algorithm.
where ϵ t in the third term of both equations is a vector of random numbers drawn from Gaussian or uniform distribution or any random vector according to the domain of the decision variables, and α is a constant number in 0, 1 and D d t is the vector from the current population towards the best search agent.
For detailed information on this algorithm and the equations for updating coefficients A and D , you can refer to the recently published paper [19].Algorithm 2 proposes the pseudocode of the IWOA using P 1 .

Solution
Algorithm.The integrated nature of the problem with several interrelated subproblems creates a complex structure to MIPDSP, which results in too many constraint groups and decision variables.To confront this complexity, a new algorithm is proposed to solve the large-scale problems with respect to the proposed procedure.The steps of the algorithm are presented in Figure 5.

Computational Experiments
This section represents the results of computational experiments designed to examine the validity of the proposed solution approach using IWOA and other metaheuristic algorithms.We performed statistical analysis using oneway analysis of variance (ANOVA) technique to compare the performance of IWOA with WOA as well as different versions of two prevalent metaheuristic algorithms, GA and PSO.

Genetic Algorithm (GA).
GA is an intelligent search algorithm that seeks optimal solutions for the NP-hard optimization problems.Furthermore, a standard and modified genetic algorithm has been used extensively to solve scheduling and routing problems [20][21][22][23].Thus, standard GA with uniform crossover and mutation operators as well as a new version of GA with multiparent crossover operator (GA-MPC) [24] are utilized as measures to evaluate the efficiency of the proposed IWOA.

Particle Swarm Optimization (PSO).
Since PSO is a swarm-based intelligent algorithm, it can be considered as a useful measure to evaluate the effectiveness of the proposed IWOA.Moreover, PSO is applied to solve several discrete optimization problems like machine scheduling, routing, and location-allocation problems [25][26][27].Thus, canonical PSO and standard PSO 2011 called SPSO [28] are used within the proposed simulation-based algorithm (Figure 2) to evaluate IWOA for randomly generated instances.

Experiment Setup.
Several experiments in different sizes including (i) small-scale, (ii) medium-scale, and (iii) largescale instances are conducted to cover all the configurations of the problem.For our experiments, metaheuristic algorithms including GA, GA-MPC, PSO, SPSO, WOA, and IWOA were coded in MATLAB (R2015b).Moreover, the linearized MINLP model was coded in the General Algebraic Modeling System (GAMS).Then, we generated several random instances for each size of the problem.Each experiment in this paper includes a random instance which is solved 30 times using each metaheuristic algorithm.Moreover, small and medium-scale instances are solved once with GAMS, CPLEX solver 12.2, restricting the run time to 500 seconds.
The experiments were performed on a Notebook PC Intel Core i7 (2.93 GH and 4 GB of RAM).We applied a design of experiments to tune the parameters of IWOA, PSO, GA, and GA-MPC which are represented in Table 3.Furthermore, we tuned the parameters of SPSO based on [28].
To evaluate the performance of the metaheuristic algorithms, we define three metrics: (i) percentage of deviation Input data input l, b, α Initialize for n = 1 to npop; generate u n i, j, f , m ∈ 0, 1 for all i, j, f , m; decode the continuous vector u n to binary variable X n ∈ 0, 1 using P 1 ; evaluate the fitness of the initial population; end; save the best solution → X * ; Main steps while t≤ maximume number of iterations; update (A, C); update p as random variable ∈ 0, 1 ; update the current search agent u n t using improved shrinking encircling mechanism; select a random search agent from u rand t ; update the current search agent using improved equations update the current search agent u n t using spiral updating mechanism; end if 3 decode the continuous vector u n t + 1 ∈ 0 1 to the binary variable X n t + 1 ∈ 0, 1 using P 1 ; implement the production and delivery scheduling considering the constraints; evaluate the fitness of the X n t + 1 ; end; update the X * if better; t = t + 1; end while; return X * ; Algorithm 2: Pseudocode for Improved WOA.
PD a is the percentage of deviation of algorithm a, Fit a is the fitness of algorithm a, and Fit * is the best solution ever found for a random instance.For small-size instances, Fit * is the optimal solution of the GAMS CPLEX solver, and for medium-and large-scale problems, it is the best solution ever found by metaheuristic algorithms or the GAMS solver.Since each experiment includes 30 observations, we used the mean value of PDs (MPD) to evaluate the metaheuristic algorithms' performances.Equation ( 47  where PD a i is the percentage of deviation for observation i of algorithm a.
In addition to graphical analysis of the results, one-way ANOVA test is conducted on the observations' MPDs for the statistical analysis of the results.

Results and Discussions.
To evaluate the complexity of the problem, a single experiment is designed to solve a medium-scale random instance with 416 decision variables (2/4/2/2) using the GAMS CPLEX solver relaxing runtime constraint.The best feasible solution obtained after a week runtime was just 9% better than the best solution found by metaheuristic algorithms in less than a minute.The results reveal the weakness of exact algorithms implemented by conventional software and the necessity of the computational methods to solve medium-and large-scale problems.

5.4.1.
Comparison of Small-Scale Instances.For small-scale problems, 30 random numerical experiments are conducted considering all of the available configurations of the problem.15 Complexity MPDs of these experiments are calculated using the optimal solutions obtained from the GAMS CPLEX solver in less than 500 seconds.The reasonable optimal solutions for the small-scale problems confirm the validity of the proposed MINLP model.Table 4 represents the comparative performance of different algorithms.Zeros in the MPD column of this table indicate the conformity of metaheuristic algorithms' solutions with the optimal solutions and insist the validity of metaheuristic algorithms as well.For small instances, while MPD of an instance is not zero, it means that there is at least one run among multiple runs of the algorithm that the solution of the metaheuristic algorithm was not equal to the optimal solution.Moreover, there was no experiment for IWOA with worse performance than other metaheuristic algorithms while there are two experiments for basic WOA that their MPDs are worse than other algorithms.Thus, we can conclude that the competence of IWOA is at least the same as GA, GA-MPC, PSO, and SPSO for small-scale problems.

Results' Comparison for Medium-Scale Instances.
To evaluate the performance of IWOA, 37 experiments are conducted for the medium-scale problems.The experiments' MPD, SD, and ACT are investigated in Table 5 for medium-scale instances.The last column of this table represents the deviation from best feasible solution and the gaps for solutions obtained from the GAMS CPLEX solver.Comparative analysis of GA and PSO with GA-MPC and SPSO indicated that MPDs of all the experiments for GA and PSO were better than their new versions.However, ACT of GA-MPC is better than GA and ACT of PSO is better than SPSO.Due to the significance of the solutions' quality in this research, traditional GA and PSO are used to evaluate IWOA for the small-scale problems.
The graphical analysis of the results, as illustrated in Figure 6, implies that the performance of IWOA is better than WOA for all the instances considering the MPD matrix.On the other hand, the performance of IWOA is better than GA for 89.2% of the experiments, and also it is better than PSO for most of the instances.Furthermore, one-way ANOVA test is conducted to compare IWOA with WOA as well as GA and PSO which had better performances.Although the graphical analysis illustrates a better performance for IWOA compared with GA and PSO, the results of ANOVA test on MPDs of these algorithms as illustrated in Table 6 and Figure 7 indicate On the other hand, MPD of IWOA and nonzero PDs of the results obtained from the GAMS CPLEX solver which is illustrated in Figure 4 reveal the worse performance of standard software and deterioration of solutions when the number of variables increases.
Figure 8 represents the comparison of SD for IWOA, WOA, and GA.Smaller SD of the experiments for IWOA confirms its improved robustness and stability of the solutions compared with WOA, GA, and PSO for the medium instances.7 represents the results of the experiments on 42 large-scale instances to compare the performance of different metaheuristic algorithms.In our experiments, GA-MPC had better performance than GA considering both MPD and ACT metrics.Comparing MPDs of GA with GA-MPC indicated better performance of GA-MPC for 76.2% of the experiments.Furthermore, comparing PSO with SPSO with respect to MPD metrics revealed the better performance of SPSO for 88.1% of the test problems.However, ACT of PSO is better than SPSO.Due to the importance of solutions' quality, we have evaluated IWOA using GA-MPC and SPSO for the large-scale problems.The graphical analysis of the results as illustrated in Figure 9, indicates that the performance of IWOA is better than WOA for all of the instances.Moreover, IWOA has better or same performance in 61.2% of the experiments compared with SPSO.However, the performance of GA-MPC is better than IWOA in 83.3% of the experiments.Furthermore, the comparative graphical analysis of IWOA and GA-MPC also reveals that the performance of IWOA has deteriorated as the number of decision variables increased.
On the other hand, one-way ANOVA test is conducted to compare the performance of IWOA with WOA as well as GA-MPC and SPSO which had better performances.The ANOVA test's results as investigated in Table 8 and Figure 10 indicate that the performance of IWOA is strongly better than that of WOA which confirms the results of the graphical analysis.Moreover, the ANOVA test's results show that IWOA has almost similar competence with SPSO while GA-MPC is more competent than other algorithms for the large-scale problems.However, IWOA reveals a  19 Complexity similar performance compared with SPSO considering MPD of the observations.Figures 11 and 12 represent almost the same performance for IWOA compared with WOA and GA-MPC with respect to solutions' SD.

Conclusion and Future Studies
In this study, we proposed a formulation for the integrated scheduling of production and distribution operations in a multifactory supply chain.Considering the make-to-order production system and direct shipment of finalized orders as well as the limited number of transportation vehicles increased the complexity of the problem.The integrated simulation-based solution approach that considers the interactions of subproblems also provides the best trade-off among the individual costs of the objective function.The advantage of our work is that it optimizes the scheduling of transportation integrated with production in the operational level of multifactory production systems.We used a conventional software to solve the small-size problems to evaluate the validity of the proposed model.We improved the newly introduced WOA to solve the medium-scale problems and computationally confirmed the efficiency of the modified algorithm (IWOA) compared to basic WOA for different scales of the problem.To evaluate IWOA, we compared the results of computational experiments with two versions of GA and PSO.The results of experiments confirmed a significant improvement in the mean value of the percentage of deviations while using IWOA rather than WOA for all sizes of the problem.However, unlike the improvements made to WOA and the better performance of IWOA for the medium-scale instances, it is better to use GA-MPC for the large-scale problems.
We considered parallel machines and immediate direct shipment of the orders in the problem of this research; other configurations of machines and shipping methods could be the trend of future studies.On the other hand, an integrated supply chain with a single objective function is considered in this paper.So, independent factories and multiobjective problem can be presented in the future.Furthermore, one can consider the application of the MIPDSP for the loading and distribution operations, and the general evaluation of IWOA as new directs of research.

Figure 3 :
Figure 3: Intelligent procedure based on simulation and metaheuristic algorithms.

Figure 8 :
Figure 8: SD comparison of IWOA with WOA and GA for medium-scale instances.

Figure 11 :
Figure 11: SD comparison of IWOA with WOA for large-scale instances.

Figure 12 :Figure 10 :
Figure 12: SD comparison of IWOA with GA-MPC for large-scale instances.

Table 1 :
Definition of sets, indexes, parameters, and variables.

Table 4 :
Average of objective values for metaheuristic algorithms and mean PDs

Table 6 :
Results of one-way ANOVA test for medium-scale instances.

Table 7 :
Mean PD % , standard deviation SD × 10 8 , and average of CPU time (s) for large-scale instances MPD SD ACT MPD SD ACT MPD SD ACT MPD SD ACT MPD SD ACT Complexity similar competence for all three algorithms.Furthermore, the ANOVA test's results imply that MPD of IWOA is meaningfully better than WOA which confirms the graphical analysis results.

Table 8 :
Results of the one-way ANOVA test for large-scale instances DF: degree of freedom; SS: sum of squares; MS: mean squares 18 Complexity 5.4.3.Results' Comparison for the Large-Scale Instances.Table