Solving a multi-objective model of job rotation minimizing the chemical exposure and cost by particle swarm optimization

Abstract

Today in developing countries, perilous chemicals are widely used that have harmful effect on human resources. National and global organizations developed different approaches to control chemicals exposure and their consequences in that one of them is job rotation. Since job rotation is categorized as nondeterministic polynomial-time hardness problem, meta-heuristic methods are used to solve it such as particle swarm optimization (PSO) algorithm. In this article, because job rotation can have different and contradictory objectives, the multi-objective PSO (MOPSO) method with parallel vector evaluated PSO approach is implemented. At first, a new MOPSO algorithm is presented that solves job rotation problems to minimize chemical exposures and is finally compared with non-dominated sorting genetic algorithm. The achievement of this algorithm reduces chemicals exposure in manufacturing processes such as casting, welding, foam injection, plastic injection, which ensure workers’ health.

Introduction

Safety and health problems are caused by excessive exposure to occupational hazards in industrial facilities such as industrial noise, heat, physical workload, and toxic chemicals and substances. To control these problems, a hierarchical approach is recommended: engineering control, administrative control, and the use of personal protection equipment. Among them, the administrative control provides a good balance between implementation cost and effectiveness. In brief, workers change their jobs during the day. In this way, hazardous jobs effect can be shared by many workers, instead of being accumulated by one worker. It is a helpful control in reducing daily hazardous jobs effect of individual workers. ¹

Usually, the workers have to endure and perform repetitive work in long time so increase fatigue and decrease efficiency. Also skill is important factor when operator assigns to work. Job rotation reduces these effects and increases the skill level of workers. ² Having employees with multiple skills is an advantage for organizations, because they not only save time but also save energy and resources. When workers doing monotonous things repeatedly, their passion to work decreases. Job rotation is expected to enable workers to work in different fields; thus, they will be able to gain more experience and knowledge to develop their self-potential. ³

As well as Deljooo and Kheyrkhah have presented a two-objective integer programming model minimizing the hearing impairment effects and back pain of employees, and results were analyzed statistically, and finally they reach best length of job rotation intervals that is 2 h. ⁴

In this research, a multi-objective problem model of job rotation is developed with two objectives minimizing chemical exposure and assigning costs, solved by particle swarm optimization (PSO) algorithm. Because the job rotation problem is a discrete and binary problem, PSO algorithm has to be transformed into a multi-objective binary PSO. At first, PSO, multi-objective PSO (MOPSO), and parallel vector evaluated PSO (VEPSO) are explained, and finally binary MOPSO with the VEPSO approach is proposed and is compared by non-dominated sorting genetic algorithm (NSGA II). ⁵

Literature review

PSO algorithms

In the last decades, unprecedented development in computational intelligence and the introduction of powerful algorithms have helped to understand the nature of problem. One of these powerful algorithms is PSO that is widely used in optimization problem. Much like other optimization algorithms such as Ants Algorithms and Genetic Algorithms, the inner workings of the PSO make use of transition rules to make searches the optimum solution in a search space. The transition rules are based on collective behavior arising out of social interaction of flocks of birds. Since its inception in 1995, PSO has found use in an ever-increasing array of complex, real-world optimization problems. Its simple representation and low number of adjustable parameters make it a popular choice for many problems. ⁶

PSO is an innovative search method, inspired by flying birds to find food, as an evolutionary algorithm ⁷ . For the first time, PSO was used in 1999 to solve multi-objective problem in informal research. ⁸ Currently, more than 25 separate studies have been presented on MOPSO.

Standard PSO is a population-based optimization technique in which a set of particles ${\left\{X_i\right\}}_{i=1}^N$ iteratively find the optimum solution in a search space. Each particle is a candidate solution equivalent to a point in d-dimensional space, so the particles candidate i can be represented as $x_i=\left(x_{i1}, x_{i2}, x_{i3}, \dots , x_{id} \right)$ . The movement of each particle depends on two important factors: ${\mathit{x}}_{\mathit{i}}^{\mathit{b}}$ the best position that the candidate i has found so far and ${\mathit{x}}^{\mathit{g}}$ the global best position found by the whole swarm (all particles). Based on these two factors, each Candidate updates its velocity and position in the iteration (n + 1) as follows:

$\begin{array}{l}{\mathit{v}}_{\mathit{i}}^{\mathit{n}+1}=\mathit{\omega }{\mathit{v}}_{\mathit{i}}^{\mathit{n}}+{\mathit{c}}_1{\mathit{r}}_1\left({\mathit{x}}_{\mathit{i}}^{\mathit{b}}-{\mathit{x}}_{\mathit{i}}^{\mathit{n}}\right)+{\mathit{c}}_2{\mathit{r}}_2\left({\mathit{x}}^{\mathit{g}}-{\mathit{x}}_{\mathit{i}}^{\mathit{n}}\right)\\ {\mathit{x}}_{\mathit{i}}^{\mathit{n}+1}={\mathit{x}}_{\mathit{i}}^{\mathit{n}}+{\mathit{v}}_{\mathit{i}}^{\mathit{n}+1}\end{array}$

1

where ω is the inertia weight, and the parameters c ₁ and ${\mathit{c}}_{2 }$ are positive constants, which balance the influence individual best and global best position. The parameters ${\mathit{r}}_1, {\mathit{r}}_{2 }\in \left(0, 1\right)$ are uniformly distributed random numbers.

PSO for multi-objective problems (MOPSO)

The success of the PSO algorithm in optimizing single-objective problems, especially in continuous spaces, led to develop multi-objective optimization. The multi-objective problems are widely used in various sciences. In such cases, objectives are optimized relative to each other, which mean there are several different solutions to these problems. Our goal in these problems is to find optimal solutions that represent the best answer among the existing objectives.

When PSO is implemented for solving multi-objective optimization problems, it has to be modified. In multi-objective optimization, there isn’t a feasible solution that minimizes all objective functions simultaneously. Therefore, Pareto optimal solutions are solutions that cannot be improved any objectives without regarding to other objectives. In mathematical terms, a feasible solution $\widehat{x}\in X$ is called efficient or Pareto optimal, if there is no other $x\in X$ such that $f\left(x\right)\le f\left(\widehat{x}\right)$ . If $\widehat{x}$ is efficient $f\left(\widehat{x}\right)$ is called nondominated point. If $x_1$ , x₂ ∈ X and $f\left(x_1\right)\le f\left(x_2\right),$ we say $x_1$ dominates $x_2$ and $f\left(x_1\right)$ dominates $f\left(x_2\right)$ . The set of all Pareto optimal solutions $\widehat{x}\in X$ is denoted by XE and called the Pareto optimal set. The set of all nondominated points $\widehat{Y}=f\left(x\right)$ ∈ Y, where $\widehat{x}\in XE$ , is denoted by YN and called the nondominated set. ⁹

According to other evolutionary algorithm, when PSO is developed to multi-objective optimization, the three main issues should be considered:

How to select particles (to be used as leaders) in order to give preference to none dominated solutions over those that are dominated.

How to retain the none dominated solutions found during the search process in order to report solutions that are none dominated with respect to all the past populations and not only with respect to the current one? Also, it is desirable that these solutions are well spread along the Pareto front.

How to maintain diversity in the swarm in order to avoid convergence to a single solution? ¹⁰

Vector evaluated PSO

One of the MOPSO methods is VEPSO that can solve multi-objective problems inspired by the concept and main ideas of the VEGA ¹¹ algorithm. In VEPSO, each group of particles optimizes an objective function and distributes its information to other groups. Its parameters are defined as follows. ¹²

The VEPSO assumes that M swarms $S_1, S_2, \dots , S_M$ of size N aim to optimize simultaneously M objective functions. Each swarm evaluated an objective function. Let $X_i^{\left[j\right]}, V_i^{\left[j\right]}, P_i^{\left[j\right]}$ be the current position, velocity, and the best previous position of the particle i, respectively, and $P_{gb}^i\left(t\right)$ is global best in the swarm j, at a given time t.

In the VEPSO, a neighboring swarm determines the search behavior of a swarm. Specifically, it proposes that the global best position $P_{gb}^s\left(t\right)$ and the best previous position of its particle $P_i^s\left(t\right)$ in the swarm s should be considered for the evaluation of the swarm j velocities. Then the velocity and position of VEPSO’s swarms update as follows:

$\begin{array}{l}{\mathit{V}}_{\mathit{i}}^{\mathit{j}}\left(\mathit{t}+1\right)={\mathit{K}}^{\mathit{j}}.\left[{\mathit{\omega }}_{\mathit{i}}^{\mathit{j}}.{\mathit{V}}_{\mathit{i}}^{\mathit{j}}\left(\mathit{t}\right)+{\mathit{c}}_1^{\mathit{j}}.{\mathit{Rand}}_1.\left({\mathit{p}}_{\mathit{i}}^{\mathit{s}}-{\mathit{s}}_{\mathit{i}}^{\mathit{j}}\left(\mathit{t}\right)\right)+{\mathit{c}}_2^{\mathit{j}}.{\mathit{Rand}}_2.\left({\mathit{p}}_{\mathit{g}\mathit{b}}^{\mathit{s}}-{\mathit{s}}_{\mathit{i}}^{\mathit{j}}\left(\mathit{t}\right)\right)\right]\\ {\mathit{x}}_{\mathit{i}}^{\mathit{j}}\left(\mathit{t}+1\right)={\mathit{x}}_{\mathit{i}}^{\mathit{j}}\left(\mathit{t}\right)+{\mathit{V}}_{\mathit{i}}^{\mathit{j}}\left(\mathit{t}+1\right)\end{array}$

2

where j is the number of swarms, ${\mathit{c}}_{1 }$ and ${\mathit{c}}_2$ are the cognitive and the social parameters, respectively, Rand₁ and Rand₂ are random numbers uniformly distributed within [0, 1], and the superscript [s] represents the ring migration defined as

$\mathit{s}=\{\begin{array}{c}M, if j=1,\\ j-1, if j=2,\dots ,M,\end{array}$

3

The inertia weighting function for the velocity of particle i is defined by the inertia weight approach.

${\mathit{w}}_i={\mathit{w}}_{\mathit{max}}-\frac{{\mathit{w}}_{\mathit{max}}-{\mathit{w}}_{\mathit{min}}}{{\mathit{t}}_{\mathit{max}}}.\mathit{t}$

4

That $t_{\mathit{max}}$ is the maximum number of iterations and t is the current number of iterations.

The role of the inertia weighting function is critical to PSO’s convergence that controls the influence of the previous history of the velocities on the current one. Also, the inertia weighting function balanced the global and local exploration abilities of the swarms. Another parameter that guarantees the convergence of the PSO algorithm is constriction factor $\mathit{K}$ :

$\mathit{K}=\frac{2}{|2-\mathit{\phi }-\sqrt{{\mathit{\phi }}^2-4\mathit{\phi }}|},\mathit{\phi }={\mathit{c}}_1+{\mathit{c}}_{2,}\mathit{\phi }>4$

5

As well as there are some constraints to system convergence such as the positions and the velocities of the particle i of the swarm j. The positions are limited by the minimum and maximum positions expressed by vectors:

$\left[{\mathit{S}}_{\mathit{i}}^{\left[\mathit{j}\right].\mathit{min}},{\mathit{S}}_{\mathit{i}}^{\left[\mathit{j}\right].\mathit{max}}\right],\left(\mathit{j}=1,2,\dots ,\mathit{M}\right),\left(\mathit{i}=1,2,\dots ,\mathit{N}\right)$

6

Also the velocities of the particle i of the swarm j are limited by:

$\left[-{\mathit{V}}_{\mathit{i}}^{\left[\mathit{j}\right].\mathit{min}},{\mathit{V}}_{\mathit{i}}^{\left[\mathit{j}\right].\mathit{max}}\right],\left(\mathit{j}=1,2,\dots ,\mathit{M}\right),\left(\mathit{i}=1,2,\dots ,\mathit{N}\right)$

7

Where the vector of maximum velocities is defined as:

${\mathit{v}}_{\mathit{i}.\mathit{l}\mathit{g}}^{\left[\mathit{j}\right],\mathit{max}}=\frac{{\mathit{s}}_{\mathit{i}.\mathit{l}\mathit{g}}^{\left[\mathit{j}\right],\mathit{max}}-{\mathit{s}}_{\mathit{i}.\mathit{l}\mathit{g}}^{\left[\mathit{j}\right].\mathit{min}}}{\mathit{N}\mathit{r}},\left(\mathit{j}=1,2,\dots ,\mathit{M}\right),\left(\mathit{i}=1,2,\dots ,\mathit{N}\right),\left(\mathit{l}\mathit{g}=1,2,\dots ,\mathit{n}\right)$

8

That $Nr$ is an important parameter in the parallel VEPSO algorithm. A low $Nr$ facilitates global exploration (searching new areas), while a high one facilitates local exploration (fine-tuning of the current search area). Not only a suitable amount for $Nr$ usually provides a balance between global and local exploration abilities, but also it reduces the number of iterations to locate the optimum solution. ¹³

Binary PSO

The PSO was originally developed for continuous-valued spaces but many problems have discrete spaces so binary PSO must be defined. In binary PSO, each value particle is 0 or 1 that can be changed from 1 to 0 or vice versa. Also the velocity of a particle defined as probability variable. The particle’s personal best and global best is updated as continuous version.

Two velocity vectors for each particle are introduced as $\overrightarrow{V_i^0}$ and $\overrightarrow{V_i^1}$ . Vector $\overrightarrow{V_i^{0 }}$ is probability of the particle that changes to zero, while $\overrightarrow{V_i^1}$ is the probability of the particle that changes to one so the velocity probability of bit j of particle i is simply defined as follows:

${\mathit{v}}_{\mathit{i},\mathit{j}}^{\mathit{C}}=\{\begin{array}{c}{\mathit{v}}_{\mathit{i},\mathit{j}}^1 \mathit{if} {\mathit{x}}_{\mathit{i},\mathit{j}}=0\\ {\mathit{v}}_{\mathit{i},\mathit{j}}^0 \mathit{if} {\mathit{x}}_{\mathit{i},\mathit{j}}=1\end{array}$

9

Also, $\overrightarrow{ V_i^0}$ and $\overrightarrow{V_i^1}$ are updated as follows:

$\begin{array}{c}\mathit{I}\mathit{f}{\mathit{P}}_{\mathit{i}\mathit{best}}^{\mathit{j}}=1\mathit{Then}{\mathit{d}}_{\mathit{i}\mathit{j},1}^1={\mathit{c}}_1{\mathit{r}}_1\mathit{and}{\mathit{d}}_{\mathit{i}\mathit{j},1}^0=-{\mathit{c}}_1{\mathit{r}}_1\\ \mathit{I}\mathit{f}{\mathit{P}}_{\mathit{i}\mathit{best}}^{\mathit{j}}=0\mathit{Then}{\mathit{d}}_{\mathit{i}\mathit{j},1}^0={\mathit{c}}_1{\mathit{r}}_1\mathit{and}{\mathit{d}}_{\mathit{i}\mathit{j},1}^1=-{\mathit{c}}_1{\mathit{r}}_1\\ \mathit{I}\mathit{f}{\mathit{P}}_{\mathit{g}\mathit{best}}^{\mathit{j}}=1\mathit{Then}{\mathit{d}}_{\mathit{i}\mathit{j},2}^1={\mathit{c}}_2{\mathit{r}}_2\mathit{and}{\mathit{d}}_{\mathit{i}\mathit{j},2}^0=-{\mathit{c}}_2{\mathit{r}}_2\\ \mathit{I}\mathit{f}{\mathit{P}}_{\mathit{g}\mathit{best}}^{\mathit{j}}=0\mathit{Then}{\mathit{d}}_{\mathit{i}\mathit{j},2}^0={\mathit{c}}_2{\mathit{r}}_2\mathit{and}{\mathit{d}}_{\mathit{i}\mathit{j},2}^1=-{\mathit{c}}_2{\mathit{r}}_2\end{array}$

10

where $d_{ij}^1, d_{ij }^0$ are two temporary values, also $r_{1 }$ and r₂ are two random variables in the range of $\left(0,1\right)$ which are updated each iteration, and c₁, c₂ are two fixed variables which are determined by the user.

$\begin{array}{l}{\mathit{v}}_{\mathit{i}\mathit{j}}^1=\mathit{w}{\mathit{v}}_{\mathit{i}\mathit{j}}^1+{\mathit{d}}_{\mathit{i}\mathit{j},1}^1+{\mathit{d}}_{\mathit{i}\mathit{j},2}^1\\ {\mathit{v}}_{\mathit{i}\mathit{j}}^0=\mathit{w}{\mathit{v}}_{\mathit{i}\mathit{j}}^0+{\mathit{d}}_{\mathit{i}\mathit{j},1}^0+{\mathit{d}}_{\mathit{i}\mathit{j},2}^0\end{array}$

11

where w is the inertia term. In this algorithm, if the bit j of the global best variable is 0 or if the bit j of the corresponding personal best variable is 0, then the velocity $\overrightarrow{V_i^0}$ is increased. At the same time, the probability of change to 1 is decreased. In addition, if the bit j in the global best variable is 1, the velocity $\overrightarrow{V_i^1}$ increases and $\overrightarrow{V_i^0}$ decreases. In this approach, previous direction of a bit uses to find next direction that can change to 1 or 0; therefore, the velocity of particles $\overrightarrow{V_i^0}$ and $\overrightarrow{V_i^1}$ updates as equation (11). Finally, sigmoid function (9) is normalized to calculate the next particles position as follows:

${\mathit{x}}_{\mathit{i}\mathit{j}}\left(\text{t+1}\right)\text{=}\{\begin{array}{c}{\overline{\mathit{x}}}_{\mathit{i}\mathit{j}}{\text{(t) if r}}_{\text{ij}}\text{ <} {\mathit{v}}_{\mathit{i}\mathit{j}}^{,}\\ {\mathit{x}}_{\mathit{i}\mathit{j}}{\text{(t) if r}}_{\text{ij}}\text{ >} {\mathit{v}}_{\mathit{i}\mathit{j}}^{,}\end{array}$

12

where ${\overline{X}}_{ij}$ is binary variable, and it is opposite of $X_{ij}$ , so that if $X_{ij}=0,$ then ${\overline{X}}_{ij}=1$ and if $X_{ij}=1,$ then ${\overline{X}}_{ij}=0$ and $r_{ij}$ is a random number between 0 and 1.

The meaning of the parameters in the velocity equation (11) is exactly like continuous PSO. The inertia maintains the previous direction of bits. Also, the velocity is the rate of change in the particle’s position. The large value of maximum velocity causes random search, while the small value of maximum velocity causes the particle to move less.

The problem definition

In this article, based on VEPSO approach, a multi-objective binary algorithm for solving MOPSO problems is proposed by making the necessary changes in the single-objective binary algorithm (“Binary PSO” section). Then, the mathematical model of the job rotation is implemented in a case study minimizing chemical exposure and costs. In this model, two objective functions are used naming ${\text{obj}}_1$ and ${\text{ obj}}_2$ . After each objective is separately solved, the optimal value is shown by ${\text{obj}}_1^{\text{*}}$ and ${\text{ obj}}_2^{\text{*}}$ . Then, an L_P-metric method is formulated as follows:

$Z=\text{min}\left[w_1.\frac{\left(ob{j}_1^{\text{*}}-ob{j}_1\right)}{ob{j}_1^{\text{*}}}+w_2.\frac{\left(ob{j}_2^{\text{*}}-ob{j}_2\right)}{ob{j}_2^{\text{*}}}\right]$

13

That $w_1$ and $w_2$ are the weight values considered by the decision maker (DM). A brief discussion on L_P-metric method is given in this section. L_P-metric method is usually discussed in multi-objective design making, which is an optimization technique combining multiple objectives into a single objective. The weighted L_P distance measure of any solution x from the ideal solution $f_i\left(x^{*}\right)$ can be minimized as follows ¹⁴ :

$L_P={\left\{{\displaystyle \sum }_{i=1}^k{W}_i{\left(f_i\left(x^{\text{*}}\right)-f_i\left(x\right)\right)}^P\right\}}^{\frac{1}{P}}$

14

where $W_i$ is a nonnegative weight that is assigned to each objective function by DM and the amount $P$ indicates the importance of each objective function deviation from its ideal value. When $P=\infty$ is considered, the L_P changes as follows:

$L_P={\displaystyle \sum }_{i=1}^k{W}_i\frac{f_i\left(x^{\text{*}}\right)-f_i\left(x\right)}{f_i\left(x^{\text{*}}\right)}$

15

It is showed when L_P method is used and $1\le \text{P}\le \infty$ and $\text{w}>$ 0 so all solutions are efficient solutions. On the other hand, if one considers constrained L_P problems, only ${\text{L}}_{\infty }$ leads to Pareto optimal solutions. ¹⁴

Proposed algorithm

As detailed in the previous chapters, the PSO algorithm is inspired by moving particles toward the target which speed and direction of movement are determined by the best position of each particle and the best position of all particles. According to formulas (1) and (3) for the PSO continuous algorithm, we have:

$\begin{array}{l}{\mathit{V}}_{\mathit{i}}\left(\mathit{t}+1\right)=\mathit{K}\text{*}\left[{\mathit{\omega }}_{\mathit{t}}\text{* }{\mathit{V}}_{\mathit{i}}\left(\mathit{t}\right)+{\mathit{c}}_1{\mathit{r}}_1\left({\mathit{p}}_{\mathit{i}}\mathit{best}\left(\mathit{t}\right)-{\mathit{X}}_{\mathit{i}}\left(\mathit{t}\right)\right)+{\mathit{c}}_2{\mathit{r}}_2\left(\mathit{g}\mathit{best}\left(\mathit{t}\right)-{\mathit{X}}_{\mathit{i}}\left(\mathit{t}\right)\right)\right]\\ {\mathit{X}}_{\mathit{i}}\left(\mathit{t}+1\right)={\mathit{X}}_{\mathit{i}}\left(\mathit{t}\right)+{\mathit{V}}_{\mathit{i}}\left(\mathit{t}+1\right)\end{array}$

16

where j is the number of swarms; $c_1$ and $c_{2 }$ are the cognitive and the social parameters, respectively, and r₁and $r_2$ are random numbers uniformly distributed within $\left[0,1\right].$ As well as $p_{ibest}$ is the best previous position of the particle i and $g_{ibest}$ is the global best in the swarm j at time t.

Factors ${\text{ ω}}_{\text{t}}$ and K are the most important parameters to the algorithm convergence. The inertia weight ${\text{ω}}_{\text{t}}$ is a linear function that decreases linearly at each event. The function of this factor is as follows:

${\mathit{w}}_{\mathit{t}}={\mathit{w}}_{\mathit{max}}-\frac{{\mathit{w}}_{\mathit{max}}-{\mathit{w}}_{\mathit{min}}}{{\mathit{t}}_{\mathit{max}}}.\mathit{t}$

17

Usually ${\omega }_{\text{max}}$ is considered 1, and ${\omega }_{\text{min}}$ is considered 0.1. Design of experiment (DOE) is used to determine the algorithm formula coefficients. During problem-solving, big inertia weight improves global searches and small inertia weight improves local searches. In the initial repetition of the problem-solving, big inertia weight causes global search and in the last repetition, small inertia weight causes local search.

In order to guarantee the convergence of the PSO algorithm, the constriction factor K is adopted:

$\mathit{K}=\frac{2}{|2-\mathit{\phi }-\sqrt{{\mathit{\phi }}^2-4\mathit{\phi }}|},\mathit{\phi }={\mathit{c}}_1+{\mathit{c}}_2,\mathit{\phi }>4$

18

For using the two above factors in the PSO discrete algorithm, the changes in formulas (11) apply as follows:

$\begin{array}{l}{\mathit{v}}_{\mathit{i}\mathit{j}}^0=\mathit{K}\text{*}\left[{\mathit{w}}_{\mathit{t}}\text{*}{\mathit{v}}_{\mathit{i}\mathit{j}}^0+{\mathit{d}}_{\mathit{i}\mathit{j},1}^0+{\mathit{d}}_{\mathit{i}\mathit{j},2}^0\right]\\ {\mathit{v}}_{\mathit{i}\mathit{j}}^1=\mathit{K}\text{*}\left[{\mathit{w}}_{\mathit{t}}\text{*}{\mathit{v}}_{\mathit{i}\mathit{j}}^1+{\mathit{d}}_{\mathit{i}\mathit{j},1}^1+{\mathit{d}}_{\mathit{i}\mathit{j},2}^1\right]\end{array}$

19

Algorithm

Algorithm assumptions

The work shift is 8 h and each shift is up to 2 h.

If the shifts exceed 8 h, the threshold limit value (TLV) will not change.

In any case, at least one person from each specialized level must be present.

The number of workers should be higher than the minimum number of workers per job.

During a day, no worker should be unemployed.

Parameters and variables of the problem

Before the model was explained, parameters and variables are defined in Table 1.

Table 1.

Variables and indicators of problem.

j=1, 2,…, J	Workgroups
t = 1,…, S	Number of the period in a day
d = 1,…, D	Number of working days in a planning period
K	The number of specialized levels of work
lk	Number of operators at level k
$l_k\left(k\right)=1,\hspace{1em}\dots , l_k$	The index of workers at level k
k = 1,…, L	Total number of workers
i = 1,…, I	Number of chemicals
Nj	The minimum number of workers in the working group j
$TL{V}_i$	Amount of exposure to chemicals i during one 8-h day in accordance with the ACGIH standard
$E_{ij}$	Exposure to chemicals i in the workgroup j
$C_{ij}$	The cost of assigning a level k to work j
D	Harmful factors index
$D_{Max}$	Maximum Chemical Contamination Index
$C_{Max}$	Maximum cost job allocation
	If the worker l with at level k at the period t from workday d is assigned to the work j
$X_{kjtd}\hspace{1em}$ 10	Otherwise

Presenting the model

The general form of the problem is formulated as follows:

$Min{D}_{Max}$

20

$Min{C}_{Max}={\displaystyle \sum }_{d=1}^D{\displaystyle \sum }_{t=1}^S{\displaystyle \sum }_{j=1}^J{\displaystyle \sum }_{k=1}^L{C}_{jk}{X}_{k,j,t,d}$

21

s.t:

${\displaystyle \sum }_{j=1}^J{X}_{k, j, t, d}\le 1\forall k, t, d$

22

${\displaystyle \sum }_{j=1}^J{\displaystyle \sum }_{t=1}^S{X}_{k, j, t, d }=S \forall d, k$

23

${\displaystyle \sum }_{k=1}^L{X}_{k,j,t,d}\ge N_j \forall j, t, d$

24

${\displaystyle \sum }_{t=1}^S{\displaystyle \sum }_{j=1}^J{\displaystyle \sum }_{i=1}^I\frac{E_{ji} . X_{k,j,t,d}}{TL{V}_i.S} < D_{\max } \forall k, d$

25

${\displaystyle \sum }_{k\in l\left(k\right)}{X}_{k,j,t,d}\ge 1 \forall j, t, d$

26

$X_{ k, j,t,d}=1 \text{\hspace{0.17em}or\hspace{0.17em}} 0$

27

Following, each objectives and constraints are explained:

The first objective function (20) shows that chemical exposures are minimized.

The second objective function (21) shows that assigning workers cost is minimized.

Constraint (22) shows that each worker is assigned in one workgroup at each period.

Constraint (23) shows that no worker can be unemployed a day.

Constraint (24) shows that the minimum number of workers is assigned in each workgroup at each period of a day.

Constraint (26) shows that at least one worker must be present at each level.

The last constraint is the variable X that can contain 0 or 1 values.

In order to calculate $D_{\max }$ , the following harmful factors index is used ¹⁵

$D=\frac{C1}{TLV1}+\frac{C2}{TLV2}+\frac{C3}{TLV3}+\dots +\frac{Cn}{TLVn}$

28

where $C_n is$ the exposure value of the chemical and TLV ( $TL{V}_{n })$ is the permitted amount of exposure to that chemical over an 8-h working day (according to ACGIH).

If D is less than 1, chemicals exposure is permitted, but if D is more than 1, chemicals exposure is unacceptable and corrective actions should be defined to reduce the chemicals exposure. ¹⁵

In the ACGIH standard, the amount is expressed in parts per million (ppm) or milligram per cubic meter (mg/m³) which according to the below formula can be converted to any other:

$x_{ppm}=y_{\frac{mg}{m3}}\text{*}\left(24.45\right)\text{*}\left(MW\right)$

29

MW is each chemical molecular weight calculated by the following website. ¹⁶ As well as following site is used to convert from ppm to mg/m³. ¹⁷ Since in each shift the chemical exposure duration is 2 h, the maximum exposure is calculated as follows:

$\frac{T=2}{8}=\frac{TL{V}_i}{TL{\text{V}}_{i^{\prime }}}$

30

The value of T is calculated as:

$T= \frac{8}{S}$

31

That S is the number of work rotations per day. As a result, we have for $TL{V}_i$ :

$TL{\text{V}}_{i^{\prime }}=TL{V}_i\text{*}S$

32

The chemical exposure level per day in a workgroup is calculated as follows:

${\displaystyle \sum }_{i=1}^I\frac{E_{j,i}}{TL{\text{V}}_{i^{\prime }}}\to {\displaystyle \sum }_{i=1}^I\frac{E_{j,i}}{TL{V}_i\text{*}S}$

33

That $E_{ji}$ is the exposure of each worker to chemical i in workgroup j.

Finally, the amount of chemical exposures of each worker is calculated as follows (Constraint (25)):

${\displaystyle \sum }_{t=1}^S{\displaystyle \sum }_{j=1}^J{\displaystyle \sum }_{i=1}^I\frac{E_{ji}\text{* }{X}_{k,j,t,d}}{TL{V}_i\text{*}S}<D_{\max }\forall k,d$

34

The second objective function (21) is minimized the costs associated with assigning workers in workgroups at each period of scheduling time.

$Min{\displaystyle \sum }_{d=1}^D{\displaystyle \sum }_{t=1}^S{\displaystyle \sum }_{j=1}^J{\displaystyle \sum }_{k=1}^L{C}_{jk}{X}_{k,j,t,d}$

35

If worker k is assigned in workgroup j at period t on day d, then the cost is equal to $C_{jk}$ that sum of these costs should be minimized.

The costs of allocating a low-level worker on a job include the cost of losses, rework, training during work, the failure of the machine, and so on, and the cost allocation with high-level worker on a job includes missed opportunity cost that worker proficiency can be used in better job.

Before solving the problem with the PSO algorithm, the solvability of mathematical model is examined by GAMS software.

Numerical study

Introduction

In this chapter, the mathematic model is solved, and the results are compared. Because of the lacks of parameters information, these parameters are determined by the DOE method. DOE helps researchers do experiments that discover phenomena and discover the truth about a system, as well as DOE determines the input variables that have the most impact on system output. ¹⁸ In accordance with the factorial method of DOE, after 32 iterations, upper and lower limits set for four parameters ( ${\text{w}}_{\text{max}}\mathrm{و,} {\text{w}}_{\text{min}},{\text{ c}}_2,{\text{c}}_1)$ of the algorithm, and the results are shown in Table 2. Also on, identify a method and solve the problem then the problem results are presented and the particle’s behavior is examined.

Table 2.

Parameter values of PSO algorithm, ${\text{w}}_{\text{max}}\mathrm{و,} {\text{w}}_{\text{min}},{\text{ c}}_2,{\text{c}}_1.$

${\mathit{w}}_{\mathit{max}}$	${\mathit{w}}_{\mathit{min}}$	C1	C2
1	0.01	3.05	1.05

PSO: particle swarm optimization.

As well as the related variables in the NSGA II are: the crossover probability is equal to 0.9 and the mutation probability is equal to 0.1. The results of MOPSO and NSGA II can be compared when the function evaluations parameter in both methods is equal. Function evaluations parameter for MOPSO and NSGA II can be determined as follows ¹⁹ :

$\text{Function\hspace{0.17em}Evaluations MOPSO}=\text{Swarm Size }\times \text{ Maximum Iteration}$

$\text{Function\hspace{0.17em}Evaluations NSGA II}=\text{Population Size }\times \text{ Maximum Generation}$

Determines the parameters of the problem

Following parameters are determined to solve the model: the number of particles is 50; the number of iterations is 50, and the dimensions of the problem are 5000. Dimensions are calculated as below:

Calculate the space dimensions of the problem

A solution is the particle position in decision variable space. In this case, there are 50 workers, 5 workgroups (welding, foam injection, sewing, patting, and assembling), and 4 periods (each period is 2 h) during a 5-day planning period (1 week’s work).

A feasible solution is a worker in a shift during a planning period assign to a task that depends on the amount of chemical exposure, as worker k in shift t is assigned to workgroups j from workday d. Each of 50 employees can be assigned to one of the five workgroups in one of the four shifts of the 5-day planning period. If a job is assigned to a worker in a shift on a day, the variable value should be 1; otherwise, it is 0.

According to Table 3, there are 20 modes in which a worker can be assigned to a work so the number of assignments for 50 workers per day is equal to 50 × 20 = 1000. As a result, the solved space dimensions of feasible solution for a 5-day planning period are 5000 that 1000 component values are 1 and others (4000 dimensions) are 0.

Table 3.

All modes that a worker assigns to one of five workgroups during 1 day.

One day
One worker
Shift 1	Shift 1	Shift 1	Shift 1

As well as seen in “Presenting the model” section, when the level of worker skill increases, the missed opportunity cost is increased too because total costs such as rework and component waste are reduced. The costs are related to level of the worker skills that have been extracted from quality costs of the case study organization.

Results and analysis

Optimal problem-solving

As shown in “Calculate the space dimensions of the problem” section, the dimensions of the mathematical model are 5000. In this section, it is shown that the optimal solution of the model at large solved space dimensions (of 100,000) is impossible (see Table 5). When one or more parameters increase, the solved space dimension is increased. For example, if there are 100 day, 100 workers, 10 work types, and 4 shifts in a day, the solved space dimension is equal to 100 × 100 × 10 × 4 = 400,000.

According to “The problem definition” section, first, each objective is solved in a single-objective problem in order to get the optimal value.

$Ob{j}_{\text{1}}\text{* = Optimal amount of the first objective }\left(\text{reduction of exposure to chemicals}\right)\text{ = 1.254.}$

$Ob{j}_{\text{2}}\text{* = Optimal amount of second objective }\left(\text{reduction of job rotation costs}\right)\text{ = 2345.}$

If w₁ = 0.4 and w₂ = 0.6 (determined by the DM), after optimizing the optimal value is Z = 0.0234, so the value of two objective functions is Obj₁ = 1.297 and Obj₂ = 2365 (Table 4).

Table 4.

Objective value separately and concurrently.

	Value z	Value Obj2	Value Obj1
Obj 1	0.0113	2384	1.254
Obj 2	0.0819	2345	1.5
LP-metric	0.0234	2365	1.297

If each object is optimally optimized alone, the value of other objective increases. When the problem is solved by the L_P-metric method, both functions are optimized simultaneously so the optimal value of this function less than optimal value of each objective functions separately. This operation can be repeated at the same time when the dimensions of the problem increase, but increasing the dimensions causes reduction in software optimization accuracy, or software is not capable of solving. These results are shown in Table 5.

Table 5.

Objective value and problem-solving time with increasing its variables.

Problem-solving time (min)	Obj2	Obj1	Number of dimensions
5	487	1.353	1000
8	944	1.353	2000
15	1884	1.351	4000
16	2360	1.298	5000
25	4703	1.485	10,000
36	7036	1.648	15,000
48	9382	1.648	20,000
60	11729	1.753	25,000
The software is not able to solve the problem in this dimension			100,000
			200,000
			300,000

Solve the model with the PSO

As seen in the previous section, the problems in large dimensions have not the optimal solution so meta-heuristic methods are used such as PSO. To compare the accuracy of solutions, the problem is solved by binary PSO and NSGA II, and solutions compared with the optimal solution. The results of 10-time ruining of the proposed binary PSO and the NSGA II with 50 times repetitive are shown in Tables 6 and 7. After 10 times, the number of closest solution to optimal solution of proposed binary PSO is 4, while only one time the NSGA II has the closest solution to optimal solution and other solutions are not feasible (do not apply in problem constraints). Also the time to solve the proposed binary PSO is greater than the NSGA II.

Table 6.

Results of 10 times the run of the NSGA II algorithm.

Best Obj2	Best Obj1	Best Z	Feasible solution
2348	1.861	0.200	No	1
2342	2.008	0.246	No	2
2360	1.861	0.203	Yes	3
2385	1.883	0.216	No	4
2373	1.936	0.230	Yes	5
2360	1.881	0.209	No	6
2365	1.500	0.088	No	7
2340	1.640	0.127	No	8
2371	1.547	0.105	No	9
2390	1.450	0.079	No	10

NSGA II: non-dominated sorting genetic algorithm.

Table 7.

Results of 10 times the run of the PSO algorithm.

Best Obj2	Best Obj1	Best Z	Feasible solution
2360	1.881	0.209	No	1
2390	1.450	0.079	No	2
2344	1.758	0.166	No	3
2360	1.861	0.203	Yes	4
2360	1.861	0.203	Yes	5
2369	1.766	0.175	No	6
2360	1.861	0.203	Yes	7
2324	1.640	0.123	No	8
2360	1.861	0.203	Yes	9
2385	1.883	0.216	No	10

PSO: particle swarm optimization.

Particle behavior survey

According to Figures 1 and 2, both algorithms converge optimal solution (Obj₂ = 2360, Obj₁ = 1.891); also according to Table 8, the number of Pareto optimal of the proposed PSO binary algorithm is greater than the NSGA II algorithm. Because the proposed algorithm has complex, accurate, and high computations, the problem-solving time is high so it works more accurate than the NSGA II algorithm.

Figure 1.

PSO algorithm converges to the optimal answer. PSO: particle swarm optimization.

Figure 2.

NSGA II algorithm converges to the optimal answer. NSGA II: non-dominated sorting genetic algorithm.

Table 8.

Number of Pareto optimal.

Time to solve (min)	Number of Pareto optimal
30	10	PSO
20	3	NSGA II

PSO: particle swarm optimization; NSGA II: non-dominated sorting genetic algorithm.

Particle movement behavior in initial and final solution

In this section, from the beginning of the problem-solving for reaching the optimal point, the particle movement behavior is examined. For example, in the initial solution with ${\text{obj}}_1=2.02$ ,، ${\text{obj}}_2=3060$ , and Z = 0.434, in the first shift of the first working day, the worker numbers {6, 8, 10, 40, 31, 22, 11, 41, 44, 48} were assigned to workgroup 1 as $x_{6,1,1,1}=1$ , and other variables for worker 6 in the first shift of the first working day are 0 as $x_{6,2,1,1}=0,x_{6,3,1,1}=0, x_{6,4,1,1}=0, x_{6,5,1,1}=0$ . Whereas in the final solution with ${\text{ obj}}_1=1.861$ , ${\text{obj}}_2=2360$ , and Z = 0.203, the worker numbers {2, 8, 10, 11, 12, 20, 21, 38, 42, 49} were assigned to workgroup 1 as $x_{2,1,1,1}=1$ , and other variables for worker 2 in the first shift of the first working day are 0 as $x_{2,2,1,1}=0, x_{2,3,1,1}=0,x_{2,4,1,1}=0, x_{2,5,1,1}=0$ . The results of assigning 50 workers to 5 working days on the first day are shown in Table 9 and other ${\text{x}}_{\text{k},\text{ j},1,1}=0$ that the number of these variables in the first shift of the first working day is = 200. As well as, in the initial solution, worker 6 assigns to workgroup 1 but in the final solution assigns to workgroup 5.

Table 9.

Assign 5 types of work to 50 workers in the first shift of the first working day.

Workgroup j	Worker k		Workgroup j	Worker k		Workgroup j	Worker k
	Initial solution	Final solution		Initial solution	Final solution		Initial solution	Final solution
1	6	2	3	18	23	5	4	4
1	8	8	3	20	24	5	9	6
1	10	10	3	27	26	5	12	7
1	11	11	3	32	32	5	15	14
1	22	12	3	35	34	5	19	15
1	31	20	3	36	35	5	21	16
1	40	21	3	37	36	5	23	17
1	41	38	3	39	41	5	24	18
1	44	42	3	49	48	5	25	19
1	48	49	4	2	5	5	26	27
2	7	13	4	3	9	5	29	28
2	14	22	4	13	29	5	34	31
2	17	25	4	16	33	5	43	44
2	30	30	4	28	39	5	47	46
2	38	37	4	33	43
2	42	40	4	45	45
2	50	47	4	46	50
3	5	1	5	1	3

Conclusion

In this case, after the multi-objective model parameters and variables have been determined, it is modeled and solved. The model includes constraints such as assigning each worker to only one type of work in each shift, the whole chemical exposure of the worker in a working day, and assigning at least one person from each level in each of the workgroups. Results of case study can be categorized as follows:

All 50 workers are assigned to existing workgroups so that chemicals exposure and workflow costs are minimized.

A job rotation program is planned assigning 50 workers to 5 workgroups in 5 working days with 4 shifts (each shift is 2 h) that there are no unemployed workers in each work shift.

The number of workers of each level in a workgroup is determined according to minimum number required workers in each of workgroups.

As well as the number of Pareto optimal of the proposed PSO binary algorithm is greater than the NSGA II algorithm. Thus, this algorithm works more accurate than the NSGA II algorithm, but it takes more time and space. Because the proposed algorithm has complex, accurate, and high computations, the problem-solving time is high so it works more accurate than the NSGA II algorithm.