diff --git a/Project.toml b/Project.toml index 1d0bc014..d85ea982 100644 --- a/Project.toml +++ b/Project.toml @@ -4,10 +4,15 @@ authors = ["Joaquim Garcia , guilhermebodin x==idx, idxs)) # findfirst can return nothing + end + return idxs +end + + +# switch input of I,J for finding cs +function non_zero_idxs_except_one_IJV(I::Vector{Int}, J::Vector{Int}, col::Int, idx::Int)::Vector{Int} + rs = PushVector{Int}() + for (i,v) in enumerate(J) + if v == col && I[i] != idx + push!(rs, I[i]) + end + end + return finish!(rs) +end + + +struct UnderDeterminedException <: Exception end + + +""" + recursive_col_search(A::AbstractArray, row::Int, col::Int, rows::AbstractVector{Int}, cols::AbstractVector{Int}) + +Given a row, col in an array A find all the connected rows and cols starting with searching col, +where "connected" rows and cols are those with non-zero values. +Search is stopped if any redundant rows or cols are added to the output arrays +(which indicates that there is a loop in the connections and the system of equations is underdetermined). + +Returns Vector{Int}, Vector{Int}, Bool where the Bool indicates if redundanct rows or columns were found +(true indicates underdertermined system). +""" +function recursive_col_search(A::AbstractArray, row::Int, col::Int, + rows::AbstractVector{Int}, cols::AbstractVector{Int}) + rs = non_zero_idxs_except_one(A[:, col], row) + if any(r in rows for r in rs) + rr = intersect(rs, rows) + @debug("Returning early from recursive_col_search due to redundant row(s)! ($rr)") + throw(UnderDeterminedException()) + end + push!(rows, rs...) + for r in rs + cs = non_zero_idxs_except_one(A[r, :], col) + if any(c in cols for c in cs) + cc = intersect(cs, cols) + @debug("Returning early from recursive_col_search due to redundant column(s)! ($cc)") + throw(UnderDeterminedException()) + end + push!(cols, cs...) + for c in cs + recursive_col_search(A, r, c, rows, cols) + end + end + return finish!(rows), finish!(cols) +end + +# a version of recursive_col_search that works with I,J,vals = findnz(V) +function recursive_col_search_IJV(I::Vector{Int}, J::Vector{Int}, row::Int, col::Int, + rows::PushVector{Int}, cols::PushVector{Int})::Tuple{Vector{Int}, Vector{Int}} + + rs = non_zero_idxs_except_one_IJV(I, J, col, row) + # rs = non_zero_idxs_except_one(A[:, col], row) + if any(r in rows for r in rs) + rr = intersect(rs, rows) + @debug("Returning early from recursive_col_search due to redundant row(s)! ($rr)") + throw(UnderDeterminedException()) + end + push!(rows, rs...) + for r in rs + cs = non_zero_idxs_except_one_IJV(J, I, r, col) + # cs = non_zero_idxs_except_one(A[r, :], col) + if any(c in cols for c in cs) + cc = intersect(cs, cols) + @debug("Returning early from recursive_col_search due to redundant column(s)! ($cc)") + throw(UnderDeterminedException()) + end + push!(cols, cs...) + for c in cs + recursive_col_search_IJV(I, J, r, c, rows, cols) + end + end + return finish!(rows), finish!(cols) +end + + +""" + find_connected_rows_cols(A::AbstractArray, row::Int, col::Int; skip_1st_col_check=false) + +Given a row, col location to start from in an array, find all the connected rows and columns +where connected is defined by following non-zero values in the array, starting with all non-zero +values in the given row (except the starting row, col) + +Return Int, Int for rows to add primal equations, and columns to add dual equations + +A = [1 0 0 1 0 0; + 0 0 1 1 0 0; + 0 0 1 0 1 0; + 0 1 1 0 0 0; + 1 0 0 0 0 1] + +find_connected_rows_cols(A, 1, 1) +([1, 2, 3, 4], [4, 3, 5, 2]) + +Example 5.2 in paper: +using BlockArrays +V = BlockArray{Float64}(undef_blocks, [3,3], [3,3,3,1,3]) +setblock!(V, Matrix{Float64}(-I, 3, 3), 1, 1) +setblock!(V, Matrix{Float64}( I, 3, 3), 1, 2) +setblock!(V, Matrix{Float64}( I, 3, 3), 1, 3) +setblock!(V, Matrix{Float64}( I, 3, 3), 2, 3) +setblock!(V, Matrix{Float64}( I, 3, 3), 2, 5) +setblock!(V, zeros(3,3), 2, 1) +setblock!(V, zeros(3,3), 2, 2) +setblock!(V, zeros(3,3), 1, 5) +setblock!(V, zeros(3,1), 1, 4) +setblock!(V, ones(3,1), 2, 4) + +find_connected_rows_cols(V, 1, 1; skip_1st_col_check=true) +([1, 4, 5, 6, 2, 3], [4, 7, 10, 11, 8, 12, 2, 5, 9, 13, 3, 6]) +""" +function find_connected_rows_cols(A::AbstractArray, row::Int, col::Int; + skip_1st_col_check=false, + finding_blocks=false, + ) + if A[row, col] == 0 + throw(@error("Linearization is undefined when the dual variable is not associated with the primal variable.")) + end + redundant_vals = false + # step 1 check if all non-zeros in A[:, col], if so the dual constraint gives linearization + if !skip_1st_col_check && length(findall(!iszero, A[:, col])) == 1 + return [], [col], redundant_vals + end + # step 2 add 1st row and any other non-zero columns + rows = PushVector{Int}() + push!(rows, row) + if finding_blocks + cols_to_check = findall(!iszero, A[row, :]) + else + cols_to_check = non_zero_idxs_except_one(A[row, :], col) + end + cols = PushVector{Int}() + push!(cols, cols_to_check...) + # step 3 recursive search to find all connections + rows_to_add, cols_to_add = Int[], Int[] + for c in cols_to_check + try + rows_to_add, cols_to_add = recursive_col_search(A, row, c, PushVector{Int}(), PushVector{Int}()) + catch e + if isa(e, UnderDeterminedException) + if finding_blocks # then we still need to add the connected rows and cols + for r in non_zero_idxs_except_one(A[:, c], row) + push!(rows, r) + push!(cols, findall(!iszero, A[r, :])...) + end + rows = unique(rows) # unique! does not work with PushVector + cols = unique(cols) + else + redundant_vals = true + end + break + else rethrow(e) + end + end + push!(rows, rows_to_add...) + push!(cols, cols_to_add...) + end + + return finish!(rows), finish!(cols), redundant_vals +end + + +function find_connected_rows_cols_cached(A, I, J, row::Int, col::Int; + skip_1st_col_check=false, + finding_blocks=false, + store_cache=true + )::Tuple{Vector{Int}, Vector{Int}, Bool} + + if (row, col, skip_1st_col_check, finding_blocks) ∉ keys(cache) + if A[row, col] == 0 + throw(@error("Linearization is undefined when the dual variable is not associated with the primal variable.")) + end + redundant_vals = false + # step 1 check if all non-zeros in A[:, col], if so the dual constraint gives linearization + if !skip_1st_col_check && length(findall(!iszero, A[:, col])) == 1 # TODO an IJV method for this? + return [], [col], redundant_vals + end + # step 2 add 1st row and any other non-zero columns + rows = Int[] + push!(rows, row) + if finding_blocks + cols_to_check = findall(!iszero, A[row, :]) # TODO an IJV method for this? + else + cols_to_check = non_zero_idxs_except_one_IJV(J, I, row, col) #non_zero_idxs_except_one(A[row, :], col) + end + cols = Int[] + push!(cols, cols_to_check...) + # step 3 recursive search to find all connections + rows_to_add, cols_to_add = Int[], Int[] + for c in cols_to_check + try + rows_to_add, cols_to_add = recursive_col_search_IJV(I, J, row, c, PushVector{Int}(), PushVector{Int}()) + catch e + if isa(e, UnderDeterminedException) + if finding_blocks # then we still need to add the connected rows and cols + for r in non_zero_idxs_except_one_IJV(I, J, c, row) #non_zero_idxs_except_one(A[:, c], row) + push!(rows, r) + push!(cols, findall(!iszero, A[r, :])...) # TODO an IJV method for this? + end + rows = unique(rows) # unique! does not work with PushVector + cols = unique(cols) + else + redundant_vals = true + end + break + else rethrow(e) + end + end + push!(rows, rows_to_add...) + push!(cols, cols_to_add...) + end + if !store_cache + # when not checking linearization conditions, storing in the cache unnecessarily consumes memory + return rows, cols, redundant_vals + end + cache[(row, col, skip_1st_col_check, finding_blocks)] = rows, cols, redundant_vals + end + return cache[(row, col, skip_1st_col_check, finding_blocks)] +end + + +function get_coef(var::MOI.VariableIndex, safts::Vector{MOI.ScalarAffineTerm{R}}) where R <: Real + coef = 0 + for saft in safts + if var == saft.variable + coef = saft.coefficient + end + end + return coef +end + + +function get_coef(var1::MOI.VariableIndex, var2::MOI.VariableIndex, sqts::Vector{MOI.ScalarQuadraticTerm{R}}) where R <: Real + coef = 0 + for sqft in sqts + if var1 == sqft.variable_1 && var2 == sqft.variable_2 + coef = sqft.coefficient + end + end + return coef +end + + +function get_coef_matrix_and_rhs_vec(m, + constraint_indices::Array{ + MOI.ConstraintIndex{ + MOI.ScalarAffineFunction{Float64}, + MOI.GreaterThan{Float64} + }, 1} + ) + nrows = length(constraint_indices) + E = spzeros(nrows, MOI.get(m, MOI.NumberOfVariables())) + f = -Inf*ones(nrows) + + for (r, ci) in enumerate(constraint_indices) + con_func = MOI.get(m, MOI.ConstraintFunction(), ci) + for term in con_func.terms + E[r, term.variable.value] = term.coefficient + end + f[r] = MOI.get(m, MOI.ConstraintSet(), ci).lower + end + return E, f +end + + +# can replace the following SingleVariable limits with an MOI function for getting bounds? +function get_coef_matrix_and_rhs_vec(m, + constraint_indices::Array{ + MOI.ConstraintIndex{ + MOI.VariableIndex, + MOI.GreaterThan{Float64} + }, 1} + ) + f = -Inf*ones(MOI.get(m, MOI.NumberOfVariables())) + + for ci in constraint_indices + var_index = MOI.get(m, MOI.ConstraintFunction(), ci) + f[var_index.value] = MOI.get(m, MOI.ConstraintSet(), ci).lower + end + return f +end + + +function get_coef_matrix_and_rhs_vec(m, + constraint_indices::Array{MOI.ConstraintIndex{ + MOI.VariableIndex, + MOI.LessThan{Float64}}, + 1} + ) + d = Inf*ones(MOI.get(m, MOI.NumberOfVariables())) + + for ci in constraint_indices + var_index = MOI.get(m, MOI.ConstraintFunction(), ci) + d[var_index.value] = MOI.get(m, MOI.ConstraintSet(), ci).upper + end + return d +end + + +function get_coef_matrix_and_rhs_vec(m, + constraint_indices::Array{ + MOI.ConstraintIndex{ + MOI.ScalarAffineFunction{Float64}, + MOI.LessThan{Float64}}, + 1} + ) + nrows = length(constraint_indices) + C = spzeros(nrows, MOI.get(m, MOI.NumberOfVariables())) + d = Inf*ones(nrows) + for (r, ci) in enumerate(constraint_indices) + con_func = MOI.get(m, MOI.ConstraintFunction(), ci) + for term in con_func.terms + C[r, term.variable.value] = term.coefficient + end + d[r] = MOI.get(m, MOI.ConstraintSet(), ci).upper + end + return C, d +end + + +function get_coef_matrix_and_rhs_vec(m, + constraint_indices::Array{ + MOI.ConstraintIndex{ + MOI.ScalarAffineFunction{Float64}, + MOI.EqualTo{Float64}}, + 1} + ) + nrows = length(constraint_indices) + V = spzeros(nrows, MOI.get(m, MOI.NumberOfVariables())) + w = spzeros(nrows) + for (r, ci) in enumerate(constraint_indices) + con_func = MOI.get(m, MOI.ConstraintFunction(), ci) + con_set = MOI.get(m, MOI.ConstraintSet(), ci) + for term in con_func.terms + V[r, term.variable.value] = term.coefficient + end + w[r] = con_set.value + end + return V, w +end + + +""" + standard_form(m; upper_var_indices=Vector{MOI.VariableIndex}()) + +Given: + + A[x:y] = b, C[x;y] ≤ d, E[x;y] ≥ f + +collect all of y bounds into: + + yl ≤ y ≤ yu + +and put all inequality constraints into equalities with slack variables s: + + Ux + V [y; s] = [b; d; f] + +NOTE: U and V are sparse arrays with columns for all variables in model s.t. that variable indices line up +""" +function standard_form(m; upper_var_indices=Vector{MOI.VariableIndex}()) + nvars = MOI.get(m, MOI.NumberOfVariables()) + con_types = MOI.get(m, MOI.ListOfConstraintTypesPresent()) + + n_equality_cons = 0 # A[x;y] = b + + if (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) in con_types + + eq_con_indices = MOI.get(m, MOI.ListOfConstraintIndices{ + MOI.ScalarAffineFunction{Float64}, + MOI.EqualTo{Float64} + }()); + + A, b = get_coef_matrix_and_rhs_vec(m, eq_con_indices) + n_equality_cons = size(A,1) + else + A, b = spzeros(0, nvars), spzeros(0) + end + # U = part of A for upper_var_indices + # V = rest of A + + #= + For each ScalarAffineFunction ≥ or ≤ Float64 we add a slack variable and make the constraint + an EqualTo{Float64} + =# + n_lessthan_cons = 0 # C[x;y] ≤ d + + if (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) in con_types + + lt_con_indices = MOI.get(m, MOI.ListOfConstraintIndices{ + MOI.ScalarAffineFunction{Float64}, + MOI.LessThan{Float64} + }()); + + C, d = BilevelJuMP.get_coef_matrix_and_rhs_vec(m, lt_con_indices) + else + C, d = spzeros(0, nvars), spzeros(0) + end + + n_greaterthan_cons = 0 # E[x;y] ≥ f + + if (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) in con_types + + gt_con_indices = MOI.get(m, MOI.ListOfConstraintIndices{ + MOI.ScalarAffineFunction{Float64}, + MOI.GreaterThan{Float64} + }()); + + E, f = BilevelJuMP.get_coef_matrix_and_rhs_vec(m, gt_con_indices) + else + E, f = spzeros(0, nvars), spzeros(0) + end + + #= + For each SingleVariable GreaterThan{Float64} or LessThan{Float64} we fill in the bounds + yl and yu + =# + + yl = -Inf*ones(MOI.get(m, MOI.NumberOfVariables())) + if (MOI.VariableIndex, MOI.GreaterThan{Float64}) in con_types + + singleVar_gt_indices = MOI.get(m, MOI.ListOfConstraintIndices{ + MOI.VariableIndex, + MOI.GreaterThan{Float64} + }()); + + yl = BilevelJuMP.get_coef_matrix_and_rhs_vec(m, singleVar_gt_indices) + end + + yu = Inf*ones(MOI.get(m, MOI.NumberOfVariables())) + if (MOI.VariableIndex, MOI.LessThan{Float64}) in con_types + + singleVar_lt_indices = MOI.get(m, MOI.ListOfConstraintIndices{ + MOI.VariableIndex, + MOI.LessThan{Float64} + }()); + + yu = get_coef_matrix_and_rhs_vec(m, singleVar_lt_indices) + end + + # remove rows from C that only apply to one variable by moving them to the bounds in yu + + rows_to_remove = Int[] + for r in 1:size(C,1) + if length(findall(!iszero, C[r, :])) == 1 # only one non-zero value in row + # then we remove the row and move it to the bounds (if the RHS is lower/higher than the upper/lower bound) + push!(rows_to_remove, r) + ci = findall(!iszero, C[r, :])[1] + upbound = d[r] / C[r,ci] # C [x;y] ≤ d + if upbound < yu[ci] + yu[ci] = upbound + end + end + end + C = C[setdiff(1:end, rows_to_remove),:] + d = d[setdiff(1:end, rows_to_remove)] + n_lessthan_cons = size(C,1) + + # remove rows from E that only apply to one variable by moving them to the bounds in yl + + rows_to_remove = Int[] + for r in 1:size(E,1) + if length(findall(!iszero, E[r, :])) == 1 # only one non-zero value in row + # then we remove the row and move it to the bounds (if the RHS is lower/higher than the upper/lower bound) + push!(rows_to_remove, r) + ci = findall(!iszero, E[r, :])[1] + lobound = f[r] / E[r,ci] + if lobound > yl[ci] + yl[ci] = lobound + end + end + end + E = E[setdiff(1:end, rows_to_remove),:] + f = f[setdiff(1:end, rows_to_remove)] + n_greaterthan_cons = size(E,1) + + + #= + Build V = | A | + | C I | + | E I | + =# + + n_vars = size(A,2) + n_rows_V = n_equality_cons + n_greaterthan_cons + n_lessthan_cons + n_cols_V = n_vars + n_greaterthan_cons + n_lessthan_cons + V = spzeros(n_rows_V, n_cols_V) + V[1:n_equality_cons, 1:n_vars] = A + GC.gc() + + V[n_equality_cons+1 : n_equality_cons+n_lessthan_cons, 1 : n_vars] = C + GC.gc() + V[n_equality_cons+1 : n_equality_cons+n_lessthan_cons, n_vars+1 : n_vars+n_lessthan_cons] = Matrix(I, n_lessthan_cons, n_lessthan_cons) + GC.gc() + + V[n_equality_cons+n_lessthan_cons+1 : n_equality_cons+n_greaterthan_cons+n_lessthan_cons, 1:n_vars] = E + GC.gc() + V[n_equality_cons+n_lessthan_cons+1 : n_equality_cons+n_greaterthan_cons+n_lessthan_cons, + n_vars+n_lessthan_cons+1 : n_vars+n_greaterthan_cons+n_lessthan_cons] = Matrix(-I, n_greaterthan_cons, n_greaterthan_cons) + GC.gc() + + # zero out the columns in V for upper level variables and build U + + rows, cols, vals = findnz(V) + Urows = PushVector{Int}() + Ucols = PushVector{Int}() + Uvals = PushVector{Float64}() + + for col in upper_var_indices + indices_to_mv = findall(i->i==col.value, cols) + push!(Urows, [rows[i] for i in indices_to_mv]...) + push!(Ucols, [cols[i] for i in indices_to_mv]...) + push!(Uvals, [vals[i] for i in indices_to_mv]...) + deleteat!(rows, indices_to_mv) + deleteat!(cols, indices_to_mv) + deleteat!(vals, indices_to_mv) + end + GC.gc() + + V = sparse(rows, cols, vals, n_rows_V, n_cols_V) + U = sparse(finish!(Urows), finish!(Ucols), finish!(Uvals), n_rows_V, n_cols_V) + # every time a new V is created for a problem we need to empty the cache used in find_connected_rows_cols_cached + empty!(cache) + w = [b; d; f] + + return U, V, w # , yu, yl, n_equality_cons, C, E + # TODO use n_equality_cons to check rows from find_connected_rows_cols for values corresponding to constraints with slack variables +end + + +""" + check_upper_objective_for_bilinear_linearization(upper, upper_to_lower_var_indices, upper_var_to_lower_ctr) + +Construct the set A_N, which is the indices of lower level variables in the upper level +objective of the form λ^T A y, where λ are the dual variables of lower level equality constraints and +y are lower level primal variables. + +Also find the quadratic upper objective terms with lower primal * lower dual and return the maps for +the upper variable index (of the lower dual variable): +- to the lower primal term and +- to the quadratic term in the UL objective + +""" +function check_upper_objective_for_bilinear_linearization(upper, upper_to_lower_var_indices, upper_var_to_lower_ctr) + nt = Threads.nthreads() + + A_N = Vector{PushVector{Int}}() + upper_dual_to_quad_term = Vector{Dict{BilevelJuMP.VI, Dict{BilevelJuMP.VI, Float64}}}() + upper_dual_to_lower_primal = Vector{Dict{BilevelJuMP.VI, Vector{BilevelJuMP.VI}}}() + lower_primal_var_to_lower_con = Vector{Dict{BilevelJuMP.VI, BilevelJuMP.CI}}() + + for _ in 1:nt + push!(A_N, PushVector{Int}()) + push!(upper_dual_to_quad_term, Dict{BilevelJuMP.VI, Dict{BilevelJuMP.VI, Float64}}()) + push!(upper_dual_to_lower_primal, Dict{BilevelJuMP.VI, Vector{BilevelJuMP.VI}}()) + push!(lower_primal_var_to_lower_con, Dict{BilevelJuMP.VI, BilevelJuMP.CI}()) + end + + # upper_primal_to_lower_primal = Dict{BilevelJuMP.VI, BilevelJuMP.VI}() + + UL_obj_type = MOI.get(upper, MOI.ObjectiveFunctionType()) + upper_obj_func_quad_terms = MOI.get( + upper, MOI.ObjectiveFunction{UL_obj_type}()).quadratic_terms + + Threads.@threads for upper_dual_var_idx in collect(eachindex(upper_var_to_lower_ctr)) + id = Threads.threadid() + for term in upper_obj_func_quad_terms + + if upper_dual_var_idx == term.variable_1 + if term.variable_2 in values(upper_to_lower_var_indices) + lower_primal_var_idx = upper_to_lower_var_indices[term.variable_2] + push!(A_N[id], lower_primal_var_idx.value) + # upper_primal_to_lower_primal[term.variable_2] = lower_primal_var_idx + # upper_dual_to_quad_term[upper_dual_var_idx] = [term] # will overwrite entry if upper_dual is multiplied by more than one lower variable in the upper objective, needs to list of pairs, muliple cases here + # upper_dual_to_quad_term is only used to find A_jn, the coefs in the upper obj of lambda_j and yn, could instead make matrix + # upper_dual_to_quad_term[upper_dual_var_idx][upper_to_lower_var_indices[term.variable_2]] = term.coefficient, which changes to variable_1 in elseif below + if !haskey(upper_dual_to_quad_term[id], upper_dual_var_idx) + upper_dual_to_quad_term[id][upper_dual_var_idx] = Dict( + upper_to_lower_var_indices[term.variable_2] => term.coefficient + + ) + else + upper_dual_to_quad_term[id][upper_dual_var_idx][upper_to_lower_var_indices[term.variable_2]] = term.coefficient + end + if !haskey(upper_dual_to_lower_primal[id], upper_dual_var_idx) + upper_dual_to_lower_primal[id][upper_dual_var_idx] = [lower_primal_var_idx] + else + push!(upper_dual_to_lower_primal[id][upper_dual_var_idx], lower_primal_var_idx) + end + lower_primal_var_to_lower_con[id][lower_primal_var_idx] = upper_var_to_lower_ctr[upper_dual_var_idx] + end + + elseif upper_dual_var_idx == term.variable_2 + if term.variable_1 in values(upper_to_lower_var_indices) + lower_primal_var_idx = upper_to_lower_var_indices[term.variable_1] + push!(A_N[id], lower_primal_var_idx.value) + # upper_primal_to_lower_primal[term.variable_1] = lower_primal_var_idx + if !haskey(upper_dual_to_quad_term[id], upper_dual_var_idx) + upper_dual_to_quad_term[id][upper_dual_var_idx] = Dict( + upper_to_lower_var_indices[term.variable_1] => term.coefficient + + ) + else + upper_dual_to_quad_term[id][upper_dual_var_idx][upper_to_lower_var_indices[term.variable_1]] = term.coefficient + end + if !haskey(upper_dual_to_lower_primal[id], upper_dual_var_idx) + upper_dual_to_lower_primal[id][upper_dual_var_idx] = [lower_primal_var_idx] + else + push!(upper_dual_to_lower_primal[id][upper_dual_var_idx], lower_primal_var_idx) + end + lower_primal_var_to_lower_con[id][lower_primal_var_idx] = upper_var_to_lower_ctr[upper_dual_var_idx] + end + end + end + end + + upper_dual_to_quad_term = mergewith(merge, upper_dual_to_quad_term...) + upper_dual_to_lower_primal = mergewith(vcat, upper_dual_to_lower_primal...) + lower_primal_var_to_lower_con = merge(lower_primal_var_to_lower_con...) + return unique(vcat(finish!.(A_N)...)), upper_dual_to_quad_term, upper_dual_to_lower_primal, lower_primal_var_to_lower_con +end + + +function get_lower_obj_coefs_of_upper_times_lower_primals( + lower, + lower_var_indices_of_upper_vars, + A_N, + lower_to_m_idxmap, + ) + AB_N = Int[] + nvars = MOI.get(lower, MOI.NumberOfVariables()) + # B is in LL objective: x^T B y + B = spzeros(nvars, nvars) + + ks = collect(keys(lower_to_m_idxmap)) + lower_var_indices = collect(filter(k -> typeof(k)==MOI.VariableIndex , ks)) + lower_only_vars = setdiff(lower_var_indices, lower_var_indices_of_upper_vars) + + # get lower objective terms for finding cost coefficients later + lower_obj_type = MOI.get(lower, MOI.ObjectiveFunctionType()) + lower_obj = MOI.get(lower, MOI.ObjectiveFunction{lower_obj_type}()) + lower_obj_quad_terms = nothing + lower_obj_terms = nothing + lower_obj_type_handled = true + + if lower_obj_type == MOI.ScalarQuadraticFunction{Float64} + lower_obj_quad_terms = lower_obj.quadratic_terms + lower_obj_terms = lower_obj.affine_terms + elseif lower_obj_type == MOI.ScalarAffineFunction{Float64} + lower_obj_terms = lower_obj.terms + else + lower_obj_type_handled = false + end + + # fill the set AB_N = {n in A_N : ∃ m ∈ M s.t. B_mn ≠ 0} and the B matrix (x^T B y in LL objective) + # NOTE B is assumed to index by upper, lower var index (m,n) indices, + # (which are not necessarily in the order entered by user). + if !isnothing(lower_obj_quad_terms) # check for values in AB_N, + for term in lower_obj_quad_terms + if term.variable_1 in lower_var_indices_of_upper_vars # UL var + if term.variable_2.value in A_N # LL var + push!(AB_N, term.variable_2.value) # AB_N is not empty + end + if term.variable_2 in lower_only_vars + B[term.variable_1.value, term.variable_2.value] = + get_coef(term.variable_1, term.variable_2, lower_obj_quad_terms) + end + elseif term.variable_2 in lower_var_indices_of_upper_vars # UL var + if term.variable_1.value in A_N # LL var + push!(AB_N, term.variable_1.value) # AB_N is not empty + end + if term.variable_1 in lower_only_vars + B[term.variable_2.value, term.variable_1.value] = + get_coef(term.variable_1, term.variable_2, lower_obj_quad_terms) + end + end + end + end + return AB_N, B, lower_obj_terms, lower_obj_type_handled +end + + +function linear_terms_for_empty_AB( + lower, + upper_var_to_lower_ctr, + bilinear_upper_dual_to_lower_primal, + V, + w, + bilinear_upper_dual_to_quad_term, + lower_obj_terms, + lower_to_m_idxmap, + lower_primal_dual_map, + lower_dual_idxmap + ) + linearizations = Vector{MOI.ScalarAffineTerm}() + con_type = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}} + I, J, vals = findnz(V) + + for (upper_var, lower_con) in upper_var_to_lower_ctr # TODO thread + j = lower_con.value + for lower_var in bilinear_upper_dual_to_lower_primal[upper_var] + n = lower_var.value + + J_j, N_n, redundant_vals = find_connected_rows_cols_cached(V, I, J, j, n, skip_1st_col_check=false) + if redundant_vals + return nothing + end + + A_jn = bilinear_upper_dual_to_quad_term[upper_var][lower_var] + V_jn = V[j,n] + + for j_prime in J_j + lower_con_index = con_type(j_prime) + lower_dual_var = lower_primal_dual_map.primal_con_dual_var[lower_con_index][1] + + push!(linearizations, + MOI.ScalarAffineTerm(A_jn / V_jn * w[j_prime], lower_dual_idxmap[lower_dual_var]) + ) + end + + # TODO assert that lower level constraints in upper_var_to_lower_ctr are linear + + num_vars = MOI.get(lower, MOI.NumberOfVariables()) + for n_prime in N_n + if n_prime > num_vars continue end # TODO do we need to add slack variables? + lower_var = MOI.VariableIndex(n_prime) + # lower primal * lower cost + lower_var_cost_coef = BilevelJuMP.get_coef(lower_var, lower_obj_terms) + push!(linearizations, + MOI.ScalarAffineTerm(-A_jn / V_jn * lower_var_cost_coef, lower_to_m_idxmap[lower_var]) + ) + # variable bound * dual variable + low_bound, upp_bound = MOIU.get_bounds(lower, Float64, lower_var) # yl[n_prime], yu[n_prime] # + lo_bound_index = MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}}(lower_var.value) + up_bound_index = MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}}(lower_var.value) + low_dual = get(lower_primal_dual_map.primal_con_dual_var, lo_bound_index, [nothing])[1] + upp_dual = get(lower_primal_dual_map.primal_con_dual_var, up_bound_index, [nothing])[1] + + # have to use opposite signs of paper for these terms (b/c Dualization sets variable bound dual variables to be non-positive?) + if low_bound != -Inf && !isnothing(low_dual) + push!(linearizations, MOI.ScalarAffineTerm(-A_jn / V_jn * low_bound, lower_dual_idxmap[low_dual])) + end + if upp_bound != Inf && !isnothing(upp_dual) # TODO add a big number in place of Inf ? + push!(linearizations, MOI.ScalarAffineTerm( A_jn / V_jn * upp_bound, lower_dual_idxmap[upp_dual])) + end + end + end + end + + return linearizations +end + + +function linear_terms_for_non_empty_AB( + lower, + upper_var_to_lower_ctr, + bilinear_upper_dual_to_lower_primal, + V, + w, + A_N, + bilinear_upper_dual_to_quad_term, + lower_obj_terms, + lower_to_m_idxmap, + lower_primal_dual_map, + lower_dual_idxmap, + store_cache, + num_blocks, + rows + )::Vector{MOI.ScalarAffineTerm} + # TODO add store_cache to empty AB functions + + linearizations = PushVector{MOI.ScalarAffineTerm}() + + con_type = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}} + I, J, _ = findnz(V) + + # only need one (j,n) pair for the linearization! (equ. 23 or 24 in https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=9729553) + # upper_var_to_lower_ctr is equivalent to set A with pairs (j,n) : A_jn ≠ 0 + + for block_int in 1:num_blocks + + # find a j in this block, i.e. a row in V or a constraint index in the lower problem + j = nothing + upper_var = nothing + for (uv, lower_con) in upper_var_to_lower_ctr + if lower_con.value in rows[block_int] + j = lower_con.value + upper_var = uv + continue + end + end + if isnothing(j) + return nothing + end + + for lower_var in bilinear_upper_dual_to_lower_primal[upper_var] # typically just one lower_var + + A_jn = bilinear_upper_dual_to_quad_term[upper_var][lower_var] + if isapprox(A_jn, 0.0; atol=1e-12) + continue + end + n = lower_var.value + # this call is same as in get_all_connected_rows_cols, hence memoization should speed things up + J_j, N_n, redundant_vals = find_connected_rows_cols_cached( + V, I, J, j, n, skip_1st_col_check=true, store_cache=store_cache + ) + if redundant_vals + return nothing + end + + V_jn = V[j,n] + p = A_jn / V_jn + for r in J_j + if isapprox(w[r], 0.0; atol=1e-12) + continue + end + lower_con_index = con_type(r) + lower_dual_var = lower_primal_dual_map.primal_con_dual_var[lower_con_index][1] + + push!(linearizations, + MOI.ScalarAffineTerm(p*w[r], lower_dual_idxmap[lower_dual_var]) + ) + end + + # TODO assert that lower level constraints in upper_var_to_lower_ctr are linear + N_n = setdiff(N_n, A_N) + num_vars = MOI.get(lower, MOI.NumberOfVariables()) + for c in N_n + if c > num_vars continue end # TODO do we need to add slack variables? + lower_var = MOI.VariableIndex(c) + # lower primal * lower cost + lower_var_cost_coef = BilevelJuMP.get_coef(lower_var, lower_obj_terms) + if !isapprox(lower_var_cost_coef, 0.0; atol=1e-12) + push!(linearizations, + MOI.ScalarAffineTerm(-p*lower_var_cost_coef, lower_to_m_idxmap[lower_var]) + ) + end + # variable bound * dual variable + low_bound, upp_bound = MOIU.get_bounds(lower, Float64, lower_var) # yl[n_prime], yu[n_prime] # + lo_bound_index = MOI.ConstraintIndex{MOI.VariableIndex, MOI.GreaterThan{Float64}}(lower_var.value) + up_bound_index = MOI.ConstraintIndex{MOI.VariableIndex, MOI.LessThan{Float64}}(lower_var.value) + low_dual = get(lower_primal_dual_map.primal_con_dual_var, lo_bound_index, [nothing])[1] + upp_dual = get(lower_primal_dual_map.primal_con_dual_var, up_bound_index, [nothing])[1] + + # have to use opposite signs of paper for these terms (b/c Dualization sets variable bound dual variables to be non-positive?) + if low_bound != -Inf && !isapprox(low_bound, 0.0; atol=1e-12) && !isnothing(low_dual) + + push!(linearizations, MOI.ScalarAffineTerm(-p * low_bound, lower_dual_idxmap[low_dual])) + end + if upp_bound != Inf && !isapprox(upp_bound, 0.0; atol=1e-12) && !isnothing(upp_dual) # TODO add a big number in place of Inf ? + push!(linearizations, MOI.ScalarAffineTerm( p * upp_bound, lower_dual_idxmap[upp_dual])) + end + end + end + end + + return linearizations +end + + +function check_condition_1(J_U::AbstractVector{Int}, U::AbstractMatrix) + # Condition 1: U_jm = 0 ∀ j ∈ J_U, ∀ m ∈ M + met_condition = true + for j in J_U, m in 1:size(U,2) + if U[j,m] ≠ 0 + met_condition = false + @warn("Condition 1 not met: at least one connected lower constraint contains upper level variables.") + break + end + end + met_condition +end + + +function check_condition_2(N_U::AbstractVector{Int}, U::AbstractMatrix, B::AbstractMatrix) + # Condition 2: B_mn = 0 ∀ m ∈ M, ∀ n ∈ N_U + met_condition = true + for m in 1:size(U,1), n in N_U + if B[m,n] ≠ 0 + met_condition = false + @warn("Condition 2 not met: at least one connected lower variable is multiplied with an upper variable in the lower objective.") + break + end + end + met_condition +end + + +function check_condition_2prime( + N_U::AbstractVector{Int}, A_N::AbstractVector{Int}, + U::AbstractMatrix, B::AbstractMatrix + ) + # Condition 2': B_mn = 0 ∀ m ∈ M, ∀ n ∈ N_U \ A_N + met_condition = true + for m in 1:size(U,2), n in setdiff(N_U, A_N) + if B[m,n] ≠ 0 + met_condition = false + @warn("Condition 2' not met: at least one connected lower variable (that is not in the upper objective with a lower dual) is multiplied with an upper variable in the lower objective.") + break + end + end + met_condition +end + + +function check_condition_3(A_N::AbstractVector{Int}, V::AbstractMatrix, lower_primal_var_to_lower_con) + # Condition 3: A_N \ n ⊆ N_n ∀ n ∈ A_n + I, J, vals = findnz(V) + met_condition = true + Threads.@threads for n in A_N + j = lower_primal_var_to_lower_con[MOI.VariableIndex(n)].value + # this call is same as in get_all_connected_rows_cols, hence memoization should speed things up + _, N_n, _ = find_connected_rows_cols_cached(V, I, J, j, n, skip_1st_col_check=true) + # NOTE not handling redundant_vals here b/c goal is to first check conditions + # (redundant_vals handled when determining linearizations) + if !(issubset(setdiff(A_N, n), N_n)) + met_condition = false + @warn("Condition 3 not met: at least one lower variable from the upper objective bilinear terms is not connected to the other lower variables in the upper bilinear terms.") + break + end + end + met_condition +end + + +function check_condition_4(A_N::AbstractVector{Int}, V::AbstractMatrix, upper_var_to_lower_ctr, bilinear_upper_dual_to_lower_primal) + # Condition 4: V_j'n = 0 ∀ j' ∈ J \ {j}, ∀ (j,n) ∈ A + met_condition = true + J = 1:size(V,1) + for (upper_var, lower_con) in upper_var_to_lower_ctr # thread? + j = lower_con.value + for lower_var in bilinear_upper_dual_to_lower_primal[upper_var] + n = lower_var.value + if !(n in A_N) continue end # TODO can we pass in a better set to loop over? + for j_prime in setdiff(J, j) + if V[j_prime, n] ≠ 0 + met_condition = false + @warn("Condition 4 not met: at least one lower variable from the upper objective bilinear terms is in more than one lower constraint.") + @goto outer_break + end + end + end + end + @label outer_break + return met_condition +end + + +function check_condition_5(A_N::AbstractVector{Int}, V::AbstractMatrix, upper_var_to_lower_ctr, bilinear_upper_dual_to_lower_primal, bilinear_upper_dual_to_quad_term) + # Condition 5: A_jn = p V_jn ∀ (j,n) ∈ A (p is proportionality constant) + met_condition = true + p = nothing + for (upper_var, lower_con) in upper_var_to_lower_ctr # thread? + j = lower_con.value + for lower_var in bilinear_upper_dual_to_lower_primal[upper_var] + n = lower_var.value + if !(n in A_N) continue end # TODO can we pass in a better set to loop over? + A_jn = bilinear_upper_dual_to_quad_term[upper_var][lower_var] + V_jn = V[j,n] + if !(isnothing(p)) && !(isapprox(p, A_jn / V_jn, atol=1e-5)) + met_condition = false + @warn("Condition 5 not met: at least one of the ratios of the upper objective bilinear coefficient to lower level constraint coefficient is not equal to the other ratios.") + end + p = A_jn / V_jn + end + end + return met_condition +end + + +""" + check_empty_AB_N_conditions(J_U, U, N_U, B) + +Check two required conditions for linearizing bilinear terms in UL of the form λj * yn when there + are no xm*yn in the LL objective (when AB_N is empty set). +""" +function check_empty_AB_N_conditions(J_U, U, N_U, B) + # Condition 1: U_jm = 0 ∀ j ∈ J_U, ∀ m ∈ M + met_condition_1 = check_condition_1(J_U, U) + + # Condition 2: B_mn = 0 ∀ m ∈ M, ∀ n ∈ N_U + met_condition_2 = check_condition_2(N_U, U, B) + + if met_condition_1 && met_condition_2 + return true + end + + return false +end + + +""" + check_non_empty_AB_N_conditions(J_U, U, N_U, A_N, B, V, lower_primal_var_to_lower_con, + upper_var_to_lower_ctr, bilinear_upper_dual_to_quad_term, bilinear_upper_dual_to_lower_primal) + +Check five required conditions for linearizing bilinear terms in UL of the form λj * yn when there + are xm*yn in the LL objective (when AB_N is not empty). +""" +function check_non_empty_AB_N_conditions(J_U, U, N_U, A_N, B, V, lower_primal_var_to_lower_con, + upper_var_to_lower_ctr, bilinear_upper_dual_to_quad_term, bilinear_upper_dual_to_lower_primal) + # Condition 1: U_jm = 0 ∀ j ∈ J_U, ∀ m ∈ M + met_condition_1 = check_condition_1(J_U, U) + + # Condition 2': B_mn = 0 ∀ m ∈ M, ∀ n ∈ N_U \ A_N + met_condition_2 = check_condition_2prime(N_U, A_N, U, B) + + # Condition 3: A_N \ n ⊆ N_n ∀ n ∈ A_n + met_condition_3 = check_condition_3(A_N, V, lower_primal_var_to_lower_con) + + # Condition 4: V_j'n = 0 ∀ j' ∈ J \ {j}, ∀ (j,n) ∈ A + met_condition_4 = check_condition_4(A_N, V, upper_var_to_lower_ctr, + bilinear_upper_dual_to_lower_primal) + # Condition 5: A_jn = p V_jn ∀ (j,n) ∈ A (p is proportionality constant) + met_condition_5 = check_condition_5(A_N, V, upper_var_to_lower_ctr, + bilinear_upper_dual_to_lower_primal, bilinear_upper_dual_to_quad_term) + + if met_condition_1 && met_condition_2 && met_condition_3 && met_condition_4 && + met_condition_5 + return true + end + + return false +end + + +""" + find_blocks(V::AbstractMatrix{<:Real}, U::AbstractMatrix{<:Real}) + +Find the blocks in the joint matrix [U V] + +For example, given: + +[U V] = | 1 0 0 1 0 0 | + | 0 1 0 0 1 0 | + | 0 1 0 0 1 0 | + +rows = [ + [1], [2,3] +] +cols = [ + [1, 4], [2,5] +] + +NOTE that rows/cols with all zeros are not included in the return values. +""" +function find_blocks(V::AbstractMatrix{<:Real}, U::AbstractMatrix{<:Real}) + + # have to create nthreads arrays and then combine them at the end (to be thread safe) + num_blocks = 0 + rows, cols = PushVector{Vector}(), PushVector{Vector}() + rows_connected = PushVector{Int}() + + # starting with first row of V, find all connected values. + # If all rows are connected then there is only one block. + # Else, search the remaining unconnected rows until all rows are accounted for. + nrows = size(V,1) + UV = U + V + for r in 1:nrows + if r in rows_connected continue end # cannot multithread this for loop b/c of this check + c = findfirst(!iszero, UV[r,:]) + if !isnothing(c) + num_blocks += 1 + rs, cs = find_connected_rows_cols(UV, r, c; skip_1st_col_check=true, finding_blocks=true) + # not worried about redundant_vals here; it is checked when forming linearizations + push!(rows, rs) + push!(cols, cs) + push!(rows_connected, rs...) + end + end + @debug("Found $(num_blocks) block(s) in V.") + return num_blocks, finish!(rows), finish!(cols) +end + + +function is_model_in_standard_form(m::MOI.ModelLike) + model_con_types = Set(MOI.get(m, MOI.ListOfConstraintTypesPresent())) + standard_form_con_types = Set([ + (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}), + (MOI.VariableIndex, MOI.GreaterThan{Float64}), + (MOI.VariableIndex, MOI.LessThan{Float64}) + ]) + issubset(model_con_types, standard_form_con_types) +end + + +function get_all_connected_rows_cols(upper_var_to_lower_ctr, bilinear_upper_dual_to_lower_primal, V, AB_N) + + # have to create nthreads arrays and then combine them at the end (to be thread safe) + J_Us = Vector{PushVector{Int}}() + N_Us = Vector{PushVector{Int}}() + for _ in 1:Threads.nthreads() + push!(J_Us, PushVector{Int}()) + push!(N_Us, PushVector{Int}()) + end + I, J, vals = findnz(V) + + Threads.@threads for upper_var in collect(eachindex(upper_var_to_lower_ctr)) # equivalent to set A with pairs (j,n) : A_jn ≠ 0 + lower_con = upper_var_to_lower_ctr[upper_var] + if !(upper_var in keys(bilinear_upper_dual_to_lower_primal)) continue end # user defined DualOf but did not use it in UL objective + j = lower_con.value + for lower_var in bilinear_upper_dual_to_lower_primal[upper_var] + n = lower_var.value + rows, cols = find_connected_rows_cols_cached(V, I, J, j, n, skip_1st_col_check=!(isempty(AB_N))) + push!(J_Us[Threads.threadid()], rows...) + push!(N_Us[Threads.threadid()], cols...) + end + end + return unique(vcat(finish!.(J_Us)...)), unique(vcat(finish!.(N_Us)...)) +end + + +""" + main_linearization( + lower, + upper, + upper_var_to_lower_ctr, + upper_to_lower_var_indices, + lower_var_indices_of_upper_vars, + lower_to_m_idxmap, + lower_primal_dual_map, + lower_dual_idxmap + ) + +The main logic for linearizing bilinear terms of lower level dual and primal variables in the upper +level model. +""" +function main_linearization( + m, + lower, + upper, + upper_var_to_lower_ctr, + upper_to_lower_var_indices, + lower_var_indices_of_upper_vars, + lower_to_m_idxmap, + upper_to_m_idxmap, + lower_primal_dual_map, + lower_dual_idxmap, + check_linearization_conditions + ) + @debug """starting main linearization at $(Dates.format(now(), "HH:MM:SS"))""" + if MOI.get(upper, MOI.ObjectiveFunctionType()) <: MOI.ScalarQuadraticFunction && + !isempty(upper_var_to_lower_ctr) + + @assert MOI.get(upper, MOI.ObjectiveSense()) == MOI.MIN_SENSE + @assert MOI.get(lower, MOI.ObjectiveSense()) == MOI.MIN_SENSE + # TODO switch signs with MAX sense + + linearize = true + # check lower constraint types and if not just equality and singlevariable bounds then linearize = false and @warn + + @debug """starting is_model_in_standard_form at $(Dates.format(now(), "HH:MM:SS"))""" + if !(is_model_in_standard_form(lower)) + linearize = false + @warn("The lower model must be in standard form to linearize bilinear terms. Skipping linearization process.") + else + @debug """starting check_upper_objective_for_bilinear_linearization at $(Dates.format(now(), "HH:MM:SS"))""" + A_N, bilinear_upper_dual_to_quad_term, bilinear_upper_dual_to_lower_primal, lower_primal_var_to_lower_con = + check_upper_objective_for_bilinear_linearization(upper, upper_to_lower_var_indices, upper_var_to_lower_ctr) + if isempty(A_N) + @debug("No bilinear products of lower level dual and primal variables found in upper level objective. Skipping linearization process.") + linearize = false + else + @debug """starting get_lower_obj_coefs_of_upper_times_lower_primals at $(Dates.format(now(), "HH:MM:SS"))""" + AB_N, B, lower_obj_terms, lower_obj_type_handled = + get_lower_obj_coefs_of_upper_times_lower_primals(lower, lower_var_indices_of_upper_vars, A_N, lower_to_m_idxmap) + if !lower_obj_type_handled + @warn("Linearizing bilinear terms does not handle lower level objective type $(MOI.get(lower, MOI.ObjectiveFunctionType())). Skipping linearization process.") + linearize = false + end + end + end + else + @debug("Upper objective must be quadratic and there must lower dual variables in the upper model to linearize bilinear terms. Skipping linearization process.") + linearize = false + end + + if linearize + @debug """starting standard_form at $(Dates.format(now(), "HH:MM:SS"))""" + U, V, w = standard_form(lower, upper_var_indices=lower_var_indices_of_upper_vars) + @debug """done standard_form at $(Dates.format(now(), "HH:MM:SS"))""" + upper_obj_func_quad_terms = MOI.get(upper, MOI.ObjectiveFunction{MOI.get(upper, MOI.ObjectiveFunctionType())}()).quadratic_terms + linearizations = nothing + m_objective = MOI.get(m, MOI.ObjectiveFunction{MOI.get(m, MOI.ObjectiveFunctionType())}()) + bilinear_upper_quad_term_to_m_quad_term = Dict{MOI.ScalarQuadraticTerm, MOI.ScalarQuadraticTerm}() + + @debug """starting loop over m_objective.quadratic_terms at $(Dates.format(now(), "HH:MM:SS"))""" + # TODO thread or speed this up some how? + for term in m_objective.quadratic_terms + mset = Set([term.variable_1, term.variable_2]) + for upper_term in upper_obj_func_quad_terms + uset = Set([ + upper_to_m_idxmap[upper_term.variable_1], + upper_to_m_idxmap[upper_term.variable_2] + ]) + if uset == mset + bilinear_upper_quad_term_to_m_quad_term[upper_term] = term + end + end + end + + # TODO check for integer x * continuous y, for now assuming continuous x conditions + if check_linearization_conditions + @debug """starting get_all_connected_rows_cols at $(Dates.format(now(), "HH:MM:SS"))""" + J_U, N_U = get_all_connected_rows_cols(upper_var_to_lower_ctr, bilinear_upper_dual_to_lower_primal, V, AB_N) + # J_U and N_U are needed only for checking linearization conditions + end + + #= + Case without x_m * y_n in LL objective for all y_n in A_N (set of bilinear UL objective terms of form λ_j * y_n) + =# + if isempty(AB_N) + @debug("set AB_N is empty") + + conditions_passed = true + if check_linearization_conditions + conditions_passed = check_empty_AB_N_conditions(J_U, U, N_U, B) + end + + if conditions_passed + linearizations = linear_terms_for_empty_AB( + lower, + upper_var_to_lower_ctr, + bilinear_upper_dual_to_lower_primal, + V, + w, + bilinear_upper_dual_to_quad_term, + lower_obj_terms, + lower_to_m_idxmap, + lower_primal_dual_map, + lower_dual_idxmap + ) + if isnothing(linearizations) + @warn("Unable to linearize bilinear terms due to underdertermined system of equations. Skipping linearization process.") + end + else + @warn("Required conditions for linearization not met. Skipping linearization process.") + end + else # AB_N is not empty + @debug("set AB_N is NOT empty") + + # TODO input flag for checking for blocks? (to save time) + + @debug """starting find_blocks at $(Dates.format(now(), "HH:MM:SS"))""" + num_blocks, rows, cols = find_blocks(V, U) + # it is necessary to find the blocks of the lower problem because the linearization + # requires one pair of (j,n) values for each block + + conditions_passed = Bool[] + + if check_linearization_conditions + nrows, ncols = size(V) + + @debug """starting condition checks at $(Dates.format(now(), "HH:MM:SS"))""" + for n in 1:num_blocks + # TODO can skip blocks that are not linked to bilinear terms? + Vblock = spzeros(nrows, ncols) + Vblock[rows[n],:] = V[rows[n],:] + check = check_non_empty_AB_N_conditions( + intersect(J_U, rows[n]), U, intersect(N_U, cols[n]), intersect(A_N, cols[n]), B, Vblock, + lower_primal_var_to_lower_con, upper_var_to_lower_ctr, + bilinear_upper_dual_to_quad_term, + bilinear_upper_dual_to_lower_primal) + push!(conditions_passed, check) + end + @debug("Done looping over V blocks") + else + push!(conditions_passed, true) + end + + if all(conditions_passed) + + @debug """starting linear_terms_for_non_empty_AB at $(Dates.format(now(), "HH:MM:SS"))""" + linearizations = linear_terms_for_non_empty_AB( + lower, + upper_var_to_lower_ctr, + bilinear_upper_dual_to_lower_primal, + V, + w, + A_N, + bilinear_upper_dual_to_quad_term, + lower_obj_terms, + lower_to_m_idxmap, + lower_primal_dual_map, + lower_dual_idxmap, + check_linearization_conditions, + num_blocks, + rows + ) + if isnothing(linearizations) + @warn("Unable to linearize bilinear terms due to underdertermined system of equations. Skipping linearization process.") + end + else + @warn("Required conditions for linearization not met. Skipping linearization process.") + end + + end + + # LINEARIZATION + + if !(isnothing(linearizations)) + + @debug """starting linearizations at $(Dates.format(now(), "HH:MM:SS"))""" + # set m's objective by replacing quadratic terms with linearizations + mobj = deepcopy(m_objective) + quadratic_terms = mobj.quadratic_terms # TODO keep quadratic terms that have not been linearized + c = mobj.constant + affine_terms = mobj.affine_terms + cvb = collect(values(bilinear_upper_quad_term_to_m_quad_term)) + new_objective = deepcopy(m_objective) + if Set(cvb) == Set(quadratic_terms) + @debug("Replacing bilinear lower dual * lower primal terms in upper objective with linear terms.") + new_objective = MOI.ScalarAffineFunction{Float64}( + append!(affine_terms, linearizations), + c + ) + MOI.set(m, MOI.ObjectiveFunction{MOI.ScalarAffineFunction}(), new_objective) + else + @warn("Unable to linearize bilinear terms due to mis-matched quadratic terms in the upper level and single level models. Skipping linearization process.") + end + end + end # if linearize + @debug """done linearizing at $(Dates.format(now(), "HH:MM:SS"))""" +end diff --git a/src/jump.jl b/src/jump.jl index e9982b1c..11ab141b 100644 --- a/src/jump.jl +++ b/src/jump.jl @@ -105,6 +105,8 @@ mutable struct BilevelModel <: AbstractBilevelModel copy_names::Bool copy_names_to_solver::Bool pass_start::Bool + linearize_bilinear_upper_terms::Bool + check_linearization_conditions::Bool # for completing the JuMP.Model API objdict::Dict{Symbol,Any} # Same that JuMP.Model's field `objdict` @@ -156,6 +158,8 @@ mutable struct BilevelModel <: AbstractBilevelModel false, false, true, + false, + true, # jump api Dict{Symbol,Any}(), ) @@ -164,11 +168,15 @@ mutable struct BilevelModel <: AbstractBilevelModel end end function BilevelModel( - optimizer_constructor; - mode::AbstractBilevelSolverMode = SOS1Mode(), - add_bridges::Bool = true, -) + optimizer_constructor; + mode::AbstractBilevelSolverMode = SOS1Mode(), + add_bridges::Bool=true, + linearize_bilinear_upper_terms=false, + check_linearization_conditions=true + ) bm = BilevelModel() + bm.linearize_bilinear_upper_terms = linearize_bilinear_upper_terms + bm.check_linearization_conditions = check_linearization_conditions set_mode(bm, mode) JuMP.set_optimizer(bm, optimizer_constructor; add_bridges = add_bridges) return bm @@ -559,29 +567,21 @@ function JuMP.optimize!( end t0 = time() - - moi_upper = JuMP.index.(collect(values(model.upper_to_lower_link))) - moi_link = convert_indices(model.link) - moi_link2 = index2(model.upper_var_to_lower_ctr_link) + + lower_var_indices_of_upper_vars = JuMP.index.(collect(values(model.upper_to_lower_link))) + upper_to_lower_var_indices = convert_indices(model.link) + upper_var_lower_ctr = index2(model.upper_var_to_lower_ctr_link) reset!(mode) # cleaup cached data # build bound for FortunyAmatMcCarlMode build_bounds!(model, mode) - 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, - ) + @debug """starting build_bilevel at $(Dates.format(now(), "HH:MM:SS"))""" + single_blm, upper_to_sblm, lower_to_sblm, lower_primal_dual_map, lower_dual_to_sblm = + build_bilevel(upper, lower, upper_to_lower_var_indices, lower_var_indices_of_upper_vars, mode, upper_var_lower_ctr, + copy_names = model.copy_names, pass_start = model.pass_start, + linearize_bilinear_upper_terms = model.linearize_bilinear_upper_terms, + check_linearization_conditions = model.check_linearization_conditions) # pass additional info (hints - not actual problem data) # for lower level dual variables (start, upper hint, lower hint) diff --git a/src/moi.jl b/src/moi.jl index bf6a51b3..37f88197 100644 --- a/src/moi.jl +++ b/src/moi.jl @@ -199,26 +199,26 @@ end function build_full_map!( mode, - upper_idxmap, - lower_idxmap, - lower_dual_idxmap, - lower_primal_dual_map, + upper_to_m_idxmap, + lower_to_m_idxmap, + lower_dual_idxmap, + lower_primal_dual_map ) return nothing end function build_full_map!( mode::FortunyAmatMcCarlMode, - upper_idxmap, - lower_idxmap, - lower_dual_idxmap, - lower_primal_dual_map, + upper_to_m_idxmap, + lower_to_m_idxmap, + lower_dual_idxmap, + lower_primal_dual_map ) _build_bound_map!( mode.cache, - upper_idxmap, - lower_idxmap, - lower_dual_idxmap, - lower_primal_dual_map, + upper_to_m_idxmap, + lower_to_m_idxmap, + lower_dual_idxmap, + lower_primal_dual_map ) return nothing end @@ -240,17 +240,17 @@ function build_full_map!( end function _build_bound_map!( mode::ComplementBoundCache, - upper_idxmap, - lower_idxmap, - lower_dual_idxmap, - lower_primal_dual_map, + upper_to_m_idxmap, + lower_to_m_idxmap, + lower_dual_idxmap, + lower_primal_dual_map ) empty!(mode.map) - for (k, v) in mode.upper - mode.map[upper_idxmap[k]] = v + for (k,v) in mode.upper + mode.map[upper_to_m_idxmap[k]] = v end - for (k, v) in mode.lower - mode.map[lower_idxmap[k]] = v + for (k,v) in mode.lower + mode.map[lower_to_m_idxmap[k]] = v end for (k, v) in mode.ldual vec = lower_primal_dual_map.primal_con_dual_var[k]#[1] # TODO check this scalar @@ -344,15 +344,15 @@ function set_with_zero(set) return MOI.copy(set) end -function build_bilevel( - upper::MOI.ModelLike, + +function build_maps( + upper::MOI.ModelLike, lower::MOI.ModelLike, - upper_to_lower_link::Dict{VI,VI}, - upper_variables::Vector{VI}, + upper_to_lower_var_indices::Dict{VI,VI}, + lower_var_indices_of_upper_vars::Vector{VI}, mode, - upper_var_to_lower_ctr::Dict{VI,CI} = Dict{VI,CI}(); - copy_names::Bool = false, - pass_start::Bool = false, + upper_var_to_lower_ctr::Dict{VI,CI}, + copy_names::Bool, ) # Start with an empty problem @@ -365,11 +365,12 @@ function build_bilevel( #= Create Lower DUAL level model =# + # TODO if linearizing bilinear terms the lower dual model must be built using standard form! + # (for now requiring that lower model is in standard form) # dualize the second level - dual_problem = Dualization.dualize( - lower; - dual_names = DualNames("dual_", "dual_"), - variable_parameters = upper_variables, + dual_problem = Dualization.dualize(lower; + dual_names = Dualization.DualNames("dual_","dual_"), + variable_parameters = lower_var_indices_of_upper_vars, ignore_objective = ignore_dual_objective(mode), consider_constrained_variables = false, ) @@ -383,9 +384,9 @@ function build_bilevel( =# # keys are from src (upper), values are from dest (m) - upper_idxmap = MOIU.default_copy_to(m, upper) + upper_to_m_idxmap = MOIU.default_copy_to(m, upper) if copy_names - pass_names(m, upper, upper_idxmap) + pass_names(m, upper, upper_to_m_idxmap) end #= @@ -394,63 +395,39 @@ function build_bilevel( handle_lower_objective_sense(lower) - # cache and delete lower objective - if !ignore_dual_objective(mode) - # get primal obj - type_primal_obj = MOI.get(lower, MOI.ObjectiveFunctionType()) - @assert type_primal_obj !== nothing - lower_primal_obj = - MOI.get(lower, MOI.ObjectiveFunction{type_primal_obj}()) - # deepcopy and delete dual obj - # MOI.set(lower, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0)) - end - - # initialize map from (original) lower model to single model (m) - lower_idxmap = MOIU.IndexMap() - for (upper_key, lower_val) in upper_to_lower_link - lower_idxmap[lower_val] = upper_idxmap[upper_key] + # initialize map to lower level model + lower_to_m_idxmap = MOIU.IndexMap() + lower_to_upper_var_indices = Dict{MOI.VariableIndex, MOI.VariableIndex}() + for (upper_key, lower_val) in upper_to_lower_var_indices + lower_to_m_idxmap[lower_val] = upper_to_m_idxmap[upper_key] + lower_to_upper_var_indices[lower_val] = upper_key end # append the second level primal - append_to(m, lower, lower_idxmap; allow_single_bounds = true) + append_to(m, lower, lower_to_m_idxmap, allow_single_bounds = true) if copy_names - pass_names(m, lower, lower_idxmap) + pass_names(m, lower, lower_to_m_idxmap) end #= Pass Dual of Lower level model =# - # initialize map to lower level model - if !ignore_dual_objective(mode) - # get dual obj - tp_dual_obj = MOI.get(lower_dual, MOI.ObjectiveFunctionType()) - @assert tp_dual_obj !== nothing - lower_dual_obj = - MOI.get(lower_dual, MOI.ObjectiveFunction{tp_dual_obj}()) - # delete dual obj - # MOI.set(lower_dual, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0)) - end - - # initialize map from *dual lower* model to single model (m) + # initialize map from lower level dual variables to m indices lower_dual_idxmap = MOIU.IndexMap() - # 1) for QP's there are dual variable that are linked to: - # 1.1) primal variables - for (lower_primal_var_key, lower_dual_quad_slack_val) in - lower_primal_dual_map.primal_var_dual_quad_slack - lower_dual_idxmap[lower_dual_quad_slack_val] = - lower_idxmap[lower_primal_var_key] + # for lower level QP's there are lower dual variables that are tied to: + # lower primal variables + for (lower_primal_var_key, lower_dual_quad_slack_val) in lower_primal_dual_map.primal_var_dual_quad_slack + lower_dual_idxmap[lower_dual_quad_slack_val] = lower_to_m_idxmap[lower_primal_var_key] end - # 1.2) and to upper level variable which are lower level parameters - for (lower_primal_param_key, lower_dual_param_val) in - lower_primal_dual_map.primal_parameter - lower_dual_idxmap[lower_dual_param_val] = - lower_idxmap[lower_primal_param_key] + # and to upper level variables (which are lower level parameters) + for (lower_primal_param_key, lower_dual_param_val) in lower_primal_dual_map.primal_parameter + lower_dual_idxmap[lower_dual_param_val] = lower_to_m_idxmap[lower_primal_param_key] end - # 2) Dual variables might appear in the upper level + # lower level dual variable -> upper level variable for (upper_var, lower_con) in upper_var_to_lower_ctr var = lower_primal_dual_map.primal_con_dual_var[lower_con][1] # TODO check this scalar - lower_dual_idxmap[var] = upper_idxmap[upper_var] + lower_dual_idxmap[var] = upper_to_m_idxmap[upper_var] end # append the second level dual @@ -458,6 +435,64 @@ function build_bilevel( if copy_names pass_names(m, lower_dual, lower_dual_idxmap) end + return m, lower_dual, lower_primal_dual_map, upper_to_m_idxmap, lower_to_m_idxmap, lower_dual_idxmap +end + + +function build_bilevel( + upper::MOI.ModelLike, + lower::MOI.ModelLike, + upper_to_lower_var_indices::Dict{VI,VI}, + lower_var_indices_of_upper_vars::Vector{VI}, + mode, + upper_var_to_lower_ctr::Dict{VI,CI} = Dict{VI,CI}(); + copy_names::Bool = false, + pass_start::Bool = false, + linearize_bilinear_upper_terms::Bool = false, + check_linearization_conditions::Bool = true + ) + + m, lower_dual, lower_primal_dual_map, upper_to_m_idxmap, lower_to_m_idxmap, lower_dual_idxmap = build_maps( + upper, + lower, + upper_to_lower_var_indices, + lower_var_indices_of_upper_vars, + mode, + upper_var_to_lower_ctr, + copy_names + ) + + # cache and delete lower objective + if !BilevelJuMP.ignore_dual_objective(mode) + # get primal obj + type_primal_obj = MOI.get(lower, MOI.ObjectiveFunctionType()) + @assert type_primal_obj !== nothing + lower_primal_obj = MOI.get(lower, MOI.ObjectiveFunction{type_primal_obj}()) + # deepcopy and delete dual obj + # MOI.set(lower, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0)) + + # get dual obj + tp_dual_obj = MOI.get(lower_dual, MOI.ObjectiveFunctionType()) + @assert tp_dual_obj !== nothing + lower_dual_obj = + MOI.get(lower_dual, MOI.ObjectiveFunction{tp_dual_obj}()) + # delete dual obj + # MOI.set(lower_dual, MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(), MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0)) + end + + #= + Check if Upper bilinear terms of the form Upper primal * Lower dual exist, and can be linearized + =# + # TODO implement for bilinear products in constraints (only checking upper objective so far) + + if linearize_bilinear_upper_terms + main_linearization( + m, lower, upper, upper_var_to_lower_ctr, upper_to_lower_var_indices, + lower_var_indices_of_upper_vars, lower_to_m_idxmap, upper_to_m_idxmap, + lower_primal_dual_map, lower_dual_idxmap, + check_linearization_conditions + ) + end #= Additional Optimiality conditions (to complete the KKT) @@ -466,10 +501,10 @@ function build_bilevel( # build map bound map for FortunyAmatMcCarlMode build_full_map!( mode, - upper_idxmap, - lower_idxmap, - lower_dual_idxmap, - lower_primal_dual_map, + upper_to_m_idxmap, + lower_to_m_idxmap, + lower_dual_idxmap, + lower_primal_dual_map ) if ignore_dual_objective(mode) @@ -479,13 +514,13 @@ function build_bilevel( if !is_equality(comp.set_w_zero) accept_vector_set(mode, comp) add_complement( - mode, - m, + mode, + m, comp, - lower_idxmap, - lower_dual_idxmap, - copy_names, - pass_start, + lower_to_m_idxmap, + lower_dual_idxmap, + copy_names, + pass_start ) else # feasible equality constraints always satisfy complementarity @@ -493,21 +528,17 @@ function build_bilevel( end else # strong duality add_strong_duality( - mode, - m, - lower_primal_obj, - lower_dual_obj, - lower_idxmap, - lower_dual_idxmap, + mode, + m, + lower_primal_obj, + lower_dual_obj, + lower_to_m_idxmap, + lower_dual_idxmap ) end add_aggregate_constraints(m, mode, copy_names) - return m, - upper_idxmap, - lower_idxmap, - lower_primal_dual_map, - lower_dual_idxmap + return m, upper_to_m_idxmap, lower_to_m_idxmap, lower_primal_dual_map, lower_dual_idxmap end function add_strong_duality( diff --git a/src/moi_utilities.jl b/src/moi_utilities.jl index 32bccfe7..44865510 100644 --- a/src/moi_utilities.jl +++ b/src/moi_utilities.jl @@ -45,7 +45,7 @@ function append_to( # The `NLPBlock` assumes that the order of variables does not change (#849) if MOI.NLPBlock() in MOI.get(src, MOI.ListOfModelAttributesSet()) error("NLP models are not supported.") - # constraint_types = MOI.get(src, MOI.ListOfConstraints()) + # constraint_types = MOI.get(src, MOI.ListOfConstraintTypesPresent()) # single_variable_types = [S for (F, S) in constraint_types # if F == MOI.VariableIndex] # vector_of_variables_types = [S for (F, S) in constraint_types diff --git a/test/Project.toml b/test/Project.toml index 9ac931ee..43ec9dc9 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -10,5 +10,6 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" [compat] Ipopt = "=1.0.2" +JuMP = "=1.0" MibS_jll = "=1.1.3" QuadraticToBinary = "=0.4.0" \ No newline at end of file diff --git a/test/bilinear_linearization.jl b/test/bilinear_linearization.jl new file mode 100644 index 00000000..1b8b7e08 --- /dev/null +++ b/test/bilinear_linearization.jl @@ -0,0 +1,362 @@ +# using Logging; global_logger(ConsoleLogger(stderr, Logging.Debug)) +#= +Testing linearization of upper objective bilinear terms of lower duals and primals and supporting +functions. + +NOTE: test/jump.jl contains a test (jump_conejo2016_linearize) that confirms same results for a +bilinear and linearized problem. +=# + + +function test_recursive_col_search() + A = ones(2,2) # redundant rows + @test_throws BilevelJuMP.UnderDeterminedException BilevelJuMP.recursive_col_search(A, 1, 1, Int[], Int[]) + + A[1,2] = 0 # redundant cols + @test_throws BilevelJuMP.UnderDeterminedException BilevelJuMP.recursive_col_search(A, 1, 1, Int[], Int[2]) +end + + +function test_find_connected_rows_cols() + A = Matrix{Float64}([ + [1 0]; + [0 1] + ]) + rows, cols, redundant_vals = BilevelJuMP.find_connected_rows_cols(A, 1, 1; + skip_1st_col_check=false, finding_blocks=false + ) + @test isempty(rows) + @test cols == [1] + @test redundant_vals == false + + A[2,1] = 1 + rows, cols, redundant_vals = BilevelJuMP.find_connected_rows_cols(A, 1, 1; + skip_1st_col_check=false, finding_blocks=true + ) + @test rows == cols + + A = ones(2,2) + rows, cols, redundant_vals = BilevelJuMP.find_connected_rows_cols(A, 1, 1; + skip_1st_col_check=false, finding_blocks=false + ) + @test redundant_vals == true +end + + +function simple_linearization(optimizer, mode = BilevelJuMP.SOS1Mode()) + + cder = 1 + clmp = 1 + ci = 1 + d1 = 1 + d2 = 2 + MOI.empty!(optimizer) + # MILP Program (bilinear terms in upper and lower objectives get linearized) + model = BilevelModel( + ()->optimizer, + mode = mode, + linearize_bilinear_upper_terms=true + ) + @variables(Upper(model), begin + 10 >= x0 >= 0 + 10 >= xe >= 0 + end) + + @variables(Lower(model), begin + 10 >= ye >= 0 + 10 >= yi >= 0 + 10 >= yder >= 0 + end) + + @objective(Lower(model), Min, cder * yder + ci * yi - xe * ye) + @constraint(Lower(model), loadbal, yi - ye + yder == d1) + + @variable(Upper(model), lambda, DualOf(loadbal)) + + @objective(Upper(model), Min, clmp * x0 + lambda * ye) + @constraint(Upper(model), x0 + ye - yi - d2 == 0) + @constraint(Upper(model), [ye, yi] in MOI.SOS1([1.0, 2.0])) + + optimize!(model) + + @expression(model, LLcost, cder * yder + ci * yi - xe * ye) + @test objective_value(model) ≈ 2.0 + @test value(LLcost) ≈ 1.0 + @test value(yi) ≈ 0.0 + @test value(ye) ≈ 0.0 + @test value(yder) ≈ 1.0 + @test value(x0) ≈ 2.0 + + @test MOI.get(model.solver, MOI.ObjectiveFunctionType()) == MathOptInterface.ScalarAffineFunction{Float64} + +end + + +""" +make some dummy models that violate the linearization conditions when the set AB_N is non-empty, +i.e. when there are xm*yn in the LL objective + + +for local test: +```julia +using BilevelJuMP, JuMP, Gurobi, Test +optimizer = Gurobi.Optimizer() +mode = BilevelJuMP.SOS1Mode() + +``` +""" +function failing_conditions_non_empty_AB_N(mode = BilevelJuMP.SOS1Mode()) + + cder = 1 + clmp = 1 + ci = 1 + d1 = 1 + d2 = 2 + + model = BilevelModel() + model.linearize_bilinear_upper_terms = true + + @variables(Upper(model), begin + 10 >= x0 >= 0 + 10 >= xe >= 0 + 10 >= xbad1 >= 0 + end) + + @variables(Lower(model), begin + 10 >= ye >= 0 + 10 >= yi >= 0 + 10 >= yder >= 0 + 10 >= ybad3 >= 0 + 10 >= ybad4 >= 0 + 10 >= ybad5 >= 0 + end) + + # condition 2: no connected variables * upper variable in LL obj. + @objective(Lower(model), Min, cder * yder + ci * yi - xe * ye + 2*xe*yi) + + # @objective(Lower(model), Min, cder * yder + ci * yi - xe * ye) + + @constraint(Lower(model), loadbal, yi - ye + yder + ybad5 == d1) + # need minus sign on ybad5 to meet condition 5 + + # # adding xbad to connected constraint results in Condition 1 not being met + @constraint(Lower(model), badcon1, yi + xbad1 == d1) + + #= + adding ye (lower variable in upper bilinear obj. terms) to another constraint violates + Condition 4 + =# + @constraint(Lower(model), badcon4, ye + ybad4 == d1) + + # conditon 3: ybad3 is not connected to other yn in bilinear upper obj. terms + @constraint(Lower(model), badcon3, ybad3 + xbad1 == d1) + + @variable(Upper(model), lambda, DualOf(loadbal)) + @variable(Upper(model), dummylambda, DualOf(badcon3)) + + # having different coef.s on ratio of obj. coef.s to constraint coef.s on the upper obj. + # bilinear terms violates Condition 5 + @objective(Upper(model), Min, clmp * x0 + lambda * ye + lambda * ybad5 - dummylambda * ybad3) + @constraint(Upper(model), x0 + ye - yi - d2 == 0) + @constraint(Upper(model), [ye, yi] in MOI.SOS1([1.0, 2.0])) + + # code from JuMP.optimize! in src/jump.jl + upper = JuMP.backend(model.upper) + lower = JuMP.backend(model.lower) + lower_var_indices_of_upper_vars = JuMP.index.(collect(values(model.upper_to_lower_link))) + upper_to_lower_var_indices = BilevelJuMP.convert_indices(model.link) + upper_var_to_lower_ctr = BilevelJuMP.index2(model.upper_var_to_lower_ctr_link) + BilevelJuMP.build_bounds!(model, mode) + copy_names = false + + @test BilevelJuMP.is_model_in_standard_form(lower) == true + + # code from build_bilevel in src/moi.jl + m, lower_dual, lower_primal_dual_map, upper_to_m_idxmap, lower_to_m_idxmap, lower_dual_idxmap = + BilevelJuMP.build_maps( + upper, + lower, + upper_to_lower_var_indices, + lower_var_indices_of_upper_vars, + mode, + upper_var_to_lower_ctr, + copy_names + ) + + # code from main_linearization in src/bilinear_linearization.jl + A_N, bilinear_upper_dual_to_quad_term, bilinear_upper_dual_to_lower_primal, lower_primal_var_to_lower_con = + BilevelJuMP.check_upper_objective_for_bilinear_linearization( + upper, upper_to_lower_var_indices, upper_var_to_lower_ctr + ) + AB_N, B, lower_obj_terms, lower_obj_type_handled = + BilevelJuMP.get_lower_obj_coefs_of_upper_times_lower_primals( + lower, lower_var_indices_of_upper_vars, A_N, lower_to_m_idxmap + ) + U, V, w = BilevelJuMP.standard_form(lower, upper_var_indices=lower_var_indices_of_upper_vars) + J_U, N_U = BilevelJuMP.get_all_connected_rows_cols(upper_var_to_lower_ctr, bilinear_upper_dual_to_lower_primal, V, AB_N) + @test !isempty(AB_N) + num_blocks, rows, cols = BilevelJuMP.find_blocks(V, U) + + @test B[2, 5] == 2 + + # have to have function calls outside of @test for codecov to consider them tested? + con1 = BilevelJuMP.check_condition_1(J_U, U) + @test !(con1) + con2 = BilevelJuMP.check_condition_2prime(N_U, A_N, U, B) + @test !(con2) + con3 = BilevelJuMP.check_condition_3(A_N, V, lower_primal_var_to_lower_con) + @test !(con3) + con4 = BilevelJuMP.check_condition_4(A_N, V, upper_var_to_lower_ctr, bilinear_upper_dual_to_lower_primal) + @test !(con4) + con5 = BilevelJuMP.check_condition_5(A_N, V, upper_var_to_lower_ctr, + bilinear_upper_dual_to_lower_primal, bilinear_upper_dual_to_quad_term) + @test !(con5) + + BilevelJuMP.main_linearization( + m, + lower, + upper, + upper_var_to_lower_ctr, + upper_to_lower_var_indices, + lower_var_indices_of_upper_vars, + lower_to_m_idxmap, + upper_to_m_idxmap, + lower_primal_dual_map, + lower_dual_idxmap, + true + ) + check_false = BilevelJuMP.check_non_empty_AB_N_conditions(J_U, U, N_U, A_N, B, V, lower_primal_var_to_lower_con, + upper_var_to_lower_ctr, bilinear_upper_dual_to_quad_term, bilinear_upper_dual_to_lower_primal) + @test check_false == false + linearizations = + BilevelJuMP.linear_terms_for_non_empty_AB( + lower, + upper_var_to_lower_ctr, + bilinear_upper_dual_to_lower_primal, + V, + w, + A_N, + bilinear_upper_dual_to_quad_term, + lower_obj_terms, + lower_to_m_idxmap, + lower_primal_dual_map, + lower_dual_idxmap, + false, + num_blocks, + rows + ) +end + + +function failing_conditions_empty_AB_N(mode = BilevelJuMP.SOS1Mode()) + + cder = 1 + clmp = 1 + ci = 1 + d1 = 1 + d2 = 2 + + model = BilevelModel() + model.linearize_bilinear_upper_terms = true + + @variables(Upper(model), begin + 10 >= x0 >= 0 + 10 >= xe >= 0 + 10 >= xbad1 >= 0 + end) + + @variables(Lower(model), begin + 10 >= ye >= 0 + 10 >= yi >= 0 + 10 >= yder >= 0 + 10 >= ybad4 >= 0 + end) + + @objective(Lower(model), Min, cder * yder + ci * yi) + @constraint(Lower(model), loadbal, yi - ye + yder + xbad1 == d1) + @constraint(Lower(model), badcon1, ye + ybad4 + xbad1 == d1) + + @variable(Upper(model), lambda, DualOf(loadbal)) + + @objective(Upper(model), Min, clmp * x0 + lambda * ye) + @constraint(Upper(model), x0 + ye - yi - d2 == 0) + @constraint(Upper(model), [ye, yi] in MOI.SOS1([1.0, 2.0])) + + + # code from JuMP.optimize! in src/jump.jl + upper = JuMP.backend(model.upper) + lower = JuMP.backend(model.lower) + lower_var_indices_of_upper_vars = JuMP.index.(collect(values(model.upper_to_lower_link))) + upper_to_lower_var_indices = BilevelJuMP.convert_indices(model.link) + upper_var_to_lower_ctr = BilevelJuMP.index2(model.upper_var_to_lower_ctr_link) + BilevelJuMP.build_bounds!(model, mode) + copy_names = false + + @test BilevelJuMP.is_model_in_standard_form(lower) == true + + # code from build_bilevel in src/moi.jl + m, lower_dual, lower_primal_dual_map, upper_to_m_idxmap, lower_to_m_idxmap, lower_dual_idxmap = + BilevelJuMP.build_maps( + upper, + lower, + upper_to_lower_var_indices, + lower_var_indices_of_upper_vars, + mode, + upper_var_to_lower_ctr, + copy_names + ) + + # code from main_linearization in src/bilinear_linearization.jl + A_N, bilinear_upper_dual_to_quad_term, bilinear_upper_dual_to_lower_primal, lower_primal_var_to_lower_con = + BilevelJuMP.check_upper_objective_for_bilinear_linearization( + upper, upper_to_lower_var_indices, upper_var_to_lower_ctr + ) + AB_N, B, lower_obj_terms, lower_obj_type_handled = + BilevelJuMP.get_lower_obj_coefs_of_upper_times_lower_primals( + lower, lower_var_indices_of_upper_vars, A_N, lower_to_m_idxmap + ) + U, V, w = BilevelJuMP.standard_form(lower, upper_var_indices=lower_var_indices_of_upper_vars) + J_U, N_U = BilevelJuMP.get_all_connected_rows_cols(upper_var_to_lower_ctr, bilinear_upper_dual_to_lower_primal, V, AB_N) + num_blocks, rows, cols = BilevelJuMP.find_blocks(V, U) + + @test isempty(AB_N) + @test num_blocks == 1 + @test rows == [[1, 2]] + @test cols == [[3, 4, 5, 6, 7]] + @test !BilevelJuMP.check_empty_AB_N_conditions(J_U, U, N_U, B) + @test_throws BilevelJuMP.UnderDeterminedException BilevelJuMP.recursive_col_search(U+V, 1, 3, Int[], Int[]) + + # to reach more code + @constraint((Lower(model)), ybad4 <= ye) + @constraint((Lower(model)), yder >= ye) + @constraint((Lower(model)), ybad4 <= 2) + @constraint((Lower(model)), yder >= 1) + U, V, w = BilevelJuMP.standard_form(lower, upper_var_indices=lower_var_indices_of_upper_vars) + +end + +function test_get_coef_matrix_and_rhs_vec() + model = BilevelModel() + + @variable(Upper(model), x) + @variable(Lower(model), y) + + @constraint(Lower(model), 2x + 3y <= 4) + @constraint(Lower(model), 5x + 6y >= 7) + lower = JuMP.backend(model.lower) + con_indices = MOI.get(lower, MOI.ListOfConstraintIndices{ + MOI.ScalarAffineFunction{Float64}, + MOI.LessThan{Float64} + }()) + A, b = BilevelJuMP.get_coef_matrix_and_rhs_vec(lower, con_indices) + @test A == [2 3.0] + @test b == [4.0] + + con_indices = MOI.get(lower, MOI.ListOfConstraintIndices{ + MOI.ScalarAffineFunction{Float64}, + MOI.GreaterThan{Float64} + }()) + A, b = BilevelJuMP.get_coef_matrix_and_rhs_vec(lower, con_indices) + @test A == [5 6.0] + @test b == [7.0] +end diff --git a/test/jump.jl b/test/jump.jl index 0e6afddd..18ec63c9 100644 --- a/test/jump.jl +++ b/test/jump.jl @@ -2336,6 +2336,45 @@ function jump_conejo2016(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Conf @test value(lambda) ≈ 15 atol=1e-3 rtol=1e-2 end +""" + jump_conejo2016_linearize(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config()) + +Test linearization of bilinear terms in upper objective of the form lower dual * lower primal. +NOTE: lower model must be in standard form +""" +function jump_conejo2016_linearize(optimizer, mode = BilevelJuMP.SOS1Mode(), config = Config()) + + atol = config.atol + start = config.start_value + + MOI.empty!(optimizer) + model = BilevelModel(()->optimizer, mode = mode, linearize_bilinear_upper_terms=true) + BilevelJuMP.set_copy_names(model) + + @variable(Upper(model), 0 <= x <= 250, start = 50) + @variable(Lower(model), 0 <= y[i=1:3] <= [300, 150, 100][i], start = [50, 150, 0][i]) + @variable(Lower(model), 300 >= s >= 0) + + @objective(Lower(model), Min, 10y[1] + 12y[2] + 15y[3]) + + @constraint(Lower(model), b, y[1] + y[2] + y[3] == 200) + @constraint(Lower(model), b2, y[1] - x + s == 0) + @variable(Upper(model), 0 <= lambda <= 20, DualOf(b), start = 15) + + @objective(Upper(model), Min, 40_000x + 8760*(10y[1]-lambda*y[1])) + + optimize!(model) + + primal_status(model) + termination_status(model) + + @test objective_value(model) ≈ -190_000 atol=1e-1 rtol=1e-2 + @test value(x) ≈ 50 atol=1e-3 rtol=1e-2 + @test value.(y) ≈ [50, 150, 0] atol=1e-3 rtol=1e-2 + @test value(lambda) ≈ 15 atol=1e-3 rtol=1e-2 + @test MOI.get(model.solver, MOI.ObjectiveFunctionType()) == MathOptInterface.ScalarAffineFunction{Float64} +end + #= Bruno Fanzeres PhD thesis Robust Strategic Bidding in Auction-Based Markets. =# diff --git a/test/runtests.jl b/test/runtests.jl index 2dd4e355..dbecd1ce 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -78,6 +78,7 @@ include("moi.jl") include("jump.jl") include("jump_unit.jl") include("jump_nlp.jl") +include("bilinear_linearization.jl") @testset "BilevelJuMP tests" begin @testset "MibS" begin @@ -280,6 +281,27 @@ include("jump_nlp.jl") jump_01_sum_agg(solver.opt) end end + for solver in solvers_nlp_sd + jump_conejo2016(solver.opt, solver.mode) + jump_conejo2016_linearize(solver.opt, solver.mode, config) + # jump_fanzeres2017(solver.opt, solver.mode) + # jump_eq_price(solver.opt, solver.mode) + end + for solver in solvers_sos_quad_bin + jump_conejo2016(solver.opt, solver.mode, config, bounds = true) # fail travis on cbc + jump_conejo2016_linearize(solver.opt, solver.mode, config) + # jump_fanzeres2017(solver.opt, solver.mode) + # jump_eq_price(solver.opt, solver.mode) # fail travis on cbc + end + for solver in solvers_fa_quad_bin + jump_conejo2016(solver.opt, solver.mode, config, bounds = true) + jump_conejo2016_linearize(solver.opt, solver.mode, config) + # jump_fanzeres2017(solver.opt, solver.mode) + # jump_eq_price(solver.opt, solver.mode) # fail with Gurobi b/c `Variable MathOptInterface.VariableIndex(17) has no upper bound`` + end + for solver in solvers_sos + jump_conejo2016_linearize(solver.opt, solver.mode, config) + end @testset "Princeton Handbook Quadratic" begin for is_min in [true, false] @@ -401,4 +423,15 @@ include("jump_nlp.jl") @time jump_conic04(solver.opt, solver.mode, config, bounds = true) end end + + @testset "bilinear_linearization" begin + test_recursive_col_search() + test_find_connected_rows_cols() + test_get_coef_matrix_and_rhs_vec() + failing_conditions_non_empty_AB_N() + failing_conditions_empty_AB_N() + for solver in solvers_sos + simple_linearization(solver.opt, solver.mode) + end + end end