# NMinimize or NMaximize for a dynamic problem

I am trying to solve a relatively simple quadratic discrete time dynamic problem using NMinimize. My code works for short time horizons (T<15) but for T longer, it ends up taking minutes and for T>=20, it never completes, even when I reduce the precision and accuracy requirements. As a reference point, I can easily and rapidly solve the same problem using Excel's Solver for 50 or more periods. I tried all the different NMinimize solution methods without success. Any insight on how to get NMinimize to solve longer problems?

ClearAll[obj, x, y]
obj = (p^(T - 1)/δ)*((d/2)*x[T]^2 + (c/2)*(r*(1 - x[T]/K))^2) +
Sum[p^t*((d/2)*x[t]^2 + (c/2)*y[t]^2), {t, 0, T - 1}];
T = 12;
x[0] = 0.05;
δ = 0.05;
p= 1./(1. + δ);
c = 1.;
r = 1.;
d = 10.;
K = 1.;
y[T] = 0;
For[t = 1, t < T + 1, t = t + 1,
x[t] = x[t - 1]*(1 + r - y[t - 1] - (r/K)*x[t - 1])];
choicevar = Table[y[i], {i, 0, T - 1}];
constraints = Flatten[Table[{y[i] >= 0, y[i] < 1}, {i, 0, T - 1}]];
eq = Prepend[constraints, obj];
AbsoluteTiming[sol = Flatten[NMinimize[eq, choicevar, Method -> {"NelderMead"},
WorkingPrecision -> 8, PrecisionGoal -> 6, AccuracyGoal -> 6]]]

• People here generally like users to post code as Mathematica code instead of images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. May 2 '19 at 17:04
• Thanks, I figured out how to get the code in. Though is there a way of getting the greek letters to display properly? May 2 '19 at 18:23
• May 2 '19 at 19:10
• Brilliant browser plugin for the greek letters. Thanks! Now I think the question is finally ready for an answer!!! May 2 '19 at 19:20
• Possibly related is this old MathGroup thread. May 2 '19 at 22:14

By forcing the recursion in obj to be done numerically at every step, instead of doing it analytically once and for all, I can compute T=50 in less than 40 seconds without even specifying any options to NMinimize:

T = 50;
Δ = 0.05;
p = 1./(1. + Δ);
c = 1.;
r = 1.;
d = 10.;
K = 1.;

obj[ylist_ /; VectorQ[ylist, NumericQ]] := Module[{xlist},
(* calculate the list of x[t]-values *)
xlist = FoldList[#1*(1 + r - #2 - (r/K)*#1) &, 0.05, ylist];
(* evaluate the obj function *)
(p^(T - 1)/Δ)*((d/2)*xlist[[T + 1]]^2 + (c/2)*(r*(1 - xlist[[T + 1]]/K))^2) +
Sum[p^t*((d/2)*xlist[[t + 1]]^2 + (c/2)*ylist[[t + 1]]^2), {t, 0, T - 1}]]

choicevar = Table[y[i], {i, 0, T - 1}];

AbsoluteTiming[
sol = NMinimize[
Prepend[Thread[0 <= choicevar < 1], obj[choicevar]],
choicevar]]


{37.7883, {8.89576, {y[0] -> 0.577873, y[1] -> 0.645035, y[2] -> 0.716978, y[3] -> 0.785202, y[4] -> 0.838979, y[5] -> 0.869077, y[6] -> 0.876428, y[7] -> 0.876755, y[8] -> 0.876754, y[9] -> 0.876754, y[10] -> 0.876755, y[11] -> 0.876755, y[12] -> 0.876755, y[13] -> 0.876754, y[14] -> 0.876755, y[15] -> 0.876755, y[16] -> 0.876755, y[17] -> 0.876756, y[18] -> 0.876753, y[19] -> 0.876756, y[20] -> 0.876754, y[21] -> 0.876755, y[22] -> 0.876756, y[23] -> 0.876754, y[24] -> 0.876754, y[25] -> 0.876756, y[26] -> 0.876753, y[27] -> 0.876756, y[28] -> 0.876755, y[29] -> 0.876755, y[30] -> 0.876756, y[31] -> 0.876752, y[32] -> 0.876758, y[33] -> 0.876754, y[34] -> 0.876755, y[35] -> 0.876754, y[36] -> 0.876755, y[37] -> 0.876758, y[38] -> 0.876752, y[39] -> 0.876754, y[40] -> 0.876762, y[41] -> 0.87675, y[42] -> 0.876755, y[43] -> 0.87676, y[44] -> 0.876748, y[45] -> 0.87676, y[46] -> 0.876752, y[47] -> 0.876762, y[48] -> 0.876745, y[49] -> 0.87676}}}

I suppose that by compiling the obj function this can be sped up a lot more. Also, using choicevar = Table[Unique[y], {i, 0, T - 1}]; instead of what you've used gives a bit of a speedup (after all, you don't need to care about the names of the optimization variables here).

Here I've put together some more speedups: avoiding a loop in obj by using only vector processing, and using Unique variables instead of indexed ones. This gives about a factor of two over the above code.

T = 50;
Δ = 0.05;
p = 1./(1. + Δ);
c = 1.;
r = 1.;
d = 10.;
K = 1.;

pt = p^Range[0, T - 1]/2;
obj[ylist_ /; VectorQ[ylist, NumericQ]] := Module[{xlist},
xlist = FoldList[#1*(1 + r - #2 - r/K #1) &, 0.05, ylist];
p^(T-1)/(2Δ)*(d*xlist[[T+1]]^2 + c*r^2*(1-xlist[[T+1]]/K)^2) +
(d*Most[xlist]^2 + c*ylist^2).pt]

choicevar = Table[Unique[y], {i, 0, T - 1}];

First@AbsoluteTiming[
sol = NMinimize[Prepend[Thread[0 <= choicevar < 1], obj[choicevar]], choicevar];]


22.4621

{sol[[1]], choicevar /. sol[[2]]}


{8.89576, {0.577873, 0.645035, 0.716978, 0.785202, 0.838979, 0.869077, 0.876428, 0.876755, 0.876754, 0.876754, 0.876755, 0.876755, 0.876755, 0.876754, 0.876755, 0.876755, 0.876755, 0.876756, 0.876753, 0.876756, 0.876754, 0.876755, 0.876756, 0.876754, 0.876754, 0.876756, 0.876753, 0.876756, 0.876755, 0.876755, 0.876756, 0.876752, 0.876758, 0.876754, 0.876755, 0.876754, 0.876755, 0.876758, 0.876752, 0.876754, 0.876762, 0.87675, 0.876755, 0.87676, 0.876748, 0.87676, 0.876752, 0.876762, 0.876745, 0.87676}}

• I suspected very strongly that forcing the recursion was the source of the problem. So it makes total sense that this is where your code makes the gains. This being said, I do not grasp precisely how your code proceeds (my lack of familiarity with what the first line does). So, I will have to study it more closely. May 2 '19 at 20:57
• I simplified the obj function a bit, using FoldList instead of a recursion. It's cleaner and a bit faster too. May 2 '19 at 21:14
• Using Unique[y] actually substantially speeds up long problems (ie Large T). But the result I get is of this form {y$113532 -> 0.5778717004, y$113533 -> 0.6450358597}. I have been looking for more than an hour on how to extract the values of y from these lists of local variables without any luck. How does one do it? May 3 '19 at 9:55
• Try choicevar /. sol[[2]] to get the solution out. See mathematica.stackexchange.com/a/18706/26598 May 3 '19 at 10:07
• Thanks again. As I push my luck and make this and very similar problems longer, I am starting to run into long runtime again. Could you illustrate what/how to compile? Cheers. May 3 '19 at 20:04