0
$\begingroup$

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]]]
$\endgroup$
  • 3
    $\begingroup$ 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. $\endgroup$ – Mariusz Iwaniuk May 2 at 17:04
  • $\begingroup$ Thanks, I figured out how to get the code in. Though is there a way of getting the greek letters to display properly? $\endgroup$ – Drr777 May 2 at 18:23
  • $\begingroup$ See this mathematica.meta.stackexchange.com/questions/1043/… $\endgroup$ – OkkesDulgerci May 2 at 19:10
  • $\begingroup$ Brilliant browser plugin for the greek letters. Thanks! Now I think the question is finally ready for an answer!!! $\endgroup$ – Drr777 May 2 at 19:20
  • $\begingroup$ Possibly related is this old MathGroup thread. $\endgroup$ – Daniel Lichtblau May 2 at 22:14
2
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Drr777 May 2 at 20:57
  • $\begingroup$ I simplified the obj function a bit, using FoldList instead of a recursion. It's cleaner and a bit faster too. $\endgroup$ – Roman May 2 at 21:14
  • $\begingroup$ 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? $\endgroup$ – Drr777 May 3 at 9:55
  • $\begingroup$ Try choicevar /. sol[[2]] to get the solution out. See mathematica.stackexchange.com/a/18706/26598 $\endgroup$ – Roman May 3 at 10:07
  • $\begingroup$ 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. $\endgroup$ – Drr777 May 3 at 20:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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