madnlp / madnlp.jl Goto Github PK

Now that we have fully migrated to NLPModels, I think it is a great time to implement an interface to Nonconvex.jl. I think we can either create an extension package MadNLPNonconvex.jl in this repository or create a NonconvexMadNLP.jl repository in the JuliaNonconvex organization. I think the latter might be better for consistency with other solver interfaces to Nonconvex.jl. Any suggestions, @frapac @mohamed82008?

Currently, we have NonlinearProgram, which provides basic model information to MadNLP. But all of the same functionality can be provided by NLPModels.jl, and it even has more features. In the future, NonlinearProgram will be replaced by NLPModels.jl , and the current MOI interface will be replaced by a thin wrapper of NLPModelsJuMP.

As described here: https://discourse.julialang.org/t/unsupported-feature-hess-with-user-defined-functions-using-jump-and-alpine/70204/2

The documentation branch is added here. The TODO list is below (will keep updated):

• Home: explain key features. Things to highlight: being a native-julia solver, interface to various modeling tools in julia ecosystem, performance (as fast and scalable as Ipopt), flexibility (abstract KKT system), runs on GPU, and exploits multi-thread in CPU.
• Update linear solver documentation.
• Add docstrings for front-end APIs
• Use DocstringExtensions.jl to document options and outputs.
• Algorithm: tutorial of interior point method
• Algorithm: tutorial of filter line search
• Algorithm: tutorial of intertia correction (inertia-based and inertia-free)
• Algorithm: tutorial of restoration phase
• Example: optimal control (using JuMP.jl)
• Example: optimal power flow (using PowerModels.jl)
• Example: gas network optimization (using Plasmo.jl)
• Example: examples of ADNLPModel.jl
• Example: hand-coded jac/hess example

Like the line-search and the barrier update, the primal-dual regularization should be modular too, to let the choice to the user to extend the regularization algorithm.

It would be good to add support for AMDGPU. We only need to add rocBLAS interface to MadNLPLapackGPU.

I believe the print of iterations shows wrong results or that the solver does not handle this problem. Objective should end at 306.5

http://www.stfmc.de/fmc/rhs/x/wschfh.shtml?ampl/tp015.mod
or
https://vanderbei.princeton.edu/ampl/nlmodels/hs/hs015.mod

Krylov.jl comes with two advantages, in my opinion:

• it implements different iterative solvers useful in optimization (and Krylov.jl is already integrated to Tulip.jl)
• the usual iterative solvers are already working on the GPU

Hi, thanks for this fantastic package. I have a question. I couldn't see any way to not pass the Hessian in and just rely on an LBFGS approximation of the (inverse) Hessian. Is there a way to do this in MadNLP now? If not, is it too much work to add?

using JuMP, MadNLP
m = Model(MadNLP.Optimizer)
@variable(m, x)
@variable(m, y)
@NLobjective(m, Min, (x-5)^2 + (y-8)^2)
@NLconstraint(m, x*y==5)
optimize!(m)
this small example errors in MadNLP. This example runs fine with Ipopt

As mentioned in #58 (comment), several major types in MadNLP include non-concrete types. Because of this, several auxiliary functions are introduced to avoid unnecessary memory allocations. Such introduction of auxiliary functions can be avoided by always using concrete types in the type definition. The following is suggested:

• Thoroughly revisit the type definitions (NonlinearProgram, Solver, SparseMatrixCOO, etc) and eliminate the abstract types in the type definition.
• Eliminate the auxiliary functions whenever there's no performance regression.
• Change Solver -> InteriorPointSolver

Does MadNLP currently support using a lazy Jacobian and Jacobian transpose? It should be possible to use an iterative solver to solve the linear system which doesn't require forming the Jacobian matrix at all.

Currently, I get an error when trying to run with MadNLPLapackGPU.QR or MadNLPLapackGPU.LU. The tril_to_full!() function is only defined for objects of type Matrix, so the function doesn't work with CuArray types on the GPU. An example script that results in this error is shown bleow.

using MadNLP, ADNLPModels, MadNLPGPU, CUDA, NLPModels

function MadNLP.jac_dense!(nls, x, jac) 
    NLPModels.increment!(nls, :neval_jac)
    copyto!(jac, zeros(size(jac)))
end

function MadNLP.hess_dense!(nls, x, w1l, hess; obj_weight = 1.0)
    NLPModels.increment!(nls, :neval_hess)
    copyto!(hess, [0 0; -20 0])
end

F(x) = [x[1] - 1.0; 10 * (x[2] - x[1]^2)]
x0 = [-1.2; 1.0]
nls = ADNLSModel(F, x0, 2)

madnlp_options = Dict{Symbol, Any}(
:kkt_system=>MadNLP.DENSE_CONDENSED_KKT_SYSTEM,
:linear_solver=>MadNLPLapackGPU,
:print_level=>MadNLP.DEBUG,
)

linear_solver_options = Dict{Symbol, Any}(
:lapackgpu_algorithm => MadNLPLapackGPU.QR,
)

TKKTGPU = MadNLP.DenseCondensedKKTSystem{Float64, CuVector{Float64}, CuMatrix{Float64}}
opt = MadNLP.Options(; madnlp_options...)
ips = MadNLP.InteriorPointSolver{TKKTGPU}(nls, opt; option_linear_solver=copy(linear_solver_options))
MadNLP.optimize!(ips)

I get an error in the attached script because the solver is trying to perform feasibility restoration. If I change the tolerance to 1e-6, there is no error. If I run the equivalent problem on the CPU using LapackCPUSolver and MadNLP.InteriorPointSolver instead of LapackGPUSolver and MadNLPGPU.CuInteriorPointSolver, then the problem solves fine and does not try to perform a feasibility restoration.

using MatrixEquations
using MadNLPGPU
using NLPModels, QuadraticModels

function MadNLP.jac_dense!(nlp::DenseLQDynamicModel{T, V, M1, M2, M3}, x, jac) where {T, V, M1<: AbstractMatrix, M2 <: AbstractMatrix, M3 <: AbstractMatrix}
    NLPModels.increment!(nlp, :neval_jac)

    J = nlp.data.A
    copyto!(jac, J)
end

function MadNLP.hess_dense!(nlp::DenseLQDynamicModel{T, V, M1, M2, M3}, x, w1l, hess; obj_weight = 1.0) where {T, V, M1<: AbstractMatrix, M2 <: AbstractMatrix, M3 <: AbstractMatrix}
    NLPModels.increment!(nlp, :neval_hess)
    H = nlp.data.H
    copyto!(hess, H)
end
function build_3D_heating_AB(dx, nx, dt)

    A = zeros(nx^3, nx^3)
    B = zeros(nx^3, 6)

    k = 400. # thermal conductivity of copper, W/(m-K)
    k2 = 400
    rho = 8960. # density of copper, kg/m^3
    specificHeat = 386. # specific heat of copper, J/(kg-K)

    conduction_constant = k * dt / rho / specificHeat / dx^2
    input_constant = k2 * dt / rho / specificHeat / dx^2

    # Set A matrix
    for i in 1:nx^3
        A[i, i] = 1 - 6 * conduction_constant
        # Set links in x direction
        if i%nx != 0 && i%nx != 1
            A[i, i - 1] = conduction_constant
            A[i, i + 1] = conduction_constant
            #y has boundaries if i%100 < 10 or i %100 >90
            #z has boundaries if i%1000 < 100 and i%1000 > 900
            #A[i, i] += -2 * conduction_constant
        elseif i%nx == 0
            A[i, i - 1] = conduction_constant
            #A[i, i] += -1 * conduction_constant
        else
            A[i, i + 1] = conduction_constant
            #A[i, i] += -1 * conduction_constant
        end

        # Set links in the y direction
        if i%(nx^2) in 1:nx
            A[i, i + nx] = conduction_constant
            #A[i, i] += -1 * conduction_constant
        elseif i%(nx^2) == 0 || i%(nx^2) > nx^2 - nx
            A[i, i - nx] = conduction_constant
            #A[i, i] += -1 * conduction_constant
        else
            A[i, i + nx] = conduction_constant
            A[i, i - nx] = conduction_constant
            #A[i, i] += -2 * conduction_constant
        end

        # Set links in the z direction
        if i <= nx^2
            A[i, i + nx^2] = conduction_constant
            #A[i, i] += -1 * conduction_constant
        elseif i > nx^3 - nx^2
            A[i, i - nx^2] = conduction_constant
            #A[i, i] += -1 * conduction_constant
        else
            A[i, i + nx^2] = conduction_constant
            A[i, i - nx^2] = conduction_constant
            #A[i, i] += -2 * conduction_constant
        end
    end

    #Set B matrix
    B[1:nx^2, 1] .= input_constant
    B[(nx^3 - nx^2):(nx^3), 2] .= input_constant
    for i in 1:nx^3
        if i%nx == 1
            B[i, 3] += input_constant
        end
        if i %nx == 0
            B[i, 4] += input_constant
        end
        if i%nx^2 in 1:nx
            B[i, 5] += input_constant
        end
        if i%nx^2 == 0 || i%nx^2 > nx^2 - nx
            B[i, 6] += input_constant
        end
    end
    return A, B
end

function set_d!(d, nx, N, Tmax, Tstart)
    fill!(d, Tstart)
    Tdiff = (Tmax - Tstart)/2
    Tmin = (Tmax - Tstart)/10
    for j in 1:(N + 1)
        for i in 1:nx^3
            x = i%nx
            y = div(i % nx^2, nx)
            z = div(i, nx^2)
            d[i, j] = Tstart + Tmin +  (1 - z/nx) * (2 * sin(3.14159 * x/nx) + 2 * sin(3.14159 * y/nx)) * Tdiff * (j / N / 3)
        end
    end
end

function build_3D_PDE(N, nx, dx, dt, Tmax, Tstart; dense::Bool = true, implicit = false)

    ns = nx^3
    nu = 6

    Q  = 10. * Matrix(LinearAlgebra.I, ns, ns)
    Qf = 10. * Matrix(LinearAlgebra.I, ns, ns)./dt
    R  = 1.0 * Matrix(LinearAlgebra.I, nu, nu)./10

    A, B = build_3D_heating_AB(dx, nx, dt)

    s0 = fill(Tstart, ns)
    sl = fill(200., ns)
    su = fill(550., ns)
    ul = fill(300., nu)
    uu = fill(500., nu)


    S = -.001 * Matrix(I, ns, nu)
    Q_scale = 1
    R_scale = 1

    if dense
        if implicit
            lqdm = DenseLQDynamicModel(s0, A, B, Q, R, N; Qf = Qf, sl = sl, su = su, ul = ul, uu = uu, S = S, implicit=implicit)
        else
            lqdm = DenseLQDynamicModel(s0, A, B, Q, R, N; Qf = Qf, sl = sl, su = su, ul = ul, uu = uu, S = S)
        end
    else
        lqdm = SparseLQDynamicModel(s0, sparse(A), sparse(B), sparse(Q), sparse(R), N; Qf = sparse(Qf), sl = sl, su = su, ul = ul, uu = uu, S = sparse(S))
    end

    d = zeros(nx^3, N + 1)

    set_d!(d, nx, N, Tmax, Tstart)

    block_Q = SparseArrays.sparse([],[],eltype(Q)[], ns * (N + 1), ns * (N + 1))

    for i in 1:N
        block_Q[(1 + (i - 1) * ns):(ns * i), (1 + (i - 1) * ns):(ns * i)] = Q
    end

    block_Q[(1 + ns * N):end, (1 + ns * N):end] = Qf

    Qd    = zeros(size(d, 1))
    Qdvec = zeros(length(d))
    dQd   = 0

    for i in 1:N
        LinearAlgebra.mul!(Qd, Q, d[:, i])
        Qdvec[(1 + ns * (i - 1)):ns * i] = Qd

        dQd += LinearAlgebra.dot(Qd, d[:, i])
    end

    LinearAlgebra.mul!(Qd, Qf, d[:, N + 1])
    Qdvec[(1 + ns * N):end] = Qd

    dQd += LinearAlgebra.dot(Qd, d[:, N + 1])


    # Add c and c0 that result from (x-d)^T Q (x-d) in the objective function
    if dense
        block_A = lqdm.blocks.A
        block_B = lqdm.blocks.B

        As0 = zeros(size(block_A, 1))
        LinearAlgebra.mul!(As0, block_A, s0)
        dQB = zeros(nu * N)
        dQB_sub_block = zeros(nu)

        for i in 1:N
            B_sub_block = block_B[(1 + ns * (i - 1)):ns * i, :]
            for j in N:-1:i
                Qd_sub_block = Qdvec[(1 + ns * j):(ns * (j + 1))]
                LinearAlgebra.mul!(dQB_sub_block, B_sub_block', Qd_sub_block)

                dQB[(1 + nu * (j - i)):nu * (j - i + 1)] .+= dQB_sub_block
            end
        end


        lqdm.data.c0 += dQd / 2
        lqdm.data.c0 += -LinearAlgebra.dot(Qdvec, As0)
        lqdm.data.c  += - dQB
    else
        uvec = zeros(nu * N)
        full_Qd = vcat(Qdvec, uvec)

        lqdm.data.c0 += dQd / 2
        lqdm.data.c  += - full_Qd
    end

    return lqdm
end

N = 250
nx = 9
lenx = .02
dt = .5
Tmax = 350.
Tstart = 300.

lqdm = build_3D_PDE(N, nx, lenx, dt, Tmax, Tstart; dense = true, implicit = false)

madnlp_options = Dict{Symbol, Any}(
    :kkt_system=>MadNLP.DENSE_CONDENSED_KKT_SYSTEM,
    :linear_solver=>LapackGPUSolver,
    :jacobian_constant=>true,
    :hessian_constant=>true,
    :lapack_algorithm=>MadNLP.CHOLESKY,
    :nlp_scaling=>false,
    :max_iter=>100
)

ips1 = MadNLPGPU.CuInteriorPointSolver(lqdm, option_dict=madnlp_options)
sol_ref = MadNLP.optimize!(ips1)```

Currently, fixed variable treatment is done by setting the associated RHS to zero and making the associated block in the augmented matrix identity. This should be changed to actually reducing the dimension of the KKT system.

Right now, instantiating a KKT system on the GPU amounts to

TKKTGPU = MadNLP.DenseKKTSystem{Float64, Vector{Float64}, Matrix{Float64}}
# Instantiate Solver with KKT on the GPU
d_ips = MadNLP.InteriorPointSolver{TKKTGPU}(lq_dense, opt; option_linear_solver=copy(madnlp_options))
MadNLP.optimize!(d_ips)

which is overly complicated.

Instead, we can have a CuInteriorPointSolver which would take care to instantiate everything properly on the GPU. Ideally, the API should be very similar to CPU's API:

# Instantiate Solver with KKT on the GPU
d_ips = MadNLP.CuInteriorPointSolver(lq_dense; option_dict=copy(madnlp_options))
MadNLP.optimize!(d_ips)

The heuristics discussed in Section 3.2 of the Ipopt paper is not implemented in MadNLP. This should be implemented for the performance improvement

@Article{wachter2006implementation,
title={On the implementation of an interior-point filter line-search algorithm for large-scale nonlinear programming},
author={W{"a}chter, Andreas and Biegler, Lorenz T},
journal={Mathematical programming},
volume={106},
number={1},
pages={25--57},
year={2006},
publisher={Springer}
}

I'm using Julia 1.7 with updated MadNLP 0.4.0:

(@v1.7) pkg> status MadNLP
      Status `~/.julia/environments/v1.7/Project.toml`
  [2621e9c9] MadNLP v0.4.0

Running @info "test" 1 in the REPL produces the following output (as expected):

julia> @info "TEST" 1
┌ Info: TEST
└   1 = 1

When I first import MadNLP and then @info "test" 1, I get the following error:

julia> import MadNLP
julia> @info "TEST" 1
ERROR: MethodError: no method matching get_level(::String)
Closest candidates are:
  get_level(::MadNLP.Logger) at ~/.julia/packages/MadNLP/OUcjZ/src/utils.jl:20
Stacktrace:
 [1] top-level scope
   @ ~/.julia/packages/MadNLP/OUcjZ/src/utils.jl:34

MadNLP seems to overwrite the default logging macros, which is not ideal when using the solver in larger projects.

At the moment, when using Bunch-Kaufman with LapackGPU, MadNLP is using cusolver to compute the factorization on the GPU, then copy the factorization back to the CPU to do the triangular solver with OpenBlas.

The reason of this additional data transfer is that the routine sytrs is not implemented properly in cusolver (there exists a function cusolverDnDsytrs in the API, but apparently it doesn't do anything, and this function was not even referenced in the cusolver's documentation).

However, it looks like they have included a new sytrs function in CUDA 11.3.1: cf the documentation (strangely, they don't mention this new function in the release note).
It would be interesting to test this new function inside MadNLP (it would require to wrap the new cusolverDnXsytrs inside CUDA.jl). If it works, we would save one expensive data transfer :-)

I have been making a link constraint in Plasmo between a variable on a master node of the primary OptiGraph and corresponding variables on subgraphs. I found that the order of the link constraint matters in determining the size of the Schur complement, and I am not sure if this should be the case. If I do @linkconstraint(graph, subgraph[:variable] == master[:variable], the Schur complement size only depends on the number of variables on the master node, but if I do @linkconstraint(graph, master[:variable] == subgraph[:variable] then the Schur complement size increases for every added subgraph variable being linked to the master variable.

I have included an example below. This difference is seen in lines 24 and 25 (one of them can be commented out when running the code). The script will print the size of the Schur complement, and it is not the same between the two formulations. Additionally, the number of subgraph variables can be increased at lines 8 and 9 which will change the Schur complement size if the link constraint is expressed as master[:variable] == subgraph[:variable].

schur_error.zip

Sometimes MadNLP seems to struggle if inertia_correction_method=MadNLP.INERTIA_BASED option is used (the default is INERTIA_AUTO, but it is set to INERTIA_BASED when direct linear solvers are used). Such a case can be reproduced with this script. Here,

MadNLP.optimize!(gas_network_stoch;inertia_correction_method=MadNLP.INERTIA_FREE)  # works well, but
MadNLP.optimize!(gas_network_stoch;inertia_correction_method=MadNLP.INERTIA_BASED) # doesn't

Potentially there is a major bug in the inertia correction procedure. Originally reported by @dlcole3.

Right now, the option DENSE_CONDENSED_KKT_SYSTEM does not reduce the equality constraints and focuses only on the inequality constraints. One has to use an elastic relaxation to reformulate the equality constraints as inequalities before condensing the system.

We should investigate a proper way to condense both the equality and the inequality constraints, without resorting to an elastic relaxation trick.

I think the current issue is caused by not having a consistent dimension in https://github.com/sshin23/MadNLP.jl/blob/master/src/interiorpointsolver.jl#L1010. It can be fixed by specifying the operation based on kkt.hess, not kkt.aug_com.

Currently MadNLP assumes that the sparsity is given in COO format. This option allows the user to tell MadNLP if the sparsity is in a certain from (dense or CSC)

Currently, InteriorPointSolver keeps views of primal vectors, such as

x_lr, x_ur, dx_lr, dx_ur, etc.

We can create a new data structure PrimalVector and keep these views inside PrimalVector to make the data more structured. Any opinion about this @frapac? I was assuming these #TODO's are inserted by you :)

On the line

In the same file

nlp.meta.x0, nlp.meta.y0

I think currently MadNLP comes with too many dependencies, which has a significant impact on the precompilation time of the package. Unfortunately, Julia does not support properly conditional dependencies (see the long awaited issue JuliaLang/Pkg.jl#1285). But, since 2020, we can define several packages with different UUID under the same git repo (see this discussion, or this PR for a practical example).

I would suggest to refactor MadNLP the following way, by keeping all the code in this repo.

• MadNLP/src would implement only the core interior point algorithm, interface to usual linear solvers (as LapackCPU, or IterativeSolvers) and the lightweights interface (NLPModels and MathOptInterface).
• we would move in a new directory MadNLP/lib all "complicated" deps. In MadNLP/lib, each subfolder would be registered as Julia package with its own UUID in Julia's registry.

For instance, we could structure MadNLP/lib as

lib/MadNLPPlasmo    # interface to Plasmo
lib/MadNLPCUDA      # interface to CUDA's Lapack
lib/MadNLPHSL         # interface to HSL
lib/MadNLPMumps    # interface to MUMPS
lib/MadNLPPardiso    # interface to Pardiso 

Leading to 5 subpackages ... For the three linear solvers HSL, Mumps and Pardiso, I am wondering if we could build something in common with JuliaSmoothOptimizers, to adopt a common interface for the linear solvers we are using. I have a look at HLS.jl, MUMPS.jl but these two packages implement only a bare-bone wrapper to the C code. I think medium term, it would be nice to share the same interface :

abstract type AbstractLinearSolver end

# dummy functions
introduce(::EmptyLinearSolver) = ""
factorize!(M::EmptyLinearSolver) = M
solve!(::EmptyLinearSolver,x) = x
is_inertia(::EmptyLinearSolver) = false
inertia(::EmptyLinearSolver) = (0,0,0)
improve!(::EmptyLinearSolver) = false
rescale!(::AbstractLinearSystemScaler) = nothing
solve_refine!(y,::AbstractIterator,x) = nothing

With DocStringExtensions, we can generate docstring for the Type's fields with the following syntax:

"""
This is an X
"""
struct X
    "This is a"
    a::String
    "This is b"
    b::Char
    "This is c"
    c
end

One shortcoming is that DocStringExtensions should be added to the dependency (not to the /doc/Project.toml, but to the main /Project.toml). But DocStringExtensions seems to be reasonably lightweight. See the discussion below:
https://discourse.julialang.org/t/documenting-elements-of-a-struct/64769

It seems to me that this might be the proper way to manage the long option documentation (which is generated manually for now). Let me know if there are any other alternatives

See: #73 (comment)

This would allows to run MadNLP on the GPU with inertia based regularization. The main showstopper is the sequential nature of the inertia computation algorithm: we may want to have a look at what is done inside Magma.

I'm a bit confused about how to interpret the result
The example https://madnlp.github.io/MadNLP.jl/dev/quickstart/ contains two variables, but

ips.x

4-element Vector{Float64}:
 -0.7921232178470456
 -1.2624298435831807
  0.9999999900052515
  0.801605892122408

has length 4. What is the definition of ips.x and which variables correspond to x[1:2]?

Another question, how do I interpret the objective value at the solution?

julia> ips.obj_val
360.37976240508465

julia> ips.obj_scale
1-element Vector{Float64}:
 1.0

in this case the scaling is 1 so I assume that the actual, unscaled, objective value is given by ips.obj_val. If the scale is not 1, do I compute the unscaled value as ips.obj_val * ips.obj_scale?

I settled the env variable pointing to the HSL source location but could not install nor build the MadNLPHSL extension.

pkg> build MadNLPHSL
    Building MadNLPHSL → `C:\Users\ilyao\.julia\scratchspaces\44cfe95a-1eb2-52ea-b672-e2afdf69b78f\78742b9e9faaba0ade5617cdaca5e34f63c5fa8e\build.log`
ERROR: Error building `MadNLPHSL`:
ERROR: LoadError: Unable to open libLLVM!
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] (::BinaryProvider.var"#open_libllvm#124")()
    @ BinaryProvider C:\Users\ilyao\.julia\packages\BinaryProvider\U2dKK\src\PlatformNames.jl:652
  [3] detect_cxx11_string_abi()
    @ BinaryProvider C:\Users\ilyao\.julia\packages\BinaryProvider\U2dKK\src\PlatformNames.jl:655
  [4] detect_compiler_abi()
    @ BinaryProvider C:\Users\ilyao\.julia\packages\BinaryProvider\U2dKK\src\PlatformNames.jl:668
  [5] top-level scope
    @ C:\Users\ilyao\.julia\packages\BinaryProvider\U2dKK\src\PlatformNames.jl:685
  [6] include(mod::Module, _path::String)
    @ Base .\Base.jl:419
  [7] include(x::String)
    @ BinaryProvider C:\Users\ilyao\.julia\packages\BinaryProvider\U2dKK\src\BinaryProvider.jl:1
  [8] top-level scope
    @ C:\Users\ilyao\.julia\packages\BinaryProvider\U2dKK\src\BinaryProvider.jl:12
  [9] include
    @ .\Base.jl:419 [inlined]
 [10] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String)
    @ Base .\loading.jl:1554
 [11] top-level scope
    @ stdin:1
in expression starting at C:\Users\ilyao\.julia\packages\BinaryProvider\U2dKK\src\PlatformNames.jl:685
in expression starting at C:\Users\ilyao\.julia\packages\BinaryProvider\U2dKK\src\BinaryProvider.jl:1
in expression starting at stdin:1
ERROR: LoadError: Failed to precompile BinaryProvider [b99e7846-7c00-51b0-8f62-c81ae34c0232] to C:\Users\ilyao\.julia\compiled\v1.8\BinaryProvider\jl_B958.tmp.
Stacktrace:
  [1] error(s::String)
    @ Base .\error.jl:35
  [2] compilecache(pkg::Base.PkgId, path::String, internal_stderr::IO, internal_stdout::IO, keep_loaded_modules::Bool)
    @ Base .\loading.jl:1705
  [3] compilecache
    @ .\loading.jl:1649 [inlined]
  [4] _require(pkg::Base.PkgId)
    @ Base .\loading.jl:1337
  [5] _require_prelocked(uuidkey::Base.PkgId)
    @ Base .\loading.jl:1200
  [6] macro expansion
    @ .\loading.jl:1180 [inlined]
  [7] macro expansion
    @ .\lock.jl:223 [inlined]
  [8] require(into::Module, mod::Symbol)
    @ Base .\loading.jl:1144
  [9] include(fname::String)
    @ Base.MainInclude .\client.jl:476
 [10] top-level scope
    @ none:5
in expression starting at C:\Users\ilyao\.julia\packages\MadNLPHSL\8XAiY\deps\build.jl:1

This is my system info:

pkg> st MadNLP MadNLPHSL
Status `D:\ilyao\Documents\research\neural_ode\transcription_neuralode\Project.toml`
  [2621e9c9] MadNLP v0.5.0
  [7fb6135f] MadNLPHSL v0.3.0

julia> versioninfo()
Julia Version 1.8.0
Commit 5544a0fab7 (2022-08-17 13:38 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: 16 × AMD Ryzen 7 4800H with Radeon Graphics
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-13.0.1 (ORCJIT, znver2)
  Threads: 1 on 16 virtual cores

Following the discussion in #33 , we figure out it would be great to support dense nonlinear programming in MadNLP. The motivation is to avoid unneeded memory transfer when we are solving the KKT system with Lapack (on the CPU or on the GPU).

Currently, MadNLP is storing the Hessian and the Jacobian in a SparseMatrixCOO matrix object, defined inside Solver. When solving one IPM iteration, the procedure is:

• evaluate the Jacobian and the Hessian in sparse format by passing to the callbacks contiguous views porting on the COO matrix. After that operations, the COO matrix is up-to-date. This operation is efficient, as the views are contiguous here.
• then, the sparse COO matrix is transferred to another matrix aug_com, storing the augmented system being solved by the linear solvers (we convert COO to CSC if we are using a sparse linear solver, and COO to dense if we are using a dense linear solver).

When solving a dense nonlinear problem, this procedure is inefficient, as we need to transfer first the dense Hessian to the SparseMatrixCOO, then the SparseMatrixCOO to another dense matrix. Instead, we could envision skipping the intermediate steps.
The problem here is we can't pass directly the KKT matrix to the Hessian and the Jacobian callbacks. Thus, it's difficult to avoid some copy operations there, as we need to copy the values of the Hessian and the Jacobian to the KKT system before solving the KKT system with the linear solver. We could abstract the problem as passing the Jacobian and Hessian matrices to the KKT system, using an Intermediate Matrix (IM):

Callbacks Matrices -> Intermediate Matrix -> KKT matrix 

with currently the Intermediate Matrix being the SparseMatrixCOO used by MadNLP. The passage between Callback Matrices and Intermediate matrix are transparent as we are using views.

I see two solutions forward, one easy, the other more complicated:

• we keep the SparseMatrixCOO as an IntermediateMatrix, and an an option to specify that the Hessian and the Jacobian are dense (so we don't need to compute the sparsity pattern). The SparseMatrixCOO stores only the upper or the lower triangular part of the KKT systems. That solution would require an efficient function to transfer the data from the Intermediate Matrix to the KKT matrix, though.
• or, instead of storing a SparseMatrixCOO, we store two dense matrices corresponding respectively to the Jacobian and the Hessian. These matrices would be passed directly to the callbacks. Then, it would remain to implement a function to build the KKT matrix in an efficient way.

I welcome any feedback on this, as I am not sure of the implementation! I have some time available in the following days, so if you want, I could give a try at implementing the support for dense nonlinear problems.

This change will allow running MadNLP with different precisions (as long as compatible derivative callbacks and linear solvers exist).

Right now, we pass the options to MadNLP as a Dict madnlp_options. The problem is that MadNLP deletes the entries in madnlp_options, requiring us to use copy to keep the initial settings in memory. Ideally, MadNLP should not delete the entries in madnlp_options

We currently have an issue when using custom KKTSystem structure, when the dimension of the reduced matrix does not match the number of variable ips.n in the original problem: https://github.com/MadNLP/MadNLP.jl/blob/master/src/interiorpointsolver.jl#L1065

If the size of the linear system is r x r, with r < ips.n, then the following check fails:

    num_pos != ips.n

even if the matrix is full rank.

We should implement the test at the KKTSystem level. We could implement a new structure to store the inertia information returned by the solver:

struct Inertia 
   n::Int
   p::Int
   q::Int
end

and pass the inertia to the KKTSystem to check if the KKT matrix is full rank:

is_inertia_correct(kkt::KKTSystem, inertia::Inertia)::Bool

Hello there,

Solving a QP using MadNLP seems to fail on my test for julia-nightly on ubuntu-latest. Here is the error

/home/runner/work/GeneralizedSmoothingSplines.jl/GeneralizedSmoothingSplines.jl/test/runtests.jl:7
  Got exception outside of a @test
  LoadError: ArgumentError: strides is invalid for SubArrays with indices of type Vector{Int64}
  Stacktrace:
    [1] substrides(strds::Tuple{Int64}, I::Tuple{Vector{Int64}})
      @ Base ./subarray.jl:366
    [2] strides(V::SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false})
      @ Base ./subarray.jl:360
    [3] stride(V::SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false}, d::Int64)
      @ Base ./subarray.jl:368
    [4] vec_pointer_stride
      @ /opt/hostedtoolcache/julia/nightly/x64/share/julia/stdlib/v1.9/LinearAlgebra/src/blas.jl:153 [inlined]
    [5] axpy!(alpha::Float64, x::Vector{Float64}, y::SubArray{Float64, 1, Vector{Float64}, Tuple{Vector{Int64}}, false})
      @ LinearAlgebra.BLAS /opt/hostedtoolcache/julia/nightly/x64/share/julia/stdlib/v1.9/LinearAlgebra/src/blas.jl:494
    [6] regular!(ips::MadNLP.InteriorPointSolver{MadNLP.SparseKKTSystem{Float64, SparseArrays.SparseMatrixCSC{Float64, Int32}}})
      @ MadNLP ~/.julia/packages/MadNLP/4CmQz/src/interiorpointsolver.jl:780
    [7] optimize!(ips::MadNLP.InteriorPointSolver{MadNLP.SparseKKTSystem{Float64, SparseArrays.SparseMatrixCSC{Float64, Int32}}})
      @ MadNLP ~/.julia/packages/MadNLP/4CmQz/src/interiorpointsolver.jl:612
    [8] optimize!(model::MadNLP.Optimizer)
      @ MadNLP ~/.julia/packages/MadNLP/4CmQz/src/Interfaces/MOI_interface.jl:1168
    [9] optimize!
      @ ~/.julia/packages/MathOptInterface/ajp5T/src/Bridges/bridge_optimizer.jl:354 [inlined]
   [10] optimize!
      @ ~/.julia/packages/MathOptInterface/ajp5T/src/MathOptInterface.jl:87 [inlined]
   [11] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{MadNLP.Optimizer}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}})

Looking at the error, it seems to refer to a interpointsolver.jl, which does not exist anymore en this repo, so it is kind of hard to debug this.

Cheers,
Mikkel

To be consistent with other JSO solvers, I suggest making the following changes:

InteriorPointSolver -> MadNLPSolver
IPMOptions -> MadNLPOptions
optimize! -> solve!

Now that the StridedOneVector PR #174 has been merged, we can focus on the time to first solve.

Using the following script:

using MadNLP
using MadNLPTests

nlp = MadNLPTests.HS15Model()
ips = MadNLP.InteriorPointSolver(nlp)
@time MadNLP.optimize!(ips)

MadNLP takes for the first solve:

Total wall-clock secs in solver (w/o fun. eval./lin. alg.)  =  8.082
Total wall-clock secs in linear solver                      =  0.590
Total wall-clock secs in NLP function evaluations           =  0.037
Total wall-clock secs                                       =  8.709

EXIT: Optimal Solution Found.
  8.147811 seconds (24.28 M allocations: 1.187 GiB, 6.16% gc time, 99.93% compilation time)

It would be ideal if we can decrease that time to improve the user experience.

The construction of an MadNLPSolver allocates a lot of memory, and repeatedly solving similar optimization problems by defining a new MadNLPSolver for each problem is thus associated with a lot of overhead.

Would it be possible to come up with either

1. An interface to update things like the initial guess (warm start), upper and lower bounds etc.
2. A constructor of MadNLPSolver that takes an existing MadNLPSolver to re-use the allocated memory

My motivation is solving MPC problems where each consecutive differs from the last only in numerical values of bounds, the structure remains the same.

At the moment MadNLP is using a classical backtracking line-search. We should implement a modular way to specify the line-search algorithm, and give the possibility to the user to use its own line-search algorithm.

When I run the NLPModels example in the docs, the number of bounds that are reported by the solver appears incorrect

Total number of variables............................:        2
                     variables with only lower bounds:        1
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        0
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        2
        inequality constraints with only lower bounds:        2
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

From the code below, I expect

Total number of variables............................:        2
                     variables with only lower bounds:        0
                variables with lower and upper bounds:        0
                     variables with only upper bounds:        1
Total number of equality constraints.................:        0
Total number of inequality constraints...............:        2
        inequality constraints with only lower bounds:        2
   inequality constraints with lower and upper bounds:        0
        inequality constraints with only upper bounds:        0

using MadNLP
using MadNLP.NLPModels


struct HS15Model <: NLPModels.AbstractNLPModel{Float64,Vector{Float64}}
    meta::NLPModels.NLPModelMeta{Float64, Vector{Float64}}
    counters::NLPModels.Counters
end

function HS15Model(x0)
    return HS15Model(
        NLPModels.NLPModelMeta(
            2,     #nvar
            ncon = 2,
            nnzj = 4,
            nnzh = 3,
            x0 = x0,
            y0 = zeros(2),
            lvar = [-Inf, -Inf],
            uvar = [0.5, Inf],
            lcon = [1.0, 0.0],
            ucon = [Inf, Inf],
            minimize = true
        ),
        NLPModels.Counters()
    )
end

function NLPModels.obj(nlp::HS15Model, x::AbstractVector)
    return 100.0 * (x[2] - x[1]^2)^2 + (1.0 - x[1])^2
end

function NLPModels.grad!(nlp::HS15Model, x::AbstractVector, g::AbstractVector)
    z = x[2] - x[1]^2
    g[1] = -400.0 * z * x[1] - 2.0 * (1.0 - x[1])
    g[2] = 200.0 * z
    return
end

function NLPModels.cons!(nlp::HS15Model, x::AbstractVector, c::AbstractVector)
    c[1] = x[1] * x[2]
    c[2] = x[1] + x[2]^2
end

function NLPModels.jac_structure!(nlp::HS15Model, I::AbstractVector{T}, J::AbstractVector{T}) where T
    copyto!(I, [1, 1, 2, 2])
    copyto!(J, [1, 2, 1, 2])
end

function NLPModels.jac_coord!(nlp::HS15Model, x::AbstractVector, J::AbstractVector)
    J[1] = x[2]    # (1, 1)
    J[2] = x[1]    # (1, 2)
    J[3] = 1.0     # (2, 1)
    J[4] = 2*x[2]  # (2, 2)
end

function NLPModels.hess_structure!(nlp::HS15Model, I::AbstractVector{T}, J::AbstractVector{T}) where T
    copyto!(I, [1, 2, 2])
    copyto!(J, [1, 1, 2])
end

function NLPModels.hess_coord!(nlp::HS15Model, x, y, H::AbstractVector; obj_weight=1.0)
    # Objective
    H[1] = obj_weight * (-400.0 * x[2] + 1200.0 * x[1]^2 + 2.0)
    H[2] = obj_weight * (-400.0 * x[1])
    H[3] = obj_weight * 200.0
    # First constraint
    H[2] += y[1] * 1.0
    # Second constraint
    H[3] += y[2] * 2.0
end

x0 = zeros(2) # initial position
nlp = HS15Model(x0)
ips = MadNLP.MadNLPSolver(nlp)
MadNLP.solve!(ips)

It would be nice if the error message below

EXIT: Invalid number in NLP function or derivative detected.

indicated which user-defined function was returning an invalid number.

This issue is used to trigger TagBot; feel free to unsubscribe.

If you haven't already, you should update your TagBot.yml to include issue comment triggers.
Please see this post on Discourse for instructions and more details.

If you'd like for me to do this for you, comment TagBot fix on this issue.
I'll open a PR within a few hours, please be patient!

Hello,
I have installed MadNLP with the command Pkg.add("MadNLP") and added using MadNLP in my code.
At the first run the following error showed up:
ERROR: LoadError: syntax: extra token "as" after end of expression
Stacktrace:
[1] top-level scope at C:\Users\umbe0\.julia\packages\MadNLP\FNmbj\src\MadNLP.jl:16
[2] include(::Module, ::String) at .\Base.jl:377
[3] top-level scope at none:2
[4] eval at .\boot.jl:331 [inlined]
[5] eval(::Expr) at .\client.jl:449
[6] top-level scope at .\none:3
in expression starting at C:\Users\umbe0\.julia\packages\MadNLP\FNmbj\src\MadNLP.jl:16
ERROR: LoadError: Failed to precompile MadNLP [2621e9c9-9eb4-46b1-8089-e8c72242dfb6] to
C:\Users\umbe0\.julia\compiled\v1.4\MadNLP\Ru3yE_ccVQJ.ji.
Stacktrace:
[1] error(::String) at .\error.jl:33
[2] compilecache(::Base.PkgId, ::String) at .\loading.jl:1272
[3] _require(::Base.PkgId) at .\loading.jl:1029
[4] require(::Base.PkgId) at .\loading.jl:927
[5] require(::Module, ::Symbol) at .\loading.jl:922

Does anyone know how to solve it?

DELETED

Currently, MadNLP assumes COO structure. But the user should be able to tell MadNLP if the sparsity structure is CSC or dense

In our team, we ran into the following issue when installing MadNLP.jl and Ipopt.jl in the same Julia session:

session> test MadNLP 
...
free(): invalid pointer

signal (6): Aborted
in expression starting at /home/fpacaud/exa/MadNLP.jl/test/MOI_interface_test.jl:12
gsignal at /lib/x86_64-linux-gnu/libc.so.6 (unknown line)
abort at /lib/x86_64-linux-gnu/libc.so.6 (unknown line)
unknown function (ip: 0x7fe0b3a01906)
unknown function (ip: 0x7fe0b3a08979)
cfree at /lib/x86_64-linux-gnu/libc.so.6 (unknown line)
dmumps_end_driver_ at /home/fpacaud/.julia/artifacts/9d7c145c183adb9f2f4df30a7531e4ef1c95cdc3/lib/libdmumps.so (unknown line)
dmumps_ at /home/fpacaud/.julia/artifacts/9d7c145c183adb9f2f4df30a7531e4ef1c95cdc3/lib/libdmumps.so (unknown line)
dmumps_f77_ at /home/fpacaud/.julia/artifacts/9d7c145c183adb9f2f4df30a7531e4ef1c95cdc3/lib/libdmumps.so (unknown line)
dmumps_c at /home/fpacaud/.julia/artifacts/9d7c145c183adb9f2f4df30a7531e4ef1c95cdc3/lib/libdmumps.so (unknown line)

As detailed in the stacktrace, the issue is caused by MUMPS. Investigating it further:

• MadNLP is working fine with MUMPS_seq_jll 5.3.5
• The crash could be reproduced if we install MadNLP alone and fix explicitly the version of MUMPS_seq_jll to 5.2.1
• the problem is that Ipopt.jl requires explicitly MUMPS_seq_jll 5.2.1: see the related discussion here jump-dev/Ipopt.jl#242

I see three solutions onward:

• we fix explicitly the version of MUMPS_seq_jll to 5.3.5 in MadNLP's Project.toml (in that case, we could not use anymore Ipopt and MadNLP together)
• we fix the error we are observing with MUMPS_seq_jll 5.2.1 (but that would be hard, I think)
• we wait till Ipopt.jl supports MUMPS_seq_jll 3.5.1: jump-dev/Ipopt.jl#246

We should allow the user to change the update strategy for the barrier term mu (fixed, monotone, adaptive, Mehrotra predictor-corrector type).

Currently, the solve! function is not properly tested, and a test need to be added

As noticed in #224, if a variable is fixed (its lower-bound equates its upper-bound) in MadNLP its multiplier is set to 0, which is wrong. The bound dual need to be properly set when terminating the solution.