mp_tsp_nc.m

marko chain project

function MP_TSP_NC(n)
%Multi Period TSP with Stochastic Demands
%n is the # of hospitals, not including DC

%Input Data
%n = 5;%Number of hospitals
% n = 2;
nState = 3^n;
p0 = 0.3;   %Probability that no demand occurs
p1 = 0.5;   %Probability of occurence of normal demand
p2 = 0.2;   %Probability of occurence of urgent demand

%In defining the states no demand:0 ,regular:1, urgent:2  
%Distance matrix
%for 2 hospitals
% C = [0 10 25;10 0 30;25 30 0];
% C = dists(nccity('XY'),nccity('XY'),2);
%for 5 hospitals
% C = [0 18 35 35 50 15;18 0 21 45 68 50;35 21 0 24 50 20;35 45 24 0 27 20;50 68 50 27 0 30;15 50 20 30 30 0];

[C,cityname,cityxy] = nccitydist(n);

%Generating all possible 3^n states
ST = zeros(nState,n); 
for i = 1:nState
    char = dec2base(i-1,3,n);
    for j = 1:n
        ST(i,j) = sscanf(char(j),'%i');
    end
end

%Generating the transition matrices 
trans_state.prob = [];

%Generating the intermediate states: the state at the end of the first day
trans_state.intermstate = [];

%Generating the reward matrices
trans_state.cost = [];

for i = 1:nState
    ind = find(ST(i,:)==1);
    trans_state(i).prob = zeros(2^length(ind),nState);
    trans_state(i).cost = zeros(2^length(ind),nState);
    if(length(ind)==0)
        trans_state(i).intermstate = zeros(1,n);
    else
        trans_state(i).intermstate = zeros(2^length(ind),n);
        for j = 2:2^length(ind)
            char = dec2bin(j-1,n);
            char = fliplr(char);%ask this
            for k = 1:length(ind)
                trans_state(i).intermstate(j,ind(k)) = sscanf(char(k),'%i');
            end
        end
        trans_state(i).intermstate = 2*trans_state(i).intermstate;
    end
end

%Generating the transition probabilities
for i = 1:nState    %for each state
    ind = find(ST(i,:)==1);
    for k = 1:2^length(ind) %for each alternative belongs to that state
        if length(find (ST(i,:)>trans_state(i).intermstate(k,:)))==0
            trans_state(i).cost(k,:) = zeros(1,nState);
        else
           c = find (ST(i,:)>trans_state(i).intermstate(k,:));
           c1=[n+1,c];
           %c1 = [3,c];
           c2 = C(c1,c1);
           [rte,TC] = tspnneighbor(c2,1);
           trans_state(i).cost(k,:) = TC*ones(1,nState);
        end
        for j = 1:nState %for each state that the process is moving to
            dummy = 1;
            for m = 1:n
                if(trans_state(i).intermstate(k,m)==0)  
                    switch(ST(j,m))
                        case(0)
                            dummy = dummy*p0;
                        case(1)
                            dummy = dummy*p1;
                        case(2)
                            dummy = dummy*p2;
                    end
                else
                    switch(ST(j,m))
                        case(0)
                            dummy = dummy*0;
                        case(1)
                            dummy = dummy*0;
                        case(2)
                            dummy = dummy*1;
                    end
                end
            end
            trans_state(i).prob(k,j) = dummy;
            if trans_state(i).prob(k,j)==0
               trans_state(i).cost(k,j) = 1000000;
            end
        end      
    end
end

q.value = []; %alternatives
q.alt = [];
q.best = [];
for i = 1:nState
    q(i).value = sum(trans_state(i).prob(:,:).*trans_state(i).cost(:,:),2);
    [q(i).best,q(i).alt] = min(q(i).value);
end

epsilon = 1E-6;
nCheap = 20;

sum_s = 1;
difference = false;

v = zeros(1,nState);
V = zeros(1,nState);
vold = zeros(1,nState);

while((sum_s > epsilon)||(difference == true))
    for nch = 1:nCheap
        vold = v;
        for i = 1:nState
            V(i) = q(i).best + sum(trans_state(i).prob(q(i).alt,:).*v);
        end
        v = V - V(end);
    end
    sum_s = sum(abs(v - vold));
    difference = false;
    dummy_p = zeros(1,nState);
    for i = 1:nState
        ind = find(ST(i,:)==1);
        dummy_V = zeros(1,2^length(ind));
        for k = 1:2^length(ind) 
            dummy_V(k) = q(i).value(k) + sum(trans_state(i).prob(k,:).*v);
        end
        [V(i),dummy_p(i)] = min(dummy_V);
    end
    policy = [];
    for t = 1:nState
        policy(t) = q(t).alt;
    end
    if(dummy_p ~= policy)
        difference = true
        for t = 1:nState
            q(t).alt = dummy_p(t);
        end
    end
end