Q:

Maximising solutions of NDSolve

I have the following system of differential equations
systeqns[gn_, gp_, sn_, sp_, K_, a_, 
  b_, n0_, p0_, tg_] := {

  n'[t] == 
   If[t <= tg, gn, sn]*n[t]*(1 - (n[t] + p[t])/K) - 
    a*n[t] + b*p[t],

  p'[t] == 
   If[t <= tg, gp, sp]*p[t]*(1 - (n[t] + p[t])/K) + 
    a*n[t] - b*p[t],

  tot'[t] == n'[t] + p'[t],

  p[0] == p0,
  n[0] == n0,
  tot[0] == p0 + n0}

sol[gn_, gp_, sn_, sp_, K_, a_, b_, 
   n0_, p0_, tg_, ts_] := 
  First[Evaluate[
    NDSolve[systeqns[gn, gp, sn, sp, 
      K, a, b, n0, p0, tg], {n, p, tot}, {t, 0, tg + ts}, 
     Method -> {"StiffnessSwitching"}, AccuracyGoal -> 1000, 
     WorkingPrecision -> 30, MaxSteps -> Infinity]]]; 

From which I calculate
avgR[gn_, gp_, sn_, sp_, K_, a_, b_, 
  n0_, p0_, tg_, ts_] := 
 Log[(tot[tg + ts]/tot[0]) /. 
   sol[gn, gp, sn, sp, K, a, b, n0, 
    p0, tg, ts]]

I define some example values
tgex = 5; tsex = 5;gnex = 2;gpex = 2/10;snex = -2;spex = -2/10;bex = 1/10000;n0ex = 1;p0ex = 1/1000;Kex = 10^9;

I am then interested in determining the value of a that maximizes avgR. Graphically, this would be around a=0.2:
vals = Table[{10^pow, 
    avgR[gnex, gpex, snex, spex, Kex, 
     10^pow, bex, n0ex, p0ex, tgex, tsex]}, {pow, -4, 0, 1/10}];
ListLogLinearPlot[vals, PlotRange -> All, Frame -> True, 
 FrameLabel -> {"Switching rate a", "Average population growth rate"},
  PlotStyle -> {Thick}, PlotTheme -> {"Classic", "LargeLabels"}, 
 Joined -> True]  

But trying to calculate this value using NMaximize returns an error:
NMaximize[{avgR[gnex, gpex, snex, spex, Kex, 
     a, bex, n0ex, p0ex, tgex, tsex], 
  a > 10^-4 && a < 1}, {a, 1/10}, Method -> "DifferentialEvolution"]
NDSolve::ndnum: Encountered non-numerical value for a derivative at t == 0.  >>

Any thoughts what I am doing wrong? 
Interpolating the values first does work btw, but that is no good for me, since I would eventually like to maximise both a and b, and this would require me calculating too big a grid :
f = Interpolation[vals];
Maximize[{f[a], 10^-4 < a < 1}, a]
{6.507434893938602167476249656, {a -> 
   0.20613498164678882000313750361}}

So I am looking for a solution to get NMaximize working...

A:

You have to suppress symbolic evaluations by replacing your definition of avgR with, e.g.
avgR[gn_, gp_, sn_, sp_, K_, a_?NumericQ, b_, n0_, p0_, tg_, ts_] := 
  Log[(tot[tg + ts]/tot[0]) /. sol[gn, gp, sn, sp, K, a, b, n0, p0, tg, ts]]

In NMaximize the 1/10 has to be removed
NMaximize[{avgR[gnex, gpex, snex, spex, Kex, a, bex, n0ex, p0ex, tgex,
tsex], a > 10^-4 && a < 1}, {a}, Method -> "DifferentialEvolution"]

{6.50728, {a -> 0.205164}

If you want to use a starting value and you are only searching for a local maximum you can use
FindMaximum[{avgR[gnex, gpex, snex, spex, Kex, a, bex, n0ex, p0ex, tgex, tsex], 
  a > 10^-4 && a < 1}, {{a, 1/10}}]

{6.50728, {a -> 0.205164}

Using FindMaximum is also much faster than using NMaximize with Method -> "DifferentialEvolution" (6.6 sec vs. 170.1 sec on my PC).  NMaximize with Method -> "NelderMead" or Method -> "SimulatedAnnealing" would be faster (13.4 sec and 14.4 sec, respectively) than "DifferentialEvolution". 
You can get rid of the messages by setting the WorkingPrecision >= the value you set for NDSolve, e.g.
NMaximize[{avgR[gnex, gpex, snex, spex, Kex, a, bex, n0ex, p0ex, tgex, tsex], 
  a > 10^-4 && a < 1}, {a}, Method -> "NelderMead", WorkingPrecision -> 30]

