(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).

```
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:

```
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:

```
p[14] = 14
my_problem(p,IpoptSolver())
```

Still ok. However:

```
p[14] = 13
my_problem(p,IpoptSolver())
```

gives

```
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:

```
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

```
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:

```
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.