A "Decoupled" Quadratic Program

Someone posted a question on OR Stack Exchange about a problem with a "nonseparable bilinear objective function" with "decoupled constraints". The problem has the following form:

\begin{alignat*}{1} \min & \sum_{i=1}^{m}\sum_{j=1}^{n}a_{ij}x_{i}y_{j}+b^{\prime}x+c^{\prime}y\\ \textrm{s.t. } & x\in X\subset\mathbb{R}^{m}\\ & y\in Y\subset\mathbb{R}^{n} \end{alignat*}

where $X$ and $Y$ are polytopes (i.e., defined by linear constraints, and bounded) in the positive orthants of their respective spaces (i.e., $x\ge 0$ and $y\ge 0$) and $a$, $b$ and $c$ are all nonnegative. So besides everything being nonnegative, the key things here are that the objective contains bilinear terms (always involving the product of an $x$ and a $y$) and each constraint involves either only $x$ or only $y$. The author was looking for a way to exploit the separability of the constraints in solving the problem.

It occurred to me that a heuristic approach might be to solve a sequence of linear programs, alternating between $X$ and $Y$. Let $h(x,y)$ be the objective function of the original problem, and let $$f(x;y) = \sum_i\sum_j a_{ij}x_i y_j + b^\prime x$$and $$g(y;x)=\sum_i\sum_j a_{ij}x_i y_j + c^\prime y$$ be respectively the portions of the objective dependent on $x$ with $y$ fixed or dependent on $y$ with $x$ fixed. Each is linear in one set of variables (the other set being fixed). So we could proceed as follows:

1. Minimize $f(x;0)$ over $X$ (a linear program) to get a starting value $x^0\in X.$
2. Minimize $g(y;x^0)$ over $Y$ to get a starting value $y^0\in Y.$ We now have an incumbent solution with objective value $h(x^0,y^0)=b^\prime x^0 + g(y^0;x^0).$ Set $t=0$.
3. Repeat the following until the first time that the incumbent does not improve.
1. Minimize $f(x;y^t)$ over $X$ to get $x^{t+1}.$ If $f(x^{t+1},y^t) + c^\prime y^t$ is less than the incumbent value, make $(x^{t+1}, y^t)$ the new incumbent, else punt.
2. Minimize $g(y;x^{t+1})$ over $Y$ to get $y^{t+1}.$ If $g(y^{t+1};x^{t+1}) + b^\prime x^{t+1}$ is less than the incumbent value, make $(x^{t+1},y^{t+1})$ the new incumbent and increment $t$, else punt.

This will generate a sequence of solutions, each a corner point of the feasible region $X\times Y,$ with monotonically decreasing objective values.

Before continuing, it is worth noting that we are guaranteed that at least one corner of $X\times Y$ will be an optimum. This is of course always true when the objective is linear, but not always true with a quadratic objective. In our case, suppose that $(x^*, y^*)$ is an arbitrary optimum for the original problem. Then $x^*$ must minimize $f(x;y^*)$ over $X$, so either $x^*$ is a corner of $X$ of (if there are multiple optimal) $x^*$ can be replaced by a corner of $X$ with the same value of $f(x;y^*)$ (and thus the same value of $h(x;y^*)$, since the difference $c^\prime y^*$ does not depend on $x$). Apply the same logic on the other variable to show that $y*$ is either a corner of $Y$ or can be replaced by one.

I'm still calling this a heuristic, because there is one more piece to the puzzle. Could the procedure stop at a corner of $X\times Y$ that is a local but not global optimum? I'm not sure. Offhand, I do not see a way to prove that it will not, but I also cannot easily conjure up an example where it would.

To test this (and to confirm that I was not hallucinating, which has been known to happen), I wrote some Java code to generate random test problems and try the procedure. The code uses CPLEX to solve the LPs. In all test cases, it terminated with a solution (at least locally optimal) in a pretty small number of iterations. 

As a benchmark, I tried solving the full QP models with CPLEX, setting its "optimality target" parameter to "OPTIMALGLOBAL", which tells CPLEX to search for a global optimum to a nonconvex problem. This causes CPLEX to turn the problem into a mixed-integer quadratic program, which shockingly takes longer to solve than an LP. In a limited number of test runs, I observed a surprisingly consistent behavior. At the root node, CPLEX immediately found a feasible solution and then immediately found a better solution that matched what my procedure produced. After than (and within a usually stingy time limit I set), CPLEX made progress on the bound but never improved on the incumbent. In two test runs, CPLEX actually reached proven optimality with that second incumbent, meaning my procedure had found a global optimum.

So perhaps my "heuristic" can actually be shown to be an exact algorithm ... or perhaps not. At least in my test runs, CPLEX working on the QP found the best solution about as fast as, or maybe slightly faster than, my procedure did. On the other hand, my procedure only requires an LP solver, so it will work with code that does not accept QPs.

My Java code (which requires CPLEX) is available here.

Addendum: User "Dusto" on the Discord Operations Research server posted a simple counterexample to global convergence. The example has no constraints other than bounds on the variables (from 0 to 10 in all cases). The $b$ and $c$ vectors are strictly positive, so the linear programs in my heuristic will start at the origin and get stuck there. The $A$ matrix is crafted so that a negative overall objective value is attainable.

Tournament Scheduling: CPLEX v. CP Optimizer (Part 3)

This is the third and final chapter in my opus about scheduling "pods" (players) in a tournament. Reading the first post will be essential to understanding what is going on here. Reading the second post (in which I formulate an integer programming model for the problem) is definitely optional.

I will stick with the notation introduced in the prior post. As before, there are $N$ pods (players) being organized into teams of two to play in $G$ games (two teams per game, one in white, one in black). Each game has two slots corresponding to the two teams in the game.

• $t$ indexes a team;
• $p$ and $q$ index pods;
• $g$ indexes a game;
• $s$ indexes a slot (and $s'$ is the other slot of the same game);
• $c$ is a jersey color (0 for white, 1 for black).

For a slot $s$, I will use $G_s$ and $C_s$ to denote respectively the game in which the slot exists and the jersey color (0 or 1) of the team in the slot. For a team $t$, $P_{t,0}$ and $P_{t,1}$ will be the two pods on the team. (The order they are listed is immaterial.)

The model variables are almost (but not exactly) what they were in the IP model.

• $x_{t}\in \lbrace 0, \dots, 2G-1\rbrace$ is the index of the slot in which team $t$ plays (where slots 0 and $G$ are the white and black teams in the first game and slots $G-1$ and $2G-1$ are the white and black teams in the final game).
  
• $y_{pqg}\in \lbrace 0,1\rbrace$ is 1 if pods $p$ and $q\neq p$ are playing in game $g$ on opposing teams.
  
• $z_{pg}\in \lbrace 0,1\rbrace$ is the color jersey worn by pod $p$ at the conclusion of game $g$ (regardless of whether or not $p$ played in game $g$).
• $w_{pg}\in \lbrace 0,1\rbrace$ is 1 if pod $p$ changes jersey color going into game $g.$
• $v_{pg} \in \lbrace 0,1\rbrace$ is 1 if pod $p$ is playing in game $g.$

The changes from the IP model are that $x$ is general integer rather than binary (with dimension 1 rather than 2), the third index of $y$ is the game rather than the slot, and $v$ is a new variable. As in the IP model, the objective is to minimize $\sum_p \sum_g w_{pg}.$

The constraints are as follows. Note that the requirement that every team play exactly once is captured by the assignment of a single slot value to each team (via $x$).

• Two different teams cannot occupy the same slot: $$\textrm{allDiff}(x).$$
• Occupying a slot implies playing in the game: $$x_t = s \implies (v_{P_{t,1},G_s} = 1 \wedge v_{P_{t,2},G_s} = 1) \quad \forall t,s.$$
• Occupying a slot determines the color of the jersey worn in the game: $$x_t = s \implies  (z_{P_{t,1},G_s} = C_s \wedge z_{P_{t,2},G_s} = C_s) \quad \forall t,s.$$
• Jersey color for a pod remains constant from game to game unless a change is recorded: $$ z_{p,g} \neq z_{p,g-1}  \iff w_{p,g} = 1 \quad \forall p, \forall g > 0.$$
• Playing in consecutive games precludes changing jerseys: $$(v_{p,g - 1} = 1 \wedge v_{p,g} = 1) \implies w_{p,g} = 0 \quad \forall p, \forall g > 0.$$
• Two pods are opponents if and only if they are playing in the same game with different colors: $$y_{pqg} = 1 \iff (v_{pg} = 1 \wedge v_{qg} = 1 \wedge z_{pg} \neq z_{qg} \quad \forall p, \forall q\neq p, \forall g.$$
• Every pair of pods faces each other exactly twice: $$\sum_g y_{pqg} = 2 \quad \forall p, \forall q\neq p.$$

Here "allDiff()" is the CPLEX notation for the "all different" constraint, a global constraint common to most constraint solvers. Given a collection of variables with the same range, the all different constraint says that no two variables in the collection can take the same value.

To put it mildly, I would not be shocked if there is a more efficient (easier to solver) way to write the problem as a CP model.

 

Tournament Scheduling: CPLEX v. CP Optimizer (Part 2)

In my previous post, I described a tournament scheduling problem and discussed my results using both an integer programming model (CPLEX) and a constraint programming model (CP Optimizer) to try to solve it. Here I will present the IP model. (This post will be gibberish unless you have read at least the start of the previous post.)

There are multiple ways to write an integer programming (IP) model for the tournament problem. My initial instinct was to think in terms of what games a pod plays in, which pods it plays with/against and so on, but I fairly quickly changed to thinking in terms of teams rather than pods. With $N=9$ pods and two pods to a team, there are $N(N-1)/2=36$ possible teams, which we can enumerate at the outset. There will be 18 games, each containing two teams, and every team will play exactly once. Each game contains two "slots", one for the team in white jerseys and the other for the team in black jerseys.

In what follows, I will (I hope) stick to the following notation:

• $t$ indexes a team;
• $p$ and $q$ index pods;
• $g$ indexes a game;
• $s$ indexes a slot (and $s'$ is the other slot of the same game);
• $c$ is a jersey color (0 for white, 1 for black).

Since I am using Java, all indexing is zero-based. Slots are numbered so that game 0 has slots 0 (white) and $G$ (black), game 1 has slots 1 (white) and $G+1$ (black), etc., where $G=N(N-1)/4$ is the number of games. I will denote by $T_p$ the set of teams $t$ containing pod $p$.

The model variables are as follows.

• $x_{ts}\in \lbrace 0,1\rbrace$ is 1 if team $t$ plays in slot $s,$ 0 if not.
• $y_{pqs}\in \lbrace 0,1\rbrace$ is 1 if pod $p$ is playing in slot $s$ and is opposed by pod $q.$
• $z_{pg}\in \lbrace 0,1\rbrace$ is the color jersey worn by pod $p$ at the conclusion of game $g$ (regardless of whether or not $p$ played in game $g$).
• $w_{pg}\in \lbrace 0,1\rbrace$ is 1 if pod $p$ changes jersey color going into game $g,$ 0 if not.

The objective is to minimize $\sum_p \sum_g w_{pg}.$ The constraints are the following.

• Every team plays exactly once. $$\sum_s x_{ts} = 1 \quad\forall t.$$
• Every schedule slot is filled exactly once. $$\sum_t x_{ts} = 1 \quad\forall s.$$
• Pods $p$ and $q$ oppose each other with $p$ in slot $s$ if and only if $p$ is playing in $s$ and $q$ is playing in $s'$. $$y_{pqs} \le \sum_{t\in T_p} x_{ts} \quad\forall p, q\neq p, s.$$ $$y_{pqs} \le \sum_{t\in T_q} x_{ts'} \quad\forall p, q\neq p, s.$$ $$y_{pqs} \ge \sum_{t\in T_p} x_{ts} + \sum_{t\in T_q} x_{ts'} - 1 \quad\forall p, q\neq p, s.$$
• Every pair of pods opposes each other exactly twice. $$\sum_s y_{pqs} = 2 \quad \forall p\neq q.$$
• No pod plays consecutive games in different jerseys.$$\sum_{t\in T_p}\left( x_{t,s-1} + x_{t,s'} \right) \le 1 \quad \forall p, \forall s\notin \lbrace 0, G\rbrace.$$
• A pod playing in a game has its jersey color determined by its team's slot. (This has the desirable side effect of preventing two teams containing the same pod from playing against each other.) $$\sum_{t\in T_p} x_{t, g+G} \le z_{pg} \le 1 - \sum_{t\in T_p} x_{t,g} \quad \forall p,g.$$(Note that for any game index $g$, slot $g$ is the white slot in game $g$ and slot $g+G$ is the black slot.)
• A pod's color for any game (playing or not) is the same as its color for the previous game, unless a change occurs. $$z_{p, g-1} - w_{pg} \le z_{pg} \le z_{p, g-1} + w_{pg} \quad\forall p, \forall g \ge 1.$$

In the next post, I will try to articulate my CP model for the problem.

Tournament Scheduling: CPLEX v. CP Optimizer (Part 1)

 A recent question on Mathematics Stack Exchange asked about scheduling $N$ "pods" (players) in a tournament under the following conditions.

• There will be $N(N-1)/4$ games played sequentially (no games in parallel).
• Games pit teams of two pods each against each other. One team wears white and one team wears black.
  
• Each pod partners with every other pod once and plays against every other pod twice.
• No pod can play in consecutive games wearing different jersey colors.

The objective is find a schedule that minimizes the total number of times pods have to change their jersey color.

The question specifies $N=9$. Note that, since there are $N(N-1)/4$ games, a necessary condition for the problem to be feasible is that either $N$ or $N-1$ must be divisible by $4$. That condition is not, however, sufficient. An obvious counterexample is when $N=4$ and there are three games to be scheduled. The first game uses all four teams, as does the second game. Since no pod can be forced to wear different jerseys in consecutive games, all four teams would go into the second game wearing the same colors as in the first game ... which means being in the same teams, violating the rule about pods being teamed together exactly once.

I coded both an integer programming model (using CPLEX) and a constraint programming model (using CP Optimizer) in Java and tried to solve the tournament problem with $N=9$. Since everything is inherently integer and the constraints are more logical than arithmetic (the only real algebra is summing up the number of jersey changes), I had two initial expectations: that CPLEX would provide tighter lower bounds than CPO (because MIP models tend to do better with bounds than CP models); and that CPO would find feasible (and possibly optimal) schedules faster than CPLEX would, since the problem is inherently logic-based (and CP models often do better than MIP models on scheduling problems). I was half right. CPLEX did provide tighter lower bounds than CPO, at least within the time limits I was willing to use, although I don't think either came remotely close to a "tight" lower bound. CPO, however, struggled massively to find feasible schedules while CPLEX did not have too much trouble doing so.

Before going any further, I need to issue three caveats. First, the longest I ran the IP model was 30 minutes and the longest I ran the CP model was an hour. Second, while I am pretty familiar with formulating IP models, I am much less familiar running CP models, so I may have missed opportunities to make the CP model more efficient. Third, while CPLEX gives the user a fair amount of indirect control over the search process (by setting branching priorities, MIP emphasis, how frequently to try certain kinds of cuts and so on), CP Optimizer may offer the user even more flexibility in how to tailor searches I am not yet familiar enough to do anything beyond trying to indicate which variables should be the first candidates for branching (and I'm not sure I got that right).

I'll end this installment with a synopsis of some runs.

• In a 30 minute run using the new setting 5 for the MIP emphasis parameter (stressing use of heuristics to improve the incumbent, at the cost of not paying too much attention to the lower bound), CPLEX found a feasible schedule with 14 jersey changes and a lower bound of 2.97 (a 79% gap). It evaluated somewhere around 18,500 nodes.
  
• In a 30 minute run, CP Optimizer branched over 957 million times but never found a feasible schedule and ended with a best bound of 0.
• Finally, I tried running CPLEX for three minutes (in which it found an incumbent with 19 jersey changes) and then used that incumbent as a starting solution for a one hour run of CP Optimizer. CPO accepted the incumbent solution, did a bit over 1.55 billion branch operations, and failed to improve on it. The best bound was again 0. 

If you want to look at my code (which will make slightly more sense after my next couple of posts, where I will describe the models), it is available from my university GitLab repository.

Taxi Dispatch

A question on OR Stack Exchange pertains to a taxi stand. The assumptions are as follows.

• There are $t$ taxis operating between points A and B, with all passenger traffic going from A to B. For concreteness, I'll refer to a "taxi stand" with a "line" of taxis, similar to what you might see at an airport or popular hotel.
• The round-trip travel time for a taxi (taking passengers from A to B and deadheading back from B to A) is deterministic with value $T$. (Note: Deterministic travel time is the author's assumption, not mine.)
• Each taxi has a capacity for $p$ passengers and is dispatched only when full.
• Passengers arrive according to a Poisson process with parameter (rate) $\lambda$.

The author was looking for the mean time between dispatches.

It has been a very long time since I taught queuing theory, and I have forgotten most of what I knew. My best guess for the mean time between dispatches, under the assumption that the system reached steady state, was $p/\lambda$, which would be the long-run average time required for $p$ passengers to arrive. To achieve steady-state, you need the average arrival rate ($\lambda$) to be less than the maximum full-blast service rate ($t p/T$). If arrivals occur faster than that, in the long term you will dispatching a taxi every $1/T$ time units while either queue explodes or there is significant balking/reneging.

To test whether my answer was plausible (under the assumption of a steady state), I decided to write a little discrete event simulation (DES). There are special purpose DES programs, but it's been way to long since I used one (or had a license), so I decided to try writing one in R. A little googling turned up at least a couple of DES packages for R, but I did not want to commit much time to learning their syntax for a one-off experiment, plus I was not sure if any of them handled batch processing. So I wrote my own from scratch.

My code is not exactly a paragon of efficiency. Given the simplicity of the problem and the fact that I would only be running it once or twice, I went for minimal programmer time as opposed to minimal run time. That said, it can run 10,000 simulated passenger arrivals in a few seconds. The program uses one data frame as what is known in DES (or at least was back in my day) an "event calendar", holding in chronological order three types of events: passengers joining the queue; taxis being dispatched; and taxis returning.

The simulation starts with all taxis in line and no passengers present. I did not bother with a separate "warm-up" period to achieve steady-state, which would be required if I were striving for more accuracy. For a typical trial run ($t=5$, $p=4$, $T=20$, $\lambda=0.7$), my guess at the mean time between dispatches worked out to 5.714 and the mean time computed from the run was 5.702. So I'm inclined to trust my answer (and to quit while I'm ahead).

If anyone is curious (and will forgive the lack of efficiency), you can take a look at my R notebook (which includes the code).