I've got a fairly straightforward constrained non-linear optimisation problem, maximising revenue, given some spend constraints (keep overall spend constant, and increase/decrease spend on each channel by upto 50%) and an objective function which is revenue = spend * roi (where roi is calculated using the log-log model coefficients for each channel).
Firstly the solution doesn't match what I believe to be the optimal solution. Also, when take the optimal spend values suggested by GEKKO, and put them into the objective function, this also doesn't match what GEKKO gives as the optimal objective function value.
If the above problem isn't articulated very well, the example code below should bring to life the problem...
Anyone know what I'm missing?
Code below:
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
m.Equation(sum(x) >= spend_min)
m.Equation(sum(x) <= spend_max)
# Log-Log model
def roi(a, b, x):
roi = a + b * m.log(x[0])
roi = m.exp(roi[0])
return roi
# Objective function
m.Maximize(m.sum(x * roi(a, b, x)))
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)}')
# THIS DOESN'T MATCH THE OBJECTIVE FUNCTION VALUE SUGGESTED BY GEKKO
opt = 0
for i in range(0, channels):
opt += x[i].value[0] * (np.exp((a[0].value[0] + b[i].value[0] * np.log(x[i].value[0]))))
print(f'optimal solution: {opt}')
# THIS IS THE EXPECETD SOLUTION, WHICH ALSO DOESN'T MATCH THE OBJECTIVE FUNCTION VALUE SUGGESTED BY GEKKO
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))
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