How to show linear combination of independent, but non-identically distributed Bernoulli's is asymptotically normal? Summary
I am curious about whether there exists theoretical justification to say a linear combination of a sufficiently large number of independent (but not identically distributed) Bernoulli random variables is asymptotically normally distributed.
I know that the sum of i.i.d. Bernoulli's (i.e.,a binomial) approaches a normal distribution, but I have two particular issues which put this in question: (1) they are unequally weighted, so it is a linear combination, and (2) the probability of success is different in most of them.
In other words, consider the random variable $Y$:
$$
Y=a_1X_1 + a_2X_2 + \dots + a_{25000}X_{25000},
$$
where $$X_i\sim \text{Bernoulli}(p_i),\: i\in \{1,2,\dots,25000\}$$
and
$$
p_i\in[0,1]\text{ and }a_i\in\mathbb{R}\:\forall i
$$
Is this $Y$ asymptotically normal?
I have run a basic test in R to evaluate this:
data <- c()
for (j in 1:500) {
  sum <- 0
  for (i in 1:500) {
    sum <- sum + runif(1,0,200)*rbinom(1, 1, runif(1))
  }
  data[j] <- sum
}
shapiro.test(data)
hist(data)

The histogram sure looks normal, and it passes the Shapiro-Wilk test, but if I am to make this assumption, I would need more explicit theoretical justification.
Background
I am evaluating the predicted expense for paying a 3rd party company to run transactions for us. I am using a set of about 25,000 payments made throughout the past two years, but we don't know the payment method they would have used under this system (card or check), and our expense depends greatly on which method they use. We know, using a different set of transactions which are already done using this 3rd party, that larger transactions are more likely to be done by check, and I have taken specific probabilities from a logistic regression, which I am imposing on this other set of payments.
This means that we have about 25,000 bernoulli random variables, the sum of which is our total expense, but our processing fee is a percentage amount of the transactions size when it is a card, so the variables are weighted unequally, and the large ones are more and more likely to be a check, which results in a flat $1 fee.

You need to impose constraints on the weights/probabilities
Without imposing some constraints on the weights, this quantity will not necessarily be asymptotically normal.  It is easy to see this by imagining that the weights increase in magnitude fast enough that a finite number of values will always "dominate" the sum.  If you want to establish a central limit theorem for this kind of quantity, it will be necessary to impose a limiting condition on the weights so that the influence of any individual value (or finite set of values) diminishes to zero in the limit.  This can be done using the Lindeberg or Lyapunov conditions.

Setup: Suppose we have two sequences $\mathbf{a} = (a_1, a_2, a_3, ...)$ and $\mathbf{p} = (p_1, p_2, p_3, ...)$ giving the weights and probabilities for your problem.  Take independent values $X_i \sim \text{Bern}(p_i)$ and define the linear combination:
$$Y_n \equiv \sum_{i=1}^n a_i X_i.$$
The mean and standard deviation of this quantity are given by:
$$\mu_n \equiv \mathbb{E}(Y_n) = \sum_{i=1}^n a_i p_i
\quad \quad \quad \quad \quad 
\sigma_n \equiv \mathbb{S}(Y_n) = \sqrt{\sum_{i=1}^n a_i^2 p_i(1-p_i)}.$$
Consequently, the standardised version of this quantity is:
$$Z_n \equiv \frac{Y_n-\mu_n}{\sigma_n} = \frac{\sum_{i=1}^n a_i (X_i - p_i)}{\sqrt{\sum_{i=1}^n a_i^2 p_i(1-p_i)}}.$$
We can see from the above that $Y_n$ is a sum of independent random variables with finite variance.  This allows us to appeal to versions of the CLT that apply to these kinds of sequences, such as the Lyapunov CLT or the Lindeberg-Feller CLT.  In the next section we show a sufficient condition to fulfil the Lyapunov condition.

Sufficient condition for CLT: Suppose we are willing to assume that as $n \rightarrow \infty$ we have:
$$\frac{\sum_{i=1}^n a_i^{2+\delta} p_i (1 - p_i)}{(\sum_{i=1}^n a_i^2 p_i (1 - p_i))^{1+\delta/2}} \rightarrow 0
\quad \quad \quad \text{for some } \delta > 0.$$
We also have:
$$\begin{align}
\mathbb{E}(|X_i - p_i|^{2+\delta})
&= (1-p_i) |- p_i|^{2+\delta} + p_i |1 - p_i|^{2+\delta} \\[6pt]
&= (1-p_i) p_i^{2+\delta} + p_i (1 - p_i)^{2+\delta} \\[6pt]
&= p_i (1-p_i) [p_i^{1+\delta} + (1 - p_i)^{1+\delta}] \\[6pt]
&\leqslant p_i (1 - p_i). \\[6pt]
\end{align}$$
We therefore have:
$$\begin{align}
\frac{1}{\sigma_n^{2+\delta}} \sum_{i=1}^n \mathbb{E}(|a_i X_i - a_i p_i|^{2+\delta})
&= \frac{1}{\sigma_n^{2+\delta}} \sum_{i=1}^n a_i^{2+\delta} \mathbb{E}(|X_i - p_i|^{2+\delta}) \\[6pt]
&= \frac{1}{\sigma_n^{2+\delta}} \sum_{i=1}^n a_i^{2+\delta} [(1-p_i) p_i^{2+\delta} + p_i (1 - p_i)^{2+\delta}] \\[6pt]
&= \frac{1}{\sigma_n^{2+\delta}} \sum_{i=1}^n (1-p_i) (a_i p_i)^{2+\delta} + \frac{1}{\sigma_n^{2+\delta}} \sum_{i=1}^n p_i (a_i(1-p_i))^{2+\delta} \\[6pt]
&\leqslant \frac{1}{\sigma_n^{2+\delta}} \sum_{i=1}^n a_i^{2+\delta} p_i (1 - p_i) \\[6pt]
&= \frac{\sum_{i=1}^n a_i^{2+\delta} p_i (1 - p_i)}{(\sum_{i=1}^n a_i^2 p_i (1 - p_i))^{1+\delta/2}} \\[12pt]
&\rightarrow 0,
\end{align}$$
which is the Lyapunov condition.  This is sufficient to ensure that the CLT holds so we then have:
$$\lim_{n \rightarrow \infty} \mathbb{P}(Z_n \leqslant z) = \Phi(z),$$
which gives the large sample approximation:
$$Y_n \overset{\text{approx}}{\sim} \text{N} \bigg( \sum_{i=1}^n a_i p_i, \sum_{i=1}^n a_i^2 p_i(1-p_i) \bigg).$$
