1. The Hohenberg-Kohn Theorems:
I think the first part of your question does indeed summarize well the Hohenberg-Kohn theorems. Apart from some minor details (such as that $V$ is unique up to an additive constant, or $\psi$ is unique up to a phase), I think I have nothing to add. Your explanation should very well agree with the ones given in the books of e.g. Dreizler.
2. $v$-Representability:
As you've pointed out already, the discussion basically takes places for a fixed interaction $U$ (say e.g. Coulomb interaction) and kinetic energy $T$ (and particle number $N$). For a fixed $U$, we are interested in all possible ground states (here we assume non-degeneracy, which is not much of a restriction for the following discussion) generated by the different external potentials $V$. Your set $\mathcal N$ is often called set of $v$-representable densities: It contains all densities which can be realized from the ground state $\psi_0$ for some external potential $V$.
Note that this is a relative notion: Given a density $n$ (what a density, or better a $N$-representable density is can be defined mathematically), it might be $v$-representable for $U$, but not for a different interaction $U^\prime$.
3. Non-interacting $v$-representability:
Of course,the Hohenberg-Kohn theorems also hold for $U=0$, i.e. non-interacting systems. Suppose the system you are interested in has for some external potential $V$ a ground state $\psi_0$ with ground state density $n_0$. This density is thus (interacting-)$v$-representable. It is additionally non-interacting $v$-representable if there is a non-interacting system $(U=0)$ for which there is an external potential $V_\mathrm{eff}$ such that the ground state $\phi_0$ gives rise to $n_0$.
Thus, if you have a density which is both $v$-representable and non-interacting $v$-representable, there are (at least) two systems which have a ground state giving the said density. Whether or not every $v$-representable density is also non-interacting $v$-representable is, as far as I know, an unsettled problem.
4. Kohn-Sham Systems:
Assuming for now that indeed every $v$-representable density is also non-interacting $v$-representable, we can discuss the Kohn-Sham systems. Suppose you knew the effective potential $V_\mathrm{eff}$ which gives rise to the ground state density $n_0$ of interest. Then you could take this potential, solve the non-interacting Schrödinger equation and then simply extract from the ground state $\phi_0$ the ground state density $n_0$, which you then may plug into the (approximated) energy functional.
However, the conditional statement is important. We know that there exists some $V_\mathrm{eff}$ with the desired properties, but we don't know a priori how it looks like. There are ways however, explained in basically any reference I gave above, to construct this potential. What you will find is that $V_\mathrm{eff}$ depends on the functional derivative of the energy functional and in particular it depends on the ground state density.
This then leads to the known fact that you have to solve the Kohn-Sham equations self-consistently. In other words, given some approximation for the energy functional, you may start with an initial guess for the ground state density $\tilde n_0$, solve the Kohn-Sham equations and get a new ground state density. Repeat the cycle until convergence (which might be a topic on its own).
Answer to questions in the comments:
The external potential operator for your system of $N$ (identical) fermions (e.g. electrons) is of the form $V=\sum\limits_{i=1}^N v(x_i)$. The function $v:\mathbb R^3\longrightarrow \mathbb R$ is the external potential, which every electron feels (e.g. the presence of the fixed ions of a lattice in some solid material). The effective potential also takes this form $V_\mathrm{eff} =\sum\limits_{i=1}^N v_{\mathrm {eff}}(x_i)$. The Hamiltonian of the non-interacting system is by definition $$H_{\mathrm{eff}}=T+V_{\mathrm{eff}}= \sum\limits_{i=1}^N h_i \quad ,$$ where $h_i =- \frac{\hbar^2}{2m}\nabla_i^2 + v_{\mathrm {eff}}(x_i)$. Without going into too much detail, we can show (with some assumptions, e.g. non-degeneracy) that the ground state is of the form of a Slater-determinant: If $$h\varphi_j=\epsilon_j \varphi_j \quad ,$$ then an eigenstate of $H_{\mathrm{eff}}$ is given by $\psi\sim \varphi_{j_1} \wedge \varphi_{j_2} \wedge \ldots \wedge \varphi_{j_N}$ with energy $E=\sum\limits_{k=1}^N \epsilon_{j_k}$. This also shows you that the ground state of this system is the one such that its corresponding energy is build from the $N$ lowest $\epsilon_j$'s.
Consequently, in order to solve the non-interacting Schrödinger equation with the Hamiltonian $H_\mathrm{eff}$, you can solve the single-particle Schrödinger equation from which you take the solutions $\varphi_j$ with the $N$ lowest energies $\epsilon_j$ to build your ground state Slater determinant, from which in turn you get the ground state density. But $v_\mathrm{eff}$ depends on the ground state density $n_0$, which we have to determine. This makes the corresponding single-particle Schrödinger equation a non-linear problem and hence we have to solve the Kohn-Sham equations self-consistently.
Best Answer
The sum of the probabilities of all mutually exclusive events must equal one. However, just because you've found an electron somewhere, it does not mean that you can't find another somewhere else (except for $N=1$).
Reference: A Primer in Density Functional Theory by Fiolhais, C., Nogueira, F., & Marques, M. A. (Eds.), Springer, chapter 1.2. Especially the discussion below equation (1.21).