pythonnonlinear-optimizationgekko

Gekko returning a clearly non-optimal solution for a simple quadratic problem?


Given the program:

from gekko import GEKKO
m = GEKKO()  

tlX = m.Var()
tlY = m.Var()

m.Equation(m.sqrt(tlX*tlX + tlY*tlY)  <= 31 )

m.Maximize(tlX * tlY)

m.solve() 

print("solution: ", tlX.value, tlY.value )

Returns:

EXIT: Optimal Solution Found.
 The final value of the objective function is   0.000000000000000E+000
[0.0] [0.0]

Which is clearly not optimal. As tlX=1, tlY=1 is closer to optimal. What am I missing?


Solution

  • The others that responded are correct that any non-zero initial guess will solve the problem. Even setting a small non-zero value for one of the variables resolves the issue.

    tlX = m.Var(value=1e-5)
    

    The problem is that the local optimizers find that the initial guess of [0,0] satisfies the Karush-Kuhn Tucker conditions. Because they are local optimizers, they don't check that the Hessian (matrix of partial 2nd derivatives) is positive definite for a minimum or negative definite for a maximum.

    KKT conditions

    A couple things to consider for global solutions are listed on this page. Here is an example with your problem.

    from gekko import GEKKO
    from hyperopt import fmin, tpe, hp
    from hyperopt import STATUS_OK, STATUS_FAIL
    
    # Define the search space for the hyperparameters
    space = {'tlX': hp.quniform('tlX', 0, 1, 0.3),
             'tlY': hp.quniform('tlY', 0, 1, 0.3)}
    
    def objective(params):
        m = GEKKO(remote=False)
        x = m.Array(m.Var,2,lb=0)
        tlX,tlY = x
        tlX.value = params['tlX']
        tlY.value = params['tlY']
        m.Equation(m.sqrt(tlX*tlX + tlY*tlY)  <= 31 )
        m.Maximize(tlX * tlY)
        m.solve(disp=False,debug=False)
        obj = m.options.objfcnval
        if m.options.APPSTATUS==1:
            s=STATUS_OK
        else:
            s=STATUS_FAIL
        m.cleanup()
        return {'loss':-obj, 'status': s, 'x':x}
    
    best = fmin(objective, space, algo=tpe.suggest, max_evals=10)
    sol = objective(best)
    print(f"Solution Status: {sol['status']}")
    print(f"Objective: {sol['loss']:.2f}")
    print(f"Solution: {sol['x']}")
    

    Here is the solution:

    Solution Status: ok
    Objective: 480.50
    Solution: [[21.920310221] [21.920310221]]