Chance of receiving all elements of a set I have a set of $n$ different elements. I will select $i$ times a subset $S_j$ of $n/k$ elements randomly. Each element can only occur once in each $S_j$, but can be part of multiple different subsets $S_j$. I want to know the probability of having all elements in the sum of the subsets.

My attempt:
This would be a possible outcome of what I want to do:
$n= 4$ | number of different items
$elements = \{ A,B,C,D \}$ | set of items
$k=2$ | nr. of elements per subset
$ i=3$ | nr. of subsets generated
$S_1=\{A,B\}$ | 1st random subset
$S_2=\{C,A\}$ | 2nd random subset
$S_3=\{B,A\}$ | 3rd random subset
$total = S_1 \cup S_2 \cup S_3 = \{A,B,C\}$ | all elements in at least 1 subset
I am looking for the probability $p(n,k,i)$ for the event "Each element from $elements$ is contained in $total$"
 A: Fix a nonnegative integer $n$, an integer $m$ with $0\leq m\leq n$, and a nonnegative integer $l$.  Let $q(n,m,l)$ be the probability that the union of random $m$-subsets $S_r$, for $r\in[l]$ with $[l]:=\{1,2,\ldots,l\}$, of the set $[n]:=\{1,2,\ldots,n\}$ is precisely $[n]$, where the random sets are chosen uniformly randomly and statistically independently.  Observe that, if $k\in\mathbb{N}$ is a divisor of $n$ and $i$ is a positive integer, then $$p(n,k,i)=q\left(n,\frac{n}{k},i\right)\,.$$  The event $E_{A,r}$ denotes the event that $A\cap S_r=\emptyset$, for $A\subseteq [n]$ and $r\in[l]$.  The notation $\mathbb{P}$ denotes the probability measure.  Also, for $A\subseteq [n]$, the union $\bigcup\limits_{r=1}^l\,E_{A,r}$ will be written as $E_A$.  Observe that $E_{A}\cap E_{B}=E_{A\cup B}$ for any $A,B\subseteq[n]$.
First, it is evident that $$\mathbb{P}\left(E_{A,r}\right)=\frac{\binom{n-|A|}{m}}{\binom{n}{m}}=\prod_{\mu=1}^m\,\left(1-\frac{|A|}{n-\mu+1}\right)\,.$$
for all $A\subseteq [n]$ and $r\in[l]$.  Thus, by statistical independence,
$$\mathbb{P}\left(E_A\right)=\mathbb{P}\left(\bigcap_{r=1}^l\,E_{A,r}\right)=\prod_{r=1}^l\,\mathbb{P}\left(E_{A,r}\right)=\prod_{\mu=1}^m\,\left(1-\frac{|A|}{n-\mu+1}\right)^l\,.$$
Finally, using the Principle of Inclusion and Exclusion, along with the fact that $\mathbb{P}\left(E_A\right)=\mathbb{P}\left(E_B\right)$ if $A,B\subseteq[n]$ are equicardinal, we obtain
$$q(n,m,l)=1-\mathbb{P}\left(\bigcup_{j=1}^n\,E_{\{j\}}\right)=\sum_{\nu=0}^n\,(-1)^{\nu}\,\binom{n}{\nu}\,\prod_{\mu=1}^m\,\left(1-\frac{\nu}{n-\mu+1}\right)^l\,,$$
under the convention that $0^0=1$. For $l>0$, the formula above is also given by
$$q(n,m,l)=\sum_{\nu=0}^{n-m}\,(-1)^{\nu}\,\binom{n}{\nu}\,\prod_{\mu=1}^m\,\left(1-\frac{\nu}{n-\mu+1}\right)^l\,.$$
I do not think that there is a nice simplification of the expression above.
For examples, $q(n,m,0)=0$ $n>0$ and $q(n,n,l)=1$ for $l>0$, as well as $q(0,0,0)=1$.  If $ml<n$ holds, then $q(n,m,l)=0$.  Also, the answer to the OP's example is $$p(4,2,3)=q(4,2,3)=\frac{19}{36}\,.$$  Anyhow, there is a good approximation when $m\ll n$ (see the hidden portion for a proof):
$$q(n,m,l)\approx\Biggl(1-\left(1-\frac{1}{n}\right)^{ml}\Biggr)^n\,.$$
Below is a comparison between the actual plot (blue circles) and its approximation (red squares) with $n=10$ and $m=2$.


  Here is a probabilistic proof of the approximation.  If $m\ll n$, then we can well approximate the selections $S_1,S_2,\ldots,S_l$ by randomly choosing elements $X_{\sigma,\rho}$ with $\sigma=1,2,\ldots,m$ and $\rho=1,2,\ldots,l$ and setting $S_r:=\left\{X_{1,r},X_{2,r},\ldots,X_{m,r}\right\}$ (as $m\ll n$, the chance that $\left|S_r\right|<m$ is slim).  For each element $j$ of $[n]$, the probability that $j$ is not in the union $\bigcup\limits_{r=1}^l\,S_r$ is roughly $\left(1-\frac{1}{n}\right)^{ml}$ (that is, each of the $ml$ random elements $X_{\sigma,\rho}$ must not be equal to $j$).  Hence, the probability that $j$ is in the union $\bigcup\limits_{r=1}^l\,S_r$ is approximately $1-\left(1-\frac{1}{n}\right)^{ml}$.  Using the statistical independence, we conclude that the probability that $[n]=\bigcup\limits_{r=1}^l\,S_r$ is well approximated by $\Big(1-\left(1-\frac{1}{n}\right)^{ml}\Big)^n$.

A: Remark, added  later. I somehow  missed the fact that  the formula
below  is the  same  as in  the  accepted post.   Maybe the  alternate
derivation  can  facilitate  understanding  of the  problem.  Text  of
original version follows.
Here is a different formula.  Building on the work by @Batominovski we
use the  notation $$q(n, m, l)$$  for the probability  in the scenario
with $n$ elements, subsets of size $m$ and $l$ samples.
Use the variable $A_\nu$ to represent the type $\nu$ from among the
$n$ available  types. The generating  function for the  $m$-subsets is
given by
$$[z^m] \prod_{\nu=1}^n (1 + zA_\nu).$$
It follows that  the generating function for the  $l$ samples is given
by
$$\left([z^m] \prod_{\nu=1}^n (1 + zA_\nu)\right)^l.$$
Now  we  are interested  in  those terms  that  contain  at least  one
instance  of the  $n$  types. Evaluating  the  generating function  at
$A_\nu=0$ and subtracting that value yields the generating function of
the  terms containing  at least  one $A_\nu.$  So we  do this  for all
$A_\nu.$ Note  however that we  have now removed terms  not containing
$A_\nu$ and $A_\mu$ two times, so we must add those back in and so on.
We have by Inclusion-Exclusion for the desired generating function
$$\sum_{S\subseteq A} (-1)^{|S|} 
\left.\left([z^m] 
\prod_{\nu=1}^n (1 + zA_\nu)\right)^l\right|_{S=0.}$$
Actually doing  the substitution  and setting the  elements of  $S$ to
zero produces  $k = |S|$ factors  that are equal to  one. Then setting
the  remaining $A_\nu$ to  one because  we are  not interested  in the
classification but rather the count  produces $n-k$ terms that are all
equal to $1+z.$ This yields for the inner term
$$[z^m] (1+z)^{n-k} = {n-k\choose m}$$
and for the probability
$$\bbox[5px,border:2px solid #00A000]{
{n\choose m}^{-l} \sum_{k=0}^n {n\choose k} (-1)^k
{n-k\choose m}^l.}$$
Sanity check. With just one sample we should get probability zero
for $m\lt n$ and probability one for $n=m.$ Setting $l=1$ we obtain
$${n\choose m}^{-1} \sum_{k=0}^n {n\choose k} (-1)^k [z^m] (1+z)^{n-k}
\\ = {n\choose m}^{-1} 
[z^m] (1+z)^n \sum_{k=0}^n {n\choose k} (-1)^k \frac{1}{(1+z)^k}
\\ = {n\choose m}^{-1} 
[z^m] (1+z)^n \left(1-\frac{1}{1+z}\right)^n
\\ = {n\choose m}^{-1} [z^m] z^n.$$
This is indeed zero when $m\lt n$ and one when $n=m.$
Here is the Maple code for comparison of the two formulae.

Q1 := (n,m,l) ->
add((-1)^nu*binomial(n, nu) *
    mul(1-nu/(n-mu+1),mu=1..m)^l, nu=0..n-m);


Q2 := (n,m,l) ->
binomial(n,m)^(-l) *
add(binomial(n,k)*(-1)^k*binomial(n-k,m)^l, k=0..n);

Addendum. We now compute a closed form for the case $l=2.$
We get for the sum term
$$\sum_{k=0}^n {n\choose k} (-1)^k
{n-k\choose m}^2.$$
Introducing 
$${n-k\choose m} =
\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{m+1}} (1+z)^{n-k} \;dz$$
we obtain for the sum
$$\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{m+1}} 
(1+z)^{n} 
\frac{1}{2\pi i}
\int_{|w|=\gamma} \frac{1}{w^{m+1}} 
(1+w)^{n} 
\\ \times \sum_{k=0}^n {n\choose k} (-1)^k
\frac{1}{(1+z)^k (1+w)^k}
\; dw \;dz
\\ = \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{m+1}} 
(1+z)^{n} 
\frac{1}{2\pi i}
\int_{|w|=\gamma} \frac{1}{w^{m+1}} 
(1+w)^{n} 
\\ \times \left(1-\frac{1}{(1+z)(1+w)}\right)^n
\; dw \;dz 
\\ = \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{m+1}} 
\frac{1}{2\pi i}
\int_{|w|=\gamma} \frac{1}{w^{m+1}} 
(z+w+wz)^n
\; dw \;dz.$$
Extracting the inner coefficient produces
$$\frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{m+1}} 
{n\choose m} (1+z)^m z^{n-m}
\;dz
\\ = {n\choose m}  \frac{1}{2\pi i}
\int_{|z|=\epsilon} \frac{1}{z^{2m-n+1}} 
(1+z)^m
\;dz
\\ = {n\choose m} {m\choose 2m-n}.$$
This yields for the probability
$$q(n, m, 2) = {n\choose m}^{-1} {m\choose n-m}.$$
