diff --git a/README.md b/README.md index 72b98d8..1ff1409 100644 --- a/README.md +++ b/README.md @@ -2,11 +2,11 @@ # Pajarito -Pajarito is a **mixed-integer convex programming** (MICP) solver package written in [Julia](http://julialang.org/). +Pavito is a **mixed-integer convex programming** (MICP) solver package written in [Julia](http://julialang.org/). MICP problems are convex except for restrictions that some variables take binary or integer values. -MICP problems are convex except for restrictions that some variables take binary or integer values. Pajarito supports both **mixed-integer conic programming**, which encodes nonlinearities using a small number of predefined cones, and more traditional **convex mixed-integer nonlinear programming**, which encodes nonlinearities with smooth functions and uses their derivatives. +Pavito solves MICP problems by constructing sequential polyhedral outer-approximations of the convex feasible set, similar to [Bonmin](https://projects.coin-or.org/Bonmin). The underlying algorithm has theoretical finite-time convergence under reasonable assumptions. Pavito accesses state-of-the-art MILP solvers and continuous convex conic solvers through the MathProgBase interface. It currently accepts mixed-integer conic problems with second-order, rotated second-order, (primal) exponential, and positive semidefinite cones. MISOCP and MISDP are two established sub-classes of MICPs that Pajarito can solve. -Pajarito solves MICP problems by constructing sequential lifted polyhedral outer-approximations of the convex feasible set to leverage the power of MILP solvers. The algorithm has theoretical finite-time convergence under reasonable assumptions. Pajarito accesses state-of-the-art MILP solvers and continuous conic or nonlinear programming (NLP) solvers through the MathProgBase interface. +For algorithms that use a derivative-based nonlinear programming (NLP) solver (e.g. Ipopt) instead of a conic solver, use [Pavito](https://github.com/JuliaOpt/Pavito.jl). Pavito is a convex mixed-integer nonlinear programming (convex MINLP) solver. As Pavito relies on gradient cuts, it can fail near points of nondifferentiability. Pajarito may be more robust than Pavito on nonsmooth problems (such as MISOCPs). ## Installation @@ -17,91 +17,82 @@ julia> Pkg.add("Pajarito") ## Usage -Pajarito has two entirely separate algorithms depending on the form in which the input is provided. The first algorithm is the **derivative-based nonlinear** (NLP) algorithm, where the approach used is analogous to that of [Bonmin](https://projects.coin-or.org/Bonmin) (and will perform similarly, with the primary advantage of being able to easily swap-in various mixed-integer and convex subproblem solvers that Bonmin does not support). The second algorithm is the **conic** algorithm, which is the new approach proposed in the publications describing Pajarito. The conic algorithm currently supports MICP with second-order, rotated second-order, (primal) exponential, and positive semidefinite cones. MISOCP and MISDP are two established sub-classes of MICPs that Pajarito can solve. - There are several convenient ways to model MICPs in Julia and access Pajarito: | | [JuMP][JuMP-url] | [Convex.jl][convex-url] | [MathProgBase][mpb-url] | |-------------|-------------------|--------------------------|--------------------------| -| NLP model | [X][JuMP-nlp-url] | | [X][mpb-nlp-url] | | Conic model | X | X | [X][mpb-conic-url] | -[mpb-nlp-url]: http://mathprogbasejl.readthedocs.io/en/latest/nlp.html [mpb-conic-url]: http://mathprogbasejl.readthedocs.io/en/latest/conic.html [JuMP-url]: https://github.com/JuliaOpt/JuMP.jl -[JuMP-nlp-url]: http://jump.readthedocs.io/en/latest/nlp.html [convex-url]: https://github.com/JuliaOpt/Convex.jl [mpb-url]: https://github.com/JuliaOpt/MathProgBase.jl -JuMP and Convex.jl are algebraic modeling interfaces, while MathProgBase is a lower-level interface for providing input in raw callback or matrix form. Convex.jl is perhaps the most user-friendly way to provide input in conic form, since it transparently handles conversion of algebraic expressions. JuMP supports general nonlinear smooth functions, e.g. by using `@NLconstraint`. JuMP also supports conic modeling, but requires cones to be explicitly specified, e.g. by using `norm(x) <= t` for second-order cone constraints and `@SDconstraint` for positive semidefinite constraints. Note that a problem provided in conic form may solve faster than an equivalent problem in NLP form because conic form naturally encodes extended formulations; however, [Hijazi et al.](http://www.optimization-online.org/DB_FILE/2011/06/3050.pdf) suggest manual reformulation techniques which achieve many of the algorithmic benefits of conic form. - -Pajarito may be accessed through MathProgBase from outside Julia by using the experimental [cmpb](https://github.com/mlubin/cmpb) interface which provides a C API to the low-level conic input format. The [ConicBenchmarkUtilities](https://github.com/mlubin/ConicBenchmarkUtilities.jl) package provides utilities to read files in the [CBF](http://cblib.zib.de/) format. +JuMP and Convex.jl are algebraic modeling interfaces, while MathProgBase is a lower-level interface for providing input in raw callback or matrix form. Convex.jl is perhaps the most user-friendly way to provide input in conic form, since it transparently handles conversion of algebraic expressions. JuMP supports conic modeling, but requires cones to be explicitly specified, e.g. by using `norm(x) <= t` for second-order cone constraints and `@SDconstraint` for positive semidefinite constraints. Pajarito may be accessed through MathProgBase from outside Julia by using the experimental [cmpb](https://github.com/mlubin/cmpb) interface which provides a C API to the low-level conic input format. The [ConicBenchmarkUtilities](https://github.com/mlubin/ConicBenchmarkUtilities.jl) package provides utilities to read files in the [CBF](http://cblib.zib.de/) format. ## MIP and continuous solvers -The algorithm implemented by Pajarito itself is relatively simple, and most of the hard work is performed by the MIP outer-approximation model solver and the continuous convex subproblem models solver. **The performance of Pajarito depends on these two types of solvers.** For best performance, use commercial solvers. - -The mixed-integer solver is specified by using the `mip_solver` option to `PajaritoSolver`, e.g. `PajaritoSolver(mip_solver=CplexSolver())`. You must first load the Julia package which provides the mixed-integer solver, e.g. `using CPLEX`. This solver is typically a mixed-integer linear solver, but if a conic problem has both second-order cones and other nonlinear cones, or if it has PSD cones, then the MIP solver can be an MISOCP solver and Pajarito can put second-order cones in the outer-approximation model. +The algorithm implemented by Pajarito itself is relatively simple, and most of the hard work is performed by the MIP solver and the continuous conic solver. **The performance of Pajarito depends on these two types of solvers.** For best performance, use commercial solvers. -The continuous convex solver is specified by using the `cont_solver` option, e.g. `PajaritoSolver(cont_solver=IpoptSolver())`. When given input in derivative-based nonlinear form, Pajarito requires a derivative-based nonlinear solver, e.g. [Ipopt](https://projects.coin-or.org/Ipopt) or [KNITRO](http://www.ziena.com/knitro.htm). When given input in conic form, the convex subproblem solver can be *either* a conic solver which supports the cones used *or* a derivative-based solver like Ipopt. If a derivative-based solver is provided in this case, then Pajarito will switch to the derivative-based algorithm by using the [ConicNonlinearBridge](https://github.com/mlubin/ConicNonlinearBridge.jl) (which does not support PSD cones). Note that using derivative-based solvers for conic problems can cause numerical instability because conic problems are not always smooth. +The mixed-integer solver is specified by using the `mip_solver` option to `PajaritoSolver`, e.g. `PajaritoSolver(mip_solver=CplexSolver())`. You must first load the Julia package which provides the mixed-integer solver, e.g. `using CPLEX`. This solver is typically a mixed-integer linear solver, but if a conic problem has both second-order cones and other nonlinear cones, or if it has PSD cones, then the MIP solver can be an MISOCP solver and Pajarito can put second-order cones in the outer approximation model. -For conic models, the predefined cones are listed in the [MathProgBase](http://mathprogbasejl.readthedocs.io/en/latest/conic.html) documentation. The following conic solvers (**O** means open-source) can be used by Pajarito to solve mixed-integer conic models with any mixture of the corresponding cones: +The continuous conic solver is specified by using the `cont_solver` option, e.g. `PajaritoSolver(cont_solver=MosekSolver())`. For conic models, the predefined cones are listed in the [MathProgBase](http://mathprogbasejl.readthedocs.io/en/latest/conic.html) documentation. The following conic solvers (**O** means open-source) can be used by Pajarito to solve mixed-integer conic models with any mixture of the corresponding cones: | | Second-order | Rotated second-order | Positive semidefinite | Primal exponential | |------------------------|--------------|----------------------|-----------------------|--------------------| -| [CSDP][csdp-url] **O** | | | X | | -| [ECOS][ecos-url] **O** | X | X | | X | -| [SCS][scs-url] **O** | X | X | X | X | -| [Mosek][mosek-url] | X | X | X | | +| [CSDP][csdp-url] **O** | | | X | | +| [ECOS][ecos-url] **O** | X | X | | X | +| [SCS][scs-url] **O** | X | X | X | X | +| [Mosek][mosek-url] | X | X | X | | [csdp-url]: https://github.com/JuliaOpt/CSDP.jl [ecos-url]: https://github.com/JuliaOpt/ECOS.jl [mosek-url]: https://github.com/JuliaOpt/Mosek.jl [scs-url]: https://github.com/JuliaOpt/SCS.jl -MIP and continuous solver parameters must be specified through their corresponding Julia interfaces. For example, to turn off the output of Ipopt solver, use `cont_solver=IpoptSolver(print_level=0)`. +MIP and continuous solver parameters must be specified through their corresponding Julia interfaces. For example, to turn off the output of Mosek solver, use `cont_solver=MosekSolver(LOG=0)`. ## Pajarito solver options -The following options can be passed to `PajaritoSolver()` to modify its behavior (see [solver.jl](https://github.com/mlubin/Pajarito.jl/blob/master/src/solver.jl) for default values; **C** means conic algorithm only): +The following options can be passed to `PajaritoSolver()` to modify its behavior (see [solver.jl](https://github.com/mlubin/Pajarito.jl/blob/master/src/solver.jl) for default values): * `log_level::Int` Verbosity flag: 0 for quiet, 1 for basic solve info, 2 for iteration info, 3 for detailed timing and cuts and solution feasibility info * `timeout::Float64` Time limit for algorithm (in seconds) * `rel_gap::Float64` Relative optimality gap termination condition * `mip_solver_drives::Bool` Let MIP solver manage convergence ("branch and cut") * `mip_solver::MathProgBase.AbstractMathProgSolver` MIP solver (MILP or MISOCP) - * `mip_subopt_solver::MathProgBase.AbstractMathProgSolver` **C** MIP solver for suboptimal solves (with appropriate options already passed) - * `mip_subopt_count::Int` **C** Number of times to use `mip_subopt_solver` between `mip_solver` solves - * `round_mip_sols::Bool` **C** Round integer variable values before solving subproblems - * `use_mip_starts::Bool` **C** Use conic subproblem feasible solutions as MIP warm-starts or heuristic solutions - * `cont_solver::MathProgBase.AbstractMathProgSolver` Continuous solver (conic or nonlinear) - * `solve_relax::Bool` **C** Solve the continuous conic relaxation to add initial subproblem cuts - * `solve_subp::Bool` **C** Solve the continuous conic subproblems to add subproblem cuts - * `dualize_relax::Bool` **C** Solve the conic dual of the continuous conic relaxation - * `dualize_subp::Bool` **C** Solve the conic duals of the continuous conic subproblems - * `soc_disagg::Bool` **C** Disaggregate SOC cones - * `soc_abslift::Bool` **C** Use SOC absolute value lifting - * `soc_in_mip::Bool` **C** Use SOC cones in the MIP model (if `mip_solver` supports MISOCP) - * `sdp_eig::Bool` **C** Use PSD cone eigenvector cuts - * `sdp_soc::Bool` **C** Use PSD cone eigenvector SOC cuts (if `mip_solver` supports MISOCP) - * `init_soc_one::Bool` **C** Use SOC initial L_1 cuts - * `init_soc_inf::Bool` **C** Use SOC initial L_inf cuts - * `init_exp::Bool` **C** Use Exp initial cuts - * `init_sdp_lin::Bool` **C** Use PSD cone initial linear cuts - * `init_sdp_soc::Bool` **C** Use PSD cone initial SOC cuts (if `mip_solver` supports MISOCP) - * `scale_subp_cuts::Bool` **C** Use scaling for subproblem cuts - * `scale_subp_factor::Float64` **C** Fixed multiplicative factor for scaled subproblem cuts - * `viol_cuts_only::Bool` **C** Only add cuts violated by current MIP solution - * `prim_cuts_only::Bool` **C** Add primal cuts, do not add subproblem cuts - * `prim_cuts_always::Bool` **C** Add primal cuts and subproblem cuts - * `prim_cuts_assist::Bool` **C** Add subproblem cuts, and add primal cuts only subproblem cuts cannot be added - * `cut_zero_tol::Float64` **C** Zero tolerance for cut coefficients - * `prim_cut_feas_tol::Float64` **C** Absolute feasibility tolerance used for primal cuts (set equal to feasibility tolerance of `mip_solver`) + * `mip_subopt_solver::MathProgBase.AbstractMathProgSolver` MIP solver for suboptimal solves (with appropriate options already passed) + * `mip_subopt_count::Int` Number of times to use `mip_subopt_solver` between `mip_solver` solves + * `round_mip_sols::Bool` Round integer variable values before solving subproblems + * `use_mip_starts::Bool` Use conic subproblem feasible solutions as MIP warm-starts or heuristic solutions + * `cont_solver::MathProgBase.AbstractMathProgSolver` Continuous conic solver + * `solve_relax::Bool` Solve the continuous conic relaxation to add initial subproblem cuts + * `solve_subp::Bool` Solve the continuous conic subproblems to add subproblem cuts + * `dualize_relax::Bool` Solve the conic dual of the continuous conic relaxation + * `dualize_subp::Bool` Solve the conic duals of the continuous conic subproblems + * `soc_disagg::Bool` Disaggregate SOC cones + * `soc_abslift::Bool` Use SOC absolute value lifting + * `soc_in_mip::Bool` Use SOC cones in the MIP model (if `mip_solver` supports MISOCP) + * `sdp_eig::Bool` Use PSD cone eigenvector cuts + * `sdp_soc::Bool` Use PSD cone eigenvector SOC cuts (if `mip_solver` supports MISOCP) + * `init_soc_one::Bool` Use SOC initial L_1 cuts + * `init_soc_inf::Bool` Use SOC initial L_inf cuts + * `init_exp::Bool` Use Exp initial cuts + * `init_sdp_lin::Bool` Use PSD cone initial linear cuts + * `init_sdp_soc::Bool` Use PSD cone initial SOC cuts (if `mip_solver` supports MISOCP) + * `scale_subp_cuts::Bool` Use scaling for subproblem cuts + * `scale_subp_factor::Float64` Fixed multiplicative factor for scaled subproblem cuts + * `viol_cuts_only::Bool` Only add cuts violated by current MIP solution + * `prim_cuts_only::Bool` Add primal cuts, do not add subproblem cuts + * `prim_cuts_always::Bool` Add primal cuts and subproblem cuts + * `prim_cuts_assist::Bool` Add subproblem cuts, and add primal cuts only subproblem cuts cannot be added + * `cut_zero_tol::Float64` Zero tolerance for cut coefficients + * `prim_cut_feas_tol::Float64` Absolute feasibility tolerance used for primal cuts (set equal to feasibility tolerance of `mip_solver`) **Pajarito may require tuning of parameters to improve convergence.** Note: - * For the conic algorithm, Pajarito usually returns a solution constructed from one of the conic solver's feasible solutions. Since the conic solver is not subject to the same feasibility tolerances as the MIP solver (which should match the absolute feasibility tolerance `prim_cut_feas_tol`), Pajarito's solution will not necessarily satisfy `prim_cut_feas_tol`. + * Pajarito usually returns a solution constructed from one of the conic solver's feasible solutions. Since the conic solver is not subject to the same feasibility tolerances as the MIP solver (which should match the absolute feasibility tolerance `prim_cut_feas_tol`), Pajarito's solution will not necessarily satisfy `prim_cut_feas_tol`. * MIP solver integrality tolerance should typically be tightened, for example to `1e-8`, for improved Pajarito performance. * `viol_cuts_only` defaults to `true` on the MIP-solver-driven algorithm and `false` on the iterative algorithm. diff --git a/REQUIRE b/REQUIRE index f81a0c7..faa8139 100644 --- a/REQUIRE +++ b/REQUIRE @@ -1,5 +1,4 @@ julia 0.6 MathProgBase 0.5 0.8 JuMP 0.15 -ConicNonlinearBridge ConicBenchmarkUtilities diff --git a/examples/cardls.jl b/examples/cardls.jl index cd9d7cc..fe18e4b 100644 --- a/examples/cardls.jl +++ b/examples/cardls.jl @@ -70,8 +70,6 @@ mip_solver = CplexSolver( ) -# Specify mixed-integer conic solver (Pajarito conic algorithm) - # using SCS # conic_solver = SCSSolver(eps=1e-6, max_iters=1000000, verbose=0) @@ -87,20 +85,6 @@ micp_solver = PajaritoSolver( ) -# Specify mixed-integer NLP solver (Pajarito nonlinear algorithm) - -using Ipopt -nlp_solver = IpoptSolver(print_level=0) - -minlp_solver = PajaritoSolver( - mip_solver_drives=mip_solver_drives, - log_level=1, - rel_gap=rel_gap, - mip_solver=mip_solver, - cont_solver=nlp_solver, -) - - #========================================================= Specify/generate data =========================================================# @@ -125,11 +109,8 @@ xB = 4 # Bound on absolute values of estimate variables (|x_j| <= xB) Solve JuMP models =========================================================# -println("\n\n****MIQP model with MINLP solver****\n") -miqp_cardls(m, d, A, b, k, rho, xB, minlp_solver) - -println("\n\n****MIQP model with conic solver****\n") +println("\n\n****MIQP model****\n") miqp_cardls(m, d, A, b, k, rho, xB, micp_solver) -println("\n\n****MISDP model with conic solver****\n") +println("\n\n****MISDP model****\n") misdp_cardls(m, d, A, b, k, rho, micp_solver) diff --git a/examples/runcbf.jl b/examples/runcbf.jl index 636adc1..36a0a8c 100644 --- a/examples/runcbf.jl +++ b/examples/runcbf.jl @@ -26,7 +26,7 @@ solver = PajaritoSolver( mip_solver=mip_solver, cont_solver=cont_solver, log_level=3, - ) +) dat = readcbfdata(ARGS[1]) (c, A, b, con_cones, var_cones, vartypes, sense, objoffset) = cbftompb(dat) diff --git a/src/Pajarito.jl b/src/Pajarito.jl index 1fab34d..c05d3f1 100644 --- a/src/Pajarito.jl +++ b/src/Pajarito.jl @@ -16,10 +16,8 @@ __precompile__() module Pajarito import MathProgBase - using ConicNonlinearBridge include("conic_dual_solver.jl") include("solver.jl") include("conic_algorithm.jl") - include("nonlinear_algorithm.jl") end diff --git a/src/conic_algorithm.jl b/src/conic_algorithm.jl index c47cdc3..07707aa 100644 --- a/src/conic_algorithm.jl +++ b/src/conic_algorithm.jl @@ -11,8 +11,6 @@ This mixed-integer conic programming algorithm is described in: (available online at http://arxiv.org/abs/1511.06710) Model MICP with JuMP.jl conic format or Convex.jl DCP format http://mathprogbasejl.readthedocs.org/en/latest/conic.html -TODO features -- implement warm-starting: use set_best_soln! =========================================================# using JuMP @@ -46,38 +44,38 @@ type PajaritoConicModel <: MathProgBase.AbstractConicModel mip_solver_drives::Bool # Let MIP solver manage convergence ("branch and cut") mip_solver::MathProgBase.AbstractMathProgSolver # MIP solver (MILP or MISOCP) mip_subopt_solver::MathProgBase.AbstractMathProgSolver # MIP solver for suboptimal solves (with appropriate options already passed) - mip_subopt_count::Int # (Conic only) Number of times to use `mip_subopt_solver` between `mip_solver` solves - round_mip_sols::Bool # (Conic only) Round integer variable values before solving subproblems - use_mip_starts::Bool # (Conic only) Use conic subproblem feasible solutions as MIP warm-starts or heuristic solutions - - cont_solver::MathProgBase.AbstractMathProgSolver # Continuous solver (conic or nonlinear) - solve_relax::Bool # (Conic only) Solve the continuous conic relaxation to add initial subproblem cuts - solve_subp::Bool # (Conic only) Solve the continuous conic subproblems to add subproblem cuts - dualize_relax::Bool # (Conic only) Solve the conic dual of the continuous conic relaxation - dualize_subp::Bool # (Conic only) Solve the conic duals of the continuous conic subproblems - - all_disagg::Bool # (Conic only) Disaggregate cuts on the nonpolyhedral cones - soc_disagg::Bool # (Conic only) Disaggregate SOC using extended formulation - soc_abslift::Bool # (Conic only) Use SOC absolute value lifting - soc_in_mip::Bool # (Conic only) Use SOC cones in the MIP model (if `mip_solver` supports MISOCP) - sdp_eig::Bool # (Conic only) Use PSD cone eigenvector cuts - sdp_soc::Bool # (Conic only) Use PSD cone eigenvector SOC cuts (if `mip_solver` supports MISOCP) - init_soc_one::Bool # (Conic only) Use SOC initial L_1 cuts - init_soc_inf::Bool # (Conic only) Use SOC initial L_inf cuts - init_exp::Bool # (Conic only) Use Exp initial cuts - init_sdp_lin::Bool # (Conic only) Use PSD cone initial linear cuts - init_sdp_soc::Bool # (Conic only) Use PSD cone initial SOC cuts (if `mip_solver` supports MISOCP) - - scale_subp_cuts::Bool # (Conic only) Use scaling for subproblem cuts - scale_subp_factor::Float64 # (Conic only) Fixed multiplicative factor for scaled subproblem cuts - scale_subp_up::Bool # (Conic only) Scale up any scaled subproblem cuts that are smaller than the equivalent separation cut - viol_cuts_only::Bool # (Conic only) Only add cuts violated by current MIP solution - sep_cuts_only::Bool # (Conic only) Add primal cuts, do not add subproblem cuts - sep_cuts_always::Bool # (Conic only) Add primal cuts and subproblem cuts - sep_cuts_assist::Bool # (Conic only) Add subproblem cuts, and add primal cuts only subproblem cuts cannot be added - - cut_zero_tol::Float64 # (Conic only) Zero tolerance for cut coefficients - mip_feas_tol::Float64 # (Conic only) Absolute feasibility tolerance used for primal cuts (set equal to feasibility tolerance of `mip_solver`) + mip_subopt_count::Int # Number of times to use `mip_subopt_solver` between `mip_solver` solves + round_mip_sols::Bool # Round integer variable values before solving subproblems + use_mip_starts::Bool # Use conic subproblem feasible solutions as MIP warm-starts or heuristic solutions + + cont_solver::MathProgBase.AbstractMathProgSolver # Continuous conic solver + solve_relax::Bool # Solve the continuous conic relaxation to add initial subproblem cuts + solve_subp::Bool # Solve the continuous conic subproblems to add subproblem cuts + dualize_relax::Bool # Solve the conic dual of the continuous conic relaxation + dualize_subp::Bool # Solve the conic duals of the continuous conic subproblems + + all_disagg::Bool # Disaggregate cuts on the nonpolyhedral cones + soc_disagg::Bool # Disaggregate SOC using extended formulation + soc_abslift::Bool # Use SOC absolute value lifting + soc_in_mip::Bool # Use SOC cones in the MIP model (if `mip_solver` supports MISOCP) + sdp_eig::Bool # Use PSD cone eigenvector cuts + sdp_soc::Bool # Use PSD cone eigenvector SOC cuts (if `mip_solver` supports MISOCP) + init_soc_one::Bool # Use SOC initial L_1 cuts + init_soc_inf::Bool # Use SOC initial L_inf cuts + init_exp::Bool # Use Exp initial cuts + init_sdp_lin::Bool # Use PSD cone initial linear cuts + init_sdp_soc::Bool # Use PSD cone initial SOC cuts (if `mip_solver` supports MISOCP) + + scale_subp_cuts::Bool # Use scaling for subproblem cuts + scale_subp_factor::Float64 # Fixed multiplicative factor for scaled subproblem cuts + scale_subp_up::Bool # Scale up any scaled subproblem cuts that are smaller than the equivalent separation cut + viol_cuts_only::Bool # Only add cuts violated by current MIP solution + sep_cuts_only::Bool # Add primal cuts, do not add subproblem cuts + sep_cuts_always::Bool # Add primal cuts and subproblem cuts + sep_cuts_assist::Bool # Add subproblem cuts, and add primal cuts only subproblem cuts cannot be added + + cut_zero_tol::Float64 # Zero tolerance for cut coefficients + mip_feas_tol::Float64 # Absolute feasibility tolerance used for primal cuts (set equal to feasibility tolerance of `mip_solver`) dump_subproblems::Bool # Save each conic subproblem in conic benchmark format (CBF) at a specified location dump_basename::String # Basename of conic subproblem CBF files: "/path/to/foo" creates files "/path/to/foo_NN.cbf" where "NN" is a counter diff --git a/src/nonlinear_algorithm.jl b/src/nonlinear_algorithm.jl deleted file mode 100644 index 28b64a1..0000000 --- a/src/nonlinear_algorithm.jl +++ /dev/null @@ -1,764 +0,0 @@ -# Copyright 2017, Chris Coey and Miles Lubin -# Copyright 2016, Los Alamos National Laboratory, LANS LLC. -# This Source Code Form is subject to the terms of the Mozilla Public -# License, v. 2.0. If a copy of the MPL was not distributed with this -# file, You can obtain one at http://mozilla.org/MPL/2.0/. - -using JuMP - -type PajaritoNonlinearModel <: MathProgBase.AbstractNonlinearModel - # SOLVER DATA: - verbose::Int # Verbosity level flag - algorithm # Choice of algorithm: "OA" or "BC" - mip_solver # Choice of MILP solver - cont_solver # Choice of Conic solver - opt_tolerance # Relative optimality tolerance - time_limit # Time limit - - solution::Vector{Float64} - status - objval::Float64 - objbound::Float64 - iterations::Int - numVar::Int - numIntVar::Int # Number of integer or binary variables - numConstr::Int - numNLConstr::Int - - A - A_lb - A_ub - lb::Vector{Float64} - ub::Vector{Float64} - l::Vector{Float64} - u::Vector{Float64} - c::Vector{Float64} - vartype::Vector{Symbol} - constrtype::Vector{Symbol} - constrlinear::Vector{Bool} - objsense::Symbol - objlinear::Bool - mip_x - d - - nlp_load_timer - total_time - - # CONSTRUCTOR: - function PajaritoNonlinearModel(verbose,algorithm,mip_solver,cont_solver,opt_tolerance,time_limit) - m = new() - m.solution = Float64[] - m.verbose = verbose - m.algorithm = algorithm - m.mip_solver = mip_solver - m.cont_solver = cont_solver - m.opt_tolerance = opt_tolerance - m.time_limit = time_limit - m.total_time = 0. - return m - end -end - -type InfeasibleNLPEvaluator <: MathProgBase.AbstractNLPEvaluator - d - numConstr::Int - numNLConstr::Int - numVar::Int - constrtype::Vector{Symbol} - constrlinear::Vector{Bool} -end - -function MathProgBase.eval_f(d::InfeasibleNLPEvaluator, x) - retval = 0.0 - # SUM UP THE SLACKS AND RETURN - for i in d.numVar+1:length(x) - retval += x[i] - end - return retval -end - -function MathProgBase.eval_grad_f(d::InfeasibleNLPEvaluator, g, x) - g[:] = [zeros(d.numVar); ones(d.numNLConstr)] -end - -function MathProgBase.eval_g(d::InfeasibleNLPEvaluator, g, x) - MathProgBase.eval_g(d.d, g, x[1:d.numVar]) - k = 1 - for i in 1:d.numConstr - d.constrlinear[i] && continue - if d.constrtype[i] == :(<=) - g[i] -= x[k+d.numVar] - else - g[i] += x[k+d.numVar] - end - k += 1 - end -end - -function MathProgBase.jac_structure(d::InfeasibleNLPEvaluator) - I, J = MathProgBase.jac_structure(d.d) - I_new = copy(I) - J_new = copy(J) - #push indices - k = 1 - for i in 1:(d.numConstr) - d.constrlinear[i] && continue - push!(I_new, i); push!(J_new, k+d.numVar); - k += 1 - end - return I_new, J_new -end - -function MathProgBase.eval_jac_g(d::InfeasibleNLPEvaluator, J, x) - MathProgBase.eval_jac_g(d.d, J, x[1:d.numVar]) - k = length(J) - d.numNLConstr + 1 - for i in 1:d.numConstr - d.constrlinear[i] && continue - if d.constrtype[i] == :(<=) - J[k] = -(1.0) - else - J[k] = 1.0 - end - k += 1 - end -end - -function MathProgBase.eval_hesslag(d::InfeasibleNLPEvaluator, H, x, σ, μ) - MathProgBase.eval_hesslag(d.d, H, x[1:d.numVar], 0.0, μ) -end - -MathProgBase.hesslag_structure(d::InfeasibleNLPEvaluator) = MathProgBase.hesslag_structure(d.d) -MathProgBase.initialize(d::InfeasibleNLPEvaluator, requested_features::Vector{Symbol}) = -MathProgBase.initialize(d.d, requested_features) -MathProgBase.features_available(d::InfeasibleNLPEvaluator) = [:Grad,:Jac,:Hess] - -# function MathProgBase.eval_jac_prod(d::InfeasibleNLPEvaluator, y, x, w) -# jac_I, jac_J = MathProgBase.jac_structure(d) -# jac_V = zeros(length(jac_I)) -# MathProgBase.eval_jac_g(d, jac_V, x) -# varidx_new = [zeros(0) for i in 1:m.numConstr] -# coef_new = [zeros(0) for i in 1:m.numConstr] -# for k in 1:length(jac_I) -# row = jac_I[k] -# push!(varidx_new[row], jac_J[k]); push!(coef_new[row], jac_V[k]) -# end -# -# for i = 1:d.numConstr -# retval = 0.0 -# for j in 1:length(varidx_new[i]) -# retval += coef_new[i][j] * w[varidx_new[i][j]] -# end -# y[i] = retval -# end -# end -# -# DONT NEED IT FOR NOW -# eval_jac_prod_t(d::AbstractNLPEvaluator, y, x, w) -# eval_hesslag_prod(d::AbstractNLPEvaluator, h, x, v, σ, μ) - -MathProgBase.isobjlinear(d::InfeasibleNLPEvaluator) = true -MathProgBase.isobjquadratic(d::InfeasibleNLPEvaluator) = true #MathProgBase.isobjquadratic(d.d) -MathProgBase.isconstrlinear(d::InfeasibleNLPEvaluator, i::Int) = MathProgBase.isconstrlinear(d.d, i) -MathProgBase.obj_expr(d::InfeasibleNLPEvaluator) = MathProgBase.obj_expr(d.d) -MathProgBase.constr_expr(d::InfeasibleNLPEvaluator, i::Int) = MathProgBase.constr_expr(d.d, i) - -function MathProgBase.loadproblem!( - m::PajaritoNonlinearModel, numVar, numConstr, l, u, lb, ub, sense, d) - - if !applicable(MathProgBase.NonlinearModel, m.cont_solver) - error("$(m.cont_solver) is not a nonlinear solver.") - end - - m.numVar = numVar - m.numConstr = numConstr - m.lb = lb - m.ub = ub - m.l = l - m.u = u - m.objsense = sense - m.d = d - m.vartype = fill(:Cont,numVar) - - MathProgBase.initialize(d, [:Grad,:Jac,:Hess]) - - m.constrtype = Array{Symbol}(numConstr) - for i = 1:numConstr - if lb[i] > -Inf && ub[i] < Inf - m.constrtype[i] = :(==) - elseif lb[i] > -Inf - m.constrtype[i] = :(>=) - else - m.constrtype[i] = :(<=) - end - end - m.solution = fill(NaN, m.numVar) -end - -function OAprintLevel(iter, mip_objval, conic_objval, optimality_gap, best_objval, primal_infeasibility, OA_infeasibility) - - if abs(conic_objval) == Inf || isnan(conic_objval) - conic_objval_str = @sprintf "%s" " " - else - conic_objval_str = @sprintf "%+.7e" conic_objval - end - if abs(optimality_gap) == Inf || isnan(optimality_gap) - optimality_gap_str = @sprintf "%s" " " - else - optimality_gap_str = @sprintf "%+.7e" optimality_gap - end - if abs(best_objval) == Inf || isnan(best_objval) - best_objval_str = @sprintf "%s" " " - else - best_objval_str = @sprintf "%+.7e" best_objval - end - - @printf "%9d %+.7e %s %s %s %+.7e %+.7e\n" iter mip_objval conic_objval_str optimality_gap_str best_objval_str primal_infeasibility OA_infeasibility -end - -function populatelinearmatrix(m::PajaritoNonlinearModel) - # set up map of linear rows - constrlinear = Array{Bool}(m.numConstr) - numlinear = 0 - constraint_to_linear = fill(-1,m.numConstr) - for i = 1:m.numConstr - constrlinear[i] = MathProgBase.isconstrlinear(m.d, i) - if constrlinear[i] - numlinear += 1 - constraint_to_linear[i] = numlinear - end - end - m.numNLConstr = m.numConstr - numlinear - - # extract sparse jacobian structure - jac_I, jac_J = MathProgBase.jac_structure(m.d) - - # evaluate jacobian at x = 0 - c = zeros(m.numVar) - x = m.solution - jac_V = zeros(length(jac_I)) - MathProgBase.eval_jac_g(m.d, jac_V, x) - MathProgBase.eval_grad_f(m.d, c, x) - m.objlinear = MathProgBase.isobjlinear(m.d) - if m.objlinear - (m.verbose > 0) && println("Objective function is linear") - m.c = c - else - (m.verbose > 0) && println("Objective function is nonlinear") - m.c = zeros(m.numVar) - end - - # Build up sparse matrix for linear constraints - A_I = Int[] - A_J = Int[] - A_V = Float64[] - - for k in 1:length(jac_I) - row = jac_I[k] - if !constrlinear[row] - continue - end - row = constraint_to_linear[row] - push!(A_I,row); push!(A_J, jac_J[k]); push!(A_V, jac_V[k]) - end - - m.A = sparse(A_I, A_J, A_V, numlinear, m.numVar) - - # g(x) might have a constant, i.e., a'x + b - # let's find b - constraint_value = zeros(m.numConstr) - MathProgBase.eval_g(m.d, constraint_value, x) - b = constraint_value[constrlinear] - m.A * x - # so linear constraints are of the form lb ≤ a'x + b ≤ ub - - # set up A_lb and A_ub vectors - m.A_lb = m.lb[constrlinear] - b - m.A_ub = m.ub[constrlinear] - b - - # Now we have linear parts - m.constrlinear = constrlinear -end - -function addCuttingPlanes!(m::PajaritoNonlinearModel, mip_model, separator, jac_I, jac_J, jac_V, grad_f, cb, mip_solution) - max_violation = -1e+5 - # EVALUATE g and jac_g AT MIP SOLUTION THAT IS INFEASIBLE - g = zeros(m.numConstr) - MathProgBase.eval_g(m.d, g, separator[1:m.numVar]) - MathProgBase.eval_jac_g(m.d, jac_V, separator[1:m.numVar]) - - # create rows corresponding to constraints in sparse format - - varidx_new = [zeros(Int, 0) for i in 1:m.numConstr] - coef_new = [zeros(0) for i in 1:m.numConstr] - - for k in 1:length(jac_I) - row = jac_I[k] - push!(varidx_new[row], jac_J[k]); push!(coef_new[row], jac_V[k]) - end - - # CREATE CONSTRAINT CUTS - for i in 1:m.numConstr - if m.constrtype[i] == :(<=) - val = g[i] - m.ub[i] - else - val = m.lb[i] - g[i] - end - lin = m.constrlinear[i] - (m.verbose > 2) && println("Constraint $i value $val linear $lin") - if !(lin) #&& (val > 1e-4) # m.ub[i] is in the constraint somehow - # CREATE SUPPORTING HYPERPLANE - (m.verbose > 2) && println("Create supporting hyperplane for constraint $i") - local new_rhs::Float64 - if m.constrtype[i] == :(<=) - new_rhs = m.ub[i] - g[i] - else - new_rhs = m.lb[i] - g[i] - end - for j = 1:length(varidx_new[i]) - new_rhs += coef_new[i][j] * separator[Int(varidx_new[i][j])] - end - (m.verbose > 2) && println("varidx $(varidx_new[i])") - (m.verbose > 2) && println("coef $(coef_new[i])") - (m.verbose > 2) && println("rhs $new_rhs") - if m.constrtype[i] == :(<=) - if cb != [] - @lazyconstraint(cb, dot(coef_new[i], m.mip_x[varidx_new[i]]) <= new_rhs) - else - @constraint(mip_model, dot(coef_new[i], m.mip_x[varidx_new[i]]) <= new_rhs) - end - viol = vecdot(coef_new[i], mip_solution[varidx_new[i]]) - new_rhs - if viol > max_violation - max_violation = viol - end - #MathProgBase.addconstr!(mip_model, varidx_new[i], coef_new[i], -Inf, new_rhs) - else - if cb != [] - @lazyconstraint(cb, dot(coef_new[i], m.mip_x[varidx_new[i]]) >= new_rhs) - else - @constraint(mip_model, dot(coef_new[i], m.mip_x[varidx_new[i]]) >= new_rhs) - end - viol = new_rhs - vecdot(coef_new[i], mip_solution[varidx_new[i]]) - if viol > max_violation - max_violation = viol - end - #MathProgBase.addconstr!(mip_model, varidx_new[i], coef_new[i], new_rhs, Inf) - end - end - end - # CREATE OBJECTIVE CUTS - if !(m.objlinear) - (m.verbose > 2) && println("Create supporting hyperplane for objective f(x) <= t") - f = MathProgBase.eval_f(m.d, separator[1:m.numVar]) - MathProgBase.eval_grad_f(m.d, grad_f, separator[1:m.numVar]) - if m.objsense == :Max - f = -f - grad_f = -grad_f - end - new_rhs = -f - varidx = zeros(Int, m.numVar+1) - for j = 1:m.numVar - varidx[j] = j - new_rhs += grad_f[j] * separator[j] - end - varidx[m.numVar+1] = m.numVar+1 - grad_f[m.numVar+1] = -(1.0) - (m.verbose > 2) && println("varidx $(varidx)") - (m.verbose > 2) && println("coef $(grad_f)") - (m.verbose > 2) && println("rhs $new_rhs") - if cb != [] - @lazyconstraint(cb, dot(grad_f, m.mip_x[varidx]) <= new_rhs) - else - @constraint(mip_model, dot(grad_f, m.mip_x[varidx]) <= new_rhs) - end - viol = vecdot(grad_f, mip_solution[varidx]) - new_rhs - if viol > max_violation - max_violation = viol - end - #MathProgBase.addconstr!(mip_model, varidx, grad_f, -Inf, new_rhs) - end - - return max_violation -end - -function loadMIPModel(m::PajaritoNonlinearModel, mip_model) - lb = [m.l; -1e6] - ub = [m.u; 1e6] - @variable(mip_model, lb[i] <= x[i=1:m.numVar+1] <= ub[i]) - numIntVar = 0 - for i = 1:m.numVar - setcategory(x[i], m.vartype[i]) - if m.vartype[i] == :Int || m.vartype[i] == :Bin - numIntVar += 1 - end - end - if numIntVar == 0 - error("No variables of type integer or binary; call the conic continuous solver directly for pure continuous problems") - end - setcategory(x[m.numVar+1], :Cont) - for i = 1:m.numConstr-m.numNLConstr - if m.A_lb[i] > -Inf && m.A_ub[i] < Inf - if m.A_lb[i] == m.A_ub[i] - @constraint(mip_model, m.A[i:i,:]*x[1:m.numVar] .== m.A_lb[i]) - else - @constraint(mip_model, m.A[i:i,:]*x[1:m.numVar] .>= m.A_lb[i]) - @constraint(mip_model, m.A[i:i,:]*x[1:m.numVar] .<= m.A_ub[i]) - end - elseif m.A_lb[i] > -Inf - @constraint(mip_model, m.A[i:i,:]*x[1:m.numVar] .>= m.A_lb[i]) - else - @constraint(mip_model, m.A[i:i,:]*x[1:m.numVar] .<= m.A_ub[i]) - end - end - c_new = [m.objsense == :Max ? -m.c : m.c; m.objlinear ? 0.0 : 1.0] - @objective(mip_model, Min, dot(c_new, x)) - - m.mip_x = x - m.numIntVar = numIntVar - #= - mip_model = MathProgBase.LinearQuadraticModel(m.mip_solver) - MathProgBase.loadproblem!(mip_model, - [m.A spzeros(size(m.A,1),1)], - [m.l; -1e4], - [m.u; Inf], - [(m.objsense == :Max)? -m.c : m.c; m.objlinear? 0.0 : 1.0], m.A_lb, m.A_ub, :Min) - MathProgBase.setvartype!(mip_model, [m.vartype; :Cont]) - =# -end - -function checkInfeasibility(m::PajaritoNonlinearModel, solution) - g = zeros(m.numConstr) - MathProgBase.eval_g(m.d, g, solution[1:m.numVar]) - max_infeas = 0.0 - for i = 1:m.numConstr - max_infeas = max(max_infeas, g[i]-m.ub[i], m.lb[i]-g[i]) - end - for i = 1:m.numVar - max_infeas = max(max_infeas, solution[i]-m.u[i], m.l[i]-solution[i]) - end - - return max_infeas -end - -function compareIntegerSolutions(m::PajaritoNonlinearModel, sol1, sol2) - int_ind = filter(i->m.vartype[i] == :Int || m.vartype[i] == :Bin, 1:m.numVar) - return round.(sol1[int_ind]) == round.(sol2[int_ind]) -end - -function MathProgBase.optimize!(m::PajaritoNonlinearModel) - start = time() - - cputime_nlp = 0.0 - cputime_mip = 0.0 - - m.nlp_load_timer = 0.0 - - # if we haven't gotten a starting point, - # (e.g., if acting as MIQP solver) assume zero is okay - if any(isnan,m.solution) - m.solution = zeros(length(m.solution)) - end - - populatelinearmatrix(m) - # solve it - # MathProgBase.optimize!(nlp_model) - - # pull out solution values - #= m.status = MathProgBase.status(nlp_model) - m.objval = MathProgBase.getobjval(nlp_model) - m.solution = MathProgBase.getsolution(nlp_model) =# - - # MODIFICATION new objective t >= f(x) - # set up cplex to solve the mixed integer linear problem - mip_model = Model(solver=m.mip_solver) - loadMIPModel(m, mip_model) - - for i in 1:m.numConstr - if !(m.constrlinear[i]) && m.constrtype[i] == :(==) - error("Nonlinear equality or two-sided constraints not accepted.") - end - end - - jac_I, jac_J = MathProgBase.jac_structure(m.d) - jac_V = zeros(length(jac_I)) - grad_f = zeros(m.numVar+1) - - ini_nlp_model = MathProgBase.NonlinearModel(m.cont_solver) - start_load = time() - MathProgBase.loadproblem!(ini_nlp_model, - m.numVar, m.numConstr, m.l, m.u, m.lb, m.ub, m.objsense, m.d) - m.nlp_load_timer += time() - start_load - - # pass in starting point - #MathProgBase.setwarmstart!(nlp_model, m.solution) - MathProgBase.setwarmstart!(ini_nlp_model, m.solution[1:m.numVar]) - - start_nlp = time() - MathProgBase.optimize!(ini_nlp_model) - cputime_nlp += time() - start_nlp - - ini_nlp_status = MathProgBase.status(ini_nlp_model) - if ini_nlp_status == :Optimal || ini_nlp_status == :Suboptimal - if m.numIntVar == 0 - m.solution = MathProgBase.getsolution(ini_nlp_model) - m.objval = MathProgBase.getobjval(ini_nlp_model) - m.status = ini_nlp_status - return - end - separator = MathProgBase.getsolution(ini_nlp_model) - addCuttingPlanes!(m, mip_model, separator, jac_I, jac_J, jac_V, grad_f, [], zeros(m.numVar+1)) - elseif ini_nlp_status == :Infeasible - (m.verbose > 0) && println("Initial NLP Relaxation Infeasible.") - m.status = :Infeasible - return - # TODO Figure out the conditions for this to hold! - elseif ini_nlp_status == :Unbounded - warn("Initial NLP Relaxation Unbounded.") - m.status = :Unbounded - return - else - warn("NLP Solver Failure.") - m.status = :Error - return - end - ini_nlp_objval = MathProgBase.getobjval(ini_nlp_model) - - (m.verbose > 0) && println("\nPajarito started...\n") - (m.verbose > 0) && println("MINLP algorithm $(m.algorithm) is chosen.") - (m.verbose > 0) && println("MINLP has $(m.numVar) variables, $(m.numConstr - m.numNLConstr) linear constraints, $(m.numNLConstr) nonlinear constraints.") - (m.verbose > 0) && @printf "Initial objective = %13.5f.\n\n" ini_nlp_objval - - m.status = :UserLimit - m.objval = Inf - iter = 0 - prev_mip_solution = fill(NaN,m.numVar) - cut_added = false - - nlp_status = :Infeasible - nlp_solution = zeros(m.numVar) - - function nonlinearcallback(cb) - if cb != [] - mip_objval = -Inf #MathProgBase.cbgetobj(cb) - mip_solution = MathProgBase.cbgetmipsolution(cb)[1:m.numVar+1] - else - m.objbound = mip_objval = getobjbound(mip_model) - mip_solution = getvalue(m.mip_x) - end - #MathProgBase.getsolution(internalmodel(mip_model)) - (m.verbose > 2) && println("MIP Solution: $mip_solution") - # solve NLP model for the MIP solution - new_u = m.u - new_l = m.l - for i in 1:m.numVar - if m.vartype[i] == :Int || m.vartype[i] == :Bin - new_u[i] = mip_solution[i] - new_l[i] = mip_solution[i] - end - end - - # set up ipopt to solve continuous relaxation - nlp_model = MathProgBase.NonlinearModel(m.cont_solver) - start_load = time() - MathProgBase.loadproblem!(nlp_model, - m.numVar, m.numConstr, new_l, new_u, m.lb, m.ub, m.objsense, m.d) - m.nlp_load_timer += time() - start_load - - # pass in starting point - #MathProgBase.setwarmstart!(nlp_model, m.solution) - MathProgBase.setwarmstart!(nlp_model, mip_solution[1:m.numVar]) - - l_inf = [new_l;zeros(m.numNLConstr)] - u_inf = [new_u;Inf*ones(m.numNLConstr)] - - d_inf = InfeasibleNLPEvaluator(m.d, m.numConstr, m.numNLConstr, m.numVar, m.constrtype, m.constrlinear) - inf_model = MathProgBase.NonlinearModel(m.cont_solver) - start_load = time() - MathProgBase.loadproblem!(inf_model, - m.numVar+m.numNLConstr, m.numConstr, l_inf, u_inf, m.lb, m.ub, :Min, d_inf) - m.nlp_load_timer += time() - start_load - - #MathProgBase.setvarUB!(nlp_model, new_u) - #MathProgBase.setvarLB!(nlp_model, new_l) - # optimize the NLP problem - start_nlp = time() - MathProgBase.optimize!(nlp_model) - cputime_nlp += time() - start_nlp - nlp_status = MathProgBase.status(nlp_model) - nlp_objval = (m.objsense == :Max ? -Inf : Inf) - #separator::Vector{Float64} - if nlp_status == :Optimal - (m.verbose > 2) && println("NLP Solved") - nlp_objval = MathProgBase.getobjval(nlp_model) - nlp_solution = MathProgBase.getsolution(nlp_model) - separator = copy(nlp_solution) - (m.verbose > 2) && println("NLP Solution: $separator") - - # KEEP TRACK OF BEST KNOWN INTEGER FEASIBLE SOLUTION - if m.objsense == :Max - if nlp_objval > -m.objval - m.objval = -nlp_objval - m.solution = separator[1:m.numVar] - end - else - if nlp_objval < m.objval - m.objval = nlp_objval - m.solution = separator[1:m.numVar] - end - end - - else - # Create the warm start solution for inf model - inf_initial_solution = zeros(m.numVar + m.numNLConstr); - inf_initial_solution[1:m.numVar] = mip_solution[1:m.numVar] - g = zeros(m.numConstr) - MathProgBase.eval_g(m.d, g, inf_initial_solution[1:m.numVar]) - k = 1 - for i in 1:m.numConstr - if !m.constrlinear[i] - if m.constrtype[i] == :(<=) - val = g[i] - m.ub[i] - else - val = m.lb[i] - g[i] - end - if val > 0 - # Because the sign of the slack changes if the constraint - # direction change - inf_initial_solution[m.numVar + k] = val - else - inf_initial_solution[m.numVar + k] = 0.0 - end - k += 1 - end - end - - MathProgBase.setwarmstart!(inf_model, inf_initial_solution) - (m.verbose > 2) && println("NLP Infeasible") - start_nlp = time() - MathProgBase.optimize!(inf_model) - cputime_nlp += time() - start_nlp - inf_model_status = MathProgBase.status(inf_model) - if inf_model_status != :Optimal && inf_model_status != :Suboptimal - warn("NLP Solver Failure.") - m.status = :Error - return - end - (m.verbose > 2) && println("INF NLP Solved") - separator = MathProgBase.getsolution(inf_model) - (m.verbose > 2) && println("INF NLP Solution: $separator") - end - # add supporting hyperplanes - cycle_indicator = (m.algorithm == "OA" ? compareIntegerSolutions(m, prev_mip_solution, mip_solution) : false) - # m.objval and mip_objval are always in minimization sense - optimality_gap = m.objval - mip_objval - primal_infeasibility = checkInfeasibility(m, mip_solution) - OA_infeasibility = 0.0 - if (optimality_gap > (abs(mip_objval) + 1e-5)*m.opt_tolerance && !cycle_indicator) || cb != [] - OA_infeasibility = addCuttingPlanes!(m, mip_model, separator, jac_I, jac_J, jac_V, grad_f, cb, mip_solution) - #(m.cut_switch > 0) && addCuttingPlanes!(m, mip_model, mip_solution, jac_I, jac_J, jac_V, grad_f, cb, mip_solution) - cut_added = true - else - if optimality_gap < (abs(mip_objval) + 1e-5)*m.opt_tolerance - (m.verbose > 1) && println("MINLP Solved") - m.status = :Optimal - m.iterations = iter - (m.verbose > 1) && println("Number of OA iterations: $iter") - else - @assert cycle_indicator - m.status = :Suboptimal - end - end - fixsense = (m.objsense == :Max) ? -1 : 1 - (m.verbose > 0) && (m.algorithm == "OA") && OAprintLevel(iter, fixsense*mip_objval, nlp_objval, optimality_gap, fixsense*m.objval, primal_infeasibility, OA_infeasibility) - (cycle_indicator && m.status != :Optimal) && warn("Mixed-integer cycling detected, terminating Pajarito...") - end - - function heuristiccallback(cb) - if nlp_status == :Optimal - for i = 1:m.numVar - setsolutionvalue(cb, m.mip_x[i], nlp_solution[i]) - end - addsolution(cb) - end - end - - if m.algorithm == "BC" - addlazycallback(mip_model, nonlinearcallback) - addheuristiccallback(mip_model, heuristiccallback) - m.status = solve(mip_model, suppress_warnings=true) - m.objbound = getobjbound(mip_model) - elseif m.algorithm == "OA" - (m.verbose > 0) && println("Iteration MIP Objective NLP Objective Optimality Gap Best Solution Primal Inf. OA Inf.") - # Check if we've been provided with an initial feasible solution - if all(isfinite,m.solution) && checkInfeasibility(m, m.solution) < 1e-5 - if m.objsense == :Max - m.objval = -MathProgBase.eval_f(m.d, m.solution) - else - m.objval = MathProgBase.eval_f(m.d, m.solution) - end - end - while (time() - start) < m.time_limit - flush(STDOUT) - cut_added = false - # WARMSTART MIP FROM UPPER BOUND - if !any(isnan,m.solution) && !isempty(m.solution) - if applicable(MathProgBase.setwarmstart!, internalmodel(mip_model), m.solution) - # Extend solution with the objective variable - MathProgBase.setwarmstart!(internalmodel(mip_model), [m.solution; m.objval]) - end - end - # solve MIP model - start_mip = time() - mip_status = solve(mip_model, suppress_warnings=true) - cputime_mip += time() - start_mip - #mip_objval = Inf - #mip_solution = zeros(m.numVar+1) - if mip_status == :Infeasible || mip_status == :InfeasibleOrUnbounded - (m.verbose > 1) && println("MIP Infeasible") - m.status = :Infeasible - break - else - (m.verbose > 1) && println("MIP Status: $mip_status") - end - mip_solution = getvalue(m.mip_x) - - nonlinearcallback([]) - - if cut_added == false - break - end - - prev_mip_solution = mip_solution - iter += 1 - end - else - error("Unspecified algorithm.") - end - - if m.objsense == :Max - m.objval = -m.objval - m.objbound = -m.objbound - end - - m.total_time = time()-start - - if m.verbose > 0 - println("\nPajarito finished...\n") - @printf "Status = %13s.\n" m.status - (m.status == :Optimal) && @printf "Optimum objective = %13.5f.\n" m.objval - (m.algorithm == "OA") && @printf "Iterations = %13d.\n" iter - @printf "Total time = %13.5f sec.\n" m.total_time - @printf "MIP total time = %13.5f sec.\n" cputime_mip - @printf "NLP total time = %13.5f sec.\n" cputime_nlp - @printf "Subprob load time = %13.5f sec.\n" m.nlp_load_timer - end -end - -MathProgBase.setwarmstart!(m::PajaritoNonlinearModel, x) = (m.solution = x) -MathProgBase.setvartype!(m::PajaritoNonlinearModel, v::Vector{Symbol}) = (m.vartype = v) -#MathProgBase.setvarUB!(m::IpoptMathProgModel, v::Vector{Float64}) = (m.u = v) -#MathProgBase.setvarLB!(m::IpoptMathProgModel, v::Vector{Float64}) = (m.l = v) - -MathProgBase.status(m::PajaritoNonlinearModel) = m.status -MathProgBase.getobjval(m::PajaritoNonlinearModel) = m.objval -MathProgBase.getobjbound(m::PajaritoNonlinearModel) = m.objbound -MathProgBase.getsolution(m::PajaritoNonlinearModel) = m.solution -MathProgBase.getsolvetime(m::PajaritoNonlinearModel) = m.total_time diff --git a/src/solver.jl b/src/solver.jl index 6a0e94f..75f7411 100644 --- a/src/solver.jl +++ b/src/solver.jl @@ -5,16 +5,13 @@ # file, You can obtain one at http://mozilla.org/MPL/2.0/. #========================================================= -This file implements the default PajaritoSolver + Pajarito solver object =========================================================# export PajaritoSolver - # Dummy solver -type UnsetSolver <: MathProgBase.AbstractMathProgSolver -end - +type UnsetSolver <: MathProgBase.AbstractMathProgSolver end # Pajarito solver type PajaritoSolver <: MathProgBase.AbstractMathProgSolver @@ -25,44 +22,43 @@ type PajaritoSolver <: MathProgBase.AbstractMathProgSolver mip_solver_drives::Bool # Let MIP solver manage convergence ("branch and cut") mip_solver::MathProgBase.AbstractMathProgSolver # MIP solver (MILP or MISOCP) mip_subopt_solver::MathProgBase.AbstractMathProgSolver # MIP solver for suboptimal solves (with appropriate options already passed) - mip_subopt_count::Int # (Conic only) Number of times to use `mip_subopt_solver` between `mip_solver` solves - round_mip_sols::Bool # (Conic only) Round integer variable values before solving subproblems - use_mip_starts::Bool # (Conic only) Use conic subproblem feasible solutions as MIP warm-starts or heuristic solutions - - cont_solver::MathProgBase.AbstractMathProgSolver # Continuous solver (conic or nonlinear) - solve_relax::Bool # (Conic only) Solve the continuous conic relaxation to add initial subproblem cuts - solve_subp::Bool # (Conic only) Solve the continuous conic subproblems to add subproblem cuts - dualize_relax::Bool # (Conic only) Solve the conic dual of the continuous conic relaxation - dualize_subp::Bool # (Conic only) Solve the conic duals of the continuous conic subproblems - - all_disagg::Bool # (Conic only) Disaggregate cuts on the nonpolyhedral cones - soc_disagg::Bool # (Conic only) Disaggregate SOC cones - soc_abslift::Bool # (Conic only) Use SOC absolute value lifting - soc_in_mip::Bool # (Conic only) Use SOC cones in the MIP model (if `mip_solver` supports MISOCP) - sdp_eig::Bool # (Conic only) Use PSD cone eigenvector cuts - sdp_soc::Bool # (Conic only) Use PSD cone eigenvector SOC cuts (if `mip_solver` supports MISOCP) - init_soc_one::Bool # (Conic only) Use SOC initial L_1 cuts - init_soc_inf::Bool # (Conic only) Use SOC initial L_inf cuts - init_exp::Bool # (Conic only) Use Exp initial cuts - init_sdp_lin::Bool # (Conic only) Use PSD cone initial linear cuts - init_sdp_soc::Bool # (Conic only) Use PSD cone initial SOC cuts (if `mip_solver` supports MISOCP) - - scale_subp_cuts::Bool # (Conic only) Use scaling for subproblem cuts - scale_subp_factor::Float64 # (Conic only) Fixed multiplicative factor for scaled subproblem cuts - scale_subp_up::Bool # (Conic only) Scale up any scaled subproblem cuts that are smaller than the equivalent separation cut - viol_cuts_only::Bool # (Conic only) Only add cuts violated by current MIP solution - prim_cuts_only::Bool # (Conic only) Add primal cuts, do not add subproblem cuts - prim_cuts_always::Bool # (Conic only) Add primal cuts and subproblem cuts - prim_cuts_assist::Bool # (Conic only) Add subproblem cuts, and add primal cuts only subproblem cuts cannot be added - - cut_zero_tol::Float64 # (Conic only) Zero tolerance for cut coefficients - prim_cut_feas_tol::Float64 # (Conic only) Absolute feasibility tolerance used for primal cuts (set equal to feasibility tolerance of `mip_solver`) + mip_subopt_count::Int # Number of times to use `mip_subopt_solver` between `mip_solver` solves + round_mip_sols::Bool # Round integer variable values before solving subproblems + use_mip_starts::Bool # Use conic subproblem feasible solutions as MIP warm-starts or heuristic solutions + + cont_solver::MathProgBase.AbstractMathProgSolver # Continuous conic solver + solve_relax::Bool # Solve the continuous conic relaxation to add initial subproblem cuts + solve_subp::Bool # Solve the continuous conic subproblems to add subproblem cuts + dualize_relax::Bool # Solve the conic dual of the continuous conic relaxation + dualize_subp::Bool # Solve the conic duals of the continuous conic subproblems + + all_disagg::Bool # Disaggregate cuts on the nonpolyhedral cones + soc_disagg::Bool # Disaggregate SOC cones + soc_abslift::Bool # Use SOC absolute value lifting + soc_in_mip::Bool # Use SOC cones in the MIP model (if `mip_solver` supports MISOCP) + sdp_eig::Bool # Use PSD cone eigenvector cuts + sdp_soc::Bool # Use PSD cone eigenvector SOC cuts (if `mip_solver` supports MISOCP) + init_soc_one::Bool # Use SOC initial L_1 cuts + init_soc_inf::Bool # Use SOC initial L_inf cuts + init_exp::Bool # Use Exp initial cuts + init_sdp_lin::Bool # Use PSD cone initial linear cuts + init_sdp_soc::Bool # Use PSD cone initial SOC cuts (if `mip_solver` supports MISOCP) + + scale_subp_cuts::Bool # Use scaling for subproblem cuts + scale_subp_factor::Float64 # Fixed multiplicative factor for scaled subproblem cuts + scale_subp_up::Bool # Scale up any scaled subproblem cuts that are smaller than the equivalent separation cut + viol_cuts_only::Bool # Only add cuts violated by current MIP solution + prim_cuts_only::Bool # Add primal cuts, do not add subproblem cuts + prim_cuts_always::Bool # Add primal cuts and subproblem cuts + prim_cuts_assist::Bool # Add subproblem cuts, and add primal cuts only subproblem cuts cannot be added + + cut_zero_tol::Float64 # Zero tolerance for cut coefficients + prim_cut_feas_tol::Float64 # Absolute feasibility tolerance used for primal cuts (set equal to feasibility tolerance of `mip_solver`) dump_subproblems::Bool # Save each conic subproblem in conic benchmark format (CBF) at a specified location dump_basename::String # Basename of conic subproblem CBF files: "/path/to/foo" creates files "/path/to/foo_NN.cbf" where "NN" is a counter end - function PajaritoSolver(; log_level = 1, timeout = Inf, @@ -109,6 +105,10 @@ function PajaritoSolver(; dump_basename = nothing, ) + if (cont_solver != UnsetSolver()) && !applicable(MathProgBase.ConicModel, cont_solver) + error("Continuous solver (cont_solver) specified is not a conic solver; if your continuous solver is a derivative-based NLP solver, try Pavito solver (Pajarito's MINLP functionality was moved to the Pavito solver package)\n") + end + if mip_solver == UnsetSolver() error("No MIP solver specified (set mip_solver)\n") end @@ -130,104 +130,71 @@ function PajaritoSolver(; PajaritoSolver(log_level, timeout, rel_gap, mip_solver_drives, deepcopy(mip_solver), deepcopy(mip_subopt_solver), mip_subopt_count, round_mip_sols, use_mip_starts, deepcopy(cont_solver), solve_relax, solve_subp, dualize_relax, dualize_subp, all_disagg, soc_disagg, soc_abslift, soc_in_mip, sdp_eig, sdp_soc, init_soc_one, init_soc_inf, init_exp, init_sdp_lin, init_sdp_soc, scale_subp_cuts, scale_subp_factor, scale_subp_up, viol_cuts_only, prim_cuts_only, prim_cuts_always, prim_cuts_assist, cut_zero_tol, prim_cut_feas_tol, dump_subproblems, dump_basename) end +# Cannot use Pajarito on an NLP model +MathProgBase.NonlinearModel(s::PajaritoSolver) = error("Pajarito solver cannot be used for NLP models (Pajarito's MINLP functionality was moved to the Pavito solver package)\n") -# Create Pajarito conic model: can solve with either conic algorithm or nonlinear algorithm wrapped with ConicNonlinearBridge +# Create Pajarito conic model function MathProgBase.ConicModel(s::PajaritoSolver) - if applicable(MathProgBase.ConicModel, s.cont_solver) || (s.cont_solver == UnsetSolver()) - if (s.solve_relax || s.solve_subp) && (s.cont_solver == UnsetSolver()) - error("Using conic relaxation (solve_relax) or subproblem solves (solve_subp), but no continuous solver specified (set cont_solver)\n") - end - - if s.soc_in_mip || s.init_sdp_soc || s.sdp_soc - # If using MISOCP outer approximation, check MIP solver handles MISOCP - if !(:SOC in MathProgBase.supportedcones(s.mip_solver)) - error("Using SOC constraints in the MIP model (soc_in_mip or init_sdp_soc or sdp_soc), but MIP solver (mip_solver) specified does not support MISOCP\n") - end - end - - if !s.all_disagg && (s.soc_disagg || s.sdp_eig || s.sdp_soc) - error("Cannot use SOC extended formulation (soc_disagg) or SDP cut disaggregation (sdp_eig) or SOC cuts for SDP cones (sdp_soc) when not also disaggregating all nonpolyhedral cone cuts (all_disagg)\n") - end - - if (s.mip_subopt_count > 0) && (s.mip_subopt_solver == UnsetSolver()) - error("Using suboptimal solves (mip_subopt_count > 0), but no suboptimal MIP solver specified (set mip_subopt_solver)\n") - end - - if s.init_soc_one && !s.soc_disagg && !s.soc_abslift - error("Cannot use SOC initial L_1 cuts (init_soc_one) if both SOC disaggregation (soc_disagg) and SOC absvalue lifting (soc_abslift) are not used\n") - end - - if s.sdp_soc && s.mip_solver_drives - warn("In the MIP-solver-driven algorithm, SOC cuts for SDP cones (sdp_soc) cannot be added from subproblems or primal solutions, but they will be added from the conic relaxation\n") - end - - if !s.solve_subp - s.prim_cuts_only = true - s.use_mip_starts = false - s.round_mip_sols = false - end - if s.prim_cuts_only - s.prim_cuts_always = true - end - if s.prim_cuts_always - s.prim_cuts_assist = true - end + if (s.solve_relax || s.solve_subp) && (s.cont_solver == UnsetSolver()) + error("Using conic relaxation (solve_relax) or subproblem solves (solve_subp), but no continuous solver specified (set cont_solver)\n") + end - if !s.all_disagg && s.prim_cuts_assist - error("Cannot use primal cuts when not disaggregating all nonpolyhedral cone cuts (all_disagg)\n") + if s.soc_in_mip || s.init_sdp_soc || s.sdp_soc + # If using MISOCP outer approximation, check MIP solver handles MISOCP + if !(:SOC in MathProgBase.supportedcones(s.mip_solver)) + error("Using SOC constraints in the MIP model (soc_in_mip or init_sdp_soc or sdp_soc), but MIP solver (mip_solver) specified does not support MISOCP\n") end - - return PajaritoConicModel(s.log_level, s.timeout, s.rel_gap, s.mip_solver_drives, s.mip_solver, s.mip_subopt_solver, s.mip_subopt_count, s.round_mip_sols, s.use_mip_starts, s.cont_solver, s.solve_relax, s.solve_subp, s.dualize_relax, s.dualize_subp, s.all_disagg, s.soc_disagg, s.soc_abslift, s.soc_in_mip, s.sdp_eig, s.sdp_soc, s.init_soc_one, s.init_soc_inf, s.init_exp, s.init_sdp_lin, s.init_sdp_soc, s.scale_subp_cuts, s.scale_subp_factor, s.scale_subp_up, s.viol_cuts_only, s.prim_cuts_only, s.prim_cuts_always, s.prim_cuts_assist, s.cut_zero_tol, s.prim_cut_feas_tol, s.dump_subproblems, s.dump_basename) - elseif applicable(MathProgBase.NonlinearModel, s.cont_solver) - return MathProgBase.ConicModel(ConicNonlinearBridge.ConicNLPWrapper(nlp_solver=s)) - else - error("Continuous solver (cont_solver) specified is not a conic or NLP solver recognized by MathProgBase\n") end -end + if !s.all_disagg && (s.soc_disagg || s.sdp_eig || s.sdp_soc) + error("Cannot use SOC extended formulation (soc_disagg) or SDP cut disaggregation (sdp_eig) or SOC cuts for SDP cones (sdp_soc) when not also disaggregating all nonpolyhedral cone cuts (all_disagg)\n") + end -# Create Pajarito nonlinear model: can solve with nonlinear algorithm only -function MathProgBase.NonlinearModel(s::PajaritoSolver) - if !applicable(MathProgBase.NonlinearModel, s.cont_solver) - error("Continuous solver (cont_solver) specified is not a NLP solver recognized by MathProgBase\n") + if (s.mip_subopt_count > 0) && (s.mip_subopt_solver == UnsetSolver()) + error("Using suboptimal solves (mip_subopt_count > 0), but no suboptimal MIP solver specified (set mip_subopt_solver)\n") end - # Translate options into old nonlinearmodel.jl fields - verbose = s.log_level - algorithm = (s.mip_solver_drives ? "BC" : "OA") - mip_solver = s.mip_solver - cont_solver = s.cont_solver - opt_tolerance = s.rel_gap - time_limit = s.timeout + if s.init_soc_one && !s.soc_disagg && !s.soc_abslift + error("Cannot use SOC initial L_1 cuts (init_soc_one) if both SOC disaggregation (soc_disagg) and SOC absvalue lifting (soc_abslift) are not used\n") + end - return PajaritoNonlinearModel(verbose, algorithm, mip_solver, cont_solver, opt_tolerance, time_limit) -end + if s.sdp_soc && s.mip_solver_drives + warn("In the MIP-solver-driven algorithm, SOC cuts for SDP cones (sdp_soc) cannot be added from subproblems or primal solutions, but they will be added from the conic relaxation\n") + end + if !s.solve_subp + s.prim_cuts_only = true + s.use_mip_starts = false + s.round_mip_sols = false + end + if s.prim_cuts_only + s.prim_cuts_always = true + end + if s.prim_cuts_always + s.prim_cuts_assist = true + end -# If input is in LinearQuadratic format, dispatch to conic or nonlinear algorithm depending on continuous solver -function MathProgBase.LinearQuadraticModel(s::PajaritoSolver) - if applicable(MathProgBase.ConicModel, s.cont_solver) || (s.cont_solver == UnsetSolver()) - MathProgBase.ConicToLPQPBridge(MathProgBase.ConicModel(s)) - else - MathProgBase.NonlinearToLPQPBridge(MathProgBase.NonlinearModel(s)) + if !s.all_disagg && s.prim_cuts_assist + error("Cannot use primal cuts when not disaggregating all nonpolyhedral cone cuts (all_disagg)\n") end + + return PajaritoConicModel(s.log_level, s.timeout, s.rel_gap, s.mip_solver_drives, s.mip_solver, s.mip_subopt_solver, s.mip_subopt_count, s.round_mip_sols, s.use_mip_starts, s.cont_solver, s.solve_relax, s.solve_subp, s.dualize_relax, s.dualize_subp, s.all_disagg, s.soc_disagg, s.soc_abslift, s.soc_in_mip, s.sdp_eig, s.sdp_soc, s.init_soc_one, s.init_soc_inf, s.init_exp, s.init_sdp_lin, s.init_sdp_soc, s.scale_subp_cuts, s.scale_subp_factor, s.scale_subp_up, s.viol_cuts_only, s.prim_cuts_only, s.prim_cuts_always, s.prim_cuts_assist, s.cut_zero_tol, s.prim_cut_feas_tol, s.dump_subproblems, s.dump_basename) end +# Create Pajaito LinearQuadratic model +MathProgBase.LinearQuadraticModel(s::PajaritoSolver) = MathProgBase.ConicToLPQPBridge(MathProgBase.ConicModel(s)) -# Return a vector of the supported cone types, for conic algorithm only +# Return a vector of the supported cone types function MathProgBase.supportedcones(s::PajaritoSolver) if s.cont_solver == UnsetSolver() # No conic solver, using primal cuts only, so support all Pajarito cones return [:Free, :Zero, :NonNeg, :NonPos, :SOC, :SOCRotated, :SDP, :ExpPrimal] - elseif applicable(MathProgBase.ConicModel, s.cont_solver) + else # Using conic solver, so supported cones are its cones (plus rotated SOC if SOC is supported) cones = MathProgBase.supportedcones(s.cont_solver) if :SOC in cones push!(cones, :SOCRotated) end return cones - else - # Solver must be NLP - error("Cannot get cones supported by continuous solver (cont_solver)\n") end end diff --git a/test/REQUIRE b/test/REQUIRE index 75b86f9..9ee360e 100644 --- a/test/REQUIRE +++ b/test/REQUIRE @@ -1,4 +1,3 @@ -GLPKMathProgInterface +GLPKMathProgInterface 0.4.0 0.4.1 ECOS SCS -Ipopt diff --git a/test/conictest.jl b/test/conictest.jl index 192896f..80050cf 100644 --- a/test/conictest.jl +++ b/test/conictest.jl @@ -4,7 +4,6 @@ # License, v. 2.0. If a copy of the MPL was not distributed with this # file, You can obtain one at http://mozilla.org/MPL/2.0/. - # Take a CBF file and conic solver and solve, redirecting output function solve_cbf(testname, probname, solver, redirect) flush(STDOUT) @@ -62,7 +61,6 @@ function solve_cbf(testname, probname, solver, redirect) return (status, time, objval, objbound, sol) end - # SOC problems for NLP and conic algorithms function run_soc(mip_solver_drives, mip_solver, cont_solver, log_level, redirect) testname = "SOC optimal" diff --git a/test/nlptest.jl b/test/nlptest.jl deleted file mode 100644 index d346821..0000000 --- a/test/nlptest.jl +++ /dev/null @@ -1,239 +0,0 @@ -# Copyright 2017, Chris Coey and Miles Lubin -# Copyright 2016, Los Alamos National Laboratory, LANS LLC. -# This Source Code Form is subject to the terms of the Mozilla Public -# License, v. 2.0. If a copy of the MPL was not distributed with this -# file, You can obtain one at http://mozilla.org/MPL/2.0/. - -# Take a JuMP model and solver and solve, redirecting output -function solve_jump(testname, m, redirect) - flush(STDOUT) - flush(STDERR) - @printf "%-30s... " testname - tic() - - if redirect - mktemp() do path,io - out = STDOUT - err = STDERR - redirect_stdout(io) - redirect_stderr(io) - - status = try - solve(m) - catch e - e - end - - flush(io) - redirect_stdout(out) - redirect_stderr(err) - end - else - status = try - solve(m) - catch e - e - end - end - flush(STDOUT) - flush(STDERR) - - if isa(status, ErrorException) - @printf ":%-16s %5.2f s\n" "ErrorException" toq() - else - @printf ":%-16s %5.2f s\n" status toq() - end - - flush(STDOUT) - flush(STDERR) - - return status -end - - -# Quadratically constrained problems compatible with MathProgBase ConicToLPQPBridge -function run_qp(mip_solver_drives, mip_solver, cont_solver, log_level, redirect) - solver=PajaritoSolver(timeout=120., mip_solver_drives=mip_solver_drives, mip_solver=mip_solver, cont_solver=cont_solver, log_level=(redirect ? 0 : 3)) - - testname = "QP optimal" - @testset "$testname" begin - m = Model(solver=solver) - - @variable(m, x >= 0, Int) - @variable(m, y >= 0) - @variable(m, 0 <= u <= 10, Int) - @variable(m, w == 1) - - @objective(m, Min, -3x - y) - - @constraint(m, 3x + 10 <= 20) - @constraint(m, y^2 <= u*w) - - status = solve_jump(testname, m, redirect) - - @test status == :Optimal - @test isapprox(getobjectivevalue(m), -12.162277, atol=TOL) - @test isapprox(getobjbound(m), -12.162277, atol=TOL) - @test isapprox(getvalue(x), 3, atol=TOL) - @test isapprox(getvalue(y), 3.162277, atol=TOL) - end - - testname = "QP maximize" - @testset "$testname" begin - m = Model(solver=solver) - - @variable(m, x >= 0, Int) - @variable(m, y >= 0) - @variable(m, 0 <= u <= 10, Int) - @variable(m, w == 1) - - @objective(m, Max, 3x + y) - - @constraint(m, 3x + 2y + 10 <= 20) - @constraint(m, x^2 <= u*w) - - status = solve_jump(testname, m, redirect) - - @test status == :Optimal - @test isapprox(getobjectivevalue(m), 9.5, atol=TOL) - @test isapprox(getobjbound(m), 9.5, atol=TOL) - @test isapprox(getvalue(x), 3, atol=TOL) - @test isapprox(getvalue(y), 0.5, atol=TOL) - end -end - -# NLP problems NOT compatible with MathProgBase ConicToLPQPBridge -function run_nlp(mip_solver_drives, mip_solver, cont_solver, log_level, redirect) - solver=PajaritoSolver(timeout=120., mip_solver_drives=mip_solver_drives, mip_solver=mip_solver, cont_solver=cont_solver, log_level=(redirect ? 0 : 3)) - - testname = "Nonconvex error" - @testset "$testname" begin - m = Model(solver=solver) - - @variable(m, x >= 0, start = 1, Int) - @variable(m, y >= 0, start = 1) - - @objective(m, Min, -3x - y) - - @constraint(m, 3x + 2y + 10 <= 20) - @NLconstraint(m, 8 <= x^2 <= 10) - - status = solve_jump(testname, m, redirect) - - @test isa(status, ErrorException) - end - - testname = "Optimal" - @testset "$testname" begin - m = Model(solver=solver) - - @variable(m, x >= 0, start = 1, Int) - @variable(m, y >= 0, start = 1) - - @objective(m, Min, -3x - y) - - @constraint(m, 3x + 2y + 10 <= 20) - @constraint(m, x >= 1) - @NLconstraint(m, x^2 <= 5) - @NLconstraint(m, exp(y) + x <= 7) - - status = solve_jump(testname, m, redirect) - - @test status == :Optimal - @test isapprox(getvalue(x), 2.0) - end - - testname = "Infeasible 1" - @testset "$testname" begin - m = Model(solver=solver) - - @variable(m, x >= 0, start = 1, Int) - @variable(m, y >= 0, start = 1) - - @objective(m, Min, -3x - y) - - @constraint(m, 3x + 2y + 10 <= 20) - @NLconstraint(m, x^2 >= 9) - @NLconstraint(m, exp(y) + x <= 2) - - status = solve_jump(testname, m, redirect) - - @test status == :Infeasible - end - - testname = "Infeasible 2" - @testset "$testname" begin - m = Model(solver=solver) - - @variable(m, x >= 0, start = 1, Int) - @variable(m, y >= 0, start = 1) - - @objective(m, Min, -3x - y) - - @constraint(m, 3x + 2y + 10 <= 20) - @constraint(m, 6x + 5y >= 30) - @NLconstraint(m, x^2 >= 8) - @NLconstraint(m, exp(y) + x <= 7) - - status = solve_jump(testname, m, redirect) - - @test status == :Infeasible - end - - testname = "Continuous error" - @testset "$testname" begin - m = Model(solver=solver) - - @variable(m, x >= 0, start = 1) - @variable(m, y >= 0, start = 1) - - @objective(m, Min, -3x - y) - - @constraint(m, 3x + 2y + 10 <= 20) - @constraint(m, x >= 1) - - @NLconstraint(m, x^2 <= 5) - @NLconstraint(m, exp(y) + x <= 7) - - status = solve_jump(testname, m, redirect) - - @test isa(status, ErrorException) - end - - testname = "Maximization" - @testset "$testname" begin - m = Model(solver=solver) - - @variable(m, x >= 0, start = 1, Int) - @variable(m, y >= 0, start = 1) - - @objective(m, Max, 3x + y) - - @constraint(m, 3x + 2y + 10 <= 20) - @NLconstraint(m, x^2 <= 9) - - status = solve_jump(testname, m, redirect) - - @test status == :Optimal - @test isapprox(getobjectivevalue(m), 9.5, atol=TOL) - end - - testname = "Nonlinear objective" - @testset "$testname" begin - m = Model(solver=solver) - - @variable(m, x >= 0, start = 1, Int) - @variable(m, y >= 0, start = 1) - - @objective(m, Max, -x^2 - y) - - @constraint(m, x + 2y >= 4) - @NLconstraint(m, x^2 <= 9) - - status = solve_jump(testname, m, redirect) - - @test status == :Optimal - @test isapprox(getobjectivevalue(m), -2.0, atol=TOL) - @test isapprox(getobjbound(m), -2.0, atol=TOL) - end -end diff --git a/test/qptest.jl b/test/qptest.jl new file mode 100644 index 0000000..0bf6f5e --- /dev/null +++ b/test/qptest.jl @@ -0,0 +1,102 @@ +# Copyright 2017, Chris Coey and Miles Lubin +# Copyright 2016, Los Alamos National Laboratory, LANS LLC. +# This Source Code Form is subject to the terms of the Mozilla Public +# License, v. 2.0. If a copy of the MPL was not distributed with this +# file, You can obtain one at http://mozilla.org/MPL/2.0/. + +# Take a JuMP model and solver and solve, redirecting output +function solve_jump(testname, m, redirect) + flush(STDOUT) + flush(STDERR) + @printf "%-30s... " testname + tic() + + if redirect + mktemp() do path,io + out = STDOUT + err = STDERR + redirect_stdout(io) + redirect_stderr(io) + + status = try + solve(m) + catch e + e + end + + flush(io) + redirect_stdout(out) + redirect_stderr(err) + end + else + status = try + solve(m) + catch e + e + end + end + flush(STDOUT) + flush(STDERR) + + if isa(status, ErrorException) + @printf ":%-16s %5.2f s\n" "ErrorException" toq() + else + @printf ":%-16s %5.2f s\n" status toq() + end + + flush(STDOUT) + flush(STDERR) + + return status +end + +# Quadratically constrained problems compatible with MathProgBase ConicToLPQPBridge +function run_qp(mip_solver_drives, mip_solver, cont_solver, log_level, redirect) + solver=PajaritoSolver(timeout=120., mip_solver_drives=mip_solver_drives, mip_solver=mip_solver, cont_solver=cont_solver, log_level=(redirect ? 0 : 3)) + + testname = "QP optimal" + @testset "$testname" begin + m = Model(solver=solver) + + @variable(m, x >= 0, Int) + @variable(m, y >= 0) + @variable(m, 0 <= u <= 10, Int) + @variable(m, w == 1) + + @objective(m, Min, -3x - y) + + @constraint(m, 3x + 10 <= 20) + @constraint(m, y^2 <= u*w) + + status = solve_jump(testname, m, redirect) + + @test status == :Optimal + @test isapprox(getobjectivevalue(m), -12.162277, atol=TOL) + @test isapprox(getobjbound(m), -12.162277, atol=TOL) + @test isapprox(getvalue(x), 3, atol=TOL) + @test isapprox(getvalue(y), 3.162277, atol=TOL) + end + + testname = "QP maximize" + @testset "$testname" begin + m = Model(solver=solver) + + @variable(m, x >= 0, Int) + @variable(m, y >= 0) + @variable(m, 0 <= u <= 10, Int) + @variable(m, w == 1) + + @objective(m, Max, 3x + y) + + @constraint(m, 3x + 2y + 10 <= 20) + @constraint(m, x^2 <= u*w) + + status = solve_jump(testname, m, redirect) + + @test status == :Optimal + @test isapprox(getobjectivevalue(m), 9.5, atol=TOL) + @test isapprox(getobjbound(m), 9.5, atol=TOL) + @test isapprox(getvalue(x), 3, atol=TOL) + @test isapprox(getvalue(y), 0.5, atol=TOL) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 8d75f82..148d1a9 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,7 +9,6 @@ import ConicBenchmarkUtilities using Pajarito using Base.Test - # Tests absolute tolerance and Pajarito printing level TOL = 1e-3 ll = 3 @@ -17,7 +16,7 @@ redirect = true # Define dictionary of solvers, using JuMP list of available solvers include(Pkg.dir("JuMP", "test", "solvers.jl")) -include("nlptest.jl") +include("qptest.jl") include("conictest.jl") solvers = Dict{String,Dict{String,MathProgBase.AbstractMathProgSolver}}() @@ -58,15 +57,6 @@ end # solvers["MILP"]["SCIP"] = solvers["MISOCP"]["SCIP"] = SCIP.SCIPSolver("display/verblevel", 0, "limits/gap", tol_gap, "numerics/feastol", tol_feas) # end -# NLP solvers -solvers["NLP"] = Dict{String,MathProgBase.AbstractMathProgSolver}() -if ipt - solvers["NLP"]["Ipopt"] = Ipopt.IpoptSolver(print_level=0) -end -# if kni -# solvers["NLP"]["Knitro"] = KNITRO.KnitroSolver(objrange=1e16, outlev=0, maxit=100000) -# end - # Conic solvers solvers["SOC"] = Dict{String,MathProgBase.AbstractMathProgSolver}() solvers["Exp+SOC"] = Dict{String,MathProgBase.AbstractMathProgSolver}() @@ -86,7 +76,6 @@ if mos solvers["Exp+SOC"]["Mosek"] = solvers["PSD+Exp"]["Mosek"] = Mosek.MosekSolver(LOG=0, MSK_DPAR_INTPNT_CO_TOL_REL_GAP=1e-9, MSK_DPAR_INTPNT_CO_TOL_PFEAS=1e-10, MSK_DPAR_INTPNT_CO_TOL_DFEAS=1e-10, MSK_DPAR_INTPNT_CO_TOL_NEAR_REL=1e3) end - println("\nSolvers:") for (stype, snames) in solvers println("\n$stype") @@ -94,8 +83,6 @@ for (stype, snames) in solvers @printf "%2d %s\n" i sname end end -println() - @testset "Algorithm - $(msd ? "MSD" : "Iter")" for msd in [false, true] alg = (msd ? "MSD" : "Iter") @@ -106,23 +93,11 @@ println() continue end - @testset "NLP models, NLP solver - $conname" for (conname, con) in solvers["NLP"] - println("\nNLP models, NLP solver: $alg, $mipname, $conname") - run_qp(msd, mip, con, ll, redirect) - run_nlp(msd, mip, con, ll, redirect) - end - @testset "LPQP models, SOC solver - $conname" for (conname, con) in solvers["SOC"] println("\nLPQP models, SOC solver: $alg, $mipname, $conname") run_qp(msd, mip, con, ll, redirect) end - @testset "Exp+SOC models, NLP solver - $conname" for (conname, con) in solvers["NLP"] - println("\nExp+SOC models, NLP solver: $alg, $mipname, $conname") - run_soc(msd, mip, con, ll, redirect) - run_expsoc(msd, mip, con, ll, redirect) - end - @testset "SOC models/solver - $conname" for (conname, con) in solvers["SOC"] println("\nSOC models/solver: $alg, $mipname, $conname") run_soc(msd, mip, con, ll, redirect)