## Wednesday, August 29, 2018

### Scheduling: easy MIP

Many scheduling problems cause headaches for MIP solvers. Here is one that is easy. The problem is from [1]:

We need to schedule tasks. Each task must be assigned to a single time slot. Not every time slot is allowed: each task has a set of time slots where it can be executed. An example small data set is:

task_availability_map = {
"T1" : [0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
"T2" : [0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
"T3" : [0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
"T4" : [0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
"T5" : [0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
"T6" : [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
"T7" : [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0],
"T8" : [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0],
"T9" : [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
"T10": [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0]
}


In addition we have a maximum of tasks that can be assigned to a single time slot (i,e, a capacity). The goal is to minimize the number of slots we need, such that all tasks are scheduled. If we can handle 3 tasks per time slot, an optimal schedule can look like:

                t6          t8          t9         t12



This solution is by no means unique.

#### MIP Model

To model this as a MIP model we can use two binary variables:$x_{i,t} = \begin{cases} 1 & \text{ if job i is assigned to time slot t}\\ 0 & \text{ otherwise}\end{cases}$ and $y_t = \begin{cases} 1 & \text{if time slot t has at least one job assigned to it}\\ 0 & \text{otherwise}\end{cases}$

A MIP model can look like: \bbox[lightcyan,10px,border:3px solid darkblue]{\begin{align}\min& \sum_t y_t\\ & \sum_t \mathit{ok}_{i,t} x_{i,t} = 1 && \forall i \\ & \sum_i \mathit{ok}_{i,t} x_{i,t} \le N && \forall t \\ & y_t \ge x_{i,t} &&\forall i,t|\mathit{ok}_{i,t}=1\\ & x_{i,t}, y_t \in \{0,1\}\end{align}}

We can re-interpret $$\mathit{ok}_{i,t}$$ as $$\mathit{ok}_{i,t}=\mathit{true}$$ if and only if $$(i,t)$$ is an allowed assignment. Then the model can be:  \bbox[lightcyan,10px,border:3px solid darkblue]{\begin{align}\min& \sum_t y_t\\ & \sum_{t| \mathit{ok}_{i,t}} x_{i,t} = 1 && \forall i \\ & \sum_{i| \mathit{ok}_{i,t}} x_{i,t} \le N && \forall t \\ & y_t \ge x_{i,t} &&\forall i,t|\mathit{ok}_{i,t}\\ & x_{i,t}, y_t \in \{0,1\}\end{align}}

There are a number of interesting issues. The first one is that we have many $$(i,t)$$ combinations that we know are not allowed. Instead of fixing those $$x_{i,t}=0$$ we can try to not even generate those variables. For this example we would reduce the number of variables from $$10 \times 17 = 170$$ to $$50$$ (the number of ones in the data table).

A GAMS model can look like:

I used the set ok to indicate if a combination $$(i,t)$$ is allowed. This makes the equations a bit simpler. Equations are the most difficult part in a model to debug, so we should go the extra mile to simplify them.

Using a time slot capacity of 3, the results can look like:

----     64 VARIABLE x.L  assignment

t6          t8          t9         t12

----     64 VARIABLE y.L  slot is used

t6  1,    t8  1,    t9  1,    t12 1

----     64 VARIABLE z.L                   =            4  objective variable


We can see we really only have 50 $$x$$ variables by looking at the GAMS column listing:

---- x  assignment

(.LO, .L, .UP, .M = 0, 0, 1, 0)
1       assign2(t4)

(.LO, .L, .UP, .M = 0, 0, 1, 0)
1       assign2(t5)

(.LO, .L, .UP, .M = 0, 0, 1, 0)
1       assign2(t6)

REMAINING 47 ENTRIES SKIPPED


There is one more issue. Some of the columns in the data matrix have only zero entries. This means some of the assign2 constraints are empty. To be more precise they are of the form $$0 \le 3$$. These are dropped by GAMS (unless they are infeasible, e.g. $$0\le -1$$, in which case GAMS would complain). You can see this in the equation listing (after increasing the limrow option):

---- assign2  =L=  assignment constraint

assign2(t1)..  x(task6,t1) =L= 3 ; (LHS = 0)

assign2(t2)..  x(task6,t2) =L= 3 ; (LHS = 0)

assign2(t3)..  x(task6,t3) =L= 3 ; (LHS = 0)

+ x(task10,t7) =L= 3 ; (LHS = 0)

(LHS = 0)

assign2(t15)..  x(task7,t15) =L= 3 ; (LHS = 0)


The equations assign2(t16) and assign(t17) are dropped.

Of course, if we want, we can make this more explicit by saying something like:

assign2(t)$sum(ok(i,t),1).. sum(ok(i,t), x(i,t)) =l= N; Here the equation is not generated if the column sum is zero. As this does not add much value, it may be better to just let GAMS take care of this. We can generalize this a little bit. The constraint assign2 is not needed if the number of entries in column $$t$$ of the parameter $$\mathit{allowed}_{i,t}$$ is less or equal to $$N$$. We can do something like: parameter cnt(t) 'count nonzero entries in column t'; cnt(t) = sum(i, allowed(i,t)); assign2(t)$(cnt(t)>N)..
sum(ok(i,t), x(i,t)) =l= N;

This saves us relatively few constraints. It is much better to focus on $$x$$ as that is 2-dimensional so we have more bang for the buck. I'll leave it to the solver's presolver to remove the constraints.

#### Binary or continuous?

If we want we can relax the binary variable $$y_t$$ to a continuous variable between 0 and 1. Intuitively that is a good idea (fewer integer variables). However note that some solvers prefer them to be binary, and they change them back:

Tried aggregator 2 times.
MIP Presolve eliminated 16 rows and 9 columns.
Aggregator did 1 substitutions.
Reduced MIP has 59 rows, 58 columns, and 163 nonzeros.
Reduced MIP has 47 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.06 sec. (0.16 ticks)
Probing time = 0.00 sec. (0.05 ticks)
Tried aggregator 1 time.
Reduced MIP has 59 rows, 58 columns, and 163 nonzeros.
Reduced MIP has 58 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.15 ticks)


#### MIP or CP (Mixed Integer Programming vs Constraint Programming)?

In the comments it was argued that a constraint programming approach is more suited for this model.  There are 2 dimensions to this argument: (1) is the problem easier to formulate as a CP problem or (2) is the problem easier to solve.

I think these two are linked. My rule of thumb is: if the problem is easy to formulate as a MIP keep it a MIP. If we need to use a lot of tricks and extra variables, try a CP formulation. An example where CP is easier is [2]. If the problem is easier to model, there is a good chance that the model also solves faster.

The poster mentioned the real problem has a few hundred tasks and time slots. Let's create some  random data for this: 200 tasks, 200 time slots, time slot capacity = 5. Each task has two random  values: a start value $$\sim U(1,200)$$ and a length $$\sim U(3,10)$$. Here $$U$$ indicates we draw from the uniform distribution. This means: we can assign a task to time slots $$\mathit{start},\dots,\mathit{start}+\mathit{length}-1$$. This gives a problem with 1,071 $$x$$ variables. It solves very quickly: less than a second using 1,022 simplex iterations and all work done in node 0. The GAMS model execution (and generation) time was about 0.1 seconds.

I don't see a good reason to consider a CP approach: the MIP model is simple and the model solves quickly.

#### Pulp?

In the post it was reported that the implementation in Pulp took a long time. I believe that when we have a properly formulated mathematical model, it should not take too much time to develop Python/Pulp code to implement this model. It can look like:

from pulp import *

#
# data
#
"T1" : [0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
"T2" : [0, 0, 0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
"T3" : [0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
"T4" : [0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
"T5" : [0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
"T6" : [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
"T7" : [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0],
"T8" : [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0],
"T9" : [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
"T10": [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0]
}

tasks = ["T%i" % i for i in range(1,11)]
slots = range(17)

N = 3 # max number of tasks in a slot

#
# model
#

x = LpVariable.dicts(
"assign",
[(i, t) for i in tasks for t in slots if task_availability_map[i][t]==1],
lowBound=0,upBound=1,cat="Integer")

y = LpVariable.dicts(
"use",slots,lowBound=0,upBound=1,cat="Integer")

prob = LpProblem("Easy MIP",LpMinimize)

prob += lpSum([y[t] for t in slots]), "Number of slots in use"

prob += lpSum([x[(i,t)] for t in slots if task_availability_map[i][t]==1]) == 1, "Assign1(%s)" % i

for t in slots:
prob += lpSum([x[(i,t)] for i in tasks if task_availability_map[i][t]==1]) <= N, "Assign2(%i)" % t

for t in slots:
prob += y[t] >= x[(i,t)], "Bound(%s,%i)" % (i,t)

#
# for debugging
#
prob.writeLP("d:\\tmp\\easy.lp")

#prob.solve()
prob.solve(PULP_CBC_CMD(msg=1))
# no solve log when running under Jupyter

print "Status: %s" % LpStatus[prob.status]
print "Number of slots needed: %.0f" % value(prob.objective)
print "Assignment:"
for t in slots:
if x[(i,t)].value() > 0.5:
print "Task %s -> slot %i" % (i,t)


The print out looks like:

Status: Optimal
Number of slots needed: 4
Assignment:


Some notes.

• The Pulp model is more unwieldy than the GAMS model due to all the explicit looping.
• To debug a Pulp model it is a good idea to print the LP file. This essentially replaces the GAMS equation and column listing we have shown earlier. We will look at this below.
• We make multiple passes over the data. If this is a real large data set we need to develop additional data structures to quickly run through the task_availability_map both row-wise and column-wise, while skipping the zero values.
• Pulp does not show the solver log, especially when running under Jupyter. This makes it important to always print the status. It would be nice if the solver log is captured and echoed to the Jupyter output cell.
• There are no built-in tools to print solution tables in an easy to read format. You need to develop your own printing solutions.
• When printing the solution we should not compare values of binary variables to exactly 1.0. For instance in [2] we see:
print "The choosen tables are out of a total of %s:"%len(possible_tables)
for table in possible_tables:
if x[table].value() == 1.0:
print table
This is not safe. A binary variable can have values like 0.999999 or 1.000001 due to tolerances. My code is a bit wild:
if x[(i,t)].value() > 0.5:
but it is as safe as can be.
• When we plug in the data set with 200 tasks by 200 slots we generated with GAMS earlier, and use a capacity of 5 tasks per slot, we can generate and solve the problem in 4 seconds.
• The performance of Pulp/CBC is somewhat erratic. I'll report on this below.
• In GAMS a variable declaration does not allocate memory for variables. Hence we can just do binary variable x(i,t). In Pulp we want to allocate only variables that are used, so the LpVariable.dicts statement needs to know which variables are actually being used.

It is interesting to see how Pulp deals with the empty columns in the data. They lead to constraints like $$0 \le 3$$. When looking at the LP file, we see:

Assign2(0): assign_('T6',_0) <= 3
Assign2(1): assign_('T6',_1) <= 3
Assign2(10): assign_('T6',_10) + assign_('T7',_10) + assign_('T8',_10) <= 3
Assign2(11): assign_('T6',_11) + assign_('T7',_11) + assign_('T8',_11) <= 3
Assign2(12): assign_('T6',_12) + assign_('T7',_12) + assign_('T8',_12) <= 3
Assign2(13): assign_('T6',_13) + assign_('T7',_13) + assign_('T8',_13) <= 3
Assign2(14): assign_('T7',_14) <= 3
_dummy: __dummy = 0
Assign2(15): __dummy <= 3
Assign2(16): __dummy <= 3
Assign2(2): assign_('T6',_2) <= 3
Assign2(3): assign_('T1',_3) + assign_('T3',_3) + assign_('T6',_3) <= 3
Assign2(4): assign_('T1',_4) + assign_('T3',_4) + assign_('T5',_4)
+ assign_('T6',_4) <= 3
Assign2(5): assign_('T1',_5) + assign_('T2',_5) + assign_('T3',_5)
+ assign_('T4',_5) + assign_('T5',_5) + assign_('T6',_5) <= 3
Assign2(6): assign_('T1',_6) + assign_('T10',_6) + assign_('T2',_6)
+ assign_('T3',_6) + assign_('T4',_6) + assign_('T5',_6) + assign_('T6',_6)
+ assign_('T9',_6) <= 3
Assign2(7): assign_('T1',_7) + assign_('T10',_7) + assign_('T3',_7)
+ assign_('T4',_7) + assign_('T5',_7) + assign_('T6',_7) + assign_('T9',_7)
<= 3
Assign2(8): assign_('T10',_8) + assign_('T6',_8) + assign_('T9',_8) <= 3
Assign2(9): assign_('T10',_9) + assign_('T6',_9) + assign_('T9',_9) <= 3

Bounds
__dummy = 0


We see these two $$0 \le 3$$ constraints are not dropped, but rather included with a dummy variable __dummy. I am not sure why an empty constraint like $$0 \le 3$$ is added to the model and passed on to the solver. I am also surprised about __dummy being to zero in two different ways: both bounds are used to fix __dummy to zero and additionally a constraint is included to achieve the same thing. For most solvers one of these measures should be enough. In practice, of course, this is not really a big problem: the presolver will take care of this, and will get rid of the equation _dummy: __dummy = 0.

We can prevent this business with the dummy variable and equation by making sure Pulp does not generate an equation when the left-hand side is empty:

tused = [False for t in slots]
for t in slots:
tused[t] = True

for t in slots:
if tused[t]:
prob += lpSum([x[(i,t)] for i in tasks if task_availability_map[i][t]==1]) <= N, "Assign2(%i)" % t


#### Printing

One thing I miss when doing modeling in programming languages like Python is an easy way to print results in a nice readable format.

GAMS can print any multi-dimensional parameter or variable using the display statement. Usually the results looks quite reasonable. Here is the output of the statement display allowed, x.l;

----     70 PARAMETER allowed

t1          t2          t3          t4          t5          t6          t7          t8          t9

task1                                            1.000       1.000       1.000       1.000       1.000
task3                                            1.000       1.000       1.000       1.000       1.000
task6        1.000       1.000       1.000       1.000       1.000       1.000       1.000       1.000       1.000

+         t10         t11         t12         t13         t14         t15

task6        1.000       1.000       1.000       1.000       1.000
task7                    1.000       1.000       1.000       1.000       1.000

----     70 VARIABLE x.L  assignment

t5          t6          t7         t13



I have heard someone saying that the display statement is the most important statement in the GAMS language. Especially during model development when I use small data sets (my favorite), I use the display statement a lot. When the data sets are large I usually use the pivot table facility in the GDX viewer, When the model is beyond its early development stage, I actually pay quite some time and effort to produce meaningful reports, often in Excel. But even when debugging problems in later stages, the display statement remains a very important tool.

Python/Pulp has no good standard printing tool for multidimensional data. We can do something simple, but the output is less clear:

 View results in Jupyter
A slightly easier way is to use: [(i,t) for (i,t) in x.keys() if x[(i,t)].value()>0.5].

Things like Panda data frames look quite nice in Jupyter, but typically we need to populate them first.

When developing and debugging models, often under pressure, life without easy data viewing is just more complicated than it should be.

#### Solver Performance

The performance of Pulp/CBC was not very predictable. Here are some results using the 200 tasks, 200 time slots random data set.

N=3 GAMS Cplex 1 68 2788 30 1
GAMS Cplex 4 68 2937 42 1.9
GAMS CBC 1 68 604336 18911 146
GAMS CBC 4 68 230016 13323 72
Pulp CBC > 2 hours
N=4 GAMS Cplex 1 54 2148 0 1
GAMS Cplex 4 54 2139 0 2.3
GAMS CBC 1 54 273966 11500 97
GAMS CBC 4 54 100467 2745 27
Pulp CBC 54 374219 13002 192
N=5 GAMS Cplex 1 49 1022 0 0.4
GAMS Cplex 4 49 1022 0 1
GAMS CBC 1 49 543 1 2.1
GAMS CBC 4 49 543 1 2.5
Pulp CBC 49 ? 0 4.5

The $$N=3$$ problem with Pulp went nowhere. The last few lines of the CBC log look like:

Cbc0010I After 690000 nodes, 17543 on tree, 69 best solution, best possible 68 (7140.39 seconds)
Cbc0010I After 691000 nodes, 17542 on tree, 69 best solution, best possible 68 (7149.65 seconds)
Cbc0010I After 692000 nodes, 17543 on tree, 69 best solution, best possible 68 (7159.53 seconds)
Cbc0010I After 693000 nodes, 17535 on tree, 69 best solution, best possible 68 (7167.80 seconds)
Cbc0010I After 694000 nodes, 17510 on tree, 69 best solution, best possible 68 (7176.48 seconds)
Cbc0010I After 695000 nodes, 17494 on tree, 69 best solution, best possible 68 (7183.40 seconds)
Cbc0010I After 696000 nodes, 17592 on tree, 69 best solution, best possible 68 (7194.70 seconds)
Cbc0010I After 697000 nodes, 17579 on tree, 69 best solution, best possible 68 (7203.96 seconds)
Cbc0010I After 698000 nodes, 17571 on tree, 69 best solution, best possible 68 (7209.88 seconds)
Cbc0010I After 699000 nodes, 17566 on tree, 69 best solution, best possible 68 (7215.16 seconds)
Cbc0027I Exiting on user event
Cbc0005I Partial search - best objective 69 (best possible 68), took 9050076 iterations and 699207 nodes (7216.40 seconds)
Cbc0032I Strong branching done 993710 times (8577199 iterations), fathomed 197962 nodes and fixed 223209 variables
Cbc0035I Maximum depth 101, 124125 variables fixed on reduced cost
Cuts at root node changed objective from 67 to 68
Probing was tried 10 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.103 seconds)
Gomory was tried 117378 times and created 50461 cuts of which 95 were active after adding rounds of cuts (286.044 seconds)
Knapsack was tried 10 times and created 2 cuts of which 0 were active after adding rounds of cuts (0.122 seconds)
Clique was tried 10 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.013 seconds)
MixedIntegerRounding2 was tried 10 times and created 3 cuts of which 0 were active after adding rounds of cuts (0.052 seconds)
FlowCover was tried 10 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.006 seconds)
TwoMirCuts was tried 117378 times and created 612 cuts of which 0 were active after adding rounds of cuts (234.192 seconds)

Result - User ctrl-cuser ctrl-c

Objective value:                69.00000000
Lower bound:                    68.000
Gap:                            0.01
Enumerated nodes:               699207
Total iterations:               9050076
Time (CPU seconds):             7216.48
Time (Wallclock seconds):       7216.48


There remain some questions about what Pulp/CBC is doing on some versions of this model.
Notes:

• Both the GAMS version and the Pulp version of CBC were based on version 2.9 of the CBC library. The Pulp version says: Build Date: Feb 12 2015. Looks like we deal with a somewhat older executable.
• Pulp uses ascii LP files to communicate the model with the solver. GAMS uses binary files. This will cause the problem to be slightly different.
• CBC needs to be properly compiled to support multiple threads. Not all CBC versions on Windows have been built with multi-thread support.
• CBC failing on this model is somewhat scary.
• With slightly different options: prob.solve(PULP_CBC_CMD(threads=1, mip=1, options=['sec','3600'],msg=1)), I saw CBC actually just hang in the middle of the run.

#### Conclusions

• Even small, simple and easy models have interesting angles. In this case we have a decision variable $$x_{i,t}$$ of which we only want to use a subset, and we may encounter some empty constraints.
• Pulp/CBC has some issues on some data sets.
• Under the hood Pulp and GAMS are different in how empty constraints are handled.

This model was more interesting than I initially thought.

#### References

1. Scheduling optimization to minimize the number of timeslots (with constraints), https://stackoverflow.com/questions/52022575/scheduling-optimization-to-minimize-the-number-of-timeslots-with-constraints
2. https://www.coin-or.org/PuLP/CaseStudies/a_set_partitioning_problem.html

## Saturday, August 25, 2018

### Best Factorization X=AB

Restated from an originally posting in [1]:

Given a matrix $$X$$ find a factorization $$X\approx AB$$ such that all elements of $$A,B$$ are between 0 and 1 and the matrix product $$AB$$ is as close as possible to $$X$$. Can we solve this as an NLP (nonlinear programming) model?

The model can simply look like: \begin{align} \min\> & ||AB-X||^2 \\ & 0 \le A,B \le 1 \end{align} or \begin{align} \min &\sum_{i,j} \left ( \sum_k a_{i,k} b_{k,j} - x_{i,j} \right )^2\\ & 0 \le a_{i,k},b_{k,j} \le 1\end{align} However, this is a rather nasty non-convex problem. As I expected to find local optima, I solved this in a loop with a different random starting points for $$A$$ and $$B$$ (each uniform between 0 and 1). With a small data set \begin{align}&i=1,\dots,8\\ &j=1,\dots,12\\ & k=1,\dots,6\end{align} and random values for $$X$$ (uniform between 0 and 6), we see:

Results with NLP solver CONOPT

----     59 PARAMETER results

obj   modelstat        time        best

iter1      129.370    localopt       0.062     129.370
iter2      129.088    localopt       0.094     129.088
iter3      126.113    localopt       0.047     126.113
iter4      127.768    localopt       0.031
iter5      128.103    localopt       0.046
iter6      128.178    localopt       0.094
iter7      126.644    localopt       0.031
iter8      129.179    localopt       0.032
iter9      128.201    localopt       0.047
iter10     130.898    localopt       0.125
iter11     125.803    localopt       0.063     125.803
iter12     127.407    localopt       0.063
iter13     127.898    localopt       0.031
iter14     127.768    localopt       0.047
iter15     129.823    localopt       0.031
iter16     125.803    localopt       0.047
iter17     129.823    localopt       0.063
iter18     128.901    localopt       0.047
iter19     126.105    localopt       0.047
iter20     128.164    localopt       0.032

Results with NLP solver IPOPT

----     59 PARAMETER results

obj   modelstat        time        best

iter1      126.113    localopt       0.344     126.113
iter2      128.129    localopt       0.266
iter3      126.113    localopt       0.313
iter4      128.103    localopt       0.218
iter5      128.901    localopt       0.266
iter6      125.803    localopt       0.344     125.803
iter7      129.330    localopt       0.375
iter8      128.816    localopt       0.360
iter9      128.103    localopt       0.406
iter10     130.276    localopt       0.390
iter11     127.255    localopt       0.328
iter12     126.113    localopt       0.281
iter13     129.002    localopt       0.297
iter14     128.149    localopt       0.235
iter15     129.823    localopt       0.266
iter16     126.545    localopt       0.328
iter17     127.094    localopt       0.204
iter18     128.164    localopt       0.344
iter19     128.201    localopt       0.313
iter20     128.723    localopt       0.281


Both solvers find a solution with an objective of 125.803. This is a small problem with $$8 \times 6 + 6 \times 12 = 120$$ nonlinear variables. This looks doable for a global solver like Baron or Antigone. The problem is that we have $$8 \times 6\times 12 = 576$$ bilinear terms $$a_{i,k}b_{k,j}$$. This will make life really difficult for the global solvers. Let's try this, while using the 125.803 solution as an initial point and with a time limit of 1000 seconds.

Results with BARON

Iteration    Open nodes         Time (s)    Lower bound      Upper bound
1             1             3.00     54.4268          125.803
*         4             3             8.00     54.4268          125.715
41            28            38.00     57.2903          125.715
76            54            69.00     58.6781          125.715
115            82           100.00     59.7597          125.715
155           111           131.00     60.2728          125.715
190           137           162.00     60.5990          125.715
229           164           192.00     60.8507          125.715
266           192           223.00     61.1276          125.715
301           216           254.00     61.3283          125.715
335           239           284.00     61.4689          125.715
373           265           315.00     61.7422          125.715
414           296           346.00     61.9213          125.715
454           321           377.00     62.1128          125.715
492           348           408.00     62.3822          125.715
528           370           438.00     62.4439          125.715
567           392           469.00     62.5986          125.715
598           412           499.00     62.6449          125.715
631           433           531.00     62.6994          125.715
664           458           561.00     62.7280          125.715
695           478           593.00     62.7932          125.715
728           501           624.00     62.8952          125.715
754           520           655.00     62.9707          125.715
782           539           685.00     63.0497          125.715
809           558           716.00     63.1053          125.715
834           577           746.00     63.1215          125.715
866           599           777.00     63.1805          125.715
894           620           808.00     63.1963          125.715
920           636           839.00     63.3036          125.715
949           655           869.00     63.3194          125.715
978           677           900.00     63.4053          125.715
1007           695           930.00     63.4625          125.715
1036           717           961.00     63.5189          125.715
1066           737           992.00     63.5608          125.715
1096           758          1023.00     63.6598          125.715
1122           778          1054.00     63.7217          125.715
1149           797          1085.00     63.8112          125.715
1160           804          1097.00     63.8213          125.715

*** Max. allowable time exceeded ***

Results with Antigone

-------------------------------------------------------------------------------
Time (s) Nodes explored Nodes remaining Best possible   Best found Relative Gap
-------------------------------------------------------------------------------

Searching for feasible solutions with 4 starting points at tree level 0 --
25              1               1    +9.689e+01   +1.258e+02   +2.298e-01
Adding 0 total cutting planes at tree level 0 ----------------------------
Adding 739 total cutting planes at tree level 0 --------------------------
Strong Branching at tree level 1 -----------------------------------------
Strong Branching at tree level 1 -----------------------------------------
561              1               2    +9.893e+01   +1.257e+02   +2.131e-01
Adding 0 total cutting planes at tree level 1 ----------------------------
Generating Edge-Concave Cuts; 0 generated so far -------------------------
Strong Branching at tree level 2 -----------------------------------------
Strong Branching at tree level 2 -----------------------------------------
684              2               3    +9.893e+01   +1.257e+02   +2.131e-01
Adding 0 total cutting planes at tree level 1 ----------------------------
Strong Branching at tree level 2 -----------------------------------------
Strong Branching at tree level 2 -----------------------------------------
Strong Branching at tree level 2 -----------------------------------------
804              3               4    +9.894e+01   +1.257e+02   +2.130e-01
Adding 0 total cutting planes at tree level 2 ----------------------------
Strong Branching at tree level 3 -----------------------------------------
Strong Branching at tree level 3 -----------------------------------------
905              4               5    +9.894e+01   +1.257e+02   +2.130e-01
921              5               6    +9.897e+01   +1.257e+02   +2.128e-01
Adding 0 total cutting planes at tree level 3 ----------------------------
962              6               7    +9.897e+01   +1.257e+02   +2.128e-01
970              7               8    +9.897e+01   +1.257e+02   +2.128e-01
979              8               9    +9.897e+01   +1.257e+02   +2.128e-01
986              9              10    +9.897e+01   +1.257e+02   +2.128e-01
993             10              11    +9.897e+01   +1.257e+02   +2.128e-01

Reached time limit of 1000 CPU s
1000             10              11    +9.897e+01   +1.257e+02   +2.128e-01


Both solvers have troubles closing the gap. Antigone seems to do better initially, but rest assured: it will slow down considerably and bounds are just not getting closer reasonably fast. But also: both solvers are able to improve the solution from 125.803 to 125.715. Interesting results.

We can make things less non-linear by considering absolute deviations:\begin{align}\min & \sum_{i,j} \left| \sum_k a_{i,k} b_{k,j} - x_{i,j}\right |\\ & 0 \le a_{i,k},b_{k,j} \le 1 \end{align} or \begin{align} \min\> & \sum_{i,j} y_{i,j} \\ & -y_{i,j} \le \sum_k a_{i,k} b_{k,j} - x_{i,j} \le y_{i,j} \\ & 0 \le a_{i,k},b_{k,j} \le 1 \\ & y_{i,j} \ge 0 \end{align} Note: when implementing this we should make sure not to duplicate the expression $$\sum_k a_{i,k} b_{k,j} - x_{i,j}$$. That means: introduce another set of variables.  An alternative least absolute deviations formulation uses variable splitting and can look like: \begin{align} \min\> & \sum_{i,j} \left ( y^+_{i,j}+y^-_{i,j}\right) \\ & y^+_{i,j}- y^-_{i,j} = \sum_k a_{i,k} b_{k,j} - x_{i,j}\\ & 0 \le a_{i,k},b_{k,j} \le 1 \\ & y^+_{i,j},y^-_{i,j}\ge 0\end{align} These least absolute deviations models do not buy us much w.r.t. performance: we are stuck with the bilinear forms $$a_{i,k} b_{k,j}$$.