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:

 

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

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.