set.seed(23111996);
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)

simul = function(N=100, arrival_rate=.2, processing_rate=1, debug=FALSE, E1 = 1, E0 = 0.5 ) {
    Arrival = cumsum(rexp(n=N, rate=arrival_rate));  
    Service = rexp(n=N, rate =processing_rate); # rep(N,x=1); rgamma(N, shape=.1,rate=.1)
    Remaining = rep(N, x=NA);
    Completion = rep(N, x=NA);
    Energy = rep(N, x=NA);
    t = 0;
    CurrentTask = NA;
    NextArrival = 1;
    while (TRUE) {
        dtA = NA;
        dtC = NA;
        if(length(Arrival[Arrival>t])>0) {
            dtA = head(Arrival[Arrival>t],n=1)-t  # temps jusqu'à la prochaine arrivée
        }
        if(!is.na(CurrentTask)) {
            dtC = Remaining[CurrentTask]; # temps jusqu'à la prochaine terminaison
        }
        if(is.na(dtA) & is.na(dtC)) {
            break;
        } 
        dt = min(dtA,dtC,na.rm=T)
        
        # Mettre à jour comme il faut:
        #   la date
        t = t + dt;
        #   les arrivées
        if((NextArrival <=N) & (Arrival[NextArrival] <= t)) { ## je met un <= et pas un == car 3-2.9!=0.1 ...
            Remaining[NextArrival] = Service[NextArrival];
            if (NextArrival == 1){
              Energy[NextArrival] = 0 + Remaining[NextArrival]*E1 + Arrival[NextArrival]*E0;
            }
            else{
              Energy[NextArrival] = Energy[NextArrival-1] + Remaining[NextArrival]*E1 + (Arrival[NextArrival] - Arrival[NextArrival-1])*E0;
            }
            NextArrival = NextArrival + 1;
            
        }
        #   le remaining 
        if(!is.na(CurrentTask)) {
            Remaining[CurrentTask] = Remaining[CurrentTask] - dt ;
            if(Remaining[CurrentTask] <= 0) {
                Completion[CurrentTask] = t;
                Remaining[CurrentTask] = NA;
            }
            CurrentTask = NA;
        }
        #   et currentTask (politique d'ordonnancement: FIFO)
        WaitingList=(1:N)[!is.na(Remaining)];
        if(length(WaitingList)>0) {
            CurrentTask = head(WaitingList,n=1);
        }
    }
    return(data.frame(
        Arrival_rate=arrival_rate,
        Processing_rate=processing_rate,
        Arrival = Arrival,
        Service = Service,
        Completion = Completion,
        Energy = Energy
        ))
        
}

#### Avec un N petit : N = 100

s1 = data.frame()

for (r in  seq(from=.1, to=.9, by=.1)) {
    s1 = rbind(s1,simul(N=100, arrival_rate=r, processing_rate=1, debug=FALSE, E1 = 1, E0 = 0.5));
}

t_resp1 = s1 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(t_resp1 = mean(Completion - Arrival),
              Response_se1 = sd(Completion - Arrival)/sqrt(n()))

s2 = data.frame()

for (r in  seq(from=.1, to=.9, by=.1)) {
    s2 = rbind(s2,simul(N=100, arrival_rate=r, processing_rate=2, debug=FALSE, E1 = 4, E0 = 0.25));
}

t_resp2 = s2 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(t_resp2 = mean(Completion - Arrival),
              Response_se2 = sd(Completion - Arrival)/sqrt(n()))

s3 = data.frame()

for (r in  seq(from=.1, to=.9, by=.1)) {
    s3 = rbind(s3,simul(N=100, arrival_rate=r, processing_rate=3, debug=FALSE, E1 =1*1/3 + 2/3*4, E0 = 0.25+0.5));
} 

t_resp3 = s3 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(t_resp3 = mean(Completion - Arrival),
              Response_se3 = sd(Completion - Arrival)/sqrt(n()))


response_plot = cbind(t_resp1,t_resp2,t_resp3)

p1_10 = ggplot(response_plot, aes(x=Arrival_rate, y=t_resp1)) + geom_line(color='red') + ylim(0,7) + geom_errorbar(aes(ymin = t_resp1 - 2*Response_se1, ymax = t_resp1 + 2*Response_se1))
 
p1_10 = p1_10 + geom_line(data=response_plot, aes(x=Arrival_rate, y=t_resp2)) + geom_errorbar(aes(ymin = t_resp2 - 2*Response_se2, ymax = t_resp2 + 2*Response_se2))

p1_10 = p1_10 + geom_line(data=response_plot, aes(x=Arrival_rate, y=t_resp3)) + geom_errorbar(aes(ymin = t_resp3 - 2*Response_se3, ymax = t_resp3 + 2*Response_se3), color='green') + scale_x_continuous(breaks = seq(from =0 , to = 1, by = 0.1))
 

energys1 = s1 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(EnergyCons1 = mean(Energy))
 
energys2 = s2 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(EnergyCons2 = mean(Energy))

energys3 = s3 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(EnergyCons3 = mean(Energy))
 
energy_plot = cbind(energys1,energys2,energys3)
 
p2_10 = ggplot(energy_plot, aes(x=Arrival_rate, y=EnergyCons1)) + geom_line(color="red") + scale_x_continuous(breaks = seq(from =0 , to = 1, by = 0.1) )
p2_10 = p2_10 + geom_line(data=energy_plot, aes(x=Arrival_rate, y=EnergyCons2)) + geom_line(data=energy_plot, aes(x=Arrival_rate, y=EnergyCons3),color = "green")
 

#### Avec un N moyen : N = 1000

s1 = data.frame()

for (r in  seq(from=.1, to=.9, by=.1)) {
    s1 = rbind(s1,simul(N=1000, arrival_rate=r, processing_rate=1, debug=FALSE, E1 = 1, E0 = 0.5));
}

t_resp1 = s1 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(t_resp1 = mean(Completion - Arrival),
              Response_se1 = sd(Completion - Arrival)/sqrt(n()))

s2 = data.frame()

for (r in  seq(from=.1, to=.9, by=.1)) {
    s2 = rbind(s2,simul(N=1000, arrival_rate=r, processing_rate=2, debug=FALSE, E1 = 4, E0 = 0.25));
}

t_resp2 = s2 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(t_resp2 = mean(Completion - Arrival),
              Response_se2 = sd(Completion - Arrival)/sqrt(n()))

s3 = data.frame()

for (r in  seq(from=.1, to=.9, by=.1)) {
    s3 = rbind(s3,simul(N=1000, arrival_rate=r, processing_rate=3, debug=FALSE, E1 =1*1/3 + 2/3*4, E0 = 0.25+0.5));
} 

t_resp3 = s3 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(t_resp3 = mean(Completion - Arrival),
              Response_se3 = sd(Completion - Arrival)/sqrt(n()))


response_plot = cbind(t_resp1,t_resp2,t_resp3)

p1_100 = ggplot(response_plot, aes(x=Arrival_rate, y=t_resp1)) + geom_line(color='red') + ylim(0,7) + geom_errorbar(aes(ymin = t_resp1 - 2*Response_se1, ymax = t_resp1 + 2*Response_se1))
 
p1_100 = p1_100 + geom_line(data=response_plot, aes(x=Arrival_rate, y=t_resp2)) + geom_errorbar(aes(ymin = t_resp2 - 2*Response_se2, ymax = t_resp2 + 2*Response_se2))

p1_100 = p1_100 + geom_line(data=response_plot, aes(x=Arrival_rate, y=t_resp3)) + geom_errorbar(aes(ymin = t_resp3 - 2*Response_se3, ymax = t_resp3 + 2*Response_se3), color='green') + scale_x_continuous(breaks = seq(from =0 , to = 1, by = 0.1))
 


energys1 = s1 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(EnergyCons1 = mean(Energy))
 
energys2 = s2 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(EnergyCons2 = mean(Energy))

energys3 = s3 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(EnergyCons3 = mean(Energy))
 
energy_plot = cbind(energys1,energys2,energys3)
 
p2_100 = ggplot(energy_plot, aes(x=Arrival_rate, y=EnergyCons1)) + geom_line(color="red") + scale_x_continuous(breaks = seq(from =0 , to = 1, by = 0.1) )
p2_100 = p2_100 + geom_line(data=energy_plot, aes(x=Arrival_rate, y=EnergyCons2)) + geom_line(data=energy_plot, aes(x=Arrival_rate, y=EnergyCons3),color = "green")
 


##### Avec un N grand : N = 10000
s1 = data.frame()

for (r in  seq(from=.1, to=.9, by=.1)) {
    s1 = rbind(s1,simul(N=10000, arrival_rate=r, processing_rate=1, debug=FALSE, E1 = 1, E0 = 0.5));
}

t_resp1 = s1 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(t_resp1 = mean(Completion - Arrival),
              Response_se1 = sd(Completion - Arrival)/sqrt(n()))

s2 = data.frame()

for (r in  seq(from=.1, to=.9, by=.1)) {
    s2 = rbind(s2,simul(N=10000, arrival_rate=r, processing_rate=2, debug=FALSE, E1 = 4, E0 = 0.25));
}

t_resp2 = s2 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(t_resp2 = mean(Completion - Arrival),
              Response_se2 = sd(Completion - Arrival)/sqrt(n()))

s3 = data.frame()

for (r in  seq(from=.1, to=.9, by=.1)) {
    s3 = rbind(s3,simul(N=10000, arrival_rate=r, processing_rate=3, debug=FALSE, E1 =1*1/3 + 2/3*4, E0 = 0.25+0.5));
} 

t_resp3 = s3 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(t_resp3 = mean(Completion - Arrival),
              Response_se3 = sd(Completion - Arrival)/sqrt(n()))


response_plot = cbind(t_resp1,t_resp2,t_resp3)

p1_1000 = ggplot(response_plot, aes(x=Arrival_rate, y=t_resp1)) + geom_line(color='red') + ylim(0,7) + geom_errorbar(aes(ymin = t_resp1 - 2*Response_se1, ymax = t_resp1 + 2*Response_se1))
 
p1_1000 = p1_1000 + geom_line(data=response_plot, aes(x=Arrival_rate, y=t_resp2)) + geom_errorbar(aes(ymin = t_resp2 - 2*Response_se2, ymax = t_resp2 + 2*Response_se2))

p1_1000 = p1_1000 + geom_line(data=response_plot, aes(x=Arrival_rate, y=t_resp3)) + geom_errorbar(aes(ymin = t_resp3 - 2*Response_se3, ymax = t_resp3 + 2*Response_se3), color='green') + scale_x_continuous(breaks = seq(from =0 , to = 1, by = 0.1))
 

energys1 = s1 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(EnergyCons1 = mean(Energy))
 
energys2 = s2 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(EnergyCons2 = mean(Energy))

energys3 = s3 %>% group_by(Arrival_rate, Processing_rate) %>%
    summarize(EnergyCons3 = mean(Energy))
 
energy_plot = cbind(energys1,energys2,energys3)
 
p2_1000 = ggplot(energy_plot, aes(x=Arrival_rate, y=EnergyCons1)) + geom_line(color="red") + scale_x_continuous(breaks = seq(from =0 , to = 1, by = 0.1) )
p2_1000 = p2_1000 + geom_line(data=energy_plot, aes(x=Arrival_rate, y=EnergyCons2)) + geom_line(data=energy_plot, aes(x=Arrival_rate, y=EnergyCons3),color = "green")

S1 est en rouge S2 en noir et S3 en vert

Simulation du temps de reponse avec N = 100

## Warning: Removed 2 rows containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_errorbar).

Simulation du temps de reponse avec N = 1000

Simulation du temps de reponse avec N = 10000

## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_errorbar).

Simulation de l’energie avec N = 100

Simulation de l’energie avec N = 1000