%-----------------------------------------------------|
%---------------FLUID SIMULATION-1--------------------|
%------------Steady State Distribution----------------|              
%-----------------------------------------------------|
%Author:H. Emre KANKAYA
%DESCRIPTION
%------------------------------------------------------

%This program simulates first order feedback fluid queues
%with multiple thresholds and plots the steady state pdf 
%and cdf of the system concerned.

%Small time epochs are used as simulation steps.
%Look at fluidsimul2.m for another (faster) simulation.

%All sorts of boundaries can be simulated.
%------------------------------------------------------
clear
a=cputime;
disp('----------------------')
disp('SIMULATION-1 STARTED...')
disp('INITIALIZATIONS')
%********************INPUT PARAMETERS***************************
simtime=2000;% simulate system for simtime secs.               
deltat=0.001;%sampling epochs                                  
deltax=0.01;%histogram intervals                               

%EXAMPLE 1
thresholds=[0,1,2]; %terminal and intermediate boundary points
%Q{1}=[-3,1,1,1;1,-3,1,1;1,1,-3,1;1,1,1,-3];%infinitesimal generator in region 1          
%Q{2}=[-3,1,1,1;1,-3,1,1;1,1,-3,1;1,1,1,-3];%infinitesimal generator in region 2          
%QT{1}=[-3,1,1,1;1,-3,1,1;1,1,-3,1;1,1,1,-3];%infinitesimal generator in region 1          
%QT{2}=[-3,1,1,1;1,-3,1,1;1,1,-3,1;1,1,1,-3];%infinitesimal generator in region 2          
%QT{3}=[-3,1,1,1;1,-3,1,1;1,1,-3,1;1,1,1,-3];%infinitesimal generator in region 1          
%R{1}=[0,2,1,-1];%rate in region 1: states 1 and 2 respectively     
%R{2}=[1,-1,2,-1];%rate in region 2: states 1 and 2 respectively     
%RT{1}=[0,2,1,0];%rate at threshold 1(AT ZERO): states 1 and 2 respectively
%RT{2}=[0,0,1,-1];%rate at threshold 2: states 1 and 2 respectively 
%RT{3}=[0,-1,0,-1];%rate at threshold 2: states 1 and 2 respectively 
thresholds=[0,1,2]; %terminal and intermediate boundary points
Q{1}=[-4,1,1,1,1;1,-4,1,1,1;1,1,-4,1,1;1,1,1,-4,1;1,1,1,1,-4];%infinitesimal generator in region 1          
Q{2}=[-4,1,1,1,1;1,-4,1,1,1;1,1,-4,1,1;1,1,1,-4,1;1,1,1,1,-4];%infinitesimal generator in region 2                    
QT{1}=Q{1};
QT{2}=Q{1};
QT{3}=Q{2};
R{1}=[1,-2,2,-4,0];%rates in region 1: states 1 and 2 respectively     
R{2}=[-1,1,3,-1,2];
RT{1}=[1,0,2,0,0];
RT{2}=[0,0,2,-4,0];
RT{3}=[-1,0,0,-1,0];
%**********************INITIALIZATIONS**************************
%SYSTEM PARAMETERS:
%**************************
N=length(thresholds)-1;%number of regions
S=length(R{1});%number of states
noofinter=ceil(thresholds(N+1)/deltax)+1;%for 3: 4 intervals, for 3.1: 5 intervals
%TEMPORARY PARAMETERS:
%**************************
ss=zeros(1,noofinter);
cstate=1;%current state: states 1->S
cregime=1;%current regime: regions 1->N, thresholds:1->N+1
coccupancy=0;%current occupancy
change=0;
t=0;%time parameter
j=1;
statemat=zeros(2,S-1);
inregion=0;%inregion=0 if occupancy is at a threshold 1 otherwise
for j=1:N
    cQ=Q{j};
    for k=1:S
        statematI{j}{k}=zeros(2,S-1);
        for i=1:S
            if i<k
                statematI{j}{k}(1,i)=-cQ(k,i)/cQ(k,k);
                statematI{j}{k}(2,i)=i;
            elseif i>k
                statematI{j}{k}(1,i-1)=-cQ(k,i)/cQ(k,k);
                statematI{j}{k}(2,i-1)=i;
            end
        end
        for i=2:S-1
            statematI{j}{k}(1,i)=statematI{j}{k}(1,i)+statematI{j}{k}(1,i-1);
        end
    end
end

for j=1:N+1
    cQ=QT{j};
    for k=1:S
        statematT{j}{k}=zeros(2,S-1);
        for i=1:S
            if i<k
                statematT{j}{k}(1,i)=-cQ(k,i)/cQ(k,k);
                statematT{j}{k}(2,i)=i;
            elseif i>k
                statematT{j}{k}(1,i-1)=-cQ(k,i)/cQ(k,k);
                statematT{j}{k}(2,i-1)=i;
            end
        end
        for i=2:S-1
            statematT{j}{k}(1,i)=statematT{j}{k}(1,i)+statematT{j}{k}(1,i-1);
        end
    end
end
statemat=statematT{cregime}{cstate};
cQ=QT{1};%current infinitesimal generator
%***************************************************************

%*************************SIMULATION*****************************
disp('SIMULATING THE SYSTEM')
while t<simtime
    %cstate
    %statemat
    %FIND THE STATE OF THE SYSTEM:
    r1=rand;%uniform random number used to simulate changes in states     
    if r1<-cQ(cstate,cstate)*deltat%state changes                                        
        r2=rand;%uniform random number used to find the next state 
        k=1;%temp variable, next state
        
        while k<S
            prob=statemat(1,k);        
            if r2<prob
                    cstate=statemat(2,k);
                    k=S+1;
            end            
            k=k+1;
        end
    end
    
%    cstate
    %FIND THE OCCUPANCY OF THE SYSTEM:
    %determine the amount of change in the occupancy
    %if occupancy is in region r then upper boundary point is thresholds(r+1)
    %lower boundary point is thresholds(r)
    %if inregion=0 then occupancy is at thresholds(cregime)
    if inregion==1
        cQ=Q{cregime};
        statemat=statematI{cregime}{cstate};
        change=deltat*R{cregime}(cstate);
        if change+coccupancy<thresholds(1+cregime)&change+coccupancy>thresholds(cregime)
            coccupancy=coccupancy+change;

        elseif change+coccupancy<=thresholds(cregime)
            if cregime==1
                coccupancy=0;
                inregion=0;
            elseif change+coccupancy>thresholds(cregime-1)
                if RT{cregime}(cstate)<0
                    coccupancy=coccupancy+change;                
                    if coccupancy==thresholds(cregime)
                        inregion=0;
                    else
                        cregime=cregime-1;
                    end
                elseif R{cregime-1}(cstate)>=0
                    coccupancy=thresholds(cregime);
                    inregion=0;
                end
            end
        elseif change+coccupancy>=thresholds(cregime+1)
            if cregime==N
                coccupancy=thresholds(cregime+1);
                inregion=0;
                cregime=N+1;
            elseif change+coccupancy<thresholds(cregime+2)   
                if RT{cregime+1}(cstate)>0
                    coccupancy=coccupancy+change;
                    cregime=cregime+1;           
                    if coccupancy==thresholds(cregime+1)
                        inregion=0;
                    end                   
                elseif R{cregime+1}(cstate)<=0
                    coccupancy=thresholds(cregime+1);
                    cregime=cregime+1;
                    inregion=0;
                end
            end
        end
    
    else%inregion==0
        cQ=QT{cregime};
        statemat=statematT{cregime}{cstate};
        
        
        change=deltat*RT{cregime}(cstate);        
        if RT{cregime}(cstate)>0
           if change+coccupancy<thresholds(cregime+1)
                inregion=1;
                coccupancy=coccupancy+change;                
           end
        elseif RT{cregime}(cstate)<0
           if change+coccupancy>thresholds(cregime-1)
                inregion=1;
                coccupancy=coccupancy+change;
                cregime=cregime-1;
           end 
        end
        
    end
    j=j+1;
    t=j*deltat;
    c=ceil(coccupancy/deltax)+1;
    ss(c)=ss(c)+1;
end


ss2=ss;
for i=2:noofinter
    ss2(i)=ss2(i-1)+ss2(i);
end

ss2=ss2/ss2(noofinter);
ss3=ss/sum(ss)/deltax;
tt=0:deltax:thresholds(N+1);
figure(4);
plot(tt,ss2)
figure(2);
plot(tt(:),ss3(:))
disp('SIMULATION TIME:')
disp(cputime-a);
disp('SIMULATION-2 ENDED...')
disp('----------------------')