Inventory problem. Computing assignment 3. Math 504, spring 2008. CSUF
by Nasser Abbasi

Define constants for the problem (give in class to use)

s=3;
S=5;

define demand function  (note the use of conditions on arguments)

demand[x_]:=(x+1)/;x≤5
demand[x_]:=(3/4)(13-x)/; 5<x≤13;
demand[x_]:=0;
points=Range[0,13]

define the and p (i, j) functions (see handout)

p[i_,j_]:=Total[Table[c*demand[v],{v,S,Length[points]}]]/;(j==0&&i<s);
p[i_,j_]:=Total[Table[c*demand[v],{v,i,Length[points]}]]/;(j==0&&i≥s);
p[i_,j_]:=c*demand[S-j]/;(i<s&& (0<j≤S));
p[i_,j_]:=c*demand[i-j]/;(i≥ s&& (0<j≤ i)); 
p[i_,j_]:=0;

define a function which determines the ' c' constant by making at area = 1 under the demand function

demandTable=Table[demand[x],{x,points[[1]],points[[-1]]}];
c=1/Total[demandTable];
Print["c=",N[c]];
demandTable=demandTable*c(*normalize it to make it probability*)

Print the demand function numerical values

t=TableForm[N[demandTable],TableHeadings->{Apply[ToString,points,1]}]

Make a histogram of the demand function

Needs["BarCharts`"]
t=Table[{x,demandTable[[x+1]],1},{x,points[[1]],points[[-1]]}];
GeneralizedBarChart[t,AxesLabel->{"x","demand(x)"},Frame->True,FrameTicks->{{Automatic,Automatic},{points,None}},FrameLabel->{"x","demand(x)","demand function"}]

and fill in the initial P matrix (ONE step matrix)

A=Table[p[i,j],{i,0,S},{j,0,S}];
A//N//MatrixForm

Part (A)

Now we investigate what happens to (A^n)as n-> ∞ we start with n=2 and try few increasing values

m=MatrixPower[A,3] //N//MatrixForm

m=MatrixPower[A,5] //N//MatrixForm

m=MatrixPower[A,12] //N//MatrixForm

We can stop now. The matrix (A^n) has converged. Let try some very large value of n to check

m=MatrixPower[A,100] //N//MatrixForm

Conclusion

As n-> ∞ the matrix A converges to a fixed value shown above. Each column has the same entries in its rows. IN addition, all entries are non-zero (this is due to the fact that there are no transiet states). In addition, the above means that starting from any initial state (Subscript[(π^0), i]) the final state ((Subscript[π, k])^((∞))) is the same for any state k

See if final state vector depends or not on the initial state

make some random initial states to use

numberOfStates=6;
n=30;
initialStates={{1,0,0,0,0,0},{.5,.5,0,0,0,0},{.4,.2,.4,0,0,0},{.1,.1,.4,.4,0,0},{.2,.3,.1,.1,.2,.1},{0,1,0,0,0,0}};
initialStates//MatrixForm

Now multiply each of the above by A raised to large power

For[i=1,i<=6,i++,
finalPi=initialStates[[i,All]].MatrixPower[A,n];
Print[MatrixForm[initialStates[[i,All]]],MatrixForm[N[finalPi]]];
]

Lets verify w = wp theory

A//MatrixForm

j=2;
Sum[finalPi[[i]]*A[[i,j]],{i,1,Length[finalPi]}]//N

finalPi[[j]]//N