# Solving a dynamic program via recursion

I'm trying to use recursion to solve a joint inventory/ dynamic pricing problem as in monahan, petruzzi and zhao 2004.

I tried to solve for y[t] and k[t] via recursion with the following code:

m = 0.5;
λ = 0.1;
k = 0;
y = 0;
Clear[y];
Clear[k];
r[t_Integer, z_] :=
1/z^m (-(1/λ)*E^(-λ*z) + 1/λ + k[t - 1]*λ*Integrate[(z - w)^m*E^(-λ*w), {w, 0, z}]) // N;
y[t_Integer] := y[t] = ArgMax[{r[t, j], j > y[t - 1]}, j] // N;
k[t_Integer] := k[t] = r[t, y[t]];
y


I get the following message:

$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >> When I solve the problem manually, it works fine: b = 2; m = 1 - 1/b; λ = 1/187; r1[z_] := 1/z^m (-(1/λ)*E^(-λ*z) + 1/λ ) // N; y1 = ArgMax[{r1[g], g >= 0}, g]; y1 r1[y1]  234.953 8.72688  r2[z_] := N[(1/z^m)*(-(1/λ)/E^(λ*z) + 1/λ + r1[y1]*λ*Integrate[(z - w)^m/E^(λ*w), {w, 0, z}])]; y2 = ArgMax[{r2[h], h >= y1}, h]; y2 r2[y2]  486.281 14.432  r3[z_] := N[(1/z^m)*(-(1/λ)/E^(λ*z) + 1/λ + r2[y2]*λ*Integrate[(z - w)^m/E^(λ*w), {w, 0, z}])]; y3 = ArgMax[{r3[i], i >= y2}, i]; y3 r3[y3]  687.782 18.982  ...  Could anybody tell me where's my mistake? This is the updated code, after editing the suggestions: m = 0.5; \[Lambda] = 1/187; r[0, z_] := 0; Clear[k, y]; k = 0; y = 0; r[t_Integer, z_] := 1/z^m (-(1/\[Lambda]) Exp[-\[Lambda] z] + 1/\[Lambda] + k[t - 1] \[Lambda] Integrate[(z - w)^m Exp[-\[Lambda] z], {w, 0, z}]) // N; y[t_Integer] := y[t] = ArgMax[{r[t, j], j > y[t - 1]}, j] // N; k[t_Integer] := k[t] = r[t, y[t]]; y k y k 234.953 8.72688 234.953 11.3039  • try assigning k,y after you Clear[k,y] .. Mar 6, 2014 at 20:27 • thanks for the reply! I do not get the mistake$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >> and y, r are calculated correctly. however in the next iteration y=y, which is not the case when manually solving the steps. Mar 6, 2014 at 21:14

I modified your code to eliminate some numerically induced imaginary fuzz.

ClearAll[y, k];

m = 0.5;
λ = 1/187.;
k = 0;
y = 0;
r[t_Integer, z_] :=
Chop @ N[1/z^m (-(1/λ)*E^(-λ*z) + 1/λ + k[t - 1]*λ*
Integrate[(z - w)^m*E^(-λ*w), {w, 0, z}])];
y[t_Integer] := y[t] = Chop @ N @ ArgMax[{Re@r[t, j], j > y[t - 1]}, j];
k[t_Integer] := k[t] = Chop @ N @ r[t, y[t]];


I was able to get

y

234.953

y

486.281

y

687.782


So, perhaps my modifications are what you are looking for.

• yes :-) sweet, thank you. Mar 6, 2014 at 21:39

Thanks again for the help. Here's the whole commented code I wrote to produce my desired output in a nice form, after updating the parameters. It also simulates random demand und uses another recursion to calculate the inventory level and the optimal prices for the given inventory level:

(*Recursive Formulation of Joint Inventory/ Dynamic-Pricing Problem \
according to Monahan, Petruzzi and Zhao 2004 for price-dependent \
stationary periodic demand Subscript[D, t](p,A)=Subscript[Ap, t]^-b, \
where A is exponentially distributed and the decision maker has a one \
time ordering decision at the beginning of the horizon*)

ClearAll[y, k, q, d];

(*The parameters between begin and end can be changed*)

(*begin*)
(*constant price-elasticity*)
b = 2;

(*Intensity of A*)
\[Lambda] = 1/187;

(*No. of periods*)
T = 15;

(*Cost per unit*)
c = 2;
(*end*)

m = 1 - 1/b;

k = 0;
y = 0;

r[t_Integer, z_] :=
Chop@N[1/z^m (-(1/\[Lambda])*Exp[-\[Lambda] z] + 1/\[Lambda] +
k[t - 1] \[Lambda] Integrate[(z - w)^m Exp[-\[Lambda] w], {w,
0, z}])];
y[t_Integer] := y[t] = Chop@N@ArgMax[{Re@r[t, j], j > y[t - 1]}, j];
k[t_Integer] := k[t] = Chop@N@Re@r[t, y[t]];

(*Optimal Inventory to stock for T periods*)
S = ((k[T] m)/c)^b;

(*Price depending on Inventory in stock at beginning of period*)
p[t_Integer] := If[q[t] == 0, 0, (y[t]/q[t])^(1 - m)];

(*Random price-dependent demand*)
d[t_] := d[t] =
If[t == 0, 0, -(1/\[Lambda]) Log[1 - RandomReal[]] p[t]^-b];

(*quantity in stock*)
q[T] = S;
q[t_Integer] :=
q[t] = If[q[t + 1] - d[t + 1] > 0, q[t + 1] - d[t + 1], 0];

(*Optimal Revenue-to-go Function*)
R[t_] := k[t] q[t]^m;

v = Table[{t, y[t], k[t], R[t], q[t], d[t], p[t]}, {t, 0, T}];
Tablevh =
Prepend[v, {"t",
"\!$$\*SuperscriptBox[SubscriptBox[\(z$$, $$t$$], $$*$$]\)",
"\!$$\*SuperscriptBox[SubscriptBox[\(r$$, $$t$$], $$*$$]\)",
"\!$$\*SubscriptBox[\(R$$, $$t$$]\)(\!$$\*SubscriptBox[\(I$$, $$t\$$]\))", "\!$$\*SubscriptBox[\(I$$, $$t$$]\)",
"\!$$\*SubscriptBox[\(D$$, $$t$$]\)(p,W)",
"\!$$\*SuperscriptBox[SubscriptBox[\(p$$, $$t$$], $$*$$]\)"}];
Print@Grid[Tablevh, Frame -> All];
Print[StringForm["The maximum expected profit is .",
AccountingForm[R[T] - c q[T]^m,
IntegerLength@Round[R[T] - c q[T]^m] + 2]]];