Nested optimization problems
This tutorial was generated using Literate.jl. Download the source as a .jl file.
This tutorial shows how to solve a nested (bilevel) optimization problem in which the upper-level problem uses the optimal value of a lower-level subproblem as an objective term. It also demonstrates how to use a cache to avoid redundant subproblem solves.
Learning intentions:
- Decompose a bilevel program and expose the lower-level value function to the upper level as a user-defined operator via
@operator - Cache the lower-level solution so that the function and gradient at the same point all share a single subproblem solve
- Use DiffOpt.jl to compute a Hessian callback by differentiating the lower-level solution with respect to the upper-level solution.
For a simpler example of writing a user-defined operator, see the User-defined Hessians tutorial.
The JuMP extension BilevelJuMP.jl can also be used to model and solve bilevel optimization problems.
Required packages
This tutorial uses the following packages:
using JuMPimport DiffOptimport IpoptFormulation
In the rest of this tutorial, our goal is to solve the bilevel optimization problem:
\[\begin{array}{r l} \min\limits_{x,z} & x_1^2 + x_2^2 + z \\ s.t. & \begin{array}{r l} z = \max\limits_{y} & x_1^2 y_1 + x_2^2 y_2 - x_1 y_1^4 - 2 x_2 y_2^4 \\ s.t. & (y_1 - 10)^2 + (y_2 - 10)^2 \le 25 \end{array} \\ & x \ge 0. \end{array}\]
This bilevel optimization problem is composed of two nested optimization problems. An upper level, involving variables $x$, and a lower level, involving variables $y$. From the perspective of the lower-level problem, the values of $x$ are fixed parameters, and so the model optimizes $y$ given those fixed parameters. Simultaneously, the upper-level problem optimizes $x$ and $z$ given the response of $y$.
Decomposition
There are a few ways to solve this problem, but we are going to use a nonlinear decomposition method. The first step is to write a function to compute the lower-level problem:
\[\begin{array}{r l} V(x_1, x_2) = \max\limits_{y} & x_1^2 y_1 + x_2^2 y_2 - x_1 y_1^4 - 2 x_2 y_2^4 \\ s.t. & (y_1 - 10)^2 + (y_2 - 10)^2 \le 25 \end{array}\]
function solve_lower_level(x...) model = Model(Ipopt.Optimizer) set_silent(model) @variable(model, y[1:2]) @objective( model, Max, x[1]^2 * y[1] + x[2]^2 * y[2] - x[1] * y[1]^4 - 2 * x[2] * y[2]^4, ) @constraint(model, (y[1] - 10)^2 + (y[2] - 10)^2 <= 25) optimize!(model) assert_is_solved_and_feasible(model) return objective_value(model), value.(y)endsolve_lower_level (generic function with 1 method)The function above takes a value of $x$ and returns the optimal lower-level objective-value and the optimal response $y$. The reason why we need both the objective and the optimal $y$ will be made clear shortly, but for now let us define:
function V(x...) f, _ = solve_lower_level(x...) return fendV (generic function with 1 method)Then, we can substitute $V$ into our full problem to create:
\[\begin{array}{r l} \min\limits_{x} & x_1^2 + x_2^2 + V(x_1, x_2) \\ s.t. & x \ge 0. \end{array}\]
This looks like a nonlinear optimization problem with a user-defined operator $V$! However, because $V$ solves an optimization problem internally, we can't use automatic differentiation to compute the first and second derivatives. Instead, we can use JuMP's ability to pass callback functions for the gradient (and later, Hessian) instead.
First up, we need to define the gradient of $V$ with respect to $x$. In general, this may be difficult to compute, but because $x$ appears only in the objective, we can just differentiate the objective function with respect to $x$, giving:
function ∇V(g::AbstractVector, x...) _, y = solve_lower_level(x...) g[1] = 2 * x[1] * y[1] - y[1]^4 g[2] = 2 * x[2] * y[2] - 2 * y[2]^4 returnend∇V (generic function with 1 method)We now have enough to define our bilevel optimization problem:
model = Model(Ipopt.Optimizer)@variable(model, x[1:2] >= 0)@operator(model, op_V, 2, V, ∇V)@objective(model, Min, x[1]^2 + x[2]^2 + op_V(x[1], x[2]))optimize!(model)assert_is_solved_and_feasible(model)solution_summary(model)solution_summary(; result = 1, verbose = false)
├ solver_name : Ipopt
├ Termination
│ ├ termination_status : LOCALLY_SOLVED
│ ├ result_count : 1
│ └ raw_status : Solve_Succeeded
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ └ objective_value : -4.18983e+05
└ Work counters
├ solve_time (sec) : 4.91825e+00
└ barrier_iterations : 8The optimal objective value is:
objective_value(model)-418983.48680640804and the optimal upper-level decision variables $x$ are:
value.(x)2-element Vector{Float64}:
154.97862343623828
180.00961424557616To compute the optimal lower-level decision variables, we need to call solve_lower_level with the optimal upper-level decision variables:
_, y = solve_lower_level(value.(x)...)y2-element Vector{Float64}:
7.072593960662295
5.946569893186169Improving performance
Our solution approach works, but it has a performance problem: every time we need to compute the value or gradient of $V$, we have to re-solve the lower-level optimization problem. This is wasteful, because we will often call the function and gradient at the same point, and so solving the problem twice with the same input repeats work unnecessarily.
We can work around this by using a cache:
mutable struct Cache x::Any f::Float64 y::Vector{Float64}endwith a function to update the cache if needed:
function _update_if_needed(cache::Cache, x...) if cache.x !== x cache.f, cache.y = solve_lower_level(x...) cache.x = x end returnend_update_if_needed (generic function with 1 method)Then, we define cached versions of our two functions which call _update_if_needed and return values from the cache.
function cached_f(cache::Cache, x...) _update_if_needed(cache, x...) return cache.fendfunction cached_∇f(cache::Cache, g::AbstractVector, x...) _update_if_needed(cache, x...) g[1] = 2 * x[1] * cache.y[1] - cache.y[1]^4 g[2] = 2 * x[2] * cache.y[2] - 2 * cache.y[2]^4 returnendcached_∇f (generic function with 1 method)Now we're ready to setup and solve the upper level optimization problem:
model = Model(Ipopt.Optimizer)@variable(model, x[1:2] >= 0)cache = Cache(Float64[], NaN, Float64[])@operator( model, op_cached_f, 2, (x...) -> cached_f(cache, x...), (g, x...) -> cached_∇f(cache, g, x...),)@objective(model, Min, x[1]^2 + x[2]^2 + op_cached_f(x[1], x[2]))optimize!(model)assert_is_solved_and_feasible(model)solution_summary(model)solution_summary(; result = 1, verbose = false)
├ solver_name : Ipopt
├ Termination
│ ├ termination_status : LOCALLY_SOLVED
│ ├ result_count : 1
│ └ raw_status : Solve_Succeeded
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ └ objective_value : -4.18983e+05
└ Work counters
├ solve_time (sec) : 1.05917e-01
└ barrier_iterations : 8and we can check we get the same objective value:
objective_value(model)-418983.48680640804and upper-level decision variable $x$:
value.(x)2-element Vector{Float64}:
154.97862343623828
180.00961424557616The problem takes 8 iterations to converge:
barrier_iterations(model)8Computing the Hessian
To improve performance we can pass a function that computes the Hessian of V with respect to the inputs x. Computing this Hessian requires some calculus.
Let $f(x, y)$ be the objective function, then:
\[\nabla^2_{xx} V(x) = \nabla^2_{xx}f(x, y^*) + \nabla^2_{xy}f(x, y^*) \cdot \frac{d y^*}{d x}\]
It is easy to compute $\nabla^2_{xx} f$ and $\nabla^2_{xy} f$ by hand. Computing $\frac{d y^*}{d x}$ (the derivative of the optimal solution $y^*$ with respect to the input $x$) is trickier. However, JuMP has a package, DiffOpt.jl which can do this for us.
function solve_lower_level_with_sensitivity(x...) model = DiffOpt.nonlinear_diff_model(Ipopt.Optimizer) set_silent(model) # Instead of directly using `x`, introduce a parameter `p` instead. @variable(model, p[i in 1:2] in Parameter(x[i])) @variable(model, y[1:2]) @objective( model, Max, p[1]^2 * y[1] + p[2]^2 * y[2] - p[1] * y[1]^4 - 2 * p[2] * y[2]^4, ) @constraint(model, (y[1] - 10)^2 + (y[2] - 10)^2 <= 25) optimize!(model) assert_is_solved_and_feasible(model) y_star = value.(y) dy_dx = zeros(2, 2) for j in 1:2 DiffOpt.set_forward_parameter.(model, p, 0.0) DiffOpt.set_forward_parameter(model, p[j], 1.0) DiffOpt.forward_differentiate!(model) dy_dx[:, j] .= DiffOpt.get_forward_variable.(model, y) end return y_star, dy_dxendfunction ∇²V(H::AbstractMatrix, x...) y, dy_dx = solve_lower_level_with_sensitivity(x...) ∇²V_xx = [2 * y[1] 0; 0 2 * y[2]] ∇²V_xy = [(2 * x[1] - 4 * y[1]^3) 0; 0 (2 * x[2] - 8 * y[2]^3)] ret = ∇²V_xx + ∇²V_xy * dy_dx for j in 1:size(H, 2), i in j:size(H, 1) H[i, j] = ret[i, j] end returnendmodel = Model(Ipopt.Optimizer)@variable(model, x[1:2] >= 0)@operator(model, op_V, 2, V, ∇V, ∇²V)@objective(model, Min, x[1]^2 + x[2]^2 + op_V(x[1], x[2]))optimize!(model)assert_is_solved_and_feasible(model)solution_summary(model)solution_summary(; result = 1, verbose = false)
├ solver_name : Ipopt
├ Termination
│ ├ termination_status : LOCALLY_SOLVED
│ ├ result_count : 1
│ └ raw_status : Solve_Succeeded
├ Solution (result = 1)
│ ├ primal_status : FEASIBLE_POINT
│ ├ dual_status : FEASIBLE_POINT
│ └ objective_value : -4.18983e+05
└ Work counters
├ solve_time (sec) : 1.45804e+01
└ barrier_iterations : 7The optimal objective value is:
objective_value(model)-418983.48680640955and the optimal upper-level decision variables $x$ are:
value.(x)2-element Vector{Float64}:
154.978623436436
180.00961424578398With the Hessian, the problem now takes 7 iterations to converge:
barrier_iterations(model)7