Some metaheuristic algorithms for solving multiple cross-functional team selection problems

MOP-TS model

The following are several variables used for the model:

• K is the number of candidates.

• G denotes the number of groups.

• S represents the number of skills in the skillset.

• Z_g stands for the number of members in group gth ∀g =1 , …, G.

• The decision variables $X\in R^{K\times G}=x_{k,g}|x_{k,g}\in \left(0,1\right);k=1\dots K,g=1\dots G$ where:

• $x_{k,g}=\left\{\begin{array}{c}1\text{if member}k\text{th}\text{selected to group}g\mathrm{th}\hfill \\ 0\mathrm{otherwise}\hfill \end{array}\right.$

• $N_{g,s}=\left\{\begin{array}{c}1\text{if skill}s\mathrm{th}\text{is required to group}g\mathrm{th}\hfill \\ 0\mathrm{otherwise}\hfill \end{array}\right.$ ∀s =1 , …, S; g =1 , …, G.

• V_k,s is the rating score for skill sth of the candidate ith ∀s =1 , …, S; k =1 , …, K.

• L_g,s is the minimum required score for skill sth for group gth.

The objective functions were defined as follows:

• To select candidates who were fluent in the required skills by group: ${\displaystyle \max \left(f_g^{\mathrm{deep}}\left(X\right)=\sum _{i=1}^K\left(\sum _{s=1}^S{x}_{i,g}\mathrm{\ast }{N}_{g,s}\mathrm{\ast }{V}_{i,s}\right)\right)\forall g=1,\dots ,G.}$

• To select candidates who knew many of the required skills by group: ${\displaystyle \max \left(f_g^{\mathrm{wide}}\left(X\right)=\sum _{i=1}^K\sum _{s=1}^{S}min\left(1,x_{i,g}\mathrm{\ast }{V}_{i,s}\mathrm{\ast }{N}_{g,s}\right)\right)\forall g=1,\dots ,G.}$

Subject to:

• No candidate could join more than one group: ${\displaystyle \sum _{g=1}^G{x}_{i,g}\le 1\forall i=1,\dots ,K.\left(C1\right)}$

• No group was over team size: ${\displaystyle \sum _{i=1}^K{x}_{i,g}=Z_g\forall g=1,\dots ,G.\left(C2\right)}$

• Selected groups must respect the minimum required skills: ${\displaystyle \sum _{k=1}^K{x}_{k,g}\mathrm{\ast }{V}_{k,s}\ge L_{g,s}\forall g=1,\dots ,G;s=1,..,S.\left(C3\right)}$

Compromise programming to MOP-TS

The idea of compromise programming (CP) (Ringuest, 1992) is based on not utilizing any preference information or depending on assumptions about the relevance of objectives. This approach does not strive to discover numerous Pareto solutions. Instead, the distance between a reference point and the feasible objective region is reduced to identify a single optimal solution. Many studies have used CP for their MOPs, such as examination timetabling (Tung Ngo et al., 2021), teacher timetabling (Ngo et al., 2021b), enrollment timetabling (Ngo et al., 2021a), learning path recommendation (Son et al., 2021a), and task assignment (Son et al., 2021b). Several methods are used to select the preferred point  (Ngo et al., 2021c) or normalize the distance function (Ngo et al., 2022). In this situation, we introduced the compromised objective function as: $\min {\left({\sum }_{i=1}^L{w}_i|\frac{F_i-z_i^{\mathrm{\ast }}}{z_i^{\mathrm{worst}}-z_i^{\mathrm{\ast }}}{|}^2\right)}^{1/2}$ where L = 2∗G and $z_i^{\mathrm{worst}}={max}_{x\in X}{f}_i\left(x\right)$.

The z^∗ and z^worst were pre-calculated as:

• $z^{\mathrm{worst}}=\left\{d^{\min }|c^{\min }\right\}$

• $z^{\mathrm{\ast }}=\left\{d^{\max }|c^{\max }\right\}$

• F = f^deep|f^wide.

Where:

• $d^{\min }=\left\{d_g^{\min }|d_g^{\min }={\sum }_{s=1}^S\left({\sum }_{i=1}^{Z_g}{N}_{g,s}\mathrm{\ast }{P}_{s,i}\right),g=1,\dots ,G\right\}$

• $c^{\min }=\left\{c_g^{\min }|c_g^{\min }={\sum }_{s=1}^{S}min\left(1,P_{s,1}\mathrm{\ast }{N}_{g,s}\right),g=1,\dots ,G\right\}$

• $d^{\max }=\left\{d_g^{\max }|d_g^{\max }={\sum }_{s=1}^S\left({\sum }_{i=1}^{Z_g}{N}_{g,s}\mathrm{\ast }{R}_{s,i}\right),g=1,\dots ,G\right\}$

• $c^{\max }=\left\{c_g^{\max }|c_g^{\max }={\sum }_{s=1}^S{\sum }_{i=1}^{Z_g}min\left(1,R_{s,1}\mathrm{\ast }{N}_{g,s}\right),g=1,\dots ,G\right\}$ and

• R_s is sorted vector of $\left(V_{1,s},V_{2,s}\dots ,V_{K,s}\right)$ by descending order.

• P_s  is sorted vector of $\left(V_{1,s},V_{2,s}\dots ,V_{K,s}\right)$ by ascending order.

Genetic algorithm

The GA scheme is illustrated in Fig. 2, where the chromosome is represented like the decision variables X, but here we used list p as List < List < Integer > > of G items. Each p_g stores the indexes of selected candidates (Fig. 3). This mechanism allows for reducing the size of the original X.

The steps of the algorithm (Fig. 2) are described as:

• Generation of the initial population: We generated population P as the set of π indiviuals.

• Computation of the fitness function: For each p in P as: ${\displaystyle p\text{.fitness}=\left\{\begin{array}{c}\left(\sum _{i=1}^{2\mathrm{\ast }G}{w}_i{\left|\frac{F_i-z_i^{\mathrm{\ast }}}{z_i^{\mathrm{worst}}-z_i^{\mathrm{\ast }}}\right|}^2\right)\mathrm{if}\mathrm{val}\left(p\right)=0.\hfill \\ 1\mathrm{otherwise}\hfill \end{array}\right.}$ Where val(p) returns 1 if the solution p violates any constraints of $\left(C1\right),\left(C2\right),$ or (C3). The function returns 0 if p is valid. The violated individuals are punished by removing their effects. These solutions are potentially replaced by better quality solutions in the following genetic operations.

• Selection: We chose the selection rate of ϕ to keep the elite individuals in the next generation.

• Crossover: μ denotes the crossover rate where the new individual was constructed as follows: Figure 4 illustrates an example of the three steps of the crossover phase.

  1. Randomly select two individuals as parents denoted by p₁, p₂

  2. Gather selected members from p₁, p₂ to list pool.

  3. Randomly select items from pool to construct the new solution for the next generation.

• Mutation: We used the mutation rate β to select individuals from the population. These individuals’ genes were modified by randomly selected candidates from the list of K candidates.

• We defined a condition that if after α generations, the system could not find any better solutions, the algorithm would stop. Otherwise, continue to the selection phase.

ACO for multiple team selection

To design the scheme of the TS-ACO, we applied a similar data structure to represent the artificial ant (solution) as GA. Figure 5 shows the flow chart of the proposed ACO algorithm.

The details of each step of the ACO is described as follows:

• Initialization of the cost matrix: We created matrix C ∈ ℝ^K×G where C_k,g represented the cost if candidate kth was chosen for group gth as following: ${\displaystyle C_{k,g}=\sum _{z=1}^2{w}_{g\mathrm{\ast }z}{\left|\frac{M_{z\left(k,g\right)}-z_g^{\mathrm{\ast }}}{z_g^{\mathrm{worst}}-z_g^{\mathrm{\ast }}}\right|}^2.}$

Where: ${\displaystyle M_1\left(k,g\right)=\sum _{s=1}^S\left(N_{g,s}\mathrm{\ast }{V}_{k,s}\right)}$

and ${\displaystyle M_2\left(k,g\right)=\sum _1^S\left(min\left(1,N_{g,s}\mathrm{\ast }{V}_{k,s}\right)\right).}$

• Initialization of the pheromone matrix: We created matrix T ∈ ℝ^K×G where $T_{k,g}=\frac{1}{C_{k,g}}$.

• Generation of ant colony: We created new population P as the set of π artificial ants.

• Update of trail and computation of fitness: the population P was constructed as: $r_z=\mathrm{rand}\left(v_g^{\prime }\right),\forall z=1\dots Z_g$ are selected as candidates for group gth, ∀g = 1…G of the ant p, ∀p ∈ P, where:

  • rand represents the cumulative distribution function.

  • v_g denotes vector column gth of matrix T.

  • $v_g^{\prime }=\left\{\frac{v_{g,1}}{{\sum }_{i=1}^K{v}_{g,k}},\frac{v_{g,2}}{{\sum }_{i=1}^K{v}_{g,k}},\dots ,\frac{v_{g,K}}{{\sum }_{i=1}^K{v}_{g,k}}\right\}.$

• Update of pheromone matrix: We selected the top Φ of the best solution from P to list E, and then the T was updated as following:

• We also defined a condition that if, after α generations, the best fitness value of the population did not change, the algorithm stopped. Otherwise, we returned to the step of generating a new ant colony.

Computational complexity

Both GA and ACO operations are stochastic. It is difficult to determine the exact complexity of these algorithms, although we can calculate the cost of each iteration (Table 3). However, it is difficult to predict the number of iterations the algorithm will need to perform in advance when GA is O(K∗H∗S) and ACO is O(K∗H∗S∗G). However, the number of iterations is not predictable. Thus, if we assign N_GA and N_ACO as the number of iterations for each algorithm, we can say that the complexity of GA and ACO are O(N_GA∗K∗H∗S) and O(N_ACO∗K∗H∗G∗S), respectively. We predicted that the computation time of GA would be better than that of ACO. We used execution time (CPU time) to illustrate this in our experiments.

Experimental design

To evaluate the performance of the proposed algorithms on a real-world dataset, we used the dataset of 500 programming contestants  (Ngo et al., 2020). Contestants tackled programming exercises and were each tagged into different types of exercises across 37 categories. Figure 6 shows the statistical numbers corresponding to the available skills in the dataset.

We displayed the dataset with a two-dimensional space using t-SNE to find a similar probability distribution across the contestants in low-dimensional space (Fig. 7). Each point represents a contestant in the original space. The locations of the points indicate that the class difference between the contestants was very distinct. This allocation affects the search results because some members often appear in most of the solutions because they dominated in several skills. A dataset of the normal distribution can have many Pareto solutions.

The quality of the metaheuristic solution depends on factors such as the distribution of the data parameters. Since the algorithms are designed based on stochastic operations, it is not easy to evaluate the algorithmic complexity. Besides using the benchmarking dataset, we also generated 50 random datasets based on different distributions to compare the solution quality based on the statistical method. The 50 datasets included 300 candidates, and 30 skills were randomized based on one of six random distributions (Hypergeometric, Poisson, Exponential, Gamma, Student, and Binomial) using the Python SciPy library. We conducted the experiments on configured computers (Table 4).

Metaheuristic algorithms operate according to user customization through parameters. They significantly affect the performance of the algorithms. For example, one can bulk order search agents to increase the likelihood of finding a better-quality solution. However, this increases the execution time. We calibrated the parameter values (Table 5) used in this experiment by re-executing the algorithm several times.

Results

The results are described in three subsections. The first section compares the designed algorithm with the designed GA by Ngo et al. (2020) for a single team selection problem. The second subsection shows the results of multiple team selection. To assess the performance of the compared algorithms, we considered the aspects of the quality of the optimal solutions, processing time, and how they dealt with different decision scenarios. The final subsection compares metaheuristics and the exact algorithm on randomly generated datasets.

(1) Single Team Selection

In Ngo et al. (2020), the author compared GA algorithms (called GA-1) and DCA, and their GA-1 showed superior results when selecting three members from the tested dataset. To evaluate the quality of the proposed algorithms, we compared our designs to GA-1 when solving the single team selection problem using the same objective function. The results in Table 6 show that CPLEX could not find a solution from 500 candidates on the tested machine due to an out-of-memory error. Even when the number of candidates was decreased to 300, the objective value obtained by CPLEX was insufficient. GA-1 was designed to solve the single team selection problem. Meanwhile, the proposed GA and ACO aim to find multiple-team selections, but they can find a quality solution as GA-1 for single team selection.

(2) Multiple Team Selection

As mentioned in Section 2.2, the original objective functions were transformed into a distance function from the actual solution point to the ideal point. This requires the use of z^∗ and z^worst. The values of their elements are shown in Table 7. These values were readily calculated based on the candidate achievement data. We pre-calculated them for different scales of top K candidates from the tested dataset. The pair values $i\in \left\{1,2\right\},\left\{3,4\right\}$ and $\left\{5,6\right\}$ of $z_i^{\ast }$ and $z_i^{\mathrm{worst}}$ represent the best and worst scores, respectively, of the groups 1st, 2nd, and 3rd.

In this experiment, we used the query shown in Fig. 8. There were three selected groups, along with the required skills required for each group. We used the indexes of skills in arrays instead of listing their names. This query was used in the next section. The heat map illustrates the minimum required scores that needed to be archived by the selected team’s skills. No required skills are displayed in white.

Figure 9 represents the results of 15 executions of each algorithm’s GA, ACO, old GA, and GA using the previous search operations (Ngo et al., 2020) to eliminate individuals that violated the constraint. GA consistently showed a better median and the best solution when compared to old GA and ACO. In terms of time execution, old GA and GA easily outperformed ACO, and, due to the mechanism of removing solutions that violated constraints, old GA ran slower than GA. This shows that this algorithm seems to reduce the diversity of the population, therefore reducing the chances and speed of acquiring good solutions. CPLEX’s fitness value and timely execution results were overwhelmed by the three algorithms above. Therefore, we concluded that our proposed scheme can adapt better than the design of Ngo et al. (2020) for multiple team selection situations and our designed algorithms are superior to CPLEX based on these statistical results.

Figure 10 illustrates the performances of the tested algorithms across different system sizes: 50, 100, 300, and 500 candidates. The ACO was configured to achieve similar search results as GA. Calculations and updates on cost and pheromone matrices are increasingly expensive, especially as the size of the search space expands. This setting causes the processing time of ACO to be a few times higher than that of GA, and this cost is increasing, corresponding to the candidate’s size. The search operations of ACO require that every ant chooses its path by calculating the probability of each member to be selected for the groups in each iteration. Meanwhile, GA does not have to scan the whole candidate for mutation or crossover. This mechanism makes ACO take longer to finish an iteration, but it allows for more search capabilities over a larger space, provided the computational resources are expanded. CPLEX is not capable of handling 500 candidates. Both CPLEX’s solution quality and computation time show that it is comprehensively inferior to the proposed algorithms.

As described in Section 3A, search agents use fitness values to assess the quality of solutions. The use of stochastic operations allows agents to explore the search space. The more population diversity is ensured, the better the population’s ability to discover. Unlike Ngo et al. (2020), we only punish the violated resolutions to their fitness values instead of removing them or repairing them. This aims to maintain the diversity of the population. Figure 11 shows the number of violated solutions over generations for the best fitness values of the proposed algorithm with 500 candidates. This shows that the number of these solutions decreases when new generations are generated. However, the mutation/random selection process still produces a certain number of these solutions.

Fig. 12 shows the value of the fitness functions and the objective functions returned through iterations of the search processes for the best-obtained fitness values to the proposed algorithms. By observing the shape of the graphs of fitness functions produced by GA and ACO, we can see that the convergence of GA is slightly better than that of ACO. ACO’s search agents need more time to complete an iteration, which leads to a costlier total processing time. GA achieved the optimal solution at the 55th iteration. Meanwhile, ACO took 101 iterations to achieve the same result. The data distribution affects the reduction of values in the distance-based fitness function. If the standard deviation is significant, the objective value has a higher impact on the fitness value calculation when solutions are projected to a specific objective function in the objective space.

(3) Different decision scenarios

Approaches to the MOP problem based on the decomposition of multi-objective functions to a ranking function have many advantages. CP is more suitable for the decision-maker who cannot indicate the preferences to trade-off the specific goals. The combination of CP and EAs allows for a seamless transition between model and algorithm design by using compromised-objective and fitness functions that are both distance-based. In contrast, it is challenging to find optimal solutions on the Pareto frontier using MOEA (Emmerich & Deutz, 2018). However, decision-makers use the weighted parameters to manipulate the optimizer for different decision scenarios. We executed three algorithms several times using 10 different sets of weight parameter values on the dataset of 300 contestants. CPLEX cannot work on the scale of 500 candidates; therefore, we only tested the proposed GA and ACO at this scale. The obtained solutions are displayed in Fig. 12 with corresponding values of the weight parameters.

The obtained solutions shown in Fig. 13 may not totally dominate each other. Therefore, it is not easy to indicate how each algorithm deals with different decision-making scenarios. To evaluate this, we measured the hypervolume HVC (Ishibuchi et al., 2019) covered by the obtained solutions. A greater HVC indicates that the algorithm can produce a better Pareto frontier. The hypervolume was computed as: ${\displaystyle HVC=\frac{\mathrm{volume}\left(\bigcup _{s\in S}\left(s,z^{\mathrm{worst}}\right)\right)}{\mathrm{volume}\left(\mathrm{cube}\left(z^{\mathrm{\ast }},z^{\mathrm{worst}}\right)\right)}}$

where:

• s denotes the solution in the Pareto solution set S that is generated by the algorithm.

• cube(a, b) denotes the oriented axis hypercube that is formulated by points a and b in the objective space.

• volume(c) denotes the volume of the hypercube c in the objective space.

To evaluate the capabilities of the proposed CP-based method, we used genetic operations designed to implement a version of the NSGA-2 algorithm (Emmerich & Deutz, 2018). The parameters used to execute the algorithm are shown in Table 8. This method has the capability to search for a Pareto front with more than 500 solutions at both scales of 300 and 500 candidates after 2,123.6 and 3,021.3 s. We used a similar cost to execute evaluated CP-based algorithms several times using different parameters. The hypervolume obtained from these solutions is displayed in Table 9. The results show that the proposed evolutionary algorithms provided superior solutions to CPLEX across different decision-making scenarios. Although the solutions generated by CP-based algorithms produced better hypervolume values than NGSA-2, it is hard to conclude that the CP-based method is better than MOEA. However, one can say that the Pareto frontier obtained by NGSA-2 may contain lower quality solutions than the proposed approach. Therefore, this approach requires significant computational overhead and is challenging to adapt to a real-world environment when the user may only need a single final solution. The user experience factor may not contribute to the effort to improve search speed.

(4) Metaheuristics vs. EA

CPLEX relaxes the original problem within bounds and solves the relaxation for mixed-integer linear programming (MILP) and the mixed-integer quadratic programming (MIQP). When CPLEX tries to solve a nonconvex MIQP to global optimality, it cannot provide the relaxation of the original problem that is not bounded. The proposed optimization problem is NP-Hard. The “branch and bound” may require exploring all possible permutations in the worst case (Chopard & Tomassini, 2018). Meanwhile, EAs (metaheuristics and general) use stochastic operations to determine possible solutions. Fortunately, the best solution can be found very early (even in the first iteration), but in the worst case, the solutions found in the first and last iteration are of the same quality. Therefore, determining the computational complexity of the two approaches from a theoretical perspective is a challenging problem and beyond the scope of this study (Darwish, 2018). Instead of giving theoretical calculations, we evaluated the performance of these algorithms using statistics.

Table 10 shows the results of three algorithms (GA, ACO, and CPLEX) when solving the problem using 50 generated datasets. In the 50 data sets, there were only four data sets that CPLEX found quality solutions as ACO or GA. The rest gave worse results. The execution time of CPLEX was also many times longer. While both GA and ACO solved the problem in less than 1 min, CPLEX took 3–5 min to execute. The above results show that the proposed metaheuristic algorithms are more efficient than the tested-exact method in terms of solution quality and processing time.