You can rephrase your question as follows: first we subtract the known vector from both and then take care of the known coordinates. So assuming the coordinates of the two points are $(\alpha,\beta)$ and $(\gamma,X)$ where $\alpha,\gamma \in \mathbb{R}^m$ are known, and $\beta \in \mathbb{R}^n$ is known, but $X$ represent the unknown coordinates constrainted to lie inside the cube $[-1,1]^n$, the integral to evaluate becomes
$$\frac{1}{2^n}\int_{[-1,1]^n} \sqrt{ |\alpha-\gamma|^2 + |\beta - X|^2 } dX$$
In the lower dimensional case this can be integrated. But an analytical expression in higher dimensions is elusive. For the case $\alpha = \gamma$ and $\beta = 0$, some bounds were obtained in an old paper of Anderssen et al. http://dx.doi.org/10.1137/0130003 For more general probability distributions there is a recent paper with some bounds by Burgstaller and Pillichshammer. http://journals.cambridge.org/action/displayAbstract?aid=6622208
Of course, one can get a fairly trivial bound by Cauchy-Schwartz
$$ \int_{[-1,1]^n} f(X) dX \leq 2^{n/2} \left( \int_{[-1,1]^n} f(X)^2 dX \right)^{1/2} $$
and that
$$ \int_{[-1,1]^n} R^2 + |\beta - X|^2 dX = 2^n (R^2 + \beta^2) + \int_{[-1,1]^n} X^2 dX$$
the last term is simply evaluated as $n 2^n / 3$, so putting it all together we have the upper bound for the expected value by
$$ \frac{1}{2^n}\int_{[-1,1]^n} \sqrt{ |\alpha-\gamma|^2 + |\beta - X|^2 } dX \leq \sqrt{ |\alpha -\gamma|^2 + \beta^2 + \frac{n}{3}}$$
which is slight improvement over the utterly trivial upper/lower bound of $\sqrt{|\alpha-\gamma|^2 + \beta^2 \pm n}$ if you just maximize/minimize each coordinates.