Costates and Hamiltonian in JuMP and MathProgbase

Specific Domains Optimization (Mathematical)
1

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 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
main700×700 33.1 KB

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

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

When it is supposed to look like
darby

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!

2

Well, there’s this comment in a post about upcoming changes in JuMP:

  • 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 ConstraintRefs as well by the way; I don’t think you should have to access the internal model:

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> 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
3

On a related note, I find it confusing that for

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, 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.

4

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.

5

@tkoolen thank you very much!

2 Likes