Clear Filters
Clear Filters

Unconstrained minimisation problem with a complicated range

27 views (last 30 days)
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?
  7 Comments

Sign in to comment.

Answers (1)

Torsten
Torsten on 10 Sep 2024
Edited: Torsten on 10 Sep 2024
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.
  8 Comments
Torsten
Torsten on 11 Sep 2024 at 17:03
Edited: Torsten on 11 Sep 2024 at 18:19
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
Torsten on 14 Sep 2024 at 14:45
Edited: Torsten on 14 Sep 2024 at 15:03
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:
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.

Sign in to comment.

Tags

Products


Release

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!