Bug with 2 warm-started optimizations

[franckgaga]

Hi,

I'm writing a ModelPredictiveControl (MPC) package based on JuMP.jl. The default solver is currently OSQP.jl. It is also possible to use my MPC package with DAQP.jl, thanks to your work and the JuMP.jl interface!

It seems that there is an issue with resolving a warm-started optimization problem. If I do:

using JuMP, DAQP
model = Model(DAQP.Optimizer)
@variable(model, -1<= x1 <=1)
@variable(model, -2<= x2 <=2)
@objective(model, Min, 0.5*(x1^2+x2^2)+x1+x2)
@constraint(model, c1, -3 <= x1 + x2 <= 3)
@constraint(model, c2, -4 <= x1 - x2 <= 4)
# one warm-started optimization : works well
set_start_value(x1, -0.5)
set_start_value(x2, -0.5)
optimize!(model); value.((x1, x2))

It returns as expected :

(-1.0, -1.0)

But if I do two optimizations:

using JuMP, DAQP
model = Model(DAQP.Optimizer)
@variable(model, -1<= x1 <=1)
@variable(model, -2<= x2 <=2)
@objective(model, Min, 0.5*(x1^2+x2^2)+x1+x2)
@constraint(model, c1, -3 <= x1 + x2 <= 3)
@constraint(model, c2, -4 <= x1 - x2 <= 4)
# one warm-started optimization : works well
set_start_value(x1, -0.5)
set_start_value(x2, -0.5)
optimize!(model)
# second warm-started optimization : crash with error
set_start_value(x1, -0.5)
set_start_value(x2, -0.5)
optimize!(model); value.((x1, x2))

it crashes with this error :

ERROR: MathOptInterface.UnsupportedAttribute{MathOptInterface.VariablePrimalStart}: Attribute MathOptInterface.VariablePrimalStart() is not supported by the model.
Stacktrace:
  [1] throw_set_error_fallback(model::Optimizer, attr::MathOptInterface.VariablePrimalStart, index::MathOptInterface.VariableIndex, value::Float64; error_if_supported::MathOptInterface.SetAttributeNotAllowed{MathOptInterface.VariablePrimalStart})
    @ MathOptInterface ~/.julia/packages/MathOptInterface/864xP/src/attributes.jl:600
  [2] throw_set_error_fallback(model::Optimizer, attr::MathOptInterface.VariablePrimalStart, index::MathOptInterface.VariableIndex, value::Float64)
    @ MathOptInterface ~/.julia/packages/MathOptInterface/864xP/src/attributes.jl:590
  [3] set(::Optimizer, ::MathOptInterface.VariablePrimalStart, ::MathOptInterface.VariableIndex, ::Float64)
    @ MathOptInterface ~/.julia/packages/MathOptInterface/864xP/src/attributes.jl:550
  [4] set(m::MathOptInterface.Utilities.CachingOptimizer{Optimizer, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, attr::MathOptInterface.VariablePrimalStart, index::MathOptInterface.VariableIndex, value::Float64)
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/864xP/src/Utilities/cachingoptimizer.jl:791
  [5] set
    @ ~/.julia/packages/MathOptInterface/864xP/src/Bridges/bridge_optimizer.jl:1235 [inlined]
  [6] set(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}, attr::MathOptInterface.VariablePrimalStart, index::MathOptInterface.VariableIndex, value::Float64)
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/864xP/src/Utilities/cachingoptimizer.jl:791
  [7] set(model::Model, attr::MathOptInterface.VariablePrimalStart, v::VariableRef, value::Float64)
    @ JuMP ~/.julia/packages/JuMP/jZvaU/src/optimizer_interface.jl:721
  [8] set_start_value
    @ ~/.julia/packages/JuMP/jZvaU/src/variables.jl:1136 [inlined]
  [9] _broadcast_getindex_evalf
    @ ./broadcast.jl:683 [inlined]
 [10] _broadcast_getindex
    @ ./broadcast.jl:656 [inlined]
 [11] getindex
    @ ./broadcast.jl:610 [inlined]
 [12] copy
    @ ./broadcast.jl:888 [inlined]
 [13] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{0}, Nothing, typeof(set_start_value), Tuple{Base.RefValue{VariableRef}, Float64}})
    @ Base.Broadcast ./broadcast.jl:873
 [14] top-level scope
    @ ~/Documents/Personnel/ModelPredictiveControl.jl/src/yo.jl:13

Thanks !

[darnstrom]

Hi Francis, First of, nice work with the MPC package! :)

Thank you for reporting this error. I have added support for VariablePrimalStart in #13. It will be merged into main when version 0.5.0 of DAQP_jll is available. Although, since DAQP is a dual solver, a primal warm start cannot directly be used (it could indirectly be used by using it to get an estimation of the active set, which I might look into in the future). But at least there should be no error now!

[franckgaga]

Thanks for the info!

My knowledge on solvers and active sets is quite limited. Out of curiosity, why the first solve accept set_start_value calls and not the second one ? I'm not sure I understrand the difference between "starting value" and "primal warm start". I know that primals are related to the decision variables and duals, to the constraints, but that's about it.

[darnstrom]

@franckgaga there is nothing in the solver itself that make a difference of "starting value" and "primal warm start". The original error that you found is more something related to the internals of JuMP: JuMP only seems to call "VariablePrimalStart" when you set the primal value to an optimizer that has been called once before (your original error was due to VariablePrimalStart not being defined in DAQP.jl) . The reason why VariablePrimalStart is only called in the second call is more a question for the JuMP developers; my guess, however, is that there are situation when one wants to do something special when a warm start is given to an already optimized model (e.g., how internal factorizations in the solver should be modified to be reused etc).

[darnstrom]

Fixed in v0.5.0 (#13)

[franckgaga]

Strangely, the update seems to resolve the error with my MWE above, but not with my MPC package and I really don't understand why. If I use DAQPv0.5.0 with my package:

using ModelPredictiveControl, ControlSystemsBase, JuMP, DAQP
mpc = LinMPC(LinModel(tf(2,[10, 1]), 0.1), optim=Model(DAQP.Optimizer))

and comment out the exception handling related to warm start in predictive_control.jl, that is :

set_start_value.(ΔŨvar, ΔŨ0)
set_objective!(mpc, ΔŨvar)
#try
optimize!(optim)
#catch err
#    if isa(err, MOI.UnsupportedAttribute{MOI.VariablePrimalStart})
#        # reset_optimizer to unset warm-start, set_start_value.(nothing) seems buggy
#        MOIU.reset_optimizer(optim) 
#        optimize!(optim)
#    else
#        rethrow(err)
#    end
#end

and run this:

mpc = LinMPC(LinModel(tf(2,[10, 1]), 0.1), optim=Model(DAQP.Optimizer))
mpc([1]); mpc([1])

The error I get in the is a bit different :

ERROR: MathOptInterface.UnsupportedAttribute{MathOptInterface.VariablePrimalStart}: Attribute MathOptInterface.VariablePrimalStart() is not supported by the model.
Stacktrace:
  [1] copy_to_check_attributes(dest::Optimizer, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}})
    @ DAQP ~/.julia/packages/DAQP/8XHFD/src/MOI_wrapper/MOI_wrapper.jl:335
  [2] copy_to(dest::Optimizer, src::MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}})
    @ DAQP ~/.julia/packages/DAQP/8XHFD/src/MOI_wrapper/MOI_wrapper.jl:289
  [3] optimize!
    @ ~/.julia/packages/MathOptInterface/864xP/src/MathOptInterface.jl:84 [inlined]
  [4] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{Optimizer, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/864xP/src/Utilities/cachingoptimizer.jl:316
  [5] optimize!
    @ ~/.julia/packages/MathOptInterface/864xP/src/Bridges/bridge_optimizer.jl:376 [inlined]
  [6] optimize!(m::MathOptInterface.Utilities.CachingOptimizer{MathOptInterface.Bridges.LazyBridgeOptimizer{MathOptInterface.Utilities.CachingOptimizer{Optimizer, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}}}, MathOptInterface.Utilities.UniversalFallback{MathOptInterface.Utilities.Model{Float64}}})
    @ MathOptInterface.Utilities ~/.julia/packages/MathOptInterface/864xP/src/Utilities/cachingoptimizer.jl:325
  [7] optimize!(model::Model; ignore_optimize_hook::Bool, _differentiation_backend::MathOptInterface.Nonlinear.SparseReverseMode, kwargs::Base.Pairs{Symbol, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
    @ JuMP ~/.julia/packages/JuMP/jZvaU/src/optimizer_interface.jl:435
  [8] optimize!
    @ ~/.julia/packages/JuMP/jZvaU/src/optimizer_interface.jl:405 [inlined]
  [9] optim_objective!(mpc::LinMPC{SteadyKalmanFilter})
    @ ModelPredictiveControl ~/Documents/Personnel/ModelPredictiveControl.jl/src/predictive_control.jl:530
 [10] moveinput!(mpc::LinMPC{SteadyKalmanFilter}, ry::Vector{Int64}, d::Vector{Float64}; R̂y::Vector{Int64}, D̂::Vector{Float64}, ym::Nothing)
    @ ModelPredictiveControl ~/Documents/Personnel/ModelPredictiveControl.jl/src/predictive_control.jl:266
 [11] moveinput!(mpc::LinMPC{SteadyKalmanFilter}, ry::Vector{Int64}, d::Vector{Float64})
    @ ModelPredictiveControl ~/Documents/Personnel/ModelPredictiveControl.jl/src/predictive_control.jl:253
 [12] #_#19
    @ ~/Documents/Personnel/ModelPredictiveControl.jl/src/predictive_control.jl:1038 [inlined]
 [13] PredictiveController
    @ ~/Documents/Personnel/ModelPredictiveControl.jl/src/predictive_control.jl:1033 [inlined]
 [14] (::LinMPC{SteadyKalmanFilter})(ry::Vector{Int64})
    @ ModelPredictiveControl ~/Documents/Personnel/ModelPredictiveControl.jl/src/predictive_control.jl:1033
 [15] top-level scope
    @ REPL[17]:1

@darnstrom, do you have any idea why ?

Additionally, do you have any idea why set_start_value.(ΔŨvar, nothing) does not unset the starting values for DAQP (I still get a Attribute MathOptInterface.VariablePrimalStart() is not supported by the model error after) ? This is the official way to unset starting values in JuMP.jl. I think it would be cleaner to use this instead of MOIU.reset_optimizer(optim) in my exception handling (since we probably loose the factorization with reset_optimizer).

Thanks for the help!

[darnstrom]

I think it is because I forgot to add a MOI.support for VariablePrimalStart (still strange that it worked for MWE). @franckgaga, could you check if the code on the branch v051 resolve this error on your end?

[franckgaga]

Sadly no, it does not seems to resolve the error.

I also created the branch v062warm that exclude the exception handling related to warm start, if you want to reproduce on your side :

(@v1.9) pkg> add ModelPredictiveControl#v062warm

after installed, you can run:

using ModelPredictiveControl, ControlSystemsBase, JuMP, DAQP
mpc = LinMPC(LinModel(tf(2,[10, 1]), 0.1), optim=Model(DAQP.Optimizer))
mpc([1]); mpc([1])

and I get the same error as above.

Thanks !

[darnstrom]

@franckgaga The code v051 seems to be working on my end now with v062warm. If you can confirm that it works as expected on your end I will merge it into main and release v0.5.1 that includes it.

[franckgaga]

Yes it is also working on my side. Thanks for your time !