Solvers fail on nonlinear problem (which has solutions)

Specific Domains Optimization (Mathematical)
1

(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 NLopts 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 SteadyStateProblems 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.

2

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.

3

I would be interested to try IntervalRootFinding.jl on this problem. Do you have it written in a more mathematical format?

4

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?

5

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.

6

Could you please give the exact code for defining the vector of ps and of xs that should be a solution?

7

Note that you can also try https://www.juliahomotopycontinuation.org/ for polynomial equations.

8

And you can try a “simple” Newton–Raphson method.

9

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?

1 Like
10

Have you considered a global optimization solver? For example, Couenne, Alpine or Baron? These are designed for problems where local solvers are problematic.

11

You could also use the multistart option of KNITRO:
https://www.artelys.com/tools/knitro_doc/2_userGuide/multistart.html

Replacing Ipopt by KNITRO and enabling the multistart option

my_problem(p, KnitroSolver(ms_enable=1)

yields the solution

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

12

@Torkel (as you already know) you can solve these kinds of problems with HomotopyContinuation.jl.

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

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 as soon as we created a new release. It would be cool to have some help with the bio background for that :slight_smile:

2 Likes
13

This is really useful, thanks everyone. I’ll try to reply to all the individual points.

14

Didn’t help in this instance, but I will keep this in mind when doing further solving using JuMP.

15

The objective is to find solutions to this system:

image
image1128Ă—456 34.6 KB

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

16

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.

17

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)

18

This is a cool idea, yes, it should be possible to use a trick like this to solve it?

19

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.

20

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.