Computation

• Imports
  

  import numpy as np
  from ortools.sat.python import cp_model as cp
  

• Input
  
  \begin{equation*} A = \begin{pmatrix} a_{11} & & \cdots & & a_{1n} \\ & \ddots & & & \\ \vdots & & a_{ij} & & \vdots \\ & & & \ddots & \\ a_{m1} & & \cdots & & a_{mn} \end{pmatrix} \in {\{0, 1\}}^{m \times n} \end{equation*}

  table = np.array(TABLE)[2:, 1:]
  A = np.where(table == '', 0, 1)
  A
  

  1	0	1	0	1	0	0	0
  0	1	1	0	1	1	0	0
  0	0	0	1	0	0	1	0
  1	0	0	0	1	0	0	1
  0	1	0	1	0	0	1	0
  1	0	1	1	0	0	0	1

  m = A.shape[0]
  n = A.shape[1]
  (m, n)
  

  (6, 8)
  

• Model
  

  model = cp.CpModel()
  

• Variables
  

  The decision variables represent the inclusion of each of the \(m\) suppliers,

  \begin{equation*} x \in {\{ 0, 1 \}}^m. \end{equation*}

  x = [model.new_bool_var(name="x_%d" % i) for i in range(m)]
  

• Constraints
  

  Every part is covered by some supplier,

  \begin{equation*} \forall j \sum_{i = 1}^{m} a_{ij} x_i \leq 1. \end{equation*}

  for j in range(n):
      model.add(sum(A[i, j] * x[i] for i in range(m)) <= 1)
  

• Objective
  

  Define the objective function

  \begin{equation*} y = \sum_{i = 1}^{m} x_i \end{equation*}

  and the objective value

  \begin{equation*} y^* = \max y. \end{equation*}

  y = sum(x[i] for i in range(m))
  model.maximize(y)
  

• Solution
  

  Find the optimum solution.

  solver = cp.CpSolver()
  status = solver.solve(model)
  assert status == cp.OPTIMAL
  y_star = int(solver.objective_value)
  print("y = %s" % y_star)
  print("x = %s" % [solver.value(x_i) for x_i in x])
  

  y = 2
  x = [0, 1, 1, 0, 0, 0]
  

• Re-parameterize
  

  Prepare for optimal solution enumeration.

  model.clear_objective()
  model.add(y == y_star)
  

• Find all optimal solutions
  

  class SolutionPrinter(cp.CpSolverSolutionCallback):
      def __init__(self):
          cp.CpSolverSolutionCallback.__init__(self)
          self.__solution_count = 0
  
      def on_solution_callback(self) -> None:
          self.__solution_count += 1
          if self.__solution_count > 1: print()
          print("Solution #%d" % self.__solution_count)
          print("x = %s" % [solver.value(x_i) for x_i in x])
          print("y = %s" % y_star)
  
  printer = SolutionPrinter()
  status = solver.solve(model, printer)
  

  Solution #1
  x = [0, 1, 1, 0, 0, 0]
  y = 2