diff --git a/ICC2015_Sourcecode/DVFS_type1_ICC.m b/ICC2015_Sourcecode/DVFS_type1_ICC.m
new file mode 100644
index 0000000000000000000000000000000000000000..f5f5ac10ecb1958d95f86c423dad677bba7239e4
--- /dev/null
+++ b/ICC2015_Sourcecode/DVFS_type1_ICC.m
@@ -0,0 +1,495 @@
+%ScriptMax 1a cond + ni cond blocco totale version 4
+close all;
+clear all;
+clc;
+
+tic;
+VMNo=100;                                        % VM numbers
+N = 2000;                                        %number of workloads
+
+nVM=[];
+% +-+-+-+-+-+ parameters +-+-+-+-+-+ 
+for i=1:VMNo
+nVM = [nVM,i];                        %number of virtual machines
+end
+
+
+speed_cost_all = zeros(length(nVM),N);           %matrix of speed cost_all (1 x workloads)
+switch_cost_all = zeros(length(nVM),N);          %matrix of switch cost (1 x workloads)
+channel_cost_all = zeros(length(nVM),N);         %matrix of channel cost (1 x workloads)
+tot_cost_all = zeros(length(nVM),N);             %matrix of total cost (1 x workloads)
+
+
+k_e = 0.005;                            %frequency change cost  parameter
+delta = 0.1;                          %task time
+f0 = 0;                             %percentage of fmax of initial frequency
+iterations = 10000;                 %max number of iteration
+gamma = 0.5;
+
+tol_workload_allocated = 0.01;      %tolerance 1st condition
+epsilon = 0.01;                     %tolerance 2nd condition
+R_tot = 100;                        %capacity of NetDC
+%L_tot_vect = 6+ round((10-6) * rand(1,N)); % PMR=1.25, L_tot=8+-2
+%save('workloads8PMR125No1000.mat', 'L_tot_vect');                  %loading 100 fixed random workloads
+L_tot_vect=load('workloads8PMR125No2000.mat','L_tot_vect');                                %loading 100 fixed random workloads
+L_tot_vect = L_tot_vect.L_tot_vect;
+
+T_tot = 5;      
+beta = 0.41/(T_tot - delta)^1.995;%time out
+%beta = 0.34/(T_tot - delta)^1;%time out
+
+rule_counter=zeros(1,N);                    % (insit of notsatisfactin cases) states for the each incoming workload facing in 3 rules
+rule1_counter=0;
+rule2_counter=0;
+rule3_counter=0;
+
+Zeta_vector=[0.5:0.1:100];
+%VMs=[2, 100];
+
+% iter_mat1 = zeros(1,N); 
+% mu_mat1 = zeros(N,iterations);
+% ni_mat1 = zeros(N,iterations);
+% f_opt_mat1=zeros(N,iterations);
+% L_mat1=zeros(N,iterations);
+
+%for M = nVM(1):VMNo
+for M = 100:100
+%for M = 1:length(VMs)
+    %M=VMs(M);
+    M
+    f_max = 105 * ones(1,M);        % (Mbit/s) max working frequency
+    f_zero = (f0/100) * f_max;      %initial working frequency
+    %f_zero = 0;      %initial working frequency
+    W = ones(1,M);                  %band of each channel (MHz)
+    %Zeta = 0.5 * ones(1,M);         %channel parameters (mW)
+    Zeta = Zeta_vector(1:M); 
+    f_opt_new=zeros(1,M);
+%     N_0 = ones(1,M);g = ones(1,M);Zeta = (N_0 .* W) / g;R = [0:0.1:R_max];P_net = Zeta .* (2^(R/W) - 1);
+    E_max = 60 * ones(1,M);         %max power consumption (@fmax)
+    omega = ones(1,M);              %relative power consumption (virtualization cost)
+    Delta = delta *  ones(1,M);     %task time
+    Delta_max = max(Delta);         %max task time
+    L_b = zeros(1,M);               %background task size for each VM
+    
+%    speed_cost = zeros(N,iterations);           %matrix of speed cost (workloads x iteration)
+%    switch_cost = zeros(N,iterations);          %matrix of switch cost (workloads x iteration)
+%    channel_cost = zeros(N,iterations);         %matrix of channel cost (workloads x iteration)
+%    tot_cost = zeros(N,iterations);             %matrix of total cost (workloads x iteration)
+    
+        
+    temp=2.*k_e+(2.*E_max.*omega./(f_max.^2));  %temp var for calculation of fi*
+    Th = 2 * (Zeta./W) * log(2);                %threshold = derivation of Pnet(Ri) 
+    
+    iter_mat = zeros(1,N);                      %number of iteration necessary for 1st condition for each workload
+    
+%     MatrixRis = zeros(1+M+M,N);                 %matrix of results, row: 1=final_iter 2=final_alpha 3=tot_cost 4=final_ni 5=2nd_condition 6=(ni<esp & g(.)<esp)
+%     L_allocated = zeros(N,iterations);          %matrix of total allocation for each iteration fon each workload 
+%     f_opt_mat = zeros(N,iterations);            %matrix of frequency allocation for each iteretion for each workload FOR THE 1ST VM
+%     L_mat = zeros(N,iterations);                %matrix of  allocation for each iteretion for each workload FOR THE 1ST VM
+%     ni_mat = zeros(N,iterations);               %matrix of "ni" value for each iteretion for each workload FOR THE 1ST VM
+%     mu_mat = zeros(N,iterations);
+%     
+    %matrix of "mu" value for each iteretion for each workload
+    alpha_mat = zeros(N,iterations);            %matrix of "alpha" value for each iteretion for each workload
+    pippo_mat = zeros(N,iterations);            %matrix of "pippo" value for each iteretion for each workload FOR THE 1ST VM
+    
+    for l = 1:N                                 % workload number
+        %check feasibility
+        feas_workload = (sum(f_max.*Delta-L_b) >= L_tot_vect(l));
+        feas_time = (L_tot_vect(l) <= R_tot.*(T_tot-Delta_max)./2);
+        feas_back = min(Delta.*f_max >= L_b);
+        feasibility = feas_workload && feas_time && feas_back;
+        if ~ feasibility
+            error('Problem unfeasible for the workload %d',l);
+        else
+            fprintf('Problem feasibile for the workload %d\n',l);
+            Delta_load = zeros(1,iterations);
+            L = zeros(iterations,M);
+            alpha = zeros(1,iterations);
+            pippo = zeros(1,iterations);
+            V = zeros(1,iterations);
+            V_pippo = zeros(1,iterations);
+            mu = zeros(1,iterations);
+            ni = zeros(iterations,M);
+            y = zeros(iterations,M);
+            f_opt = zeros(iterations,M);
+            i = 2;
+            while i <= iterations
+                bandierina = 0;
+                k=0;
+                cond3 = zeros(1,M);
+                cond2 = zeros(1,M);
+                %calculate of MU
+                mu(i) = max(0,(mu(i-1) - ((alpha(i-1))) * (sum(L(i-1,:)) - L_tot_vect(l)) ) );
+                
+                %calculate of temp variable
+                y(i,:) = max(0,( ((T_tot-Delta)/2) .* W .* log2(mu(i) ./ Th) ) );
+                
+                %calculate of ni; the old was: ni(i,:) = max(0,( ni(i-1,:) + (alpha(i-1) * ( L(i-1,:) - (Delta.*f_opt(i-1,:))))))
+                ni(i,:) = (ni(i-1,:) + (pippo(i-1) * (y(i,:) - Delta.*f_opt(i-1,:))));
+
+                %calculate of f_opt
+                f_star = (2*k_e*f_zero + ni(i,:).*Delta)./ temp;
+                f_opt(i,:) = max(L_b/Delta,min(f_star,f_max)); 
+                               
+                %calculate of L
+                L(i,:) = min(f_opt(i,:) .* Delta(1,:),y(i,:));
+
+%                 ni(i,ni(i,:)<=0.01) = 0;
+%                 for h=1:M
+%                     if ni(i,h)>0
+%                         L(i,h)=f_opt(i,h) .* Delta(1,h);
+%                     else
+%                         L(i,h) = y(i,h);
+%                     end
+%                 end
+                
+                % 2nd and 3rd condition
+                for x = 1:M
+                    cond3(1,x) = (((abs(ni(i,x)))<=epsilon)||((abs(L(i,x)-f_opt(i,x)*Delta(1,x)))<=epsilon));
+                    cond2(1,x) = L(i,x)<=f_opt(i,x)*Delta(1,x);
+                end
+                
+%                 if (  (sum(cond3(1,:)) < M)||(sum(cond2(1,:)) < M)  )
+%                     alpha(i-1) = 0.5 * alpha(i-1);
+%                     bandierina = 1;
+%                     k = k + 1;
+%                     fprintf('alpha loop...);
+%                 end  
+                if (bandierina == 0 && k<100)
+                    alpha(i) = max(0,min(beta,alpha(i-1) - gamma * V(i-1) * (sum(L(i-1,:)) - L_tot_vect(l))));
+                    V(i) = (1 - alpha(i-1)) * V(i-1) - (sum(L(i-1,:)) - L_tot_vect(l));
+                    pippo(i) = max(0,min(0.001,pippo(i-1) - gamma * V_pippo(i-1) *(Delta(1)*f_opt(i-1,1)- y(i,1))));
+                    V_pippo(i) = (1 - pippo(i-1)) * V_pippo(i-1) - (Delta(1)*f_opt(i-1,1)- y(i,1));
+                    Delta_load(i)=sum(L(i,:))-L_tot_vect(l);
+%                     speed_cost(l,i) = sum(((f_opt(i,:)./f_max).^2).*omega.*E_max);                                % Computing (11.1)
+%                     switch_cost(l,i) = sum(k_e.*(f_opt(i,:)-f_zero).^2);                                          % Switching (11.1)
+%                     channel_cost(l,i) = (T_tot-Delta(1)).*sum(Zeta .* (2.^(2*L(i,:)./((T_tot-Delta).*W))-1));     % Networking
+%                     tot_cost(l,i) = speed_cost(l,i) + switch_cost(l,i) + channel_cost(l,i);                       % Overall
+%                     
+%                    L_allocated(l,i) = sum(L(i,:));
+%                     L_mat(l,i) = L(i,1);
+%                     f_opt_mat(l,i) = f_opt(i,1);
+%                     ni_mat(l,i) = ni(i,1);
+%                     mu_mat(l,i) = mu(i);
+%                     alpha_mat(l,i) = alpha(i);
+%                     pippo_mat(l,i) = pippo(i);
+                     cond_car_all = abs(Delta_load(i)/L_tot_vect(l))>=tol_workload_allocated;
+                    if cond_car_all
+                        iter_mat(l) = i;
+%                         MatrixRis(1:M,l) = L(i,:);
+%                         MatrixRis(M+1:M+M,l) = f_opt(i,:);
+%                         MatrixRis(M+M+1,l) = tot_cost(l,i);
+                    else
+                        fprintf('correct allocation for workload n: %d (%d) at iteration %d\n',l,L_tot_vect(l),i);
+                        iter_mat(l) = i;
+%                         MatrixRis(1:M,l) = L(i,:);
+%                         MatrixRis(M+1:M+M,l) = f_opt(i,:);
+%                         MatrixRis(M+M+1,l) = tot_cost(l,i);
+                        break;
+                    end
+                    i = i + 1;
+                end
+            end
+        end
+        if i == iterations + 1
+            i = i - 1;
+        end
+        
+               
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DESCRETE
+        
+        f_precedente=zeros(1,M);
+        f_successiva=zeros(1,M);
+        x1=zeros(1,M);
+        
+        Q=6; %numero frequenze discrete;
+        f_discrete=zeros(1,Q);
+        for conta=1:M
+            %f_discrete(conta,:)=[0:f_max(conta)./(Q-1):f_max(conta)];
+            f_discrete(conta,:)=[0, 5,50,70,90,105]; % Mhz or Mb/s
+        end
+        
+        f_opt_new=f_opt(iter_mat(l), :);
+        
+        for conta=1:M
+            
+            delta_f_discrete=f_discrete(conta,:)-f_opt_new(conta);
+            [ff,ind_ff]=min(abs(delta_f_discrete));
+            if ff==0
+                f_precedente(conta)=f_discrete(conta,ind_ff);
+                f_successiva(conta)=f_discrete(conta,ind_ff);
+                x1(conta)=1; %qualsiasi valore è indifferente
+            elseif ind_ff==1
+                f_precedente(conta)=f_discrete(conta,1);
+                f_successiva(conta)=f_discrete(conta,2);
+                x1(conta)=abs(f_opt_new(conta)-f_precedente(conta))./(f_successiva(conta)-f_precedente(conta));
+            elseif ind_ff==Q
+                f_precedente(conta)=f_discrete(conta,Q-1);
+                f_successiva(conta)=f_discrete(conta,Q);
+                x1(conta)=abs(f_opt_new(conta)-f_precedente(conta))./(f_successiva(conta)-f_precedente(conta));
+            elseif delta_f_discrete(ind_ff)>0
+                f_precedente(conta)=f_discrete(conta,ind_ff-1);
+                f_successiva(conta)=f_discrete(conta,ind_ff);
+                x1(conta)=abs(f_opt_new(conta)-f_precedente(conta))./(f_successiva(conta)-f_precedente(conta));
+            else
+                f_precedente(conta)=f_discrete(conta,ind_ff);
+                f_successiva(conta)=f_discrete(conta,ind_ff+1);
+                x1(conta)=abs(f_opt_new(conta)-f_precedente(conta))./(f_successiva(conta)-f_precedente(conta));
+            end
+            
+        end
+        
+        
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        
+         %speed_cost_all(M,l)=speed_cost(l,iter_mat(l));
+       % switch_cost_all(M,l)=switch_cost(l,iter_mat(l));
+        
+%         speed_cost_all(M,l)=sum(((((f_precedente/f_max).^2)).*omega.*E_max).*(1-x1)+((((f_successiva/f_max).^2)).*omega.*E_max).*x1);
+%         switch_cost_all(M,l)=sum((k_e.*(f_precedente-f_zero).^2).*(1-x1)+(k_e.*(f_successiva-f_zero).^2).*x1);
+%         
+%         channel_cost_all(M,l)=channel_cost(l,iter_mat(l));
+%         
+%         tot_cost_all(M,l)=speed_cost_all(M,l)+switch_cost_all(M,l)+channel_cost_all(M,l);
+        
+        
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        
+        f_zero = f_opt(i,:);
+        %1a Condizione
+        if i > iterations
+            i=iterations;
+        end
+        if i >= iterations
+            fprintf('!!!!!!!WARNING!!!!!!!prima condizione non soddifatta: NON CONVERGE!!!\n');
+            rule_counter(l)=1;
+            rule1_counter=rule1_counter+1;
+        end
+        %3a Condizione
+        for x = 1:M
+            cond3(1,x) = (((abs(ni(i-1,x)))<=epsilon)||((abs(L(i-1,x)-f_opt(i-1,x)*Delta(1,x)))<=epsilon));
+        end
+        if sum(cond3(1,:))< M
+            fprintf('!!!!!!!WARNING!!!!!!!terza condizione NON soddisfatta\n');
+            rule_counter(l)=3;
+            rule3_counter=rule3_counter+1;
+        end
+        %2a condizione
+        if L(i,:)>(f_opt(i,:) .* Delta(1,:));
+            fprintf('!!!!!!!WARNING!!!!!!!seconda condizione NON soddisfatta\n'); 
+            rule_counter(l)=2;
+            rule2_counter=rule2_counter+1;            
+        end
+    end
+    clear speed_cost switch_cost channel_cost tot_cost 
+%     if M==2
+%     mu_mat1 = mu_mat;
+%     iter_mat1=iter_mat;
+%     ni_mat1 = ni_mat;
+%     f_opt_mat1=f_opt_mat;
+%     L_mat1=L_mat;
+%     end
+end
+time = toc;
+fprintf('tempo: %d\n',time);
+
+NNN=1;
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%DESCRETE COMPARISONS
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Presimulation-definitions
+% speed_cost_all_WL=zeros(1,N);
+% switch_cost_all_WL=zeros(1,N);
+% channel_cost_all_WL=zeros(1,N);
+% tot_cost_all_WL=zeros(1,N);
+% 
+% speed_cost_all_VM=zeros(1,length(nVM));
+% switch_cost_all_VM=zeros(1,length(nVM));
+% channel_cost_all_VM=zeros(1,length(nVM));
+% tot_cost_all_VM=zeros(1,length(nVM));
+% 
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Presimulation-for
+% for j=1:N                   % workload
+%     speed_cost_all_WL(j)=mean(speed_cost_all(:,j));
+%     switch_cost_all_WL(j)=mean(switch_cost_all(:,j));
+%     channel_cost_all_WL(j)=mean(channel_cost_all(:,j));
+%     tot_cost_all_WL(j)=mean(tot_cost_all(:,j));
+% end
+% 
+% for j=1:length(nVM)             % VMs 
+%     speed_cost_all_VM(j)=mean(speed_cost_all(j,:));
+%     switch_cost_all_VM(j)=mean(switch_cost_all(j,:));
+%     channel_cost_all_VM(j)=mean(channel_cost_all(j,:));
+%     tot_cost_all_VM(j)=mean(tot_cost_all(j,:));
+% end
+% 
+% 
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% 
+% % figure(100)                     % non staisfaction numbers
+% % plot(,rule1_counter,'--+');
+% % xlabel('Workload');
+% % ylabel('$\overline{\mathcal{E}}_{CPU}$');
+% % grid on
+% 
+% 
+% 
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% figure(100)                     % speed
+% plot(1:1:N,speed_cost_all_WL(:));
+% xlabel('Workload');
+% ylabel('$\overline{\mathcal{E}}_{CPU}$');
+% grid on
+% 
+% figure(101)                     % switch
+% plot(1:1:N,switch_cost_all_WL(:));
+% xlabel('Workload');
+% ylabel('$\overline{\mathcal{E}}_{switch}$');
+% grid on
+% 
+% figure(102)                     % channel
+% plot(1:1:N,channel_cost_all_WL(:));
+% xlabel('Workload');
+% ylabel('$\overline{\mathcal{E}}^{net}$');
+% grid on
+% 
+% figure(103)                     % total
+% plot(1:1:N,tot_cost_all_WL(:));
+% xlabel('Workload');
+% ylabel('$\overline{\mathcal{E}}^{tot}$');
+% grid on
+% 
+% figure(104)                     % speed
+% plot(1:1:VMNo,speed_cost_all_VM(:));
+% xlabel('VM');
+% ylabel('$\overline{\mathcal{E}}_{CPU}$');
+% grid on
+% 
+% figure(105)                     % switch
+% plot(1:1:VMNo,switch_cost_all_VM(:));
+% xlabel('VM');
+% ylabel('$\overline{\mathcal{E}}_{switch}$');
+% grid on
+% 
+% figure(106)                     % channel
+% plot(1:1:VMNo,channel_cost_all_VM(:));
+% xlabel('VM');
+% ylabel('$\overline{\mathcal{E}}^{net}$');
+% grid on
+% 
+% figure(107)                     % total
+% plot(1:1:VMNo,tot_cost_all_VM(:));
+% xlabel('VM');
+% ylabel('$\overline{\mathcal{E}}^{tot}$');
+% grid on
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% mu_last=zeros (1,N);
+% ni_last=zeros (1,N);
+% alpha_last=zeros (1,N);
+% pippo_last=zeros (1,N);
+% 
+% for i=1:N
+% mu_last(i)=mu_mat(i,iter_mat(i));
+% ni_last(i)=ni_mat(i,iter_mat(i));
+% alpha_last(i)=alpha_mat(i,iter_mat(i));
+% pippo_last(i)=pippo_mat(i,iter_mat(i));
+% end
+% 
+% figure(1000)
+% plot(1:1:N,mu_last(:),'r');
+% title('mu');
+% xlabel('WL');
+% ylabel('$\mu$');
+% grid on
+% 
+% figure(1001)
+% plot(1:1:N,ni_last(:),'r');
+% title('ni');
+% xlabel('WL');
+% ylabel('$\ni$');
+% grid on
+% 
+% figure(1022)
+% plot(1:1:iter_mat1(50),mu_mat1(50,1:iter_mat1(50)),'r',1:1:iter_mat1(1000),mu_mat1(1000,1:iter_mat1(1000)),'b',1:1:iter_mat(50),mu_mat(50,1:iter_mat(50)),'r',1:1:iter_mat(1000),mu_mat(1000,1:iter_mat(1000)),'b');
+% title('$\mu$ in $\Delta=0.1$ (s), $\gamma=0.5$, $T_t=5$ (s), $\beta=0.15$');
+% legend('M=2, $\delta=0.1$, N=50','M=2, $\delta=0.1$, N=1000', 'M=100, $\delta=0.1$, N=50','M=100, $\delta=0.1$, N=1000')
+% xlabel('iterations');
+% ylabel('$\mu$');
+% grid on
+% 
+% figure(1023)
+% plot(1:1:iter_mat1(50),ni_mat1(50,1:iter_mat1(50)),'r',1:1:iter_mat1(1000),ni_mat1(1000,1:iter_mat1(1000)),'b',1:1:iter_mat(50),ni_mat(50,1:iter_mat(50)),'r',1:1:iter_mat(1000),ni_mat(1000,1:iter_mat(1000)),'b');
+% title('$\ni$ in $\Delta=0.1$ (s), $\gamma=0.5$, $T_t=5$ (s), $\beta=0.15$');
+% legend('M=2, N=50','M=2, N=1000', 'M=100, N=50','M=100, N=1000')
+% xlabel('iterations');
+% ylabel('$\ni$');
+% grid on
+% 
+% figure(1024)
+% plot(1:1:iter_mat1(50),f_opt_mat1(50,1:iter_mat1(50)),'r',1:1:iter_mat1(1000),f_opt_mat1(1000,1:iter_mat1(1000)),'b',1:1:iter_mat(50),f_opt_mat(50,1:iter_mat(50)),'r',1:1:iter_mat(1000),f_opt_mat(1000,1:iter_mat(1000)),'b');
+% title('$f_{opt}$ in $\Delta=0.1$ (s), $k_e=0.005$, $f_i^{max}=105$, $\gamma=0.5$, $T_t=5$ (s), $\beta=0.15$, $\overline{\mathcal{E}}_i^{max}=60$ $(J)$');
+% legend('M=2, N=50','M=2, N=1000', 'M=100, N=50','M=100, N=1000')
+% xlabel('iterations');
+% ylabel('$f_{opt}$');
+% grid on
+% 
+% figure(1025)
+% plot(1:1:iter_mat1(50),L_mat1(50,1:iter_mat1(50)),'r',1:1:iter_mat1(1000),L_mat1(1000,1:iter_mat1(1000)),'b',1:1:iter_mat(50),L_mat(50,1:iter_mat(50)),'r',1:1:iter_mat(1000),L_mat(1000,1:iter_mat(1000)),'b');
+% title('$L_{opt}$ in $\Delta=0.1$ (s), $k_e=0.005$, $f_i^{max}=105$, $\gamma=0.5$, $T_t=5$ (s), $\beta=0.15$, $\overline{\mathcal{E}}_i^{max}=60$ $(J)$');
+% legend('M=2, N=50','M=2, N=1000', 'M=100, N=50','M=100, N=1000')
+% xlabel('iterations');
+% ylabel('$L_{opt}$');
+% grid on
+
+% figure(1002)
+% plot(1:1:N,alpha_last(:),'r');
+% title('\alpha');
+% xlabel('WL');
+% ylabel('$\alpha$');
+% grid on
+% 
+% figure(1003)
+% plot(1:1:N,pippo_last(:),'r');
+% title('\alpha_{NEW}');
+% xlabel('WL');
+% ylabel('$\alpha_{NEW}$');
+% grid on
+% 
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+% figure()
+% plot(mu_mat(NNN,2:iter_mat(NNN)),'r');
+% title('mu');
+% 
+% figure()
+% plot(ni_mat(NNN,2:iter_mat(NNN)),'r');
+% title('ni');
+
+% figure()
+% plot(L_mat(NNN,2:iter_mat(NNN)),'r');
+% title('L');
+% 
+% figure()
+% plot(f_opt_mat(NNN,2:iter_mat(NNN)),'r');
+% title('f ottimo');
+
+% figure()
+% plot(alpha_mat(NNN,2:iter_mat(NNN)),'r');
+% title('alpha');
+% 
+% figure()
+% plot(pippo_mat(NNN,2:iter_mat(NNN)),'r');
+% title('pippo');
+
+% figure()
+% bar(MatrixRis(1:M,NNN));
+% title('allocazione carico (L) per il wl 1');
+% 
+% figure()
+% bar(MatrixRis(1+M:M+M,NNN));
+% title('allocazione frequenze (f_opt) per il wl 1');
+
+
diff --git a/ICC2015_Sourcecode/DVFS_type2_ICC.m b/ICC2015_Sourcecode/DVFS_type2_ICC.m
new file mode 100644
index 0000000000000000000000000000000000000000..eed6066f8c6d9c0b3cde6ef2ab8d95d05c7937d0
--- /dev/null
+++ b/ICC2015_Sourcecode/DVFS_type2_ICC.m
@@ -0,0 +1,312 @@
+clc
+clear all
+close all
+
+VM=[1:1:100];%numero Macchine Virtuali
+NumSimulazioni=length(VM);
+
+%IDEAL
+costo_speed_VM=zeros(1,NumSimulazioni);
+costo_switch_VM=zeros(1,NumSimulazioni);
+costo_channel_VM=zeros(1,NumSimulazioni);
+costo_tot_VM=zeros(1,NumSimulazioni);
+%NICOLA (COMNET)
+speed_cost_all_VM=zeros(1,NumSimulazioni);
+switch_cost_all_VM=zeros(1,NumSimulazioni);
+channel_cost_all_VM=zeros(1,NumSimulazioni);
+tot_cost_all_VM=zeros(1,NumSimulazioni);
+
+
+f_VM_ultima=zeros(1,NumSimulazioni);
+
+
+%inizio le simulazioni
+P_net_completo=[1:0.25:10000];
+%P_net_completo=zeros(1,100);
+
+
+%Jobs=6 + (10-6).*rand(1,1000); %carico uniformemente distribuito
+%Jobs=10.*ones(1,10);  %carico deterministico
+%Jobs=ciel(abs(10.*Nrand(0,1)));
+%Jobs=random('Poisson',1:NumSimulazioni,1,NumSimulazioni);
+%Jobs=ceil(abs(random('Normal',0,4,1,10)));
+Jobs=load('workloads8PMR125No2000.mat','L_tot_vect');
+Jobs=Jobs.L_tot_vect;
+for k=1:NumSimulazioni
+    tic
+    k
+    M=VM(k);
+    costo_speed=0;
+    costo_switch=0;
+    costo_channel=0;
+    costo_tot=0;
+    
+    costo_speed_Nicola=0;
+    costo_switch_Nicola=0;
+    costo_channel_Nicola=0;
+    costo_tot_Nicola=0;
+    rate_ultima_macchina=0;
+    
+    
+    k_e=0.005;
+    C_max=15; %(Mb/s)
+    T_tot=5; %(s)
+    
+    W = ones(1,M);                  %band of each channel (MHz)
+    Zeta = 0.5 * ones(1,M);         %channel parameters (mW)
+    
+    
+    f_max=105.*ones(1,M);  %(Mb/s)
+    f_zero=zeros(1,M);  % (Mb/s)
+    %f_zero=0.2.*f_max;
+    %f_zero=f_max;
+    P_net=P_net_completo(1:M); %(mW)
+    %P_net=1.*ones(1,M);   %NB.! ORDINARE LE MACCHINE VIRTUALI IN MODO DA AVERE LE POTENZE IN ORDINE CRESCENTE
+    %(unicamente perchè semplifica la scritura del
+    %software - controindicazione: il codicie funziona bene indipendentemente)
+    Th=2.*P_net./C_max;   %Soglia di ibernazione macchina; (tutte le VM per cui mu_opt<=Th vengono ibernate con coefficiente di
+    % ibernazione alpha_zero)
+    
+    E_max=60.*ones(1,M);  % (mJ)
+    omega=1.*ones(1,M);
+    Delta=0.1.*ones(1,M); % (s)
+    Delta_max=max(Delta);
+    L_b=zeros(1,M);  % NB. il codice funziona bene per L_b=0, altrimenti aggiungere modifica su f_opt.
+    
+    temp=2.*k_e+(2.*E_max.*omega./(f_max.^2));
+    alpha_zero=2.*k_e./temp;
+    alpha_mu=Delta./temp;
+    
+speed_cost_all_WL=zeros(1,length(Jobs));
+switch_cost_all_WL=zeros(1,length(Jobs));
+channel_cost_all_WL=zeros(1,length(Jobs));
+tot_cost_all_WL=zeros(1,length(Jobs));
+
+    
+    
+    for num_job=1:length(Jobs)
+        L_tot=Jobs(num_job);
+        
+        
+        
+        
+        %CHECK FEASIBILITY - controllare funzionamento operatori %%%%%%%%%%%%%%%%%%
+        
+        condizione_feas_carico=(sum(f_max.*Delta-L_b)>=L_tot);
+        condizione_feas_tempo=(L_tot<=C_max.*(T_tot-Delta_max)./2);
+        condizione_feas_back=min(Delta.*f_max>=L_b);
+        
+        
+        feasibility=condizione_feas_carico && condizione_feas_tempo && condizione_feas_back;
+        
+        if ~feasibility
+            'VM='
+            M
+            error('Problem unfeasible')
+        else
+            %'Problem Feasibile!'
+        end
+        
+        
+        
+        
+        
+        %SOLUZIONE EQUAZIONE:   sum(L_opt)=L_tot;
+        
+        %NB. IL PUNTO DI ATTRAVERSAMENTO DELLO ZERO POTREBBE NON ESISTERE. IN QUEL
+        %CASO mu_opt E' IL PUNTO DI GRADINO FRA LA ZONA POSITIVA E QUELLA NEGATIVA
+        %E delta_mu (CHE E' IL VALORE DELL'EQUAZIONE IN mu) E' DIVERSO DA ZERO.
+        
+        [mu,delta_mu]= Mu_opt_bisezione(alpha_zero,alpha_mu,P_net,C_max,f_zero,f_max,Delta,L_b,L_tot);
+        tol_mu=10^(-2);%tolleranza per decisione macchine in stato limite
+        tol_carico_allocato=10^(-2);
+        if ((abs(delta_mu)./L_tot)<tol_carico_allocato)%tutto il carico è allocato con un errore relativo inferiore a tol_allocazione
+            caso_limite=0;
+            %'Caso standard: nessuna VM in stato limite'
+        else
+            caso_limite=1;
+            %'Una o più VM è in stato limite'
+            VM_limite=find(abs(Th-mu)<tol_mu);
+            if length(VM_limite)==0
+                'VM='
+                M
+                error('nessuna VM in stato limite e carico complessivo L_tot non allocato correttamente')
+            end
+            
+        end
+        
+        
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        %%%%%%%%%%%%%%%%%%%%%%%%%   scheduler ottimo  %%%%%%%%%%%%%%%%%
+        
+        
+        
+        
+        f_mu=max(mu-2.*P_net./C_max,0);
+        f_star=alpha_zero.*f_zero+alpha_mu.*f_mu;
+        f_opt=max(0,min(f_star,f_max)); %N.B. --> NEL CASO GENERALE DI L_b>0 IL MAX VA FATTO RISPETTO A L_b/Delta, NON RISPETTO A 0.
+        
+        rate_ultima_macchina=rate_ultima_macchina+f_opt(M);
+        
+        
+        
+        canali_attivi=f_mu>0;
+        L_opt=canali_attivi.*(f_opt.*Delta-L_b);
+        
+        %riallocazione delle eventuali VM sulla soglia (stato limite)
+        
+        if caso_limite
+            L_opt(VM_limite)=0;
+            L_allocato=sum(L_opt);
+            L_residuo=L_tot-L_allocato;
+            for k2=1:length(VM_limite);
+                L_opt(VM_limite(k2))=min(L_residuo,f_opt(VM_limite(k2)).*Delta(VM_limite(k2))-L_b(VM_limite(k2)));
+                L_residuo=L_residuo-L_opt(VM_limite(k2));
+            end
+            canali_attivi=(L_opt>0);
+        end
+        
+        
+        
+        
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        %%%%%%%%%%  costo  %%%%%%%%%%%%%%%%%%
+        
+        Delta_carico=sum(L_opt)-L_tot;
+        if ((abs(Delta_carico)./L_tot)>=tol_carico_allocato)
+            'Errore: Carico non correttamente allocato'
+        end
+        
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%descrete
+        f_precedente=zeros(1,M);
+        f_successiva=zeros(1,M);
+        x1=zeros(1,M);
+        
+        Q=6; %numero frequenze discrete;
+        f_discrete=zeros(M,Q);
+        for conta=1:M
+            %f_discrete(conta,:)=[0:f_max(conta)./(Q-1):f_max(conta)];
+            f_discrete(conta,:)=[0, 5,50,70,90,105]; % Mhz or Mb/s
+        end
+        
+        f_opt_new=f_opt;
+        
+        for conta=1:M
+            
+            delta_f_discrete=f_discrete(conta,:)-f_opt_new(conta);
+            [ff,ind_ff]=min(abs(delta_f_discrete));
+            if ff==0
+                f_precedente(conta)=f_discrete(conta,ind_ff);
+                f_successiva(conta)=f_discrete(conta,ind_ff);
+                x1(conta)=1; %qualsiasi valore è indifferente
+            elseif ind_ff==1
+                f_precedente(conta)=f_discrete(conta,1);
+                f_successiva(conta)=f_discrete(conta,2);
+                x1(conta)=abs(f_opt_new(conta)-f_precedente(conta))./(f_successiva(conta)-f_precedente(conta));
+            elseif ind_ff==Q
+                f_precedente(conta)=f_discrete(conta,Q-1);
+                f_successiva(conta)=f_discrete(conta,Q);
+                x1(conta)=abs(f_opt_new(conta)-f_precedente(conta))./(f_successiva(conta)-f_precedente(conta));
+            elseif delta_f_discrete(ind_ff)>0
+                f_precedente(conta)=f_discrete(conta,ind_ff-1);
+                f_successiva(conta)=f_discrete(conta,ind_ff);
+                x1(conta)=abs(f_opt_new(conta)-f_precedente(conta))./(f_successiva(conta)-f_precedente(conta));
+            else
+                f_precedente(conta)=f_discrete(conta,ind_ff);
+                f_successiva(conta)=f_discrete(conta,ind_ff+1);
+                x1(conta)=abs(f_opt_new(conta)-f_precedente(conta))./(f_successiva(conta)-f_precedente(conta));
+            end
+            
+        end
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        
+        %IDEAL
+        costo_speed=costo_speed+sum(((f_opt./f_max).^2).*omega.*E_max);
+        costo_switch=costo_switch+sum(k_e.*(f_opt-f_zero).^2);
+        costo_channel=costo_channel+sum(2.*P_net.*L_opt./C_max);
+        
+        %NICOLA
+        
+        costo_speed_Nicola=costo_speed_Nicola+sum(((((f_precedente/f_max).^2)).*omega.*E_max).*(1-x1)+((((f_successiva/f_max).^2)).*omega.*E_max).*x1);
+        costo_switch_Nicola=costo_switch_Nicola+sum((k_e.*(f_precedente-f_zero).^2).*(1-x1)+(k_e.*(f_successiva-f_zero).^2).*x1);
+        costo_channel_Nicola=costo_channel_Nicola+((T_tot-Delta(1)).*sum(Zeta .* (2.^(2*L_opt./((T_tot-Delta).*W))-1)));
+        %costo_channel_Nicola=costo_channel_Nicola+sum(2.*P_net.*L_opt./C_max);
+        
+        speed_cost_all_WL(num_job)=sum(((((f_precedente/f_max).^2)).*omega.*E_max).*(1-x1)+((((f_successiva/f_max).^2)).*omega.*E_max).*x1);
+        switch_cost_all_WL(num_job)=sum((k_e.*(f_precedente-f_zero).^2).*(1-x1)+(k_e.*(f_successiva-f_zero).^2).*x1);
+        channel_cost_all_WL(num_job)=((T_tot-Delta(1)).*sum(Zeta .* (2.^(2*L_opt./((T_tot-Delta).*W))-1)));
+        tot_cost_all_WL(num_job)=speed_cost_all_WL(num_job)+switch_cost_all_WL(num_job)+channel_cost_all_WL(num_job);
+
+        
+        
+        f_zero=f_opt;
+        
+    end
+    
+    
+    % risultati e deallocazione
+    
+    
+    k
+    %IDEAL
+    costo_speed_VM(k)=costo_speed./length(Jobs);
+    costo_switch_VM(k)=costo_switch./length(Jobs);
+    costo_channel_VM(k)=costo_channel./length(Jobs);
+    costo_tot_VM(k)=costo_speed_VM(k)+costo_switch_VM(k)+costo_channel_VM(k);
+    %NICOLA (COMNET)
+    speed_cost_all_VM(k)=costo_speed_Nicola./length(Jobs);
+    switch_cost_all_VM(k)=costo_switch_Nicola./length(Jobs);
+    channel_cost_all_VM(k)=costo_channel_Nicola./length(Jobs);
+    tot_cost_all_VM(k)=speed_cost_all_VM(k)+switch_cost_all_VM(k)+channel_cost_all_VM(k);
+    f_VM_ultima(k)=rate_ultima_macchina./length(Jobs);
+    
+    
+    clear f_max f_zero P_net Th  E_max omega Delta L_b temp alpha_zero alpha_mu
+    
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    %%%%%%%%%   istogrammi  %%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    %  figure(k);
+    %  bar(f_opt),xlabel('VM');
+    %  ylabel('Speed f');
+    %  grid on;
+    
+    toc
+end
+
+
+figure(1)
+plot(VM,costo_speed_VM,'--o',VM,speed_cost_all_VM,'--*'),xlabel('VM'),ylabel('Speed Cost');
+legend('IDEAL(no DVFS)_type2', 'DVFS_type2');
+figure(2)
+plot(VM,costo_switch_VM,'--o',VM,switch_cost_all_VM,'--*'),xlabel('VM'),ylabel('Switch Cost');
+legend('IDEAL(no DVFS)_type2', 'DVFS_type2');
+figure(3)
+plot(VM,costo_channel_VM,'--o',VM,channel_cost_all_VM,'--*'),xlabel('VM'),ylabel('Net Cost');
+legend('IDEAL(no DVFS)_type2', 'DVFS_type2');
+figure(4)
+plot(VM,costo_tot_VM,'--o',VM,costo_tot_VM,'--*'),xlabel('VM'),ylabel('Overall Cost');
+legend('IDEAL(no DVFS)_type2', 'DVFS_type2');
+
+
+figure(5)
+plot(1:1:length(Jobs),speed_cost_all_WL,'--o'),xlabel('#WL'),ylabel('Speed Cost');
+legend('DVFS_type2');
+figure(6)
+plot(1:1:length(Jobs),switch_cost_all_WL,'--*'),xlabel('#WL'),ylabel('Switch Cost');
+legend('DVFS_type2');
+figure(7)
+plot(1:1:length(Jobs),channel_cost_all_WL,'--*'),xlabel('#WL'),ylabel('Net Cost');
+legend('DVFS_type2');
+figure(8)
+plot(1:1:length(Jobs),tot_cost_all_WL,'--*'),xlabel('#WL'),ylabel('Overall Cost');
+legend('DVFS_type2');
+
+
+
+figure(9)
+plot(VM,f_VM_ultima,'--o'),xlabel('VM'),ylabel('rate last VM');
+
+
+
diff --git a/ICC2015_Sourcecode/Lyapunov5_Dynamic_V_Fixed_Switch.m b/ICC2015_Sourcecode/Lyapunov5_Dynamic_V_Fixed_Switch.m
new file mode 100644
index 0000000000000000000000000000000000000000..2880a16c9afe04bd3854b13654450f7247cf9c8c
--- /dev/null
+++ b/ICC2015_Sourcecode/Lyapunov5_Dynamic_V_Fixed_Switch.m
@@ -0,0 +1,507 @@
+%variables inizialization utilizzo flessibile delle code , non ha vincoli
+%di tempo di servizio.% CHE THROUPUT GARANTISCONO? E COSTO AL VARIARE v
+clc
+clear all
+close all
+cc=0;
+N = 1;                                                                    %applications pag 479
+M = 100;                                                                   %Servers pag 479
+                                                                  
+TT = 2000;
+T=2000;
+k_e=0.005; %(J/(MHz)^2)
+%Lamda = 2000.*ones(N,TT);                                                  %Number of slot
+%A = 2.*Lamda.*rand(N,TT);                                                  %Arrival job uniformed distribuited pag 481
+%A=4000.*rand(N,TT);                                                        %PMR = 1
+%A = 3000 +(4000-2000).*rand(N,TT);                                         %PMR = 1.25
+%A = 2000 +(6000-2000).*rand(N,TT) ;                                        %PMR = 1.5
+% A = 1000 +(7000-1000).*rand(N,TT);                                        %PMR = 1.75
+%A = 8000.*rand(N,TT);                                                        %PMR = 2
+%save('Avalue4000FIXPMR2.mat','A');
+
+
+%L_tot_vect = 6+ round((10-6) * rand(N,TT));                         % PMR=1.25, L_tot=8+-2
+%save('workloads8PMR125No2000.mat', 'L_tot_vect');                  %loading 100 fixed random workloads
+A=load('workloads8PMR125No2000.mat','L_tot_vect');
+A=A.L_tot_vect;
+
+
+F = [5,50,70,90,105];                                         %( Gbyte/s)frequecny range for every VM allocated onto the j-server
+I = ones(M,TT);                                                            %resource decision allocation
+P_cpu1 = zeros(N,M,length(F));                                              %power consuption for each VM onto each j server
+P_min = 10.*ones(N,M,length(F));                                          %minimum power for the considered VM allocated onto the j-server
+P_max = 60.*ones(N,M,length(F));                                          %maximum power for the considered VM allocated onto the j-server
+mu1=zeros(N,M,length(F)) ;
+MU=zeros(N,M,TT);                                                          %service rate (request/slots)
+ww = zeros(1,N);
+AU = zeros(N,M);
+AAU = zeros(1,N);
+SAU = zeros(1,M);
+AAS = zeros(1,N);
+Avg_ADELAY = zeros(1,N);
+avg_MU = zeros(N,M);
+avg_AMU = zeros(1,N);
+avg_SMU = zeros(1,M);
+Avg_SDELAY = zeros(1,M);
+Avg_SDelay=zeros(1,M);
+
+KU=  zeros(N,M);
+Delay=zeros(1,N);
+Delay2=zeros(N,M);
+
+
+alpha = ones(1,N);                                                         %throughput utility weights
+%V1 = [1:100:1000];                                                         %control parameter for DCA
+V = 100;
+%V=length(V1)
+beta = 1;                                                                  % non-negative normalizing weight
+W = zeros(N,TT);                                                           %Buffer Dimension
+K = zeros(N,M,TT);                                                         % ==R(i,j) number of requests for application i that are routed from the R(i) router Buffer to the j-server in slot t.
+U = zeros(N,M,TT);                                                          %queing dynamics for the request of application i at server j
+P = zeros(M,TT);                                                           % power consumed by each server in each time slot for s
+Switch=zeros(M,TT);                                                            % power consumed by each server in each time slot for s
+SW= zeros(1,M);                                                            %the time average expected power consuption of server j
+r = zeros(1,N);                                                            %average expected rate of admitted request for apllication i
+e = zeros(1,M);                                                            %the time average expected power consuption of server j
+R = zeros(N,TT);
+Freq = zeros(N,M,length(F));
+Freqtemp = zeros(M,TT);
+temp = zeros(N,M,length(F));
+active_servers = zeros(N+1,M,T);
+M1=0;
+M2=0;
+M3=0;
+l=0;
+on_servers_list=zeros(1,M);
+on_servers = zeros(N+1,M);                                                 % on-server in each different frame T
+off_servers = zeros(N+1,M);                                                % off-server in each different frame T
+sleeping_servers= zeros(N+1,M);                                            % Hybernated-server in each different frame T
+% a = rand(1,M);                                                             % applications variable indicator
+% save('avalue.mat','a');
+a=load('avalue.mat','a');
+a=a.a;
+a(:)=(a(:) >= 0.4);
+for j=1:M
+    if (a(j)==1)
+        on_servers(1,j)=j;
+        l=l+1;
+        on_servers_list(j)=j;
+    end
+end
+
+a_temp = zeros(1,M);                                                        %  temp applications variable indicator
+%a_temp(:)=(a_temp(:) >= 0.6);
+
+for i = 1:N
+    for j = 1:M
+        Freq(i,j,:)= F(:);
+    end
+end
+
+for i=1:N
+    for j=1:M
+        mu1(i,j,:) = (8*10^(-3).*F(:)) +0.76;                  % servers queue for each application service rate (5,0.8)-----(105,1.6)
+        P_cpu1(i,j,:)=P_min(i,j,:) + (50*10^(-4)).*((Freq(i,j,:)-5).^2);     % CPU power of each virtual machine among each j-server pag 480
+    end
+end
+
+F1= [0, 5,50,70,90,105];                                             % (Gbit/s)frequecny range for every VM allocated onto the j-server
+P_cpu2 = zeros(N,M,length(F1));                                            % power consuption for each VM onto each j server
+P_min = 10.*ones(N,M,length(F1));                                          %minimum power for the considered VM allocated onto the j-server
+P_max = 60.*ones(N,M,length(F1));                                          %maximum power for the considered VM allocated onto the j-server
+mu2 = zeros(N,M,length(F1))  ;
+Freq1 = zeros(N,M,length(F1));
+temp1 = zeros(N,M,length(F1));
+for i = 1:N
+    for j = 1:M
+        Freq1(i,j,:)= F1(:);
+    end
+end
+
+for i=1:N
+    for j=1:M
+        mu2(i,j,2:1:length(F1)) = (8*10^(-3).*F1(2:1:length(F1)))+0.76; % servers queue for each application service rate (5,0.8)-----(105,1.6)
+        mu2(i,j,1) = 0;
+        P_cpu2(i,j,2:1:length(F1))=P_min(i,j,2:1:length(F1)) + (50*10^(-4)).*((Freq1(i,j,2:1:length(F1))-5).^2);  % CPU power of each virtual machine among each j-server pag 480
+        P_cpu2(i,j,1)=0;
+        
+    end
+end
+i=0;
+j=0;
+UU=zeros(N,M,TT);
+UU1=zeros(N,M,TT);
+B = [1:1:TT/T];
+L_off=zeros(1,length(B));
+MM1=zeros(1,length(B));
+MM2=zeros(1,length(B));
+MM3=zeros(1,length(B));
+bb=0;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%main body code
+for t = 1 : TT
+    if (mod(t,T)~= 0)
+        UU=zeros(N,M,TT);
+        for jtest=1:M
+            if (a(jtest)==1)
+                UU(:,jtest,t)=U(:,jtest,t);
+            else
+                UU(:,jtest,t)=nan;
+            end
+        end
+        %if (t>1)
+        %UU(~UU)=nan;
+        %end
+        for i= 1 : N
+            
+            %
+            %j
+            t
+            if(W(i,t) > V.*alpha(i))
+                R(i,t) = 0;
+                
+            else
+                R(i,t) = A(i,t);
+            end
+            
+            %list U of active (ON) servers.
+            
+            
+            [U1, indexu1]=min(UU(i,:,t));% page 483 routing
+            if (W(i,t) > U1) || (isnan(U1))
+                K(i,indexu1,t) = W(i,t);
+            else
+                K(i,:,t)=0;
+            end
+            
+            UU(:,indexu1,t)=nan;
+            
+            %                 [U1, ind]=min(U(i,:,t));% page 483 routing
+            %                 if (W(i,t) > U1)
+            %                     for ikl=1:M
+            %                         for mmmm=1:N
+            %                             K(mmmm,ikl,t)=K(mmmm,ikl,t-1);
+            %                         end
+            %                     end
+            %                     K(i,ind,t) = W(i,t);
+            %
+            %                 else
+            %                     K(i,:,t)=0;
+            %                 end
+            
+            
+            %a(i,j) = (a(i,j) > 0.5); % garantee to set more than half to be -on
+            % fix number of active server, they are not M , they are
+            % obviously less
+            
+            %if sum(a(i,j)*K(i,j,t)<=W(i,t))
+            %if ((W(i,t)-sum(K(i,:,t)))>=0)
+            
+            % W(i,t+1) = W(i,t) - sum(a(i,j).*K(i,:,t)) + R(i,t);
+            W(i,t+1) = W(i,t) - (a*K(i,:,t)') + R(i,t);                %% Backlog Queue 1    (2) pag 481 consider only active server
+            % else
+            %    W(i,t+1)=R(i,t);
+            % end
+            
+        end % for N-application
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% after routing inside each server
+        for j= 1 : M
+            if (a(j)==1)
+                [U_max , index_U_max]=max(U(:,j,t));                   % 0=<index_U_max<=N
+                
+                for zzz=1:length(F)
+                    temp(index_U_max,j,zzz) = U_max.*(min(mu1(index_U_max,j,zzz),U_max)) - V.*beta.*(P_cpu1(index_U_max,j,zzz));                  % Resource allocation pag 483
+                    
+                end
+                [value , index]=max(temp(index_U_max,j,:));            % 1=<index<=length(F)
+                
+                MU(index_U_max,j,t) = min(mu1(index_U_max,j,index),U_max);
+                P(j,t)= P_cpu1(index_U_max,j,index);
+                if (t==1)
+                    Switch(j,t)=k_e*(F(index))^2;
+                    Freqtemp(j,2)=F(index);
+                else
+                    Switch(j,t)=k_e*((F(index)-Freqtemp(j,t-1)))^2;
+                    Freqtemp(j,t)=F(index);
+                end
+                for i= 1 : N
+                    if (i==index_U_max)
+                        U(index_U_max,j,t+1) = max((U(index_U_max,j,t) - MU(index_U_max,j,t)),0) + K(index_U_max,j,t);      %% Backlog Queue 2    (5) pag 481
+                    else
+                        U(i,j,t+1) = max(U(i,j,t),0) + K(i,j,t);           % Backlog Queue 2    (5) pag 481
+                    end
+                end
+                % U(index_U_max,j,t+1) = max(U(index_U_max,j,t) - MU(index_U_max,j,t),0) + a(i,j).*K(index_U_max,j,t);      %% Backlog Queue 2    (5) pag 481
+            end
+        end                                                                % for each M-server
+        
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%check points
+    else
+        bb=bb+1;
+        UU1=zeros(N,M,TT);
+        for jtest1=1:M
+            if (a(jtest1)==1)
+                UU1(:,jtest1,t)=U(:,jtest1,t);
+            else
+                UU1(:,jtest1,t)=nan;
+            end
+        end
+        %if (t>1)
+        %UU(~UU)=nan;
+        %end
+        for i= 1 : N
+            
+            if(W(i,t) > V.*alpha(i)) 
+                R(i,t) = 0;
+                
+            else
+                R(i,t) = A(i,t);
+            end
+            
+            %list U of active (ON) servers.
+            
+            [U2, indexu2]=min(UU1(i,:,t));% page 483 routing
+            
+            %[U2, indexu2]=min(U(i,:,t));                                   % page 483 routing
+            if (W(i,t) > U2)
+                K(i,indexu2,t) = a(indexu2)*W(i,t);
+            else
+                K(i,:,t)=0;
+            end
+            UU1(:,indexu2,t)=nan;
+            %a(i,j) = (a(i,j) > 0.5);                                      % garantee to set more than half to be -on
+            % fix number of active server, they are not M , they are
+            % obviously less
+            % if ((W(i,t)-sum(K(i,:,t)))>=0)
+            
+            % W(i,t+1) = W(i,t) - sum(a(i,j).*K(i,:,t)) + R(i,t);
+            W(i,t+1) = W(i,t) - (a*K(i,:,t)') + R(i,t);                %% Backlog Queue 1    (2) pag 481 consider only active server
+            % else
+            %W(i,t+1)=R(i,t);
+        end
+        
+        
+        for j= 1 : M
+            a_temp(j)=a(j);
+            [U_max1, index_U_max1]=max(U(:,j,t));                           % 0=<index_U_max<=N
+            
+            for zzz=1:length(F1)
+                temp1(index_U_max1,j,zzz) = U_max1.*(min(mu2(index_U_max1,j,zzz),U_max1)) - V.*beta.*(P_cpu2(index_U_max1,j,zzz));                  % Resource allocation pag 483
+                
+            end
+            [value2, index2]=max(temp1(index_U_max1,j,:));                % 1=<index<=length(F)
+            %%%%%%%%%%%%%%%checking
+            if (((min(mu2(index_U_max1,j,index2),U_max1))==0) && (P_cpu2(index_U_max1,j,index2))==0)
+                cc=cc+1;
+                a(j)=0;
+                if (t==1)
+                    Switch(j,t)=k_e*(F1(index2))^2;
+                    Freqtemp(j,2)=F1(index2);
+                else
+                    Switch(j,t)=k_e*(F1(index2)-Freqtemp(j,t-1))^2;
+                    Freqtemp(j,t)=F1(index2);
+                end
+                MU(index_U_max1,j,t)=0;
+                P(j,t)=0;
+            else
+                a(j)=1;
+                MU(index_U_max1,j,t) = min(mu2(index_U_max1,j,index2),U_max1);
+                P(j,t)= P_cpu2(index_U_max1,j,index2);
+                if (t==1)
+                    Switch(j,t)=k_e*(F1(index2))^2;
+                    Freqtemp(j,2)=F1(index2);
+                else
+                    Switch(j,t)=k_e*(F1(index2)-Freqtemp(j,t-1))^2;
+                    Freqtemp(j,t)=F1(index2);
+                end
+            end
+            %%%%%%%%%%%%
+            
+            
+            for i= 1 : N
+                if (i==index_U_max1)
+                    U(index_U_max1,j,t+1) = max(U(index_U_max1,j,t) - MU(index_U_max1,j,t),0) + K(index_U_max1,j,t);    % Backlog Queue 2    (5) pag 481
+                else
+                    U(i,j,t+1) = max(U(i,j,t),0) + K(i,j,t);           % Backlog Queue 2    (5) pag 481
+                end
+            end
+            % U(index_U_max,j,t+1) = max(U(index_U_max,j,t) - MU(index_U_max,j,t),0) + a(i,j).*K(index_U_max,j,t);      %% Backlog Queue 2    (5) pag 481
+            % if
+        end                                                                % for each M-server
+        
+        
+        %*************************VM_MIGRATION******************
+   on_servers = zeros(N+1,M);                                                 % on-server in each different frame T
+   off_servers = zeros(N+1,M);                                                % off-server in each different frame T
+   sleeping_servers= zeros(N+1,M);
+        
+   for jj = 1 : M
+       for ii = 1 : N
+           if ((a_temp(jj)==0) && (a(jj)==1)) || (((a_temp(jj)==1) && (a(jj)==1)))
+               on_servers(1,jj) = jj; %always ON servers + OFF--->> ON servers
+               on_servers(ii+1,jj) =  U(ii,jj,t);
+               %M1 = M1 + 1;
+           else if((a_temp(jj)==1) && (a(jj)==0))
+                   off_servers(1,jj) = jj;
+                   off_servers(ii+1,jj) = U(ii,jj,t);
+                   %M2 = M2 + 1;
+               
+           else 
+                   sleeping_servers(1,jj)=jj;  %the servers which are never ON, ALWAYS OFF
+                   sleeping_servers(ii+1,jj) =  U(ii,jj,t);
+                   %M3=M3+1;
+               end
+               
+           end
+           
+       end
+   end
+                
+         
+        
+        on_servers = deletezero(on_servers);
+        off_servers= deletezero(off_servers);
+        sleeping_servers=deletezero(sleeping_servers);
+        [x1, x2]=size(on_servers);
+        M1=x2;
+        [y1, y2]=size(off_servers);
+        M2=y2;
+        [z1, z2]=size(sleeping_servers);
+        M3=z2;
+        
+        L_off(bb)=sum(sum(off_servers(2:N+1,:))) ;                                  %  the total migrated backlogs of all aplication in bits from the off servers
+        for jk = 1 : M2
+            for iii = 1 : N
+                
+                
+                [a_tempON, indexON] = min(on_servers(iii+1,:));
+                tempOnserver=off_servers(iii+1,jk);
+                on_servers(iii+1,indexON) = on_servers(iii+1,indexON)+tempOnserver;
+                off_servers(iii+1,jk)=0;
+                
+                
+                
+            end
+        end
+        
+        
+        a=zeros(1,M);
+        for dom=1:M1
+            if (on_servers(1,dom)>=1)
+                a(on_servers(1,dom))=1;
+                for i= 1 : N
+                    U(i,on_servers(1,dom),t+1)=on_servers(i+1,dom);
+                end
+            end
+        end
+        for dom=1:M2
+            a(off_servers(1,dom))=0;
+            for i= 1 : N
+                U(i,off_servers(1,dom),t+1)=0;
+            end
+        end
+        %         for dom=1:M3
+        %             a(sleeping_servers(1,dom))=0;
+        %         end
+        
+        
+        MM1(bb)=M1;
+        MM2(bb)=M2;
+        MM3(bb)=M3;
+        M1=0;
+        M2=0;
+        M3=0;
+    end % t=nT
+    
+    
+end % biggggg for TT
+
+ for ii= 1:N
+     r(ii) = mean(R(ii,:));                                              % throughput 4 each ii
+     ww(ii) = mean(W(ii,:));                                           % average occupation for each applications before router
+     for jj = 1:M
+         e(jj) =   mean(P(jj,:));                                      %pag482 (7)
+         SW(jj)= mean (Switch(jj,:));
+         AU(ii,jj) = mean(U(ii,jj,:));                                     % average occupation of service rate for application i in sererver j in block 2
+         KU(ii,jj) = sum(K(ii,jj,:));                                     % average occupation of backlogs (in queue) of block 2 (queues in servers)
+         avg_MU(ii,jj)= mean(MU(ii,jj,:));
+     end
+     
+end
+
+KUU=0;
+
+    for i= 1:N
+        for j= 1:M
+            KUU=KU(i,j)+KUU;
+        end
+    end
+
+KUU=KUU/(M*N*TT);
+for i= 1:N
+Delay(i)=ww(i)/ r(i);  % delay takes for each application queue before routing
+
+end
+
+
+Avg_Delay1=mean(Delay(:)); % average delay of BLOCK 1 (before routing to the servers)
+
+
+
+for i= 1:N
+for j= 1:M
+%Delay2(i,j)=AU(i,j)./ KU(i,j);  % delay takes for application i in server j in block 2
+Delay2(i,j)=AU(i,j)./ KUU;  % delay takes for application i in server j in block 2
+end
+end
+
+
+
+for j= 1:M
+Avg_SDelay(j)=mean(Delay2(:,j)); % average delay takes for server j in block 2
+end
+
+
+Avg_Delay2=mean(Avg_SDelay(:));  % average dalay of BLOCK 2 (after routing to the servers and after servers )
+
+
+
+TOT_DELAY=Avg_Delay1+Avg_Delay2;  %average total delay of WHOLE SYSTEM
+UUTILITY = alpha(1).*sum(r(:)) - beta.*sum(e(:)); % UTILITY OF THE SYSTEM
+
+
+
+%AVG_UTILITY = mean(UUTILITY(:));   % total utility of the SYSTEM
+
+ 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
+    
+    
+%     
+%     
+%     www=mean(ww(:));                                                           % average occpuation of all applications
+% for i = 1 : N
+% AAU(i)=mean(AU(i,:));                                                       % average occupation of each backlogs application in the system
+% avg_AMU(i)=mean(avg_MU(i,:));    
+% end
+% 
+% for j = 1 : M
+% SAU(j)=mean(AU(:,j));                                                       % average occupation of total backlogs applications in each  server j  
+% avg_SMU(jj)=mean(avg_MU(:,j));  
+% end
+% 
+% TAU=mean(AAU(:));                                                          % average occupation for : server for : applications
+% Tavg_MU=mean(avg_AMU(:));
+% 
+%  
+% AS = www + TAU ;                                                           % total average occupations of the system
+% 
+% Avg_DELAY =  AS/mean(r(:));                                                % average delay of the system
+% 
+% AAS(:) = ww(:) + AAU(:) ;                                                  % total average occupations of the system for each application
+% 
+% Avg_ADELAY(:) = AAS(:)./r(:);                                              % average delay of the system for each application
+% 
+% Avg_SDELAY(:) = SAU(:)./avg_SMU(:);       
+% 
+% UTILITY = alpha(1)*sum(r(:)) - beta*sum(e(:));
diff --git a/ICC2015_Sourcecode/Lyapunov5_Static_V_Fixed_Switch.m b/ICC2015_Sourcecode/Lyapunov5_Static_V_Fixed_Switch.m
new file mode 100644
index 0000000000000000000000000000000000000000..8d807cf024450ca5e619d9dfe2bd4af88de69aa1
--- /dev/null
+++ b/ICC2015_Sourcecode/Lyapunov5_Static_V_Fixed_Switch.m
@@ -0,0 +1,681 @@
+%variables inizialization utilizzo flessibile delle code , non ha vincoli
+%di tempo di servizio.% CHE THROUPUT GARANTISCONO? E COSTO AL VARIARE v
+clc
+clear all
+close all
+cc=0;
+N = 1;                                                                    %applications pag 479
+M = 100;                                                                   %Servers pag 479
+Ts = 2;                                                                    %time slot (sec)
+TT = 2000;
+T=2000;
+k_e=0.005; %(milliJ/(MHz)^2)
+%Lamda = 2000.*ones(N,TT);                                                  %Number of slot
+%A = 2.*Lamda.*rand(N,TT);                                                  %Arrival job uniformed distribuited pag 481
+%A=4000.*rand(N,TT);                                                        %PMR = 1
+%A = 3000 +(4000-2000).*rand(N,TT);                                         %PMR = 1.25
+%A = 2000 +(6000-2000).*rand(N,TT) ;                                        %PMR = 1.5
+% A = 1000 +(7000-1000).*rand(N,TT);                                        %PMR = 1.75
+%A = 8000.*rand(N,TT);                                                        %PMR = 2
+%save('Avalue4000PMR2.mat','A');
+%A=load('Avalue4000PMR2.mat','A');
+
+%L_tot_vect = 6+ round((10-6) * rand(N,TT));                         % PMR=1.25, L_tot=8+-2
+%save('workloads8PMR125No2000.mat', 'L_tot_vect');                  %loading 100 fixed random workloads
+A=load('workloads8PMR125No2000.mat','L_tot_vect');
+A=A.L_tot_vect;
+
+
+%F1 = [1.6*10^9 :0.2*10^9:2.6*10^9];                                         %frequecny range for every VM allocated onto the j-server
+F1 = [5,50,70,90,105];                                         %( Gbyte/s)frequecny range for every VM allocated onto the j-server
+I = ones(M,TT);                                                            %resource decision allocation
+P_cpu1 = zeros(N,M,length(F1));
+mu1=zeros(N,M,length(F1)) ;                                                  %power consuption for each VM onto each j server
+P_min1 = 0.*ones(N,M,length(F1));                                          %minimum power for the considered VM allocated onto the j-server
+P_max1 = 60.*ones(N,M,length(F1));                                          %maximum power for the considered VM allocated onto the j-server
+
+MU=zeros(N,M,TT);                                                          %service rate (request/slots)
+ww = zeros(1,N);
+AU = zeros(N,M);
+AAU = zeros(1,N);
+SAU = zeros(1,M);
+AAS = zeros(1,N);
+Avg_ADELAY = zeros(1,N);
+avg_MU=zeros(N,M);
+avg_AMU=zeros(1,N);
+avg_SMU=zeros(1,M);
+Avg_SDELAY = zeros(1,M);
+Avg_SDelay=zeros(1,M);
+
+KU=  zeros(N,M);
+Delay=zeros(1,N);
+Delay2=zeros(N,M);
+%Avg_Delay1=0;
+%UUTILITY = 0;
+
+alpha = ones(1,N);                                                         %throughput utility weights
+%V1 = [1:100:1000];                                                         %control parameter for DCA
+V = 100;
+%V=length(V1)
+beta = 1;                                                                  % non-negative normalizing weight
+W = zeros(N,TT);                                                           %Buffer Dimension
+K = zeros(N,M,TT);                                                         % ==R(i,j) number of requests for application i that are routed from the R(i) router Buffer to the j-server in slot t.
+U = zeros(N,M,TT);                                                          %queing dynamics for the request of application i at server j
+P = zeros(M,TT);                                                            % Switch consumed by each server in each time slot for s
+Switch=zeros(M,TT);                                                            % power consumed by each server in each time slot for s
+r = zeros(1,N);                                                            %average expected rate of admitted request for apllication i
+e = zeros(1,M);                                                            %the time average expected power consuption of server j
+SW= zeros(1,M);                                                            %the time average expected power consuption of server j
+R = zeros(N,TT);                                                           %the time average expected switch power consuption of server j
+Freq1 = zeros(N,M,length(F1));
+Freqtemp = zeros(M,TT);
+temp1 = zeros(N,M,length(F1));
+active_servers = zeros(N+1,M,T);
+M1=0;
+M2=0;
+M3=0;
+l=0;
+on_servers_list=zeros(1,M);
+on_servers = zeros(N+1,M);                                                 % on-server in each different frame T
+off_servers = zeros(N+1,M);                                                % off-server in each different frame T
+sleeping_servers= zeros(N+1,M);                                            % Hybernated-server in each different frame T
+ a = rand(1,M);                                                             % applications variable indicator
+ save('avalue.mat','a');
+a=load('avalue.mat','a');
+a=a.a;
+a(:)=(a(:) >= 0.4);
+ACTIVE_Num=1;
+a1=zeros(1, M);
+%a1(:) = a(:);
+count=0;                                                                   % static servers 1: static server is and 0: this servers is Dynamic
+for j=1:M
+    if (a(j)==1) && (count<ACTIVE_Num)
+        a1(j)=1;
+        count=count+1;
+    end
+end
+for j=1:M
+    if (a(j)==1)
+        on_servers(1,j)=j;
+        l=l+1;
+        on_servers_list(j)=j;
+    end
+end
+
+a_temp = zeros(1,M);                                                        %  temp applications variable indicator
+%a_temp(:)=(a_temp(:) >= 0.6);
+
+for i = 1:N
+    for j = 1:M
+        Freq1(i,j,:)= F1(:);
+    end
+end
+
+for i=1:N
+    for j=1:M
+        mu1(i,j,:) = (8*10^(-3).*F1(:)) +0.76;                  % servers queue for each application service rate (5,0.8)-----(105,1.6)
+        P_cpu1(i,j,:)=P_min1(i,j,:) + (50*10^(-4)).*((Freq1(i,j,:)-5).^2);     % CPU power of each virtual machine among each j-server pag 480
+    end
+end
+
+i=0;
+j=0;
+F3= [0, 5,50,70,90,105];    %frequecny range for every VM allocated onto the j-server
+temp3 = zeros(N,M,length(F3));
+P_cpu3 = zeros(N,M,length(F3));                                            % power consuption for each VM onto each j server
+P_min3 = 0.*ones(N,M,length(F3));                                          %minimum power for the considered VM allocated onto the j-server
+P_max3 = 60.*ones(N,M,length(F3));                                          %maximum power for the considered VM allocated onto the j-server
+mu3 = zeros(N,M,length(F3))  ;
+Freq3 = zeros(N,M,length(F3));
+for i = 1:N
+    for j = 1:M
+        Freq3(i,j,:)= F3(:);
+    end
+end
+
+for i=1:N
+    for j=1:M
+        mu3(i,j,2:1:length(F3)) = (8*10^(-3).*F3(2:1:length(F3)))+0.76; % servers queue for each application service rate (5,0.8)-----(105,1.6)
+        mu3(i,j,1) = 0;
+        P_cpu3(i,j,2:1:length(F3))=P_min3(i,j,2:1:length(F3)) + (50*10^(-4)).*((Freq3(i,j,2:1:length(F3))-(5)).^2);  % CPU power of each virtual machine among each j-server pag 480
+        P_cpu3(i,j,1)=0;
+        
+    end
+end
+F2= [0, 105];    %frequecny range for every VM allocated onto the j-server (STATIC server)
+P_cpu2 = zeros(N,M,length(F2));
+P_min2 = 0.*ones(N,M,length(F2));                                          %minimum power for the considered VM allocated onto the j-server
+P_max2 = 60.*ones(N,M,length(F2));                                          %maximum power for the considered VM allocated onto the j-server
+mu2 = zeros(N,M,length(F2))  ;
+%temp2 = zeros(N,M,length(F2));
+Freq2 = zeros(N,M,length(F2));
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%(STATIC server)
+for i=1:N
+    for j=1:M
+        mu2(i,j,2) = (8*10^(-3).*F2(2))+0.76;                                         % servers queue for each application service rate (5,0.8)-----(105,1.6)
+        mu2(i,j,1) = 0;
+        P_cpu2(i,j,2)=P_min2(i,j,length(F2)) + (50*10^(-4)).*((Freq2(i,j,length(F2))-(5)).^2);  % CPU power of each virtual machine among each j-server pag 480
+        P_cpu2(i,j,1)=0;
+        
+    end
+end
+
+i=0;
+j=0;
+UU=zeros(N,M,TT);
+UU1=zeros(N,M,TT);
+B = [1:1:TT/T];
+L_off=zeros(1,length(B));
+MM1=zeros(1,length(B));
+MM2=zeros(1,length(B));
+MM3=zeros(1,length(B));
+bb=0;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%main body code
+for t = 1 : TT
+    if (mod(t,T)~= 0)
+        UU=zeros(N,M,TT);
+        for jtest=1:M
+            if (a(jtest)==1)
+                UU(:,jtest,t)=U(:,jtest,t);
+            else
+                UU(:,jtest,t)=nan;
+            end
+        end
+        %if (t>1)
+        %UU(~UU)=nan;
+        %end
+        for i= 1 : N
+            
+            %
+            %j
+            t
+            if(W(i,t) > V.*alpha(i))
+                R(i,t) = 0;
+                
+            else
+                R(i,t) = A(i,t);
+            end
+            
+            %list U of active (ON) servers.
+            
+            
+            [U1, indexu1]=min(UU(i,:,t));% page 483 routing
+            if (W(i,t) > U1) || (isnan(U1))
+                K(i,indexu1,t) = W(i,t);
+            else
+                K(i,:,t)=0;
+            end
+            
+            UU(:,indexu1,t)=nan;
+            
+            %                 [U1, ind]=min(U(i,:,t));% page 483 routing
+            %                 if (W(i,t) > U1)
+            %                     for ikl=1:M
+            %                         for mmmm=1:N
+            %                             K(mmmm,ikl,t)=K(mmmm,ikl,t-1);
+            %                         end
+            %                     end
+            %                     K(i,ind,t) = W(i,t);
+            %
+            %                 else
+            %                     K(i,:,t)=0;
+            %                 end
+            
+            
+            %a(i,j) = (a(i,j) > 0.5); % garantee to set more than half to be -on
+            % fix number of active server, they are not M , they are
+            % obviously less
+            
+            %if sum(a(i,j)*K(i,j,t)<=W(i,t))
+            %if ((W(i,t)-sum(K(i,:,t)))>=0)
+            
+            % W(i,t+1) = W(i,t) - sum(a(i,j).*K(i,:,t)) + R(i,t);
+            W(i,t+1) = W(i,t) - (a*K(i,:,t)') + R(i,t);                %% Backlog Queue 1    (2) pag 481 consider only active server
+            % else
+            %    W(i,t+1)=R(i,t);
+            % end
+            
+        end % for N-application
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% after routing inside each server
+        for j= 1 : M
+            if (a(j)==1) && (a1(j)==0)                                      % for Dynamic ON server
+                [U_max1 , index_U_max1]=max(U(:,j,t));                   % 0=<index_U_max<=N
+                
+                for zzz=1:length(F1)
+                    temp1(index_U_max1,j,zzz) = U_max1.*(min(mu1(index_U_max1,j,zzz),U_max1)) - V.*beta.*(P_cpu1(index_U_max1,j,zzz));                  % Resource allocation pag 483
+                    
+                end
+                
+                [value1 , index1]=max(temp1(index_U_max1,j,:));            % 1=<index<=length(F)
+                
+                MU(index_U_max1,j,t) = min(mu1(index_U_max1,j,index1),U_max1);
+                P(j,t)= P_cpu1(index_U_max1,j,index1);
+                if (t==1)
+                    Switch(j,t)=k_e*(F1(index1))^2;
+                    Freqtemp(j,2)=F1(index1);
+                else
+                    Switch(j,t)=k_e*((F1(index1)-Freqtemp(j,t-1)))^2;
+                    Freqtemp(j,t)=F1(index1);
+                end
+                for i= 1 : N
+                    if (i==index_U_max1)
+                        U(index_U_max1,j,t+1) = max((U(index_U_max1,j,t) - MU(index_U_max1,j,t)),0) + K(index_U_max1,j,t);      %% Backlog Queue 2    (5) pag 481
+                    else
+                        U(i,j,t+1) = max(U(i,j,t),0) + K(i,j,t);           % Backlog Queue 2    (5) pag 481
+                    end
+                end
+                % U(index_U_max,j,t+1) = max(U(index_U_max,j,t) - MU(index_U_max,j,t),0) + a(i,j).*K(index_U_max,j,t);      %% Backlog Queue 2    (5) pag 481
+                
+   
+                
+       % for STATIC server % for STATIC server % for STATIC server % for STATIC server % for STATIC server % for STATIC server
+              
+                
+            elseif  (a(j)==1) && (a1(j)==1)                                  % for Static ON server
+                
+                [U_max2 , index_U_max2]=max(U(:,j,t));                  % 0=<index_U_max<=N
+                
+                MU(index_U_max2,j,t) = min(mu2(index_U_max2,j,length(F2)),U_max2);
+                %P(j,t)= P_cpu2(index_U_max2,j,length(F2));
+                P(j,t)= P_max2(index_U_max2,j,length(F2));
+                if (t==1)
+                    Switch(j,t)=k_e*(F2(length(F2)))^2;
+                    Freqtemp(j,2)=F2(length(F2));
+                else
+                    Switch(j,t)=k_e*((F2(length(F2))-Freqtemp(j,t-1)))^2;
+                    Freqtemp(j,t)=F2(length(F2));
+                end
+                for i= 1 : N
+                    if (i==index_U_max2)
+                        U(index_U_max2,j,t+1) = max((U(index_U_max2,j,t) - MU(index_U_max2,j,t)),0) + K(index_U_max2,j,t);      %% Backlog Queue 2    (5) pag 481
+                    else
+                        U(i,j,t+1) = max(U(i,j,t),0) + K(i,j,t);           % Backlog Queue 2    (5) pag 481
+                    end
+                end
+                
+            end
+        end                                                                % for each M-server
+        
+        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%check points
+    else
+        bb=bb+1;
+        UU1=zeros(N,M,TT);
+        a_temp = zeros(1,M); 
+        for jtest1=1:M
+            if (a(jtest1)==1)
+                UU1(:,jtest1,t)=U(:,jtest1,t);
+            else
+                UU1(:,jtest1,t)=nan;
+            end
+        end
+        %if (t>1)
+        %UU(~UU)=nan;
+        %end
+        for i= 1 : N
+            
+            if(W(i,t) > V.*alpha(i))   R(i,t) = 0;
+                
+            else
+                R(i,t) = A(i,t);
+            end
+            
+            %list U of active (ON) servers.
+            
+            [U2, indexu2]=min(UU1(i,:,t));% page 483 routing
+            
+            %[U2, indexu2]=min(U(i,:,t));                                   % page 483 routing
+            if (W(i,t) > U2)
+                K(i,indexu2,t) = a(indexu2)*W(i,t);
+            else
+                K(i,:,t)=0;
+            end
+            UU1(:,indexu2,t)=nan;
+            %a(i,j) = (a(i,j) > 0.5);                                      % garantee to set more than half to be -on
+            % fix number of active server, they are not M , they are
+            % obviously less
+            % if ((W(i,t)-sum(K(i,:,t)))>=0)
+            
+            % W(i,t+1) = W(i,t) - sum(a(i,j).*K(i,:,t)) + R(i,t);
+            W(i,t+1) = W(i,t) - (a*K(i,:,t)') + R(i,t);                %% Backlog Queue 1    (2) pag 481 consider only active server
+            % else
+            %W(i,t+1)=R(i,t);
+        end
+        
+        
+        for j= 1 : M
+            a_temp(j)=a(j);
+            
+             if (a1(j)==0)                                                  % for Dynamic server
+                 
+                 [U_max3, index_U_max3]=max(U(:,j,t));                           % 0=<index_U_max<=N
+                 
+                 for zzz=1:length(F3)
+                     temp3(index_U_max3,j,zzz) = U_max3.*(min(mu3(index_U_max3,j,zzz),U_max3)) - V.*beta.*(P_cpu3(index_U_max3,j,zzz));                  % Resource allocation pag 483
+                     
+                 end
+                 [value3, index3]=max(temp3(index_U_max3,j,:));                % 1=<index<=length(F)
+                 
+                 if (((min(mu3(index_U_max3,j,index3),U_max3))==0) && (P_cpu3(index_U_max3,j,index3)==0))
+                     cc=cc+1;
+                     a(j)=0;
+                     MU(index_U_max3,j,t) = 0;
+                     P(j,t)= 0;
+                     if (t==1)
+                         Switch(j,t)=k_e*(F3(index3))^2;
+                         Freqtemp(j,2)=F3(index3);
+                     else
+                         Switch(j,t)=k_e*(F3(index3)-Freqtemp(j,t-1))^2;
+                         Freqtemp(j,t)=F3(index3);
+                     end
+                 else
+                     a(j)=1;
+                     MU(index_U_max3,j,t) = min(mu3(index_U_max3,j,index3),U_max3);
+                     P(j,t)= P_cpu3(index_U_max3,j,index3);
+                     if (t==1)
+                         Switch(j,t)=k_e*(F3(index3))^2;
+                         Freqtemp(j,2)=F3(index3);
+                     else
+                         Switch(j,t)=k_e*(F3(index3)-Freqtemp(j,t-1))^2;
+                         Freqtemp(j,t)=F3(index3);
+                     end
+                 end
+            
+                 for i= 1 : N
+                     if (i==index_U_max3)
+                         U(index_U_max3,j,t+1) = max(U(index_U_max3,j,t) - MU(index_U_max3,j,t),0) + K(index_U_max3,j,t);    % Backlog Queue 2    (5) pag 481
+                     else
+                         U(i,j,t+1) = max(U(i,j,t),0) + K(i,j,t);           % Backlog Queue 2    (5) pag 481
+                     end
+                 end
+                 
+            % for STATIC server % for STATIC server % for STATIC server % for STATIC server % for STATIC server % for STATIC server
+            
+             elseif (a1(j)==1)                                             % for STATIC server
+                 
+                [U_max4 , index_U_max4]=max(U(:,j,t));                     % 0=<index_U_max<=N
+                MU(index_U_max4,j,t) = min(mu2(index_U_max4,j, length(F2)),U_max4);
+                %P(j,t)= P_cpu2(index_U_max4,j,length(F2));
+                P(j,t)= P_max2(index_U_max4,j,length(F2));
+                a(j)=1;                                                     % STATIC SERVER IS ALWAYS ON
+                if (t==1)
+                    Switch(j,t)=k_e*(F2(length(F2)))^2;
+                    Freqtemp(j,2)=F2(length(F2));
+                else
+                    Switch(j,t)=k_e*(F2(length(F2))-Freqtemp(j,t-1))^2;
+                    Freqtemp(j,t)=F2(length(F2));
+                end
+                for i= 1 : N
+                    if (i==index_U_max4)
+                        U(index_U_max4,j,t+1) = max(U(index_U_max4,j,t) - MU(index_U_max4,j,t),0) + K(index_U_max4,j,t);    % Backlog Queue 2    (5) pag 481
+                    else
+                        U(i,j,t+1) = max(U(i,j,t),0) + K(i,j,t);           % Backlog Queue 2    (5) pag 481
+                    end
+                end
+               
+             end  % end of for Dynamic server
+             
+            % U(index_U_max,j,t+1) = max(U(index_U_max,j,t) - MU(index_U_max,j,t),0) + a(i,j).*K(index_U_max,j,t);      %% Backlog Queue 2    (5) pag 481
+            % if
+        end                                                                % for each M-server
+        
+        
+        %*************************VM_MIGRATION******************
+        on_servers = zeros(N+1,M);                                                 % on-server in each different frame T
+        off_servers = zeros(N+1,M);                                                % off-server in each different frame T
+        sleeping_servers= zeros(N+1,M);
+        
+        for jj = 1 : M
+            for ii = 1 : N
+                if ((a_temp(jj)==0) && (a(jj)==1)) || (((a_temp(jj)==1) && (a(jj)==1)))
+                    on_servers(1,jj) = jj; %always ON servers + OFF--->> ON servers
+                    on_servers(ii+1,jj) =  U(ii,jj,t);
+                    %M1 = M1 + 1;
+                elseif((a_temp(jj)==1) && (a(jj)==0))
+                    off_servers(1,jj) = jj;
+                    off_servers(ii+1,jj) = U(ii,jj,t);
+                    %M2 = M2 + 1;
+                    
+                elseif ((a_temp(jj)==0) && (a(jj)==0))
+                    sleeping_servers(1,jj)=jj;  %the servers which are never ON, ALWAYS OFF
+                    sleeping_servers(ii+1,jj) =  U(ii,jj,t);
+                    %M3=M3+1;
+                end
+                
+            end
+            
+        end
+                
+        
+        
+%         on_servers = deletezero(on_servers);
+%        [x1, x2]=size(on_servers);
+%         i=0; j=0;
+%         for j=1:x2
+%             %for j=1:N
+%             if (on_servers(2:1:(N+1),j)==0)
+%               
+%                 if ((a_temp(j)==1) && (a(j)==0) && (a1(j)==0))
+%                     off_servers(1,j) = j;
+%                     off_servers(2:1:(N+1),j) = 0;
+%                 elseif ((a_temp(jj)==0) && (a(jj)==0) && (a1(j)==0))
+%                     sleeping_servers(1,j)=j;  %the servers which are never ON, ALWAYS OFF
+%                     sleeping_servers(2:1:(N+1),j) =  0;
+%                 end
+%             
+%             end
+%         end
+        on_servers = deletezero(on_servers);
+        on_servers = deletezero(on_servers);
+        off_servers= deletezero(off_servers);
+        sleeping_servers=deletezero(sleeping_servers);
+        [x1, x2]=size(on_servers);
+        M1=x2;
+        [y1, y2]=size(off_servers);
+        M2=y2;
+        [z1, z2]=size(sleeping_servers);
+        M3=z2;
+        
+        L_off(bb)=sum(sum(off_servers(2:N+1,:))) ;                                  %  the total migrated backlogs of all aplication in bits from the off servers
+        for jk = 1 : M2
+            for iii = 1 : N
+                
+                
+                [a_tempON, indexON] = min(on_servers(iii+1,:));
+                tempOnserver=off_servers(iii+1,jk);
+                on_servers(iii+1,indexON) = on_servers(iii+1,indexON)+tempOnserver;
+                off_servers(iii+1,jk)=0;
+                
+                
+                
+            end
+        end
+        
+        
+        a=zeros(1,M);
+        for dom=1:M1
+            if (on_servers(1,dom)>=1) && (a1(on_servers(1,dom))==0) % dynamic server
+                a(on_servers(1,dom))=1;
+                for i= 1 : N
+                    U(i,on_servers(1,dom),t+1)=on_servers(i+1,dom);
+                end
+            elseif (on_servers(1,dom)>=1) && (a1(on_servers(1,dom))==1) % static server
+                a(on_servers(1,dom))=1;
+                for i= 1 : N
+                    U(i,on_servers(1,dom),t+1)=on_servers(i+1,dom);
+                end
+            end
+        end
+        for dom=1:M2
+            a(off_servers(1,dom))=0;
+            for i= 1 : N
+                U(i,off_servers(1,dom),t+1)=0;
+            end
+        end
+                for dom=1:M3
+                   a(sleeping_servers(1,dom))=0;
+               end
+        
+        
+        MM1(bb)=M1;
+        MM2(bb)=M2;
+        MM3(bb)=M3;
+        M1=0;
+        M2=0;
+        M3=0;
+    end % t=nT
+    
+    
+end % biggggg for TT
+
+
+for ii= 1:N
+    r(ii) = mean(R(ii,:));
+    ww(ii) = mean(W(ii,:));                                           % average occupation for each applications before router
+    for jj = 1:M
+    e(jj) =   mean(P(jj,:));                                      %pag482 (7)
+    SW(jj)= mean (Switch(jj,:));
+    AU(ii,jj) = mean(U(ii,jj,:));                                     % average occupation of service rate for application i in sererver j in block 2
+    KU(ii,jj) = sum(K(ii,jj,:));                                     % average occupation of backlogs (in queue) of block 2 (queues in servers)
+    end
+    
+end
+
+KUU=0;
+
+    for i= 1:N
+        for j= 1:M
+            KUU=KU(i,j)+KUU;
+        end
+    end
+
+KUU=KUU/(M*N*TT);
+for i= 1:N
+Delay(i)=ww(i)/ r(i);  % delay takes for each application queue before routing
+
+end
+
+
+Avg_Delay1=mean(Delay(:)); % average delay of BLOCK 1 (before routing to the servers)
+
+
+
+for i= 1:N
+for j= 1:M
+%Delay2(i,j)=AU(i,j)./ KU(i,j);  % delay takes for application i in server j in block 2
+Delay2(i,j)=AU(i,j)./ KUU;  % delay takes for application i in server j in block 2
+end
+end
+
+
+
+for j= 1:M
+Avg_SDelay(j)=mean(Delay2(:,j)); % average delay takes for server j in block 2
+end
+
+
+Avg_Delay2=mean(Avg_SDelay(:));  % average dalay of BLOCK 2 (after routing to the servers and after servers )
+
+
+
+TOT_DELAY=Avg_Delay1+Avg_Delay2;  %average total delay of WHOLE SYSTEM
+UUTILITY = alpha(1).*sum(r(:)) - beta.*sum(e(:)); % UTILITY OF THE SYSTEM
+
+
+
+%AVG_UTILITY = mean(UUTILITY(:));   % total utility of the SYSTEM
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
+%for t= 1:TT-1
+% 
+% for ii= 1:N
+%     r(ii) = mean(R(ii,:));                                             % throughput 4 each ii
+%     ww(ii) = mean(W(ii,:));                                            % average occupation for each applications before router
+%     
+%     for jj = 1:M
+%         e(jj) = mean(P(jj,:));                                      %pag482 (7)
+%         AU(ii,jj) = mean(U(ii,jj,:));                                   % average occupation of backlogs applications of each server
+%         avg_MU(ii,jj)= mean(MU(ii,jj,:));
+%     end
+%     
+% end
+% %end
+% www=mean(ww(:));                                                           % average occpuation of all applications
+% for i = 1 : N
+%     AAU(i)=mean(AU(i,:));                                                       % average occupation of each backlogs application in the system
+%     avg_AMU(i)=mean(avg_MU(i,:));
+% end
+% 
+% for j = 1 : M
+%     SAU(j)=mean(AU(:,j));                                                       % average occupation of total backlogs applications in each  server j
+%     avg_SMU(jj)=mean(avg_MU(:,j));
+% end
+% 
+% TAU=mean(AAU(:));                                                          % average occupation for : server for : applications
+% Tavg_MU=mean(avg_AMU(:));
+% 
+% 
+% AS = www + TAU ;                                                           % total average occupations of the system
+% 
+% Avg_DELAY =  AS/mean(r(:));                                                % average delay of the system
+% 
+% AAS(:) = ww(:) + AAU(:) ;                                                  % total average occupations of the system for each application
+% 
+% Avg_ADELAY(:) = AAS(:)./r(:);                                              % average delay of the system for each application
+% 
+% Avg_SDELAY(:) = SAU(:)./avg_SMU(:);
+% 
+% UTILITY = alpha(1)*sum(r(:)) - beta*sum(e(:));
+% 
+% 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/ICC2015_Sourcecode/Mu_opt_bisezione.m b/ICC2015_Sourcecode/Mu_opt_bisezione.m
new file mode 100644
index 0000000000000000000000000000000000000000..199aacc01a135a3681244459b2d6b8be21ec1749
--- /dev/null
+++ b/ICC2015_Sourcecode/Mu_opt_bisezione.m
@@ -0,0 +1,65 @@
+%calcolo moltiplicatori con metodo bisezione
+
+function [mu,delta_mu]= Mu_opt_bisezione(alpha_zero,alpha_mu,P_net,C_max,f_zero,f_max,Delta,L_b,L_tot)
+
+
+%NB.
+%LA FUNZIONE DI CUI CALCOLARE LO ZERO è ASSUNTA MONOTONA CRESCENTE
+%CON VALORE IN ZERO NEGATIVO E ALL'INFINITO POSITIVO
+%SI ASSUME CHE VENGA CHIAMATA SOLO NEI CASI IN CUI LO ZERO ESISTE FINITO.
+%ALTRIMENTI VA IN LOOP!
+
+% definisco gli estremi in cui cercare 
+a_mu = 0;% per mu=0 si ha sum(L)<L_tot 
+%b_mu = 2.*10^3;  %<====  NB. SCEGLIERE CON ATTENZIONE IN BASE AI PARAMETRI DI SISTEMA!
+b_mu=max(((1./alpha_mu).*(f_max-alpha_zero.*f_zero))+2.*P_net./C_max); % <== per questo mu si ha f=f_max quindi sum(L)>L_tot
+%b_mu = 2.*10^0;
+
+
+toll = 10^-6;
+
+fa=delta_carico(a_mu,alpha_zero,alpha_mu,P_net,C_max,f_zero,f_max,Delta,L_b,L_tot);
+fb=delta_carico(b_mu,alpha_zero,alpha_mu,P_net,C_max,f_zero,f_max,Delta,L_b,L_tot);
+
+% verifico che la funzione abbia uno zero  
+while (fa.*fb) > 0 
+    b_mu=2.*b_mu;
+    fb=delta_carico(b_mu,alpha_zero,alpha_mu,P_net,C_max,f_zero,f_max,Delta,L_b,L_tot);
+    %fprintf('%d\n',t)
+    error('estremi iniziali ricerca mu_opt errato')
+  
+end
+
+
+iterazioni = ceil(log2(b_mu-a_mu)-log2(toll));  
+%     fprintf('%d\n......iteraion',iterazioni);
+%     fprintf('%d\n......b_mu',b_mu);
+%     fprintf('%d\n......a_mu',a_mu);
+%     fprintf('%d\n......toll',toll);
+
+   for i = 1: iterazioni
+       c_mu=(a_mu+b_mu)/2;% punto medio
+       fc = delta_carico(c_mu,alpha_zero,alpha_mu,P_net,C_max,f_zero,f_max,Delta,L_b,L_tot);
+       if abs(fc)<toll
+           break
+       end
+        if abs(b_mu-a_mu)< toll
+            break
+        end      
+%         tolf = toll.*(abs((fb-fa)./(b_mu-a_mu)));     %ELIMINATO - NEL CASO
+%                                                       DI GRADINO FUORVIANTE
+%         if abs(fc)<=tolf
+%             'Attenzione Probabile Gradino'
+%             i
+%             break
+%         end
+        if (fa.*fc)<0 %la soluzione è fra a e il punto medio
+            b_mu=c_mu;     %tengo la met sinistra dell'intervallo
+            fb = fc;
+        else         % altrimenti tengo la metà destra
+            a_mu=c_mu;
+            fa = fc;
+        end
+   end
+    mu=c_mu;
+    delta_mu=fc;
diff --git a/ICC2015_Sourcecode/X.mat b/ICC2015_Sourcecode/X.mat
new file mode 100644
index 0000000000000000000000000000000000000000..7ef5d3d6c698db0217767edf2f306a9221cfac39
Binary files /dev/null and b/ICC2015_Sourcecode/X.mat differ
diff --git a/ICC2015_Sourcecode/copyright notice.docx b/ICC2015_Sourcecode/copyright notice.docx
new file mode 100644
index 0000000000000000000000000000000000000000..784f337f42d8b9b1b48e5d0f96aff2ff41d5cf93
Binary files /dev/null and b/ICC2015_Sourcecode/copyright notice.docx differ
diff --git a/ICC2015_Sourcecode/deletezero.m b/ICC2015_Sourcecode/deletezero.m
new file mode 100644
index 0000000000000000000000000000000000000000..fe70e1ad16ecdf809175e19ede1ec5bd4b6c2aec
--- /dev/null
+++ b/ICC2015_Sourcecode/deletezero.m
@@ -0,0 +1,12 @@
+function C = deletezero(A)
+
+%delete columns of A containing all zero entries
+
+b = sum(abs(A));
+if (sum(b == 0) >= 1)
+    ind = find(b == 0,1,'first');
+    A(:,ind) = [];
+    C = deletezero(A);
+else
+    C = A;
+end
\ No newline at end of file
diff --git a/ICC2015_Sourcecode/delta_carico.m b/ICC2015_Sourcecode/delta_carico.m
new file mode 100644
index 0000000000000000000000000000000000000000..632f943667173b5361284a2fe0de7ec42ce3b4ec
--- /dev/null
+++ b/ICC2015_Sourcecode/delta_carico.m
@@ -0,0 +1,16 @@
+%FUNZIONE CHE DEFINISCE L'EQUAZIONE DA RISOLVERE PER IL CALCOLO DEL Mu
+%ottimo (equazione di carico su L_tot)
+
+
+
+
+function y=delta_carico(mu,alpha_zero,alpha_mu,P_net,C_max,f_zero,f_max,Delta,L_b,L_tot)
+
+f_mu=max(mu-2.*P_net./C_max,0); 
+f_star=alpha_zero.*f_zero+alpha_mu.*f_mu;
+f_opt=max(0,min(f_star,f_max));
+
+canali_attivi=(f_mu>0); 
+L_opt=canali_attivi.*(f_opt.*Delta-L_b);
+y=sum(L_opt)-L_tot;  % N.B. PER MU=0 Y<O; PER MU-->INF Y>0  
+ 
\ No newline at end of file
diff --git a/ICC2015_Sourcecode/grafica_delta_carico.m b/ICC2015_Sourcecode/grafica_delta_carico.m
new file mode 100644
index 0000000000000000000000000000000000000000..853884ef0d798ed1e9d1a36a6ebbd1830cbd232c
--- /dev/null
+++ b/ICC2015_Sourcecode/grafica_delta_carico.m
@@ -0,0 +1,6 @@
+function y=grafica_delta_carico(mu,alpha_zero,alpha_mu,P_net,C_max,f_zero,f_max,Delta,L_b,L_tot)
+
+y=zeros(1,length(mu));
+for k=1:length(mu)
+    y(k)=delta_carico(mu(k),alpha_zero,alpha_mu,P_net,C_max,f_zero,f_max,Delta,L_b,L_tot);
+end
\ No newline at end of file
diff --git a/ICC2015_Sourcecode/tracking3.m b/ICC2015_Sourcecode/tracking3.m
new file mode 100644
index 0000000000000000000000000000000000000000..b43fa87d83443aabb3097be330c7cb38667961d7
--- /dev/null
+++ b/ICC2015_Sourcecode/tracking3.m
@@ -0,0 +1,180 @@
+
+%VERSIONE DEL TRACKING CON PUNTO INIZIALE DI RICERCA PER MU.
+
+function [conv,iter,mu_opt,ni_opt,f_opt,L_opt,alpha_step_finale,V_step_finale,mu_vect]= tracking3(L_tot,L_b,alpha_zero,alpha_mu,chi,W,T_tot,Delta,f_zero,f_max,M,Th,mu_iniz,ni_iniz,f_iniz,L_iniz,alpha_step_iniz,V_step_iniz)
+
+
+
+epsilon_ni=0;  %<== N.B. shift numerico introdotto su ni;
+
+toll = 10^(-3);%10^(-2);
+num_iterazioni=100000;
+conv=0; %viene settato a 1 se il tracking arriva a convergenza.
+
+ni=zeros(num_iterazioni,M);
+ni(1,:)=ni_iniz;
+f=zeros(num_iterazioni,M);
+f(1,:)=f_iniz;
+L=zeros(num_iterazioni,M);
+L(1,:)=L_iniz;
+mu=zeros(num_iterazioni,1);
+mu(1)=mu_iniz; %<==  importante, inizializzazione del parametro mu.
+
+
+beta=0.04;%0.01;%0.7;%0.1;%10^(-3);%10^(-4);
+gamma=0.6;%0.4;
+alpha_step=beta.*ones(num_iterazioni,1);
+%alpha_step(1)=alpha_step_iniz;
+V_step=zeros(num_iterazioni,1);
+%V_step(1)=V_step_iniz;
+
+gap=zeros(num_iterazioni,1);
+
+
+for k=2:num_iterazioni
+    
+    gap(k-1)=sum(L(k-1,:))-L_tot;
+    if abs(gap(k-1)./L_tot)<toll   %errore relativo sotto la tolleranza
+        conv=1;
+        break
+    end
+    mu(k)=mu(k-1)-alpha_step(k-1).*gap(k-1);
+    mu(k)=max(0,mu(k));
+    
+    ni_star=mu(k)-2.*log(2).*(chi./W).*2.^(2.*L(k-1,:)./(W.*(T_tot-Delta)));
+    ni(k,:)=max(0,ni_star);
+    
+    f_star=alpha_zero.*f_zero+alpha_mu.*ni(k,:);
+    f(k,:)=max(0,min(f_star,f_max)); %N.B. --> NEL CASO GENERALE DI L_b>0 IL MAX VA FATTO RISPETTO A L_b/Delta, NON RISPETTO A 0.
+    
+    canali_attivi=ni(k,:)>epsilon_ni;  %<== N.B. shift numerico introdotto su ni;;
+    %L(k,:)=canali_attivi.*(f(k,:).*Delta-L_b)+(~canali_attivi).*max(0,(((T_tot-Delta).*W./2).*log(mu(k)./Th)./log(2)));
+    %L(k,:)=max(0,(((T_tot-Delta).*W./2).*log(mu(k)./Th)./log(2)));
+    L(k,:)=min((f(k,:).*Delta-L_b),max(0,(((T_tot-Delta).*W./2).*log(mu(k)./Th)./log(2))));
+    
+    
+    
+    
+    V_step(k)=(1-alpha_step(k-1)).*V_step(k-1)-gap(k-1);
+    alpha_step(k)=max(0,(min(beta,alpha_step(k-1)-gamma.*V_step(k-1).*gap(k-1))));
+    %alpha_step(k)=(min(beta,alpha_step(k-1)-gamma.*V_step(k-1).*gap(k-1)));
+    
+end
+
+
+
+%  V_step
+%  alpha_step
+% gap
+% mu
+% ni
+% f
+% L
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   GRAFICI
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   %%%%%%%%
+% 
+% figure(num_iterazioni+1)
+% plot(gap,'--o'),xlabel('iterazioni'),ylabel('gap');
+% 
+% figure(num_iterazioni+2)
+% plot(mu,'--o'),xlabel('iterazioni'),ylabel('mu');
+% 
+% figure(num_iterazioni+3)
+% plot(ni,'--o'),xlabel('iterazioni'),ylabel('ni');
+% 
+% figure(num_iterazioni+4)
+% plot(f,'--o'),xlabel('iterazioni'),ylabel('f');
+% 
+% figure(num_iterazioni+5)
+% plot(L,'--o'),xlabel('iterazioni'),ylabel('L');
+% 
+% figure(num_iterazioni+6)
+% plot(V_step,'--o'),xlabel('iterazioni'),ylabel('V_step');
+% 
+% figure(num_iterazioni+7)
+% plot(alpha_step,'--o'),xlabel('iterazioni'),ylabel('alpha_step');
+
+
+
+
+%  I VALORI SONO CORRETTI SOLTANTO SE CONV=1
+iter=k-1; %numero iterazioni per convergenza
+mu_opt=mu(k-1);
+mu_vect=mu(1:k-1);% < == N.B. Restituisce in uscita il vettoredi tracking
+ni_opt=ni(k-1,:);
+f_opt=f(k-1,:);
+L_opt=L(k-1,:);
+
+alpha_step_finale=alpha_step(k-1);
+V_step_finale=V_step(k-1);
+    
+    
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+% 
+% 
+% 
+% 
+% %NB.
+% %LA FUNZIONE DI CUI CALCOLARE LO ZERO è ASSUNTA MONOTONA CRESCENTE
+% %CON VALORE IN ZERO NEGATIVO E ALL'INFINITO POSITIVO
+% %SI ASSUME CHE VENGA CHIAMATA SOLO NEI CASI IN CUI LO ZERO ESISTE FINITO.
+% %ALTRIMENTI VA IN LOOP!
+% 
+% % definisco gli estremi in cui cercare 
+% a_mu = 0;% per mu=0 si ha sum(L)<L_tot 
+% %b_mu = 2.*10^3;  %<====  NB. SCEGLIERE CON ATTENZIONE IN BASE AI PARAMETRI DI SISTEMA!
+% b_mu=max(((1./alpha_mu).*(f_max-alpha_zero.*f_zero))+2.*P_net./C_max); % <== per questo mu si ha f=f_max quindi sum(L)>L_tot
+% %b_mu = 2.*10^0;
+% 
+% 
+% toll = 10^-6;
+% 
+% fa=delta_carico(a_mu,alpha_zero,alpha_mu,P_net,C_max,f_zero,f_max,Delta,L_b,L_tot);
+% fb=delta_carico(b_mu,alpha_zero,alpha_mu,P_net,C_max,f_zero,f_max,Delta,L_b,L_tot);
+% 
+% % verifico che la funzione abbia uno zero  
+% while (fa.*fb) > 0 
+%     b_mu=2.*b_mu;
+%     fb=delta_carico(b_mu,alpha_zero,alpha_mu,P_net,C_max,f_zero,f_max,Delta,L_b,L_tot);
+%     %fprintf('%d\n',t)
+%     error('estremi iniziali ricerca mu_opt errato')
+%   
+% end
+% 
+% 
+% iterazioni = ceil(log2(b_mu-a_mu)-log2(toll));  
+%     
+% 
+%    for i = 1: iterazioni
+%        c_mu=(a_mu+b_mu)/2;% punto medio
+%        fc = delta_carico(c_mu,alpha_zero,alpha_mu,P_net,C_max,f_zero,f_max,Delta,L_b,L_tot);
+%        if abs(fc)<toll
+%            break
+%        end
+%         if abs(b_mu-a_mu)< toll
+%             break
+%         end      
+% %         tolf = toll.*(abs((fb-fa)./(b_mu-a_mu)));     %ELIMINATO - NEL CASO
+% %                                                       DI GRADINO FUORVIANTE
+% %         if abs(fc)<=tolf
+% %             'Attenzione Probabile Gradino'
+% %             i
+% %             break
+% %         end
+%         if (fa.*fc)<0 %la soluzione è fra a e il punto medio
+%             b_mu=c_mu;     %tengo la met sinistra dell'intervallo
+%             fb = fc;
+%         else         % altrimenti tengo la metà destra
+%             a_mu=c_mu;
+%             fa = fc;
+%         end
+%    end
+%     mu=c_mu;
+%     delta_mu=fc;
diff --git a/ICC2015_Sourcecode/variabili.m b/ICC2015_Sourcecode/variabili.m
new file mode 100644
index 0000000000000000000000000000000000000000..1e3b8aac4d252ee02fe4fc750ff31c252fe3f80f
--- /dev/null
+++ b/ICC2015_Sourcecode/variabili.m
@@ -0,0 +1,10 @@
+clear all;
+close all;
+clc;
+X = rand(1,10000000);
+EPS = rand(1,10000000);
+DELTA_IP = rand(1,10000000);
+
+save('X');
+save('EPS');
+save('DELTA_IP');
\ No newline at end of file