神秘的资源分配问题产生次优解决方案

Question

考虑下面这个顺序资源分配问题的符号和目标：

花费/成本时间步长 i 通道 j：C_{i, j}

总资源：B

地平线：H

时间步之间的最大分配跳跃：1 / 0.4

最小分配：min_total_spend_channel_i

目标投资回报率：目标值

计划内的支出：P = (C_{1, 1}, C_{1, 2}.....C_{H, 1}, C_{H, 2})

我们有一个特定日期销售额的时变预测器。对于某个时间步，它采用以下形式：

sales_{i, j} = coeff_{i, j} * spend_{i, j}^{saturation_{{i, j}}}

其中 i 表示时间步长，j 表示通道。

以下成立：0 <= saturation_{i ,j} <= 1 , notice that this curve exhibits diminishing marginal returns and is concave when minimizing thus convex when maximizing.

地平线上的总销售额可以表示为：

totalsales(P) = \sum_{i}\sum_{j} coeff_{i, j} * spend_{i, j}^{saturation_{{i, j}}}

此优化任务的目的是在保持目标 ROI（投资回报率）的同时最大化支出。我们也有一些我们希望我们的解决方案遵守的约束，例如每个渠道的最低支出、总预算约束、一天与另一天之间的预算之间的最大跳转。

我们构建以下优化问题：

这是一个玩具示例，而不是我正在研究的真实示例。

但是，在谈到这一点之前，我想让它开始工作。 由于我的实际问题更复杂，尤其是销售预测器，我一直在尝试寻找可以解决受约束的非线性非凸优化问题的 python 优化器。

浏览完这个论坛后，我偶然发现了 mystic 并尝试了一下：

代码如下：

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)

结果远非最佳，惩罚被严重违反。我浏览了文档，但找不到明显的解释。

请注意，存在一个解决方案，打印出的结果 penaltyvalue 应该接近于零，并且整个预算 B 都花光了。事实并非如此，消除惩罚使得至少可以将预算花费到任意常数，但一旦惩罚开始发挥作用，我们就离解决方案很远了。

Answer 1

在（6）的两边都乘以分母后，这个问题在我看来是线性的。线性求解器会更合适。

警告：我不知道

totalSales(P)

是什么。我认为它可能是外生的（因为它只取决于外生的

P

），但这是一个猜测。