# Optimization formulation has slow performance

I am formulating and solving a scheduling problem. The problem consists of scheduling items on a single machine, where the only time element is the transition from item $$i$$ to item $$j$$. For example, it takes 3 units of time to go from item 1 to item 2, and 2 units of time to go from item 2 to item 1. The formulation I follow is using a binary variable $$y_{ij}$$ to represent item $$i$$ immediately preceding item $$j$$.

The constraints and variables are as follows:

Variables: $$y_{ij}, \forall i \in \text{Items},j\in\text{Items},i \neq j\\ s_i, \forall i \in \text{Items}$$ Constraints:
An item will only be behind one other item $$\sum_{j\neq i}y_{ij} \le 1, \forall i \in \text{Items}$$ An item will only be in front of one other item $$\sum_{i\neq j}y_{ij} \le 1, \forall j \in \text{Items}$$ An item cannot be both behind and in front of another item $$y_{ij} + y_{ji} \le 1, \forall i,j \in \text{Items}$$ Only one item will not be before another item $$\sum_{i,j,i \neq j} y_{ij} \ge |\text{Items}| - 1$$ Tracking start times $$s_i - s_j + t_{ij} \le M ( 1 - y_{ij} )$$ Objective: $$\min \sum_{i} s_i$$

Using Pulp with the default solver, I can solve an instances of up to 7 items. I have tried 10 but killed it after waiting >2 minutes. Any feedback on improving the formulation would be greatly appreciated.

• Isn't this the traveling salesman problem, where each item is a city? – RobPratt Mar 20 at 22:41
• Yes! But Im looking to formulate it along the lines of a scheduling formulation. While this is the traveling salesman problem, I was looking to add more scheduling elements, such as tardiness and more machines. I think these could still be formulated as a form of the traveling salesman, but I was under the impression other, more efficient formulations exist. – Bob Jeans Mar 21 at 3:18

## 3 Answers

since you deal with scheduling on top of MIP you could try CPOptimizer scheduling.

For instance you could start with https://github.com/AlexFleischerParis/howtowithopl/blob/master/tspcpo.mod

using CP;
int     n       = ...;
range   Cities  = 1..n;

int realCity[i in 1..n+1]=(i<=n)?i:1;

// Edges -- sparse set
tuple       edge        {int i; int j;}
setof(edge) Edges       = {<i,j> | ordered i,j in 1..n};
setof(edge) Edges2       = {<i,j> | i,j in 1..n+1};  // node n+1 is node 1

int         dist[Edges] = ...;
int         dist2[<i,j> in Edges2]=(realCity[i]==realCity[j])?0:
((realCity[i]<realCity[j])?dist[<realCity[i],realCity[j]>]:dist[<realCity[j],realCity[i]>]);

dvar interval itvs[1..n+1] size 1;

dvar sequence seq in all(i in 1..n+1) itvs[i];

execute
{

cp.param.TimeLimit=60;
var f = cp.factory;
cp.setSearchPhases(f.searchPhase(seq));
}

tuple triplet { int c1; int c2; int d; };
{triplet} Dist = {
<i-1,j-1,dist2[<i ,j >]>
|  i,j in 1..n+1};

minimize endOf(itvs[n+1]) - (n+1);
subject to
{
startOf(itvs[1])==0; // break sym
noOverlap(seq,Dist,true);   // nooverlap with a distance matrix
last(seq, itvs[n+1]); // last node
}

int  x[<i,j> in Edges]=prev(seq,itvs[i],itvs[j])+prev(seq,itvs[j],itvs[i]);
int  isPrevFromNPlus1[i in 1..n]=prev(seq,itvs[i],itvs[n+1]);
int l=first({i | i in 1..n : isPrevFromNPlus1[i]==1});
edge el=<1,l>;

execute
{
isPrevFromNPlus1;
x;
x[el]=1;
}

// Let us check here that the constraints of the IP model are ok
assert forall (j in Cities)
as:sum (<i,j> in Edges) x[<i,j>] + sum (<j,k> in Edges) x[<j,k>] == 2;

// Let us compute here the objective the IP way
int cost=sum (<i,j> in Edges) dist[<i,j>]*x[<i,j>];

execute
{
writeln(cost);
}


which gives

In

noOverlap(seq,Dist,true);


Dist is the transition time matrix

NB: CPOptimizer also available in python through docplex.

Such a problem is difficult to model and solve following a MILP approach, as you observed. Boolean modeling approaches are tedious to write and don't scale well for this kind of problem.

Your problem can be modeled compactly by following a list-based modeling approach instead of the classical Boolean modeling approach, as you described in your question. This is a modeling approach offered by LocalSolver, which is different from traditional MILP solvers. Note that LocalSolver is commercial software. Nevertheless, it is free for faculty and students.

Since you're interested in scheduling the jobs on several machines, and not only one as described in your question, you may be interested in the job shop scheduling problem, more than in the traveling salesman problem. Below is the code snippet to model the job shop scheduling problem by using the LocalSolver modeling language, namely LSP:

function model() {
// Integer decisions: start time of each activity
// start[j][m] is the start time of the activity of job j which is processed on machine m
start[j in 0..nbJobs-1][m in 0..nbMachines-1] <- int(0, maxStart);
end[j in 0..nbJobs-1][m in 0..nbMachines-1] <- start[j][m] + processingTime[j][m];

// Precedence constraints between the activities of a job
for [j in 0..nbJobs-1][k in 1..nbMachines-1]
constraint start[j][machineOrder[j][k]] >= end[j][machineOrder[j][k-1]];

// Sequence of activities on each machine
jobsOrder[m in 0..nbMachines-1] <- list(nbJobs);

for [m in 0..nbMachines-1] {
// Each job has an activity scheduled on each machine
constraint count(jobsOrder[m]) == nbJobs;

// Disjunctive resource constraints between the activities on a machine
constraint and(0..nbJobs-2, i => start[jobsOrder[m][i+1]][m] >= end[jobsOrder[m][i]][m]);
}

// Minimize the makespan: end of the last activity of the last job
makespan <- max[j in 0..nbJobs-1] (end[j][machineOrder[j][nbMachines-1]]);
minimize makespan;
}


You can write the same model in Python, Java, C#, or C++. Just have a look at the job shop scheduling example here.

Having followed this modeling approach, LocalSolver finds quality solutions in minutes on a standard computer for instances with thousands of jobs to schedule. LocalSolver embeds and combines neighborhood search, constraint programming, and mixed-integer linear programming techniques under the hood to get such results.

You did not define all your symbols, so I am going to infer from the model that $$t_{ij}$$ is the transition time from item $$i$$ to item $$j$$ (a parameter) and $$s_i$$ is the time that item $$i$$ is processed (a variable). Assuming that is correct, your objective function looks unlikely to be correct. The sum of the start/end times has no intrinsic meaning. It is more likely that you want to minimize the makespan (the time it takes to process all items). Normally that is accomplished by adding a continuous variable $$z$$ for the makespan, minimizing $$z$$, and defining $$z$$ via the constraints $$z \ge s_i\quad\forall i.$$ In your case, however, since there is no processing time (only setup time), we can make things even simpler. The makespan will be the sum of the setup times, so you can eliminate the $$s$$ variables and the big-M constraints and just minimize $$\sum_{i,j} t_{ij} y_{ij}$$ subject to your first four constraints. Eliminating the big-M constraints should speed things up a fair bit.

• Thank you for the fast response! And yes you are correct, t_ij is the time to go from item_i to item_j, and s_i is the starting time for item_i and is a continuous variable. Im working on understanding and implementing scheduling formulations, so the intent would be to have processing time in the future, and because of that I want to keep the bigM formulation. I am trying to minimize the makespan, and minimizing your variable z should be the same as minimizing the sum of s_i? Would that improve my performance? – Bob Jeans Mar 20 at 22:00
• No, minimizing $z$ would minimize the makespan, which in your case is the largest $s_i$, not the sum of the $s_i$. The sum of the $s_i$, as I said, has no useful interpretation. For instance, if you have three items processed at times 0, 3 and 7, the makespan (with no processing times) would be 7; the sum of the $s_i$ would be 10. – prubin Mar 21 at 23:06