Mixed-integer Linear Programming with double sum constraints and 3D optimization variable

Hello,

I am trying to implement the below MILP problem. I did not add all of the constraints to not to complicate my question. P and W are integer decision variables (Mx1).

s.t.

Below is my code to implement. I tried to use MATLAB's optimization language but I couldn't use it correcly so I used an inefficient for loop. Commented alternative MATLAB script on top of each constraint.

data = someFun2ImportData('data.csv');
len = length(data);
% Pairwise distance matrix
dist_matrix = squareform(pdist([data(:,2:3)],@lldistmile));

P = optimvar('P',len,'Type','integer','LowerBound',0,'UpperBound',1);
W = optimvar('w', len,len,len, 'Type', 'integer', 'LowerBound',0,'UpperBound',1);

beta = 150;
t = 5;

% Constraints
for i=1:len
    const1 = sum(W(:,:,i)) == 1;
    ndesign.Constraints.const1 = const1;
    
    % const2 = sum(W,2) <= P;
    const2 = sum(W(:,i,:)) <= P(i);
    ndesign.Constraints.const2 = const2;
    
    %const3 = sum(W*dist_matrix/beta,[1 2])*60 <= t;
    for k=1:len
        const3 = sum(W(i,k,i).*dist_matrix(i,k))/beta*60 <= t;
        sum(W(i,k,i).*dist_matrix(i,k))/beta*60
        ndesign.Constraints.const3 = const3;
    end
    

end

ndesign = optimproblem; % Create optimization problem

objfun = sum(P); % Objective Function

ndesign.Objective = objfun; % Assign objfun to the problem

opts = optimoptions('intlinprog','Display','off','PlotFcn',@optimplotmilp);
[sol,fval,exitflag,output] = solve(networkdesign,'options',opts);

I need your help to implement this problem.

Answer 1

[Matt J]

data = someFun2ImportData('data.csv');
len = length(data);
% Pairwise distance matrix
dist_matrix = squareform(pdist([data(:,2:3)],@lldistmile));


P = optimvar('P',[len,1],'Type','integer','LowerBound',0,'UpperBound',1);
W = optimvar('w', [len,len,len,] , 'Type', 'integer', 'LowerBound',0,'UpperBound',1);

beta = 150;
t = 5;

ndesign = optimproblem('Objective',sum(P)); % Create optimization problem

 e=ones(1,len);
 ndesign.Constraints.const1 = ( sum(reshape(W,[],len),1)==1  );
 ndesign.Constraints.const2 = ( squeeze(sum(W,2))<=P*e  );
 ndesign.Constraints.const3 = ( dist_matrix(:).'*reshape(W,[],len)<=beta*t   );
 
opts = optimoptions('intlinprog','Display','off','PlotFcn',@optimplotmilp);
[sol,fval,exitflag,output] = solve(ndesign,'options',opts);

[Emin Burak Onat]

Thanks a lot Matt! You were really helpful and fast!