Let $\nabla F= \langle \partial_x F, \partial_{x_1} F, \partial_y F, \partial_{y_1} F \rangle$. The method of Lagrange in this case requires the introduction of two multipliers. We should solve:
$ \nabla d = \lambda_1\nabla g+ \lambda_2\nabla h $
subject to constraints 1 and 2 as listed in your post.
The reason for this is as follows: if an extrema exists then curves $t \mapsto \alpha(t) $ which pass through the extremal point must make $\eta=d \circ \alpha$ extreme at the corresponding point in the domain. Suppose $t=0$ gives $\alpha(0)$ the extremal point (we can shift the parameter to make this happen, nothing is lost in this convenience).
The curve $\alpha$ lies on the intersection of the level sets given by 1 and 2. We have
$ g (\alpha(t)) =0 \qquad h( \alpha(t))=0 $
The chain-rule yields
$ \nabla g (\alpha(t)) \cdot \alpha'(t)=0 \qquad \nabla h (\alpha(t)) \cdot \alpha'(t)=0 $
Likewise, since $\eta$ is extremal at $t=0$ the chain rule gives
$ \nabla d (\alpha(0)) \cdot \alpha'(0)=0 $
Summarizing, the tangent vector field $\alpha'$ is orthogonal to the gradient fields of $g$ and $h$ where they can be compared and at $\alpha(0)$ the tangent $\alpha'(0)$ is orthogonal to $\nabla d$. The point $\alpha(0)$ is special in that we obtain orthogonality with respect to $\nabla d, \nabla g$ and $\nabla h$.
At first glance this would not appear to connect $\nabla d, \nabla g$ and $\nabla h$ in any particular way. However, there is not just one curve on the constraint surface. Provided the constraints 1. and 2. are nondegenerate the level set they define is two-dimensional and there will be a two-dimensional plane of tangent vectors which are found orthogonal to $\nabla d, \nabla g$ and $\nabla h$. But, this means that $\nabla d, \nabla g$ and $\nabla h$ are linearly dependent since $\mathbb{R}^4$ should be the direct sum of the tangent and normal space. For these reasons we introduce multipliers to ascertain the location of the max/min solution.
Notice the method is based on the existence of extreme solutions. For the continuous function $d$ these are known to exist if the constraint is a compact surface. Sometimes the method still "works" for non-compact constraints, but beware the limit of the method.