Я пытаюсь оптимизировать бинарную проблему для моего веб-сайта.
Данные содержат примерно 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, но не смог заставить его работать, так как я застрял, переводя свою проблему в код. Может кто-нибудь мне помочь?
Используя комментарии @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. Я посмотрю дальше, чтобы найти решение!
Я не смог заставить ваш код работать, но повторил его. Результаты в другом комментарии в этой теме. Не могли бы вы проверить его и понять, почему он все еще не работает для меня? Целевая функция всегда возвращает 0.
Спасибо за ответ. Я понял, что одно ограничение было указано неправильно - "model.AddMultiplicationEquality(y_x_i_w_i, [x_i_w_i, y])" (это правильная версия, раньше была ошибка). Я предоставил исправленный код выше. Теперь выбираются элементы 2,3,4. Имеет ли смысл этот выбор?
Я пробовал следующее:
# 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 и, таким образом, может быть оценено правильно.
Надеюсь, это поможет любому, у кого есть похожая проблема. Кроме того, большое спасибо Джони и Бартенду за то, что указали мне правильное направление!
Прежде чем писать код, вы должны правильно сформулировать задачу оптимизации на бумаге. Для начала: введите двоичную переменную 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, см. этот пост для более подробной информации.