Lets define the following action :
\begin{equation}\tag{1}
S' = \int_{t_1}^{t_2} L' \, dt,
\end{equation}
where
\begin{equation}\tag{2}
L' = L - \lambda \, \frac{d}{dt} (q^k \, p_k).
\end{equation}
An arbitrary variation $\delta q^k$ gives the following :
\begin{align}
\delta S' &= \int_{t_1}^{t_2} \Big( \frac{\partial L}{\partial q^i} \; \delta q^i + \frac{\partial L}{\partial \dot{q}^i} \; \delta \dot{q}^i \Big) \, dt - \lambda \, \delta (q^k \, p_k) \Big|_{t_1}^{t_2} \\[12pt]
&= \int_{t_1}^{t_2} \Big[ \frac{\partial L}{\partial q^i} - \frac{d}{dt} \Big( \frac{\partial L}{\partial \dot{q}^i} \Big) \Big] \, \delta q^i \; dt + \int_{t_1}^{t_2} \frac{d}{dt} \Big( \frac{\partial L}{\partial \dot{q}^i} \; \delta q^i \Big) \, dt - \lambda \, \delta (q^k \, p_k) \Big|_{t_1}^{t_2} \\[12pt]
&= (\text{usual Euler-Lagrange variation of } S) + \big(p_i \; \delta q^i - \lambda \, \delta (q^k \; p_k)\big)\Big|_{t_1}^{t_2} \\[12pt]
&= (\ldots) + \big( (1 - \lambda) \, p_k \; \delta q^k - \lambda \, q^k \, \delta p_k \big) \Big|_{t_1}^{t_2}. \tag{3}
\end{align}
The first term cancels from the Euler-Lagrange equations. To find a proper quantity that is fixed at the endpoints $t_1$ and $t_2$, the last term need to be written like this (look at the position of the indices. The bar is here on purpose, since there are two solutions) :
\begin{equation}\tag{4}
(1 - \lambda) \, p_k \; \delta q^k - \lambda \, q^k \, \delta p_k = -\; Q^k \; \delta P_k \quad \text{or} \quad \bar{P}_k \; \delta \bar{Q}^k.
\end{equation}
We could then fix $P_k$ or $\bar{Q}^k$ at $t_1$ and $t_2$ : $\delta P_k(t_1) = 0$ and so on. The two solutions are apparently a consequence of the symmetry between the variables $q^k$ and $p_k$ in the product $q^k \, p_k$. From the indices position, we can guess the following :
\begin{align}
Q^k &= \mathcal{Q}(q, p) \, q^k,
& P_k &= \mathcal{P}(q, p) \, p_k, \tag{5} \\[18pt]
\bar{Q}^k &= \bar{\mathcal{Q}}(q, p) \, q^k,
& \bar{P}_k &= \bar{\mathcal{P}}(q, p) \, p_k. \tag{6}
\end{align}
From dimensional ground and indices position, we could expect that $\mathcal{P}(q, p) \propto (q^i \, p_i)^a$ and $\bar{\mathcal{Q}}(q, p) \propto (q^i \, p_i)^b$, where $a$ and $b$ are unknown exponents to be found as functions of $\lambda$. The variation of $P_k$ and $\bar{Q}^k$, substituted into the equation above, give us some simple algebraic equations that are easy to solve :
\begin{equation}
-\; a \, \mathcal{Q} \, (q^i \, p_i)^a \, p_k \; \delta q^k - (1 + a) \, \mathcal{Q} \, (q^i \, p_i)^a \, q^k \, \delta p_k = (1 - \lambda) \, p_k \; \delta q^k - \lambda \, q^k \, \delta p_k,
\end{equation}
or :
\begin{equation}
(1 + b) \, \bar{\mathcal{P}} \, (q^i \, p_i)^a \, p_k \; \delta q^k + b \, \bar{\mathcal{P}} \, (q^i \, p_i)^a \, q^k \, \delta p_k = (1 - \lambda) \, p_k \; \delta q^k - \lambda \, q^k \, \delta p_k.
\end{equation}
The complete algebraic solutions are then $a = \lambda - 1$, $b = - \lambda$ and
\begin{align}
\mathcal{Q}(q, p) &= (q^i \, p_i)^{1 - \lambda},
& \mathcal{P}(q, p) &= (q^i \, p_i)^{\lambda - 1}, \\[18pt]
\bar{\mathcal{Q}}(q, p) &= (q^i \, p_i)^{- \lambda},
& \bar{\mathcal{P}}(q, p) &= (q^i \, p_i)^{\lambda}.
\end{align}
Finally, the quantities that are fixed at the endpoints are these :
\begin{equation}\tag{7}
P_i = (q^k \, p_k)^{\lambda - 1} \, p_i,
\end{equation}
or :
\begin{equation}\tag{8}
\bar{Q}^i = (q^k \, p_k)^{- \lambda} \, q^i.
\end{equation}
Now I'm still puzzled as why there are two solutions. Also, since $q^i \, p_i$ may turn out negative, these quantities may be ill-defined for some values of $\lambda$ ($\lambda = \tfrac{1}{2}$, for example).
Some special cases :
If $\lambda = 0$, the natural solution is given by the second one : $\bar{Q}^i = q^i$ (Dirichlet conditions). Why is there another solution (first one above) :
\begin{equation}\tag{9}
P_i = \frac{p_i}{q^k \, p_k} \; ?
\end{equation}
If $\lambda = 1$, then the natural solution is given by the first one : $P_i = p_i$. Why is there another solution (second one above) :
\begin{equation}\tag{10}
\bar{Q}^i = \frac{q^i}{q^k \, p_k} \; ?
\end{equation}
Someone confirms the calculations above, and have an interpretation of the dual solutions ? And what about the sign problem, for some specific values of $\lambda$ ?
EDIT 1 : It is funny to notice that the $Q^i$ and $P_i$ (or $\bar{Q}^i$ and $\bar{P}_i$) defined above (expressions (5) or (6)) are canonical transformations of the $q^i$ and $p_i$ :
\begin{align}\tag{11}
\{ Q^i, \, Q^j \} &= 0,
&\{ P_i, \, P_j \} &= 0,
&\{ Q^i, \, P_j \} &= \delta^i_j.
\end{align}
EDIT 2 : The expression $q^k \, p_k$ is invariant under the canonical transformation
\begin{align}\tag{12}
Q^k &= (q^i p_i)^{-1} \, q^k,
& P_k &= (q^i p_i) \, p_k, \\[12pt]
q^k &= (Q^i P_i) \, Q^k,
& p_k &= (Q^i P_i)^{-1} \, P_k, \tag{13}
\end{align}
so $q^k \, p_k \equiv Q^k \, P_k$, and it's easy to verify the identity
\begin{equation}\tag{14}
- \; q^k \; \delta p_k \equiv P_k \; \delta Q^k.
\end{equation}
So it is not surprising after all that there's a second solution to the problem. This $Q^k$ may be interpreted as another quantity that is fixed at the boundaries, in the case of $\lambda = 1$ ($p_k$ and $Q^k$ are both "natural" quantities that are fixed in this case).