使用两个参数的 Scipy Optimizer 约束

Question

我试图通过找到一个人会使用的最佳N 单位来最大化效用函数。限制之一是他们的资金有限，m。所以我试图设置一个约束，其中长度为 3 的数组，N 乘以价格，P 也是长度为 3 的数组，不能大于m。

如下示例：

P = np.array([3,4,5])
N = np.array([1,2,1])
m = 50
sum(P*N) > m

对于这个优化，P 是基于之前的优化给出的。现在这是我的实际代码：

cons_c = [{'type':'ineq', 'fun': lambda N: 10 - sum(np.round(N)*P)},{'type':'ineq', 'fun': lambda N: 24 - sum(N*T)}]
bnds = [(0.,None) for x in range(len(N))]
optimized_c = scipy.optimize.minimize(utility_c, N, (P,Q,T), method='SLSQP', bounds=bnds, constraints=cons_c)

功能：

def utility_c(N,P,Q,T):

    print "N: {0}".format(N)
    print "P: {0}".format(P)
    print "Q: {0}".format(Q)
    print "T: {0}".format(T)
    N = np.round(N)
    m = 10 - sum(N*P)
    b = sum(N*Q)
    t = 24 - sum(N*T)
    print "m in C: {0}".format(m)
    print "b: {0}".format(b)
    print "t: {0}".format(t)
    # if m < 0 or t < 0:
    #     return 0
    return 1/ ((b**0.3)*(t**0.7))+(5*(m**0.5))

问题是我仍然否定m！很明显，我没有正确通过约束。我猜这是因为P 没有正确使用？

输出：

N: [ 1.  1.  1.]
P: [  5.  14.   4.]
Q: [ 1.  3.  1.]
T: [ 1.    1.    1.01]
m in C: -13.0

我的尝试：

我也尝试过在 args 中传递 P，如下所示：

cons_c = [{'type':'ineq', 'fun': lambda N,P: 10 - sum(np.round(N)*P), 'args':P},{'type':'ineq', 'fun': lambda N: 24 - sum(N*T)}]

但它告诉我`Lambda 想要 2 个参数并收到 4 个

** 更新：**

在'args' 中使用(F,) 不允许程序在不引发错误的情况下运行，但约束仍然无法成立。

另外，nan 在 m 定义为负值之后返回，这当然会使整个 scipy 优化变得异常。

** 完整项目代码：**

import scipy.optimize
import numpy as np
import sys

def solve_utility(P,Q,T):
    """
    Here we are given the pricing already (P,Q,T), but solve for the quantities each type
    would purchase in order to maximize their utility (N).
    """

    def utility_a(N,P,Q,T):
        N = np.round(N)
        m = 50 - sum(N*P)
        b = sum(N*Q)
        t = 8 - sum(N*T)
        return 1/ ((b**0.5)*(t**0.5))+(5*(m**0.5))

    def utility_b(N,P,Q,T):
        N = np.round(N)
        m = 50 - sum(N*P)
        b = sum(N*Q)
        t = 8 - sum(N*T)
        return 1/ ((b**0.7)*(t**0.3))+(5*(m**0.5))

    def utility_c(N,P,Q,T):
        N = np.round(N)
        print "N: {0}".format(N)
        print "P: {0}".format(P)
        print "Q: {0}".format(Q)
        print "T: {0}".format(T)
        m = 10 - sum(N*P)
        b = sum(N*Q)
        t = 24 - sum(N*T)
        print "m in C: {0}".format(m)
        print "b: {0}".format(b)
        print "t: {0}".format(t)
        return 1/ ((b**0.3)*(t**0.7))+(5*(m**0.5))



    # Establishing constraints so no negative money or time:
    N = np.array([2,2,1])
    cons_a = [{'type':'ineq', 'fun': lambda N, P: 50 - sum(np.round(N)*P), 'args':(P,)},{'type':'ineq', 'fun': lambda N: 8 - sum(N*T)}]
    cons_b = [{'type':'ineq', 'fun': lambda N, P: 50 - sum(np.round(N)*P), 'args':(P,)},{'type':'ineq', 'fun': lambda N: 8 - sum(N*T)}]
    cons_c = [{'type':'ineq', 'fun': lambda N, P: 10 - sum(np.round(N)*P), 'args':(P,)},{'type':'ineq', 'fun': lambda N: 24 - sum(N*T)}]

    maxes = P/50
    bnds = [(0.,None) for x in range(len(N))]
    b = [()]

    optimized_a = scipy.optimize.minimize(utility_a, N, (P,Q,T), method='SLSQP', constraints=cons_a)
    optimized_b = scipy.optimize.minimize(utility_b, N, (P,Q,T), method='SLSQP', constraints=cons_b)
    optimized_c = scipy.optimize.minimize(utility_c, N, (P,Q,T), method='SLSQP', constraints=cons_c)

    if not optimized_a.success:
        print "Solving Utilities A didn't work..."
        return None
    if not optimized_b.success:
        print "Solving Utilities B didn't work..."
        return None
    if not optimized_c.success:
        print "Solving Utilities C didn't work..."
        return None
    else:
        print "returning N: {0}".format(np.array([optimized_a.x,optimized_b.x,optimized_c.x]))
        return np.array([optimized_a.x,optimized_b.x,optimized_c.x])


# solve_utility(P,Q,T,N)


def solve_profits():
    """ 
    Here we build the best pricing strategy to maximize solve_profits
    """

    P = np.array([   3,      10.67,       2.30]) # Pricing
    Q = np.array([   1,       4,       1]) # Quantity of beer for each unit
    T = np.array([   1,       1,    4]) # Time cost per unit
    N = np.array([   1,       0,       1]) # Quantities of unit taken by customer

    def profit(X):
        P,Q,T = X[0:3], X[3:6], X[6:9]
        Q[1] = round(Q[1]) # needs to be an integer
        N = solve_utility(P,Q,T)
        print "N: {0}".format(N)
        N = np.sum(N,axis=1)
        # print "P: {0}".format(P)
        # print "Q: {0}".format(Q)
        # print "T: {0}".format(T)
        denom = sum(N*P*Q) - sum(Q*N)
        return 1/ (sum(N*P*Q) - sum(Q*N))

    cons = [{'type':'ineq', 'fun': lambda X: X[8] - X[6] - 0.01 }, # The time expense for a coupon must be 0.01 greater than regular
            {'type':'ineq', 'fun': lambda X: X[4] - 2 }, # Packs must contain at least 2 beers
            {'type':'eq',   'fun': lambda X: X[3] - 1}, # Quantity has to be 1 for single beer
            {'type':'eq',   'fun': lambda X: X[5] - 1}, # same with coupons
            {'type':'ineq', 'fun': lambda X: X[6] - 1}, # Time cost must be at least 1
            {'type':'ineq', 'fun': lambda X: X[7] - 1}, 
            {'type':'ineq', 'fun': lambda X: X[8] - 1},
            ] 
    X = np.concatenate([P,Q,T])
    optimized = scipy.optimize.minimize(profit, X, method='L-BFGS-B', constraints=cons)

    if not optimized.success:
        print "Solving Profits didn't work..."
    else:
        return optimized.x, N

X, N = solve_profits()
print "X: {0} N {1}".format(X,N)
P,Q,T = X[0:3], X[3:6], X[6:9]
rev = sum(P * Q * N)
cost = sum(Q * N)
profit = (rev-cost)*50
print "N: {0}".format(N)
print "P: {0}".format(P)
print "Q: {0}".format(Q)
print "T: {0}".format(T)
print "profit = {0}".format(profit)

Answer 1

如果您隔离 for optimize_a 并运行，您会看到它抛出的错误是错误 8 - 这是正导数错误。

BFGS 和 SLSQP 都是梯度搜索方法，这意味着它们会根据您的初始猜测，评估梯度及其导数，并寻找最佳的踏入方向，始终下坡并在值变化时停止低于您设置的容差或达到最小值。

该错误表明（至少在您最初的猜测中），该问题没有强导数。一般来说，SQLSP 最适用于可以表述为平方和的问题。也许尝试更现实的初始猜测会有所帮助。我肯定会丢弃大部分代码并首先使用 optimize_a 运行一个最小的示例，一旦你开始工作，其余的就可以跟随了。

也许基于非梯度的求解器会起作用，或者根据问题大小和您对参数的实际限制，全局优化可能是可行的。

如果你没有很好的衍生品可以遵循，scipy 优化就不是很好