NMinimize or NMaximize for a dynamic problem

Question

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. — Mariusz Iwaniuk, May 02 '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? — Drr777, May 02 '19 at 18:23
See this https://mathematica.meta.stackexchange.com/questions/1043/additional-useful-buttons-for-our-m-se-editor — OkkesDulgerci, May 02 '19 at 19:10
Brilliant browser plugin for the greek letters. Thanks! Now I think the question is finally ready for an answer!!! — Drr777, May 02 '19 at 19:20

Roman · Accepted Answer · 2019-05-04T15:17:59.633

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. — Drr777, May 02 '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. — Roman, May 02 '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? — Drr777, May 03 '19 at 09:55
Try choicevar /. sol[[2]] to get the solution out. See https://mathematica.stackexchange.com/a/18706/26598 — Roman, May 03 '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. — Drr777, May 03 '19 at 20:04
Please make this compilation into a separate question (you can reference this answer). I am not a specialist in compilation, and in this case some tricks on vector access will be necessary. — Roman, May 03 '19 at 20:11
Thanks for the continued support. I am making progress and building bigger problems. Could you tell me how I can impose an endogenous constraint on NMinimize that relies on values inside the module? The simplest example for the problem above would be that all elements of xlist are bound above or below or both. More general would be to impose a constraint that ties the alphas of choicevar to an element of xlist. E.g. it could be alpha(t)<=x(t). — Drr777, May 05 '19 at 09:00
Inside the obj function you could have statements that terminate the function with a super-bad result so that NMinimize will avoid this region. For example, if you want all elements of xlist to be nonnegative, you could do If[Min[xlist]<0, Return[1000]] inside of the obj function, so that this function gives a terrible result (1000) whenever one of the $x$-elements becomes negative. — Roman, May 05 '19 at 09:11

NMinimize or NMaximize for a dynamic problem

1 Answers1

Linked