How can I determine the amount of times a certain value can be achieved by summing values in a matrix?

I have a vector of up to 50 whole numbers, all of which are length values for wall panels. The wall panels come in sections of 144 inches, so my goal is to apply the lengths from the matrix and determine the number of full wall panels I would need. For example:

LengthVal = [50 50 60 70]; %User Input

With the values above, I know I would ned two 144 in panels: the 70 and 60 on one with the 50 and 50 on the other. I want to do this programetrically with a much larger number of length values. I assume the solution will require several for loops, but I am unsure about how to implement that. Are there any short cuts that you may know to perform a similar function?

Edit: Essentially, I want to get something with similar functionality to this website.

Thanks for any help!

[Walter Roberson]

Your Title suggests that you also want to produce all of the solutions, not just a minimal solution ? For example you would also want the solution [50] [50] [60] [70] (throwing away the rest of each panel) ??

[Kevin Johnston]

Apologies, I should have been more clear. Yes, I am just looking for the minimal solution.

[Cris LaPierre]

What defines a minimal solution? Why is [50 50] and [60 70] preferred over say [50 60] and [50 70]?

[Walter Roberson]

Is the requirement the minimum number of pieces used and any solution that achieves that is acceptable? Is the requirement the minimum total waste? Is the requirement that assuming that you are using the minimum number of pieces, that you want to maximize the waste width -- because the larger the waste pieces the more likely that you could use one of them for a later job?

Answer 1

[Matt J]

This solution (EDITED) requires the Optimization Toolbox:

LengthVal = [50,50 ,60,70]; %User Input


N=numel(LengthVal);

X=optimvar('X',[N,N],'LowerBound',0,'UpperBound',1,'type','integer');
y=optimvar('y',[1,N],'LowerBound',0,'UpperBound',1,'type','integer');

prob=optimproblem('Objective',sum(y));

prob.Constraints.row=sum(X,2)==1;   %assign exactly 1 panel to each wall
prob.Constraints.capacity=LengthVal*X<=144; %respect panel capacity
prob.Constraints.ylb=sum(X,1)/N<=1.1*y;
prob.Constraints.yub=sum(X,1)>=y;

sol=solve(prob);

[~,~,groups]=unique((1:N)*round(sol.X.'),'stable') %the result, as a grouping vector

groups =

     1
     2
     1
     2

[Bruno Luong]

Neat!!!

[Matt J]

Thanks, but could you explain to me why it works? I now realize that the constraint,

prob.Constraints.row=sum(X,2)==1; 

ensures that all feasible solutions will have the same objective function value, N. No solution is more optimal than any other, and yet the solver finds the one we want....

[Matt J]

I think I fixed it (see my edited answer), but I still find it an interesting mystery that the original draft below ever worked!

LengthVal = [50 50 60 70]; %User Input

N=numel(LengthVal);

X=optimvar('X',[N,N],'LowerBound',0,'UpperBound',1,'type','integer');

prob=optimproblem('Objective',sum(X(:)));

prob.Constraints.row=sum(X,2)==1;   %each wall must be completely covered by one of the panels
prob.Constraints.col=LengthVal*X<=144; %length limit each panel can accomodate

sol=solve(prob);

[~,~,groups]=unique((1:N)*round(sol.X.'),'stable') %the result, as a grouping vector

[Bruno Luong]

Yeah now you have mentionned it.

I have impression that the "right" solution is found because that INTLINPROG translates to the dual problem, which can be understood for a human being as "try to pack as much paper walls into full-length one" and miraclly it gives the right result.

Not a proof but this code shows that when we set algorithm to primal-simplex then it sucks, the result is wrong:

LengthVal = [50 50 60 70]; %User Input
N=numel(LengthVal);
X=optimvar('X',[N,N],'LowerBound',0,'UpperBound',1,'type','integer');
prob=optimproblem('Objective',sum(X(:)));
prob.Constraints.row=sum(X,2)==1;   %each wall must be completely covered by one of the panels
prob.Constraints.col=LengthVal*X<=144; %length limit each panel can accomodate

problem = prob2struct(prob);

problem.options = optimoptions('intlinprog','RootLPAlgorithm','dual-simplex');
x = intlinprog(problem);
X = round(reshape(x,[N,N]));
[~,~,groups]=unique((1:N)*X.','stable') % Seems OK

problem.options = optimoptions('intlinprog','RootLPAlgorithm','primal-simplex');
x = intlinprog(problem);
X = round(reshape(x,[N,N]));
[~,~,groups]=unique((1:N)*X.','stable') % fails

Now the REAL question is "does the dual-simplex really works"? EDIT: NO

Not sure but it seems giving reasonable solution and pretty fast (your last modification makes the algo turn foreever for N=50).

But when the outcome depends on detail of the algorithm and we don't access, it's hard to make any claim.

[Bruno Luong]

Even the dual doesn't give the best result as showed in this example:

LengthVal = [36    55   112    68    64    22   129    24    35    37];
N=numel(LengthVal);

X=optimvar('X',[N,N],'LowerBound',0,'UpperBound',1,'type','integer');
prob=optimproblem('Objective',sum(X(:)));
prob.Constraints.row=sum(X,2)==1;   %each wall must be completely covered by one of the panels
prob.Constraints.col=LengthVal*X<=144; %length limit each panel can accomodate
sol=solve(prob);
X1 = round(sol.X);
[~,~,groups1]=unique((1:N)*X1','stable'); %the result, as a grouping vector

X=optimvar('X',[N,N],'LowerBound',0,'UpperBound',1,'type','integer');
y=optimvar('y',[1,N],'LowerBound',0,'UpperBound',1,'type','integer');
prob=optimproblem('Objective',sum(y));
prob.Constraints.row=sum(X,2)==1;   %assign exactly 1 panel to each wall
prob.Constraints.capacity=LengthVal*X<=144; %respect panel capacity
prob.Constraints.ylb=sum(X,1)/N<=y;
sol=solve(prob);
X2 = round(sol.X);
[~,~,groups2]=unique((1:N)*X2','stable'); %the result, as a grouping vector

np1 = sum(sum(X1,1)>0) % 6
np2 = sum(sum(X2,1)>0) % 5
clf
subplot(1,2,1); imagesc(X1);
subplot(1,2,2); imagesc(X2);

[Matt J]

I was missing a constraint

prob.Constraints.yub=sum(X,1)>=y;

[Bruno Luong]

Sorry I can't see the need yub.

Alternative is just create a bunch of constraints (at least this formulation has better lisibility for other persons who read):

X=optimvar('X',[N,N],'LowerBound',0,'UpperBound',1,'type','integer');
y=optimvar('y',[1,N],'LowerBound',0,'UpperBound',1,'type','integer');
prob=optimproblem('Objective',sum(y));
prob.Constraints.row=sum(X,2)==1;   %assign exactly 1 panel to each wall
prob.Constraints.capacity=LengthVal*X<=144; %respect panel capacity
for j=1:N
    for i=1:N
        prob.Constraints.(sprintf('Bnd_i%dj%d',i,j))=X(i,j)<=y(1,j);
    end
end
sol=solve(prob);
X3 = round(sol.X);

[Matt J]

We want y=any(X,1).

ylb forces y(i) to 1 when any(X(:,1))=1

yub forces y(i) to 0 when X(:,i) contains only zeros.

[Bruno Luong]

It looks redundant with

y=optimvar('y',[1,N],'LowerBound',0,'UpperBound',1,'type','integer')

[Matt J]

Consider this hypothetical X,y pair with N=5 and LengthVal(:)=1. It satisfies all the constraints except yub, but y(2) is not equal to any(X(:,2)) as we desire.

X =

     1     0     0     0     0
     1     0     0     0     0
     1     0     0     0     0
     1     0     0     0     0
     1     0     0     0     0


y =

     1     1     0     0     0

[Bruno Luong]

Yes but the objective is minimizing sum(y) therefore INTLINPROG must return

y = [1 0 0 0 0]

which also satisfies the constraints, or do I miss something?

[Matt J]

Ok, so you are saying that yub will inherently be enforced by the minimization. We can be sure it will be satisfied at the minimum even without enforcing it explciitly...

[Bruno Luong]

Exactly.

[Matt J]

Very well, but I wonder if it is good to have as an explicit constraint anyway to help the algorithm narrow the search...

[Bruno Luong]

In my experience no.

By removing the redundancy of the constraints, we'll help the minimizer for not doing extra calculation for when it needs to update the active set.

Answer 2

[Bruno Luong]

I put here the code with limiting CPU time so one can run up to 50 samples. Matt should get credit for his original idea.

fulllength = 144;
% Random LengthVal
LengthVal = randi(fulllength,1,50)
N=numel(LengthVal);

X=optimvar('X',[N,N],'LowerBound',0,'UpperBound',1,'type','integer');
y=optimvar('y',[1,N],'LowerBound',0,'UpperBound',1,'type','integer');
prob=optimproblem('Objective',sum(y));
prob.Constraints.row=sum(X,2)==1;   %assign exactly 1 panel to each wall
prob.Constraints.capacity=LengthVal*X<=fulllength; %respect panel capacity
for j=1:N
    for i=1:N
        prob.Constraints.(sprintf('Bnd_i%dj%d',i,j))=X(i,j)<=y(1,j);
    end
end
problem = prob2struct(prob);
problem.options = optimoptions('intlinprog','MaxTime',60); % limits the computation to 60s
x = intlinprog(problem);

% Post process data
X = reshape(x(1:N^2),[N,N]);
X = round(X);
X(:,all(X==0,1)) = [];
ng = size(X,2);
FillLengths=LengthVal*X;
L = X.*LengthVal(:);
N = cumsum(X,1).*X;
[~,g,n] = find(N);
data = accumarray([g,n],L(X>0),[],[],NaN);

% Graphic output
figure(1);
clf

subplot(1,2,1);
imagesc(X);
xlabel(sprintf('group (1-%d)', ng));
ylabel('Wall paper number');

subplot(1,2,2);
bar(1:ng,data,'stacked');
xlim([0.5 ng+0.5]);
ylim([0 fulllength]);
xlabel(sprintf('group (1-%d)', ng));
ylabel('stack lengths');

Examples of results:

Answer 3

[Matt J]

Below is a more advanced version of my original answer, which is showing greater empirical succes. It incorporates initialization and restart strategies to iterative reduce problem dimension, and hence the run time. On 3 consecutive sets of randomized inputs with N=50,

LengthVal = (randi(10,1,50)*10);

it was able to complete the optimization within a few minutes. However, for some cases, it does still fall into a protracted branch and bound phase, so I have chosen an empirical cap for the MaxNodes option parameter:

N=numel(LengthVal);
e=ones(N,1);

sol=init(LengthVal);
M=sum(sol.y);

tic;
while 1
    
    X=optimvar('X',[N,M],'LowerBound',0,'UpperBound',1,'type','integer');
    y=optimvar('y',[1,M],'LowerBound',0,'UpperBound',1,'type','integer');
    
    prob=optimproblem('Objective',sum(y));
    
    prob.Constraints.rowLB=sum(X,2)>=0.5;   %assign exactly 1 panel to each wall
    prob.Constraints.rowUB=sum(X,2)<=1.5;
    prob.Constraints.capacity=LengthVal*X<=144; %respect panel capacity
%     prob.Constraints.ylb=sum(X,1)/N<=1.1*y;
    prob.Constraints.ylb=X/2<=e*y;
    %prob.Constraints.yub=sum(X,1)+0.5>=y;
    prob.Constraints.mono=sum(X(:,1:M-1),1)>=sum(X(:,2:M),1);
    
    of=@(x,o,s) customFcn(x,o,s,M);
    
    cut=logical(sol.y);
    sol.X=sol.X(:,cut); 
    [~,is]=sort(sum(sol.X,1),'descend');
    sol.X=sol.X(:,is);
    sol.y=sol.y(cut);
    
    
    [sol,fval]=solve(prob,sol,...
                     'options',optimoptions(@intlinprog,'OutputFcn',of,...
                     'MaxNodes',M*20000));

    
    fval=round(fval),
    if fval==M, break; end
    
    M=fval;
    

end
toc;
    
[~,~,groups]=unique((1:M)*round(sol.X.'),'stable'); %the result, as a grouping vector
 grous=groups(:).',

function stop = customFcn(x,optimValues,state,M)
   stop=false;
   if ~isempty(optimValues.fval) && ~isempty(x)
     stop=optimValues.fval<M;
   end
end

function sol=init(LengthVal)

 N=numel(LengthVal);
 sol.X=zeros(N);
 j=1; accum=0;

 
  for i=1:N
      
     accum=accum+LengthVal(i); 
     if accum<=144
      sol.X(i,j)=1;
     else
        accum=LengthVal(i); j=j+1;
        sol.X(i,j)=1;
     end
  end
  
  sol.y=double(any(sol.X,1));
  
end