6.7. The Job-Shop Problem: and now with Local Search!

1. C++ code:
   1. jobshop_ls.h
   2. jobshop_ls1.cc
   3. jobshop_ls2.cc
   4. jobshop_ls3.cc
2. Data files:
   1. first_example_jssp.txt
   2. abz9

We have seen in the previous section how to implement Local Search on our dummy example. This time, we apply Local Search on a real problem and present the real thing: the MakeNextNeighbor() method, the delta and deltadelta Assignments[1].

To solve the Job-Shop Problem, we’ll define two basic LocalSearchOperators. First, we’ll apply them separately and then we’ll combine them to get better results. In doing so, we will discover that Local Search is very sensitive to the initial solution used to start it and that the search is path-dependent.

6.7.1. LocalSearchOperators: the real thing!

Until now, we only have redefined the MakeOneNeighbor()[2] method of a LocalSearchOperator but as we have seen in the sub-section The basic Local Search algorithm and the callback hooks for the SearchMonitors , the real method called by the MakeOneNeighbor DecisionBuilder is MakeNextNeighbor().

Before we dissect MakeNextNeighbor(), we quickly explain again what the deltas are.

6.7.1.1. Deltas and DeltaDeltas

The idea behind the Deltas and DeltaDeltas is really simple: efficiency. Only the modified part of the solution is broadcast:

• Delta: the difference between the initial solution that defines the neighborhood and the current candidate solution.
• DeltaDelta: the difference between the current candidate solution and the previous candidate solution.

Delta and DeltaDelta are just Assignments only containing the changes.

6.7.1.2. MakeNextNeighbor()

The signature of the MakeNextNeighbor() method is:

bool MakeNextNeighbor(Assignment* delta, Assignment* deltadelta)

This method constructs the delta and deltadelta corresponding to the new candidate solution and returns true. If the neighborhood has been exhausted, i.e. the LocalSearchOperator cannot find another candidate solution, this method returns false.

When you write your own MakeNextNeighbor() method, you have to provide the new delta but you can skip the deltadelta if you prefer. This deltadelta can be convenient when you define your filters and you can gain some efficiency over the sole use of deltas.

To help you construct these deltas, we provide an inner mechanism that constructs automatically these deltas when you use the following self-explanatory setters:

• for IntVarLocalSearchOperators only:
  • SetValue(int64 index, int64 value);
• for SequenceVarLocalSearchOperators only:
  • SetForwardSequence(int64 index, const std::vector<int>& value);
  • SetBackwardSequence(int64 index, const std::vector<int>& value);
• for both:
  • Activate(int64 index);
  • Deactivate(int64 index).

If you only use these methods to change the current solution, you then can automatically construct the deltas by calling the ApplyChanges() method and revert these changes by calling the RevertChanges() method.

We recommend to use the following template to define your MakeNextNeighbor() method:

virtual bool MakeNextNeighbor(Assignment* delta,
                              Assignment* deltadelta) {
  CHECK_NOTNULL(delta);
  while (true) {
    RevertChanges(true);
    if (NEIGHBORHOOD EXHAUSTED) {  // CONDITION TO BE DEFINED
      return false;
    }

    // CONSTRUCT NEW CANDIDATE SOLUTION
    ...
    if (ApplyChanges(delta, deltadelta)) {
      return true;
    }
  }
  return false;
}

Currently, ApplyChanges() always returns true but this might change in the future and then you might have to revert the changes, hence the while() loop.

We also provide several getters:

• for IntVarLocalSearchOperators only:
  • int64 Value(int64 index);
  • IntVar* Var(int64 index);
  • int64 OldValue(int64 index);
• for SequenceVarLocalSearchOperators only:
  • const std::vector<int>& Sequence(int64 index);
  • SequenceVar* Var(int64 index);
  • const std::vector<int>& OldSequence(int64 index);
• for both:
  • bool IsIncremental();
  • bool Activated(int64 index);

Why would I want to use MakeNextNeighbor() instead of MakeOneNeighbor()?

One reason is efficiency: you skip one callback. But the real reason is that you might need other methods than the ones that are provided to construct your candidate solution. In this case, you have no other choice than to reimplement the MakeNextNeighbor() method.

6.7.2. The initial solution

We let the CP solver construct the initial solution for us. What about reusing the DecisionBuilder defined in the section An implementation of the disjunctive model and grab its first feasible solution?

// This decision builder will rank all tasks on all machines.
DecisionBuilder* const sequence_phase =
    solver.MakePhase(all_sequences, Solver::SEQUENCE_DEFAULT);

// After the ranking of tasks, the schedule is still loose and any
// task can be postponed at will. Fix the objective variable to its
// minimum value.
DecisionBuilder* const obj_phase =
    solver.MakePhase(objective_var,
                     Solver::CHOOSE_FIRST_UNBOUND,
                     Solver::ASSIGN_MIN_VALUE);

// Initial solution for the Local Search.
Assignment* const first_solution = solver.MakeAssignment();
first_solution->Add(all_sequences);
first_solution->AddObjective(objective_var);

// Store the first solution.
DecisionBuilder* const store_db =
                       solver.MakeStoreAssignment(first_solution);

// The main decision builder (ranks all tasks, then fixes the
// objective_variable).
DecisionBuilder* const first_solution_phase =
    solver.Compose(sequence_phase, obj_phase, store_db);

LOG(INFO) << "Looking for the first solution";
const bool first_solution_found = solver.Solve(first_solution_phase);
if (first_solution_found) {
  LOG(INFO) << "Solution found with makespan = "
            << first_solution->ObjectiveValue();
} else {
  LOG(INFO) << "No initial solution found!";
  return;
}

If you have some troubles to follow, go back to the section An implementation of the disjunctive model to understand the sequence_phase and obj_phase DecisionBuilders. Here, we simply add a StoreAssignment DecisionBuilder at the leaf of the search tree to collect the solutions with the first_solution_phase DecisionBuilder. Our initial solution will be stored in the first_solution Assignment.

Next, we define a first LocalSearchOperator.

6.7.3. Exchanging two IntervalVars on a SequenceVar

You’ll find the code in the file jobshop_ls1.cc and the SwapIntervals operator in the file jobshop_ls.

The idea of exchanging two IntervalVars on a SequenceVar is very common and the corresponding operator is often referred to as the 2-opt-, 2-exchange- or swap- operator.

We implement a basic version that systematically exchanges all IntervalVars for all SequenceVars one after the other in the order given by the std::vectors. We use three indices:

• int current_var_: the index of the processed SequenceVar;
• int current_first_: the index of the first IntervalVar variable to swap;
• int current_second_: the index of the second IntervalVar variable to swap.

We proceed sequentially with the first SequenceVar (current_var_ = 0) and exchange the first and second IntervalVars, then the first and the third IntervalVars and so on until exhaustion of all possibilities. Here is the code that increments these indices to create each candidate solution:

bool Increment() {
  const SequenceVar* const var = Var(current_var_);
  if (++current_second_ >= var->size()) {
    if (++current_first_ >= var->size() - 1) {
      current_var_++;
      current_first_ = 0;
    }
    current_second_ = current_first_ + 1;
  }
  return current_var_ < Size();
}

This Increment() method returns a bool that indicates when the neighborhood is exhausted, i.e. it returns false when there are no more candidate to construct. Size() and Var() are helper methods defined in the SequenceVarLocalSearchOperator class. We start with current_var_, current_first_ and current_second_ all set to 0. Pay attention to the fact that current_first_ and current_second_ are also updated inside the if conditions.

We are now ready to define the OnStart() and MakeNextNeighbor() methods.

The OnStart() method is straightforward:

virtual void OnStart() {
  current_var_ = 0;
  current_first_ = 0;
  current_second_ = 0;
}

For the MakeNextNeighbor() method, we use our template:

virtual bool MakeNextNeighbor(Assignment* delta,
                              Assignment* deltadelta) {
  CHECK_NOTNULL(delta);
  while (true) {
    RevertChanges(true);
    if (!Increment()) {
      return false;
    }

    std::vector<int> sequence = Sequence(current_var_);
    const int tmp = sequence[current_first_];
    sequence[current_first_] = sequence[current_second_];
    sequence[current_second_] = tmp;
    SetForwardSequence(current_var_, sequence);

    if (ApplyChanges(delta, deltadelta)) {
      return true;
    }
  }
  return false;
}

If Increment() returns false, we have exhausted the neighborhood and MakeNextNeighbor() must return false. Sequence() and SetForwardSequence() are two helper methods from the SequenceVarLocalSearchOperator class that allow us to use the ApplyChanges() method to construct the deltas.

And that’s it! Our LocalSearchOperator operator is completed. Let’s test it!

First, we need our LocalSearchOperator:

LocalSearchOperator* const swap_operator =
solver.RevAlloc(new SwapIntervals(all_sequences));

Then we need a complementary DecisionBuilder to construct feasible candidate solutions. We don’t want to spent too much time on the completion of our solutions. We will use the CHOOSE_RANDOM_RANK_FORWARD strategy:

DecisionBuilder* const random_sequence_phase =
               solver.MakePhase(all_sequences,
                                Solver::CHOOSE_RANDOM_RANK_FORWARD);

DecisionBuilder* const complementary_ls_db =
          solver.MakeSolveOnce(solver.Compose(random_sequence_phase,
                                              obj_phase));

If we run the program jobshop_ls1 with our instance problem (file first_example_jssp.txt), we get the optimal solution. Always a good sign. With the instance in abz9 however, we only get a solution with a cost of 1051 in 51,295 seconds:

Time (in s.)	Value	Candidates	Solutions
51,295	1051	31172	26

Not very satisfactory: 1051 is really far from the optimal value of 679. Let’s try to generalize our operator. Instead of just swapping two IntervalVars, we’ll shuffle an arbitrary number of IntervalVars per SequenceVar in the next subsection.

6.7.4. Exchanging an arbitrary number of contiguous IntervalVars on a SequenceVar

You’ll find the code in the file jobshop_ls2.cc and the ShuffleIntervals operator in the file jobshop_ls.

After having implemented the SwapIntervals operator, the only real difficulty that remains is to implement a permutation. This is not an easy task but we’ll elude this difficulty by only exchanging contiguous IntervalVars and using the std::next_permutation() function. You can find the declaration of this function in the header algorithm. Its customizable version reads like:

template <class BidirectionalIterator, class Compare>
bool next_permutation (BidirectionalIterator first,
                       BidirectionalIterator last, Compare comp);

We accept the default values for the BidirectionalIterator and the Compare classes. It will rearrange the elements in the range [first,last) into the next lexicographically greater permutation. An example will clarify this jargon:

No	Permutations
1	0 1 2
2	0 2 1
3	1 0 2
4	1 2 0
5	2 0 1
6	2 1 0

We have generated the permutations of 0,1,2 with std::next_permutation(). There are 3! = 6 permutations (the first permutation is given to std::next_permutation() and is not generated by it) and you can see that the permutations are ordered by value, i.e. 0 1 2 is smaller than 0 2 1 that itself is smaller than 1 0 2, etc[3].

As usual with the std, the last element is not involved in the permutation. There is only one more detail we have to pay attention to. We ask the user to provide the length of the permutation with the gflags flag FLAGS_shuffle_length. First, we have to test if this length makes sense but we also have to adapt it to each SequenceVar variable.

Without delay, we present the constructor of the ShuffleIntervals LocalSearchOperator:

ShuffleIntervals(const std::vector<operations_research::SequenceVar*>& \
                 vars,
                 int max_length) :
  SequenceVarLocalSearchOperator(vars, size),
  max_length_(max_length),
  current_var_(-1),
  current_first_(-1),
  current_length_(-1) {}

vars is the array of SequenceVars. max_length is the length of the sequence of IntervalVars to shuffle. Because you can have less IntervalVars for a given SequenceVar, we have named it max_length.

The indices are very similar to the ones of the SwapIntervals operator:

• current_var_: the index of the processed SequenceVar;
• current_first_: the index of the first IntervalVar variable to shuffle;
• current_length_: the length of the current sub-array of indices to shuffle. It must be smaller or equal to the number of IntervalVars in the SequenceVar.

Here is the code to increment the next permutation:

bool Increment() {
  if (!std::next_permutation(current_permutation_.begin(),
                             current_permutation_.end())) {
    // No permutation anymore -> update indices
    if (++current_first_ >
                      Var(current_var_)->size() - current_length_) {
      if (++current_var_ >= Size()) {
        return false;
      }
      current_first_ = 0;
      current_length_ = std::min(Var(current_var_)->size(),
                                 max_length_);
      current_permutation_.resize(current_length_);
    }
    // Reset first permutation in case we have to increase
    // the permutation.
    for (int i = 0; i < current_length_; ++i) {
      current_permutation_[i] = i;
    }
    // Start with the next permutation, not the identity
    // just constructed.
    if(!std::next_permutation(current_permutation_.begin(),
                              current_permutation_.end())) {
      LOG(FATAL) << "Should never happen!";
    }
  }
  return true;
}

Thanks to the std::next_permutation() function, this is a breeze! The OnStart() method is again straightforward:

virtual void OnStart() {
  current_var_ = 0;
  current_first_ = 0;
  current_length_ = std::min(Var(current_var_)->size(), max_length_);
  current_permutation_.resize(current_length_);
  for (int i = 0; i < current_length_; ++i) {
    current_permutation_[i] = i;
  }
}

We just have to pay attention to resize() the std::vector current_permutation_ of indices and we start with the same permutation: [0, 1, 2, 3, ...].

We again use our template for the MakeNextNeighbor() method:

virtual bool MakeNextNeighbor(Assignment* delta,
                              Assignment* deltadelta) {
  CHECK_NOTNULL(delta);
  while (true) {
    RevertChanges(true);
    if (!Increment()) {
      return false;
    }
    std::vector<int> sequence = Sequence(current_var_);
    std::vector<int> sequence_backup(current_length_);
    for (int i = 0; i < current_length_; ++i) {
      sequence_backup[i] = sequence[i + current_first_];
    }
    for (int i = 0; i < current_length_; ++i) {
      sequence[i + current_first_] =
          sequence_backup[current_permutation_[i]];
    }
    SetForwardSequence(current_var_, sequence);
    if (ApplyChanges(delta, deltadelta)) {
      return true;
    }
  }
  return false;
}

If Increment() returns false, we have exhausted the neighborhood and MakeNextNeighbor() must return false. After the call to Increment(), we simply copy the indices according to the new generated permutation and call the helper method SetForwardSequence() to update the current SequenceVar variable. ApplyChanges() constructs the deltas for us.

File jobshop_ls2.cc is exactly the same as file jobshop_ls1.cc except that we use the ShuffleIntervals operator instead of the SwapIntervals operator.

We again obtain the optimal solution on our instance problem (file first_example_jssp.txt whether shuffle_length=2 or shuffle_length=3). What about the abz9 instance? The next table summarize some tests with different values for the suffle_length parameter:

suffle_length	Time (in s.)	Value	Candidates	Solutions
2	12,301	1016	4302	32
3	21,312	1087	7505	15
4	170,087	1034	70854	33
5	584,173	1055	268478	27

These results are typical for a Local Search operator. There certainly are several lessons to be drawn from these results, but let’s focus on one of the most basic and important ones. The path taken to find the local optimum is crucial. Even if the neighborhoods (theoretically) constructed with suffle_length set to 2 are all contained in the neighborhoods constructed with suffle_length set to 3, we don’t reach the same local optimum. This is very important to understand. The paths taken in both cases are different. The (practical) construction of the neighbourhoods is dynamic and path-dependent. Good (meta-)heuristics are path-aware: these heuristics take the path (and thus the history of the search) into account. Moreover, bigger neighbourhoods (shuffle_length = 3) aren’t necessarily better than smaller ones (shuffle_length = 2). We obtain a better solution quicker with shuffle_length=2 than with suffle_length=3.

The best solution obtained so far has a value of 1016. Can we do better? That’s the topic of next sub-section!

6.7.5. Can we do better?

You’ll find the code in the file jobshop_ls3.cc.

You should know by now that whenever we ask this question in this manual, the answer is yes. To find a better solution, we’ll first investigate how important the initial solution is and then we’ll enlarge our definition of a neighborhood by combining our two LocalSearchOperators.

6.7.5.1. The initial solution

Local search is strongly dependent on the initial solution. Investing time in finding a good solution is a good idea. We’ll use... Local Search to find an initial solution to get the real Local Search started! The idea is that maybe we can find an even better solution in the vicinity of this initial solution. We don’t want to spend too much time to find it though and we’ll limit ourselves to a custom-made SearchLimit. To define this SearchLimit, we construct a callback:

class LSInitialSolLimit : public ResultCallback<bool> {
  public:
    LSInitialSolLimit(Solver * solver, int64 global_time_limit,
                      int solution_nbr_tolerance) :
      solver_(solver), global_time_limit_(global_time_limit),
      solution_nbr_tolerance_(solution_nbr_tolerance),
      time_at_beginning_(solver_->wall_time()),
      solutions_at_beginning_(solver_->solutions()),
      solutions_since_last_check_(0) {}

    //  Returns true if limit is reached, false otherwise.
    virtual bool Run() {
      bool limit_reached = false;

      //  Test if time limit is reached.
      if ((solver_->wall_time() - time_at_beginning_)
                                             > global_time_limit_) {
        limit_reached = true;
        //  Test if we continue despite time limit reached.
        if (solver_->solutions() - solutions_since_last_check_
                                       >= solution_nbr_tolerance_) {
          //  We continue because we produce enough new solutions.
          limit_reached = false;
        }
      }
      solutions_since_last_check_ = solver_->solutions();

      return limit_reached;
    }

  private:
    Solver * solver_;
    int64 global_time_limit_;
    int solution_nbr_tolerance_;

    int64 time_at_beginning_;
    int solutions_at_beginning_;
    int solutions_since_last_check_;
};

The main method in this callback is the virtual bool Run() method. This method returns true if our limit has been reached and false otherwise. The time limit in ms is given by global_time_limit. If the Search is still producing a certain amount solution_nbr_tolerance of solutions, we let the search continue.

To initialize our first Local Search that finds our initial solution, we use the same code as in the file jobshop_ls2.cc (we call this first solution first_solution).

To find an initial solution, we use Local Search and start form the first_solution found. We only use a ShuffleIntervals operator with a shuffle length of 2. This time, we limit this Local Search with our custom limit:

SearchLimit * initial_search_limit = solver.MakeCustomLimit(
              new LSInitialSolLimit(&solver,
                                    FLAGS_initial_time_limit_in_ms,
                                    FLAGS_solutions_nbr_tolerance));

FLAGS_initial_time_limit_in_ms and FLAGS_solutions_nbr_tolerance are the two gflags flags we use in the constructor of the callback LSInitialSolLimit described above to limit the search.

The initial solution is stored in an Assignment initial_solution.

Now, we are ready to prepare the Local Search with our two LocalSearchOperators combined.

6.7.5.2. Combining the two LocalSearchOperators

Often, one LocalSearchOperator isn’t enough to define a good neighborhood. Finding a good definition of a neighborhood is an art and is really difficult. One way to diversify a neighborhood is to combine several basic LocalSearchOperators. Here, we combine SwapIntervals and ShuffleIntervals:

std::vector<LocalSearchOperator*> operators;
LocalSearchOperator* const swap_operator =
          solver.RevAlloc(new SwapIntervals(all_sequences));
operators.push_back(swap_operator);

LocalSearchOperator* const shuffle_operator =
        solver.RevAlloc(new ShuffleIntervals(all_sequences,
                                             FLAGS_shuffle_length));
operators.push_back(shuffle_operator);

LocalSearchOperator* const ls_concat =
                       solver.ConcatenateOperators(operators, true);

The ConcatenateOperators() method takes an std::vector of LocalSearchOperator and a bool that indicates if we want to restart the operators one after the other in the order given by this vector once a solution has been found.

The rest of the code is similar to that in the file jobshop_ls2.cc.

6.7.5.3. Results

If we solve our problem instance (file first_example_jssp.txt), we still get the optimal solution. No surprise here. What about the abz9 instance?

With our default value of

• time_limit_in_ms = 0, thus no time limit;
• shuffle_length = 4;
• initial_time_limit_in_ms = 20000, thus a time of 20 seconds to find an initial solution with Local Search and the ShuffleIntervals operator with a shuffle length of 2 and;
• solutions_nbr_tolerance = 1,

we are not able to improve our best solution so far!

As we said, Local Search is very sensitive to the initial solution chosen. In the next table, we start with different initial solutions:

Initial time limit	Initial sol. obj.	Time	Value	Candidates	Solutions
1,000	1114	81,603	983	49745	35
2,000	1103	103,139	936	70944	59
3,000	1093	104,572	931	70035	60
4,000	1089	102,860	931	68359	60
5,000	1073	84,555	931	63949	60
6,000	1057	42,235	1012	29957	32
7,000	1042	36,935	1012	26515	32
...	...	...	...	...	...
>= 13,000	1016	19,229	1016	13017	32

The first column lists the times allowed to find the initial solution with the ShuffleIntervals operator (with its shuffle length set to 2) and the second column collects the objective values of this initial solution. The more time given to the first Local Search, the better the objective values. The next four columns are the same as before.

You might think that starting from a better solution would give better results but it is no necessarily the case. Our best result, 931 is obtained when we start from solutions with an average objective value. When we start with better solutions, like the one with an objective value of 1016, we completely miss the 931 solution!

This 931 solution seems to be a local optimum for our Local Search and it seems we can not escape it. In the chapter Meta-heuristics: several previous problems, we’ll see how some meta-heuristics escape this local minimum. For now, we turn our attention to another preoccupation: if you read the Candidates column and compare it with the Solutions column, you can see that our algorithm produces lots of candidates and very few solutions. This is normal. Remember that every time a candidate (a neighbor) is produced, the CP solver takes the time to verify if this candidate is a feasible solution. This is costly. In the next section, we’ll see a mechanism to shortcut this verification and command the solver to disregard some candidates without the need for the solver to test them explicitly.

Footnotes

[1]	But not how to implement incremental LocalSearchOperators. This is a very advanced feature.

[2]	Remember our footnote about the PathOperator classes? They use a special customized MakeNeighbor() method instead of the MakeOneNeighbor() method.

[3]

This explanation is not rigorous but it is simple and you can fill the gaps. What happens if you start with 1 0 2? The std::next_permutation() function simply “returns” 1 2 0 (oops, there goes our rigour again!). If you give it 2 1 0, this function returns false but there is a side effect as the array will be ordered! Thus in our case, we’ll get 0 1 2!