Линейное программирование с scipy.optimize.linprog - переменные коэффициенты

Question

Попытка оптимизировать с помощью scipy.optimize.linprog функцию стоимости, где коэффициенты стоимости являются функцией переменных; например

Стоимость = c1 * x1 + c2 * x2 (x1,x2 - переменные)

например

если х1 = 1, с1 = 0,5

если х1 = 2, с1 = 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

Значения fi () рассчитываются для каждого возможного значения xi

Я надеюсь, что это проясняет модель.

Answer 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 (of 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)