% cost and decision matrix for inv control problem pg 26 prob 3 in Denardo

clc
clear
clc
%set up the model parameters
DD=0;
%**********************
DD= [2 4 1 ];
maxhold=2; %2 or 4
maxprod=3;
startinv=1;
endinv=0;
prevprod =1;
%************************
cc=0;
ccc=0;
cccc=0;
xx=0;
xxx=0;
xxxx=0;
a=0;
sizeDD=size(DD);
nstate= maxhold+1;
nstage= sizeDD(1,2)+1;
cc=zeros(nstage*nstate, nstage*nstate);
xx=zeros(nstage*nstate, nstage*nstate);
countr=1;
countc=nstate+1;
m=1;
D=DD;
sizeD=size(D);

% these steps are for the first month
k=1;
i=startinv;
x=0;
temps=0;
for j=0:1:maxhold
        x=j-i+D(1,1);
        
        
        if x>=0 & x<=maxprod
            temps(m,1)=i;
        temps(m,2)=j;
        temps(m,3)=x;
            temps(m,4)=k;
        end
        m=m+1;
end

for k=1:1:1
    
    D=D(1,k)

for i= 0:1:maxhold
    for j=0:1:maxhold
        x=j-i+D;
        
        
        if x>=0 && x<=maxprod
          
            if x-prevprod==0
                if k==1
                    c=3*x;
                elseif k==2
                    c=5*x;
                elseif k==3
                    c=3*x;
                end
            elseif x-prevprod ==1 || x-prevprod ==-1
               if k==1
                    c=3*x+2;
                elseif k==2
                    c=5*x+2;
                elseif k==3
                    c=3*x+2;
               end  
             elseif x-prevprod ==2 || x-prevprod ==-2
               if k==1
                    c=3*x+4;
                elseif k==2
                    c=5*x+4;
                elseif k==3
                    c=3*x+4;
               end  
             elseif x-prevprod ==3 || x-prevprod ==-3
               if k==1
                    c=3*x+6;
                elseif k==2
                    c=5*x+6;
                elseif k==3
                    c=3*x+6;
               end  
            end
            if c==0 % arc exists
                c=0.0001;
                end
        
        elseif x<0
            c=100; % infeasible arcs
        elseif x>maxprod
            c=100; % infeasible arcs
        end
       
        a(i+1,j+1)=c;
        xx(countr,countc)=x;
        cc(countr,countc)=c;
        countc=countc+1;
    
    end
    countr=countr+1;
    countc=countc-nstate;
   
end
countc=countc+nstate;
     D= DD;       
    a
end

% these steps are for month 2 and onwards
% month t+1 requires action data from month t for a given starting
% inventory of month t+1. This data is stored in tempss.
tempss=temps;
temps=0;
sizetempss=size(tempss);
m=1;
for k=2:1:sizeD(1,2)
    
    D=D(1,k)

for i= 0:1:maxhold
    for j=0:1:maxhold
        x=j-i+D;
        
        
        if x>=0 && x<=maxprod
            temps(m,1)=i;
        temps(m,2)=j;
        temps(m,3)=x;
        temps(m,4)=k;
           for u=1:1:sizetempss(1,1)
               if i==tempss(u,2);
                   prevprod=tempss(u,3);
               end
           end
               
            if x-prevprod==0
                if k==1
                    c=3*x;
                elseif k==2
                    c=5*x;
                elseif k==3
                    c=3*x;
                end
            elseif x-prevprod ==1 || x-prevprod ==-1
               if k==1
                    c=3*x+2;
                elseif k==2
                    c=5*x+2;
                elseif k==3
                    c=3*x+2;
               end  
             elseif x-prevprod ==2 || x-prevprod ==-2
               if k==1
                    c=3*x+4;
                elseif k==2
                    c=5*x+4;
                elseif k==3
                    c=3*x+4;
               end  
             elseif x-prevprod ==3 || x-prevprod ==-3
               if k==1
                    c=3*x+6;
                elseif k==2
                    c=5*x+6;
                elseif k==3
                    c=3*x+6;
               end  
            end
            if c==0 % arc exists
                c=0.0001;
                end
        
        elseif x<0
            c=100; % infeasible arcs
        elseif x>maxprod
            c=100; % infeasible arcs
        end
       
        a(i+1,j+1)=c;
        xx(countr,countc)=x;
        cc(countr,countc)=c;
        countc=countc+1;
         m=m+1;
    end
    countr=countr+1;
    countc=countc-nstate;
   
end
countc=countc+nstate;
     D= DD;       
    a
    tempss=0;
    tempss=temps;
temps=0;
sizetempss=size(tempss);
m=1;
end
cc;
xx;
 

 % the following code will trim c and X  to have one source and one
 % sink only
 %for inv control problem
 %one source
 c=cc;
 X=xx;
sizec=size(c);
sizex=size(X);
for i=1:1:sizec(1,2)
c2(1,i)=c(startinv+1,i);
X2(1,i)=X(startinv+1,i);
end
rowc=2;
 for i=nstate+1:1:sizec(1,1)
    for j= 1:1:sizec(1,2)
    c2(rowc,j)=c(i,j);
    X2(rowc,j)=X(i,j);
    end
 rowc=rowc+1;
 end
 
%one sink
sizec2=size(c2);
sizeX2=size(X2);
colc=1;
for i=1:1:(nstate*(nstage-1))
     for j= 1:1:sizec2(1,1)
     c3(j,colc)=c2(j,i);
     X3(j,colc)=X2(j,i);
     end
colc=colc+1;
end
for i=1:1:sizec2(1,1)
c3(i,1+(nstate*(nstage-1)))=c2(i,((nstage*nstate)-nstate+endinv+1));
X3(i,1+(nstate*(nstage-1)))=X2(i,((nstage*nstate)-nstate+endinv+1));
end

%the following will trim the zeros in the first few colums and the last few
%rows.Keep only first column and last row with all zeros
cc=c3; 
xx=X3;
cccc=0;
 sizecc=size(cc);
 ccc= cc(1:sizecc(1,1)-nstate+1,:);
 cccc= ccc(:, nstate:sizecc(1,2));
 cccc;
 
 xxxx=0;
 sizexx=size(xx);
 xxx= xx(1:sizexx(1,1)-nstate+1,:);
 xxxx= xxx(:, nstate:sizexx(1,2));
 xxxx;
 
 c=cccc % c has the correct rows and columns with one source and sink
 X=xxxx % X has the correct rows and columns with one source and sink
 
 
 % these files are written in the same folder where this code is saved.
 xlswrite('costinv1.xlsx',c)
 xlswrite('decisioninv1.xlsx',X)
 
