Hybrid Particle Swarm and Whale Optimization Algorithm for Multi-Visit and Multi-Period Dynamic Workforce Scheduling and Routing Problems

Abstract

This paper focuses on the dynamic workforce scheduling and routing problem for the maintenance work of harvesters in a sugarcane harvesting operation. Technician teams categorized as mechanical, hydraulic, and electrical teams are assumed to have different skills at different levels to perform services. The jobs are skill-constrained and have time windows. During a working day, a repair request from a sugarcane harvester may arrive, and as time passes, the harvester’s position may shift to other sugarcane fields. We formulated this problem as a multi-visit and multi-period dynamic workforce scheduling and routing problem (MMDWSRP) and our study is the first to address the workforce scheduling and routing problem (WSRP). A mixed-integer programming formulation and a hybrid particle swarm and whale optimization algorithm (HPSWOA) were firstly developed to solve the problem, with the objective of minimizing the total cost, including technician labor cost, penalty for late service, overtime, travel, and subcontracting costs. The HPSWOA was developed for route planning and maintenance work for each mechanical harvester to be provided by technician teams. The proposed algorithm (HPSWOA) was validated against Lingo computational software using numerical experiments in respect of static problems. It was also tested against the current practice, the traditional whale optimization algorithm (WOA), and traditional particle swarm optimization (PSO) in respect of dynamic problems. The computational results show that the HPSWOA yielded a solution with significantly better quality. The HPSWO was also tested against the traditional genetic algorithm (GA), bat algorithm (BA), WOA, and PSO to solve the well-known CEC 2017 benchmark functions. The computational results show that the HPSWOA achieved more superior performance in most cases compared to the GA, BA, WOA, and PSO algorithms.

1. Introduction

Sugarcane is one of the major economic agricultural crops of Thailand. In 2021, Thailand was ranked the world’s second largest sugar exporter, and the fourth largest sugar producer [1]. However, maintaining the maximum milling capacity is still a challenge in Thai sugar mills. There is a constant need to improve inbound logistics activities, especially the harvesting process, to obtain a steady supply of sugarcane for the sugar mills.

In sugarcane inbound logistics, the sugarcane harvest is an important process since the supply of sugarcane must match the capacity of the mill [2]. In the past, manual harvesting was intensively used but it is not presently practical due to the challenging nature of the work, leading to labor shortages and high labor costs. Thus, mechanical harvesters have been introduced and widely used to overcome the labor scarcity, to control the high labor wages, and importantly, to obtain a daily steady sugarcane supply that matches the mill’s maximum capacity. Improvement of sugarcane delivery from fields to the sugar mill at the right time is thus crucial and is required during the crushing season [3].

Often, the challenge of using mechanical harvesters is their low utilization. In actual practice, operation of a mechanical harvester is a complicated activity since there are various teams used to operate the sugarcane harvester. The teams involved include (1) transporting trailers or trucks (so-called forwarders) moving by the side of the mechanical harvester to transport the chopped sugarcane to the mill, (2) teams providing maintenance and fuel supply services, and (3) teams supporting emergency services (e.g., fire trucks). These three components are usually referred to as harvest fronts [4]. It is apparent that harvesters are the most crucial machines in the mechanical harvesting system, since sugar production relies on machine performance, which is directly related to field efficiency. Therefore, time loss, especially due to delays in machine repair, results in a significant reduction in sugarcane harvest capacity and the mill’s fluctuating production capacity, as well as increased production costs.

In the current practice in Thailand, time loss due to delayed machine maintenance is high due to inappropriate field service scheduling, which is complicated due to the increasing number of harvesters, number and skills of technicians, and types of maintenance services for the harvesters (preventive maintenance (PM), break-down maintenance (BM), and corrective maintenance (CM)), each of which includes three sub-systems comprising mechanical, hydraulic, and electrical systems [5]. Normally, each sub-system requires specific types of technicians to perform maintenance tasks. The scheduling is complicated when a harvester requires more than one service operation (PM, CM, or BM) simultaneously. In this case, if any service operation is not completed, the harvester is still not able to work. More importantly, in order to avoid losses in harvester productivity, a subcontractor is required when harvesters require servicing but a technician team is not available. Time loss also occurs because mechanical harvesting is a dynamic process. Harvesters and technicians usually move and operate in different locations during each time period. More importantly, uncertainties related to harvesters usually occur in the planning horizon (e.g., harvester breakdowns). These aforementioned factors cause delays in servicing, resulting in reduced utilization of the harvesters. Therefore, the potential approach to maintaining the harvesting capacity is to keep the harvesters in operation for the longest period by minimizing the time loss due to these issues.

Given the mechanical harvest plan, in this paper, the allocation and scheduling of multi-heterogeneous technician teams, with consideration of transportation for the teams and the operations of sugarcane harvesters, were addressed to reduce time loss. This is because, in actual practice, harvesters and technicians usually move and operate in different locations during each time period, and uncertainties related to harvesters usually occur in the planning horizon (e.g., harvester breakdowns). All these factors, involving limited resources but more dynamic operations, were taken into consideration in this paper. In other words, a more realistic and challenging dynamic workforce scheduling and routing problem (DWSRP) was addressed. A decision in terms of the allocation and scheduling of different types of service staff to perform operations at several locations was taken such that the delay in harvester servicing and travel distance of technicians were minimized. In this problem, the objective was to minimize the total cost, including technician labor cost, penalty for late service, overtime, travel, and subcontracting costs.

Since there are several types of teams, each of which includes many teams, and some harvesters might require more than one type of service during each time period, this problem was formulated as a multi-visit and multi-period dynamic workforce scheduling and routing problem (MMDWSRP) (see the systematic diagram shown in Figure 1). The problem consisted of two sub-problems: (1) technician team allocation and scheduling and (2) routing of the technicians. In the first sub-problem, workforce allocation of technician teams was performed through consideration of clusterization in certain geographical areas, eligibility, and skills. Then, their sequences and routings were determined in the second sub-problem, with consideration of time windows, maximum duration of work, and dynamic requests of harvesters during a day. In general, only a few studies in the literature have formulated the workforce scheduling and routing problem (WSRP) as a dynamic routing problem. To the best of our knowledge, the DWSRP has been implemented in a home health-care routing problem (DHHC) [6,7], a telecommunication problem [8,9], and a technician routing and scheduling problem (DTRSP) [10].

The main contribution of this paper is the novelty in considering dynamics to be part of the problem. In other words, we believe our work is the first of its kind in addressing the WSRP problem, which combines a dynamic problem and multiple customer visits by technician teams. In such a problem, the service skills of staff (mechanical, hydraulic, and electrical), types of maintenance (BM, CM, and PM), time windows, eligibility, maximum duration of work, dynamic requests of harvesters, and the shifts in the harvester’s position during a day may be considered simultaneously. To solve a realistic-sized dynamic problem, a novel swarm-based meta-heuristic optimization algorithm, the whale optimization algorithm (WOA), was proposed. The WOA has been used intensively to solve complex optimization problems [11] such as allocation [12], scheduling [13,14,15], routing [16,17], and manufacturing [18,19]. Even though the WOA has been used extensively in various problems, it has never been used to solve WSRPs. This study is believed to be the first to employ the WOA and the hybrid particle swarm and whale optimization algorithm (HPSWOA) to solve the MMDWRSP. The proposed algorithm was validated against Lingo computational software using numerical experiments in static problems. It was also tested against the current practice, the traditional whale optimization algorithm (WOA), and traditional particle swarm optimization (PSO) in respect of dynamic problems. The computational results demonstrate that the HPSWOA yielded a solution with significantly better quality. The HPSWO was also tested against the traditional genetic algorithm (GA), bat algorithm (BA), WOA, and PSO to solve the well-known CEC 2017 benchmark functions. The computational results show that the HPSWOA achieved a more superior performance in most cases compared to the GA, BA, WOA, and PSO algorithms.

The remainder of this paper is organized as follows. In Section 2, a brief literature review of the WSRP is presented. Section 3 presents the mathematical model. In Section 4, the HPSWOA algorithm for MMDWSRP problems is developed. Computational results are discussed in Section 5. The last section delineates conclusions and future research directions.

2. Literature Review

Presently, a significant amount of research has been undertaken in respect of WSRP, with consideration of several factors including staff skills, planning horizon, time windows, transportation modes, teaming, and the dynamics of the problem. These factors may result in different WSRP characteristics. The most studied and uncomplicated WSRP involved a single-period planning horizon without team building. Studies of this problem type included those of Xu and Chiu [20], Dohn et al. [21], Pillac et al. [22], Pinheiro et al. [23], and Mathlouthi et al. [24], while the studies of Cordeau et al. [25], Kovacs et al. [26], Anoshkina et al. [27], and Çakırgil et al. [28] focused on a WSRP with a single period including team building. To solve a WSRP with a single period, both mathematical and heuristic approaches (e.g., variable neighborhood search, adaptive large neighborhood search, and decomposition solution methods) were used to solve the problem.

During the past decade, a multiple-period WSRP with team building has been investigated by several researchers. For example, in the study of Blakeley et al. [29], Zamorano and Stolletz [30], Pekel [31], and Punyakum et al. [32], technicians with different skills were considered. A tabu search approach was applied to solve the problem. Later, Zamorano and Stolletz [30] studied a technician routing and scheduling problem (TRSP) with consideration of the skills of technicians and time windows in a multiple-period horizon in order to decide the daily team formation, and daily team assignments to tasks and to routes. A branch-and-price approach was applied to obtain the operating costs. Recently, Pekel [31] focused on a multi-period TRSP with the objective of minimizing operating costs. An improved particle swarm optimization (IPSO) approach, which is a hybrid of particle swarm optimization (PSO) and the neighborhood operator, was used to solve the problem. Later, Punyakum et al. [32] concentrated on a multi-visit and multi-period WSRP (MMWSRP) in the field service operation of a sugarcane mill company. A hybrid differential evolution and particle swarm optimization (HDEPSO) algorithm was used to obtain the solution, with the objective of minimizing operating costs. The HDEPSO was tested against the PSO and differential evolution. It was found that the HDPSO resulted in a solution of much higher quality.

For a dynamic problem, prior to planning, information about customers may be unknown or only partially known. For example, in the study of Lesaint et al. [8], a workforce scheduling problem was considered for British Telecom, in which technicians were assigned to tasks that arrived dynamically during the day. The objective of the problem was minimizing cost. The simulated annealing approach was employed to solve the problem. Borenstein et al. [9] modified Lesaint et al.’s problem [8] with the objective of maximizing the number of completed tasks. The decomposition method was applied to solve the problem. In 2018, Pillac et al. [10] studied a dynamic technician routing and scheduling problem (DTRSP) in which new requests appeared over time. Time windows, tools, spare parts, and skills were taken into account in designing staff routing to respond to requests for service, and an adaptive large neighborhood search approach was applied to solve the problem, with the objective of minimizing the total working time and total distance. Later, Ouertani et al. [6] studied a dynamic home health-care routing problem (DHHC) in which new calls from patients could arrive after the departure of the caregiver, and the caregiver’s route had to be replanned to include these new patients. The objective of DHHC was to minimize the travelling cost. The hypermutation genetic algorithm method was used to obtain the solution. Later, Demirbilek et al. [33] studied the home healthcare routing and scheduling of multiple nurses in a dynamic environment. A scenario-based method for several nurses was used to determine the maximum number of patient visits for a group of nurses during the planning horizon.

In our study, a more realistic dynamic workforce scheduling and routing problem (DWSRP) was investigated, which could be formulated as a multi-visit and multi-period dynamic workforce scheduling and routing problem (MMDWSRP) with the objective of minimizing the total cost, consisting of labor cost, traveling cost, penalty for late service, overtime, and subcontracting cost. To the best of our knowledge, the MMDWSRP has not been addressed in the literature. Among the studies reported in the literature, the only one that closely matches this study is the MMWSRP (Punyakum et al. [32]). That study does not involve the dynamic WSRP, employs a different methodology, and does not consider technician labor cost in the objective function.

Recently, the whale optimization algorithm (WOA) has been efficiently applied to solve complicated optimization problems [11]. The WOA was developed by Mirjalili and Lewis [34]. It imitates the hunting manipulation of humpback whales. The WOA with a stochastic nature is simple, flexible, and fast in convergence speed; it has gained a great deal of attention among researchers in various disciplines [35]. The WOA is widely used in several fields such as manufacturing, allocation, scheduling, and routing problems. For example, in 2019, Saleh et al. [12] determined the optimal allocation of distributed generation and capacitors in radial distribution systems while reducing single- and multi-objective functions (network power losses, voltage deviation, and total operating cost). The weighted sum method was used to create the multi-objective function. The WOA was used to obtain the solution. Later, Liu et al. [15] used a hybrid WOA improved with Lévy flight and differential evolution to solve job shop scheduling problems. The proposed algorithm achieved superior performance when solving small and large benchmark instances compared to nine other algorithms and was not inferior when solving hard instances. In the same year, Tanvir et al. [19] studied stainless steel 304 in respect of turning operations. The effects of machining parameters such as the feed rate and depth of cut on surface roughness, cutting forces, cutting speed, peak tool temperature, power, material removal rate, and heat rate were studied. A hybrid WOA and Grey relational analysis was used to decide the machining conditions for optimizing the production unit cost and production quality. Later, Dewi and Utama [16] considered a green vehicle routing problem (GVRP) with the objective of minimizing the distribution cost, considering the vehicle usage cost, fuel consumption, and carbon emissions. A hybrid whale optimization algorithm developed by combining the WOA algorithm, TS algorithm, and local searching was used to obtain the minimum distribution cost.

Hybrid algorithms are more frequently employed. A hybrid algorithm combines two or more methods and is intended to perform better than each algorithm individually [36]. Wang and Shoup [37] investigated a hybrid optimization technique that combines the best elements of four conventional techniques for optimization, including the Hooke and Jeeves pattern search approach, the simplex method, the pattern search method, and the random search method. Their method was tested against 11 conventional test functions, and the results show that it performed better than eight of the selected search methodologies. Khajehzadeh et al. [38] combined the adaptive gravitational search algorithm and pattern search method to solve the multi-objective optimization of reinforced concrete retaining walls. The embedded CO₂ emissions and total cost were the objective functions. Later, Koessler and Almomani [39] studied a hybrid particle swarm optimization and pattern search algorithm (PSOPS) to solve benchmark functions and the optimization of basin networks. When compared against other hybrid techniques for a basin network optimization problem, the PSOPS algorithm outperformed them. Recently, Khajehzadeh et al. [40] proposed the hybrid tunicate swarm algorithm and pattern search (HTSAPS). In benchmark test functions and model application, the HTSAPS was also compared to the tunicate swarm algorithm, sine-cosine algorithm, and grey wolf optimization method. The numerical studies demonstrated that the HTSAPS offered better optimal solutions and outperformed the other techniques. In 2004, Victoire et al. [41] concentrated on the economic dispatch with valve-point effect problem. A hybrid particle swarm optimization and sequential quadratic programming approach was developed to solve the problem with the objective of minimizing fuel costs. Later, Eslami et al. [42] considered damping controller design for power system oscillations. The damping ratio and damping factor were included in the objective functions. The hybrid genetic algorithm and sequential quadratic programming were used to solve the problem. Recently, Abdel-Mawgoud et al. [43] studied integrating a battery energy storage system into distribution networks. A hybrid arithmetic optimization algorithm and a sine cosine algorithm were developed to solve the problem with the objective of minimizing the whole active power loss of the network.

From the literature, it was found that the most commonly used solution methods for WSRPs are mathematical programming and metaheuristics. However, the WOA has not been used for solving WSRPs. In this study, the WOA, PSO, and the hybrid particle swarm and whale optimization algorithm were proposed to solve the MMDWRSP. The main contribution of this paper involves considering the WSRP by combining a dynamic problem and multiple customer visits by technician teams, which has not previously been proposed in the literature.

3. Mathematical Formulation

In this study, there were several types of teams and each of which included various teams. Some mechanical harvesters might require more than one type of service during each time period; this problem was therefore developed as a multi-visit and multi-period dynamic workforce scheduling and routing problem (MMDWSRP), where new maintenance service requests in respect of sugarcane mechanical harvesters appear while the technician teams are performing maintenance tasks on their routes. In this section, the mathematical model of the MMDWSRP is developed with static data. It is extended from the mathematical model of Punyakum et al. [32], which was developed as a multi-visit and multi-period workforce scheduling and routing problem (MMWSRP).

In the MMDWSRP, $J$ harvesting jobs require $S$ sub-systems, namely mechanical, hydraulic, and electrical, for maintenance. Job $m$ might require more than one sub-system’s services at once throughout each period $t$. For each sub-system, $s\in S$ represents a set of sub-systems (mechanical, hydraulic, and electrical), $R_s$ represents a set of teams in sub-system $s$, and $L_s$ represents a set of skill types in sub-system $s$ (1 = basic, 2 = medium, 3 = expert). In a day, each technician team may provide services in respect of more than one job (e.g., sugarcane mechanical harvester) dispersed in several sugarcane fields, and the jobs may require servicing by more than one technician team for different sub-systems. Since sugarcane fields are far apart, jobs may have to wait a long time for servicing by the technician teams. A technician team usually works in more than one sugar field in a day, and thus usually moves from one field to another. In order to prevent losses in harvester productivity when no technician team is available, a subcontractor is needed. The following are the assumptions for mathematical formulation of the MMDWSRP: (i) only one technician team in a sub-system is allowed to service a job; (ii) each technician team always has the required spare parts and equipment; (iii) before starting the servicing, technicians are allocated to a team; (iv) only one size of transport vehicle is considered; (v) no splitting of maintenance services is allowed; (vi) a technician team will go back to the company only when they finish the assigned servicing; (vii) degree of dynamism = 0. Therefore, modifications in the MMDWSRP were made from the MMWSRP model developed by Punyakum et al. [32] as follows: (1) the labor cost of technicians was added to the objective function; (2) since this problem is a dynamic WSRP, constraints 13 and 14 in this mathematical model were revised from constraints 16–19 [32] of the MMWSRP mathematical model to ensure that the number of technicians and the technician teams’ skill levels and types are sufficient for servicing; and (3) constraint 15 was added to the MMDWSRP mathematical model to determine the number of teams required to service the mechanical harvesters. The following is a list of the notation used in this paper.

Indices:

$m,n$	Jobs and depot indices ($m,n=1,2,\dots ,$ J)
$s$	Sub-systems; mechanical, hydraulic, and electrical sub-systems ($s=1,2,\dots ,S$)
$r$	Technician teams ($r=1,2,\dots ,R_s$)
$t$	Working periods ($t=1,2,\dots ,T$)
$l$	Technician team skills ($l=1,2,\dots ,L_s$)

Parameters:

$J$	Set of jobs and depot
$J\prime$	Set of jobs ($m,n=2,3,\dots ,J$)
$S$	Set of sub-systems
$R_s$	Set of technician teams in sub-system $s$
$T$	Set of working periods
$L_s$	Set of skill types in sub-system $s$
$LC{c}_{r,s}$	Labor cost of technician team $r$ in sub-system $s$ (units)
$TC{c}_{m,n}$	Technician transportation cost moving from jobs $m$ to $n$ (units)
$LT{c}_m$	Penalty cost from a job $m$’s that was finished late (unit per minutes)
$OV{c}_{r,s}$	Overtime cost of technician team $r$ giving service sub-system $s$ (unit per minutes)
$ST{c}_{m,s}$	Subcontracting cost of sub-system $s$ of jobs $m$ (unit)
$e_t$	Starting time of period $t$
$f_t$	Completion time of period $t$
$a_m$	Starting time of job $m$
$b_m$	Completion time of job $m$
$t{v}_{m,n}$	Travel time from job $m$ to job $n$ (minutes)
$p{t}_{m,s}$	Service time of job $m$ in sub-system $s$ (minutes)
$w_{r,s,l}$	Proficiency of technician team of skill $l$ in sub-system $s$
$u_{m,s,l}$	For job $m$, a technician team with the attribute $w_{r,s,l}$ is required
$o_s$	The quantity of technician teams in sub-system $s$
$q_{m,s}$	The quantity of technicians needed to complete job $m$ for sub-system $s$
$k_{r,s}$	The quantity of technician team $r$ in sub-system $s$
$d_{m,s}$	= 1 if job $m$ needs a team of sub-system $s$ = 0, otherwise
$H$	Large number
$LTm$	Maximum allowed delay time
$OVm$	Maximum allowed overtime

Decision variables:

$L{T}_m$	Amount of service delay time of job $m$ (minutes)
$O{V}_{s,r,t}$	Overtime of team $r$ in sub-system $s$ in period $t$ (minutes)
$c_{m,s,r,t}$	Starting time of team $r$ in sub-system $s$ for job $m$ in period $t$
$X_{m,n,s,r,t}$	= 1, if team $r$ of sub-system $s$ leaves job $m$ for job $n$ in period $t$; = 0, otherwise
$Y_{m,s,r,t}$	= 1, if team $r$ of sub-system $s$ services job $m$ in period $t$; = 0, otherwise
$Z_{s,r,t}$	= 1, if using technician team $r$ of sub-system $s$ in period $t$; = 0, otherwise
$S{T}_{m,s}$	= 1, if job $m$ requires a service from a subcontract for subsystem $s$; = 0, otherwise

Objective function:

$$\begin{gathered} Minimize:{\displaystyle \sum }_{s\in S}^{}{\displaystyle \sum }_{r\in R_s}^{}{\displaystyle \sum }_{t\in T}^{}LC{c}_{r,s}·Z_{s,r,t}+{\displaystyle \sum }_{m\in J}^{}{\displaystyle \sum }_{n\in J}^{}{\displaystyle \sum }_{s\in S}^{}{\displaystyle \sum }_{r\in R_s}^{}{\displaystyle \sum }_{t\in T}^{}{X}_{m,n,s,r,t}·TC{c}_{m,n}+{\displaystyle \sum }_{m\in J}^{}L{T}_m·LT{c}_m \\ +{\displaystyle \sum }_{s\in S}^{}{\displaystyle \sum }_{r\in R_s}^{}{\displaystyle \sum }_{t\in T}^{}O{V}_{s,r,t}·OV{c}_{r,s}+{\displaystyle \sum }_{m\in J}^{}{\displaystyle \sum }_{s\in S}^{}S{T}_{m,s}·ST{c}_{m,s} \end{gathered}$$

(1)

Subject to:

${\displaystyle {\displaystyle \sum }_{r\in R_s}^{}{\displaystyle \sum }_{t\in T}^{}}{Y}_{m,s,r,t}+S{T}_{m,s}=d_{m,s}$	$\forall m\in J\prime ,s\in S$	(2)
${\displaystyle {\displaystyle \sum }_{n\in J\prime }^{}}{X}_{m,n,s,r,t}=Y_{m,s,r,t}$	$$\begin{gathered} \forall m\in J\prime ,s\in S, \\ r\in R_s,t\in T;m\ne n \end{gathered}$$	(3)
${\displaystyle {\displaystyle \sum }_{n\in J\prime }^{}}(X_{m,n,s,r,t}-X_{n,m,s,r,t})=0$	$$\begin{gathered} \forall m\in J\prime ,s\in S, \\ r\in R_s,t\in T;m\ne n \end{gathered}$$	(4)
$(c_{m,s,r,t}+p{t}_{m,s}+t{v}_{m,n})·Y_{m,s,r,t}-H\left(1-X_{n,m,s,r,t}\right)\le s_{m,s,r,t}$	$$\begin{gathered} \forall m\in J, n\in J\prime , \\ s\in S,r\in R_s,t\in T; \\ m\ne n \end{gathered}$$	(5)
$c_{m,s,r,t}\ge a_m·Y_{m,s,r,t}$	$$\begin{gathered} \forall m\in J\prime ,s\in S, \\ r\in R_s,t\in T \end{gathered}$$	(6)
$c_{m,s,r,t}-\left(b_m-p{t}_{m,s}\right)·Y_{m,s,r,t}\le L{T}_m$	$$\begin{gathered} \forall m\in J\prime ,s\in S, \\ r\in R_s,t\in T \end{gathered}$$	(7)
$L{T}_m\le LTm$	$\forall m\in J\prime$	(8)
$c_{n,s,r,t}\ge (e_t+t{v}_{m,n})·X_{m,n,s,r,t}$	$$\begin{gathered} \forall n\in J\prime ,s\in S, \\ r\in R_s,t\in T,m=1 \end{gathered}$$	(9)
$c_{m,s,r,t}-O{V}_{s,r,t}+(t{v}_{m,n}+p{t}_{m,s})·X_{m,n,s,r,t}\le f_t$	$$\begin{gathered} \forall m\in J\prime , s\in S, \\ r\in R_s,t\in T,n=1 \end{gathered}$$	(10)
$O{V}_{s,r,t}\le OVm$	$s\in S,r\in R_s,t\in T$	(11)
${\displaystyle {\displaystyle \sum }_{n=J\prime }^{}{\displaystyle \sum }_{r\in R_s}^{}}{X}_{m,n,s,r,t}\le o_s$	$\forall s\in S,t\in T,m=1$	(12)
$q_{m,s}·Y_{m,s,r,t}\le k_{r,s}$	$$\begin{gathered} \forall m\in J\prime ,s\in S, \\ r\in R_s,t\in T \end{gathered}$$	(13)
$u_{m,s,l}·Y_{m,s,r,t}\le w_{r,s,l}$	$$\begin{gathered} \forall m\in J\prime ,s\in S, \\ r\in R_s,t\in T,l\in L_s \end{gathered}$$	(14)
$Z_{s,r,t}\ge Y_{m,s,r,t}$	$$\begin{gathered} \forall m\in J\prime ,s\in S, \\ r\in R_s,t\in T \end{gathered}$$	(15)
$X_{m,n,s,r,t}=\left\{0, 1\right\}$	$$\begin{gathered} \forall m,n\in J,s\in S, \\ r\in R_s,t\in T \end{gathered}$$	(16)
$Y_{m,s,r,t}=\left\{0, 1\right\}$	$$\begin{gathered} \forall m\in J,s\in S, \\ r\in R_s,t\in T \end{gathered}$$	(17)
$Z_{s,r,t}=\left\{0, 1\right\}$	$$\begin{gathered} \forall s\in S,r\in R_s, \\ t\in T \end{gathered}$$	(18)
$S{T}_{m,s}=\left\{0, 1\right\}$	$\forall m\in J,s\in S$	(19)

Objective function (1) is used to minimize the total cost, including labor cost, travel cost of technician teams, penalty for late services, overtime, and subcontracting cost. The mathematical formulation constraints are described in

Table 1. Mathematical formulation constraints.

Constraint	Description
(2)	Ensures that job $m$ in sub-system $s$ is serviced either by a team or a subcontractor.
(3)	Indicates that job $m$ is serviced by team $r$ in sub-system $s$.
(4)	Ensures that, after service completion, the team transport vehicles had to leave the job.
(5)	Indicates that the start time for servicing job $n$ is equal to the sum of the completion time of job $m$ and the travel time of team $r$.
(6)	Ensures that a team is able to start servicing only when the job is ready.
(7)	Ensures that if a team cannot finish the service in time, a delay occurs.
(8)	Ensures that the service delay must not exceed an allowed delay time.
(9)	Specifies the time when a team can start a service.
(10)	Indicates team overtime for late arrival at the depot.
(11)	Ensures that the amount of overtime used by each team in a given period must not exceed the allowed overtime.
(12)	Ensures that the number of teams servicing the jobs must not exceed the number of teams that are actually available.
(13)	Ensures that the number of technicians in team $r$ in sub-system $s$ is sufficient for servicing job $m$.
(14)	Ensures that the technician teams’ skill levels and types are adequate for servicing job $m$.
(15)	Determines technician team $r$ of the sub-system $s$ required during period $t$.
(16)–(19)	The binary variable constraints.

.

4. The Proposed Method

The hybrid particle swarm and whale optimization algorithm (HPSWOA) is described below.

4.1. Initial Solution

The proposed algorithm starts by generating an initial population size or the number of vectors of each sub-system (mechanical, hydraulic, and electrical). The vector consists of a job vector and a technician team vector. Firstly, each job vector will be randomly generated, equal to the number of jobs. Decryption arranges each vector’s rank order value in ascending order to obtain the sequence of jobs to be completed. An additional vector is the technician team vector, which will also be randomly generated and equal to the number of technician teams. Secondly, each technician team vector will be randomly generated, equal to the number of jobs. The technician team vectors are used to assign a technician team to a job. Figure 2 shows an illustrative example of the construction of a vector of a sub-system. There are five jobs and two technician teams. For example, the job vector is 0.51-0.83-0.64-0.11-0.25. The job ID is 1-2-3-4-5. After the job vector is sorted in ascending order, the vector obtained is 0.11-0.25-0.51-0.64-0.83. Therefore, the job ID is 4-5-1-3-2. The first technician team vector is 0.43–0.24. The technician team ID is 1–2. After the technician team vector is sorted in ascending order, the technician team vector obtained is: 0.24–0.43. Therefore, the technician team rank is 2–1. To establish the sequence for all technician teams, the procedure is as follows. The jobs are assigned by rank of the job ID to the technician team rank. If the first team’s skill and the number of team members do not match with the job, the job will be assigned to the next technician team, and so on until all jobs have been assigned to the teams. Jobs that cannot be assigned to the technician teams will be outsourced to subcontractors. The objective calculation as shown in Equation (1) is used to minimize the total cost.

4.2. Particle Swarm Optimization (PSO)

Eberhart and Kennedy [44] introduced PSO. The simulation of social and psychological expression in fish and birds served as the inspiration for this algorithm. In looking for food, a living swarm (such as a school of fish or a flock of birds) will scatter and then regroup in search of the best feeding area. The mathematical formulae that produce updates in respect of position and velocity throughout the course of one iteration are given by Equations (20) and (21):

$$X_{i,j}\left(t+1\right)=X_{i,j}\left(t\right)+V_{i,j}\left(t+1\right)$$

(20)

$$V_{i,j}\left(t+1\right)=w{V}_{i,j}\left(t\right)+C_1{r}_1\left(Pbes{t}_{i,j}\left(t\right)-X_{i,j}\left(t\right)\right)+C_2{r}_2\left(Gbes{t}_j\left(t\right)-X_{i,j}\left(t\right)\right)$$

(21)

where $X_{i,j}$ and $V_{i,j}$ are the position and velocity of the $j$th member of the $i$th vector (or particle), $t$ is the current iteration, $Gbes{t}_j$ is the global best position of all vectors, and $Pbes{t}_{i,j}$ is the personal best position of the $j$th member of the $i$th vector. $C_1$ and $C_2$ are personal and global best position weights, respectively; $r_1$ and $r_2$ are random numbers from 0 to 1; and $w$ is the inertial weight.

4.3. Whale Optimization Algorithm (WOA)

This algorithm is built on the copying of natural phenomena in a way that mimics animal behavior. The algorithm was developed by Mirjalili and Lewis [34] and inspired by the foraging behavior of humpback whales. Humpback whales have a unique prey trait when they know the location of their prey. Humpback whales create a spiral bubble around their prey that continuously encircles them, while swimming in a spiral and rising to the top to eat their prey.

4.3.1. Encircling Prey

When humpback whales know the location of their prey, they surround them. This determines that the best solution or the closest is the location of the target prey. The other solution then updates their position with reference to the prey’s position. This is similar to humpback whales swarming around their prey, which is represented by Equations (22) and (23).

$$\stackrel{\rightharpoonup }{X}\left(t+1\right)={\stackrel{\rightharpoonup }{X}}_{best}\left(t\right)-\stackrel{\rightharpoonup }{A}·\stackrel{\rightharpoonup }{D}$$

(22)

$$\stackrel{\rightharpoonup }{D}=\left|\stackrel{\rightharpoonup }{C}·{\stackrel{\rightharpoonup }{X}}_{best}\left(t\right)-\stackrel{\rightharpoonup }{X}\left(t\right)\right|$$

(23)

$$\stackrel{\rightharpoonup }{A}=2\stackrel{\rightharpoonup }{a}·\stackrel{\rightharpoonup }{r}-\stackrel{\rightharpoonup }{a}$$

(24)

$$\stackrel{\rightharpoonup }{C}=2\stackrel{\rightharpoonup }{r}$$

(25)

where $\stackrel{\rightharpoonup }{A}$ and $\stackrel{\rightharpoonup }{D}$ are coefficient vectors, $t$ is the current iteration, $\stackrel{\rightharpoonup }{X}$ is a position vector, ${\stackrel{\rightharpoonup }{X}}_{best}$ is the position vector of the best current solution (leader whale), $\stackrel{\rightharpoonup }{r}$ is a random number between 0 and 1, and $\stackrel{\rightharpoonup }{a}$ is a variable that decreases linearly from 2 to 0 during the course of the iteration.

4.3.2. Bubble-Net Attacking Method

This behavior simulates the spiral motion of humpback whales surrounding prey and is represented by Equations (26) and (27).

$$\stackrel{\rightharpoonup }{X}\left(t+1\right)=\stackrel{\rightharpoonup }{D\prime }{e}^{b·l}\cos 2\pi l+{\stackrel{\rightharpoonup }{X}}_{best}\left(t\right)$$

(26)

$$\stackrel{\rightharpoonup }{D\prime }=\left|{\stackrel{\rightharpoonup }{X}}_{best}\left(t\right)-\stackrel{\rightharpoonup }{X}\left(t\right)\right|$$

(27)

where $l$ is a random number from −1 to 1 and $b$ is a constant that defines the logarithmic shape.

4.3.3. Searching for Prey

To allow humpback whales to explore more, random searching for prey may be represented by Equations (28) and (29).

$$\stackrel{\rightharpoonup }{X}\left(t+1\right)={\stackrel{\rightharpoonup }{X}}_{rand}\left(t\right)-\stackrel{\rightharpoonup }{A}·\stackrel{\rightharpoonup }{D}$$

(28)

$$\stackrel{\rightharpoonup }{D}=\left|\stackrel{\rightharpoonup }{C}·{\stackrel{\rightharpoonup }{X}}_{rand}\left(t\right)-\stackrel{\rightharpoonup }{X}\left(t\right)\right|$$

(29)

where ${\stackrel{\rightharpoonup }{X}}_{rand}$ is a random search agent from the current population.

4.3.4. Updating Position

Humpback whales exhibit swimming behavior either by swimming towards prey (updating position) or moving around prey (spiral updating position). To simulate this behavior, the WOA assigns a probability equal to 50%, as shown in Equation (26). Updating the position of a humpback whale may be determined by using either Equation (22) ($\left|\stackrel{\rightharpoonup }{A}\right|<1$) or Equation (28) ($\left|\stackrel{\rightharpoonup }{A}\right|\ge 1$).

$$\stackrel{\rightharpoonup }{X}\left(t+1\right)=\left\{\begin{array}{l}{\stackrel{\rightharpoonup }{X}}_{best}\left(t\right)-\stackrel{\rightharpoonup }{A}·\stackrel{\rightharpoonup }{D} if p<0.5\\ \stackrel{\rightharpoonup }{D\prime }{e}^{b·l}\cos 2\pi l+{\stackrel{\rightharpoonup }{X}}_{best}\left(t\right) if p\ge 0.5\end{array}\right.$$

(30)

where $p$ is a random number from 0 to 1.

4.4. Self-Adaptive Control Parameter

In using the PSO to determine an optimum solution, low inertial weight will lead to local trapping, while high inertial weight delays the optimization’s convergence. Jia et al. [45] proposed a self-adaptive parameter control technique by considering each individual’s fitness data to deal with local trapping and slow convergence. This technique was used in this work for better searching as shown in Equation (31).

$$w_i\left(t+1\right)=\left\{\begin{array}{c}{w}_{min}+\left(w_i\left(t\right)-w_{min}\right)\times \frac{f\left(X_i\left(t+1\right)\right)-f_{min}}{f_{avg}-f_{min}} if(r_3<\tau and f\left(X_i\left(t+1\right)\right)<f_{avg}) \\ r_4 if(r_3<\tau and f\left(X_i\left(t+1\right)\right)\ge f_{avg})\\ w_i\left(t\right) if\left(othewise\right)\end{array}\right.$$

(31)

where $\tau$ indicates probabilities to adjust factor $w_i$, $w_{max}$ and $w_{min}$ are the maximum and minimum values of the inertial weights, $r_3$ is a random number from 0 to 1, and $r_4$ is a random number from $w_{min}$ to $w_{max}$.

4.5. Hybrid Particle Swarm and Whale Optimization Algorithm (HPSWOA)

This algorithm was proposed by Trivedi et al. [46]. Due to its constant inertial weight, PSO’s only limitation is its coverage of a narrow search space; hence, it is not suitable for tackling higher-order or complex problems. The HPSWOA can be used to solve this problem because it extracts the quality characteristics of both the WOA and PSO. Since the WOA uses a logarithmic spiral function, it covers a larger area in an uncertain search space during the exploration phase. Therefore, we combined the WOA and PSO to increase their ability to find the solution. Additionally, the WOA accelerates the vector direction toward the optimal value and reduces the computing time. In order to achieve the desired optimal solution, HPSWOA combines the WOA and PSO by having the WOA’s vector ($\stackrel{\rightharpoonup }{X}(t+1)$) replace PSO’s personal best ($Pbes{t}_{i,j}(t)$). The HPSWOA updates position and velocity using Equations (32) and (33). The HPSWOA’s self-adaptation can regulate inertial weight values using Equation (31). The pseudo-code for the HPSWOA for a multi-visit and multi-period dynamic workforce scheduling and routing problem (MMDWSRP) is shown in Algorithms 1 and 2. The flowchart of the HPSWOA is shown in Figure 3.

$$X_{i,j}\left(t+1\right)=X_{i,j}\left(t\right)+V_{i,j}\left(t+1\right)$$

(32)

$$V_{i,j}\left(t+1\right)=w_i\left(t\right)V_{i,j}\left(t\right)+C_1{r}_1\left({\stackrel{\rightharpoonup }{X}}_{i,j}\left(t+1\right)-X_{i,j}\left(t\right)\right)+C_2{r}_2\left(Gbes{t}_j\left(t\right)-X_{i,j}\left(t\right)\right)$$

(33)

4.6. Current Practice Method (CP)

In the CP, jobs will be assigned to technician teams on a first-come first-served (FCFS) basis. If the number of team members and the chosen team’s skill do not meet the work requirement, the next technician team will be chosen instead. The overall details of the CP method are given in Algorithms 1 and 3.

Algorithm 1. Re-schedule for multi-visit and multi-period dynamic workforce scheduling and routing problem (MMDWSRP).

Input: MMDWSRP data, Total cost = 0, $\aleph =\varnothing$ Output: Total cost For each period     Set technician teams’ position at depot     Add job for each sub-system that knows the repair time at start period time to $\aleph$.     While job in period is still in need of completion do         Plan $\aleph$ by selected method         Assign the technician teams by the selected method plan         If new job or harvester relocation then          Update Time          If job in the plan completes before updated time then           Remove the job in $\aleph$ and add to period plan           Update technician teams’ position          End if          Add new job to $\aleph$ or update job’s position         End if       End of while       Calculate period plan fitness by Equation (1)       Total cost += period plan fitness End for

Algorithm 2. Hybrid particle swarm and whale optimization algorithm (HPSWOA).

Input: MMDWSRP data, HPSWOA parameters ($\tau ,b,C_1,C_2,w_{min},w_{max}$), maximum iteration, NP, $\aleph$ Output: plan $\aleph$ for technician teams Randomly generate a set of WOA vectors $i$ ($i$ = 1…NP) (Section 4.1) Randomly generate a set of PSO vectors $i$ ($i$ = 1…NP) (Section 4.1)   While maximum iteration not reached do       For $i$ = 1 to NP         Update position WOA vectors $i$ (Section 4.3)         Update position and velocity PSO vectors $i$ (Section 4.5)         Fitness evaluation PSO vectors $i$ (Section 4.1)         If PSO vectors $i$ fitness < Gbest fitness then           Assign particles $i$ to Gbest           Gbest fitness = PSO vectors $i$ fitness         End if         Update inertial weight PSO vectors $i$ (Section 4.4)       End for       If Gbest fitness < leader whale fitness then         Assign Gbest to leader whale         leader whale fitness = Gbest fitness       End if End of while

Algorithm 3. Current practice method.

Input: MMDWSRP data, $\aleph$ Output: plan $\aleph$ for technician teams Sort jobs in $\aleph$ with ready time in ascending order Sort technician teams with technician team cost in ascending order For sorted jobs     For sorted technician teams       If technician team can service then         Assign the job to the technician team service       Else then         Assign the job to new technician team service       End if     End for End for

5. Computational Results

Problems were solved to illustrate the efficiency and effectiveness of the HPSWOA in respect of the CEC 2017 benchmark functions [47] and the MMDWSRP. In doing so, the self-adaptive control inertial weight probability was set as $\tau$ = 0.1, $w_{min}=0.4,$ and $w_{max}=0.9$ [46]. The control parameters of traditional PSO and the WOA were set as $w$ = 0.729844, $C_1=C_2$ = 1.49618 [48], and $b$ = 0.5 [34]. The numerical experiments in respect of the MMDWSRP consisted of 8 static problem sets (Section 5.2) and 18 dynamic problem sets (Section 5.3) with vector sizes = 100, max iterations = 100, and 10 runs per instance. The proposed algorithms were run using Python on a 2.38 GHz PC, with 8 GB of RAM, for testing and evaluation.

5.1. Optimization of CEC 2017 Benchmark Functions

In this section, the proposed algorithm (HPSWOA) is evaluated using benchmark functions of the well-known CEC 2017. F1 to F9 of the CEC 2017 benchmark functions were applied. Their specifications are detailed in

Table 2. Descriptions of the CEC 2017 benchmark functions.

No.	Group	Description	${\mathit{F}}_{\mathit{i}}^{*}$
F1	Unimodal functions	Shifted and rotated bent cigar function	100
F2		Shifted and rotated Zakharov function	200
F3	Simple multimodal functions	Shifted and rotated Rosenbrock’s function	300
F4		Shifted and rotated Rastrigin’s function	400
F5		Shifted and rotated expanded Scaffer’s F6 function	500
F6		Shifted and rotated Lunacek bi-Rastrigin function	600
F7		Shifted and rotated non-continuous Rastrigin’s function	700
F8		Shifted and rotated Levy function	800
F9		Shifted and rotated Schwefel’s function	900

. From this table, F1 and F2 are unimodal functions, which were employed to evaluate the optimization accuracy and the convergence speed of the proposed algorithms (i.e., the traditional GA, BA, WOA, PSO, and HPSWOA), while F3 to F9 are simple multimodal functions. The proposed algorithm (HPSWOA) was also tested against the traditional genetic algorithm (GA) [49], bat algorithm (BA) [50], WOA, and PSO to solve the CEC 2017 benchmark functions. The computational results for the CEC 2017 benchmark functions, with Dim = 10 in terms of mean and standard deviation of the error values, are shown in

Table 3. Comparative results for the CEC 2017 benchmark functions.

No.	GA		BA		WOA		PSO		HPSWOA
	Mean	Std.	Mean	Std.	Mean	Std.	Mean	Std.	Mean	Std.
F1	9.3355 × 103	5.2101 × 103	1.1542 × 1013	6.3627 × 1012	5.1413 × 103	3.3779 × 103	2.8359 × 103	1.8965 × 103	3.5739 × 103	2.5432 × 103
F2	1.4435 × 103	9.1837 × 102	3.5692 × 107	2.9560 × 107	6.3598 × 10−6	4.1574 × 10−6	0.0000 × 100	0.0000 × 100	0.0000 × 100	0.0000 × 100
F3	5.6102 × 100	4.8250 × 10−1	6.8654 × 105	5.4167 × 105	1.9155 × 100	4.9560 × 10−1	9.8217 × 100	1.4493 × 100	1.1292 × 100	8.7680 × 10−1
F4	9.5142 × 100	2.0899 × 100	5.9447 × 105	2.6542 × 104	2.3524 × 101	5.2721 × 100	1.2314 × 102	2.5796 × 101	1.4057 × 101	4.9011 × 100
F5	2.1875 × 100	1.3174 × 100	6.2376 × 105	5.7474 × 104	1.1509 × 101	8.3052 × 100	5.2498 × 101	1.3300 × 101	1.6095 × 100	8.0240 × 10−1
F6	2.6048 × 101	3.3328 × 100	7.4935 × 105	2.2487 × 104	4.9409 × 101	1.1489 × 101	1.1152 × 102	1.4160 × 10−1	2.0654 × 101	5.4430 × 100
F7	6.5510 × 100	3.4600 × 100	8.4397 × 105	1.6496 × 104	2.3388 × 101	6.4430 × 100	8.3493 × 102	0.0000 × 100	2.1755 × 101	4.4055 × 100
F8	4.3192 × 100	2.8240 × 100	4.7410 × 106	1.7545 × 105	1.6726 × 10−1	1.1135 × 10−1	8.7776 × 102	1.1230 × 101	1.7789 × 10−4	1.2279 × 10−4
F9	3.7047 × 102	2.4671 × 102	3.6634 × 106	3.6975 × 105	6.8409 × 102	2.2967 × 102	1.8222 × 103	4.7987 × 102	2.5487 × 102	1.2286 × 102

. From this table, it can be seen that the HPSWOA and PSO achieved the optimum results in F2. Furthermore, the HPSWOA yielded the best results for F2 to F6 and F8 to F9 (seven out of nine functions). Therefore, the HPSWOA achieved more superior performance in most cases compared to other algorithms (the traditional GA, BA, WOA, and PSO).

5.2. Static Problem

We generated eight static problems (

Table 4. Static problem results for Lingo, WOA, PSO, and HPSWOA.

Ins.	No. of Job	No. of Day	Lingo		WOA			PSO			HPSWOA
			Best	CPU Time (min)	Best	Avg	CPU Time (s)	Best	Avg	CPU Time (s)	Best	Avg	CPU Time (s)
1	10	3	16,847.00	0.03	16,847.00	17,256.60	4.76	16,847.00	17,689.10	7.01	16,847.00	17,084.00	7.26
2	10	7	23,192.00	0.07	23,192.00	23,192.00	2.90	23,192.00	23,192.00	2.84	23,192.00	23,192.00	2.78
3	20	3	20,192.00	6.83	20,192.00	20,665.20	5.00	20,192.00	20,693.60	5.82	20,192.00	20,546.10	6.08
4	20	7	40,042.00	37.33	40,042.00	40,042.00	2.77	40,042.00	40,042.00	2.80	40,042.00	40,042.00	2.87
5	20	15	44,348.00	19.61	44,348.00	44,348.00	3.27	44,348.00	44,348.00	3.38	44,348.00	44,348.00	2.82
6	30	3	31,093.00 1	1440.00	30,697.00	32,056.90	14.68	30,938.00	32,028.00	23.46	30,068.00	30,938.10	23.29
7	30	7	34,853.00 1	1440.00	34,833.00	35,698.40	6.11	35,049.00	36,268.60	9.37	34,833.00	35,446.80	10.00
8	30	15	56,124.00	10.03	56,124.00	56,124.00	3.03	56,124.00	56,124.00	2.84	56,124.00	56,124.00	2.78

¹ Not the optimal solution.

) to compare the proposed methods with Lingo computational software. The penalty cost for each customer was generated in the value interval 15 to 25 for each problem. Job service times and time windows were set in the intervals 30 to 360 and 180 to 600 min, respectively. Overtime cost for each technician team was generated in the value interval 6 to 10 for each problem. The technician teams begin working at the 0th minute and continue until the 480th minute has passed. The waiting time for a job must not exceed 30 min. Overtime hours are limited to 120 min. The labor cost (technician team cost) and outsourcing service were priced in the intervals 1600 to 2000 and 300 to 7200, respectively. The travel time and travel cost were set in the intervals 5 to 100 min and 20 to 400, respectively. In addition, the job requirement for skills and the number of technicians were 2 skills (1 to 3 levels) and 1 to 2, respectively. We tested the traditional PSO, WOA, and HPSWOA performances against the performance of Lingo using 8 test instances. The results for the PSO, WOA, HPSWOA, and Lingo are shown in Table 4. The computational tests compared the best total cost. For instances 1 to 5 and 8, the WOA, PSO, and HPSWOA found that the best total cost equals the optimal solution. For instances 6 and 7, there were a high number of variables, and the MMDWSRP was NP-hard. In addition, the sugar mill only accepts a computational time that is less than 1 day, implying that the planning time for scheduling must not exceed one day. The optimal solution by Lingo could not be reached in one day (1440 min). A statistical test was performed to check if the proposed method was significantly different using the paired t-test, and the results of the test are shown in

Table 5. Statistical test results of the best total cost obtained from Table 4.

	WOA	PSO	HPSWOA
Lingo	0.757	0.444	0.492

. In Table 5, the statistical test revealed that the results of Lingo were not significantly different from those of the WOA, PSO, and HPSWOA. Thus, PSO, the WOA, and the HPSWOA may be used in place of Lingo to solve static problems. The next section outlines the use of PSO, the WOA, and the HPSWOA to solve dynamic problems.

5.3. Dynamic Problem

In actual practice, the harvesters and technicians usually move and operate in different locations during each time period, and there are uncertainties such as breakdowns of the harvesters, together with limited resources, so it is more realistic to consider the problem as a dynamic one. Lund et al. [51] proposed the degree of dynamism (DoD) as a scale from 0 (before planning, all information is readily available) to 1 (before planning, all necessary information is not known) to quantify the proportion of dynamic requests to all requests. In a dynamic problem, we take parameters from the static problem and add the DoD to create the dynamic problem. We generated 18 dynamic problems as shown in

Table 6. The dynamic problem results for CP, the WOA, PSO, and the HPSWOA.

Ins.	No. of Job	No. of Day	DoD	CP	WOA			PSO			HPSWOA
				Best	Best	Avg	CPU Time (s)	Best	Avg	CPU Time (s)	Best	Avg	CPU Time (s)
1	60	7	0.1	60,742.77	56,867.55	58,751.37	43.15	57,888.63	59,843.40	72.82	54,930.89	56,321.80	55.77
2	60	7	0.2	66,165.01	61,926.64	63,313.08	47.46	62,064.82	63,910.11	103.30	60,443.12	61,737.93	75.73
3	60	7	0.3	68,547.37	62,147.38	63,655.08	68.76	63,026.41	64,310.59	83.25	60,583.23	62,002.82	50.22
4	100	7	0.1	86,060.89	85,507.40	87,915.87	245.82	85,040.03	87,839.67	448.50	80,137.68	84,612.61	354.72
5	100	7	0.2	90,143.77	88,425.89	89,749.63	389.28	89,820.82	92,619.14	431.87	85,672.18	87,408.17	429.48
6	100	7	0.3	98,786.54	89,394.79	91,258.24	465.19	90,127.95	93,335.72	517.32	84,678.02	87,069.23	524.65
7	100	15	0.1	115,122.91	102,471.15	103,113.43	56.29	102,673.57	103,419.85	72.03	101,024.25	101,588.70	78.62
8	100	15	0.2	113,841.56	107,658.95	108,713.87	49.96	107,658.95	109,082.14	63.19	106,202.02	108,143.62	55.69
9	100	15	0.3	120,072.12	111,284.17	113,801.47	41.05	111,281.22	113,817.79	62.25	110,977.17	113,055.55	49.14
10	150	7	0.1	157,205.32	138,068.72	141,631.88	510.00	139,019.63	142,985.52	576.20	127,254.70	132,442.84	512.08
11	150	7	0.2	172,360.12	141,335.98	145,410.89	761.75	142,160.74	146,246.04	1317.51	137,330.67	141,279.73	1045.58
12	150	7	0.3	175,128.80	145,251.82	148,517.14	593.44	144,339.02	147,761.72	868.92	136,238.08	140,712.74	837.83
13	150	15	0.1	148,614.33	140,222.92	142,179.03	227.98	141,804.39	146,990.98	324.40	130,367.17	133,773.83	289.68
14	150	15	0.2	153,608.22	140,301.37	142,183.72	261.65	140,290.30	143,741.14	356.78	131,435.74	138,424.81	306.14
15	150	15	0.3	154,851.33	141,343.62	144,553.90	302.06	142,864.80	145,484.64	505.76	134,668.31	141,013.72	439.32
16	150	30	0.1	187,507.74	177,029.34	177,213.14	20.15	177,029.34	177,550.58	30.99	176,457.84	176,572.68	15.99
17	150	30	0.2	195,290.08	186,768.27	188,166.19	12.10	187,051.56	188,778.25	24.12	185,870.04	187,509.73	20.90
18	150	30	0.3	206,263.06	198,529.40	199,915.41	6.89	198,241.63	199,534.03	15.39	197,074.98	198,647.00	6.86

. For the dynamic problems, the number of jobs, number of days, and DoD were set in intervals of 60 to 150 jobs, 7 to 30 days, and 0.1 to 0.3, respectively. The scatter plot of the average computational time in Figure 4 was based on the data from Table 6. From Figure 4, it can be seen that the average computational time of the WOA was the lowest, followed by HPSWOA and PSO. From

Table 7. Statistical test results for the average computational time obtained from Table 6.

	PSO	HPSWOA
WOA	0.008	0.010
PSO		0.016

, it was found that the p-values were all less than 0.05, indicating that the WOA, PSO, and HPSWOA were significantly different in their average computational times. Therefore, it can be concluded that the WOA improves the average computational time of PSO in HPSWOA. Statistical testing was performed to check if the proposed method was significantly different using the paired t-test, and the results of the test are shown in

Table 8. Statistical test results for the best total cost obtained from Table 6.

	WOA	PSO	HPSWOA
CP	0.001	0.001	0.001
WOA		0.018	0.001
PSO			0.001

. In Table 8, the p-values were all less than 0.05, implying that current practice (CP), the WOA, PSO, and the HPSWOA were significantly different in determining the best total cost. The scatter plot of the best total cost (Figure 5) was based on the data from Table 6. From Figure 5, it can be seen that the HPSWOA yielded the best total cost.

Table 9. Average percentage change of the objective function from the best solution obtained from the HPSWOA.

	CP	WOA	PSO
HPSWOA	11.06%	3.47%	3.91%

shows that the average percentage changes of the best total cost obtained from the HPSWOA were 11.06%, 3.47%, and 3.91% when compared with CP, the WOA, and PSO, respectively.

6. Conclusions

In this paper, a new dynamic workforce scheduling and routing problem for the maintenance work of harvesters in a sugarcane harvest operation was introduced. In this problem, technician teams categorized as mechanical, hydraulic, and electrical were assumed to have different skills at different levels to perform services. The jobs were skill-constrained and had time windows. During a working day, a repair request from a sugarcane harvester may arrive, and as time passes, the harvester’s position may shift to other sugarcane fields. This problem was formulated as a multi-visit and multi-period dynamic workforce scheduling and routing problem (MMDWSRP), and this study is the first to address this problem in respect of the WSRP. The objective function of this research is used to minimize the total cost, including technician labor cost, travel cost, late service penalty, overtime, and subcontracting costs. To obtain the harvester maintenance plan, an integer linear programming formulation was developed to solve small problems. For large, practical problems, a hybrid particle swarm and whale optimization algorithm (HPSWOA) was proposed, using an approach based on two well-known methods, namely the WOA and PSO, to solve a multi-visit and multi-period dynamic workforce scheduling and routing problem (MMDWSRP). This approach can be used for maintenance jobs and route planning for each mechanical harvester to be requested by technician teams of relevant sub-systems. The method can be used to determine routes for each individual technician team so that the total cost, consisting of labor cost, traveling cost of technician teams, overtime, penalty, and subcontract costs, is minimized. The proposed algorithms were validated against Lingo computational software using numerical experiments in respect of static problems. Based on the statistical t-test, the WOA, PSO, and HPSWOA were shown to be as efficient as the Lingo computational software.

The proposed method was also tested against the current practice, the traditional WOA, and the traditional PSO, to solve dynamic problems. The experimental results show that the HPSWOA, which extended WOA’s vector to PSO’s personal best vector for a high search exploration advantage, yielded more superior performance than all the other algorithms in terms of solution quality. The total cost of the proposed method was reduced by 11.06%, 3.47%, and 3.91% as compared with the current practice, the traditional WOA, and the traditional PSO, respectively. The HPSWO was also tested against the traditional genetic algorithm (GA), bat algorithm (BA), WOA, and PSO to solve the well-known CEC 2017 benchmark functions. The computational results show that the HPSWOA achieved a superior performance in most cases compared to the GA, BA, WOA, and PSO algorithms.

For future work, more practical limitations can be included in the research. These might include a dynamic workforce constraint, where the number of workforces differs for each day; a route blocking constraint if the maintenance is delayed to the next day resulting a changed route; and maintenance service time as a fuzzy variable of the problem. Additionally, technology adoption such as through Internet of Things (IoT) technology, sensors, Cloud technology, and mobile applications can be developed to create a more intelligent maintenance system for mechanical harvesters. We believe that these can be added to the capability of our study to describe real-life problems and will be a useful expansion.