Modeling and Solving Non-Convex Quadratic Problems

Tutorial
Senior Optimization Scientist

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. LocalSolver can compute lower bounds and prove optimality using global, exact optimization techniques. This article presents the results of LocalSolver Optimizer with default parameter setting on the non-convex quadratic models presented in Erwin Kalvelagen’s blog post. In particular, LocalSolver 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 LocalSolver, 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 LocalSolver mathematical features.

import localsolver
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 localsolver.LocalSolver() as ls:
        model = ls.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()
        ls.solve()

For 50 random points, this problem shows no difficulty for LocalSolver 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 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 LocalSolver 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 localsolver.LocalSolver() as ls:
        model = ls.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()
        ls.solve()

Despite the objective becoming cubic instead of quadratic, LocalSolver 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, LocalSolver 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 LocalSolver, you can write the objective straightforwardly using the min operator.

import localsolver

if __name__ == '__main__':
    with localsolver.LocalSolver() as ls:
        model = ls.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()
        
        ls.param.time_limit=60
        ls.solve()

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

You can find below the visual of the solution. 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, LocalSolver 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 LocalSolver after 1 minute and 1 hour.

LocalSolver 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 be glad to exchange with you to understand your problems and needs in optimization.

Share

Recent posts

Categories