Large sparse transportation model with CVXPY,CVXR

In [1] I was trying out different formulations of a large, sparse (but very easy) transportation model using different modeling tools and solvers. The conclusion was:

• Exploiting sparsity leads to the best performance
• GAMS/Cplex is about 10 times as fast as Python/PuLP/CBC
• GAMS/Cplex is about 1000 times as fast as R/OMPR/GLPK

Don't overinterpret these numbers: this is a single model, which may not be at all representable for your models. Let's see how CVXPY[2] and CVXR[3] are doing. To recap we want to solve a large (but easy) transportation model:

Dense Transportation Model
\[\begin{align}\min&\sum_{i,j}\color{darkblue}c_{i,j}\cdot\color{darkred}x_{i,j}\\ & \sum_j \color{darkred}x_{i,j} \le \color{darkblue}a_i && \forall i \\& \sum_i \color{darkred}x_{i,j} \ge \color{darkblue}b_j && \forall j \\ & \color{darkred}x_{i,j} \ge 0 \end{align}\]

The additional wrinkle is that only a small fraction of the links \(i \rightarrow j\) exist.

There are (at least) three ways to model this:

1. Use a large cost coefficient for forbidden links.
2. Fix \(\color{darkred}x_{i,j}=0\) for non-existing links.
3. Exploit sparsity directly by only summing over existing links.

CVXPY/CVXR are matrix-based modeling tools. That gives some extra things to think about. A transportation model in matrix-notation can look like:

Transportation Model in Matrix Notation
\[\begin{align}\min\>&{\bf tr}(\color{darkblue}C^T\color{darkred}X)\\&\color{darkred}X\color{darkblue}e \le \color{darkblue}a \\ & \color{darkred}X^t\color{darkblue}e \ge \color{darkblue}b \\ & \color{darkred}X \ge 0 \end{align}\]

A small time arbitrage model

In [1] the following question was posed:

We have monthly prices for a product and I have predictions for the next 12 months. The idea is to exploit this and buy when prices are low and sell when they are high. We have limits on how much we can buy and sell and there are also inventory capacity constraints. What plan is the most profitable?

This is an interesting little problem. We will use R to experiment with this.

Data

The data set is as follows:

# data
price = c(12, 11, 12, 13, 16, 17, 18, 17, 18, 16, 17, 13)
capacity = 25
max_units_buy = 4
max_units_sell = 8

Model

The main vehicle we will use is the well-known inventory balance equation:\[\mathit{inv}_t = \mathit{inv}_{t-1} + \mathit{buy}_t - \mathit{sell}_t\] I assume that the initial inventory is \[\mathit{inv}_0 = 0\]

The optimization model for this problem can look like:

Model 1

\begin{align}\max&\sum_t \color{darkblue}{\mathit{price}}_t \cdot(\color{darkred}{\mathit{sell}}_t-\color{darkred}{\mathit{buy}}_t)\\ &\color{darkred}{\mathit{inv}}_t = \color{darkred}{\mathit{inv}}_{t-1} + \color{darkred}{\mathit{buy}}_t - \color{darkred}{\mathit{sell}}_t\\ &\color{darkred}{\mathit{inv}}_0 = 0\\ &\color{darkred}{\mathit{inv}}_t \in \{0,1,\dots,\color{darkblue}{\mathit{capacity}}\}\\ &\color{darkred}{\mathit{buy}}_t \in \{0,1,\dots,\color{darkblue}{\mathit{maxbuy}}\}\\ &\color{darkred}{\mathit{sell}}_t \in \{0,1,\dots,\color{darkblue}{\mathit{maxsell}}\}\\ \end{align}

I used here integer decision variables: the assumption is we cannot deal with fractional values. Obviously, we can easily relax this and work with continuous variables.

Implimentation

To try this out, I am using CVXR [2] with the Glpk solver [3].

CVXR is using a matrix oriented notation. It is interesting to see how this can be applied to our inventory balance equation. The main trick is that we use a Lag operator implemented as a matrix-vector multiplication: \[ L \cdot \mathit{inv}\]

The Lag operator. We know that multiplication by an identity matrix gives us the same vector: \[ v = I \cdot v \] If we shift the diagonal down by one, we get what we want. Let's do a little experiment:

> n <- 5
> v <- 1:n
> v
[1] 1 2 3 4 5
> # matrix-vector multiplication with identity matrix
> diag(n) %*% v
     [,1]
[1,]    1
[2,]    2
[3,]    3
[4,]    4
[5,]    5
> # lag operator
> L = cbind(rbind(0,diag(n-1)),0)
> L
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    1    0    0    0    0
[3,]    0    1    0    0    0
[4,]    0    0    1    0    0
[5,]    0    0    0    1    0
> # matrix-vector multiplication with Lag operator
> L %*% v
     [,1]
[1,]    0
[2,]    1
[3,]    2
[4,]    3
[5,]    4
> 

We see that \(L\) is indeed the identity matrix but shifted one position down. When we multiply \(Lv\) we see that the result vector is also shifted one down. The newly inserted element at the first position is zero, which is exactly what we want.

Note: if we want a non-zero initial inventory we need to add a vector \(\mathit{initinv}\) with the first element being the initial inventory and the other elements all zero. Then we can write the inventory balance equation in matrix notation as \[\mathit{inv} = \mathit{initinv} + L \cdot \mathit{inv} + \mathit{buy} - \mathit{sell}\] In the R code below, we assumed this was not needed.

With this, we are ready to solve our problem:

> library(CVXR)
> 
> # data
> price = c(12, 11, 12, 13, 16, 17, 18, 17, 18, 16, 17, 13)
> capacity = 25
> max_units_buy = 4
> max_units_sell = 8
> 
> # number of time periods
> NT <- length(price)
> 
> # Decision variables
> inv = Variable(NT,integer=T)
> buy = Variable(NT,integer=T)
> sell = Variable(NT,integer=T)
> 
> # Lag operator
> L = cbind(rbind(0,diag(NT-1)),0)
> 
> # optimization model
> problem <- Problem(Maximize(sum(price*(sell-buy))),
+                    list(inv == L %*% inv + buy - sell,
+                         inv >= 0, inv <= capacity,
+                         buy >= 0, buy <= max_units_buy,
+                         sell >= 0, sell <= max_units_sell))
> result <- solve(problem,verbose=T)
GLPK Simplex Optimizer, v4.47
84 rows, 36 columns, 119 non-zeros
*     0: obj =  0.000000000e+000  infeas = 0.000e+000 (12)
*    35: obj = -1.040000000e+002  infeas = 0.000e+000 (0)
OPTIMAL SOLUTION FOUND
GLPK Integer Optimizer, v4.47
84 rows, 36 columns, 119 non-zeros
36 integer variables, none of which are binary
Integer optimization begins...
+    35: mip =     not found yet >=              -inf        (1; 0)
+    35: >>>>> -1.040000000e+002 >= -1.040000000e+002   0.0% (1; 0)
+    35: mip = -1.040000000e+002 >=     tree is empty   0.0% (0; 1)
INTEGER OPTIMAL SOLUTION FOUND
> cat("status:",result$status)
status: optimal
> cat("objective:",result$value)
objective: 104
> # print results
> data.frame(price,
+            buy=result$getValue(buy),
+            sell=result$getValue(sell),
+            inv=result$getValue(inv))
   price buy sell inv
1     12   4    0   4
2     11   4    0   8
3     12   4    0  12
4     13   4    0  16
5     16   4    0  20
6     17   0    8  12
7     18   0    8   4
8     17   4    0   8
9     18   0    8   0
10    16   4    0   4
11    17   0    4   0
12    13   0    0   0
> 

We see indeed we are buying when things are cheap and selling when prices are high. This gives us a net profit of 104. The inventory capacity is never binding here.

A complication

In a subsequent post, a non-trivial wrinkle was added to the problem [4].

price = c(12, 11, 12, 13, 16, 17, 18, 17, 18, 16, 17, 13)
capacity = 25
max_units_buy_30 = 4 # when inventory level is lower then 30% it is possible to buy 0 to 4 units
max_units_buy_65 = 3 # when inventory level is between 30% and 65% it is possible to buy 0 to 3 units
max_units_buy_100 = 2 # when inventory level is between 65% and 100% it is possible to buy 0 to 2 units
max_units_sell_30 = 4 # when inventory level is lower then 30% it is possible to sell 0 to 4 units
max_units_sell_70 = 6 # when inventory level is between 30% and 70% it is possible to sell 0 to 6 units
max_units_sell_100 = 8 # when inventory level is between 70% and 100% it is possible to sell 0 to 8 units

This requires a little bit of work.

Time

The first thing to consider is how time is handled in the model. In Model 1, we just used a time index \(t\) and did not really say much about it. Here, we need to be a bit more precise.

The variables buy and sell are what we call flow variables. Buying and selling takes place during time period \(t\).  Inventory, however, is measured at a specific point in time. In our case, this is at the end of period \(t\).  This is called a stock variable.

Looking at the picture we should put limits on \(\mathit{buy}_t\) and \(\mathit{sell}_t\) based on the inventory level at the end of the previous period: \(\mathit{inv}_{t-1}\).

Segmenting inventory levels

From the question, we see that we have different intervals:

• For buying, we have: 0%-30%, 30%-65% and 65%-100%
• For selling, we have: 0%-30%, 30%-70% and 70%-100%

I prefer to deal with one set of intervals, so we combine this into:

• 0%-30%, 30%-65%, 65%-70% and 70%-100%

The idea is to introduce a binary variable \(\delta_{k,t}\) indicating in which segment \(k\) we are. Obviously, we can only be in one segment, so: \[\begin{align} &\sum_k \delta_{k,t} = 1 && \forall t\\ & \delta_{k,t} \in \{0,1\}\end{align}\] So depending on which \(\delta_{k,t}\) is turned on, we have different lower- and upperbounds on the inventory levels.

This means we can write:\[\sum_k \mathit{invlb}_k \cdot \delta_{k,t}\le \mathit{inv}_{t-1} \le \sum_k \mathit{invub}_k \cdot \delta_{k,t}\] Notice that we put the bounds on the inventory level at the end of the previous period.

Limiting buying and selling

Here we need to link the values for \(\delta_{k,t}\) to bounds on buying and selling. This is very similar to what we did with respect to inventory: \[\begin{align}\mathit{buy}_t \le \sum_k \mathit{buyub}_k \cdot \delta_{k,t} && \forall t \\\mathit{sell}_t \le \sum_k \mathit{sellub}_k \cdot \delta_{k,t} && \forall t\end{align}\]

Model

We now have all the pieces to write down the complete model.

Model 2

\begin{align}\max&\sum_t \color{darkblue}{\mathit{price}}_t \cdot(\color{darkred}{\mathit{sell}}_t-\color{darkred}{\mathit{buy}}_t)\\ &\color{darkred}{\mathit{inv}}_t = \color{darkred}{\mathit{inv}}_{t-1} + \color{darkred}{\mathit{buy}}_t - \color{darkred}{\mathit{sell}}_t&& \forall t\\ &\color{darkred}{\mathit{inv}}_0 = 0\\ &\color{darkred}{\mathit{inv}}_{t-1} \ge \sum_k \color{darkblue}{\mathit{invlb}}_k \cdot \color{darkred}\delta_{k,t} && \forall t \\ &\color{darkred}{\mathit{inv}}_{t-1} \le \sum_k \color{darkblue}{\mathit{invub}}_k \cdot \color{darkred}\delta_{k,t} && \forall t \\ &\color{darkred}{\mathit{buy}}_t \le \sum_k \color{darkblue}{\mathit{buyub}}_k \cdot \color{darkred}\delta_{k,t} && \forall t \\ &\color{darkred}{\mathit{sell}}_t \le \sum_k \color{darkblue}{\mathit{sellub}}_k \cdot \color{darkred}\delta_{k,t} && \forall t \\ &\sum_k \color{darkred}\delta_{k,t} = 1 && \forall t\\ &\color{darkred}{\mathit{inv}}_t \in \{0,1,\dots,\color{darkblue}{\mathit{capacity}}\}\\ &\color{darkred}{\mathit{buy}}_t \in \{0,1,\dots\}\\ &\color{darkred}{\mathit{sell}}_t \in \{0,1,\dots\}\\ &\color{darkred}\delta_{k,t} \in \{0,1\} \end{align}

Implementation

The R implementation can look like:

> library(CVXR)
> 
> # data
> price = c(12, 11, 12, 13, 16, 17, 18, 17, 18, 16, 17, 13)
> capacity = 25
> max_units_buy = 4
> max_units_sell = 8
> 
> # capacity segments
> s <- c(0,0.3,0.65,0.7,1)
> 
> # corresponding lower and upper bounds
> invlb <- s[1:(length(s)-1)] * capacity
> invlb
[1]  0.00  7.50 16.25 17.50
> invub <- s[2:length(s)] * capacity
> invub
[1]  7.50 16.25 17.50 25.00
> 
> buyub <- c(4,3,2,2)
> sellub <- c(4,6,6,8)
> 
> # number of time periods
> NT <- length(price)
> NT
[1] 12
> 
> # number of capacity segments
> NS <- length(s)-1
> NS
[1] 4
> 
> # Decision variables
> inv = Variable(NT,integer=T)
> buy = Variable(NT,integer=T)
> sell = Variable(NT,integer=T)
> delta = Variable(NS,NT,boolean=T)
> 
> # Lag operator
> L = cbind(rbind(0,diag(NT-1)),0)
> 
> # optimization model
> problem <- Problem(Maximize(sum(price*(sell-buy))),
+                    list(inv == L %*% inv + buy - sell,
+                         sum_entries(delta,axis=2)==1, 
+                         L %*% inv >= t(delta) %*% invlb,
+                         L %*% inv <= t(delta) %*% invub,
+                         buy <= t(delta) %*% buyub,
+                         sell <= t(delta) %*% sellub,
+                         inv >= 0, inv <= capacity,
+                         buy >= 0, sell >= 0))
> result <- solve(problem,verbose=T)
GLPK Simplex Optimizer, v4.47
120 rows, 84 columns, 369 non-zeros
      0: obj =  0.000000000e+000  infeas = 1.200e+001 (24)
*    23: obj =  0.000000000e+000  infeas = 0.000e+000 (24)
*    85: obj = -9.875986758e+001  infeas = 0.000e+000 (2)
OPTIMAL SOLUTION FOUND
GLPK Integer Optimizer, v4.47
120 rows, 84 columns, 369 non-zeros
84 integer variables, 48 of which are binary
Integer optimization begins...
+    85: mip =     not found yet >=              -inf        (1; 0)
+   123: >>>>> -8.800000000e+001 >= -9.100000000e+001   3.4% (17; 0)
+   126: >>>>> -9.000000000e+001 >= -9.100000000e+001   1.1% (9; 11)
+   142: mip = -9.000000000e+001 >=     tree is empty   0.0% (0; 35)
INTEGER OPTIMAL SOLUTION FOUND
> cat("status:",result$status)
status: optimal
> cat("objective:",result$value)
objective: 90
> # print results
> data.frame(price,
+            buy=result$getValue(buy),
+            sell=result$getValue(sell),
+            inv=result$getValue(inv),
+            maxbuy=t(result$getValue(delta)) %*% buyub,
+            maxsell=t(result$getValue(delta)) %*% sellub)
   price buy sell inv maxbuy maxsell
1     12   3    0   3      4       4
2     11   4    0   7      4       4
3     12   4    0  11      4       4
4     13   3    0  14      3       6
5     16   3    0  17      3       6
6     17   1    0  18      2       6
7     18   0    8  10      2       8
8     17   0    6   4      3       6
9     18   0    4   0      4       4
10    16   4    0   4      4       4
11    17   0    4   0      4       4
12    13   0    0   0      4       4
> # display the optimal values for delta
> print(result$getValue(delta))
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
[1,]    1    1    1    0    0    0    0    0    1     1     1     1
[2,]    0    0    0    1    1    0    0    1    0     0     0     0
[3,]    0    0    0    0    0    1    0    0    0     0     0     0
[4,]    0    0    0    0    0    0    1    0    0     0     0     0
> 

The results are interesting. Compared to the first model, the profit has decreased from 140 to 90.

One question that comes to mind is: why only buy 3 in the first period? Well, if we buy 4 and 4 units in periods 1 and 2, we can no longer buy 4 units in period 3. These kinds of decisions are very difficult to make by hand: it is just too complex.

Conclusion

This is a nice little example, where some formal modeling really helps. It also shows both the strengths and weaknesses of a CVX-like modeling tool: it can make matrix- and vector-based models very easy to express. But when we stray away from this a bit and want to implement modeling concepts like lags, things become somewhat more complicated. 

References

1. R optimisation buy sell, https://stackoverflow.com/questions/61873468/r-optimisation-buy-sell
2. CVXR, Convex Optimization in R, https://cvxr.rbind.io/
3. Package Rglpk, https://cran.r-project.org/web/packages/Rglpk/Rglpk.pdf
4. R optimisation max buy/sell dependent on stock level, https://stackoverflow.com/questions/61900256/r-optimisation-max-buy-sell-dependent-on-stock-level

Elementwise vs matrix multiplication

Introduction

In this post I want to demonstrate some interesting issues in how CVXR and CVXPY deal with elementwise vs matrix multiplication and compare that to the host languages (R and Python, respectively). CVXR and CVXPY deviate a bit from what their host languages do. When implementing optimization models with these tools, it is important to be aware of these details.

Elementwise and matrix multiplication 

There are two often used methods to perform the multiplication of matrices (and vectors). The first is simply elementwise multiplication: \[c_{i,j} = a_{i,j} \cdot b_{i,j}\] In mathematics, this is sometimes referred to as the Hadamard product. The notation is typically: \[C = A \circ B\] but sometimes we see: \[C = A \odot B\] This product only works when \(A\) and \(B\) have the same shape. I.e. \[\begin{matrix} C&=&A&\circ&B\\ (m \times n) &&(m \times n)&&(m \times n)\end{matrix}\]

The standard matrix multiplication \(C=A\cdot B\) or \[c_{i,j} = \sum_k a_{i,k} b_{k,j}\] has a different rule for conformance:  \[\begin{matrix} C&=&A&\cdot&B\\ (m \times n) &&(m \times k)&&(k \times n)\end{matrix}\] Most of the time, the dot operator is dropped and we write \(C = A B\).

In optimization modeling both forms are used a lot.

R, CVXR

R has two operators for multiplication:

• * for elementwise multiplication
• %*% for matrix multiplication

A vector is considered as a column vector (i.e. an \(n\)-vector is like a \(n \times 1\) matrix).  There is one special thing in R, as shown here:

> m <- 2
> n <- 3
> A <- matrix(c(1,2,3,4,5,6),m,n)
> A
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> 
> x <- c(1,2,3)
>
 
#
# matrix multiplication
# 
> A %*% x
     [,1]
[1,]   22
[2,]   28
> 

#
# elementwise multiplication
#
> A * A
     [,1] [,2] [,3]
[1,]    1    9   25
[2,]    4   16   36

#
# but this also works
#
> A * x
     [,1] [,2] [,3]
[1,]    1    9   10
[2,]    4    4   18
>

The last multiplication is surprising as \(A\) is a \(2 \times 3\) matrix and \(x\) is different in shape. Well, R may extend and recycle vectors to make them as large as needed. In this case \(x\) is duplicated and then considered as a  \(2 \times 3\) matrix. More or less like:

> x
[1] 1 2 3
> matrix(x,2,3)
     [,1] [,2] [,3]
[1,]    1    3    2
[2,]    2    1    3

The modeling tool CVXR follows R notation and implements both * and %*% (for elementwise and matrix multiplication). However, CVXR is not implementing the extending and recycling of vectors that are too small.

Concept: recycling

Just to emphasize the concept of recycling. If an operation requires two vectors of the same length, R may make the shorter vector longer by recycling (duplicating). Here is an example:

> a <- 1:10
> b <- 1:2
> c <- a + b
> c
 [1]  2  4  4  6  6  8  8 10 10 12

In this example the vector a has elements 1 through 10. The vector b is too short, so it is recycled. When added to a, b is functionally equal to rep(c(1,2),5). When multiples of b do not exactly fit a, we effectively have a fractional duplication number. E.g. when we use b <- 1:3, we get a message:
Warning message:
In a + b : longer object length is not a multiple of shorter object length

This recycling trick is somewhat unique to R (I don't know about other languages doing this).

Example

In [2] the product: \[b_{i,j} = a_{i,j} \cdot x_i \] is implemented in R with the elementwise product:

> m <- 2
> n <- 3
> A <- matrix(c(1,2,3,4,5,6),m,n)
> A
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> x <- c(1,2)
> B <- A*x
> B
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    4    8   12

This again uses recycling of \(x\). CVXR is not doing this automatically. We can see this here:

> library(CVXR)
> x <- Variable(m)
> B <- Variable(m,n)
> e <- rep(1,n)
> problem <- Problem(Minimize(0),
+                    list(x == c(1,2), 
+                         B == A * x ))
Error in sum_shapes(lapply(object@args, function(arg) { : 
  Incompatible dimensions

Here \(x\) is now a CVXR variable. As the vector \(x\) is not recycled, we end up with two different shapes and elementwise multiplication is refused. So how do we do something like this in CVXR?

The recycling operation:

> x
[1] 1 2
> matrix(x,m,n)
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    2

can be expressed in matrix notation as: \[x \cdot e^T\] where \(e\) is a (column) vector of ones. This is sometimes called an outerproduct. I.e. we can write our assignment as \[B = A \circ (x \cdot e^T)\] In a CVXR model this can look like:

> library(CVXR)
> x <- Variable(m)
> B <- Variable(m,n)
> e <- rep(1,n)
> problem <- Problem(Minimize(0),
+                    list(x == c(1,2), 
+                         B == A * (x %*% t(e))))
> sol <- solve(problem)
> sol$status
[1] "optimal"
> sol$getValue(x)
     [,1]
[1,]    1
[2,]    2
> sol$getValue(B)
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    4    8   12

Note: the constraint B == A * (x %*% t(e)) has both a matrix multiplication and a elementwise multiplication. This is rather funky.

Conclusion: if your matrix operations rely on recycling, you will need to rework things a bit to have this work correctly in CVXR. CVXR does not do recycling.

Python, CVXPY

In the previous section we saw that there are subtle differences between R's and CVXR's elementwise multiplication semantics.

Let's now look at Python and CVXPY.

Since Python 3.5 we have two multiplication operators:

• * for elementwise multiplication 
• @ for matrix multiplication 

CVXPY has different rules:

• *, @  and matmul for matrix multiplication
• multiply for elementwise multiplication

Example

import numpy as np
import cvxpy as cp

#
# In Python/numpy * indicates elementwise multiplication
#
A = np.array([[1,2],[3,4]])
B = np.array([[1,1],[2,2]])
C = A*B
print(A)
# output:
# [[ 1 2]
#  [ 3 4]]
print(C)
# output:
# [[ 1 2]
#  [ 6 8]]

#
# In CVXPY * indicates matrix multiplication
#
A = cp.Variable((2,2))
C = cp.Variable((2,2))
prob = cp.Problem(cp.Minimize(0),
                  [A == [[1,2],[3,4]],
                   C == A*B])
prob.solve(verbose=True)
print(A.value)
# output:
# [[1. 3.]
#  [2. 4.]]
print(C.value)
# [[ 7.  7.]
#  [10. 10.]]

Here we see some differences between Python/Numpy and CVXPY. First the interpretation of Python lists (the values for A) is different. And secondly: the semantics of * are different. This may cause some confusion.

References

1. CVXR, Convex Optimization in R, https://cvxr.rbind.io/
2. CVXR Elementwise multiplication of matrix with vector, https://stackoverflow.com/questions/59224555/cvxr-elementwise-multiplication-of-matrix-with-vector