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

Question

我正在尝试求解已存在初始解的最小化问题，并且 objective 函数基于此初始解。

我有一些行 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) == 1 用于所有 i in i_idx。我的 objective 是，从最初的 y_line 开始，每个售出的资源都为我们带来了节省，每个新购买的资源都让我们付出了代价，并且解决方案（新线路）y 有运营的处理成本。我定义了 objective 如下：

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)

但是，如果我求解此模型（完整代码见下文），则 objective 值 (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

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 模型我不知道如何解决它。

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

Minimization model with initial starting values

python

or-tools