Can Julia solve a nonlinear system of equation which nests a linear programming problem?

General Usage
,
1

I am trying to find a “real-valued” solution of a fixed point problem which nests a constrained linear programming problem by using NLsolve.jl and Gurobi.jl. Concretely, the outer problem is a nonlinear system of equations to find a fixed point given the inner solution by NLsolve.jl, whereas the inner problem solves a linear programming given the outer solution by Gurobi.jl. Technically, it is a bilevel problem. Here is the example code.

using Gurobi, NLsolve
N = 2
function inner_linear_programming(outer_sol)
  utility = zeros(N,N)
  # coefficients of utility are arbitrary but continuous in outer_sol 
  utility[1,1] = 2 * outer_sol[1,1];
  utility[1,2] = 2 * outer_sol[2,1]-1;
  utility[2,1] = 2 * outer_sol[1,2];
  utility[2,2] = 2 * outer_sol[1,1]-1;
  # vectorize utility term for Gurobi format
  f_utility = zeros(N*N)
  for i in 1:N
    for j in 1:(N)
        f_utility[(i-1)*N+j] = utility[i,j];
    end
  end
  # feasibility constraints mu:
  A_feas_mu_const = [0.0 1.0 1.0 0.0]
  b_feas_mu_const = 1.0
  # feasibility constraints rho(measure const):
  A_feas_rho_const = [1.0 1.0 0.0 0.0;
                     0.0 0.0 1.0 1.0]
  b_feas_rho_const = ones(N)
  A_feas_mu_rho = vcat(A_feas_mu_const, A_feas_rho_const)
  b_feas_mu_rho = vcat(b_feas_mu_const, b_feas_rho_const)
  # solve inner problem by plain-Gurobi
  env = Gurobi.Env()
  setparam!(env, "Presolve", 0)
  # For Gurobi format
  @time model = gurobi_model(env;
    name = "lp_01", # minimizing problem
    f = -f_utility, # vectorized coefficients in LP objective function
    Aeq = A_feas_mu_rho, # LHS of equality constraints matrix
    beq = b_feas_mu_rho, # RHS of equality constraints (vector)
    lb = zeros(N*N), # inequality consraints vectorized
    ub = ones(N*N) # inequality consraints vectorized
    )
  @time Gurobi.optimize(model)
  sol = get_solution(model)
  opt_x = convert(Array, transpose(reshape(sol, N, N)))
  return opt_x
end

function f!(F, x) # outer fixed point problem
    # solve 0 = inner_linear_programming(x) - x
    F = inner_linear_programming(x) - x
    @show inner_linear_programming(x)
    @show x
end
res=NLsolve.nlsolve(f!, [0.0 1.0; 1.0 0.0])
res.zero # solution found by NLsolve (not converged)
@show inner_linear_programming(res.zero) - res.zero
# it should be zeros if it is evaluated at fixed point

The problem has a real-valued fixed point solution theoretically by Kakutani’s fixed point theorem, allowing arbitrary (continuous) specification of coefficients of objective function. However, NLsolve.jl (root-finding)solver cannot find the “real-valued” solution and the inner linear programming problem oscillates between integer solutions. If the true solution is integer, NLsolve.jl can find it. Other minimization solvers such as NLopt.jl and Optim.jl work. However, these solve the inner and outer problem iteratively so oscillate between integer solutions similarly. This type of linear programming is likely to have an integer solution if solved iteratively, so I ideally need to solve whole problem simultaneously as a root-finding problem.

My question is whether Julia can solve nonlinear system of equations nesting linear programming solver such as Gurobi “simultaneously” to find a “real-valued” solution.

Thanks in advance.

2

These problems can be very tricky. I would try to reformulate, eg find a contraction for the fixed point part and iterate (at least initially, until I am close to the root), or rewrite the whole thing using a primal/dual approach.

If you otherwise have full control of the model, it may be cleaner to just switch to a utility/production function which gives you a clean interior solution you can characterize with FOC.

1 Like
3

Thanks for considering this tricky problem.
The problem does not have any contraction property even though the existence of solution is verified theoretically. This is a serious problem. Also, the inner constrained linear programming problem cannot be characterized as some convenient FOC of outer solutions.
So, it needs “joint” optimization procedure explained above.

4

I think the first step here is not a julia one. You need to figure out the best formulation of this in terms of optimization lingo/etc. and then you can search to see if a julia package exists to implement that formulation.

One trick at the top level is that large systems can be cast as nonlinear least squares, making it a nested optimization problem rather than a system with optimization inside. Otherwise, while these nested problems often have an MPEC form, I think that yours is kind of the other way around (i.e. usually the MPEC is an optimization with the nested problem written as FOCs complementarity conditions). I believe the optimization crew might call the generalization of this stuff EMP? Extended Mathematical Programming - Wikipedia

If you can figure out how this stuff might fit into those frameworks, it is possible that Knitro and others could implement it. stackexchange is often useful for getting help on optimization formulations.

2 Likes
5

Thanks for considering this tricky problem.

I think that yours is kind of the other way around (i.e. usually the MPEC is an optimization with the nested problem written as FOCs complementarity conditions). I believe the optimization crew might call the generalization of this stuff EMP?

Yes, this is the MPEC. Actually, I have tried bilevel reformulations and used GAMS and AMPL/KNITRO with/without calling EMP options so far. Then, the model became a nonlinear optimization problem with many constraints and locally found a unique real-valued solution in some cases from any starting value for both solvers. So, the reformulation was correct. However, I found that the outer objective value evaluated at their fixed point solution was not zero so it is not a fixed point, particularly the inner linear programming problem at their fixed point solution generated a different integer solution not real-valued one. It implies that the procedure is not adequate for finding a fixed point solution in my case.

This is the reason why I want to solve this tricky problem ``jointly" in Julia, which is easy to implement flexible combinations of many useful commercial solvers.