# Advanced-Use Example 2

This example is also provided here as a Jupyter Notebook.

### An $\alpha$BB Example for a QCQP.

In this example, we will demonstrate the use of a user-defined lower-bounding problem that uses $\alpha$BB convex relaxations. In this example, we wish to solve the *nonconvex* quadratically constrained quadratic program (QCQP):

\[\begin{aligned} & \min_{\mathbf x \in \mathbb{IR}^{2}} \frac{1}{2} \mathbf x^{\rm T} \mathbf Q_{f} \mathbf x + \mathbf c_{f}^{\rm T} \mathbf x \\ {\rm s.t.} \; \; & g_{1}(\mathbf x) = \frac{1}{2} \mathbf x^{\rm T} \mathbf Q_{g_{1}}\mathbf x + \mathbf c_{g_{1}}^{\rm T} \mathbf x \leq 0 \\ & g_{2}(\mathbf x) = \frac{1}{2} \mathbf x^{\rm T} \mathbf Q_{g_{2}} \mathbf x + \mathbf c_{g_{2}}^{\rm T} \mathbf x \leq 0 \\ \end{aligned}\]

with $\mathbf Q_i\in\mathbb R^{2\times 2}$ not positive semidefinite for any $i$.

`using JuMP, EAGO, Ipopt, LinearAlgebra`

## Custom Function Definitions

For convenience, we'll define the following function that returns all the problem data $\mathbf Q_i$ and $\mathbf c_i$.

```
function QCQP_setup()
Qf = [3.0, 3/2; 3/2, -5.0]
cf = [3.0; 2.0]
Qg1 = [-2.0, 5.0; 5.0, -2.0]
cg1 = [1.0; 3.0]
Qg2 = [-6.0, 3.0; 3.0, 2.0]
cg2 = [2.0; 1.0]
return Qf, cf, Qg1, cg1, Qg2, cg2
end
```

The next function we'll define will take as input data for a particular quadratic function and the interval bounds on the decision variables, and construct an $\alpha$BB convex relaxation of that function. Since we're solving a QCQP, we'll use the `eigvals`

function to directly compute the eigenvalues of the input $\mathbf Q_i$ matrix.

```
function αBB_relax(Q::Matrix{T}, c::Vector{T}, xL::Vector{T}, xU::Vector{T}, x::Real...) where {T<:Float64}
α = max(0.0, -minimum(eigvals(Q))/2)
y = [x[1]; x[2]]
cv = 1/2*y'*Q*y + c'*y + α*(xL - y)'*(xU - y)
return cv
end
```

## Redefine the Lower-Bounding Problem

The following code first defines our EAGO extension (custom version) struct and then it redefines the lower-bounding problem as our own version. That is, when we call this customized version of EAGO to solve the problem, it'll deploy this version of the lower-bounding problem instead of the default version.

```
import EAGO: Optimizer, GlobalOptimizer
struct αBB_Convex <: EAGO.ExtensionType end
import EAGO: lower_problem!
function EAGO.lower_problem!(t::αBB_Convex, opt::GlobalOptimizer)
# Get active node
n = opt._current_node
# Get bounds on active node for calculating relaxations
xL = n.lower_variable_bounds[1:2]
xU = n.upper_variable_bounds[1:2]
# Get the problem data
Qf, cf, Qg1, cg1, Qg2, cg2 = QCQP_setup()
# Define the JuMP model and declare the solver
mL = JuMP.Model(JuMP.optimizer_with_attributes(Ipopt.Optimizer,
"tol" => 1.0e-6,
"print_level" => 0))
@variable(mL, xL[i] <= x[i=1:2] <= xU[i])
# Define the function closures for the user-defined relaxations
fcv(x...) = αBB_relax(Qf, cf, xL, xU, x...)
g1cv(x...) = αBB_relax(Qg1, cg1, xL, xU, x...)
g2cv(x...) = αBB_relax(Qg2, cg2, xL, xU, x...)
# Register the user-defined functions
# Note: If the gradients and Hessians are directly available, they could
# be passed as arguments to the register function to speed things up.
JuMP.register(mL, :fcv, 2, fcv, autodiff=true)
JuMP.register(mL, :g1cv, 2, g1cv, autodiff=true)
JuMP.register(mL, :g2cv, 2, g2cv, autodiff=true)
# Declare the objective function and constraints
@NLobjective(mL, Min, fcv(x[1], x[2]))
@NLconstraint(mL, g1cv(x[1], x[2]) <= 0.0)
@NLconstraint(mL, g2cv(x[1], x[2]) <= 0.0)
# Solve the relaxed problem
JuMP.optimize!(mL)
# Get primal status, termination status, determine if a global solution was obtained
tstatus = MOI.get(mL, MOI.TerminationStatus())
pstatus = MOI.get(mL, MOI.PrimalStatus())
solution = JuMP.value.(x)
# Interpret status codes for branch-and-bound
if EAGO.local_problem_status(tstatus, pstatus) == EAGO.LRS_FEASIBLE
opt._lower_objective_value = JuMP.objective_value(mL)
opt._lower_solution[1:length(solution)] = solution
opt._lower_feasibility = true
opt._cut_add_flag = false
else
opt._lower_feasibility = false
opt._lower_objective_value = -Inf
opt._cut_add_flag = false
end
return
end
```

By default, EAGO solves the epigraph reformulation of your original problem, which increases the original problem dimensionality by +1 with the introduction of an auxiliary variable. When defining custom routines (such as the lower-bounding problem here) that are intended to work nicely with default EAGO routines (such as preprocessing), the user must account for the *new* dimensionality of the problem. In the code above, we wish to access the information of the specific B&B node and define an optimization problem based on that information. However, in this example, the node has information for 3 variables (the original 2 plus 1 for the auxiliary variable appended to the original variable vector) as $(x_{1}, x_{2}, \eta)$. The lower-bounding problem was defined to optimize the relaxed problem with respect to the original 2 decision variables. When storing the results of this subproblem to the current B&B node, it is important to take care to store the information at the appropriate indices and not inadvertently redefine the problem dimensionality (i.e., by simply storing the optimization solution as the `lower_solution`

of the current node). For problems that are defined to only branch on a subset of the original variables, the optimizer has a member `_sol_to_branch_map`

that carries the mapping between the indices of the original variables to those of the variables being branched on. Visit our quasiconvex example to see how this is done.

## (Optional) Turn Off Processing Routines

Turn off preprocessing routines if you don't want to use them as defined in EAGO.

```
import EAGO: preprocess!
function EAGO.preprocess!(t::αBB_Convex, x::GlobalOptimizer)
x._preprocess_feasibility = true
return
end
```

Turn off postprocessing routines if you don't want to use them as defined in EAGO.

```
import EAGO: postprocess!
function EAGO.postprocess!(t::αBB_Convex, x::GlobalOptimizer)
x._postprocess_feasibility = true
return
end
```

## Construct the JuMP Model and Optimize

Now, we'll tell EAGO to use our custom/extended solver, set up the main JuMP model, and solve it with our custom solver.

```
factory = () -> EAGO.Optimizer(SubSolvers(; t = αBB_Convex() ))
m = JuMP.Model(optimizer_with_attributes(factory,
"relative_tolerance" => 1e-3,
"verbosity" => 1,
"output_iterations" => 1,
"branch_variable" => Bool[true; true],
))
Qf, cf, Qg1, cg1, Qg2, cg2 = QCQP_setup() # Get QCQP data
xL = [-3.0; -5.0] # Lower bounds on x
xU = [1.0; 2.0] # Upper bounds on x
@variable(m, xL[i] <= x[i=1:2] <= xU[i])
# Define objective and constraints
@objective(m, Min, 1/2*x'*Qf*x + cf'*x)
@constraint(m, 1/2*x'*Qg1*x + cg1'*x <= 0.0)
@constraint(m, 1/2*x'*Qg2*x + cg2'*x <= 0.0)
# Solve the problem
@time optimize!(m)
```

## Retrieve Results

We then recover the objective value, the solution values, and termination status codes using standard JuMP syntax.

```
println("x1* = ", JuMP.value(x[1]), " x2* = ",
JuMP.value(x[2])," f* = ",JuMP.objective_value(m))
TermStatus = JuMP.termination_status(m)
PrimStatus = JuMP.primal_status(m)
println("Algorithm terminated with a status of $TermStatus and result code of $PrimStatus")
```