Summary
Formula of area of a triangle
Assuming our triangle has corner points \(({\color{darkred}x}_1,{\color{darkred}y}_1)\), \(({\color{darkred}x}_2,{\color{darkred}y}_2)\) and \(({\color{darkred}x}_3,{\color{darkred}y}_3)\), the determinant formula for the area is: \[{\color{darkred}A} = 0.5 \> {\mathbf{abs}} \> {\mathbf{det}}\begin{bmatrix}{\color{darkred}x}_1 & {\color{darkred}y}_1 & 1\\{\color{darkred}x}_2 & {\color{darkred}y}_2 & 1 \\ {\color{darkred}x}_3 & {\color{darkred}y}_3 & 1\end{bmatrix}\] This can be written as: \[{\color{darkred}A}=0.5 \cdot \left|{\color{darkred}x}_1({\color{darkred}y}_2-{\color{darkred}y}_3) + {\color{darkred}x}_2({\color{darkred}y}_3-{\color{darkred}y}_1) + {\color{darkred}x}_3({\color{darkred}y}_1-{\color{darkred}y}_2) \right| \] The absolute value can be linearized easily. What remains is a non-convex quadratic expression.
Constraints forcing points to be inside the triangle
Here I used a concept called barycentric coordinates [2]. If we can form weights \({\color{darkred}\lambda}_k\) such that \[\begin{align}& {\color{darkblue}p} = \sum_k {\color{darkred}\lambda}_k \cdot {\color{darkred}t}_k \\ & {\color{darkred}\lambda}_k \ge 0 \\ & \sum_k {\color{darkred}\lambda}_k = 1 \end{align}\] then point \({\color{darkblue}p}\) is inside the triangle formed by corner points \({\color{darkred}t}_k\). This again introduces a non-convex quadratic expression as both \({\color{darkred}\lambda}_k\), and \({\color{darkred}t}_k\) are endogenous (i.e. variables). Note that \({\color{darkblue}p}\) and \({\color{darkred}t}_k\) are two-dimensional, while \({\color{darkred}\lambda}_k\) are scalars.
Complete model
Putting these two things together, we arrive at our optimization model:
| Smallest triangle model |
|---|
| \[\begin{align}\min\> & {\color{darkred}A}^+ + {\color{darkred}A}^- \\ & {\color{darkred}A}^+ - {\color{darkred}A}^- = 0.5 \cdot \left[{\color{darkred}t}_{1,x}({\color{darkred}t}_{2,y}-{\color{darkred}t}_{3,y}) + {\color{darkred}t}_{2,x}({\color{darkred}t}_{3,y}-{\color{darkred}t}_{1,y}) + {\color{darkred}t}_{3,x}({\color{darkred}t}_{1,y}-{\color{darkred}t}_{2,y}) \right] \\ & {\color{darkblue}p}_{i,c} = \sum_k {\color{darkred}\lambda}_{i,k} \cdot {\color{darkred}t}_{k,c} && \forall i,c\\ & \sum_k {\color{darkred}\lambda}_{i,k} = 1 && \forall i \\ & {\color{darkred}A}^+, {\color{darkred}A}^- \ge 0 \\ & {\color{darkred}\lambda}_{i,k} \ge 0 \\ & {\color{darkred}t}_{k,c}\> {\mathbf{free}}\\ & c \in \{x,y\}\\ & k\in \{1,2,3\}\\ & i \in \{1,\dots,n\}\end{align}\] |
We can add a constraint that orders the corner points of the triangle: \[{\color{darkred}t}_{1,x} \le {\color{darkred}t}_{2,x} \le {\color{darkred}t}_{3,x}\] This reduces the symmetry of the model a bit and can possible help global solvers.
GAMS Model
$onText
Given are n data points (2d).
Find the smallest triangle that contains all points.
$offText
option nlp=scip, reslim=1000;
*---------------------------------------------------------------------
* data: points
*---------------------------------------------------------------------
set
i 'points' /point1*point25/
c 'coordinates' /x,y/
;
parameter p(i,c) 'data points';
p(i,c) = uniform(0,100);
*---------------------------------------------------------------------
* find smallest triangle to contain all points
* original formulation
*---------------------------------------------------------------------
sets
k 'corner points of triangle' /corner1*corner3/
pm 'plusmin -- used in linearizing abs()' /'+','-'/
;
* shorthands to make our area calculation easier
singleton sets
x1(k,c) /'corner1'.'x'/
x2(k,c) /'corner2'.'x'/
x3(k,c) /'corner3'.'x'/
y1(k,c) /'corner1'.'y'/
y2(k,c) /'corner2'.'y'/
y3(k,c) /'corner3'.'y'/
;
variable
t(k,c) 'triangle'
z 'objective'
;
positive variable
area(pm) 'area (using variable splitting)'
lambda(i,k) 'barycentric coordinates'
;
equations
calcArea 'calculate area given its three corner points'
calcLambda(i,c) 'solve for barycentric coordinates'
sumLambda(i) 'lambdas need to add up to one'
obj 'objective'
order 'order corner points by their x coordinate'
;
calcArea.. area('+')-area('-') =e= 0.5*[t(x1)*(t(y2)-t(y3)) + t(x2)*(t(y3)-t(y1)) + t(x3)*(t(y1)-t(y2))];
calcLambda(i,c).. p(i,c) =e= sum(k, lambda(i,k)*t(k,c));
sumLambda(i).. sum(k, lambda(i,k)) =e= 1;
obj.. z =e= sum(pm,area(pm));
order(k-1).. t(k,'x') =g= t(k-1,'x');
* some reasonable bounds
t.lo(k,c) = -1000;
t.up(k,c) = +1000;
model m1 /all/;
solve m1 minimizing z using nlp;
*---------------------------------------------------------------------
* reporting
*---------------------------------------------------------------------
display p,t.l,area.l,lambda.l;
Looking at the picture, it may make sense to preprocess the data and only keep the points that are part of the convex hull. Furthermore, and may be less interesting, the area of the triangle can not be smaller than the area of the convex hull.
Smallest case
For the \(n=3\) case, the problem is exceedingly simple. The optimal corner points are obviously the same as the three data points. The resulting non-convex NLP model is very small:
MODEL STATISTICS BLOCKS OF EQUATIONS 5 SINGLE EQUATIONS 13 BLOCKS OF VARIABLES 4 SINGLE VARIABLES 18 NON ZERO ELEMENTS 60 NON LINEAR N-Z 42 CODE LENGTH 84 CONSTANT POOL 16
The data and and solution are listed here:
---- 79 PARAMETER p data points x y point1 17.175 84.327 point2 55.038 30.114 point3 29.221 22.405 ---- 79 VARIABLE t.L triangle x y corner1 17.175 84.327 corner2 29.221 22.405 corner3 55.038 30.114 ---- 79 VARIABLE area.L area (using variable splitting) + 845.721 ---- 79 VARIABLE lambda.L barycentric coordinates corner1 corner2 corner3 point1 1.000 point2 1.000 point3 1.000
We can see that the solution indeed is optimal: the corner points of the triangle are indentical to our data points.
However none of the solvers I tried, could prove optimality for this small case in 1,000 seconds. They all had troubles in moving the lower bound away from 0. Some fragments of the progress logs are listed below.
SCIP log
time | node | left |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr| dualbound | primalbound | gap | compl. 998s| 1087k| 19335 | 28046k| 25.8 | 433M | 56 | 48 | 12 | 234 |8946k| 2 |2323 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 998s| 1087k| 19323 | 28048k| 25.8 | 433M | 56 | 48 | 12 | 0 |8947k| 0 |2323 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 998s| 1088k| 19325 | 28051k| 25.8 | 433M | 56 | 48 | 12 | 262 |8948k| 5 |2323 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 998s| 1088k| 19337 | 28053k| 25.8 | 433M | 56 | 48 | 12 | 252 |8949k| 3 |2323 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 998s| 1088k| 19341 | 28056k| 25.8 | 433M | 56 | 48 | 12 | 272 |8950k| 3 |2323 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 998s| 1088k| 19337 | 28059k| 25.8 | 433M | 56 | 48 | 12 | 267 |8951k| 2 |2323 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 998s| 1088k| 19335 | 28061k| 25.8 | 433M | 56 | 48 | 12 | 257 |8951k| 2 |2323 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 999s| 1088k| 19343 | 28064k| 25.8 | 433M | 56 | 48 | 12 | 221 |8952k| 3 |2323 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 999s| 1088k| 19343 | 28067k| 25.8 | 433M | 56 | 48 | 12 | 255 |8953k| 5 |2324 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 999s| 1088k| 19353 | 28069k| 25.8 | 433M | 56 | 48 | 12 | 240 |8954k| 3 |2324 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 999s| 1088k| 19355 | 28072k| 25.8 | 433M | 56 | 48 | 12 | 239 |8955k| 2 |2324 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 999s| 1088k| 19355 | 28075k| 25.8 | 433M | 56 | 48 | 12 | 219 |8956k| 1 |2324 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 999s| 1089k| 19351 | 28076k| 25.8 | 434M | 56 | 48 | 12 | 202 |8957k| 4 |2324 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 999s| 1089k| 19373 | 28078k| 25.8 | 434M | 56 | 48 | 12 | 237 |8957k| 2 |2324 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 999s| 1089k| 19367 | 28080k| 25.8 | 434M | 56 | 48 | 12 | 259 |8958k| 4 |2324 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% time | node | left |LP iter|LP it/n|mem/heur|mdpt |vars |cons |rows |cuts |sepa|confs|strbr| dualbound | primalbound | gap | compl. 999s| 1089k| 19355 | 28081k| 25.8 | 434M | 56 | 48 | 12 | 231 |8959k| 3 |2325 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 999s| 1089k| 19357 | 28083k| 25.8 | 434M | 56 | 48 | 12 | 257 |8960k| 2 |2325 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% 999s| 1089k| 19365 | 28085k| 25.8 | 434M | 56 | 48 | 12 | 207 |8961k| 2 |2325 | 0 | 0.000000e+00 | 8.457215e+02 | Inf | 42.07% SCIP Status : solving was interrupted [time limit reached] Solving Time (sec) : 1000.00 Solving Nodes : 1089501 Primal Bound : +8.45721492842725e+02 (3 solutions) Dual Bound : +0.00000000000000e+00 Gap : infinite
The primal bound represents the best solution found so far, and the dual bound means the best possible solution. The last column has an indication of how far we are in the completion of the run. That is an estimate.
Baron log
Doing local search Preprocessing found feasible solution with value 845.721 Solving bounding LP Starting multi-start local search Done with local search =========================================================================== Iteration Time (s) Mem Lower bound Upper bound Progress 1 0.17 12MB 0.00000 845.721 NA 46639 29.98 12MB 0.00000 845.721 7.66% 83130+ 59.91 12MB 0.00000 845.721 9.98% 108649 89.88 12MB 0.00000 845.721 11.00% 135767 119.81 12MB 0.00000 845.721 12.04% 160770 149.72 12MB 0.00000 845.721 13.22% 183770 179.69 12MB 0.00000 845.721 14.57% 203933 209.64 12MB 0.00000 845.721 15.47% 222412 239.53 12MB 0.00000 845.721 15.98% 241749 269.39 12MB 0.00000 845.721 16.90% 260570 299.28 12MB 0.00000 845.721 17.73% 277359 329.22 12MB 0.00000 845.721 18.28% 293639 359.16 12MB 0.00000 845.721 18.63% 310561 389.03 12MB 0.00000 845.721 19.02% 326955+ 418.97 12MB 0.00000 845.721 19.68% 341310 448.86 12MB 0.00000 845.721 19.92% 356540 478.78 12MB 0.00000 845.721 20.08% 371601 508.66 12MB 0.00000 845.721 20.38% 386201 538.64 12MB 0.00000 845.721 20.68% 400581 568.58 12MB 0.00000 845.721 21.40% 414407 598.48 12MB 0.00000 845.721 21.74% 427627 628.44 12MB 0.00000 845.721 21.96% 440293 658.41 12MB 0.00000 845.721 22.18% 451183 688.34 12MB 0.00000 845.721 22.46% 463509 718.30 12MB 0.00000 845.721 22.86% 475163 748.22 12MB 0.00000 845.721 23.30% 484121 778.16 12MB 0.00000 845.721 23.53% 495774 808.09 12MB 0.00000 845.721 23.73% 507127 837.98 12MB 0.00000 845.721 24.06% 517685 867.92 12MB 0.00000 845.721 24.34% 528924 897.78 12MB 0.00000 845.721 24.48% 540003 927.70 12MB 0.00000 845.721 24.92% 551003 957.62 12MB 0.00000 845.721 25.21% 561686 987.58 12MB 0.00000 845.721 25.46% 566020 1000.00 12MB 0.00000 845.721 25.60% *** Max. allowable time exceeded *** Wall clock time: 1002.74 Total CPU time used: 1000.00 Total no. of BaR iterations: 566020 Best solution found at node: -1 Max. no. of nodes in memory: 590 All done
Gurobi log
Continuous model is non-convex -- solving as a MIP
Presolve time: 0.00s
Presolved: 108 rows, 42 columns, 291 nonzeros
Presolved model has 24 bilinear constraint(s)
Variable types: 42 continuous, 0 integer (0 binary)
Found heuristic solution: objective 845.7214933
Root relaxation: objective 0.000000e+00, 34 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
0 0 0.00000 0 24 845.72149 0.00000 100% - 0s
0 0 0.00000 0 24 845.72149 0.00000 100% - 0s
H 0 0 845.7214932 0.00000 100% - 0s
0 2 0.00000 0 24 845.72149 0.00000 100% - 0s
87509 2820 infeasible 56 845.72149 0.00000 100% 7.2 5s
181089 4571 0.00000 55 7 845.72149 0.00000 100% 7.5 10s
279430 5733 0.00000 34 10 845.72149 0.00000 100% 7.2 15s
365052 7563 infeasible 57 845.72149 0.00000 100% 7.1 20s
448542 7794 0.00000 58 9 845.72149 0.00000 100% 7.5 25s
542861 9533 infeasible 75 845.72149 0.00000 100% 7.5 30s
634169 11687 0.00000 76 11 845.72149 0.00000 100% 7.6 35s
723679 13899 0.00000 39 8 845.72149 0.00000 100% 7.6 40s
812858 16683 infeasible 42 845.72149 0.00000 100% 7.6 45s
900001 19199 infeasible 65 845.72149 0.00000 100% 7.6 50s
981975 21796 0.00000 48 19 845.72149 0.00000 100% 7.7 55s
1060351 22782 0.00000 66 6 845.72149 0.00000 100% 7.8 60s
1146405 23023 0.00000 49 14 845.72149 0.00000 100% 7.9 65s
1231170 23371 0.00000 40 13 845.72149 0.00000 100% 8.0 70s
1311008 23228 0.00000 53 15 845.72149 0.00000 100% 8.2 75s
1392856 23692 0.00000 65 10 845.72149 0.00000 100% 8.3 80s
1484084 24094 0.00000 82 13 845.72149 0.00000 100% 8.3 85s
1586221 25844 0.00000 56 7 845.72149 0.00000 100% 8.2 90s
1683920 26664 0.00000 51 14 845.72149 0.00000 100% 8.2 95s
1781376 27253 infeasible 60 845.72149 0.00000 100% 8.2 100s
1879069 27233 0.00000 57 22 845.72149 0.00000 100% 8.2 105s
1978438 27317 0.00000 62 5 845.72149 0.00000 100% 8.2 110s
2079788 27113 infeasible 46 845.72149 0.00000 100% 8.1 115s
2176369 27299 0.00000 53 7 845.72149 0.00000 100% 8.2 120s
2273400 27973 infeasible 118 845.72149 0.00000 100% 8.2 125s
2382555 28076 infeasible 52 845.72149 0.00000 100% 8.1 130s
2490061 28289 infeasible 51 845.72149 0.00000 100% 8.0 135s
2596000 29069 0.00000 53 16 845.72149 0.00000 100% 7.9 140s
2704736 28963 infeasible 47 845.72149 0.00000 100% 7.8 145s
2814008 28685 0.00000 70 13 845.72149 0.00000 100% 7.7 150s
2923779 29203 infeasible 80 845.72149 0.00000 100% 7.7 155s
3023688 29277 infeasible 67 845.72149 0.00000 100% 7.6 160s
3106302 29297 infeasible 45 845.72149 0.00000 100% 7.7 165s
3185810 30783 infeasible 42 845.72149 0.00000 100% 7.7 170s
3262367 30697 0.00000 43 22 845.72149 0.00000 100% 7.7 175s
3336095 30796 0.00000 72 13 845.72149 0.00000 100% 7.7 180s
3407884 30514 0.00000 50 10 845.72149 0.00000 100% 7.7 185s
3489921 30840 0.00000 62 14 845.72149 0.00000 100% 7.7 190s
3579425 30649 0.00000 52 12 845.72149 0.00000 100% 7.7 195s
3662157 31234 infeasible 74 845.72149 0.00000 100% 7.7 200s
3772556 31551 infeasible 73 845.72149 0.00000 100% 7.6 205s
3879636 31731 infeasible 61 845.72149 0.00000 100% 7.6 210s
3983022 32340 infeasible 76 845.72149 0.00000 100% 7.6 215s
4085295 32582 0.00000 60 6 845.72149 0.00000 100% 7.5 220s
4195244 32612 0.00000 85 6 845.72149 0.00000 100% 7.5 225s
4301486 32927 0.00000 68 8 845.72149 0.00000 100% 7.5 230s
4406655 33963 infeasible 68 845.72149 0.00000 100% 7.4 235s
4515096 36747 infeasible 49 845.72149 0.00000 100% 7.4 240s
4616212 39197 infeasible 89 845.72149 0.00000 100% 7.4 245s
4713644 39205 infeasible 65 845.72149 0.00000 100% 7.4 250s
4811609 39036 infeasible 93 845.72149 0.00000 100% 7.4 255s
4911613 39133 0.00000 46 12 845.72149 0.00000 100% 7.4 260s
5004637 39457 infeasible 78 845.72149 0.00000 100% 7.4 265s
5100653 39188 0.00000 57 8 845.72149 0.00000 100% 7.4 270s
5196726 39046 infeasible 59 845.72149 0.00000 100% 7.4 275s
5294558 39395 infeasible 50 845.72149 0.00000 100% 7.5 280s
5392429 39475 0.00000 44 12 845.72149 0.00000 100% 7.5 285s
5490479 39240 0.00000 63 6 845.72149 0.00000 100% 7.5 290s
5584131 39506 infeasible 47 845.72149 0.00000 100% 7.5 295s
5675931 39754 0.00000 44 6 845.72149 0.00000 100% 7.5 300s
5766067 40631 infeasible 39 845.72149 0.00000 100% 7.5 305s
5858195 40882 infeasible 62 845.72149 0.00000 100% 7.5 310s
5957180 41961 infeasible 67 845.72149 0.00000 100% 7.5 315s
6044532 42208 0.00000 38 13 845.72149 0.00000 100% 7.6 320s
6139051 42361 infeasible 71 845.72149 0.00000 100% 7.6 325s
6240644 42187 0.00000 49 13 845.72149 0.00000 100% 7.6 330s
6332834 42892 infeasible 81 845.72149 0.00000 100% 7.6 335s
6431141 43094 infeasible 94 845.72149 0.00000 100% 7.6 340s
6529490 42822 0.00000 68 10 845.72149 0.00000 100% 7.6 345s
6633781 42809 0.00000 64 5 845.72149 0.00000 100% 7.6 350s
6731081 42397 infeasible 46 845.72149 0.00000 100% 7.6 355s
6820057 42950 0.00000 71 12 845.72149 0.00000 100% 7.6 360s
6923260 42686 0.00000 65 10 845.72149 0.00000 100% 7.6 365s
7020177 42544 infeasible 67 845.72149 0.00000 100% 7.6 370s
7109967 42336 infeasible 79 845.72149 0.00000 100% 7.6 375s
7204990 43282 0.00000 56 11 845.72149 0.00000 100% 7.6 380s
7310428 43391 0.00000 53 9 845.72149 0.00000 100% 7.6 385s
7411947 44120 cutoff 67 845.72149 0.00000 100% 7.6 390s
7514017 44406 infeasible 68 845.72149 0.00000 100% 7.6 395s
7624230 44197 0.00000 52 15 845.72149 0.00000 100% 7.6 400s
7724371 44350 0.00000 63 16 845.72149 0.00000 100% 7.6 405s
7824546 43992 0.00000 82 9 845.72149 0.00000 100% 7.5 410s
7931695 44272 0.00000 70 5 845.72149 0.00000 100% 7.5 415s
8032771 45062 0.00000 53 8 845.72149 0.00000 100% 7.5 420s
8133197 45314 infeasible 89 845.72149 0.00000 100% 7.5 425s
8239816 44949 infeasible 73 845.72149 0.00000 100% 7.5 430s
8341105 45211 infeasible 80 845.72149 0.00000 100% 7.5 435s
8444474 46196 0.00000 70 6 845.72149 0.00000 100% 7.5 440s
8540717 46132 0.00000 73 5 845.72149 0.00000 100% 7.5 445s
8633853 47775 infeasible 51 845.72149 0.00000 100% 7.5 450s
8718813 47992 infeasible 55 845.72149 0.00000 100% 7.5 455s
8801173 48442 infeasible 76 845.72149 0.00000 100% 7.5 460s
8888173 48821 infeasible 82 845.72149 0.00000 100% 7.5 465s
8969838 48709 0.00000 76 7 845.72149 0.00000 100% 7.5 470s
9054901 48735 0.00000 86 12 845.72149 0.00000 100% 7.5 475s
9139983 49013 infeasible 58 845.72149 0.00000 100% 7.5 480s
9222824 49141 0.00000 58 12 845.72149 0.00000 100% 7.5 485s
9305478 49141 0.00000 53 5 845.72149 0.00000 100% 7.5 490s
9386036 49122 0.00000 55 14 845.72149 0.00000 100% 7.5 495s
9470421 48960 infeasible 54 845.72149 0.00000 100% 7.6 500s
9554389 49140 infeasible 43 845.72149 0.00000 100% 7.6 505s
9652375 49287 0.00000 47 5 845.72149 0.00000 100% 7.5 510s
9761322 49050 0.00000 68 11 845.72149 0.00000 100% 7.5 515s
9867290 48930 0.00000 53 18 845.72149 0.00000 100% 7.5 520s
9965913 49688 0.00000 28 10 845.72149 0.00000 100% 7.5 525s
10071315 51026 infeasible 53 845.72149 0.00000 100% 7.5 530s
10171544 50502 infeasible 85 845.72149 0.00000 100% 7.5 535s
10274791 51927 0.00000 120 19 845.72149 0.00000 100% 7.5 540s
10384316 52084 cutoff 70 845.72149 0.00000 100% 7.5 545s
10494572 51820 0.00000 99 17 845.72149 0.00000 100% 7.5 550s
10602563 51755 infeasible 117 845.72149 0.00000 100% 7.5 555s
10710904 51585 infeasible 73 845.72149 0.00000 100% 7.5 560s
10811637 51431 0.00000 49 16 845.72149 0.00000 100% 7.5 565s
10915146 51471 infeasible 59 845.72149 0.00000 100% 7.4 570s
11014799 51749 infeasible 87 845.72149 0.00000 100% 7.5 575s
11116945 51772 0.00000 44 9 845.72149 0.00000 100% 7.5 580s
11221881 51849 0.00000 64 8 845.72149 0.00000 100% 7.5 585s
11321489 51975 infeasible 81 845.72149 0.00000 100% 7.5 590s
11424641 51778 0.00000 73 12 845.72149 0.00000 100% 7.4 595s
11525725 51688 0.00000 66 5 845.72149 0.00000 100% 7.5 600s
11626520 51958 0.00000 79 8 845.72149 0.00000 100% 7.5 605s
11737524 51917 0.00000 84 12 845.72149 0.00000 100% 7.4 610s
11842191 53086 infeasible 59 845.72149 0.00000 100% 7.4 615s
11944591 53012 infeasible 46 845.72149 0.00000 100% 7.4 620s
12050709 52948 0.00000 78 9 845.72149 0.00000 100% 7.4 625s
12150122 52943 0.00000 67 20 845.72149 0.00000 100% 7.4 630s
12256127 52960 infeasible 80 845.72149 0.00000 100% 7.4 635s
12358157 52755 0.00000 53 13 845.72149 0.00000 100% 7.4 640s
12460085 52935 0.00000 74 13 845.72149 0.00000 100% 7.4 645s
12560196 52859 infeasible 66 845.72149 0.00000 100% 7.4 650s
12669140 52870 infeasible 61 845.72149 0.00000 100% 7.4 655s
12777770 53088 cutoff 58 845.72149 0.00000 100% 7.4 660s
12889503 53282 infeasible 67 845.72149 0.00000 100% 7.4 665s
12998803 53170 0.00000 66 9 845.72149 0.00000 100% 7.4 670s
13105332 53408 infeasible 44 845.72149 0.00000 100% 7.4 675s
13203483 53394 infeasible 68 845.72149 0.00000 100% 7.4 680s
13305462 53413 infeasible 45 845.72149 0.00000 100% 7.4 685s
13409735 54002 0.00000 54 12 845.72149 0.00000 100% 7.4 690s
13517804 53935 0.00000 85 7 845.72149 0.00000 100% 7.4 695s
13619933 53675 0.00000 47 13 845.72149 0.00000 100% 7.4 700s
13728120 54287 infeasible 59 845.72149 0.00000 100% 7.3 705s
13830286 54754 infeasible 63 845.72149 0.00000 100% 7.3 710s
13938426 54802 infeasible 65 845.72149 0.00000 100% 7.3 715s
14037753 55496 0.00000 72 6 845.72149 0.00000 100% 7.4 720s
14136366 55286 0.00000 62 9 845.72149 0.00000 100% 7.4 725s
14236512 55419 0.00000 57 8 845.72149 0.00000 100% 7.4 730s
14336714 55893 infeasible 88 845.72149 0.00000 100% 7.4 735s
14438374 55833 0.00000 77 5 845.72149 0.00000 100% 7.4 740s
14541236 55950 infeasible 61 845.72149 0.00000 100% 7.4 745s
14643803 56430 infeasible 89 845.72149 0.00000 100% 7.4 750s
14745667 56477 infeasible 67 845.72149 0.00000 100% 7.4 755s
14849383 56235 0.00000 47 13 845.72149 0.00000 100% 7.4 760s
14953386 56491 0.00000 96 6 845.72149 0.00000 100% 7.4 765s
15058140 56442 0.00000 100 8 845.72149 0.00000 100% 7.4 770s
15165635 56486 infeasible 63 845.72149 0.00000 100% 7.4 775s
15267411 56360 infeasible 43 845.72149 0.00000 100% 7.4 780s
15377210 56653 infeasible 47 845.72149 0.00000 100% 7.3 785s
15482977 57402 0.00000 115 13 845.72149 0.00000 100% 7.3 790s
15587571 56840 infeasible 85 845.72149 0.00000 100% 7.3 795s
15685438 57174 0.00000 85 7 845.72149 0.00000 100% 7.3 800s
15780069 57188 infeasible 84 845.72149 0.00000 100% 7.3 805s
15879876 57387 0.00000 78 12 845.72149 0.00000 100% 7.3 810s
15983035 57256 infeasible 112 845.72149 0.00000 100% 7.3 815s
16088941 57276 infeasible 95 845.72149 0.00000 100% 7.3 820s
16194954 56854 infeasible 83 845.72149 0.00000 100% 7.3 825s
16314096 56933 0.00000 63 5 845.72149 0.00000 100% 7.3 830s
16420921 56936 infeasible 66 845.72149 0.00000 100% 7.3 835s
16526928 57212 infeasible 80 845.72149 0.00000 100% 7.3 840s
16632943 57436 cutoff 71 845.72149 0.00000 100% 7.3 845s
16740271 57705 infeasible 60 845.72149 0.00000 100% 7.3 850s
16848014 58106 infeasible 70 845.72149 0.00000 100% 7.3 855s
16956660 58089 infeasible 108 845.72149 0.00000 100% 7.3 860s
17061684 58186 0.00000 98 13 845.72149 0.00000 100% 7.3 865s
17165491 57648 infeasible 62 845.72149 0.00000 100% 7.3 870s
17263196 57489 0.00000 69 10 845.72149 0.00000 100% 7.3 875s
17363860 57819 0.00000 72 3 845.72149 0.00000 100% 7.3 880s
17466234 57672 0.00000 74 7 845.72149 0.00000 100% 7.3 885s
17568389 57655 infeasible 76 845.72149 0.00000 100% 7.3 890s
17668292 57462 infeasible 54 845.72149 0.00000 100% 7.3 895s
17764449 57740 infeasible 61 845.72149 0.00000 100% 7.3 900s
17857454 57815 0.00000 76 16 845.72149 0.00000 100% 7.3 905s
17954845 57880 0.00000 40 8 845.72149 0.00000 100% 7.3 910s
18055547 57558 infeasible 66 845.72149 0.00000 100% 7.3 915s
18168860 57680 infeasible 63 845.72149 0.00000 100% 7.3 920s
18276707 57757 infeasible 65 845.72149 0.00000 100% 7.3 925s
18377132 57949 0.00000 57 16 845.72149 0.00000 100% 7.3 930s
18477251 58047 0.00000 61 8 845.72149 0.00000 100% 7.3 935s
18574743 57844 infeasible 59 845.72149 0.00000 100% 7.3 940s
18680752 57956 infeasible 61 845.72149 0.00000 100% 7.3 945s
18783011 57797 infeasible 60 845.72149 0.00000 100% 7.3 950s
18903457 57699 0.00000 47 17 845.72149 0.00000 100% 7.3 955s
19009899 57713 infeasible 74 845.72149 0.00000 100% 7.3 960s
19106418 57638 0.00000 45 20 845.72149 0.00000 100% 7.3 965s
19202306 57961 infeasible 70 845.72149 0.00000 100% 7.3 970s
19301058 58150 0.00000 49 17 845.72149 0.00000 100% 7.3 975s
19398515 58371 infeasible 66 845.72149 0.00000 100% 7.3 980s
19494042 58225 infeasible 73 845.72149 0.00000 100% 7.3 985s
19590342 58223 0.00000 56 7 845.72149 0.00000 100% 7.3 990s
19686931 58659 infeasible 67 845.72149 0.00000 100% 7.3 995s
19779988 58357 infeasible 57 845.72149 0.00000 100% 7.3 1000s
Explored 19780322 nodes (145136295 simplex iterations) in 1000.00 seconds (806.64 work units)
Thread count was 4 (of 72 available processors)
Solution count 1: 845.721
Time limit reached
Best objective 8.457214931859e+02, best bound 0.000000000000e+00, gap 100.0000%
Antigone log
Before Pre-processing:
17 Variables
17 Continuous
13 Equations
After Pre-processing:
16 Variables
16 Continuous
31 Equations
5 Linear
26 Nonconvex nonlinear
60 Nonlinear Terms
60 Bilinear/Quadratic
64 Possible Reformulation Linearization Technique (RLT) equations
18 RLT Equations Added Outright to Formulation
Constituent Libraries:
CPLEX Solving relaxations
CONOPT Finding feasible points
LAPACK Addressing linear systems
Boost Bounding Intervals
-------------------------------------------------------------------------------
Time (s) Nodes explored Nodes remaining Best possible Best found Relative Gap
-------------------------------------------------------------------------------
0 1 1 +0.000e+00 +1.240e+03 --
5 433 379 +0.000e+00 +8.457e+02 --
10 843 702 +0.000e+00 +8.457e+02 --
15 1347 1108 +0.000e+00 +8.457e+02 --
20 1825 1486 +0.000e+00 +8.457e+02 --
25 2387 1894 +0.000e+00 +8.457e+02 --
30 2990 2336 +0.000e+00 +8.457e+02 --
35 3697 2832 +0.000e+00 +8.457e+02 --
40 4438 3361 +0.000e+00 +8.457e+02 --
45 5086 3773 +0.000e+00 +8.457e+02 --
50 5830 4297 +0.000e+00 +8.457e+02 --
55 6688 4786 +0.000e+00 +8.457e+02 --
60 7211 5208 +0.000e+00 +8.457e+02 --
65 7883 5540 +0.000e+00 +8.457e+02 --
70 8602 6102 +0.000e+00 +8.457e+02 --
75 9508 6670 +0.000e+00 +8.457e+02 --
80 10473 7206 +0.000e+00 +8.457e+02 --
85 11018 7613 +0.000e+00 +8.457e+02 --
90 11636 8067 +0.000e+00 +8.457e+02 --
95 12379 8369 +0.000e+00 +8.457e+02 --
100 13051 8867 +0.000e+00 +8.457e+02 --
105 13795 9370 +0.000e+00 +8.457e+02 --
110 14819 9922 +0.000e+00 +8.457e+02 --
115 15700 10385 +0.000e+00 +8.457e+02 --
120 16594 10889 +0.000e+00 +8.457e+02 --
125 17116 11255 +0.000e+00 +8.457e+02 --
130 17758 11702 +0.000e+00 +8.457e+02 --
135 18329 12114 +0.000e+00 +8.457e+02 --
140 19193 12401 +0.000e+00 +8.457e+02 --
145 19894 12929 +0.000e+00 +8.457e+02 --
150 20660 13440 +0.000e+00 +8.457e+02 --
155 21367 13927 +0.000e+00 +8.457e+02 --
160 22378 14606 +0.000e+00 +8.457e+02 --
165 23197 14983 +0.000e+00 +8.457e+02 --
170 24116 15410 +0.000e+00 +8.457e+02 --
175 25021 15898 +0.000e+00 +8.457e+02 --
180 25633 16196 +0.000e+00 +8.457e+02 --
185 26072 16523 +0.000e+00 +8.457e+02 --
190 26607 16849 +0.000e+00 +8.457e+02 --
195 27320 17301 +0.000e+00 +8.457e+02 --
200 27946 17721 +0.000e+00 +8.457e+02 --
205 28466 18098 +0.000e+00 +8.457e+02 --
210 29406 18276 +0.000e+00 +8.457e+02 --
215 30095 18787 +0.000e+00 +8.457e+02 --
220 30626 19108 +0.000e+00 +8.457e+02 --
225 31459 19587 +0.000e+00 +8.457e+02 --
230 32137 20010 +0.000e+00 +8.457e+02 --
235 32828 20394 +0.000e+00 +8.457e+02 --
240 33823 21065 +0.000e+00 +8.457e+02 --
245 34913 21483 +0.000e+00 +8.457e+02 --
250 35789 21843 +0.000e+00 +8.457e+02 --
255 36828 22172 +0.000e+00 +8.457e+02 --
260 37917 22540 +0.000e+00 +8.457e+02 --
265 38605 22856 +0.000e+00 +8.457e+02 --
270 39081 23131 +0.000e+00 +8.457e+02 --
275 39616 23418 +0.000e+00 +8.457e+02 --
280 40207 23790 +0.000e+00 +8.457e+02 --
285 40831 24234 +0.000e+00 +8.457e+02 --
290 41596 24664 +0.000e+00 +8.457e+02 --
295 42243 25099 +0.000e+00 +8.457e+02 --
300 42746 25481 +0.000e+00 +8.457e+02 --
305 43307 25886 +0.000e+00 +8.457e+02 --
310 44249 26106 +0.000e+00 +8.457e+02 --
315 45119 26582 +0.000e+00 +8.457e+02 --
320 45788 27037 +0.000e+00 +8.457e+02 --
325 46548 27358 +0.000e+00 +8.457e+02 --
330 47445 27875 +0.000e+00 +8.457e+02 --
335 48144 28225 +0.000e+00 +8.457e+02 --
340 49027 28754 +0.000e+00 +8.457e+02 --
345 49770 29128 +0.000e+00 +8.457e+02 --
350 50888 29861 +0.000e+00 +8.457e+02 --
355 52067 30355 +0.000e+00 +8.457e+02 --
360 53094 30614 +0.000e+00 +8.457e+02 --
365 54119 30950 +0.000e+00 +8.457e+02 --
370 55342 31356 +0.000e+00 +8.457e+02 --
375 56501 31597 +0.000e+00 +8.457e+02 --
380 57249 31963 +0.000e+00 +8.457e+02 --
385 57778 32215 +0.000e+00 +8.457e+02 --
390 58335 32567 +0.000e+00 +8.457e+02 --
395 59041 32807 +0.000e+00 +8.457e+02 --
400 59728 33196 +0.000e+00 +8.457e+02 --
405 60379 33594 +0.000e+00 +8.457e+02 --
410 61072 33963 +0.000e+00 +8.457e+02 --
415 61930 34437 +0.000e+00 +8.457e+02 --
420 62622 34877 +0.000e+00 +8.457e+02 --
425 63207 35278 +0.000e+00 +8.457e+02 --
430 63779 35683 +0.000e+00 +8.457e+02 --
435 64456 36194 +0.000e+00 +8.457e+02 --
440 65477 36428 +0.000e+00 +8.457e+02 --
445 66509 36986 +0.000e+00 +8.457e+02 --
450 67297 37415 +0.000e+00 +8.457e+02 --
455 68090 37786 +0.000e+00 +8.457e+02 --
460 68948 38110 +0.000e+00 +8.457e+02 --
465 69802 38599 +0.000e+00 +8.457e+02 --
470 70692 38968 +0.000e+00 +8.457e+02 --
475 71378 39369 +0.000e+00 +8.457e+02 --
480 72273 39905 +0.000e+00 +8.457e+02 --
485 73018 40298 +0.000e+00 +8.457e+02 --
490 73961 40756 +0.000e+00 +8.457e+02 --
495 75039 41424 +0.000e+00 +8.457e+02 --
500 76196 42090 +0.000e+00 +8.457e+02 --
505 77301 42170 +0.000e+00 +8.457e+02 --
510 78362 42470 +0.000e+00 +8.457e+02 --
515 79504 42730 +0.000e+00 +8.457e+02 --
520 80809 43201 +0.000e+00 +8.457e+02 --
525 81915 43466 +0.000e+00 +8.457e+02 --
530 82721 43861 +0.000e+00 +8.457e+02 --
535 83358 44003 +0.000e+00 +8.457e+02 --
540 83897 44219 +0.000e+00 +8.457e+02 --
545 84406 44557 +0.000e+00 +8.457e+02 --
550 85043 44791 +0.000e+00 +8.457e+02 --
555 85811 45070 +0.000e+00 +8.457e+02 --
560 86460 45465 +0.000e+00 +8.457e+02 --
565 87057 45832 +0.000e+00 +8.457e+02 --
570 87828 46153 +0.000e+00 +8.457e+02 --
575 88499 46548 +0.000e+00 +8.457e+02 --
580 89462 47010 +0.000e+00 +8.457e+02 --
585 90263 47576 +0.000e+00 +8.457e+02 --
590 91003 47997 +0.000e+00 +8.457e+02 --
595 91682 48408 +0.000e+00 +8.457e+02 --
600 92325 48863 +0.000e+00 +8.457e+02 --
605 92957 49268 +0.000e+00 +8.457e+02 --
610 93766 49738 +0.000e+00 +8.457e+02 --
615 94749 49965 +0.000e+00 +8.457e+02 --
620 95742 50399 +0.000e+00 +8.457e+02 --
625 96694 50847 +0.000e+00 +8.457e+02 --
630 97431 51270 +0.000e+00 +8.457e+02 --
635 98201 51608 +0.000e+00 +8.457e+02 --
640 99008 51882 +0.000e+00 +8.457e+02 --
645 99887 52171 +0.000e+00 +8.457e+02 --
650 100738 52691 +0.000e+00 +8.457e+02 --
655 101821 53028 +0.000e+00 +8.457e+02 --
660 102468 53323 +0.000e+00 +8.457e+02 --
665 103272 53792 +0.000e+00 +8.457e+02 --
670 104029 54237 +0.000e+00 +8.457e+02 --
675 104831 54726 +0.000e+00 +8.457e+02 --
680 105636 54997 +0.000e+00 +8.457e+02 --
685 106768 55411 +0.000e+00 +8.457e+02 --
690 107980 56031 +0.000e+00 +8.457e+02 --
695 109089 56669 +0.000e+00 +8.457e+02 --
700 110305 57001 +0.000e+00 +8.457e+02 --
705 111414 57065 +0.000e+00 +8.457e+02 --
710 112274 57323 +0.000e+00 +8.457e+02 --
715 112983 57580 +0.000e+00 +8.457e+02 --
720 114078 57791 +0.000e+00 +8.457e+02 --
725 115340 58157 +0.000e+00 +8.457e+02 --
730 116389 58318 +0.000e+00 +8.457e+02 --
735 117148 58603 +0.000e+00 +8.457e+02 --
740 117879 58836 +0.000e+00 +8.457e+02 --
745 118530 58918 +0.000e+00 +8.457e+02 --
750 119007 59167 +0.000e+00 +8.457e+02 --
755 119550 59442 +0.000e+00 +8.457e+02 --
760 120089 59725 +0.000e+00 +8.457e+02 --
765 120721 59866 +0.000e+00 +8.457e+02 --
770 121444 60090 +0.000e+00 +8.457e+02 --
775 122090 60387 +0.000e+00 +8.457e+02 --
780 122749 60693 +0.000e+00 +8.457e+02 --
785 123223 61036 +0.000e+00 +8.457e+02 --
790 123870 61339 +0.000e+00 +8.457e+02 --
795 124649 61611 +0.000e+00 +8.457e+02 --
800 125328 61941 +0.000e+00 +8.457e+02 --
805 126257 62287 +0.000e+00 +8.457e+02 --
810 127197 62692 +0.000e+00 +8.457e+02 --
815 127980 63130 +0.000e+00 +8.457e+02 --
820 128837 63618 +0.000e+00 +8.457e+02 --
825 129481 63912 +0.000e+00 +8.457e+02 --
830 130108 64249 +0.000e+00 +8.457e+02 --
835 130672 64643 +0.000e+00 +8.457e+02 --
840 131323 65072 +0.000e+00 +8.457e+02 --
845 131905 65372 +0.000e+00 +8.457e+02 --
850 132518 65732 +0.000e+00 +8.457e+02 --
855 133400 66162 +0.000e+00 +8.457e+02 --
860 134488 66280 +0.000e+00 +8.457e+02 --
865 135411 66728 +0.000e+00 +8.457e+02 --
870 136514 67209 +0.000e+00 +8.457e+02 --
875 137480 67582 +0.000e+00 +8.457e+02 --
880 138321 67995 +0.000e+00 +8.457e+02 --
885 139135 68352 +0.000e+00 +8.457e+02 --
890 140030 68540 +0.000e+00 +8.457e+02 --
895 141029 68847 +0.000e+00 +8.457e+02 --
900 142007 69231 +0.000e+00 +8.457e+02 --
905 142902 69762 +0.000e+00 +8.457e+02 --
910 144160 69988 +0.000e+00 +8.457e+02 --
915 144823 70250 +0.000e+00 +8.457e+02 --
920 145536 70583 +0.000e+00 +8.457e+02 --
925 146533 71159 +0.000e+00 +8.457e+02 --
930 147549 71433 +0.000e+00 +8.457e+02 --
935 148371 71859 +0.000e+00 +8.457e+02 --
940 149116 72119 +0.000e+00 +8.457e+02 --
945 150064 72519 +0.000e+00 +8.457e+02 --
950 151389 72954 +0.000e+00 +8.457e+02 --
955 152552 73388 +0.000e+00 +8.457e+02 --
960 153719 73969 +0.000e+00 +8.457e+02 --
965 155020 74494 +0.000e+00 +8.457e+02 --
970 156280 74434 +0.000e+00 +8.457e+02 --
975 157462 74682 +0.000e+00 +8.457e+02 --
980 158482 75013 +0.000e+00 +8.457e+02 --
985 159801 75193 +0.000e+00 +8.457e+02 --
990 161233 75520 +0.000e+00 +8.457e+02 --
995 162321 75697 +0.000e+00 +8.457e+02 --
Reached time limit of 1000 CPU s
1000 163199 75922 +0.000e+00 +8.457e+02 --
-------------------------------------------------------------------------------
Termination Status : Found feasible point; reached time limit
Best Feasible Point: +8.457215e+02
Best Possible Point: +0.000000e+00
Relative Gap: +1.000000e+00
Algorithm analysis :
163199 Nodes explored
75922 Nodes remaining
23 Maximum tree depth
24394 Cutting Planes
24394 Edge-Concave
1000.00 Total time (wall s)
0.00 Pre-processing
406.37 Solving MILP relaxations
8.71 Searching for feasible solutions
580.09 Variable bounds tightening
3.57 OBBT
577.47 FBBT (26.1 EC; 195.8 RLT; 21.0 Factoring)
527.50 Branching
0.59 Reliability branching
This is an excellent small benchmark model. I think this is the most stubborn, tiny, quadratic model I have ever seen.
Question: is it worthwhile to run this model for a day or so and see what happens? I have my doubts.
Can we do something about this?
Instead of just giving up, let's work on this problem a bit further. We don't know where the real underlying problem is (is it the minimizing area or the point inclusion constraint).
The first experiment is unlikely to make difference. But, let's make sure. We can use a different method to implement the absolute value. We can form \[\begin{align}\min \>&{\color{darkred}A}\\ & -{\color{darkred}A} \le 0.5 \cdot \left[{\color{darkred}t}_{1,x}({\color{darkred}t}_{2,y}-{\color{darkred}t}_{3,y}) + {\color{darkred}t}_{2,x}({\color{darkred}t}_{3,y}-{\color{darkred}t}_{1,y}) + {\color{darkred}t}_{3,x}({\color{darkred}t}_{1,y}-{\color{darkred}t}_{2,y}) \right] \le {\color{darkred}A} \end{align}\] I tried this out, and, as expected, it made no difference: the solvers have difficulties moving the lower bound.
One other thing we can try, is to use a bit of a different definition for the size of the triangle. In the previous attempt we used the area, which is the most natural measure for size. But we could opt for something else, e.g. circumference \[\min\>\sum_{p,q|p\lt q} \sqrt{({\color{darkred}t}_{p,x}-{\color{darkred}t}_{q,x})^2+({\color{darkred}t}_{p,y}-{\color{darkred}t}_{q,y})^2}\] where we sum over all combinations of two corner points (there are just 3 of them). This actually gives much better performance: the lower bound starts to move quickly.
We can even take a step further, and remove the square root. Now the objective is just a (convex) sum of squares: \[\min\>\sum_{p,q|p\lt q} \left[ ({\color{darkred}t}_{p,x}-{\color{darkred}t}_{q,x})^2+({\color{darkred}t}_{p,y}-{\color{darkred}t}_{q,y})^2\right]\] This improves performance even further.
I tried all four models:
- Original: minimize area using variable splitting
- Minimize area using bounding
- Minimize circumference
- Minimize sum of squares
with SCIP on a data set with 50 random points. This gave:
---- 173 PARAMETER results combined results area time m1 area var split 10460.682 1000.000 m2 area bnd 10460.682 1000.000 m3 circumference 10470.137 133.000 m4 sos 10519.844 117.000
The first two models were stopped on the time limit of 1,000 seconds. The last two models were proven optimal (or near optimal and stopped on very small gap — GAMS, by default, uses a gap tolerance of 0.01%).
The triangles are very close.
GAMS Model
$onText
Given are n data points (2d). Find the smallest triangle that contains all points. Three different definitions of size. $offText
option nlp=scip, reslim=1000; *--------------------------------------------------------------------- * data: points *--------------------------------------------------------------------- sets i 'points' /point1*point50/ c 'coordinates' /x,y/ ; parameter p(i,c) 'data points'; p(i,c) = uniform(0,100);
*--------------------------------------------------------------------- * constraints: * points inside triangle using barymetric coordinates * order corner points by x coordinate *---------------------------------------------------------------------
set k 'corner points of triangle' /corner1*corner3/ ;
positive variable lambda(i,k) 'barycentric coordinates' ;
free variable t(k,c) 'triangle' ; * some reasonable bounds t.lo(k,c) = -1000; t.up(k,c) = +1000;
equations calcLambda(i,c) 'solve for barycentric coordinates' sumLambda(i) 'lambdas need to add up to one' order 'order corner points by their x coordinate' ;
calcLambda(i,c).. p(i,c) =e= sum(k, lambda(i,k)*t(k,c)); sumLambda(i).. sum(k, lambda(i,k)) =e= 1; order(k-1).. t(k,'x') =g= t(k-1,'x');
model cons 'constraints' /calcLambda,sumLambda,order/;
*--------------------------------------------------------------------- * objective 1: area (variable splitting) *---------------------------------------------------------------------
parameter results(*,*) 'combined results' triangle(*,k,c) 'solution' ;
set pm 'plusmin -- used in linearizing abs()' /'+','-'/ ; * shorthands to make our area calculation easier singleton sets x1(k,c) /'corner1'.'x'/ x2(k,c) /'corner2'.'x'/ x3(k,c) /'corner3'.'x'/ y1(k,c) /'corner1'.'y'/ y2(k,c) /'corner2'.'y'/ y3(k,c) /'corner3'.'y'/ ; variable z 'objective' ; positive variable area(pm) 'area (using variable splitting)' ; equations calcArea 'calculate area given its three corner points' obj 'objective' ; calcArea.. area('+')-area('-') =e= 0.5*[t(x1)*(t(y2)-t(y3)) + t(x2)*(t(y3)-t(y1)) + t(x3)*(t(y1)-t(y2))]; obj.. z =e= sum(pm,area(pm)); model m1 'area -- var splitting' /calcArea,obj,cons/; solve m1 minimizing z using nlp;
results('m1 area var split','area') = abs(0.5*[t.l(x1)*(t.l(y2)-t.l(y3)) + t.l(x2)*(t.l(y3)-t.l(y1)) + t.l(x3)*(t.l(y1)-t.l(y2))]); results('m1 area var split','time') = m1.resusd; triangle('m1',k,c) = t.l(k,c); display results,triangle;
*--------------------------------------------------------------------- * objective 2: area (bounding) *---------------------------------------------------------------------
variable a 'area' absa 'abs area' ; absa.lo = 0;
equations calcArea2 'calculate area given its three corner points' bnd1 'objective bound' bnd2 'objective bound' ;
calcArea2.. a =e= 0.5*[t(x1)*(t(y2)-t(y3)) + t(x2)*(t(y3)-t(y1)) + t(x3)*(t(y1)-t(y2))]; bnd1.. a =l= absa; bnd2.. a =g= -absa;
model m2 'area -- bounding' /calcArea2,bnd1,bnd2,cons/; solve m2 minimizing absa using nlp;
results('m2 area bnd','area') = abs(0.5*[t.l(x1)*(t.l(y2)-t.l(y3)) + t.l(x2)*(t.l(y3)-t.l(y1)) + t.l(x3)*(t.l(y1)-t.l(y2))]); results('m2 area bnd','time') = m2.resusd; triangle('m2',k,c) = t.l(k,c); display results,triangle;
*--------------------------------------------------------------------- * objective 3: circumference *---------------------------------------------------------------------
alias(k,kk);
equations circumference 'calculate perimeter of triangle' ;
circumference.. z =e= sum((k,kk)$(ord(k)<ord(kk)), sqrt(sum(c,sqr(t(k,c)-t(kk,c)))));
model m3 'area -- bounding' /circumference,cons/; solve m3 minimizing z using nlp;
results('m3 circumference','area') = abs(0.5*[t.l(x1)*(t.l(y2)-t.l(y3)) + t.l(x2)*(t.l(y3)-t.l(y1)) + t.l(x3)*(t.l(y1)-t.l(y2))]); results('m3 circumference','time') = m3.resusd; triangle('m3',k,c) = t.l(k,c); display results,triangle;
*--------------------------------------------------------------------- * objective 4: sum of squares *---------------------------------------------------------------------
alias(k,kk);
equations sos 'drop square root' ;
sos.. z =e= sum((k,kk)$(ord(k)<ord(kk)), sum(c,sqr(t(k,c)-t(kk,c))));
model m4 'area -- bounding' /sos,cons/; solve m4 minimizing z using nlp;
results('m4 sos','area') = abs(0.5*[t.l(x1)*(t.l(y2)-t.l(y3)) + t.l(x2)*(t.l(y3)-t.l(y1)) + t.l(x3)*(t.l(y1)-t.l(y2))]); results('m4 sos','time') = m4.resusd; triangle('m4',k,c) = t.l(k,c); display results,triangle;
*--------------------------------------------------------------------- * visualization *---------------------------------------------------------------------
$set html plot.html $set data data.js
file fdata /%data%/; put fdata;
* points put "px=["; loop(i,put p(i,'x'):0:3,","); put "];"/; put "py=["; loop(i,put p(i,'y'):0:3,","); put "];"/; * triangles put "txm1=["; loop(k,put triangle('m1',k,'x'):0:3,","); put "];"/; put "tym1=["; loop(k,put triangle('m1',k,'y'):0:3,","); put "];"/; put "txm3=["; loop(k,put triangle('m3',k,'x'):0:3,","); put "];"/; put "tym3=["; loop(k,put triangle('m3',k,'y'):0:3,","); put "];"/; put "txm4=["; loop(k,put triangle('m4',k,'x'):0:3,","); put "];"/; put "tym4=["; loop(k,put triangle('m4',k,'y'):0:3,","); put "];"/;
set mall /m1,m3,m4/;
parameter box; box(c,'lo') = smin((mall,k),triangle(mall,k,c)); box(c,'up') = smax((mall,k),triangle(mall,k,c)); box(c,'range') = box(c,'up')-box(c,'lo'); box(c,'lo2') = box(c,'lo') - 0.1*box(c,'range'); box(c,'up2') = box(c,'up') + 0.1*box(c,'range'); display box;
put "cmin=["; loop(c, put box(c,'lo2'):0:3,","); put "];"/; put "cmax=["; loop(c, put box(c,'up2'):0:3,","); put "];"/;
$onecho > %html% <html> <script src="https://cdn.plot.ly/plotly-3.4.0.min.js" charset="utf-8"></script> <script src="%data%" charset="utf-8"></script> <h1>Smallest Encompassing Triangle</h1> <div id="plotDiv"></div> <script> var data = { x: px, y: py, mode: 'markers', type: 'scatter', name: 'data points' };
txm1[3]=txm1[0]; tym1[3]=tym1[0]; var minarea = { x: txm1, y: tym1, mode: 'lines', type: 'scatter', name: 'min area' };
txm3[3]=txm3[0]; tym3[3]=tym3[0]; var mincircumf = { x: txm3, y: tym3, mode: 'lines', type: 'scatter', name: 'min circumference' };
txm4[3]=txm4[0]; tym4[3]=tym4[0]; var minsos = { x: txm4, y: tym4, mode: 'lines', type: 'scatter', name: 'min sum of squares' };
var layout = { autosize: false, width: 700, height: 500, showlegend: true } var trc = [data,minarea,mincircumf,minsos]; var options = {staticPlot: true, displayModeBar: false}; Plotly.newPlot('plotDiv', trc, layout, options); </script> </html> $offecho
executetool 'win32.ShellExecute "%html%"'; |
Conclusions
Sometimes when a model cannot be solved within acceptable time limits, it is a good idea to solve something different. Here we redefine what the size of a triangle means. Using this freedom, we can solve the previously extremely difficult models quite easily.
This was a very interesting exercise.
References
- Tiny non-convex quadratic model brings solvers to their knees, https://yetanothermathprogrammingconsultant.blogspot.com/2023/01/tiny-non-convex-quadratic-model-brings.html
- Barycentric coordinate system, https://en.wikipedia.org/wiki/Barycentric_coordinate_system
This is another example where a trivially true geometric fact is quite hard to derive algebraically for a specific implementation of the model. It would be a nice exercise to prove algebraically that the objective value of the model is lower bounded by the area of any triangle in the point set. Once you have such an algebraic proof it can be turned into a presolver or cut generation technique to strengthen the nonlinear solvers.
ReplyDeleteInteresting! Though, once changed the definition of area, I it is strange to claim optimal solutions, given the definitions of areas are different, and therefore, the objectives are different. But indeed, if there is a "flexible" definition of area, improvements can be made "outside the box". I wonder if there are different definitions of area that are provably lower bounds of the "actual" area. This may improve lower bounds.
ReplyDelete