[R] optim() not finding optimal values

[dave fournier]

 If you are going to make this program available for general
  use you want to take every precaution to make it bulletproof.

  This is a fairly informative data set. The model will undoubtedly
  be used on far less informative data.  While the model looks
  pretty simple it is  very challenging from a numerical point of view.
  I took a moment to code it up in AD Model Builder. The true minimum is
  1619.480495 So I think Ravi has finally arrived pretty close to the
answer.
  One way of judging the difficulty of a model is to look at the
  eigenvalues of the Hessian at the minimum. They are

       3.122884668e-09 1.410866202e-08  1866282.520 1.330233652e+13

  so the condition number is around 1.e+21. One begins to see why these
  models are challenging.  The model as formulated represents the state
  of the art in fisheries models circa 1985.
  A lot of progress has been made since that time.
  Using B_t for the biomass and C_t for the catch the equation
  in the code is

            B_{t+1} = B_t + r *B_t*(1-B_t/K) -C_t  (1)
  First  notice that
  for (1) to make sense the following conditions must be satisfied

       B_t > 0 for all t
       r > 0
       K>0
  Strictly speaking it is not necessary that B_t<=K but if B_t>K and r
  is large then B_{t+1} could be <0.  So formulation (1) gives
  Murphys law a good chance.  How to fix it. Notice that (1) is really
  a rough approximation to the solution of a differential equation

      B'(t) =  r *B(t)*(1-B(t)/K) -C  (2)

  where in (2) C is a constant catch rate.  To fix (1) we use
  a semi-implicit differencing scheme. Because it is useful to
  allow smaller step sizes than one we denote them by d.

       B_{t+d} = B_t + d* r *B_t*(1-B_{t+d}/K) -d*C_t*B_{t+d}/B_t  (1)

  The idea is that the quantity  1-x with x>0 will be replaced by
  1/(1+x).  Expanding 2 and solving for B_{t+d} yields
      B_{t+d} = (1+d*r) B_t / (1+d*r*B_t/K +d*C_t/B_t)  (3)

   So long as r>0, K>0 C_t>0 then starting from an initial value
   B_0 > 0 ensures that B_t> 0 for all t>0.  We can let
   d=1/nsteps where nsteps is the number of steps in the
   approximate integration for each year
   which can be increased until the solution is judged to be close
   enough to the exact solution from (2)

   Notice that in (3) as C_t --> infinity  B_{t+d} --> 0
   So that you can never catch more fish than you have.

   I coded up this version of the model in AD Model Builder and
   fit it to the data. It is now much more resistant to bad
   starting values for the parameters etc.

   If anyone wants the tpl file for the model in ADMB they can
   contact me off list.