# An Improved Simulated Annealing-Based Decision Model for the Hybrid Flow Shop Scheduling of Aviation Ordnance Handling

Authors: Xianglei Meng; Nengjian Wang; Jue Liu; Qinhui Liu
Journal: Computational Intelligence and Neuroscience
 (2022)

---

## Abstract

Aviation ordnance handling is critical to the firepower projection of the time-critical cyclic flight operation on aircraft carriers. The complexity of the problem depends on the supply and demand features of ordnance. This paper examines the scheduling of aviation ordnance handling of an operational aircraft carrier under the framework of hybrid flow shop scheduling (HFS) and derives a method based on the simulated annealing (SA) algorithm to get the HFS problem’s solution. The proposed method achieves the minimum possible flow time by optimizing the ordnance assignment through different stages. The traditional SA algorithm depends heavily on the heuristic scheme and consumes too much time to compute the optimal solution. To solve the problem, this paper improves the SA by embedding a task-based encoding method and a matrix perturbation method. The improved SA remains independent of the heuristic scheme and effectively propagates the local search process. Since the performance of SA is also influenced by its embedded parameters, orthogonal tests were carried out to carefully compare and select these parameters. Finally, different ordnance loading plans were simulated to reveal the advantage of the improved SA. The simulation results show that the improved SA (ISA) can generate better and faster solution than the traditional SA. This research provides a practical solution to stochastic HFS problems.


---

## Body

## 1. Introduction

This study focuses on the ordnance dispatching scheduling problem observed onboard the aircraft carrier flight operation, which plays an important role in the air wing firepower projection in its sortie generation [1]. The ordnance handling process involves many stages, equipment, and hundreds of personnel operating in a limited work space [2]; finding an optimal dispatch scheduling for a given ordnance load plan, plus the time critical nature of cyclic flight window requires, is a challenging problem. Traditionally, ordnance dispatching scheduling is made by a human operator’s hand in a spreadsheet with experience, which is always nonoptimal, or even leading to delays that left aircraft launching without firepower. Thus, robust optimal scheduling is essential to conduct the ordnance handling procedure. However, such problem has seldomly been studied, which can be casted in the hybrid flow shop (HFS) scheduling framework.The HFS scheduling problem [3] can be regarded as the combination of the flow shop scheduling (FSS) [4] problem and parallel machine scheduling (PMS) problem, where the former is to decide the job sequences through the shop and the latter is to allocate jobs to machines, given the processing times of each job on each machine according to one or several given criteria, aiming to minimize the makespan [5–7]. For an n jobs m stages problem, there are a total of (n!)m possible schedules, which proves to be NP-hard [8]. If the numbers n of jobs and m of stages are very small, the optimal schedule may be determined by exhaustion, such as branch and bound (B&B) [9] or integer programming techniques [10]. However, these approaches are not applicable to HFS problems with numerous jobs and stages, due to their enormous computing time and memory occupation.Thus, for scheduling different HFS configurations, a large number of approximation and heuristics algorithms have been proposed [11, 12]. The computational complexity of HFS propelled scholars to develop many heuristics algorithms to obtain good enough solutions in a short time for medium-to-large problems, such as different scheduling rules [13], but the heuristic methods are too problem-specific, it often cannot be applied to generalized problem. For the past decades, many general schemes on improving the performance of simple heuristics have been successfully developed, most of which are named as metaheuristics, such as genetic algorithm (GA) [14], ant colony optimization (ACO) [15], tabu search (TS) [16], neural networks (NN) [17], artificial immune systems (AIS) [18], and simulated annealing (SA) [19]. They inhere with higher level of abilities in searching the vast solution space, which have better performance than the simple heuristic methods.Since different heuristics work effectively for different problems, when it encounters the flow shop scheduling problem, Maaroju [20] tested all the metaheuristic methods and found that the genetic algorithm and simulated annealing outperformed others, for hill climbing, swarm intelligence, and neural networks yielded only marginal improvements. However, the computation time for the GA is larger than that for SA. Thus, the SA-based algorithm is chosen in dealing with the ordnance handling problem under the HFS framework. Simulated annealing origins from the metallurgy technology, where a material cools down from high temperature to get minimum energy state. In the algorithm, the current state s and neighbor states s’ are considered, and the algorithm decides the state transition probability from s to s’ based on current system energy (known as temperature). This process continues until a good enough state has been found or the computation threshold has been reached. Such mechanism guarantees approximating to global optimum without getting stuck in local minimum for solving large complex optimization problems. However, the traditional SA algorithm has several defects [21], which include heuristic-dependent, parameter-specific, and long computation time; thus, the performance of the algorithm is yet to be improved. To overcome the above defects, this paper presents a Monte Carlo [22] perturbation method, which directly perturbs the solution matrix in each iteration of SA cooling, eliminating the dependence on any heuristic method, whereas SA performance also depends on cooling parameters; this paper carefully plans the calibration of these parameters to accelerate the computation process by adding double thresholds and setting the memory method of the SA.The organizations of this paper are as follows: Section2 introduces the ordnance handling process in detail; Section 3 discusses the methodology of the SA and improves the SA by embedding a new decoding method and matrix perturbation method; Section 4 evaluates the improved SA algorithm through computational experiments; Section 5 summarizes the research findings and gives the directions of future research.


## 2. Ordnance Handling Procedure

The ordnance handling procedure is specified in a daily loading plan, which lists the amount and types of weapons (throughout this paper, ordnance and weapon are used interchangeably) to be loaded onto the corresponding aircraft. Figure1 shows the layout of aircraft carrier decks, where the construction and transfer of ordnance origin from the magazines located in lower decks to the awaiting aircraft on flight deck following a series of stages. In stage I, the ordnances are retrieved from magazines by bomb skids and delivered to lower-stage elevators. In stage II, the ordnances are lifted to the hangar deck by lower-stage elevators. In stage III, the ordnances are transferred to the staging area of the hangar deck, assembled in that area, and moved to upper-stage elevators. In stage IV, the ordnances are transferred to the flight deck by upper-stage elevators. In stage V, the ordnances are moved directly to and loaded on the aircraft waiting on the flight deck. The flowchart of this procedure is shown in Figure 2. For a common aircraft carrier layout (Figure 1), there are at least 4 magazines in the delivering stage, 8 elevators in the lower-lifting stage, 2 assemble centers in the assembling stage, 4 elevators in the upper-lifting stage, and around 10 aircraft spots in the loading stage. According to Gupta [8], the two-stage flow shop problem with one stage containing a single machine can be NP-hard. Thus, the ordnance handling problem is far from trivial, especially for ordnance officers in making timely decisions of the flight operations.Figure 1
The layout of the ordnance handling routes aboard an aircraft carrier.Figure 2
Flowchart of ordnance handling process. Note:M = magazine; LE = lower-stage elevator; AD = assembly department; UE = upper-stage elevator; P = aircraft parking spot.This paper examines the ordnance handling problem in the HFS framework. The definition of hybrid flow shop system is as follows: in a factory, the set ofn jobs J=1,2,…,j,…,n is going to be processed through m stages M=1,2,…,i,…,m in sequence, while each stage i contains Mi=1,2,…,k,…,mi identical machines, and the processing time of job j on machine k is pjk≥0. The objective is always to minimize makespan. Similarly, in the ordnance handling problem, n batches of weapons are considered as jobs; they also shall be processed in the same order through m stages by facility k (as shown in Figure 2) with processing time pjk, and the objective is to decide the weapons’ sequences and the allocations of weapons to facilities to get the minimum flow time. This is a combinatorial optimization problem with n!m possible schedules, which is considered as NP-hard so that it is difficult to find the optimal solution in polynomial time. For a simple case of 10 batches of weapons in our problem, there can be 10!5=6.3×1032 different schedules for the ordnance officer to choose, which is beyond human mind’s reach in conducting the time-critical flight operations.The ordnance configuration of aircraft on the carrier depends on the specific mission [23]. It is assumed that each transfer equipment or personnel only transfers one type of ordnance at a time [24]. For each type of ordnance, the number loaded in one skid is denoted as rw Thus, the batches of ordnances needed to complete the task of all aircraft can be determined by(1)task=11,num12,num2...w,numw21,num12,num2...w,numw...W,numW...a1,num12,num2...w,numw...W,numW,w=1,2,…W,where a is the number of aircraft to be loaded, w is the type of ordnance, and numw is the actual number of ordnance types w. (2)ξ=∑i=1a∑j=1wnumwrw,i=1,2,…a,j=1,2,…w,where rw is the number of type w ordnances in one skid, a is the number of aircraft, and w is the number of ordnance types to be loaded on the aircraft. That is, an ordinance should be transferred to the required aircraft, once being retrieved from the magazine.


### 2.1. Stage I: Weapons Retrieving

Multiple magazines are located in the bow and aft of the carrier, and the ordnances can be transferred by multiple elevators. The ordnances are firstly retrieved from magazines by skids with the setup timeT0. From the same magazine, the ordnances should be retrieved with an interval no shorter than tint. Then, the skids deliver the ordnances to lower-stage elevators. The time consumed to transfer ordnances from magazines to lower-stage elevators can be expressed as(3)TML=t11t12⋯t1Lt21t22⋯t2L⋮⋮⋮⋮tM1tM2⋯tML,where M is a magazine and L is a lower-stage elevator.


### 2.2. Stage II: Weapons Buildup

The ordnances are loaded onto lower-stage elevators and lifted vertically with constant speed to the hangar deck. The time consumed in this stage (lifting time) is denoted asTL.


### 2.3. Stage III: Weapons Assembling

The ordnances are preassembled in the staging area of the hangar deck, with a sufficient lead time to meet the short turnaround time of the flight schedule. The assembling timeTKass varies with the types of ordnances. Note that the assembling time of the staff fluctuates in the real world. Therefore, the interval of assembling time was set to [−Tf1, Tf1]. The real assembling time is denoted as TKass+t1, where t1 is a random number within [−Tf1, Tf1].


### 2.4. Stage IV: Weapons Striking Up

The ordnances are transferred to the flight deck by upper-stage elevators. The time consumed in this stage (transport time) is denoted asTU.


### 2.5. Stage V: Weapons Loading

Some ordnance crew members on the flight deck transport the ordnances from the upper-stage elevators to the aircraft. The time consumed in this stage can be expressed as(4)TUA=t11t12⋯t1At21t22⋯t2A⋮⋮⋮⋮tU1tU2⋯tUA,where U is an upper-stage elevator and A is an aircraft.The other ordnance crew members load the ordnances onto the aircraft. It is assumed that the different groups consume the same time to load the same ordnance and different types of ordnances need different time to be loaded. The time needed to load each type of ordnance is denoted asTKload. The interval of the loading time was set as [−Tf2, Tf2]. The real loading time is denoted as TKload+t2, where t2 is a random number within [−Tf2, Tf2].The ordnance crew members can be shared across groups. The loading cannot proceed unless ordnance crew members are available. The walking time for interstation transfer between different aircraft can be expressed as(5)TAA=t11t12⋯t1At21t22⋯t2A⋮⋮⋮⋮tA1tA2⋯tAA.Referring to the standard three-field notation for scheduling problems, our problem can be described as follows:FH5,PMkk2=1,RMkk=34|prmu,Mj5,block|Cmax. Specifically, FH5 is a five-stage HFS problem: stage I involves M(1) identical magazines that store ordnances; stage II involves M(2) lower-stage elevators to transport the ordnances; stage III has M(3) identical assembling personnel to assemble the ordnances; stage IV has RM(4) independent upper-stage elevators to transport the ordnances; and stage V has RM(5) independent aircraft to be loaded. Note that prmu indicates that the ordnances are handled in the same order in every stage; Mj5 (eligibility constraint) means that the handling of ordnance j is limited to the aircraft set M in stage V; block indicates that the capacity of buffer between stages is constrained, for instance, the weapons have to wait in the current stage till enough room is released for the next stage of handling.In total, the completion time for batchi of ordnances can be calculated by(6)Ci=TML+TL+TKass+t1+TU+TUA+TAA+TKload+t2+Twating,i=1,2,…ξ,where Twating is the whole waiting time in the transporting process, for an ordinance cannot be handled unless machines or ordnance crew members are available.The general objective of ordnance handling is to complete all transporting operations as efficiently as possible within the specified time and to generate a reasonable schedule for ordnance handling. Therefore, for our ordnance handling problem, the minimization of makespan is set as goal, so that sufficient ordnances can be loaded to the awaiting aircraft to fly in the next flying window.(7)Object=minCmax.


## 2.1. Stage I: Weapons Retrieving

Multiple magazines are located in the bow and aft of the carrier, and the ordnances can be transferred by multiple elevators. The ordnances are firstly retrieved from magazines by skids with the setup timeT0. From the same magazine, the ordnances should be retrieved with an interval no shorter than tint. Then, the skids deliver the ordnances to lower-stage elevators. The time consumed to transfer ordnances from magazines to lower-stage elevators can be expressed as(3)TML=t11t12⋯t1Lt21t22⋯t2L⋮⋮⋮⋮tM1tM2⋯tML,where M is a magazine and L is a lower-stage elevator.


## 2.2. Stage II: Weapons Buildup

The ordnances are loaded onto lower-stage elevators and lifted vertically with constant speed to the hangar deck. The time consumed in this stage (lifting time) is denoted asTL.


## 2.3. Stage III: Weapons Assembling

The ordnances are preassembled in the staging area of the hangar deck, with a sufficient lead time to meet the short turnaround time of the flight schedule. The assembling timeTKass varies with the types of ordnances. Note that the assembling time of the staff fluctuates in the real world. Therefore, the interval of assembling time was set to [−Tf1, Tf1]. The real assembling time is denoted as TKass+t1, where t1 is a random number within [−Tf1, Tf1].


## 2.4. Stage IV: Weapons Striking Up

The ordnances are transferred to the flight deck by upper-stage elevators. The time consumed in this stage (transport time) is denoted asTU.


## 2.5. Stage V: Weapons Loading

Some ordnance crew members on the flight deck transport the ordnances from the upper-stage elevators to the aircraft. The time consumed in this stage can be expressed as(4)TUA=t11t12⋯t1At21t22⋯t2A⋮⋮⋮⋮tU1tU2⋯tUA,where U is an upper-stage elevator and A is an aircraft.The other ordnance crew members load the ordnances onto the aircraft. It is assumed that the different groups consume the same time to load the same ordnance and different types of ordnances need different time to be loaded. The time needed to load each type of ordnance is denoted asTKload. The interval of the loading time was set as [−Tf2, Tf2]. The real loading time is denoted as TKload+t2, where t2 is a random number within [−Tf2, Tf2].The ordnance crew members can be shared across groups. The loading cannot proceed unless ordnance crew members are available. The walking time for interstation transfer between different aircraft can be expressed as(5)TAA=t11t12⋯t1At21t22⋯t2A⋮⋮⋮⋮tA1tA2⋯tAA.Referring to the standard three-field notation for scheduling problems, our problem can be described as follows:FH5,PMkk2=1,RMkk=34|prmu,Mj5,block|Cmax. Specifically, FH5 is a five-stage HFS problem: stage I involves M(1) identical magazines that store ordnances; stage II involves M(2) lower-stage elevators to transport the ordnances; stage III has M(3) identical assembling personnel to assemble the ordnances; stage IV has RM(4) independent upper-stage elevators to transport the ordnances; and stage V has RM(5) independent aircraft to be loaded. Note that prmu indicates that the ordnances are handled in the same order in every stage; Mj5 (eligibility constraint) means that the handling of ordnance j is limited to the aircraft set M in stage V; block indicates that the capacity of buffer between stages is constrained, for instance, the weapons have to wait in the current stage till enough room is released for the next stage of handling.In total, the completion time for batchi of ordnances can be calculated by(6)Ci=TML+TL+TKass+t1+TU+TUA+TAA+TKload+t2+Twating,i=1,2,…ξ,where Twating is the whole waiting time in the transporting process, for an ordinance cannot be handled unless machines or ordnance crew members are available.The general objective of ordnance handling is to complete all transporting operations as efficiently as possible within the specified time and to generate a reasonable schedule for ordnance handling. Therefore, for our ordnance handling problem, the minimization of makespan is set as goal, so that sufficient ordnances can be loaded to the awaiting aircraft to fly in the next flying window.(7)Object=minCmax.


## 3. Improved Simulated Annealing Algorithm

The SA is a technique capable of searching for good solutions to various combinatorial problems in material science and physics. The pseudocode of the algorithm is as follows.Algorithm 1: Standard simulated annealing.
Input: A initialS solution and a cost value function F(x).Output: AS′ solution that minimizes the cost value function F(x).T ← initializing Temperature//a method for assigning an initial temperaturewhileT > B not freezing do//a definition of “frozen,”fori = 1 to C do while not at equilibrium do//a definition of “equilibrium,”S’← new permutation of S.IfF(S′) < F(S) or Random value < eFS’−FS/DT then S ← S’//a selection criterionendT ← reduced temperature//a way of calculating the next temperatureendThe SA includes three major functions: state generation, state acceptance, and temperature update. The first is to make perturbations of the given initial solution, in order to search for the optimal solution effectively in the vast solution space. The second decides whether to accept a newly generated solution with a certain probability in case of trapping in the local minimum. The third offers a cooling scheme that mimics the physical annealing process to get a stable state for the problem. The SA performance can be augmented by adjusting various parameters and operators [25], such as initial temperature, descent gradient of temperature, and termination rule, which have to be adjusted manually. This paper mainly improves the search efficiency (timeliness) of the SA, without sacrificing the optimization quality. Thus, encoding scheme, initial solution, neighborhood search structure (NSS), type of cooling schedule, and two criteria (internal cycle termination and external cycle termination) are modified as described in the following subsections.


### 3.1. Task Encoding

This paper presents a task-based encoding method for the ordnance handling problem, where ordnances are coupled with aircraft. Each job of the HFS is defined as an operationoxy, where x is the type of ordnance and y is the aircraft to be loaded. Then, a solution can be expressed as(8)S=o1xyo2xym11…m21…mn1mn1⋮⋮⋮oixymn1⋯mnn,where the first column is a permutation of the task sequence, the following columns are the corresponding stages, and m1 to mn are machines assigned randomly to execute the tasks.


### 3.2. Initial Solution

As mentioned above, the initial solution has a great impact on the final solution of the SA. According to the heuristic methods mentioned in [26], we generate the initial solution by the following rule: first, the permutation of tasks is determined by assigning each type of weapon x to its corresponding aircraft y, defined as operation oixy, according to the ordnance loading plan P. Then, the tasks are assigned to the earliest available machine. If there are more than one earliest available machines, one of them will be chosen randomly. The initial solution will be generated as Si×m with i tasks through m stages.For example, Table1 shows the ordnance loading plan of the aircraft waiting to operate in the next fly window. There are three aircraft and four types of ordnances. Aircraft 1 requires 2 skids of type 1 ordnances and 1 skid of type 2 ordnances.Table 1
Ordnance loading plan.
Aircraft type of ordnance123121—21——3—1—4——1Here, the operations are permutated by aircraft number (equal priority) in ascending order:oixy=o111,o211,o321,o412,o532,o643. First, a random permutation of jobs (tasks) is generated. Then, the optimal available elevator is assigned to each operation. The final plan can be written as(9)311111522222133133444241211112622223.The meaning of such matrix can be explained as follows: taken the third row for illustration, operationo111 is extracted out of the third magazine, transferred by the third lower-stage elevator to the first staging area, moved by the third upper-stage elevator to the flight deck, and loaded to the corresponding aircraft by the No. 3 ordnance crew.


### 3.3. Neighborhood Perturbation

The NSS can generate a new solution by slightly modifying the current candidate solution. Traditionally, many different NSSs are adopted in each iteration of SA computation, such as swap, shift, and reversion [27], to realize the task permutation in the first stage. Then, the same heuristic rules are applied in the following stages to generate a new solution S’. The quality of the solution depends on the selected heuristic rules. If the rules are too greedy, the algorithm may fall into the local optimum, being unable to get convergence to the final optimal result.The quality of the SA solution is highly sensitive to the selection of candidate solutions. Therefore, the perturbation scheme is crucial to the good performance of the SA algorithm. To ease the dependence of neighborhood search on heuristics, this paper proposes a Monte Carlo perturbation technique, which directly perturbs the initial solution matrix. The initial solution is changed significantly in one step, eliminating the effect by heuristic methods. The matrix perturbation is described as follows.Algorithm 2: Matrix perturbation.
(1)
Randomly generate a rectangleRa×b, with width a∈0,i and length b∈0,m;(2)
PutR in the solution matrix S and mark the four vertices oixy,oi+ax,y,oi+ax+b,y+b,oix+b,y+b;(3)
The subsequence of theR area is reversed, oixy⟶oi+ax,y and oix+b,y+b⟶oi+ax+b,y+b;(4)
Get the new solutionS′=oi+axy…oi+ax+b,y+b⋮⋱⋮oixy⋯oix+b,y+b.The following is an example of the matrix perturbation process: for instance, a6×3 matrix, a rectangle Ra×b randomly generated a size of 4×3 matrix.which covers the matrix of 421632213121.Then, matrixR is reversed as R’ =  121213632421.Then, the new solution can be obtained asS’ =  312121213632421533.Hence, the matrix perturbation process is completed:312421632213121533⟶312121213632421533.


### 3.4. Cooling Schedule

The SA behavior can be regulated by the temperature and its descent gradient. To avoid the local optimum trap, inferior solutions may be accepted depending on the falling temperature, under the mechanism of cooling schedule. Here, the exponential cooling rate is adopted:(10)Tl=T0−TfN+1Nl+1+T0−T0−TfN+1N;l=1,2,…,N,where T0 and Tf are initial temperature and final temperature, respectively, and N is the number of temperatures between T0 and Tf.The SA needs to accept the new state through probability judgment, in order to avoid the local minimum. When the initial temperature is sufficiently high, the cooling is slow enough (i.e., each temperature is held for a sufficiently long time), and the final temperature approaches zero; the SA will converge to the global optimal solution with the probability of 1. However, it is very difficult to fulfil the global convergence condition. Besides, the current state may be worse than some intermediate states in the search trajectory, owing to the probability acceptance mechanism. Thus, the SA algorithm often converges to an approximate optimal solution, or a solution poorer than the best intermediate solution. The search efficiency is inevitably affected. To preserve the best-known state and improve search efficiency, this section makes the following improvements to the SA:(1)
Memorize the best intermediate solution in the search process and update it immediately. The improvement of memory turns the SA into an intelligent algorithm.(2)
Set up two thresholds, internal cycle threshold and external cycle threshold, to reduce the computing load while maintaining optimality. The internal cycle threshold refers to the number of cycles that the new solution of continuous disturbance does not generate a better solution at a certain temperature, while the external cycle threshold refers to the number of cycles that the new solution generated by continuous cooling does not generate a better solution. The two thresholds are determined as follows.First, determine whether the number of internal cycles reaches the threshold; if yes, lower the temperature by one step; otherwise, judge if it conforms to Markov chain. If not, reconduct the process of state generation, state acceptance, and algorithm termination; otherwise, lower the temperature by one step. Second, determine whether the number of external cycles reaches the threshold; if yes, terminate the algorithm and obtain the final solution; otherwise, judge whether the algorithm meets the termination conditions. If yes, terminate the algorithm and obtain the final solution; otherwise, reconduct the process of state generation, state acceptance, and algorithm termination. Terminate the algorithm once the number of iterationsi surpasses the prior fixed constant MaxIter. In our experiments, MaxIter was set to 104.The flow of the improved SA is shown in Figure3.Figure 3
Flow of the improved SA.


## 3.1. Task Encoding

This paper presents a task-based encoding method for the ordnance handling problem, where ordnances are coupled with aircraft. Each job of the HFS is defined as an operationoxy, where x is the type of ordnance and y is the aircraft to be loaded. Then, a solution can be expressed as(8)S=o1xyo2xym11…m21…mn1mn1⋮⋮⋮oixymn1⋯mnn,where the first column is a permutation of the task sequence, the following columns are the corresponding stages, and m1 to mn are machines assigned randomly to execute the tasks.


## 3.2. Initial Solution

As mentioned above, the initial solution has a great impact on the final solution of the SA. According to the heuristic methods mentioned in [26], we generate the initial solution by the following rule: first, the permutation of tasks is determined by assigning each type of weapon x to its corresponding aircraft y, defined as operation oixy, according to the ordnance loading plan P. Then, the tasks are assigned to the earliest available machine. If there are more than one earliest available machines, one of them will be chosen randomly. The initial solution will be generated as Si×m with i tasks through m stages.For example, Table1 shows the ordnance loading plan of the aircraft waiting to operate in the next fly window. There are three aircraft and four types of ordnances. Aircraft 1 requires 2 skids of type 1 ordnances and 1 skid of type 2 ordnances.Table 1
Ordnance loading plan.
Aircraft type of ordnance123121—21——3—1—4——1Here, the operations are permutated by aircraft number (equal priority) in ascending order:oixy=o111,o211,o321,o412,o532,o643. First, a random permutation of jobs (tasks) is generated. Then, the optimal available elevator is assigned to each operation. The final plan can be written as(9)311111522222133133444241211112622223.The meaning of such matrix can be explained as follows: taken the third row for illustration, operationo111 is extracted out of the third magazine, transferred by the third lower-stage elevator to the first staging area, moved by the third upper-stage elevator to the flight deck, and loaded to the corresponding aircraft by the No. 3 ordnance crew.


## 3.3. Neighborhood Perturbation

The NSS can generate a new solution by slightly modifying the current candidate solution. Traditionally, many different NSSs are adopted in each iteration of SA computation, such as swap, shift, and reversion [27], to realize the task permutation in the first stage. Then, the same heuristic rules are applied in the following stages to generate a new solution S’. The quality of the solution depends on the selected heuristic rules. If the rules are too greedy, the algorithm may fall into the local optimum, being unable to get convergence to the final optimal result.The quality of the SA solution is highly sensitive to the selection of candidate solutions. Therefore, the perturbation scheme is crucial to the good performance of the SA algorithm. To ease the dependence of neighborhood search on heuristics, this paper proposes a Monte Carlo perturbation technique, which directly perturbs the initial solution matrix. The initial solution is changed significantly in one step, eliminating the effect by heuristic methods. The matrix perturbation is described as follows.Algorithm 2: Matrix perturbation.
(1)
Randomly generate a rectangleRa×b, with width a∈0,i and length b∈0,m;(2)
PutR in the solution matrix S and mark the four vertices oixy,oi+ax,y,oi+ax+b,y+b,oix+b,y+b;(3)
The subsequence of theR area is reversed, oixy⟶oi+ax,y and oix+b,y+b⟶oi+ax+b,y+b;(4)
Get the new solutionS′=oi+axy…oi+ax+b,y+b⋮⋱⋮oixy⋯oix+b,y+b.The following is an example of the matrix perturbation process: for instance, a6×3 matrix, a rectangle Ra×b randomly generated a size of 4×3 matrix.which covers the matrix of 421632213121.Then, matrixR is reversed as R’ =  121213632421.Then, the new solution can be obtained asS’ =  312121213632421533.Hence, the matrix perturbation process is completed:312421632213121533⟶312121213632421533.


## 3.4. Cooling Schedule

The SA behavior can be regulated by the temperature and its descent gradient. To avoid the local optimum trap, inferior solutions may be accepted depending on the falling temperature, under the mechanism of cooling schedule. Here, the exponential cooling rate is adopted:(10)Tl=T0−TfN+1Nl+1+T0−T0−TfN+1N;l=1,2,…,N,where T0 and Tf are initial temperature and final temperature, respectively, and N is the number of temperatures between T0 and Tf.The SA needs to accept the new state through probability judgment, in order to avoid the local minimum. When the initial temperature is sufficiently high, the cooling is slow enough (i.e., each temperature is held for a sufficiently long time), and the final temperature approaches zero; the SA will converge to the global optimal solution with the probability of 1. However, it is very difficult to fulfil the global convergence condition. Besides, the current state may be worse than some intermediate states in the search trajectory, owing to the probability acceptance mechanism. Thus, the SA algorithm often converges to an approximate optimal solution, or a solution poorer than the best intermediate solution. The search efficiency is inevitably affected. To preserve the best-known state and improve search efficiency, this section makes the following improvements to the SA:(1)
Memorize the best intermediate solution in the search process and update it immediately. The improvement of memory turns the SA into an intelligent algorithm.(2)
Set up two thresholds, internal cycle threshold and external cycle threshold, to reduce the computing load while maintaining optimality. The internal cycle threshold refers to the number of cycles that the new solution of continuous disturbance does not generate a better solution at a certain temperature, while the external cycle threshold refers to the number of cycles that the new solution generated by continuous cooling does not generate a better solution. The two thresholds are determined as follows.First, determine whether the number of internal cycles reaches the threshold; if yes, lower the temperature by one step; otherwise, judge if it conforms to Markov chain. If not, reconduct the process of state generation, state acceptance, and algorithm termination; otherwise, lower the temperature by one step. Second, determine whether the number of external cycles reaches the threshold; if yes, terminate the algorithm and obtain the final solution; otherwise, judge whether the algorithm meets the termination conditions. If yes, terminate the algorithm and obtain the final solution; otherwise, reconduct the process of state generation, state acceptance, and algorithm termination. Terminate the algorithm once the number of iterationsi surpasses the prior fixed constant MaxIter. In our experiments, MaxIter was set to 104.The flow of the improved SA is shown in Figure3.Figure 3
Flow of the improved SA.


## 4. Experiments

To test the effectiveness of SA-based algorithm, we first evaluate the control factors of the SA and suggest a good parameter setting. Then, the solution quality and efficiency of the ISA were verified through several experiments. The algorithms are implemented in our previously published carrier-based flight operations simulation [28], which is written by C++ and ran on Microsoft Windows operating system with 4 GB RAM and dual core CPU.


### 4.1. Parameter Tuning

The efficiency of the SA algorithm is greatly affected by the design of parameters and operators. The full factorial design tests all possible combinations. Such an approach becomes too laborious in the face of numerous factors. Taguchi utilized orthogonal arrays to examine lots of decision variables in a few tests [29] and measured the importance of each factor by its influence on algorithm performance, using the signal-to-noise ratio: 10log10objective2. Following Taguchi’s method, the SA control factors were configured as follows: initial solution, initial temperature, cooling rate, and number of neighborhood searches in every temperature. Table 2 shows the different levels of these factors.Table 2
The SA control factors.
FactorsLevelsTypesInitial solution (A)3(A1) randomness (A2) task-based (A3) ascending orderInitial temperature (B)4B1 = 50, B2 = 100, B3 = 500, B4 = 1000Cooling rate (C)4C1 = 0.85, C2 = 0.9, C3 = 0.95, C4 = 0.99Markov chain (d)4D1 = 100, D2 = 500, D3 = 1000, D4 = 2000Hence, the SA has one 3-level factor and three 4-level factors. The best design among the orthogonal arrays is L16. Thus, additional transform was performed to fit L16 (Table3).Table 3
Orthogonal array of L16 of our algorithm test.
ExperimentABCD112332241233434412115131462131731138233291142102323113341123124131421142244153222161443The relative percentage deviation (RPD) was also adopted to measure the performances:(11)RPD=Algsol−MinsolMinsol⋅100%,where the best solution obtained for one instance is denoted as Minsol, while the objective value is marked as Algsol. Table 4 lists the S/N ratios and RPD values of each level of the factor value. The results show that A(3), B(4), C(4), and D(4) are the best levels of the factors.Table 4
Results of the orthogonal test.
Factor with levelCompletion timeMean S/N ratioMean RPDRandom initial (A1)−73.67792.399Descending initial (A2)−73.66342.216Ascending initial∗ (A3)−73.6504∗2.060∗T0 = 50 (B1)−73.66562.239T0 = 100 (B2)−73.66382.223T0 = 500 (B3)−73.66692.252T0 = 1000∗ (B4)−73.6629∗2.212∗Alpha = 0.85 (C1)−73.72662.962Alpha = 0.9 (C2)−73.67912.397Alpha = 0.95 (C3)−73.65472.113Alpha = 0.99∗ (C4)−73.5987∗1.453∗L = 100 (D1)−73.76023.360L = 500 (D2)−73.67162.310L = 1000 (D3)−73.63441.869L = 2000∗ (D4)−73.5930∗1.387∗


### 4.2. Experimental Settings

The size of the test instances was set ton=15,30,45,60 tasks, which corresponds to a common mission of strike sorties of 5, 10, 15, and 20 aircraft, respectively. The processing time for jobs on each machine was generated by triangle distribution with mean time according to [30]. There are resource constraints of five types of weapons, four lower-stage elevators, 10 assembling crews, four upper-stage elevators, and six loading crew members, see Tables 5–9 for comprehensive data.Table 5
Ordnance loading plan.
Ordnance typeSpot taskA1A2A3A4A5A6A7A8A9A1014(2)24(2)32(1)2(1)2(1)2(1)2(1)2(1)2(1)2(1)42(1)2(1)2(1)2(1)2(1)2(1)2(1)2(1)56(2)6(2)6(2)6(2)62(2)2(2)2(2)2(2)Table 6
Ordnance assembling time (min).
Ordnance type123456Loading time/skid13877116Table 7
Ordnance transport time on flight deck.
Time (s)Upper-stageelevator numberSpotA1A2A3A4A5A6A7A8A9A1012451991483061107214213402342210660190926667961593452863381337262166127787491267208444339832022818913411880217161Table 8
Ordnance loading time (s).
Ordnance type123456Time (skid/s)660600480480780360Table 9
Ordnance crew walking time (s).
SpotSpotA1A2A3A4A5A6A7A8A9A10A1047.53946168.2222.6266.9331366446.9639578.8A247.50142.9181.2225.6289.4324.9404.6596.4536.1A3168.2142.9097.8131.5189.8220302.5492432.9A4222.6181.297.8044.4108.5143.7224.3416.5356.3A5266.9225.6131.544.4064.199.3180.1372.3312.1A6331289.4189.9108.564.1036.1116308.2278.2A7366324.9220143.799.336.1082.6274214.1A8446.9404.6302.5224.3180.111682.60192.2132A9639596.4492416.5372.3308.1623274.0073192.1666060.3A10578.8536.1432.9356.3312.1278.2214.113260.30


### 4.3. Makespan Analysis

The test of problem uses an ordnance loading plan, ranging from 5 to 20 aircraft. At first, the experiment tries to solve a standard small size problem with optimal solution and preliminarily demonstrates the adaptability and feasibility of the ISA. Then, the proposed ISA was adopted to solve larger size problems and compared with the other methods to reveal its superiority.To compare the ISA with the SA, control factors were configured as those in the preceding section. The permutation of machines was arranged in ascending order. The initial temperatureT0 = 103, final temperature ε = 0.5, cooling rate α = 0.99, and length of Markov chain L = 2000.The base case includes 15 tasks. As the temperature declined (Figure4(a)), the ISA converged to the optimal solution in 2.705 s, faster than the SA, which converged to the optimal solution in 2.965 s (Figure 5(a)). Although the initial scheduling time of the ISA was higher than that of SA, its faster convergence procedure suggests the good performance of the ISA. Note that the computing time of the ISA was 85.3637 s, much shorter than that (143.7861 s) of the SA.Figure 4
Optimal schedule of 15 tasks (5 aircraft) derived by the ISA. (a) Convergence curve. (b) Gantt chart.
(a)(b)Figure 5
Optimal schedule of 15 tasks (5 aircraft) derived by the SA. (a) Convergence curve. (b) Gantt chart.
(a)(b)Next, the ISA performance on larger problems was tested, whereas more tasks bring a greater computing effort for the heuristic. Table10 lists the average scheduling time and computing time of different tasks, with each task running for 30 times. Figure 6 shows the corresponding plot. The results showed that the ISA achieved a shorter scheduling time and a faster computing speed than the SA. Comparing with the SA, the average scheduling time of the ISA reduces from 300 to 600 seconds in each ordnance turnaround cycle for 30 to 60 tasks and the average computing time saves from 110 to 190 seconds. Note that there are usually ten or more cycles in a general flight day; the scheduling time saved by implementing ISA equals one more group sortie generation, which in turn enhanced the firepower capacity of carrier air wing.Table 10
Comparison between the ISA and the SA.
Task numbersAverage scheduling time (s)Average computing time (s)ISASAISASA152965310060150304200450015026045610065002003706080008600260450Figure 6
Performance of the ISA and the SA on medium and large problems. (a) Scheduling time. (b) Computing time.
(a)(b)The results also show that the bottleneck of ordnance handling is the loading process, where the number of loading crews heavily influences the aircraft turnaround time. When there is high intensity of surge operations, more loading crews should be arranged to handle ordnances.


## 4.1. Parameter Tuning

The efficiency of the SA algorithm is greatly affected by the design of parameters and operators. The full factorial design tests all possible combinations. Such an approach becomes too laborious in the face of numerous factors. Taguchi utilized orthogonal arrays to examine lots of decision variables in a few tests [29] and measured the importance of each factor by its influence on algorithm performance, using the signal-to-noise ratio: 10log10objective2. Following Taguchi’s method, the SA control factors were configured as follows: initial solution, initial temperature, cooling rate, and number of neighborhood searches in every temperature. Table 2 shows the different levels of these factors.Table 2
The SA control factors.
FactorsLevelsTypesInitial solution (A)3(A1) randomness (A2) task-based (A3) ascending orderInitial temperature (B)4B1 = 50, B2 = 100, B3 = 500, B4 = 1000Cooling rate (C)4C1 = 0.85, C2 = 0.9, C3 = 0.95, C4 = 0.99Markov chain (d)4D1 = 100, D2 = 500, D3 = 1000, D4 = 2000Hence, the SA has one 3-level factor and three 4-level factors. The best design among the orthogonal arrays is L16. Thus, additional transform was performed to fit L16 (Table3).Table 3
Orthogonal array of L16 of our algorithm test.
ExperimentABCD112332241233434412115131462131731138233291142102323113341123124131421142244153222161443The relative percentage deviation (RPD) was also adopted to measure the performances:(11)RPD=Algsol−MinsolMinsol⋅100%,where the best solution obtained for one instance is denoted as Minsol, while the objective value is marked as Algsol. Table 4 lists the S/N ratios and RPD values of each level of the factor value. The results show that A(3), B(4), C(4), and D(4) are the best levels of the factors.Table 4
Results of the orthogonal test.
Factor with levelCompletion timeMean S/N ratioMean RPDRandom initial (A1)−73.67792.399Descending initial (A2)−73.66342.216Ascending initial∗ (A3)−73.6504∗2.060∗T0 = 50 (B1)−73.66562.239T0 = 100 (B2)−73.66382.223T0 = 500 (B3)−73.66692.252T0 = 1000∗ (B4)−73.6629∗2.212∗Alpha = 0.85 (C1)−73.72662.962Alpha = 0.9 (C2)−73.67912.397Alpha = 0.95 (C3)−73.65472.113Alpha = 0.99∗ (C4)−73.5987∗1.453∗L = 100 (D1)−73.76023.360L = 500 (D2)−73.67162.310L = 1000 (D3)−73.63441.869L = 2000∗ (D4)−73.5930∗1.387∗


## 4.2. Experimental Settings

The size of the test instances was set ton=15,30,45,60 tasks, which corresponds to a common mission of strike sorties of 5, 10, 15, and 20 aircraft, respectively. The processing time for jobs on each machine was generated by triangle distribution with mean time according to [30]. There are resource constraints of five types of weapons, four lower-stage elevators, 10 assembling crews, four upper-stage elevators, and six loading crew members, see Tables 5–9 for comprehensive data.Table 5
Ordnance loading plan.
Ordnance typeSpot taskA1A2A3A4A5A6A7A8A9A1014(2)24(2)32(1)2(1)2(1)2(1)2(1)2(1)2(1)2(1)42(1)2(1)2(1)2(1)2(1)2(1)2(1)2(1)56(2)6(2)6(2)6(2)62(2)2(2)2(2)2(2)Table 6
Ordnance assembling time (min).
Ordnance type123456Loading time/skid13877116Table 7
Ordnance transport time on flight deck.
Time (s)Upper-stageelevator numberSpotA1A2A3A4A5A6A7A8A9A1012451991483061107214213402342210660190926667961593452863381337262166127787491267208444339832022818913411880217161Table 8
Ordnance loading time (s).
Ordnance type123456Time (skid/s)660600480480780360Table 9
Ordnance crew walking time (s).
SpotSpotA1A2A3A4A5A6A7A8A9A10A1047.53946168.2222.6266.9331366446.9639578.8A247.50142.9181.2225.6289.4324.9404.6596.4536.1A3168.2142.9097.8131.5189.8220302.5492432.9A4222.6181.297.8044.4108.5143.7224.3416.5356.3A5266.9225.6131.544.4064.199.3180.1372.3312.1A6331289.4189.9108.564.1036.1116308.2278.2A7366324.9220143.799.336.1082.6274214.1A8446.9404.6302.5224.3180.111682.60192.2132A9639596.4492416.5372.3308.1623274.0073192.1666060.3A10578.8536.1432.9356.3312.1278.2214.113260.30


## 4.3. Makespan Analysis

The test of problem uses an ordnance loading plan, ranging from 5 to 20 aircraft. At first, the experiment tries to solve a standard small size problem with optimal solution and preliminarily demonstrates the adaptability and feasibility of the ISA. Then, the proposed ISA was adopted to solve larger size problems and compared with the other methods to reveal its superiority.To compare the ISA with the SA, control factors were configured as those in the preceding section. The permutation of machines was arranged in ascending order. The initial temperatureT0 = 103, final temperature ε = 0.5, cooling rate α = 0.99, and length of Markov chain L = 2000.The base case includes 15 tasks. As the temperature declined (Figure4(a)), the ISA converged to the optimal solution in 2.705 s, faster than the SA, which converged to the optimal solution in 2.965 s (Figure 5(a)). Although the initial scheduling time of the ISA was higher than that of SA, its faster convergence procedure suggests the good performance of the ISA. Note that the computing time of the ISA was 85.3637 s, much shorter than that (143.7861 s) of the SA.Figure 4
Optimal schedule of 15 tasks (5 aircraft) derived by the ISA. (a) Convergence curve. (b) Gantt chart.
(a)(b)Figure 5
Optimal schedule of 15 tasks (5 aircraft) derived by the SA. (a) Convergence curve. (b) Gantt chart.
(a)(b)Next, the ISA performance on larger problems was tested, whereas more tasks bring a greater computing effort for the heuristic. Table10 lists the average scheduling time and computing time of different tasks, with each task running for 30 times. Figure 6 shows the corresponding plot. The results showed that the ISA achieved a shorter scheduling time and a faster computing speed than the SA. Comparing with the SA, the average scheduling time of the ISA reduces from 300 to 600 seconds in each ordnance turnaround cycle for 30 to 60 tasks and the average computing time saves from 110 to 190 seconds. Note that there are usually ten or more cycles in a general flight day; the scheduling time saved by implementing ISA equals one more group sortie generation, which in turn enhanced the firepower capacity of carrier air wing.Table 10
Comparison between the ISA and the SA.
Task numbersAverage scheduling time (s)Average computing time (s)ISASAISASA152965310060150304200450015026045610065002003706080008600260450Figure 6
Performance of the ISA and the SA on medium and large problems. (a) Scheduling time. (b) Computing time.
(a)(b)The results also show that the bottleneck of ordnance handling is the loading process, where the number of loading crews heavily influences the aircraft turnaround time. When there is high intensity of surge operations, more loading crews should be arranged to handle ordnances.


## 5. Conclusions

This paper treats the aviation ordnance scheduling problem under the HFS framework with multistages, independent parallel machines, and several processing constraints. The simulated annealing algorithm was modified with dual threshold selection to generate faster and better schedules using the proposed matrix perturbation method that keeps the SA independent of the heuristic schemes. The influencing parameters of the improved algorithm are carefully tuned by Taguchi’s method. The experimental results demonstrate the effectiveness of ISA, which provides a practical solution to a broad application in dealing with stochastic hybrid flow shop scheduling problems.



---

*Source: 1843675-2022-02-02.xml*