cooperative logistics - profit maximization

I need help about this code dealing with a problem of cooperative logistics. I have to maximize profits, that is, demonstrate that with the three carriers together you get a profit greater than the sum of the profits given by the three carriers taken individually. MatLab takes over 30 minutes and gives infeasible. Thank you so much!

clc
clear all

%%%%%%%%%%%%%%%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%

Nr=20; % number of trips 
Nc=3; % number of carriers
Vr=22; % number of vehicles

alpha=1; % extra charge on revenue;
cE= 0.2; 
beta=0.8;
mu=80;

% dummy quantity
M=1000000;

% data about trips 20
d= [85, 80, 123, 111, 153, 112, 150, 107, 181, 135, 109, 106, 113, 148, 164, 127, 148, 141, 108, 118];
delta=[200, 200, 200, 200, 200, 200, 200, 200, 300, 300, 300, 300, 300, 300, 300, 400, 400, 400, 400, 400];
p=[25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25, 25];

% data about carriers, INVENTATI
I=[3, 3, 3];
cc=[0.8, 1.0, 0.9];


% data about vehicles, 22
T= [8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8];
CAP= [30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30, 30];
Eempty= [0.3, 0.2, 0.3, 0.3, 0.2, 0.3, 0.3, 0.2, 0.3, 0.2, 0.2, 0.3, 0.2, 0.2, 0.2, 0.2, 0.3, 0.2, 0.2, 0.3, 0.3, 0.2];
Efull= [0.3, 0.4, 0.3, 0.3, 0.3, 0.4, 0.3, 0.4, 0.3, 0.4, 0.4, 0.4, 0.4, 0.4, 0.3, 0.3, 0.3, 0.4, 0.4, 0.3, 0.4, 0.4]; % it is already Efull-Eempty
IsE=[0, 0, 0, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0]; %1 if the vehicle is electric
Owner= [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3];

% repositioing distances between trips 20x20
eps= [14	16	6	19	26	27	21	10	16	18	29	7	30	27	28	10	8	11 17 10;
    27	9	20	11	9	27	22	18	27	14	29	22	5	25	9	21	17	18	11	16;
    25	11	10	30	19	14	7	9	19	26	27	22	29	7	12	29	24	7	9	30;
    13	16	6	19	7	7	24	5	27	28	6	18	27	9	17	30	24	14	18	10;
    22	23	21	13	16	16	25	21	14	8	14	13	18	17	11	7	16	29	6	6;
    5	21	9	17	26	30	27	11	24	28	16	6	25	22	27	11	18	30	13	7;
    13	11	15	24	20	30	26	11	7	21	7	17	29	5	11	12	22	22	8	5;
    28	7	14	20	7	25	20	12	15	5	27	10	5	27	20	26	22	25	14	28;
    16	19	11	12	8	25	24	9	26	26	17	15	8	5	20	19	7	21	28	24;
    29	8	13	18	17	13	15	12	21	11	10	25	17	29	22	11	28	27	14	13;
    10	22	26	16	22	19	7	28	10	14	12	28	27	10	12	6	23	30	28	26;
    18	11	15	25	20	7	21	17	21	9	16	5	11	5	18	16	17	26	21	19;
    14	15	25	19	11	16	8	25	21	9	11	8	13	5	10	7	5	28	30	8;
    26	15	11	19	21	23	9	24	9	26	16	17	9	27	25	27	20	21	15	13;
    27	27	14	19	28	9	22	24	18	11	22	15	15	26	27	27	27	19	17	19;
    25	11	17	24	6	30	15	23	21	30	16	5	18	25	12	20	9	11	25	27;
    29	9	17	11	24	8	29	28	21	16	20	15	6	28	8	22	7	13	26	22;
    23	30	28	15	11	20	13	10	25	18	29	13	12	11	16	20	24	11	6	30;
    28	23	12	29	16	29	26	22	28	17	20	8	30	22	8	25	12	7	13	9;
    7	7	22	11	23	9	17	25	27	26	21	29	23	14	8	6	6	26	24	30];

% 1 if the pair of trips can be combined INVENTATO 20x20
O=[0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
    0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
    0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
    0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
    0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
    0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
    0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
    0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
    0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
    0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1;
    1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0;
    1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0];
%O=[0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 1;
%1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0];

% distances between the depot of vehicle v (column) and 
% the origin of trip n (row) 20x22 dati mancanti metto 1000 da 1 a 10
dO= [3 8 4 2 2 4 3 9 3 4 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    9 3 5 7 8 2 5 9 1 4 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    8 9 5 1 3 9 5 6 8 5 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    5 8 2 2 9 2 8 2 7 7 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    1 6 9 5 6 5 6 3 1 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    1 5 7 7 4 8 1 5 2 3 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    1 3 1 5 8 2	6 1 5 2 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    5 8 7 6 3 7	7 9 1 6 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 3 8 8 2 7 7 7 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 7 7 8 2 8 3 9 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 6 7 2 3 3 6 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 6 8 5 4 2 8 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1 2 6 5 6 7 7 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 9 1 2 7 8	8 5 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 8 9 8 9 8 9 5 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 8 6 7 9 8;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1 5 5 7 9;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 1 8 6 4;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 4 1 3 8;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 8 9 7 3];

% distances between the destination of trip n (row)
% and the depot of vehicle v (column) 20x22 dati mancanti metto 1000 da 1 a 10

dD=[7 5 1 3 9 1 5 2 1 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    7 2 7 3 3 1 5 8 1 1 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    1 7 2 5 3 5 8 4 7 8 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    7 7 4 8 3 7 3 1 7 1 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    5 1 1 6 7 2 4 9 9 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    6 7 1 1 5 2 4 3 9 5 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    8 7 2 9 5 6 6 5 8 9 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    6 8 9 7 6 8 2 5 8 8 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 7 7 6 3 9 1 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 9 7 3 6 9 6 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 3 4 6	6 5 9 9 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 7	7 1 7 2 1 7 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 6 7 9 8 6 2 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 4 1 5 4 1 3 3 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 6	6 6 4 6 9 5 1000 1000 1000 1000 1000;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1 5 6 7 2;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 9 9 6 8 2;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 2 9 4 6 2;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 5 5 4 9 1;
    1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 1000 4 8 6 9 9];

% operative cost 3x22
c=[0.7, 0.8, 0.8, 0.6, 0.9, 0.8, 0.6, 0.9, 0.6, 0.8, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100;	
    100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 0.9, 0.8, 0.8, 1, 0.7, 0.7, 0.8, 100, 100, 100, 100, 100;
    100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 0.9, 0.9, 0.6, 0.7, 0.7];

% 1 if the trip (row) is owned by the carrier (column) INVENTATO 20x3 
z=[1 0 0;
    1 0 0;	
    1 0 0;
    1 0 0;
    1 0 0;
    1 0 0;	
    1 0 0;
    1 0 0;
    0 1 0;
    0 1 0;	
    0 1 0;
    0 1 0;
    0 1 0
    0 1 0;
    0 1 0;
    0 0 1;
    0 0 1;
    0 0 1;	
    0 0 1;
    0 0 1];

% initial profit of carriers INVENTATO dati ottenuti dal prob1 per ciascun
% carrier
S0=  [3380  2382 1561];
%S0 = [1000 1000 1000];

TOA=zeros(Nr,Vr);

% computation of TOA(n,v) for single trips;
for n=1:Nr
    for v=1:Vr
        TOA(n,v) = dO(n,v)+d(n)+dD(n,v);
    end
end

E=zeros(Nr,Vr);
% computation of cost E;
for n=1:Nr
    for v=1:Vr
        E(n,v)= TOA(n,v)*Eempty(v) + (Efull(v)*d(n)*p(n)/CAP(v));
    end
end

TOAC=zeros(Nr,Nr,Vr);
% computation of TOAC(n,k,v) for pairs of trips;
for n=1:Nr
    for k=1:Nr
        for v=1:Vr
            TOAC(n,k,v) = dO(n,v)+d(n)+d(k)+eps(n,k)+dD(k,v);
        end
    end
end

EC=zeros(Nr,Nr,Vr);
% computation of cost E for combined trips;
for n=1:Nr
    for k=1:Nr
        for v=1:Vr
            EC(n,k,v)= TOAC(n,k,v)*Eempty(v) + (Efull(v)*(d(n)*p(n)+d(k)*p(k))/CAP(v));
        end
    end
end


% % %%%%%%%%%%%%%%%% MAX %%%%%%%%%%%%%%%%%%%%%
% 
%TO BE NOTED: I have not included the carrier in decision variables 
% because here we are not includeing vehicles in the data structrue related
% to carriers but we are considering the whole list of vehicles
% than we associate vehicles with carriers by using the parameters Owner(v)


% decision variables with lower and upper bounds
x = optimvar('x',Nr,Vr,'Type','integer','LowerBound',0, 'UpperBound',1);
y = optimvar('y',Nr,Nr,Vr,'Type','integer','LowerBound',0, 'UpperBound',1);

% problema definition
prob = optimproblem('ObjectiveSense','max');

% objective function computation

obj=0;
% %profit S1 of carriers;

for r=1:Nc
    for n=1:Nr
        for v=1:Vr
            if (Owner(v) == r)
                obj=obj+(alpha* d(n)*I(r) - TOA(n,v)*c(r,v)-(1-IsE(v))*beta*cE*E(n,v))*x(n,v);
            end
        end
    end
end
% 
%profit S2 of carriers;

for r=1:Nc
    for n=1:Nr
        for k=1:Nr 
            if (k ~= n)
                for v=1:Vr
                    if (Owner(v) == r)
                        obj= obj+ (alpha*(d(n)+d(k))*I(r)-TOAC(n,k,v)*c(r,v)-(1-IsE(v))*beta*cE*EC(n,k,v))*y(n,k,v);
                    end
                end
            end
        end
    end
end

%profit S3 of carriers;

for r=1:Nc
    %compensation profits and costs for single trips
    for n=1:Nr
        for rr=1:Nc
            if (rr ~= r)
                for v=1:Vr
                    if (Owner(v) == rr)
                        obj=obj+ delta(n)*z(n,r)*x(n,v);
                    else
                        if (Owner(v) == r)
                            obj=obj- delta(n)*z(n,rr)*x(n,v);
                        end
                    end
                end
            end
        end
    end
    %compensation profits and costs for combined trips
    for n=1:Nr
        for k=1:Nr
            for rr=1:Nc
                if (rr ~= r)
                    for v=1:Vr
                        if (Owner(v) == rr)
                            obj = obj+ y(n,k,v)*(delta(n)*z(n,r)+delta(k)*z(k,r));
                        else
                            if (Owner(v) == r)
                                obj = obj- y(n,k,v)*(delta(n)*z(n,rr)+delta(k)*z(k,rr));
                            end   
                        end
                    end
                end
            end
        end
    end
end

%

prob.Objective = obj;

%%% CONSTRAINTS

% profits higher than initial profits

constr_profit = optimineq(Nc,1);
% 
for r=1:Nc
    S=0;
    %S1 computation
    for n=1:Nr
        for v=1:Vr
            if (Owner(v) == r)
                S=S+(alpha* d(n)*I(r) - TOA(n,v)*c(r,v)-(1-IsE(v))*beta*cE*E(n,v))*x(n,v);
            end
        end
    end
    
    
    %S2 computation
    for n=1:Nr
        for k=1:Nr 
            if (k ~= n)
                for v=1:Vr
                    if (Owner(v) == r)
                        S= S + (alpha*(d(n)+d(k))*I(r)-TOAC(n,k,v)*c(r,v)-(1-IsE(v))*beta*cE*EC(n,k,v))*y(n,k,v);
                    end
                end
            end
        end
    end
    
    
    % S3 computation;
    
    %compensation profits and costs for single trips
    for n=1:Nr
        for rr=1:Nc
            if (rr ~= r)
                for v=1:Vr
                    if (Owner(v) == rr)
                        S=S+ delta(n)*z(n,r)*x(n,v);
                    else
                        if (Owner(v) == r)
                            S=S- delta(n)*z(n,rr)*x(n,v);
                        end
                    end
                end
            end
        end
    end
    %compensation profits and costs for combined trips
    for n=1:Nr
        for k=1:Nr
            for rr=1:Nc
                if (rr ~= r)
                    for v=1:Vr
                        if (Owner(v) == rr)
                            S = S+ y(n,k,v)*(delta(n)*z(n,r)+delta(k)*z(k,r));
                        else
                            if (Owner(v) == r)
                                S = S- y(n,k,v)*(delta(n)*z(n,rr)+delta(k)*z(k,rr));
                            end   
                        end
                    end
                end
            end
        end
    end
    constr_profit(r)=S>=S0(r);
end

prob.Constraints.constr_profit = constr_profit;


% all demand is executed (14);

constr_demand_sat=optimeq(1,1);
% 
for n=1:Nr
    temp=0;
    for v=1:Vr
        temp=temp+x(n,v);
        for k=1:Nr
            if (k ~= n)
                temp=temp+y(n,k,v)+y(k,n,v);
            end
        end
    end
    constr_demand_sat(n)=temp==1;
end

prob.Constraints.constr_demand_sat = constr_demand_sat;

% Every vehicle must be assigned to at most one trip;

constr_assign=optimineq(Vr,1);

for v=1:Vr
    temp=0;
    for n=1:Nr
        temp=temp+x(n,v);
        for k=1:Nr
            if (k ~= n)
                temp=temp+y(n,k,v);
            end
        end
        
    end
    constr_assign(v)=temp<=1;
end

prob.Constraints.constr_assign = constr_assign;

% duration of the batteries of electric trucks;

constr_batteries_single=optimineq(Nr,Vr);

for n=1:Nr
    for v=1:Vr
        if (IsE(v) == 1)
            constr_batteries_single(n,v)=x(n,v)*TOA(n,v)<=T(v)*mu;
        else
            constr_batteries_single(n,v)=x(n,v)*TOA(n,v)<=M;
        end
    end  
end

prob.Constraints.constr_batteries_single = constr_batteries_single;

constr_batteries_combined=optimineq(Nr,Nr,Vr);

for n=1:Nr
    for k=1:Nr
        if (k ~= n)
            for v=1:Vr
                if (IsE(v) == 1)
                    constr_batteries_combined(n,k,v)=y(n,k,v)*TOAC(n,k,v)<=T(v)*mu;
                else
                    constr_batteries_combined(n,k,v)=y(n,k,v)*TOAC(n,k,v)<=M;
                end
            end
        end
    end
end

prob.Constraints.constr_batteries_combined = constr_batteries_combined;

% constraints on trips combination;

constr_combined=optimineq(Nr,Nr,Vr,1);

for n=1:Nr
    for k=1:Nr
        for v=1:Vr
            constr_combined(n,k,v)=y(n,k,v)<=O(n,k);
        end
    end
end

prob.Constraints.constr_combined = constr_combined;

%  
% % problem solution
show(prob)
[sol,val] = solve(prob);

% retrieve of optimal values
opt_variables_x=sol.x;
opt_variables_y=sol.y;

opt_obj=val;

% profits higher than initial profits


% %SS=zeros(Nc,1);
% %for r=1:Nc
%   SS(r)=0;
%   %S1 computation
%   for n=1:Nr
%         for v=1:Vr
%             if (Owner(v) == r)
%               SS(r)=SS(r)+(alpha* d(n)*I(r) - TOA(n,v)*c(r,v)-(1-IsE(v))*beta*cE*E(n,v))*sol.x(n,v);
%             end
%         end
%   end
% 
% 
%   %S2 computation
%   for n=1:Nr
%     for k=1:Nr 
%       if (k ~= n)
%          for v=1:Vr
%             if (Owner(v) == r)
%                  SS(r)= SS(r) + (alpha*(d(n)+d(k))*I(r)-TOAC(n,k,v)*c(r,v)-(1-IsE(v))*beta*cE*EC(n,k,v))*sol.y(n,k,v);
%             end
%          end
%       end
%     end
%   end
% 
% 
% % S3 computation;
% 
%   %compensation profits and costs for single trips
%   for n=1:Nr
%     for rr=1:Nc
%       if (rr ~= r)
%         for v=1:Vr
%           if (Owner(v) == rr)
%               SS(r)=SS(r)+ delta(n)*z(n,r)*sol.x(n,v);
%           else
%             if (Owner(v) == r)
%               SS(r)=SS(r)- delta(n)*z(n,rr)*sol.x(n,v);
%             end
%           end
%         end
%       end
%     end
%   end
%   %compensation profits and costs for combined trips
%   for n=1:Nr
%     for k=1:Nr
%       for rr=1:Nc
%         if (rr ~= r)
%           for v=1:Vr
%             if (Owner(v) == rr)
%               SS(r) = SS(r)+ sol.y(n,k,v)*(delta(n)*z(n,r)+delta(k)*z(k,r));
%             else
%                if (Owner(v) == r)
%                  SS(r) = SS(r)- sol.y(n,k,v)*(delta(n)*z(n,rr)+delta(k)*z(k,rr));
%                end   
%             end
%           end
%          end
%       end
%     end
%   end
% 
%   y2=0;
%   for i=1:Nr
%       for v=1:Vr
%           y2=y2+opt_variables_y(i,2,v)+opt_variables_y(2,i,v);
%       end;
%   end;
% 
% end

----

[Torsten]

I suggest reducing the size of the problem first to find the problematic constraints leading to infeasibility.

What do you get if you use

solvers(prob)

before invoking

[sol,val] = solve(prob);

?