diff --git a/src/nlp.jl b/src/nlp.jl index 5bc50040517..5263b566533 100644 --- a/src/nlp.jl +++ b/src/nlp.jl @@ -9,10 +9,11 @@ type NLPData nlobj nlconstr::Vector{NonlinearConstraint} nlconstrlist::ReverseDiffSparse.ExprList + nlconstrDuals::Vector{Float64} evaluator end -NLPData() = NLPData(nothing, NonlinearConstraint[], ExprList(), nothing) +NLPData() = NLPData(nothing, NonlinearConstraint[], ExprList(), Float64[], nothing) Base.copy(::NLPData) = error("Copying nonlinear problems not yet implemented") @@ -22,6 +23,15 @@ function initNLP(m::Model) end end +function getDual(c::ConstraintRef{NonlinearConstraint}) + initNLP(c.m) + nldata::NLPData = c.m.nlpdata + if length(nldata.nlconstrDuals) != length(nldata.nlconstr) + error("Dual solution not available. Check that the model was properly solved.") + end + return nldata.nlconstrDuals[c.idx] +end + type JuMPNLPEvaluator <: MathProgBase.AbstractNLPEvaluator m::Model A::SparseMatrixCSC{Float64,Int} # linear constraint matrix @@ -501,8 +511,18 @@ function solvenlp(m::Model; suppress_warnings=false) m.objVal = MathProgBase.getobjval(m.internalModel) m.colVal = MathProgBase.getsolution(m.internalModel) - if stat != :Optimal && !suppress_warnings - warn("Not solved to optimality, status: $stat") + if stat != :Optimal + suppress_warnings || warn("Not solved to optimality, status: $stat") + else # stat == :Optimal + if applicable(MathProgBase.getconstrduals, m.internalModel) && applicable(MathProgBase.getreducedcosts, m.internalModel) + nlduals = MathProgBase.getconstrduals(m.internalModel) + m.linconstrDuals = nlduals[1:length(m.linconstr)] + # quadratic duals currently not available, formulate as nonlinear constraint if needed + nldata.nlconstrDuals = nlduals[length(m.linconstr)+length(m.quadconstr)+1:end] + m.redCosts = MathProgBase.getreducedcosts(m.internalModel) + else + suppress_warnings || Base.warn_once("Nonlinear solver does not provide dual solutions") + end end m.internalModelLoaded = true diff --git a/test/nonlinear.jl b/test/nonlinear.jl index ae94130e397..bba8ea0d871 100644 --- a/test/nonlinear.jl +++ b/test/nonlinear.jl @@ -313,6 +313,86 @@ context("With solver $(typeof(nlp_solver))") do @fact getValue(zexpr[1]) => roughly(-(1/4)*log(1/4), 1e-4) end; end; end +facts("[nonlinear] Test nonlinear duals") do +for nlp_solver in nlp_solvers +applicable(MathProgBase.getconstrduals, MathProgBase.model(nlp_solver)) || continue +context("With solver $(typeof(nlp_solver))") do + modA = Model(solver=nlp_solver) + @defVar(modA, x >= 0) + @defVar(modA, y <= 5) + @defVar(modA, 2 <= z <= 4) + @defVar(modA, 0 <= r[i=3:6] <= i) + @setNLObjective(modA, Min, -((x + y)/2.0 + 3.0)/3.0 - z - r[3]) + @addConstraint(modA, cons1, x+y >= 2) + @addConstraint(modA, cons2, sum{r[i],i=3:5} <= (2 - x)/2.0) + @addNLConstraint(modA, cons3, 7.0*y <= z + r[6]/1.9) + + # Solution + @fact solve(modA) => :Optimal + @fact getObjectiveValue(modA) => roughly(-5.8446115, 1e-6) + @fact getValue(x) => roughly(0.9774436, 1e-6) + @fact getValue(y) => roughly(1.0225563, 1e-6) + @fact getValue(z) => roughly(4.0, 1e-6) + @fact getValue(r)[3] => roughly(0.5112781, 1e-6) + @fact getValue(r)[4] => roughly(0.0, 1e-6) + @fact getValue(r)[5] => roughly(0.0, 1e-6) + @fact getValue(r)[6] => roughly(6.0, 1e-6) + + # Reduced costs + @fact getDual(x) => roughly( 0.0, 1e-6) + @fact getDual(y) => roughly( 0.0, 1e-6) + @fact getDual(z) => roughly(-1.0714286, 1e-6) + @fact getDual(r)[3] => roughly( 0.0, 1e-6) + @fact getDual(r)[4] => roughly(1.0, 1e-6) + @fact getDual(r)[5] => roughly(1.0, 1e-6) + @fact getDual(r)[6] => roughly(-0.03759398, 1e-6) + + # Row duals + @fact getDual(cons1) => roughly( 0.333333, 1e-6) + @fact getDual(cons2) => roughly(-1.0, 1e-6) + @fact getDual(cons3) => roughly(-0.0714286, 1e-6) +end; end; end + +facts("[nonlinear] Test nonlinear duals (Max)") do +for nlp_solver in nlp_solvers +applicable(MathProgBase.getconstrduals, MathProgBase.model(nlp_solver)) || continue +context("With solver $(typeof(nlp_solver))") do + modA = Model(solver=nlp_solver) + @defVar(modA, x >= 0) + @defVar(modA, y <= 5) + @defVar(modA, 2 <= z <= 4) + @defVar(modA, 0 <= r[i=3:6] <= i) + @setNLObjective(modA, Max, ((x + y)/2.0 + 3.0)/3.0 + z + r[3]) + @addConstraint(modA, cons1, x+y >= 2) + @addConstraint(modA, cons2, sum{r[i],i=3:5} <= (2 - x)/2.0) + @addNLConstraint(modA, cons3, 7.0*y <= z + r[6]/1.9) + + # Solution + @fact solve(modA) => :Optimal + @fact getObjectiveValue(modA) => roughly(5.8446115, 1e-6) + @fact getValue(x) => roughly(0.9774436, 1e-6) + @fact getValue(y) => roughly(1.0225563, 1e-6) + @fact getValue(z) => roughly(4.0, 1e-6) + @fact getValue(r)[3] => roughly(0.5112781, 1e-6) + @fact getValue(r)[4] => roughly(0.0, 1e-6) + @fact getValue(r)[5] => roughly(0.0, 1e-6) + @fact getValue(r)[6] => roughly(6.0, 1e-6) + + # Reduced costs + @fact getDual(x) => roughly( 0.0, 1e-6) + @fact getDual(y) => roughly( 0.0, 1e-6) + @fact getDual(z) => roughly(1.0714286, 1e-6) + @fact getDual(r)[3] => roughly( 0.0, 1e-6) + @fact getDual(r)[4] => roughly(-1.0, 1e-6) + @fact getDual(r)[5] => roughly(-1.0, 1e-6) + @fact getDual(r)[6] => roughly(0.03759398, 1e-6) + + # Row duals + @fact getDual(cons1) => roughly(-0.333333, 1e-6) + @fact getDual(cons2) => roughly(1.0, 1e-6) + @fact getDual(cons3) => roughly(0.0714286, 1e-6) +end; end; end + ############################################################################# # Test that output is produced in correct MPB form