An assignment problem with a wrinkle

In [1], the following problem is posed:

Consider two arrays \(\color{darkblue}a_i\) (length \(\color{darkblue}m\)) and \(\color{darkblue}b_j\) (length \(\color{darkblue}n\)) with \(\color{darkblue}m \lt \color{darkblue}n\). Assign all values \(\color{darkblue}a_i\) to a \(\color{darkblue}b_j\) such that:

• Each \(\color{darkblue}b_j\) can have 0 or 1 \(\color{darkblue}a_i\) assigned to it.
• The assignments need to maintain the original order of \(\color{darkblue}a_i\). I.e. if \(\color{darkblue}a_i \rightarrow \color{darkblue}b_j\) then \(\color{darkblue}a_{i+1}\) must be assigned to a slot in \(\color{darkblue}b\) that is beyond slot \(j\). In the picture below that means that arrows cannot cross.
• Do this while minimizing the sum of the products.

Mixed-integer programming model

This problem can be viewed as an assignment problem with a side constraint:

MIP Model
\[\begin{align}\min& \sum_{i,j}\color{darkred}x_{i,j}\cdot\color{darkblue}a_i\cdot\color{darkblue}b_j \\ &\sum_j \color{darkred}x_{i,j}=1 &&\forall i\\ & \sum_i \color{darkred}x_{i,j}\le 1 &&\forall j\\ & \color{darkred}v_i = \sum_j j \cdot \color{darkred} x_{i,j}\\ & \color{darkred}v_i \ge \color{darkred}v_{i-1}+1\\ & \color{darkred}x_{i,j} \in \{0,1\} \\ & \color{darkred}v_i \ge 1\end{align}\]

Here \(\color{darkred}x_{i,j}\) indicates the assignment \(\color{darkblue}a_i \rightarrow \color{darkblue}b_j\). The variable \(\color{darkred}v_i\) represents the position in \(\color{darkblue}b\) to which \(\color{darkblue}a_i\) is assigned.

The output of this model can look like: 

----     40 SET i  

i1,    i2,    i3


----     40 SET j  

j1,    j2,    j3,    j4,    j5,    j6


----     40 PARAMETER a  

i1 1.000,    i2 2.000,    i3 3.000


----     40 PARAMETER b  

j1  4.000,    j2  9.000,    j3  5.000,    j4  3.000,    j5  2.000,    j6 10.000


----     40 VARIABLE z.L                   =       16.000  objective

----     40 VARIABLE x.L  assignment

            j1          j4          j5

i1       1.000
i2                   1.000
i3                               1.000


----     40 VARIABLE v.L  position of a(i) in b

i1 1.000,    i2 4.000,    i3 5.000

Notes:

• We can tighten the bounds on variable \(\color{darkred}v_i\) to \(\color{darkred}v_i \in [i,\color{darkblue}n-\color{darkblue}m+i]\).
• The variable \(\color{darkred}v_i\) is automatically integer-valued. We can declare it as an integer or continuous variable.
• If declared as a continuous variable, Cplex may make it an integer variable. We can see this in the log: Reduced MIP has 17737 binaries, 99 generals, 0 SOSs, and 0 indicators.
• A larger problem with \(\color{darkblue}m=50, \color{darkblue}n=500\) solves in about 50 seconds.

It is interesting to see what happens when we feed this model an even larger data set. I generated some random data: \[\begin{align} & \color{darkblue}m=100\\ &\color{darkblue}n=1000 \\ & \color{darkblue}a_i \sim U(0,100) \\ & \color{darkblue}b_j \sim U(0,100)\end{align}\] This is not a small problem: we will have 100,000 binary variables. This model has no longer a network structure, so compared to a pure assignment problem we see much poorer performance:

	Assignment problem as LP	Assignment problem as MIP	Full MIP model
Rows/Columns	1,100/100,000	1,100/100,000	1,300/100,100
Objective	11,524.178	11,524.178	14,371.455
Time (seconds)	1	7	1,853
Nodes	-	0	3,419
Iterations	2,672	2,727	30,7429

The assignment problem was formed by dropping the constraints involving \(\color{darkred}v_i\). Further, the problem was solved as an LP. It is noted that Cplex behaves differently when we use continuous variables \(\color{darkred}x_{i,j}\) or declare them as binary variables, When they are binary variables, the presolver works differently, and as a result, the solution of the relaxed problem is no longer optimal: it has to do some work to find the optimal integer solution. In this case, this was done in MIP preprocessing: no branching was needed. When solved as a pure LP, the optimal solution is integer-valued automatically. There is a  difference in solution time: the LP solves in 1 second, while the MIP version takes 7 seconds. 

Let's see if we can shave something off this 1,800 seconds solution time. The first experiment I did was to provide an initial solution. I took the 100 smallest values of \(\color{darkblue}b_j\) and assigned the 100 \(\color{darkblue}a_i\)'s to these while maintaining the order. This is very simple, so this took no time. The MIP was solved with a MIPSTART option. The results are:

 

	Assignment problem as LP	Assignment problem as MIP	Full MIP model	Using initial solution
Rows/Columns	1,100/100,000	1,100/100,000	1,300/100,100	1,300/100,100
Objective	11,524.178	11,524.178	14,371.455	initial: 17,739.702 optimal: 14,371.455
Time (seconds)	1	7	1,853	1,458
Nodes	-	0	3,419	3,546
Iterations	2,672	2,727	307,429	343,969

This model solves faster. We don't see this in the node and iteration count, mainly because Cplex spends much time in preprocessing and restarts. Our initial objective is not that great. But this leads to the question, can we select the say 200 smallest values of \(\color{darkblue}b_j\) and solve this smaller problem first. Then do a second solve with use all values of \(\color{darkblue}b_j\), again using a MIPSTART. So the algorithm becomes:

1. Select 100 smallest values of \(\color{darkblue}b_j\). Assign \(\color{darkblue}a_i\) in order. This gives an initial solution of 17,739.702. 
2. Pick 200 of of the smallest values of \(\color{darkblue}b_j\). Solve our MIP problem on this subset. This gives an optimal solution of 14,451.861 in  118 seconds.
3. Use this as an initial solution for the whole problem. This now solves in 972 seconds. Total solution time is 118+972=1,090 seconds.

Adding this to our table gives:

	Assignment problem as LP	Assignment problem as MIP	Full MIP model	Using initial solution	Using two models
Rows/Columns	1,100/100,000	1,100/100,000	1,300/100,100	1,300/100,100	subset: 500/20,100 full: 1,300/100,100
Objective	11,524.178	11,524.178	14,371.455	initial: 17,739.702 optimal: 14,371.455	initial: 17,739.702 subset: 14,451.861 full: 14,371.455
Time (seconds)	1	7	1,853	1,458	subset: 118 full: 972
Nodes	-	0	3,419	3,546	subset: 3,205 full: 3,500
Iterations	2,672	2,727	307,429	343,969	subset: 276,796 full: 456,537

Conclusions

• Sometimes it makes sense to replace one big solve with a number of solves. Each one will start from and improve a previously found solution. 
• In production environments, a setup like this has another advantage: if a single model fails, you still get good solutions.
• In the above table, the statistics for the number of nodes and Simplex iterations are not very indicative of the total effort. As MIP solvers become more complex, the best performance measure is solution time.
• An assignment problem (or other network problems) should be solved as an LP, not as a MIP.

References

1. Minimize sum of product of two uneven consecutive arrays, https://stackoverflow.com/questions/65884107/minimize-sum-of-product-of-two-uneven-consecutive-arrays
2. Assignment problem with a wrinkle formulated as a network problem, https://yetanothermathprogrammingconsultant.blogspot.com/2021/01/assignment-problem-with-wrinkle.html. Write-up on how to formulate this problem as a shortest path problem.

Appendix: GAMS model

This is the model that I used to try out the three-tier approach.

$ontext     Solve ordered assignment problem in three stages:     1. Find a trivial initial solution (100 smallest values in b)     2. Find a better solution using 200 smallest values     3. Solve the full problem $offtext *------------------------------------------------------- * data *------------------------------------------------------- sets   i /i1*i100/   j /j1*j1000/ ; parameter   A(i)   B(j) ; a(i) = uniform(0,100); b(j) = uniform(0,100); display a,b; *------------------------------------------------------- * optimization model *------------------------------------------------------- set subset(j) 'subset of b'; binary variable    x(i,j) 'assign' ; variable    v(i) 'position of a(i) in b'    z    'objective' ; * tight bounds on v v.lo(i) = ord(i); v.up(i) = card(j) - (card(i)-ord(i)); equations     assign1(i) 'assignment constraint'     assign2(j) 'assignment constraint'     calcv      'caculate v'     order      'ordering constraint'     obj        'objective' ; obj..     z =e= sum((i,subset(j)),x(i,j)*a(i)*b(j)); assign1(i)..     sum(subset(j), x(i,j)) =e= 1; assign2(subset(j))..     sum(i, x(i,j)) =l= 1; calcv(i)..     v(i) =e= sum(subset(j), ord(j)*x(i,j)); order(i-1)..     v(i) =g= v(i-1)+1; model m /all/; option optcr=0, threads=8; *------------------------------------------------------- * set initial values using trivial solution of * 100 smallest values in b. *------------------------------------------------------- * find 100 smallest values in b subset(j) = no; scalars    bmin   'minimum value in remaining b'    k      'loop variable' ; for(k=1 to 100,    bmin = smin(j$(not subset(j)), b(j));    loop(j$(b(j)=bmin),       subset(j) = yes;       break;    ); ); * let x reflect this initial subset loop(i,    loop(j$subset(j),       x.l(i,j) = 1;       subset(j) = 0;       break;    ); ); z.l = sum((i,j),x.l(i,j)*a(i)*b(j)); v.l(i) = sum(j, ord(j)*x.l(i,j)); option v:0; display "***********  initial",z.l,v.l; *------------------------------------------------------- * limit search to 200 values *------------------------------------------------------- for(k=1 to 200,    bmin = smin(j$(not subset(j)), b(j));    loop(j$(b(j)=bmin),       subset(j) = yes;       break;    ); ); m.optfile=1; solve m minimizing z using mip; display "***********  solve 1",z.l,v.l; *------------------------------------------------------- * solve full model *------------------------------------------------------- subset(j) = yes; solve m minimizing z using mip; display "***********  solve 2",z.l,v.l; $onecho > cplex.opt mipstart 1 $offecho