Solvers fail on nonlinear problem (which has solutions)

jump

#21

This sounds promising, I run into a problem with KNITRO installation, but got help which should fix it. I hope to test this tomorrow.


#22

Just a second thanks to everyone, and a note why I love the community around Julia. I go to bed and the next morning I have loads of really good ideas on how to fix my problem, looking forward to looking deeper into this the following few days.


#23

arctan was for removing a discontinuity. x^2 or exp transformations are for positivity constraints. In fact, we may want to bake that into DiffEqBiological somehow.


#24

The problem with this is that it might create a lot of “extra” roots that have some x_i < 0.


#25

You could try https://github.com/rdeits/CouenneNL.jl, which just downloads Couenne binaries instead of building them from source and hides the AmplNLWriter layer. Note that Couenne is not the fastest global solver, but it’s free.


#26

Indeed, doing a simple Newton method, I always end up finding roots that have some component negative:


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]

using ForwardDiff

function newton(f, x, n=1000)

    for i in 1:n

        J = ForwardDiff.jacobian(f, x)

        δ = J \ (-f(x))

        x += δ
    end
    
    return x
end

x0 = 100 * rand(10)
 
root = newton(x -> f(x, p), x0)

@show maximum(abs.(f(root, p)))
root

gives

maximum(abs.(f(root, p))) = 4.547473508864641e-13
10-element Array{Float64,1}:
  0.11482204904345317   
  0.00032484595992817946
  0.0406348879740191    
  6.852413274079044     
 33.75115705986466      
 -0.06847078576137997   
 16.87400939428231      
  0.1015090297557086    
  2.062941412150575     
 -1.6629414121505748    

#27

OK, just had to run it a few times from random initial conditions until I happened upon a non-negative solution.


#28

How do you prove that there is a unique non-negative root?


#29

This worked well, for those problems that Ipopt also can solve it is slower (Ipopt~0.025s and CouenneNL~0.25s). However it do work well, thanks.


#30

Yes, using your approach is much better. This should probably be fast enough (or rather it is, I will see what happens when I increase complexity of the problem, but I think it should work well, especially if you are about to release an improved version).

I’d be more than happy to help you with the blog post.


#31

I don’t have any proof really, but I know that the system neither can cross into negative values, nor approach infinity (due to the nature of the model and both of these alternatives are biologically weird). Also, the system should be quite well behaved and I do not expect any strange cases. There will be at least one root.

In principle there could be more than one root, however, after playing around with this model for a while I am quite certain that this is not the case (and for this specific application I am satisfied with finding a single solution).