python-3.x 如何在ortools /cp模型中用线性规划求解二元优化问题？

e3bfsja2 于 2023-02-14 发布在 Python

关注(0)|答案(3)|浏览(248)

我正在为我的网站优化一个二进制问题。
数据包含大约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
   }
}

我用下式计算整个数据集的期望值
瓦尔=（w1 p1 + w2 p2 +...+ wn pn）/总和（w1 + w2 +... wn）
与
总和（w1 + w2 +... wn）= 23665（考虑所有项目）
到目前为止一切顺利，但现在来了棘手的部分。并不是所有的项目都是想要的，也就是说，他们的价值较低和/或有一个高权重，稀释了池，我可以借鉴。
通过“阻塞”或删除最多3个项目，我只能从剩余项目中提取，这样做可以最大化加速值函数。哪些项目要删除？由于价格随着时间的推移而变化，我必须定期检查要删除的项目。
我已经开始简单地删除具有最高权重和最低价格的项目，但我敢肯定，这只代表了一个局部最优，并会有一个更优的战略。
在查看了一些网站后，似乎混合整数线性规划（MILP）或特别是BILP（binary...）可以解决我的问题，我找到了https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.milp.html，但无法使其工作，因为我坚持将我的问题转换为代码。有人能帮助我吗？

python-3.x

来源：https://stackoverflow.com/questions/75038973/how-to-solve-a-binary-optimization-problem-with-linear-programming-in-ortools

3条答案

按热度按时间

6xfqseft1#

我尝试了以下方法：

# 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，尽管它应该返回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

我目前没有解决办法，但会在未来几天内尝试解决它。

赞(0）回复(0）举报 2023-02-14

ddarikpa2#

利用@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）回复(0）举报 2023-02-14

qojgxg4l3#

我的解决方案是：

- 1.使结果可跟踪**

在或tools中的cpmodel中，您可以通过使用

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.创建模型**

这样一来，我发现我不必遵循Jonis的建议，因为ortool的cpmodel可以自己处理二进制变量，下面的代码很适合我：

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，因此可以正确计算。
希望这对任何有类似问题的人都有帮助。另外，非常感谢乔尼和巴滕杜为我指明了正确的方向!

赞(0）回复(0）举报 2023-02-14

我来回答

python-3.x 如何在ortools /cp模型中用线性规划求解二元优化问题？

3条答案

相关问题

热门标签

最新问答