# Solvers fail on nonlinear problem (which has solutions)
**URL:** https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051
**Category:** Optimization (Mathematical)
**Tags:** jump
**Created:** [January 24, 2019, 5:52pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051 "2019-01-24T17:52:43Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@Torkel](https://discourse.julialang.org/u/Torkel)
#### Post date: [January 24, 2019, 5:52pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/1 "2019-01-24T17:52:43Z")
(quite long text, but I wanted to put down everything I had on the problem)
I have been having this set of non-linear equations which I have been trying to solve for a while now. I initially tried `NLsolve.jl`, however, it does not allow for specifying variables as non-negative. I have instead phrased it as an optimisation problem and used JuMP.
I have defined this function which solves the problem for me (to easily allow for different parameter value and solvers (I can also define a desired initial point for the solver, if I want).
```julia
function my_problem(p,solver;xstart=fill(1.,10))
m = JuMP.Model(solver=solver)
JuMP.@variable m x1 >= 0.0 start = xstart[1]
JuMP.@variable m x2 >= 0.0 start = xstart[2]
JuMP.@variable m x3 >= 0.0 start = xstart[3]
JuMP.@variable m x4 >= 0.0 start = xstart[4]
JuMP.@variable m x6 >= 0.0 start = xstart[6]
JuMP.@variable m x5 >= 0.0 start = xstart[5]
JuMP.@variable m x7 >= 0.0 start = xstart[7]
JuMP.@variable m x8 >= 0.0 start = xstart[8]
JuMP.@variable m x9 >= 0.0 start = xstart[9]
JuMP.@variable m x10 >= 0.0 start = xstart[10]
JuMP.@NLconstraint m -1 * p[17] * x1 + -2 * p[1] * (x1 ^ 2 / 2) + 2 * p[2] * x2 + p[21] * p[18] * ((1 + p[19] * x7) / (p[20] + x7)) == 0
JuMP.@NLconstraint m -1 * p[17] * x2 + p[1] * (x1 ^ 2 / 2) + -1 * p[2] * x2 + -1 * p[4] * x2 * x4 + p[9] * x3 + p[14] * x3 + -1 * p[6] * x2 * x7 + p[11] * x8 == 0
JuMP.@NLconstraint m -1 * p[17] * x3 + p[4] * x2 * x4 + -1 * p[9] * x3 + -1 * p[5] * x3 * x4 + p[10] * x5 + -1 * p[14] * x3 + p[15] * x5 + p[7] * x8 * x4 + -1 * p[12] * x3 * x7 == 0
JuMP.@NLconstraint m -1 * p[17] * x4 + -1 * p[4] * x2 * x4 + p[9] * x3 + -1 * p[5] * x3 * x4 + p[10] * x5 + -1 * p[7] * x8 * x4 + p[12] * x3 * x7 + p[16] * x9 + p[22] * p[18] * ((1 + p[19] * x7) / (p[20] + x7)) == 0
JuMP.@NLconstraint m -1 * p[17] * x5 + p[5] * x3 * x4 + -1 * p[10] * x5 + -1 * p[15] * x5 == 0
JuMP.@NLconstraint m -1 * p[17] * x6 + p[14] * x3 + p[15] * x5 + -1 * p[8] * x6 * x10 + p[13] * x9 == 0
JuMP.@NLconstraint m -1 * p[17] * x7 + -1 * p[6] * x2 * x7 + p[11] * x8 + p[7] * x8 * x4 + -1 * p[12] * x3 * x7 + p[18] * ((1 + p[19] * x7) / (p[20] + x7)) == 0
JuMP.@NLconstraint m -1 * p[17] * x8 + p[6] * x2 * x7 + -1 * p[11] * x8 + -1 * p[7] * x8 * x4 + p[12] * x3 * x7 == 0
JuMP.@NLconstraint m -1 * p[17] * x9 + p[8] * x6 * x10 + -1 * p[13] * x9 + -1 * p[16] * x9 == 0
JuMP.@constraint m x10 + x9 == p[23]
JuMP.@objective m Max 1
JuMP.solve(m)
return [getvalue(x1), getvalue(x2), getvalue(x3), getvalue(x4), getvalue(x5), getvalue(x6), getvalue(x7), getvalue(x8), getvalue(x9), getvalue(x10)]
end
```
This problem is a biochemical reaction network, which yields a system of ODEs. I am trying to find a fixed point to these ODEs by solving these equations. These systems will always have at least one such point for X\>=0 (If it is well defined. I am very certain this system is properly defined). For this specific system I know that there exist exactly one fixed point (solution to my system). Depending on parameters this may either be a stable fixed point, or an unstable one around which the system has a limit cycle.
Now I try to solve for some specific parameter values:
```julia
using JuMP, NLopt, Ipopt
p = [3600, 18, 18, 3600, 3600, 3600, 1800, 3600, 18, 18, 18, 1800, 18, 36, 11, 180, 0.7, 0.4, 30, 0.2, 4, 4.5, 0.4]
my_problem(p,IpoptSolver())
```
This works fine, but when I start to changing the parameters things starts to behave funny:
```julia
p[14] = 14
my_problem(p,IpoptSolver())
```
Still ok. However:
```julia
p[14] = 13
my_problem(p,IpoptSolver())
```
gives
```julia
This is Ipopt version 3.12.10, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
Number of nonzeros in equality constraint Jacobian...: 45
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 51
Total number of variables............................: 10
variables with only lower bounds: 10
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 10
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 -1.0000000e+00 6.94e+03 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 -1.0000000e+00 6.90e+03 3.02e+01 -1.0 3.52e+02 - 1.57e-02 2.81e-03h 1
2 -1.0000000e+00 6.90e+03 6.42e+03 -1.0 6.43e+01 - 1.53e-01 1.54e-04h 1
3r-1.0000000e+00 6.90e+03 1.00e+03 2.0 0.00e+00 - 0.00e+00 4.68e-07R 3
4r-1.0000000e+00 2.06e+03 9.94e+02 2.0 9.14e+04 - 5.33e-03 1.02e-03f 1
5r-1.0000000e+00 1.41e+03 9.85e+02 0.6 3.54e+03 - 1.09e-02 4.77e-03f 1
6r-1.0000000e+00 1.26e+03 6.61e+02 0.6 4.47e+02 - 3.64e-01 1.27e-02f 1
7r-1.0000000e+00 1.03e+03 5.63e+02 0.6 6.19e+00 - 1.58e-01 1.06e-01f 1
8r-1.0000000e+00 8.14e+02 3.64e+02 0.6 3.53e+00 - 6.88e-01 1.45e-01f 1
9r-1.0000000e+00 2.17e+02 3.82e+02 0.6 1.40e+00 - 9.36e-01 5.50e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10r-1.0000000e+00 1.15e+02 1.15e+03 0.6 2.51e+00 - 3.76e-01 9.02e-01f 1
11r-1.0000000e+00 2.98e+01 3.25e+02 0.6 1.25e+00 - 4.49e-01 7.89e-01f 1
12r-1.0000000e+00 1.44e+01 1.95e+02 0.6 4.46e-01 - 9.17e-01 1.00e+00f 1
13r-1.0000000e+00 4.56e+00 2.43e+01 -0.1 9.01e-02 2.0 9.14e-01 9.37e-01f 1
14r-1.0000000e+00 2.16e+01 1.91e+02 -0.8 3.18e+01 - 6.63e-02 4.73e-02f 1
15r-1.0000000e+00 2.11e+01 1.59e+02 -0.8 2.57e+00 - 4.97e-01 1.14e-01f 1
16r-1.0000000e+00 7.81e+01 2.79e+02 -0.8 3.72e+00 - 6.60e-01 2.75e-01f 1
17r-1.0000000e+00 6.54e+01 1.39e+03 -0.8 1.85e+00 - 1.90e-02 1.44e-01f 1
18r-1.0000000e+00 3.31e+01 7.05e+02 -0.8 4.48e-01 - 5.28e-01 5.57e-01f 1
19r-1.0000000e+00 2.54e+01 3.86e+02 -0.8 9.60e-01 - 2.06e-01 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20r-1.0000000e+00 1.44e+00 2.06e+01 -0.8 2.43e-01 1.5 1.00e+00 9.54e-01h 1
21r-1.0000000e+00 2.22e+00 1.72e+01 -0.8 1.06e-01 1.0 1.00e+00 1.00e+00f 1
22r-1.0000000e+00 2.08e+00 2.19e+00 -0.8 6.81e-02 0.6 1.00e+00 1.00e+00h 1
23r-1.0000000e+00 9.90e+01 1.01e+02 -1.5 2.57e+00 - 3.28e-01 1.00e+00f 1
24r-1.0000000e+00 7.41e+01 6.02e+02 -1.5 6.90e+00 - 1.54e-01 1.00e+00f 1
25r-1.0000000e+00 4.15e+01 9.82e+01 -1.5 4.34e+00 - 7.49e-01 1.00e+00f 1
26r-1.0000000e+00 3.65e+01 8.57e+01 -1.5 6.46e+00 - 5.21e-01 1.03e-01f 1
27r-1.0000000e+00 2.45e+01 3.63e+01 -1.5 3.96e+00 - 1.00e+00 1.00e+00f 1
28r-1.0000000e+00 1.97e+00 9.21e-01 -1.5 2.55e-01 - 1.00e+00 1.00e+00h 1
29r-1.0000000e+00 1.39e+00 1.24e+00 -2.3 1.97e-02 - 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30r-1.0000000e+00 6.27e+01 5.55e+01 -2.3 8.25e+00 - 9.86e-01 6.78e-01f 1
31r-1.0000000e+00 4.46e+01 3.18e+01 -2.3 7.64e+00 - 1.00e+00 3.76e-01f 1
32r-1.0000000e+00 3.49e+01 1.47e+02 -2.3 2.71e-01 - 1.00e+00 2.16e-01h 1
33r-1.0000000e+00 7.22e-01 3.36e-01 -2.3 1.40e-01 - 1.00e+00 1.00e+00h 1
34r-1.0000000e+00 1.63e-03 2.77e-04 -2.3 8.60e-03 - 1.00e+00 1.00e+00h 1
35r-1.0000000e+00 5.73e-04 1.49e+00 -5.1 2.11e-02 - 1.00e+00 1.00e+00f 1
36r-1.0000000e+00 4.53e-04 6.17e+02 -5.1 5.84e-04 - 1.00e+00 2.06e-01h 1
37r-1.0000000e+00 3.86e-08 5.01e-05 -5.1 1.92e-05 - 1.00e+00 1.00e+00f 1
Number of Iterations....: 37
(scaled) (unscaled)
Objective...............: -1.0000000000000000e+00 -1.0000000000000000e+00
Dual infeasibility......: 0.0000000000000000e+00 0.0000000000000000e+00
Constraint violation....: 5.3435900520436833e-10 3.8573773508687736e-08
Complementarity.........: 0.0000000000000000e+00 0.0000000000000000e+00
Overall NLP error.......: 5.3435900520436833e-10 3.8573773508687736e-08
Number of objective function evaluations = 41
Number of objective gradient evaluations = 5
Number of equality constraint evaluations = 41
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 39
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 37
Total CPU secs in IPOPT (w/o function evaluations) = 0.020
Total CPU secs in NLP function evaluations = 0.008
EXIT: Optimal Solution Found.
This is Ipopt version 3.12.10, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
Number of nonzeros in equality constraint Jacobian...: 45
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 51
Total number of variables............................: 10
variables with only lower bounds: 10
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 10
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 -1.0000000e+00 6.94e+03 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 -1.0000000e+00 6.90e+03 3.05e+01 -1.0 3.55e+02 - 1.56e-02 2.79e-03h 1
2 -1.0000000e+00 6.90e+03 6.44e+03 -1.0 6.44e+01 - 1.53e-01 1.54e-04h 1
3 -1.0000000e+00 4.34e+03 1.02e+06 -1.0 4.14e+00 - 1.00e+00 1.78e-01h 1
4 -1.0000000e+00 4.34e+03 8.84e+08 -1.0 2.39e+00 - 1.00e+00 2.11e-05h 1
5 -1.0000000e+00 3.61e+03 5.92e+08 -1.0 5.47e+00 - 3.30e-01 9.87e-02h 1
6 -1.0000000e+00 3.61e+03 5.90e+06 -1.0 2.99e+00 - 9.90e-01 1.61e-03h 1
7 -1.0000000e+00 3.61e+03 6.77e+10 -1.0 1.72e+00 - 1.00e+00 1.56e-05h 1
8r-1.0000000e+00 3.61e+03 1.00e+03 1.9 0.00e+00 - 0.00e+00 3.94e-07R 2
9r-1.0000000e+00 2.03e+03 2.58e+03 1.9 6.85e+04 - 4.07e-03 1.01e-03f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10r-1.0000000e+00 1.04e+03 2.75e+03 1.2 3.93e+03 - 9.62e-03 3.83e-03f 1
11r-1.0000000e+00 5.89e+02 6.01e+02 1.2 5.13e+02 - 5.07e-01 1.25e-02f 1
12r-1.0000000e+00 3.65e+02 1.21e+03 1.2 5.31e+00 - 1.57e-01 9.84e-02f 1
13r-1.0000000e+00 3.13e+02 2.59e+03 1.2 3.64e+00 - 7.29e-01 1.39e-01f 1
14r-1.0000000e+00 2.53e+02 2.03e+03 1.2 2.15e+00 - 1.13e-01 1.97e-01f 1
15r-1.0000000e+00 8.62e+01 1.15e+03 1.2 1.58e+00 - 6.20e-01 8.04e-01f 1
16r-1.0000000e+00 5.30e+01 1.79e+02 1.2 9.74e-01 - 4.24e-01 4.67e-01f 1
17r-1.0000000e+00 2.08e+01 1.32e+02 0.5 6.72e-01 - 5.84e-01 7.12e-01f 1
18r-1.0000000e+00 8.45e+00 1.40e+02 -0.2 2.01e+00 - 4.52e-01 6.89e-01f 1
19r-1.0000000e+00 4.80e+01 1.14e+02 -0.2 2.95e+00 - 9.94e-01 7.37e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20r-1.0000000e+00 4.88e+00 3.69e+02 -0.2 3.12e-01 - 4.10e-01 1.00e+00f 1
21r-1.0000000e+00 4.11e+01 4.96e+02 -0.2 3.81e-01 - 4.92e-01 5.00e-01f 2
22r-1.0000000e+00 7.81e+01 2.80e+02 -0.2 3.35e+00 0.0 8.55e-02 9.13e-02f 2
23r-1.0000000e+00 9.50e+01 4.14e+02 -0.2 6.69e+00 0.4 2.18e-02 7.47e-02f 2
24r-1.0000000e+00 8.12e+01 5.84e+02 -0.2 4.74e-01 0.9 1.00e+00 1.79e-01h 1
25r-1.0000000e+00 5.21e+00 5.25e+01 -0.2 2.70e-01 - 9.24e-01 1.00e+00f 1
26r-1.0000000e+00 1.25e+01 2.60e+02 -0.2 6.17e-01 - 6.50e-01 1.00e+00f 1
27r-1.0000000e+00 4.54e+01 1.98e+01 -0.2 1.51e+00 - 1.00e+00 1.00e+00F 1
28r-1.0000000e+00 1.67e+00 1.39e+01 -0.9 4.50e-01 - 1.00e+00 1.00e+00f 1
29r-1.0000000e+00 9.21e+00 1.37e+01 -0.9 2.25e-01 0.4 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30r-1.0000000e+00 1.98e+00 1.78e+00 -0.9 7.01e-02 0.8 1.00e+00 1.00e+00h 1
31r-1.0000000e+00 1.15e+01 7.51e+00 -1.6 3.51e-01 0.3 9.58e-01 1.00e+00f 1
32r-1.0000000e+00 1.77e+01 7.85e+00 -1.6 3.10e-01 0.7 1.00e+00 1.00e+00h 1
33r-1.0000000e+00 2.42e+01 1.06e+02 -1.6 1.34e+00 - 2.87e-01 6.65e-02h 1
34r-1.0000000e+00 6.89e+01 3.19e+02 -1.6 9.69e-01 - 1.00e+00 1.95e-01f 1
35r-1.0000000e+00 1.46e+02 5.39e+02 -1.6 3.48e-01 - 4.30e-01 5.97e-01h 1
36r-1.0000000e+00 1.23e+02 1.31e+04 -1.6 1.71e+00 0.3 1.00e+00 2.07e-01f 1
37r-1.0000000e+00 7.08e+01 4.05e+03 -1.6 6.27e-01 - 2.63e-01 3.86e-01f 1
38r-1.0000000e+00 6.10e+01 1.70e+04 -1.6 6.18e-02 - 6.51e-01 1.38e-01h 1
39r-1.0000000e+00 2.00e+00 1.02e+03 -1.6 4.87e-02 - 9.28e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40r-1.0000000e+00 2.00e+00 1.33e+01 -1.6 8.44e-03 - 1.00e+00 1.00e+00h 1
41r-1.0000000e+00 2.00e+00 6.63e+00 -1.6 2.51e-03 - 1.00e+00 1.00e+00h 1
42r-1.0000000e+00 2.00e+00 5.25e-02 -1.6 1.30e-04 - 1.00e+00 1.00e+00h 1
43r-1.0000000e+00 2.00e+00 2.01e+00 -3.7 7.45e-03 - 9.99e-01 9.80e-01f 1
44r-1.0000000e+00 2.00e+00 2.23e+00 -3.7 1.44e-03 - 1.00e+00 1.00e+00h 1
45r-1.0000000e+00 2.00e+00 2.98e-03 -3.7 1.67e-05 - 1.00e+00 1.00e+00h 1
46r-1.0000000e+00 2.00e+00 5.99e+01 -5.5 8.41e-05 - 1.00e+00 9.05e-01f 1
47r-1.0000000e+00 2.00e+00 1.60e-02 -5.5 1.97e-06 - 1.00e+00 1.00e+00f 1
48r-1.0000000e+00 2.00e+00 9.35e-05 -5.5 8.94e-07 - 1.00e+00 1.00e+00f 1
49r-1.0000000e+00 2.00e+00 2.90e+01 -8.3 3.62e-07 - 1.00e+00 9.54e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50r-1.0000000e+00 2.00e+00 3.56e+00 -8.3 4.55e-06 - 8.74e-01 1.00e+00f 1
51r-1.0000000e+00 2.00e+00 2.58e+00 -8.3 2.21e-05 - 1.00e+00 4.35e-01h 2
52r-1.0000000e+00 2.00e+00 1.61e-05 -8.3 6.33e-06 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 52
(scaled) (unscaled)
Objective...............: -1.0000000000000000e+00 -1.0000000000000000e+00
Dual infeasibility......: 4.1973005524903554e-04 4.1973005524903554e-04
Constraint violation....: 3.7041736792201733e-02 1.9999986103698808e+00
Complementarity.........: 5.5547830709406970e-09 5.5547830709406970e-09
Overall NLP error.......: 3.7041736792201733e-02 1.9999986103698808e+00
Number of objective function evaluations = 71
Number of objective gradient evaluations = 10
Number of equality constraint evaluations = 71
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 55
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 53
Total CPU secs in IPOPT (w/o function evaluations) = 0.024
Total CPU secs in NLP function evaluations = 0.000
EXIT: Converged to a point of local infeasibility. Problem may be infeasible.
┌ Warning: Not solved to optimality, status: Infeasible
└ @ JuMP ~/.julia/packages/JuMP/PbnIJ/src/nlp.jl:1283
******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
Ipopt is released as open source code under the Eclipse Public License (EPL).
For more information visit http://projects.coin-or.org/Ipopt
******************************************************************************
This is Ipopt version 3.12.10, running with linear solver mumps.
NOTE: Other linear solvers might be more efficient (see Ipopt documentation).
Number of nonzeros in equality constraint Jacobian...: 45
Number of nonzeros in inequality constraint Jacobian.: 0
Number of nonzeros in Lagrangian Hessian.............: 51
Total number of variables............................: 10
variables with only lower bounds: 10
variables with lower and upper bounds: 0
variables with only upper bounds: 0
Total number of equality constraints.................: 10
Total number of inequality constraints...............: 0
inequality constraints with only lower bounds: 0
inequality constraints with lower and upper bounds: 0
inequality constraints with only upper bounds: 0
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
0 -1.0000000e+00 6.94e+03 1.00e+00 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0
1 -1.0000000e+00 6.90e+03 3.05e+01 -1.0 3.55e+02 - 1.56e-02 2.79e-03h 1
2 -1.0000000e+00 6.90e+03 6.44e+03 -1.0 6.44e+01 - 1.53e-01 1.54e-04h 1
3 -1.0000000e+00 4.34e+03 1.02e+06 -1.0 4.14e+00 - 1.00e+00 1.78e-01h 1
4 -1.0000000e+00 4.34e+03 8.84e+08 -1.0 2.39e+00 - 1.00e+00 2.11e-05h 1
5 -1.0000000e+00 3.61e+03 5.92e+08 -1.0 5.47e+00 - 3.30e-01 9.87e-02h 1
6 -1.0000000e+00 3.61e+03 5.90e+06 -1.0 2.99e+00 - 9.90e-01 1.61e-03h 1
7 -1.0000000e+00 3.61e+03 6.77e+10 -1.0 1.72e+00 - 1.00e+00 1.56e-05h 1
8r-1.0000000e+00 3.61e+03 1.00e+03 1.9 0.00e+00 - 0.00e+00 3.94e-07R 2
9r-1.0000000e+00 2.03e+03 2.58e+03 1.9 6.85e+04 - 4.07e-03 1.01e-03f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
10r-1.0000000e+00 1.04e+03 2.75e+03 1.2 3.93e+03 - 9.62e-03 3.83e-03f 1
11r-1.0000000e+00 5.89e+02 6.01e+02 1.2 5.13e+02 - 5.07e-01 1.25e-02f 1
12r-1.0000000e+00 3.65e+02 1.21e+03 1.2 5.31e+00 - 1.57e-01 9.84e-02f 1
13r-1.0000000e+00 3.13e+02 2.59e+03 1.2 3.64e+00 - 7.29e-01 1.39e-01f 1
14r-1.0000000e+00 2.53e+02 2.03e+03 1.2 2.15e+00 - 1.13e-01 1.97e-01f 1
15r-1.0000000e+00 8.62e+01 1.15e+03 1.2 1.58e+00 - 6.20e-01 8.04e-01f 1
16r-1.0000000e+00 5.30e+01 1.79e+02 1.2 9.74e-01 - 4.24e-01 4.67e-01f 1
17r-1.0000000e+00 2.08e+01 1.32e+02 0.5 6.72e-01 - 5.84e-01 7.12e-01f 1
18r-1.0000000e+00 8.45e+00 1.40e+02 -0.2 2.01e+00 - 4.52e-01 6.89e-01f 1
19r-1.0000000e+00 4.80e+01 1.14e+02 -0.2 2.95e+00 - 9.94e-01 7.37e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
20r-1.0000000e+00 4.88e+00 3.69e+02 -0.2 3.12e-01 - 4.10e-01 1.00e+00f 1
21r-1.0000000e+00 4.11e+01 4.96e+02 -0.2 3.81e-01 - 4.92e-01 5.00e-01f 2
22r-1.0000000e+00 7.81e+01 2.80e+02 -0.2 3.35e+00 0.0 8.55e-02 9.13e-02f 2
23r-1.0000000e+00 9.50e+01 4.14e+02 -0.2 6.69e+00 0.4 2.18e-02 7.47e-02f 2
24r-1.0000000e+00 8.12e+01 5.84e+02 -0.2 4.74e-01 0.9 1.00e+00 1.79e-01h 1
25r-1.0000000e+00 5.21e+00 5.25e+01 -0.2 2.70e-01 - 9.24e-01 1.00e+00f 1
26r-1.0000000e+00 1.25e+01 2.60e+02 -0.2 6.17e-01 - 6.50e-01 1.00e+00f 1
27r-1.0000000e+00 4.54e+01 1.98e+01 -0.2 1.51e+00 - 1.00e+00 1.00e+00F 1
28r-1.0000000e+00 1.67e+00 1.39e+01 -0.9 4.50e-01 - 1.00e+00 1.00e+00f 1
29r-1.0000000e+00 9.21e+00 1.37e+01 -0.9 2.25e-01 0.4 1.00e+00 1.00e+00f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
30r-1.0000000e+00 1.98e+00 1.78e+00 -0.9 7.01e-02 0.8 1.00e+00 1.00e+00h 1
31r-1.0000000e+00 1.15e+01 7.51e+00 -1.6 3.51e-01 0.3 9.58e-01 1.00e+00f 1
32r-1.0000000e+00 1.77e+01 7.85e+00 -1.6 3.10e-01 0.7 1.00e+00 1.00e+00h 1
33r-1.0000000e+00 2.42e+01 1.06e+02 -1.6 1.34e+00 - 2.87e-01 6.65e-02h 1
34r-1.0000000e+00 6.89e+01 3.19e+02 -1.6 9.69e-01 - 1.00e+00 1.95e-01f 1
35r-1.0000000e+00 1.46e+02 5.39e+02 -1.6 3.48e-01 - 4.30e-01 5.97e-01h 1
36r-1.0000000e+00 1.23e+02 1.31e+04 -1.6 1.71e+00 0.3 1.00e+00 2.07e-01f 1
37r-1.0000000e+00 7.08e+01 4.05e+03 -1.6 6.27e-01 - 2.63e-01 3.86e-01f 1
38r-1.0000000e+00 6.10e+01 1.70e+04 -1.6 6.18e-02 - 6.51e-01 1.38e-01h 1
39r-1.0000000e+00 2.00e+00 1.02e+03 -1.6 4.87e-02 - 9.28e-01 1.00e+00h 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
40r-1.0000000e+00 2.00e+00 1.33e+01 -1.6 8.44e-03 - 1.00e+00 1.00e+00h 1
41r-1.0000000e+00 2.00e+00 6.63e+00 -1.6 2.51e-03 - 1.00e+00 1.00e+00h 1
42r-1.0000000e+00 2.00e+00 5.25e-02 -1.6 1.30e-04 - 1.00e+00 1.00e+00h 1
43r-1.0000000e+00 2.00e+00 2.01e+00 -3.7 7.45e-03 - 9.99e-01 9.80e-01f 1
44r-1.0000000e+00 2.00e+00 2.23e+00 -3.7 1.44e-03 - 1.00e+00 1.00e+00h 1
45r-1.0000000e+00 2.00e+00 2.98e-03 -3.7 1.67e-05 - 1.00e+00 1.00e+00h 1
46r-1.0000000e+00 2.00e+00 5.99e+01 -5.5 8.41e-05 - 1.00e+00 9.05e-01f 1
47r-1.0000000e+00 2.00e+00 1.60e-02 -5.5 1.97e-06 - 1.00e+00 1.00e+00f 1
48r-1.0000000e+00 2.00e+00 9.35e-05 -5.5 8.94e-07 - 1.00e+00 1.00e+00f 1
49r-1.0000000e+00 2.00e+00 2.90e+01 -8.3 3.62e-07 - 1.00e+00 9.54e-01f 1
iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls
50r-1.0000000e+00 2.00e+00 3.56e+00 -8.3 4.55e-06 - 8.74e-01 1.00e+00f 1
51r-1.0000000e+00 2.00e+00 2.58e+00 -8.3 2.21e-05 - 1.00e+00 4.35e-01h 2
52r-1.0000000e+00 2.00e+00 1.61e-05 -8.3 6.33e-06 - 1.00e+00 1.00e+00h 1
Number of Iterations....: 52
(scaled) (unscaled)
Objective...............: -1.0000000000000000e+00 -1.0000000000000000e+00
Dual infeasibility......: 4.1973005524903554e-04 4.1973005524903554e-04
Constraint violation....: 3.7041736792201733e-02 1.9999986103698808e+00
Complementarity.........: 5.5547830709406970e-09 5.5547830709406970e-09
Overall NLP error.......: 3.7041736792201733e-02 1.9999986103698808e+00
Number of objective function evaluations = 71
Number of objective gradient evaluations = 10
Number of equality constraint evaluations = 71
Number of inequality constraint evaluations = 0
Number of equality constraint Jacobian evaluations = 55
Number of inequality constraint Jacobian evaluations = 0
Number of Lagrangian Hessian evaluations = 53
Total CPU secs in IPOPT (w/o function evaluations) = 0.172
Total CPU secs in NLP function evaluations = 0.056
EXIT: Converged to a point of local infeasibility. Problem may be infeasible.
┌ Warning: Not solved to optimality, status: Infeasible
└ @ JuMP ~/.julia/packages/JuMP/PbnIJ/src/nlp.jl:1283
```
(One run of the function yields many such sections, each ending with `EXIT: Converged to a point of local infeasibility. Problem may be infeasible.`)
However, it still have a solution (it is `[0.102, 0.416, 1.79, 0.0269, 5.84, 22.0, 0.112, 7.87, 0.399, 0.000999]`). I can even confirm that this is indeed a solution using this function:
```julia
function my_problem_func(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,p)
output = fill(0.,10)
output[1] = -1 * p[17] * x1 + -2 * p[1] * (x1 ^ 2 / 2) + 2 * p[2] * x2 + p[21] * p[18] * ((1 + p[19] * x7) / (p[20] + x7))everything
output[2] = -1 * p[17] * x2 + p[1] * (x1 ^ 2 / 2) + -1 * p[2] * x2 + -1 * p[4] * x2 * x4 + p[9] * x3 + p[14] * x3 + -1 * p[6] * x2 * x7 + p[11] * x8
output[3] = -1 * p[17] * x3 + p[4] * x2 * x4 + -1 * p[9] * x3 + -1 * p[5] * x3 * x4 + p[10] * x5 + -1 * p[14] * x3 + p[15] * x5 + p[7] * x8 * x4 + -1 * p[12] * x3 * x7
output[4] = -1 * p[17] * x4 + -1 * p[4] * x2 * x4 + p[9] * x3 + -1 * p[5] * x3 * x4 + p[10] * x5 + -1 * p[7] * x8 * x4 + p[12] * x3 * x7 + p[16] * x9 + p[22] * p[18] * ((1 + p[19] * x7) / (p[20] + x7))
output[5] = -1 * p[17] * x5 + p[5] * x3 * x4 + -1 * p[10] * x5 + -1 * p[15] * x5
output[6] = -1 * p[17] * x6 + p[14] * x3 + p[15] * x5 + -1 * p[8] * x6 * x10 + p[13] * x9
output[7] = -1 * p[17] * x7 + -1 * p[6] * x2 * x7 + p[11] * x8 + p[7] * x8 * x4 + -1 * p[12] * x3 * x7 + p[18] * ((1 + p[19] * x7) / (p[20] + x7))
output[8] = -1 * p[17] * x8 + p[6] * x2 * x7 + -1 * p[11] * x8 + -1 * p[7] * x8 * x4 + p[12] * x3 * x7
output[9] = -1 * p[17] * x9 + p[8] * x6 * x10 + -1 * p[13] * x9 + -1 * p[16] * x9
output[10] = x10 + x9 - p[23]
output
end
```
If I input this solution as starting point I can solve the problem
```julia
my_problem(p,IpoptSolver(),xstart=[0.102, 0.416, 1.79, 0.0269, 5.84, 22.0, 0.112, 7.87, 0.399, 0.000999])
```
however I want to solve for many parameters, and several similar models, so it is hard to start making guesses manually.
I have tried altering all kinds of parameters in `Ipopt`, but to no help (setting `tol` really high lead to a problem in the “restoration”, but otherwise I just got the same error.
I have tried all kinds of solvers I could find, the only other one which seems to attempt on the problem is `NLopt`s `LD_SLSQP`. However it simply fails to return a correct answer:
```julia
sol = my_problem(p,NLoptSolver(algorithm=:LD_SLSQP))
```
returns `sol = [0.337, 0.00, 5.59, 0.000831, 3.36, 0.00172, 0.000185, 2.87, 0.607, 0.00]`, which is simply wrong.
I am no expert with this type of work, but I have really been trying all I know, with no success. I have also tried `DifferentialEquations.jl`, however simulating the system does not work out when there is an oscillation, and using `SteadyStateProblem`s may find negative solutions. I have also been using homotopy continuations, which can solve single instances, but it is way to slow when I want to solve 100 or even 1000 related systems.
In principle it feels like quite a straightforward problem, 10 variables, mostly 2nd order polynomials and a single fraction. Only one solution and it is guaranteed to be there. If anyone with more experience could help me that would be immensely useful. Thanks.
---
### Author:  [@odow](https://discourse.julialang.org/u/odow)
#### Post date: [January 24, 2019, 6:51pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/2 "2019-01-24T18:51:18Z")
Nonlinear problems are hard to solve, so it’s not surprising that you run into issues.
Even though you don’t care about the objective, try providing an arbitrary one (like `@objective(m, Min, sum(x))`. This might help.
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [January 24, 2019, 7:01pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/3 "2019-01-24T19:01:04Z")
I would be interested to try `IntervalRootFinding.jl` on this problem. Do you have it written in a more mathematical format?
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [January 24, 2019, 8:55pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/4 "2019-01-24T20:55:22Z")
> [@Torkel](#):
>
> [0.102, 0.416, 1.79, 0.0269, 5.84, 22.0, 0.112, 7.87, 0.399, 0.000999]
What is `everything` in the third line of `my_problem_func` at the end?
The solution you give does not seem to give 0 when I put it in `my_problem_func`?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 24, 2019, 8:55pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/5 "2019-01-24T20:55:31Z")
> [@Torkel](#):
>
> I am no expert with this type of work, but I have really been trying all I know, with no success. I have also tried `DifferentialEquations.jl` , however simulating the system does not work out when there is an oscillation
What integrator and tolerance? Can you share the code?
You may want be encoding those domain constraints through transformations of the equation. That would make it a lot easier on the methods.
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [January 24, 2019, 8:56pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/6 "2019-01-24T20:56:24Z")
Could you please give the exact code for defining the vector of `p`s and of `x`s that should be a solution?
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [January 24, 2019, 10:34pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/7 "2019-01-24T22:34:19Z")
Note that you can also try [https://www.juliahomotopycontinuation.org/](https://www.juliahomotopycontinuation.org/) for polynomial equations.
---
### Author:  [@dpsanders](https://discourse.julialang.org/u/dpsanders)
#### Post date: [January 24, 2019, 10:35pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/8 "2019-01-24T22:35:15Z")
And you can try a “simple” Newton–Raphson method.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [January 25, 2019, 3:27am UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/9 "2019-01-25T03:27:04Z")
> [@Torkel](#):
>
> I initially tried `NLsolve.jl` , however, it does not allow for specifying variables as non-negative.
You can always just plug `x^2` in for a variable that you want to be non-negative, to allow the solver to use an unconstrained `x`, no?
---
### Author:  [@ccoffrin](https://discourse.julialang.org/u/ccoffrin)
#### Post date: [January 25, 2019, 3:38am UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/10 "2019-01-25T03:38:25Z")
Have you considered a global optimization solver? For example, Couenne, Alpine or Baron? These are designed for problems where local solvers are problematic.
---
### Author:  [@frapac](https://discourse.julialang.org/u/frapac)
#### Post date: [January 25, 2019, 10:17am UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/11 "2019-01-25T10:17:54Z")
You could also use the multistart option of KNITRO:
[https://www.artelys.com/tools/knitro\_doc/2\_userGuide/multistart.html](https://www.artelys.com/tools/knitro_doc/2_userGuide/multistart.html)
Replacing Ipopt by KNITRO and enabling the multistart option
```julia
my_problem(p, KnitroSolver(ms_enable=1)
```
yields the solution
```julia
[0.101743, 0.415984, 1.79022, 0.026919, 5.84132, 22.0398, 0.112034, 7.87432, 0.399001, 0.000999219]
```
which seems close to the solution you are looking for.
If you do not have KNITRO available, you can request a free trial license.
---
### Author:  [@saschatimme](https://discourse.julialang.org/u/saschatimme)
#### Post date: [January 25, 2019, 1:43pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/12 "2019-01-25T13:43:12Z")
@Torkel (as you already know) you can solve these kinds of problems with [HomotopyContinuation.jl](https://www.juliahomotopycontinuation.org).
But if you want to solve the system for many parameters the fastest (and numerically best) method is to first solve a “generic” system and then using a parameter homotopy to your specific parameters.
I constructed an IJulia notebook to demonstrate this:
[https://gist.github.com/saschatimme/38fc2bbe641cfd7f3985c0169f6efed6](https://gist.github.com/saschatimme/38fc2bbe641cfd7f3985c0169f6efed6)
Note that the generic solve works fastest on a soon to be released version due to some undesired behaviour in the current release.
I would love to turn this also in a blog post on [https://www.juliahomotopycontinuation.org/blog](https://www.juliahomotopycontinuation.org/blog) as soon as we created a new release. It would be cool to have some help with the bio background for that 🙂
---
### Author:  [@Torkel](https://discourse.julialang.org/u/Torkel)
#### Post date: [January 25, 2019, 3:32pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/13 "2019-01-25T15:32:53Z")
This is _really_ useful, thanks everyone. I’ll try to reply to all the individual points.
---
### Author:  [@Torkel](https://discourse.julialang.org/u/Torkel)
#### Post date: [January 25, 2019, 3:33pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/14 "2019-01-25T15:33:47Z")
Didn’t help in this instance, but I will keep this in mind when doing further solving using JuMP.
---
### Author:  [@Torkel](https://discourse.julialang.org/u/Torkel)
#### Post date: [January 25, 2019, 3:54pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/15 "2019-01-25T15:54:39Z")
The objective is to find solutions to this system:
However, although this is 10 equations it is slight degenerate (correct word?), and these actually exists a one dimensional line of solutions. So I have removed the last equation and replaced it with:
_vPp+phos = pTot_
The parameter values are then in the following ones
[kBw kDw kD kB1 kB2 kB3 kB4 kB5 kD1 kD2 kD3 kD4 kD5 kK1 kK2 kP kDeg v0 F K λW λV pTot]
[3600, 18, 18, 3600, 3600, 3600, 1800, 3600, 18, 18, 18, 1800, 18, 13, 11, 180, 0.7, 0.4, 30, 0.2, 4, 4.5, 0.4]
The `everything ` in `my_problem_func` is a simple typo. Your right the solution is not exact, my fault. I only included the first few decimals, which considering some of the large parameters amplifies and make it seem non zero. I will add a few more decimals to make it more correct.
I used HomotopyContinuation to obtain the real solution (and simulations using DifferentialEquations), my problem was that when I wanted to do a lot of solving that method was to slow, so I figured I try something else (with the disadvantage that I could not find all solutions, but looking at your interval root finding that one might be able to do that).
---
### Author:  [@saschatimme](https://discourse.julialang.org/u/saschatimme)
#### Post date: [January 25, 2019, 4:25pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/16 "2019-01-25T16:25:53Z")
> [@Torkel](#):
>
> used HomotopyContinuation to obtain the real solution (and simulations using DifferentialEquations), my problem was that when I wanted to do a lot of solving that method was to slow,
Is the approach I described above fast enough for you? Otherwise it is also possible (with a little bit code) to speed up things even more by 1) keeping the internal data structures around for reuse and 2) stopping the computation as soon as you find the desired solution.
---
### Author:  [@Torkel](https://discourse.julialang.org/u/Torkel)
#### Post date: [January 25, 2019, 4:34pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/17 "2019-01-25T16:34:21Z")
Haven’t tried it yet (department free beer event happened…), but it looks very interesting. I am basically doing a bifurcation analysis, and I would then be able to solve the problem for a specific parameter, and then use that solution to create a homotopy with yields the solution as that parameter changes? If so that would probably be the best way.
(HomotopyContinuation is great because it finds all states without me having to give extra input, in this specific case I knew there were only one solution so figured I could speed it up by using another method, but this way I might be able to keep to Homotopy Continuation)
---
### Author:  [@Torkel](https://discourse.julialang.org/u/Torkel)
#### Post date: [January 25, 2019, 4:35pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/18 "2019-01-25T16:35:05Z")
This is a cool idea, yes, it should be possible to use a trick like this to solve it?
---
### Author:  [@Torkel](https://discourse.julialang.org/u/Torkel)
#### Post date: [January 25, 2019, 4:37pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/19 "2019-01-25T16:37:27Z")
You mean one could use a transformation with arctan to ensure that non-negative values are not considered (I think remember you suggested using arctan for a very different problem).
I’ll dig up the code I used.
---
### Author:  [@Torkel](https://discourse.julialang.org/u/Torkel)
#### Post date: [January 25, 2019, 4:38pm UTC](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051/20 "2019-01-25T16:38:44Z")
Looked at Couenne, but couldn’t get it to work, will make another attempt. had some issue with the installation of Alpine, but I hope I should be able to get it to work.
[Next page](https://discourse.julialang.org/t/solvers-fail-on-nonlinear-problem-which-has-solutions/20051.md?page=2)