Modeling and Solving Non-Convex Quadratic Problems using Hexaly

Tutorial

Frédéric Senior Optimization Scientist July 11, 2022

A recent post on Erwin Kalvelagen’s blog, Yet Another Math Programming Consultant, shows results obtained by state-of-the-art solvers on some non-convex quadratic problems. Hexaly can compute lower bounds and prove optimality using global, exact optimization techniques. This article presents the results of Hexaly Optimizer with default parameter setting on the non-convex quadratic models presented in Erwin Kalvelagen’s blog post. In particular, Hexaly is easier to use and faster than Cplex and Gurobi.

Largest Empty Rectangle

We try to find the largest empty rectangle in a space filled with random points. To model this problem with Hexaly, you only need 4 decision variables: the coordinates of two opposite corners of the rectangle. You then ensure that each point is either left, right, below, or above the rectangle. Below is the Python code to solve this problem. You can notice the easy formulation of the constraint thanks to the or operator available among the Hexaly mathematical modeling features.

import hexaly
import random

if __name__ == '__main__':
    random.seed(2)
    n_points = 50
    points = []
    
    for i in range(n_points):
        points.append((random.uniform(0, 1), random.uniform(0, 1)))

    with hexaly.Hexaly() as hx:
        model = hx.model
        x1, y1 = model.float(0,1), model.float(0,1)
        x2, y2 = model.float(0,1), model.float(0,1)
        
        model.constraint(x1 <= x2)
        model.constraint(y1 <= y2)
        
        for x,y in points:
            model.constraint(model.or_(x1 >= x, x2 <= x, y1 >= y, y2 <= y))
       
        model.maximize((x2 - x1) * (y2 - y1))
                
        model.close()
        hx.solve()

For 50 random points, this problem shows no difficulty for Hexaly and is solved optimally in less than 1 second. On the other hand, such models require specific parameters to be solved by Cplex and Gurobi. Moreover, the resulting MIQP models are much harder to write using big-M or indicator constraints; please look at Erwin Kalvelagen’s blog post for details.

Adding new points does not increase the number of decisions in the model but increases the number of constraints. The same problem with 1000 random points instead of 50 is proven optimal within a minute. For larger-scale instances, high-quality feasible solutions are provided in minutes, but high-quality bounds are harder to prove.

Largest Empty Box

If we move on to the 3D version of the same problem, the formulation in Hexaly remains almost the same.

if __name__ == '__main__':
    random.seed(7)
    n_points = 50
    points = []
    
    for i in range(n_points):
        points.append((random.uniform(0, 1), random.uniform(0, 1), random.uniform(0, 1)))

    with hexaly.Hexaly() as ls:
        model = hx.model
        x1, y1, z1 = model.float(0,1), model.float(0,1), model.float(0,1)
        x2, y2, z2 = model.float(0,1), model.float(0,1), model.float(0,1)
        
        model.constraint(x1 <= x2)
        model.constraint(y1 <= y2)
        model.constraint(z1 <= z2)
        
        for x,y,z in points:
            model.constraint(model.or_(x1 >= x, x2 <= x, y1 >= y, y2 <= y, z1 >= z, z2 <= z))
       
        model.maximize((x2 - x1) * (y2 - y1) * (z2 - z1))
                
        model.close()
        hx.solve()

Despite the objective becoming cubic instead of quadratic, Hexaly can still handle this problem. In this example, with 50 random points, optimality is proved in 1 second. Erwin Kalvelagen’s blog post explained that Cplex could not solve this problem because the MIQP reformulation induces non-convex constraints.

In the same way as for the rectangle, one can increase the number of random points in the box. Within 1 minute, we obtain the optimality for instances of 350 points and a gap below 5% for instances up to 500 points. For larger instances, with 1000 points, Hexaly can still find good solutions, but bounds are harder to obtain.

Hostile Brothers

The second non-convex quadratic model considered on the blog is the hostile brothers’ problem. In this case, we need to find a way to decide the position of each brother’s house on a plot of land, such that the minimum distance between two brothers is maximized. Using Hexaly, you can write the objective straightforwardly using the min operator.

import hexaly

if __name__ == '__main__':
    with hexaly.Hexaly() as ls:
        model = hx.model
        points = []
        N = 50
        for i in range(N):
            points.append((model.float(0, 1), model.float(0, 1)))
        
        dists = []
        for i in range(N):
            for j in range(i+1, N):
                dx = points[i][0] - points[j][0]
                dy = points[i][1] - points[j][1]
                dists.append(dx * dx + dy * dy)
        
        model.maximize(model.min(dists))
        model.close()
        
        hx.param.time_limit=60
        hx.solve()

Within a time limit of 1 minute, Hexaly obtains a solution with an objective value of 0.02423.

You can find the visual of the solution below. The pattern is not perfectly regular, but the distribution is already quite good. As stated in Erwin Kalvelagen’s post, Gurobi cannot find a quality, feasible solution despite 2 hours of running times, while Cplex cannot solve this problem again.

For a larger number of points, Hexaly can still come up with solutions that keep improving with time. The two figures below are the results for the hostile brothers’ problem with 500 points, obtained with Hexaly after 1 minute and 1 hour.

Hexaly is the solver of choice for tackling non-convex quadratic problems. Are you interested in trying it out? Get free trial licenses here. For more information, don’t hesitate to get in touch with us. We will gladly exchange with you to understand your problems and needs in mathematical optimization.