clear all
close all
clc
% Modelo con solo ruteo +1 vehículos, +1 residuo, para +1 periodo
%Enfatizando en la restricción de carga acumulada para evitar subtours
% CON VIDA ÚTIL

% Condiciones parámetros - De prueba

% 1. Número de farmacias: Peq: 5:25 - Gra: 40:50 [Ajus P: 5:20 - Gra: 30:50]
I = 6; % 2 farmacias

% 2. Número de centros de tratamiento: 1 de cada tipo - total: 3
J = 8; % 3 centros de tratamiento

% 3. Número de depósitos: 1 (siempre)
d = 1; 

% Número de nodos (d+i+j)
n = 8; 

% 4. Número de vehículos: Máx 3 categorías dif con hasta 2 vehículos c/u
K = 2;
M = 50;
g = 2; %Nodos entre recolección-disposición

% 5. Horizonte temporal: Periodos diarios con ciclos de 5,7 y 20 días
T = 3;

% 6. Tipos de residuo: 2 peligroso y no peligroso
R = 2;

% 7. Capacidad del vehículo: debe tener coherencia con las unds a mover
Q_k = [120 80];

% 8. Costos: costo fijo de transporte es 5 veces mayor al costo del
% transporte y el costo de inventario se mantiene en el mismo rango que
% este último

% 8.1. Costo de transporte: tamaño nxn, penalización (valor más alto) del
% depósito a los centros de tratamiento, valor cero desde los Ct al D,
% valores entre 10 - 100
% Dimensión de la matriz - número de nodos - n
rng(50);
CT = randi([10,100],n,n,K,T); % Generar matriz
CT(1,n-1:end) = 1000;
CT(n-2:end,1) = 0;

% 8.2. Costo inventario: tamaño IxI
%I = 5; % Dimensión de la matriz - número de nodos
rng(45)
L_max = 2; % Número de niveles (ajustable si cambia la estructura)
F = I-1;
CI = randi([50,500],F,L_max,T); % Generar matriz

% 8.3. Costo vehículo (fijo)
% Dimensión de la matriz - número de nodos
rng(34)
CV = randi([20,250],K,T); % Generar matriz/ vector columna

% 9. Inventario a recolectar en farmacias: Niveles de inventario de desde 8
% a 24 unidades, la tasa de generación suele ser de 1 a 5. Encontrar una
% combinación entre estos rangos
% Z_ur_nor = [0 10 2 -1 -11; 0 3 2 -3 -2; 0 1 8 -1 -8; 0 7 6 -7 -6];
% Cargar la matriz desde el archivo .mat
load('IIP_3.mat');  % Ahora A está disponible
load('CIIP_3.mat');
% Z_urt = reshape(Z_urt_P,[n,R,T]);}

% whos('IP_1')
% whos('CIP_1')

% Variables
X_uvkt = optimvar('X_uvkt',[n, n, K, T],'Type','integer','LowerBound',0,'UpperBound',1); % Binaria para el ruteo
w_urkt = optimvar('w_urkt',[n,R,K,T], 'Type', 'integer', 'LowerBound', 0); % Carga acumulada
y_kt = optimvar('y_kt',[K,T], 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1); % Binaria activación vehículos

disp(size(X_uvkt))
disp(size(w_urkt))
disp(size(y_kt))

% Modelo
Z = optimproblem("ObjectiveSense","min"); % Modelo de minimización
Z.Objective = sum(X_uvkt.*CT,"all") + sum(y_kt.*CV,'all') + sum(INVFCilt.*CI,'all'); % FO considerando costo de transporte y costo fijo de los vehículos

% Restricciones


% R1 - Activa los arcos de i (farmacias) a u (todos) - Aplicando que para
% cada vehículo puede ser 1 (se activa) o cero (no se usa)
for t = 1:T
    for k = 1:K;
        for i = 2:I;
        suma_Xu = 0;
        for u = 1:n;
            suma_Xu = suma_Xu + X_uvkt(i,u,k,t);
        end
        Z.Constraints.(['R1_ikt',num2str(i),'_',num2str(k),'_',num2str(t)]) = suma_Xu <= 1;
        end
    end
end
disp('R1_ikt')

% R2 - Activa los arcos de u (todos) a i (farmacias)
for t = 1:T
    for k = 1:K;
        for i = 2:I;
        suma_Xu = 0;
        for u = 1:n;
            suma_Xu = suma_Xu + X_uvkt(u,i,k,t);
        end
        Z.Constraints.(['R2_ikt',num2str(i),'_',num2str(k),'_',num2str(t)]) = suma_Xu <= 1;
        end
    end
end
disp('R2_ikt')

% R3: Garantiza que para cada nodo la suma de arcos activos con todos los
% vehículos puede ser mínimo 1
for t = 1:T
    for u = 1:n;
        suma_Xvk = 0;
        for v = 1:n;
            for k = 1:K;
                if u ~= v;
                suma_Xvk = suma_Xvk + X_uvkt(v,u,k,t);
                end
            end
        end
    Z.Constraints.(['R3A_ut',num2str(u),num2str(t)]) = suma_Xvk >= 1;
    end
end
disp('R3A_ut')

% R4 - Activa los arcos de j (centros de tratamiento) a u (todos)
for t = 1:T
    for k = 1:K;
        for j = I+1:J;
        suma_Xu = 0;
        for u = 1:n;
            suma_Xu = suma_Xu + X_uvkt(j,u,k,t);
        end
        Z.Constraints.(['R4_jkt',num2str(j),'_',num2str(k),'_',num2str(t)]) = suma_Xu <= 1;
        end
    end
end
disp('R4_jkt')

% R5 - Activa los arcos de u (todos) a j (centros de tratamiento)
for t = 1:T
    for k = 1:K;
        for j = I+1:J;
        suma_Xu = 0;
        for u = 1:n;
            suma_Xu = suma_Xu + X_uvkt(u,j,k,t);
        end
        Z.Constraints.(['R5_jkt',num2str(j),'_',num2str(k),'_',num2str(t)]) = suma_Xu <= 1;
        end
    end
end
disp('R5_jkt')

% R6: Conservación/balance de flujo para cada vehículo
for t = 1:T
    for v = 1:n;
        for k = 1:K;
            suma_Xuv = 0;
            suma_Xvu = 0;
            for u = 1:n;
                if u ~= v;
                    suma_Xuv = suma_Xuv + X_uvkt(u,v,k,t);
                    suma_Xvu = suma_Xvu + X_uvkt(v,u,k,t);
                end
            end
            Z.Constraints.(['R6_vkt',num2str(v),num2str(k),num2str(t)]) = suma_Xuv - suma_Xvu == 0;
        end
    end
end
disp('R6_vkt')

% R7: Cantidad de vehículos que salen del depósito 
for t = 1:T
    for k = 1:K;
        suma_Xi = 0; 
        for i = 2:I;
            suma_Xi = suma_Xi + X_uvkt(1,i,k,t);
        end
        Z.Constraints.(['R7A_kt',num2str(k),num2str(t)]) = suma_Xi == y_kt(k,t);
    end
end 
disp('R7A_kt')

for t = 1:T
    for k = 1:K;
        suma_Xj = 0; 
        for j = I+1:J;
            suma_Xj = suma_Xj + X_uvkt(1,j,k,t);
        end
        Z.Constraints.(['R7B_kt',num2str(k),num2str(t)]) = suma_Xj == 0;
    end
end 
disp('R7B_kt')

% R8: Cantidad de vehículos que regresan al depósito (activados)
for t = 1:T
    for k = 1:K;
        suma_Xj = 0; 
        for j = I+1:J;
            suma_Xj = suma_Xj + X_uvkt(j,1,k,t);
        end
        Z.Constraints.(['R8A_kt',num2str(k),num2str(t)]) = suma_Xj == y_kt(k,t);
    end
end
disp('RA8_kt')

for t = 1:T
    for k = 1:K;
        suma_Xi = 0; 
        for i = 1:I;
            suma_Xi = suma_Xi + X_uvkt(i,1,k,t);
        end
        Z.Constraints.(['R8B_kt',num2str(k),num2str(t)]) = suma_Xi == 0;
    end
end
disp('R8B_kt')


% R9: Límite superior de la suma de vehículos activos
% No está influyendo de la forma que debería en el modelo
for t = 1:T
    suma_yk = 0;
    for k = 1:K;
        suma_yk = suma_yk + y_kt(k,t);
    end
    Z.Constraints.(['R9A_t',num2str(t)]) = suma_yk <= K;
    Z.Constraints.(['R9B_t',num2str(t)]) = suma_yk >= 1;
end
disp('R9A_t')
disp('R9B_t')


% R10: Uso/activación de vehículos
for t = 1:T
    for k = 1:K;
        suma_Xuv = 0;
        for u = 1:n;
            for v = 1:n;
                if u ~= v;
                suma_Xuv = suma_Xuv + X_uvkt(u,v,k,t);
                end
            end
        end
        % Z.Constraints.(['R10A_kt',num2str(k),num2str(t)]) = suma_Xuv <= Q_k(k)*y_kt(k,t); %min
        Z.Constraints.(['R10A_kt',num2str(k),num2str(t)]) = suma_Xuv <= 2*n*y_kt(k,t); %min
        Z.Constraints.(['R10B_kt',num2str(k),num2str(t)]) = suma_Xuv >= y_kt(k,t); %max
    end
end
disp('R10A_kt')
disp('R10B_kt')

% R11: Puente (solo 1 paso de i a j)
for t = 1:T
    for k = 1:K;
        suma_Xij = 0;
        for i = 2:I;
            for j = I+1:J;
            suma_Xij =  suma_Xij + X_uvkt(i,j,k,t);
            end
        end
        Z.Constraints.(['R11_kt',num2str(k),num2str(t)]) = suma_Xij == y_kt(k,t);
    end
end
disp('R11_kt')

% R12: Carga transportada no excede capacidad del vehículo - Ver final
for t = 1:T
    for k = 1:K;
        suma_wr = 0;
        for r = 1:R;
            suma_wr = suma_wr + w_urkt(1,r,k,t);
        end
        Z.Constraints.(['R12AA_ukt',num2str(k),num2str(t)]) = suma_wr == 0;
    end
end
disp('R12AA_ukt')

% Límite superior carga acumulada
for t = 1:T
    for u = 1:n;
        for k = 1:K;
            suma_wr = 0;
            for r = 1:R;
                suma_wr = suma_wr + w_urkt(u,r,k,t);
            end
            Z.Constraints.(['R12B_ukt',num2str(u),num2str(k),num2str(t)]) = suma_wr <= Q_k(k)*y_kt(k,t);
        end 
    end
end
disp('R12B_ukt')

% % 3. También es necesario asegurar que la carga acumulada respete la capacidad del vehículo
% % Mejorar la restricción R12B para considerar correctamente la carga por tipo de residuo
% for t = 1:T
%     for k = 1:K
%         % Suma total de todos los residuos que recolecta el vehículo k en el período t
%         suma_residuos_total = 0;
%         for i = 2:I
%             for r = 1:R
%                 suma_Xi = 0;
%                 for v = 1:n
%                     if i ~= v
%                         suma_Xi = suma_Xi + X_uvkt(i,v,k,t);
%                     end
%                 end
%                 suma_residuos_total = suma_residuos_total + suma_Xi * Z_urt(i,r,t);
%             end
%         end
% 
%         % La cantidad total recolectada no debe exceder la capacidad del vehículo
%         Z.Constraints.(['R12C_kt',num2str(k),num2str(t)]) = suma_residuos_total <= Q_k(k);
%     end
% end
% disp('R12C_kt');

% % R13: Emparejamiento nodos de recolección y disposición
% for t = 1:T
%     for k = 1:K;
%         for i = 2:I;
%             suma_Xj = 0;
%             suma_Xgj = 0;
%             for j = I+1:J;
%                 if i ~= j;
%                     if max(Q_k) < sum(Z_urt(:,:,t))
%                     if j ~= i+g;
%                         suma_Xj = suma_Xj + X_uvkt(i,j,k,t);
%                         suma_Xgj = suma_Xgj + X_uvkt(i+g,j,k,t);
%                     end
%                     else 
%                         suma_Xj = suma_Xj + X_uvkt(i,j,k,t);
%                         suma_Xgj = suma_Xgj + X_uvkt(i+g,j,k,t);
%                    end
%                 end
%             end 
%         Z.Constraints.(['R13_ikt',num2str(i),num2str(k),num2str(t)]) = suma_Xj == suma_Xgj;
%         end
%     end 
% end
% disp('R13_ikt')

% R14: Evita subtours - carga positiva
for t = 1:T
    for i = 2:I;
        for v = 2:n;
            for k = 1:K;
                for r = 1:R;
                    if i ~= v;
                        Z.Constraints.(['R14_ivrkt',num2str(i),num2str(v),num2str(r),num2str(k),num2str(t)]) = ...
                            w_urkt(v,r,k,t) >= w_urkt(i,r,k,t) + Z_urt(v,r,t)*X_uvkt(i,v,k,t) - Q_k(k)*(1 - X_uvkt(i,v,k,t));
                    end
                end
            end
        end
    end
end
disp('R14_ivrkt')

% R15: Evita subtours - carga negativa
for t = 1:T
    for u = 2:n;
        for j = I+1:J;
            for k = 1:K;
                for r = 1:R;
                    if j ~= u;
                        Z.Constraints.(['R15_ujrkt',num2str(u),num2str(j),num2str(r),num2str(k)]) = ...
                            w_urkt(u,r,k,t) >= w_urkt(j,r,k,t) - Z_urt(u,r,t)*X_uvkt(u,j,k,t) - Q_k(k)*(1 - X_uvkt(u,j,k,t));
                    end
                end
            end
        end
    end
end
disp('R15_ujrkt')

% Asegura que toda la demanda se satisfaga
for t = 1:T
    for u = 2:I
        for r = 1:R;
            if Z_urt(u,r,t) > 0
                suma_Zrec = 0;
                for k = 1:K
                    suma_Xv = 0;
                    for v = 1:n
                        suma_Xv = suma_Xv + X_uvkt(u,v,k,t);
                    end
                    suma_Zrec = suma_Zrec + suma_Xv*Z_urt(u,r,t);
                end
                Z.Constraints.(['R16_urt',num2str(u),num2str(r),num2str(t)]) = suma_Zrec == Z_urt(u,r,t);
            end
        end
    end
end
disp('R16_urt')

% R17: Restricción para asegurar que el residuo tipo 1 va al nodo n-1 y el residuo tipo 2 va al nodo n

% Para residuo tipo 1 dirigido al nodo n-1 (penúltimo nodo)
for t = 1:T
    for k = 1:K
        for i = 2:I
            % Si hay residuo tipo 1 en el nodo i que es recogido por el vehículo k
            if Z_urt(i,1,t) > 0
                % Debe existir un arco que conecte algún nodo con el nodo n-1 para el vehículo k
                suma_X_to_n1 = 0;
                for u = 1:n
                    if u ~= (n-1) % Cualquier nodo diferente a n-1
                        suma_X_to_n1 = suma_X_to_n1 + X_uvkt(u,n-1,k,t);
                    end
                end
                
                % Si el vehículo k recoge en el nodo i, debe visitar el nodo n-1
                suma_X_from_i = 0;
                for v = 1:n
                    if i ~= v
                        suma_X_from_i = suma_X_from_i + X_uvkt(i,v,k,t);
                    end
                end
                
                Z.Constraints.(['R17A_ikt',num2str(i),num2str(k),num2str(t)]) = ...
                    suma_X_to_n1 >= suma_X_from_i * (Z_urt(i,1,t) > 0);
            end
        end
    end
end
disp('R17A_ikt');

% Para residuo tipo 2 dirigido al nodo n (último nodo)
for t = 1:T
    for k = 1:K
        for i = 2:I
            % Si hay residuo tipo 2 en el nodo i que es recogido por el vehículo k
            if Z_urt(i,2,t) > 0
                % Debe existir un arco que conecte algún nodo con el nodo n para el vehículo k
                suma_X_to_n = 0;
                for u = 1:n
                    if u ~= n % Cualquier nodo diferente a n
                        suma_X_to_n = suma_X_to_n + X_uvkt(u,n,k,t);
                    end
                end
                
                % Si el vehículo k recoge en el nodo i, debe visitar el nodo n
                suma_X_from_i = 0;
                for v = 1:n
                    if i ~= v
                        suma_X_from_i = suma_X_from_i + X_uvkt(i,v,k,t);
                    end
                end
                
                Z.Constraints.(['R17B_ikt',num2str(i),num2str(k),num2str(t)]) = ...
                    suma_X_to_n >= suma_X_from_i * (Z_urt(i,2,t) > 0);
            end
        end
    end
end
disp('R17B_ikt');

% Restricción adicional: La descarga de residuo tipo 1 debe ser en n-1 y tipo 2 en n
for t = 1:T
    for k = 1:K
        for r = 1:R
            % Cantidad total de residuo tipo r recogido en todos los nodos por el vehículo k
            suma_recogida = 0;
            for i = 2:I
                suma_X_from_i = 0;
                for v = 1:n
                    if i ~= v
                        suma_X_from_i = suma_X_from_i + X_uvkt(i,v,k,t);
                    end 
                end
                suma_recogida = suma_recogida + suma_X_from_i * Z_urt(i,r,t);
            end
            
            % Si r=1, todo va a n-1; si r=2, todo va a n
            if r == 1
                Z.Constraints.(['R17C_rkt',num2str(r),num2str(k),num2str(t)]) = ...
                    w_urkt(n-1,r,k,t) == suma_recogida;
            else % r == 2
                Z.Constraints.(['R17C_rkt',num2str(r),num2str(k),num2str(t)]) = ...
                    w_urkt(n,r,k,t) == suma_recogida;
            end
        end
    end
end
disp('R17C_rkt');

% Modificar la restricción R12A para asegurar que la carga refleje correctamente la demanda recogida:
for t = 1:T
    for u = 2:n
        for k = 1:K
            for r = 1:R;
                suma_Xv = 0;
                for v = 1:n
                    suma_Xv = suma_Xv + X_uvkt(u,v,k,t);
                end
                % Si el vehículo k visita el nodo u, debe manejar su demanda
                Z.Constraints.(['R12A_ukt',num2str(u),num2str(r),num2str(k),num2str(t)]) = suma_Xv*Z_urt(u,r,t) <= w_urkt(u,r,k,t);
            end
        end
    end
end
disp('R12A_ukt')

% 1. Agregar restricción que verifica si la demanda total excede la capacidad de un vehículo
% y activa automáticamente más vehículos según sea necesario

% R18: Restricción para asegurar que se activan suficientes vehículos para cubrir la demanda total
for t = 1:T
    % Calcular la demanda total para el período t (en volumen)
    demanda_total = 0;
    for i = 2:I
        for r = 1:R
            if Z_urt(i,r,t) > 0
                demanda_total = demanda_total + Z_urt(i,r,t);
            end
        end
    end
    
    % Vector para almacenar la contribución de capacidad por cada vehículo activo
    capacidad_acumulada = optimexpr(K);

for k = 1:K
    capacidad_acumulada(k) = Q_k(k) * y_kt(k,t);
end
    
    % La suma de capacidades de vehículos activos debe ser suficiente para cubrir la demanda total
    Z.Constraints.(['R18_t',num2str(t)]) = sum(capacidad_acumulada) >= demanda_total;
end
disp('R18_t');

% Solver
sol = solve(Z,'Options',optimoptions('intlinprog','ConstraintTolerance',1e-06, 'Display', 'iter'));
disp(sol.X_uvkt);
disp(sol.w_urkt);
disp(sol.y_kt);
% disp(Z.Constraints);