% Value iteration using discounetd cost for machine replacement problem
% Stopping criteria using epsilon
% Algorithm loops over all states.
clc
clear
V=0;
Vminnext=0;
Vmin=0;
%State 0,1,2,3 is 1,2,3,4 for Matlab
C=[	0	4000	6000
	1000	4000	6000
	3000	4000	6000
	10000	10000	6000
];

T1=[0	0.875	0.0625	0.0625
0	1	0	0
1	0	0	0
];

T2=[0	0.75	0.125	0.125
0	1	0	0
1	0	0	0
];

T3=[0	0	0.5	0.5
0	1	0	0
1	0	0	0
];

T4=[0	0	0	1
0	1	0	0
1	0	0	0
];

Aph=0.9;
E=0.01 %epsilon

%***************************************
ite=0; %iteration number
V=zeros(4,3);
Vmin=zeros(4,1);
Vminnext=zeros(4,1);
VV(1,:)=V(:,1)';
%***************************************
ite=1;
for a =1 :1:3
V(1,a)= C(1,a)+ (Aph*(T1(a,:)*Vmin));
end
temp1=min(V(1,:));

for a =1 :1:3
V(2,a)= C(2,a)+ (Aph*(T2(a,:)*Vmin));
end
temp2=min(V(2,:));

for a =1 :1:3
V(3,a)= C(3,a)+ (Aph*(T3(a,:)*Vmin));
end
temp3=min(V(3,:));

for a =1 :1:3
V(4,a)= C(4,a)+ (Aph*(T4(a,:)*Vmin));
end
temp4=min(V(4,:));

Vminnext(1,1)=temp1;
Vminnext(2,1)=temp2;
Vminnext(3,1)=temp3;
Vminnext(4,1)=temp4;
VV(2,:)=Vminnext(:,1)';
%***************************************
n=3;
                
while ite<100000 
ite=ite+1;
Vmin=Vminnext;
for a =1 :1:3
V(1,a)= C(1,a)+ (Aph*(T1(a,:)*Vmin));
end
temp1=min(V(1,:));

for a =1 :1:3
V(2,a)= C(2,a)+ (Aph*(T2(a,:)*Vmin));
end
temp2=min(V(2,:));

for a =1 :1:3
V(3,a)= C(3,a)+ (Aph*(T3(a,:)*Vmin));
end
temp3=min(V(3,:));

for a =1 :1:3
V(4,a)= C(4,a)+ (Aph*(T4(a,:)*Vmin));
end
temp4=min(V(4,:));
Vminnext(1,1)=temp1;
Vminnext(2,1)=temp2;
Vminnext(3,1)=temp3;
Vminnext(4,1)=temp4;

VV(n,:)=Vminnext(:,1)';
n=n+1;

if abs(Vminnext(1,1)-Vmin(1,1))<E  % check stopping criteria
    if abs(Vminnext(2,1)-Vmin(2,1))<E
        if abs(Vminnext(3,1)-Vmin(3,1))<E
            if abs(Vminnext(4,1)-Vmin(4,1))<E
                Actualite=ite
                ite=100002;

            end
        end
    end
end
end
%******************************************
% finding arg min for action
for a =1 :1:3
if temp1-V(1,a)<0.001 & temp1-V(1,a)>-0.001
Act1=a;
end
end
for a =1 :1:3
if temp2-V(2,a)<0.001 & temp2-V(2,a)>-0.001
Act2=a;
end
end
for a =1 :1:3
if temp3-V(3,a)<0.001 & temp3-V(3,a)>-0.001
Act3=a;
end
end
for a =1 :1:3
if temp4-V(4,a)<0.001 & temp4-V(4,a)>-0.001
Act4=a;
end
end

Policy(1,1)=Act1;
Policy(1,2)=Act2;
Policy(1,3)=Act3;
Policy(1,4)=Act4;
Policy

plot(VV)