clc,clear
%%=========初始化
tic
tanker=[10,8,8,6,6,5];
day=30;
G=10;%油罐總容量（萬噸）
g(1)=0.5*G;h(1)=0;
u=2.5*1e-2;v=1*1e-2;%油罐加油／出油速度
busy=0;%碼頭忙的時間
sum_H=0;% 總運輸油量
K=0;%油輪來的數量
N=2400*day;
T=1;
%%==========時間推進，步長0.01小時
for t=1:N
    %%=======有油輪到來
    if t==T
       
    k=ceil(6*rand(1));H=tanker(k);
    T=floor(-1000*log(rand(1)))+T;
    h(t)=0.95*H+h(t); K=K+1;sum_H=sum_H+0.95*H;
    end
    %%=====設置標幟旗
    if g(t)>=G, flag1=0; end
    if g(t)<=0.8*G, flag1=1; end
    if g(t)>=5, flag2=1; end
    if g(t)<=0.5, flag2=0; end
    if mod(t,2400)>=1800  %時間從6 a.m.算起
        flag3=0;
    else 
        flag3=1;
    end
    %%==========判斷狀態，修改變量
    if h(t)~=0 && flag1==1 && flag3==1
        h(t+1)=h(t)-u; % 油輪卸油
        g(t)=g(t)+u; % 油罐加油
        busy=busy+1; % 碼頭忙
    else
        h(t+1)=h(t);
    end
    if h(t)<=0
        h(t+1)=0;% 卸完油
    end
    if flag2==1
        g(t+1)=g(t)-v; %供油
    else
        g(t+1)=g(t);
    end
end
toc
subplot(211)
plot(g),grid,ylabel('油罐存儲量')
subplot(212)
plot(h),grid,ylabel('油輪碼頭的存儲量')
%%======平均運行時間137s／cpu：2.2GHz,512M內存
% conclusion
busy_rate=busy/N
mean_oil=mean(g)
day_supply=(sum_H-g(N)-h(N))/day
C_waiting=(busy/100)/K