Probability Mechanism Based Particle Swarm Optimization Algorithm and Its Application in Resource-Constrained Project Scheduling Problems

In this paper, a new probability mechanism based particle swarm optimization (PMPSO) algorithm is proposed to solve combinatorial optimization problems. Based on the idea of traditional PSO, the algorithm generates new particles based on the optimal particles in the population and the historical optimal particles in the individual changes. In our algorithm, new particles are generated by a specially designed probability selection mechanism. We adjust the probability of each child element in the new particle generation based on the difference between the best particles and the elements of each particle. To this end, we redefine the speed, position, and arithmetic symbols in the PMPSO algorithm. To test the performance of PMPSO, we used PMPSO to solve resource-constrained project scheduling problems. Experimental results validated the efficacy of the algorithm.


Introduction
Particle Swarm Optimization (PSO) is an evolutionary computational technique proposed by Kennedy and Eberhart [1] in 1995 for continuous optimization problems.Since the PSO algorithm only needs to adjust a small number of parameters and it is easy to implement, theories related to PSO algorithm has been greatly developed in the last two decades [2].In the literature, there are many successful PSO applications for single-or multiobjective optimization [3,4].
Although PSO was developed for continuous optimization problems initially, some work focused on its discrete version (discrete PSO and DPSO).Kennedy and Eberhart proposed a discrete binary version of particle swarm algorithm in 1997 [5], where the trajectories of the particles are changes in the probability that a coordinate will take on a zero or one value.Clerc [6] gave out a brief outline of the PSO method for solving Traveling Salesman Problems (TSP).Hendtlass used the PSO algorithm to solve small-size TSP problems and improved its performance by adding a memory capacity to each particle [7].Pang et al. [8] used fuzzy matrixes to represent the position and velocity of the particles in the PSO and the operators in the original PSO formulas were redefined.Wang et al. [9] redefined the PSO operators by introducing the concepts of "Swap operator" and "Swap sequence".They solved the TSP problems in another way.We refer to all the DPSO algorithms above as actual DPSO algorithms, since they are from the original PSO.These studies have chosen some of the benchmark questions in the Traveling Salesman's Problems Library (TSPLIB), where the city size does not exceed 52.
Some researchers combined the DPSO with other algorithms to obtain better results or enlarge the problems scale.Zhong et al. introduced a new parameter named mutation factor to the DPSO for solving the TSP [10], and the particle is not a permutation of numbers but a set of edges.Yuan et al. [11] proposed a novel algorithm based on particle optimization algorithm (PSO) and Chaotic Optimization Algorithm (COA) to solve TSP problems.Shi et al. used an uncertain searching strategy and a crossover eliminated technique to design a novel particle swarm optimization algorithm for the traveling salesman problem [12].Machado et al. [13] presented a new hybrid model, based on particle swarm optimization, genetic algorithms, and fast local search, for solving traveling salesman problem.Fang et al. [14] combined particle swarm optimization with simulated annealing for TSP.Except traveling salesman problem, other problems were taken into account by many scholars.Anghinol and Paolucci [15] presented a new discrete particle swarm optimization approach for the single-machine total weighted tardiness scheduling problem.Jia et al. [16] improved particle swarm optimization (PSO) algorithm with rank-prioritybased representation employing a double justification skill for solving the resource-constrained project scheduling problem (RCPSP).Although most of the above hybrid algorithms have achieved good results, it is difficult to say that the discretization of PSO is successful because they all combine other algorithms.
For problems such as RCPSP and TSP, an excellent solution corresponds to a good arrangement consisting of good elements.In order to obtain these excellent elements, we propose a new DPSO algorithm PMPSO (probability-based particle swarm optimization algorithm).PMPSO is effective in jumping out of local optimum, computationally efficient, and easy to implement.We use this algorithm to solve the RCPSP problems, and the results proved the efficacy of the method.

Probability Mechanism Based Particle Swarm Optimization (PMPSO) Algorithm
In the PSO algorithm, a particle swarm is composed of a specific number of particles.At each iteration, all the particles move in the N-dimensional problem space with a certain velocity to find the global optima.The position of particle  ( = 1 . . .) is composed of D vectors denoted   , where    denotes vector at iteration .The velocity of particle  is denoted as of a D-dimensional vectors V  , where V   denotes the velocity vector at iteration .The velocity and position of each particle is adjusted according the following formulas: where    is the local best position that the i-th particle had reached at iteration k and    is the global best position that all the particles had reached at iteration k.  1 and  2 are random numbers between 0 and 1. w is the inertia weight. 1 and  2 are cognitive confidence coefficients, which are all constant numbers.

The Representation of Particles.
When solving continuous optimization problems, the viable solution is usually directly viewed as a particle   which can be updated according to (1) and (2).However, in discrete problems, the update operation of a viable solution is mostly meaningless unless we can represent a viable solution in a feasible way.Therefore, when the PSO algorithm is applied to the discrete domain, the representation of the particle becomes a key issue.In other words, the key point is to establish a link between particles and the solution of the problem.Most discrete problems can be regarded as sequencing problems.For example, the scheduling problem to a project comprised of five activities {1, 2, 3, 4, 5} can be considered as a sequencing problem to the five elements.Usually, the solution is in a form of permutation, e.g., (2,3,1,4,5).A permutation is formed by some fundamental elements, for example, the permutation of (1,2,3,4,5) is formed by the fundamental elements of (1, 2), (2,3), (3,4), and (4, 5).We call each fundamental element as subpermutation.A superior permutation must be consisted of some reasonably good subpermutations.To find these good subpermutations according to (1), we describe the permutation by a two-dimensional array called adjacency matrix, in which the elements correspond to the subpermutations are set to 1, the other elements are set to 0.Here are two examples of adjacency matrixes representing permutations: (1, 2, 3, 4, 5) can be represented by In the adjacency matrix, when the element  , =1, it states that the number "i" is in front of the number "j" and they are adjacent in the permutation.In our algorithm, the position of a particle is a (0, 1) adjacency matrix, which corresponds to a permutation combination.

Notation and Definitions of Operators. The operators in
(1) and ( 2) are redefined because the position and velocity are represented in the form of a matrix in discrete problems.
where    denotes the position of a particle at iteration k and    denotes the transition probability velocity vector at iteration k.    is the local best position that the i-th particle had reached at iteration k, and    is the global best position that all the particles had reached at iteration k.    ,    and    are all (0,1) adjacency matrixes corresponding to their permutations. 1 and  2 are matrixes with the same dimensions as   , whose elements are random numbers between 0 and 1. w is called inertia weight and its value is between 0 and 1.  1 and  2 are cognitive confidence coefficients, which are constant numbers.The symbol "⊖" denotes subtraction between matrixes and the negative elements in the subtraction result are set to 0. The symbol "⊕" denotes addition between matrixes and all elements greater than 1 in the addition result are set to 1.The symbol "b" denotes element-wise multiplication between matrixes.Suppose A and B are two matrixes with the same dimension, then the result of AbB is the same dimension matrix, where the elements are the products of corresponding elements in A and B. The symbol "⊗" is used to denote the modified multiplication.Let c be a real number, then c⊗A means all the elements of the matrix A are multiplied by c, and the elements with a value greater than 1 are set to 1.
The reason why one permutation combination is distinguished from another is its subpermutation.Therefore, a good permutation combination must contain subpermutations that are good and not available in other permutations.Based on the above ideas, we have designed a probability coefficient selection method.The basic idea of the method is that new permutations are generated based on probability selection, and those good subpermutations should have high probabilities of being selected.We try to find these reasonably good subpermutations through the PSO algorithm and make them have higher probabilities of being selected when generating new permutations.According to (5), the transition probability velocity vector  +1  is combined by three parts: .    is used to govern the extent to which the old velocity    determines the new velocity  +1  .Conventionally, w is assigned a value of 0.9.Moreover, c1 and c2 are the learning factors used to control the effects of the local experience and global experience on the new velocity, respectively.Conventionally, they are set to 1.    ⊖    is used to find the unique subpermutations in the local best position    but not in the individual    .   ⊖    is used to find the unique subpermutations in the global best position    but not in the individual    .In the following, we give an example to illustrate the operators in (2), ( 5), ( 6), (7), and (8).Let    be the matrix of representing the permutation (1, 2, 3, 4, 5), and    be matrix of ( ( representing the permutation (2, 3, 1, 4, 5).Then and this result denotes that the subpermutations (3, 1) and (1,4) are unique for    comparing with    .If the random matrix R 1 is ( ( ( 0.2 0.5 0.8 0.7 0.2 0.4 0.6 0.9 0.3 0.7 0.3 0.1 0.8 0.9 0.2 0.4 0.5 0.9 0.4 0.6 0.9 0.2 0.1 0.8 0.7 then the 0-1 matrix is transformed to ( ( ( 0 0 0 0.7 0 0 0 0 0 0 0.3 0 0 0 0 0 0 0 0 0 by The elements in the matrix corresponding to the subpermutations (3,1) and (1,4) are highlighted, which increases the probabilities that they are selected in the new permutation.
New particles are obtained by selection based on probability coefficients.According the main idea of PSO, we need to update the selection probability coefficients considering the current position and velocity.Therefore, we must convert    to the same form as probability coefficients matrix.In (6),    is defined as probability matrix corresponding to    , which stems from    according to   function.The function of   can be freely constructed according to the following rules: ensure the elements in    corresponding to the elements of 1 in    have big probabilities and the value of the elements in    that cannot be used to make a selection is 0.Here we give an example.Let    be ( ( ( 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 In (7), the selection probability coefficients matrix is updated considering the current position and the velocity.Suppose In (8),   ( +1  ) was defined as probability selecting operation, which selects the right elements and set them to one in  +1 according to the given probability coefficients in  +1  .The selection operation must obey the rule that if an element in  +1  has a larger probability coefficient value, the element with the same index in  +1  has a greater probability of being selected and set to 1. Then we get the new adjacency matrix and the new permutation is obtained by decoding it.

Selection Mechanism Based on Probability Coefficients.
A solution formed by means of probability coefficients selection can also break the local optimum and without losing the opportunity to find better solutions.We illustrate the selection process as follows.
Process 1 (process of probability coefficients selection).
Step 1. Set all the elements of adjacency matrix  +1  to zero.
Step 2. Select a row according to the row number of the row of the maximum element in the probability matrix of  +1  .
Step 3. The corresponding element in  +1  is used as the probability selection coefficient, and one element in the row is randomly selected accordingly.The selection rule is that elements with larger probability coefficients have a greater probability of being selected.
Step 4. Let the selected element be equal to 1, and set the corresponding so-called contradictory elements in  +1  to zero.There are three types of contradictory elements: (1) elements with the same row number as the selected element; (2) elements with the same column number as the selected element; (3) elements which can lead to incompletely permutation.
Step 5.If the stopping criterion is not satisfied, go to Step 2.
In the above example in Section 2.
We then use the same example to explain the implementation of Process 1. Firstly,  +1  is set to ( ( ( 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 According to Process 1, we select the row of (0.8007,0.0003, 0.0000,0.1386,0.0008)as coefficients, because the biggest value of the matrix is in this row.We randomly select a column according to the coefficients.According to the selection rule, each element could be selected except the third element, which coefficients equal to zero.However, the elements which have bigger probability coefficients have bigger probabilities to be selected.Hence, we assume that the first element was selected.Then in this example,  +1  is ( ( ( 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 ) .

Description of the Algorithm. PMPSO algorithm employs
the probability coefficients selection mechanism to improve the exploration of the solution space and it assures the convergence through the PSO evolutionary mechanism.We summarize the procedures of the proposed algorithm as follows.
Step 1. Initialization Step 1.1.Set the number of the particles and the max number of iterations.
Step 1.2.Each particle gets a random position matrix   corresponding to a permutation.
Step 1.3.Decode   and calculate the fitness of the position.
Step 1.4.Initialize the local best   =  0  and the global best   is the best of all the   .
Step 2. Calculate the position and velocity of each particle.
Step 3. Calculate the fitness of each position.If the fitness of the new position is better than that of the local best position of the particle, update the best position with the new position.If the fitness of the local best position of some particles is better than that of the global position, then update the global best with the local best position.
Step 4. If the stopping criterion is not satisfied, go to Step 2.
Step 5. Output the best solution and the fitness.

Problem Statement
Resource-constrained project scheduling problems (RCPSP) were investigated in many situations such as building construction, facility and equipment maintenance, and new product development.Due to their academic appeal and wide application, many researchers have always shown interest in these issues.The problem of RCPSP can be stated as follows.A project consists of a set  = {0, 1, . . ., J+1} of activities, where activities 0 and J+1 are dummy activities that represent the events of the start and finish time of the project respectively.The RCPSP is based on the following assumptions: (1) the durations of activities composing the project are deterministic and known; (2) all activities must obey their precedence constraints; (3) resources can be available in limited amounts and renewable from period to period;(4) activities cannot be interrupted when in progress; (5) the goal of optimization is to minimize project duration.The duration of activity  is denoted   , and resource type  requested by activity  is denoted r jk .For the dummy activities, for any resource type , we have d 0 = d J+1 = 0, and  0, =  +1, = 0.The availability of each unit of resource type  is   ( = 1 . . .).A project schedule can be represented by the start time of each activity denoted (s 0 , s 1 , . . ., s j+1 ) with s 0 = 0 and the finish time of each activity denoted by (f 0 , f 1 , . . ., f J+1 ) with f 0 = 0. Note that f J+1 = s J+1 and f j = s j + d j for all .  is the set of tasks being executed in the time period of .A schedule is feasible if it satisfies precedence relationships and all these constraints during each execution time period.The mathematical model of the problem can be formulated as follows: Equation ( 28) is the objective function that minimizes the completion time of the last activity of the project, that is, minimizes the span of the project.Equation ( 29) means each task's start time must be subject to the finish time of its preceding activities.It also can be expressed as activity  cannot start until all of the predecessor activities have finished, where   stands for the set of preceding activities of activity .Equation (30) ensure that the total number of used resources of type  cannot exceed its available maximum, denoted as   .
It was proved that the optimization of the RCPSP is NPhard in the strong sense [17].Mingozzi et al. [18] developed an exact method for solving the problems with sizes up to 30 activities.Donrndorf et al. [19] also reported satisfactory results using a branch-and-bound algorithm for the problems with 30 activities.For large-scale RCPSP instances, most researchers are employing metaheuristics algorithm to solve them, including tabu algorithms [20], genetic algorithms [21,22], ant colony optimization [23], simulated annealing [24], and variable neighborhood search [25].However, the heuristic methods based on heuristic rules are often variable effectiveness on different cases and may be trapped within local optima [21][22][23].Since the advent of the PSO algorithm, there have been many attempts to apply PSO to solve RCPSP.In contrast to other heuristic methods, PSO has the advantages of computational efficiency, great capability of escaping local optima, and easy implementation [26,27].Zhang et al. [28] took the priorities of the scheduling activities as particles to develop a DPSO algorithm for solving the RCPSP and the size of the problem investigated was 25.Xie et al. [29] combined the chaos algorithm with PSO to solve the typical multiple resources-constrained project scheduling problem.Ni et al. [30] proposed the hybrid particle swarm optimization based on differential evolution for solving the RCPSP with fifteen activities.In the following, we use the proposed PMPSO algorithm to solve the RCPSP problems.

RCPSP Solution Based on PMPSO Algorithm
To solve a RCPSP problem using PMPSO, we must have the information of the project, such as the duration of each activity, each resource required for each activity, the availability of each resource, and precedence relationships of the project activities.Then we can use PMPSO algorithm to solve the problem.The flow chart of the solution is given as Figure 1.

4.1.
Encoding and Decoding of the Particles.We directly take the set of activity number sequence as the source of particles, which has the advantage of simplicity and effectiveness.Then the permutation is transformed to a 0-1 adjacency matrix, which is the Encoding of particle.The 0-1 adjacency matrix is considered to be the position of the particle.According to the 0-1 adjacency matrix, we can derive the permutation, which called decoding of the particle.

Generation of Feasible Sequences.
A feasible sequence S = { 1 ,  2 , . . .,   } is a sequence in which the predecessor activity of any activity in the sequence must be arranged to the left of the activity.A feasible sequence is the basis of solving the RCPSP, regardless the algorithm.Researchers used various methods to get the feasible sequence.Montoya-Torres et al. [31] transformed sequences to feasible sequences using SSGS (Serial Schedule Generation Schemes) when they solved the RCPSP problem with genetic algorithm.Peng and Wang [32] adopted the priority list in the genetic representation combing the SGS to generate the feasible sequence.Xie et al. [33] allotted the priorities to the activities and then applied the bubble method to generate the feasible sequence.Zhou et al. [33] eliminated the improper sequence with a logic judgment procedure to get the feasible sequence.In all these articles, the feasible scheduling sequences of the activities were acquired indirectly.In order to improve the validity of the calculation, we have designed a method called Modulation Weighting Mechanism for Generating feasible scheduling sequences (MWMG), which can directly obtain a feasible solution.First, the precedence relationships of the activities must be given as a 0-1 matrix R, in which the element   =1 means that the activity  is a predecessor of activity .Secondly, an activity must be scheduled after all his predecessor activities.If the activity  as the predecessor of the activity  is scheduled, we can deactivate the coupling relationship between the activities  and  by setting the element   =0.Finally, when we generate the scheduling scheme according the probability weight matrix, the weight matrix should be modulated by the relation 0-1 matrix  to ensure that all the predecessor relations are obeyed.Because the activity "0" is dummy activity which represents the event of the project start, we begin to schedule from the activity "0".The generation process of the scheduling sequence is illustrated as follows.
Process 3 (process of generating feasible sequence).
Step 1.A 0-1 matrix representing the scheduling scheme is initialized by setting all elements to zero.Set the first column of the weight matrix to zero to let the project begin from activity "0".Set the first row of the relationship matrix to zero to release subsequent activities as activity "0" so that these activities can be scheduled.
Step 2. The value of the element in the modulation vector depends on the value of the above row vector.When the element value of the row vector is zero, the corresponding element value in the modulation vector is 1; otherwise, the element value is 0, where the element in the modulation vector equals to 1 means that the activity denoted by the column number can be scheduled, because all of its predecessor activities have been scheduled.
Step 3. A row is selected from the weight matrix as a weight vector and is modulated by multiplying each element in the weight vector by its corresponding element in the modulation vector.At the beginning, the selected row is the first row.
Step 4. Select the column number according the modulated weight vector randomly.Then, the corresponding element value in the 0-1 matrix of the scheduling scheme is set to 1 according to the row number and the column number.
Step 5. Renew the weight matrix according the rules in Section 2.3.
Step 6. Update the row number by the column number.
Step 7. Set the row of the relationship matrix to zero to release the subsequent activities of the activity represented as the column number so that they can be scheduled.
Step 8.If the stopping criterion is not satisfied, go to Step 3.
Step 9. Decode the 0-1 matrix of the scheduling scheme to get the feasible sequence.
Through the above generation mechanism, not only the feasible scheme can be directly obtained, but also the excellent subpermutations have better chance to enhance their probabilities.

Generation of a Scheduling Scheme.
We adopt the method of Serial Scheduling Generation Scheme (SSGS) to obtain a feasible scheduling scheme.According to the sequence, we schedule the activities successively.The following gives the generation process of the scheduling scheme.
Process 4 (process of generating a scheduling scheme).
Step 1. Select an activity  from the feasible sequence in order, which duration is d i , and calculate the early start time of the activity.
Step 2. Schedule the start time of the activity as early as possible, denoted t is , and calculate its finish time denoted t if , where t if = t if + d i .
Step 3. Calculate the availability of the resources required for the activity  during the time period in which the activity is performed.

4.
If the availability of the resources is satisfied, go to Step 5; otherwise, let t is = t is + 1 and go to Step 2.
Step 5.If there are still unscheduled activities, go to Step 1; otherwise go to Step 6.
Step 6.Output the result.

Total Process of Solving a RCPSP Problem Using PMPSO
Algorithm.According to Processes 1-4 we summarize the general process of solving a RCPSP problem using PMPSO algorithm as follows.
Process 5 (process of solving the RCPSP problem).
Step 1. Input the information of the project to be scheduled and initialize the parameters of the PMPSO algorithm.
Step 2. Initialize the probability coefficients matrixes according to the rules mentioned in Process 1.
Step 3. Generate 0-1 matrixes of the scheduling scheme using the method in Process 3.
Step 4. Decode each 0-1 matrix to get a scheduling scheme and calculate the fitness function value according to Process 4.
Step 5. Update the best position of the swarm and the best positions of each particle experienced.
Step 7. If the stopping criterion of the PMPSO algorithm is not satisfied, go to Step 3; otherwise go to Step 8.
Step 8. Output the result.

Computational Experiments
We conduct the experiments for PMPSO in this section.The test set comes from J30, J60, J90 and J120 of the PSPLIB, provided on the website (http://www.om-db.wi.tum.de/psplib/)[34].J30, J60 and J90 all contain 480 instances, while J120 contains 600 instances.Each instance set has four resource types.The currently best known solution in PSPLIB as reported on June 17, 2015.For the J30 set, these solutions are all optimal.The test environment is as follows: CPU: INTEL-2.4G,RAM: 4GB DDR 333, HD: 250G 7200 rpm, OS: Win 10.
We set  to 0.98, both c 1 and c 2 to 1, and the maximal number of iteration to 5000.Then we randomly select 10 instances from each test set and run the scheduler 100 times on each problem instance.Table 1 shows the statistical results of these solutions.
From Table 1, we observed that PMPSO can get the best known solution for all of the instances randomly selected from J30, and the maximum deviation of the mean of all solutions from the optimal value is only 0.17%.For other test sets, PMPSO can get most of the best known solutions, and the max deviation from the best known solution is 1.69%.As we can see from Table 1, the best known solutions of these instances were found by various researchers using various methods.In this paper, we use PMPSO to get some of the best known solutions and some solutions deviating from the best known solutions slightly.We also execute the scheduling program for each instance of each problem set.Table 2 shows the results of relative deviations from the currently best known solution.
Tables 1 and 2 summarize the experimental results of PMPSO.We observed that the larger the problem, the greater the deviation from the best value or the current best known solution, and as the size of the problem increases, the average number of best solutions obtained decreases.For most instances, we can get the best known solutions.Moreover, as it is shown in the Table 2, using the algorithm we can get a very high percentage of hits for the problem size within 30.We note that the best known solutions in PSPLIB are found by different people using kinds of algorithms.Therefore, PMPSO is an algorithm that adapts to problems of various scales and has unique advantages especially for smallscale problems.
To observe the performance of PMPSO, we also illustrate the evolution curves of the instances of J305-2 and J1201-5, which are showed in Figures 2 and 3, respectively.
Figure 2 shows that evolution curve drops sharply and then converges at the best value.This shows that the PMPSO algorithm has the advantage of fast convergence, which is the most advantage of the original PSO algorithm.It also demonstrates the outstanding performance of solving small-scale problems.
For large-scale problems, PMPSO is doing well too, which can be observed from Figure 3.For the instance of J1201-5, which has 120 activities, PMPSO can get its best known solution at 3200 iterations.
In Figure 3, before the curve levels off to the minimal value, it experiences several platforms.This shows that the PMPSO algorithm can effectively avoid converging on local optimum and reduce solution time greatly.

Conclusions
This paper proposes a probability mechanism based particle swarm optimization algorithm.The design of the algorithm includes the following features: a special decoding scheme, a mechanism for finding good subpermutations, a new particle swarm iteration formula, and a selection mechanism based on probability coefficients.When applying the PMPSO algorithm to discrete optimization problems, we only need to convert the solution to a 0-1 matrix form, and then use our mutated PSO formula.The experimental results showed that state-of-the-art solutions can be obtained using the PMPSO algorithm to solve some benchmark problems.In our future work, we will establish a better probability selection mechanism to improve the selection process.We will implement dynamic adjustment of parameters to improve the performance of PMPSO.In addition, we can combine PMPSO with other algorithms, such as local search, metaheuristic, to get better results.In the experiments of this paper, PMPSO is capable of solving the RCPSP efficiently.Although these experiments are made based on PSPLIB, PMPSO should perform equally well for the other types of sequencing problems, such as production scheduling, route planning, and the assignment problem.

Figure 1 :
Figure 1: Flow chart of the solution.
Secondly, we divide all the elements by 1000, except for the elements corresponding to the elements of one in    .Finally, we set the diagonal elements to be 0, because the diagonal elements in    cannot equal one and the corresponding elements in    cannot be used for the selection operation.