Euler-Lagrange with derivative of multipliers Let $\Omega\subset \mathbb{R}^2$ be some simply connected domain and consider the functional
\begin{align}
V\left[u(x,y),v(x,y),w(x,y)\right]=\int_\Omega dxdy\:\left(u_x v_y-u_y v_x\right)w,
\end{align}
We wish to find the functions that maximize $V$ constrained to
\begin{align}
(u_x)^2+(v_x)^2+(w_x)^2&=1\\
(u_y)^2+(v_y)^2+(w_y)^2&=1.
\end{align}
If we introduce the multipliers $\lambda_1$ and $\lambda_2$ for the constraints, we have the Lagrangian
\begin{align}
\mathcal{L}=\left(u_x v_y-u_y v_x\right)w-\frac{\lambda_1}{2}\Big[(u_x)^2+(v_x)^2+(w_x)^2-1 \Big]-\frac{\lambda_2}{2}\Big[(u_y)^2+(v_y)^2+(w_y)^2-1 \Big]
\end{align}
that leads to the equations
\begin{align}
0&=u_x v_y-u_y v_x+\left(\lambda_1 w_x \right)_x+\left(\lambda_2 w_y \right)_y\\
0&=-\left(wv_y-\lambda_1 u_x\right)_x+\left(wv_x+\lambda_2 u_y\right)_y\\
0&=\left(wu_y+\lambda_1 v_x\right)_x-\left(w u_x-\lambda_2 v_y\right)_y
\end{align}
The question is, how to solve for the multipliers $\lambda_1$ and $\lambda_2$ ?
In the typical classroom examples the constraints do not depend on the derivatives of the dependent variables, thus the derivatives of the multipliers don't appear in the equations and one can solve for them algebraically. Here, however, that's not the case.
(We impose just the boundary condition $w\left(\partial \Omega\right)=0$, leaving $u$ and $v$ free. This should not be relevant for the question, which is only about how to consistently eliminate the multipliers.)
 A: 5 coupled non-linear PDEs seem a too difficult task to solve. In this answer, we will only discuss the boundary conditions for the Lagrange multipliers $\lambda_1$ and $\lambda_2$.

*

*The De Donder/Weyl polymomenta are
$$ {\bf p}_u~=~\begin{pmatrix} v_y w -\lambda_1 u_x \cr 
-v_x w -\lambda_2 u_y \end{pmatrix}, $$
$$ {\bf p}_v~=~\begin{pmatrix} -u_y w -\lambda_1 v_x \cr 
u_x w -\lambda_2 v_y \end{pmatrix}, $$
$$ {\bf p}_w~=~\begin{pmatrix} -\lambda_1 w_x \cr 
 -\lambda_2 w_y \end{pmatrix}, $$
$$ {\bf p}_{\lambda_1}~=~\begin{pmatrix} 0 \cr 
 0 \end{pmatrix}, $$
$$ {\bf p}_{\lambda_2}~=~\begin{pmatrix} 0 \cr 
 0 \end{pmatrix}. $$


*OP lists the following boundary conditions:
$$ u,v \text{ free and }  \left.w \right|_{\partial \Omega}~=~0. $$
We assume that this also means that the derivatives of $u,v$ are free.


*The infinitesimal variation of OP's ${\cal L}$ reads
$$\delta {\cal L} ~=~\sum_{\phi=1}^5\left\{ \underbrace{\left(\frac{\partial {\cal L}}{\partial\phi} -\nabla\cdot{\bf p}_{\phi} \right)}_{\text{EL expression}} 
\delta \phi+ \nabla\cdot\left({\bf p}_{\phi}\delta \phi \right)\right\},$$
where $\phi$ runs over the 5 variables $u,v,w,\lambda_1,\lambda_2$.


*In order for the variational principle for OP's ${\cal L}$ to be well-posed, we must demand that $$  {\bf p}_u \text{ and } {\bf p}_v \text{ are parallel to the boundary  at the boundary } \partial \Omega,$$
cf. the divergence theorem.


*Generically this implies that
$$ \left.\lambda_1\right|_{\partial \Omega}~=~0 \quad\text{ and } \quad \left.\lambda_2\right|_{\partial \Omega}~=~0, $$
which in turn are also sufficient conditions for well-posedness.
