Tag: optimization

I’d like to start by pointing out that Hector Pefo also posted a nice solution to this problem. There were also diagrams of solutions posted on Twitter, such as those from @DimEarly and from @chattinthehatt. The solution I’ll describe here is similar to Hector’s, but I’ll explain how to model it as an integer linear program and solve it that way. Of course both solutions produce the same answers (because they’re both correct!) but it’s nice to see different approaches. So without further ado…

Integer linear program

Let’s suppose there are $N$ shapes and we’d like to use $C$ colors. Define:
\[
x_{ij} = \begin{cases}1&\text{if shape }i\text{ is colored using color }j\\0&\text{otherwise}\end{cases}
\]The first constraint is that each shape must be assigned to exactly one color. We can express this algebraically as follows:
\[
\sum_{j=1}^C x_{ij} = 1\qquad\text{for }i=1,\dots,N
\]Now suppose that shape $j$ has area $a_j$. The constraint that each color has an equal share of the area means that the sum of the areas of the shapes of each color must be $A_\text{tot}/C$, where $A_\text{tot}$ is the total area of the Ostomachion. Algebraically, we can write this as:
\[
\sum_{i=1}^N a_i x_{ij} = A_\text{tot}/C\qquad\text{for }j=1,\dots,C
\]Finally, we must deal with the edge constraints. We say that $(p,q)$ is an “edge” if shape $p$ shares an edge with shape $q$. Shapes that share an edge must be different colors. We can encode this algebraically as:
\[
x_{pj} + x_{qj} \le 1\qquad\text{for all edges }(p,q)\text{ and for }j=1,\dots,C
\]And that’s it! We described the feasible set of solutions as the set of matrices $x \in \{0,1\}^{N\times C}$ that satisfy the three sets of linear constraints above. This sort of optimization model is a modified version of a Knapsack problem. More generally, it’s an Integer Linear Program (ILP). Such problems are known to be NP-hard, which roughly means that in the worst-case, there is no efficient way to solve them short of trying all possible colorings. However, modern ILP solvers are quite sophisticated and can usually solve these types of problems very efficiently.

Finding all solutions

When using an ILP model, most commercial solvers are designed to find a single solution. For this problem, we are interested in finding all solutions. One way to get around the limitation is to use Barvinok’s algorithm. One can also use variants of branch-and-cut. Yet another way, which is not as systematic but much easier to explain and implement, is to use a randomized approach.

Since all our variables are binary (0 or 1), the convex hull of the feasible set cannot contain any interior solution points. All solutions will be on the boundary. This means that any solution can be found by specifying the right objective function for the optimization. As stated, our optimization is just a feasibility problem and there is no objective specified. For this problem, we start by letting $r_{ij}$ be random numbers (we can draw them i.i.d. from a standard normal distribution, for example). Then, we simply add the objective function:
\[
\text{Minimize}\qquad \sum_{i=1}^N\sum_{j=1}^C r_{ij} x_{ij}
\]and solve the problem. If we do this repeatedly with different realizations of the $\{r_{ij}\}$, we will eventually enumerate all possible solutions of the optimization problem.

Solutions

I implemented the ILP described earlier using Julia and JuMP and I solved everything using the Gurobi solver. Although the randomized approach for finding all solutions works well, Gurobi has built-in functionality that allows one to directly find all solutions. In addition to doing this, I also included logic that removes duplicate solutions arising from different re-colorings of the same solution. e.g. if you take a solution and just swap red for green, I counted that as the same solution. Here are the solutions I found. I included the time it took Gurobi (on my laptop) to find all solutions for each case:

For $C=3$ colors, there are $2$ distinct solutions (took 1.6 ms to solve)

For $C=4$ colors, there are $9$ distinct solutions (took 32 ms to solve):

For $C=6$ colors, there are $33$ distinct solutions (took 700 ms to solve):

Interestingly, using fewer colors makes it easier to satisfy the equal-area constraint, but it makes it harder to satisfy the edge constraints. Ultimately, there are only solutions for the cases above ($C=3,4,6$). There are no solutions for any other values of $C$.

If you are interested in seeing my code, here is my IJulia notebook.

Let’s give names to the variables so that the unit conversions don’t become cumbersome. Define:
\begin{align}
x_\text{light} &= \text{position of the traffic light (2 km)} \\
x_\text{end} &= \text{position of the finish line (4 km)} \\
v_\text{max} &= \text{speed limit (100 km/h)} \\
a_\text{max} &= \text{maximum acceleration (2 m/sec}^2\text{)} \\
\tau_y &= \text{time that a yellow light lasts (5 sec)} \\
\tau_r &= \text{time that a red light lasts (20 sec)}
\end{align}

It is implicitly assumed in the problem that we can never run a red light. So even though we don’t know when the light will turn yellow, we can’t take any chances; we must make sure that even the unlikeliest most ill-timed yellow light can never cause us to run a red.

There are many different approaches to solving this problem but since I’m currently teaching a class on optimization modeling, I figured I would go that route. The idea behind a modeling approach is to imagine that the position, velocity, and acceleration of the car are decision variables that must satisfy several constraints.

Basic variables and constraints

Let’s discretize the time interval and suppose we make decisions at every $\Delta t$ seconds. So if $\Delta t = 1$, we make decisions every second. Then, define the following decision variables:
\begin{align}
x_t & \ge 0 && \text{position at time }t\\
0 \le v_t &\le v_\text{max} && \text{speed at time }t\\
-a_\text{max} \le a_t &\le a_\text{max} && \text{acceleration at time }t
\end{align}where the time indices range over $t=0,1,\dots,T$ where $T$ is sufficiently large. Of course, we cannot choose positions and velocities arbitrarily; they must be consistent. By approximating the acceleration as constant over each $\Delta t$ step, we may write that:
\begin{align}
v_{t+1} &= v_t + a_t \Delta t && \text{(change in speed with constant accel.)}\\
x_{t+1} &= x_t + v_t\Delta t + \tfrac{1}{2}a_t \Delta t^2 && \text{(change in pos. with constant accel.)}
\end{align}We also have initial conditions for position and speed:
\[
x_0 = 0 \qquad\text{and}\qquad v_0 = v_\text{max}
\]Note that we constrained the speed to be nonnegative. It turns out this doesn’t matter and there is nothing to be gained by allowing backward movement. We get the same answer either way!

Contingency planning

Since we can never risk running a red light, we must ensure that at every instant, one of two things must happen. Either:

1. We can reach $x_\text{light}$ within $\tau_y$ sec. (we run the yellow) or:
2. We never pass $x_\text{light}$ within $\tau_y+\tau_r$ sec. (we can avoid running a red).

These constraints can be encoded algebraically. For each $t$, We’ll define a new set of positions, velocities, and accelerations that last for $\tau_y+\tau_r$ seconds. These are contingency variables that allow us to ensure nothing can go wrong if the light turns yellow at time $t$. So we define the variables:
\begin{align}
\hat x_{t,k} & \ge 0 && \text{(contingency position)}\\
0 \le \hat v_{t,k} &\le v_\text{max} && \text{(contingency velocity)}\\
-a_\text{max} \le \hat a_{t,k} &\le a_\text{max} && \text{(contingency acceleration)}
\end{align}These new contingency variables are defined over $0\le t \le T$ and $0 \le k \le \tau_y+\tau_r$. Just like the ordinary position and velocity, they must satisfy contingency constraints:
\begin{align}
\hat v_{t,k+1} &= \hat v_{t,k} + \hat a_{t,k} \Delta t && \text{(change in speed with const. accel.)}\\
\hat x_{t,k+1} &= \hat x_{t,k} + \hat v_{t,k} \Delta t + \tfrac{1}{2}\hat a_{t,k} \Delta t^2 && \text{(change in pos. with const. accel.)}
\end{align}We also have initial conditions for position and speed:
\[
\hat x_{t,0} = x_t \qquad\text{and}\qquad \hat v_{t,0} = v_t
\]Now that we have the contingency trajectories set up, we must encode our binary constraint. Either we can run the yellow, or we are safe from running the red. Let’s define binary variables:
\begin{align}
z^y_t \in \{0,1\} & \qquad \text{(can we run the yellow if light turns yellow at }t\text{?)}\\
z^r_t \in \{0,1\} & \qquad \text{(are we safe from running the red if light turns yellow at }t\text{?)}
\end{align}The constraints that we can run the yellow and not run a red are:
\begin{align}
\hat x_{t,\tau_y} &\ge x_\text{light} z^y_t &&\quad\text{for all }t\text{, and}\\
\hat x_{t,\tau_y+\tau_r} &\le x_\text{light} + (1-z^r_t)x_\text{end} &&\quad\text{for all }t\text{, respectively.}
\end{align}together with the logic constraint $z^y_t + z^r_t \ge 1$ for all $t$, which ensures that we can run the yellow OR we are safe from running the red. And that’s it!

Objective function

The work above describes the feasible set of position, velocity, and acceleration profiles that guarantee we never break the law. Among these profiles, we should choose the one that gets us to the finish line as fast as possible on average. I struggled to find a simple objective function that captures this notion precisely. As a compromise, I simply minimized the duration of the trip assuming the light never turns red. This should be a decent approximation to the true solution because yellow lights are relatively rare. Most of the time when the light turns yellow, it doesn’t actually affect us; it’s only a concern if the light turns yellow when we are close to it, which won’t happen all that often.

To characterize the finish time, I defined the binary variables:
\[
z^e_t \in \{0,1\} \qquad \text{(are we finished at time }t\text{?)}
\]together with the constraint:
\[
x_t \ge x_\text{end} z^e_t\qquad\text{for all }t.
\]Now, our binary variable will be $1$ whenever we have crossed the finish line. So by maximizing $\sum_{t=0}^T z^e_t$, we get to the end as soon as possible.

Solution

From a modeling standpoint, all our constraints are linear, which is nice. We do have some binary variables, which makes this a mixed-integer linear program (MILP). Nevertheless, the problem is relatively easy to solve using a standard solver. For this problem, I coded the model in Julia together with the JuMP package and the Gurobi solver. Some details:

• Solving with $\Delta t = 1$. After presolve, the problem had: 3492 constraints, 6216 continuous variables, and 108 binary variables. Solving took about 0.3 seconds on a laptop.
• Solving with $\Delta t = 0.5$. After presolve, the problem had: 15499 constraints, 29602 continuous variables, and 216 binary variables. Solving took about 1.3 seconds on a laptop.
• Solving with $\Delta t = 0.25$. After presolve, the problem had: 63925 constraints, 125250 continuous variables, and 436 binary variables. Solving took about 12 seconds on a laptop.

No matter the time discretization, all solutions looked the same. Here is a plot of the solution for $\Delta t = 0.25$, zoomed in near the point where the car approaches the traffic light.

The plot looks constant outside this range. So, the optimal thing to do is maintain a speed of 100 km/h for 65 seconds, then decelerate as hard as possible for 2.5 seconds, then accelerate as hard as possible for 2.5 seconds, and then continue on at 100 km/h. The reason for the dip in velocity is that we would otherwise be unable to stop in time if the light were to turn yellow, so we must be risk-averse. After 67.5 seconds, we can safely accelerate, knowing that if the light turns yellow we will make it within the next 5 seconds.

If you’re interested in seeing my code in full detail, you can view/download my notebook here.