File content as of revision 1:ed464543ca57:

#include "mbed.h"

Serial pc(USBTX, USBRX); // tx, rx
DigitalOut myled(LED2);
DigitalOut led_a(LED3);


#define u1          8
#define M1          1
#define u2          3
#define M2          3
#define u3          15
#define M3          1
#define pQ3Q1       0.1         // probability of sending Q3 output to the Q1 input
#define pQ3Out      0.9        // probability of delivering Q3 output as system output
#define tf          50
#define lambda      10
#define pQ1I_Q2I    0.5
#define sz_a        1000

int Q1[sz_a]= {0};              // initial Q1
int Q2[sz_a]= {0};              // initial Q2
int Q3[sz_a]= {0};              // initial Q3
float t = 0;               // initial simulation time
float t_events_Q1[sz_a] = {0};     // time stamps for the k events
float t_events_Q2[sz_a] = {0};     // time stamps for the k events
float t_events_Q3[sz_a]= {0};     // time stamps for the k events
int k=1;                 // event number

double SI = 0.1;

double Q1I1 = tf;
double Q1I2 = 0;
double Q1O = 0;

double Q2I1 = 0;
double Q2O = 0;

double Q3I1 = 0;
double Q3I2 = 0;
double Q3O = tf;


double SO = 0;

float Q1_lambda [sz_a]= {lambda*0.5};
float Q2_lambda [sz_a]= {lambda*0.5};
float Q3_lambda [sz_a]= {1/(u1*M1) + 1/(u2*M2)};

float t_Q1=0;
float t_Q2=0;
float t_Q3=0;

float eventtime_Q1=0;
float eventtime_Q2=0;
float eventtime_Q3=0;
float eventtime=0;

double min_2(double x, double y);
double min_3(double x, double y, double z);
int max_2(int x, int y);
double max_3(double x, double y, double z);
double var_rand_q1(int x);
double var_rand_s(int x);
double var_rand_q2(int x);
double var_rand_q3(int x);

int Q1k;
int Q2k;
int Q3k;



uint32_t x_i_s=7; //semilla aleatoria
uint32_t x_i_q1=11;
uint32_t x_i_q2= 17;
uint32_t x_i_q3=3;
int i=0;
int main()
{
      
    while( t <= tf ){
        myled=0;
    //exprnd(1/freq) simulate an exponential pdf with parameter = average time = 1 / arriving rate    
    
    // System input     
    
    SI = var_rand_s(lambda); 
    SI= SI/100;
    
    
    // Q1
    Q1I1 = Q3O/pQ3Q1;
    Q1I2 = SI/0.5;
    printf("Sl: %lf \n", SI);
    Q1_lambda [k]= (1/Q1I1 + 1/Q1I2);
    
    //printf("Q1I1: %lf \n", Q1I1);
    //printf("Q1I2: %lf \n", Q1I2);
    //printf("Q1_lambda: %f \n", Q1_lambda[k]);
     
    if(Q1[k] > 0){
        int q1O= u1*M1; // service time in the queue
        Q1O= var_rand_q1(q1O);
        Q1O= Q1O/100;
    }
    else{
        Q1O= min_2(Q1I1,Q1I2)+1.0;
    }
    printf("Q1O: %lf \n", Q1O);   
    
    // Update queue Q1 
    t_events_Q1[k]=(t+min_3(Q1I1,Q1I2,Q1O));
    //printf("Q1K :%d \n",Q1[k]);
    if ( (Q1O < Q1I1) && (Q1O < Q1I2) ){
        Q1k=Q1[k]- 1;  
        Q1[k+1]=max_2(0,Q1k); // Q1 is reduced by one
    
    }
    else        {
        Q1[k+1]=Q1[k]+1; // Q1 is increased by one due to Q1 input (only one of the two inputs is considered)
    }
    
    //printf("Q1K :%d \n",Q1[k+1]);
    eventtime_Q1 = min_3(Q1I1, Q1I2, Q1O);
    t_Q1 = t;
    t_Q1 = t_Q1 + eventtime_Q1;
        
    // Q2
    Q2I1 = SI/pQ1I_Q2I;
    //printf("Q2I1 :%lf \n",Q2I1);
    Q2_lambda[k] = (1.0/Q2I1);
    //printf("Q2_lambda :%f \n",Q2_lambda[k]);
    
    if(Q2[k] > 0){
        int q2O = u2*M2; // service time in the queue
        Q2O= var_rand_q2(q2O);
        Q2O= Q2O/100;
    }
    else{
        Q2O = Q2I1 + 1.0;
    }
    printf("Q20 :%f \n",Q2O);
    
    // Update queue Q2 
    t_events_Q2[k]=(t+min_2(Q2I1, Q2O));
    //printf("Q2[k]: %d\n",Q2[k]);
    if ( Q2O < Q2I1 ){
        Q2k=Q2[k]- 1;
    //    printf("Q2K %d: \n",Q2k);
        Q2[k+1]=max_2(0,Q2k); // Q2 is reduced by one
    }
    else{
        Q2[k+1]= (Q2[k])+1; // Q2 is increased by one due to Q2 input
    }
    //printf("Q2[k]: %d\n",Q2[k+1]);
    eventtime_Q2 = min_2(Q2I1, Q2O);
    t_Q2 = t;
    t_Q2 = t_Q2 + eventtime_Q2;
    
    
    //Q3
    Q3I1 = Q1O;
    Q3I2 = Q2O;
    //printf("Q3I1: %lf \n",Q3I1);
    //printf("Q3I2: %lf \n",Q3I2);
    Q3_lambda [k]= 1/Q3I1 + 1/Q3I2;
    //printf("Q3_lambda: %f \n",Q3_lambda[k]);
    
    if(Q3[k] > 0){
        int q3O = u3*M3; // service time in the queue
        Q3O= var_rand_q3(q3O);
    //    printf("Q3O_ant %lf", Q3O);   
        Q3O= Q3O/100; 
    }
    else{
        Q3O = min_2(Q1O, Q2O) + 1.0;
    }
    
    printf("Q3O: %lf \n",Q3O);
    // Update queue Q3 
    t_events_Q3[k]= t+min_3(Q3I1, Q3I2, Q3O);    
    //printf("Q3[k]: %d \n", Q3[k]);
    if ( (Q3O < Q3I1) && (Q3O < Q3I2) ){
        Q3k=Q3[k]- 1;
    //    printf("Q3k: %f \n", Q3k); 
        Q3[k+1]=max_2(0,Q3k); // Q3 is reduced by one
    }
    else        {
        Q3[k+1]= Q3[k] + 1; // Q3 is increased by one due to Q3 input (only one of the two inputs is considered)

    }
    //printf("Q3[k]: %d \n",Q3[k+1]);
    eventtime_Q3 = min_3(Q3I1, Q3I2, Q3O);
    t_Q3 = t;
    t_Q3 = t_Q3 + eventtime_Q3;
    
    // System output
    SO = Q3O/pQ3Out;
    
    eventtime = max_3(eventtime_Q1, eventtime_Q2, eventtime_Q3);
    t = t + eventtime;
    k=k+1;
    
    myled=1;
    wait(0.02);      
    }          
    
    
    led_a=0;
    wait(.5);
    pc.printf("k=");
    pc.printf("%d\r\n",k);
    wait(.5);
    
    pc.printf("Q1\n");
    for (i=0; i<k; i++){
        pc.printf("%d, ",Q1[i]);
        
        }
    
    pc.printf("\n Q2\n");  
    for (i=0; i<k; i++){
        pc.printf("%d, ",Q2[i]);
        
        }  
        
    pc.printf("\n Q3\n");  
    for (i=0; i<k; i++){
        pc.printf("%d, ",Q3[i]);
        
        }  
    
    pc.printf("\n Q1_lambda\n");  
    for (i=0; i<k; i++){
        pc.printf("%f, ",Q1_lambda[i]);
        
        }  
    
    pc.printf("\n Q2_lambda\n");  
    for (i=0; i<k; i++){
        pc.printf("%f, ",Q2_lambda[i]);
        
        }  
    pc.printf("\n Q3_lambda\n");    
    for (i=0; i<k; i++){
        pc.printf("%f, ",Q3_lambda[i]);
        
        }  
    
    pc.printf("\n t_events_Q1\n");    
    for (i=0; i<k; i++){
        pc.printf("%f, ",t_events_Q1[i]);
        
        }  
    pc.printf("\n t_events_Q2\n");    
    for (i=0; i<k; i++){
        pc.printf("%f, ",t_events_Q2[i]);
        
        } 
        
    pc.printf("\n t_events_Q3\n");    
    for (i=0; i<k; i++){
        pc.printf("%f, ",t_events_Q3[i]);
        
        } 
            
    
    led_a=1;
}


//-----------------------Función para obtener el minimo de 2 números ---------------------------------------------//
double min_2(double x, double y){
     double min=0;
     if (x<y){
         min= x;
         }
     else {
        min=y;
        }
        return min;
    
    }
    
///-----------------------Función para obtener el minimo de 3 números ---------------------------------------------//
double min_3(double x, double y, double z){
     double min=0;
     if ((x<y)&&(x<z)){
         min= x;
         }
     else if ((y<x)&&(y<z)){
         min= y;
         }
     else {
         min=z;
         }
         return min;
}

//-----------------------Función para obtener el miáximo de 2 números ---------------------------------------------//
int max_2(int x, int y){
    int max=0;
    if(x>y){
        max=x;
        }     
    else {
        max=y;
        }
        return max;
    }


//-----------------------Función para obtener el máximo de 3 números ---------------------------------------------//
double max_3(double x, double y, double z){
     double max=0;
     if ((x>y)&&(x>z)){
         max= x;
         }
     else if ((y>x)&&(y>z)){
         max= y;
         }
     else {
         max=z;
         }
         return max; 
    }

double var_rand_s(int x){
    int32_t a_s = 40014;
    int32_t m_s=2147283563;
    double exprnd_s;
    int32_t x_i1_s=(a_s*x_i_s) % m_s;
    x_i_s=x_i1_s; 
    double U_s= float(x_i1_s)/float(m_s); //Uniforme [0,1] 

     
    if(U_s>=0.2 && U_s<=0.9){
        exprnd_s=-1.0*float(x)*log(U_s);
        }
    else {
        x_i1_s=(a_s*x_i_s) % m_s;
        U_s= float(x_i1_s)/float(m_s);
        }
   // printf("exprnd_s: %lf", exprnd_s);
    return exprnd_s;
    }
    
double var_rand_q1(int x){
    int32_t a_q1 = 40014;
    int32_t m_q1=2147283563;
    double exprnd_q1;
    int32_t x_i1_q1=(a_q1*x_i_q1) % m_q1;
    x_i_q1=x_i1_q1; 
    double U_q1= float(x_i1_q1)/float(m_q1); //Uniforme [0,1]    
    
    
    if(U_q1>=0.4 && U_q1<=0.9){
        exprnd_q1=-1.0*float(x)*log(U_q1);
        }
    else {
        x_i1_q1=(a_q1*x_i_q1) % m_q1;
        U_q1= float(x_i1_q1)/float(m_q1);
        }
        return exprnd_q1;
    }
    
    double var_rand_q2(int x){
    int32_t a_q2 = 40014;
    int32_t m_q2=2147283563;
    double exprnd_q2;
    int32_t x_i1_q2=(a_q2*x_i_q2) % m_q2;
    x_i_q2=x_i1_q2; 
    double U_q2= float(x_i1_q2)/float(m_q2); //Uniforme [0,1]
    
    if(U_q2>=0.4 && U_q2<=0.9){
        exprnd_q2=-1.0*float(x)*log(U_q2);
        }
    else {
        x_i1_q2=(a_q2*x_i_q2) % m_q2;
        U_q2= float(x_i1_q2)/float(m_q2);
        }
        return exprnd_q2;
    }
    
    double var_rand_q3(int x){
    int32_t a_q3 = 40014;
    int32_t m_q3=2147283563;
    double exprnd_q3;
    int32_t x_i1_q3=(a_q3*x_i_q3) % m_q3;
    x_i_q3=x_i1_q3; 
    double U_q3= float(x_i1_q3)/float(m_q3); //Uniforme [0,1]
    
    if(U_q3>=0.4 && U_q3<=0.9){
        exprnd_q3=-1.0*float(x)*log(U_q3);
        }
    else {
        x_i1_q3=(a_q3*x_i_q3) % m_q3;
        U_q3= float(x_i1_q3)/float(m_q3);
        }
    
    return exprnd_q3;
    }