python - 具有初始起始值的最小化模型

Question

我正在尝试解决已经存在初始解决方案并且目标函数基于此初始解决方案的最小化问题。

我有某种线y_line，它是资源和站点的初始映射：

y_line = np.array([[1, 0, 0],
                   [0, 1, 0],
                   [0, 0, 1]])

此外，我有一个用于从生产线销售的储蓄阵列S，一个用于购买新EC产品和用于加工的阵列P

S = np.array([[-260., -260., -260.],
              [-30.,  -30.,  -30.],
              [360.,  360.,  360.]], dtype=int)
EC = np.array([[1000, 1000, 1000],
               [2000, 2000, 2000],
               [5000, 5000, 5000]], dtype=int)
P = np.array([[720.,  720.,  720.],
              [1440., 1440., 1440.],
              [3600., 3600., 3600.]], dtype=int)

仅使用一个简化的约束：每个工作站i必须至少有一个资源j-> sum(y[i, j] for j in j_idx) == 1for all i in i_idx。我的目标是，从最初售出的每一个资源都会y_line为我们带来节省，每一次新购买的资源都会给我们带来成本，而解决方案（新生产线）y则需要运营处理成本。我将目标定义如下：

y_delta = y - y_line  # delta between new line (y) and old line (y_line)
y_delta_plus = np.zeros(y.shape, dtype=object)  # 1
y_delta_minus = np.zeros(y.shape, dtype=object)  # 2

# I -> new bought resources
y_delta_plus[y_delta >= 0] = y_delta[y_delta >= 0]
# II -> sold resources
y_delta_minus[y_delta <= 0] = y_delta[y_delta <= 0]

c_i = y_delta_plus * EC  # invest
c_s = y_delta_minus * S  # savings
c_p = y * P  # processing cost
c_y = np.sum(c_s + c_i + c_p)

但是，如果我解决了这个模型（完整代码见下文），那么目标值 ( 5760) 与我的健全性检查计算 ( ) 不匹配12430。是否可以为 设置初始值y[i, j]？还是有其他功能可以实现这一点？

from ortools.linear_solver import pywraplp
import numpy as np


y_line = np.array([[1, 0, 0],
                   [0, 1, 0],
                   [0, 0, 1]])
S = np.array([[-260., -260., -260.],
              [-30.,  -30.,  -30.],
              [360.,  360.,  360.]], dtype=int)
EC = np.array([[1000, 1000, 1000],
               [2000, 2000, 2000],
               [5000, 5000, 5000]], dtype=int)
P = np.array([[720.,  720.,  720.],
              [1440., 1440., 1440.],
              [3600., 3600., 3600.]], dtype=int)

solver = pywraplp.Solver('stack', pywraplp.Solver.SAT_INTEGER_PROGRAMMING)

y = np.zeros_like(y_line, dtype=object)

i_idx = range(y_line.shape[0])
j_idx = range(y_line.shape[1])

for i in i_idx:
    for j in j_idx:
        y[i, j] = solver.IntVar(0, 1, 'y[%i_%i]' % (i, j))

for i in i_idx:
    solver.Add(
        sum(y[i, j] for j in j_idx) == 1
    )


def objective(y, y_line):
    y_delta = y - y_line  # delta between new line (y) and old line (y_line)
    y_delta_plus = np.zeros(y.shape, dtype=object)  # 1
    y_delta_minus = np.zeros(y.shape, dtype=object)  # 2

    # I -> new bought resources
    y_delta_plus[y_delta >= 0] = y_delta[y_delta >= 0]
    # II -> sold resources
    y_delta_minus[y_delta <= 0] = y_delta[y_delta <= 0]

    c_i = y_delta_plus * EC  # invest
    c_s = y_delta_minus * S  # savings
    c_p = y * P  # processing
    return np.sum(c_s + c_i + c_p)

c_y = objective(y=y, y_line=y_line)
solver.Minimize(
    c_y
)

# [START solve]
print("Number of constraints:", solver.NumConstraints())
print("Number of variables:", solver.NumVariables())
status = solver.Solve()
# [END solve]

y_new = np.zeros_like(y)
for i in range(y_line.shape[0]):
    for j in range(y_line.shape[1]):
        if y[i, j].solution_value() > 0:
            y_new[i, j] = y[i, j].solution_value()

print(f"Objective sat: {solver.Objective().Value()}")
print(y_new)

# Number of constraints: 3
# Number of variables: 9
# Objective sat: 5760.0
# [[1.0 0 0]
#  [1.0 0 0]
#  [1.0 0 0]]

# %%
c_y_test = objective(y=y_new, y_line=y_line)
c_y_test # -> 12430.0

标签： pythonor-tools

Answer 1

模型可以解决。然而，不是用这种方法，我首先选择了。使用pywraplp模型它不起作用，但是cp_model可以使用预定义的变量（如@sascha 提到的）来解决它。数组y_line、S和EC与P上述相同。庄严的约束也是一样的。然而，我可以使用以下方法解决的“过滤”：

for i in range(len(y_cp.flatten())):
    model.AddElement(i, y_delta.flatten().tolist(), y_cp.flatten().tolist()[i] - y_line.flatten().tolist()[i])

for i in i_idx:
    for j in j_idx:
        model.AddMaxEquality(y_delta_plus[i, j], [y_delta[i, j], model.NewConstant(0)])
        model.AddMinEquality(y_delta_minus[i, j], [y_delta[i, j], model.NewConstant(0)])

model.Minimize(
    np.sum(y_delta_plus * EC) + np.sum(y_delta_minus * S) + np.sum(y_cp * P)
)

求解和完整性检查产生：

solver_cp = cp_model.CpSolver()
solver_cp.Solve(model)

y_new_cp = np.zeros_like(y_cp)
for i in i_idx:
    for j in j_idx:
        if solver_cp.Value(y_cp[i, j]) > 0:
            y_new_cp[i, j] = solver_cp.Value(y_cp[i, j])

print(f"Objective cp: {solver_cp.ObjectiveValue()}")
print(y_new_cp)

# Objective cp: 5760.0
# [[1 0 0]
#  [0 1 0]
#  [1 0 0]]

c_y_test = objective(y=y_new_cp, y_line=y_line)
c_y_test # -> 5760 -> Correct

可以解决它cp_model并匹配健全性检查。使用该pywraplp模型，我无法弄清楚如何解决它。