Preface

After writing Global Optimization Algorithms – Theory and Applications [205] during my time as PhD student a long time ago, I now want to write a more direct guide to optimization and metaheuristics. Currently, this book is in an early stage of development and work-in-progress, so expect many changes. It is available as pdf, html, epub, and azw3.

The text tries to introduce optimization in an accessible way for an audience of undergraduate and graduate students without background in the field. It tries to provide an intuition about how optimization algorithms work in practice, what things to look for when solving a problem, or how to get from a simple, working, proof-of-concept approach to an efficient solution for a given problem. We follow a learning-by-doing approach by trying to solve one practical optimization problem as example theme throughout the book. All algorithms are directly implemented and applied to that problem after we introduce them. This allows us to discuss their strengths and weaknesses based on actual results. We try to improve the algorithms step-by-step, moving from very simple approaches, which do not work well, to efficient metaheuristics.

@book{aitoa,
  author    = {Thomas Weise},
  title     = {An Introduction to Optimization Algorithms},
  year      = {2018--2020},
  publisher = {Institute of Applied Optimization ({IAO}),
               School of Artificial Intelligence and Big Data,
               Hefei University},
  address   = {Hefei, Anhui, China},
  url       = {http://thomasweise.github.io/aitoa/},
  edition   = {2020-12-26}
}

We use concrete examples and algorithm implementations written in Java. While originally designed for educational purposes, the code is general and may applicable to real research experiments, too. All of it is freely available in the repository thomasWeise/aitoa-code on GitHub under the MIT License. Often, we will just look at certain portions of the code, maybe parts of a class where we omit methods or member variables, or even just snippets from functions. Each source code listing is accompanied by a (src) link in the caption linking to the current full version of the file in the GitHub repository. If you discover an error in any of the examples, please file an issue.

This book is written using our automated book writing environment, which integrates GitHub, Travis CI, and docker-hub. Many of the charts and tables in the book have been generated with R scripts, whose source code is available in the aitoaEvaluate on GitHub under the MIT License, too. The experimental results are available in the repository aitoa-data. The text of the book itself is actively written and available in the repository thomasWeise/aitoa on GitHub. There, you can also submit issues, such as change requests, suggestions, errors, typos, or you can inform me that something is unclear, so that I can improve the book.

This book is released under the Attribution-NonCommercial-ShareAlike 4.0 International license (CC BY‑NC‑SA 4.0), see http://creativecommons.org/licenses/by-nc-sa/4.0/ for a summary.

1 Introduction

Today, algorithms influence a bigger and bigger part in our daily life and the economy. They support us by suggesting good decisions in a variety of fields, ranging from engineering, timetabling and scheduling, product design, over travel and logistic planning to even product or movie recommendations. They will be the most important element of the transition of our industry to smarter manufacturing and intelligent production, where they can automate a variety of tasks, as illustrated in Figure 1.1.

Optimization and Operations Research provide us with algorithms that propose good solutions to such a wide range of questions. Usually, it is applied in scenarios where we can choose from many possible options. The goal is that the algorithms propose solutions which minimize (at least) one resource requirement, be it costs, energy, space, etc. If they can do this well, they also offer another important advantage: Solutions that minimize resource consumption are often not only cheaper from an immediate economic perspective, but also better for the environment, i.e., with respect to ecological considerations.

Thus, we already know three reasons why optimization will be a key technology for the next century, which silently does its job behind the scenes:

1. Any form of intelligent production or smart manufacturing needs automated decisions. Since these decisions should be intelligent, they can only come from a process which involves optimization in one way or another.
2. In global and local competition in all branches of industry and all service sectors those institutions who can reduce their resource consumption and costs while improving product quality and production efficiency will have the edge. One key technology for achieving this is better planning via optimization.
3. Our world suffers from both depleting resources and too much pollution. Optimization can give us more while needing less. It often inherently leads to more environmentally friendly processes.

But how can algorithms help us to find solutions for hard problems in a variety of different fields? What does variety even mean? How general are these algorithms? And how can they help us to make good decisions? And how can they help us to save resources?

In this book, we will try to answer all of these questions. We will explore quite a lot of different optimization algorithms. We will look at their actual implementations and we will apply them to example problems to see what their strengths and weaknesses are.

1.1 Examples

Let us first get a feeling about typical use cases of optimization.

1.1.1 Example: Layout of Factories

When we think about intelligent production, then there both dynamic and static aspects, as well as all sorts of nuances in between. The question how a factory should look like is a rather static, but quite important aspect. Let us assume we own a company and buy a plot of land to construct a new factory. Of course we know which products we will produce in this factory. We therefore also know the set of facilities that we need to construct, i.e., the workshops, storage depots, and maybe an administrative building. What we need to decide is where to place them on our land, as illustrated in Figure 1.2.

Let us assume we have n locations on our plot of land that we can use for our n facilities. In some locations, there might be already buildings, in others, we may need to construct them anew. For each facility and location pair, a different construction cost may arise. A location with an existing shed might be a good solution to put a warehouse in but may need to be demolished if we want to put the administration department there.

But there also are costs arising from the relative distances between the facilities that we wish to place. Maybe there is a lot of material flow between two workshops. Some finished products and raw material may need to be transported between a workshop and the storage depot. Between the administration building and the workshops, on the other hand, there will usually be no material flow. Of course, the distance between two facilities will depend on the locations we pick for them.
For each pair of facilities that we place on the map, flow costs will arise as a function (e.g., the product) of the amount of material to be transported between them and the distance of their locations.

The total cost of an assignment of facilities to locations is therefore the sum of the resulting base costs and flow costs. Our goal would be to find the assignment with the smallest possible total cost.

This scenario is called quadratic assignment problem (QAP) [31]. It has been subject to research since the 1950s [20]. QAPs appear in wide variety of scenarios such as the location of facilities on a plot of land or the placement of work stations on the factory floor. But even if we need to place components on a circuit board in a way that minimizes the total wire length, we basically have a QAP, too [190]! Despite being relatively simple to understand, the QAP is hard to solve [176].

1.1.2 Example: Route Planning for a Logistics Company

Another, more dynamic application area for optimization is logistics. Let us look at a typical real-world scenario from this field [211,212]: the situation of a logistics company that fulfills delivery tasks for its clients. A client can order one or multiple containers to be delivered to her location within a certain time window. She will fill the containers with goods, which are then to be transported to a destination location, again within a certain time window. The logistics company may receive many such customer orders per day, maybe several hundreds or even thousands. The company may have multiple depots, where containers and trucks are stored. For each order, it needs to decide which container(s) to use and how to get them to the customer, as sketched in Figure 1.3. The trucks it owns may have different capacities and could, e.g., carry either one or two containers. Besides using trucks, which can travel freely on the map, it may also be possible to utilize trains. Trains have higher capacities and can carry many containers. Different from trucks, they must follow specific schedules. They arrive and depart at fixed times to/from fixed locations. For each possible vehicle, different costs could occur. Containers can be exchanged between different vehicles at locations such as parking lots, depots, or train stations.

The company could have the goals to fulfill all transportation requests at the lowest cost. Actually, it might seek to maximize its profit, which could even mean to outsource some tasks to other companies. The goal of optimization then would be to find the assignment of containers to delivery orders and vehicles and of vehicles to routes, which maximizes the profit. And it should do so within a limited, feasible time.

Of course, there is a wide variety of possible logistics planning tasks. Besides our real-world example above, a classical task is the Traveling Salesman Problem (TSP) [8,94,133], where the goal is to find the shortest round-trip tour throughn cities, as sketched in Figure 1.4. Many other scenarios can be modeled as such logistics questions, too: If a robot arm needs to several drill holes into a circuit board, finding the shortest tour means solving a TSP and will speed up the production process, for instance [91].

1.1.3 Example: Packing, Cutting Stock, and Knapsack

Let’s say that your family is moving to a new home in another city. This means that you need to transport all of your belongings from your old to your new place, your PC, your clothes, maybe some furniture, a washing machine, and a fridge, as sketched in Figure 1.5. You cannot pack everything into your car at once, so you will have to drive back and forth a couple of times. But how often will you have to drive? Packing problems [68,180] aim to package sets of objects into containers as efficient as possible, i.e., in such a way that we need as few containers as possible. Your car can be thought of as a container and whenever it is filled, you drive to the new flat. If you need to fill the container four times, then you have to drive back and forth four times.

Such bin packing problems exist in many variants and are very related to cutting stock problems [68]. They can be one-dimensional [56], for example if we want to transport dense/heavy objects with a truck where the maximum load weight is limiting factor while there is enough space capacity. This is similar to having a company which puts network cables into people’s homes and therefore bulk purchases reels with 100m of cables each. Of course, each home needs a different required total length of cables and we want to cut our cables such that we need as few reels as possible.

A two-dimensional variant [136] could correspond to printing a set of (rectangular) images of different sizes on (rectangular) paper. Assume that more than one image fits on a sheet of paper but we have too many images for one piece of paper. We can cut the paper after printing to separate the single images. We then would like to arrange the images such that we need as few sheets of paper as possible.

The three-dimensional variant then corresponds to our moving homes scenario. Of course, there are many more different variants – the objects we want to pack could be circular, rectangular, or have an arbitrary shape. We may also have a limited number of containers and thus may not be able to pack all objects, in which case we would like to only package those that give us the most profit (arriving at a task called knapsack problem [141]).

1.1.4 Example: Job Shop Scheduling Problem

Another typical optimization task arises in manufacturing, namely the assignment (scheduling) of tasks (jobs) to machines in order to optimize a given performance criterion (objective). Scheduling [165,166] is one of the most active areas of operational research for more than six decades.

In the Job Shop Scheduling Problem (JSSP) [33,43,66,88,132,134], we have a factory (shop) with several machines. We receive a set of customer orders for products which we have to produce. We know the exact sequence in which each product/order needs to pass through the machines and how long it will need at each machine. Each production job has one sub-job (operation) for each machine on which it needs to be processed. These operations must be performed in the right sequence. Of course, no machine can process more than one operation at the same time. While we must obey these constraints, we can decide about the time at which each of the operations should begin. Often, we are looking for the starting times that lead to the earliest completion of all jobs, i.e., the shortest makespan.

Such a scenario is sketched in Figure 1.6, where four orders for different types of shoe should be produced. The resulting jobs pass through different workshops (or machines, if you want) in different order. Some, like the green sneakers, only need to be processed by a subset of the workshops.

This general scenario encompasses many simpler problems. For example, if we only produce one single product, then all jobs would pass through the same machines in the same order. Customers may be able to order different quantities of the product, so the operations of the different jobs for the same machine may need different amounts of time. This is the so-called Flow Shop Scheduling Problem (FSSP) – and it has been defined back in 1954 [121]!

Clearly, since the JSSP allows for an arbitrary machine order per job, being able to solve the JSSP would also enable us to solve the FSSP, where the machine order is fixed. We will introduce the JSSP in detail in Section 2.2.2 and use it as the main example in this book on which we will step-by-step exercise different optimization methods.

1.1.5 Summary

The four examples we have discussed so far are, actually, quite related. They all fit into the broad areas of operational research, which and smart manufacturing [52,106]. The goal of smart manufacturing is to optimize development, production, and logistics in the industry. Therefore, computer control is applied to achieve high levels of adaptability in the multi-phase process of creating a product from raw material. The manufacturing processes and maybe even whole supply chains are networked. The need for flexibility and a large degree of automation require automatic intelligent decisions. The key technology necessary to propose such decisions are optimization algorithms. In a perfect world, the whole production process as well as the warehousing, packaging, and logistics of final and intermediate products would take place in an optimized manner. No time or resources would be wasted as production gets cleaner, faster, and cheaper while the quality increases.

1.2 Metaheuristics: Why do we need them?

The main topic of this book will be metaheuristic optimization (although I will eventually also discuss some other methods (remember: work in progress). So why do we need metaheuristic algorithms? Why should you read this book?

1.2.1 Good Solutions within Acceptable Time

The first and foremost reason is that they can provide us good solutions within reasonable time. It is easy to understand that there are some problems which are harder to solve than others. Everyone of us already knows this from the mathematics classes in school. Of course, the example problems discussed before cannot be attacked as easily as solving a single linear equation. They require algorithms, they require computer science.

Ever since primary school, we have learned many problems and types of equations that we can solve. Unfortunately, theoretical computer science shows that for many problems, the time we need to find the best-possible solution can grow exponentially with the number of involved variables in the worst case. The number of involved variables here could be the number of cities in a TSP, the number of jobs or machines in a JSSP, or the number of objects to pack in a, well, packing problem. A big group of such complicated problems are called \mathcal{NP}‑hard [37,134]. Unless some unlikely breakthrough happens [44,124], there will be many problems that we cannot always solve exactly within reasonable time. Each and every one of the example problems discussed belongs to this type!

As sketched in Figure 1.7, the exponential function rises very quickly. One idea would be to buy more computers for bigger problems and to simply parallelize the computation. Well, parallelization can provide a linear speed-up at best, but we are dealing with problems where the runtime requirements may double every time we add a new decision variable. And no: Quantum computers are not the answer. Most likely, they cannot even solve these problems qualitatively faster either [1].

So what can we do to solve such problems? The exponential time requirement occurs if we make guarantees about the solution quality, especially about its optimality, over all possible scenarios. What we can do, therefore, is that we can trade-in the guarantee of finding the best possible solution for lower runtime requirements. We can use algorithms from which we hope that they find a good approximation of the optimum, i.e., a solution which is very good with respect to the objective function, but which do not guarantee that their result will be the best possible solution. We may sometimes be lucky and even find the optimum, while in other cases, we may get a solution which is close enough. And we will get this within acceptable time limits.

In Figure 1.8 we illustrate this idea on the example of the Traveling Salesman Problem [8,94,133] briefly mentioned in Section 1.1.2. The goal of solving the TSP is to find the shortest round trip tour through n cities. The TSP is \mathcal{NP}‑hard [79,94]. Today, it is possible to solve many large instances of this problem to optimality by using sophisticated exact algorithms [45,46]. Yet, finding the shortest possible tour for a particular TSP might still take many years if you are unlucky. Finding just one tour is, however, very easy: I can write down the cities in any particular order. Of course, I can visit the cities in an arbitrary order. That is an entirely valid solution, and I can obtain it basically in 0 time. This tour would probably be very bad, very long, and generally not a good idea.

In the real world, we need something in between. We need a solution which is as good as possible as fast as possible. Heuristic and metaheuristic algorithms offer different trade-offs of solution quality and runtime. Different from exact algorithms, they do not guarantee to find the optimal solution and often make no guarantee about the solution quality at all. Still, they often allow us to get very good solutions for computationally hard problems in short time. They may often still discover them (just not always, not guaranteed).

1.2.2 Good Solutions within Acceptable Time

Saying that we need a good algorithm to solve a given problem is very easy. Developing a good algorithm to solve a given problem is not, as any graduate student in the field can probably confirm. Before, I stated that great exact algorithms for the TSP exist [45,46], that can solve many TSPs quickly (although not all). There are years and years of research in these algorithms. Even the top heuristic and metaheuristic algorithm for the TSP today result from many years of targeted research [105,153,219] and their implementation from the algorithm specification alone can take months [215]. Unfortunately, if you do not have plain TSP, but one with some additional constraints – say, time windows to visit certain cities – the optimized, state-of-the-art TSP solvers are no longer applicable. And in a real-world application scenario, you do not have years to develop an algorithm. What you need are simple, versatile, general algorithm concepts that you can easily adapt to your problem at hand. Something that can be turned into a working prototype within a few weeks.

Metaheuristics are the answer. They are general algorithm concepts into which we can plug problem-specific modules. General metaheuristics are usually fairly easy to implement and deliver acceptable results. Once a sufficiently well-performing prototype has been obtained, we could go and integrate it into the software ecosystem of the customer. We also can try to improve its performance using different ideas … and years and years of blissful research, if we are lucky enough to find someone paying for it.

2 The Structure of Optimization

2.1 Introduction

From the examples that we have seen, we know that optimization problems come in different shapes and forms. Without training, it is not directly clear how to identify, define, understand, or solve them. The goal of this chapter is to bring some order into this mess. We will approach an optimization task step-by-step by formalizing its components, which will then allow us to apply efficient algorithms to it. This structure of optimization is a blueprint that can be applied in many different scenarios as basis to apply different optimization algorithms.

First, let us clarify what optimization problems actually are.

Definition 1. An optimization problem is a situation which requires deciding for one choice from a set of possible alternatives in order to reach a predefined/required benefit at minimal costs.

Definition 1 presents an economical point of view on optimization in a rather informal manner. We can refine it to the more mathematical formulation given in Definition 2.

Definition 2. The goal of solving an optimization problem is finding an input value y^{\star}\in\mathbb{Y} from a set \mathbb{Y} of allowed values for which a function f:\mathbb{Y}\mapsto\mathbb{R} takes on the smallest value.

From these definitions, we can already deduce a set of necessary components that make up such an optimization problem. We will look at them from the perspective of a programmer:

1. First, there is the problem instance data \mathcal{I}, i.e., the concrete situation which defines the framework conditions for the solutions we try to find. This input data of the optimization algorithm is discussed in Section 2.2.
2. The second component is the data structure \mathbb{Y} representing possible solutions to the problem. This is the output of the optimization software and is discussed in Section 2.3.
3. Finally, the objective function f:\mathbb{Y}\mapsto\mathbb{R} rates the quality of the candidate solutions y\in\mathbb{Y}. This function embodies the goal that we try to achieve, e.g., (minimize) costs. It is discussed in Section 2.4.

If we want to solve a Traveling Salesmen Problem (see Section 1.1.2), then the instance data includes the names of the cities that we want to visit and a map with the information of all the roads between them (or simply a distance matrix). The candidate solution data structure could simply be a city list containing each city exactly once and prescribing the visiting order. The objective function would take such a city list as input and compute the overall tour length. It would be subject to minimization.

Usually, in order to actually practically implement an optimization approach, there often will also be

1. a search space \mathbb{X}, i.e., a simpler data structure for internal use, which can more efficiently be processed by an optimization algorithm than \mathbb{Y} (Section 2.6),
2. a representation mapping \gamma:\mathbb{X}\mapsto\mathbb{Y}, which translates points x\in\mathbb{X} from the search space \mathbb{X} to candidate solutions y\in\mathbb{Y} in the solution space \mathbb{Y} (Section 2.6),
3. search operators \textrm{searchOp}:\mathbb{X}^n\mapsto\mathbb{X}, which allow for the iterative exploration of the search space \mathbb{X} (Section 2.7), and
4. a termination criterion, which tells the optimization process when to stop (Section 2.8).

At first glance, this looks a bit complicated – but rest assured, it won’t be. We will explore all of these structural elements that make up an optimization problem in this chapter, based on a concrete example of the Job Shop Scheduling Problem (JSSP) from Section 1.1.4 [33,66,88,132,134]. This example should give a reasonable idea about how the structural elements and formal definitions involved in optimization can be realized in practice. While any actual optimization problem can require very different data structures and operations from what we will discuss here, the general approach and ideas that we will discuss on specific examples should carry over to many scenarios.

At this point, I would like to make explicitly clear that the goal of this book is NOT to solve the JSSP particularly well. Our goal is to have an easy-to-understand yet practical introduction to optimization. This means that sometimes I will intentionally and knowingly choose an easy-to-understand approach, algorithm, or data structure over a better but more complicated one. Also, our aim is to nurture the general ability to come up with a solution approach to a new optimization problem within a reasonably short time, i.e., without being able to conduct research over several years. That being said, the algorithms and approaches discussed here are not necessarily inefficient. While having much room for improvement, we eventually reach approaches that find decent solutions (see, e.g., Section 3.5.5).

2.2 Problem Instance Data

2.2.1 Definitions

We implicitly distinguish optimization problems (see Definition 2) from problem instances. While an optimization problem is the general blueprint of the tasks, e.g., the goal of scheduling production jobs to machines, the problem instance is a concrete scenario of the task, e.g., a concrete lists of tasks, requirements, and machines.

Definition 3. A concrete instantiation of all information that are relevant from the perspective of solving an optimization problems is called a problem instance \mathcal{I}.

The problem instance is the input of the optimization algorithms. A problem instance is related to an optimization problem in the same way an object/instance is related to its class in an object-oriented programming language like Java or a struct in C. The class defines which member variables exists and what their valid ranges are. An instance of the class is a piece of memory which holds concrete values for each member variable.

2.2.2 Example: Job Shop Scheduling

2.2.2.1 JSSP Instance Structure

So how can we characterize a JSSP instance \mathcal{I}? In the a basic and yet general scenario [66,88,132,134], our factory has m\in\mathbb{N}_1 machines.¹ At each point in time, a machine can either work on exactly one job or do nothing (be idle). There are n\in\mathbb{N}_1 jobs that we need to schedule to these machines. For the sake of simplicity and for agreement between our notation here, the Java source code, and the example instances that we will use, we reference jobs and machines with zero-based indices from 0\ldots(n-1) and 0\ldots(m-1), respectively.

Each of the n jobs is composed of m sub-jobs – the operations – one for each machine. Each job may need to pass through these machines in a different order. The operation j of job i must be executed on machine M_{i,j}\in 0\ldots(m-1) and doing so needs T_{i,j}\in\mathbb{N}_0 time units for completion.

This definition at first seems strange, but upon closer inspection is quite versatile. Assume that we have a factory that produces exactly one product, but different customers may order different quantities. Here, we would have JSSP instances where all jobs need to be processed by exactly the same machines in exactly the same sequence. In this case M_{i_1,j}=M_{i_2,j} would hold for all jobs i_1 and i_2 and all operation indices j. The jobs would pass through all machines in the same order but may have different processing times (due to the different quantities).

We may also have scenarios where customers can order different types of products, say the same liquid soap, but either in bottles or big cannisters. Then, different machines may be needed for different orders. This is similar to the situation illustrated in Figure 1.6, where a certain job i does not need to be executed on a machine j'. We then can simply set the required time T_{i,j} to 0 for the operation j with M_{i,j}=j'.

In other words, the JSSP instance structure described here already encompasses a wide variety of real-world production situations. This means that if we can build an algorithm which can solve this general type of JSSP well, it can also automatically solve the above-mentioned special cases.

2.2.2.2 Sources for JSSP Instances

In order to practically play around with optimization algorithms, we need some concrete instances of the JSSP. Luckily, the optimization community provides benchmark instances for many different optimization problems. Such common, well-known instances are important, because they allow researchers to compare their algorithms. The eight classical and most commonly used sets of benchmark instances are published in [3,9,57,66,74,135,191,224]. Their data can be found (sometimes partially) in several repositories in the internet, such as

• the OR‑Library managed by Beasley [18],
• the comprehensive set of JSSP instances provided by van Hoorn [199,201], where also state-of-the-art results are listed,
• Oleg Shylo’s Page [185], which, too, contains up-to-date experimental results,
• Éric Taillard’s Page, or, finally,
• my own repository jsspInstancesAndResults [207], where I collect all the above problem instances and many results from existing works.

We will try to solve JSSP instances obtained from these collections. They will serve as illustrative example of how to approach optimization problems. In order to keep the example and analysis simple, we will focus on only four instances, namely

1. instance abz7 by Adams et al. [3] with 20 jobs and 15 machines
2. instance la24 by Lawrence [135] with 15 jobs and 10 machines,
3. instance swv15 by Storer et al. [191] with 50 jobs and 10 machines, and
4. instance yn4 by Yamada and Nakano [224] with 20 jobs and 20 machines.

These instances are contained in text files available at http://people.brunel.ac.uk/~mastjjb/jeb/orlib/files/jobshop1.txt, http://raw.githubusercontent.com/thomasWeise/jsspInstancesAndResults/master/data-raw/instance-data/instance_data.txt, and in http://jobshop.jjvh.nl/. Of course, if we really want to solve a new type of problem, we will usually use many benchmark problem instances to get a good understand about the performance of our algorithm(s). Only for the sake of clarity of presentation, we will here limit ourselves to these above four problems.

2.2.2.3 File Format and demo Instance

For the sake of simplicity, we created one additional, smaller instance to describe the format of these files, as illustrated in Figure 2.1.

In the simple text format used in OR‑Library, several problem instances can be contained in one file. Each problem instance \mathcal{I} is starts and ends with a line of several + characters. The next line is a short description or title of the instance. In the third line, the number n of jobs is specified, followed by the number m of machines. The actual IDs or indexes of machines and jobs are 0-based, similar to array indexes in Java. The JSSP instance definition is completed by n lines of text, each of which specifying the operations of one job i\in0\ldots(n-1). Each operation j is specified as a pair of two numbers, the ID M_{i,j} of the machine that is to be used (violet), from the interval 0\ldots(m-1), followed by the number of time units T_{i,j} the job will take on that machine. The order of the operations defines exactly the order in which the job needs to be passed through the machines. Of course, each machine can only process at most one job at a time.

In our demo instance illustrated in Figure 2.1, this means that we have n=4 jobs and m=5 machines. Job 0 first needs to be processed by machine 0 for 10 time units, it then goes to machine 1 for 20 time units, then to machine 2 for 20 time units, then to machine 3 for 40 time units, and finally to machine 4 for 10 time units. This job will thus take at least 100 time units to be completed, if it can be scheduled without any delay or waiting period, i.e., if all of its operations can directly be processed by their corresponding machines. Job 3 first needs to be processed by machine 4 for 50 time units, then by machine 3 for 30 time units, then by machine 2 for 15 time units, then by machine 0 for 20 time units, and finally by machine 1 for 15 time units. It would not be allowed to first send Job 3 to any machine different from machine 4 and after being processed by machine 4, it must be processed by machine 3 – althoug it may be possible that it has to wait for some time, if machine 3 would already be busy processing another job. In the ideal case, job 3 could be completed after 130 time units.

2.2.2.4 A Java Class for JSSP Instances

This structure of a JSSP instance can be represented by the simple Java class given in Listing 2.1.

public class JSSPInstance {
  public int m;
  public int n;
  public int[][] jobs;
}

Here, the two-dimensional array jobs directly receives the data from operation lines in the text files, i.e., each row stands for a job and contains machine IDs and processing times in an alternating sequence. The actual source file of the class JSSPInstance accompanying our book also contains additional code, e.g., for reading such data from the text file, which we have omitted here as it is unimportant for the understanding of the scenario.

2.3 The Solution Space

2.3.1 Definitions

As stated in Definition 1, an optimization problem asks us to make a choice between different possible solutions. We call them candidate solutions.

Definition 4. A candidate solution y is one potential solution of an optimization problem.

Definition 5. The solution space \mathbb{Y} of an optimization problem is the set of all of its candidate solutions y\in\mathbb{Y}.

Basically, the input of an optimization algorithm is the problem instance \mathcal{I} and the output would be (at least) one candidate solution y\in\mathbb{Y}. This candidate solution is the choice that the optimization process proposes to the human operator. It therefore holds all the data that the human operator needs to take action, in a for that the human operator can understand, interpret, and execute. During the optimization process, many such candidate solutions may be created and compared to find and return the best of them.

From the programmer’s perspective, the solution space is again a data structure, e.g., a class in Java. An instance of this data structure is the candidate solution.

2.3.2 Example: Job Shop Scheduling

What would be a candidate solution to a JSSP instance as defined in Section 2.2.2? Recall from Section 1.1.4 that our goal is to complete the jobs, i.e., the production tasks, as soon as possible. Hence, a candidate solution should tell us what to do, i.e., how to process the jobs on the machines.

2.3.2.1 Idea: Gantt Chart

This is basically what Gantt chart [126,223] are for, as illustrated in Figure 2.2. A Gantt chart defines what each of our m machines has to do at each point in time. The operations of each job are assigned to time windows on their corresponding machines.

The Gantt chart contains one row for each machine. It is to be read from left to right, where the x-axis represents the time units that have passed since the beginning of the job processing. Each colored bar in the row of a given machine stands for a job and denotes the time window during which the job is processed. The bar representing operation j of job i is painted in the row of machine M_{i,j} and its length equals the time requirement T_{i,j}.

The chart given in Figure 2.2, for instance, defines that job 0 starts at time unit 0 on machine 0 and is processed there for ten time units. Then the machine idles until the 70th time unit, at which point it begins to process job 1 for another ten time units. After 15 more time units of idling, job 3 will arrive and be processed for 20 time units. Finally, machine 0 works on job 2 (coming from machine 3) for ten time units starting at time unit 150.

Machine 1 starts its day with an idle period until job 2 arrives from machine 2 at time unit 30 and is processed for 20 time units. It then processes jobs 1 and 0 consecutively and finishes with job 3 after another idle period. And so on.

If we wanted to create a Java class to represent the complete information from a Gantt diagram, it could look like Listing 2.2. Here, for each of the m machines, we create one integer array of length 3n. Such an array stores three numbers for each of the n operations to be executed on the machine: the job ID, the start time, and the end time.

public class JSSPCandidateSolution {
  public int[][] schedule;
}

Of course, we would not strictly need a class for that, as we could as well use the integer array int[][] directly.

Also the third number, i.e., the end time, is not strictly necessary, as it can be computed based on the instance data as start+T_{i,j'} for job i on machine j after searching j' such that M_{i,j'}=j. Another form of representing a solution would be to just map each operation to a starting time, leading to m*n integer values per candidate solution [200]. But the presented structure – illustrated on an example in Figure 2.3 – is handy and easier to understand. It allows the human operator to directly see what is going on, to directly tell each machine or worker what to do and when to do it, without needing to look up any additional information from the problem instance data.

2.3.2.2 Size of the Solution Space

We choose the set of all Gantt charts for m machines and n jobs as our solution space \mathbb{Y}. Now it is not directly clear how many such Gantt charts exist, i.e., how big \mathbb{Y} is. If we allow arbitrary useless waiting times between jobs, then we could create arbitrarily many different valid Gantt charts for any problem instance. Let us therefore assume that no time is wasted by waiting unnecessarily.

There are n!=\prod_{i=1}^{n} i possible ways to arrange n jobs on one machine. n!, called the factorial of n, is the number of different permutations (or orderings) of n objects. If we have three jobs a, b, and c, then there are 3!=1*2*3=6 possible permutations, namely (a,b,c), (a,c,b), (b,a,c), (b, c, a), (c, a, b), and (c, b, a). Each permutation would equal one possible sequence in which we can process the jobs on one machine. If we have three jobs and one machine, then six is the number of possible different Gantt charts that do not waste time.

If we would have n=3 jobs and m=2 machines, we then would have (3!)^2=36 possible Gantt charts, as for each of the 6 possible sequence of jobs on the first machines, there would be 6 possible arrangements on the second machine. For m=2 machines, it is then (n!)^3, and so on. In the general case, we obtain Equation 2.1 for the size \left|\mathbb{Y}\right| of the solution space \mathbb{Y}.

\left|\mathbb{Y}\right| = (n!)^{m} \qquad(2.1)

However, the fact that we can generate (n!)^{m} possible Gantt charts without useless delay for a JSSP with n jobs and m machines does not mean that all of them are actual feasible solutions.

2.3.2.3 The Feasibility of the Solutions

Definition 6. A constraint is a rule imposed on the solution space \mathbb{Y} which can either be fulfilled or violated by a candidate solution y\in\mathbb{Y}.

Definition 7. A candidate solution y\in\mathbb{Y} is feasible if and only if it fulfills all constraints.

Definition 8. A candidate solution y\in\mathbb{Y} is infeasible if it is not feasible, i.e., if it violates at least one constraint.

In order to be a feasible solution for a JSSP instance, a Gantt chart must indeed fulfill a couple of constraints:

1. all operations of all jobs must be assigned to their respective machines and properly be completed,
2. only the jobs and machines specified by the problem instance must occur in the chart,
3. a operation will must be assigned a time window on its corresponding machine which is exactly as long as the operation needs on that machine,
4. the operations cannot intersect or overlap, each machine can only carry out one job at a time, and
5. the precedence constraints of the operations must be honored.

While the first four constraints are rather trivial, the latter one proofs problematic. Imagine a JSSP with n=2 jobs and m=2 machines. There are (2!)^2=(1*2)^2=4 possible Gantt charts. Assume that the first job needs to first be processed by machine 0 and then by machine 1, while the second job first needs to go to machine 1 and then to machine 0. A Gantt chart which assigns the first job first to machine 1 and the second job first to machine 0 cannot be executed in practice, i.e., is infeasible, as such an assignment does not honor the precedence constraints of the jobs. Instead, it contains a deadlock.

The third schedule in the first column of Figure 2.4 illustrates exactly this case. Machine 0 should begin by doing job 1. Job 1 can only start on machine 0 after it has been finished on machine 1. At machine 1, we should begin with job 0. Before job 0 can be put on machine 1, it must go through machine 0. So job 1 cannot go to machine 0 until it has passed through machine 1, but in order to be executed on machine 1, job 0 needs to be finished there first. Job 0 cannot begin on machine 1 until it has been passed through machine 0, but it cannot be executed there, because job 1 needs to be finished there first. A cyclic blockage has appeared: no job can be executed on any machine if we follow this schedule. This is called a deadlock. No jobs overlap in the schedule. All operations are assigned to proper machines and receive the right processing times. Still, the schedule is infeasible, because it cannot be executed or written down without breaking the precedence constraint.

Hence, there are only three out of four possible Gantt charts that work for this problem instance. For a problem instance where the jobs need to pass through all machines in the same sequence, however, all possible Gantt charts will work, as also illustrated in Figure 2.4. The number of actually feasible Gantt charts in \mathbb{Y} thus is different for different problem instances.

This is very annoying. The potential existence of infeasible solutions means that we cannot just pick a good element from \mathbb{Y} (according to whatever good means), we also must be sure that it is actually feasible. An optimization algorithm which might sometimes return infeasible solutions will not be acceptable.

2.3.2.4 Summary

We illustrate some examples for the number \left|\mathbb{Y}\right| of schedules which do not waste time uselessly for different values of n and m in Table 2.1. Since we use instances for testing our JSSP algorithms, we have added their settings as well and listed them in column name. Of course, there are infinitely many JSSP instances for a given setting of n and m and our instances always only mark single examples for them.

We find that even small problems with m=5 machines and n=5 jobs already have billions of possible solutions. The four more realistic problem instances which we will try to solve here already have more solutions that what we could ever enumerate, list, or store with any conceivable hardware or computer. As we cannot simply test all possible solutions and pick the best one, we will need some more sophisticated algorithms to solve these problems. This is what we will discuss in the following.

The number \#\text{feasible} of possible feasible Gantt charts can be different, depending on the problem instance. For one setting of m and n, we are interested in the minimum \min(\#\text{feasible}) of this number, i.e., the smallest value that \#\text{feasible} can take on over all possible instances with n jobs and m machines. It is not so easy to find a formula for this minimum, so we won’t do this here. Instead, in Table 2.1, we provided the corresponding numbers for a few selected instances. We find that, if we are unlucky, most of the possible Gantt charts for a problem instance might be infeasible, as \min(\#\text{feasible}) can be much smaller than \left|\mathbb{Y}\right|.

2.4 Objective Function

We now know the most important input and output data for an optimization algorithm: the problem instances \mathcal{I} and candidate solutions y\in\mathbb{Y}, respectively. But we do not just want to produce some output, not just want to find any candidate solution – we want to find the good ones. For this, we need a measure rating the solution quality.

2.4.1 Definitions

Definition 9. An objective function f:\mathbb{Y}\mapsto\mathbb{R} rates the quality of a candidate solution y\in\mathbb{Y} from the solution space \mathbb{Y} as real number.

Definition 10. An objective value f(y) of the candidate solution y\in\mathbb{Y} is the value that the objective function f takes on for y.

Without loss of generality, we assume that all objective functions are subject to minimization, meaning that smaller objective values are better. In this case, a candidate solution y_1\in\mathbb{Y} is better than another candidate solution y_2\in\mathbb{Y} if and only if f(y_1)<f(y_2). If f(y_1)>f(y_2), then y_2 would be better and for f(y_1)=f(y_2), there would be no benefit of either solution over the other, at least from the perspective of the optimization criterion f. The minimization scenario fits to situations where f represents a cost, a time requirement, or, in general, any number of required resources.

Maximization problems, i.e., where the candidate solution with the higher objective value is better, are problems where the objective function represents profits, gains, or any other form of positive output or result of a scenario. Maximization and minimization problems can be converted to each other by simply negating the objective function. In other words, if f is the objective function of a minimization problem, we can solve the maximization problem with -f and get the same result, and vice versa.

From the perspective of a Java programmer, the general concept of objective functions can be represented by the interface given in Listing 2.3. The evaluate function of this interface accepts one instance of the solution space class Y and returns a double value. doubles are floating point numbers in Java, i.e., represent a subset of the real numbers. We keep the interface generic, so that we can implement it for arbitrary solution spaces. Any actual objective function would then be an implementation of that interface.

public interface IObjectiveFunction<Y> {
  double evaluate(Y y);
}

2.4.2 Example: Job Shop Scheduling

As stated in Section 1.1.4, our goal is to complete the production jobs as soon as possible. This means that we want to minimize the makespan, the time when the last job finishes. Obviously, the smaller this value, the earlier we are done with all jobs, the better is the plan. As illustrated in Figure 2.5, the makespan is the time index of the right-most edge of any of the machine rows/schedules in the Gantt chart. In the figure, this happens to be the end time 230 of the last operation of job 0, executed on machine 4.

Our objective function f is thus equivalent to the makespan and subject to minimization. Based on our candidate solution data structure from Listing 2.2, we can easily compute f. We simply have to look at the last number in each of the integer arrays stored in the member schedule, as it represents the end time of the last job processed by a machine. We then return the largest of these numbers. We implement the interface IObjectiveFunction in class JSSPMakespanObjectiveFunction accordingly in Listing 2.4.

public class JSSPMakespanObjectiveFunction
    implements IObjectiveFunction<JSSPCandidateSolution> {
  public double evaluate(JSSPCandidateSolution y) {
    int makespan = 0;
// look at the schedule for each machine
    for (int[] machine : y.schedule) {
// the end time of the last job on the machine is the last number
// in the array, as array machine consists of "flattened" tuples
// of the form ((job, start, end), (job, start, end), ...)
      int end = machine[machine.length - 1];
      if (end > makespan) {
        makespan = end; // remember biggest end time
      }
    }
    return makespan;
  }
}

With this objective function f, subject to minimization, we have defined that a Gantt chart y_1 is better than another Gantt chart y_2 if and only if f(y_1)<f(y_2).²

2.5 Global Optima and the Lower Bound of the Objective Function

We now know the three key-components of an optimization problem. We are looking for a candidate solution y^{\star}\in\mathbb{Y} that has the best objective value f(y^{\star}) for a given problem instance \mathcal{I}. But what is the meaning best?

2.5.1 Definitions

Assume that we have a single objective function f:\mathbb{Y}\mapsto\mathbb{R} defined over a solution space \mathbb{Y}. This objective function is our primary guide during the search and we are looking for its global optima.

Definition 11. If a candidate solution y^{\star}\in\mathbb{Y} is a global optimum for an optimization problem defined over the solution space \mathbb{Y}, then there is no other candidate solution in \mathbb{Y} which is better.

Definition 12. For every global optimum y^{\star}\in\mathbb{Y} of single-objective optimization problem with solution space \mathbb{Y} and objective function f:\mathbb{Y}\mapsto\mathbb{R} subject to minimization, it holds that f(y) \geq f(y^{\star}) \forall y\in \mathbb{Y}.

Notice that Definition 12 does not state that the objective value of y^{\star} needs to be better than the objective value of all other possible solutions. The reason is that there may be more than one global optimum, in which case all of them have the same objective value. Thus, a global optimum is not defined as a candidate solutions better than all other solutions, but as a solution for which no better alternative exists.

The real-world meaning of a globally optimal is nothing else than superlative [30]. If we solve a JSSP for a factory, our goal is to find the shortest makespan. If we try to pack the factory’s products into containers, we look for the packing that needs the least amount of containers. Thus, optimization means searching for such superlatives, as illustrated in Figure 2.6. Vice versa, whenever we are looking for the cheapest, fastest, strongest, best, biggest or smallest thing, then we have an optimization problem at hand.

For example, for the JSSP there exists a simple and fast algorithm that can find the optimal schedules for problem instances with exactly m=2 machines and if all n jobs need to be processed by the two machines in exactly the same order [121]. If our application always falls into a special case of the problem, we may be lucky to find an efficient way to always solve it to optimality. The general version of the JSSP, however, is \mathcal{NP}‑hard [37,134], meaning that we cannot expect to solve it to global optimality in reasonable time. Developing a good (meta-)heuristic algorithm, which cannot provide guaranteed optimality but will give close-to-optimal solutions in practice, is a good choice.

2.5.2 Bounds of the Objective Function

If we apply an approximation algorithm, then we do not have the guarantee that the solution we get is optimal. We often do not even know if the best solution we currently have is optimal or not. In some cases, we be able to compute a lower bound \textrm{lb}(f) for the objective value of an optimal solution, i.e., we know It is not possible that any solution can have a quality better than \textrm{lb}(f), but we may not know whether a solution actually exists that has quality \textrm{lb}(f). This is not directly useful for solving the problem, but it can tell us whether our method for solving the problem is good. For instance, if we have developed an algorithm for approximately solving a given problem and we know that the qualities of the solutions we get are close to a the lower bound, then we know that our algorithm is good. We then know that improving the result quality of the algorithm may be hard, maybe even impossible, and probably not worthwhile. However, if we cannot produce solutions as good as or close to the lower quality bound, this does not necessarily mean that our algorithm is bad.

It should be noted that it is not necessary to know the bounds of objective values. Lower bounds are a nice to have feature allowing us to better understand the performance of our algorithms.

2.5.3 Example: Job Shop Scheduling

We have already defined our solution space \mathbb{Y} for the JSSP in Listing 2.2 and the objective function f in Listing 2.3. A Gantt chart with the shortest possible makespan is then a global optimum. There may be multiple globally optimal solutions, which then would all have the same makespan.

When facing a JSSP instance \mathcal{I}, we do not know whether a given Gantt chart is the globally optimal solution or not, because we do not know the shortest possible makespan. There is no direct way in which we can compute it. But we can, at least, compute some lower bound \textrm{lb}(f) for the best possible makespan.

For instance, we know that a job i needs at least as long to complete as the sum \sum_{j=0}^{m-1} T_{i,j} over the processing times of all of its operations. It is clear that no schedule can complete faster then the longest job. Furthermore, we know that the makespan of the optimal schedule also cannot be shorter than the latest finishing time of any machine j. This finishing time is at least as big as the sum b_{j} of the runtimes of all the operations assigned to this machine. But it may also include a least initial idle time a_{j}: If the operations for machine j never come first in their job, then for each job, we need to sum up the runtimes of the operations coming before the one on machine j. The least initial idle time a_{j} is then the smallest of these sums. Similarly, there is a least idle time c_{j} at the end if these operations never come last in their job. As lower bound for the fastest schedule that could theoretically exist, we therefore get:

\textrm{lb}(f) = \max\left\{\max_{i}\left\{ \sum_{j=0}^{m-1} T_{i,j}\right\} \;,\; \max_{j} \left\{ a_{j}+b_{j} +c_{j}\right\}\right\} \qquad(2.2)

More details are given in Section 6.1.1 and [66]. Often, we may not have such lower bounds, but it does never hurt to think about them, because it will provide us with some more understanding about the nature of the problem we are trying to solve.

Even if we have a lower bound for the objective function, we can usually not know whether any solution of that quality actually exists. In other words, we do not know whether it is actually possible to find a schedule whose makespan equals the lower bound. There simply may not be any way to arrange the jobs such that no operation stalls any other operation. This is why the value \textrm{lb}(f) is called lower bound: We know no solution can be better than this, but we do not know whether a solution with such minimal makespan exists.

However, if our algorithms produce solutions with a quality close to \textrm{lb}(f), we know that we are doing well. Also, if we would actually find a solution with that makespan, then we would know that we have perfectly solved the problem. The lower bounds for the makespans of our example problems are illustrated in Table 2.2.

Figure 2.7 illustrates the globally optimal solution for our small demo JSSP instance defined in Figure 2.1 (we will get to how to find such a solution later). Here we were lucky: The objective value of this solution happens to be the same as the lower bound for the makespan. Upon closer inspection, the limiting machine is the one at index 3.

We will find this by again looking at Figure 2.1. Regardless with which job we would start here, it would need to initially wait at least a_{3}=30 time units. The reason is that no first operation of any job starts at machine 3. Job 0 would get to machine 3 the earliest after 50 time units, job 1 after 30, job 2 after 62, and job 3 after again 50 time units. Also, no job in the demo instance finishes at machine 3. Job 0, for instance, needs to be processed by machine 4 for 10 time units after it has passed through machine 3. Job 1 requires 80 more time units after finishing at machine 3, job 2 also 10 time units, and job 3 again 50 time units. In other words, machine 3 needs to wait at least 30 time units before it can commence its work and will remain idle for at least 10 time units after processing the last sub job. In between, it will need to work for exactly 140 time units, the total sum of the running time of all operations assigned to it. This means that no schedule can complete faster than 30+140+10=180 time units. Thus, Figure 2.7 illustrates the optimal solution for the demo instance.

Then, all the jobs together on the machine will consume b_{3}=150 time units if we can execute them without further delay. Finally, it regardless with which job we finish on this machine, it will lead to a further waiting time of c_{3}=10 time units. This leads to a lower bound \textrm{lb}(f) of 180 and since we found the illustrated candidate solution with exactly this makespan, we have solved this (very easy) JSSP instance.

public interface IObjectiveFunction<Y> {
  double evaluate(Y y);
  default double lowerBound() {
    return Double.NEGATIVE_INFINITY;
  }
}

We can extend our interface for objective functions in Listing 2.5 to now also allow us to implement a function lowerBound which returns, well, the lower bound. If we have no idea how to compute that for a given problem instance, this function can simply return -\infty.

2.6 The Search Space and Representation Mapping

The solution space \mathbb{Y} is the data structure that makes sense from the perspective of the user, the decision maker, who will be supplied with one instance of this structure (a candidate solution) at the end of the optimization procedure. But \mathbb{Y} it not necessarily is the space that is most suitable for searching inside.

We have already seen that there are several constraints that apply to the Gantt charts. For every problem instance, different solutions may be feasible. Besides the constraints, the space of Gantt charts also looks kind of unordered, unstructured, and messy. It would be nice to have a compact, clear, and easy-to-understand representation of the candidate solutions.

2.6.1 Definitions

Definition 13. The search space \mathbb{X} is a representation of the solution space \mathbb{Y} suitable for exploration by an algorithm.

Definition 14. The elements x\in\mathbb{X} of the search space \mathbb{X} are called points in the search space.

Definition 15. The representation mapping \gamma:\mathbb{X}\mapsto\mathbb{Y} is a left-total relation which maps the points x\in\mathbb{X} of the search space \mathbb{X} to the candidate solutions y\in\mathbb{Y} in the solution space \mathbb{Y}.

For applying an optimization algorithm, we usually choose a data structure \mathbb{X} which we can understand intuitively. Ideally, it should be possible to define concepts such as distances, similarity, or neighborhoods on this data structure. Spaces that are especially suitable for searching in include, for instances:

1. subsets of s-dimensional real vectors, i.e., \mathbb{R}^s,
2. the set \mathbb{P}(s) of sequences/permutations of s objects, and
3. a number of s yes-no decisions, which can be represented as bit strings of length s and spans the space \{0,1\}^s.

For such spaces, we can relatively easily define good search methods and can rely on a large amount of existing research work and literature. If we are lucky, then our solution space \mathbb{Y} is already similar to one of these well-known and well-researched data structures. Then, we can set \mathbb{X}=\mathbb{Y} and use the identity mapping \gamma(x)=x\forall x\in\mathbb{X} as representation mapping. In other cases, we will often prefer to map \mathbb{Y} to something similar to these spaces and define \gamma accordingly.

The mapping \gamma does not need to be injective, as it may map two points x_1 and x_2 to the same candidate solution even though they are different (x_1\neq x_2). Then, there exists some redundancy in the search space. We would normally like to avoid redundancy, as it tends to slow down the optimization process [128]. Being injective is therefore a good feature for \gamma.

The mapping \gamma also does not necessarily need to be surjective, i.e., there can be candidate solutions y\in\mathbb{Y} for which no x\in\mathbb{X} with \gamma(x)=y exists. However, such solutions then can never be discovered. If the optimal solution would reside in the set of such solutions to which no point in the search space can be mapped, then, well, it could not be found by the optimization process. Being surjective is therefore a good feature for \gamma.

public interface IRepresentationMapping<X, Y> {
  void map(Random random, X x, Y y);
}

The interface given in Listing 2.6 provides a function map which maps one point x in the search space class X to a candidate solution instance y of the solution space class Y. We define the interface as generic, because we here do not make any assumption about the nature of X and Y. This interface therefore truly corresponds to the general definition \gamma:\mathbb{X}\mapsto\mathbb{Y} of the representation mapping. Side note: An implementation of map will overwrite whatever contents were stored in object y in the process, i.e., we assume that Y is a class whose instances can be modified.

2.6.2 Example: Job Shop Scheduling

In our JSSP example, we have developed the class JSSPCandidateSolution given in Listing 2.2 to represent the data of a Gantt chart (candidate solution). It can easily be interpreted by the user and we have defined a suitable objective function for it in Listing 2.4. Yet, it is not that clear how we can efficiently create such solutions, especially feasible ones, let alone how to search in the space of Gantt charts.³ What we would like to have is a search space \mathbb{X}, which can represent the possible candidate solutions of the problem in a more machine-tangible, algorithm-friendly way. While comprehensive overviews about different such search spaces for the JSSP can be found in [2,40,217,225], we here develop only one single idea which I find particularly appealing.

2.6.2.1 Idea: 1-dimensional Encoding

Imagine you would like to construct a Gantt chart as candidate solution for a given JSSP instance. How would you do that? Well, we know that each of the n jobs has m operations, one for each machine. We could simply begin by choosing one job and placing its first operation on the machine to which it belongs, i.e., write it into the Gantt chart. Then we again pick a job, take the first not-yet-scheduled operation of this job, and add it to the end of the row of its corresponding machine in the Gantt chart. Of course, we cannot pick a job whose operations all have already be assigned. We can continue doing this until all jobs are assigned – and we will get a valid solution.

This solution is defined by the order in which we chose the jobs. Such an order can be described as a simple, linear string of job IDs, i.e., of integer numbers. If we process such a string from the beginning to the end and step-by-step assign the jobs, we get a feasible Gantt chart as result.

The encoding and corresponding representation mapping can best be described by an example. In the demo instance, we have m=5 machines and n=4 jobs. Each job has m=5 operations that must be distributed to the machines. We use a string of length m*n=20 denoting the priority of the operations. We know the order of the operations per job as part of the problem instance data \mathcal{I}. We therefore do not need to encode it. This means that we just include each job’s id m=5 times in the string. This was the original idea: The encoding represents the order in which we assign the n jobs, and each job must be picked m times. Our search space is thus somehow similar to the set \mathbb{P}(n*m) of permutations of n*m objects mentioned earlier, but instead of permutations, we have permutations with repetitions.

A point x\in\mathbb{X} in the search space \mathbb{X} for the demo JSSP instance would thus be an integer string of length 20. As example, we chose x=(0, 2, 1, 0, 3, 1, 0, 1, 2, 3, 2, 1, 1, 2, 3, 0, 2, 0, 3, 3) in Figure 2.8. The representation mapping starts with an empty Gantt chart. This string is interpreted from left to right, as illustrated in the figure. The first value is 0, which means that job 0 is assigned to a machine first. From the instance data, we know that job 0 first must be executed for 10 time units on machine 0. The job is thus inserted on machine 0 in the chart. Since machine 0 is initially idle, it can be placed at time index 0. We also know that this operation can definitely be executed, i.e., won’t cause a deadlock, because it is the first operation of the job.

The next number in the string is 2, so job 2 is next. This job needs to go for 30 time units to machine 2, which also is initially idle. So it can be inserted into the candidate solution directly as well (and cannot cause any deadlock either).

Then job 1 is next in x, and from the instance data we can see that it will go to machine 1 for 20 time units. This machine is idle as well, so the job can start immediately.

We now encounter job 0 again in the integer string. Since we have already performed the first operation of job 0, we now would like to schedule its second operation. According to the instance data, the second operation takes place on machine 1 and will take 20 time units. We know that completing the first operation took 10 time units. We also know that machine 1 first has to process job 1 for 20 time units. The earliest possible time at which we can begin with the second operation of job 0 is thus at time unit 20, namely the bigger one of the above two values. This means that job 0 has to wait for 10 time units after completing its first operation and then can be processed by machine 1. No deadlock can occur, as we made sure that the first operation of job 0 has been scheduled before the second one.

We now encounter job 3 in the integer string, and we know that job 3 first goes to machine 4, which currently is idle. It can thus directly be placed on machine 4, which it will occupy for 50 time units.

Then we again encounter job 1 in the integer string. Job 1 should, in its second operation, go to machine 0. Its first operation to 20 time units on machine 1, while machine 0 was occupied for 10 time units by job 0. We can thus start the second operation of job 1 directly at time index 20.

Further processing of y leads us to job 0 again, which means we will need to schedule its third operation, which will need 20 time units on machine 2. Machine 2 is occupied by job 2 from time unit 0 to 30 and becomes idle thereafter. The second operation of job 0 finishes on time index 40 at machine 1. Hence, we can begin with the third operation at time index 40 at machine 2, which had to idle for 10 time units.

We continue this iterative processing until reaching the end of the string x. We now have constructed the complete Gantt chart y illustrated in Figure 2.8. Whenever we assign a operation i>0 of any given job to a machine, then we already had assigned all operations at smaller indices first. No deadlock can occur and y must therefore be feasible.

public class JSSPRepresentationMapping implements
    IRepresentationMapping<int[], JSSPCandidateSolution> {
  public void map(Random random, int[] x,
      JSSPCandidateSolution y) {
// create variables machineState, machineTime of length m and
// jobState, jobTime of length n, filled with 0 [omitted brevity]
// iterate over the jobs in the solution
    for (int nextJob : x) {
// get the definition of the steps that we need to take for
// nextJob from the instance data stored in this.m_jobs
      int[] jobSteps = this.mJobs[nextJob];
// jobState tells us the index in this list for the next step to
// do, but since the list contains machine/time pairs, we
// multiply by 2 (by left-shifting by 1)
      int jobStep = (jobState[nextJob]++) << 1;

// so we know the machine where the job needs to go next
      int machine = jobSteps[jobStep]; // get machine

// start time is maximum of the next time when the machine
// becomes idle and the time we have already spent on the job
      int start =
          Math.max(machineTime[machine], jobTime[nextJob]);
// the end time is simply the start time plus the time the job
// needs to spend on the machine
      int end = start + jobSteps[jobStep + 1]; // end time
// it holds for both the machine (it will become idle after end)
// and the job (it can go to the next station after end)
      jobTime[nextJob] = machineTime[machine] = end;

// update the schedule with the data we have just computed
      int[] schedule = y.schedule[machine];
      schedule[machineState[machine]++] = nextJob;
      schedule[machineState[machine]++] = start;
      schedule[machineState[machine]++] = end;
    }
  }
}

In Listing 2.7, we illustrate how such a mapping can be implemented. It basically is a function translating an instance of int[] to JSSPCandidateSolution. This is done by keeping track of time that has passed for each machine and each job, as well as by remembering the next operation for each job and the position in the schedule of each machine.

2.6.2.2 Advantages of a very simple Encoding

This is a very nice and natural way to represent Gantt charts with a much simpler data structure. As a result, it has been discovered by several researchers independently, the earliest being Gen et al. [81], Bierwirth [24,25], and Shi et al. [182], all in the 1990s.

But what do we gain by using this search space and representation mapping? First, well, we now have a very simple data structure \mathbb{X} to represent our candidate solutions. Second, we also have very simple rules for validating a point x in the search space: If it contains the numbers 0\ldots(n-1) each exactly m times, it represents a feasible candidate solution.

Third, the candidate solution corresponding to a valid point from the search space will always be feasible [24]. The mapping \gamma will ensure that the order of the operations per job is always observed. We do not need to worry about the issue of deadlocks mentioned in Section 2.3.2.3. We know from Table 2.1, that the vast majority of the possible Gantt charts for a given problem may be infeasible – and now we do no longer need to worry about that. Our mapping also makes sure of the more trivial constraints, such as that each machine will process at most one job at a time and that all operations are eventually processed.

Finally, we also could modify our representation mapping \gamma to adapt to more complicated and constraint versions of the JSSP if need be: For example, imagine that it would take a job- and machine-dependent time requirement for carrying a job from one machine to another, then we could facilitate this by changing \gamma so that it adds this time to the starting time of the job. If there was a job-dependent setup time for each machine [5], which could be different if job 1 follows job 0 instead of job 2, then this could be facilitated easily as well. If our operations would be assigned to machine types instead of machines and there could be more than one machine per machine type, then the representation mapping could assign the operations to the next machine of their type which becomes idle. Our representation also trivially covers the situation where each job may have more than m operations, i.e., where a job may need to cycle back and pass one machine twice. It is also suitable to simple scenarios, such as the Flow Shop Problem, where all jobs pass through the machines in the same, pre-determined order [66,80,217].

Many such different problem flavors can now be reduced to investigating the same space \mathbb{X} using the same optimization algorithms, just with different representation mappings \gamma and/or objective functions f. Additionally, it became very easy to indirectly create and modify candidate solutions by sampling points from the search space and moving to similar points, as we will see in the following chapters.

2.6.2.3 Size of the Search Space

It is relatively easy to compute the size \left|\mathbb{X}\right| of our proposed search space \mathbb{X} [182]. We do not need to make any assumptions regarding no useless waiting time, as in Section 2.3.2.2, since this is not possible by default. Each element x\in\mathbb{X} is a permutation of a multiset where each of the n elements occurs exactly m times. This means that the size of the search space can be computed as given in Equation 2.3.

\left|\mathbb{X}\right| = \frac{\left(m*n\right)!}{ \left(m!\right)^{n} } \qquad(2.3)

We give some example values for this search space size \left|\mathbb{X}\right| in Table 2.3. From the table, we can immediately see that the number of points in the search space, too, grows very quickly with both the number of jobs n and the number of machines m of an JSSP instance \mathcal{I}.

For our demo JSSP instance with n=4 jobs and m=5 machines, we already have about 12 billion different points in the search space and 7 million possible, non-wasteful candidate solutions.

We now find the drawback of our encoding: There is some redundancy in our mapping, \gamma here is not injective. If we would exchange the first three numbers in the example string in Figure 2.8, we would obtain the same Gantt chart, as jobs 0, 1, and 2 start at different machines.

As said before, we should avoid redundancy in the search space. However, here we will stick with our proposed mapping because it is very simple, it solves the problem of feasibility of candidate solutions, and it allows us to relatively easily introduce and discuss many different approaches, algorithms, and sub-algorithms.

2.7 Search Operations

One of the most important design choices of a metaheuristic optimization algorithm are the search operators employed.

2.7.1 Definitions

Definition 16. An k‑ary search operator \textrm{searchOp}:\mathbb{X}^k\mapsto\mathbb{X} is a left-total relation which accepts k points in the search space \mathbb{X} as input and returns one point in the search space as output.

Special cases of search operators are

• nullary operators (k=0, see Listing 2.8), which sample a new point from the search space without using any information from an existing points,
• unary operators (k=1, see Listing 2.9), which sample a new point from the search space based on the information of one existing point,
• binary operators (k=2, see Listing 2.10), which sample a new point from the search space by combining information from two existing points, and
• ternary operators (k=3), which sample a new point from the search space by combining information from three existing points.

public interface INullarySearchOperator<X>
    extends ISetupPrintable {
  void apply(X dest, Random random);
}

public interface IUnarySearchOperator<X>
    extends ISetupPrintable {
  void apply(X x, X dest, Random random);
}

public interface IBinarySearchOperator<X>
    extends ISetupPrintable {
  void apply(X x0, X x1, X dest, Random random);
}

Whether, which, and how such such operators are used depends on the nature of the optimization algorithms and will be discussed later on.

Search operators are often randomized, which means invoking the same operator with the same input multiple times may yield different results. This is why Listings 2.8-2.10 all accept an instance of java.util.Random, a pseudorandom number generator.

Operators that take existing points in the search space as input tend to sample new points which, in some sort, are similar to their inputs. They allow us to define proximity-based relationships over the search space, such as the common concept of neighborhoods.

Definition 17. A unary operator \textrm{searchOp}:\mathbb{X}\mapsto\mathbb{X} defines a neighborhood relationship over a search space where a point x_1\in\mathbb{X} is called a neighbor of a point x_2\in\mathbb{X} are called neighbors if and only if x_1 could be the result of an application of \textrm{searchOp} to x_2.

2.7.2 Example: Job Shop Scheduling

We will step-by-step introduce the concepts of nullary, unary, and binary search operators in the subsections of Section 3 on metaheuristics as they come. This makes more sense from a didactic perspective.

2.8 The Termination Criterion and the Problem of Measuring Time

We have seen that the search space for even small instances of the JSSP can already be quite large. We simply cannot enumerate all of them, as it would take too long. This raises the question: If we cannot look at all possible solutions, how can we find the global optimum? We may also ask: If we cannot look at all possible solutions, how can we know whether a given candidate solution is the global optimum or not? In many optimization scenarios, we can use theoretical bounds to solve that issue, but a priori, these questions are valid and their answer is simply: No. No, without any further theoretical investigation of the optimization problem, we don’t know if the best solution we know so far is the global optimum or not. This leads us to another problem: If we do not know whether we found the best-possible solution or not, how do we know if we can stop the optimization process or should continue trying to solve the problem? There are two basic answers: Either when the time is up or when we found a reasonably-good solution.

2.8.1 Definitions

Definition 18. The termination criterion is a function of the state of the optimization process which becomes true if the optimization process should stop (and then remains true) and remains false as long as it can continue.

public interface ITerminationCriterion {
  boolean shouldTerminate();
}

With a termination criterion defined as implementation of the interface given in Listing 2.11, we can embed any combination of time or solution quality limits. We could, for instance, define a goal objective value g good enough so that we can stop the optimization procedure as soon as a candidate solution y\in\mathbb{Y} has been discovered with f(y)\leq g, i.e., which is at least as good as the goal. Alternatively – or in addition – we may define a maximum amount of time the user is willing to wait for an answer, i.e., a computational budget after which we simply need to stop. Discussions of both approaches from the perspective of measuring algorithm performance are given in Sections 4.2, 4.3.

2.8.2 Example: Job Shop Scheduling

In our example domain, the JSSP, we can assume that the human operator will input the instance data \mathcal{I} into the computer. Then she may go drink a coffee and expect the results to be ready upon her return. While she does so, can we solve the problem? Unfortunately, probably not. As said, for finding the best possible solution, if we are unlucky, we would need in invest a runtime growing exponentially with the problem size, i.e., m and n [37,134]. So can we guarantee to find a solution which is, say, 1% worse, until she finishes her drink? Well, it was shown that there is no algorithm which can guarantee us to find a solution only 25% worse than the optimum within a runtime polynomial in the problem size [117,222] in 1997. Since 2011, we know that any algorithm guaranteeing to provide schedules that are only be a constant factor (be it 25% or 1’000’000) worse than the optimum may need the dreaded exponential runtime [142]. So whatever algorithm we will develop for the JSSP, defining a some limit solution quality based on the lower bound of the objective value at which we can stop makes little sense.

Hence, we should rely on the simple practical concern: The operator drinks a coffee. A termination criterion granting three minutes of runtime seems to be reasonable to me here. We should look for the algorithm implementation that can give us the best solution quality within that time window.

Of course, there may also be other constraints based on the application scenario, e.g., whether a proposed schedule can be implemented/completed within the working hours of a single day. We might let the algorithm run longer than three minutes until such a solution was discovered. But, as said before, if a very odd scenario occurs, it might take a long time to discover such a solution, if ever. The operator may also need to be given the ability to manually stop the process and extract the best-so-far solution if needed. For our benchmark instances, however, this is not relevant and we can limit ourselves to the runtime-based termination criterion.

2.9 Solving Optimization Problems

Thank you for sticking with me during this long and a bit dry introduction chapter. Why did we go through all of this long discussion? We did not even solve the JSSP yet…

Well, in the following you will see that we now are actually only a few steps away from getting good solutions for the JSSP. Or any optimization problem. Because we now have actually exercise a the basic process that we need to go through whenever we want to solve a new optimization task.

1. The first thing to do is to understand the scenario information, i.e., the input data \mathcal{I} that our program will receive.
2. The second step is to understand what our users will consider as a solution – a Gantt chart, for instance. Then we need to define a data structure \mathbb{Y} which can hold all the information of such a candidate solution.
3. Once we have the data structure \mathbb{Y} representing a complete candidate solution, we need to know when a solution is good. We will define the objective function f, which returns one number (say the makespan) for a given candidate solution.
4. If we want to apply any of the optimization algorithms introduced in the following chapters, then we also to know when to stop. As already discussed, we usually cannot solve instances of a new problem to optimality within feasible time and often don’t know whether the current-best solution is optimal or not. Hence, a termination criterion usually arises from practical constraints, such as the acceptable runtime.

All the above points need to be tackled in close collaboration with the user. The user may be the person who will eventually, well, use the software we build or at least a domain expert. The following steps then are our own responsibility:

5. In the future, we will need to generate many candidate solutions quickly, and these better be feasible. Can this be done easily using the data structure \mathbb{Y}? If yes, then we are good. If not, then we should think about whether we can define an alternative search space \mathbb{X}, a simpler data structure. Creating and modifying instances of such a simple data structure \mathbb{X} is much easier than \mathbb{Y}. Of course, defining such a data structure \mathbb{X} makes only sense if we can also define a mapping \gamma from \mathbb{X} to \mathbb{Y}.
6. We select optimization algorithms and plug in the representation and objective function. We may need to implement some other algorithmic modules, such as search operations. In the following chapters, we discuss a variety of methods for this.
7. We test, benchmark, and compare several algorithms to pick those with the best and most reliable performance (see Section 4).

3 Metaheuristic Optimization Algorithms

Optimization problems are solved by optimization algorithms. We can roughly divide these into exact and heuristic methods.

An exact algorithm guarantees to find the optimal solution if sufficient runtime is granted. This required runtime might, in the worst case, exceed what we can afford, in particular for \mathcal{NP}‑hard problems, such as the JSSP. Alternatively, many exact methods can be halted before completing their run and they can then still provide an approximate solution (without the guarantee that it is optimal).

For heuristic algorithms, this directly is the basic premise; They give us some approximate solution relatively quickly. They either do not make any guarantees at all how good it will be or, sometimes, provide some bound guarantee (like: This solution will not cost more than two times of the optimal cost.) Simple heuristics are usually tailor-made for specific problems, like the TSP or JSSP.

Definition 19. A metaheuristic is a general algorithm that can produce approximate solutions for a class of different optimization problems.

Metaheuristics [42,82,83,205] are the most important class of algorithms that we explore in this book. These algorithms have the advantage that we can easily adapt them to new optimization problems. As long as we can construct the elements discussed in Section 2 for a problem, we can attack it with a metaheuristic. We will introduce several such general algorithms in this book. We explore them by again using the Job Shop Scheduling Problem (JSSP) from Section 1.1.4 as example.

3.1 Common Characteristics

Before we delve into our first algorithms, let us first take a look on some things that all metaheuristics have in common.

3.1.1 Anytime Algorithms

Definition 20. An anytime algorithm is an algorithm which can provide an approximate result during almost any time of its execution.

All metaheuristics – and many other optimization and machine learning methods – are anytime algorithms [27]. The idea behind anytime algorithms is that they start with (potentially bad) guess about what a good solution would be. During their course, they try to improve their approximation quality, by trying to produce better and better candidate solutions. At any point in time, we can extract the current best guess about the optimum (and stop the optimization process if we want to). This fits to the optimization situation that we have discussed in Section 2.8: We often cannot find out whether the best solution we currently have is the globally optimal solution for the given problem instance or not, so we simply continue trying to improve upon it until a termination criterion tells us to quit, e.g., until the time is up.

3.1.2 Return the Best-So-Far Candidate Solution

This one is actually quite simple, yet often ignored and misunderstood by novices to the subject: Regardless what the optimization algorithm does, it will never NEVER forget the best-so-far candidate solution. Often, this is not explicitly written in the formal definition of the algorithms, but there always exists a special variable somewhere storing that solution. This variable is updated each time a better solution is found. Its value is returned when the algorithm stops.

3.1.3 Randomization

Often, metaheuristics make randomized choices. In cases where it is not clear whether doing A or doing B is better, it makes sense to simply flip a coin and do A if it is heads and B if it is tails. That our search operator interfaces in Listings 2.8-2.10 all accept a pseudorandom number generator as parameter is one manifestation of this issue. Random number generators are objects which provide functions that can return numbers from certain ranges, say from [0,1) or an integer interval. Whenever we call such a function, it may return any value from the allowed range, but we do not know which one it will be. Also, the returned value should be independent from those returned before, i.e., from known the past random numbers, we should not be able to guess the next one. By using such random number generators, we can let an algorithm make random choices, randomly pick elements from a set, or change a variable’s value in some unpredictable way.

3.1.4 Black-Box Optimization

The concept of general metaheuristics, the idea to attack a very wide class of optimization problems with one basic algorithm design, can be realized when following a black-box approach. If we want to have one algorithm that can be applied to all the examples given in in the introduction, then this can best be done if we hide all details of the problems under the hood of the structural elements introduced in Section 2. For a black-box metaheuristic, it does not matter how the objective function f works. The only thing that matters is that gives a rating of a candidate solution y\in\mathbb{Y} and that smaller ratings are better.

There are different degrees of how general black-box metaheuristics can be. For the most general algorithms, it does not matter what exactly the search operators do or even what data structure is used as search space \mathbb{X}. For them, it only matters that these operators can be used to get to new points in the search space (which can be mapped to candidate solutions y via a representation mapping \gamma whose nature is also unimportant for the metaheuristic). Then, even the nature of the candidate solutions y\in\mathbb{Y} and the solution space \mathbb{Y} play no big role for black-box optimization methods, as they only work on and explore the search space \mathbb{X}.

Then there are also black-box metaheuristics that demand a special type of search space, e.g., a specific subset of the n-dimensional real numbers (\mathbb{X}\subset\mathbb{R}^n), bit strings of a given length, or permutations of the first n natural numbers. These algorithms can still be general, as they may make very few assumptions about the nature of the objective function f defined over the solution space \mathbb{Y}.

The solution space is relevant for the human operator using the algorithm only, the search space is what the algorithm works on. Of course, in many cases, \mathbb{X}=\mathbb{Y}.

Thus, a black-box metaheuristic is a general algorithm into which we can plug representations, objective functions, and often also search operations as needed by a specific application. This is illustrated in Figure 3.1. Black-box optimization is the highest level of abstraction on which we can work when trying to solve complex problems.

3.1.5 Putting it Together: A simple API

Before, I promised that we will implement all the algorithms discussed in this book.

If we would be dealing with classical algorithms, things would be somewhat easier: A classical algorithm has clearly defined input and output data structures. Dijkstra’s shortest path algorithm [60], for instance, gets fed with a graph of weighted edges and a start node and will return the paths of minimum weight to all other nodes (or a specified target node). In Machine Learning, the situation is quite similar: We would have a lot of specialized algorithms for clearly defined situations. The input and output data would usually adhere to some basic, fixed structures. If you implement k‑means clustering [77,104], for instance, you have real vectors coming in and k real vectors going out of your algorithm and that’s that. Deep Learning [87] basically takes, as input, a set of labeled real vectors (plus a network structure) and, as output, produces the vector of weights for the network. However, we have to deal with the black-box concept, meaning that our algorithms will be very variable in terms of the data structures we can feed to them. Matter in fact: Any of the three scenarios above can be modeled as optimization problem. Any of them can be tackled with (most of) the metaheuristics in this book as well!

This is challenging from a programming perspective, especially when we try to tackle this in an educational setting, where stuff should not be overly complicated. What we want is to not just implement general algorithms, but also be able to execute them and obtain their results in a convenient fashion. Ideally, we do not want to be bothered with too much book keeping or the creation of log files and such and such. It should be possible to implement the most general type of black-box methods as well as problem-specific optimization methods and to run them in a uniform environment.

I therefore try to define a simple API for black-box optimization that combines all of our considerations far. The goal is to make the implementation of metaheuristics as simple as possible. We do this by clearly dividing between the optimization algorithms for solving a problem on one side and the structural components of the problem on the other side. The algorithms that we will implement will be general black-box methods. At the same time, we will develop the components that we need to plug into them to solve JSSPs as educational example.

We therefore first need to consider what an optimization process needs as input. Obviously, in the most common case, these are all the items we have discussed in the previous section, ranging from the termination criterion over the search operators (which we will discuss later) and the representation mapping to the objective function. Let us therefore define an interface that can provide all these components with corresponding getter methods. We call this interface IBlackBoxProcess<X,Y> from which an excerpt is given in Listing 3.1. The interface is generic, meaning it allows us to provide a search space \mathbb{X} as type parameter X and a solution space \mathbb{Y} via the type parameter Y.

public interface IBlackBoxProcess<X, Y> extends
    IObjectiveFunction<X>, ITerminationCriterion, Closeable {
  Random getRandom();
  ISpace<X> getSearchSpace();
  double getBestF();
  double getGoalF();
  void getBestX(X dest);
  void getBestY(Y dest);
  long getConsumedFEs();
  long getLastImprovementFE();
  long getMaxFEs();
  long getConsumedTime();
  long getLastImprovementTime();
  long getMaxTime();
}

Actually, such an interface does not need to expose the representation mapping \gamma and objective function f as separate components to an optimization algorithm. It is sufficient if the interface directly implements an evaluate that takes, as input, an element x\in\mathbb{X}, internally performs the representation mapping y=\gamma(x), then invokes the objective function f(y), and returns its result. This evaluate method could then even be implemented such that it remembers the best-so-far-solution. We then no longer need to keep track of it in the optimization itself.

Of course, we can also implement logging of the search progress inside of evaluate, which would make this functionality available to all of our experiments in a transparent fashion. Furthermore, we could also keep track of the total number of calls to the objective function as well as of the consumed runtime. This, in turn, can be used to implement the termination criterion.

All in all, this interface allows us to create transparent implementations that

1. provide a random number generator to the algorithm,
2. wrap an objective function f together with a representation mapping \gamma to allow us to evaluate a point in the search space x\in\mathbb{X} in a single step, effectively performing f(\gamma(x)),
3. keep track of the elapsed runtime and FEs as well as when the last improvement was made by updating said information when necessary during the invocations of the wrapped objective function,
4. keep track of the best points in the search space and solution space so far as well as their associated objective value in special variables by updating them whenever the wrapped objective function discovers an improvement (taking care of the issue from Section 3.1.2 automatically),
5. represent a termination criterion based on the above information (e.g., maximum FEs, maximum runtime, reaching a goal objective value), and
6. log the improvements that the algorithm makes to a text file, so that we can use them to make tables and draw diagrams.

Along with the interface class IBlackBoxProcess, we also provide a builder for instantiating it. The actual implementation behind this interface does not matter here. It is clear what it does, and the actual code is simple and not contributing to the understand of the algorithms or processes. Thus, you do not need to bother with it, just the assumption that an object implementing IBlackBoxProcess has the abilities listed above shall suffice here.

When instantiating this interface via the builder, we can provide the termination criterion, representation mapping, and objective function. Additionally, we also need to provide the functionality to instantiate and copy the data structures making up the spaces \mathbb{X} and \mathbb{Y}. On one hand, these are needed for the internal book-keeping so that evaluate can internally make a copy of the best-so-far elements in \mathbb{X} and \mathbb{Y}). On the other hand, the black-box optimization algorithms that we will implement also must be able to make such copies. Since we do not make any assumption about the Java classes corresponding to the spaces \mathbb{X} and \mathbb{Y}, we add another easy-to-implement and very simple interface, namely ISpace, see Listing 3.2. It can be implemented for each problem type and then provides the required functionality.

public interface ISpace<Z> {
  Z create();
  void copy(Z from, Z to);
}

Equipped with this, defining an interface for black-box metaheuristics becomes easy: The optimization algorithms themselves then are implementations of the generic interface IMetaheuristic<X,Y> given in Listing 3.3. As you can see, this interface only really needs a single method, solve. This method will use the functionality provided by the IBlackBoxProcess handed to it as parameter process. The implemented algorithm can generate new points in the search space X and send them to the evaluate method of process. This is the core behavior of every basic metaheuristic and in the rest of this chapter, we will learn how different algorithms realize it.

public interface IMetaheuristic<X, Y> extends ISetupPrintable {
  void solve(IBlackBoxProcess<X, Y> process);
}

Notice that the interface IMetaheuristic is, again, generic, allowing us to specify a search space \mathbb{X} as type parameter X and a solution space \mathbb{Y} via the type parameter Y. Whether an implementation of this interface is generic too or whether it ties down X or Y to concrete types will then depend on the algorithms we try to realize. The most general black-box metaheuristics may be able to deal with any search- and solution space, as long they are provided with the right operators. Still, we could also implement an algorithm specified to numerical problems where \mathbb{X}\subset\mathbb{R}^n, by tying down X to double[] in the algorithm class specification.

3.1.6 Example: Job Shop Scheduling

What we need to provide for our JSSP example are implementations of our ISpace interface for both the search and the solution space, which are given in Listing 3.4 and Listing 3.5, respectively. These classes implement the methods that an IBlackBoxProcess implementation needs under the hood to, e.g., copy and store candidate solutions and points in the search space.

public class JSSPSearchSpace implements ISpace<int[]> {
  public int[] create() {
    return new int[this.mLength];
  }
  public void copy(int[] from, int[] to) {
    System.arraycopy(from, 0, to, 0, this.mLength);
  }
}

public class JSSPSolutionSpace
    implements ISpace<JSSPCandidateSolution> {
  public JSSPCandidateSolution create() {
    return new JSSPCandidateSolution(this.instance.m,
        this.instance.n);
  }
  public void copy(JSSPCandidateSolution from,
      JSSPCandidateSolution to) {
    int n = this.instance.n * 3;
    int i = 0;
    for (int[] s : from.schedule) {
      System.arraycopy(s, 0, to.schedule[i++], 0, n);
    }
  }
}

With the exception of the search operators, which we will introduce when they are needed, we have already discussed how the other components needed to solve a JSSP can be realized in Section 2.3.2.1, Section 2.6.2, Section 2.4.2, and Section 2.8.2.

3.2 Random Sampling

If we have our optimization problem and its components properly defined according to Section 2, then we already have the proper tools to solve the problem. We know

• how a solution can internally be represented as point x in the search space \mathbb{X} (Section 2.6),
• how we can map such a point x\in\mathbb{X} to a candidate solution y in the solution space \mathbb{Y} (Section 2.3) via the representation mapping \gamma:\mathbb{X}\mapsto\mathbb{Y} (Section 2.6), and
• how to rate a candidate solution y\in\mathbb{Y} with the objective function f (Section 2.4).

The only question left for us to answer is how to create a point x in the search space. We can then apply \gamma(x) and get a candidate solution y whose quality we can assess via f(y). Let us look at the problem as a black box (Section 3.1.4). In other words, we do not really know structures and features make a candidate solution good. Hence, we do not know how to create a good solution either. Then the best we can do is just create the solutions randomly.

3.2.1 Ingredient: Nullary Search Operation for the JSSP

For this purpose, we need to implement the nullary search operation from Listing 2.8. We create a new search operator which needs no input and returns a point in the search space. Recall that our representation (Section 2.6.2) requires that each index i\in 0\ldots(n-1) of the n must occur exactly m times in the integer array of length m*n, where m is the number of machines in the JSSP instance. In Listing 3.6, we achieve this by first creating the sequence (n-1,n-2,\ldots,0) and then copying it m times in the destination array dest. We then randomly shuffle dest by applying the Fisher-Yates shuffle algorithm [76,129], which simply brings the array into an entirely random order.

public class JSSPNullaryOperator
    implements INullarySearchOperator<int[]> {
  public void apply(int[] dest, Random random) {
// create first sequence of jobs: n-1, n-2, ..., 0
    for (int i = this.n; (--i) >= 0;) {
      dest[i] = i;
    }
// copy this m-1 times: n-1, n-2, ..., 0, n-1, ... 0, n-1, ...
    for (int i = dest.length; (i -= this.n) > 0;) {
      System.arraycopy(dest, 0, dest, i, this.n);
    }
// now randomly shuffle the array: create a random sequence
    RandomUtils.shuffle(random, dest, 0, dest.length);
  }
}

The apply method of our implemented operator creates one random point in the JSSP search space. We can then pass this point through the representation mapping that we already implemented in Listing 2.7 and get a Gantt diagram. Easily we then obtain the quality, i.e., makespan, of this candidate solution as the right-most edge of any an job assignment in the diagram, as defined in Section 2.4.2.

3.2.2 Single Random Sample

3.2.2.1 The Algorithm

Now that we have all ingredients ready, we can test the idea. In Listing 3.7, we implement this algorithm (here called 1rs) which creates exactly one random point x in the search space. It then takes this point and passes it to the evaluation function of our black-box process, which will perform the representation mapping  y=\gamma(x) and compute the objective value f(y). Internally, we implemented this function in such a way that it automatically remembers the best candidate solution it ever has evaluated. Thus, we do not need to take care of this in our algorithm, which makes the implementation so short.

public class SingleRandomSample<X, Y>
    extends Metaheuristic0<X, Y> {
  public void solve(IBlackBoxProcess<X, Y> process) {
// Allocate data structure for holding 1 point from search space.
    X x = process.getSearchSpace().create();

// Apply the nullary operator: Fill data structure with a random
// but valid point from the search space.
    this.nullary.apply(x, process.getRandom());

// Evaluate the point: process.evaluate automatically applies the
// representation mapping and calls objective function. The
// objective value is ignored here (not stored anywhere), but
// "process" will remember the best solution. Thus, whoever
// called this "solve" function can obtain the result.
    process.evaluate(x);
  }
}

3.2.2.2 Results on the JSSP

Of course, since the algorithm is randomized, it may give us a different result every time we run it. In order to understand what kind of solution qualities we can expect, we hence have to run it a couple of times and compute result statistics. We therefore execute our program 101 times and the results are summarized in Table 3.1, which describes them using simple statistics whose meanings are discussed in-depth in Section 4.4.

What we can find in Table 3.1 is that the makespan of the best solution that any of the 101 runs has delivered for each of the four JSSP instances is roughly between 60% and 100% longer than the lower bound. The arithmetic mean and median of the solution qualities are even between 10% and 20% worse. In the Gantt charts of the median solutions depicted in Figure 3.2, we can find big gaps between the operations.

This is completely reasonable. After all, we just create a single random solution. We can hardly assume that doing all jobs of a JSSP in a random order would be good idea.

But we also notice more. The median time t(med) until the runs stop improving is approximately 0s. The reason is that we only perform one single objective function evaluation per run, i.e., 1 FE. Creating, mapping, and evaluating a solution can be very fast, it can happen within milliseconds. However, we had originally planned to use up to three minutes for optimization. Hence, almost all of our time budget remains unused. At the same time, we already know that that there is a 10-20% difference between the best and the median solution quality among the 101 random solutions we created. The standard deviation sd of the solution quality also is always above 100 time units of makespan. So why don’t we try to make use of this variance and the high speed of solution creation?

3.2.3 Random Sampling Algorithm

3.2.3.1 The Algorithm

Random sampling algorithm, also called random search, repeats creating random solutions until the computational budget is exhausted [188]. In our corresponding Java implementation given in Listing 3.8, we therefore only needed to add a loop around the code from the single random sampling algorithm from Listing 3.7.

public class RandomSampling<X, Y>
    extends Metaheuristic0<X, Y> {
  public void solve(IBlackBoxProcess<X, Y> process) {
// Allocate data structure for holding 1 point from search space.
    X x = process.getSearchSpace().create();
    Random random = process.getRandom();// get random gen

    do { // Repeat until budget is exhausted.
      this.nullary.apply(x, random); // Create random point in X.
// Evaluate the point: process.evaluate applies the
// representation mapping and calls objective function. It
// remembers the best solution, so the caller can obtain it.
      process.evaluate(x);
    } while (!process.shouldTerminate()); // do until time is up
  }
}

The algorithm can be described as follows:

1. Set the best-so-far objective value z to +\infty and the best-so-far candidate solution y to NULL.
2. Create random point x' in search space \mathbb{X} by using the nullary search operator.
3. Map the point x' to a candidate solution y' by applying the representation mapping y'=\gamma(x').
4. Compute the objective value z' of y' by invoking the objective function z'=f(y').
5. If z' is better than best-so-far-objective value z, i.e., z'<z, then
   1. store z' in z and
   2. store y' in y.
6. If the termination criterion is not met, return to step 2.
7. Return the best-so-far objective value z and the best solution y to the user.

In actual program code, steps 3 to 5 can again be encapsulate by a wrapper around the objective function. This reduces a lot of potential programming mistakes and makes the code much shorter. This is what we did with the implementations of the black-box process interface IBlackBoxProcess given in Listing 3.1.

3.2.3.2 Results on the JSSP

Let us now compare the performance of this iterated random sampling with our initial method. Table 3.2 shows us that the iterated random sampling algorithm is better in virtually all relevant aspects than the single random sampling method. Its best, mean, and median result quality are significantly better. Since creating random points in the search space is so fast that we can sample many more than 101 candidate solutions, even the median and mean result quality of the rs algorithm are better than the best quality obtainable with 1rs. Matter of fact, each run of our rs algorithm can create and test several million candidate solutions within the three minute time window, i.e., perform several million FEs. By remembering the best of millions of created solutions, what we effectively do is we exploit the variance of the quality of the random solutions. As a result, the standard deviation of the results becomes lower. This means that this algorithm has a more reliable performance, we are more likely to get results close to the mean or median performance when we use rs compared to 1rs.

Actually, if we would have infinite time for each run (instead of three minutes), then each run would eventually guess the optimal solutions. The variance would become zero and the mean, median, and best solution would all converge to this global optimum. Alas, we only have three minutes, so we are still far from this goal.

In Figure 3.3, we now again plot the solutions of median quality, i.e., those which are in the middle of the results, quality-wise. The improved performance becomes visible when comparing Figure 3.3 with Figure 3.2. The spacing between the jobs on the machines has significantly reduced. Still, the schedules clearly have a lot of unused time, visible as white space between the operations on the machines. We are also still relatively far away from the lower bounds of the objective function, so there is lots of room for improvement.

3.2.3.3 Progress over Time and the Law of Diminishing Returns

Another new feature of our rs algorithm is that it truly is an Anytime Algorithm (Section 3.1.1). It begins with an entirely random solution and tries to find better solutions as time goes by. Let us take a look at Figure 3.4, which illustrates how the solution quality of the runs improves over time. At first glance, this figure looks quite nice. For each of the four problem instances we investigate, our algorithm steadily and nicely improves the solution quality. Each single line (one per run) keeps slowly going down, which means that the makespan (objective value) of its best-so-far solution decreases steadily.

However, upon closer inspection, we notice that the time axes in the plots are logarithmically scaled. The first of the equally-spaces axis tick marks is at 1s, the second one at 10s, the third one at 100s, and so on. The progress curves plotted over these logarithmically scaled axes seem to progress more or less like straight linear lines, maybe even slower. A linear progress over a logarithmic time scale could mean, for instance, that we may make the same improvements in the time intervals 1s\ldots 9s, 10s\ldots 99s, 100s\ldots 999s, and so on. In other words, the algorithm improves the solution quality slower and slower.

This is the first time we witness a manifestation of a very basic law in optimization. When trying to solve a problem, we need to invest resources, be it software development effort, research effort, computational budget, or expenditure for hardware, etc. If you invest a certain amount a of one of these resources, you may be lucky to improve the solution quality that you can get by, say, b units. Investing 2a of the resources, however, will rarely lead to an improvement by 2b units. Instead, the improvements will become smaller and smaller the more you invest. This is exactly the Law of Diminishing Returns [177] known from the field of economics.

And this makes a lot of sense here. On one hand, the maximal possible improvement of the solution quality is bounded by the global optimum – once we have obtained it, we cannot improve the quality further, even if we invest infinitely much of an resource. On the other hand, in most practical problems, the amount of solutions that have a certain quality gets the smaller the closer said quality is to the optimal one. This is actually what we see in Figure 3.4: The chance of randomly guessing a solution of quality F becomes the smaller the better (smaller) F is.

From the diagrams we can also see that random sampling is not a good method to solve the JSSP. It will not matter very much if we have three minutes, six minutes, or one hour. In the end, the improvements we would get by investing more time would probably become smaller and the amount of time we need to invest to get any improvement would keep to increase. The fact that random sampling can be parallelized perfectly does not help much here, as we would need to provide an exponentially increasing number of processors to keep improving the solution quality.

3.2.4 Summary

With random sampling, we now have a very primitive way to tackle optimization problems. In each step, the algorithm generates a new, entirely random candidate solution. It remembers the best solution that it encounters and, after its computational budget is exhausted, returns it to the user.

Obviously, this algorithm cannot be very efficient. But we already also learned one method to improve the result quality and reliability of optimization methods: restarts. Restarting an optimization can be beneficial if the following conditions are met:

1. If we have a budget limitation, then most of the improvements made by the algorithm must happen early during the run. If the algorithm already uses its budget well can keep improving even close to its end, then it makes no sense to stop and restart. The budget must be large enough so that multiple runs of the algorithm can complete or at least deliver reasonable results.
2. Different runs of the algorithm must generate results of different quality. A restarted algorithm is still the same algorithm. It just exploits this variance, i.e., we will get something close to the best result of multiple runs. If the different runs deliver bad results anyway, doing multiple runs will not solve the problem.

Above we said that random sampling is not a very efficient algorithm. This is true in most reasonable scenarios. In problems where information about existing good solutions does not help us in any way to find new good solutions, we cannot really do better than random sampling. In most reasonable problems that one may try to solve, however, such information is helpful. Random sampling then is also a basic yardstick: An optimization algorithm that does not significantly outperform random sampling is useless.

3.3 Hill Climbing

Our first algorithm, random sampling, is not very efficient. It does not make any use of the information it sees during the optimization process. Each search step consists of creating an entirely new, entirely random candidate solution. Each search step is thus independent of all prior steps. If the problem that we try to solve is entirely without structure, then this is already the best we can do. But our JSSP problem is not without structure. For example, we can assume that if we swap the last two jobs in a schedule, the makespan of the resulting new schedule should still be somewhat similar to the one of the original plan. Actually, most reasonable optimization problems have a lot of structure. We therefore should try to somehow make use of the information gained from sampling candidate solutions.

Local search algorithms [112,205] offer one idea for how to do that. They remember one point x^m in the search space \mathbb{X}. In every step, a local search algorithm investigates g\geq1 points x_i which are derived from and are similar to x^m. From the joined set of these g points and x^m, only one point is chosen as the (new or old) x^m. All the other points are discarded.

Definition 21. Causality means that small changes in the features of an object (or candidate solution) also lead to small changes in its behavior (or objective value).

Local search exploits a property of many optimization problems which is called causality [168,169,208,216]. If we have two points in the search space that only differ a little bit, then they likely map to similar schedules, which, in turn, likely have similar makespans. This means that if we have a good candidate solution, then there may exist similar solutions which are better. We hope to find one of them and then continue trying to do the same from there.

3.3.1 Ingredient: Unary Search Operation for the JSSP

For now, let us limit ourselves to local searches creating g=1 new point in each iteration. The question arises how we can create a candidate solution which is similar to, but also slightly different from, one that we already have? Our search algorithms are working in the search space \mathbb{X}. So we need one operation which accepts an existing point x\in\mathbb{X} and produces a slightly modified copy of it as result. In other words, we need to implement a unary search operator!

On a JSSP with m machines and n jobs, our representation \mathbb{X} encodes a schedule as an integer array of length m*n containing each of the job IDs (from 0\ldots(n-1)) exactly m times. The sequence in which these job IDs occur then defines the order in which the jobs are assigned to the machines, which is realized by the representation mapping \gamma (see Listing 2.7).

One idea to create a slightly modified copy of such a point x in the search space would be to simply swap two of the jobs in it. Such a 1swap operator can be implemented as follows:

1. Make a copy x' of the input point x from the search space.
2. Pick a random index i from 0\ldots(m*n-1).
3. Pick a random index j from 0\ldots(m*n-1).
4. If the values at indexes i and j in x' are the same, then go back to point 3.
5. Swap the values at indexes i and j in x'.
6. Return the now modified copy x' of x.

Point 4 is important since swapping the same values makes no sense, as we would then get x'=x. Then, also the mappings \gamma(x) and \gamma(x') would be the same, i.e., we would actually not make a move and just waste time.

public class JSSPUnaryOperator1Swap
    implements IUnarySearchOperator<int[]> {
  public void apply(int[] x, int[] dest,
      Random random) {
// copy the source point in search space to the dest
    System.arraycopy(x, 0, dest, 0, x.length);

// choose the index of the first operation to swap
    int i = random.nextInt(dest.length);
    int jobI = dest[i]; // remember job id

    for (;;) { // try to find a location j with a different job
      int j = random.nextInt(dest.length);
      int jobJ = dest[j];
      if (jobI != jobJ) { // we found two locations with two
        dest[i] = jobJ; // different values
        dest[j] = jobI; // then we swap the values
        return; // and are done
      }
    }
  }
}

We implemented this operator in Listing 3.9. Notice that the operator is randomized, i.e., applying it twice to the same point in the search space will likely yield different results.

In Figure 3.5, we illustrate the application of this operator to one point x in the search space for our demo JSSP instance. It swaps the two jobs at index i=10 and j=15 of x. In the new, modified copy x', the jobs 3 and 0 at these indices have thus traded places. The impact of this modification becomes visible when we map both x and x' to the solution space using the representation mapping \gamma. The 3 which has been moved forward now means that job 3 will be scheduled before job 1 on machine 2. As a result, the last two operations of job 3 can now finish earlier on machines 0 and 1, respectively. However, time is wasted on machine 2, as we first need to wait for the first two operations of job 3 to finish before we can execute it there. Also, job 1 finishes now later on that machine, which also delays its last operation to be executed on machine 4. This pushes back the last operation of job 0 (on machine 4) as well. The new candidate solution \gamma(x') thus has a longer makespan of f(\gamma(x'))=195 compared to the original solution with f(\gamma(x))=180.

In other words, our application of 1swap in Figure 3.5 has led us to a worse solution. This will happen most of the time. As soon as we have a good solution, the solutions similar to it tend to be worse in average and the number of even better solutions in the neighborhood tends to get smaller. However, if we would have been at x' instead, an application of 1swap could well have resulted in x. In summary, we can hope that the chance to find a really good solution by iteratively sampling the neighborhoods of good solutions is higher compared to trying to randomly guessing them (as rs does)   even if most of our samples are worse.

3.3.2 Stochastic Hill Climbing Algorithm

3.3.2.1 The Algorithm

Stochastic Hill Climbing [173,187,205] is the simplest implementation of local search. It is also sometimes called localized random search [188]. It proceeds as follows:

1. Create one random point x in the search space \mathbb{X} using the nullary search operator.
2. Map the point x to a candidate solution y by applying the representation mapping y=\gamma(x).
3. Compute the objective value by invoking the objective function z=f(y).
4. Repeat until the termination criterion is met:
   1. Apply the unary search operator to x to get a slightly modified copy x' of it.
   2. Map the point x' to a candidate solution y' by applying the representation mapping y'=\gamma(x').
   3. Compute the objective value z' by invoking the objective function z'=f(y').
   4. If z'<z, then store x' in x, store y' in y, and store z' in z.
5. Return the best encountered objective value z and the best encountered solution y to the user.

This algorithm is implemented in Listing 3.10 and we will refer to it as hc.

public class HillClimber<X, Y>
    extends Metaheuristic1<X, Y> {
  public void solve(IBlackBoxProcess<X, Y> process) {
// initialize local variables xCur, xBest, random
// Create starting point: a random point in the search space.
    this.nullary.apply(xBest, random); // xBest = random point
    double fBest = process.evaluate(xBest); // map & evaluate

    while (!process.shouldTerminate()) {
// Create a slightly modified copy of xBest and store in xCur.
      this.unary.apply(xBest, xCur, random);
// Map xCur from X to Y and evaluate candidate solution.
      double fCur = process.evaluate(xCur);
      if (fCur < fBest) { // we found a better solution
// Remember best objective value and copy xCur to xBest.
        fBest = fCur;
        process.getSearchSpace().copy(xCur, xBest);
      } // Otherwise, i.e., fCur >= fBest: Just forget xCur.
    } // Repeat until computational budget is exhausted.
  } // `process` has remembered the best candidate solution.
}

If you are wondering what would happen if we would accept the new solution x' also if z'=z, i.e., replace the z'<z with an z'\leq z in point 4.d of the algorithm definition: This algorithm is called (1+1) EA and discussed later in Section 3.4.6.1.

3.3.2.2 Results on the JSSP

We now plug our unary operator 1swap into our hc algorithm and apply it to the JSSP. We will refer to this setup as hc_1swap and present its results with those of rs in Table 3.3.

The hill climber outperforms random sampling in almost all aspects. It produces better mean, median, and best solutions. Actually, its median and mean solutions are better than the best solutions discovered by rs. Furthermore, it finds its solutions much much faster. The median time med(t) consumed until the algorithm converges is not more than one seconds. The median number of consumed FEs med(FEs) to find the best solutions per run is between 7000 and 105’000, i.e., between one 50^th and one 2500^th of the number of FEs needed by rs.

It may be interesting to know that this simple hc_1swap algorithm can already achieve some remotely acceptable performance, even though it is very far from being useful. For instance, on instance abz7, it delivers better best and mean results than all four Genetic Algorithms (GAs) presented in [122]. On la24, only one of the four (GA‑PR) has a better best result and all lose in terms of mean result. On this instance, hc_1swap finds a better best solution than all six GAs in [2] and better mean results than five of them. In Section 3.4, we will later introduce (better-performing!) Evolutionary Algorithms, to which GAs belong.

The Gantt charts of the median solutions of hc_1swap are illustrated in Figure 3.6. They are more compact than those discovered by rs and illustrated in Figure 3.3.

Figure 3.7 shows how both hc_1swap and rs progress over time. In Figure 3.4, we plotted every individual run. This time, we plot the median (see Section 4.4.2) of achieved quality at each time step as thick line. In the background, we plot the whole range of the values as semi-transparent region. The more runs fall into one region, the stronger we plot the color: In the outer borders of the lightest color shade mark the range between the best and worst results of any run at a given time. The next stronger-shaded region contains about 95% of the runs. Then follow by 68% of the runs while the strongest-shaded region holds half of the runs.

It should be noted that I designed the experiments in such a way that there were 101 different random seeds per instance. For each instance, all algorithms use the same random seeds, i.e., the hill climber and random sampling start with the same initial solutions. Still, the shaded regions of the two algorithms separate almost immediately.

We already knew from Table 3.3 that hc_1swap converges very quickly. After initial phases with quick progress, it stops making any further progress, usually before 1000 milliseconds have been consumed. This fits well to the values med(t) given in Table 3.3. With the exception of instance la24, where two runs of the hill climber performed exceptionally bad, there is much space between the runs of rs and hc_1swap. We can also see again that there is more variance in the end results of hc_1swap compared to those of rs, as they are spread wider in the vertical direction.

3.3.3 Stochastic Hill Climbing with Restarts

Upon close inspection of the results, we notice that we are again in the same situation as with the 1rs algorithm: There is some variance between the results and most of the action takes place in a short time compared to our total computational budget (1 second vs. 3 minutes). Back in Section 3.2.3 we made use of this situation by simply repeating 1rs until the computational budget was exhausted, which we called the rs algorithm. Now the situation is a bit different, however. 1rs creates exactly one solution and is finished, whereas our hill climber does not actually finish. It keeps creating modified copies of the current solution, only that these eventually do not mark improvements anymore. Then, the algorithm has converged into a local optimum.

Definition 22. A local optimum is a point x^{\times} in the search space which maps to a better candidate solution than any other points in its neighborhood (see Definition 17).

Definition 23. An optimization process has prematurely converged if it has not yet discovered the global optimum but can no longer improve its approximation quality. [208,216]

Due to the black-box nature of our basic hill climber algorithm, it is not really possible to know when the complete neighborhood of the current solution has already been tested. We thus cannot know whether or not the algorithm is trapped in a local optimum and has prematurely converged. However, we can try to guess it: If there has not been any improvement for a high number L of steps, then the current point x in the search space is probably a local optimum. If that happens, we just restart at a new random point in the search space. Of course, we will remember the best-so-far candidate solution in a special variable y_{b} over all restarts and return it to the user in the end.

3.3.3.1 The Algorithm

1. Set counter C of unsuccessful search steps to 0.
2. Set the best-so-far objective value z_{b} to +\infty and the best-so-far candidate solution y_{b} to NULL.
3. Create a random point x in the search space \mathbb{X} using the nullary search operator.
4. Map the point x to a candidate solution y by applying the representation mapping y=\gamma(x).
5. Compute the objective value by invoking the objective function z=f(y).
6. If z<z_{b}, then store z in z_{b} and store y in y_{b}.
7. Repeat until the termination criterion is met:

   1. Apply the unary search operator to x to get the slightly modified copy x' of it.

   2. Map the point x' to a candidate solution y' by applying the representation mapping y'=\gamma(x').

   3. Compute the objective value z' by invoking the objective function z'=f(y').

   4. If z'<z, then

      1. store z' in z and store x' in x and
      2. set C to 0.
      3. If z'<z_{b}, then store z' in z_{b} and store y' in y_{b}.

      otherwise to step 7d, i.e., if z'\geq z, then

      4. increment C by 1.
      5. If C\geq L, then perform a restart by going back to step 3.
8. Return best encountered objective value z_{b} and the best encountered solution y_{b} to the user.

public class HillClimberWithRestarts<X, Y>
    extends Metaheuristic1<X, Y> {
  public void solve(IBlackBoxProcess<X, Y> process) {
// omitted for brevity initialize local variables xCur, xBest,
// random, failsBeforeRestart, and failCounter=0
    while (!(process.shouldTerminate())) { // outer loop: restart
      this.nullary.apply(xBest, random); // start=random solution
      double fBest = process.evaluate(xBest); // evaluate it
      long failCounter = 0L; // initialize counters

      while (!(process.shouldTerminate())) { // inner loop
        this.unary.apply(xBest, xCur, random); // try to improve
        double fCur = process.evaluate(xCur); // evaluate

        if (fCur < fBest) { // we found a better solution
          fBest = fCur; // remember best quality
          process.getSearchSpace().copy(xCur, xBest); // copy
          failCounter = 0L; // reset number of unsuccessful steps
        } else { // ok, we did not find an improvement
          if ((++failCounter) >= this.failsBeforeRestart) {
            break; // jump back to outer loop for restart
          } // increase fail counter
        } // failure
      } // inner loop
    } // outer loop
  } // process has stored best-so-far result
}

Now this algorithm – implemented in Listing 3.11 – is a bit more elaborate. Basically, we embed the original hill climber into a loop. This inner hill climber will stop after a certain number L of unsuccessful search steps, which then leads to a new round in the outer loop. In combination with the 1swap operator, we refer to this algorithm as hcr_L_1swap, where L is to be replaced with the actual value of the parameter L.

3.3.3.2 The Right Setup

We now realize that we do not know which value of L is good. If we pick it too low, then the algorithm will restart before it actually converges to a local optimum, i.e., stop while it could still improve. If we pick it too high, we waste runtime and do fewer restarts than what we could do.

If we do not know which value for a parameter is reasonable, we can always do an experiment to investigate. Since the order of magnitude of the proper value for L is not yet clear, it makes sense to test exponentially increasing numbers. Here, we test the powers of two from 2^7=128 to 2^{18}=262'144. For each value, we plot the scaled median result quality over the 101 runs in Figure 3.8. In this diagram, the horizontal axis is logarithmically scaled.

From the plot, we can confirm our expectations: Small numbers of L perform bad and high numbers of L cannot really improve above the basic hill climber. Actually, if we would set L to a number larger than the overall budget, then we would obtain exactly the original hill climber, as it would never perform any restart. For different problem instances, different values of L perform good, but L\approx2^{14}=16'384 seems to be a reasonable choice for three of the four instances.

3.3.3.3 Results on the JSSP

The performance indicators of three settings of our hill climber with restarts in comparison with the plain hill climber are listed in Table 3.4. We know that L=2^{14} seems a reasonable setting. Additionally, we also list the adjacent setups, i.e., give the results for L\in\{2^{13},2^{14},2^{15}\}.

Table 3.4 shows us that the restarted algorithms hcr_L_1swap almost always provide better best, mean, and median solutions than hc_1swap. Only the overall best result of hcr_8192_1swap on abz7 is worse than for hc_1swap – on all other instances and for all other quality metrics, hc_1swap loses.

The standard deviations of the end results of the variants with restarts are also always smaller, meaning that these algorithms perform more reliably. Their median time until they converge is now higher, which means that we make better use of our computational budget.

As a side note, the median and mean result of the three listed setups of our very basic hcr_L_1swap algorithms for instance la24 are already better than the best result (982) delivered by the Grey Wolf Optimization algorithm proposed in [118]. In other words, even with such a simple algorithm we can already achieve results not very far from recent publications…

The median solutions discovered by hcr_16384_1swap, illustrated in Figure 3.9, again show less wasted time. The scheduled jobs again move a bit closer together.

From the progress diagrams plotted in Figure 3.10, we can see that the algorithm version with restart initially behave very similar to the original hill climber. Its median quality is almost exactly the same during the first 100 ms. This makes sense, because until its first restart, hcr_16384_1swap is identical to hc_1swap. However, when hc_1swap has converged and stops making improvements, hcr_16384_1swap still continues to make progress.

Of course, our way of finding the right value for the restart parameter L was rather crude. Most likely, the right choice could be determined per instance based on the number m of machines and the number n of jobs. But even with such a coarse way of algorithm configuration, we managed to get rather good results.

3.3.3.4 Drawbacks of the Idea of Restarts

With our restart method, we could significantly improve the results of the hill climber. It directly addressed the problem of premature convergence, but it tried to find a remedy for its symptoms, not for its cause.

Basically, a restarted algorithm is still the same algorithm – we just restart it again and again. If there are many more local optima than global optima, every restart will probably end again in a local optimum. If there are many more bad'' local optima thangood’’ local optima, every restart will probably end in a ``bad’’ local optimum. While restarts improve the chance to find better solutions, they cannot solve the intrinsic shortcomings of an algorithm.

Another problem is: With every restart we throw away all accumulated knowledge and information of the current run. Restarts are therefore also somewhat wasteful.

3.3.4 Hill Climbing with a Different Unary Operator

3.3.4.1 Small vs. Large Neighborhoods – and Uniform vs. Non-Uniform Sampling

One of issues limiting the performance of our restarted hill climber is the design of the unary operator. 1swap will swap two jobs in an encoded solution. Since the solutions are encoded as integer arrays of length m*n, there are m*n possible choices when picking the index of the first job to be swapped. We swap only with different jobs and each job appears m times in the encoding. This leaves m*(n-1) choices for the second swap index, because we will only use a second index that points to a different job ID. If we think about the size of the neighborhood spanned by 1swap, we can also ignore equivalent swaps: Exchanging the jobs at indexes (10,5) and (5,10), for example, would result in the same outcome. In total, from any given point in the search space, 1swap may reach 0.5*(m*n)*[m*(n-1)]=0.5*m^2(n^2-n) different other points. Some of these points might still actually encode the same candidate solutions, i.e., identical schedules. In other words, the neighborhood spanned by our 1swap operator equals only a tiny fraction of the big search space (remember Table 2.3).

This has two implications:

1. The chance of premature convergence for a hill climber applying this operator is relatively high, since the neighborhoods are relatively small. If the neighborhood spanned by the operator was larger, it would contain more, potentially better solutions. Hence, it would take longer for the optimization process to reach a point where no improving move can be discovered anymore.
2. Assume that there is no better solution in the 1swap neighborhood of the current best point in the search space. There might still be a much better, similar solution. Finding it could, for instance, require swapping three or four jobs – but the hc_1swap algorithm will never find it, because it can only swap two jobs. If the search operator would permit such moves, then even the plain hill climber may discover this better solution.

So let us try to think about how we could define a new unary operator which can access a larger neighborhood. As we should always do, we first consider the extreme cases.

On the one hand, we have 1swap which samples from a relatively small neighborhood. Because the neighborhood is small, the stochastic hill climber will eventually have visited all of the solutions it contains. If none of them is better than the current best solution, it will not be able to depart from it.

The other extreme could be to use our nullary operator as unary operator: It would return an entirely random point from the search space \mathbb{X} and ignore its input. Then, each point x\in\mathbb{X} would have the whole \mathbb{X} as the neighborhood. Using such a unary operator would turn the hill climber into random sampling (and we do not want that).

From this thought experiment we know that unary operators which indiscriminately sample from very large neighborhoods are probably not very good ideas, as they are too random. They also make less use of the causality of the search space, as they make large steps and their produced outputs are very different from their inputs.

Using an operator that creates larger neighborhoods than 1swap, which are still smaller than \mathbb{X} would be one idea. For example, we could always swap three jobs instead of two. The more jobs we swap in each application, the larger the neighborhood gets. Then we will be less likely to get trapped in local optima (as there will be fewer local optima). But we will also make less and less use of the causality property, i.e., the solutions we derive from the current best one will be more and more different from it. Where should we draw the line? How many jobs should we swap?

Well, there is one more aspect of the operators that we did not think about yet. An operator does not just span a neighborhood, but it also defines a probability distribution over it. So far, our 1swap unary operator samples uniformly from the neighborhood of its input. In other words, all of the 0.5*m^2(n^2-n) new points that it could create in one step have exactly the same probability, the same chance to be chosen.

But we do not need to do it like that. We could construct an operator that often creates outputs very similar to its input (like 1swap), but also, from time to time, samples points a bit farther away in the search space. This operator could have a huge neighborhood – but sample it non-uniformly.

3.3.4.2 Second Unary Search Operator for the JSSP

We define the nswap operator for the JSSP as follows and implement it in Listing 3.12:

public class JSSPUnaryOperatorNSwap
    implements IUnarySearchOperator<int[]> {
  public void apply(int[] x, int[] dest,
      Random random) {
// copy the source point in search space to the dest
    System.arraycopy(x, 0, dest, 0, x.length);

// choose the index of the first operation to swap
    int i = random.nextInt(dest.length);
    int first = dest[i];
    int last = first; // last stores the job id to "swap in"

    boolean hasNext;
    do { // we repeat a geometrically distributed number of times
      hasNext = random.nextBoolean();
      inner: for (;;) { // find a location with a different job
        int j = random.nextInt(dest.length);
        int jobJ = dest[j];
        if ((last != jobJ) && // don't swap job with itself
            (hasNext || (first != jobJ))) { // also not at end
          dest[i] = jobJ; // overwrite job at index i with jobJ
          i = j; // remember index j: we will overwrite it next
          last = jobJ; // but not with the same value jobJ...
          break inner;
        }
      }
    } while (hasNext); // Bernoulli process

    dest[i] = first; // write back first id to last copied index
  }
}

1. Make a copy x' of the input point x from the search space.
2. Pick a random index i from 0\ldots(m*n-1).
3. Store the job-id at index i in the variable f for holding the very first job, i.e., set f=x'_{i}.
4. Set the job-id variable l for holding the last-swapped-job to x'_{i} as well.
5. Repeat
   1. Decide whether we should continue the loop after the current iteration (TRUE) or not (FALSE) with equal probability and remember this decision in variable n.
   2. Pick a random index j from 0\ldots(m*n-1).
   3. If l=x'_{j}, go back to point b.
   4. If f=x_{j} and we will not do another iteration (n=FALSE), go back to point b.
   5. Store the job-id at index j in the variable l.
   6. Copy the job-id at index j to index i, i.e., set x'_{i}=x'_{j}.
   7. Set i=j.
6. If we should do another iteration (n=TRUE), go back to point 5.
7. Store the first-swapped job-id f in x'_{i}.
8. Return the now modified copy x' of x.

Here, the idea is that we will perform at least one iteration of the loop (point 5). If we would do exactly one iteration, then we would pick two indices i and j where different job-ids are stored, as l must be different from f (point c and d). We would then swap the jobs at these indices (points f, g, and 7). In the case of exactly one iteration of the main loop, this operator behaves the same as 1swap. This takes place with a probability of 0.5 (point a).

If we do two iterations, i.e., pick true the first time we arrive at point a and false the second time, then we swap three job-ids instead. Let us say we picked indices \alpha at point 2, \beta at point b, and \gamma when arriving the second time at b. We will store the job-id originally stored at index \beta at index \alpha, the job originally stored at index \gamma at index \beta, and the job-id from index \gamma to index \alpha. Condition c prevents index \beta from referencing the same job-id as index \alpha and index \gamma from referencing the same job-id as what was originally stored at index \beta. Condition d only applies in the last iteration and prevents \gamma from referencing the original job-id at \alpha.

This three-job swap will take place with probability 0.5*0.5=0.25. Similarly, a four-job-swap will happen with half of that probability, and so on. In other words, we have something like a Bernoulli process, where we decide whether or not to do another iteration by flipping a fair coin, where each choice has probability 0.5. The number of iterations will therefore be geometrically distributed with an expectation of two job swaps. Of course, we only have m different job-ids in a finite-length array x', so this is only an approximation. Generally, this operator will most often apply small changes and sometimes bigger steps. The bigger the search step, the less likely will it be produced. The operator therefore can make use of the causality while – at least theoretically – being able to escape from any local optimum.

3.3.4.3 Results on the JSSP

Let us now compare the end results that our hill climbers can achieve using either the 1swap or the new nswap operator after three minutes of runtime on my laptop computer in Table 3.5.

From Table 3.5, we find that hc_nswap performs almost always better than hc_1swap. Only on instance abz7, hc_1swap finds the better best solution. For all other instances, hc_nswap has better best, mean, and median results. It also converges much later and often performs 7 to 15 million function evaluations and consumes 14% to 25% of the three minute budget before it cannot improve anymore. Still, the hill climber hcr_16384_1swap using the 1swap operator with restarts tends to outperform hc_nswap.

Figure 3.11 illustrates the progress of the hill climbers with the 1swap and nswap operators. Initially, both algorithms behave very similar in median. This may be because half of the time, hc_nswap also performs single-swap moves. However, while hc_1swap stops improving even before one second has elapsed, hc_nswap can still improve after ten seconds. Still, these late improvements tend to be small and occur infrequently. It may be that the algorithm arrives in local optima from which it can only escape with complicated multi-swap moves, which are harder to discover by chance.

3.3.5 Combining Bigger Neighborhood with Restarts

Both restarts and the idea of allowing bigger search steps with small probability are intended to decrease the chance of premature convergence, while the latter one also can investigate more solutions similar to the current best one. We have seen that both measures work separately. The fact that hc_nswap improves more and more slowly towards the end of the computational budget means that it could be interesting to try to combine both ideas, restarts and larger neighborhoods.

We plug the nswap operator into the hill climber with restarts and obtain algorithm hcr_L_nswap. We perform the same experiment to find the right setting for the restart limit L as for the hcr_L_1swap algorithm and illustrate the results in Figure 3.12.

The sweet spot for the number of unsuccessful FEs before a restart has increased compared to before. This makes sense, because we already know that nswap can keep improving longer.

3.3.5.1 Results on the JSSP

From Table 3.6, where we print the results of hcr_32768_nswap and hcr_65536_nswap, we can find that the algorithm version with restarts performs better in average than the one without. However, it does not always find the best solution, as can be seen on instance swv15, where hc_nswap finds a schedule of length 3602. The differences between hcr_16384_1swap and hcr_L_nswap, however, are quite small. If we compare the progress over time of hcr_16384_1swap and hcr_65536_nswap, then the latter seems to have a slight edge over the former – but only by about half of a percent. This small difference is almost indistinguishable in the progress diagram Figure 3.13.

3.3.5.2 Testing for Significance

Still, even small improvements can have a big economical impact. Saving 0.5% of 10’000’000 RMB is still 50’000 RMB. The problem is knowing whether such small improvements are true improvements or artifacts of the randomness in the algorithm.

In order to understand the latter situation, consider the following thought experiment. Assume you have a completely unbiased, uniform source of true random real numbers from the interval [0,1). You draw 500 such numbers, i.e., have a list A containing 500 numbers, each from [0,1). Now you repeat the experiment and get another such list B. Since the numbers stem from a random source, we can expect that A\neq B. If we compute the medians A and B, they are likely to be different as well. Actually, I just did exactly this in the R programming language and got median(A)=0.5101432 and median(B)=0.5329007. Does this mean that the generator producing the numbers in A creates somehow smaller numbers than the generator from which the numbers in B stem? Obviously not, because we sampled the numbers from the same source. Also, every time I would repeat this experiment, I would get different results.

So how do we know whether or not the sources of A and B are truly different? Well, we cannot really know for sure. But we can we can make a statement which is wrong with at most a given probability. This is called statistical testing, and we discuss it in detail in Section 4.5. Thus, in order to see whether the observed small performance difference of the hcr setups is indeed real or just random jitter, we compare their sets of 101 end results on each of the problem instances. For this purpose, we use the Mann-Whitney U test, as prescribed in Section 4.5.4.

From Table 3.7 we know that if we would claim hcr_16384_1swap tends to produce results with large makspan than hcr_65536_nswap on abz7, then, from our experimental data, we can estimate our chance to be wrong to be about 30%. In other words, making that claim would be quite a gamble and we can conclude that here, the differences we observed in the experiment are not statistically significant (marked with ? in the table). However, if we would claim the same for swv15, our chance to be wrong is about 1\!\cdot\!10^{-10}, i.e., very small. So on swv15, we find that hcr_65536_nswap very likely performs better.

In summary, Table 3.7 tells us that it seems that the hill climber with restarts using the nswap operator can outperform the one using the 1swap operator on swv15 and yn4. For abz7, there certainly is no significant difference. The p-value of 3.57\!\cdot\!10^{-3} on la24 is fairly small but still above the Bonferroni-corrected \alpha'=7.14\!\cdot\!10^{-4}, so we cannot make any statement here.

We found that the nswap operator tends to be better then the 1swap operator when plugged into the hill climber. We found that restarts are also beneficial to improve the performance of our hill climber. While we did not find that combining these two methods will yield a convincing improvement, we can at least conclude that preferring hcr_65536_nswap over hcr_16384_1swap seems reasonable. It might sometimes give better performance while it is unlikely to perform worse.

3.3.6 Summary

In this section, we have learned about our first reasonable optimization method. The stochastic hill climbing algorithm always remembers the best-so-far point in the search space. In each step, it applies the unary operator to obtain a similar but slightly different point. If it is better, then it becomes the new best-so-far point. Otherwise, it is forgotten.

The performance of hill climbing depends very much on the unary search operator. If the operator samples from a very small neighborhood only, like our 1swap operator does, then the hill climber might quickly get trapped in a local optimum. A local optimum here is a point in the search space which is surrounded by a neighborhood that does not contain any better solution. If this is the case, the two conditions for doing efficient restarts may be fulfilled: quick convergence and variance of result quality.

The question when to restart then arises, as we usually cannot find out if we are actually trapped in a local optimum or whether the improving move (application of the unary operator) just has not been discovered yet. The most primitive solution is to simply set a limit L for the maximum number of moves without improvement that are permitted.

Our hcr_L_1swap was born. We configured L in a small experiment and found that L=16384 seemed to be reasonable. The setup hcr_16384_1swap performed much better than hc_1swap. It should be noted that our experiment used for configuration was not very thorough, but it should suffice at this stage. We can also note that it showed that different settings of L are better for different instances. This is probably related to the corresponding search space size – but we will not investigate this any further here.

A second idea to improve the hill climber was to use a unary operator spanning a larger neighborhood, but which still most often sampled solutions similar to current one. The nswap operator gave better results than than the 1swap operator in the basic hill climber. The take-away message is that different search operators may (well, obviously) deliver different performance and thus, testing some different operators can always be a good idea.

Finally, we tried to combine our two improvements, restarts and better operator, into the hcr_L_nswap algorithm. Here we learned the lesson that performance improvements do not necessarily add up. If we have a method that can deliver an improvement of 10% of solution quality and combine it with another one delivering 15%, we may not get an overall 25% improvement. Indeed, our hcr_65536_nswap algorithm only performed a little bit better than hcr_16384_1swap.

From this chapter, we also learned one more lesson: Many optimization algorithms have parameters. Our hill climber had two: the unary operator and the restart limit L. Configuring these parameters well can lead to significant improvements.

3.4 Evolutionary Algorithms

We now already have one functional, basic optimization method – the hill climber. Different from the random sampling approach, it makes use of some knowledge gathered during the optimization process, namely the best-so-far point in the search space. However, only using this single point led to the danger of premature convergence, which we tried to battle with two approaches, namely restarts and the search operator nswap, spanning a larger neighborhood from which we sampled in a non-uniform way. These concepts can be transfered rather easily to many different kinds of optimization problems. Now we will look at a third concept to prevent premature convergence: Instead of just remembering and utilizing only one single point from the search space in each iteration of our algorithm, we will now work on an array of points!

3.4.1 Evolutionary Algorithm without Recombination

Today, there exists a wide variety of different Evolutionary Algorithms (EAs) [17,41,54,84,145,146,188,205]. We will here discuss a simple yet efficient variant: the (\mu+\lambda) EA without recombination.⁴ This algorithm always remembers the best \mu\in\mathbb{N}_1 points in the search space found so far. In each step, it derives \lambda\in\mathbb{N}_1 new points from them by applying the unary search operator.

3.4.1.1 The Algorithm (without Recombination)

The basic (\mu+\lambda) Evolutionary Algorithm works as follows:

1. I\in\mathbb{X}\times\mathbb{R} be a data structure that can store one point x in the search space and one objective value z.
2. Allocate an array P of length \mu+\lambda of instances of I.
3. For index i ranging from 0 to \mu+\lambda-1 do
   1. Create a random point from the search space using the nullary search operator and store it in P_{i}.x.
   2. Apply the representation mapping y=\gamma(P_{i}.x) to get the corresponding candidate solution y.
   3. Compute the objective objective value of y and store it at index i as well, i.e., P_{i}.z=f(y).
4. Repeat until the termination criterion is met:
   4. Sort the array P in ascending order according to the objective values, i.e., such that the records r with better associated objective value r.z are located at smaller indices.
   5. Shuffle the first \mu elements of P randomly.
   6. Set the source index p=-1.
   7. For index i ranging from \mu to \mu+\lambda-1 do
      1. Set the source index p to p=(p+1)\textrm{~mod~}\mu.
      2. Apply unary search operator to the point stored at index p and store result at index i, i.e., set P_{i}.x=\textrm{searchOp}_1(P_{p}.x).
      3. Apply the representation mapping y=\gamma(P_{i}.x) to get the corresponding candidate solution y.
      4. Compute the objective objective value of y and store it at index i as well, i.e., P_{i}.z=f(y).
5. Return the candidate solution corresponding to the best record in P (i.e., the best-ever encountered solution) to the user.

public class EA<X, Y> extends Metaheuristic2<X, Y> {
  public void solve(IBlackBoxProcess<X, Y> process) {
// omitted: initialize local variables random, searchSpace, and
// the array P of length mu+lambda
// first generation: fill population with random solutions
    for (int i = P.length; (--i) >= 0;) {
      X x = searchSpace.create(); // allocate point
      this.nullary.apply(x, random); // fill with random data
      P[i] = new Record<>(x, process.evaluate(x)); // evaluate
    }

    for (;;) { // main loop: one iteration = one generation
// sort the population: mu best records at front are selected
      Arrays.sort(P, Record.BY_QUALITY);
// shuffle the first mu solutions to ensure fairness
      RandomUtils.shuffle(random, P, 0, this.mu);
      int p1 = -1; // index to iterate over first parent

// overwrite the worse lambda solutions with new offsprings
      for (int index = P.length; (--index) >= this.mu;) {
        if (process.shouldTerminate()) { // we return
          return; // best solution is stored in process
        }

        Record<X> dest = P[index];
        p1 = (p1 + 1) % this.mu; // step the parent 1 index
        Record<X> sel = P[p1];
// create modified copy of parent using unary operator
          this.unary.apply(sel.x, dest.x, random);
// map to solution/schedule and evaluate quality
        dest.quality = process.evaluate(dest.x);
      } // the end of the offspring generation
    } // the end of the main loop
  }
}

This algorithm is implemented in Listing 3.13. There, we make use of instances of the utility class Record<X>, which holds one point x in the search space along with their corresponding objective values stored in the field quality.

Basically, the algorithm starts out by creating and evaluating \mu+\lambda random candidate solutions (point 3).

Definition 24. Each iteration of the main loop of an Evolutionary Algorithm is called a generation.

Definition 25. The array of solutions under investigation in an EA is called population.

In each generation, the \mu best points in the population P are retained and the other \lambda solutions are overwritten with newly sampled points.

Definition 26. The selection step in an Evolutionary Algorithm picks the set of points in the search space from which new points should be derived. This usually involves choosing a smaller number \mu\in\mathbb{N}_1 of points from a larger array P. [17,26,36,146,205]

Selection can be done by sorting the array P (point 4d). This way, the best \mu solutions end up at the front of the array on the indices from 0 to \mu-1. The worse \lambda solutions are at index \mu to \mu+\lambda-1. These are overwritten by sampling points from the neighborhood of the \mu selected solutions by applying the unary search operator (which, in the context of EAs, is often called mutation operator).

Definition 27. The selected points in an Evolutionary Algorithm are called parents.

Definition 28. The points in an Evolutionary Algorithm that are sampled from the neighborhood of the parents by applying search operators are called offspring.

Definition 29. The reproduction step in an Evolutionary Algorithm uses the selected \mu\in\mathbb{N}_1 points from the search space to derive \lambda\in\mathbb{N}_1 new points.

For each new point to be created during the reproduction step, we apply a search operator to one of the selected \mu points. The index p in steps 4f to 4g identifies the point to be used as source for sampling the next new solution. By incrementing p before each application of the search operator, we try to make sure that each of the selected points is used approximately equally often to create new solutions. Of course, \mu and \lambda can be different (often \lambda>\mu), so if we would just keep increasing p for \lambda times, it could exceed \mu. We thus perform a modulo division with \mu in step 4g.i, i.e., set p to the remainder of the division with \mu, which makes sure that p will be in 0\ldots(\mu-1).

If \mu\neq\lambda, then the best solutions in P tend to be used more often, since they may survive selection several times and often be at the front of P. This means that, in our algorithm, they would be used more often as input to the search operator. To make our algorithm more fair, we randomly shuffle the selected \mu points (point 4e). This does not change the fact that they have been selected.

Since our algorithm will never prefer a worse solution over a better one, it will also never lose the best-so-far solution. We therefore can simply pick the best element from the population once the algorithm has converged. It should be mentioned that there are selection methods in EAs which might reject the best-so-far solution. In this case, we would need to remember it in a special variable like we did in case of the hill climber with restarts. Here, we do not consider such methods, as we want to investigate a plain EA.

3.4.1.2 The Right Setup

After implementing the (\mu+\lambda) EA as discussed above, we already have all the ingredients ready to apply to the JSSP. We need to decide which values for \mu and \lambda we want to use. The configuration of EAs is a whole research area itself. The question arises which values for \mu and \lambda are reasonable. Without investigating whether this is the best idea, let us set \mu=\lambda here, so we only have two parameters to worry about: \mu and the unary search operator. We already have two unary search operators. Let us call our algorithms of this type ea_mu_unary, where mu will stand for the value of \mu and \lambda and unary can be either 1swap or nswap. We no therefore do a similar experiment as in Section 3.3.3.2 in order to find the right parameters.

In Figure 3.14, we illustrate this experiment. Regarding \mu and \lambda, we observe the same situation as with the restarts parameters hill climber. There is a sweet spot somewhere between small and large population sizes. For small values of \mu, the algorithm may end up in a local optimum, whereas for large values, it may not be able to perform sufficiently many generations to arrive at a good solution. Nevertheless, compared to the hill climber with restart and with the exception for instance swv15, the EA performs quite stable: There are only relatively little differences in result quality for both unary operators and over many scales of \mu. This generally a nice feature, as we would like to have algorithms that are not too sensitive to parameter settings.

The setting \mu=\lambda=16'384 seems to work well for instances abz7, la25, and yn4. Interestingly, instance swv15 behaves different: here, the setting \mu=\lambda=1024 works best. It is quite common in optimization that different problem instances may require different setups to achieve the best performance. Here we see this quite pronounced. This is generally a feature that we do not like. We would ideally like to have algorithms where a good parameter setting performs well on many instances as opposed to such where each instance requires a totally different setup. Luckily, this here only concerns one of our example problem instances.

Regarding the choice of the unary search operator: With the exception of problem swv15, both operators provide the same median result quality. In the other setups, if one of the two is better, it is most of the time nswap.

Therefore, we will consider the two setups ea_16384_nswap and ea_1024_nswap when evaluating the performance of our Evolutionary Algorithms.

3.4.1.3 Results on the JSSP

Let us now compare the results of the best two EA setups with those that we obtained with the hill climber.

Table 3.8 shows us that we can improve the best, mean, and median solution quality that we can get within three minutes of runtime when using our either of the two EA setups instead of the hill climber. The exception is case la24, where the hill climber already came close to the lower bound of the makespan and has a better best solution than ea_16384_nswap. Here, the best solution encountered now has a makespan which is only 0.7% longer than what is theoretically possible. Nevertheless, we find quite a tangible improvement in case swv15 on ea_1024_nswap.

Our ea_16384_nswap outperforms the four Evolutionary Algorithms from [122] both in terms of mean and best result quality on abz7 and la24. It does the same for HIMGA-Mutation, the worst of the four HIMGA variants introduced in [130], for abz7, la24, and yn4. It obtains better results than the PABC from [175] on swv15. On la24, both in terms of mean and best result, it outperforms also all six EAs from [2], both variants of the EA in [157], and the LSGA from [158]. The best solution quality for abz7 delivered by ea_16384_nswap is better than the best result found by the old Fast Simulated Annealing algorithm which was improved in [4].

The Gantt charts of the median solutions of ea_16384_nswap are illustrated in Figure 3.15. They appear only a bit denser than those in Figure 3.9.

More interesting are the progress diagrams of the ea_16384_nswap, ea_1024_nswap, hcr_65536_nswap, and hc_nswap algorithms given in Figure 3.16. Here we find big visual differences between the way the EAs and hill climbers proceed. The EAs spend the first 10 to 1000 ms to discover some basins of attraction of local optima before speeding up. The larger the population, the longer it takes them until this happens: The difference between ea_16384_nswap, is very obvious ea_1024_nswap in this respect. It is interesting to notice that the two problems where the EAs visually outperform the hill climber the most, swv15 and yn4, are also those with the largest search spaces (see Table 2.3). la24, however, which already can almost be solved by the hill climber and where there are the smallest differences in performance, is the smallest instance. The population used by the EA seemingly guards against premature convergence and allows it to keep progressing for a longer time.

We also notice that – for the first time in our endeavor to solve the JSSP – runtime is the limiting factor. If we would have 20 minutes instead of three, then we could not expect much improvement from the hill climbers. Even with restarts, they already improve very slowly towards the end of the computational budget. Tuning their parameter, e.g., increasing the time until a restart is performance, would probably not help then either. However, we can clearly see that ea_16384_nswap has not fully converged on neither abz7, swv15, nor on yn4 after the three minutes. Its median curve is still clearly pointing downwards. And even if it had converged: If we had a larger computational budget, we could increase the population size. The reason why the performance for larger populations in Figure 3.14 gets worse is very likely the limited time budget.

3.4.1.4 Exploration versus Exploitation

Naturally, when discussing EAs, we may ask why the population is helpful in the search. While we now have some anecdotal evidence that this is the case, we may also think about this more philosophically. Let us therefore do a thought experiment. If we would set \mu=1 and \lambda=1, then the EA would always remember the best solution we had so far and, in each generation, derive one new solution from it. If it is better, it will replace the remembered solution. Such an EA is actually a hill climber.

Now imagine what would happen if we would set \mu to infinity instead. We would not even complete one single generation. Instead, if \mu\rightarrow\infty, it would also take infinitely long to finish creating the first population of random solutions. This does not even require infinity \mu – \mu just needs to be large enough so that the complete computational budget (in our case, three minutes) is consumed before creating the initial, random candidate solutions is completed. In other words, the EA would then equal random sampling.

The parameter \mu basically allows us to tune between these two behaviors [214]! If we pick it small, our algorithm becomes more greedy. It will spend more time investigating (exploiting) the neighborhood of the current best solutions. It will trace down local optima faster but be trapped more easily in local optima as well.

If we set \mu to a larger value, we will keep more not-that-great solutions in its population. The algorithm spends more time exploring the neighborhoods of solutions which do not look that good, but from which we might eventually reach better results. The convergence is slower, but we are less likely to get trapped in a local optimum.

The question on which of the two to focus is known as the dilemma of Exploration versus Exploitation [48,70,205,208,216]. To make matters worse, theorists have proofed that there are scenarios where only a small population can perform well, while there are other scenarios where only a large population works well [45]. In other words, if we apply an EA, we always need to do at least some rudimentary tuning of \mu and \lambda.

3.4.2 Ingredient: Binary Search Operator

On one hand, keeping a population of the \mu>1 best solutions as starting points for further exploration helps us to avoid premature convergence. On the other hand, it also represents more information. The hill climber only used the information in current-best solution as guide for the search (and the hill climber with restarts used, additionally, the number of steps performed since the last improvement). Now we have a set of \mu selected points from the search space. These points have, well, been selected. At least after some time has passed in our optimization process, being selected means being good. If you compare the Gantt charts of the median solutions of ea_16384_nswap (Figure 3.15) and hcr_16384_1swap (Figure 3.9), you can see some good solutions for the same problem instances. These solutions differ in some details. Wouldn’t it be nice if we could take two good solutions and derive a solution in between, a new solution which is similar to both of its parents?

This is the idea of the binary search operator (also often referred to as recombination or crossover operator). By defining such an operator, we hope that we can merge the good characteristics of two selected solutions to obtain one new (ideally better) solution [53,109]. If we are lucky and that works, then ideally such good characteristics could aggregate over time [84,147].

How can we define a binary search operator for our JSSP representation? One possible idea would be to create a new encoded solution x' by processing both input points {x}1 and {x}2 from front to back and schedule their not-yet scheduled job IDs into x' similar to what we do in our representation mapping.

1. Allocate a data structure x' to hold the new point in the search space that we want to sample.
2. Set the index i where the next operation should be stored in x' to i=0.
3. Repeat
   1. Randomly choose one of the input points {x}1 or {x}2 with equal probability as source x.
   2. Select the first (at the lowest index) operation in x that is not marked yet and store it in variable J.
   3. Set x'_{i}=J.
   4. Increase i by one (i=i+1).
   5. If i=n*m, then all operations have been assigned. We exit and returning x'.
   6. Mark the first unmarked occurrence of J as already assigned in {x}1.
   7. Mark the first unmarked occurrence of J as already assigned in {x}2.

This can be implemented efficiently by keeping indices of the first unmarked element for both {x}1 and {x}2, which we do in Listing 3.14.

public class JSSPBinaryOperatorSequence
    implements IBinarySearchOperator<int[]> {
  public void apply(int[] x0, int[] x1,
      int[] dest, Random random) {
// omitted: initialization of arrays doneX0 and doneX1 (that
// remember the already-assigned operations from x0 and x1) of
// length=m*n to all false; and indices desti, x0i, x10 to 0
    for (;;) { // repeat until dest is filled, i.e., desti=length
// randomly chose a source point and pick next operation from it
      int add = random.nextBoolean() ? x0[x0i] : x1[x1i];
      dest[desti++] = add; // we picked a operation and added it
      if (desti >= length) { // if desti==length, we are finished
        return; // in this case, desti is filled and we can exit
      }

      for (int i = x0i;; i++) { // mark operation as done in x0
        if ((x0[i] == add) && (!doneX0[i])) { // find added job
          doneX0[i] = true;// found it and marked it
          break; // quit operation finding loop
        }
      }
      while (doneX0[x0i]) { // now we move the index x0i to the
        x0i++; // next, not-yet completed operation in x0
      }

      for (int i = x1i;; i++) { // mark operation as done in x1
        if ((x1[i] == add) && (!doneX1[i])) { // find added job
          doneX1[i] = true; // found it and marked it
          break; // quit operation finding loop
        }
      }
      while (doneX1[x1i]) { // now we move the index x1i to the
        x1i++; // next, not-yet completed operation in x0
      }
    } // loop back to main loop and to add next operation
  } // end of function
}

As we discussed in Section 2.6.2, our representation mapping processes the elements x\in\mathbb{X} from the front to the back and assigns the jobs to machines according to the order in which their IDs appear. It is a natural idea to design a binary operator that works in a similar way. Our sequence recombination processes two points from the search space {x}1 and {x}2 from their beginning to the end. At each step randomly picks one of them to extract the next operation, which is is then stored in the output x' and marked as done in both {x}1 and {x}2.

If it would, by chance, always choose {x}1 as source, then it would produce exactly {x}1 as output. If it would always pick {x}2 as source, then it would also return {x}2. If it would pick {x}1 for the first half of the times and then always pick {x}2, it would basically copy the first half of {x}1 and then assign the rest of the operations in exactly the order in which they appear in {x}2.

For illustration purposes, one example application of this operator is sketched in Figure 3.17. As input, we chose to points {x}1 and {x}2 from the search space for our demo instance. They encode two different corresponding Gantt charts, {y}1 and {y}2, with makespans of 202 and 182 time units, respectively.

Our operator begins by randomly choosing {x}1 as the source of the first operation for the new point x'. The first job ID in {x}1 is 2, which is placed as first operation into x'. We also mark the first occurrence of 2 in {x}2, which happens to be at position 4, as already scheduled. Then, the operator again randomly picks {x}1 as source for the next operation. The first not-yet marked element in {x}1 is now at the second 0, so it is placed into x' and marked as scheduled in {x}2, where the fifth element is thus crossed out. As next source, the operator, again, chooses&bsnp;{x}1. The first unmarked operation in {x}1 is 3 at position 3, which is added to x' and leads to the first element of {x}2 being marked. Finally, for picking the next operation, {x}2 is chosen. The first unmarked operation there has ID 1 and is located at index 2. It is inserted at index 4 into x'. It also occurs at index 4 in {x}1, which is thus marked. This process is repeated again and again, until x' is constructed completely, at which point all the elements of {x}1 and {x}2 are marked.

The application of our binary operator yields a new point x' which corresponds to the Gantt chart y' with makespan 192. This new candidate solution clearly inherits some characteristics from either of its parents.

3.4.3 Evolutionary Algorithm with Recombination

We now want to utilize this new operator in our EA. The algorithm now has two ways to create new offspring solutions: either via the unary operator (mutation, in EA-speak) or via the binary operator (recombination in EA-speak). We modify the original EA as follows.

3.4.3.1 The Algorithm (with Recombination)

We introduce a new parameter cr\in[0,1], the so-called crossover rate. It is used whenever we want to derive a new points in the search space from existing ones. It denotes the probability that we apply the binary operator (while we will otherwise apply the unary operator, i.e., with probability 1-cr). The basic (\mu+\lambda) Evolutionary Algorithm with recombination works as follows:

1. I\in\mathbb{X}\times\mathbb{R} be a data structure that can store one point x in the search space and one objective value z.
2. Allocate an array P of length \mu+\lambda of instances of I.
3. For index i ranging from 0 to \mu+\lambda-1 do
   1. Create a random point from the search space using the nullary search operator and store it in P_{i}.x.
   2. Apply the representation mapping y=\gamma(P_{i}.x) to get the corresponding candidate solution y.
   3. Compute the objective objective value of y and store it at index i as well, i.e., P_{i}.z=f(y).
4. Repeat until the termination criterion is met:
   4. Sort the array P in ascending order according to the objective values, i.e., such that the records r with better associated objective value r.z are located at smaller indices.
   5. Shuffle the first \mu elements of P randomly.
   6. Set the first source index p1=-1.
   7. For index i ranging from \mu to \mu+\lambda-1 do
      1. Set the first source index p1 to p1=(p1+1)\textrm{~mod~}\mu.
      2. Draw a random number c uniformly distributed in [0,1).
      3. If c is less than the crossover rate cr, then we apply the binary operator: A. Randomly choose a second index p2 from 0\ldots(\mu-1) such that p2\neq p1. B. Apply binary search operator to the points stored at index p1 and p2 and store result at index i, i.e., set P_{i}.x=\textrm{searchOp}_2(P_{p1}.x, P_{p2}.x).
      4. otherwise to step 4g.iii, i.e., if c\geq cr, then we apply the unary operator: C. Apply unary search operator to the point stored at index p1 and store result at index i, i.e., set P_{i}.x=\textrm{searchOp}_1(P_{p1}.x).
      5. Apply the representation mapping y=\gamma(P_{i}.x) to get the corresponding candidate solution y.
      6. Compute the objective objective value of y and store it at index i as well, i.e., P_{i}.z=f(y).
5. Return the candidate solution corresponding to the best record in P to the user.

public class EA<X, Y> extends Metaheuristic2<X, Y> {
  public void solve(IBlackBoxProcess<X, Y> process) {
// omitted: initialize local variables random, searchSpace, and
// array P of length mu+lambda
// first generation: fill population with random solutions
    for (int i = P.length; (--i) >= 0;) {
      X x = searchSpace.create(); // allocate point
      this.nullary.apply(x, random); // fill with random data
      P[i] = new Record<>(x, process.evaluate(x)); // evaluate
    }

    for (;;) { // main loop: one iteration = one generation
// sort the population: mu best records at front are selected
      Arrays.sort(P, Record.BY_QUALITY);
// shuffle the first mu solutions to ensure fairness
      RandomUtils.shuffle(random, P, 0, this.mu);
      int p1 = -1; // index to iterate over first parent

// overwrite the worse lambda solutions with new offsprings
      for (int index = P.length; (--index) >= this.mu;) {
        if (process.shouldTerminate()) { // we return
          return; // best solution is stored in process
        }

        Record<X> dest = P[index];
        p1 = (p1 + 1) % this.mu; // step the parent 1 index
        Record<X> sel = P[p1];
        if (random.nextDouble() <= this.cr) { // crossover!
          do { // find a second, different record
            p2 = random.nextInt(this.mu);
          } while (p2 == p1); // repeat until p1 != p2
// perform recombination of the two selected records
          this.binary.apply(sel.x, P[p2].x, dest.x, random);
        } else {
// create modified copy of parent using unary operator
          this.unary.apply(sel.x, dest.x, random);
        }
// map to solution/schedule and evaluate quality
        dest.quality = process.evaluate(dest.x);
      } // the end of the offspring generation
    } // the end of the main loop
  }
}

This algorithm, implemented in Listing 3.15 only differs from the variant in Section 3.4.1.1 by choosing whether to use the unary or binary operator to sample new points from the search space (steps 4g.iii.A, B, and C). If cr is the probability to apply the binary operator and we draw a random number c which is uniformly distributed in [0,1), then the probability that c<cr is exactly cr (see point iii).

3.4.3.2 The Right Setup

Unfortunately, with cr\in[0,1], a new algorithm parameter has emerged. It is not really clear whether a large or a small crossover rate is good. We already tested cr=0: The EA without recombination. Similar to our other small tuning experiments, let us compare the performance of different settings of cr. We investigate the crossover rates cr\in\{0, 0.05, 0.3, 0.98\}. We will stick with nswap as unary operator and keep \mu=\lambda.

From Figure 3.18, we immediately find that the large crossover rate cr=0.98 is performing much worse than using no crossover at all (cr=0) on all instances. Smaller rates cr\in\{0.05,0.3\} tend to be sometimes better and sometimes worse than cr=0, but there is no big improvement. On instance swv15, the binary operator does not really help. On instance la24, cr=0.05 performs best on mid-sized populations, while there is no distinguishable difference between cr=0.05 and cr=0 for large populations. On abz7 and yn4, cr=0.05 always seems to be a good choice. On these three instances, the population size \mu=\lambda=8192 in combination with cr=0.05 looks promising. We will call this setup ea_8192_5%_nswap in the following, where 5% stands for the crossover rate of 0.05 using the binary sequence operator.

3.4.3.3 Results on the JSSP

We can now investigate whether our results have somewhat improved.

The results in Table 3.9 show that a moderate crossover rate of 0.05 can indeed improve our algorithm’s performance – a little bit. For swv15, we already know that recombination was not helpful and this is confirmed in the table. On the three other instances, ea_8192_5%_nswap has better mean and median results than ea_1024_nswap, ea_8192_nswap, and ea_16384_nswap. On abz7, it can also slightly improve on the best result we got so far.

Interestingly, on swv15, the ea_1024_nswap stops improving in median after about 87 seconds, whereas all other EAs keep finding improvements even very close to the end of the 180 seconds budget. This probably means that they would also perform better than ea_1024_nswap if only we had more time. There might just not be enough time for any potential benefits of the binary operator to kick in. This could also be a valuable lesson: it does not help if the algorithm gives better results if it needs too much time. Any statement about an achieved result quality is only valid if it also contains a statement about the required computational budget. If we would have let the algorithms longer, maybe the setups using the binary operator would have given more saliently better results … but these would then be useless in our real-world scenario, since we only have 3 minutes of runtime.

By the way: It is very important to always test the cr=0 rate! Only by doing this, we can find whether our binary operator is designed properly. It is a common fallacy to assume that an operator which we have designed to combine good characteristics from different solutions will actually do that. If the algorithm setups with cr=0 would be better than those that use the binary operator, it would be a clear indication that we are doing something wrong. So we need to carefully analyze whether the small improvements that our binary operator can provide are actually significant. We therefore apply the same statistical approach as already used in Section 3.3.5.2 and later discussed in detail in Section 4.5.

In Table 3.10, we find that that EA using the binary sequence operator to generate 5% of the offspring solutions leads to significantly better results on abz7. It never performs significantly worse, not even on swv15, and the p-values are below \alpha but above \alpha' on the other two instances. Overall, this is not a very convincing result. Not enough for us to claim that our particular recombination operator is a very good invention – but enough to convince us that using it may sometimes be good and won’t hurt too much otherwise.

If we look at Figure 3.19, we can confirm that using the binary sequence operator at the low 5% rate does make some difference in how the median solution quality over time changes, although only a small one. On abz7, ea_8192_5%_nswap improves faster than the setup without recombination. On la24 and yn4, there may be a small advantage during the phase when the algorithm improves the fastest, but this could also be caused by the randomness of the search. On swv15, there is a similarly small disadvantage of ea_8192_5%_nswap.

In summary, it seems that using our binary operator is reasonable. Different from what we may have hoped for (and which would have been very nice for this book…), it does not improve the results by much. We now could try to design a different recombination operator in the hope to get better results, similar to what we did with the unary operator by moving from 1swap to nswap. We will not do this here – the interested reader is invited to do that by herself as an exercise.

As the end of this section, let me point out that binary search operators are a hot and important research topic right now. On one of the most well-known classical optimization problems, the Traveling Salesman Problem mentioned already back in the introduction in Section 1.1.2, they are part of the most efficient algorithms [153,178,197,219]. It also has theoretically been proven that a binary operator can speed-up optimization on some problems [61].

3.4.4 Ingredient: Diversity Preservation

In Section 3.4.1.4, we asked why a population is helpful for optimization. Our answer was that there are two opposing goals in metaheuristic optimization: On the one hand, we want to get results quickly and, hence, want that the algorithms quickly trace down to the bottom of the basins around optima. On the other hand, we want to get good results, i.e., better local optima, preferably global optima. A smaller population is good for the former and forsters exploitation. A larger population is good for the latter, as it invests more time on exploration.

Well. Not necessarily. Imagine we discover a good local optimum, a solution better than everything else we have in the population. Great. It will survive selection and we will derive offspring solutions from it. Since it is a local optimum, these will probably be worse. They might also encode the same solution as the parent, which is entirely possible in our JSSP scenario. But even if they are worse, they maybe just good enough to survive the next round of selection. Then, their (better) parent will also survive. We will thus get more offspring from this parent. But also offsprings from its surviving offsprings. And some of these may again be the same as the parent. If this process keeps continuing, the population may slowly be filling with copies of that very good local optimum. The larger our population, the longer it will take, of course. But unless we somehow encounter a different, similarly good or even better solution, it will probably happen eventually.

What does this mean? Recombination of two identical points in the search space should yield the very same point as output, i.e., the binary operator will become useless. This would leave only the unary operator as possible source of randomness. We then practical have one point in the search space to which only the unary operator is applied. Our EA has become a weird hill climber.

If we would want that, then we would have implemented a hill climber, i.e., a local search, instead. In order to enable global exploration and to allow for the binary search operators to work, it makes sense to try to preserve the diversity in the population [193].

Now there exist quite a few ideas how to do that [48,184,189] and we also discuss some concepts later in Section 5.1.2.4. Many of them are focused on penalizing candidate solutions which are too similar to others in the selection step. Similarity could be measured via computing the distance in the search space, the distance in the solution space, the difference of the objective values. The penalty could be achieved by using a so-called fitness as basis for selection instead of the objective value. In our original EA, the fitness would be the same as the objective value. In a diversity-preserving EA, we could add a penalty value to this base fitness for each solution based on the distance to the other solutions.

3.4.5 Evolutionary Algorithm with Clearing in the Objective Space

Let us now test whether a diversity preserving strategy can be helpful in an EA. We will only investigate one very simple approach: Avoiding objective value duplicates [78,184]. In the rest of this section, we will call this method clearing, as it can be viewed as the strictest possible variant of the clearing [161,162] applied in the objective space.

Put simply, we will ensure that all records that survive selection have different objective values. If two good solutions have the same objective value, we will discard one of them. This way, we will ensure that our population remains diverse. No single candidate solution can take over the population.

3.4.5.1 The Algorithm (with Recombination and Clearing)

We can easily extend our (\mu+\lambda) EA with recombination from Section 3.4.3.1 to remove duplicates of the objective value. We need to consider that a full population of \mu+\lambda records may contain less than \mu different objective values. Thus, in the selection step, we may obtain 1\leq u \leq \mu elements, where u can be different in each generation. If u=1, we cannot apply the binary operator regardless of the crossover rate cr.

1. I\in\mathbb{X}\times\mathbb{R} be a data structure that can store one point x in the search space and one objective value z.
2. Allocate an array P of length \mu+\lambda of instances of I.
3. For index i ranging from 0 to \mu+\lambda-1 do
   1. Create a random point from the search space using the nullary search operator and store it in P_{i}.x.
   2. Apply the representation mapping y=\gamma(P_{i}.x) to get the corresponding candidate solution y.
   3. Compute the objective objective value of y and store it at index i as well, i.e., P_{i}.z=f(y).
4. Repeat until the termination criterion is met:
   4. Sort the array P in ascending order according to the objective values, i.e., such that the records r with better associated objective value r.z are located at smaller indices.
   5. Iterate over P from front to end and delete all records with an objective value already seen in this iteration. The number of remaining records be w. Set the number u of selected records to u=\min\{w,\mu\}.
      
   6. Shuffle the first u elements of P randomly.
   7. Set the first source index p1=-1.
   8. For index i ranging from u to  \mu+\lambda-1 do
      1. Set the first source index p1 to p=(p1+1)\textrm{~mod~}u.
      2. Draw a random number c uniformly distributed in [0,1).
      3. If u>1 and c<cr, then we apply the binary operator: A. Randomly choose another index p2 from 0\ldots(u-1) such that p2\neq p1. B. Apply binary search operator to the points stored at index p1 and p2 and store result at index i, i.e., set P_{i}.x=\textrm{searchOp}_2(P_{p1}.x, P_{p2}.x).
      4. otherwise to *step 4h.iii, i.e., if c\geq cr or u=1, we apply the unary search operator to the point stored at index p1 and store result at index i, i.e., set P_{i}.x=\textrm{searchOp}_1(P_{p1}.x).
      5. Apply the representation mapping y=\gamma(P_{i}.x) to get the corresponding candidate solution y.
      6. Compute the objective objective value of y and store it at index i as well, i.e., P_{i}.z=f(y).
5. Return the candidate solution corresponding to the best record in P to the user.

public class EAWithClearing<X, Y>
    extends Metaheuristic2<X, Y> {
  public void solve(IBlackBoxProcess<X, Y> process) {
// Omitted: Initialize local variables random, searchSpace, set
// arrays P of length mu+lambda, and array T to null. Fill P with
// random solutions + evaluate.
    for (int i = P.length; (--i) >= 0;) {
      X x = searchSpace.create();
      this.nullary.apply(x, random);
      P[i] = new Record<>(x, process.evaluate(x));
      if (process.shouldTerminate()) { // we return
        return; // best solution is stored in process
      }
    }

    while (!process.shouldTerminate()) { // main loop
      RandomUtils.shuffle(random, P, 0, P.length); // make fair
      int u = Utils.qualityBasedClearing(P, this.mu);
// Now we have 1 <= u <= mu unique solutions.
      RandomUtils.shuffle(random, P, 0, u); // for fairness
      int p1 = -1; // index to iterate over first parent
// Overwrite the worse (mu + lambda - u) solutions.
      for (int index = P.length; (--index) >= u;) {
// Omitted: Quit loop if process.shouldTerminate()
        Record<X> dest = P[index]; // offspring
        p1 = (p1 + 1) % u; // parent 1 index
        Record<X> sel = P[p1]; // parent 1
        if ((u >= 2) && (random.nextDouble() <= this.cr)) {
          do { // find a second, different record
            p2 = random.nextInt(u);
          } while (p2 == p1); // Of course, can't be p1.
          this.binary.apply(sel.x, P[p2].x, dest.x, random);
        } else { // Otherwise: Mutation.
          this.unary.apply(sel.x, dest.x, random);
        }
        dest.quality = process.evaluate(dest.x);
      } // the end of the offspring generation
    } // the end of the main loop
  }
}

public static int qualityBasedClearing(Record<?>[] array,
    int max) {

  Arrays.sort(array, Record.BY_QUALITY);// best -> first

  int unique = 0;
  double lastQuality = Double.NEGATIVE_INFINITY; // impossible

  for (int index = 0; index < array.length; index++) {
    Record<?> current = array[index];
    double currentQuality = current.quality;
    if (currentQuality > lastQuality) { // unique so-far
      if (index > unique) { // need to move forward?
        Record<?> other = array[unique];
        array[unique] = current; // swap with first non-unique
        array[index] = other;
      }
      lastQuality = currentQuality; // update new quality
      if ((++unique) >= max) { // are we finished?
        return unique; // then quit: unique == max
      }
    }
  }

  return unique; // return number of unique: 1<=unique<=max
}

This algorithm, implemented in Listing 3.16 and using the routine given in Listing 3.17 differs from the variant in Section 3.4.3.1 mainly in step 4e. There, the sorted population P is processed from beginning to end. Whenever an objective value is found in a record which has already been encountered during this processing step, the record is removed. Since P is sorted, this means that the record at (zero-based) index k is deleted if and only if k>0 and P_{k}.z=P_{k-1}.z. As a result, the number u of selected records with unique objective value may be less than \mu (while always being greater or equal to 1). Therefore, we need to adjust the parts of the algorithm where parent solutions are selected for generating offsprings. Also, we generate \mu+\lambda-u offspring, to again obtain a total of \mu+\lambda elements.

In the actual implementation in Listing 3.16, we do not delete the records but move them to the end of the list, so we can re-use them later. We also stop processing P as soon as we have \mu unique records, as it does not really matter whether un-selected records are unique. This is slightly more efficient, but would be harder to write in pseudo-code.

We will name setups of this algorithm in the same manner as those of the original EA, except that we start the names with the prefix eac_ instead of ea_.

3.4.5.2 The Right Setup

With the simple diversity-preservation mechanism in place, we may wonder which population sizes are good. It is easy to see that findings from Section 3.4.1.2, where we found that \mu=\lambda=16'384 is reasonable, may longer hold: We know from Section 2.5.3 of the makespan for any solution on the JSSP instance abz7 is 656. It can be doubted whether it is even possible to generate 16'384 schedules with different makespans. If not, then the number u of selected records would always be less than \mu, which would make choosing a large \mu useless.

Since this time smaller population sizes may be interesting, we investigate all powers of 2 for \mu=\lambda from 4 to 65’536. In Figure 3.20, we find that the eac_mu_5%_nswap behave entirely different from those of ea_mu_5%_nswap. If clearing is applied, the smallest investigated setting, \mu=\lambda=4, seems to be the right choice. This setup has the best performance on abz7 and swv15, while being only a tiny bit worse than the best choices on la24 and yn4. Larger populations lead to worse results at the end of the computational budget of three minutes.

Why could this be? One possible reason could be that maybe not very many different candidate solutions are required. If only few are needed and we can maintain sufficiently many in a small population, then the advantage of the small population is that we can do many more iterations within the same computational budget.

This very small population size means that an EA with clearing in the objective space should be configured quite similar to a hill climber!

3.4.5.3 Results on the JSSP

We can now investigate whether our results have somewhat improved.

Table 3.11 shows that our EA with recombination and clearing in the objective space outperforms the EA without clearing on every single instance in terms of best, mean, and median result quality. Especially on instance swv15, the new algorithm performs much better than ea_8192_5%_nswap. Most remarkable is that we can even solve instance la24 to optimality once. The median solutions of our new algorithm variant are illustrated in Figure 3.21. Compared to Figure 3.15, the result on swv15 has improved: especially the upper-left and lower-right corners of the Gantt chart, where only few jobs are scheduled, have visibly become smaller.

We plot the discovered optimal solution for la24 in Figure 3.23. Comparing it with the median solution for la24 in Figure 3.21, time was saved, e.g., by arranging the jobs in the top-left corner in a tighter pattern.

From the progress charts plotted in Figure 3.22, we can confirm that eac_4_5%_nswap indeed behaves more similar to the hill climber hcr_65536_nswap than to the other EA ea_8192_5%_nswap. This is due to its small population size. However, unlike the hill climber, its progress curve keeps going down for a longer time. After only 0.1 seconds, it keeps producing better results in median.

We can confirm that even the simple and rather crude pruning in the objective space can significantly improve the performance of our EA. We did not test any sophisticated diversity preservation method. We also did not re-evaluate which crossover rate cr and which unary operator (1swap or nswap) works best in this scenario. This means that we might be able to squeeze out more performance, but we will leave it at this. The small populations working so well make us curious, however.

3.4.6 (1+1) EA

The (1+1) EA is a special case of the (\mu+\lambda) EA with \mu=1 and \lambda=1. Since the number of parents is \mu=1, it does not apply the binary recombination operator. We explicitly discuss it here because it is usually defined slightly differently from what we did in Section 3.4.1.1 and implemented as Listing 3.13.

3.4.6.1 The Algorithm

The (1+1) EA works like a hill climber (see Section 3.3), with the difference that the new solution replaces the current one if it is better or equally good, instead of just replacing it if it is better.

1. Create one random point x in the search space \mathbb{X} using the nullary search operator.
2. Map the point x to a candidate solution y by applying the representation mapping y=\gamma(x).
3. Compute the objective value by invoking the objective function z=f(y).
4. Repeat until the termination criterion is met:
   1. Apply the unary search operator to x to get a slightly modified copy x' of it.
   2. Map the point x' to a candidate solution y' by applying the representation mapping y'=\gamma(x').
   3. Compute the objective value z' by invoking the objective function z'=f(y').
   4. If z'\leq z, then store x' in x, store y' in y, and store z' in z.
5. Return the best encountered objective value z and the best encountered solution y to the user.

This algorithm is implemented in Listing 3.18. The only difference to the hill climber in Section 3.3.2.1 and Listing 3.10 is the z'\leq z in point 4.d instead of z'<z: the new solution does not need to be strictly better to win over the old one, it is sufficient if it is not worse. In the (\mu+\lambda) EA implementation that we had before, we could also set \mu=1 and \lambda=1. The algorithm would be slightly different from the (1+1) EA: In the (1+1) EA as defined here, the new solution y' will always win against the current solution y if z'=z, whereas in our shuffle-and-sort method of the population as implemented in Listing 3.13, either one could win with probability 0.5.

public class EA1p1<X, Y> extends Metaheuristic1<X, Y> {
  public void solve(IBlackBoxProcess<X, Y> process) {
// initialize local variables xCur, xBest, random
    X xCur = process.getSearchSpace().create();
    X xBest = process.getSearchSpace().create();
    Random random = process.getRandom();// get random gen

// create starting point: a random point in the search space
    this.nullary.apply(xBest, random); // xBest=random point
    double fBest = process.evaluate(xBest); // map & evaluate

    while (!process.shouldTerminate()) {
// create a slightly modified copy of xBest and store in xCur
      this.unary.apply(xBest, xCur, random);
// map xCur from X to Y and evaluate candidate solution
      double fCur = process.evaluate(xCur);
      if (fCur <= fBest) { // we found a not-worse solution
// remember best objective value and copy xCur to xBest
        fBest = fCur;
        process.getSearchSpace().copy(xCur, xBest);
      } // otherwise, i.e., fCur > fBest: just forget xCur
    } // until time is up
  } // process will have remembered the best candidate solution
}

3.4.6.2 Results on the JSSP

We can now apply our (1+1) EA to the four JSSP instances either using the 1swap operator (ea_1+1_1swap) or the nswap operator (ea_1+1_nswap). Astonishingly, Table 3.12 reveals that it performs better than our best EA so far, namely eac_4_5%_nswap. Both (1+1) EA setups also perform much better than our hill climber. The log-scaled Figure 3.25 shows that the two EAs without population have better median solution almost always during the runs. And the Gantt charts of the median solutions of ea_1+1_1swap, illustrated in Figure 3.24, again appear denser.

3.4.6.3 Discussion

The big advantage of EAs is that their population guards against premature convergence. This can bring better end results, but comes at a trade-off of slower convergence [214]. The fact that the EA without population performs best in our experimental setup is annoying for the author. It does not mean that this is always the case, though. But why is it the case here:

The answer consists, most likely, of two parts: First, as we know, our search space \mathbb{X} is much larger than the solution space \mathbb{Y}, i.e., |\mathbb{X}|\gg|\mathbb{Y}. If we have an integer string x_1\in\mathbb{X} representing a Gantt chart y_1=\gamma(x_1), then we can swap two jobs and get a new string x_2\in\mathbb{X} with x_2\neq x_2, but this does not necessarily mean that this new string maps to a different solution. It could well be that \gamma(x_1)=\gamma(x_2). Then, we have made a neutral move. Of course, our move could also be neutral if \gamma(x_1)\neq\gamma(x_2) but f(\gamma(x_1))=f(\gamma(x_2)) Our (1+1) EA can drift along a network of points in the search space which all map to solutions with the same makespan. This means that the (1+1) EA can explore far beyond the neighborhood visible to the hill climber. By drifting over the network, it may eventually reach a point which maps to a better solution. This reasoning is supported by the fact that our (1+1) EA perform much better than the hill climbers. You can find a more detailed discussion on neutrality in Section 5.4.1.

The second reason is probably that the three minutes of runtime are simply not enough for the EAs with the really big populations to reap the benefit of being resilient against premature convergence.

3.4.7 Summary

In this chapter, we have introduced Evolutionary Algorithms as methods for global optimization. We have learned the two key concepts that distinguish them from local search: the use of a population of solutions and of a binary search operator. We found that even the former alone can already outperform the simple hill climber with restarts. While we were a bit unlucky with our choice of the binary operator, our idea did work at least a bit.

We then noticed that populations are only useful if they are diverse. If all the elements in the population are very similar, then the binary operator stops working and the EA has converged. From our experiments with the hill climber, we already know premature convergence to a local optimum as something that should be avoided.
It therefore makes sense to try adding a method for enforcing diversity as another ingredient into the algorithm. Our experiments with a very crude diversity enhancing method – only allowing one solution per unique objective value in the population – confirmed that this can lead to better results.

We also found that the strict criterion to only accept solutions which are strictly better than the current one, as practiced in our hill climbers Section 3.3 may not be a good idea. Accepting solutions which are not worse than the current one allows for drift in the search space, which can yield better results.

3.5 Simulated Annealing

So far, we have only discussed two variants of local search: the hill climbing algorithm and the (1+1) EA. A hill climbing algorithm is likely to get stuck at local optima, which may vary in quality. We found that we can utilize this variance of the result quality by restarting the optimization process when it could not improve any more in Section 3.3.3. Such a restart is costly, as it forces the local search to start completely from scratch (while we, of course, remember the best-ever solution in a variable hidden from the algorithm). The (1+1) EA may additionally utilize drift over solutions of the same quality, but there is no reason to assume that it cannot get stuck in local optima as well.

Another way to look at this is the following: A schedule which is a local optimum probably is somewhat similar to what the globally optimal schedule would look like. It must, obviously, also be somewhat different. This difference is shaped such that it cannot be conquered by the unary search operator that we use, because otherwise, the basic hill climber could already move from the local to the global optimum. If we do a restart, we also dispose of the similarities to the global optimum that we have already discovered. We will subsequently spend time to re-discover them in the hope that this will happen in a way that allows us to eventually reach the global optimum itself. But maybe there is a less-costly way? Maybe we can escape from a local optimum without discarding the entirety good solution characteristics we already have discovered?

3.5.1 Idea: Accepting Worse Solutions with Decreasing Probability

Simulated Annealing (SA) [47,116,125,164] is a local search which provides another approach to escape local optima [188,205]. The algorithm is inspired by the idea of simulating the thermodynamic process of annealing using statistical mechanics, hence the naming [144]. Instead of restarting the algorithm when reaching a local optimum, it tries to preserve the parts of the current solution by permitting search steps towards worsening objective values. This algorithm therefore introduces three principles:

1. Worse candidate solutions are sometimes accepted, too.
2. The probability P of accepting them is decreases with increasing differences \Delta E of their objective values to the current solution.
3. The probability also decreases with the number of performed search steps.

These three principles are injected into the main loop of the hill climber. This is realized as follows.

Let us assume that x\in\mathbb{X} is the current point that our local search maintains. x'\in\mathbb{X} is the newly sampled point, i.e., the result of the application of the unary search operator to x. Then, \Delta E be the difference between the objective value corresponding to x' and x. In other words, if \gamma is the representation mapping and f the objective function, then:

\Delta E = f(\gamma(x')) - f(\gamma(x)) \qquad(3.1)

Clearly, if we try to minimize the objective function f, then \Delta E < 0 means that x' is better than x since f(\gamma(x')) < f(\gamma(x)). If \Delta E>0, on the other hand, the new solution is worse. The probability P to overwrite x with x' then be

P = \left\{\begin{array}{rl} 1 & \text{if~}\Delta E \leq 0\\ e^{-\frac{\Delta E}{T}} & \text{if~}\Delta E >0 \land T > 0\\ 0 & \text{otherwise~}(\Delta E > 0 \land T=0) \end{array} \right. \qquad(3.2)

In other words, if the new point x' is actually better (or, at least, not worse) than the current point x, i.e., \Delta E \leq 0, then we will definitely accept it. (This is exactly the acceptance criterion used in the (1+1) EA.) However, the new solution is not necessarily rejected otherwise: If the new point x' is worse (\Delta E > 0), then the acceptance probability

1. gets smaller the larger \Delta E is and
2. gets smaller the smaller the so-called temperature T\geq 0 is.

Both the temperature T>0 and the objective value difference \Delta E>0 enter Equation 3.2 in the exponential term and the two above points follow from e^{-a}<e^{-b}\forall a>b. We also have e^{-a}\in(0,1)\forall a>0, so it can be used as probability value.

The temperature will be changed automatically such that it decreases and approaches zero with a rising number \tau of algorithm iterations, i.e., the performed objective function evaluations. The optimization process is initially hot and T is high. Then, even significantly worse solutions may be accepted Over time, the process cools down and T decreases. The search slowly accepts fewer and fewer worse solutions and more likely such which are only a bit worse. Eventually, at temperature T=0, the algorithm only accepts better solutions. In other words, T is actually a monotonously decreasing function T(\tau) called the temperature schedule. It holds that \lim_{\tau\rightarrow+\infty} T(\tau) = 0.

3.5.2 Ingredient: Temperature Schedule

The temperature schedule T(\tau) determines how the temperature changes over time (where time is measured in algorithm steps \tau). It begins with an start temperature T_s at \tau=1. This temperature is the highest, which means that the algorithm is more likely to accept worse solutions. It will then behave a bit similar to a random walk and put more emphasis on exploring the search space than on improving the objective value. As time goes by and \tau increases, T(\tau) decreases and may even reach 0 eventually. Once T gets small enough, then Simulated Annealing will behave exactly like a hill climber and only accepts a new solution if it is better than the current solution. This means the algorithm tunes itself from an initial exploration phase to strict exploitation.

Consider the following perspective: An Evolutionary Algorithm allows us to pick a behavior in between a hill climber and a random sampling algorithm by choosing a small or large population size. The Simulated Annealing algorithm allows for a smooth transition of a random search behavior towards a hill climbing behavior over time.

public abstract class TemperatureSchedule
    implements ISetupPrintable {
  public double startTemperature;

  public abstract double temperature(long tau);
}

The ingredient needed for this tuning, the temperature schedule, can be expressed as a class implementing exactly one simple function that translates an iteration index \tau to a temperature T(\tau), as defined in Listing 3.19.

The two most common temperature schedule implementations may be the exponential and the logarithmic schedule.

3.5.2.1 Exponential Temperature Schedule

In an exponential temperature schedule, the temperature decreases exponentially with time (as the name implies). It follows Equation 3.3 and is implemented in Listing 3.20. Besides the start temperature T_s, it has a parameter \epsilon\in(0,1) which tunes the speed of the temperature decrease. Higher values of \epsilon lead to a faster temperature decline.

T(\tau) = T_s * (1 - \epsilon) ^ {\tau- 1} \qquad(3.3)

public static class Exponential
    extends TemperatureSchedule {
  public double epsilon;

  public double temperature(long tau) {
    return (this.startTemperature
        * Math.pow((1d - this.epsilon), (tau - 1L)));
  }
}

3.5.2.2 Logarithmic Temperature Schedule

The logarithmic temperature schedule will prevent the temperature from becoming very small for a longer time. Compared to the exponential schedule, it will thus longer retain a higher probability to accept worse solutions. It obeys Equation 3.4 and is implemented in Listing 3.21.. It, too, has the parameters \epsilon\in(0,\infty) and T_s. Larger values of \epsilon again lead to a faster temperature decline.

T(\tau) = \frac{T_s}{\ln{\left(\epsilon(\tau-1)+e\right)}} \qquad(3.4)

public static class Logarithmic
    extends TemperatureSchedule {
  public double epsilon;
  public double temperature(long tau) {
    if (tau >= Long.MAX_VALUE) {
      return 0d;
    }
    return (this.startTemperature
        / Math.log(((tau - 1L) * this.epsilon) + Math.E));
  }
}

3.5.3 The Algorithm

Now that we have the blueprints for temperature schedules, we can completely define our SA algorithm and implement it in Listing 3.22.

1. Create random point x in the search space \mathbb{X} by using the nullary search operator.
2. Map the point x to a candidate solution y by applying the representation mapping y=\gamma(x).
3. Compute the objective value by invoking the objective function z=f(y).
4. Store y in y_{b} and z in z_{b}, which we will use to preserve the best-so-far results.
5. Set the iteration counter \tau to \tau=1.
6. Repeat until the termination criterion is met:
   1. Set \tau=\tau+1.
   2. Apply the unary search operator to x to get the slightly modified copy x' of it.
   3. Map the point x' to a candidate solution y' by applying the representation mapping y'=\gamma(x').
   4. Compute the objective value z' by invoking the objective function z'=f(y').
   5. If z'\leq z, then
      1. Store x' in x and store z' in z.
      2. If z'\leq z_{b}, then store y' in y_{b} and store z' in z_{b}.
      3. Perform next iteration by going to step 6.
   6. Compute the temperature T according to the temperature schedule, i.e., set T=T(\tau).
   7. If T\leq 0 the perform next iteration by going to step 6.
   8. Set \Delta E = z' - z (see Equation 3.1).
   9. Compute P=e^{-\frac{\Delta E}{T}} (see Equation 3.2).
   10. Draw a random number r uniformly distributed in [0,1).
   11. If r\leq P, then store x' in x, store z' in z, and perform next iteration by going to step 6.
7. Return best encountered objective value z_{b} and the best encountered solution  z_{b} to the user.

public class SimulatedAnnealing<X, Y>
    extends Metaheuristic1<X, Y> {
  public TemperatureSchedule schedule;

  public void solve(IBlackBoxProcess<X, Y> process) {
// init local variables xNew, xCur, random
// create starting point: a random point in the search space
    this.nullary.apply(xCur, random); // put random point in xCur
    double fCur = process.evaluate(xCur); // map & evaluate
    long tau = 1L; // initialize step counter to 1

    do { // repeat until budget exhausted
// create a slightly modified copy of xCur and store in xNew
      this.unary.apply(xCur, xNew, random);
      ++tau; // increase step counter
// map xNew from X to Y and evaluate candidate solution
      double fNew = process.evaluate(xNew);
      if ((fNew <= fCur) || // accept if better solution OR
          (random.nextDouble() < // probability is e^(-dE/T)
              Math.exp((fCur - fNew) / // -dE == -(fNew-fCur)
                  this.schedule.temperature(tau)))) {
// accepted: remember objective value and copy xNew to xCur
        fCur = fNew;
        process.getSearchSpace().copy(xNew, xCur);
      } // otherwise fNew > fCur and not accepted
    } while (!process.shouldTerminate()); // until time is up
  } // process will have remembered the best candidate solution
}

There exist a several proofs [89,156] showing that, with a slow-enough cooling schedule, the probability that Simulated Annealing will find the globally optimal solution approaches 1. However, the runtime one would need to invest to actually cash in on this promise exceeds the time needed to enumerate all possible solutions [156]. In Section 1.2.1 we discussed that we are using metaheuristics because for many problems, we can only guarantee to find the global optimum if we invest a runtime growing exponentially with the problem scale (i.e., proportional to the size of the solution space). So while we have a proof that SA will eventually find a globally optimal solution, this proof is not applicable in any practical scenario and we instead use SA as what it is: a metaheuristic that will hopefully give us good approximate solutions in reasonable time.

3.5.4 The Right Setup

Our algorithm has four parameters:

• the start temperature T_s,
• the parameter \epsilon,
• the type of temperature schedule to use (here, logarithmic or exponential), and
• the unary search operator (in our case, we could use 1swap or nswap).

We will only consider 1swap as choice for the unary operator and focus on the exponential temperature schedule. We have two more parameters to set: T_s and \epsilon and thus refer to the settings of this algorithm with the naming scheme sa_exp_Ts_epsilon_1swap.

At first glance, it seems entirely unclear how what to do with these parameters. However, we may get some ideas about their rough ranges if we consider Simulated Annealing as an improved hill climber. Then, we can get some leads from our prior experiments with that algorithm.

In Table 3.13, we print the standard deviation sd of the final result qualities that our hc_1swap algorithm achieved on our four JSSP instances. This tells us something about how far the different local optima at which hc_1swap can get stuck are apart in terms of objective value. The value of sd ranges from 28 on abz7 to 137 on swv15. The median standard deviation over all four instances is about 50. Thus, accepting a solution which is worse by 50 units of makespan, i.e., with \Delta E\approx 50, should be possible at the beginning of the optimization process.

How likely should accepting such a value be? Unfortunately, we are again stuck at making an arbitrary choice – but at least we can make a choice from within a well-defined region: Probabilities must be in [0,1] and can be understood relatively intuitively, whereas it was completely unclear in what range reasonable temperatures would be located. Let us choose that the probability P_{50} to accept a candidate solution that is 50 makespan units worse than the current one, i.e., has \Delta E = 50, should be P_{50}=0.1 at the beginning of the search. In other words, there should be a 10% chance to accept such a solution at \tau=1. At \tau=1, T(\tau)=T_s for both temperature schedules. Of course, at \tau=1 we do not really use the probability formula to decide whether or not to accept a solution, but we can expect that the temperature at \tau=2 would still be very similar to T_s and we are using rough and rounded approximations anyway, so let’s not make our life unnecessarily complicated here. We can now solve Equation 3.2 for T_s:

\begin{array}{rl} P_{50} =& e^{-\frac{\Delta E}{T(\tau)}}\\ 0.1 =& e^{-\frac{50}{T_s}}\\ \ln{0.1}=& -\frac{50}{T_s}\\ T_s=&-\frac{50}{\ln{0.1}}\\ T_s\approx&21.7 \end{array}

A starting temperature T_s of approximately T_s=20 seems to be suitable in our scenario. Of course, we just got there using very simplified and coarse estimates. If we would have chosen P_{50}=0.5, we would have gotten T_s\approx 70 and if we additionally went with the maximum standard deviation 137 instead of the median one, we would obtain T_s\approx200. But at least we have a first understanding of the range where we will probably find good values for T_s.

But what about the \epsilon parameters? In order to get an idea for how to set it, we first need to know a proper end temperature T_e, i.e., the temperature which should be reached by the end of the run. It cannot be 0, because while both temperature schedules do approach zero for \tau\rightarrow\infty, they will not actually become 0 for any finite number \tau of iterations.

So we are stuck with the task to pick a suitably small value for T_e. Maybe here our previous findings from back when we tried to restart the hill climber can come in handy. In Section 3.3.3.2, we learned that it makes sense to restart the hill climber after L=16'384 unsuccessful search steps. So maybe a terminal state for our Simulated Annealing could be a scenario where the probability P_e of accepting a candidate solution which is \Delta E=1 makespan unit worse than the current one should be P_e=1/16'384? This would mean that the chance to accept a candidate solution being marginally worse than the current one would be about as large as making a complete restart in hcr_16384_1swap. Of course, this is again an almost arbitrary choice, but it at least looks like a reasonable terminal state for our Simulated Annealing algorithm. We can now solve Equation 3.2 again to get the end temperature T_e:

\begin{array}{rl} P_e =& e^{-\frac{\Delta E}{T(\tau)}}\\ 1/16'384 =& e^{-\frac{1}{T_e}}\\ \ln{(1/16'384)}=& -\frac{1}{T_e}\\ T_e=&-\frac{1}{\ln{(1/16'384)}}\\ T_e\approx&0.103 \end{array}

We choose a final temperature of T_e=0.1. But when should it be reached? In Table 3.13, we print the total number of function evaluations (FEs) that our hc_1swap algorithm performed on the different problem instances. We find that it generated and evaluated between 22 million on swv15 and 71 million on la24 candidate solutions.⁵ The overall median is at about 30 million FEs within the 3 minute computational budget. From this, we can conclude that after about 30 million FEs, we should reach approximately T_e=0.1. We can solve Equation 3.3 for \epsilon to configure the exponential schedule:

\begin{array}{rl} T(\tau) =& T_s * (1 - \epsilon) ^ {\tau- 1} \\ T(30'000'000) =& 20 * (1 - \epsilon) ^ {30'000'000 - 1} \\ 0.1 =& 20 * (1 - \epsilon) ^ {29'999'999} \\ 0.1/20 =& (1 - \epsilon) ^ {29'999'999} \\ 0.005^{1/29'999'999} =& 1 - \epsilon \\ \epsilon =& 1 - 0.005^{1/29'999'999}\\ \epsilon \approx& 1.776*10^{-7} \end{array}

We can conclude, for an exponential temperature schedule, settings for \epsilon somewhat between 1!\cdot\!10^{-7} and 2\!\cdot\!10^{-7} seem to be a reasonable choice if the start temperature T_s is set to 20.

In Figure 3.26, we illustrate the behavior of the exponential temperature schedule for starting temperature T_s=20 and the six values 5\!\cdot\!10^{-8}, 1\!\cdot\!10^{-7}, 1.5\!\cdot\!10^{-7}, 2\!\cdot\!10^{-7}, 4\!\cdot\!10^{-7}, and 8\!\cdot\!10^{-7}. The sub-figure on top shows how the temperature declines over the performed objective value evaluations. Starting at T_s=20 it reaches close to zero for \epsilon\geq 1.5\!\cdot\!10^{-7} after about \tau=30'000'000 FEs. For the smaller \epsilon values, the temperature would need longer to decline, while for larger values, it declines quicklier. The next three sub-figures show how the probability to accept candidate solutions which are worse by \Delta E units of the objective value decreases with \tau for \Delta E\in\{1, 3, 10\}. This decrease is, of course, the direct result of the temperature decrease. Solutions with larger \Delta E clearly have a lower probability of being accepted. The larger \epsilon, the earlier and faster does the acceptance probability decrease.

In Figure 3.27, we illustrate the normalized median result quality that can be obtained by Simulated Annealing with starting temperature T_s=20, exponential schedule, and 1swap operator for different values of the parameter \epsilon, including those from Figure 3.26. Interestingly, it turns out that \epsilon=2 is the best choice for the instances abz7, swv15, and yn4, which is quite close to what we could expect from our calculation. Smaller values will make the temperature decrease more slowly and lead to too much exploration and too little exploitation, as we already know from Figure 3.26. They work better on instance la24, which is no surprise: From Table 3.13, we know that on this instance, we can conduct more than twice as many objective value evaluations than on the others within the three minute budget: On la24, we can do enough steps to let the temperature decrease sufficiently even for smaller \epsilon.

3.5.5 Results on the JSSP

We can now evaluate the performance of our Simulated Annealing approach on the JSSP. We choose the setup with exponential temperature schedule, the 1swap operator, the starting temperature T_s=20, and the parameter \epsilon=2\!\cdot\!10^{-7}. For simplicity, we refer to it as sa_exp_20_2_1swap.