One might be able to get your second approach, with integrator, to work, however I would think adding feed-forward might be an easier solution. Namely when you reach your reference goal your system has to satisfy $C\,x_r + D\,u_r = r$ and $\dot{x}_r = A\, x_r + B\,u_r = 0$, or when writing it as one system of linear equations
$$
\begin{bmatrix}
A & B \\ C & D
\end{bmatrix}
\begin{bmatrix}
x_r \\ u_r
\end{bmatrix} =
\begin{bmatrix}
0_{12\times1} \\ r
\end{bmatrix}.
$$
Assuming that that matrix is invertible, then solving for $x_r$ and $u_r$ yields
$$
\begin{bmatrix}
x_r \\ u_r
\end{bmatrix} =
\begin{bmatrix}
A & B \\ C & D
\end{bmatrix}^{-1}
\begin{bmatrix}
0_{12\times1} \\ r
\end{bmatrix}.
$$
By defining new state space coordinates and input $\hat{x} = x - x_r$ and $\hat{u} = u - u_r$ respectively, then its dynamics can be written as
$$
\dot{\hat{x}} = \dot{x} - \dot{x}_r = A\, x + B\, u = A (\hat{x} + x_r) + B (\hat{u} + u_r) = A\, \hat{x} + B\, \hat{u} + \underbrace{A\, x_r + B\,u_r}_{0}.
$$
So $x$ can be driven to $x_r$ ($\hat{x} = 0$) by using the control law $\hat{u} = -K\, \hat{x}$ (for all $K$ such that $A-B\,K$ has all its eigenvalues in the open left half plane). This implies that the actual input $u$ should be equal to
$$
u = \hat{u} + u_r = -K (x - x_r) + u_r.
$$
When substituting in the solution for $x_r$ and $u_r$ yields
$$
u = -K\,x +
\begin{bmatrix}
K & I_{4\times4}
\end{bmatrix}
\begin{bmatrix}
A & B \\ C & D
\end{bmatrix}^{-1}
\begin{bmatrix}
0_{12\times4} \\ I_{4\times4}
\end{bmatrix} r.
$$
Since all these matrices are constant, then by defining a new matrix as
$$
F =
\begin{bmatrix}
K & I_{4\times4}
\end{bmatrix}
\begin{bmatrix}
A & B \\ C & D
\end{bmatrix}^{-1}
\begin{bmatrix}
0_{12\times4} \\ I_{4\times4}
\end{bmatrix}
$$
you can see that with your first attempt you where not that far off from this feed-forward method in terms of control law.
It can be noted that this method assumes that $r$ is constant, however when the reference changes over time, then your tracking performance might not be great. And after some playing with this myself it even seemed that choosing a $K$, such that the real parts of the poles of $A-B\,K$ are more negative, negatively affected the tracking performance a high frequencies.
It might also be worth noting that one can define a new state space model which captures the resulting reference-output dynamics as follows:
$$
\left\{
\begin{array}{l}
\dot{x} = \underbrace{(A - B\,K)}_{A_r}\,x + \underbrace{B\,F}_{B_r}\,r \\
y = \underbrace{(C - D\,K)}_{C_r}\,x + \underbrace{D\,F}_{D_r}\,r
\end{array}\right.
$$
I did a quick search online and found this document which basically contains the same derivation of the feed-forward method. It also contains the integral method, however it stops once they defined the extended state space and only suggested a control law $u = K\,x - K_I\,e$, but not how to find such $K$ and $K_I$.
I was determined to also find a solution to the integrator approach myself. But ran into difficulties using only the state space extension of $\dot{e} = r - y$. After some trial and error I decided on an extend state space of
$$
z =
\begin{bmatrix}
x \\ e \\ r \\ u
\end{bmatrix}
$$
By defining a new input to this system as $v = \dot{u}$ then the dynamics can be formulate as
$$
\dot{z} = \underbrace{
\begin{bmatrix}
A & 0 & 0 & B \\
C & 0 & -I & D \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0
\end{bmatrix}}_{\hat{A}} z + \underbrace{
\begin{bmatrix}
0 \\
0 \\
0 \\
I
\end{bmatrix}}_{\hat{B}} v.
$$
As time goes to infinity one would like $y-r$, so $\dot{e}$, to go to zero, however one would also like the system to be internally stable, this would happen if $\dot{x}$ and $\dot{u}$ would go to zero as well. So a cost function for an optimal controller could be defined as
$$
\begin{array}{rl}
J = & \int \left(\dot{e}^\top Q_e\, \dot{e} + \dot{x}^\top Q_x\, \dot{x} + \dot{u}^\top Q_u\, \dot{u}\right) dt \\
= & \int \dot{z}^\top \underbrace{
\begin{bmatrix}
Q_x & 0 & 0 & 0 \\
0 & Q_e & 0 & 0 \\
0 & 0 & 0 & 0 \\
0 & 0 & 0 & Q_u
\end{bmatrix}}_{\hat{Q}} \dot{z}\, dt
\end{array}
$$
where $Q_e$, $Q_x$ and $Q_u$ are positive definite matrices. Substituting in the expression for $\dot{z}$ gives
$$
J = \int \left(z^\top Q\, z + v^\top R\, v + 2\, z^\top N\, v\right) dt,
$$
$$
Q = \hat{A}^\top \hat{Q}\, \hat{A},
$$
$$
R = \hat{B}^\top \hat{Q}\, \hat{B} = Q_u,
$$
$$
N = \hat{A}^\top \hat{Q}\, \hat{B} = 0.
$$
This just in the standard LQR form. However I had some trouble solving this consistently, because often due to numerically errors $Q$ was not positive semi-definite and $(\hat{A},\hat{B})$ has a boundary stable mode that is not controllable. One way I was able to get this to work with some success was by adding an identity matrix times a very small positive value $\epsilon$ to $Q$ to increase all its eigenvalues by $\epsilon$ and by adding very weak damping to $r$ by setting the third matrix on the diagonal of $\hat{A}$ equal to an identity matrix times a very small negative value. If you did manage to find a controller $v = -K\,z$ with
$$
K =
\begin{bmatrix}
-K_x & -K_e & -K_r & -K_u
\end{bmatrix},
$$
then the closed-loop could also be written as
$$
\begin{bmatrix}
\dot{x} \\ \dot{e} \\ \dot{u}
\end{bmatrix} =
\begin{bmatrix}
A & 0 & B \\
C & 0 & D \\
K_x & K_e & K_u
\end{bmatrix}
\begin{bmatrix}
x \\ e \\ u
\end{bmatrix} +
\begin{bmatrix}
0 \\ -I \\ K_r
\end{bmatrix} r,
$$
$$
y =
\begin{bmatrix}
C & 0 & D
\end{bmatrix}
\begin{bmatrix}
x \\ e \\ u
\end{bmatrix}.
$$
When looking at the relevant CARE I believe it will always be the case that $K_e = 0$. If this is true then the closed-loop can be reduced even further, by removing the state $e$. Doing this yields
$$
\begin{bmatrix}
\dot{x} \\ \dot{u}
\end{bmatrix} =
\begin{bmatrix}
A & B \\
K_x & K_u
\end{bmatrix}
\begin{bmatrix}
x \\ u
\end{bmatrix} +
\begin{bmatrix}
0 \\ K_r
\end{bmatrix} r,
$$
$$
y =
\begin{bmatrix}
C & D
\end{bmatrix}
\begin{bmatrix}
x \\ u
\end{bmatrix}.
$$
After some more searching I also came across what is called Linear-Quadratic-Integral. This indeed does defines tries to find an input of the form $u = -K\,x - F\,e$ which tries to minimize the following cost function
$$
J = \int \left(\begin{bmatrix}x \\ e\end{bmatrix}^\top Q\, \begin{bmatrix}x \\ e\end{bmatrix} + u^\top R\, u + 2\,\begin{bmatrix}x \\ e\end{bmatrix}^\top N\,u\right) dt.
$$
According to the function definition of lqi in MATLAB its corresponding ARE it does indeed tries to solve the same as your first ARE
$$
\begin{bmatrix}A & 0 \\ -C & 0\end{bmatrix}^\top P + P
\begin{bmatrix}A & 0 \\ -C & 0\end{bmatrix} - \left(P
\begin{bmatrix}B & -D\end{bmatrix} + N\right) R^{-1} \left(
\begin{bmatrix}B & -D\end{bmatrix}^\top P + N^\top\right) + Q = 0.
$$
After some testing I found that this method seems to be identical to mine when you choose $Q$ for lqi the same as your $\mathcal{Q}$, so placing $Q_x$ and $Q_e$ on its diagonal, $R = Q_u$ and $N = 0$. This was somewhat unexpected since I used those same matrices to minimize $\dot{x}$, $\dot{e}$ and $\dot{u}$ instead of $x$, $e$ and $u$. The advantage of this method over mine would be that it does not seem to run into the same troubles as mine did and the corresponding ARE it tries to solve is of smaller dimension which should also make it a little faster to solve.
Best Answer
The Hamiltonian is convex in $p$ because it is affine (linear + constant) with respect this variable. The Hamiltonian is always affine with respect to the dual variables, regardless of the problem at hand.
Edit 1:
The supremum w.r.t. $u$ is inconsequential on the relationship w.r.t. $p$. Saying "the supremum w.r.t. $u$" is equivalent of saying "for fixed $u$".
Edit 2:
You don't need to take the supremum. You can igrone the supremum. Just take the partial derivative w.r.t. $p$. Is this derivative constant with respect to $p$? If so, then the function is affine with respect to $p$. And if a function is affine with respect to a variable, then it is convex with respect to that variable. It is a direct implication.
In the case of the Hamiltonian, the partial derivative with respect to the costates is always the right hand side of the equations of motion which are independent of the costates. Therefore the Hamiltonian is affine with respect to the costates. And therefore it is convex with respect to the costates. Always, no exception.