## Overview

The goal of this notebook is to illustrate and compare different LCP solutions to Ben Moll's [vanilla example](http://www.princeton.edu/~moll/HACTproject/option_simple.pdf).

Here is a summary of the methods tried and relevant information.

| Method | Successful? | Time | 
|--------|----------|--------|
| `PATHSolver.jl` | Yes | 4e-6 seconds (`1e-12` tolerance, second run) |
| `Complementarity.jl` | No (JuMP) | N/A (Note: Calls to `PATHSolver` for solution) | 
| `NLsolve.jl` | Yes | 60e-6 seconds (`ftol = 1e-12`, second run) | 
| `VariationalInequality.jl` | No (JuMP) | N/A | 
| `JuMP` + Ipopt | Yes | 55e-6 seconds | 

#### Setup

Let's pull parameters from the [MATLAB code](https://github.com/econtoolkit/continuous_time_methods/blob/master/matlab/tests/simple_optimal_stopping_diffusion_test.m).


```julia
using Parameters, LinearAlgebra

StoppingProblem = @with_kw (μ_bar = -0.01, 
                            σ_bar = 0.01, 
                            S_bar = 10.0, 
                            γ = 0.5, # u(x) = x^γ
                            ρ = 0.05, # discount rate
                            x_min = 0.01,
                            x_max = 1.0,
                            M = 1000) # num of grid points 
```




    #3 (generic function with 2 methods)



Let's also copy the discretization method from [PTW](https://github.com/jlperla/PerlaTonettiWaugh.jl/commits/master/src/utils/diffusionoperators.jl). This is for a reflecting barrier at each boundary.


```julia
# Diffusion operators with regular grids
function diffusionoperators(x::AbstractRange; ξ = 0.0) # ξ = 0.0 corresponds to the unrescaled case from PTW
    Δ = step(x)
    P = length(x)

    dl_1 = zeros(P-1)
    d_1 = -ones(P)
    du_1 = ones(P-1)
    d_1[end] = d_1[end] + du_1[end] * (1-ξ*Δ)
    L_1_plus = Tridiagonal(dl_1, d_1, du_1)/Δ

    dl_m1 = -ones(P-1)
    d_m1 = ones(P)
    du_m1 = zeros(P-1)
    d_m1[1] = d_m1[1] + dl_m1[1] * (1+ξ*Δ)
    L_1_minus = Tridiagonal(dl_m1, d_m1, du_m1)/Δ

    dl_2 = ones(P-1)
    d_2 = -2 * ones(P)
    d_2[1] = -2 + (1+ξ*Δ)
    d_2[end] = -2 + (1-ξ*Δ)
    du_2 = ones(P-1)
    L_2 = Tridiagonal(dl_2, d_2, du_2)/(Δ^2)

    return (x, L_1_minus, L_1_plus, L_2)
end
```




    diffusionoperators (generic function with 1 method)



### Approach 1: Simplest LCP with PATHSolver.jl
-------

Here, we consider the case with constrant drift $\mu(x) = \bar{\mu}$, stopping value $S(x) \equiv \bar{S}$, and variance $\sigma(x) = \bar{\sigma}$.

First, here's a function to take model parameters and return LCP objects.


```julia
function LCP_objects(sp)
    # setup work 
    @unpack μ_bar, σ_bar, S_bar, γ, ρ, x_min, x_max, M = sp 
    grid = range(x_min, x_max, length = M)
    tmp, L_1_minus, L_1_plus, L_2 = diffusionoperators(grid) 
    μ = μ_bar 
    S = x -> S_bar 
    u = x -> x^γ # L36 from MATLAB
    σ = σ_bar 
    # construct objects
    A = μ * L_1_minus + 1/2 * σ^2 * L_2 # stochastic generator. use L_1_minus because drift is always negative
    B = ρ*I - A 
    S_vec = S.(grid)
    -u.(grid)
    q_vec = -u.(grid) + B*S_vec
    # return 
    return (A = A, B = B, S = S_vec, q = q_vec)
end
```




    LCP_objects (generic function with 1 method)



Let's generate these objects for the default case.


```julia
sp = StoppingProblem()
@unpack A, B, S, q = LCP_objects(sp);
```

This suggests the following linear complementarity problem:

$$ 
\begin{align}
z'(\mathbf{B}z + q) &= 0 \\ 
z &\geq 0 \\ 
\mathbf{B}z + q &\geq 0
\end{align}
$$

Let's give it to PATHSolver.


```julia
using PATHSolver
```


```julia
f = z -> B*z + q
n = sp.M
lb = zeros(n)
ub = 1000*ones(n) # a reasonable guess?
options(convergence_tolerance = 1e-12, output = :no, time_limit = 600) # 10 minute budget
@time exit_code, sol_z, sol_f = solveLCP(f, lb, ub)
```

      1.377233 seconds (3.40 M allocations: 197.704 MiB, 4.10% gc time)
    Reading options file path.opt
     > convergence_tolerance 1.0e-12
    Read of options file complete.
    
    1000 row/cols, 2998 non-zeros, 0.30% dense.
    
    Path 4.7.03 (Thu Jan 24 15:44:12 2013)
    Written by Todd Munson, Steven Dirkse, and Michael Ferris
    
     ** EXIT - other error.





    (:StationaryPointFound, [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0  …  7.6867, 7.69559, 7.70406, 7.71202, 7.71937, 7.72598, 7.73171, 7.73638, 7.73977, 7.74162], [0.4, 0.395162, 0.390538, 0.386101, 0.381831, 0.37771, 0.373723, 0.369858, 0.366105, 0.362454  …  -3.58602e-14, 1.84297e-14, -7.07212e-14, 2.39808e-14, -3.58602e-14, 5.09592e-14, 8.43769e-14, -1.4877e-14, 7.63833e-14, 5.68434e-14])




```julia
exit_code # this lets us handle logic like converged(exit_code) || error("foobar")
```




    :StationaryPointFound



Now, we can unpack the solutions and plot them.


```julia
using Plots
v = sol_z + S 
plot(v, title = "Value Function for Vanilla Model", legend = false, xlabel = "x", ylabel = "v(x)")
```




    

    



This is the kind of thing we expect --- below some threshold you get the stopping value $S \equiv 10$, and above it you get a smoothly increasing value.

**Note:** `PATHSolver.jl` is the backend used by `Complementarity.jl`, so there's no need to test that package separately (and, in fact, I couldn't get the macro/JuMP setup to work.)

### Approach 2: Same Problem with NLsolve
------


```julia
using NLsolve
```

    ┌ Info: Recompiling stale cache file /Users/Jozefiak/.julia/compiled/v1.1/NLsolve/KFCNP.ji for NLsolve [2774e3e8-f4cf-5e23-947b-6d7e65073b56]
    └ @ Base loading.jl:1184



```julia
function f!(F, z)
    F .= B*z + q # in place assignment
end
```




    f! (generic function with 1 method)




```julia
@time r = mcpsolve(f!, lb, ub,
                zeros(n), # initial condition
                reformulation = :smooth, # uses a so-called "Fischer function" to smooth out the problem
                autodiff = :forward,
                ftol = 1e-12);
```

      5.725109 seconds (12.87 M allocations: 1.253 GiB, 9.23% gc time)


As before, we can plot the resulting value function.


```julia
v = r.zero + S 
plot(v, title = "Value Function for Vanilla Model", legend = false, xlabel = "x", ylabel = "v(x)")
```




    

    



Reassuringly, they look about the same.

### Approach 3: JuMP
--------


```julia
using JuMP, Ipopt
```


```julia
#m = Model(solver=IpoptSolver(tol=1e-12))
m = Model(with_optimizer(Ipopt.Optimizer, tol=1e-12))

```




$$ \begin{alignat*}{1}\text{feasibility}\\
\text{Subject to} \quad\end{alignat*}
 $$




```julia
@variable(m, z[1:sp.M])
```




    1000-element Array{VariableRef,1}:
     z[1]   
     z[2]   
     z[3]   
     z[4]   
     z[5]   
     z[6]   
     z[7]   
     z[8]   
     z[9]   
     z[10]  
     z[11]  
     z[12]  
     z[13]  
     ⋮      
     z[989] 
     z[990] 
     z[991] 
     z[992] 
     z[993] 
     z[994] 
     z[995] 
     z[996] 
     z[997] 
     z[998] 
     z[999] 
     z[1000]




```julia
@constraint(m, z .>= 0) # first nonnegativity constraint
```




    1000-element Array{ConstraintRef{Model,MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}},ScalarShape},1}:
     z[1] ≥ 0.0   
     z[2] ≥ 0.0   
     z[3] ≥ 0.0   
     z[4] ≥ 0.0   
     z[5] ≥ 0.0   
     z[6] ≥ 0.0   
     z[7] ≥ 0.0   
     z[8] ≥ 0.0   
     z[9] ≥ 0.0   
     z[10] ≥ 0.0  
     z[11] ≥ 0.0  
     z[12] ≥ 0.0  
     z[13] ≥ 0.0  
     ⋮            
     z[989] ≥ 0.0 
     z[990] ≥ 0.0 
     z[991] ≥ 0.0 
     z[992] ≥ 0.0 
     z[993] ≥ 0.0 
     z[994] ≥ 0.0 
     z[995] ≥ 0.0 
     z[996] ≥ 0.0 
     z[997] ≥ 0.0 
     z[998] ≥ 0.0 
     z[999] ≥ 0.0 
     z[1000] ≥ 0.0




```julia
@constraint(m, B*z + q .>= 0) # second nonlinearity complaint
```




    1000-element Array{ConstraintRef{Model,MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}},ScalarShape},1}:
     50.963223140495856 z[1] - 50.91322314049586 z[2] ≥ -0.4                                                
     -61.00413223140495 z[1] + 111.96735537190081 z[2] - 50.91322314049586 z[3] ≥ -0.3951620727457512       
     -61.00413223140495 z[2] + 111.96735537190081 z[3] - 50.91322314049586 z[4] ≥ -0.39053776001763396      
     -61.00413223140495 z[3] + 111.96735537190081 z[4] - 50.91322314049586 z[5] ≥ -0.3861010405098138       
     -61.00413223140495 z[4] + 111.96735537190081 z[5] - 50.91322314049586 z[6] ≥ -0.3818307825026454       
     -61.00413223140495 z[5] + 111.96735537190081 z[6] - 50.91322314049586 z[7] ≥ -0.3777095467547562       
     -61.00413223140495 z[6] + 111.96735537190081 z[7] - 50.91322314049586 z[8] ≥ -0.37372274177066866      
     -61.00413223140495 z[7] + 111.96735537190081 z[8] - 50.91322314049586 z[9] ≥ -0.3698580123983526       
     -61.00413223140495 z[8] + 111.96735537190081 z[9] - 50.91322314049586 z[10] ≥ -0.3661047875094203      
     -61.00413223140495 z[9] + 111.96735537190081 z[10] - 50.91322314049586 z[11] ≥ -0.3624539389190383     
     -61.00413223140495 z[10] + 111.96735537190081 z[11] - 50.91322314049586 z[12] ≥ -0.3588975198308767    
     -61.00413223140495 z[11] + 111.96735537190081 z[12] - 50.91322314049586 z[13] ≥ -0.35542856125476224   
     -61.00413223140495 z[12] + 111.96735537190081 z[13] - 50.91322314049586 z[14] ≥ -0.35204091142529215   
     ⋮                                                                                                      
     -61.00413223140495 z[988] + 111.96735537190081 z[989] - 50.91322314049586 z[990] ≥ 0.494534614329171   
     -61.00413223140495 z[989] + 111.96735537190081 z[990] - 50.91322314049586 z[991] ≥ 0.495032708050275   
     -61.00413223140495 z[990] + 111.96735537190081 z[991] - 50.91322314049586 z[992] ≥ 0.49553055256021883 
     -61.00413223140495 z[991] + 111.96735537190081 z[992] - 50.91322314049586 z[993] ≥ 0.49602814823269203 
     -61.00413223140495 z[992] + 111.96735537190081 z[993] - 50.91322314049586 z[994] ≥ 0.4965254954404511  
     -61.00413223140495 z[993] + 111.96735537190081 z[994] - 50.91322314049586 z[995] ≥ 0.4970225945553226  
     -61.00413223140495 z[994] + 111.96735537190081 z[995] - 50.91322314049586 z[996] ≥ 0.4975194459482072  
     -61.00413223140495 z[995] + 111.96735537190081 z[996] - 50.91322314049586 z[997] ≥ 0.49801604998908167 
     -61.00413223140495 z[996] + 111.96735537190081 z[997] - 50.91322314049586 z[998] ≥ 0.49851240704700306 
     -61.00413223140495 z[997] + 111.96735537190081 z[998] - 50.91322314049586 z[999] ≥ 0.49900851749011177 
     -61.00413223140495 z[998] + 111.96735537190081 z[999] - 50.91322314049586 z[1000] ≥ 0.49950438168563416
     -61.00413223140495 z[999] + 61.05413223140495 z[1000] ≥ 0.5                                            




```julia
@objective(m, Min, z' * B * z + z' * q); # objective for problem
```


```julia
#@time res = solve(m);
@time res = optimize!(m);
```

    
    ******************************************************************************
    This program contains Ipopt, a library for large-scale nonlinear optimization.
     Ipopt is released as open source code under the Eclipse Public License (EPL).
             For more information visit http://projects.coin-or.org/Ipopt
    ******************************************************************************
    
    This is Ipopt version 3.12.10, running with linear solver mumps.
    NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
    
    Number of nonzeros in equality constraint Jacobian...:        0
    Number of nonzeros in inequality constraint Jacobian.:     3998
    Number of nonzeros in Lagrangian Hessian.............:     1999
    
    Total number of variables............................:     1000
                         variables with only lower bounds:        0
                    variables with lower and upper bounds:        0
                         variables with only upper bounds:        0
    Total number of equality constraints.................:        0
    Total number of inequality constraints...............:     2000
            inequality constraints with only lower bounds:     2000
       inequality constraints with lower and upper bounds:        0
            inequality constraints with only upper bounds:        0
    
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
       0  0.0000000e+00 5.00e-01 1.17e+00  -1.0 0.00e+00    -  0.00e+00 0.00e+00   0
       1 -5.7605060e-01 4.92e-01 1.14e+00  -1.0 4.87e-01    -  2.56e-02 2.30e-02h  1
       2 -9.1960923e-01 4.91e-01 1.61e+00  -1.0 3.47e+00    -  4.81e-03 3.13e-03h  1
       3 -2.0520392e+00 4.89e-01 1.87e+00  -1.0 5.76e+00    -  4.28e-03 3.92e-03h  1
       4 -4.4155948e+00 4.86e-01 1.87e+00  -1.0 6.86e+00    -  5.02e-03 4.84e-03h  1
       5 -8.4489764e+00 4.83e-01 1.86e+00  -1.0 7.37e+00    -  5.87e-03 6.20e-03f  1
       6 -1.4465608e+01 4.79e-01 1.85e+00  -1.0 7.59e+00    -  7.03e-03 7.92e-03f  1
       7 -2.2954073e+01 4.73e-01 2.07e+00  -1.0 7.66e+00    -  8.50e-03 1.03e-02f  1
       8 -3.4644384e+01 4.66e-01 3.63e+00  -1.0 7.65e+00    -  1.04e-02 1.38e-02f  1
       9 -5.0729506e+01 4.57e-01 5.77e+00  -1.0 7.58e+00    -  1.27e-02 1.92e-02f  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
      10 -7.3194051e+01 4.43e-01 9.14e+00  -1.0 7.46e+00    -  1.53e-02 2.81e-02f  1
      11 -1.0134764e+02 4.25e-01 1.39e+01  -1.0 7.25e+00    -  1.64e-02 3.86e-02f  1
      12 -1.2617900e+02 4.08e-01 2.45e+01  -1.0 6.98e+00    -  5.72e-03 3.86e-02f  1
      13 -1.4445614e+02 3.94e-01 2.76e+01  -1.0 6.71e+00    -  5.71e-03 3.22e-02f  1
      14 -2.0576257e+02 3.35e-01 3.67e+01  -1.0 6.49e+00    -  1.99e-02 1.44e-01f  1
      15 -2.3885991e+02 2.09e-01 7.60e+01  -1.0 5.55e+00    -  1.74e-02 3.61e-01f  1
      16  3.8931952e+01 0.00e+00 2.45e+01  -1.0 3.55e+00    -  1.97e-01 1.00e+00h  1
      17  5.0712046e+01 0.00e+00 6.53e+00  -1.0 1.44e-01    -  5.20e-01 1.00e+00f  1
      18  6.7813057e+01 0.00e+00 2.31e+00  -1.0 3.11e-01    -  6.03e-01 1.00e+00f  1
      19  8.9753484e+01 0.00e+00 4.56e-01  -1.0 4.41e-01    -  7.32e-01 1.00e+00f  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
      20  2.9826750e+01 0.00e+00 2.00e-07  -1.7 3.88e-01    -  1.00e+00 1.00e+00f  1
      21  5.6140821e+00 0.00e+00 2.83e-08  -2.5 3.34e-01    -  1.00e+00 1.00e+00f  1
      22  6.8595691e-01 0.00e+00 1.50e-09  -3.8 1.78e-01    -  1.00e+00 1.00e+00f  1
      23  2.1807995e-01 0.00e+00 1.50e-09  -3.8 7.51e-02    -  1.00e+00 1.00e+00f  1
      24  1.5867396e-02 0.00e+00 3.08e-05  -5.7 3.65e-02    -  1.00e+00 9.99e-01f  1
      25  3.5324511e-03 0.00e+00 1.84e-11  -5.7 1.41e-02    -  1.00e+00 1.00e+00f  1
      26  2.0275883e-03 0.00e+00 1.84e-11  -5.7 4.99e-03    -  1.00e+00 1.00e+00f  1
      27  3.7912734e-04 0.00e+00 1.96e-04  -8.6 2.05e-03    -  9.00e-01 8.16e-01f  1
      28  6.9968770e-05 0.00e+00 1.42e-04  -8.6 1.50e-03    -  8.28e-01 7.77e-01f  1
      29  1.3505373e-06 0.00e+00 1.97e-04  -8.6 1.10e-03    -  8.25e-01 7.37e-01f  1
    iter    objective    inf_pr   inf_du lg(mu)  ||d||  lg(rg) alpha_du alpha_pr  ls
      30 -1.8173662e-05 0.00e+00 1.07e-04  -8.6 8.48e-04    -  8.40e-01 7.68e-01f  1
      31 -2.4259210e-05 0.00e+00 5.12e-13  -8.6 1.07e-03    -  1.00e+00 1.00e+00f  1
      32 -2.4336890e-05 0.00e+00 4.55e-13  -8.6 2.31e-04    -  1.00e+00 1.00e+00h  1
      33 -2.6819044e-05 0.00e+00 7.17e-07 -12.9 7.26e-05    -  9.77e-01 9.89e-01h  1
      34 -2.6847078e-05 1.51e-13 5.12e-13 -12.9 4.33e-06    -  1.00e+00 1.00e+00f  1
    
    Number of Iterations....: 34
    
                                       (scaled)                 (unscaled)
    Objective...............:  -2.6847077606362291e-05   -2.6847077606362291e-05
    Dual infeasibility......:   5.1159076974727213e-13    5.1159076974727213e-13
    Constraint violation....:   1.3516965324811281e-13    1.5134588600728054e-13
    Complementarity.........:   8.2060540906518042e-13    8.2060540906518042e-13
    Overall NLP error.......:   8.2060540906518042e-13    8.2060540906518042e-13
    
    
    Number of objective function evaluations             = 35
    Number of objective gradient evaluations             = 35
    Number of equality constraint evaluations            = 0
    Number of inequality constraint evaluations          = 35
    Number of equality constraint Jacobian evaluations   = 0
    Number of inequality constraint Jacobian evaluations = 1
    Number of Lagrangian Hessian evaluations             = 1
    Total CPU secs in IPOPT (w/o function evaluations)   =      2.200
    Total CPU secs in NLP function evaluations           =      0.748
    
    EXIT: Optimal Solution Found.
     11.105514 seconds (24.43 M allocations: 1.191 GiB, 8.91% gc time)



```julia
#plot(getvalue(z) + S, legend = false, title = "Value Function for Vanilla Model")
plot(value.(z) + S, legend = false, title = "Value Function for Vanilla Model")
```




    

    




### Benchmarking (including commercial)
------


```julia
using Suppressor
```


```julia
function LCP_1()
    sp = StoppingProblem()
    @unpack A, B, S, q = LCP_objects(sp)
    f = z -> B*z + q
    n = sp.M
    lb = zeros(n)
    ub = 1000*ones(n) # a reasonable guess?
    @suppress options(convergence_tolerance = 1e-9, output = :no, time_limit = 600) # 10 minute budget
    @suppress solveLCP(f, lb, ub);
end
```




    LCP_1 (generic function with 1 method)




```julia
function LCP_2()
    sp = StoppingProblem()
    @unpack A, B, S, q = LCP_objects(sp)
    f = z -> B*z + q
    n = sp.M
    lb = zeros(n)
    ub = 1000*ones(n) # a reasonable guess?
    function f!(F, z)
        F .= B*z + q # in place assignment
    end
    @suppress mcpsolve(f!, lb, ub,
                zeros(n), # initial condition
                reformulation = :smooth, # uses a so-called "Fischer function" to smooth out the problem
                autodiff = :forward,
                ftol = 1e-9)
end
```




    LCP_2 (generic function with 1 method)




```julia
function LCP_3()
    sp = StoppingProblem()
    @unpack A, B, S, q = LCP_objects(sp)
    f = z -> B*z + q 
    n = sp.M
    #m = Model(solver = IpoptSolver(tol=1e-9))
    m = Model(with_optimizer(Ipopt.Optimizer, tol=1e-9))
    @variable(m, z[1:n])
    @constraint(m, z .>= 0)
    @constraint(m, B*z + q .>= 0)
    @objective(m, Min, z' * B * z + z' * q)
    #@suppress solve(m);
    @suppress optimize!(m);
end
```




    LCP_3 (generic function with 1 method)




```julia
using Gurobi, Mosek, OSQP, MathOptInterfaceMosek
```

    ┌ Info: Recompiling stale cache file /Users/Jozefiak/.julia/compiled/v1.1/Mosek/vS2ln.ji for Mosek [6405355b-0ac2-5fba-af84-adbd65488c0e]
    └ @ Base loading.jl:1184
    ┌ Info: Recompiling stale cache file /Users/Jozefiak/.julia/compiled/v1.1/MathOptInterfaceMosek/sqIeN.ji for MathOptInterfaceMosek [0087ddc6-3964-5e57-817f-9937aefb0357]
    └ @ Base loading.jl:1184



```julia
function LCP_4()
    sp = StoppingProblem()
    @unpack A, B, S, q = LCP_objects(sp)
    f = z -> B*z + q 
    n = sp.M
    #m = Model(solver = GurobiSolver())
    m = Model(with_optimizer(Gurobi.Optimizer))
    @variable(m, z[1:n])
    @constraint(m, z .>= 0)
    @constraint(m, B*z + q .>= 0)
    @objective(m, Min, z' * B * z + z' * q)
    #@suppress solve(m);
    @suppress optimize!(m);
end
```




    LCP_4 (generic function with 1 method)




```julia
function LCP_5()
    sp = StoppingProblem()
    @unpack A, B, S, q = LCP_objects(sp)
    f = z -> B*z + q 
    n = sp.M
    #m = Model(solver = MosekSolver(MSK_DPAR_BASIS_TOL_X=1e-9))
    m = Model(with_optimizer(MosekOptimizer, MSK_DPAR_BASIS_TOL_X=1e-9))
    @variable(m, z[1:n])
    @constraint(m, z .>= 0)
    @constraint(m, B*z + q .>= 0)
    @objective(m, Min, z' * B * z + z' * q)
    #@suppress solve(m);
    @suppress optimize!(m);
end
```




    LCP_5 (generic function with 1 method)




```julia
function LCP_6()
    sp = StoppingProblem()
    @unpack A, B, S, q = LCP_objects(sp)
    f = z -> B*z + q 
    n = sp.M
    #m = Model(solver = OSQPMathProgBaseInterface.OSQPSolver(tol=1e-9))
    m = Model(with_optimizer(OSQP.MathOptInterfaceOSQP.Optimizer))
    @variable(m, z[1:n])
    @constraint(m, z .>= 0)
    @constraint(m, B*z + q .>= 0)
    @objective(m, Min, z' * B * z + z' * q)
    #@suppress solve(m);
    @suppress optimize!(m)
end
```




    LCP_6 (generic function with 1 method)




```julia
using BenchmarkTools
@btime LCP_1();
```

      16.726 ms (989 allocations: 25.29 MiB)



```julia
@btime LCP_2();
```

      609.998 ms (34321 allocations: 634.90 MiB)



```julia
@btime LCP_3();
```

      231.883 ms (245017 allocations: 13.61 MiB)



```julia
# This was not in the notebook originally
#@btime LCP_4(); 
```


```julia
@btime LCP_5();
```


```julia
@btime LCP_6();
```

      184.438 ms (261042 allocations: 15.51 MiB)


d95f090783b533c1118dd51ec9dd3b898d148300
ipynb
Jupyter Notebook
continuous_time_methods/LCP.ipynb
snowdj/computing_and_datascience
3deef184cc389a681cf305296a7f4197831f4469
MIT
2018-07-27T00:56:06.000Z
2021-08-15T00:20:26.000Z
continuous_time_methods/LCP.ipynb
snowdj/computing_and_datascience
3deef184cc389a681cf305296a7f4197831f4469
MIT
2017-12-22T18:26:34.000Z
2019-05-13T22:53:47.000Z
continuous_time_methods/LCP.ipynb
snowdj/computing_and_datascience
3deef184cc389a681cf305296a7f4197831f4469
MIT
2018-01-05T19:05:21.000Z
2021-11-19T17:01:09.000Z
Qwen/Qwen-72B
1. YES
2. YES