Here are two methods you can use.
Numerical derivative
As @gunes explains, you can estimate the gradient of $F$ using numerical differentiation methods, such as finite differences. In particular, the partial derivative $\frac{\partial F}{\partial x_i}(x)$ can be approximated as
$$\frac{\partial F}{\partial x_i}(x) \approx \frac{F(x + \delta e_i) - F(x - \delta e_i)}{2\delta},$$
where $\delta>0$ is sufficiently small and $e_i=(0,\dots,0,1,0,\dots,0)$ is a vector that is zeros in all coordinates except in the $i$th dimension. If $F$ is a function of $n$ variables, doing this $n$ times will give you an approximation of the gradient of $F$ at $x$. This then lets you compute the gradient of the loss function $L$.
Now you can use gradient descent to train the neural network $NN$. You will use gradient descent to minimize the loss function. This requires knowledge of the gradient of $F$, which you can compute using the method above.
How well this will work is likely to depend on how "locally smooth" $F$ is.
Approximation with a neural network
A different option is to train a neural network $\tilde{F}$ to approximate the function $F$. Then, $\tilde{F}$ will be differentiable, so you can minimize the loss
$$\tilde{L}_{\pi} = || Y_{true} - \tilde{F}(X_{true}, NN_{\pi}(X_{true}) ) ||^2 $$
directly using gradient descent.
How do you train $\tilde{F}$? Well, you pick a random input $x$ and some plausible $\pi$, compute $y=F(x, NN_{\pi}(x))$, and then add $(w,y)$ to the training set for $\tilde{F}$, where $w=(x, NN_{\pi}(x))$. Repeat many times until you have a sufficiently large training set for $\tilde{F}$, then use standard methods for training neural nets to train $\tilde{F}$ on this training set.
A tricky bit is that the there is a circular dependency here: the training set for $\tilde{F}$ depends on $\pi$, but $\pi$ is obtained by training $NN$ on a training set and you have to have $\tilde{F}$ to do that. So, it might be better to use optimization to jointly optimize both $\pi$ and $\tilde{F}$. In particular, we might construct the following loss function:
$$L^* = \| Y - \tilde{F}_{\rho}(X, NN_{\pi}(X) )\|^2 + \lambda \cdot \| \tilde{F}_{\rho}(X, NN_{\pi}(X)) - F(X, NN_{\pi}(X))) \|^2,$$
where here $(X,Y)$ are chosen randomly from the training distribution (i.e., the training set). Then we can use gradient descent to simultaneously find parameters $\pi,\rho$
that minimize $\mathbb{E}[L^*]$. Note that, for a fixed $X,Y$, you can compute the gradient of $L^*$ with respect to $\rho,\pi$ without difficulty.
Here $\lambda>0$ is some hyperparameter that you can set using cross-validation.
How well this will work likely depends on how well-behaved $F$ is and how easy it is to approximate, given some input-output samples for it.