There is no maximum unless $\beta+\lambda\equiv1$ and $\gamma\equiv\lambda$. The terms containing $p_s$ add up to
$\iint_{\mathbb R_2}p_s(1-\beta-\lambda)\mathrm dx\mathrm dy\;,$
so you can make the integral arbitrarily large by choosing $p_s$ arbitrarily large (positive or negative depending on the sign of $1-\beta-\lambda$). Likewise, the terms containing $p_t$ add up to
$\iint_{\mathbb R_2}p_t(\lambda-\gamma)\mathrm dx\mathrm dy\;.$
This fixes $\beta$ and $\gamma$ in terms of $\lambda$, and then $p_s$ and $p_t$ no longer occur in the integral. We can drop the constant terms containing $\beta$ and $\gamma$, leaving us with
$I=\iint_{\mathbb R^2}\left[- \alpha|\vec{p}|+ \lambda(\vec{\nabla}\cdot\vec{p})\right]\mathrm dx\mathrm dy\;.$
Let's restrict ourselves to solutions $\vec{p}$ that decay sufficiently rapidly at infinity that we can integrate by parts and omit the boundary term:
$I=\iint_{\mathbb R^2}\left[- \alpha|\vec{p}|- (\vec{\nabla}\lambda)\cdot\vec{p}\right]\mathrm dx\mathrm dy\;.$
So we want to find $\vec{p}$ such that its magnitude (weighted by $\alpha$) is small but its component along $-\vec{\nabla}\lambda$ is large. This will work best if we choose $\vec{p}$ to always point along $-\vec{\nabla}\lambda$, i.e. $\vec{p}=-\mu\vec{\nabla}\lambda$, with as yet undetermined $\mu(x,y)$. With $\vec{g}:=\vec{\nabla}\lambda$, this gives
$ \begin{eqnarray} I &=& \iint_{\mathbb R^2}\left[- \alpha|-\mu\vec{g}|- \vec{g}\cdot(-\mu\vec{g})\right]\mathrm dx\mathrm dy \\ &=& \iint_{\mathbb R^2}\left[- \alpha|\mu||\vec{g}|+\mu|\vec{g}|^2\right]\mathrm dx\mathrm dy \\ &=& \iint_{\mathbb R^2}|\vec{g}|\left[- \alpha|\mu|+\mu|\vec{g}|\right]\mathrm dx\mathrm dy \;. \end{eqnarray}$
Thus, we must have $\alpha\ge|\vec{\nabla}\lambda|$ everywhere; else we could make the integral arbitrarily large by choosing $\vec{p}$ to vanish wherever that inequality holds and to make the integral arbitrarily large where it doesn't.
Given that inequality, the integral is forced to be non-positive, and thus it is maximized for $\vec{p}\equiv0$.
The first conclusion, $\alpha\ge|\vec{\nabla}\lambda|$, holds independent of the assumption about the behaviour at infinity, since we can choose $\vec{p}$ to have that behaviour if we want. The second conclusion, $\vec{p}\equiv0$, however, holds only under that assumption, since the boundary term might compensate a negative value of the integral. I'm not sure how to treat the general case, without decay at infinity.
[Edit:] Quite generally speaking, the parts of the integral that depend on $\vec{p}$ all scale linearly if you scale $\vec{p}$ by a positive scale factor, so there can't be any maximum other than $\vec{p}=0$ -- all you can do by imposing constraints on $\alpha$, $\beta$, $\gamma$ and $\lambda$ is to ensure that $\vec{p}=0$ is indeed a maximum.