Optimizing an objective function with a vector containing constants, variables, and functions of variables using JuMP

Specific Domains Optimization (Mathematical)
1

I’m working on a set of optimization problems where I need to maximize something that looks a lot like a quadratic form but isn’t. The following example is one such problem:

Maximize f^TAf

where f = [1, x, x^2] and A is a 3x3 PSD matrix. I believe JuMP paired with a semidefinite solver like SCS will provide everything I need for the task, but I can’t seem to get Julia to accept my program. I feel like I should be able to do something like:

model = Model(SCS.Optimizer);
@variable(model, x)
f = [1, x, x^2]
@constraint(model, c1, x <= 1);
@constraint(model, c2, x >= -1);
@objective(model, Max, f'*A*f)
optimize!(model)

but an error is thrown when I define my objective. The particular error is:

ERROR: *(::GenericQuadExpr,::GenericQuadExpr) is not defined. Are you trying to build a nonlinear problem? Make sure you use @NLconstraint/@NLobjective. If you are using an `@NL` macro and you encountered this error message, it is because you are attempting to use another unsupported function which calls this method internally.

Any suggestions will be appreciated!

2

Hi @Hyrum_Hansen, welcome to the forum.

SCS supports only linear and quadratic functions. Your problem is not convex because it contains cubic and quartic terms.

In addition, even if you could write your f, if A is PSD, then you cannot maximize the objective because it is a convex function. Do you mean Min instead?

You could reformulate your problem to this:

using JuMP, SCS
A = rand(3, 3)
A += A'
model = Model(SCS.Optimizer)
@variable(model, x)
@variable(model, y)
@constraint(model, y >= x^2)
f = [1, x, y]
@objective(model, Min, f'*A*f)
optimize!(model)

but it might not be what you’re trying to achieve.

3

Apologies, I forgot some critical information (but thank you for your quick response!). x is bounded from above and below; -1<=x<=1. Still, the problem is nonconvex, but with a finite search space there should be a solver that can handle it right? Something like Gloptipoly in Matlab I’m thinking.

4

SCS requires the problem to be convex.

You can solve the problem using Gurobi: jump-dev/Gurobi.jl · JuMP

julia> using JuMP, Gurobi

julia> A = rand(3, 3)
3×3 Matrix{Float64}:
 0.760386  0.407191  0.723833
 0.604206  0.364107  0.623264
 0.356807  0.60037   0.521208

julia> A += A'
3×3 Matrix{Float64}:
 1.52077  1.0114    1.08064
 1.0114   0.728214  1.22363
 1.08064  1.22363   1.04242

julia> model = Model(Gurobi.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Gurobi

julia> set_attribute(model, "NonConvex", 2)
Set parameter NonConvex to value 2

julia> @variable(model, -1 <= x <= 1)
x

julia> @variable(model, y)
y

julia> @constraint(model, y == x^2)
-x² + y = 0

julia> f = [1, x, y]
3-element Vector{AffExpr}:
 1
 x
 y

julia> @objective(model, Max, f' * A * f)
0.7282140530433079 x² + 2.447268930048988 x*y + 1.0424158567197273 y² + 2.0227956982742628 x + 2.16127995587574 y + 1.520771522270953

julia> optimize!(model)
Set parameter NonConvex to value 2
Gurobi Optimizer version 10.0.0 build v10.0.0rc2 (mac64[x86])

CPU model: Intel(R) Core(TM) i5-8259U CPU @ 2.30GHz
Thread count: 4 physical cores, 8 logical processors, using up to 8 threads

Optimize a model with 0 rows, 2 columns and 0 nonzeros
Model fingerprint: 0x6fc2d173
Model has 3 quadratic objective terms
Model has 1 quadratic constraint
Coefficient statistics:
  Matrix range     [0e+00, 0e+00]
  QMatrix range    [1e+00, 1e+00]
  QLMatrix range   [1e+00, 1e+00]
  Objective range  [2e+00, 2e+00]
  QObjective range [1e+00, 5e+00]
  Bounds range     [1e+00, 1e+00]
  RHS range        [0e+00, 0e+00]

Continuous model is non-convex -- solving as a MIP

Presolve time: 0.00s
Presolved: 7 rows, 5 columns, 15 nonzeros
Presolved model has 3 bilinear constraint(s)
Variable types: 5 continuous, 0 integer (0 binary)

Root relaxation: objective 9.922746e+00, 2 iterations, 0.00 seconds (0.00 work units)

    Nodes    |    Current Node    |     Objective Bounds      |     Work
 Expl Unexpl |  Obj  Depth IntInf | Incumbent    BestBd   Gap | It/Node Time

*    0     0               0       9.9227460    9.92275  0.00%     -    0s

Explored 1 nodes (2 simplex iterations) in 0.01 seconds (0.00 work units)
Thread count was 8 (of 8 available processors)

Solution count 1: 9.92275 

Optimal solution found (tolerance 1.00e-04)
Best objective 9.922746016233e+00, best bound 9.922746016233e+00, gap 0.0000%

User-callback calls 100, time in user-callback 0.00 sec

julia> value(x)
1.0
1 Like