User Tools

Site Tools


jump

Differences

This shows you the differences between two versions of the page.

Link to this comparison view

jump [2018/12/27 13:27] (current)
Line 1: Line 1:
 +
 +~~CLOSETOC~~
 +
 +~~TOC 1-3 wide~~
 +
 +---
 +
 +^  [[roots|Roots and Solving]] ​ ^  [[Optim|Compact Optimization]] ​ ^  Linear and Squared Optimization ​ ^
 +
 +---
 +
 +
 +```juliarepl
 +julia> pkgchk.( [ "​julia"​ => v"​1.0.3"​ , "​JuMP"​ => v"​0.18.0",​ "​Clp"​ => v"​0.3.1",​ "​DataStructures"​ => v"​0.7.4",​ "​Ipopt"​ => v"​0.3.0"​ ]; )
 +```
 +
 +
 +# JuMP --- Linear and Quadratic Constrained Problems
 +
 +JuMP 0.18 under Julia 0.6.2 is best suited to solving linear and quadratic problems, subjects to linear constraints.
 +
 +In principle, it can also do general non-linear problems, too.  However, JuMP (1) it seems to be too much tool for the job; (2) suffers from changing interfaces; and (3) just seems unreliable and not yet ready. ​ We recommend [IntervalOptimisation.jl](https://​github.com/​JuliaIntervals/​IntervalOptimisation.jl),​ [Optim.jl](https://​github.com/​JuliaNLSolvers/​Optim.jl),​ or [NLopt.jl](https://​github.com/​JuliaOpt/​NLopt.jl) for unconstrained non-linear optimization,​ instead.
 +
 +
 +## A Linear Optimization Program
 +
 +This is the introductory problem in the documentation.
 +
 +```juliarepl
 +julia> using JuMP                           ## the frontend which interprets user's commands
 +
 +julia> using Clp                            ## the backend solver
 +
 +julia> mymodel = Model(solver = ClpSolver())
 +Feasibility problem with:
 + * 0 linear constraints
 + * 0 variables
 +Solver is ClpMathProg
 +
 +julia> @variable(mymodel,​ 0 <= x <= 2 )      ## first variable, with box constraint
 +x
 +
 +julia> @variable(mymodel,​ 0 <= y <= 30 )     ## second variable, with box constraint
 +y
 +
 +julia> @constraint(mymodel,​ 1x + 5y <= 3.0 ) ## a non-box constraint
 +x + 5 y ≤ 3
 +
 +julia> @objective(mymodel,​ Max, 5*x + 3*y )  ## the objective function
 +5 x + 3 y
 +
 +julia> print(mymodel) ​                       ## nice summary
 +Max 5 x + 3 y
 +Subject to
 + x + 5 y ≤ 3
 + 0 ≤ x ≤ 2
 + 0 ≤ y ≤ 30
 +
 +julia> status = solve(mymodel) ​              ## now find the solution
 +:Optimal
 +
 +julia> ( getobjectivevalue(mymodel),​ getvalue(x),​ getvalue(y) )
 +(10.6, 2.0, 0.2)
 +
 +```
 +
 +* Array optimizing variables are more flexible than suggested here, e.g., `@variable(m,​ x[1:M,1:N] >= 0 )`, `@variable(m,​ x[i=1:​10,​j=1:​10;​ isodd(i+j)] >= 0)`, or `@variable(m,​ x[i=1:10], start=(i/​2))`.
 +
 +
 +
 +
 +### Interrogating the Model
 +
 +
 +```julia
 +using JuMP, Clp, DataStructures
 +
 +mymodel = Model(solver = ClpSolver());​
 +
 +@variable(mymodel,​ 0 <= x <= 2 );
 +
 +@variable(mymodel,​ 0 <= y <= 30 );
 +
 +@constraint(mymodel,​ 1x + 5y <= 3.0 );   ## cone constraints:​ @constraint(m,​ norm(A*x) <= 2w - 1)
 +
 +@objective(mymodel,​ Max, 5*x + 3*y );
 +
 +status = solve(mymodel); ​    ## solution could be optimal, unbounded, infeasible, userlimit, error, or notsolved
 +
 +function describeModel( mymodel )
 + println(mymodel)
 + jumpfuns= OrderedDict{Function,​ String}(
 +     MathProgBase.numvar ​ => " number of variables associated with the Model m ",
 +     MathProgBase.numvar => "​number of variables ",
 +     MathProgBase.numlinconstr => "​number of linear constraints ",
 +     MathProgBase.numquadconstr => "​number of quadratic constraints ",
 +     MathProgBase.numconstr => "total number of constraints ",
 +     JuMP.numsocconstr => "​number of second order cone constraints ",
 +     JuMP.numsosconstr => "​number of sos constraints ",
 +     JuMP.numsdconstr => "​number of semi-definite constraints ",
 +     JuMP.numnlconstr => "​number of nonlinear constraints ",
 +     getsolvetime => "solve time reported by the solver if it is implemented.",​
 +     getnodecount => "​number of explored branch-and-bound nodes, if it is implemented.",​
 +     getobjbound => "best known bound on the optimal objective value --- e.g., when branch-and-bound method stops early",​
 +     getobjectivevalue => "​objective function value",​
 +     getobjectivesense => "min or max",
 +     getobjectivebound => "best known bound on the objective",​
 +     getobjgap => "final relative optimality gap as optimization terminated. That is, it returns |b−f||f|, where b is the best bound and f is the best feasible objective value.",​
 +     getsimplexiter => "​cumulative number of simplex iterations during the optimization process. In particular, for a MIP ittotal simplex iterations for all nodes.",​
 +     getbarrieriter => "​cumulative number of barrier iterations during the optimization process."​
 + );
 + for k in keys(jumpfuns)
 +     funresult= try k(mymodel) catch "not implemented"​ end
 +     println( "| ", k, " ​ |  ", funresult, " ​ | ", jumpfuns[k],​ " |" )
 + end#for
 +end##​function##​
 +
 +describeModel( mymodel )
 +```
 +
 +provides
 +
 +```text
 +Max 5 x + 3 y
 +Subject to
 + x + 5 y ≤ 3
 + 0 ≤ x ≤ 2
 + 0 ≤ y ≤ 30
 +```
 +
 +^ **Variable** ​ ^  **Value** ^ **Explanation** ​ ^
 +| MathProgBase.SolverInterface.numvar ​ |  2  | number of variables |
 +| MathProgBase.SolverInterface.numlinconstr ​ |  1  | number of linear constraints |
 +| MathProgBase.SolverInterface.numquadconstr ​ |  0  | number of quadratic constraints |
 +| MathProgBase.SolverInterface.numconstr ​ |  1  | total number of constraints |
 +| JuMP.numsocconstr ​ |  0  | number of second order cone constraints |
 +| JuMP.numsosconstr ​ |  0  | number of sos constraints |
 +| JuMP.numsdconstr ​ |  0  | number of semi-definite constraints |
 +| JuMP.numnlconstr ​ |  0  | number of nonlinear constraints |
 +| MathProgBase.SolverInterface.getsolvetime ​ |  not implemented ​ | solve time reported by the solver if it is implemented. |
 +| MathProgBase.SolverInterface.getnodecount ​ |  not implemented ​ | number of explored branch-and-bound nodes, if it is implemented. |
 +| MathProgBase.SolverInterface.getobjbound ​ |  not implemented ​ | best known bound on the optimal objective value --- e.g., when branch-and-bound method stops early |
 +| JuMP.getobjectivevalue ​ |  10.6  | objective function value |
 +| JuMP.getobjectivesense ​ |  Max  | min or max |
 +| JuMP.getobjectivebound ​ |  NaN  | best known bound on the objective |
 +| MathProgBase.SolverInterface.getobjgap ​ |  not implemented ​ | final relative optimality gap as optimization terminated. |
 +| MathProgBase.SolverInterface.getsimplexiter ​ |  not implemented ​ | cumulative number of simplex iterations during the optimization process. |
 +| MathProgBase.SolverInterface.getbarrieriter ​ |  not implemented ​ | cumulative number of barrier iterations during the optimization process. |
 +
 +
 +There are other model variables, like getrawsolver (returns an object that may be used to access a solver-specific API"), but they are less useful. ​ At the moment, the `getvalue()` functions strangely do not require the model input.
 +
 +
 +
 +## Caveat: The Chosen Solver Must Support The Optimization Problem
 +
 +JuMP is not capable of choosing the best problem for the solver. ​ This task is left to the user:
 +
 +```juliarepl
 +julia> using JuMP, Clp
 +
 +julia> mymodel = Model(solver = ClpSolver());​
 +
 +julia> @variable(mymodel,​ 0 <= x <= 2 );
 +
 +julia> @variable(mymodel,​ 0 <= y <= 30 );
 +
 +julia> @objective(mymodel,​ Max, 5*x^2 + 3*y^2+x*y );
 +
 +julia> status = solve(mymodel) ​              ## now find the solution
 +ERROR: Solver does not support quadratic objectives
 +Stacktrace:
 + [1] setquadobjterms!(::​Clp.ClpMathProgSolverInterface.ClpMathProgModel,​ ::​Array{Int32,​1},​ ::​Array{Int32,​1},​ ::​Array{Float64,​1}) at /​Users/​ivo/​.julia/​v0.6/​MathProgBase/​src/​SolverInterface/​LinearQuadratic.jl:​71
 + [2] addQuadratics(::​JuMP.Model) at /​Users/​ivo/​.julia/​v0.6/​JuMP/​src/​solvers.jl:​435
 + [3] #​build#​119(::​Bool,​ ::Bool, ::​JuMP.ProblemTraits,​ ::Function, ::​JuMP.Model) at /​Users/​ivo/​.julia/​v0.6/​JuMP/​src/​solvers.jl:​373
 + [4] (::​JuMP.#​kw##​build)(::​Array{Any,​1},​ ::​JuMP.#​build,​ ::​JuMP.Model) at ./<​missing>:​0
 + [5] #​solve#​116(::​Bool,​ ::Bool, ::Bool, ::​Array{Any,​1},​ ::Function, ::​JuMP.Model) at /​Users/​ivo/​.julia/​v0.6/​JuMP/​src/​solvers.jl:​168
 + [6] solve(::​JuMP.Model) at /​Users/​ivo/​.julia/​v0.6/​JuMP/​src/​solvers.jl:​150
 +```
 +
 +
 +Sometimes, for a non-linear `NLobjective`,​ the error message can be more descriptive
 +
 +```julianoeval
 +julia> @NLobjective(m,​ Min, (1-x)^2 + 100(y-x^2)^2)
 +
 +ERROR: No solver was provided. JuMP has classified this model as NLP. Julia
 +packages which provide solvers for this class of problems include Ipopt,
 +KNITRO, and Mosek. The solver must be specified by using either the "​solver="​
 +keyword argument to "​Model()"​ or the "​setsolver()"​ method.
 +
 +```
 +
 +
 +
 +
 +## Non-Linear Optimization
 +
 +Again, we are not (yet) recommending JuMP for unconstrained non-linear optimization. ​ However, it can work.
 +
 +
 +### A Working Example
 +
 +Here is a working example:
 +
 +```juliarepl
 +julia> using JuMP, Ipopt
 +
 +julia> m = Model( solver = IpoptSolver());​
 +
 +julia> @variable(m,​ 0 <= x[1:2] <= 4);
 +
 +julia> @constraint(m,​ sum(x) == 1.5);
 +
 +julia> @constraint(m,​ x[1] <= x[2]);
 +
 +julia> @NLconstraint(m,​ x[1]^2+x[2]^2 == 2);
 +
 +julia> @objective(m,​ Min, x[1]*x[2]);
 +
 +julia> solve(m);
 +
 +******************************************************************************
 +This program contains Ipopt, a library for large-scale nonlinear optimization.
 + Ipopt is released as open source code under the Eclipse Public License (EPL).
 +         For more information visit http://​projects.coin-or.org/​Ipopt
 +******************************************************************************
 +
 +This is Ipopt version 3.12.8, running with linear solver mumps.
 +NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
 +
 +Number of nonzeros in equality constraint Jacobian...: ​       4
 +Number of nonzeros in inequality constraint Jacobian.: ​       2
 +Number of nonzeros in Lagrangian Hessian.............: ​       3
 +
 +Total number of variables............................: ​       2
 +                     ​variables with only lower bounds: ​       0
 +                variables with lower and upper bounds: ​       2
 +                     ​variables with only upper bounds: ​       0
 +Total number of equality constraints.................: ​       2
 +Total number of inequality constraints...............: ​       1
 +        inequality constraints with only lower bounds: ​       0
 +   ​inequality constraints with lower and upper bounds: ​       0
 +        inequality constraints with only upper bounds: ​       1
 +
 +iter    objective ​   inf_pr ​  ​inf_du lg(mu) ​ ||d||  lg(rg) alpha_du alpha_pr ​ ls
 +   ​0 ​ 9.9999800e-05 2.00e+00 1.00e+00 ​ -1.0 0.00e+00 ​   -  0.00e+00 0.00e+00 ​  0
 +   ​1 ​ 5.9107754e-01 8.12e-01 5.58e+08 ​ -1.0 7.96e-01 ​   -  1.40e-02 1.00e+00f ​ 1
 +   ​2 ​ 1.1193824e-02 3.33e-01 1.05e+10 ​ -1.0 1.19e+01 ​   -  6.14e-03 1.22e-01f ​ 1
 +   ​3 ​ 1.1624711e-01 1.75e-02 1.26e+09 ​ -1.0 1.84e-01 ​   -  9.87e-01 1.00e+00f ​ 1
 +   ​4 ​ 1.2495708e-01 8.58e-05 1.03e+07 ​ -1.0 1.31e-02 ​   -  1.00e+00 1.00e+00h ​ 1
 +   ​5 ​ 1.2500000e-01 2.11e-09 5.07e+02 ​ -1.0 6.49e-05 ​   -  1.00e+00 1.00e+00h ​ 1
 +
 +Number of Iterations....:​ 5
 +
 +                                   ​(scaled) ​                ​(unscaled)
 +Objective...............: ​  ​1.2499999894745066e-01 ​   1.2499999894745066e-01
 +Dual infeasibility......: ​  ​6.6613381477509392e-16 ​   6.6613381477509392e-16
 +Constraint violation....: ​  ​2.1050983178838578e-09 ​   2.1050983178838578e-09
 +Complementarity.........: ​  ​0.0000000000000000e+00 ​   0.0000000000000000e+00
 +Overall NLP error.......: ​  ​2.1050983178838578e-09 ​   2.1050983178838578e-09
 +
 +
 +Number of objective function evaluations ​            = 6
 +Number of objective gradient evaluations ​            = 6
 +Number of equality constraint evaluations ​           = 6
 +Number of inequality constraint evaluations ​         = 6
 +Number of equality constraint Jacobian evaluations ​  = 6
 +Number of inequality constraint Jacobian evaluations = 6
 +Number of Lagrangian Hessian evaluations ​            = 5
 +Total CPU secs in IPOPT (w/o function evaluations) ​  ​= ​     0.116
 +Total CPU secs in NLP function evaluations ​          ​= ​     0.024
 +
 +EXIT: Optimal Solution Found.
 +
 +julia> getvalue(x),​ getobjectivevalue(m)
 +([0.0885622,​ 1.41144], 0.12499999894745066)
 +```
 +
 +
 +
 +### Non-Working Examples
 +
 +#### Bad Syntax Now
 +
 +This example used to work, but no longer does.  It complains about the model syntax.
 +
 +```julianoeval
 +julia> using JuMP, Ipopt
 +
 +julia> function testme(n)
 +           ​myf(a...) = sum(collect(a).^2) ​   ## the simplest quadratic function is to be optimized
 +           m= JuMP.Model(solver=Ipopt.IpoptSolver())
 +           ​JuMP.register(m,​ :myf, n, myf, autodiff=true) ​  ## important: myf needs to be registered
 +           ​@JuMP.variable(m,​ x[1:(n-1)] >= 0.5)
 +           ​JuMP.setNLobjective(m,​ :Min, Expr(:call, :myf, [x[i] for i=1:n]...))
 +           ​JuMP.solve(m)
 +           ​return [getvalue(x[i]) for i=1:n]
 +       end##
 +testme (generic function with 1 method)
 +
 +julia> testme(3)
 +
 +******************************************************************************
 +This program contains Ipopt, a library for large-scale nonlinear optimization.
 + Ipopt is released as open source code under the Eclipse Public License (EPL).
 +         For more information visit http://​projects.coin-or.org/​Ipopt
 +******************************************************************************
 +
 +This is Ipopt version 3.12.8, running with linear solver mumps.
 +NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
 +
 +Number of nonzeros in equality constraint Jacobian...: ​       0
 +Number of nonzeros in inequality constraint Jacobian.: ​       0
 +Number of nonzeros in Lagrangian Hessian.............: ​       0
 +
 +Total number of variables............................: ​       3
 +                     ​variables with only lower bounds: ​       3
 +                variables with lower and upper bounds: ​       0
 +                     ​variables with only upper bounds: ​       0
 +Total number of equality constraints.................: ​       0
 +Total number of inequality constraints...............: ​       0
 +        inequality constraints with only lower bounds: ​       0
 +   ​inequality constraints with lower and upper bounds: ​       0
 +        inequality constraints with only upper bounds: ​       0
 +
 +iter    objective ​   inf_pr ​  ​inf_du lg(mu) ​ ||d||  lg(rg) alpha_du alpha_pr ​ ls
 +   ​0 ​ 7.8029997e-01 0.00e+00 2.00e-02 ​  0.0 0.00e+00 ​   -  0.00e+00 0.00e+00 ​  0
 +   ​1 ​ 7.5030000e-01 0.00e+00 9.70e-03 ​ -8.0 1.01e-02 ​   -  1.00e+00 9.80e-01f ​ 1
 +   ​2 ​ 7.5000291e-01 0.00e+00 1.11e-15 -10.0 9.90e-05 ​   -  1.00e+00 1.00e+00f ​ 1
 +   ​3 ​ 7.4999997e-01 0.00e+00 1.11e-16 -11.0 9.80e-07 ​   -  1.00e+00 1.00e+00f ​ 1
 +
 +Number of Iterations....:​ 3
 +
 +                                   ​(scaled) ​                ​(unscaled)
 +Objective...............: ​  ​7.4999997003576602e-01 ​   7.4999997003576602e-01
 +Dual infeasibility......: ​  ​1.1102230246251565e-16 ​   1.1102230246251565e-16
 +Constraint violation....: ​  ​0.0000000000000000e+00 ​   0.0000000000000000e+00
 +Complementarity.........: ​  ​1.1921907667178425e-11 ​   1.1921907667178425e-11
 +Overall NLP error.......: ​  ​1.1921907667178425e-11 ​   1.1921907667178425e-11
 +
 +
 +Number of objective function evaluations ​            = 4
 +Number of objective gradient evaluations ​            = 4
 +Number of equality constraint evaluations ​           = 0
 +Number of inequality constraint evaluations ​         = 0
 +Number of equality constraint Jacobian evaluations ​  = 0
 +Number of inequality constraint Jacobian evaluations = 0
 +Number of Lagrangian Hessian evaluations ​            = 0
 +Total CPU secs in IPOPT (w/o function evaluations) ​  ​= ​     0.101
 +Total CPU secs in NLP function evaluations ​          ​= ​     0.023
 +
 +EXIT: Optimal Solution Found.
 +3-element Array{Float64,​1}:​
 + 0.5
 + 0.5
 + 0.5
 +```
 +
 +
 +#### Incorrect Answer
 +
 +The following problem feeds a kinky unconstrained function to JuMP.  JuMP returns an incorrect answer.
 +
 +```julianoeval
 +using JuMP, Ipopt
 +
 +nearneghyperbola( a... ) = prod( -1.0./​(abs.( collect(a) .+ 3) + 1e-2) )
 +
 +function learn_nl_jump(userfun,​ n)
 +    mymdl= JuMP.Model( solver=Ipopt.IpoptSolver() )
 +    JuMP.register( mymdl, :userfun, n, userfun, autodiff=true )
 +    @JuMP.variable( mymdl, -100 <= x[1:n] <= 100 )
 +    JuMP.setNLobjective( mymdl, :Min, Expr(:call, :userfun, [x[i] for i=1:n]...) )
 +    JuMP.solve( mymdl )
 +    return [ getvalue(x[i]) for i=1:n ]        ## why does getvalue not refer to mymdl??
 +end
 +
 +xopt= learn_nl_jump(nearneghyperbola,​ 3)
 +
 +println( "x^* =: ", xopt, "​\n"​)
 +println( "​minimum f(x^*)= ", nearneghyperbola(xopt...) )
 +println( "f(neg 3s)= ", nearneghyperbola(-3,​ -3, -3) )
 +
 +```
 +
 +The obvious solution is `( -3, -3, -3 )` for n=3.  However, while it would have been ok for JuMP to report failure, it is more worrisome that JuMP failed to recognize that it had gone off the deep end:
 +
 +```text
 +EXIT: Solved To Acceptable Level.
 +x^* =: [-99.9897, -99.9897, 94.0]
 +
 +minimum f(x^*)= -1.0955764466765225e-6
 +f(neg 3s)= -1.0e6
 +```
 +
 +
 +# Backmatter
 +
 +## Useful Packages on Julia Repository
 +
 +## Notes
 +
 +### Non-Linear Optimization
 +
 +* See also [IntervalOptimisation.jl](https://​github.com/​JuliaIntervals/​IntervalOptimisation.jl).
 +
 +* See also [Optim.jl](https://​github.com/​JuliaNLSolvers/​Optim.jl)
 +
 +* See also [NLopt.jl](https://​github.com/​JuliaOpt/​NLopt.jl)
 +
 +## References
  
jump.txt · Last modified: 2018/12/27 13:27 (external edit)