Thursday, January 22, 2015

The Hostile Brothers Problem (2)

In http://yetanothermathprogrammingconsultant.blogspot.com/2015/01/the-hostile-brothers-problem-1.html we discussed how to find local optimal solutions for a simple geometric problem. Using a local solver, we use:

image 

where model m2 uses the second (smooth) variant of:

image[5]

Multi-start Algorithm

A simple multi-start algorithm would just slap a loop around this. Indeed we can find slightly better solutions with this:

 

obj

min distance

time

local solve conopt

0.01694

0.13016

1

multi start sequential (N=50, conopt)

0.01744

0.13208

93

multi start sequential (N=50, snopt)

0.01744

0.13205

39

multi start parallel (N=500, 4 threads, snopt)

0.01748

0.13219

113

As the problems are independent of each other it is easy to solve them in parallel. Here we used a scheme where we solve four problems in parallel. If a problem is done, a new one is started, so we keep those four threads busy (no starvation). The code managing all this is written in GAMS using a few loops. Being able to write little algorithms in a modeling language can be an important feature. Using the parallel algorithm we can solve 500 problems on 4 worker threads in under 2 minutes.

Local Solver

Interestingly, the heuristic solver “Localsolver”  likes the first formulation, model m1, better. This is more compact but non-differentiable. This solver does not care about the objective function not being smooth.

LocalSolver      24.4.1 r50296 Released Dec 20, 2014 WEI x86 64bit/MS Windows

LocalSolver 5.0 (Win64, build 20141219)
Copyright (C) 2014 Innovation 24, Aix-Marseille University, CNRS.
All rights reserved (see Terms and Conditions for details).

Searching for equations that define variables... found 0
Instance:      150 Variables
                 0 Constraints
             19583 Expressions
             41929 Operands

    time |           iterations |       objective
      0s |                    0 | 0.00000000e+000
      1s |                 4274 | 5.16875195e-003
      2s |                14315 | 8.25206698e-003
      3s |                24635 | 9.05091978e-003
      4s |                35315 | 9.66455121e-003
      5s |                45914 | 9.75280995e-003
      6s |                56555 | 1.02390051e-002
      7s |                67534 | 1.05097217e-002
      8s |                78520 | 1.05893800e-002
      9s |                89524 | 1.06696606e-002
     10s |               100346 | 1.08827077e-002
    . . . 
    time |           iterations |       objective 
    495s |              5307542 | 1.58464650e-002
    496s |              5318409 | 1.58464650e-002
    497s |              5329301 | 1.58464650e-002
    498s |              5340270 | 1.58464650e-002
    499s |              5350938 | 1.58464650e-002
    500s |              5361814 | 1.58464650e-002

Stopped due to time limit.
Feasible solution found.
Objective Value: 1.58464650e-002

For this particular problem the multi-start + local NLP solver seems more efficient.

Next: more experiments.

 

Update: LocalSolver statistics

The counts for the expressions and operands are not immediately clear. To investigate we use n=3 points. The log file shows:

Instance:        6 Variables
                 0 Constraints
                35 Expressions
                61 Operands

With an option we can write out an .LSP file with localsolver function code :

.LSP file My comment

function model() {
    n0 <- 0.0;
    n1 <- 1.0;
    n2 <- float(n0, n1);
    n3 <- -1;
    n4 <- prod(n2, n3);
    n5 <- float(n0, n1);
    n6 <- sum(n4, n5);
    n7 <- prod(n6, n6);
    n8 <- float(n0, n1);
    n9 <- prod(n8, n3);
    n10 <- float(n0, n1);
    n11 <- sum(n9, n10);
    n12 <- prod(n11, n11);
    n13 <- sum(n7, n12);
    n14 <- prod(n2, n3);
    n15 <- float(n0, n1);
    n16 <- sum(n14, n15);
    n17 <- prod(n16, n16);
    n18 <- prod(n8, n3);
    n19 <- float(n0, n1);
    n20 <- sum(n18, n19);
    n21 <- prod(n20, n20);
    n22 <- sum(n17, n21);
    n23 <- prod(n15, n3);
    n24 <- sum(n23, n5);
    n25 <- prod(n24, n24);
    n26 <- prod(n19, n3);
    n27 <- sum(n26, n10);
    n28 <- prod(n27, n27);
    n29 <- sum(n25, n28);
    n30 <- min(n13, n22, n29);
    n31 <- -1.0;
    n32 <- prod(n30, n31, n31);
    n33 <- sum(n32);
    n34 <- 1;
    maximize n33;
}

 
n0,n1: Bounds on each variable
 
n2 := x1 in [0,1]
 
n4 := -x1
n5 := x2 in [0,1]  
n6 :=  (-x1) + x2 
n7 :=  (x2-x1)^2
n8 := y1 in [0,1]
n9 := -y1
n10 := y2 in [0,1]
n11 := (-y1)+y2
n12 := (y2-y1)^2
n13:= (x2-x1)^2+(y2-y1)^2
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
n30: mindist := min(n13,n22,n29)
n31-n33: Related to dummy objective function row

n33:= mindist 
n34: not used?
maximize mindist

Indeed the generated code is very low level. This certainly does not look optimal. Lets say there is room for improvement. I suspect that a smarter code will increase the number of iterations per second only (the search space is related to the number of these float variables which stays the same). Update: see comment below where some experiments show no performance degradation for this code.