GitHub

Quadratic Program (QP) Example in Julia

We will show how to solve the Basic QP example problem both natively in Clarabel.jl and also by solving with Clarabel.jl within either JuMP or Convex.jl.

Clarabel.jl native interface

using Clarabel, SparseArrays

P = sparse([3. 0.;0. 2.].*2)

q = [-1., -4.]

A = sparse([1. -2.;    #<-- LHS of equality constraint
            1.  0.;    #<-- LHS of inequality constraint (upper bound)
            0.  1.;    #<-- LHS of inequality constraint (upper bound)
           -1.  0.;    #<-- LHS of inequality constraint (lower bound)
            0. -1.;    #<-- LHS of inequality constraint (lower bound)
            ])

b = [0.;        #<-- RHS of equality constraint
     ones(4)]   #<-- RHS of inequality constraints

cones =
    [Clarabel.ZeroConeT(1),           #<--- for the equality constraint
     Clarabel.NonnegativeConeT(4)]    #<--- for the inequality constraints

settings = Clarabel.Settings()

solver   = Clarabel.Solver()

Clarabel.setup!(solver, P, q, A, b, cones, settings)

result = Clarabel.solve!(solver)
-------------------------------------------------------------
           Clarabel.jl v0.8.1  -  Clever Acronym
                   (c) Paul Goulart
                University of Oxford, 2022
-------------------------------------------------------------

problem:
  variables     = 2
  constraints   = 5
  nnz(P)        = 2
  nnz(A)        = 6
  cones (total) = 2
    : Zero        = 1,  numel = 1
    : Nonnegative = 1,  numel = 4

settings:
  linear algebra: direct / qdldl, precision: Float64
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
  static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step
---------------------------------------------------------------------------------------------
  0  -3.2878e-01  -1.8973e+01  1.86e+01  9.76e-17  9.02e-17  1.00e+00  1.77e+00   ------
  1  -5.7453e-01  -1.8040e+00  1.23e+00  4.81e-17  3.23e-17  6.08e-02  1.15e-01  9.44e-01
  2  -6.4203e-01  -7.5345e-01  1.11e-01  1.63e-17  3.25e-17  5.21e-03  8.97e-03  9.90e-01
  3  -6.4286e-01  -6.4541e-01  2.55e-03  1.08e-16  3.26e-17  1.19e-04  2.00e-04  9.90e-01
  4  -6.4286e-01  -6.4288e-01  2.63e-05  4.89e-17  7.99e-18  1.23e-06  2.06e-06  9.90e-01
  5  -6.4286e-01  -6.4286e-01  2.63e-07  8.82e-15  1.81e-14  1.23e-08  2.06e-08  9.90e-01
  6  -6.4286e-01  -6.4286e-01  2.63e-09  2.30e-16  4.02e-16  1.23e-10  2.06e-10  9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time =  603μs
>>> Clarabel - Results
Status: SOLVED
Iterations: 6
Objective: -0.6429
Solve time:  603μs

result.x
2-element Vector{Float64}:
 0.4285714281988429
 0.2142857140994215

Using JuMP

We can solve the same problem a little more easily by using Clarabel.jl as the backend solver within JuMP. Here is the same problem again:

using Clarabel, JuMP

model = JuMP.Model(Clarabel.Optimizer)
set_optimizer_attribute(model, "verbose", true)

@variable(model, x[1:2])
@constraint(model, x[1] == 2x[2])
@constraint(model,  -1 .<= x .<= 1)
@objective(model, Min, 3x[1]^2 + 2x[2]^2 - x[1] - 4x[2])

optimize!(model)
-------------------------------------------------------------
           Clarabel.jl v0.8.1  -  Clever Acronym
                   (c) Paul Goulart
                University of Oxford, 2022
-------------------------------------------------------------

problem:
  variables     = 2
  constraints   = 5
  nnz(P)        = 2
  nnz(A)        = 6
  cones (total) = 2
    : Zero        = 1,  numel = 1
    : Nonnegative = 1,  numel = 4

settings:
  linear algebra: direct / qdldl, precision: Float64
  max iter = 200, time limit = Inf,  max step = 0.990
  tol_feas = 1.0e-08, tol_gap_abs = 1.0e-08, tol_gap_rel = 1.0e-08,
  static reg : on, ϵ1 = 1.0e-08, ϵ2 = 4.9e-32
  dynamic reg: on, ϵ = 1.0e-13, δ = 2.0e-07
  iter refine: on, reltol = 1.0e-13, abstol = 1.0e-12,
               max iter = 10, stop ratio = 5.0
  equilibrate: on, min_scale = 1.0e-04, max_scale = 1.0e+04
               max iter = 10

iter    pcost        dcost       gap       pres      dres      k/t        μ       step
---------------------------------------------------------------------------------------------
  0  -3.2878e-01  -1.8973e+01  1.86e+01  9.48e-17  1.50e-16  1.00e+00  1.77e+00   ------
  1  -5.7453e-01  -1.8040e+00  1.23e+00  4.89e-17  6.32e-17  6.08e-02  1.15e-01  9.44e-01
  2  -6.4203e-01  -7.5345e-01  1.11e-01  5.91e-17  3.34e-17  5.21e-03  8.97e-03  9.90e-01
  3  -6.4286e-01  -6.4541e-01  2.55e-03  1.14e-16  0.00e+00  1.19e-04  2.00e-04  9.90e-01
  4  -6.4286e-01  -6.4288e-01  2.63e-05  4.71e-17  3.20e-17  1.23e-06  2.06e-06  9.90e-01
  5  -6.4286e-01  -6.4286e-01  2.63e-07  8.82e-15  1.81e-14  1.23e-08  2.06e-08  9.90e-01
  6  -6.4286e-01  -6.4286e-01  2.63e-09  1.96e-16  3.17e-16  1.23e-10  2.06e-10  9.90e-01
---------------------------------------------------------------------------------------------
Terminated with status = solved
solve time =  576μs

Here is the solution

JuMP.value.(x)
2-element Vector{Float64}:
 0.4285714281988431
 0.21428571409942154

and the solver termination status again

JuMP.termination_status(model)
OPTIMAL::TerminationStatusCode = 1

Using Convex.jl

We can likewise solve the same problem a using Clarabel.jl as the backend solver within Convex.jl. Here is the same problem one more time:

using Clarabel, Convex

x = Variable(2)
problem = minimize(3square(x[1]) + 2square(x[2]) - x[1] - 4x[2])
problem.constraints = [x[1] == 2x[2]]
problem.constraints += [x >= -1; x <= 1]
solve!(problem, Clarabel.Optimizer; silent_solver = false)
Problem statistics
  problem is DCP         : true
  number of variables    : 1 (2 scalar elements)
  number of constraints  : 3 (5 scalar elements)
  number of coefficients : 12
  number of atoms        : 23

Solution summary
  termination status : OPTIMAL
  primal status      : FEASIBLE_POINT
  dual status        : FEASIBLE_POINT
  objective value    : -0.6429

Expression graph
  minimize
   └─ + (convex; real)
      ├─ * (convex; positive)
      │  ├─ [3;;]
      │  └─ qol_elem (convex; positive)
      │     ├─ …
      │     └─ …
      ├─ * (convex; positive)
      │  ├─ [2;;]
      │  └─ qol_elem (convex; positive)
      │     ├─ …
      │     └─ …
      ├─ Convex.NegateAtom (affine; real)
      │  └─ index (affine; real)
      │     └─ …
      ⋮
  subject to
   ├─ == constraint (affine)
   │  └─ + (affine; real)
   │     ├─ index (affine; real)
   │     │  └─ …
   │     └─ Convex.NegateAtom (affine; real)
   │        └─ …
   ├─ ≥ constraint (affine)
   │  └─ + (affine; real)
   │     ├─ 2-element real variable (id: 348…882)
   │     └─ Convex.NegateAtom (constant; positive)
   │        └─ …
   └─ ≤ constraint (affine)
      └─ + (affine; real)
         ├─ 2-element real variable (id: 348…882)
         └─ Convex.NegateAtom (constant; negative)
            └─ …

Here is our solution

evaluate(x)
2-element Vector{Float64}:
 0.4285731865024157
 0.21428659325133423

and the solver termination status again

problem.status
OPTIMAL::TerminationStatusCode = 1
Warning

Note that in the Clarabel.jl output that follows the call to solve! using Convex.jl, the problem posed to the solver has been converted to a second-order cone program with a linear objective. You can see this because now nnz(P) == 0 (there is no quadratic term in the objective) and the solver reports two second order cone constraints.

Although the solution will be the same, the required number of iterations and solve time are slightly higher. When solving problems with quadratic objectives in Clarabel.jl, it is generally preferable to use either the native Clarabel.jl interface or JuMP, both of which handle quadratic terms in the objective directly.

This page was generated using Literate.jl.