Sugar cane distribution optimization problem is a critical problem in Tucum´an, Argentina, with geographical, environmental and economic characteristics that makes it different from similar transportation or supply chain problems. Motivated by this problem, we introduce in this paper a new heuristic algorithm to solve it. The problem is represented as a non-linear variant of the Generalized Assignment Problem, and the algorithm introduced is a Binary Particle Swarm Optimization implemented on a MultiAgent System, written in SARL language. The new algorithm was tested on a previously developed benchmark suite, and it can efficiently solve the problem, finding the optimal solution for all test instances of different types.
In a previous work presented by the authors [Majorel Padilla et al., 2015], a non-linear variant of the well known Generalized Assignment Problem (GAP) was introduced, as well as a benchmark suite with its optimal results, for using as a reference for future algorithms. Moreover, a model of the problem was presented, using the organizational approach from MultiAgent Systems (MAS).
This variant was designed to solve the sugar cane distribution optimization problem in Tucum´an, a province in the northwest of Argentina. Being one of the largest industries in the province, the distribution of sugar canes to sugar mills lacks an optimized scheduling of the complete sugar cane supply chain. This fact not only reduces profits to both sugar cane farmers and sugar mill owners, but also has an important impact on highway traffic, road safety and environmental pollution.
There have been several efforts to optimize similar sugar cane supply chain optimization problems throughout the world [Giles et al., 2009, Le Gal et al., 2009, Yosnual & Supsomboon, 2004]. However, because of the substantial differences between the supply chain processes among countries, the solutions proposed in those cases can not be applied directly.
Considering this background, this paper proposes a new heuristic algorithm for the aforementioned problem. The resulting algorithm is a Particle Swarm Optimization (PSO), in its Binary version, implemented on a MultiAgent System. The benchmark suite developed for the problem is used to test the new algorithm, and the results of these tests are exposed.
The new algorithm was written using SARL [Rodriguez et al., 2014]. SARL is the first general purpose language for agent oriented programming, based on the Java language, fully independent of any other platform or any other agent architecture. Through its Capacity-Role-Interaction-Organization metamodel, SARL provides a well defined set of basic concepts (agent oriented first-class abstraction) needed to implement MAS. It is an effort to improve user experience in developing complex systems, thanks to its simple and extensible syntax.
The rest of the paper is structured as follows: in Section 2 a brief recapitulation of the GAP is presented, as well as the specific variant which constitutes the problem to solve. This section also describes some related works between MAS and PSO. Next, Section 3 explains the basics of PSO and details about the MAS implementation are introduced. Finally, in Section 4 the testing methodology is described and the results presented, with final conclusions and future work in Section 5. 2
In this section we will give a brief overview of the GAP, and the variant introduced to solve the proposed problem, as well as related works of MAS and PSO which were used as basis for the proposed algorithm. 2.1
The Generalized Assignment Problem is defined in [Chu & Beasley, 1997] as a combinatorial optimization problem, which consists in finding the assignment of tasks into agents that minimizes cost, with the constraint that each task is assigned to exactly one agent, and subject to agent’s capacity. For the rest of this paper, we follow the notation defined for GAP in [Ozbakir et al., 2010]: n represents the number of agents, m is the number of tasks; when task i is assigned to agent j, this assignation is expressed as xij , and it has an associated cost cij or profit pij . Besides, each agent has a maximum capacity bj , and task i consumes aij of agent j resources.
GAP is an NP-hard problem [Fisher et al., 1986], and finding if a valid solution exists is NPcomplete [Martello & Toth, 1990]. Many deterministic and heuristic algorithms can be found in literature [Cattrysse & Van Wassenhove, 1992, Osman, 1995, Ramalhinho & Serra, 2008, Yagiura et al., 1998], and also many variants [Krumke & Thielen, 2013, Laguna et al., 1995, C. Rainwater & Romeijn, 2009].
Nevertheless, none of these variants matches all the features required by the particular problem we were trying to solve, so we proposed a new variant that models sugar mills as GAP agents and sugar cane farms as GAP tasks, and called this variant Variable Profit GAP (VPGAP). In VPGAP, the assignment xij represents when a sugar cane farm i is assigned to sugar mill j, it has an associated profit pij , and this sugar cane farm provides aij tons of sugar cane to the mill for processing. The sugar mill’s maximum capacity is modeled as bj , and it also has the following distinct features:
The agents (sugar mills) have a non-linear efficiency, hj , related to the amount of resources (sugar cane) consumed by the tasks (sugar cane farms) associated to it.
Agent’s peak efficiency is about 85% of its maximum capacity. Its efficiency is constantly increasing until this value, and constantly decreasing after it.
n m maximize Z(X) = ∑ ∑ pij xij hj i=1 j=1 ( n ∑ aij xij i=1
) dj n ∑ aij xij bj i=1 m ∑ xij = 1 i = 1; : : : ; n j=1
j = 1; : : : ; m xij = { 1; if task i is assigned to agent j 0; in other case. (1) (2) (3) (4) where Z(X) is the system’s global efficiency, and hj represents agent’s own efficiency, as a percentage of its maximum capacity. The function used in this work is a Gaussian with a mean of 0:85 and a standard deviation of 0:25, in the interval [0.30, 1.00] and 0 outside this interval. This simple function accomplishes with the constraints from above. 2.2
Over the past three decades, MultiAgent Systems emerged as a new computational paradigm, and have been used in a wide range of complex problems, from robotics to distributed problem solving [Wooldridge, 2009]. The paradigm’s main concept is the agent, which is defined as “a physical or virtual entity with a high degree of autonomy, independent, capable of cooperating, competing, communicating, acting and taking control of its own behavior within its own goals” [Weiss, 2013]. These entities operate within an environment, which they are able to perceive and to modify. Besides, agents can group themselves into societies, in which they can cooperate with each other to reach common or individual goals which could not be achieved individually. Therefore, it is commonly expressed that the solution emerges as result of the interaction among agents.
However, only a few examples can be found in the literature of MAS in conjunction with PSO. One of these examples is [Zhao et al., 2005], whose MultiAgent based Particle Swarm Optimization (MAPSO) was used to find optimal solution to the distribution of reactive energy in power systems, with excellent results. A different example is found in [Pugh & Martinoli, 2007], where a multi-robot search algorithm is used to find optima in a multi- dimensional function space, inspired by PSO, and with several advantages over similar methods. In [Wang et al., 2010] a MAS is used for remote house control, while PSO is used to find the optimal set of control points of the system. Finally, a new algorithm named HMAPSO (Hybrid MultiAgent-based Particle Swarm Optimization) was introduced in [Kumar et al., 2011] to find a solution for the economic problem of energy distribution. 3
In this section, the proposed algorithm is presented. First we describe the general PSO algorithm, its components and behaviors, and then Binary PSO algorithm is described. Finally, the mixing between Binary PSO and MAS is introduced. 3.1
Particle Swarm Optimization algorithm [Kennedy & Eberhart, 1995] is a bio- inspired algorithm, based on a meta-heuristic inspired in social behavior of certain species, such as flock of birds’ flight, or fish bank movement, and it has grown in popularity thanks to its quick convergence and its simple implementation. In the literature, PSO is used in traditional optimization [Hu et al., 2004], in multi-objective optimization [Zhang et al., 2003], and in dynamic optimization [Blackwell, 2007], among others.
PSO is a population based algorithm, but it is not an evolutionary algorithm, because it lacks a selection mechanism. In PSO, every particle (or agent) represents a possible solution to the problem, it is part of a swarm of particles with similar characteristics, and it has to make decisions regarding its behavior. These decisions are made taking into account two major components: an individual component of each particle, representing results obtained in past explorations, and a social component, representing the influence of the rest of the particles of the swarm. With these two components, the particle can determine where it will move to reach a new solution of the problem solutions space.
It is important to note that in PSO the population is fixed-sized, the solution’s space is predefined (the particles are not allowed to move outside this boundary), and the particles are constantly searching for new solutions. Moreover, PSO makes little assumptions about the problem to be solved, and is capable of searching in very large solution’s spaces. As a disadvantage, PSO can not guarantee that the optimal solution to the problem can be found.
In its more general form, each PSO particle is composed by the following elements: a vector Xi representing the actual position of the particle; a vector pBesti holding the particle’s best solution found at the moment; a velocity vector Vi holding the particle’s future direction; a f itness xi value, indicating how good the actual solution is; and a f itness pBesti value, indicating how good the particle’s best solution found at the moment is.
The population is initialized with random positions and velocities for each particle of the swarm. Then, in an iterative process, particles start to move in the search space until a new position is found. In this new position, the value f itness xi is calculated, and if this value is better than the best fitness found until then, pBesti and f itness pBesti are updated.
In each iteration, the particle’s velocity is updated following equation 5, which is comprised of three terms: the first one represents inertia, the particle’s tendency to keep moving in the same direction [Shi & Eberhart, 1998]; second one represents history, the particle’s tendency to stay closer to positions that produced good results previously; and the third one represents the influence of the rest of the swarm, or the imitation of the best particle of the swarm (variable gBesti is the value of the best solution found at the moment by the rest of the swarm particles).
Vi(t + 1) = w
Vi(t) + c1 rand1 (pBesti + c2 rand2 (gBesti
Xi(t)) Xi(t))
As the vectors change in every iteration, the variable t is used to indicate time. The value w is named inertial factor. c1 and c2 are constants known as speedup constants and they usually have the same value. rand1 and rand2 are random numbers between 0 and 1.
When every particle has updated its velocity, it sums up with the actual position in order to determine the new position for the next iteration, as stated in equation 6.
Xi(t + 1) = Xi(t) + Vi(t + 1)
Many articles can be found in the literature related to the correct choice of w, c1 and c2 values, which are critical to the algorithm’s convergence [Clerc & Kennedy, 2002, Van Den Bergh, 2002]. In [Shi & Eberhart, 1998] the authors showed that when w is small, the algorithm makes a more localized search and it is strongly dependent of the initialization. While w increases, the search widens, but it is more difficult to find the optimal value. So, it is recommended that the value of w decreases over time. But later, in [Shi & Eberhart, 1999] and in [Zhang et al., 2003], authors concluded that whether the value of w should increase or decrease over time depends on the problem to be solved. 3.2
In [Kennedy & Eberhart, 1997], Kennedy proposed a binary version of the PSO algorithm, with the goal of applying PSO to combinatorial optimization. This modification affects each particle’s codification, and the position updating equation. In Binary PSO, the velocity is used as an input to a sigmoid function indicating the probability that the position takes the value 1. Then, a random number is generated and the new position is determined by equation 7.
Xi(t + 1) = { 1; if rand() < sig(Vi(t + 1)) 0; other:
Result: best solution found Pop = CreatePopulation(N); while termination condition not reached do for i = 1 to size(Pop) do
Evaluate particle Xi of Pop; if tness(Xi) is better than tness(pBesti) then pBesti = Xi; fitness(pBesti) = fitness(xi); end end for i = 1 to size(Pop) do
Choose gBesti; Calculate new Vi; for d = 1 to n do if rand() < sig(Vid) then
Xid = 1 end else end
Xid = 0 end end end
Algorithm 1: Binary PSO algorithm pseudo-code. 3.3
With the approach explained in previous sections, there exists a one-to-one relationship between a PSO particle and an MAS agent. At the same time, each one of these agents is identified as a possible solution for the optimization problem. In our proposed MAS model, agents should have one of two different roles: Swarm or Particle. In Figure 1 a simple scheme shows the interaction between the two agents and their methods and required capacities.
The Swarm agent has a unique instance in simulation time, and is responsible for the creation of all of the instances of Particle agents defined in the population. Besides, Swarm agent receives from every Particle agent the actual position values, Xi.
Each Particle agent receives from the Swarm agent its corresponding initialization values, and with these values it calculates its new velocity Vi, its new position Xi and its new f itness xi value. This f itness xi value is also passed to the Swarm agent so it can determine which of the Particle agents has reached the best solution so far.
Thanks to the interoperability of SARL with Java (the possibility to create Java objects from SARL, and invoke its methods from SARL-defined agents in a completely transparent way), the code migration from [Majorel Padilla et al., 2015], written for Janus platform, was simple and straightforward. Additionally, SARL is fully distributed and allows for parallel execution of agents’ behavior and automatic detection of execution kernels, features that can be used for distributed computation of simulation.
The code implemented is open-source and can be found in https://bitbucket.org/gitia/zafra-pso-sarl. 4
In Section 2.1 the problem to be solved was presented, as a variant to GAP. A benchmarking suite is provided for this problem, and is available through http://www.gitia.org/projects/vpgap/. The benchmarking instances of this suite are identical to other standard suites for this kind of problems [Chu & Beasley, 1997, Yagiura et al., 1998], but with a smaller size (because the problem is much more difficult, as stated before in Section 2.1). The benchmarks were made with instances of types A, B, C and D, being D-type ones the hardest to solve, because cij and aij values are inversely proportional. In summary, 20 different instances were used, of low and medium complexity.
Benchmarks were executed on a single computer with a Intel Core i7-5500U @2.4 GHz microprocessor, 6 GB of RAM and a 64-bit Operating System. Ten independent executions of the PSO algorithm were made for each one of the instances, and the results are summarized in Table 1.
In Table 1, columns labeled Type, n and m indicate type of instance used, the number of agents and the number of tasks used, respectively. Column labeled Opt shows the optimal value as a reference, obtained previously for that particular instance in [Majorel Padilla et al., 2015]. The next three columns provide the number of particles used in the swarm, the number of iterations used by the simulations, and the PSO parameters used. In column labeled Parameters, three values are separated by semicolons: the range for inertial factor w, and the values for constants c1 and c2.
The best result obtained by the PSO algorithm for every instance is under the column labeled Found. It is remarkable that for all instances this value is exactly the same as the optimal value used as reference. The mean time required by the ten independent executions of the PSO algorithm to find the optimal value is also presented. Although these execution times may seem large enough for the simulation, the main objective of the tests was to prove that the algorithm can find the optimal solutions, without any timing restrictions. For example, the stopping criteria used was the number of iterations, for letting the particles to continue the searching of the solution’s space even though they have reached the optimal value.
Finally, the last column of the table, labeled Con dence, shows how many times the PSO algorithm could find the optimal value. This values should ideally be 100%, but as PSO algorithm can not guarantee its convergence, there is a probability that particles end up stuck at a local maximum. Nevertheless, for almost every instance this value is 80% or higher, with the exception of a D-Type instance with 4 agents and 12 tasks, for whom the optimal value could be found only in half of the tests. 5
This paper introduced a new heuristic algorithm, a Particle Swarm Optimization implemented on a MultiAgent System, in order to solve a non-linear variant of the Generalized Assignment Problem. This variant was developed with the goal of optimizing sugar cane distribution in Tucum´an, Argentina.
The algorithm was tested with 200 independent tests, part of a benchmark suite previously developed specially for this problem, and it could reach the optimal solution 91% of the times. We are comfortable to conclude that the algorithm can solve the problem, finding the optimal solution for all test instances (of different types). Moreover, results were significantly better than previous ones reported in [Majorel Padilla et al., 2015]. The execution times were larger than expected, but this is partially because we did not fully use the inherent parallel capabilities of MAS yet, and this is left as future work.
We wish to expand the results presented here in two different ways: first using larger test instances, for modeling real situations, and using the well known MAS scalability to run the simulations in a parallel environment.
This research was supported by Universidad Tecnol´ogica Nacional under project UTI-UTN 1318.