Untitled

%% Set
clc

N = 8;
M = 5;

%% Parametri
BigM = 1000;

p = [5 9 6 1 7 2 8 3;
     1 1 2 2 6 6 9 2;
     3 7 7 1 9 2 10 6;
     9 6 6 5 1 6 10 6;
     1 3 1 9 10 1 6 8];

%% Variabili Decisionali
x = optimvar('x', N, N, M, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);
y = optimvar('y', N, 'Type', 'integer', 'LowerBound', 0, 'UpperBound', 1);
s = optimvar('s', N, M, 'LowerBound', 0);
c = optimvar('c', N, M, 'LowerBound', 0);
C = optimvar('C', 1, 'LowerBound', 0);

%% Funzione Obiettivo
prob = optimproblem('Objective', C, 'ObjectiveSense', 'min');

%% Vincoli 

% 1° e 2° gruppo sequenziamento dei job
k = 1;
vincolo1 = optimconstr(N*(N-1)*M);
vincolo2 = optimconstr(N*(N-1)*M);

for i = 1:N
    for j = 1:N
        if i ~= j
            for m = 1:M
                if m == 2 || m == 4
                    vincolo1(k) = s(i, m) >= c(j, m) - BigM * x(i, j, m) - BigM * y(i);
                    vincolo2(k) = s(j, m) >= c(i, m) - BigM * (1 - x(i, j, m)) - BigM * y(i);
                elseif m == 5
                    vincolo1(k) = s(i, m) >= c(j, m) - BigM * x(i, j, m);
                    vincolo2(k) = s(j, m) >= c(i, m) - BigM * (1 - x(i, j, m));
                else
                    vincolo1(k) = s(i, m) >= c(j, m) - BigM * x(i, j, m) - BigM * (1 - y(i));
                    vincolo2(k) = s(j, m) >= c(i, m) - BigM * (1 - x(i, j, m)) - BigM * (1 - y(i));
                end
                k = k + 1;
            end
        end
    end
end

prob.Constraints.constr1 = vincolo1;
prob.Constraints.constr2 = vincolo2;

% 3° gruppo 
k = 1;
vincolo3 = optimconstr(N * M);

for i = 1:N
    for m = 1:M
        if m == 2 || m == 4
            vincolo3(k) = c(i, m) >= s(i, m) + p(m, i) - BigM * y(i);
        elseif m == 5
            vincolo3(k) = c(i, m) >= s(i, m) + p(m, i);
        else
            vincolo3(k) = c(i, m) >= s(i, m) + p(m, i) - BigM * (1 - y(i));
        end
        k = k + 1;
    end
end

prob.Constraints.constr3 = vincolo3;

% 4° gruppo
k = 1;
vincolo4 = optimconstr(4 * N);

for i = 1:N
    vincolo4(k) = s(i, 3) >= c(i, 1) - BigM * (1 - y(i)); % precedenza di operazioni dei job i-esimi sulle macchine 
    k = k + 1;
    vincolo4(k) = s(i, 4) >= c(i, 2) - BigM * y(i);
    k = k + 1;
    vincolo4(k) = s(i, 5) >= c(i, 4) - BigM * y(i);
    k = k + 1;
    vincolo4(k) = s(i, 5) >= c(i, 3) - BigM * (1 - y(i));
    k = k + 1;
end

prob.Constraints.constr4 = vincolo4;

% 5° gruppo
k = 1;
vincolo5 = optimconstr(N * M);

for i = 1:N
    for m = 1:M
        vincolo5(k) = C >= c(i, m); 
        k = k + 1;
    end
end

prob.Constraints.constr5 = vincolo5;

%% Risolvere il problema
opts = optimoptions('intlinprog', 'Display', 'iter');
sol = solve(prob, 'Options', opts);

disp(sol);