Enhanced memetic search for reducing energy consumption in fuzzy flexible job shops

2.1Fuzzy processing times

We model processing times as triangular fuzzy numbers (TFNs), a particular type of fuzzy numbers [67], with an interval [a1,a3] of possible values and a modal value a2. A TFN can be represented as a triplet a=(a1,a2,a3) and its membership function is given by the following expression:

For our problem, we need two arithmetic operations between TFNs, the sum and the maximum, as well as the scalar multiplication with a real number d. All these operations can be obtained using the extension principle but, unfortunately, the set of TFNs is not closed under the maximum and for this reason we rely on an approximation thereof. This way, the sum of two TFNs a and b is given by:

the scalar multiplication is given by:

and the maximum is approximated using interpolation on its three defining points as:

When the result of the actual maximum operation is a TFN, then it coincides with the approximated value from Eq. (4). When this is not the case, the approximation maintains the support and modal value of the non-approximated maximum. Maintaining the support is in fact a very important property in scheduling since a good model of uncertainty should contemplate all possible scenarios. If the support were not maintained under the maximum, the resulting fuzzy schedule might not cover starting or completion times for operations which are actually attainable in a real scenario, resulting in poor and incomplete solutions. The interested reader is referred to [57] and references therein for further arguments supporting the use of this approximation, widely used in fuzzy scheduling.

Another important characteristic to be taken into account when working with TFNs (or fuzzy numbers in general) is that the subtraction is not the inverse of the addition. Indeed, following the Extension Principle, the subtraction of two TFNs a and b is defined as:

Clearly (a+b)-b≠a. This is so because the definition assumes that a and b are non-interactive variables, that is, they can be assigned values within their support independently of each other [67]. However, when both TFNs a and b are linked, applying the subtraction results in an over-imprecise quantity.

In scheduling, it may be the case that we have two TFNs a and b which are indeed linked in the sense that a is the result of adding b to a third TFN c, that is, there exists c such that a=b+c, and we may be interested in recovering c from a and b. In this case, we can compute the inverse of the sum as follows:

so c=a⊖b. One could say that we want to “remove” b from a to obtain c, so we shall refer to the operation ⊖ as operation remove. The use of such an inverse of the fuzzy addition in scheduling dates back to [68].

Finally, since there exists no natural relation of total order for TFNs, we need to resort to some ranking method. Here, we will use a ranking based on the expected value of TFNs, so

where 𝔼⁢[a]=a1+2⁢a2+a34 is the expression of the expected value of a fuzzy number in the particular case of a TFN [69]. Besides its wide use in the fuzzy scheduling literature, the choice of this ranking method is supported by some empirical analysis that suggests that using it contributes to the robustness of the obtained solutions [57].

2.2Fuzzy schedules

A schedule is a solution to the problem, a pair (𝝉,𝐬) consisting of both a resource assignment 𝝉 and starting time assignment 𝐬 for all operations. A solution is said to be feasible if all constraints hold. For an operation o∈𝒪, let τo be the resource assigned to o in this solution and let so and co=so+po⁢τo be its starting and completion times respectively. Then, for an arbitrary pair of operations u,v∈𝒪 precedence constraints within a job hold if and only if ∀i⁢cui⩽svi when χu=χv,ηu<ηv and capacity constraints hold if and only if ∀i⁢cui⩽svi∨∀i⁢cvi⩽sui when τu=τv for i∈{1,2,3}.

Notice that the starting time assignment 𝐬 induces a global operation processing order 𝝈 and a resource operation processing order 𝜹r for every r∈ℛ. The position of operation o in 𝝈 is denoted by σo and the position in which operation o is executed in resource τo is denoted by δo.

For a given a feasible solution (𝝉,𝐬), the fuzzy makespan is defined as:

where Cj denotes the completion time of job j, that is, Cj=co for the operation o such that χo=j,ηo=nj.

2.3Total energy objective function

In the fuzzy flexible job shop problem we consider two types of energy: active energy and passive energy. Passive energy 𝑃𝐸r is intrinsic to each resource r∈ℛ and is consumed whenever the resource is on. It is the product of the passive power consumption of the resource 𝑃𝑃r and the time r is on. We consider that all resources are turned on at the same time (at instant 0) and turned off when all operations are completed, i.e., resources are on for the entire makespan. Thus,

Active energy consumption 𝐴𝐸𝑜𝑟 occurs when an operation o is processed in a resource r. It depends on the power 𝐴𝑃𝑜𝑟 required to execute operation o in resource r, so:

We adopt an additive model so the total energy consumption Er of a resource r∈ℛ is the result of adding the passive and the active energy as follows:

Given these definitions, the total energy consumption is obtained as the sum of the total energy consumption of each resource:

This energy model corresponds to the first high-level approach to energy consumption as explained in Section 1. However, it is formulated in a somewhat different way to how it is done in the deterministic setting. In deterministic scheduling, the possible states of the resources are usually disjoint, so they can contribute to either the processing energy cost (the energy consumed when they are executing) or the non-processing energy cost (the energy consumed when they are idle). To compute the latter, it is necessary to calculate the idle time periods for each resource, i.e., the time span between the end of an operation and the start of the following one in the same resource. When durations are real numbers, this is easily done by subtracting the completion time of the first operation from the starting time of the second. However, a direct translation of this approach to the fuzzy job shop means subtracting non-interactive quantities, thus introducing some undesired and artificial uncertainty, as explained above in Section 2.1. For this reason, inspired by energy models used to reduce energy consumption in data centers, we propose instead to make use of overlapping states, resulting in an additive model where the two involved energies are aggregated. Notice that in the deterministic case both approaches, the standard one using disjoint states and the one adopted here using overlapping states are in fact equivalent. In the fuzzy setting, this alternative model has the advantage of avoiding added artificial uncertainty. This is so because it decomposes the total energy calculation in a way that does not imply subtracting interactive fuzzy numbers. The use of overlapping resource states with the goal of avoiding computing the length of idle periods under uncertainty is also underlying the definition of fuzzy total energy proposed in [65], although here the authors assume that resources are turned off individually.

2.4MILP model

In this section we present a MILP model based on [70]. As mentioned in Section 2.2, a solution to the problem is given by a pair (𝝉,𝐬) consisting of both a resource assignment 𝝉 and a starting time assignment 𝐬 for all operations. This suggests defining two types of binary decision variables. For every operation o∈𝒪 and every resource r∈ℛo, let

and for every pair of operations u,v∈𝒪, let

It is trivial that the variables X𝑜𝑟 give us enough information to build 𝝉, and Y𝑢𝑣 to build a global operation order 𝝈. However, for the solution to be complete we need the actual starting times, represented by another set of variables soi, i∈{1,2,3}, o∈𝒪, corresponding to the i-th component of the starting time of operation so, a TFNs. The resulting MILP model is:

subject to

Equation (16) ensures that operations are assigned to exactly one resource and Eq. (17) ensures that operations have at most one immediate predecessor in the resource. These two equations conform the two basic constraints of the problem. In order to define the starting time of operations, Eq. (18) ensures that the starting time of an operation is greater or equal than the addition of the starting and processing times of its predecessor in the job, so precedence constraints within jobs hold. Analogously, Eqs (19) and (20) prevents overlapping of operations in the same resource. Finally, Eqs (22)–(26) model the objective function, that is, the total energy consumption.

3.1Evolutionary algorithm

The evolutionary algorithm maintains a population of solutions so, at each generation, the individuals in the population are combined to obtain a new population that will replace the old one. A solution (𝝉,𝐬) is represented by the tuple (𝝉,𝝈), i.e., the resource assignments and the global processing order induced by 𝐬. To decode an individual, each operation is assigned the earliest starting time such that the order defined by 𝝈 is not altered. Individuals in the population are randomly matched, giving everyone an equal chance to reproduce, and each pair is combined by means of a crossover operator that generates two offspring. Here we use the extension of the Generalized Order Crossover (GOX) proposed in [71]. After crossover, the generated offspring are improved using tabu search. Then, a tournament among the parents and their two offspring is used to select the two individuals that will pass on to the next population. Because tournament is used as the replacement operator, the crossover operator is applied unconditionally. The evolution process finishes after ev_st_it generations without improvement. So far, this search procedure essentially corresponds to the simple memetic algorithm sma from [41]. We now propose to enhance this framework with two new additional components that improve its overall performance: a heuristic seeding mechanism and an immigration operator.

3.2Tabu search algorithm

Tabu search is a local search algorithm that keeps a memory structure, called tabu list, where it stores a trace of the recently visited search space. In particular, to avoid undoing recently made moves, we store in the tabu list the inverse of the moves performed to obtain the neighbors. Our tabu list has a dynamic size, similar to the one introduced in [75], so the size of the list can vary between a lower bound tb_lst_ubr and an upper bound tb_lst_lbr. When the selected neighbor is worse (resp. better) than the current solution and the upper (resp. lower) bound has not been reached, the list’s size increases (resp. decreases) in one unit. If the selected neighbor is the best solution found so far, the list is cleared; this is similar to restarting the search from this solution. We also incorporate an aspiration criterion, so a tabu move can be executed if it improves the best solution found up to this moment. In the rare situation that all neighbors are tabu, we choose the best one, clear the tabu list and slightly change its bounds by picking a random number within a given range.

4.1Comparison with state of the art

In Table 3 we can see the results obtained with 30 runs of our metaheuristic ema including a comparison with the state-of-the-art method which, to the best of our knowledge, is sma. Each row except the last one corresponds to an instance. The instance name constitutes the first column. The next columns report for each instance the best and mean expected value of the total energy consumption together with the average runtime obtained with ema. Next to each energy and runtime value we find its relative difference with respect to the same value for sma. The last row contains the average relative differences across all instances.

We can see that, overall, the mean and best energy values obtained by ema improve 0.29% and 0.48% respectively with a reduction in CPU time of 17.87%. We can conclude, supported by running a Wilcoxon signed rank test (after rejecting normality with a Shapiro–Wilk test), that our proposal ema is significantly faster than sma while obtaining comparable results quality-wise.

A more detailed look at the results shows that for those instances with the lowest flexibility ema performs worse both in quality and in time. This can be explained because instances with low flexibility are closer to the traditional job shop, i.e., there exists little variation for executing operations in different resources, and thus their search space is smaller. Compared to the parameter setting for sma presented at [41] we have decreased the initial size of the population ev_pop_size, since the immigration operator in ema introduces new individuals during the search to increase diversity. Less flexible instances do not seem to need the extra diversity we are seeking to obtain with the immigration operator and instead benefit from a larger initial population that is intensified from the start of the search. However, for instances with higher flexibility, sma improves both in quality and in time. Besides, for larger instances that also have a higher ratio n/m, sma obtains a bigger improvement in quality than for instances with a lower ratio. For example, instances 18a and 15a improve 1.65% and 1.67% in quality respectively whereas 12a and 09a only improve 0.87%. However, the improvement in time is the opposite, 12a and 09a improve 56.2% and 46% while 18a and 15a only improve 33% and 36.7%. This suggests that the ratio is an important difficulty indicator, as is in the traditional job shop. This behaviour is represented in Fig. 2, where above the name of each instance we indicate the relative improvement in quality and below, the relative improvement in time.

In [41] it is shown that it is precisely in those instances with large size and high flexibility where exact methods (in particular, a constraint programming solver) struggle to find a solution and metaheuristics become more relevant. In consequence, the improvements obtained with ema are very significant. Moreover, we have compared instances in terms of relative error, but it is important to keep in mind that a little relative worsening in a small quantity is not as significant as a big improvement in a large quantity.

4.2Synergies study

To check the synergies between the evolutionary and local search components of ema, we run individually each component for the same CPU time as the hybridized version. Finding a fair way of doing this experiment is not easy because each algorithm has different requirements and we have to find a set of parameters that make good use of the total CPU time.

Although we could have stopped the execution of each component after the average running time of ema, by doing this the search could be stopped in an improvement phase or be unnecessarily extended when it is trapped in an unpromising one. For this reason, each algorithm has been individually configured in such a way that the average execution time is similar to that of ema but giving each run independence to last longer or shorter, as necessary. To do this, in the evolutionary component we have increased the initial population size to 30000, we add 50 individuals every 5 iterations without improving and we stop the search after 50 iterations without an improvement. In the case of the tabu search, we make 15000 restarts and stop each one after 500 iterations without an improvement. All other parameters are kept as in the main experiment.

Table 4 reports the relative difference between the results obtained with each component of ema separately and the results obtained with the complete algorithm, already shown in Table 3. These relative difference values are also depicted as a bar chart in Fig. 3. Only relative differences are shown because neither the evolutionary algorithm nor the tabu search on their own can match the results obtained by ema and relative differences provide a better insight into how much the results have been altered by using only one of the search strategies.

Figure 3.

Synergy study.

In the case of the tabu search, it gives worse results because local search is very dependent on the starting point, so it is crucial to provide it with a good initial solution. Choosing these initial solutions at random, even with thousands of runs, there is a very low probability of finding a good one. In ema the evolutionary component is not so much responsible for finding good solutions, but for finding good starting points for the tabu search. This is something evolutionary algorithms excel at because they explore different areas of the search space. On the other hand, the evolutionary algorithm lacks a proper intensification mechanism, so even if it is able to explore different areas, it fails at performing a deeper search in the most promising ones. This can be clearly seen as in smaller instances with little flexibility 07a and 10a the pure evolutionary algorithm can compete with the tabu search, but as instances get harder differences become bigger. In smaller search spaces a shallow exploration yields acceptable results but as the search space increases, the need of a proper intensification becomes evident. We can also conclude that the tabu search clearly dominates the evolutionary module and its relative difference with respect to the memetic ema is pretty constant across instances.

Figure 4.

Effect of the filter.

Figure 5.

Comparison between neighborhood functions.

4.3Filtering mechanism

In our next experiment we evaluate empirically how many neighbors are filtered using the mechanism described in Section 3.2.2. First of all, it is important to note that the used lower bound is interesting because it is very fast to calculate, as we do not need to propagate the changes made by the neighbors. It is also important to remember that this mechanism leads to the same selection as using the exact value, so there is not any difference as far as the search process is concerned, only in the time taken by the method. In other words, neither the diversification nor the intensification in the algorithm are altered, even if neighbor selection becomes more efficient. In Table 5 we can see for each instance the total number of generated neighbors (column 2) together with the number of neighbors exactly evaluated by ema to obtain their real energy value (column 3) and the percentage of discarded neighbors with respect to the total (column 3). The last row provides the average across all instances of this percentage of neighbors discarded by the filtering mechanism, which is also illustrated in Fig. 4.

The results clearly show that the filtering mechanism has a strong effect, since it enables ema to discard over 95% of the neighbors. This plays a crucial role in the speed of the algorithm. It is also interesting to notice that the effect is greater as flexibility increases. Larger flexibility means a larger search space and therefore ema generates more neighbors, corresponding to alternative resource assignments. However, many of these reassignments are likely to be bad because they correspond either with slow resources with high execution time or inefficient resources with high energy consumption or even both. In any case, these alternatives must be checked and discarded and this is what the filtering mechanism achieves with undoubted speed and effectiveness.

In the light of these results, one may wonder if some of the neighborhoods are considering numerous bad-quality movements, thus increasing the number of generated neighbors without a real effect on the search, and if similar results could be achieved without using them. This leads to our last experiment.

4.4Neighborhoods comparison

The goal of this last experiment is to assess the performance of the different neighborhoods defined in Section 3.2.1. To this end, we compare the results of running ema with the tabu search using any possible combination of N𝑀𝐶𝐸𝑇, N𝑀𝐶𝑂𝑅𝑅 and N𝑂𝑃𝐸𝑅𝑅 as neighborhood function. As it was the case in the synergies study, comparing different neighborhoods can be difficult because they may have different characteristics. In this work, we have opted for changing only the neighborhoods while keeping all other parameters of the search unaltered. Given that ema uses an adaptive stopping criterion in the tabu search, if one of the neighborhoods takes longer to converge, it will be given longer running time.

To evaluate the results obtained with the different neighborhood combinations, we take as reference a lower bound of the optimal energy for each instance obtained with a constraint programming solver in [41]. Figure 5 depicts the relative difference between this lower bound and the average and best energy value obtained by ema with the different combinations of neighborhoods. If we pay attention to the average results, N𝑀𝐶𝐸𝑇 is the best standalone neighborhood. Notice however that this advantage is mainly due to those instances with low flexibility. Excluding these instances from the comparison, N𝑀𝐶𝑂𝑅𝑅 yields better results. This behavior can be explained because N𝑀𝐶𝐸𝑇 only optimizes passive energy as it does not change resource assignment. Therefore, on instances with low flexibility and few possible resource assignments, N𝑀𝐶𝐸𝑇 can obtain pretty good results starting from the assignments introduced in the initial population using heuristic seeding. On the other hand, although N𝑀𝐶𝑂𝑅𝑅 is designed to reduce passive energy consumption, it has a side effect of reducing active energy. This is because it moves operations from one resource to another and the tabu search takes the best non-forbidden neighbor. In the case of N𝑂𝑃𝐸𝑅𝑅, it can also reduce the passive energy of a solution, but when used on its own it has a very limited effect, due to the fact that it only considers the most efficient resource for a given operation. In fact, it is clear from the results that for instances with low flexibility including this neighborhood has a negative effect. Again, this can be explained because with low flexibility the algorithm does not need to focus on moving operations to other resources, so those movements are essentially introducing noise in the search. We can also conclude that the union of N𝑀𝐶𝐸𝑇∪N𝑀𝐶𝑂𝑅𝑅 yields competitive results, near those obtained with the complete neighborhood N, with the addition of N𝑂𝑃𝐸𝑅𝑅 becoming more relevant as flexibility increases.