## 2 Discrete-Event-Simulation-Based Method
### An Improved NEH Algorithm  
The basic steps of the classical NEH algorithm can be described as follows:  
(1) Sort the n workpieces in descending order of the average processing time of each job on m machines to obtain an initial sequence of these n workpieces to be processed. (2) Select the first k \(-\) 1 jobs from the initial sequence as the current sequence (k starts at 3). Insert the ith job of the current sequence (i from k to n) into each of the k possible positions of the current sequence to obtain a sequence of k partial jobs and calculate the maximum completion time corresponding to these k partial sequences. The idea of generating a sequence of k partial jobs in this step is shown in Figure 4.  
The improvement method of the NEH algorithm generally lies in two points: the priority rule and the tie-breaking mechanism. In the reported literature, NEHIJP1 is the most effective NEH-based algorithm, and PRLJP1 and TBLJP1 contained in NEHIJP1 have been proven to be the best priority rule and tie-breaking mechanism, respectively. Although we tried several new tie-breaking mechanisms, their performance never exceeded that of TBLJP1. As a result, we simply proposed a new priority rule and used TBLJP1. The new improved NEH algorithm presented in this paper is represented by NEHL (PRL combined with TBLJP1). Valid priority rules such as PRO [30], PRD [34], and PRLJP1 [8] are all based on a metric to determine the initial job orders. These metrics can be regarded as the sum of two parts: the mean value, which describes the overall level of the data, and the standard deviation or skewness, which describes the degree of data chaos. According to these priority rules and their computational results, the higher the data mean and the more chaotic the data, the higher the job order should be. Based on this conclusion, kurtosis is introduced in this paper to further evaluate the degree of data chaos for the first time. Skewness and kurtosis are two indices used to describe the degree of data chaos in two directions. Our additional computational results show that a better solution quality can be obtained by multiplying the two together rather than adding them. Skewness and kurtosis are only significant for sufficiently large samples. For the PFSP with fewer machines, the statistical significance of skewness and kurtosis is not pronounced enough. Therefore, we introduce a coefficient represented by q before the product of skewness and kurtosis and set q = 1 when the number of machines is greater than 10; otherwise, we set q = 0.5. Then, all jobs can be ordered by the non-increasing sum of  
\[f=AVG_{i}+STD_{i}+q\times abs(SKE_{i})\times KUR_{i} \tag{1}\]  
where \(AVG_{i}\) is the average processing time of job \(i\), \(STD_{i}\) stands for the standard deviation of job \(i\), abs (\(SKE_{i}\)) represents the absolute value of the skewness of job \(i\), and \(KUR_{i}\) stands  
Figure 4: The idea of creating a sequence of k partial jobs. for the kurtosis of job \(i\). The parameters are defined as follows, and the pseudo-code for the proposed algorithm is shown in Figure 5.  
\[AVG_{i}=\frac{1}{\text{m}}\sum_{k=1}^{m}t_{i,k} \tag{2}\]  
\[STD_{i}=\sqrt{\frac{1}{\text{m}-1}\sum_{k=1}^{m}(t_{i,k}-AVG_{i})^{2}} \tag{3}\]  
\[q=\left\{ \begin{array}{cc}0.5&\text{m}\leq 10\\ 1&\text{m}>10\end{array} \right.