Unconstrained minimisation problem with a complicated range

I have 3N objects with properties p1,p2,p3. I need to organise these objects into groups of N objects such that the sum of property p1 is the same. I think I can do with the functional:

f(p1)=(sum(p1,1)-sum(p1,2)).^2+(sum(p1,1)-sum(p1,3)).^2+(sum(p1,3)-sum(p1,2)).^2

|Where sum(p1,i) denotes the sum over the subset of p1's for N objects. The same for other properties, so the objective functionals will simply add(I think). Is there a way of doing this? I guess, if I can do it for one property, I can do it for three?

[Torsten]

Will it always be possible to find an arrangement that makes f(p1) = 0 or might there be cases where you can only find a "best" arrangement p1* that minimizes f, but f(p1*) > 0 will hold ?

[Mat]

I don't think there will be any arrangement for which f(p1)=0.

[David Goodmanson]

Hi Mat, does p1 take on integer values whose sum must be identical for all three groups, or does the problem concern something else?

[Matt J]

Where sum(p1,i) denotes...

It already denotes something in Matlab, so please don't use it to denote something else. What are the independent variables here? Are p1,p2 and p3 each vectors of length N? Are they the design variables or are they functions of something else?

[Mat]

Matt J, I have 3*N objects and I want them arranged in three groups in such a way such that the functional is minimised.

[Steven Lord]

So you're trying to solve a bin packing problem?

[Mat]

@Steven Lord Is that what it's called?

Answer 1

[Torsten]

I'll formulate the problem for one property - a generalization to three properties is obvious.

I'll assume that each of the 3*N objects is assigned a number p1j which stands for the amount of property 1 in object j ( 1<= j <= 3*N).

This can be formulated as a mixed-integer linear optimization problem:

min: d1 + d2 + d3

s.t.

-d1 <= sum_{i=1}^{i=N} sum_{j=1}^{3*N} x_ij * p1j - sum_{i=N+1}^{i=2*N} sum_{j=1}^{3*N} x_ij * p1j

d1 >= sum_{i=1}^{i=N} sum_{j=1}^{3*N} x_ij * p1j - sum_{i=N+1}^{i=2*N} sum_{j=1}^{3*N} x_ij * p1j

-d2 <= sum_{i=1}^{i=N} sum_{j=1}^{3*N} x_ij * p1j - sum_{i=2*N+1}^{i=3*N} sum_{j=1}^{3*N} x_ij * p1j

d2 >= sum_{i=1}^{i=N} sum_{j=1}^{3*N} x_ij * p1j - sum_{i=2*N+1}^{i=3*N} sum_{j=1}^{3*N} x_ij * p1j

d3 <= sum_{i=N+1}^{i=2*N} sum_{j=1}^{3*N} x_ij * p1j - sum_{i=2*N+1}^{i=3*N} sum_{j=1}^{3*N} x_ij * p1j

-d3 >= sum_{i=N+1}^{i=2*N} sum_{j=1}^{3*N} x_ij * p1j - sum_{i=2*N+1}^{i=3*N} sum_{j=1}^{3*N} x_ij * p1j

sum_{i=1}^{3*N} x_ij = 1 for 1 <= j <= 3*N

sum_{j=1}^{3*N} x_ij = 1 for 1 <= i <= 3*N

0 <= x_ij <= 1 integer for 1 <= i,j <= 3*N

d1, d2, d3 >= 0

You can use MATLAB's "intlinprog" to solve.

[Mat]

The numbers in general aren't integers, and as part of the solution requires where certain objects are placed where.

[Torsten]

The 3Nx3N - matrix x_ij defines a permutation of the 3N objects.

The columns of row 1 to row N where the 1s are placed by the optimizer in the solution denote the first group of assembled objects, the columns of row N+1 to row 2N where the 1s are placed by the optimizer in the solution denote the second group of assembled objects and the column numbers of row 2N+1 to row 3N where the 1s are placed by the optimizer in the solution denote the third group of assembled objects.

The first 6 inequalities measure the difference of p1 between the three groups.

The objective function as the sum of absolute values of the differences of p1 between the three groups (not squared differences as in your functional to keep the problem linear) is to be minimized.

There is no need for p1j to be integers.

Maybe there is an easier way to achieve optimum homogenity of the three groups - I don't know.

If you use the problem-based approach as @Matt J suggests, I needed 16 lines of code to setup and solve the problem as formulated by the equalities and inequalities above. So it's worth to start reading the documentation.

[Matt J]

You can use MATLAB's "intlinprog" to solve.

But I wouldn't use intlinprog directly. You should just use the Problem-Based set-up tools to define the problem and let solve pick the optimizer for you.

[Mat]

This was solved using chatGPT.

[Torsten]

So problem solved ?

Can you give us an indication about the size of N ?

[Mat]

Apparently. N=72

[Torsten]

I think N = 72 will be too large for the way I suggested.

N = 6 cannot finish with MATLAB Online (i.e. within 235 seconds).

You should search for approximate (heuristic) algorithms. And you should search for the name of this problem so that you can find information on how it can be efficiently tackled (I don't think it can be classified as "bin packing problem" as @Steven Lord suggests).

But test it first on your PC:

rng("default")
N = 5;
p1 = rand(3*N,1);
prob = optimproblem;
X = optimvar('X',3*N,3*N,'Type','integer','LowerBound',0,'UpperBound',1);
d = optimvar('d',3,1,'LowerBound',0);
prob.Objective = sum(d);
y = X*p1;
prob.Constraints.c11 = sum(y(1:N))-sum(y(N+1:2*N)) <= d(1);
prob.Constraints.c12 = sum(y(1:N))-sum(y(N+1:2*N)) >= -d(1);
prob.Constraints.c21 = sum(y(1:N))-sum(y(2*N+1:3*N)) <= d(2);
prob.Constraints.c22 = sum(y(1:N))-sum(y(2*N+1:3*N)) >= -d(2);
prob.Constraints.c31 = sum(y(N+1:2*N))-sum(y(2*N+1:3*N)) <= d(3);
prob.Constraints.c32 = sum(y(N+1:2*N))-sum(y(2*N+1:3*N)) >= -d(3);
prob.Constraints.scr = sum(X,1) == ones(1,3*N);
prob.Constraints.scc = sum(X,2) == ones(3*N,1);
sol = solve(prob)

Solving problem using intlinprog.
Running HiGHS 1.6.0: Copyright (c) 2023 HiGHS under MIT licence terms
Presolving model
36 rows, 228 cols, 1356 nonzeros
36 rows, 228 cols, 1296 nonzeros

Solving MIP model with:
   36 rows
   228 cols (225 binary, 0 integer, 0 implied int., 3 continuous)
   1296 nonzeros

        Nodes      |    B&B Tree     |            Objective Bounds              |  Dynamic Constraints |       Work      
     Proc. InQueue |  Leaves   Expl. | BestBound       BestSol              Gap |   Cuts   InLp Confl. | LpIters     Time

         0       0         0   0.00%   0               inf                  inf        0      0      0         0     0.0s
 R       0       0         0   0.00%   0               4.56709472       100.00%        0      0      0        29     0.0s
 C       0       0         0   0.00%   0               1.582483596      100.00%      248      8      0        64     0.0s
 L       0       0         0   0.00%   0               0.2066885925     100.00%     1121     28      0       294     0.2s

Symmetry detection completed in 0.0s
Found 14 generators

 L       0       0         0   0.00%   0               0.0551451725     100.00%     1121      5      0       752     0.3s
 B      63      17        18   2.97%   0               0.0426948352     100.00%     1138      5    167      1731     0.3s
 B      75      17        24   3.01%   0               0.0426948352     100.00%     1145      5    241      1813     0.4s
 B     152      26        56  55.40%   0               0.040381773      100.00%     1233     12    519      2452     0.4s
 B     183      27        70  55.45%   0               0.040381773      100.00%     1238     12    602      2653     0.4s
 B     351      42       141  58.99%   0               0.00327539493    100.00%     1353     13   1510      4038     0.5s
 B     355      42       142  58.99%   0               0.00327539493    100.00%     1355     13   1519      4093     0.5s
 T     356      42       143  58.99%   0               0.00327539493    100.00%     1363     13   1559      4116     0.5s
 B    2589      71      1177  79.35%   0               0.00327539493    100.00%     1534     25   7638     28242     2.3s
 L    2747      76      1253  79.95%   0               0.00327539492    100.00%     1519     25   8090     29729     2.7s

Solving report
  Status            Optimal
  Primal bound      0.00327539491502
  Dual bound        0.00327539491502
  Gap               0% (tolerance: 0.01%)
  Solution status   feasible
                    0.00327539491502 (objective)
                    0 (bound viol.)
                    1.83012701771e-11 (int. viol.)
                    0 (row viol.)
  Timing            4.84 (total)
                    0.00 (presolve)
                    0.00 (postsolve)
  Nodes             5496
  LP iterations     59098 (total)
                    4633 (strong br.)
                    3347 (separation)
                    3021 (heuristics)

Optimal solution found.

Intlinprog stopped because the objective value is within a gap tolerance of the optimal value, options.AbsoluteGapTolerance = 1e-06. The intcon variables are integer within
tolerance, options.ConstraintTolerance = 1e-06.

sol = struct with fields:
    X: [15x15 double]
    d: [3x1 double]

Xsol = sol.X

Xsol = 15x15

0 0 0 0 1.0000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0000 0 0 0 0 -0.0000 0 0 0 0 1.0000 0 0 -0.0000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0000 0 0 0 0 0 0.0000 0 0 0 0 0 0 0 1.0000 0 0 0 0 1.0000 0 0 0.0000 0 0 0 0 0 0 -0.0000 0 0 0 0 0 0 0 1.0000 0 0 0 -0.0000 -0.0000 0 0 0 0 0 0 0 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.0000 0 0 0 0 0 0 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

dsol = sol.d

dsol = 3x1

0.0012 0.0004 0.0016

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

ysol = Xsol*p1;
sum(ysol(1:N))

ans = 3.2034

sum(ysol(N+1:2*N))

ans = 3.2022

sum(ysol(2*N+1:3*N))

ans = 3.2039

[row,~] = find(abs(Xsol.'-1)<1e-3);
sort(row(1:N))

ans = 5x1

5 9 11 12 14

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

p1(sort(row(1:N)))

ans = 5x1

0.6324 0.9575 0.1576 0.9706 0.4854

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

sort(row(N+1:2*N))

ans = 5x1

1 3 4 8 15

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

p1(sort(row(N+1:2*N)))

ans = 5x1

0.8147 0.1270 0.9134 0.5469 0.8003

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

sort(row(2*N+1:3*N))

ans = 5x1

2 6 7 10 13

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

p1(sort(row(2*N+1:3*N)))

ans = 5x1

0.9058 0.0975 0.2785 0.9649 0.9572

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[Torsten]

The name of the problem is "Balanced number partitioning".

If you choose n = 3*N, m = 3 and k = N, you can check heuristic methods to solve your problem here:

https://en.wikipedia.org/wiki/Balanced_number_partitioning

Especially have a look at this publication:

Babel, Luitpold; Kellerer, Hans; Kotov, Vladimir (1998-02-01). "The k-partitioning problem". Mathematical Methods of Operations Research. 47 (1): 59–82. doi:10.1007/BF01193837. ISSN 1432-5217. S2CID 5594197.