NMaximize not always returning solution

Question

This question is related to a previous question and answer here:

NMinimize or NMaximize for a dynamic problem

but tackles a different problem.

I am trying to solve a simple linear quadratic optimal depletion problem.
Choose $q(t)$ from $t=0$ to $ T$ so as to maximize the objective when $R(t+1)=R(t)-q(t).$

My code is here:

 T = 100; p = 1.; δ = 0.05127; ρ = 1./(1. + δ); c = 50.; R0 = 1.;
obj[qlist_ /; VectorQ[qlist, NumericQ]] :=
 Module[{Rlist, obj},
  Rlist = FoldList[(#1 - #2) &, R0, qlist];

If[Min[Rlist] < 0, Return[-1000]]; (* Constraint to keep R(t)>=0 for all t *)

(* This is the objective function *)
Sum[ρ^t*(p*qlist[[t + 1]] - 
      0.5 c ((qlist[[t + 1]])^2)/Rlist[[t + 1]]), {t, 0, T - 1}]
  ]

choicevar = Table[Unique[q], {i, 0, T - 1}];
sol = NMaximize[Prepend[Thread[0 <= choicevar], obj[choicevar]], 
   choicevar];
sol[[1]]

I have an analytical solution to this class of problems and this code works perfectly when $R0=100$.
But for $R0=1$ as above, it fails.
The problem also wont solve if I set $T=250$ and $R0=100$.

The failures are of two kinds.

1) failures to converge to desired accuracy, or
2) violation of the constraint $(R(t)=>0)$, Returning a solution -1000.

The problems that will not solve all seem to have solutions (in the later periods) where q gets close to zero (the real solution may be or order 10^-3, for instance) so I suspect the issue has to do with accuracy or precision but I tried variations and allowing up to 10000 iterations and it does not work.

Can anyone please help me unlock the power of NMaximize?
At this point, Excel Solver seems to be a better solution but I refuse to believe that! Thanks.

Answer 1

Updated)

In general,we need to carefully set the feasible region of control.

I mean,

sol = NMaximize[Prepend[Thread[0 <= choicevar <= 0.020], obj[choicevar]],choicevar];

so just check the result.

T = 100; p = 1.; \[Delta] = 0.05127; \[Rho] = 
 1./(1. + \[Delta]); c = 50.; R0 = 1.;
obj[qlist_ /; VectorQ[qlist, NumericQ]] := 
 Module[{Rlist, obj}, Rlist = FoldList[(#1 - #2) &, R0, qlist];
  If[Min[Rlist] < 0, 
   Return[-1000]];(*Constraint to keep R\[GreaterEqual]0*)(*This is \
the objective function*)
  Sum[\[Rho]^
     t*(p*qlist[[t + 1]] - 
      0.5 c ((qlist[[t + 1]])^2)/Rlist[[t + 1]]), {t, 0, T - 1}]]

choicevar = Table[Unique[q], {i, 0, T - 1}];
sol = NMaximize[
   Prepend[Thread[0 <= choicevar <= 0.02], obj[choicevar]], choicevar];
sol[[1]]

=>

0.150353

ListLinePlot@Values@Last@sol

If you can not solve the problem, please try to make sure that the conditions you imposed. Even if not in the text, you often need additional one.

---

Perhaps I think the function approximation using the pseudospectral method is the most efficient, but I have never done it...

so,I used heuristics.

First,Let's sampling and guess the optimal solution.

Table[try = RandomVariate[UniformDistribution[{0, 0.014}], {T, 10}];
  data = {try[[All, #]], obj[try[[All, #]]]} & /@ Range[1, 10];
  SortBy[data, -#[[2]] &][[1]], {1000}];
Mean@sample[[All, 1]] // ListLinePlot

It seems that we can approximate it with $a*Exp(-ax)$

Define it and maximize the obj for its parameter $a$

f[x_, a_] := a*Exp[-a*x];

run[a_] := Block[{candsol}, candsol = Table[f[i, a], {i, 100}];
  obj[candsol]];

Table[{i, run[i]}, {i, 0, 1, 0.01}] // SortBy[#, -#[[2]] &] & // First

{0.02, 0.148254}

Plot[f[x, 0.02], {x, 0, 100}];

May it help you solve your problem.

Answer 2

This is what I think the problem boils down to:

obj2[qlist_ /; VectorQ[qlist, NumericQ]] := Module[{Rlist, obj},
   Rlist = FoldList[(#1 - #2) &, R0, qlist];
   (*This is the objective function*)
   ρ^Range[0, Length@qlist - 1].(p*qlist - 
      0.5 c*qlist^2/Most@Rlist)];
grad[qlist_ /; VectorQ[qlist, NumericQ]] := Module[{Rlist, obj, dQ},
   Rlist = FoldList[(#1 - #2) &, R0, qlist];
   dQ = IdentityMatrix[Length@qlist];
   (*This is the jacobian function*)
   ρ^Range[0, Length@qlist - 1].(p*dQ - 
      c*DiagonalMatrix@qlist/Most@Rlist + 
      0.5 c*qlist^2 LowerTriangularize[
         ConstantArray[-1, {Length@qlist, Length@qlist}], -1]/
        Most[Rlist]^2)];

sol = FindMaximum[{obj2[choicevar], 
     And @@ Append[Thread[0 <= choicevar], Total@choicevar <= 1]}, 
    Transpose@{choicevar, 
      RandomVariate[DirichletDistribution[ConstantArray[1, T + 1]]]}, 
    Gradient -> grad[choicevar]]; // AbsoluteTiming
sol[[1]]
(*
  {0.719767, Null}
  0.150353
*)

I computed & coded the gradient grad[] by hand, because it was easier to derive an efficient numerical code that way.

FindMaximum on the OP's functions took over a minute. NMaximize was much slower, naturally, but since the OP mentions that the function is "strictly concave", FindMaximum seems the better tool. If you want to try NMaximize, then the DirichletDistribution is a convenient way to get initial points for this particular domain (NMaximize has trouble finding them automatically):

Method -> {Automatic, (* or whatever *)
  "InitialPoints" -> 
    RandomVariate[DirichletDistribution[ConstantArray[1, T + 1]], 10], 
  "SearchPoints" -> 10}