%function [inventorycost]=Qsolve(v,l)
%v,l分別是路線的集合，數組及顧客的個數
%計算路線總長度
%所有參數部分


v=[3 2 1];
l=2;
customer=4;
demand=[155 98 474
        178 191 365
        78 88 365
        50 65 200
        38 133 0
         22 76 0
        174 193 0];
    depot=[38 133
           22 76
           174 193];
       UD=[3 8 6 7];
       
c=25;
cm=1;
h=0.5;
A=20;
hs=2;


%計算路線總長度
leng=0;
for i=2:l
leng=distance(demand,v(i),v(i+1))+leng;
end

leng=leng+distance(demand,v(1)+customer,v(2))+distance(demand,v(l+1),v(1)+customer);
leng
%計算路線總需求
D=sum(demand(v(2:l+1),3));
%計算路線上提前期的平均需求
u=0;
for i=2:l+1
   u=UD(v(i))+u;
end



R0=0;
BR0=(u-R0)^2/(2*u);
Q0=sqrt((2*D*(cm*leng+c+hs*BR0))/h)
R1=u-(u*h*Q0)/(hs*D)
BR1=(u-R1)^2/(2*u);
Q1=sqrt((2*D*(cm*leng+c+hs*BR1))/h) 
n=1;
while abs(Q1-Q0)>0.01
    n=n+1
   R0=R1;
   BR0=(u-R0)^2/(2*u);
   Q0=sqrt((2*D*(cm*leng+c+hs*BR0))/h);
   R1=u-(u*h*Q0)/(hs*D);
   BR1=(u-R1)^2/(2*u);
   Q1=sqrt((2*D*(cm*leng+c+hs*BR1))/h);
end
%計算庫存成本
D
u                           
Q1
R1
inventorycost=((c+cm*leng)*D)/Q1+(Q1/2+R1-u/2)*h+(D/Q1)*A+(hs*BR1*D)/Q1;