# NMaximize not always returning solution

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.

At this point, Excel Solver seems to be a better solution but I refuse to believe that! Thanks.

• Isn't the step size for candidate solution in the NMaximize too large? you may want to use another solver (AMPL, GAMS,or Excel,of course) that you have. – Xminer May 23 '19 at 1:10
• A the moment, I have reasons for wanting a Mathematica based solution. I have a hard time believing that Mathematica cannot at least be on par with Excel Solver. I don't expect it to be as good as GAMS, but that's beside the point. If Mathematica can perform, it is best for me to make it work. – Drr777 May 23 '19 at 1:31
• ok,hope someone post the good answer :) by the way,the solution itself exists. try  obj[RandomVariate[UniformDistribution[{0, 0.010}], T]] this returns positive value. – Xminer May 23 '19 at 1:37
• As I posted, I have an analytical solution to this particular problem. Call it a benchmark if you will. – Drr777 May 23 '19 at 1:42
• Are you asking for "generic solver for optimal control problem" ? – Xminer May 23 '19 at 8:04

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]]]},
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}

• Good ideas @Michael E2. I am trying to make this flexible so I went for setting random initial points rather than the gradient. The DirichletDistribution functions better but I substituted it for a uniform distribution for a more general set up for other problems. With a scaling factor, it works better than a hard constraint on choicevariables. Also, Total@choicevar<=1 is slower than Min[Rlist<...] in obj[]. NMinimize on original problem, with Uniform solves in under 30s. Much slower than with gradient, but much less time than coding the gradient! FindMin is not as fast with this setup. – Drr777 May 24 '19 at 23:44

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}];