# Package in julia for solving optimization problems with cubic objective function
**URL:** https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055
**Category:** Optimization (Mathematical)
**Created:** [April 17, 2024, 9:47am UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055 "2024-04-17T09:47:46Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@heiwie](https://discourse.julialang.org/u/heiwie)
#### Post date: [April 17, 2024, 9:47am UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/1 "2024-04-17T09:47:46Z")
Hello, I am trying to solve an optimization problem which looks like this:
```julia
@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:
````julia
@constraint(model, sum(probs[18][m][n] * (beta_18[m, n])^2 for n in 1:num_nodes) <= A^2)```
````
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [April 17, 2024, 9:49am UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/2 "2024-04-17T09:49:37Z")
How large are your variables? Can you post a complete example?
---
### Author:  [@heiwie](https://discourse.julialang.org/u/heiwie)
#### Post date: [April 17, 2024, 9:57am UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/3 "2024-04-17T09:57:39Z")
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.
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [April 17, 2024, 10:08am UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/4 "2024-04-17T10:08:41Z")
Sorry I meant the dimension of the variables ^^
---
### Author:  [@heiwie](https://discourse.julialang.org/u/heiwie)
#### Post date: [April 17, 2024, 10:12am UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/5 "2024-04-17T10:12:22Z")
Oh ok! Each beta is a scalar.
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [April 17, 2024, 10:14am UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/6 "2024-04-17T10:14:18Z")
I meant how big is `num_nodes` 😉 aka the number of variables total
---
### Author:  [@heiwie](https://discourse.julialang.org/u/heiwie)
#### Post date: [April 17, 2024, 10:17am UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/7 "2024-04-17T10:17:38Z")
Ah I see, num\_nodes equals 25.
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [April 17, 2024, 3:42pm UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/8 "2024-04-17T15:42:22Z")
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?
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [April 17, 2024, 3:45pm UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/9 "2024-04-17T15:45:04Z")
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](https://jump.dev/JuMP.jl/stable/installation/#Supported-solvers) and keeping only the ones that support “NLP”
---
### Author:  [@heiwie](https://discourse.julialang.org/u/heiwie)
#### Post date: [April 17, 2024, 4:34pm UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/10 "2024-04-17T16:34:46Z")
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.
---
### Author:  [@heiwie](https://discourse.julialang.org/u/heiwie)
#### Post date: [April 17, 2024, 4:36pm UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/11 "2024-04-17T16:36:00Z")
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?
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [April 17, 2024, 4:50pm UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/12 "2024-04-17T16:50:07Z")
> [@heiwie](#):
>
> Say that the problem had a min objective function instead, would the objective function still be convex in this case?
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).
---
### Author:  [@heiwie](https://discourse.julialang.org/u/heiwie)
#### Post date: [April 17, 2024, 4:59pm UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/13 "2024-04-17T16:59:40Z")
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!
---
### Author:  [@mohamed82008](https://discourse.julialang.org/u/mohamed82008)
#### Post date: [April 17, 2024, 9:32pm UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/14 "2024-04-17T21:32:49Z")
Try [GitHub - PSORLab/EAGO.jl: A development environment for robust and global optimization](https://github.com/PSORLab/EAGO.jl) or [GitHub - lanl-ansi/Alpine.jl: A Julia/JuMP-based Global Optimization Solver for Non-convex Programs](https://github.com/lanl-ansi/Alpine.jl).
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [April 18, 2024, 12:46am UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/15 "2024-04-18T00:46:09Z")
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](https://www.gurobi.com/documentation/11.0/refman/quadratic_constraints.html) with the [NonConvex](https://www.gurobi.com/documentation/11.0/refman/nonconvex.html) settings. In JuMP, this would be set as
> [@Calling Mathematica into Julia to symbolically solve a system of non-linear equations](https://discourse.julialang.org/t/calling-mathematica-into-julia-to-symbolically-solve-a-system-of-non-linear-equations/80518/4):
>
> ```julia
> model = JuMP.Model(Gurobi.Optimizer)
> JuMP.set_optimizer_attribute(model, "NonConvex", 2)
>
> ```
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.
---
### Author:  [@heiwie](https://discourse.julialang.org/u/heiwie)
#### Post date: [April 19, 2024, 8:26am UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/16 "2024-04-19T08:26:15Z")
> [@jd-foster](#):
>
> 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](https://www.gurobi.com/documentation/11.0/refman/quadratic_constraints.html) with the [NonConvex](https://www.gurobi.com/documentation/11.0/refman/nonconvex.html) settings. In JuMP, this would be set as
Is that only for Gurobi 11.0 or higher? Because when I try to apply this to a different program with the following constraint:
```julia
@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:
```julia
@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.
---
### Author:  [@jd-foster](https://discourse.julialang.org/u/jd-foster)
#### Post date: [April 19, 2024, 10:42am UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/17 "2024-04-19T10:42:36Z")
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.
---
### Author:  [@mtanneau](https://discourse.julialang.org/u/mtanneau)
#### Post date: [April 20, 2024, 2:41pm UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/18 "2024-04-20T14:41:18Z")
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](https://jump.dev/JuMP.jl/stable/tutorials/conic/tips_and_tricks/#Second-Order-Cone)
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](https://jump.dev/JuMP.jl/stable/tutorials/conic/tips_and_tricks/#Second-Order-Cone).
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.
---
### Author:  [@heiwie](https://discourse.julialang.org/u/heiwie)
#### Post date: [April 21, 2024, 9:28am UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/19 "2024-04-21T09:28:42Z")
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().
---
### Author:  [@mbesancon](https://discourse.julialang.org/u/mbesancon)
#### Post date: [April 21, 2024, 5:04pm UTC](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055/20 "2024-04-21T17:04:14Z")
> How would such a constraint look like in JuMP?
Something like this:
```julia
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:
```julia
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 🙂
[Next page](https://discourse.julialang.org/t/package-in-julia-for-solving-optimization-problems-with-cubic-objective-function/113055.md?page=2)