clear all;
clc;

tStart = tic;
filename='Input.csv';
%sheet='Table';

N1=120;%number of customers type 1 in system
N2=120;%number of customers type 2 in system

AI = getenv('PBS_ARRAY_INDEX');
NAI = str2num (AI);


% N= NAI + 1;
% cellReferenceB = sprintf('B%d', N);
% lam1=xlsread (filename,sheet,cellReferenceB,'basic');
% cellReferenceC = sprintf('C%d', N);
% lam2=xlsread (filename,sheet,cellReferenceC,'basic');
% cellReferenceD = sprintf('D%d', N);
% mu1=xlsread (filename,sheet,cellReferenceD,'basic');%mu1 service rate for queue1
% cellReferenceE = sprintf('E%d', N);
% mu2=xlsread (filename,sheet,cellReferenceE,'basic');%mu2 service rate for queue2
% cellReferenceF = sprintf('F%d', N);
% theta1=xlsread (filename,sheet,cellReferenceF,'basic');
% cellReferenceG = sprintf('G%d', N);
% theta2=xlsread (filename,sheet,cellReferenceG,'basic');
% cellReferenceH = sprintf('H%d', N);
% gamma1=xlsread (filename,sheet,cellReferenceH,'basic');
% cellReferenceI = sprintf('I%d', N);
% gamma2=xlsread (filename,sheet,cellReferenceI,'basic');
input = csvread('Input.csv',NAI,1,[NAI 1 NAI 8]);
lam1 = input(1);
lam2 = input(2);
mu1 = input(3);
mu2 = input(4);
theta1 = input(5);
theta2 = input(6);
gamma1 = input(7);
gamma2= input(8);

h1 = 0; % Holding cost at station 1
h2 = 0; % Holding cost  at station 2
% 
% lam1 = 0.3;%arrival rate to queue1
% lam2 = 0.3;%arrival rate to queue1
% 
% mu1 = 3; %service rate for queue1
% mu2 = 1;%service rate for queue2
% 
% theta1 = 2;%reneging rate at queue1
% theta2 = 1;%reneging rate at queue2

W(1:N1+1,1:N2+1)=0;
V(1:N1+1,1:N2+1)=0;%value matrix
A(1:N1+1,1:N2+1)=0;%action decision

%equality_percision = 0.001;
optimality_percision = 0.01;

nu=lam1+lam2+mu1+mu2+N1*theta1+N2*theta2;
iter=20000;

for step=1:iter
    for i=1:N1+1
        for j=1:N2+1
            if i == 1 && j == 1
                b1 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(i+1,j) + lam2*W(i,j+1)...
                    + mu1*W(i,j) + mu2*W(i,j)...
                    + (i-1)*theta1*W(i,j) + (j-1)*theta2*W(i,j) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                b2 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(i+1,j) + lam2*W(i,j+1)...
                    + mu1*W(i,j) + mu2*W(i,j)...
                    + (i-1)*theta1*W(i,j) + (j-1)*theta2*W(i,j) + (N1-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));               
                V(i,j) = b1;
                A(i,j) = 1;
            end
            
            if i == 1 && 1 < j && j < N2+1
                b1 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(i+1,j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i,j) + mu2*W(i,j)...
                    + (i-1)*theta1*W(i,j) + (j-1)*theta2*W(i,j-1) + (N1-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                b2 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(i+1,j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i,j) + mu2*W(i,j-1)...
                    + (i-1)*theta1*W(i,j) + (j-1)*theta2*W(i,j-1) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                V(i,j) = b2;
                A(i,j) = 2;
            end
            
            if i == 1 && j == N2+1
                b1 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(i+1,j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i,j) + mu2*W(i,j)...
                    + (i-1)*theta1*W(i,j) + (j-1)*theta2*W(i,j-1) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                b2 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(i+1,j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i,j) + mu2*W(i,j-1)...
                    + (i-1)*theta1*W(i,j) + (j-1)*theta2*W(i,j-1) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                V(i,j) = b2;
                A(i,j) = 2;
            end
            
            if 1 < i && i< N1+1 && j == 1
                b1 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(min(i+1,N1+1),j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i-1,j) + mu2*W(i,j)...
                    + (i-1)*theta1*W(i-1,j) + (j-1)*theta2*W(i,j) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                b2 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(min(i+1,N1+1),j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i,j) + mu2*W(i,j)...
                    + (i-1)*theta1*W(i-1,j) + (j-1)*theta2*W(i,j) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                V(i,j) = b1;
                A(i,j) = 1;
            end
            
            if i == N1+1 && j == 1
                b1 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(min(i+1,N1+1),j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i-1,j) + mu2*W(i,j)...
                    + (i-1)*theta1*W(i-1,j) + (j-1)*theta2*W(i,j) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                b2 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(min(i+1,N1+1),j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i,j) + mu2*W(i,j)...
                    + (i-1)*theta1*W(i-1,j) + (j-1)*theta2*W(i,j) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                V(i,j) = b1;
                A(i,j) = 1;
            end
            
            if 1 < i && i < N1+1 && 1 < j && j< N2+1
                b1 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(min(i+1,N1+1),j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i-1,j) + mu2*W(i,j)...
                    + (i-1)*theta1*W(i-1,j) + (j-1)*theta2*W(i,j-1) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                b2 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(min(i+1,N1+1),j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i,j) + mu2*W(i,j-1)...
                    + (i-1)*theta1*W(i-1,j) + (j-1)*theta2*W(i,j-1) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                V(i,j) = b2;
                A(i,j) = 2;
            end
            
            if 1 < i && i < N1+1 && j == N2+1
                b1 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(min(i+1,N1+1),j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i-1,j) + mu2*W(i,j)...
                    + (i-1)*theta1*W(i-1,j) + (j-1)*theta2*W(i,j-1)+ (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                b2 = (1/(nu))*(h1*(i-1) + h2*(j-1)+gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(min(i+1,N1+1),j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i,j)+ mu2*W(i,j-1)...
                    + (i-1)*theta1*W(i-1,j) + (j-1)*theta2*W(i,j-1) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                V(i,j) = b2;
                A(i,j) = 2;
            end
            
            if i == N1+1 && 1 < j && j < N2+1
                b1 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(min(i+1,N1+1),j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i-1,j) + mu2*W(i,j)...
                    + (i-1)*theta1*W(i-1,j) + (j-1)*theta2*W(i,j-1) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                b2 = (1/(nu))*(h1*(i-1) + h2*(j-1)+gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(min(i+1,N1+1),j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i,j)+ mu2*W(i,j-1)...
                    + (i-1)*theta1*W(i-1,j) + (j-1)*theta2*W(i,j-1) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                V(i,j) = b2;
                A(i,j) = 2;
            end
            
            if i == N1+1 && j == N2+1
                b1 = (1/(nu))*(h1*(i-1) + h2*(j-1) + gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(min(i+1,N1+1),j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i-1,j) + mu2*W(i,j)...
                    + (i-1)*theta1*W(i-1,j) + (j-1)*theta2*W(i,j-1) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                b2 = (1/(nu))*(h1*(i-1) + h2*(j-1)+gamma1*theta1*(i-1) + gamma2*theta2*(j-1)...
                    + lam1*W(min(i+1,N1+1),j) + lam2*W(i,min(j+1,N2+1))...
                    + mu1*W(i,j) + mu2*W(i,j-1)...
                    + (i-1)*theta1*W(i-1,j) + (j-1)*theta2*W(i,j-1) + (N2-(i-1))*theta1*W(i,j) + (N2-(j-1))*theta2*W(i,j));
                V(i,j) = b2;
                A(i,j) = 2;
            end
        end
    end
    
    VW = V-W;
    
    alpha = min(min(min(VW)));
    beta = max(max(max(VW)));
    
    eps = (beta-alpha)/ alpha;
    if eps < optimality_percision, break, end
    W=V;
    
end

step;
% mmax;
% mmin;
cost=(alpha+beta)/2; % total average cost during a period = g
Avgcost = cost*nu % average cost per unit time = g/ (avg lenght of a period = 1/nu)

m=zeros(1,3);

m(1,1)= Avgcost;
m(1,2) = NAI;
m(1,3) = A(4,5);

dlmwrite('outputcost.csv',m, '-append');

%m1 = [full(sparse(1,1,1,1,N2+1)).*NAI;A]; %add a row with the NAI
%dlmwrite('output.csv',m1, '-append','roffset', 1);

tEnd = toc(tStart);
fprintf('%d minutes and %f seconds\n',floor(tEnd/60),rem(tEnd,60));