A question regarding formulating a PDE for a minimization problem. Question : I'd like to formulate a pde for the following minimization problem.

Let $\Omega$ be a convex, closed, compact set in $\mathbb{R}^d$ with a smooth boundary.
Given a data $(x_i,d_i)$, $x_i \in \Omega^{\mathrm{o}} $ ,$d_i \in \mathbb{R}$, $i = 1,2,3...N$, $N>d$ and $\sum\limits_{i=1}^N d_i = 0$. Also given that, there are always $d$ vectors in $\{x_i\}$ which are linearly independent.
Let $A = \int_{\Omega}dx$
I want to find a continuous function $f:\Omega \to \mathbb{R}$, such that,


*

*$\int_{\Omega}f(x)dx = 0$ and

*$C(f)$ is minimum, where $$C(f) = \frac{A^{\frac{1}{d+1}}}{N}\left\{\sum\limits_{i=1}^N |f(x_i)-d_i|^{d+1}\right\}^{\frac{1}{d+1}} +\|f\|_{L^{d+1}}+ A^{\frac{1}{d+1}} \||\nabla f|\|_{L^{d+1}}$$


The minimum exists, and is unique and is atleast Holder continuous with $\alpha = \frac{1}{d+1}$, due to Sobolev embedding theorem and Morrey's inequality
Reference :  This Q&A from math.stackexchange.
Motivation: I am interested in a stronger result, that is the minimum $f_{min}$ is atleast Lipschitz. We havent leveraged the fact that it is a minimum, while arriving at the weaker result that it is Holder continuous.
 A: Update : My attempt at solution and a problem arising with that.
If I drop the exponents $\frac{1}{d+1}$ and take it outside of one big integral of all three terms, and derive Euler-Lagrange equation I get the following.
$$\frac{A}{N^{d+1}}\sum\limits_i\{ \mathrm sgn(f\delta(x-x_i) - d_i\delta(x-x_i)) |f\delta(x-x_i) - d_i\delta(x-x_i)|^d\delta(x-x_i)\} +  \mathrm sgn(f)|f|^d - A\nabla.(|\nabla f|^{d-1} \nabla f)  = 0                                 $$
Problem with this equation is, the first term contains a product of delta functions ($\delta(x-x_i)\delta(x-x_i)$), the second term is continuous as the solution is proven to be holder continuous, and the third terms contains derivatives of $f$, while RHS is constant and zero. In the midst of these things, I wonder how this equation can have continuous solution, leave alone Lipschitz.
So I wonder, are distributions even allowed in Euler Lagrange equations?  
A: Okay, the main problem is indeed your first term and the Dirac-delta it produces. (There also seems to be a mistake in your derivation, as the product of two deltas doesn't make any mathematical sense)
The key to understanding the EL-equation is calculating the first variation, that is $\frac{d}{ds}|_{s=0} C(f+s\phi)$ for any $\phi$ with $\int_\Omega \phi dx = 0$ (This is so that $f+s\phi$ still is admissible). For a minimizer this has to be zero.
I will do so for each term independently, starting with the easiest. (I will also use the variant of the problem without the exponents $\frac{1}{d+1}$, the basic idea is the same, but the terms will look far uglier)
$$\frac{d}{ds}|_{s=0} \| f+s\phi\|_{L^{d+1}} = \frac{d}{ds}|_{s=0} \int_\Omega |f+s\phi|^{d+1} \stackrel{*}{=} \int_\omega \frac{d}{ds}|_{s=0} |f+s\phi|^{d+1}$$
$$= \int \phi d (f+s\phi) |f+s\phi|^{d-1}|_{s=0} = \int_\Omega \phi d f |f|^{d-1}$$
where * can be justified by the dominated convergence theorem. 
In the same way we have (skipping the pre-factor for now, as all of this is linear)
$$\frac{d}{ds}|_{s=0} \| \nabla(f+s\phi)\|_{L^{d+1}} = \frac{d}{ds}|_{s=0} \int_\Omega |\nabla f+s \nabla \phi|^{d+1} = \int_\Omega \nabla \phi d \nabla f |\nabla f|^{d-1}.$$
For the first term we can try the same (again skipping factors and the sum)
$$\frac{d}{ds}|_{s=0} |f(x_i)+s\phi(x_i)-d_i|^{d+1} =  d \phi(x_i) f(x_i) |f(x_i)|^{d-1} $$
Putting everything together, you want $f$ to satisfy (dividing by the common factor d)
$$0=\frac{A}{N^{d+1}} \sum_i \phi(x_i) f(x_i) |f(x_i)|^{d-1} + \int_\Omega \phi  f |f|^{d-1} +A\int_\Omega \nabla \phi d \nabla f |\nabla f|^{d-1}$$ 
for all sufficiently smooth $\phi$ with mean zero. This is the so called weak formulation of you PDE. To turn in into an actual PDE, you need to rewrite it into the form "something times $\phi$" or in other words find the corresponding distribution. For the first term note that $\phi(x_i) = \delta_{x_i} (\phi)$. For the last term, we can try integration by parts (or rather Green's identities). Here things get a bit complicated, since this produces boundary terms. Those vanish if $\phi$ vanishes at the boundary (which it would have, if we had fixed boundary conditions), however instead the mean vanishes. Still if we disregard that then we end up with
$$0=\frac{A}{N^{d+1}} \sum_i f(x_i) |f(x_i)|^{d-1} \delta_{x_i} + f |f|^{d-1} -A \nabla (d \nabla f |\nabla f|^{d-1})$$
as a sort of Euler-Lagrange equation. Which is however not what I would call a nice equation.
There is however a good workaround, especially since you say that your goal is the regularity of solutions. Whenever $f$ minimises your original problem, f also is a minimum of the simpler problem:
$$ \text{minimise } D(f) := \|u\|_{L^{d+1}}^{d+1} + A^{\frac{1}{d+1}}\| \nabla u \|_{L^{d+1}}^{d+1}$$
with the added constraint that $u(x_i) = c_i$. To get back to the original problem you could then try to optimise the $c_i$, which is a finite dimensional problem and thus somewhat easier. Also, if you understand the structure of solutions to this problem, you understand your original problem, as both would need to have the same regularity.
By the way, as a tangent, note that this problem is rather close to the notion of $p$-capacity (for p=d+1), from the calculus of variations, where $cap_p(C)$ is defined by minimising $\|\nabla u\|_{L^{p}}$, with $u|_{\partial \Omega} = 0$ and $u|_C \geq 1$. So there might be some useful results from that direction as well, although capacity is mostly studied for p < d.
Now finally to get to the example I mentioned in my comment: Assume $\Omega = B_1(0)$, $x_1 =0$, $c_1=1$ (and in a slight disregard of the other conditions, no other $x_i$). Then by symmetry of the problem, the solution should also be symmetric, so we can try the ansatz $u(x) = v(|x|)$. Then we have (I'm skipping the details here, as this answer is getting longer and longer)
$$D(u) = \int_0^1 r^{d-1} |v(r)|^{d+1} + A \int_0^1 r^{d-1} |v'(r)|^{d+1} $$
which has the (1d) EL-equation
$$0= r^{d-1} v(r)|v(r)|^{d-1} - A\frac{d}{dr} \left( r^{d-1} v'(r) |v'(r)|^{d-1}\right)$$
which to me at least looks like it can be solved explicitly. I would assume that solutions of your general problem will look similar around each $x_i$.
