-1

尝试使用 scipy.optimize.linprog 优化成本函数，其中成本系数是变量的函数；例如

Cost = c1 * x1 + c2 * x2 # (x1,x2 are the variables)

例如

if x1 = 1, c1 = 0.5

if x1 = 2, c1 = 1.25

等等

*只是为了澄清*

我们正在寻找最小的变量成本；xi; i=1,2,3,... xi 是正整数。

然而，每个 xi 的成本系数是 xi 值的函数。成本是x1*f1(x1) + x2*f2(x2) + ... + c0

fi - 是一个“比率”表；例如- f1(0) = 0; f1(1) = 2.00; f1(2) = 3.00，等等。

xi 受约束，不能为负数，不能超过 qi =>

0 <= xi <= qi

为每个可能的 xi 值计算 fi() 值

我希望它能澄清模型。

4

1 回答 1

1

这是一些原型代码，向您展示如何解决您的问题（关于公式和性能；前者在代码中可见）。

该实现使用 cvxpy 进行建模（仅限凸编程）并基于混合整数方法。

代码

    import numpy as np
    from cvxpy import *

    """
    x0 == 0 -> f(x) = 0
    x0 == 1 -> f(x) = 1
    ...
    x1 == 0 -> f(x) = 1
    x1 == 1 -> f(x) = 4
    ...
    """
    rate_table = np.array([[0, 1, 3, 5], [1, 4, 5, 6], [1.3, 1.7, 2.25, 3.0]])
    bounds_x = (0, 3)  # inclusive; bounds are needed for linearization!

    # Vars
    # ----
    n_vars = len(rate_table)
    n_values_per_var = [len(x) for x in rate_table]

    I = Bool(n_vars, n_values_per_var[0])  # simplified assumption: rate-table sizes equal
    X = Int(n_vars)
    X_ = Variable(n_vars, n_values_per_var[0])  # X_ = mul_elemwise(I*X) broadcasted

    # Constraints
    # -----------
    constraints = []

    # X is bounded
    constraints.append(X >= bounds_x[0])
    constraints.append(X <= bounds_x[1])

    # only one value in rate-table active (often formulated with SOS-type-1 constraints)
    for i in range(n_vars):
        constraints.append(sum_entries(I[i, :]) <= 1)

    # linearization of product of BIN * INT (INT needs to be bounded!)
    # based on Erwin's answer here:
    # https://www.or-exchange.org/questions/10775/how-to-linearize-product-of-binary-integer-and-integer-variables
    for i in range(n_values_per_var[0]):
        constraints.append(bounds_x[0] * I[:, i] <= X_[:, i])
        constraints.append(X_[:, i] <= bounds_x[1] * I[:, i])
        constraints.append(X - bounds_x[1]*(1-I[:, i]) <= X_[:, i])
        constraints.append(X_[:, i] <= X - bounds_x[0]*(1-I[:, i]))

    # Fix chosings -> if table-entry x used -> integer needs to be x
    # assumptions:
    # - table defined for each int
    help_vec = np.arange(n_values_per_var[0])
    constraints.append(I * help_vec == X)

    # ONLY FOR DEBUGGING -> make simple max each X solution infeasible
    constraints.append(sum_entries(mul_elemwise([1, 3, 2], square(X))) <= 15)

    # Objective
    # ---------
    objective = Maximize(sum_entries(mul_elemwise(rate_table, X_)))

    # Problem & Solve
    # ---------------
    problem = Problem(objective, constraints)
    problem.solve()  # choose other solver if needed, e.g. commercial ones like Gurobi, Cplex
    print('Max-objective: ', problem.value)
    print('X:\n' + str(X.value))

输出

('Max-objective: ', 20.70000000000001)
X:
[[ 3.]
 [ 1.]
 [ 1.]]

主意

  • 转变目标max: x0*f(x0) + x1*f(x1) + ...
    • 进入：x0*f(x0==0) + x0*f(x0==1) + ... + x1*f(x1==0) + x1*f(x1==1)+ ...
  • 引入二元变量来制定：
    • f(x0==0) as I[0,0]*table[0,0]
    • f(x1==2) as I[1,2]*table[0,2]
  • 添加约束以限制上述内容I为每个变量仅具有一个非零条目x_i（只有一个扩展的目标组件将处于活动状态）
  • 线性化乘积x0*f(x0==0) == x0*I[0,0]*table(0,0)（整数 * 二进制 * 常数）
  • 修复表查找：使用索引为 x（x0）的表条目应该会导致x0 == x
    • I * help_vec == X)假设表格中没有空白，这可以表述为help_vec == vector(lower_bound, ..., upper_bound)

cvxpy自动（通过构造）证明，我们的公式是凸的，这是大多数求解器所需要的（并且通常不容易识别）。

只是为了好玩：更大的问题和商业解决者

输入生成：

def gen_random_growing_table(size):
    return np.cumsum(np.random.randint(1, 10, size))
SIZE = 100
VARS = 100
rate_table = np.array([gen_random_growing_table(SIZE) for v in range(VARS)])
bounds_x = (0, SIZE-1)  # inclusive; bounds are needed for linearization!
...
...
constraints.append(sum_entries(square(X)) <= 150)

输出：

Explored 19484 nodes (182729 simplex iterations) in 129.83 seconds
Thread count was 4 (of 4 available processors)

Optimal solution found (tolerance 1.00e-04)
Warning: max constraint violation (1.5231e-05) exceeds tolerance
Best objective -1.594000000000e+03, best bound -1.594000000000e+03, gap 0.0%
('Max-objective: ', 1594.0000000000005)
于 2016-08-22T22:44:12.667 回答