Bien, el problema principal es efectivamente su primer término y el delta de Dirac que produce. (También parece haber un error en tu derivación, ya que el producto de dos deltas no tiene ningún sentido matemático)
La clave para entender la ecuación EL es calcular la primera variación, es decir $\frac{d}{ds}|_{s=0} C(f+s\phi)$ para cualquier $\phi$ con $\int_\Omega \phi dx = 0$ (Esto es para que $f+s\phi$ sigue siendo admisible). Para un minimizador esto tiene que ser cero.
Lo haré para cada término de forma independiente, empezando por el más fácil. (También utilizaré la variante del problema sin los exponentes $\frac{1}{d+1}$ (la idea básica es la misma, pero los términos serán mucho más feos)
$$\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}$$
donde * puede justificarse por el teorema de convergencia dominada.
De la misma manera tenemos (omitiendo el factor previo por ahora, ya que todo esto es lineal)
$$\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}.$$
Para el primer término podemos intentar lo mismo (de nuevo omitiendo los factores y la suma)
$$\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} $$
Si se combina todo, se quiere $f$ para satisfacer (dividiendo por el factor común 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}$$ para todos los casos suficientemente suaves $\phi$ con media cero. Esta es la llamada formulación débil de la EDP. Para convertirla en una EDP real, hay que reescribirla en la forma "algo por $\phi$ "o, en otras palabras, encontrar la distribución correspondiente. Para el primer término hay que tener en cuenta que $\phi(x_i) = \delta_{x_i} (\phi)$ . Para el último término, podemos probar la integración por partes (o más bien las identidades de Green). Aquí las cosas se complican un poco, ya que esto produce términos de frontera. Estos desaparecen si $\phi$ desaparece en la frontera (lo que ocurriría si tuviéramos condiciones de frontera fijas), pero en cambio la media desaparece. Aun así, si no tenemos en cuenta esto, acabamos con
$$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})$$ como una especie de ecuación de Euler-Lagrange. Que sin embargo no es lo que yo llamaría una ecuación agradable.
Sin embargo, hay una buena solución, sobre todo porque dices que tu objetivo es la regularidad de las soluciones. Siempre que $f$ minimiza su problema original, f también es un mínimo del problema más simple:
$$ \text{minimise } D(f) := \|u\|_{L^{d+1}}^{d+1} + A^{\frac{1}{d+1}}\| \nabla u \|_{L^{d+1}}^{d+1}$$ con la restricción añadida de que $u(x_i) = c_i$ . Para volver al problema original se podría intentar optimizar el $c_i$ que es un problema de dimensión finita y, por tanto, algo más fácil. Además, si entiendes la estructura de las soluciones de este problema, entiendes tu problema original, ya que ambos tendrían que tener la misma regularidad.
Por cierto, como tangente, nótese que este problema está bastante cerca de la noción de $p$ -(para p=d+1), del cálculo de variaciones, donde $cap_p(C)$ se define minimizando $\|\nabla u\|_{L^{p}}$ con $u|_{\partial \Omega} = 0$ y $u|_C \geq 1$ . Así que podría haber algunos resultados útiles desde esa dirección también, aunque la capacidad se estudia principalmente para p < d.
Ahora, finalmente, para llegar al ejemplo que mencioné en mi comentario: Supongamos que $\Omega = B_1(0)$ , $x_1 =0$ , $c_1=1$ (y en un ligero descuido de las otras condiciones, no hay otra $x_i$ ). Entonces, por simetría del problema, la solución debe ser también simétrica, por lo que podemos probar el ansatz $u(x) = v(|x|)$ . Entonces tenemos (me salto los detalles aquí, ya que esta respuesta es cada vez más larga) $$D(u) = \int_0^1 r^{d-1} |v(r)|^{d+1} + A \int_0^1 r^{d-1} |v'(r)|^{d+1} $$ que tiene la ecuación EL (1d) $$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)$$ que al menos a mí me parece que se puede resolver explícitamente. Yo asumiría que las soluciones de su problema general serán similares alrededor de cada $x_i$ .
