We have a number of audio tracks, and a number of groups called “sides” (each corresponding to one side of a vinyl record). For each track, we know its name and its length in seconds. Each side also has a length, which is the sum of the lengths of the tracks on the side. Typically one side will be longer than all the others. This the maximum side length, and it varies depending on how we arrange the tracks. What arrangement of the tracks results in the SMALLEST maximum side length? The problem reduces to Multiway number partitioning.
User @erwin-kalvelagen proposed using a MIP solver (CPLEX). However his GAMS code generates the input song lengths randomly, whereas I need to input specific song lengths. Also, his code outputs the computed side lengths, but not which song goes on which side. I managed to adapt his code to solve both of these problems. The corrected code is below, along with its NEOS output.
set
i 'songs' /song1*song23/
j 'sides' /side1*side8/
;
parameter len(i) 'length of song in seconds' /
song1 323,
song2 289,
song3 329,
song4 424,
song5 351,
song6 369,
song7 298,
song8 296,
song9 358,
song10 358,
song11 342,
song12 275,
song13 300,
song14 280,
song15 283,
song16 348,
song17 347,
song18 357,
song19 243,
song20 244,
song21 255,
song22 308,
song23 360
/;
binary variable x(i,j) 'assignment';
variable z 'maximum length of side';
equations
maximum(j) 'bound on z'
assign(i) 'assignment constraints'
;
maximum(j).. z =g= sum(i, len(i)*x(i,j));
assign(i).. sum(j, x(i,j)) =e= 1;
model m /all/;
option threads=0, mip=cplex;
solve m minimizing z using mip;
parameter slen(j) 'solution: length of sides';
slen(j) = sum(i, len(i)*x.l(i,j));
option decimals = 0;
option x:0:0:1;
display len,slen,z.l,x.l;
Resulting output from CPLEX at NEOS (abbreviated):
Proven optimal solution
MIP Solution: 936.000000 (160060 iterations, 15252 nodes)
Final Solve: 936.000000 (0 iterations)
Best possible: 936.000000
Absolute gap: -0.000000
Relative gap: -0.000000
[snip]
---- 56 PARAMETER slen solution: length of sides
side1 936, side2 936, side3 936, side4 793, side5 931, side6 936
side7 933, side8 936
---- 56 VARIABLE z.L = 936 maximum length of si
de
---- 56 VARIABLE x.L assignment
song1 .side1 1
song2 .side3 1
song3 .side7 1
song4 .side4 1
song5 .side8 1
song6 .side4 1
song7 .side2 1
song8 .side6 1
song9 .side2 1
song10.side1 1
song11.side8 1
song12.side5 1
song13.side3 1
song14.side2 1
song15.side6 1
song16.side5 1
song17.side3 1
song18.side6 1
song19.side8 1
song20.side7 1
song21.side1 1
song22.side5 1
song23.side7 1
Probably the easiest is to develop a small MIP (Mixed-Integer Programming) model. This can be solved with readily available MIP solvers.
The very simple mathematical model can look like:
Indices:
i : song
j : side
Data:
length[i]
Decision variables:
x[i,j]=1 if song i is assigned to side j
0 otherwise
z : maximum length
Optimization Model:
minimize z
subject to the constraints:
z >= sum(i, x[i,j]*length[i]) ∀j
sum(j, x[i,j]) = 1 ∀i
x[i,j] ∈ {0,1}
Of course, this must be transcribed into code to be fed to a solver. As the model is very simple, this is quite easy. The exact syntax depends on the chosen language, modeling tool and solver.
To give an idea about performance: a problem with n=200 songs (with random lengths) and m=10 sides took 9 seconds to solve to proven optimality (using a good MIP solver). I also see good performance with Google's OR-Tools CP-Sat solver.
Test data and results
---- 27 PARAMETER len random data: length of song
song1 200.000, song2 282.000, song3 246.000, song4 216.000, song5 215.000, song6 207.000
song7 222.000, song8 283.000, song9 188.000, song10 240.000, song11 300.000, song12 250.000
song13 299.000, song14 272.000, song15 195.000, song16 257.000, song17 199.000, song18 210.000
song19 260.000, song20 232.000, song21 223.000, song22 222.000, song23 195.000, song24 198.000
song25 251.000, song26 280.000, song27 207.000, song28 260.000, song29 273.000, song30 216.000
song31 193.000, song32 240.000, song33 199.000, song34 285.000, song35 212.000, song36 214.000
song37 251.000, song38 267.000, song39 256.000, song40 236.000, song41 230.000, song42 194.000
song43 218.000, song44 185.000, song45 220.000, song46 202.000, song47 258.000, song48 247.000
song49 273.000, song50 216.000, song51 259.000, song52 271.000, song53 255.000, song54 214.000
song55 190.000, song56 192.000, song57 257.000, song58 245.000, song59 183.000, song60 275.000
song61 188.000, song62 201.000, song63 243.000, song64 270.000, song65 201.000, song66 184.000
song67 250.000, song68 255.000, song69 227.000, song70 223.000, song71 209.000, song72 209.000
song73 195.000, song74 292.000, song75 225.000, song76 274.000, song77 216.000, song78 195.000
song79 270.000, song80 188.000, song81 204.000, song82 180.000, song83 212.000, song84 240.000
song85 198.000, song86 201.000, song87 220.000, song88 218.000, song89 218.000, song90 296.000
song91 300.000, song92 224.000, song93 225.000, song94 273.000, song95 227.000, song96 290.000
song97 194.000, song98 268.000, song99 186.000, song100 249.000, song101 186.000, song102 180.000
song103 228.000, song104 242.000, song105 256.000, song106 207.000, song107 227.000, song108 213.000
song109 198.000, song110 293.000, song111 231.000, song112 196.000, song113 226.000, song114 225.000
song115 212.000, song116 294.000, song117 202.000, song118 215.000, song119 189.000, song120 228.000
song121 192.000, song122 226.000, song123 219.000, song124 203.000, song125 193.000, song126 252.000
song127 241.000, song128 185.000, song129 274.000, song130 294.000, song131 252.000, song132 253.000
song133 223.000, song134 251.000, song135 262.000, song136 241.000, song137 199.000, song138 259.000
song139 243.000, song140 195.000, song141 299.000, song142 207.000, song143 261.000, song144 273.000
song145 292.000, song146 204.000, song147 215.000, song148 203.000, song149 209.000, song150 258.000
song151 268.000, song152 190.000, song153 198.000, song154 232.000, song155 202.000, song156 263.000
song157 272.000, song158 198.000, song159 227.000, song160 264.000, song161 282.000, song162 254.000
song163 298.000, song164 183.000, song165 202.000, song166 190.000, song167 245.000, song168 195.000
song169 268.000, song170 193.000, song171 239.000, song172 276.000, song173 239.000, song174 244.000
song175 181.000, song176 245.000, song177 234.000, song178 298.000, song179 202.000, song180 199.000
song181 182.000, song182 201.000, song183 187.000, song184 182.000, song185 281.000, song186 252.000
song187 183.000, song188 203.000, song189 295.000, song190 220.000, song191 251.000, song192 211.000
song193 257.000, song194 198.000, song195 235.000, song196 227.000, song197 277.000, song198 245.000
song199 227.000, song200 247.000
---- 27 PARAMETER slen solution: length of sides
side1 4625.000, side2 4624.000, side3 4625.000, side4 4624.000, side5 4620.000, side6 4625.000
side7 4624.000, side8 4624.000, side9 4625.000, side10 4625.000
---- 27 VARIABLE z.L = 4625.000 maximum length of side
This was solved with GAMS and Cplex (commercial). The implementation looks like:
set
i 'songs' /song1*song200/
j 'sides' /side1*side10/
;
parameter len(i) 'random data: length of song';
len(i) = uniformint(3*60,5*60);
option len:0;
display len;
binary variable x(i,j) 'assignment';
variable z 'maximum length of side';
equations
maximum(j) 'bound on z'
assign(i) 'assignment constraints'
;
maximum(j).. z =g= sum(i, len(i)*x(i,j));
assign(i).. sum(j, x(i,j)) =e= 1;
model m /all/;
option threads=0, mip=cplex;
solve m minimizing z using mip;
parameter slen(j) 'solution: length of sides';
slen(j) = sum(i, len(i)*x.l(i,j));
display len,slen,z.l;
abort$(sum(i,len(i))<>sum(j,slen(j))) "solution does not add up";
You can submit this model to https://neos-server.org/neos/index.html if you don't have access to these tools.
Here is another implementation using ORTools CP-Sat solver (open source). Again, I just transcribed the mathematical model, which is really a trivial task (both for GAMS and for OR-tools). For the CP model syntax see the guides here: https://developers.google.com/optimization/cp.
from ortools.sat.python import cp_model
'''
Optimization Model
------------------
min z
z >= sum(i, x[i,j]*length[i]) ∀j
sum(j, x[i,j]) = 1 ∀i
x[i,j] ∈ {0,1}
For CP model the assumptions are:
All variables and all coefficients are integers
Note: CP model has also max() function. Here we just use the MIP
formulation.
Documentation on ORTools is here: https://developers.google.com/optimization
'''
lengths = [
200, 282, 246, 216, 215, 207, 222, 283, 188, 240, 300, 250, 299, 272, 195, 257,
199, 210, 260, 232, 223, 222, 195, 198, 251, 280, 207, 260, 273, 216, 193, 240,
199, 285, 212, 214, 251, 267, 256, 236, 230, 194, 218, 185, 220, 202, 258, 247,
273, 216, 259, 271, 255, 214, 190, 192, 257, 245, 183, 275, 188, 201, 243, 270,
201, 184, 250, 255, 227, 223, 209, 209, 195, 292, 225, 274, 216, 195, 270, 188,
204, 180, 212, 240, 198, 201, 220, 218, 218, 296, 300, 224, 225, 273, 227, 290,
194, 268, 186, 249, 186, 180, 228, 242, 256, 207, 227, 213, 198, 293, 231, 196,
226, 225, 212, 294, 202, 215, 189, 228, 192, 226, 219, 203, 193, 252, 241, 185,
274, 294, 252, 253, 223, 251, 262, 241, 199, 259, 243, 195, 299, 207, 261, 273,
292, 204, 215, 203, 209, 258, 268, 190, 198, 232, 202, 263, 272, 198, 227, 264,
282, 254, 298, 183, 202, 190, 245, 195, 268, 193, 239, 276, 239, 244, 181, 245,
234, 298, 202, 199, 182, 201, 187, 182, 281, 252, 183, 203, 295, 220, 251, 211,
257, 198, 235, 227, 277, 245, 227, 247
]
N = len(lengths) # number of songs
M = 10 # number of sides
I = range(N) # songs
J = range(M) # sides
#
# Specify Optimization Model
# This is a direct transcription of the mathematical model.
#
model = cp_model.CpModel()
x = {(i,j):model.NewBoolVar(f'x[{i},{j}]') for i in I for j in J}
z = model.NewIntVar(0,sum(lengths),'z')
for j in J:
model.Add(z >= sum([x[(i,j)]*lengths[i] for i in I]))
for i in I:
model.Add(sum([x[(i,j)] for j in J]) == 1)
model.Minimize(z)
#
# Solve
# for better performance on large instances use multiple worker threads
#
solver = cp_model.CpSolver()
status = solver.Solve(model)
print(f'Status: {solver.StatusName(status)}')
print(f'Objective value: {solver.ObjectiveValue()}')
#
# reporting
#
for j in J:
print(f'side: {j}, length: {sum([solver.Value(x[(i,j)])*lengths[i] for i in I])}')
The results should be:
Status: OPTIMAL
Objective value: 4625.0
side: 0, length: 4625
side: 1, length: 4621
side: 2, length: 4624
side: 3, length: 4622
side: 4, length: 4625
side: 5, length: 4624
side: 6, length: 4625
side: 7, length: 4625
side: 8, length: 4625
side: 9, length: 4625
Are any improvements possible? I would experiment with some additional symmetry-breaking constraints. They can make a lot of difference in performance.