You might be able to use something like the linear combination of unitaries protocol from (Childs, 2017). Summarizing Lemmas 6-7 in that paper, you will need any unitary $V$ that can prepare the following state on $m = \lceil \log_2 N \rceil$ qubits:
$$
V|0^m\rangle = \frac{1}{\sqrt{c}} \sum_{i=1}^N \sqrt{c_i} |i\rangle \tag{1}
$$
where $c = \sum_i c_i$, and you need the ability to implement the unitary
\begin{equation}
U = \sum_{i=1}^N |i\rangle \langle i| \otimes U_i \tag{2}
\end{equation}
which you can think of as a controlled gate generalized to qudits - see answers to this question for more discussion on such gates. Then assuming that your unitaries $U_i$ each act on $n$ qubits you can show that
$$
(V^\dagger \otimes I) U (V \otimes I)|0^m\rangle |0^n\rangle = \frac{1}{c} |0^m\rangle |\Psi\rangle + |\Phi^\perp\rangle \tag{3}
$$
where $|\Psi\rangle$ is the state you're trying to prepare. The $|\Phi^\perp\rangle$ appears when you implement $V^\dagger$ - since we only specify the first column in defining $V$, applying $V^\dagger$ gives us a "garbage" output associated with every other basis state in the first register besides $|0^m\rangle$.
Once you've prepared the state in Eq. $(3)$, you can just measure the first $m$ qubits and with probability $(\lVert |\Psi\rangle\rVert/c)^2 $ you will observe $0^m$, which means you've prepared $|\Psi\rangle$ in the second register. Alternatively, you can gain a square-root speedup by first performing amplitude amplification on the first term of the state given in Eq. $(3)$.