# Reducing number of suppliers for product portfolio

I have the following matrix of suppliers who are able to make a certain product, against all products in my portfolio.

What is the best way of finding the solution to "the least suppliers necessary to deliver the whole portfolio" - and "which suppliers are necessary to deliver the whole portfolio"?

Ideally looking for a solution in R, since I have manipulated the data beforehand in order to get to this one-hot coded matrix. But generally I'm trying to understand first how to approach this.

The full dataset is obviously larger - I can see myself that in this case only Sup1 and Sup7 would be required :-)

I received some great answers below. Unfortunately I cannot install the "pulp" module for python on my laptop (restrictions from work - I know it doesn't seem to make sense). I am trying to convert the below python script into R (ompr by @dirks user:2798441). However I struggle with the syntax. I receive an error message with below code:

require(ompr)
require(ompr.roi)
require(dplyr)
require(ROI)
require(ROI.plugin.glpk)

test <-rbind(c(),
c("Prd1", 1, NA, NA, NA, NA, NA, NA, 1, NA, NA, NA),
c("Prd2",1, NA, NA, NA, NA, 1, NA, NA, NA, 1, NA),
c("Prd3",NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA),
c("Prd4",1, 1, 1, 1, NA, NA, NA, NA, NA, NA, NA),
c("Prd5",NA, NA, NA, NA, NA, NA, 1, NA, NA, NA, NA),
c("Prd6",1, NA, NA, NA, NA, 1, NA, NA, NA, 1, NA),
c("Prd7",1, NA, NA, NA, 1, NA, NA, NA, 1, NA, NA),
c("Prd8",NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA),
c("Prd9",NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA))
dims <- dim(test)
test <- as.numeric(test)
dim(test) <- dims
test.df <- data.frame(test)
colnames(test.df) <- c("Product","Sup1", "Sup2", "Sup3", "Sup4", "Sup5", "Sup6", "Sup7", "Sup8", "Sup9", "Sup10", "Sup11")
test.df[,1] <- c("Prd1","Prd2","Prd3","Prd4","Prd5","Prd6","Prd7","Prd8","Prd9")

n <- nrow(test.df)
m <- ncol(test.df)
M <- 100
set_n <- range(0,n)
set_m <- range(0,m)

model <- MIPModel() %>%
add_variable(b[i,j], i = set_n, j = set_m, type = "binary") %>%
add_variable(x[j], j = set_m, type = "binary") %>%
set_objective(sum_expr(b[i,j] * 2, i = set_n, j = set_m) - sum_expr(x[j], j = set_m)) %>%
add_constraint((sum_expr(b[i,j], j = set_m)) <= 1) %>%
add_constraint(M * x[j] >= sum_expr(b[i,j], i = set_n) - 1 + 0.001) %>%
add_constraint(M * ( 1-  x[j]) >= ( 1 - sum_expr(b[i,j], i = set_n) - 0.001)) %>%
solve_model(with_ROI(solver = "symphony", verbosity = 1)) %>%
get_solution(x[i, j]) %>%
filter(value > 0) %>%
arrange(i)


If someone has some kind of experience with ompr I would appreciate a nudge in the right direction.

• Is there a parameter regarding the cost every supplier charges for delivering, or that doesn't really matter in this case? – dhasson Jan 30 '20 at 15:26
• In this first stage the idea is merely to reduce the number of suppliers. Once I figure this out, I will work on a model where a price for each product/supplier is considered. I just saw some really promising answer below, and will try to implement and test them asap. – Roman Jan 31 '20 at 14:16

You can model this as a set covering problem (or hitting set problem -- different terminology, same mathematical model). The (binary) decision variables would determine which suppliers you select; the constraints would be "select at least one that can handle this product" for each product. The objective would be to minimize the number of selections (sum of the binary variables). Note that this model does not take into account what the suppliers charge, nor their quality (other than it presumably must be adequate to be included in the matrix), nor does it consider any supplier capacity limits. It also leaves unanswered how to allocate demand among multiple suppliers if you wind up with more than covering some product.

As for R, you can create and solve the model in the ompr package (available on CRAN), but you will also need to install one of the solvers it supports.

Let $$b_{i,j}$$ denote product $$i$$ is taken from supplier $$j$$. and $$x_j$$ denote whether supplier $$j$$ is selected once or not. $$i \in n$$ (product) $$j \in m$$ (supplier)

Now:

\begin{align}\max&\quad \sum_{i,j} b_{i,j}\times 2 - \sum_j x_j\\\text{s.t.}&\quad\sum_j b_{i,j} \leq 1\\&\quad M \times x_j \geq \sum_i b_{i,j} -1 + \delta\\&\quad M \times (1- x_j) \geq 1 - \sum_i b_{i,j} - \delta\\&\quad b_{i,j} =0 \; \forall A[i,j] = 0\end{align}

$$\delta$$ is a small value.

The objective function tries to maximize the number of product picked while keeping less number of suppliers. $$A[i,j]$$ is your original matix.

$$b_{i,j}$$ is multiplied by 2 because when $$x_j = 1$$ and $$\sum_i b_{i,j} = 1$$ solver will end up ignoring it as $$x_j - \sum_i b_{i,j} = 0$$ at that time.

Here is the code but in python, you can convert it to R::

import numpy as np
import pulp as pl

p = [[1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0],
[1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0],
[1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0],
[1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]

p = np.array(p)

n = 9
m = 11
set_n = range(n)
set_m = range(m)

prob = pl.LpProblem("Problem", pl.LpMaximize)
b = {(i, j): pl.LpVariable("b_{0}_{1}".format(i, j), 0, 1, pl.LpBinary) for i in set_n for j in set_m}
x = {(j): pl.LpVariable("x_{0}".format(j), 0, 1, pl.LpBinary) for j in set_m}

prob += pl.lpSum(b[i, j] * 2 for i in set_n for j in set_m) - pl.lpSum(x[j] for j in set_m)

M = 100
for i in set_n:
prob += pl.lpSum(b[i, j] for j in set_m) <= 1

for j in set_m:
prob += M * x[j] >= (pl.lpSum(b[i, j] for i in set_n) - 1 + 0.001)
prob += M * (1 - x[j]) >= (1 - pl.lpSum(b[i, j] for i in set_n) - 0.001)

for i in set_n:
for j in set_m:
if p[i, j] == 0:
prob += b[i, j] == 0

# print(prob)
prob.solve()
print(prob.status)
for v in prob.variables():
if v.varValue == 1.0:
print(v.name, "=", v.varValue)


Please try the below code. I modified the way you import i, j vectors. Also, check the initial matrix on how products are linked with suppliers, there are suppliers with no product and vice versa. Then you should fix some b[i,j] values accordingly when NA in initial matrix 0 in b[i,j].

 model <- MIPModel() %>%
add_variable(b[i,j], i = 1:n, j = 1:m, type = "binary") %>%
add_variable(x[j], j = 1:m, type = "binary") %>%
set_objective(sum_expr(b[i,j] * 2, i = 1:n, j = 1:m) - sum_expr(x[j], j = 1:m)) %>%
add_constraint((sum_expr(b[i,j], j = 1:m)) <= 1, i = 1:n) %>%
add_constraint(M * x[j] >= sum_expr(b[i,j], i = 1:n) - 1 + 0.001, j = 1:m) %>%
add_constraint(M * ( 1-  x[j]) >= ( 1 - sum_expr(b[i,j], i = 1:n) - 0.001), j = 1:m) %>%

solve_model(with_ROI(solver = "symphony", verbosity = 1))

get_solution(model,b[i,j])
get_solution(model,x[j])