Package in julia for solving optimization problems with cubic objective function

Specific Domains Optimization (Mathematical)
1

Hello, I am trying to solve an optimization problem which looks like this:

@objective(model, Max, sum(beta_2[n] * node_info[2][1][n][1] * prod_levels[2][1][n] * uncon_probs[2][1][n] for n in 1:num_nodes) 
                        + sum(sum(beta_2[m] * beta_6[m, n] * static_beta^3 * node_info[6][m][n][1] * prod_levels[6][m][n] * uncon_probs[6][m][n] for m in 1:num_nodes) for n in 1:num_nodes)
                        + sum(sum(beta_2[get_n_index_first_split(m)] * beta_6[get_n_index_first_split(m), get_n_index_second_split(m)] * beta_18[m, n] * static_beta^14 
                        * uncon_probs[18][m][n] * node_info[18][m][n][1] * prod_levels[18][m][n] for m in 1:num_nodes^2) for n in 1:num_nodes))

The details of the problem are unimportant, however as you can see the objective function is cubic since the decision variables (beta_2, beta_6 and beta_18) are multiplied together. Is there any package in Julia which can solve such optimization problems to a global optima?

I should also mention that the problem also contains quadratic constraints:

@constraint(model, sum(probs[18][m][n] * (beta_18[m, n])^2 for n in 1:num_nodes) <= A^2)```
2

How large are your variables? Can you post a complete example?

3

The decision variables, i.e. the betas, must be larger than zero but there isnt really an upper limit. Generally they will lie in the range [0, 10], but in theory there is nothing stopping them from being very large. If you are referring to the other variables in the objective function, then the unconditional probabilities will naturally be below 1 and the “node_info” represents prices of a scenario tree and will generally lie in the interval 30-100. The static_beta is 1.00038.

4

Sorry I meant the dimension of the variables ^^

5

Oh ok! Each beta is a scalar.

6

I meant how big is num_nodes :wink: aka the number of variables total

7

Ah I see, num_nodes equals 25.

1 Like
8

So if the decision variables are larger than zero, does that mean the objective function is convex even though it is cubic?
And the quadratic constraint defines a convex feasible set?

9

If that’s right, then basically any solver worth its salt will converge to the optimal solution. JuMP.jl + Ipopt.jl is probably a good solution, but you can have more options by checking out the list Installation Guide · JuMP and keeping only the ones that support “NLP”

10

I am afraid I am a little shaky on convexity for cubic functions. Say that the problem had a min objective function instead, would the objective function still be convex in this case? The constraints should define a convex set, yes.

11

Thanks! I am familiar with Ipopt, and have tried to run a few of the other NLP packages such as MadNLP, but have found that it converges much slower on the solution. Is it generally known that Ipopt is “best practice” in most cases in terms of runtime, or is it perhaps more dependent on the particular problem you aim to solve?

12

Oh crap I didn’t check but you’re maximizing, not minimizing. Forget anything I said about convexity: the objective is convex indeed but it’s not very useful information. It only means the optimum will be found at a boundary of the domain.

That means there are few (if any) packages who can guarantee global optimality. And you should probably look into solvers designed for global optimization over low-dimensional spaces, rather than local ones like Ipopt. I don’t know much about best practices in that case. Basically in the convex case, best practices don’t matter as much cause everyone will find the right solution (after varying runtime).

1 Like
13

I see, I actually also define a very similar optimization problem in which I minimize the same objective function, so it is good to know that it should work with Ipopt or something similar in this case. Thank you!

1 Like
14

Try GitHub - PSORLab/EAGO.jl: A development environment for robust and global optimization or GitHub - lanl-ansi/Alpine.jl: A Julia/JuMP-based Global Optimization Solver for Non-convex Programs.

15

Just to weigh in here, if you have access to Gurobi as a solver, then you can solve problems to global optimality provided they are quadratic constraints with the NonConvex settings. In JuMP, this would be set as

Of course, the issue is having cubic / trilinear terms. A reformulation to quadratic / bilinear form is possible, such as rewriting z = x_1 x_2 x_3 as y = x_1 x_2 and z = y x_3.

16

Is that only for Gurobi 11.0 or higher? Because when I try to apply this to a different program with the following constraint:

@constraint(model, B[t] * (weight[t-1, m] - weight[t, node_index]) + eta[t, m, n] 
                >= A * sum(probabilities[t+1][node_index][v] * (eta[t+1, node_index, v])^2 for v in 1:num_children_nodes[t+2])^(1/2) - node_info[t][m][n][1] * x[t, m, n])

and objective function:

@objective(model, Min, -node_info[1][1][1][1] * x[1, 1, 1] + B[1] * weight[1, 1] + A * sum(probabilities[2][1][v] * (eta[2, 1, v])^2 for v in 1:25)^(1/2))

Then it says that it does not support this type of problem. Please note here that x, weight and eta are the decision variables. I have access to Gurobi 10.0.

17

Gurobi 9 or later should work.

The issue is probably to do with the square root ^(1/2) in the constraint. The objective and constraints must be linear or quadratic for the solver to accept the problem.

1 Like
18

About that different program: I took a quick glance at the constraint you mentioned, and if I read correctly it has the form

(w_{t-1} - w_{t}) + \eta_{0} \geq A \sqrt{\sum_{i} p_{i} \eta_{i}^{2}} - x

(I removed all constants and dropped some indices to make the constraint more readable).

Assuming A \geq 0 and p_{i} \geq 0 (here p_{i} corresponds to probabilities[2][1][v] in your post), then this constraint can be written equivalently as

(w_{t-1} - w_{t}) + \eta_{0} + x \geq A \| \sqrt{p} \cdot \eta \|_{2}

where the left-hand side is linear, and the right-hand side is the euclidean norm of (\sqrt{p_{1}} \eta_{1}, ..., \sqrt{p_{k}} \eta_{k}) multiplied by A \geq 0.

Assuming that A \geq 0 and p \geq 0, then this constraint can be represented by introducing an additional variable t such that

(w_{t-1} - w_{t}) + \eta_{0} + x \geq t,

and a second-order cone constraint

t \geq \| \sqrt{A} \cdot \sqrt{p} \cdot \eta \|_{2}

Note that you can use this same t variable to replace the same sum(...)^(1/2) in your objective function. This is valid because you are minimizing and the objective function is convex (again under the assumption that A, p \geq 0).

A small modelling note: JuMP allows you to build second-order constraints directly.
If you use Gurobi, you can write a constraint t \geq \|x \|_{2} either as (t, x) \in K_{SOC} (this is the conic form that I linked above), or as a quadratic constraint t^{2} \geq \sum_{i} x_{i}^{2}, t \geq 0. Gurobi will automatically detect that the latter is convex, and handle it accordingly.

1 Like
19

Very interesting, @mtanneau! How would such a constraint look like in JuMP? I have understood that second-order cone constraints require a somewhat different formulation using SecondOrderCone().

20

How would such a constraint look like in JuMP?

Something like this:

julia> m = Model()
julia> @variable(m, x[1:3])
julia> @constraint(m, x in SecondOrderCone())

Or more generically with z, x such that z \geq \|x\|_2:

m = Model()
@variable(m, z >= 0)
@variable(m, x[1:3])

@constraint(m, [z; x] in SecondOrderCone())

For the original non-convex problem, an alternative is to optimize it with the SCIP wrapper, which is an exact global optimizer. The only reformulation you will need is a linear objective:
\max_{z,x} z \;\;\; \text{s.t.}\;\;\; z \leq \text{my_objective}(x).

As a disclaimer, I work on the solver and am a bit biased :slight_smile:

1 Like