Consider the notation and objective below for this sequential resource allocation problem:
Spend/Cost timestep i channel j: C_{i, j}
Total resource: B
Horizon: H
Max allocation jump between time-steps: 1 / 0.4
Minimum allocation : min_total_spend_channel_i
target ROI: targetvalue
planned spend over horizon: P = (C_{1, 1}, C_{1, 2}.....C_{H, 1}, C_{H, 2})
We have a time-variant predictor for the amount of sales on a specific day. For a certain timestep it takes the following form:
sales_{i, j} = coeff_{i, j} * spend_{i, j}^{saturation_{{i, j}}}
where as i denotes the timestep and j denotes the channel.
The following holds: 0 <= saturation_{i ,j} <= 1 , notice that this curve exhibits diminishing marginal returns and is concave when minimizing thus convex when maximizing.
The total sales over the horizon can then be expressed as:
totalsales(P) = \sum_{i}\sum_{j} coeff_{i, j} * spend_{i, j}^{saturation_{{i, j}}}
The aim of this optimization task is to maximize the spend whilst maintaining a target ROI(return on investment). We also have some constraints we want our solution to adhere to such as minimum spend per channel, total budget constraint, max jump between the budgeting between one day and another.
We construct the following optimization problem:
This is a toy example and not the real life example i am working on.
However, before touching on that i want to get this working. Since my real problem is more complex and especially the sales-predictor i been trying to find a python optimizer that can solve constrained nonlinear nonconvex optimization problems.
After browsing through this forum i stumbled upon mystic and gave it a go:
Here is the code:
import numpy as np
from typing import List
import mystic
import pandas as pd
from mystic.symbolic import simplify
from mystic.symbolic import generate_constraint, generate_solvers
from mystic.monitors import VerboseMonitor
from mystic.penalty import quadratic_inequality
from mystic.penalty import linear_inequality
def adstock_geometric(x: List, theta: float):
x_decayed = np.zeros_like(x)
x_decayed[0] = x[0]
for xi in range(1, len(x_decayed)):
x_decayed[xi] = x[xi] + theta * x_decayed[xi - 1]
return x_decayed
def getobjective2_(x):
adstocks = {'adstock-channel_1': np.array(0.08223851),
'adstock-channel_2': np.array(0.28487006)}
channels = ['channel_1', 'channel_2']
horizon = 12
coeffs = pd.DataFrame({'channel_2-saturationcoefficients': [0.73651 for _ in range(horizon)],
'timevar-channel_2': [0.27527, 0.14287, 0.12847, 0.02112,
0.03544, 0.05230, 0.21811, 0.27527,
0.14287, 0.12847, 0.02112, 0.03544],
'channel_1-saturationcoefficients': [0.70259 for _ in range(horizon)],
'timevar-channel_1': [0.13827, 0.24760, 0.28849,
0.33613, 0.34027, 0.28988,
0.21095, 0.13827, 0.24760,
0.28849, 0.33613, 0.34027]})
incomingcostchannel_1 = np.array([0.10303531, 0.07283391, 0.14465764, 0.12315725, 0.07384496, 0.07709148,
0.06425854, 0.09207129, 0.04114219, 0.04066119, 0.09510273, 0.06785441,
0.08381945, 0.05025105, 0.0617348, 0.08135422, 0.04006317, 0.0348892,
0.05374626, 0.05904413])
incomingcoostchannel_2 = np.array([0.04199663, 0.05350016, 0.05427114, 0.05519786, 0.05449486, 0.05733955,
0.0572036, 0.05479, 0.06165305, 0.06347637, 0.05837555, 0.05841943,
0.05922741, 0.05826971, 0.05652555, 0.06605174, 0.0603985, 0.0617391,
0.05811569, 0.05972648])
def saturate(cost, channel):
return np.power(cost, coeffs[f'{channel}-saturationcoefficients'].values)
tbret = 0
for idx, channel in enumerate(channels):
try:
adstockchannel = adstocks[f'adstock-{channel}']
except KeyError:
adstockchannel = 0.0
costchannel = x[(idx * horizon): (idx * horizon) + horizon]
if channel == "channel_1":
incomingcost = incomingcostchannel_1
else:
incomingcost = incomingcoostchannel_2
zerosoverhorizon = np.zeros(horizon)
incomingcostoverhorizon = np.append(incomingcost, zerosoverhorizon)
incomingcostoverhorizonadstocked = adstock_geometric(x=incomingcostoverhorizon,
theta=adstockchannel)[-horizon:]
costchannel += incomingcostoverhorizonadstocked
costadstocked = adstock_geometric(x=costchannel,
theta=adstockchannel)
costmediaprocessed = saturate(cost=costadstocked, channel=channel)
costready = costmediaprocessed * coeffs[f'timevar-{channel}']
tbret += sum(costready)
return -tbret
def generatemysticconstraints(horizon, totalbudget, channels, minspendtotalconstraints,
maxjumps, maxcovers):
xs = [f"x{idx}" for idx in range(horizon * len(channels))]
constraints = ""
totalbudgetconstraint = f"{totalbudget} >= {'+'.join(xs)}"
constraints += totalbudgetconstraint
breakpoints = {channel: idx * horizon for idx, channel in enumerate(channels)}
for channel, cutoff in breakpoints.items():
fromx = cutoff
tox = cutoff + horizon
minspendchannel = minspendtotalconstraints[channel]
maxjump = maxjumps[channel]
maxcover = maxcovers[channel]
minspendconstrainttotal = f"{'+'.join(xs[fromx: tox])} >= {minspendchannel}"
constraints += '\n'
constraints += minspendconstrainttotal
for idx in range(fromx + 1, tox):
maxjumpconstraint_ = f"({xs[idx - 1]} * {maxjump}) >= {xs[idx]}"
constraints += '\n'
constraints += maxjumpconstraint_
for idx in range(fromx, tox - 1):
maxcoverconstraint_ = f"({xs[idx + 1]} * {maxcover}) >= {xs[idx]}"
constraints += '\n'
constraints += maxcoverconstraint_
return constraints
def run():
scalefactor = 11621.64675797
scalefactorresponse = 118257.501641172
totalbudget = 150_000 / scalefactor
horizon = 12
minspendtotalconstraints = {'channel_2': 0.0,
'channel_1': 0.0}
targetvalue = 4.0
channels = ['channel_1', 'channel_2']
objective = getobjective2_
strategy = "targetroas"
maxjumps = {'channel_1': 2.5,
'channel_2': 1.5}
maxcovers = {'channel_1': 2.5,
'channel_2': 1.5}
x0 = [0.017209265103744407, 0.017209265103744407, 0.017209265103744407, 0.017209265103744407,
0.017209265103744407, 0.017209265103744407, 0.017209265103744407, 0.017209265103744407,
0.017209265103744407, 0.017209265103744407, 0.017209265103744407, 0.017209265103744407,
0.017209265103744407, 0.017209265103744407, 0.017209265103744407, 0.017209265103744407,
0.017209265103744407, 0.017209265103744407, 0.017209265103744407, 0.017209265103744407,
0.017209265103744407, 0.017209265103744407, 0.017209265103744407, 0.017209265103744407]
boundsextrapolation = [(0.014964578351082093, 5.602713613312933), (0.014964578351082093, 5.602713613312933),
(0.014964578351082093, 5.602713613312933), (0.014964578351082093, 5.602713613312933),
(0.014964578351082093, 5.602713613312933), (0.014964578351082093, 5.602713613312933),
(0.014964578351082093, 5.602713613312933), (0.014964578351082093, 5.602713613312933),
(0.014964578351082093, 5.602713613312933), (0.014964578351082093, 5.602713613312933),
(0.014964578351082093, 5.602713613312933), (0.014964578351082093, 5.602713613312933),
(0.014964578351082093, 1.158821015684604), (0.014964578351082093, 1.158821015684604),
(0.014964578351082093, 1.158821015684604), (0.014964578351082093, 1.158821015684604),
(0.014964578351082093, 1.158821015684604), (0.014964578351082093, 1.158821015684604),
(0.014964578351082093, 1.158821015684604), (0.014964578351082093, 1.158821015684604),
(0.014964578351082093, 1.158821015684604),
(0.014964578351082093, 1.158821015684604), (0.014964578351082093, 1.158821015684604),
(0.014964578351082093, 1.158821015684604)]
constraints = generatemysticconstraints(horizon=horizon,
totalbudget=totalbudget,
channels=channels,
minspendtotalconstraints=minspendtotalconstraints,
maxjumps=maxjumps,
maxcovers=maxcovers)
def targetpenalty(x, target):
if strategy == "targetroas":
return target * (sum(x) * scalefactor) - (-objective(x) * scalefactorresponse)
elif strategy == "targetcpa":
return (sum(x) * scalefactor) - target * (-objective(x) * scalefactorresponse)
else:
raise NotImplementedError
@linear_inequality(condition=targetpenalty, kwds={'target': targetvalue})
def penalty(x):
return 0.0
constraintsimplified = simplify(constraints,
all=True)
from mystic.constraints import and_
generatedconstraints = generate_constraint(generate_solvers(constraintsimplified), join=and_)
stepmon = VerboseMonitor(1)
cost = lambda x: -sum(x)
bounds = boundsextrapolation
optimum = mystic.solvers.diffev2(cost,
x0=x0,
constraints=generatedconstraints,
bounds=bounds,
penalty=penalty,
full_output=True,
stepmon=stepmon,
disp=True,
tightrange=True,
npop=400,
gtol=200,
maxiter=10000)
optimalvalueslst = optimum[0]
totalspent = np.sum(optimalvalueslst) * scalefactor
sales = -objective(optimalvalueslst) * scalefactorresponse
estimatedcpa = totalspent / sales
estimatedroas = sales / totalspent
penaltyvalue = targetvalue * sum(optimalvalueslst) - (-objective(optimalvalueslst))
print(f"total spent: {totalspent}")
print(f"estimated cpa: {estimatedcpa}")
print(f"estimated roas: {estimatedroas}")
print(f"penalty value: {penaltyvalue}")
return {'x': optimalvalueslst}
if __name__ == '__main__':
vals = run()
print(vals)
The result is far from optimal, the penalty gets heavily violated. I browsed through the docs but couldnt find an obvious explanation.
Note that an solution exists and that the resulting penaltyvalue printed should be close to zero and the whole budget B spent. Which is not the case, removing the penalty makes it possible to at least spend the budget up to an arbitrary constant, but as soon as the penalty gets into play we get very very far of from the solution.
After multiplying both sides of (6) by the denominator, the problem looks linear to me. A linear solver would be more appropriate.
Caveat: I don't know what
totalSales(P)is. I think it may be exogenous (as it depends only onPwhich is exogenous), but that is a guess.