Sorry to be late to the party. Unfortunately, you are not using L-BFGS-B in Optim, you’re using a much less efficient barrier method. In NLSolvers.jl I have a projected newton method, but since you’re already using finite differences for your gradients, this might be a very bad idea. It shouldn’t be too bad to adjust it to work with a BFGS-approximation though.
So I think this is three things: you’re using a suboptimal solver (“my fault”), you’re using finite differences (“your problem’s fault”), and your problem is not differentiable (it’s a simulation). I think the last two are actually worse for the barrier method than for l-bfgs-b, but let that be a conjecture.
If I look at a result I get
julia> result.minimizer
30-element Vector{Float64}:
1.62658413314516e-10
0.9999999998063619
0.999999999921528
0.9999999999434931
0.9999999999526344
0.9999999999586152
0.9999999999615918
0.9999999999635981
0.9999999999650679
0.999999999966145
0.9999999999668981
0.9999999999676399
0.9999999999681027
0.999999999968442
0.9999999999688338
0.9999999999689775
0.9999999999690062
0.9999999999691217
0.9999999999690876
0.9999999999689795
0.9999999999687171
0.9999999999682864
0.9999999999677559
0.9999999999668858
0.9999999999658892
0.9999999999643274
0.9999999999623169
0.9999999999596437
0.9999999999558721
0.999999999950147
this is not good for the barrier method. My guess is that the solution is at the boundary (one is 0 the rest are 1) and L-BFGS-B quickly lands exactly on a face of the constraint “box”. I 100% intend to implement L-BFGS-B in NLSolvers.jl but I think I might do projected BFGS first.