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:
- Use a large cost coefficient for forbidden links.
- Fix \(\color{darkred}x_{i,j}=0\) for non-existing links.
- 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}\] |
#------------------------------------------------- # dense transportation model # matrix formulation #------------------------------------------------- # X has same shape as our cost matrix X = cp.Variable(cost.shape,"X") # two different e vectors ei = np.ones(capacity.shape) ej = np.ones(demand.shape) prob = cp.Problem( cp.Minimize(cp.trace(cost.T@X)), [X.T@ei >= demand, X@ej <= capacity, X>=0]) prob.solve(verbose=True) print("status:",prob.status) print("objective:",prob.value)
Here @ indicates matrix multiplication. Arguably, a more readable model is using elementwise multiplication and summations:
Here cp.multiply() indicates elementwise multiplication.
#------------------------------------------------- # dense transportation model # sum formulation #------------------------------------------------- # X has same shape as our cost matrix X = cp.Variable(cost.shape,"X") prob = cp.Problem( cp.Minimize(cp.sum(cp.multiply(cost,X))), [cp.sum(X,axis=0) >= demand, cp.sum(X,axis=1) <= capacity, X >= 0]) prob.solve(verbose=True) print("status:",prob.status) print("objective:",prob.value)
Here cp.multiply() indicates elementwise multiplication.
The same model in R, using cvxr can look like:
X <- Variable(rows=n,cols=m,name='X',nonneg=T) prob <- Problem(Minimize(sum_entries(cost*X)), constraints = list( sum_entries(X,axis=1) >= dem, sum_entries(X,axis=2) <= cap)) res <- solve(prob,verbose=T) res$status res$value
Solvers
I use the default solvers for this exercise, and default settings. For CVXPY the default solver is ECOS[4], and CVXR uses OSQP[5] by default. These are not standard Simplex-based solvers.
Formulation 1: Dense Transportation Model with Large Cost Coefficients.
Here we apply the costs:
largecost = cost + (1-link)*1.0e6
where link is a 0/1 constant matrix with values: 0 if the link does not exist, and 1 if the link exists.
CVXPY uses the conic algorithm ECOS for solving LPs. This solver cannot solve the problem using default settings:
98 -4.660e+04 -4.660e+04 +5e+05 1e-04 7e-15 8e-02 2e+00 0.0232 6e-01 1 2 2 | 0 0 99 -4.655e+04 -4.655e+04 +4e+05 1e-04 9e-15 1e-01 2e+00 0.0272 7e-01 1 2 2 | 0 0 100 -4.654e+04 -4.654e+04 +4e+05 1e-04 1e-14 1e-01 2e+00 0.0407 7e-01 1 2 1 | 0 0 Maximum number of iterations reached, recovering best iterate (88) and stopping. RAN OUT OF ITERATIONS (reached feastol=1.2e-05, reltol=4.7e-02, abstol=1.1e+03). Runtime: 13.731298 seconds.
We can add a larger iteration limit. This gives:
257 +1.202e+04 +1.202e+04 +2e-01 3e-08 9e-17 4e-07 7e-07 0.9890 7e-02 1 1 1 | 0 0 258 +1.202e+04 +1.202e+04 +5e-03 1e-09 6e-17 1e-08 2e-08 0.9713 1e-04 2 1 1 | 0 0 259 +1.202e+04 +1.202e+04 +5e-05 1e-11 8e-17 1e-10 2e-10 0.9890 1e-04 2 1 1 | 0 0 OPTIMAL (within feastol=1.4e-11, reltol=4.4e-09, abstol=5.3e-05). Runtime: 34.450007 seconds. ------------------------------------------------------------------------------- Summary ------------------------------------------------------------------------------- (CVXPY) Aug 10 12:44:50 PM: Problem status: optimal (CVXPY) Aug 10 12:44:50 PM: Optimal value: 1.202e+04 (CVXPY) Aug 10 12:44:50 PM: Compilation took 5.605e+00 seconds (CVXPY) Aug 10 12:44:50 PM: Solver (including time spent in interface) took 3.446e+01 seconds status: optimal objective: 12015.21812055369
9600 2.5164e+05 5.37e-04 6.25e+00 1.01e-01 6.13e+01s 9800 2.5184e+05 5.37e-04 5.97e+00 1.01e-01 6.25e+01s 10000 2.5167e+05 5.37e-04 5.98e+00 1.01e-01 6.38e+01s status: solved inaccurate number of iterations: 10000 optimal objective: 251671.7459 run time: 6.38e+01s optimal rho estimate: 1.53e-01 Warning messages: 1: In size(const) : NAs introduced by coercion to integer range 2: In size(const) : NAs introduced by coercion to integer range
optimal objective: 94120.8210we are actually not close to the optimal solution.
This is a disappointing result.
Formulation 2: Dense Transportation Model with Fixed Variables.
The second formulation is using:
x[i,j] = 0 if link[i,j] = 0
I don't really know how to do this in a vectorized way in CVXPY/CVXR. Very strange that this is not supported, as this is an operation that we see quite often. We could form a long series of scalar equalities. I chose to formulate:
x[i,j]*(1-link[i,j]) = 0
or in matrix notation: \[\color{darkred}X \circ (\color{darkred}e\cdot\color{darkred}e^T-\color{darkblue}{\mathit link}) = 0\] where \(\circ\) indicates the Hadamar or element-wise product. In Python notation this becomes: cp.multiply(1-link,X)==0.
When running this with CVXPY/ECOS we see success!
C:\tmp>py trnsport_cvx.py CHECK: sum capacity = 10000.000000000016, sum demand = 10000.000000000002 =============================================================================== CVXPY v1.2.1 =============================================================================== (CVXPY) Aug 10 01:16:27 AM: Your problem has 250000 variables, 4 constraints, and 0 parameters. (CVXPY) Aug 10 01:16:27 AM: It is compliant with the following grammars: DCP, DQCP (CVXPY) Aug 10 01:16:27 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.) (CVXPY) Aug 10 01:16:27 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution. ------------------------------------------------------------------------------- Compilation ------------------------------------------------------------------------------- (CVXPY) Aug 10 01:16:27 AM: Compiling problem (target solver=ECOS). (CVXPY) Aug 10 01:16:27 AM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> ECOS (CVXPY) Aug 10 01:16:27 AM: Applying reduction Dcp2Cone (CVXPY) Aug 10 01:16:27 AM: Applying reduction CvxAttr2Constr (CVXPY) Aug 10 01:16:27 AM: Applying reduction ConeMatrixStuffing (CVXPY) Aug 10 01:16:27 AM: Applying reduction ECOS (CVXPY) Aug 10 01:16:27 AM: Finished problem compilation (took 4.100e-01 seconds). ------------------------------------------------------------------------------- Numerical solver ------------------------------------------------------------------------------- (CVXPY) Aug 10 01:16:27 AM: Invoking solver ECOS to obtain a solution. ECOS 2.0.10 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS It pcost dcost gap pres dres k/t mu step sigma IR | BT 0 +5.471e+04 +5.471e+04 +2e+06 5e-01 5e-01 1e+00 1e+01 --- --- 1 1 - | - - 1 +3.188e+04 +3.188e+04 +1e+06 3e-01 2e-01 2e+00 5e+00 0.6991 3e-01 1 1 1 | 0 0 2 +2.670e+04 +2.670e+04 +1e+06 2e-01 2e-01 3e+00 4e+00 0.3511 6e-01 1 1 0 | 0 0 3 +1.558e+04 +1.558e+04 +4e+05 8e-02 6e-02 2e+00 2e+00 0.9518 3e-01 1 1 1 | 0 0 4 +1.381e+04 +1.381e+04 +2e+05 3e-02 2e-02 8e-01 7e-01 0.7207 2e-01 1 1 0 | 0 0 5 +1.318e+04 +1.318e+04 +1e+05 2e-02 2e-02 4e-01 4e-01 0.6765 5e-01 1 1 0 | 0 0 6 +1.261e+04 +1.261e+04 +5e+04 8e-03 7e-03 2e-01 2e-01 0.6174 1e-01 1 1 1 | 0 0 7 +1.235e+04 +1.235e+04 +3e+04 4e-03 3e-03 9e-02 1e-01 0.6661 3e-01 1 1 0 | 0 0 8 +1.218e+04 +1.218e+04 +1e+04 2e-03 1e-03 3e-02 5e-02 0.7733 3e-01 1 1 1 | 0 0 9 +1.211e+04 +1.211e+04 +6e+03 8e-04 7e-04 2e-02 2e-02 0.7199 3e-01 1 1 1 | 0 0 10 +1.206e+04 +1.206e+04 +3e+03 4e-04 4e-04 8e-03 1e-02 0.7449 3e-01 1 1 1 | 0 0 11 +1.204e+04 +1.204e+04 +1e+03 2e-04 2e-04 3e-03 5e-03 0.6062 9e-02 1 1 1 | 0 0 12 +1.202e+04 +1.202e+04 +5e+02 7e-05 7e-05 1e-03 2e-03 0.6782 1e-01 1 1 1 | 0 0 13 +1.202e+04 +1.202e+04 +2e+02 2e-05 2e-05 4e-04 6e-04 0.7878 1e-01 1 1 1 | 0 0 14 +1.202e+04 +1.202e+04 +9e+01 1e-05 1e-05 2e-04 4e-04 0.4928 2e-01 2 1 1 | 0 0 15 +1.202e+04 +1.202e+04 +3e+01 3e-06 3e-06 6e-05 1e-04 0.8573 2e-01 2 1 1 | 0 0 16 +1.202e+04 +1.202e+04 +2e+01 2e-06 2e-06 4e-05 6e-05 0.6635 4e-01 3 1 2 | 0 0 17 +1.202e+04 +1.202e+04 +2e+00 2e-07 2e-07 5e-06 7e-06 0.8947 2e-02 4 1 2 | 0 0 18 +1.202e+04 +1.202e+04 +2e-01 3e-08 2e-08 5e-07 8e-07 0.9343 4e-02 6 2 3 | 0 0 19 +1.202e+04 +1.202e+04 +2e-03 3e-10 3e-10 6e-09 9e-09 0.9890 1e-03 1 1 1 | 0 0 20 +1.202e+04 +1.202e+04 +3e-05 3e-12 3e-12 6e-11 1e-10 0.9890 1e-04 2 0 0 | 0 0 OPTIMAL (within feastol=3.4e-12, reltol=2.2e-09, abstol=2.6e-05). Runtime: 4.309761 seconds. ------------------------------------------------------------------------------- Summary ------------------------------------------------------------------------------- (CVXPY) Aug 10 01:16:31 AM: Problem status: optimal (CVXPY) Aug 10 01:16:31 AM: Optimal value: 1.202e+04 (CVXPY) Aug 10 01:16:31 AM: Compilation took 4.100e-01 seconds (CVXPY) Aug 10 01:16:31 AM: Solver (including time spent in interface) took 4.322e+00 seconds status: optimal objective: 12015.21812552263 C:\tmp>
The first thing we notice is that counting is difficult: Your problem has 250000 variables, 4 constraints. We seem to count here individual variables and constraint blocks. I guarantee you there are more than 4 constraints. That sounds strange. I mean, why in the world would you do that? Choose the same units!
Now let's see if we also have some success with the CVXR/OSQP setup. Unfortunately, this still does not solve:
-----------------------------------------------------------------
OSQP v0.6.0 - Operator Splitting QP Solver
(c) Bartolomeo Stellato, Goran Banjac
University of Oxford - Stanford University 2019
-----------------------------------------------------------------
problem: variables n = 250000, constraints m = 501000
nnz(P) + nnz(A) = 950334
settings: linear system solver = qdldl,
eps_abs = 1.0e-05, eps_rel = 1.0e-05,
eps_prim_inf = 1.0e-04, eps_dual_inf = 1.0e-04,
rho = 1.00e-01 (adaptive),
sigma = 1.00e-06, alpha = 1.60, max_iter = 10000
check_termination: on (interval 25),
scaling: on, scaled_termination: off
warm start: on, polish: on, time_limit: off
iter objective pri res dua res rho time
1 -5.4672e+05 4.33e+01 4.22e+01 1.00e-01 8.30e-01s
200 1.2118e+04 5.96e-02 1.12e-01 1.00e-01 3.08e+00s
400 1.2047e+04 1.95e-02 9.76e-02 1.00e-01 5.46e+00s
600 1.2031e+04 7.57e-03 3.55e-02 1.00e-01 7.73e+00s
800 1.2043e+04 4.77e-03 2.27e-02 1.00e-01 9.95e+00s
1000 1.2023e+04 1.63e-02 1.40e-02 1.87e-02 1.23e+01s
1200 1.2055e+04 7.39e-03 8.88e-03 1.87e-02 1.45e+01s
1400 1.2034e+04 6.94e-03 7.34e-03 1.87e-02 1.67e+01s
1600 1.2048e+04 3.92e-03 4.51e-03 1.87e-02 1.89e+01s
1800 1.2047e+04 3.34e-03 4.63e-03 1.87e-02 2.12e+01s
2000 1.2038e+04 3.35e-03 2.90e-03 1.87e-02 2.34e+01s
2200 1.2041e+04 2.34e-03 2.50e-03 1.87e-02 2.56e+01s
2400 1.2040e+04 2.32e-03 2.47e-03 1.87e-02 2.78e+01s
2600 1.2046e+04 1.75e-03 1.82e-03 1.87e-02 3.00e+01s
2800 1.2044e+04 1.86e-03 1.45e-03 1.87e-02 3.22e+01s
3000 1.2044e+04 1.39e-03 1.47e-03 1.87e-02 3.45e+01s
3200 1.2041e+04 1.33e-03 1.34e-03 1.87e-02 3.67e+01s
3400 1.2045e+04 1.02e-03 1.62e-03 1.87e-02 3.89e+01s
3600 1.2043e+04 9.00e-04 1.20e-03 1.87e-02 4.12e+01s
3800 1.2043e+04 1.04e-03 1.12e-03 1.87e-02 4.34e+01s
4000 1.2040e+04 1.32e-03 1.17e-03 1.87e-02 4.56e+01s
4200 1.2044e+04 7.46e-04 9.03e-04 1.87e-02 4.78e+01s
4400 1.2043e+04 6.78e-04 8.80e-04 1.87e-02 5.00e+01s
4600 1.2042e+04 6.31e-04 1.19e-03 1.87e-02 5.22e+01s
4800 1.2042e+04 7.80e-04 8.06e-04 1.87e-02 5.44e+01s
5000 1.2042e+04 6.57e-04 8.05e-04 1.87e-02 5.67e+01s
5200 1.2044e+04 6.17e-04 9.06e-04 1.87e-02 5.89e+01s
5400 1.2042e+04 6.12e-04 6.06e-04 1.87e-02 6.11e+01s
5600 1.2041e+04 6.23e-04 6.60e-04 1.87e-02 6.33e+01s
5800 1.2041e+04 3.80e-04 7.29e-04 1.87e-02 6.55e+01s
6000 1.2043e+04 5.39e-04 6.00e-04 1.87e-02 6.78e+01s
6200 1.2043e+04 4.63e-04 6.70e-04 1.87e-02 7.00e+01s
6400 1.2042e+04 4.04e-04 6.56e-04 1.87e-02 7.22e+01s
6600 1.2041e+04 5.45e-04 4.48e-04 1.87e-02 7.44e+01s
6800 1.2042e+04 6.47e-04 5.94e-04 1.87e-02 7.66e+01s
7000 1.2043e+04 3.81e-04 5.65e-04 1.87e-02 7.88e+01s
7200 1.2043e+04 7.66e-04 4.04e-04 1.87e-02 8.11e+01s
7400 1.2042e+04 2.75e-04 5.74e-04 1.87e-02 8.33e+01s
7600 1.2043e+04 5.49e-04 4.90e-04 1.87e-02 8.55e+01s
7800 1.2042e+04 3.88e-04 3.64e-04 1.87e-02 8.77e+01s
8000 1.2042e+04 3.03e-04 5.45e-04 1.87e-02 8.99e+01s
8200 1.2043e+04 3.03e-04 5.13e-04 1.87e-02 9.21e+01s
8400 1.2043e+04 3.23e-04 3.85e-04 1.87e-02 9.43e+01s
8600 1.2042e+04 3.14e-04 4.93e-04 1.87e-02 9.65e+01s
8800 1.2043e+04 2.98e-04 3.68e-04 1.87e-02 9.87e+01s
9000 1.2042e+04 3.10e-04 2.46e-04 1.87e-02 1.01e+02s
9200 1.2043e+04 2.24e-04 4.20e-04 1.87e-02 1.03e+02s
9400 1.2042e+04 2.69e-04 3.78e-04 1.87e-02 1.05e+02s
9600 1.2042e+04 2.80e-04 2.56e-04 1.87e-02 1.08e+02s
9800 1.2043e+04 2.09e-04 4.05e-04 1.87e-02 1.10e+02s
10000 1.2042e+04 2.78e-04 3.05e-04 1.87e-02 1.12e+02s
status: solved inaccurate
number of iterations: 10000
optimal objective: 12042.4795
run time: 1.12e+02s
optimal rho estimate: 9.00e-03
Warning messages: 1: In size(const) : NAs introduced by coercion to integer range 2: In size(const) : NAs introduced by coercion to integer range 3: In size(const) : NAs introduced by coercion to integer range
Formulation 3: exploit sparsity
This is a more difficult exercise. CVXPY and CVXR do not have a conditional generation facility. Or in other words, we cannot generate a sparse variable. In practical modeling, this is rather essential. So how can we generate only variables that are related to existing links? Well, we can do everything ourselves and basically do all the work that CVXPY/CVXR is supposed to do. We can also fudge things a bit. We still generate all variables, but we multiply each variable in each constraint by 0 if it should not exist. This is only a partial solution, as we still generate all variables. However, a variable with all zeros in its column should be removed by the solver. We hope that is happening, but we have little feedback that is actually happening.
In CVXPY I did:
#---------------------------------------------------- # model 3 # exploit sparsity #---------------------------------------------------- # X has same shape as our cost matrix X = cp.Variable(cost.shape,"X") prob = cp.Problem( cp.Minimize(cp.sum(cp.multiply(cp.multiply(cost,plink),X))), [cp.sum(cp.multiply(plink,X),axis=0) >= demand, cp.sum(cp.multiply(plink,X),axis=1) <= capacity, cp.multiply(plink,X)>=0] ) prob.solve(verbose=True) print("status:",prob.status) print("objective:",prob.value)
This actually seems to help a bit:
C:\tmp>py trnsport_cvx.py CHECK: sum capacity = 10000.000000000016, sum demand = 10000.000000000002 =============================================================================== CVXPY v1.2.1 =============================================================================== (CVXPY) Aug 10 01:57:25 AM: Your problem has 250000 variables, 3 constraints, and 0 parameters. (CVXPY) Aug 10 01:57:25 AM: It is compliant with the following grammars: DCP, DQCP (CVXPY) Aug 10 01:57:25 AM: (If you need to solve this problem multiple times, but with different data, consider using parameters.) (CVXPY) Aug 10 01:57:25 AM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution. ------------------------------------------------------------------------------- Compilation ------------------------------------------------------------------------------- (CVXPY) Aug 10 01:57:25 AM: Compiling problem (target solver=ECOS). (CVXPY) Aug 10 01:57:25 AM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> ECOS (CVXPY) Aug 10 01:57:25 AM: Applying reduction Dcp2Cone (CVXPY) Aug 10 01:57:25 AM: Applying reduction CvxAttr2Constr (CVXPY) Aug 10 01:57:25 AM: Applying reduction ConeMatrixStuffing (CVXPY) Aug 10 01:57:25 AM: Applying reduction ECOS (CVXPY) Aug 10 01:57:25 AM: Finished problem compilation (took 1.830e-01 seconds). ------------------------------------------------------------------------------- Numerical solver ------------------------------------------------------------------------------- (CVXPY) Aug 10 01:57:25 AM: Invoking solver ECOS to obtain a solution. ECOS 2.0.10 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS It pcost dcost gap pres dres k/t mu step sigma IR | BT 0 +5.471e+04 +5.471e+04 +2e+06 5e-01 3e-01 1e+00 1e+01 --- --- 1 1 - | - - 1 +3.188e+04 +3.188e+04 +1e+06 3e-01 1e-01 2e+00 5e+00 0.6991 3e-01 1 1 1 | 0 0 2 +2.670e+04 +2.670e+04 +1e+06 2e-01 1e-01 3e+00 4e+00 0.3511 6e-01 1 1 0 | 0 0 3 +1.558e+04 +1.558e+04 +4e+05 8e-02 4e-02 2e+00 2e+00 0.9519 3e-01 1 1 0 | 0 0 4 +1.381e+04 +1.381e+04 +2e+05 3e-02 2e-02 8e-01 7e-01 0.7207 2e-01 1 0 0 | 0 0 5 +1.318e+04 +1.318e+04 +1e+05 2e-02 1e-02 4e-01 4e-01 0.6765 5e-01 1 0 0 | 0 0 6 +1.261e+04 +1.261e+04 +5e+04 8e-03 5e-03 2e-01 2e-01 0.6174 1e-01 1 0 0 | 0 0 7 +1.235e+04 +1.235e+04 +3e+04 4e-03 3e-03 9e-02 1e-01 0.6661 3e-01 1 0 0 | 0 0 8 +1.218e+04 +1.218e+04 +1e+04 2e-03 1e-03 3e-02 5e-02 0.7733 3e-01 1 0 1 | 0 0 9 +1.211e+04 +1.211e+04 +6e+03 8e-04 6e-04 2e-02 2e-02 0.7198 3e-01 1 1 1 | 0 0 10 +1.206e+04 +1.206e+04 +3e+03 4e-04 3e-04 8e-03 1e-02 0.7449 3e-01 1 1 1 | 0 0 11 +1.204e+04 +1.204e+04 +1e+03 2e-04 1e-04 3e-03 5e-03 0.6063 9e-02 1 1 1 | 0 0 12 +1.202e+04 +1.202e+04 +5e+02 7e-05 5e-05 1e-03 2e-03 0.6782 1e-01 1 1 1 | 0 0 13 +1.202e+04 +1.202e+04 +2e+02 2e-05 2e-05 4e-04 6e-04 0.7878 1e-01 1 1 1 | 0 0 14 +1.202e+04 +1.202e+04 +9e+01 1e-05 9e-06 2e-04 4e-04 0.4929 2e-01 1 1 1 | 0 0 15 +1.202e+04 +1.202e+04 +3e+01 3e-06 3e-06 6e-05 1e-04 0.8568 2e-01 1 1 1 | 0 0 16 +1.202e+04 +1.202e+04 +2e+01 2e-06 2e-06 4e-05 6e-05 0.6638 4e-01 1 1 1 | 0 0 17 +1.202e+04 +1.202e+04 +2e+00 2e-07 2e-07 5e-06 7e-06 0.8945 2e-02 1 1 1 | 0 0 18 +1.202e+04 +1.202e+04 +2e-01 3e-08 2e-08 5e-07 8e-07 0.9341 4e-02 2 1 1 | 0 0 19 +1.202e+04 +1.202e+04 +2e-03 3e-10 2e-10 6e-09 9e-09 0.9890 1e-03 3 0 0 | 0 0 20 +1.202e+04 +1.202e+04 +3e-05 3e-12 3e-12 6e-11 1e-10 0.9890 1e-04 5 0 0 | 0 0 OPTIMAL (within feastol=3.4e-12, reltol=2.2e-09, abstol=2.6e-05). Runtime: 2.398611 seconds. ------------------------------------------------------------------------------- Summary ------------------------------------------------------------------------------- (CVXPY) Aug 10 01:57:27 AM: Problem status: optimal (CVXPY) Aug 10 01:57:27 AM: Optimal value: 1.202e+04 (CVXPY) Aug 10 01:57:27 AM: Compilation took 1.830e-01 seconds (CVXPY) Aug 10 01:57:27 AM: Solver (including time spent in interface) took 2.408e+00 seconds status: optimal objective: 12015.218125523525 C:\tmp>
I also tried this with CVXR/OSQP, but the result was still a non-optimal solution: solved inaccurate. After increasing the iteration limit the status became solved. However, the objective value was not really close to the optimal objective.
Conclusion
The conclusion is not encouraging for OSQP. That solver has difficulties finding the optimal solution.
| CVXPY/ECOS | CVXR/OSQP | |||||
|---|---|---|---|---|---|---|
| Large Cost | Fixed Variables | Sparse Formulation | Large Cost | Fixed Variables | Sparse Formulation | |
| Variables | 250,000 | 250,000 | 250,000 | 250,000 | 250,000 | 250,000 |
| Constraints | na | na | na | na | na | na |
| Status | Optimal | Optimal | Optimal | solved | solved | solved |
| Objective | 12015.218 | 12015.218 | 12015.218 | 94120.8210 | 12042.4629 | 12042.4477 |
| Compilation Time | 5.6 | 0.4 | 0.1 | |||
| Solver Time | 34.5 | 4.4 | 2.3 | 69.3 | 169 | 93.8 |
| Iterations | 259 | 20 | 20 | 10775 | 14450 | 16100 |
CVXPY and CVXR do not allow variables with more than 2 dimensions. Also, it does not support sparse variables. ECOS delivers correct optimal solutions, but OSQP gives suboptimal solutions when solving this very large but very easy LP. Note that this is with default tolerances eps_rel and eps_abs. We can tighten these to get better solutions, but at a considerable computational cost. OSQP is a bit of a special solver based on a so-called first-order method (allowing relatively large problems to be approximated quickly using limited amounts of memory). It may not be the most suited as a default solver.
It is interesting to see that our haphazard sparse formulation works in the sense that it gives the correct solution, but in less time (compilation time and solver time).
The sparse formulation in CVXPY/ECOS is about 10 times as slow as GAMS/CPLEX [1].
References
- Transportation model with some non-existing links, https://yetanothermathprogrammingconsultant.blogspot.com/2022/07/transportation-model-with-some-non.html
- Welcome to CVXPY 1.2, https://www.cvxpy.org/
- Convex Optimization in R, https://cvxr.rbind.io/
- Alexander Domahidi, Eric Chu, Stephen Boyd, ECOS: An SOCP solver for embedded systems, Proceedings European Control Conference, pages 3071-3076, Zurich, July 2013.
- OSQP, https://osqp.org/
Appendix: Python code
import numpy as np import pickle import pandas as pd import cvxpy as cp # # get the GAMS data # stored in a pickle file # datafile = r"c:\tmp\assignment.pkl" data = pickle.load(open(datafile,'rb')) capacity = data['capacity'] demand = data['demand'] cost = data['cost'] plink = data['plink'] # # we should have a balanced problem # print(f"CHECK: sum capacity = {sum(capacity)}, sum demand = {sum(demand)}") #------------------------------------------------- # problem 1: large cost #------------------------------------------------- largecost = cost+(1-plink)*1.0e6 #------------------------------------------------- # dense transportation model # matrix formulation #------------------------------------------------- # X has same shape as our cost matrix X = cp.Variable(cost.shape,"X") # two different e vectors ei = np.ones(capacity.shape) ej = np.ones(demand.shape) prob = cp.Problem( cp.Minimize(cp.trace(largecost.T@X)), [X.T@ei >= demand, X@ej <= capacity, X>=0]) prob.solve(verbose=True,max_iters=1000) print("status:",prob.status) print("objective:",prob.value) #------------------------------------------------- # dense transportation model # sum formulation #------------------------------------------------- # X has same shape as our cost matrix X = cp.Variable(cost.shape,"X",nonneg=True) prob = cp.Problem( cp.Minimize(cp.sum(cp.multiply(largecost,X))), [cp.sum(X,axis=0) >= demand, cp.sum(X,axis=1) <= capacity]) prob.solve(verbose=True,max_iters=1000) print("status:",prob.status) print("objective:",prob.value) #---------------------------------------------------- # model 2 # fix variables #---------------------------------------------------- # X has same shape as our cost matrix X = cp.Variable(cost.shape,"X",nonneg=True) prob = cp.Problem( cp.Minimize(cp.sum(cp.multiply(cost,X))), [cp.sum(X,axis=0) >= demand, cp.sum(X,axis=1) <= capacity, cp.multiply(1-plink,X)==0]) prob.solve(verbose=True) print("status:",prob.status) print("objective:",prob.value) #---------------------------------------------------- # model 3 # exploit sparsity (somewhat) #---------------------------------------------------- # X has same shape as our cost matrix X = cp.Variable(cost.shape,"X") prob = cp.Problem( cp.Minimize(cp.sum(cp.multiply(cp.multiply(cost,plink),X))), [cp.sum(cp.multiply(plink,X),axis=0) >= demand, cp.sum(cp.multiply(plink,X),axis=1) <= capacity, cp.multiply(plink,X)>=0] ) prob.solve(verbose=True) print("status:",prob.status) print("objective:",prob.value)
Appendix: R code
#----------------------------------------------------------- # libraries #----------------------------------------------------------- library(CVXR) #----------------------------------------------------------- # read GAMS data # capacity, demand, c, plink #----------------------------------------------------------- load(file="c:/tmp/gamsdata.Rdata") #----------------------------------------------------------- # convert from dataframes to matrices/vectors #----------------------------------------------------------- convSup <- 1:nrow(capacity) names(convSup) <- capacity$i convDem <- 1:nrow(demand) names(convDem) <- demand$j cost<-matrix(data=0,nrow=length(convSup),ncol=length(convDem)) # loop is fast enough for (r in 1:nrow(c)) { cost[convSup[c$i[r]],convDem[c$j[r]]] = c$value[r] } islink<-matrix(data=0,nrow=length(convSup),ncol=length(convDem)) # loop is fast enough for (r in 1:nrow(plink)) { islink[convSup[plink$i[r]],convDem[plink$j[r]]] = plink$value[r] } cap <- capacity$value dem <- demand$value paste(sum(cap),sum(dem)) n <- length(cap) m <- length(dem) #----------------------------------------------------------- # Model 1: large cost #----------------------------------------------------------- largecost = cost + (1-islink)*1e6 X <- Variable(rows=n,cols=m,name='X',nonneg=T) prob <- Problem(Minimize(sum_entries(largecost*X)), constraints = list( sum_entries(X,axis=1) >= dem, sum_entries(X,axis=2) <= cap)) res <- solve(prob,verbose=T,max_iter=1000000) res$status res$value #----------------------------------------------------------- # Model 2: Fixed Variables #----------------------------------------------------------- X <- Variable(rows=n,cols=m,name='X',nonneg=T) prob <- Problem(Minimize(sum_entries(cost*X)), constraints = list( sum_entries(X,axis=1) >= dem, sum_entries(X,axis=2) <= cap, X*(1-islink)==0)) res <- solve(prob,verbose=T,max_iter=1000000) res$status res$value #----------------------------------------------------------- # Model 3: Exploit sparsity #----------------------------------------------------------- X <- Variable(rows=n,cols=m,name='X') prob <- Problem(Minimize(sum_entries(cost*islink*X)), constraints = list( sum_entries(islink*X,axis=1) >= dem, sum_entries(islink*X,axis=2) <= cap, islink*X>=0)) res <- solve(prob,verbose=T,max_iter=1000000) res$status res$value
No comments:
Post a Comment