Suppose, we want to locate some given facilities $\left \{ (i,j) \ |\ (i,j) \in \text[{1,\cdots, N}]\right \}$ in a specific area. Each facility has a predefined dimension with a length $l_{i}$ and width $w_{i}$, (they are not necessarily square). The overall dimensions of the area are $L$ and $W$. The objective is to find the best fit shape of the facilities in the area. As far as I know, this problem can be formulated as mixed-integer programming by defining the non-overlap constraints as follows:
\begin{align} &x_i+w_i \le x_j & \text{or} \\ &x_j+w_j \le x_i & \text{or} \\ &y_i+h_i \le y_j & \text{or} \\ &y_j+h_j \le y_i \end{align}
or other formulations like $\text{ABS}$ models in the facility location literature.
The important thing is that to achieve a good feasible solution by the above-mentioned formulations, we should define the facility shape in the fixed orientation as prior, for example horizontal. What I am interested to know is about the orientation of the facility shape and rotation of them by formulation. I was wondering if, how we can change the above formulation to capture the rotation of the facility shape automatically by using mixed-integer programming?