# Solving nonlinear constrained optimization problem
**URL:** https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255
**Category:** Optimization (Mathematical)
**Tags:** question, jump
**Created:** [March 6, 2024, 5:58pm UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255 "2024-03-06T17:58:20Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@rtapia](https://discourse.julialang.org/u/rtapia)
#### Post date: [March 6, 2024, 5:58pm UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/1 "2024-03-06T17:58:20Z")
Hello,
I am trying to solve the following nonconvex problem in Julia using Optim.jl:
min x’_P_x
s.t: 1 -x’\*x \<=0
where P is a positive definite matrix. I have defined the following
```julia
using JuMP, Optim
n = 1500;
A = 10*rand(Float64,(n,n));
P = A'*A;
m = Model(Optim.Optimizer)
set_optimizer_attribute(m, "method", BFGS())
@variable(m, x[1:n]);
@objective(m,Min, x'*P*x);
f = 1 - x'*x;
@constraint(m, f1, f <= 0)
optimize!(m)
```
But I get a a stackoverflowerror, so I am not sure what I am missing. I know there are multiple ways of defining the problem in Optim. But the only constrained example I see is the boxconstraint: [Interior point Newton - Optim.jl (julianlsolvers.github.io)](https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/)
Thanks
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 6, 2024, 6:10pm UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/2 "2024-03-06T18:10:08Z")
> [@rtapia](#):
>
> I am trying to solve the following nonconvex problem in Julia using [Optim.jl](https://juliahub.com/ui/Packages/General/Optim):
> min x’_P_x
> s.t: 1 -x’\*x \<=0
Isn’t this analytically solvable? According to the [min–max theorem](https://en.wikipedia.org/wiki/Min-max_theorem), your minimum will be the smallest eigenvalue of P, achieved when x is a corresponding eigenvector (normalized to unit length).
If you are just using this as a test case for JuMP, I don’t think the BFGS algorithm in Optim.jl supports arbitrary nonlinear constraints? You may need some other engine, like Nonconvex.jl or NLopt.jl or Ipopt.jl or …
---
### Author:  [@rtapia](https://discourse.julialang.org/u/rtapia)
#### Post date: [March 6, 2024, 10:39pm UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/3 "2024-03-06T22:39:43Z")
Yes this is just a test case, according to [Optim.jl](https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/) optim can handle generic nonlinear constraints using IPNewton.
I have since corrected the code to:
```julia
set_optimizer_attribute(m, "method", IPNewton());
x0=ones(n,1);
@variable(m, x[1:n]);
#@variable(m, x[i=1:n], start = x0[i]);
set_start_value.(x,1);
```
I will try Nonconvex.jl, I have tried NLopt and Ipopt. Just getting familiar with the available packages in Julia.
---
### Author:  [@rtapia](https://discourse.julialang.org/u/rtapia)
#### Post date: [March 7, 2024, 2:35am UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/4 "2024-03-07T02:35:24Z")
Hello,
I have formed the following problem using Optim.jl with a nonconvex constraint:
min x’_P_x
1-x’\*x \<=0
According to: [Optim.jl](https://julianlsolvers.github.io/Optim.jl/stable/#examples/generated/ipnewton_basics/)
Optim should be able to solve it using IPNewton?
the code is as follows:
```julia
using JuMP, Optim, LineSearches
n = 10;
A = 10*rand(Float64,(n,n));
P = A'*A;
# cost function
f0(x) = x'*P*x;
# gradient of cost function
function f0_grad!(gf0,x)
for j=1:n
gf0[j] = 2*P[j,:]'*x;
end
end
# Hessian is a constant
function f0_hess!(h,x)
for j=1:n
for i=1:n
h[i,j] = 2*P[i,j];
end
end
end
# nonlinear constraint
con_c!(c,x) = (c[1] = 1 - x'*x);
# Jacobian of constraint
function con_jacobian!(J,x)
for j=1:n
J[j] = -2*x[j];
end
end
# Hessian of Constraint
function con_h!(h2,x,λ)
#h = zeros(n,n);
for j=1:n
h2[j,j] = -2*λ[1];
end
end
#upper and lower bounds
lx = -10*ones(n); ux = 10*ones(n);
lc = [-Inf]; uc = [-1];
#initial condition
x0 = 1.1*ones(n);
df = TwiceDifferentiable(f0, f0_grad!, f0_hess!, x0)
dfc = TwiceDifferentiableConstraints(con_c!, con_jacobian!, con_h!, lx, ux, lc, uc)
res = optimize(df, dfc, x0, IPNewton())
```
The solution states “success” but the solution is incorrect. The solution of x should correspond to the eigenvector (normalized) associated with the lowest eigenvalue.
Thanks,
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 7, 2024, 2:44am UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/5 "2024-03-07T02:44:12Z")
How far off is the minimum from the minimal eigenvalue? I notice that you don’t set a tolerance or any other convergence criteria — maybe you just need to run for more iterations / lower the tolerance?
---
### Author:  [@rtapia](https://discourse.julialang.org/u/rtapia)
#### Post date: [March 7, 2024, 3:11am UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/6 "2024-03-07T03:11:07Z")
Its far off, I increased the number iterations and lowered the tolerance, but it seems to make it worse. The value should be small, but I get a solution ~e3.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [March 7, 2024, 3:37am UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/7 "2024-03-07T03:37:26Z")
You could do:
```julia
using JuMP, Ipopt
function main(A)
n = size(A, 1)
model = Model(Ipopt.Optimizer)
x0 = zeros(n)
_, i = findmin(vec(sum(A; dims = 2)))
x0[i] = 1.0
y0 = A * x0
@variable(model, x[i in 1:n], start = x0[i]);
@variable(model, y[i in 1:n], start = y0[i]);
@objective(model, Min, sum(y.^2))
@constraint(model, A * x .== y)
@constraint(model, 1 <= sum(x.^2))
optimize!(model)
@assert is_solved_and_feasible(model)
return value.(x)
end
n = 1_500
A = 10 * rand(Float64, (n, n))
main(A)
```
But note that Ipopt assumes the problem is convex, so you are not guaranteed to find a global minima.
---
### Author:  [@rtapia](https://discourse.julialang.org/u/rtapia)
#### Post date: [March 7, 2024, 4:46am UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/8 "2024-03-07T04:46:02Z")
Well that is working much better, I modified it to use @NLconstraint, that should technically work. I am trying to see if increasing tolerance could have an impact on the solution, but I am not sure what is the correct way to implement.
```julia
using JuMP, Ipopt, LinearAlgebra
x1 = zeros(n);
function main(A)
n = size(A, 1)
model = Model(Ipopt.Optimizer)
set_attribute(model,"max_iter", 10000)
model.xtol_rel = 1e-6;
set_optimizer_attribute(model,"acceptable_constr_viol_tol", 1e-10)
set_optimizer_attribute(model,"constr_viol_tol", 1e-6)
x0 = zeros(n)
x0 = ones(n)
y0 = A * x0
@variable(model, x[i in 1:n], start = x0[i]);
@variable(model, y[i in 1:n], start = y0[i]);
@objective(model, Min, sum(y.^2))
f = x'*x;
@constraint(model, A * x .== y)
@NLconstraint(model, 1 <= f)
JuMP.optimize!(model)
@assert is_solved_and_feasible(model)
global x1 =value.(x);
return JuMP.objective_value(model)
end
n = 10;
A = 10 * rand(Float64, (n, n));
A = A'*A;
main(A)
V = eigvals(A);
print("solution should be = ",V[1]);
```
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [March 7, 2024, 6:17am UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/9 "2024-03-07T06:17:20Z")
> [@rtapia](#):
>
> I modified it to use @NLconstraint
Please don’t 😆 the `@NL` interface is now legacy (see [ANN: JuMP v1.15 is released](https://discourse.julialang.org/t/ann-jump-v1-15-is-released/103875)) and it should not be encouraged for new uses. Also, it means the constraint is added as a generic nonlinear constraint. If you use `@constraint` JuMP and Ipopt can detect it as quadratic.
> global x1 =value.(x);
Also, don’t use `global`. Return the value from the function with
```Julia
return value.(x), JuMP.objective_value(model)
```
and then call it as
```julia
x1, obj = main(A)
```
> `model.xtol_rel = 1e-6;`
This I not valid JuMP syntax.
You can see the full list of Ipopt options here: [Ipopt: Ipopt Options](https://coin-or.github.io/Ipopt/OPTIONS.html#OPTIONS_REF)
The way to set an option is `set_attribute(model, "tol", 1e-6)`.
---
### Author:  [@kellertuer](https://discourse.julialang.org/u/kellertuer)
#### Post date: [March 7, 2024, 7:12am UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/10 "2024-03-07T07:12:16Z")
Just a guess, but you are mainly tackling the Rayleigh quotient here (up to scaling, and your constraints).
Your constraints basically say you want t norm (squared) larger than 1. So the minimum will lie _on_ the boundary (you could also take norme _equal_) to one.
If you try to solve the (unconstraint) Rayleigh problem with any Newton based iteration, one can compute that the update will boil down to multiplying your (nearly) eigenvector with a (positive) scalar. Might that be a cause here?
---
### Author:  [@rtapia](https://discourse.julialang.org/u/rtapia)
#### Post date: [March 7, 2024, 5:16pm UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/11 "2024-03-07T17:16:06Z")
Clearly a noob here 🤣, thanks for the corrections. I learned based off older versions.
@kellertuer I have solved this with a Primal-Dual interior point solver
---
### Author:  [@rtapia](https://discourse.julialang.org/u/rtapia)
#### Post date: [March 7, 2024, 9:16pm UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/12 "2024-03-07T21:16:01Z")
for anyone interested I got the exact solution (or near exact) using Nonconvex.jl
```julia
using Nonconvex, NonconvexIpopt, LinearAlgebra, Printf
n = 10;
A = 10*rand(Float64,(n,n));
P = A'*A;
f(x) = x'*P*x;
f1(x) = 1- x'*x;
lx = -500*ones(n); ux = 500*ones(n);
model = Model(f)
addvar!(model, lx,ux)
add_ineq_constraint!(model, x -> f1(x))
x0 = 1.5*ones(n);
#alg = NLoptAlg(:LD_MMA)
alg = IpoptAlg()
#alg = MMA02()
options = IpoptOptions(tol=1e-10, first_order = false)
#options = MMAOptions()
res = optimize(model, alg, x0, options =options, convcriteria = KKTCriteria())
V = eigvals(P);
@printf("Answer should be = %0.4f\n",V[1])
@printf("Answer is = %0.4f\n", res.minimum)
```
---
### Author:  [@sshin23](https://discourse.julialang.org/u/sshin23)
#### Post date: [March 9, 2024, 6:20am UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/13 "2024-03-09T06:20:27Z")
If you want something scalable, try using [MadNLP](https://github.com/MadNLP/MadNLP.jl).
```julia
using MadNLP, LinearAlgebra
struct MyProblem{T, VT, MT} <: MadNLP.AbstractNLPModel{T, VT}
Q::MT
counters::MadNLP.NLPModels.Counters
meta::MadNLP.NLPModelMeta{T,VT}
end
MyProblem(Q) = MyProblem(
Q,
MadNLP.NLPModels.Counters(),
MadNLP.NLPModelMeta(
size(Q,1);
x0 = fill!(similar(Q,size(Q,1)), 10),
ncon = 1,
ucon = fill!(similar(Q,1), 0)
)
)
MadNLP.obj(m::MyProblem, x::AbstractVector) = dot(x,m.Q,x)
MadNLP.grad!(m::MyProblem, x::AbstractVector, y::AbstractVector) = (mul!(y,m.Q,x); y.*=2)
MadNLP.cons!(m::MyProblem, x::AbstractVector, y::AbstractVector) = fill!(y, 1 - dot(x,x))
MadNLP.jac_dense!(m::MyProblem, x::AbstractVector, j::AbstractMatrix) = (j .= x'; j.*=-2)
function MadNLP.hess_dense!(m::MyProblem, x::AbstractVector, y::AbstractVector, h::AbstractMatrix; obj_weight= one(eltype(x)))
copyto!(h, m.Q)
h .*= obj_weight
view(h,1: n+1: n^2) .-= y[1]
end
n = 100;
A = 10*rand(Float64,(n,n));
P = A'*A;
m = MyProblem(P)
madnlp(m; lapack_algorithm=MadNLP.CHOLESKY, max_iter= 10000)
# got an Nvidia gpu?
using MadNLPGPU, CUDA
MadNLP.obj(m::MyProblem, x::CuVector) = dot(x,m.Q * x)
function MadNLP.hess_dense!(m::MyProblem, x::CuVector, y::CuVector, h::CuMatrix; obj_weight= one(eltype(x)))
copyto!(h, m.Q)
h .*= obj_weight
CUDA.@allowscalar view(h,1: n+1: n^2) .-= y[1]
end
m = MyProblem(CuArray(P))
madnlp(m; lapack_algorithm=MadNLP.CHOLESKY, max_iter= 10000)
```
---
### Author:  [@rtapia](https://discourse.julialang.org/u/rtapia)
#### Post date: [March 9, 2024, 8:57pm UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/14 "2024-03-09T20:57:33Z")
I’ve had issues with MadNLP, get the following error. When I search madnlp function it does not have the option to specify the algorithm.
---
### Author:  [@sshin23](https://discourse.julialang.org/u/sshin23)
#### Post date: [March 10, 2024, 2:38am UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/15 "2024-03-10T02:38:12Z")
Could you try with a more recent version of MadNLP? This seems to be using `v0.6.0` or older. The most recent version is `v0.8.0`.
---
### Author:  [@rtapia](https://discourse.julialang.org/u/rtapia)
#### Post date: [March 10, 2024, 3:31am UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/16 "2024-03-10T03:31:13Z")
I get the following error when trying to update:
---
### Author:  [@sshin23](https://discourse.julialang.org/u/sshin23)
#### Post date: [March 10, 2024, 3:47am UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/17 "2024-03-10T03:47:18Z")
Could you try to install MadNLP in a separate project environment? Most likely some other package is blocking the upgrade.
```julia
cd to/some/tmp/directory
julia --project=.
]add MadNLP
```
---
### Author:  [@rtapia](https://discourse.julialang.org/u/rtapia)
#### Post date: [March 10, 2024, 5:49am UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/18 "2024-03-10T05:49:38Z")
That worked, but the output is incorrect, answer should be the smallest eigenvalue of P.
I am a little confused on the line:
MadNLP.grad!(m::MyProblem, x::AbstractVector, y::AbstractVector) = (mul!(y,m.Q,x); y.\*=2)
is this stating the gradient is [Px; 2Px] ? or just 2Px?
---
### Author:  [@sshin23](https://discourse.julialang.org/u/sshin23)
#### Post date: [March 10, 2024, 8:59pm UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/19 "2024-03-10T20:59:17Z")
It first put `m.Q*x` into `y`, and then multiply the entries of `y` by 2, which has the same effect of
`y .= 2 .* m.Q *x`, but just reduces memory allocation.
MadNLP is a local solver, like Ipopt and many other NLP solvers. So, it doesn’t aim to find the global minimum.
Just curious, what’s the motivation for solving this as an optimization problem rather than using more specialized methods (e.g., power iteration)?
---
### Author:  [@rtapia](https://discourse.julialang.org/u/rtapia)
#### Post date: [March 10, 2024, 9:14pm UTC](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255/20 "2024-03-10T21:14:11Z")
Hello,
The objective value was smaller than the possible minimum. I am just testing, then I will scale to a larger optimization problem that is nonconvex. The two problems are not related but would like to seek a “reliable” solver.
Thanks
[Next page](https://discourse.julialang.org/t/solving-nonlinear-constrained-optimization-problem/111255.md?page=2)