Unexpected results using solver pao.pyomo.FA

[etoenges]

I applied the pao.pyomo.FA solver to the following problem:

from pao.pyomo import *

M = pe.ConcreteModel()
M.x = pe.Var(bounds=(0, None), domain=pe.Reals)

M.L = SubModel(fixed=M.x)
M.L.y = pe.Var(bounds=(0, 10), domain=pe.Reals) #This is the relevant code line

M.L.ymax = 3
M.xmin = 2

def ul_obj_rule(M):
    """upper-level objective"""
    mo = M.model()
    return(mo.x + mo.L.y)

def ll_obj_rule(M):
    mo = M.model()
    return(mo.x + mo.L.y)

def ul_constraint_rule(M):
    mod = M.model()
    return(M.x >= mod.xmin)

def ll_constraint_rule(M):
    mod = M.model()
    return(mod.L.y <= mod.x + M.ymax)

M.obj = pe.Objective(rule=ul_obj_rule, sense=pe.minimize)
M.L.obj = pe.Objective(rule=ll_obj_rule, sense=pe.maximize)
M.c1 = pe.Constraint(rule=ul_constraint_rule)
M.L.c1 = pe.Constraint(rule=ll_constraint_rule)

solver = Solver('pao.pyomo.FA')
results = solver.solve(M, tee=True)

Since x is the upper-level variable and y is the lower-level variable, the solution should be x=2 and y=5 according to the objectives and constraints (with an objective value of 7).

Though, the solver sets y to 1e-4 with an objective value close to 2, so the upper-level seems to determine not only x, but also y. If I set the upper-bound of y to "None" in the "relevant code line", the expected objective value of 7 is found. For an upper bound of 1e5, the value for y is set to 1 (instead of 5).

Running the same code with the pao.pyomo.PCCG solver or pao.pyomo.REG, this behaviour doesn't occur.

Since my model contains only linear equations, real upper- and lower-level variables and is a bilevel problem, I assume that all requirements to apply the pao.pyomo.FA solver are fulfilled. Can anyone please help me understand why the described (unexpected) results are found?

[dpiendl]

Hello,
I'm not a contributor to PAO, so I unfortunately cannot help you why this behavior occurs. However, I also solved this problem using
solver = Solver("pao.pyomo.FA", mip_solver="cplex") or
solver = Solver("pao.pyomo.FA", mip_solver="cbc")
and both yield the correct result of x=2 and y=5. So I guess the problem occurs somewhere along where FA passes the MILP to GLPK, which is the default MIP solver for FA.