Advanced Metaheuristic Method for Decision-Making in a Dynamic Job Shop Scheduling Environment

: As a well-known NP-hard problem, the dynamic job shop scheduling problem has signiﬁ-cant practical value, so this paper proposes an Improved Heuristic Kalman Algorithm to solve this problem. In Improved Heuristic Kalman Algorithm, the cellular neighbor network is introduced, together with the boundary handling function, and the best position of each individual is recorded for constructing the cellular neighbor network. The encoding method is introduced based on the relative position index so that the Improved Heuristic Kalman Algorithm can be applied to solve the dynamic job shop scheduling problem. Solving the benchmark example of dynamic job shop scheduling problem and comparing it with the original Heuristic Kalman Algorithm and Genetic Algorithm-Mixed, the results show that Improved Heuristic Kalman Algorithm is effective for solving the dynamic job shop scheduling problem. The convergence rate of the Improved Heuristic Kalman Algorithm is reduced signiﬁcantly, which is beneﬁcial to avoid the algorithm from falling into the local optimum. For all 15 benchmark instances, Improved Heuristic Kalman Algorithm and Heuristic Kalman Algorithm have obtained the best solution obtained by Genetic Algorithm-Mixed. Moreover, for 9 out of 15 benchmark instances, they achieved signiﬁcantly better solutions than Genetic Algorithm-Mixed. They have better robustness and reasonable running time (less than 30 s even for large size problems), which means that they are very suitable for solving the dynamic job shop scheduling problem. According to the dynamic job shop scheduling problem applicability, the integration-communication protocol was presented, which enables the transfer and use of the Improved Heuristic Kalman Algorithm optimization results in the conventional Simio simulation environment. The results of the integration-communication protocol proved the numerical and graphical matching of the optimization results and, thus, the correctness of the data transfer, ensuring high-level usability of the decision-making method in a real-world environment.


Introduction
The use of evolutionary computational methods has been used in decision-making processes in production systems for many years. Heuristic and metaheuristic methods enable obtaining satisfactorily good optimization results for a wide variety of NP-hard decision problems, where a decision problem H is NP-hard when for every problem L in NP, there is a polynomial-time many-one reduction from L to H. The globalized market and digitally supported industry, regardless of the type of production, from the most basic job shop to mass-personalized production, aim to optimize production systems. The paper presents a new metaheuristic method based on the Heuristic Kalman Algorithm [1] which enables optimal scheduling of Dynamic Job Shop Scheduling Problem (DJSSP), and the exchange and use of optimization results in a conventional simulation environment. Within comprehensively from the various production parameters. The research work shows the evaluation of parameters: Machine utilization, machine process time and orders' flow time. Validation of optimization and simulation results, however, can confirm or refute the importance of proper DJSSP planning and scheduling.
Mathematics 2021, 9,909 3 of 23 protocol of optimization results` integration into the conventional Simio simulation environment, as presented in Figure 1. With the data transfer between the decisive IHKA algorithm and the simulation model in Simio, it is possible to evaluate the DJSSP optimization problem comprehensively from the various production parameters. The research work shows the evaluation of parameters: Machine utilization, machine process time and orders` flow time. Validation of optimization and simulation results, however, can confirm or refute the importance of proper DJSSP planning and scheduling. This research paper is organized as follows: in the second section, a problem description of DJSSP with corresponding mathematical model representation is given. This section is followed by the third section, where a newly proposed metaheuristic algorithm IHKA is presented. Next, the experimental design section including the application of IHKA to DJSSP in terms of the evaluation process and simulation model structure proposal is presented. In section five detailed numerical and graphical results of algorithm comparisons are made, followed by the evaluation of the simulation model approach in terms of transformation of optimization results from IHKA to conventional simulation modeling environment. Finally, discussion and conclusions are presented, evaluating results, suggestions, and opportunities for future work.

Dynamic Job Shop Scheduling Problem
A Dynamic Job Shop Scheduling Problem (DJSSP) is a combinatorial optimization problem in which we have job orders given in the introduction and ′ job orders that arrive after the start of scheduling the number of job orders. All job orders ( and ′) must be executed on available machines, considering the limitations of the DJSSP optimization problem [19]:  All machines from a set of are available at time 0;  An individual operation can only be performed on one machine at a time;  A single machine can only perform one operation at a time;  The operation performed on machine can only be interrupted when there is a dynamic event of machine breakdown appearing;  The next operation can be performed only when the previous one is completed; This research paper is organized as follows: in the second section, a problem description of DJSSP with corresponding mathematical model representation is given. This section is followed by the third section, where a newly proposed metaheuristic algorithm IHKA is presented. Next, the experimental design section including the application of IHKA to DJSSP in terms of the evaluation process and simulation model structure proposal is presented. In section five detailed numerical and graphical results of algorithm comparisons are made, followed by the evaluation of the simulation model approach in terms of transformation of optimization results from IHKA to conventional simulation modeling environment. Finally, discussion and conclusions are presented, evaluating results, suggestions, and opportunities for future work.

Dynamic Job Shop Scheduling Problem
A Dynamic Job Shop Scheduling Problem (DJSSP) is a combinatorial optimization problem in which we have n job orders given in the introduction and n job orders that arrive after the start of scheduling the n number of job orders. All job orders (n and n ) must be executed on m available machines, considering the limitations of the DJSSP optimization problem [19]: The operation performed on machine i can only be interrupted when there is a dynamic event of machine breakdown appearing; • The next operation can be performed only when the previous one is completed; • The process times of the operation and the assigned machine i are known in advance. During individual operation, the processing time may change due to a dynamic event of changing the original operation processing time;

•
The original operation processing time; • Setup times do not depend on the operation execution sequence and the machine on which the operation will be performed, but are included in the processing time of the operation; • Transport time between machines is 0. Problem formulation was made using the following notations for decision variables, datasets and parameters: competition time of job j on machine i c ij competition time of new job j' on machine i y ij starting time of operation (i, j) y ij starting time of new operation (i, j ) rp the start time of the rescheduled job order tm i the start that the machine i will be idled at the start of rescheduling job order In the mathematical description of the single-objective optimization problem, the DJSSP focuses on the minimizing job's makespan (completion time of all jobs) C max , where c ij is the completion time of job j on one machine i, described by Equation (1).
Two constraints (Equations (2) and (3)) ensure that C max is greater than or equal to the completion time of job j and new job j on machine i.i. C max ≥ c ij , where i = 1, . . . , m and j = 1, . . . , n C max ≥ c ij , where i = 1, . . . , m and j = 1, . . . , n The makespan of individual operation of orders n and new orders n is determined by Equations (4) and (5).
c ij = y ij + p ij , where i = 1, . . . , m and j = 1, . . . , n Equations (6) and (7) define the optimization problem constraint condition that the next operation can be performed only when the previous one is completed.
In addition, some constrains, presented by Equations (8)- (11), must be made according to requirement that only one job can be processed on a machine at any time. where conditions: z ijk = 1, if J j precedes J k on machine M i (z ijk = 0, otherwise), z ij k = 1, if J j precedes J k on machine M i (z ij k = 0, otherwise), t ij = 1, if J j will be processed on machine M i after rescheduling (t ij = 0, otherwise) and t ij = 1, if J j will be processed in machine M i after rescheduling (t ij = 0, otherwise). The start time of the rescheduled job orders and the start time of the machine i at the start of the rescheduled job orders is defined by Equations (12) and (13).
y ij ≥ (tm i + rp) * t ij , where i = 1, . . . , m and j = 1, . . . , n As stated in the literature [19] dynamic scheduling can be divided into reactive, predictive-reactive, or pro-active scheduling. This paper uses the predictive-reactive scheduling approach incorporated in the proposed metaheuristic algorithm. For a more real-world usable method, the continuous rescheduling approach was added, where rescheduling is performed each time a dynamic event, new job n , arrives. The continuous rescheduling approach gives advantages when the optimization results of DJSSP are transferred into the simulation modeling environment.

Metaheuristic Algorithm
A metaheuristic algorithm is introduced in order to solve the DJSSP. Following the introduction of the original metaheuristic algorithm, an improved version is proposed.

Heuristic Kalman Algorithm
The Heuristic Kalman Algorithm (HKA) is a Kalman filtering-based heuristic approach, which has been proposed by Toscano and Lyonnet for solving continuous and non-convex optimization problems [1,20]. The HKA is very easy to use, for only a few parameters (only three) need to be set by the user. In the HKA, the search heuristic is entirely different from other population-based stochastic optimization algorithms, which is explicitly considered in the optimization problem as a measurement process designed to give an estimate of the optimum. Through the measurement process, the HKA develops a specific procedure based on the Kalman estimator to improve the quality of the estimate obtained. In the practical implementation, the HKA needs to initialize parameters first, such as the Gaussian distribution, the user-defined parameters and the stopping rule. During the optimization process of the HKA, first, the Gaussian distribution parametrized by a given mean vector with a given variance-covariance matrix generates the solutions, then the evaluation process of the generated solutions, followed by the measurement procedure, and, finally, the optimum estimator of the parameters of the random generator is introduced for the next generation. Figure 2 shows the flowchart of the HKA [1,20].

Improved Heuristic Kalman Algorithm
Due to its very strong convergence, the HKA can easily fall into the local minimum [18,21]. Therefore, the HKA needs to be improved. By introducing the cellular neighbor network, this paper proposes an Improved Heuristic Kalman Algorithm (IHKA). In the IHKA, a function is introduced that handles the boundary constraint after the solutions are generated by the Gaussian distribution, then the best position of each individual is stored and updated after the evaluation process of the generated solutions, and, finally, the cellular neighbor network is introduced during the measurement process to select the best samples under consideration. The general pseudo-code of the IHKA is shown in Algorithm 1 and the flowchart of the IHKA is shows in Figure 3 [21]. atics 2021, 9, 909 6 of 23

Improved Heuristic Kalman Algorithm
Due to its very strong convergence, the HKA can easily fall into the local minimum [18,21]. Therefore, the HKA needs to be improved. By introducing the cellular neighbor network, this paper proposes an Improved Heuristic Kalman Algorithm (IHKA). In the IHKA, a function is introduced that handles the boundary constraint after the solutions are generated by the Gaussian distribution, then the best position of each individual is stored and updated after the evaluation process of the generated solutions, and, finally, the cellular neighbor network is introduced during the measurement process to select the best samples under consideration. The general pseudo-code of the IHKA is shown in Algorithm 1 and the flowchart of the IHKA is shows in Figure 3 [21].

Initializing the Variables
In the IHKA, the initialization parameter is followed by the initialization variable. In the process of initializing variables, the mean vector, the variance vector, the best position of each individual and the best position of the population are initialized.
The mean vector and the variance vector are the parameters of the Gaussian distribution generator. To cover the entire search space, they, denoted m and S, are initialized as follows [1]: where lu(1, j) (respectively, lu(2, j)) is the lower (respectively, upper) boundary of the j-th dimension of the problem. The best position of each individual is stored to construct the cellular neighbor network. In the initialization process, the best position and fitness of each individual are initialized. The best position of the individual is initialized as a zero vector. Its fitness is initialized to infinity (respectively, negative infinity) for the minimization (respectively, maximization) problem. For minimization problem, they, denoted Pbx and Pb, are initialized as follows: where zeros (respectively, in f ) is a function that produces a matrix of all zeros (respectively, infinity). The best position of the population is used to handling the constraints. When a certain dimension of an individual violates the constraints, it may be replaced by the corresponding dimension of the best position of the population. It needs to be used during the first processing of handling the constraints, so it is randomly initialized in the search space. They, denoted Gbx and Gb, are initialized as follows: where rand is a function to generate random numbers in the interval (0, 1) and f is a function to evaluate the fitness of the individual. Step 0: Initialization. Set the number of rows of the cellular neighbor network r, the number of columns of the cellular neighbor network c, the neighbor type of the cellular neighbor network nt, the rewiring probability of the cellular neighbor network p, the minimum maximum depth of the cellular neighbor network d, the k nearest neighbors of the cellular neighbor network k (Note that k is set to N ξ − 1), the size of the population N = r * c, the number of dimensions of the practical problem Nd, the number of best samples under consideration N ξ , the slowdown coefficient α and the maximum number of iterations MaxIter.
Initialize the mean vector m and the variance vector S by Equations (14) and (15), respectively. Initialize the best position Pbx and fitness Pb of each individual by Equations (16) and (17). Initialize the best position Gbx and fitness Gb of the population by Equations (18) and (19). Initialize the cellular neighbor network by Algorithm 2: cnn = createcnn(r, c, nt, p, d, k) where createcnn is a custom function that generates a cellular neighbor network.
Step 1: Iteration. for ite = 1 : MaxIter Step 1.1: Random generator. Generate a population x with N individuals by a Gaussian distribution parametrized by m and S: where mvnrnd(.) is a function that generates random vectors from the multivariate normal distribution and diag(.) is a function that generates diagonal matrices or diagonals of a matrix.
Step 1.2: Handling the constraints of the problem by Equation (20).
Step 1.3: Evaluate fitness. Calculate the individual fitness.
Step 1.4: Update the individual best position by Equation (21).
Step 1.5: Update the global best position.
Step 1.6: Measurement process. The N ξ best samples under consideration are selected from the k nearest neighbors of the global best position in the cellular neighbor network, which is composed of the best positions of individuals: where Gbi is the global best position index of the current population and knni is a function to find the k nearest neighbors (containing itself) in the cellular neighbor network. Note that neighbors of each individual can be obtained in advance to reduce the running time.
Compute the measurement ξ and the variance vector V by Equations (22) and (23).

Initializing the Variables
In the IHKA, the initialization parameter is followed by the initialization variable. In the process of initializing variables, the mean vector, the variance vector, the best position of each individual and the best position of the population are initialized.
The mean vector and the variance vector are the parameters of the Gaussian distribution generator. To cover the entire search space, they, denoted and , are initialized as follows [1]: where (1, ) (respectively, (2, )) is the lower (respectively, upper) boundary of the -th dimension of the problem.

Handling the Constraints
To handle the boundary constraints, a handling constraints function is introduced in the IHKA. If a dimension value violates the boundary constraint, it will be replaced. According to a large number of preliminary experiments, it has a 50% probability that it will be replaced by a random number. The probability of being replaced by the corresponding dimension value of the global best position is 25%. Otherwise, it is replaced by a boundary value. If it is greater than the upper boundary value, it will be replaced by the upper boundary value, otherwise, it will be replaced by the lower boundary value. Note that, if the boundary value cannot be taken, it will be replaced with a value that is infinitely close to the boundary value. The rule for handling the constraints is used as follows: where x(i, j) represents the j-th dimension of the i-th individual who violated the boundary constraints in the population x, and ra is a random number generated by the random function.

Updating the Individual Best
The best position in the history of each individual is stored in the IHKA. They are used to construct the cellular neighbor network to ensure the convergence of the IHKA. If its fitness is not worse than the original, the newly generated individual will replace the individual at the corresponding position in the cellular neighbor network. The rule for updating the individual best is used as follows: where f (i) represents the fitness of the newly generated individual i.

Introducing the Cellular Neighbor Network
During the measurement process, the cellular neighbor network is introduced to select the best samples under consideration in the IHKA. The cellular neighbor network was proposed by Li et al. for a cellular neighbors-based structure can achieve good performance in human society or a network [22]. In the cellular neighbor network, each individual is composed of individuals at corresponding positions in the population. Each individual in the cellular neighbor network uses the best position explored by the corresponding position individual in the population during the optimization process. In a two-dimensional space, there are two well-known structures in the neighborhood of a cell, which are the Von Neumann neighborhood [23] and the Moore neighborhood [24]. Their neighborhood structure is shown in Figures 4a and 4b, respectively. We set the minimum maximum depth of the cellular neighbor network to ensure that each node in the network is connected. The construction of the cellular neighbor network is described in detail in Algorithm 2. The instances of cellular neighbor networks constructed by the proposed algorithm are shown in Figures 4c and 4d, respectively. Algorithm 2. The general pseudo-code for creating cellular neighbor network.
Input: the number of rows of the cellular neighbor network r, the number of columns of the cellular neighbor network c, the neighbor type of the cellular neighbor network nt, the rewiring probability of the cellular neighbor network p, the minimum maximum depth of the cellular neighbor network d, the k nearest neighbors of the cellular neighbor network k.
Output: the cellular neighbor network.Initialize the network.
Step 0: Initialize a regular cellular neighbor network based on the type of cellular neighborhood.
Step 1: Create the network. Construct a cellular neighbor network based on the rewiring probability of the cellular neighbor network p. After cutting, the minimum maximum depths of the two nodes selected for deletion are checked, if the requirements are not met, the disconnection is canceled. Note that this check does not affect the addition of new random connections.
Step 2: Obtain k nearest neighbors. Obtain the k nearest neighbors (containing itself) of each node in the network.
Their neighborhood structure is shown in Figure 4a and 4b, respectively. We set the minimum maximum depth of the cellular neighbor network to ensure that each node in the network is connected. The construction of the cellular neighbor network is described in detail in Algorithm 2. The instances of cellular neighbor networks constructed by the proposed algorithm are shown in Figure 4c and 4d, respectively. Input: the number of rows of the cellular neighbor network , the number of columns of the cellular neighbor network , the neighbor type of the cellular neighbor network , the rewiring probability of the cellular neighbor network , the minimum maximum depth of the cellular neighbor network , the nearest neighbors of the cellular neighbor network .

Computing the Measurement
In the measurement process, N ξ best samples under consideration are selected, followed by the calculation of the measurement and the variance vector. They, denoted ξ and V are define as follows [1]:

Computing the Posterior Estimation
In the optimal estimation process, the mean vector and the variance vector are calculated for the next step. They are the parameters of the Gaussian distribution generator of the next-generation population. They, denoted m_pe and S_pe, are define as follows [1]: where L, W and τ are unknown vector which have to be determined to ensure optimal estimation, a is the slowdown factor, mean (.) is a function that calculates the mean value and the symbol./ (respectively, . * ) stands for a component-wise divide (respectively, product).

Experiment Design
In order to verify the performance of the proposed metaheuristic algorithm, this paper designs an experiment. The experiment contains the numerical experiment of the algorithm and the simulation experiment of the optimization results of the algorithm for the benchmark instances of DJSSP.

Applying the Algorithm on DJSSP
The HKA was originally proposed to solve continuous problems, so we introduced the relative position indexing [25] to encode. The instance of the solution encoding is shown in Figure 5. Each dimension of the individual is initialized randomly to an open interval of 0 and 1 (see Figure 5 S 0 ). The position of the individual is transferred into the discrete domain by using the relative position index (see Figure 5 S 1 ). That is, the value of each dimension of the individual represents the priority, the smaller the value, the higher the priority it will be processed with. For instance, in Figure 5 S 0 , the value of the third (respectively, sixth) position is the smallest (respectively, largest), so its processing priority is the highest (respectively, lowest), that is, it is ranked first (respectively, last) in the processing priority. Then the processing sequence of the jobs can be obtained from the DJSSP Table (see Figure 5 S 2 ). That is, sort the "Jobs" column in the DJSSP Table  (not including the rows of the machine breakdowns). In Figure 5 S 2 , the jobs processing sequence is as follows: (2, 2), (1, 1), (4, 1), (2, 1), (1, 2), (4, 2), (3,2), (3,1), where (., .) means job number and machine number. For instance, (2, 2) means job 2 is processed on machine 2. In Figure 5 S 2 , the k-th occurrence of a job number represents the k-th operation of the job. Finally, the final solution can be obtained (see Figure 5 S 3 ). In Figure 5 Table) and machine number. For instance, (3,2), means that the third row in the DJSSP Table  (job 2) is processed on machine 2.

Evaluating the Algorithm
In order to evaluate the proposed algorithm, we introduced the benchmark instances of DJSSP and set the parameters of the algorithm.

Selecting the Benchmark Instance
The benchmark instances of DJSSP were selected in order to verify the performance of the proposed algorithm to solve the DJSSP. All selected benchmark instances of DJSSP are from [19], which were generated based on Lawrence static job shop benchmark problem instances. In all the benchmark instances of DJSSP, each instance considers dynamic events such as machine breakdown, new job arrival and change in the processing time. In [19], without considering the arrival of new jobs, 15 benchmark instances of DJSSP are divided into three categories: Small (less than 60), medium (greater than, or equal to 60 and less than 100) and large (greater than, or equal to 100) according to the number of operations. For the benchmark instances of DJSSP see Table 1. The detailed data of all the selected benchmark instances of DJSSP can be obtained from [19]. In the benchmark instances of DJSSP, it focuses on the minimizing job's makespan (completion time) C max , described by Equation (1), that is, the optimal fitness. Figure 5 2 , the -th occurrence of a job number represents the -th operation of the job. Finally, the final solution can be obtained (see Figure 5 3 ). In Figure 5 3 , the operations processing sequence is as follows: (3,2), (1, 1), (7, 1), (4, 1), (2, 2), (8, 2), (5, 2), (6, 1), where (., .) means number (corresponding to the number in the "No" column in the DJSSP Table) and machine number. For instance, (3,2), means that the third row in the DJSSP Table (job 2) is processed on machine 2.

Evaluating the Algorithm
In order to evaluate the proposed algorithm, we introduced the benchmark instances of DJSSP and set the parameters of the algorithm.

Selecting the Benchmark Instance
The benchmark instances of DJSSP were selected in order to verify the performance of the proposed algorithm to solve the DJSSP. All selected benchmark instances of DJSSP are from [19], which were generated based on Lawrence static job shop benchmark

Defining the Parameters
In this paper, the MaxIter of the IHKA and HKA was set according to a large number of preliminary experiments, and the others were set according to the previous literature. The parameters r, c, k, N, N ξ and α were set according to [18,21]. The parameters p and d were set according to [22] and [26], respectively. In the large number of preliminary experiments, the MaxIter was set to a different value, and the algorithm was independently run many times to solve each of the selected 15 benchmark instances, and the convergence process of the algorithm was observed. If it converged after the maximum iteration of each independent run, the MaxIter corresponding to the benchmark instance was determined. According to the large number of preliminary experiments, the MaxIter was set to a different value corresponding to the size type of each benchmark instance. The parameters of the IHKA and HKA are shown in Table 2. The parameters defined in Table 2 will be used in all the experiments.

Transferring Optimization Results to Simulation Evaluation
Given the basic purpose of the DJSSP optimization problem, in which dynamic events such as machine breakdowns (MB), operation process time changes (PTC) and new job arrivals (NJA), occur dynamically, this optimization problem shows the most realistic problem of job shop scheduling using the IHKA, and the high ability to schedule the job and operations of the DJSSP optimization problem has been proven. In the next experimental phase, it is necessary to transfer the optimization results into the real production system environment. To transfer the optimization results of the IHKA to the real-world environment, the simulation modeling method in the commercially available Simio software environment is presented the simulation model of benchmark dataset 6 × 5 is shown in Figure 6. Established communication protocol was introduced when transferring the optimization results of the IHKA to the commercial simulation environment [27].

Transferring Optimization Results to Simulation Evaluation
Given the basic purpose of the DJSSP optimization problem, in which dynamic events such as machine breakdowns (MB), operation process time changes (PTC) and new job arrivals (NJA), occur dynamically, this optimization problem shows the most realistic problem of job shop scheduling using the IHKA, and the high ability to schedule the job and operations of the DJSSP optimization problem has been proven. In the next experimental phase, it is necessary to transfer the optimization results into the real production system environment. To transfer the optimization results of the IHKA to the real-world environment, the simulation modeling method in the commercially available Simio software environment is presented the simulation model of benchmark dataset 6×5 is shown in Figure 6. Established communication protocol was introduced when transferring the optimization results of the IHKA to the commercial simulation environment [27]. The simulation model includes five machines (from M1 to M5), a warehouse of input material (jobs` arrival point) and a warehouse of finished products (jobs` finished point), and the transport of semi-finished orders between stations is provided by a forklift. Table  3 shows the optimization results of the IHKA transmitted via the MATLAB-Excel-Simio communication protocol to a data-driven simulation model. The 6×5 DJSSP dataset used contains six initial jobs (from J1 to J6) and three additional job arrivals (J7 to J9). Four jobs are static (J1, J2, J5 and J6), and in addition five jobs are dynamic (J3, J4, J7, J8, J9). In dynamic job J3, two dynamic events occur: Interruption due to machine breakdown in operation J3,1, which is 9 min long, and change in process time in operation J3,3, from the original process time of 19 min to 23 min.  The simulation model includes five machines (from M 1 to M 5 ), a warehouse of input material (jobs' arrival point) and a warehouse of finished products (jobs' finished point), and the transport of semi-finished orders between stations is provided by a forklift. Table  3 shows the optimization results of the IHKA transmitted via the MATLAB-Excel-Simio communication protocol to a data-driven simulation model. The 6 × 5 DJSSP dataset used contains six initial jobs (from J 1 to J 6 ) and three additional job arrivals (J 7 to J 9 ). Four jobs are static (J 1 , J 2 , J 5 and J 6 ), and in addition five jobs are dynamic (J 3 , J 4 , J 7 , J 8 , J 9 ). In dynamic job J 3 , two dynamic events occur: Interruption due to machine breakdown in operation J 3,1 , which is 9 min long, and change in process time in operation J 3,3 , from the original process time of 19 min to 23 min. Similar dynamic events occur in job J 4 : an interruption is performed due to a machine breakdown in operation J 4,1 , with duration of 5 min and a change in process time in operation J 4,4 from the original process time of 93 min to 48 min. The other three dynamic jobs are J 7 , J 8 and J 9 represent new job arrivals. In a dynamic machine breakdown event, the simulation model uses a data table with the event time-based interruption function assigned to a specific operation and machine. The time of execution of the breakdown event is defined by the Time to Repair function, in which the repair time value is specified: For J 3,1 the repair time is 9 min, and for J 4,1 the time is 5 min, respectively. This value is added to the initial value of the process time of the operation, and, as a sum, represents the process time of the entire operation. The dynamic event of the change of the operation process time is modeled with a data-driven table of process times, in which the values of the initial and updated process times are given. The simulation model, when new job arrivals appear, obtains data on these via MatLab-Excel-Simio communication, dynamically, according to the performed optimization with the IHKA. The arrival of a new job in the existing order queue is assigned by the new job arrival variable. Despite the fact that the presented simulation model is limited to displaying the 6 × 5 DJSSP dataset, the proposed build-in decision logic is adaptable, and it is possible to change the simulation model continuously according to the characteristics of the devalued production system.

Result
The corresponding experimental results were obtained after the experimental design, numerical experiments and simulation experiments were executed in the MATLAB and Simio software.

Results Comparison to Other Algorithms
The original HKA and the Genetic Algorithm-Mixed (GAM) [19] were selected as a comparison to prove the performance of the proposed IHKA to solve the DJSSP. In order to prove the robustness, the IHKA and HKA were run independently, 30 times for each benchmark instance of DJSSP. All the algorithms in this paper are implemented in MATLAB language. All experiments were simulated in version R2020b of MATLAB. The computer used in all experiments was a PC with x64-processor Intel(R) Core(TM) i7-8550U CPU @ 1.80 GHz 1.99 GHz and 16 GB RAM in the environment of the Windows 10 version 20 H2. Figure 7 shows the IHKA and HKA average convergence of the best solutions for the benchmark instances of DJSSP. The convergence curve in the Figure is the average of 30 independent runs of each algorithm. It can be seen from the Figure that the IHKA reduced the convergence rate significantly, which is beneficial to avoid the algorithm from falling into the local optimum during the optimization process. In the IHKA, the influence of different cellular neighborhood structures on the convergence process is not very obvious. In the IHKA, the average convergence process of 30 independent runs of two cellular neighborhood structures, the Von Neumann neighborhood and the Moore neighborhood, are almost the same. The Figure can also clarify that both HKA and IHKA can converge to the best solution for all benchmark instances of DJSSP. For most of the benchmark instances of DJSSP, such as 10 × 8, 10 × 9, 20 × 5, 10 × 10, 12 × 10, 13 × 10, 20 × 7 and 15 × 10, both IHKA and HKA can converge significantly to the best solution obtained by GAM. Therefore, the convergence performance of IHKA and HKA was clarified to be very good. Moreover, we can also see that, in some instances, such as 10 × 8, 10 × 9, 12 × 10, 13 × 10, 20 × 7 and 15 × 10, the IHKA algorithm still shows a trend of convergence at termination. It means that better solutions may be obtained by increasing the number of iterations of the algorithm for each of these instances. The computational statistics of the IHKA and HKA on the fitness for the benchmark instances of DJSSP are shown in Table 4. In Table 4, for each benchmark instance, if the algorithm has the best statistical performance, then the row is bolded, otherwise the minimum value is bolded. It can be clarified from the Table 4 that both the IHKA and HKA performed very well. For most of the benchmark instances of DJSSP, such as 6 × 8, 10 × 8, 10 × 9, 20 × 5, 10 × 10, 12 × 10, 13 × 10, 20 × 7 and 15 × 10, the best solutions obtained by IHKA and HKA were better than those obtained by GAM. As an example, the best solution and the Gantt chart of the 6 × 5 DJSSP with IHKA VonNeumann are given in Table 5 and Figure 8, respectively. In Figure 8, the yellow is used to represent the period of machine breakdown. For 30 independent runs of all the benchmark instances of DJSSP, the IHKA VonNeumann had a 100% success rate (SR) for obtaining the best solution at least better than that obtained by GAM. Except for the 10 × 9 DJSSP, the success rate of The computational statistics of the IHKA and HKA on the fitness for the benchmark instances of DJSSP are shown in Table 4. In Table 4, for each benchmark instance, if the algorithm has the best statistical performance, then the row is bolded, otherwise the minimum value is bolded. It can be clarified from the Table 4 that both the IHKA and HKA performed very well. For most of the benchmark instances of DJSSP, such as 6 × 8, 10 × 8, 10 × 9, 20 × 5, 10 × 10, 12 × 10, 13 × 10, 20 × 7 and 15 × 10, the best solutions obtained by IHKA and HKA were better than those obtained by GAM. As an example, the best solution and the Gantt chart of the 6 × 5 DJSSP with IHKA VonNeumann are given in Table 5 and Figure 8, respectively. In Figure 8, the yellow is used to represent the period of machine breakdown. For 30 independent runs of all the benchmark instances of DJSSP, the IHKA VonNeumann had a 100% success rate (SR) for obtaining the best solution at least better than that obtained by GAM. Except for the 10 × 9 DJSSP, the success rate of IHKA Moore and HKA was also 100%. In 30 independent runs of each instance, the IHKA obtained the best solution for most instances, which was 12 (considering the best solution obtained by IHKA VonNeumann and IHKA Moore) out of 15 instances. It shows that the IHKA can escape the local minimum, which is better than the HKA that achieved 11 out of 15 instances. However, both IHKA VonNeumann and IHKA Moore are not as good as HKA, respectively, they were 10 out of 15 instances and 9 out of 15 instances. The IHKA VonNeumann performed better than the IHKA Moore, that is, it is better for the Von Neumann neighborhood structure to avoid the IHKA from falling into the local minimum than the Moore neighborhood structure.  Table 5. The best solution of the 6 × 5 DJSSP with IHKA VonNeumann.

Simulation Modeling Results
The main purpose of using the simulation model was the data transfer between the numerical results of the optimization IHKA and the real, dynamically changing production system. Given the globally highly competitive environment and the need for a flexible and dynamically fast-responsive production system, a fashion that enables the connection between the decision method and the real environment is crucial. Table 7 shows the job's sequence according to the assigned machine and the used variable Job Priority Number, which determines the sequence of performing operations on an individual machine specifically. The variable used allows the simulation model to use the optimization data generated by the IHKA. The specifically assigned value of the Job Priority Number variable for an individual operation also defines the dynamic events type (MB, PTC, or NJA) and characteristics that run during the simulation model execution. Given the limitations of decision algorithms related to the lack of connections with the real-world environment, the proposed simulation model allows the validation of key production system parameters. Table 8 presents the numerical values of the machine utilization parameter and the total process time of an individual machine. The values in the Table prove that the IHKA allows a high degree of scheduling justification from the point of view of high machine utilization, as the average value is 78.9%, which is very high, considering that the data describe job shop production. However, with a high theoretical machine utilization, as shown in the M 3 machine, its performance needs to be evaluated in detail, as such a high theoretical utilization could lead to a bottleneck in the production system. A detailed analysis confirmed, from the graphical results in Figure 7, that the job scheduling performed with the IHKA was able to distribute the optimal schedule, and, thus, allow the smooth operation of the production system. Numerical values of the machine process time parameter prove the adequacy of communication between the IHKA optimization algorithm and the simulation model, as the sums of process times of the simulation model results match the process times assumed in the datasets and terminated by the IHKA. A detailed analysis of the simulation results evaluates the order flow time. The flow time of the order depends on the time of performing an individual operation on the machine and a dynamic event that can extend or shorten its duration. Table 9 shows the flow times of orders, and identifies the dynamic events that occurred during model execution. The numerical results prove that the longest order flow times are those in which the dynamic event of an unexpected machine breakdown occurred. This time affected the total order flow time and, consequently, the total average order flow time. It can be found that, in the event of a machine breakdown, it is important to reduce its time to a minimum, if possible, or to adjust the time of performing upcoming operations as shown in the job J 4 . The result of the simulation modeling is a graphical representation of the scheduling of the working tasks of the DJSSP dataset 6 × 5. The Gantt chart in Figure 9 corresponds exactly to the Gantt chart in Figure 8, proving that the presented method of using the optimization results of the IHKA enable communication between the decision algorithm and the real-world environment staged in the simulation environment.  The result of the simulation modeling is a graphical representation of the scheduling of the working tasks of the DJSSP dataset 6 × 5. The Gantt chart in Figure 9 corresponds exactly to the Gantt chart in Figure 8, proving that the presented method of using the optimization results of the IHKA enable communication between the decision algorithm and the real-world environment staged in the simulation environment.

Discussion
The IHKA proposed in this paper was applied to solve the 15 benchmark instances of DJSSP. The performance of IHKA is clarified by comparing it with the original HKA and GAM.
First, the results show that both IHKA and HKA performed very well in solving the 15 benchmark instances of DJSSP, and their results were better than GAM [19]. In 30 independent runs of each benchmark instance of DJSSP, they all showed strong robustness. Their robustness of fitness and running time performed very well. Moreover, their running time was very short (less than 30 s even for large size problems), that is, each algorithm can obtain the best solution to the problem in a reasonable time for each instance of DJSSP. It shows that these two algorithms are suitable to solve DJSSP.
The results show that the convergence rate of IHKA was significantly reduced, which helped to avoid the algorithm falling into the local optimum. That is, the introduced cellular neighbor network played an important role. Furthermore, two different cellular neighborhood structures, Von Neumann neighborhood and Moore neighborhood, are used by IHKA. The results show that different cellular neighborhood structures affect the performance of the algorithm. In the 15 benchmark instances of DJSSP, the neighborhood structure of Moore was not as good as Von Neumann. It shows that the cellular neighborhood structure of the cellular neighbor network introduced by IHKA is worthy of further exploration.
According to the proposed integration-communication protocol, Figure 1

Discussion
The IHKA proposed in this paper was applied to solve the 15 benchmark instances of DJSSP. The performance of IHKA is clarified by comparing it with the original HKA and GAM.
First, the results show that both IHKA and HKA performed very well in solving the 15 benchmark instances of DJSSP, and their results were better than GAM [19]. In 30 independent runs of each benchmark instance of DJSSP, they all showed strong robustness. Their robustness of fitness and running time performed very well. Moreover, their running time was very short (less than 30 s even for large size problems), that is, each algorithm can obtain the best solution to the problem in a reasonable time for each instance of DJSSP. It shows that these two algorithms are suitable to solve DJSSP.
The results show that the convergence rate of IHKA was significantly reduced, which helped to avoid the algorithm falling into the local optimum. That is, the introduced cellular neighbor network played an important role. Furthermore, two different cellular neighborhood structures, Von Neumann neighborhood and Moore neighborhood, are used by IHKA. The results show that different cellular neighborhood structures affect the performance of the algorithm. In the 15 benchmark instances of DJSSP, the neighborhood structure of Moore was not as good as Von Neumann. It shows that the cellular neighborhood structure of the cellular neighbor network introduced by IHKA is worthy of further exploration.
According to the proposed integration-communication protocol, Figure 1, a simulation model was built in the conventional Simio simulation environment, made for evaluation of the 6 × 5 DJSSP dataset. The integration-communication protocol allows the use of the scheduling results of the IHKA in the decision logic of the simulation model. The results prove the correspondence of numerical (machine process time, orders' flow time and makespan) and graphical optimization results (Gantt chart of IHKA and simulation model) with the results of the simulation model. These results prove how important the use of evolutionary computation methods is in solving NP-hard optimization, since the production parameters (machine utilization, orders' flow time) proved the highly optimized degree of DJSSP scheduling with the IHKA, although dynamic events appeared. The proven successful integration-communication protocol affects significantly the high ability to use the proposed decision-making logic in a real-world production environment, which distinguishes the presented research work, along with a new own IHKA DJSSP decision method, from the research work of other researchers. The use of a conventional simulation tool enables the smooth adaptation of the simulation model structure according to the individually analyzed system, while the proposed decision logic and integrationcommunication protocol remain sustainably applicable.
Given that the presented research work addresses the DJSSP optimization problem, one of its theoretical limitations is the pre-known sequence of machines assigned to perform individual operations. The transition from single-objective optimization (current state) to DFJSSP, and, thus, to multi-objective decision-making, would enable the elimination of this limitation, which has now been partially eliminated within the simulation model.

Conclusions
This paper proposes IHKA by introducing the cellular neighbor network. The cellular neighbor network consists of the best positions experienced by individuals in the population. A boundary constraint handling function is also introduced to enhance the diversity of the population. An encoding method based on relative position index was introduced, so that IHKA can be applied to solve the DJSSP. Both IHKA and HKA were used to solve the 15 benchmark instances of DJSSP, and their performances clarified. The introduction of the cellular neighbor network reduced the convergence rate of HKA. Moreover, the use of multiple cellular neighborhood structures can prevent IHKA from falling into the local optimum. However, further research is needed to improve the performance of IHKA using a cellular neighbor structure.
In conclusion, we find that the newly developed decision-making method allows a satisfactorily good solution for the DJSSP optimization problem. It complements its application usability with an efficient integration-communication protocol, which transfers the optimization results to a conventional simulation environment. A comprehensive treatment of the DJSSP optimization problem allows the application of the proposed methods directly in a real-world production environment.
In the following, it would be useful to evaluate the performance of methods on datasets obtained from real-world production systems, and, thus, examine the capabilities of dynamic decision-making in real-time. Given the limitations of the DJSSP optimization problem, the trend of further research will focus on extending the proposed single-criterion methods to multi-criteria methods (DFJSSP consideration), since then IHKA will also be able to determine machine sequence according to a possible set of machines to perform an individual operation. This development step will contribute to an even higher level of applicability of the proposed method, as, in the real-world environment, decisions are often made as to which machine is the most suitable for and most efficient at performing a specific operation.