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

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.

  • $\begingroup$ 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. $\endgroup$
    – Xminer
    May 23, 2019 at 1:10
  • $\begingroup$ 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. $\endgroup$
    – Drr777
    May 23, 2019 at 1:31
  • $\begingroup$ 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. $\endgroup$
    – Xminer
    May 23, 2019 at 1:37
  • $\begingroup$ As I posted, I have an analytical solution to this particular problem. Call it a benchmark if you will. $\endgroup$
    – Drr777
    May 23, 2019 at 1:42
  • $\begingroup$ Are you asking for "generic solver for optimal control problem" ? $\endgroup$
    – Xminer
    May 23, 2019 at 8:04

2 Answers 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]/

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

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}
  • $\begingroup$ 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. $\endgroup$
    – Drr777
    May 24, 2019 at 23:44


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*)
     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];




enter image description here

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

enter image description here

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

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

{0.02, 0.148254}

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

enter image description here

May it help you solve your problem.

  • $\begingroup$ This is nice analysis @Xminer, but my problem is not specifically the solution to the depletion problem above. I already know that solution. It is really about trying to setup a quick mathematica applet that might allow for relatively reliable solutions for this class of problems. The corresponding Hamiltonians are almost always strictly concave in the state and control and typically can be shown to have a unique interior solution. So they are very regular problems that I would think Mathematica would handle without much of a fuss! $\endgroup$
    – Drr777
    May 23, 2019 at 3:35
  • $\begingroup$ This makes a lot of sense. Will need to investigate how long this will work before I bug down again. As I always maintained. It would surprise me if mma really could not handle problems that excel solver routinely works with. $\endgroup$
    – Drr777
    May 24, 2019 at 2:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.