Monte Carlo simulation based optimization and analysis of
3.2 Determining optimal safety stock level in multi-echelon supply chains
3.2.6 Further improvement of PSO algorithm by memory- memory-based gradient search
As we described earlier, the classical PSO algorithms provide unsatisfactory results, especially in the optimization of problems with nonlinear constraints.
In the previous section, we provided a novel gradient-based method to over-come this diculty. However, classical gradient calculation cannot be applied to stochastic and uncertain systems. In these situations stochastic techniques like Monte-Carlo (MC) simulation can be applied to determine the gradient.
These techniques require additional function evaluations. We have devel-oped a more economic, memory based algorithm where instead of generating and evaluating new simulated samples the stored and shared former func-tion evaluafunc-tions of the particles are sampled to estimate the gradients by local weighted least squares regression. The performance of the resulted fully informed, regional gradient based PSO is veried by several benchmark problems.
As we have seen previously, the equations for calculating velocities were modied. We will discuss these modications in more details, and provide further improvements to propose a memory-based approach. So the modied equation for velocities is the following:
vj(k+ 1) = w·vj(k) +c1·rand()·(xpBest,j−xj(k)) + (3.16) c2·rand()·(xgBest−xj(k)) +
c3 ·gj(f(x(k))) wheregj(f(x(k))) = ∂f∂x(x(k))
j(k) represents the partial derivatives of the objective function, and is the weight for the gradient term. It should be noted that this concept can be interpreted as inserting a gradient-descent update step into the iterations of classical PSO,x(k+ 1) =x(k)−η5f(x(k))where the learning rate is equal to η=c3dt.
The above algorithm can be applied only to continuously dierentiable objective functions5f(x(k)). The simples approach to calculate the gradient is HGPSO [119] is the numerical approximation of the gradient.
∂
∂xi(x(k)) = f(x(k) +Ei)−f(x(k))
(3.17)
The main drawback of this approach is that the step size is dicult to design and the whole approach is selective to noise and uncertainties. It is interesting to note that PSO itself can also be interpreted as a gradient based search algorithm where point dierences are used as approximation of the regional gradient. The normalized gradient evaluated as the point dierence method is
e = f(xj)−f(xi)
k(xj−xi)k (3.18)
This point?dierence estimate can be considered as regional gradient for the local region of xi and xj. Hence, the velocity of PSO can be interpreted as a weighted combination of a point-dierence estimated global regional gra-dient(xgBest−xj(k))and a point-dierence estimated ner regional gradient (xpBest−xj(k)).
Our aim is to further improve the optimization by providing robust yet accurate estimation of gradients. To obtain a robust estimate a so called regional gradient should be calculated. When the function is dierentiable the gradient for a region S it is calculated as
5f(x)∗ = 1 volume(S)
Z
x∈S
5f(x)dx (3.19)
where S represents the local region where the gradient is calculated. How-ever, when heuristic optimization algorithm should be applied the objective
function is mostly not continuously dierentiable or not explicitly given due to limited knowledge. Therefore the gradient is calculated as follows:
5f(x)∗ = R
x∈Sf(x)dx R
x∈Sdx (3.20)
An interesting example for this approach is how regional gradient is calcu-lated in the Evolutionary-Gradient-Search (EGS) procedure proposed by R.
Solomon [143]. In EGS at each iteration λ test candidates are generated by applying random "mutations" of x(k).
vi =x(k) +zi (3.21) where zi is a Gaussian distributed variable with zero mean and standard deviationσ/√
n. For n1 these test points will be distributed on a hyper-sphere with radius σ. By using information given by all candidates the procedure calculates the gradient and a unit vectore(k)that points into the direction of the estimated gradient:
These techniques require additional function evaluations. It is important to note that this concept discards all information related to the evaluation of vi.
We have developed a more economic, memory based algorithm where instead of generating and evaluating new simulated samples the stored and shared former function evaluations of the particles are sampled to estimate the gradients by local weighted least squares regression. This idea is partly similar to the concept of the fully-informed particle swarm (FIPS) algorithm proposed by Mendes et al. [110]. FIPS can be also considered as a hybrid method for estimating the gradient by a point dierence of the weighted regional gradient estimate(Pj−x(k))based on thelBestsolutions and adding an additional gradient related term to the velocity adaptation:
νj(k+ 1) =· · ·+c3(νj(k−1) +ϕ(Pj(k)−x(k))) (3.24)
FIPS utilizes only the current xlBest,i(k) values so it does not have a memory. The main concept of our work is the eective utilization of the previous function evaluations. So instead of generating new and new samples and loosing information from previous generations the whole trajectories of the particles are utilized.
The weighted regression problem that gives a robust estimates of the gra-dients is formulated by arranging these former function evaluations{x(k), f(x(k))}
into indexed data pairs{vi, f(vi)} and calculating the following dierences:
5fi(k) = f({vi)−f({x(k)) (3.26)
where λrepresents the number of particles. The weighted least squares esti-mate is calculated as
gj(k) = 4XT(k)Wj(k)4X(k)−1
4XT(k)WjT(k)4f(k)) (3.29) where the Wj(K) weighting matrix is a diagonal matrix representing the region of the jth particle.
Similarly to EGS a Gaussian distributed weighting is used:
βj,i(k) = 1
(2π)n/2exp 12(vi−x(k))T P−1
(vi−x(k)) (3.30) where vi is the ith row of the matrix 4X, j represents the jth particle, wj,i(k) = Pβj,i(k)
i=1βj,i(k) is the normalized probability of the sample, P
= Iσ is a diagonal matrix where parameter σ represents the size of the region used to calculate the gradients. By using the information given by all previous states of the particles it is possible to calculate a unit vector that points into the direction of the estimated (global) gradient. The resulted algorithm is given in Algorithm 2.
Algorithm 2 Pseudo code of the improved PSO algorithm
1: Initialize particles
2: while not terminate do
3: for all particle do
4: Calculate tness value
5: if tness < pBest then
6: pBest= tness
7: end if
8: end for
9: Choose the best particle as the gBest
10: for all particle do
11: Calculate normalized distance base weights of previous function eval-uations particles by equation 3.30
12: Calculate regional gradients by equation 3.29
13: Calculate particle velocity by equation 3.16
14: Update particle position
15: Store particle position and related cost function in a database, {vi, f(vi)}
16: end for
17: end while Results
We tested the novel algorithm using several functions, including deterministic and stochastic ones as well. Fig. 3.11 presents four of them and Table 3.2 contains the mathematical representation of the analyzed functions.
During our tests, the found global best value from the population, the mean of the best values of each individual, the iteration number what al-gorithm performed before termination and the iteration number when the global best was found were registered. In all tests, we applied Monte Carlo method to evaluate the performance of our approach eectively, thus the weight of the gradient part in the objective value calculation for the indi-viduals was adjusted by 0.1 in the domain of [0,1.25], and for each gradient weight, 500 MC simulations were performed. The result of our tests for the gBest values is presented in Fig. 3.12. It can be deducted from the gure that the best global optima was found using 0.7 as the weight of the gra-dient part (since it is a minimization task), while the standard deviation is slightly smaller than in the classic PSO which is presented by the rst data point in Fig. 3.12, where the weight is equal to 0 (w-grad = 0). To prove the eectiveness of our method, we present the details of our tests in Table 3.3.
Figure 3.11. Surface of the tness function called "dropwave" (a) the
"griewangks" function (b), the stochastic function we used (c), and a stochas-tic version of "griewangks" (noise added) (d).
Including two dierent deterministic functions (dropwave, griewanks) and two stochastic one (griewankgs with noise and another stochastic function), highlighting the best results in each cases.
Using our method, PSO nds better solution, i.e. the objective value of the best individual and the mean of all objective values in the population is decreased while the number of iterations until the nal solution found is decreased also. It yields stronger convergence during the iterations of the algorithm, thus the novel method increases the eciency of PSO. Obviously
Function Equation
a stochastic function f(x, y) = sin(120/x) + cos(60/y) Table 3.2. Mathematical equations of the analyzed functions.
the proper setup for the parameters of the PSO and the weight of the gradient is highly problem-dependent, however, during our tests, we found 0.7 as a generally applicable weight for the gradient, and 10 percent of the domain for σ, which setting in most cases improves the eciency of PSO. We propose a simple ne-tuning technique in the following to setup the parameters of the algorithm.
1. Set all parameters to zero, i.e. c0 =c1 =c2 =c3.
2. Tunec3 according to the learning method of classic gradient methods, i.e. increase c3 gradually. If oscillation occurs, divide it by 10. Find a stable setting.
3. Set the momentum, i.e. c0, which is typically 0.1 or 0.2 in the literature.
Increase it gradually, until some improvement is achieved. Find a stable setting.
4. Tunec1, i.e. increase it gradually until some improvement is achieved.
5. Finally, setc2 = 1.25−c3
These techniques propose a reliable method for tuning the parameters, however, our tests showed clearly that c3 = 0.7 is a generally good choice, and with c1 = 0.5 and c0 = 0.6the algorithm operates stable.
Figure 3.12. Histograms for the "gBest" values using the function called
"griewangks" with noise. In the title of the subgures, mean represents the mean value of the histogram, std is the standard deviation and w-grad is the weight of the gradient part in the objective value calculation of the individuals.