## 4 Design of Experiment for Parameter Tuning  
In this section, we present a Design of Experiments (DOE) approach [41] for parameter settings of the VBIH algorithm. In order to carry out experiments, we generate random instances with the method proposed in [9]. In other words, random instances are generated for each combination of \(n\in\{100,\,200,\,300,\,400,\,500,\,600,\,700,\,800\}\) and \(m\in\{20,\,40,\,60\}\). Five instances are generated for each job and machine combination. Ultimately, we obtained \(1200\) instances in total. We consider three parameters in the DOE approach. These are maximum block size (\(bMax\)), temperature adjustment parameter (\(\tau P\)), and the decision of whether or not to implement the local search to the partial solution after removal of a block of jobs. We have taken the maximum block size with seven levels as \(bMax\in(2,3,4,5,6,7,8)\); the temperature adjustment parameter with ten levels as \(\tau P\in(0.1,\,0.2,\,0.3,\,0.4,\,0.5)\); and the decision on the local search to partial solutions as \(pL\in(1,2)\). Note that \(pL=1\) means that the local search is applied to partial solutions whereas \(pL=2\) does not apply the local search to partial solutions. In the design of VBIH algorithm, there are \(7\times 5\times 2=70\) algorithm configurations, i.e., treatments. The VBIH algorithm is coded in C++ programming language on Microsoft Visual Studio 2013, and a full factorial design of experiments is carried out for each algorithm on a Core i5, 3.40 GHz, 8 GB RAM computer. Each instance is run for 70 treatments with a maximum CPU time \(T_{max}=10\times n\times m\) milliseconds. Note that it took 18 days to run the full factorial design. We calculate the relative percent deviation (RPD) for each instance-treatment pair as follows:  
\[RPD=(\frac{\mathit{CMAX}_{i}-\mathit{CMAX}_{min}}{\mathit{CMAX}_{min}})\times 100 \tag{18}\]  
where \(\mathit{CMAX}_{i}\) is the makespan value generated by the VBIH algorithm in each treatment and \(\mathit{CMAX}_{min}\) is the minimum makespan value found amongst 70 treatments. For each job size-treatment pair, the average RPD value is calculated by taking the average of five instances in each job size. Then, the response variable (ARPD) of each treatment is obtained by averaging these RPD values of all job sizes. After determining the ARPD values for each treatment as mentioned above, the main effects plots of the parameters are analyzed and given in Figure 3.  
As it can be seen from Figure 3, the following parameters have better ARPD values than the others: \(bMax=2\), \(\tau P=0.5\), and \(pL=1\). Furthermore, in order to see whether or not there is an interaction effect between parameters, an ANOVA analysis is also given in Table 2.  
Table 2 indicates that \(bMax\), \(tP\), and \(pL\) were statistically significant since higher magnitude of \(F\) values and \(p\)-values of parameter interaction effects are less than the significance level \(\alpha=0.05\). High magnitude of \(F\) value for \(pL\) also suggest that applying local search to partial solutions has a significant impact on the solution quality as mentioned in [5]. In terms of interaction effects, it can be observed that \(bMax\times tP\) interaction is not significant because the \(p\)-value is much higher than the significance level \(\alpha=0.05\). However, \(bMax\times pL\) and \(tP\times pL\) interactions were significant since their p values are less than the significance level \(\alpha=0.05\). The interaction effects plot for \(bMax\times pL\) is given in Figure 4.