Restrictions problem

Question

I tried to program a simple linear optimization model and don't know where my mistake lies. The function should be a sum over "t" and "z" from the product of two constants ("e" and "H") and the decision variable "X". The first constraint should say that every decision variable "X" (for all "t" and "z") can only take the value 0 or 1. The second constraint should say that every decision variable "X" (for all "t") can only select one option z. Therefore the sum of all "X" over "z" (for all "t") has to be equaling 1. The third constraint should say the same as the second only the other way around.

NMinimize[{Sum[
   Subscript[e, t]*Subscript[H, z]*Subscript[X, t, z] , {t, 1, 
    34}, {z, 1, 34}u], 
  ForAll[{t, z}, Subscript[X, t, z] \[Element] Integers, 
    0 <= Subscript[X, t, z] <= 1] && 
   ForAll[t, Sum[Subscript[X, t, z], {z, 1, 34}] == 1, {t, 1, 34} ] &&
    ForAll[z, 
    Sum[Subscript[X, t, z], {t, 1, 34}] == 1, {z, 1, 
     34}]}, {Subscript[X, t, z]}]

Answer 1

It is dangerous to try to describe a problem using non-working code, but if I have understood it correctly, the specific problem is to minimize $\mathbb{e' X h}$ given fixed $n$-vectors $\mathbb{e}$ and $\mathbb{h}$ subject to the constraints that the elements of $\mathbb{X}$ be in $\{0,1\}$ and that all the rows and all the columns of $\mathbb{X}$ sum to unity.

Perhaps we can take the general problem to be that of asking Mathematica to solve any such binary integer linear program where a linear function of 0-1 variables is to be optimized subject to linear constraints (including both equalities and inequalities).

Let's first solve the specific problem. The constraints on $\mathbb{X}$ compel each row and each column to contain exactly one $1$. Thus, the solution establishes a pairing between $\mathbb{e}$ and $\mathbb{h}$ to minimize the sum of products of pairs. Simply by renaming the indexes, we may assume both $\mathbb{e}$ and $\mathbb{h}$ are sorted in ascending order, in which case the obvious solution (which is unique when all elements in each vector are distinct) matches the elements in reverse order: the largest in $\mathbb{e}$ with the smallest in $\mathbb{x}$, etc.

Now that we know what the solution is, we can use this to check how well Mathematica does. All we have to do is set up the problem in a standard format. This shows how to solve the general problem.

1. The linear constraints (row sums and column sums) must each be given in the form $\mathbb{m}_i \cdot \mathbb{x} = b_i$. The vectors $\mathbb{m}_i$ are assembled into a matrix $\mathbb{M}$ and the constants on the right hand sides are collected correspondingly into a vector $\mathbb{b}$. The constraints are equivalent to the matrix equation $\mathbb{M x = b}$. One way to accomplish this is to work out first which indexes in $\mathbb{x}$ are involved in the row and column sums, then use these to construct $\mathbb{M}$ as a SparseArray object. I am going to do this in a way that may be inefficient and the opposite of terse, to make it clear what is being done and to make it easy to check its correctness. As a check, we can see the regular patterns formed by the row and column sums in a plot of $\mathbb{M}$.

   n = 34;
   rowSums = Range[#, # + (n - 1) n, n] & /@ Range[n];
   colSums = Range[n (# - 1) + 1, n # ] & /@ Range[n];
   sums = rowSums~Join~colSums;
   rules = Flatten[Table[{i, #} -> 1 & /@ sums[[i]], {i, 1, Length[sums]}]];
   m = SparseArray[rules];
   b = ConstantArray[{1, 0}, 2 n]; (* There are 2n sum-to-unity constraints *)
   ArrayPlot[m]
   

   

   (The elements of b are pairs {1,0} where the first entry is $b_i$ and the second entry, $0$, is a code indicating that the constraint is an equality rather than an inequality. See the help page for more information.)

2. Although it is not necessary in this case, it helps to stipulate that the coefficients of $\mathbb{x}$ must each lie between $0$ and $1$ and that they must be integers: the solution will be found much faster. These interval constraints (inequalities) are specified by an array of their endpoint values ({0,1}), one per variable:

   lu = ConstantArray[{0, 1}, n^2];
   

3. To see how this all works, generate some coefficients for $\mathbb{e}$ and $\mathbb{h}$:

   e = Sort[RandomReal[{-1, 1}, n]]; h = Sort[RandomReal[{-1, 1}, n]];
   {e}\[Transpose] . {h} // ReliefPlot
   

   

   The objective function, $\mathbb{e X h}$, must be written in the form $\mathbb{c} \cdot \mathbb{x}$ where $\mathbb{x}$ is a vector consisting of the entries in $\mathbb{X}$. Assuming $\mathbb{X}$ will be unraveled in column-first order,

   c = Flatten[{e}\[Transpose] . {h}];
   

Computing a solution is anticlimactic and fast:

AbsoluteTiming[x = LinearProgramming[c, m, b, lu, Integers];]

{0.0156000, Null}

A natural way to inspect the solution is to ravel $\mathbb{x}$ back into the matrix form $\mathbb{X}$, column by column:

ArrayPlot[Partition[x, n]\[Transpose]]

And there it is: the elements of $\mathbb{e}$ (representing rows) are matched in reverse order with the elements of $\mathbb{h}$ (representing columns).