support for duals in nonlinear problems, ref #456

jump-dev · Jul 8, 2015 · 88cf9ea · 88cf9ea
1 parent 3debb48
commit 88cf9ea
Show file tree

Hide file tree

Showing 2 changed files with 103 additions and 3 deletions.
diff --git 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
@@ -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