Consider variable bounds clean

[LukasBarner]

Sorry for the mess on #189, I'm really not experienced in this at all :D

This is a clean version including all old tests now for both consider_constrained_variables true and false.
Still some failing tests for "Bilevel Conic JuMP NLP" and "Sum aggregation_group".
Since some tests in these testsets failed also on master (for me), "Bilevel Conic JuMP NLP" does not seem to be a problem. However I'm not sure about "Sum aggregation_group"...

I guess it could make sense to limit aggregation_group when working with consider_constrained_variables ?

Please let me know if there is anything in terms of style etc. so I can improve.

[joaquimg]

Good start!

• We need this to work with vector sets.
• Aggregation groups should just work (we can leave it for a little later)
• I don't think we should test everything twice in the final version (we can leave it there while we finish the code.
• Bilevel Conic requires vector sets

[joaquimg]

these two functions look the same

[LukasBarner]

Yes, this was an old relict as I wasn't sure about constants etc in the beginning. Now, one is deleted and the other is defined for union of lesser/greater. Equalities still throw an error.

[joaquimg]

I dont think any of this is needed

[LukasBarner]

true, this threw occasional errors earlier. Found out that this was related to missing dual variables for which start values could not be passed.
This is now resolved, since hints etc. are no longer passed if the dual to bounds does not exist in the first place...

[joaquimg]

this should be a field in config.
Also, I am not sure this should be run in all tests.

[LukasBarner]

Left this for now, as I'm not sure which tests to select best etc.
Personally, I don't think that running tests twice adds a huge amount of time, as almost everything should be compiled already and the problems are fairly small.
On my machine, running tests like this adds about 6 seconds (less than 7%)

[LukasBarner]

I could also do the subset, but I'm not sure which ones to select...
Do you have a preference?

[joaquimg]

if the overhead is only 7% for running twice, this would indeed be fine

[joaquimg]

I dont think this should be run in all tests that would lead to run everything twice.
I think we should select a subset.

[LukasBarner]

See above

[joaquimg]

any reason for not supporting other sets like vector sets?
you can follow the code in get_canonical_complement, as you have been doing

[LukasBarner]

left this out because I did not really feel certain.
Now added, can you check if it is correct?
Also not sure about the todo (see code for the new function)

[joaquimg]

Note that CI will only run if you fix conflicts with master.

[LukasBarner]

can you tell me if I resolved the conflicts correctly?

[joaquimg]

same function twice?

[LukasBarner]

sorry, thought I already replied to these when I made the commits ...
now a conditional

[joaquimg]

why a function if a simple conditional would suffice?

[LukasBarner]

done, see above

[joaquimg]

SCALAR_SETS ?

[LukasBarner]

Would prefer this to throw an error in case of equalities (i.e. if something unexpected happens).
Could also go with SCALAR_SETS for code readability in case you prefer and do @Assert inside the function?

[joaquimg]

you need MOI.copy just like in get_canonical_complement

[LukasBarner]

Now done, but only for the dual_func. This is the only necessary one, right?
I'm asking because the copy is done for sets etc. in other functions (but that should be handled by set_with_zero, right?

PS: This was a big one, so I'm a bit surprised it didn't affect testing...

[joaquimg]

do what?

[LukasBarner]

There was the following todo in the code of get_canonical_complement (note that I commented one additional line out):

    # dim = MOI.dimension(set)
    # vector sets have no constant
    # for i in 1:dim
    #     func.constant[i] = Dualization.set_dot(i, set, T) *
    #         Dualization.get_scalar_term(primal_model, ci, i)
    # end
    # todo - set dot on function

[LukasBarner]

Quick heads up: I just noticed that there are still some problems with setting correct start values. However, considering #185, this PR is still an overall improvement.
In comparison to before, ProductMode works for constrained variables now (probably also a few others).
I can still confirm issues with at least SOS1Mode and FortunyAmatMcCarlMode. There are a few things that need additional adjustment for the variable complements here.
However, there are some more things that do not work as expected even when not using the newly added variable complements...
I have a larger sample problem for which I can compute a feasible starting point (by finding a solution to the lower level with fixed upper level variables). When passing it in ProductMode, I get a feasible starting point to the single level model. However, when I use FortunyAmatMcCarlMode, starting points are not accepted (not an issue of different solver tolerances etc., violations are orders of magnitude larger).
We could either take care of investigating the start value thing (for both newly added and existing functionality) in this PR, or move either both or one to a new issue/PR. I think I would prefer a new issue/PR for both, what do you think @joaquimg ?

[joaquimg]

We have to figure out what the issues are with SOS1Mode and FortunyAmatMcCarlMode.
They may be more fundamental things that are working in ProductMode just by chance.

Start values should be passed, this should not be hard.

Start values are frequently ignored by MIP solvers (used in FortunyAmatMcCarlMode). How do you know they are not accepted? MIP solvers usually require starting points to absolutely all variables involved.

If things don't work for existing functionality, we should have issues to discuss them separately.

This PR should not be blocked by things that were not previously working, but it should be blocked by things that previously worked.

[LukasBarner]

Quick disclaimer: As I said, I'm pretty new to this and not 100% sure about everything. In case it is nonsense, please excuse the effort of reading through all of it...

First a few replies to the text above:

Start values should be passed, this should not be hard.

100% agree with the first part (especially since it will be a heavily used functionality), but not so sure that this is a trivial issue (more in the second part...)

Start values are frequently ignored by MIP solvers (used in FortunyAmatMcCarlMode). How do you know they are not accepted?

For smaller examples, start values are sometimes accepted. Using Gurobi, I get notifications when start values violate constraints (i.e. constraint X is violated by 60.08). This is usually the case for larger models.

MIP solvers usually require starting points to absolutely all variables involved.

My code in general does set all start values correctly, this can be verified when switching to ProductMode with Ipopt, where initial point primal infeasibility of the single level model constraints is about e-7. Note that this only holds if consider_constrained_variables = true, otherwise the duals to variable bounds are not being set correctly (more on this also below). Also, I do not believe that ProductMode is working "by chance", things there should be largely correct...

Now a few remarks to the existing functionality and why I'm not sure starting values are a trivial issue.

Currently, dual start values to lower level variable bounds are not set. Some exceptions exist, such as FortunyAmatMcCarlMode, where duals are sometimes set to zero (which tbh I'm not sure is actually correct: If the constraint was tight, why do we set the dual variable start to 0? ):

BilevelJuMP.jl/src/moi.jl

Lines 1181 to 1185 in 565c0ef

	bin = MOI.add_variable(m)
	if pass_start && has_start && is_tight
	MOI.set(m, MOI.VariablePrimalStart(), bin, 1.0)
	MOI.set(m, MOI.VariablePrimalStart(), dual, 0.0)
	else

As mentioned in #185, model.ctr_lower does not contain variable bounds currently (see the issue for an MWE demonstrating this). I believe, that in the following lines, this leads to their dual starts not being passed at all (exception above):

BilevelJuMP.jl/src/jump.jl

Lines 515 to 524 in 2d25399

	for (idx, info) in model.ctr_info
	if haskey(model.ctr_lower, idx)
	ctr = model.ctr_lower[idx]
	# this fails for vector-constrained variables due dualization 0.3.5
	# because of constrained variables that change the dual
	pre_duals = lower_primal_dual_map.primal_con_dual_var[JuMP.index(ctr)] # vector
	duals = map(x->lower_dual_to_sblm[x], pre_duals)
	pass_dual_info(single_blm, duals, info)
	end
	end

Furhter, I believe another reason why I have failing tests on larger models with start values in FortunyAmatMcCarlMode could be:

BilevelJuMP.jl/src/moi.jl

Lines 1158 to 1167 in 565c0ef

	if pass_start
	val = MOIU.eval_variables(
	x -> nothing_to_nan(MOI.get(m, MOI.VariablePrimalStart(), x)),
	f_dest,
	)
	if !isnan(val)
	is_tight = abs(val) < 1e-8
	has_start = true
	end
	end

Potentially, numerical imprecision and solver tolerances may lead to deviations larger than 1e-8, falsely classifying a constraint as not tight. I know this to be true in my case, but I'm not sure its the only issue.
Another way to set the start values of binaries could be to take a look at the dual start values and check which deviation from zero of either the function or the dual variable is bigger.
A problem though is that the dual starts are not yet available in this part of the code...

These lines are actually also a problem for the newly introduced variable complements:

BilevelJuMP.jl/src/moi.jl

Lines 1158 to 1167 in 565c0ef

	if pass_start
	val = MOIU.eval_variables(
	x -> nothing_to_nan(MOI.get(m, MOI.VariablePrimalStart(), x)),
	f_dest,
	)
	if !isnan(val)
	is_tight = abs(val) < 1e-8
	has_start = true
	end
	end

Since duals do not yet have start values set, val may evaluate to NaN for the new variable complements despite all starts being set correctly...

For these reasons, I think it might make sense to already pass start values of duals in the 'build_bilevel' function before starting to set up complementary slackness (and to make sure that all starts are actually passed of course :D).

As you can see, there are a few issues that can be discussed, but I'm afraid it would blow up this PR...
This is why I would prefer to make start values a new issue once this PR is ready.
What do you think @joaquimg ?

PS: I have not taken a deeper dive into SOS1, but I figure similar issues are going on there...

[LukasBarner]

@joaquimg, did you have a chance to take a look the comment above yet?
I would like to take a deeper look as I need the start values, but I'm not sure which way to go (new PR with a more structured list of todos etc., or integrating into this one).

[LukasBarner]

We have to figure out what the issues are with SOS1Mode and FortunyAmatMcCarlMode.
They may be more fundamental things that are working in ProductMode just by chance.

You were right here! There is something more fundamental...
I already caught one problem when handling scalar set constants, let me check if there is another.
Will correct soon...

[LukasBarner]

Okay, think I have figured this one out...
Also have everything that is written below on a new branch (which I could either make a new start value PR, or include to this one, whichever you prefer).
The codes mentioned below are also tested on a larger model and work as expected there.
I have split things up into stuff related to this PR and things that are also valid for the current master.
Sry for the long text...

Things that are fundamentally wrong (but easy to fix)

Obvious mistake in this PR:

First the general mistake, these lines I have written are obvious nonsense:
https://github.com/LukasBarner/BilevelJuMP.jl/blob/422b4060f28006abdfb69d3921475670b28d8743/src/moi.jl#L317-L321

if MOI.constant(dual_set) > 0
    dual_func.constant = dual_func.constant - MOI.constant(dual_set)  
elseif MOI.constant(dual_set) < 0
    dual_func.constant = dual_func.constant + MOI.constant(dual_set)
end

They should instead be:

dual_func.constant = dual_func.constant - MOI.constant(dual_set)

Passing of dual starts:

Currently, dual starts are passed after the bilevel problem is built:

BilevelJuMP.jl/src/jump.jl

Lines 571 to 598 in 565c0ef

	single_blm,
	upper_to_sblm,
	lower_to_sblm,
	lower_primal_dual_map,
	lower_dual_to_sblm = build_bilevel(
	upper,
	lower,
	moi_link,
	moi_upper,
	mode,
	moi_link2;
	copy_names = model.copy_names,
	pass_start = model.pass_start,
	)
	# pass additional info (hints - not actual problem data)
	# for lower level dual variables (start, upper hint, lower hint)
	for (idx, info) in model.ctr_info
	if haskey(model.ctr_lower, idx)
	ctr = model.ctr_lower[idx]
	# this fails for vector-constrained variables due dualization 0.3.5
	# because of constrained variables that change the dual
	pre_duals =
	lower_primal_dual_map.primal_con_dual_var[JuMP.index(ctr)] # vector
	duals = map(x -> lower_dual_to_sblm[x], pre_duals)
	pass_dual_info(single_blm, duals, info)
	end
	end

This is fine if we do not have the new "variable complements". However, the new "variable complements" also have dual variables inside the ConstraintFunction, and some start value related functionalities (i.e. setting starting values of binary variables) rely upon evaluating starts:

BilevelJuMP.jl/src/moi.jl

Lines 1158 to 1167 in 565c0ef

	if pass_start
	val = MOIU.eval_variables(
	x -> nothing_to_nan(MOI.get(m, MOI.VariablePrimalStart(), x)),
	f_dest,
	)
	if !isnan(val)
	is_tight = abs(val) < 1e-8
	has_start = true
	end
	end

Since the dual start values are not yet passed for the variable complements, this tanks the whole warmstart procedure for MIP solvers that require fully feasibly start values (IMO usually the case).
To make warmstarts work here, we have to move passing dual starts into the build_bilevel(...) function. Luckily, this is just copy-pasting and works fine.
We need to pass the bilevel model as an argument to build_bilevel and then copy this piece

BilevelJuMP.jl/src/jump.jl

Lines 586 to 609 in 565c0ef

	# pass additional info (hints - not actual problem data)
	# for lower level dual variables (start, upper hint, lower hint)
	for (idx, info) in model.ctr_info
	if haskey(model.ctr_lower, idx)
	ctr = model.ctr_lower[idx]
	# this fails for vector-constrained variables due dualization 0.3.5
	# because of constrained variables that change the dual
	pre_duals =
	lower_primal_dual_map.primal_con_dual_var[JuMP.index(ctr)] # vector
	duals = map(x -> lower_dual_to_sblm[x], pre_duals)
	pass_dual_info(single_blm, duals, info)
	end
	end
	# pass lower & upper level primal variables info (upper, lower)
	for (idx, info) in model.var_info
	if haskey(model.var_lower, idx)
	var = lower_to_sblm[JuMP.index(model.var_lower[idx])]
	elseif haskey(model.var_upper, idx)
	var = upper_to_sblm[JuMP.index(model.var_upper[idx])]
	else
	continue
	end
	pass_primal_info(single_blm, var, info)
	end

and paste it here:

BilevelJuMP.jl/src/moi.jl

Lines 460 to 464 in 565c0ef

	end
	#=
	Additional Optimiality conditions (to complete the KKT)
	=#

prior to the complement part (Additional Optimiality conditions).

Start values for SOS1 and FortunyAmatMcCarl

Now to the issue with start values for SOS1 and FortunyAmatMcCarl:

SOS1 Starts:

Currently, no start values are set for the newly added slack variables here:

BilevelJuMP.jl/src/moi.jl

Lines 765 to 803 in 565c0ef

	function add_complement(
	mode::SOS1Mode{T},
	m,
	comp::Complement,
	idxmap_primal,
	idxmap_dual,
	copy_names::Bool,
	pass_start::Bool,
	) where {T}
	f = comp.func_w_cte
	s = comp.set_w_zero
	v = comp.variable
	if comp.is_vec
	error("Vector constraint is not supported by SOS1 mode")
	end
	slack, slack_in_set = MOI.add_constrained_variable(m, s)
	f_dest = MOIU.map_indices.(Ref(idxmap_primal), f)
	new_f = MOIU.operate(-, T, f_dest, slack)
	equality = MOIU.normalize_and_add_constraint(m, new_f, MOI.EqualTo(zero(T)))
	dual = idxmap_dual[v]
	c1 = MOI.add_constraint(
	m,
	MOI.VectorOfVariables([slack, dual]),
	MOI.SOS1([1.0, 2.0]),
	)
	if copy_names
	nm = MOI.get(m, MOI.VariableName(), dual)
	MOI.set(m, MOI.VariableName(), slack, "slk_($(nm))")
	# MOI.set(m, MOI.ConstraintName(), slack_in_set, "bound_slk_($(nm))")
	MOI.set(m, MOI.ConstraintName(), equality, "eq_slk_($(nm))")
	MOI.set(m, MOI.ConstraintName(), c1, "compl_sos1_($(nm))")
	end
	return slack, slack_in_set, equality, c1
	end

This also holds if we do not consider constrained variables during dualization. However, it can easily be fixed by something like this inside the add_complement(...) function:

slack_start = MOIU.eval_variables(
    x -> nothing_to_nan(MOI.get(m, MOI.VariablePrimalStart(), x)),
    f_dest,
)
MOI.set(m, MOI.VariablePrimalStart(), slack, slack_start)

Now, SOS1 start values work as expected!

FortunyAmatMcCarl Starts:

Here, binary starts (and duals to zero) were mixed up:

BilevelJuMP.jl/src/moi.jl

Lines 1182 to 1187 in 565c0ef

	if pass_start && has_start && is_tight
	MOI.set(m, MOI.VariablePrimalStart(), bin, 1.0)
	MOI.set(m, MOI.VariablePrimalStart(), dual, 0.0)
	else
	MOI.set(m, MOI.VariablePrimalStart(), bin, 0.0)
	end

Which I believe should instead be:

if pass_start && has_start && is_tight
    MOI.set(m, MOI.VariablePrimalStart(), bin, 0.0)
else
    MOI.set(m, MOI.VariablePrimalStart(), bin, 1.0)
    MOI.set(m, MOI.VariablePrimalStart(), dual, 0.0)
end

Since MIP solvers usually require a complete set of starts, I believe that user provided duals should not be over-written, so

if pass_start && has_start && is_tight
    MOI.set(m, MOI.VariablePrimalStart(), bin, 0.0)
else
    MOI.set(m, MOI.VariablePrimalStart(), bin, 1.0)
end

might even be a better solution...

Finally, this line can make problems when providing warmstarts:

BilevelJuMP.jl/src/moi.jl

Line 1164 in 565c0ef

is_tight = abs(val) < 1e-8

a) For larger models, it is likely that some numerical imprecision occurs.
b) If warmstarts are obtained by solving the lower level for fixed upper level variables, it is not guaranteed that solver tolerances are tighter than 1e-8 (Gurobi for example uses 1-6 by default).

One workaround would be to make the tolerance an argument that can be specified by the user.
But since MIP solvers usually require a complete set of warmstarts, there is an easier solution. We can just compare val to the dual variable (this should be more numerically stable and eliminates user defined tolerances).
Something like this does work (only if dual starts are already passed in the build_bilevel(...) function though):

is_tight = abs(val) < abs(nothing_to_nan(MOI.get(m, MOI.VariablePrimalStart(), dual)))

[joaquimg]

Your comments seem reasonable indeed.
You can either update this PR or propose new ones, I am very keen to see them.
Remember the key point:
Small PR's are much easier to review and, consequently, quicker to integrate. We don't want to add bugs by accident.
Large PR's are much harder and slower to review since I have to go through each line and think about possible issues.

[LukasBarner]

Ok, decided to keep them in this PR. This should now also fix #185.
I understand that this makes reviews more difficult, but there are synergies between the consider_constrained_variables case and the start values.
I think warm starts (especially for the integer methods) are extremely useful and will be frequently used, so fixing them should have a high priority. In case this review will take significantly longer to get ready (or fail due to other problems), I could also do an "interim PR" that ignores some of the consider_constrained_variables functionality to make start value fixes available quicker. What do you think?