# Costates and Hamiltonian in JuMP and MathProgbase

**URL:** <https://discourse.julialang.org/t/costates-and-hamiltonian-in-jump-and-mathprogbase/6160>\
**Category:** Optimization (Mathematical)\
**Created:** [September 29, 2017, 11:09pm UTC](https://discourse.julialang.org/t/costates-and-hamiltonian-in-jump-and-mathprogbase/6160 "2017-09-29T23:09:44Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![huckl3b3rry87](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/huckl3b3rry87/32/3163_2.png) [@huckl3b3rry87](https://discourse.julialang.org/u/huckl3b3rry87)\
**Post date:** [September 29, 2017, 11:09pm UTC](https://discourse.julialang.org/t/costates-and-hamiltonian-in-jump-and-mathprogbase/6160/1 "2017-09-29T23:09:44Z")

</div>

**Part # 1**  
I was interested in getting the costates. The costates are also known as the lagrange multipliers or dual variables. I looked through the MathProgBase docs and found [this function](http://mathprogbasejl.readthedocs.io/en/latest/nlp.html#getconstrduals) called `getconstduals` does that give me my lagrange multipliers?

It does not seems like it. If I am solving a the Bryson Denham problem, I can match the analytical solution for the controls and states

 ![main](https://global.discourse-cdn.com/julialang/original/3X/5/7/575465987e8adcc7f95a324759e8ae54e27b3fd3.png)

But, if I use try to pull out the costates from my model

```julia
m=JuMP.internalmodel(n.mdl);
c=MathProgBase.getconstrduals(m)

```

and then find the indecies that match my constraints (for the first state x1), and do some scaling, I get a picture that looks like this  
 ![Figure_1-1](https://global.discourse-cdn.com/julialang/original/3X/9/f/9feb4945793465d6334f7a2daa0557c7f195d0a8.png)

When it is supposed to look like  
 ![darby](https://global.discourse-cdn.com/julialang/original/3X/4/7/4718272b1263ca30476e706d611f52d546008bf9.png)

So, they are flipped. Similarly, for the other costate it is the same trend, but flipped.

Does anyone know what is going on here?

**Part # 2**  
Also, I was interested in getting the Hamiltonian, but I did not have much luck on this one. Any thoughts/suggestions?

Thanks!

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [September 30, 2017, 5:36am UTC](https://discourse.julialang.org/t/costates-and-hamiltonian-in-jump-and-mathprogbase/6160/2 "2017-09-30T05:36:06Z")

</div>

Well, there’s this comment in [a post about upcoming changes in JuMP](https://discourse.julialang.org/t/mathoptinterface-and-upcoming-breaking-changes-in-jump-0-19/4874):

> - Duality. MOI aims to define consistent conventions for duals across problem classes. In some cases, the signs may differ from the current MPB conventions.

Another possibility (just guessing), but maybe the sign change could be due to an internal rewrite of the constraints in a canonical form. If you have, say, `@constraint(m, x == 1)` in your code, it could be that while setting up the internal model, the constraint is rewritten as `1 - x == 0`, with the Lagrange multiplier (dual variable) acting on `1 - x`, while you were maybe expecting the Lagrange multiplier to act on `x - 1`.

JuMP defines `getdual` methods for `ConstraintRef`s as well by the way; I don’t think you should have to access the internal model:

```julia
julia> methods(getdual)
# 7 methods for generic function "getdual":
getdual(c::JuMP.ConstraintRef{JuMP.Model,JuMP.GenericRangeConstraint{JuMP.NonlinearExprData}}) in JuMP at /Users/twan/code/julia/RigidBodyDynamics/v0.6/JuMP/src/nlp.jl:56
getdual(c::JuMP.ConstraintRef{JuMP.Model,JuMP.SDConstraint}) in JuMP at /Users/twan/code/julia/RigidBodyDynamics/v0.6/JuMP/src/JuMP.jl:715
getdual(c::JuMP.ConstraintRef{JuMP.Model,JuMP.GenericSOCConstraint{JuMP.GenericNormExpr{2,Float64,JuMP.Variable}}}) in JuMP at /Users/twan/code/julia/RigidBodyDynamics/v0.6/JuMP/src/JuMP.jl:676
getdual(c::JuMP.ConstraintRef{JuMP.Model,JuMP.GenericRangeConstraint{JuMP.GenericAffExpr{Float64,JuMP.Variable}}}) in JuMP at /Users/twan/code/julia/RigidBodyDynamics/v0.6/JuMP/src/JuMP.jl:586
getdual(v::JuMP.Variable) in JuMP at /Users/twan/code/julia/RigidBodyDynamics/v0.6/JuMP/src/JuMP.jl:466
getdual(x::AbstractArray) in JuMP at /Users/twan/code/julia/RigidBodyDynamics/v0.6/JuMP/src/JuMPContainer.jl:109
getdual(x::JuMP.JuMPContainer) in JuMP at /Users/twan/code/julia/RigidBodyDynamics/v0.6/JuMP/src/JuMPContainer.jl:114

```

So you can do e.g.:

```julia
julia> m = Model(solver = GurobiSolver(OutputFlag = 0));

julia> @variable(m, x);

julia> @constraint(m, c, x == 1)
x = 1

julia> @objective(m, :Min, x)
x

julia> solve(m)
:Optimal

julia> getdual(c)
1.0

```

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [September 30, 2017, 2:40pm UTC](https://discourse.julialang.org/t/costates-and-hamiltonian-in-jump-and-mathprogbase/6160/3 "2017-09-30T14:40:14Z")

</div>

On a related note, I find it confusing that for

```julia
using JuMP, Gurobi
m = Model(solver = GurobiSolver(OutputFlag = false))
@variable(m, x)
@objective(m, :Min, x)
@constraint(m, c, x >= 1)
solve(m)
getdual(c)

```

we get 1.0, but for the same code with the `@objective` line replaced by `@objective(m, :Max, -x)` (really, the same problem), we get -1.0.

The Duality chapter in [Boyd and Vandenberghe](https://web.stanford.edu/~boyd/cvxbook/), starts with the canonical form (5.1), i.e. a minimization. I would expect that for a maximization problem, JuMP would automatically turn maximizing `-x` into minimizing `x` to fit (5.1), and then return the duals for that problem (so that the dual would be the same in both cases above). Maybe this is part of the ‘consistent conventions for duals’ mentioned in the post.

---

<div class="post-metadata">

**Author:** ![tkoolen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tkoolen/32/1603_2.png) [@tkoolen](https://discourse.julialang.org/u/tkoolen)\
**Post date:** [September 30, 2017, 2:51pm UTC](https://discourse.julialang.org/t/costates-and-hamiltonian-in-jump-and-mathprogbase/6160/4 "2017-09-30T14:51:10Z")

</div>

As for the Hamiltonian thing, that’s a concept from optimal control, not really from generic optimization, so you’ll have to code that up yourself.

---

<div class="post-metadata">

**Author:** ![huckl3b3rry87](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/huckl3b3rry87/32/3163_2.png) [@huckl3b3rry87](https://discourse.julialang.org/u/huckl3b3rry87)\
**Post date:** [October 1, 2017, 2:53am UTC](https://discourse.julialang.org/t/costates-and-hamiltonian-in-jump-and-mathprogbase/6160/5 "2017-10-01T02:53:42Z")

</div>

@tkoolen thank you very much!