Как решить задачу бинарной оптимизации линейным программированием в ortools/cpmodel?

Я пытаюсь оптимизировать бинарную проблему для моего веб-сайта.

Данные содержат примерно 75 предметов, и каждый из них имеет вес (от 50 до 1000) и цену. Вот фрагмент данных:

{"weighting":{
      "0":500,
      "1":50,
      "2":50,
      "3":50,
      "4":250,
      "5":1000
   },
   "price":{
      "0":4,
      "1":78,
      "2":75,
      "3":170,
      "4":5,
      "5":4
   }
}

Я рассчитываю ожидаемое значение всего набора данных с помощью

exp_val = (w1 p1 + w2 p2 + ... + wn pn) / сумма (w1 + w2 + ... wn)

с

сумма (w1 + w2 + ... wn) = 23665 (учитывая все элементы)

Пока все хорошо, но теперь начинается сложная часть. Не все предметы желательны, то есть они стоят меньше и/или имеют большой вес, который разбавляет пул, из которого я могу черпать.

«Блокируя» или удаляя до 3 элементов, я могу извлечь только из оставшихся элементов и тем самым максимизировать функцию ускоренного значения. Вопрос: какие элементы удалить? Поскольку цены со временем меняются, мне приходится регулярно проверять, что нужно удалить.

Я начал с того, что просто удалил товары с наибольшим весом и наименьшей ценой, но я уверен, что это представляет собой только локальный оптимум, и будет более оптимальная стратегия.

После проверки некоторых веб-сайтов кажется, что смешанно-целочисленное линейное программирование (MILP) или, в частности, BILP (двоичное...) может решить мою проблему, и я нашел https://docs.scipy.org/doc/scipy/ reference/generated/scipy.optimize.milp.html, но не смог заставить его работать, так как я застрял, переводя свою проблему в код. Может кто-нибудь мне помочь?

Прежде чем писать код, вы должны правильно сформулировать задачу оптимизации на бумаге. Для начала: введите двоичную переменную x_i, которая равна 1, если элемент i выбран, и 0 в противном случае. Затем вы хотите максимизировать сумму (y * x_i * w_i * p_i) с учетом ограничения 1/sum (x_i * w_i) = y. Обратите внимание, что второе ограничение эквивалентно 1 = sum(y * x_i * w_i), и вы можете линеаризовать произведения y * x_i, см. этот пост для более подробной информации.

joni 07.01.2023 10:41
Стоит ли изучать PHP в 2023-2024 годах?
Стоит ли изучать PHP в 2023-2024 годах?
Привет всем, сегодня я хочу высказать свои соображения по поводу вопроса, который я уже много раз получал в своем сообществе: "Стоит ли изучать PHP в...
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
Поведение ключевого слова "this" в стрелочной функции в сравнении с нормальной функцией
В JavaScript одним из самых запутанных понятий является поведение ключевого слова "this" в стрелочной и обычной функциях.
Приемы CSS-макетирования - floats и Flexbox
Приемы CSS-макетирования - floats и Flexbox
Здравствуйте, друзья-студенты! Готовы совершенствовать свои навыки веб-дизайна? Сегодня в нашем путешествии мы рассмотрим приемы CSS-верстки - в...
Тестирование функциональных ngrx-эффектов в Angular 16 с помощью Jest
В системе управления состояниями ngrx, совместимой с Angular 16, появились функциональные эффекты. Это здорово и делает код определенно легче для...
Концепция локализации и ее применение в приложениях React ⚡️
Концепция локализации и ее применение в приложениях React ⚡️
Локализация - это процесс адаптации приложения к различным языкам и культурным требованиям. Это позволяет пользователям получить опыт, соответствующий...
Пользовательский скаляр GraphQL
Пользовательский скаляр GraphQL
Листовые узлы системы типов GraphQL называются скалярами. Достигнув скалярного типа, невозможно спуститься дальше по иерархии типов. Скалярный тип...
1
1
94
3
Перейти к ответу Данный вопрос помечен как решенный

Ответы 3

Используя комментарии @joni, я перевел их в код. Умножение переменных решения не является линейным, и вам придется линеаризовать продукты, добавив еще одну промежуточную непрерывную переменную решения.

Однако, если вы используете решатель Google-ortool CP-SAT (который я использовал), он может обрабатывать несколько нелинейных операций, таких как умножение, деление и т. д. переменных решения, которые чисто линейный решатель не поддерживает.

Список кодов

from ortools.sat.python import cp_model as cp

model = cp.CpModel()

data = {"weighting":{
      "0":500,
      "1":50,
      "2":50,
      "3":50,
      "4":250,
      "5":1000
   },
   "price":{
      "0":4,
      "1":78,
      "2":75,
      "3":170,
      "4":5,
      "5":4
   }
}

num_items = len(data["weighting"])

dv_select_items = {i : model.NewBoolVar("item_" + i) for i in data["weighting"]}

# constraint : keep only 3 items
model.Add(sum(dv_select_items[i] for i in dv_select_items) == 3)


y = model.NewIntVar(0, 1000000, "")

x_i_w_i_p_i = model.NewIntVar(0, 1000000, "") # x_i * w_i * p_i
model.Add(x_i_w_i_p_i == sum(dv_select_items[i] * data["weighting"][i] * data["price"][i] for i in dv_select_items))

y_x_i_w_i_p_i = model.NewIntVar(0, 1000000000000, "") # y * x_i * w_i * p_i
model.AddMultiplicationEquality(y_x_i_w_i_p_i, [x_i_w_i_p_i, y])

# 1 = sum(y * x_i * w_i)  # constraint
x_i_w_i = model.NewIntVar(0, 1000000, "")
model.Add(x_i_w_i == sum(dv_select_items[i] * data["weighting"][i] for i in dv_select_items))

y_x_i_w_i = model.NewIntVar(0, 1000000000000, "")
model.AddMultiplicationEquality(y_x_i_w_i, [x_i_w_i, y])

remaining = model.NewIntVar(0, 10000, "")
model.Add(1000000 == y_x_i_w_i + remaining)
# 1 = sum(y * x_i * w_i) 
# with above y will be fractional. Now, CP-SAT solver does not support
# fractional values. With using integer values (for y) there is no
# gurantee that 1 = sum(y * x_i * w_i), hence adding "remaining" part
# but we should always try to make "remaining" as close to zero

model.Maximize(y_x_i_w_i_p_i)
solver = cp.CpSolver()
solver.Solve(model)

model.Minimize(remaining)
model.Add(y_x_i_w_i_p_i >= int(solver.ObjectiveValue()))
solver = cp.CpSolver()
solver.Solve(model) 


# inspect the solution
objective_function_value = solver.ObjectiveValue()

[(i, solver.Value(dv_select_items[i])) for i in dv_select_items]
# 2st, 3th and 4th item gets selected

Спасибо, что приложили усилия. Тем не менее, это удаляет самые прибыльные элементы. Для данных примера ожидаемое значение равно 0,18. Если убрать пункты 2, 4 и 6, ожидаемое значение снизится до 0,11. Я посмотрю дальше, чтобы найти решение!

PoE Academy 08.02.2023 16:38

Я не смог заставить ваш код работать, но повторил его. Результаты в другом комментарии в этой теме. Не могли бы вы проверить его и понять, почему он все еще не работает для меня? Целевая функция всегда возвращает 0.

PoE Academy 08.02.2023 22:22

Спасибо за ответ. Я понял, что одно ограничение было указано неправильно - "model.AddMultiplicationEquality(y_x_i_w_i, [x_i_w_i, y])" (это правильная версия, раньше была ошибка). Я предоставил исправленный код выше. Теперь выбираются элементы 2,3,4. Имеет ли смысл этот выбор?

Bhartendu Awasthi 09.02.2023 18:58

Я пробовал следующее:

# new data!
data = {"w_default": {
    "0": 500,
    "1": 250,
    "2": 500,
    "3": 1000,
    "4": 1000,
    "5": 500
},
    "chaos": {
        "0": 8,
        "1": 5,
        "2": 5,
        "3": 4,
        "4": 4,
        "5": 2
    }
}

UB = 10000000
model = cp.CpModel()

num_items = len(data["w_default"])

# create boolean coefficients
dv_select_items = {i: model.NewBoolVar("item_" + str(i)) for i in data["w_default"]}

# constraint: remove exactly 3 items
# TODO: remove exactly 3 items because right now the objective function does not work yet (even without removing any items)
model.Add(sum(dv_select_items[i] for i in dv_select_items) == num_items)

##### numerator equation #####
# x_i * w_i * p_i
xi_wi_pi = model.NewIntVar(0, UB, "xi_wi_pi")
model.Add(xi_wi_pi == sum(dv_select_items[i] * data["w_default"][i] * data["chaos"][i] for i in dv_select_items))

##### denominator equation #####
xi_wi = model.NewIntVar(0, UB, "xi_wi")
model.Add(xi_wi == sum(dv_select_items[i] * data["w_default"][i] for i in dv_select_items))

y_xi_wi_pi = model.NewIntVar(0, UB, "y_xi_wi_pi")
model.AddDivisionEquality(y_xi_wi_pi, xi_wi_pi, xi_wi)

# set target
model.Maximize(xi_wi_pi)

solver = cp.CpSolver()
solver.Solve(model)

# inspect the solution
objective_function_value = solver.ObjectiveValue()
solution_info = solver.SolutionInfo()

result = [(i, solver.Value(dv_select_items[i])) for i in dv_select_items]

Кажется, что проблема заключается в этой части:

##### denominator equation #####
xi_wi = model.NewIntVar(0, UB, "xi_wi")
model.Add(xi_wi == sum(dv_select_items[i] * data["w_default"][i] for i in dv_select_items))

y_xi_wi_pi = model.NewIntVar(0, UB, "y_xi_wi_pi")
model.AddDivisionEquality(y_xi_wi_pi, xi_wi_pi, xi_wi)

Как бы я это ни перефразировал (попробовал предложение @Bhartendu Awasthi с

model.AddMultiplicationEquality(y_x_i_w_i, [y_x_i_w_i, y])

сначала но целевая функция всегда возвращает 0, хотя так и должно быть (ничего не удаляя)

xi_wi_pi = 16750 
xi_wi = 3750
objective_function_value = 16750 / 3750 = 4,466

Если исключить пункты 3, 4 и 5, результат должен быть

objective_function_value = 7750 / 1250 = 6,2

На данный момент у меня нет решения, но я постараюсь решить его в ближайшие дни.

Ответ принят как подходящий

Вот мое решение:

1. Сделайте результаты отслеживаемыми

В cpmodel от ortools вы можете сделать это, введя класс с

from ortools.sat.python import cp_model as cp 

class VarArraySolutionPrinter(cp.CpSolverSolutionCallback):
"""Print intermediate solutions."""

def __init__(self, variables):
    cp.CpSolverSolutionCallback.__init__(self)
    self.__variables = variables
    self.__solution_count = 0

def on_solution_callback(self):
    self.__solution_count += 1
    for v in self.__variables:
        print('%s=%i' % (v, self.Value(v)), end=' ')
    print()

def solution_count(self):
    return self.__solution_count

и чем структурировать свой код следующим образом

# ############################################
# content of your objective function goes here
# ############################################

# solve the model
solver = cp.CpSolver()
status = solver.Solve(model)

# https://developers.google.com/optimization/cp/cp_solver?hl=de
# debugging
solution_printer = VarArraySolutionPrinter([objective, xi_wi_pi, xi_wi])
solver.parameters.enumerate_all_solutions = True

# inspect the solution
objective_function_value = solver.ObjectiveValue()
solution_info = solver.SolutionInfo()
status = solver.Solve(model, solution_printer)

обратите внимание, что в

solution_printer = VarArraySolutionPrinter([objective, xi_wi_pi, xi_wi])

вы хотите добавить имена переменных, т.е. используя третий аргумент (str) при создании переменной:

xi_wi = model.NewIntVar(0, 100, "xi_wi")

2. Создайте модель

После этого я обнаружил, что мне не нужно следовать совету Джониса, потому что cpmodel из ortool может сам обрабатывать двоичные переменные. Этот код работает для меня:

from ortools.sat.python import cp_model as cp

# w_default = weighting
# chaos = price
data = {"w_default":{
  "0":500,
  "1":50,
  "2":50,
  "3":50,
  "4":250,
  "5":1000
},
"chaos":{
  "0":4,
  "1":78,
  "2":75,
  "3":170,
  "4":5,
  "5":4
   }
}

model = cp.CpModel()
num_items = len(data["w_default"])

# create boolean coefficients
dv_select_items = {i: model.NewBoolVar("item_" + str(i)) for i in data["w_default"]}

# constraint: remove exactly 3 items
model.Add(sum(dv_select_items[i] for i in dv_select_items) == num_items - 3)

# ##### numerator equation #####
constant = 1000
# x_i * w_i * p_i // sum of weightings * prices = 200.000 -> UB 500.000 to give some space?
xi_wi_pi = model.NewIntVar(50000 * constant, 500000 * constant, "xi_wi_pi")
model.Add(xi_wi_pi == sum(
    dv_select_items[i] * data["w_default"][i] * data["chaos"][i] * constant for i in dv_select_items))

##### denominator equation #####
# xi_wi // sum of weightings 23665: 20665 with 3 blocked
lb_weights = int(tot_weight * 0.75)
xi_wi = model.NewIntVar(lb_weights, tot_weight, "xi_wi")
model.Add(xi_wi == sum(dv_select_items[i] * data["w_default"][i] for i in dv_select_items))

objective = model.NewIntVar(0, 100 * constant, "objective")
model.AddDivisionEquality(objective, xi_wi_pi, xi_wi)

# set target
model.Maximize(objective)

# solve the model
solver = cp.CpSolver()
status = solver.Solve(model)

# https://developers.google.com/optimization/cp/cp_solver?hl=de
# debugging
solution_printer = VarArraySolutionPrinter([objective, xi_wi_pi, xi_wi])
solver.parameters.enumerate_all_solutions = True

# inspect the solution
objective_function_value = solver.ObjectiveValue()
solution_info = solver.SolutionInfo()
status = solver.Solve(model, solution_printer)

3. Масштабируйте числитель, потому что AddDivisionEquality округляется до целых чисел.

Вам может быть интересно, почему я добавил в код константу (без нее она не работает). Это потому что

model.AddDivisionEquality(objective, xi_wi_pi, xi_wi)

всегда округляет результат до целого числа, и поскольку результаты находятся в диапазоне 8.something, целевая функция всегда возвращает 8. Однако, если вы умножите числитель на 1000, 8,3456 теперь станет 8345, а 8,4334 станет 8433 и, таким образом, может быть оценено правильно.

Надеюсь, это поможет любому, у кого есть похожая проблема. Кроме того, большое спасибо Джони и Бартенду за то, что указали мне правильное направление!

Другие вопросы по теме