I am working on an optimization problem which I am stuck on towards the end.

Essentially, I have two probability density functions in $\mathbb{R}^2$, call them $q(x,y)$ and $p(x,y)$, now I define the objective functional to be:

$C(p,q) = \int p(x,y) \ln\left[\frac{p(x,y)}{q(x,y)}\right] dxdy \cdots (1)$

Now assume we already know what $q(x,y)$ is, ie, it is pre-determined, so my aim is to minimise $(1)$ with respect to the unknown density $p(x,y)$ subject to the following constraints:

$\int_{-\infty}^{\infty} \int_{X^x_d}^{\infty} p(x,y) dx dy = PoD^x$

$\int_{-\infty}^{\infty} \int_{X^y_d}^{\infty} p(x,y) dy dx = PoD^y$

$\int_{-\infty}^{\infty} \int_{-\infty}^{\infty} p(x,y) dx dy = 1$

$p(x,y) \ge 0 \ \forall (x,y) \in \mathbb{R}^2$

Now if we define:

$ I_{[X^x_d, \infty)} = \begin{cases} 1 & \text{if} x \ge X^x_d \\ 0 & \text{if} x < X^x_d \end{cases}$

$ I_{[X^y_d, \infty)} = \begin{cases} 1 & \text{if} y \ge X^y_d \\ 0 & \text{if} y < X^y_d \end{cases}$

Then our constraints become:

$\int \int p(x,y) I_{[X^x_d, \infty)} dx dy = PoD^x \cdots (2)$

$\int \int p(x,y) I_{[X^y_d, \infty)} dy dx = PoD^y \cdots (3)$

$\int \int p(x,y) dx dy = 1 \cdots (4)$

So I set up the Lagrangian as follows:

$L(p,q) = \int \int p(x,y) \ln[p(x,y)] dxdy - \int \int p(x,y) \ln[q(x,y)] dxdy + \lambda_1 \left[ \int \int p(x,y) I_{[X^x_d, \infty)} dx dy - PoD^x\right] + \lambda_2 \left[\int \int p(x,y) I_{[X^y_d, \infty)} dydx - PoD^y \right] + \mu \left[\int \int p(x,y) dxdy -1\right]$

which can further be simplified:

$L(p,q) = \int \int p(x,y)\left[\ln[p(x,y)] - \ln[q(x,y)]\right] dxdy + \int \int p(x,y) \left[\lambda_1 I_{[X^x_d, \infty)} + \lambda_2 I_{[X^y_d, \infty)} + \mu \right] dxdy - \lambda_1 PoD^x - \lambda_2 PoD^y - \mu$

Now to solve this for $p(x,y)$ I was told to use the calculus of variations and the optimal solution would be:

$\widehat{p(x,y)} = q(x,y) \exp\{-\left[1+\hat{\mu} + \left(\hat{\lambda_1} I_{[X^x_d, \infty)} \right) + \left(\hat{\lambda_2} I_{[X^y_d, \infty)} \right) \right] \}$

Now I am not very familiar with calculus of variations, I only really know how to optimize functionals with respect to algebraic constraints, however I am not sure what to do when the constraints are (multiple) integrals, any help would be appreciated :)