How to use correctly SLSQP algoritm with non-linear constraints?

Question

I need to find the rectangle with max area inside an ellipse (which may be tilted). The goal is to gerealize this problem to N dimension, so when we set N=2 we find our rectangle iside an ellipse.

So mathematically speaking I want to maximize the function

with the constraint

I use SLSQP from scipy.optimiize.minimize() to find the optimal , but I dont know how to give the constraint correctly:

here's my code :

you will need these two functions to plot the ellipse and the rectangle

import numpy as np
import matplotlib.pyplot as plt

def ellipse(A, ax=plt):
    eigenvalues, eigenvectors = np.linalg.eig(A)
    theta = np.linspace(0, 2*np.pi, 1000);
    ellipsis = (1/np.sqrt(eigenvalues[None,:]) * eigenvectors) @ [np.sin(theta), np.cos(theta)]
    ax.plot(ellipsis[0,:], ellipsis[1,:])
    
def rectangle(x, y, ax=plt):
    ax.plot([-x,-x,x,x,-x],
            [-y,y,y,-y,-y])

import scipy
from scipy.optimize import NonlinearConstraint

def f(X):
    return -2**len(X) * np.prod(X) # we minimize -f

def g(X):
    A =np.array([[1, -0.7],
                 [-0.7, 4]])
    return X.T @ A @ X - 1

contraintes = {'type' : 'ineq', 'fun' : g} # automatically g >= 0 
res = scipy.optimize.minimize(f, x0, method="SLSQP", constraints=contraintes)
x, y = res.x

#================== or ===================
# contraintes = NonlinearConstraint(g, 0, np.inf)
# res = scipy.optimize.minimize(f, x0, method="SLSQP", constraints=contraintes)
# x, y = res.x

A =np.array([[1, -0.7],
             [-0.7, 4]])

fig, ax = plt.subplots(1,1, figsize=(5,5))
ellipse(A,ax=ax)
rectangle(x, y, ax=ax)

here's the output:

as you can see the rectangle is not inside. How do we give the constraints correctly to optimize function? I checked the scipy page and it didn't help me to understand.

Answer 1

First of all, please note that scipy.optimize.minimize expects the inequality constraints to be non-negative, i.e. g(x) >= 0, so your code comment # automatically g >= 0 is wrong. You need to explicitly transform your constraint by multiplying it with -1.0. And speaking of your constraints, please also note that the constraint x.T @ A @ x - 1 <= 0 only ensures that one vertex (x1,x2) and its mirror point (-x1,-x2) (due to symmetry) of your rectangle lies inside the ellipsis but not the whole rectangle.

For this end, you need to impose a additional constraint such that (-x1, x2) and (x1, -x2) lie within the ellipsis too:

import numpy as np
from scipy.optimize import minimize

def f(x):
    return -2**len(x)*np.prod(x)

def g(x, A):
    xx = [-1.0, 1.0] * x
    ineq1 = x.T @ A @ x - 1
    ineq2 = xx.T @ A @ xx - 1
    return np.array((ineq1, ineq2))

A = np.array([[1, -0.7], [-0.7, 4]])

# constraints
cons = [{'type': 'ineq', 'fun': lambda x: -1.0*g(x, A)}]

# optimize
res = minimize(f, x0=0.1*np.ones(2), method="SLSQP", constraints=cons)