As he regularly does, Jason Cole recently delighted the readers of his blog with an interesting problem that I will quickly sum up.

The problem

You have been assigned the task to use a given amount of concrete $A$ to build a speed bump over a section of a road of length $2\\,\ell$ . But being a regular user of the road, you are trying to do your job while minimizing the jolt felt by a driver passing above it.

We will consider the road as a 2D curve and the car as a point-like particle on this curve. We are then searching for a function $x \mapsto y(x)$ , with domain and values $[-\ell,\ell] \to \mathbb{R}$ . We will assume that the road outside of this interval lies at $y=0$ ; the continuity conditions are then $y(-\ell) = y(\ell) = 0$ , and we will also impose the nullity of the first derivatives of $y$ at these two points, which we will note with a parenthesized subscript $y_{(x)}(-\ell) = y_{(x)}(\ell) = 0$ .

The author chose to minimize the integral of the square of the vertical acceleration felt by the car passing over the bump, $\int_0^T{a_y(t)^2\mathrm{d}t}$ .

Analytical results

Going quickly over the details because they already are in the original post, conservation of energy and geometrical considerations give

which, derived with respect to time, yields

The quantity that we want to minimize is then

The author then proceeds to two approximations. The first one is the limit of "infinite velocity", which seems more appropriate for some drivers than some others, but is actually justified for everyone if we compute the quantity $\frac{{v_0}^2}{2\\,g\\,y} \approx 35 \gg 1$ . This approximation allows to set $g$ to $0$ :

The second one is that of small angles, consisting of neglecting $y_{(x)}$ in front of $1$ . This ones strips the criterion even further:

such that the problem becomes tractable. Indeed, applying the Euler-Lagrange equation on the Lagrangian $\mathcal{L}(y,y_{(x)},y_{(x,x)}) = {v_0}^3{y_{(x,x)}}^2-\lambda y$ yields the optimal solution of this twice simplified problem:

He finally programs a stochastic optimizer with the exact equations, and looks how this solution changes. He didn't observe many changes, and concluded that the approximations were not too coarse.

Numerical optimization of the exact problem

I decided to challenge this observation, and proceeded to submit the problem to Bocop, an optimal control solver developed by the Inria team COMMANDS. Bocop takes as input a system of ordinary differential equations involving one or several control variables, a criterion to optimize, and possibly constraints, and determines the temporal profile of the control variables that optimizes the criterion while satisfying the constraints.

For example, it is able to determine how a model micro-swimmer should twist in order to move fast. Or how a rocket should adapt its thrust during flight to go the highest.

Bocop is designed to models dynamical systems (with time derivatives), and even though the equation of motion of the car involves time derivatives, it seemed more natural here to remove the time and only reason with the independent variable $x$ . Since we go numeric, let's make the problem dimensionless with the change of variables $\hat{x} = \frac{x}{\ell}$ , $\hat{y} = y\frac{\ell}{A}$ , $\hat{A} = \frac{A}{\ell^2}$ , $\hat{g} = g\frac{\ell}{{v_0}^2}$ . The criterion becomes

and from now on, I will drop all the hats and work only in reduced variables.

The system can be modelled with a system of four first-order ordinary differential equations: $y(x)$ the profile of the road, $y_{(x)}(x)$ its first derivative, $c(x)$ the criterion and $Y(x) = \int_{-1}^x{y(u)\mathrm{d}u}$ the quantity of concrete used. The control variable is $y_{(x,x)}$ , the second derivative of the profile. The dynamics is straightforward:

and is implemented in the dynamics.tpp file

Tdouble yxx = control[0];
Tdouble y = state[0];
Tdouble yx = state[1];
Tdouble c = state[2];
Tdouble Y = state[3];
double g = constants[0];
double A = constants[1];

Tdouble vx = sqrt((1-2*A*g*y)/(1+A*A*yx*yx));
Tdouble ay = yxx/(1+A*A*yx*yx) - A*g*yx*yx/(1-2*A*g*y);

state_dynamics[0] = yx;             // dy/dx
state_dynamics[1] = yxx;            // dyx/dx
state_dynamics[2] = vx*vx*vx*ay*ay; // dc/dx
state_dynamics[3] = y;              // dY/dx

The criterion is the final state of $c$ :

criterion = final_state[2];

The boundary conditions take care of the continuity conditions, initialize the criterion, and enforce the area constraint on the profile:

double A = constants[1];
boundary_conditions[0] = initial_state[0];      // y(-1)  = 0
boundary_conditions[1] =   final_state[0];      // y(1)   = 0
boundary_conditions[2] = initial_state[1];      // yx(-1) = 0
boundary_conditions[3] =   final_state[1];      // yx(1)  = 0
boundary_conditions[4] = initial_state[2];      // c(-1)  = 0
boundary_conditions[5] = initial_state[3];      // Y(-1)  = 0
boundary_conditions[6] =   final_state[3] - 1;  // Y(1)   = 1

Finally we can run the solver, let's first check if we can reach the analytical solution of the simple problem (Tdouble vx = 1, ay = yxx;):

It seems to work, and the objective value $22.5$ is consistent with what we could have expected from the analytical result. Now let's try the exact problem, with different values of $A$ and $g$ :

Note that "realistic" values of $A$ and $g$ could be, respectively, $0.1$ and $0.15$ . In this regime, the two approximations (infinite speed, and small angles) are indeed very reasonable. However, for higher bumps, the exact solution starts to diverge noticeably from the approximation. Could we guess the analytical form of the solution? That seems difficult in the general case, but we could try it in the limit of infinite speed ($g=0$ ). Actually, the objective value for the solution calculated by Bocop for the condition $A=1$ , $g=0$ is $6.30$ , which is suspiciously close to an interesting number, and the profile itself looks like four arcs of circles. Coincidence?

Hypothesis for the infinite speed limit

Let's do the maths in the infinite speed limit, for the function

The objective is reduced to $\mathrm{Obj}\_{g=0}[y] = \int_{-1}^{1}{\frac{{y_{(x,x)}}^2}{(1+A^2{y_{(x)}}^2)^{7/2}}\mathrm{d}x}$ .

Mathematica helps and indeed $\mathrm{Obj}_{A=1,g=0}[\mathrm{circles}] = 2\\,\pi$ . To know if it is actually the optimal profile, we need to plug it into the Euler-Lagrange equation. Spoilers, it is not.

This profile is not optimal, but we can at least compare it with $\mathrm{approx}(x) = \frac{15}{16}(1-x^2)^2$ , the analytical solution of the approximated problem. Using the same criterion, we get $\mathrm{Opt}_{A=1,g=0}[\mathrm{approx}] = 10.39$ , which proves that the four circles are a better profile for the limit of infinite speed, at $A=1$ . We can also compare these profiles as a function of $A$ :

It was expected that the polynomial would win at low $A$ since low $A$ means small angles. As we have seen, the solution in circles is better at high $A$ . What was more difficult to understand for me was that the circles become infinitely bad at low $A$ (like $\frac{16}{15\\,a^4}$ ). I think that the reason for this behaviour is the two points of infinite slope at $x= \pm\frac{1}{2}$ , which never disappear even when the profile gets flattened by a decreasing $A$ , and create two discontinuity points where most of the acceleration is taken.

Can we find a better solution? Is the optimal solution tractable in the limit of infinite speed? What if we wanted to minimize the maximum of the acceleration over the period (the infinite norm) instead of the L2 norm? Perhaps for another time!