I know quite a lot of Variational principles (VP) yielding systems of classical mechanics. By a VP, I mean something like
$$\delta{\cal L}[U]=0$$
where ${\cal L}$ is a functional and the field belongs to some suitable space and its variations as well. For instance, in $n$-dimensional gas dynamics, the Euler system (conservation of mass, momentum and energy) is associated with the Lagrangian
$${\cal L}(\alpha,s)=\int\int\left(\frac{|m|^2}{2\rho}+f(\rho,s)\right)dxdt,$$
where $\alpha$ is a closed $n$-form whose coordinates are $(\rho,m)$ (the mass density and linear momentum) while $s$ is a function, that is a $0$-form (the entropy). The VP writes
$$\left.\frac{d}{d\epsilon}\right|_{\epsilon=0}{\cal L}[\phi_\epsilon^*(\alpha,s)]=0$$
for every flow of smooth, compactly supported vector field. Notice that the conservation of mass is encoded in the closedness of $\alpha$ : $\partial_t\rho+{\rm div}_xm=0$. See D.S. Sur le principe variationnel des équations de la mécanique des fluides parfaits. *M2AN - Modél. Math. et Anal. Num.*, **27** (1993), pp 739-758.

My question is whether such a VP is available for the relativistic version of the Euler system, where now we do not distinguish between conservation of mass and energy. This system writes \begin{eqnarray*} \partial_t\left(\frac{\rho c^2+p}{c^2-|v|^2}\,-\,\frac p{c^2}\right)+{\rm div}_x\left(\frac{\rho c^2+p}{c^2-|v|^2}\,v\right) & = & 0, \\ \partial_t\left(\frac{\rho c^2+p}{c^2-|v|^2}\,v\right)+{\rm Div}_x\left(\frac{\rho c^2+p}{c^2-|v|^2}\,v\otimes v\right)+\nabla_x p & = & 0. \end{eqnarray*}

Feel free to assume $p=a^2\rho$ for some positive constant $a$.

**Added after Igor's answer.** I am particularly interested in formulations that avoid Lagrange multipliers.