GEKKO returned non-optimal solution

user_16802832 0
Viewed 62

I want to use GEKKO to solve the following optimization problem:

Minimize x'Qx + 1e-10 * sum_{i=1}^n x_i^0.1
subject to 1' x = 1 and x >= 0

However, the following code returns sol = [0., 0., 0. ,0. ,1.] and Objective: 1.99419 as a solution. Which is far from optimal, I'll explain why below.

import numpy as np
from gekko import GEKKO
n = 5
m = GEKKO(remote=False)
m.options.SOLVER = 1
m.options.IMODE = 3

x = [m.Var(lb=0, ub=1) for _ in range(n)]

m.Equation(m.sum(x) == 1)
np.random.seed(0)
Q = np.random.uniform(-1, 1, size=(n, n))
Q = np.dot(Q.T, Q)

## Add h_i^p
c, p = 1e-10, 0.1
for i in range(n):
    m.Obj(c * x[i] ** p)
    for j in range(n):
        m.Obj(x[i] * Q[i, j] * x[j])

m.solve(disp=True)
sol = np.array(x).flatten()

This is clearly wrong since if we only optimize the quadratic part (x'Qx) using below code, and put the solution to the initial objective, we get a much smaller objective value (Objective: 0.02489503). The 1e-10 * sum_{i=1}^n x_i^p is esentially ignored since it is very small.

m1 = GEKKO(remote=False)
m1.options.SOLVER = 1
m1.options.OTOL = 1e-10
x1 = [m1.Var(lb=0, ub=1) for _ in range(n)]
m1.Equation(m1.sum(x1) == 1)
m1.qobj(b=np.zeros(n), A=2 * Q, x=x1, otype='min')
m1.solve(disp=True)
sol = np.array(x1).flatten()

Is there any way to resolve this? Thank you!

gekko mathematical-optimization python
2 Answers

Gekko solves nonlinear programming optimization problems with gradient-based methods: interior point and active set SQP. It looks like there is a problem with the objective function. Use matrix operations in Numpy to simplify the objective definition.

## Create Objective
c, p = 1e-10, 0.1
obj = np.dot(np.dot(x,Q),x) + c*m.sum([xi**p for xi in x])
m.Minimize(obj)

Here is the modified script that solves with Gekko. Increase MAX_ITER if the default limit of 250 is reached.

import numpy as np
from gekko import GEKKO
n = 5
m = GEKKO(remote=False)
m.options.SOLVER = 3
m.options.IMODE = 3

x = m.Array(m.Var,n,value=0.1, lb=1e-6, ub=1)

m.Equation(m.sum(x) == 1)

np.random.seed(0)
Q = np.random.uniform(-1, 1, size=(n, n))
Q = np.dot(Q.T, Q)
print(Q)

## Create Objective
c, p = 1e-10, 0.1
obj = np.dot(np.dot(x,Q),x) + c*m.sum([xi**p for xi in x])
m.Minimize(obj)

# adjust solver tolerance
m.options.RTOL=1e-10
m.options.OTOL=1e-10
m.options.MAX_ITER = 1000

m.solve(disp=True)
sol = np.array(x).flatten()
print('x: ', sol)
print('obj: ', m.options.OBJFCNVAL)

This gives an optimal solution that is also global because it is a Quadratic Programming (QP) problem (convex optimization). Using a nonlinear programming (SQP) solver for QP problems gives a solution with the IPOPT solver:

x:  [[0.36315827507] [0.081993130341] [1e-06] [0.086231281612] [0.46861632269]]
obj:  0.024895918696
user_23669410

As far as I could see, gekko looks like it's built for machine learning, which focuses on local optimization opposed to global optimization, and typically most libraries will not be able to guarantee you optimal solutions.

If you really want optimal solutions, than for this case I would suggest looking into interval arithmetic. There are packages such as mpmath which can offer this, though I have yet to see optimizers using it in my brief time searching.

The TL;DR on how interval arithmetic works is you feed in a range of inputs and get back a range of outputs. For example, you can test if 1 is in the range of possible outputs for x1 + x2 + x3 + x4, and you can see the minimum/maximum potential values for your objective function. In this way, you can progressively split your intervals in half, keeping only intervals for which your constraints are potentially satisfied and for which your objective function's maximum potential is at least the largest minimum potential. This allows you to achieve guaranteed convergence to global optimums at the cost of a lot more computation.

user_77415910
Related