Constrained non-linear optimisation in GEKKO, objective function not matching expected solution

There was a problem with how the objective function is defined. The fix to the problem is to use a[i], b[i], x[i] one at a time in the functions. Gekko allows more than one objective function definition. When maximizing, it converts to a minimization problem with min = -max. The maximized objective is therefore -m.options.OBJFCNVAL.

# Log-Log model
def roi(ai, bi, xi):
    return m.exp(ai + bi * m.log(xi))

# Objective function
for i,xi in enumerate(x):
    m.Maximize(xi * roi(a[i], b[i], x[i]))

It also helps to define a new variable sum_x with upper spend_max and lower spend_min bounds.

# Constraints
sum_x = m.Var(lb=spend_min,ub=spend_max)
m.Equation(sum_x==m.sum(x))

Here is the complete script.

from gekko import GEKKO
import numpy as np
import pandas as pd

df1 = pd.DataFrame({'channel': ['fb', 'tv', 'ppc'],
                    'value': [5000, 5000, 5000],                    
                    'alpha': [1.00, 1.00, 1.00],
                    'beta': [0.03, 0.02, 0.01]
})

total_spend = df1['value'].sum()

m = GEKKO()

# Constraint parameters
spend_inc = 0.00
spend_dec = 0.00
channel_inc = 0.50
channel_dec = 0.50

channels = len(df1)

# Initialise decision variables
x = m.Array(m.Var, (channels), integer=True)
i = 0
for xi in x:
    xi.value = df1['value'][i]
    xi.lower = df1['value'][i] * (1 - channel_dec)
    xi.upper = df1['value'][i] * (1 + channel_inc)
    i += 1

# Initialise alpha
a = m.Array(m.Param, (channels))
i = 0
for ai in a:
    ai.value = df1['alpha'][i]
    i += 1
    
# Initialise beta
b = m.Array(m.Param, (channels))
i = 0
for bi in b:
    bi.value = df1['beta'][i]
    i += 1
    
# Initial global contraints
spend_min = total_spend * (1 - spend_dec)
spend_max = total_spend * (1 + spend_dec)
    
# Print out variabales
print(f'spend min: {spend_min}')
print(f'spend max: {spend_max}')
print('')
for i in range(0, channels):
    print(f'x{i+1} value: {str(x[i].value)}')    
    print(f'x{i+1} lower: {str(x[i].lower)}')
    print(f'x{i+1} upper: {str(x[i].upper)}')
    print(f'x{i+1} alpha: {str(a[i].value)}')
    print(f'x{i+1} beta: {str(b[i].value)}')
    print('')

# Constraints
sum_x = m.Var(lb=spend_min,ub=spend_max)
m.Equation(sum_x==m.sum(x))

# Log-Log model
def roi(ai, bi, xi):
    return m.exp(ai + bi * m.log(xi))

# Objective function
for i,xi in enumerate(x):
    m.Maximize(xi * roi(a[i], b[i], x[i]))
m.options.IMODE = 3
m.solve()

for i in range(0, channels):
    print(f'x{i+1}: {str(x[i].value[0])}')
print('')
print(f'optimal solution: {str(-m.options.objfcnval)}')

opt = 0
for i in range(0, channels):
    opt += x[i].value[0] * (np.exp((a[i].value[0] \
                                    + b[i].value[0] \
                                    * np.log(x[i].value[0]))))
print(f'optimal solution: {opt}')

# THIS IS THE EXPECTED SOLUTION
# IT NOW MATCHES THE OBJECTIVE FUNCTION VALUE SUGGESTED BY GEKKO
sol = 7500 * np.exp(1.00 + 0.03 * np.log(7500)) \
      + 5000 * np.exp(1.00 + 0.02 * np.log(5000)) \
      + 2500 * np.exp(1.00 + 0.01 * np.log(2500))
print(f'expected optimal solution: {sol}')

The solution is shown below.

EXIT: Optimal Solution Found.
 
 The solution was found.
 
 The final value of the objective function is   -50108.7613900549     
 
 ---------------------------------------------------
 Solver         :  IPOPT (v3.12)
 Solution time  :   9.499999985564500E-003 sec
 Objective      :   -50108.7611889392     
 Successful solution
 ---------------------------------------------------
 
x1: 7500.0
x2: 4999.9999497
x3: 2500.0

optimal solution: 50108.761189
optimal solution: 50108.7611890521
expected optimal solution: 50108.76135441651