Linear Programming problem for finding supplier

Question

I have a linear programming problem in front of me that I am searching to solve in R (any package). It seems like an unbalanced assignment problem to me with relaxed constraints on columns (repeated assignment is possible) or more general possibly a transportation problem:

I got 188 materials that I can get from 7 suppliers. Not all suppliers can provide all of the materials but each just some subset of the material demand.

What I know is

1. the demand for each of the material
2. how much I can order for the month at each supplier
3. which supplier can provide which material in which quantity.

One can imagine the sample problem as a table where the marginals are known, the x means that I can order at that supplier(s), the quantity could be split into smaller chunks, costs let us say are constant for all material/supplier combination e.g. 1 (or could be 33):

   Sup1  Sup2 Sup3  (Demand-Quantity of the material to get)
m1  x               0.12
m2  x    x          1.8
m3            x     1.9
m4       x          2.1
m5 x     x    x     0.1
   30    50   12
   (Capacity of suppliers 1 to 3)

Hence the problem is to minimise the order (cost) given I get all materials and not cross over the capacity of each suppliers. I need to get, given the demand, capacity constraints, a cost matrix and demand satisfaction matrix (which supplier provides which material), which quantities of each material to get at each supplier.

A concrete case, demand and capacity in tons (cost are constant for demonstration, all = 1):

demand <- c(10.994, 0.322, 0.046, 1.449, 0.253, 0.368, 5.221, 2.208, 0.989, 
            11.983, 0.046, 0.092, 0.115, 0.023, 1.633, 0.253, 0.276, 0.483, 
            0.966, 8.763, 0.276, 1.173, 17.963, 0.943, 4.232, 0.575, 23.138, 
            0.207, 10.465, 0.138, 0.92, 11.247, 6.486, 4.807, 16.606, 0.966, 
            4.554, 1.15, 0.644, 0.046, 0.046, 0.023, 0.989, 1.288, 0.966, 
            13.34, 5.796, 12.213, 10.097, 5.865, 5.704, 38.916, 7.291, 18.193, 
            0.184, 4.439, 58.305, 0.115, 1.38, 3.22, 8.234, 5.29, 2.369, 
            0.736, 2.668, 4.968, 3.496, 0.552, 0.598, 4.922, 0.322, 0.161, 
            0.483, 0.644, 0.736, 0.069, 0.644, 14.789, 3.243, 1.242, 3.381, 
            41.423, 0.644, 56.603, 3.841, 9.154, 0.23, 2.07, 40.871, 0.575, 
            0.276, 0.276, 4.255, 0.092, 0.253, 0.161, 0.115, 0.23, 0.092, 
            0.138, 0.874, 0.184, 0.483, 0.115, 0.161, 0.069, 0.184, 1.633, 
            0.184, 0.023, 0.621, 0.046, 0.184, 2.415, 4.301, 0.138, 0.667, 
            3.818, 0.138, 0.276, 0.023, 0.874, 1.702, 0.253, 0.299, 0.552, 
            7.153, 0.966, 1.702, 1.104, 0.713, 4.14, 0.575, 0.322, 0.115, 
            0.621, 1.932, 0.759, 0.69, 0.207, 0.759, 0.966, 0.023, 3.358, 
            1.357, 0.207, 1.242, 0.138, 1.265, 7.958, 0.529, 2.53, 0.552, 
            0.552, 0.046, 0.207, 0.069, 0.368, 0.092, 0.023, 0.966, 1.633, 
            3.151, 2.53, 9.2, 0.598, 0.138, 1.173, 0.345, 2.001, 0.506, 0.598, 
            1.104, 4.278, 0.345, 1.104, 0.253, 0.736, 0.138, 0.276, 4.761, 
            0.851, 2.622, 0.345, 0.92, 1.357, 0.092, 0.368)
n_demand <- length(demand)
capacity = c(55, 70, 60, 51, 10, 90, 10, 10, 10)
n_capacity <- length(capacity)
material_supplier_mat <- structure(c(1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 
0L, 0L, 0L, 0L, 1L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 
1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 
0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 
0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 
0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 
1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 
1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 
1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 
0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 
0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 
1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 
1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 
0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 
1L, 0L, 0L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 1L, 1L, 1L, 1L, 
1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 
1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 
1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 
0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 
0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 
1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 
1L, 1L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 
1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 
1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 0L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 
0L, 1L, 1L, 1L, 1L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 
0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 
0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 
1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 
0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 
1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 
0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 
0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 
1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 
0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 
0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 
0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 
1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 
0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 
1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 
0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 
0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 
0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 
0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 
0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 
1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 
1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 
1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 
0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 
0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 
1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 
0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 
0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 
0L, 0L, 0L, 1L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 1L, 
1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 
1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 
1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 
0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 
0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 
1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 
0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 
0L, 1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 
0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 
0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 
1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L), .Dim = c(9L, 
188L))

I tried using ompr like follows

library(ompr)
library(ROI)
library(ROI.plugin.glpk)
library(ompr.roi)
model <- MIPModel() %>%
  add_variable(x[i], i = seq_len(n_demand), type = 'continuous') %>%
  add_variable(b[i], i = seq_len(n_demand), type = 'binary') %>%
  set_objective(sum_expr(x[i], i = seq_len(n_demand)), sense = 'min')
for (el in seq_len(n_demand)) {
  model %<>% 
    set_bounds(x[el], lb = demand[el]) 
}
for(el in seq_len(n_capacity)) {
column_index <- which(material_supplier_mat[el, ] == 1L) # the sum of these columns/material quantities that need to be taken into consideration 
  cat("column_index:" , column_index, "\n")
  nr_columns <- sum(material_supplier_mat[el, ] == 1L) # for the constraint to allow up to 'nr_columns' to be tried to be taken into account7
  cat("nr_columns:" , nr_columns, rep.int("\n", 2L))
model %<>%
    add_constraint(sum_expr(demand[i] * b[i]) <= capacity[el], i = column_index) %>%
    add_constraint(sum_expr(b[i]) <= nr_columns, i = column_index)
}
result <- model %>% 
  solve_model(with_ROI(solver = "glpk", verbose = TRUE))

But I don't think it is the right formulation and solution, at the end I was hoping to get the objective value as well as decision of the material quantities purchase for each supplier (for instance get m2 0.8t at S1 and 1.0t at S2).

Answer 1

The problem you mentioned sounds like a variant of the lot-sizing problem in which there are multiple sources to produce one of the materials. Let's say, the variable $x_{i,s,t}$ represents the amount of material $i$ that is produced by supplier $s$ in period $t$. Also, the binary variable $y_{i,t}$ is being introduced to ensure that if any material $i$ will be started in the period $t$ is $1$, otherwise is $0$.

\begin{align} \min&\quad Z=\sum_{i=1}^{I} \sum_{s=1}^{S} \sum_{t=1}^{T} c_{i,s,t}\cdot x_{i,s,t}\\\text{s.t.}&\quad \sum_{s=1}^{S} x_{i,s,t} \geq \text{demand}_{i,t} \quad ,\forall i=1, \ldots, I, \forall t=1, \ldots, T & \text{(1)}\\&\quad x_{i,s,t} \leq y_{i,t}\cdot\text{lotsize}_{i} \quad ,\forall i=1, \ldots, I, \forall s=1, \ldots, S, \forall t=1, \ldots, T & \text{(2)}\\&\quad \sum_{i=1}^{I} x_{i,s,t}\cdot a_{i,t} \leq \text{cap}_{s} \quad ,\forall s=1, \ldots, S, \forall t=1, \ldots, T & \text{(3)}\\&\quad x_{i,t,s} \in ℝ , y_{i,t} \in \{0,1\} \end{align}

where, $c_{i,s,t}$ is corresponding order cost, $\text{demand}_{i,t}$ is demand of material $i$ in period $t$, $\text{lotsize}_{i}$ is the minimum lot/batch size that should be ordered, $a_{i,t}$ is the required capacity to produce materials and $\text{cap}_{s}$ is the capacity of supplier $s$. Indeed one can define a constraint $y_{i,t} \geq x_{i,s,t}/M$ to set a lower bound on the variable $y$.

Also, to extend the model into the standard material requirement planning (MRP), you can define some extra constraints to capture inventory level based on whose lead time.

Answer 2

Your code now runs, but I'm fairly certain the model is wrong.

To get the values of x as a vector, you can execute get_solution(result, x[i])[, "value"]. To get it as a data frame (with columns for variable name, index and value), just use get_solution(result, x[i]). The same approach works for b. The expression result$objective_value will get you the optimal objective value.

For your model, the solution is trivial. Each x[i] is equal to demand[i] because you set that as a lower bound. Every b[i] is 0.

The description of your problem is a bit unclear (I have no idea what b[] represents), but you might want to have a look at the "transportation problem" (which applies to more things than just transportation of goods). An easy introduction is here.

Answer 3

As @prubin corretly pointed out, the essential problem at hand is the transportation problem with some additional constraints on non-available quantities for relevant material-supplier pairs.

I could model the data using ompr as follows:

  model <- MIPModel() %>%
    # quantity to get from the suppliers
    add_variable(x[i, j], i = 1:n_capacity, j = 1:n_demand, type = 'continuous') %>%
    # objective function
    set_objective(sum_over(c[i, j] * x[i, j], i = 1:n_capacity, j = 1:n_demand), 'min') %>%
    # constraints
    add_constraint(sum_over(x[i, j], i = 1:n_capacity) >= demand[j], j = 1:n_demand) %>%
    add_constraint(sum_over(x[i, j], j = 1:n_demand) <= capacity[i], i = 1:n_capacity) %>%
    add_constraint(x[i, j] >= 0, i = 1:n_capacity, j = 1:n_demand)
zero_constraint_on_quantity <- which(c == 0, arr.ind = TRUE)
  for(el in seq_len(NROW(zero_constraint_on_quantity))) {
    i_val <- zero_constraint_on_quantity[el, 'row']
    j_val <- zero_constraint_on_quantity[el, 'col']
    model %<>%
      add_constraint(x[i, j] == 0, i = i_val, j = j_val)
  }

where c is an aditional cost matrix, with cost for every material-supplier pair.

Additionaly, I find it interessting that an implementation of an elastic filter (out of scope for this question) can calculate over/under utilisation of the capacity, which is quite handy for further SCM planning.