Yet Another Math Programming Consultant: Largest Empty Rotated Square: lots of trigonometry

Here, I delve into the problem of finding the largest empty rotated square. Given $n$ data points, find the largest square rotated by an angle $\theta$, such that the square does not contain any of the points. In [1], I discussed two special cases:

an axis-aligned square (angle is 0°),
a diamond shape, which is a square rotated by 45°.

Here we focus on arbitrary angles. A picture can help a bit: this is the largest square of all the squares with a 25° angle, which does not cover any of the given data points and is complete inside the outer box:

There are two cases I want to consider:

The angle is given. In modeling terms, the angle is exogenous.
Let the model find the best angle. I.e., the angle is endogenous.

Model 1. Exogenous Angle

Here, we assume we know the rotation angle ${\color{darkblue}\theta}$ in advance. The idea behind the approach is as follows. We need to check whether a point is outside the rotated square. That is a bit difficult. So, instead, we rotate all data points by $-{\color{darkblue}\theta}$. Now the square is axis-aligned and we can use the techniques from [1]. Well, almost: a complication is that the square must be inside the outer box.

We use the same data set as in [1]:

Instead of rotating our square, we rotate the data. That will result in something like:

Here, the angle is 25°.

Now, we can simply compare the new points against an axis-aligned square. This rotation by $-{\color{darkblue}\theta}$ is implemented as follows: \[\begin{align} & {\color{darkblue}p}'_{i,x} := {\color{darkblue}p}_{i,x} \cos({\color{darkblue}\theta}) + {\color{darkblue}p}_{i,y} \sin({\color{darkblue}\theta}) \\ & {\color{darkblue}p}'_{i,y} := -{\color{darkblue}p}_{i,x} \sin({\color{darkblue}\theta}) + {\color{darkblue}p}_{i,y} \cos({\color{darkblue}\theta}) \end{align}\] These formulas can be easily derived from the standard rotation matrix [2] and using $\sin(-x)=-\sin(x)$ and $\cos(-x)=\cos(x)$. Now we can just try to find a large standard square, as documented in [1]. I want to emphasize that the calculation of ${\color{darkblue}p}'_{i,c}, c \in \{x,y\}$ is just a data-manipulation step; it is an assignment and not a model equation.

One complication is that we need to ensure the square is within our outer box. This is now much more complicated. My approach is as follows. We need to construct all corner points (that is easy), and then rotate them back to our original coordinate system. Let's denote by ${\color{darkred}q}_{j,c}'$ the coordinates of the four corner points in the transformed space (so $j=1,\dots,4$). Then rotate them back to ${\color{darkred}q}_{j,c}$. Now the constraints are simple again: \[ 0 \le {\color{darkred}q}_{j,c} \le {\color{darkblue}{\mathit{maxSize}}}\].

A complete model can look like:

Exogenous Angle Model
\[\begin{align} \max\> & {\color{darkred}s} \\&{\color{darkblue}p}'_{i,x} \le {\color{darkred}x}' + {\color{darkblue}M} {\color{darkred}\delta}_{i,1} && \forall i \\&{\color{darkblue}p}'_{i,x} \ge {\color{darkred}x}' + {\color{darkred}s} - {\color{darkblue}M} {\color{darkred}\delta}_{i,2} && \forall i \\&{\color{darkblue}p}'_{i,y} \le {\color{darkred}y}' + {\color{darkblue}M} {\color{darkred}\delta}_{i,3} && \forall i \\&{\color{darkblue}p}'_{i,y} \ge {\color{darkred}y}' + {\color{darkred}s} - {\color{darkblue}M} {\color{darkred}\delta}_{i,4} && \forall i \\ & \sum_k {\color{darkred}\delta}_{i,k} = 3 && \forall i \\&{\color{darkred}q}_{j,x}' = {\color{darkred}x}' + {\color{darkred}s} {\color{darkblue}\Delta}_{j,x} && \forall j \\&{\color{darkred}q}_{j,y}' = {\color{darkred}y}' + {\color{darkred}s} {\color{darkblue}\Delta}_{j,y} && \forall j \\&{\color{darkred}q}_{j,x} ={\color{darkred}q}_{j,x}' \cos({\color{darkblue}\theta}) - {\color{darkred}q}_{j,y}' \sin({\color{darkblue}\theta}) && \forall j \\&{\color{darkred}q}_{j,y} ={\color{darkred}q}_{j,x}' \sin({\color{darkblue}\theta}) + {\color{darkred}q}_{j,y}' \cos({\color{darkblue}\theta}) && \forall j \\&{\color{darkred}q}_{j,c} \in [0, {\color{darkblue}{\mathit{maxSize}}}] \\&{\color{darkred}s} \in [0,{\color{darkblue}{\mathit{maxSize}}}] \\&{\color{darkred}\delta}_{i,k} \in \{0,1\} \\&{\color{darkblue}p}'_{i,x} := {\color{darkblue}p}_{i,x} \cos({\color{darkblue}\theta}) + {\color{darkblue}p}_{i,y} \sin({\color{darkblue}\theta}) \\&{\color{darkblue}p}'_{i,y} := -{\color{darkblue}p}_{i,x} \sin({\color{darkblue}\theta}) + {\color{darkblue}p}_{i,y} \cos({\color{darkblue}\theta}) \\&{\color{darkblue}\Delta} = \begin{pmatrix} 0 & 0 \\ 0 & 1 \\ 1 & 0 \\ 1 & 1\end{pmatrix} \\ & i\in \{1,\dots,n\} \\ & c \in \{x,y\} \\& k \in \{1,\dots,4\}\\& j \in \{1,\dots,4\} \end{align}\]

Exogenous Angle Model

\[\begin{align} \max\> & {\color{darkred}s}
\\&{\color{darkblue}p}'_{i,x} \le {\color{darkred}x}' + {\color{darkblue}M} {\color{darkred}\delta}_{i,1} && \forall i
\\&{\color{darkblue}p}'_{i,x} \ge {\color{darkred}x}' + {\color{darkred}s} - {\color{darkblue}M} {\color{darkred}\delta}_{i,2} && \forall i
\\&{\color{darkblue}p}'_{i,y} \le {\color{darkred}y}' + {\color{darkblue}M} {\color{darkred}\delta}_{i,3} && \forall i
\\&{\color{darkblue}p}'_{i,y} \ge {\color{darkred}y}' + {\color{darkred}s} - {\color{darkblue}M} {\color{darkred}\delta}_{i,4} && \forall i \\ & \sum_k {\color{darkred}\delta}_{i,k} = 3 && \forall i
\\&{\color{darkred}q}_{j,x}' = {\color{darkred}x}' + {\color{darkred}s} {\color{darkblue}\Delta}_{j,x} && \forall j
\\&{\color{darkred}q}_{j,y}' = {\color{darkred}y}' + {\color{darkred}s} {\color{darkblue}\Delta}_{j,y} && \forall j
\\&{\color{darkred}q}_{j,x} ={\color{darkred}q}_{j,x}' \cos({\color{darkblue}\theta}) - {\color{darkred}q}_{j,y}' \sin({\color{darkblue}\theta}) && \forall j
\\&{\color{darkred}q}_{j,y} ={\color{darkred}q}_{j,x}' \sin({\color{darkblue}\theta}) + {\color{darkred}q}_{j,y}' \cos({\color{darkblue}\theta}) && \forall j
\\&{\color{darkred}q}_{j,c} \in [0, {\color{darkblue}{\mathit{maxSize}}}]
\\&{\color{darkred}s} \in [0,{\color{darkblue}{\mathit{maxSize}}}]
\\&{\color{darkred}\delta}_{i,k} \in \{0,1\}
\\&{\color{darkblue}p}'_{i,x} := {\color{darkblue}p}_{i,x} \cos({\color{darkblue}\theta}) + {\color{darkblue}p}_{i,y} \sin({\color{darkblue}\theta})
\\&{\color{darkblue}p}'_{i,y} := -{\color{darkblue}p}_{i,x} \sin({\color{darkblue}\theta}) + {\color{darkblue}p}_{i,y} \cos({\color{darkblue}\theta})
\\&{\color{darkblue}\Delta} = \begin{pmatrix} 0 & 0 \\ 0 & 1 \\ 1 & 0 \\ 1 & 1\end{pmatrix} \\ & i\in \{1,\dots,n\} \\ & c \in \{x,y\} \\& k \in \{1,\dots,4\}\\& j \in \{1,\dots,4\} \end{align}\]

Notes:

$\sin({\color{darkblue}\theta})$ and $\cos({\color{darkblue}\theta})$ are constants. Unbelievably, this is actually a linear MIP model!
${\color{darkred}x}',{\color{darkred}y}'$ denotes the lower-left corner of the axis-aligned square (like in [1]).
${\color{darkred}q}_{j,c}', j=1,\dots,4$ are the four corner points (in transformed space). We map them back to ${\color{darkred}q}_{j,c}$, which form a rotated square.
We can tighten the $\color{darkblue}M$ values considerable by making them dependent on the point, i.e., ${\color{darkblue}M}_{i,k}$.
Below, we plot the solution both in the transformed space and in the original coordinate system.
This is a very easy MIP:

Proven optimal solution
MIP Solution: 3.085938 (657 iterations, 0 nodes)

The solution in transformed coordinates looks like:

When we transform this back to our original coordinate system, we see:

This solution is slightly better (i.e., we have a larger optimal square) than the simple axis-aligned square in [1].

Obviously, we already have two good test cases for this model from [1]: 0° and 45°.

Next thing to try: make $\theta$ endogenous.

GAMS Model

$onText

Largest Empty Rotated Square.

This version: exogenous angle

This gives us a MIP model.

$offText

*---------------------------------------------------------------

* options

*---------------------------------------------------------------

* 0: no plot

* 1: produce HTML plot

$set htmlplot 1

* this model turned to linear

option mip=cplex;

option seed=12345;

*---------------------------------------------------------------

* data

*---------------------------------------------------------------

scalar maxsize 'size of the box' / 10 /;

sets

i 'points' /point1*point100/

a 'square attributes' /x,y,s/

c(a) 'coordinates' /x,y/

;

parameter p(i,c) 'points';

p(i,c) = uniform(0,maxsize);

display p;

*---------------------------------------------------------------

* MIP MODEL

*---------------------------------------------------------------

scalars

angle 'degrees' / 25 /

theta 'angle in radians'

M 'big-M'

;

M = 2 * maxSize;

theta = angle*pi/180;

display theta;

* transform data points

parameter q(i,c) 'transformed data points: rotate by -theta';

q(i,'x') = p(i,'x')*cos(theta) + p(i,'y')*sin(theta);

q(i,'y') = -p(i,'x')*sin(theta) + p(i,'y')*cos(theta);

display q;

* transform corner points (for picture)

* not used in the model, but useful for visualization

set k 'corner points box' /0-0, 0-1, 1-0, 1-1/;

parameter box(k,c) 'box corners';

box('0-0','x') = 0; box('0-0','y') = 0;

box('0-1','x') = 0; box('0-1','y') = maxsize;

box('1-0','x') = maxsize; box('1-0','y') = 0;

box('1-1','x') = maxsize; box('1-1','y') = maxsize;

display box;

parameter tbox(k,c) 'box corners (transformed)';

tbox(k,'x') = box(k,'x')*cos(theta) + box(k,'y')*sin(theta);

tbox(k,'y') = -box(k,'x')*sin(theta) + box(k,'y')*cos(theta);

display tbox;

set

case 'for comparison' /less,greater/

corners /00,01,10,11/

;

table offset(corners,c) 'unit offset of corner points from lower-left'

x y

00 0 0

01 0 1

10 1 0

11 1 1

;

variable

sq(a) 'square location and size (transformed coordinates)'

sq2(corners,c) 'corner points (original coordinates)'

;

sq2.lo(corners,c) = 0;

sq2.up(corners,c) = maxsize;

binary variable d(i,c,case) 'd=1 means: relax no-overlap constraint';

variable z 'objective';

Equations

emptysq1(i,c) 'no-overlap constraint: less-than version'

emptysq2(i,c) 'no-overlap constraint: greater-than version'

sumd(i) '(at least) one no-overlap constraint must hold'

transformx "transform x',y' -> x back to original coordinates"

transformy "transform x',y' -> y back to original coordinates"

objsq 'objective'

;

* this is like the original axis-aligned square model

objsq.. z =e= sq('s');

emptysq1(i,c).. q(i,c) =l= sq(c) + d(i,c,'less')*M;

emptysq2(i,c).. q(i,c) =g= sq(c) + sq('s') - d(i,c,'greater')*M;

sumd(i).. sum((c,case),d(i,c,case)) =e= 3;

* transform back to original coordinates

transformx(corners).. sq2(corners,'x') =e= (sq('x')+offset(corners,'x')*sq('s'))*cos(theta) - (sq('y')+offset(corners,'y')*sq('s'))*sin(theta);

transformy(corners).. sq2(corners,'y') =e= (sq('x')+offset(corners,'x')*sq('s'))*sin(theta) + (sq('y')+offset(corners,'y')*sq('s'))*cos(theta);

model emptysquare /all/;

solve emptysquare maximizing z using mip;

parameter square(*,*) 'optimal values';

square(corners,c) = sq2.l(corners,c);

square('side','-') = sq.l('s');

square('area','-') = sqr(square('side','-'));

display square;

*---------------------------------------------------------------------

* visualization

*---------------------------------------------------------------------

$set html plot.html

$set data data.js

$if %htmlplot%==0 $goto skipplot

file fdata /%data%/; put fdata;

* points

put "p=["/;

loop(i,

put " {i:'",i.tl:0,"',x:",p(i,'x'):0:4,",y:",p(i,'y'):0:4,"},"/;

);

put "]"/;

put "q=["/;

loop(i,

put " {i:'",i.tl:0,"',x:",q(i,'x'):0:4,",y:",q(i,'y'):0:4,"},"/;

);

put "]"/;

put "box={x:0,y:0,size:",maxsize:0:4,"}"/;

put "angle=",angle:0:4/;

put "tbox=["/;

loop(k,

put " {k:'",k.tl:0,"',x:",tbox(k,'x'):0:4,",y:",tbox(k,'y'):0:4,"},"/;

);

put "]"/;

put "tr_sq={x:",sq.l('x'):0:4,",y:",sq.l('y'):0:4,",side:",sq.l('s'):0:4,",area:",(sq.l('s')**2):0:4,"}"/;

put "cp=["/;

loop(corners,

put " {k:'",corners.tl:0,"',x:",square(corners,'x'):0:4,",y:",square(corners,'y'):0:4,"},"/;

);

put "]"/;

putclose;

$onecho > %html%

<html>

<style>

.brdr {

border: 1px solid black;

border-collapse: collapse;

}

th, td {

padding-left: 4px;

padding-right: 4px;

}

</style>

<table>

</td><td>

</td></tr></table>

<h1>Solution in Transformed Coordinates</h1>

<tr><th colspan=2 class="brdr">Square</td></tr>

</table>

</td><td>

</td></tr></table>

<h1>Solution in Original Coordinates</h1>

<tr><th colspan=3 class="brdr">Rotated Square</td></tr>

<tr><td class="brdr">Corner 1</td><td class="brdr" id="c00_x"></td><td class="brdr" id="c00_y"></td></tr>

<tr><td class="brdr">Corner 2</td><td class="brdr" id="c01_x"></td><td class="brdr" id="c01_y"></td></tr>

<tr><td class="brdr">Corner 3</td><td class="brdr" id="c10_x"></td><td class="brdr" id="c10_y"></td></tr>

<tr><td class="brdr">Corner 4</td><td class="brdr" id="c11_x"></td><td class="brdr" id="c11_y"></td></tr>

</table>

</td><td>

</td></tr></table>

px = p.map((x) => x['x'])

py = p.map((x) => x['y'])

qx = q.map((x) => x['x'])

qy = q.map((x) => x['y'])

document.getElementById("tr_x").innerHTML = tr_sq['x'];

document.getElementById("tr_y").innerHTML = tr_sq['y'];

document.getElementById("tr_s").innerHTML = tr_sq['side'];

document.getElementById("tr_area").innerHTML = tr_sq['area'];

document.getElementById("c00_x").innerHTML = cp[0]['x'];

document.getElementById("c00_y").innerHTML = cp[0]['y'];

document.getElementById("c01_x").innerHTML = cp[1]['x'];

document.getElementById("c01_y").innerHTML = cp[1]['y'];

document.getElementById("c10_x").innerHTML = cp[2]['x'];

document.getElementById("c10_y").innerHTML = cp[2]['y'];

document.getElementById("c11_x").innerHTML = cp[3]['x'];

document.getElementById("c11_y").innerHTML = cp[3]['y'];

var path2 = `M ${cp[0]['x']} ${cp[0]['y']} L ${cp[1]['x']} ${cp[1]['y']} L ${cp[3]['x']} ${cp[3]['y']} L ${cp[2]['x']} ${cp[2]['y']} Z`;

var data1 = {

x: px,

y: py,

mode: 'markers',

type: 'scatter',

color: 'darkblue'

};

var data2 = {

x: qx,

y: qy,

mode: 'markers',

type: 'scatter',

color: 'darkblue'

};

var layout1 = {

autosize: false,

width: 650,

height: 600,

showlegend: true,

title: {text:'Original Data'},

shapes: [{type:'rect',

xref:'x',

yref:'y',

x0:box['x'],

y0:box['y'],

x1:box['x']+box['size'],

y1:box['y']+box['size'],

line:{color:'black',width:2,opacity:0.5},

fillcolor:'rgba(0, 250, 250, 0.2)',

name:'box',

showlegend:true

}

]

}

var path = `M ${tbox[0]['x']} ${tbox[0]['y']} L ${tbox[1]['x']} ${tbox[1]['y']} L ${tbox[3]['x']} ${tbox[3]['y']} L ${tbox[2]['x']} ${tbox[2]['y']} Z`;

var layout2 = {

autosize: false,

width: 650,

height: 600,

showlegend: true,

title: {text:'Transformed Data'},

shapes: [{type:'path',

path:path,

xref:'x',

yref:'y',

line:{color:'black',width:2,opacity:0.5},

fillcolor:'rgba(0, 250, 250, 0.2)',

name:'transformed box',

showlegend:true

}

]

}

var layout3 = {

autosize: false,

width: 650,

height: 600,

showlegend: true,

title: {text:'Transformed Solution'},

shapes: [

{type:'path',

path:path,

xref:'x',

yref:'y',

line:{color:'black',width:2,opacity:0.5},

fillcolor:'rgba(0, 250, 250, 0.2)',

name:'transformed box',

showlegend:true

{type:'rect',

xref:'x',

yref:'y',

x0:tr_sq['x'],

y0:tr_sq['y'],

x1:tr_sq['x']+tr_sq['side'],

y1:tr_sq['y']+tr_sq['side'],

line:{color:'darkred',width:2,opacity:0.5},

fillcolor:'rgba(0, 150, 150, 0.5)',

name:'optimal square',

showlegend:true

]

}

var layout4 = {

autosize: false,

width: 650,

height: 600,

showlegend: true,

title: {text:'Solution'},

shapes: [{type:'rect',

xref:'x',

yref:'y',

x0:box['x'],

y0:box['y'],

x1:box['x']+box['size'],

y1:box['y']+box['size'],

line:{color:'black',width:2,opacity:0.5},

fillcolor:'rgba(0, 250, 250, 0.2)',

name:'box',

showlegend:true

{type:'path',

path:path2,

xref:'x',

yref:'y',

line:{color:'darkred',width:2,opacity:0.5},

fillcolor:'rgba(0, 150, 150, 0.5)',

name:'optimal square',

showlegend:true

]

}

var trc1 = [data1];

var trc2 = [data2];

var options1 = {staticPlot: true, displayModeBar: false, responsive: false};

Plotly.newPlot('plotDiv0', trc1, layout1, options1);

Plotly.newPlot('plotDiv1', trc2, layout2, options1);

Plotly.newPlot('plotDiv2', trc2, layout3, options1);

Plotly.newPlot('plotDiv3', trc1, layout4, options1);

</script>

$offecho

executetool 'win32.ShellExecute "%html%"';

$label skipplot

Model 2. Endogenous Angle

This is essentially the same as the previous model, except:

${\color{darkred}{\theta}}$ is now a variable. We allow it to vary between 0 and 90°, or rather between 0 and $\pi/2$ radians.
${\color{darkred}p}'_{i,c}$ is now a variable. The assignments become constraints: \[\begin{align} & {\color{darkred}p}'_{i,x} = {\color{darkblue}p}_{i,x} \cos({\color{darkred}\theta}) + {\color{darkblue}p}_{i,y} \sin({\color{darkred}\theta}) \\ & {\color{darkred}p}'_{i,y} = -{\color{darkblue}p}_{i,x} \sin({\color{darkred}\theta}) + {\color{darkblue}p}_{i,y} \cos({\color{darkred}\theta}) \end{align}\]

These are simple and quick edits, so I won't reproduce the complete model. The model is now a nonconvex MINLP model. Not all global MINLP solvers support the trigonometric functions like $\sin(.)$ and $\cos(.)$: Baron and Antigone are examples of that. We can try another global MINLP solver such as SCIP, Gurobi or Xpress. SCIP is struggling a bit with this model, but Gurobi and Xpress have no problem with it and solve it to (proven) global optimality very quickly (less than 10 seconds):

Explored 6106 nodes (93842 simplex iterations) in 5.16 seconds (5.34 work units)

Optimal solution found (tolerance 1.00e-04)
Best objective 3.141077560753e+00, best bound 3.141077560753e+00, gap 0.0000%

The solution looks like:

----    136 PARAMETER square  optimal values

                      x           y           -

00                4.071       4.434
01                2.659       7.240
10                6.877       5.847
11                5.464       8.652
side                                      3.141
area                                      9.866
angle (rad)                               0.466
angle (deg)                              26.728

Although the angle is very close to our earlier run with 25°, the area is significantly better. The fact that a degree or two can make a significant difference in the objective probably means that running a number of the exogenous angle models with different rotation values is not a good idea. Similarly, I have doubts about using linear interpolation between different angles.

It is noted that calculating better values for the big-M constants is now more difficult than when we know the angle in advance (the exogenous angle case).

A great trick (see the comments) is to replace $\sin({\color{darkred}\theta})$ and $\cos({\color{darkred}\theta})$ by variables ${\color{darkred}{\alpha}}$ and ${\color{darkred}{\beta}}$ (both with bounds $[0,1]$), and introducing the constraint \[{\color{darkred}{\alpha}}^2 +{\color{darkred}{\beta}}^2 = 1\] This will make the model a non-convex MIQCP. We can recover the angle afterwards, e.g. by \[{\color{darkred}\theta} := \arctan({\color{darkred}{\alpha}}/{\color{darkred}{\beta}})\] To make this safe, we would have to change the lower bound ${\color{darkred}{\beta}}\ge 0.00001$. That is better anyway, as it reduces some symmetry: the solutions with angles 0° and 90° are the same. So it is better to exclude the 90° solution.

Not only does this reformulation allow more solvers to be used, but it also seems a bit faster:

Explored 3276 nodes (41404 simplex iterations) in 2.88 seconds (2.55 work units)

Optimal solution found (tolerance 1.00e-04)
Best objective 3.141077907762e+00, best bound 3.141077907762e+00, gap 0.0000%

GAMS Models

$onText

Largest Empty Rotated Square.

This version: endogenous angle

This gives us a MINLP or MIQCP model.

$offText

*---------------------------------------------------------------

* options

*---------------------------------------------------------------

* 0: no plot

* 1: produce HTML plot

$set htmlplot 1

* baron, antigone don't have sin(),cos() functions, so we use xpress and gurobi for this example

option minlp=xpress,miqcp=xpress;

*option minlp=gurobi,miqcp=gurobi;

option seed=12345;

*---------------------------------------------------------------

* data

*---------------------------------------------------------------

scalar maxsize 'size of the outer box' / 10 /;

sets

i 'points' /point1*point100/

a 'square attributes' /x,y,s/

c(a) 'coordinates' /x,y/

;

parameter p(i,c) 'points';

p(i,c) = uniform(0,maxsize);

display p;

*---------------------------------------------------------------

* MINLP model with angles

*---------------------------------------------------------------

scalars M 'big-M';

M = 2 * maxSize;

set

case 'for comparison' /less,greater/

corners 'enumerate four corners' /00,01,10,11/

;

table offset(corners,c) 'unit offset of corner points from lower-left'

x y

00 0 0

01 0 1

10 1 0

11 1 1

;

variable

q(i,c) 'transformed data points: rotate by -theta'

sq(a) 'square location and size (transformed coordinates)'

sq2(corners,c) 'corner points (original coordinates)'

theta 'angle of rotation'

;

q.lo(i,c) = -2*maxsize;

q.up(i,c) = 2*maxsize;

theta.lo = 0;

theta.up = 90*pi/180 - 0.00001;

sq.lo(a) = -2*maxsize;

sq.up(a) = 2*maxsize;

sq2.lo(corners,c) = 0;

sq2.up(corners,c) = maxsize;

binary variable d(i,c,case) 'd=1 means: relax no-overlap constraint';

variable z 'objective';

Equations

transfdatax(i) 'transform data points: x coordinate'

transfdatay(i) 'transform data points: y coordinate'

emptysq1(i,c) 'no-overlap constraint: less-than version'

emptysq2(i,c) 'no-overlap constraint: greater-than version'

sumd(i) '(at least) one no-overlap constraint must hold'

transformx(corners) "transform x',y' -> x back to original coordinates"

transformy(corners) "transform x',y' -> y back to original coordinates"

objsq 'objective'

;

* transform data points

transfdatax(i).. q(i,'x') =e= p(i,'x')*cos(theta) + p(i,'y')*sin(theta);

transfdatay(i).. q(i,'y') =e= -p(i,'x')*sin(theta) + p(i,'y')*cos(theta);

* this is like the original axis-aligned square model

objsq.. z =e= sq('s');

emptysq1(i,c).. q(i,c) =l= sq(c) + d(i,c,'less')*M;

emptysq2(i,c).. q(i,c) =g= sq(c) + sq('s') - d(i,c,'greater')*M;

sumd(i).. sum((c,case),d(i,c,case)) =e= 3;

* transform back to original coordinates

transformx(corners).. sq2(corners,'x') =e= (sq('x')+offset(corners,'x')*sq('s'))*cos(theta) - (sq('y')+offset(corners,'y')*sq('s'))*sin(theta);

transformy(corners).. sq2(corners,'y') =e= (sq('x')+offset(corners,'x')*sq('s'))*sin(theta) + (sq('y')+offset(corners,'y')*sq('s'))*cos(theta);

* to simulate earlier runs

*theta.fx = 0;

*theta.fx = 25*pi/180;

*theta.fx = 45*pi/180;

model emptysquare1 /all/;

solve emptysquare1 maximizing z using minlp;

display theta.l, sq.l, sq2.l;

parameter m1results(*,*) 'optimal values model 1';

m1results(corners,c) = sq2.l(corners,c);

m1results('side','-') = sq.l('s');

m1results('area','-') = sqr(m1results('side','-'));

m1results('angle (rad)','-') = theta.l;

m1results('angle (deg)','-') = theta.l * 180 / pi;

display m1results;

*---------------------------------------------------------------

* MIQCP model without angles

*---------------------------------------------------------------

* reset variable levels so we don't have a good starting point for the MIQCP model

* just in case the solver would use these

sq.l(a) = 0;

sq2.l(corners,c) = 0;

d.l(i,c,case) = 0;

q.l(i,c) = 0;

variables

sintheta 'sine of angle'

costheta 'cosine of angle'

;

* note: angle is between 0 and 0.5pi, so cos and sin are nonnegative

costheta.lo = 1e-6;

sintheta.lo = 0;

costheta.up = 1;

sintheta.up = 1;

equations

transfdata2x(i) 'transform data points: x coordinate'

transfdata2y(i) 'transform data points: y coordinate'

transform2x(corners) "transform x',y' -> x back to original coordinates"

transform2y(corners) "transform x',y' -> y back to original coordinates"

sincos 'unit circle constraint'

;

* transform data points

transfdata2x(i).. q(i,'x') =e= p(i,'x')*costheta + p(i,'y')*sintheta;

transfdata2y(i).. q(i,'y') =e= -p(i,'x')*sintheta + p(i,'y')*costheta;

* transform back to original coordinates

transform2x(corners).. sq2(corners,'x') =e= (sq('x')+offset(corners,'x')*sq('s'))*costheta - (sq('y')+offset(corners,'y')*sq('s'))*sintheta;

transform2y(corners).. sq2(corners,'y') =e= (sq('x')+offset(corners,'x')*sq('s'))*sintheta + (sq('y')+offset(corners,'y')*sq('s'))*costheta;

sincos.. sqr(costheta) + sqr(sintheta) =e= 1;

model emptysquare2 /all-transfdatax-transfdatay-transformx-transformy/;

solve emptysquare2 maximizing z using miqcp;

scalar angle;

* recover the angle from the sine and cosine values

angle = arctan(sintheta.l/costheta.l);

display angle,sintheta.l, costheta.l, sq.l, sq2.l;

parameter m2results(*,*) 'optimal values model 2';

m2results(corners,c) = sq2.l(corners,c);

m2results('side','-') = sq.l('s');

m2results('area','-') = sqr(m2results('side','-'));

m2results('angle (rad)','-') = angle;

m2results('angle (deg)','-') = angle * 180 / pi;

display m1results,m2results;

parameter perf(*,*) 'performance comparison';

$macro modelperf(m,s) \

perf(s,'#vars') = m.numvar;\

perf(s,'#dvars') = m.numdvar;\

perf(s,'#equs') = m.numequ;\

perf(s,'obj') = m.objval;\

perf(s,'time') = m.resusd;\

perf(s,'its') = m.iterusd;\

perf(s,'nodes') = m.nodusd;

modelperf(emptysquare1,'model 1')

modelperf(emptysquare2,'model 2')

display perf;

*---------------------------------------------------------------

* visualization

*---------------------------------------------------------------

$set html plot.html

$set data data.js

$if %htmlplot%==0 $goto skipplot

file fdata /%data%/; put fdata;

* points

put "p=["/;

loop(i,

put " {i:'",i.tl:0,"',x:",p(i,'x'):0:4,",y:",p(i,'y'):0:4,"},"/;

);

put "]"/;

put "box={x:0,y:0,size:",maxsize:0:4,"}"/;

put "angle=",m1results('angle (rad)','-'):0:4/;

put "angledeg=",m1results('angle (deg)','-'):0:4/;

put "side=",m1results('side','-'):0:4/;

put "area=",m1results('area','-'):0:4/;

put "cp=["/;

loop(corners,

put " {k:'",corners.tl:0,"',x:",m1results(corners,'x'):0:4,",y:",m1results(corners,'y'):0:4,"},"/;

);

put "]"/;

putclose;

$onecho > %html%

<html>

<style>

.brdr {

border: 1px solid black;

border-collapse: collapse;

}

th, td {

padding-left: 4px;

padding-right: 4px;

}

</style>

<h1>Solution in Original Coordinates</h1>

<tr><th colspan=3 class="brdr">Rotated Square</td></tr>

<tr><td class="brdr">Corner 1</td><td class="brdr" id="c00_x"></td><td class="brdr" id="c00_y"></td></tr>

<tr><td class="brdr">Corner 2</td><td class="brdr" id="c01_x"></td><td class="brdr" id="c01_y"></td></tr>

<tr><td class="brdr">Corner 3</td><td class="brdr" id="c10_x"></td><td class="brdr" id="c10_y"></td></tr>

<tr><td class="brdr">Corner 4</td><td class="brdr" id="c11_x"></td><td class="brdr" id="c11_y"></td></tr>

<tr><td class="brdr">Angle</td><td class="brdr" id="sol_angle"></td><td class="brdr"></td></tr>

</table>

</td><td>

</td></tr></table>

px = p.map((x) => x['x'])

py = p.map((x) => x['y'])

document.getElementById("c00_x").innerHTML = cp[0]['x'];

document.getElementById("c00_y").innerHTML = cp[0]['y'];

document.getElementById("c01_x").innerHTML = cp[1]['x'];

document.getElementById("c01_y").innerHTML = cp[1]['y'];

document.getElementById("c10_x").innerHTML = cp[2]['x'];

document.getElementById("c10_y").innerHTML = cp[2]['y'];

document.getElementById("c11_x").innerHTML = cp[3]['x'];

document.getElementById("c11_y").innerHTML = cp[3]['y'];

document.getElementById("sol_angle").innerHTML = angledeg+'°';

document.getElementById("sol_side").innerHTML = side;

document.getElementById("sol_area").innerHTML = area;

var path2 = `M ${cp[0]['x']} ${cp[0]['y']} L ${cp[1]['x']} ${cp[1]['y']} L ${cp[3]['x']} ${cp[3]['y']} L ${cp[2]['x']} ${cp[2]['y']} Z`;

var data1 = {

x: px,

y: py,

mode: 'markers',

type: 'scatter',

color: 'darkblue'

};

var layout4 = {

autosize: false,

width: 650,

height: 600,

showlegend: true,

title: {text:'Solution'},

shapes: [{type:'rect',

xref:'x',

yref:'y',

x0:box['x'],

y0:box['y'],

x1:box['x']+box['size'],

y1:box['y']+box['size'],

line:{color:'black',width:2,opacity:0.5},

fillcolor:'rgba(0, 250, 250, 0.2)',

name:'box',

showlegend:true

{type:'path',

path:path2,

xref:'x',

yref:'y',

line:{color:'darkred',width:2,opacity:0.5},

fillcolor:'rgba(0, 150, 150, 0.5)',

name:'optimal square',

showlegend:true

]

}

var trc1 = [data1];

var options1 = {staticPlot: true, displayModeBar: false, responsive: false};

Plotly.newPlot('plotDiv3', trc1, layout4, options1);

</script>

$offecho

executetool 'win32.ShellExecute "%html%"';

$label skipplot

Conclusions

Finding the largest empty square rotated by a given angle ${\color{darkblue}\theta}$ leads, maybe somewhat surprisingly, to a linear MIP model. This solves very quickly.

Also, surprisingly, when we make ${\color{darkred}\theta}$ endogenous, i.e., let the model find the best angle, the resulting nonconvex MINLP model is also not very difficult. I somewhat expected a range of problems, but this model solves them rather quickly (within seconds) with standard commercial solvers. Gurobi and Xpress seem to do exceptionally well on this problem. We see that off-the-shelf global solvers can tackle non-trivial problems and deliver proven global solutions. Of course, we should always be aware that there are problems where this is not the case [3].

Note that there are algorithmic approaches to this problem [4]. In this post, we focused only on formal mathematical optimization.

References

2 comments:

Andrea NardiMay 18, 2026 at 4:13 AM
I have not tried myself, but what if we do the variable re-write $a = sin(\theta)$ and $b = cos(\theta)$, with the additional constraint of $a^2 + b^2 = 1$? This would remove the trigonometric functions and the variable $\theta$, and I foresee it becoming a (non-convex) quadratic integer linear program, but again, I have not tried myself.

Saturday, May 16, 2026

Largest Empty Rotated Square: lots of trigonometry

Model 1. Exogenous Angle

Model 2. Endogenous Angle

Conclusions

References

2 comments: