Un Peu de Math avec CPLEX
There is one thing I didn't disclose in my previous post (Un Peu de Math) on the Anal Let me start with a statement of the optimization problem to be solved. A permutation x is a function such that x(i) is in the set {0,1,...,n1} for all i in that set. Furthermore, x(i) is different from x(j) whenever i is different from j. Given a permutation x we define three functions of x, where a stands for the absolute value of a: u(x) = sum_{i} x(i)  i v(x) = sum_{i} (x(i)  (n1 i) t(x) = min( u(x), v(x) ). The problem is to find permutations that maximize t, for various values of n. How did I model this problem as a mixed integer programming (MIP) problem? Issue was that absolute value isn't a linear function. I could have tried some of the methods described in Paul Rubin's excellent post, but I decided to use a different approach. A classical way to linearize functions of integer variables is to replace the integer decision variables by a set of binary variables, one per possible value of the original integer variables. I therefore used one binary variable x_{i,j} equal to 1 if x(i) = j. Then the objective function was straightforward to express. The resulting MIP model was maximize t subject to t <= sum_{i,j} j  i x_{i,j} t <= sum_{i,j} j  (ni1) x_{i,j} sum_{i} x_{i,j} = 1 for all j sum_{j} x_{i,j} = 1 for all i x_{i,j} binary The first two constraints define the objective function, and the last two express that x is a permutation. This model worked well for small value of n but running times grew quickly with n. I needed to do better. A mathematical analysis (see Un Peu de Math for details) of the problem gave me a nice upper bound for t: t(x) <= 6m^{2} + 3mr + r(r1)/2 where m is the quotient and let r the remainder of the Euclidean division of n by 4: n = 4m + r, 0 <= r < 4, Adding this bound as a constraint on t in the model didn't help. as running times increased. This was probably due to the fact that the extra constraints makes the dual degenerate, as explained in this nice post. I therefore discarded the constraint. A side effect of the computation of the above bound was that the permutations that are likely to maximize t have a specific structure. If we map the variables x_{i,j} on the nxn square, then we are looking for permutations that set to 1 variables that are as close as possible to the sides of the square. More formally, let's define concentric rings as in Fig. 1. Then we are looking for permutations that set 4 variables to 1 in each of the first m rings (see Un Peu de Math for details)
Figure 1. Concentric rings.
I could have set a constraint per ring to enforce that but it didn't seem appealing for two reasons: (a) writing them requires convoluted index manipulation (yes, I'm lazy), and (b) it adds many constraints. Moreover, CPLEX was likely to select these variables anyway given these are the ones that are likely to maximize the objective. I rather thought about consequences of the structure of the permutations we are looking for. Let's consider a variable x_{i,j }in the lower left corner, i.e. with i < m and j < m, as in Fig 2. If we set it to 1, then no other variable can be set to 1 in the ith column or in the jth row. If i=j it means that there will be at most 3 variables set in ring R_{i}. because it is both in the ith row and the ith column of R_{i}. If i < j it means that no variable can be set to 1 in the lower side of R_{j}. therefore only 3 variables could be set to 1 in it. For these reasons, x_{i,j} must be set to 0. Figure 2. The black square represents x_{i,j} = 1 Similarly, we can sett variables to 0 in the other three corners of the square. As a result we set to 0 all the variables x_{i,j} such that: (0 <= i < m or n  m <= i < n) and (0 <= j < m or n  m <= j < n), as depicted in Fig. 3. Figure 3. Setting some variables to 0 Adding these constraints helped a lot, and the MIP model were solved at the root node for many values of n. However, for many values of n as well, the models still seemed hard to solve. Looking at the hard ones I soon discovered that these were the models for which the bound 6m^{2} + 3mr + r(r1)/2 was odd. It was quite easy to obtain permutations that were only 1 apart from the bound, but there seemed to be none for which the bound could be reached. I then resumed my mathematical analysis and was able to prove that t was always even (see Un Peu de Math for details). I then added a constraint: to make t even: t = 2 y y integer With these, all models seemed easy to solve. At least they were solved at the root node for all the values I tried (up to n=1000). The final model is given below. maximize t subject to t <= sum_{i,j} j  i x_{i,j} t <= sum_{i,j} j  (ni1) x_{i,j} t = 2 y sum_{i} x_{i,j} = 1 for all j sum_{j} x_{i,j} = 1 for all i x_{i,j} = 0 for all i,j s.t. (0 <= i < m or n  m <= i < n) and (0 <= j < m or n  m <= j < n) y integer x_{i,j} binary This comforted me in thinking that it was always possible to reach the bound if it was even, or to be 1 apart from it when it was odd. Looking at the structure of the permutations found by my model, I then looked for generic permutations by selecting m cells in each of the areas A,B,C,D of the nxn square, and r cells in the ring R of the nxn square. These regions are depicted on Fig. 4. Figure 4. Areas used for constructing optimal permutations A is the set of cells (i, x(i)) with 0 <= i < m, m <= x(i) < nm B is the set of cells (i, x(i)) with m <= i < n m, nm <= x(i) < n C is the set of cells (i, x(i)) with nm <= i < n, m <= x(i) < nm D is the set of cells (i, x(i)) with m <= i < nm,0 <= x(i) < m R is the ring R_{m} From it I was able to derive a generic form displayed in figures 3 to 6 in my previous post, which concluded my proof. There are two lessons here.
First, I don't think I could have solved that math Second, one way to improve MIP models is to better understand the mathematical structure or the optimization problem to be solved. It is not necessary to establish formal proofs as I did in this case, but mathematical insights can usually be translated into effective model modifications.
