Denote $s_n=(q^n)_\infty=\prod_{k\geqslant n}(1-q^k)$. Note that $$\sum_{k\geqslant n}q^ks_{k+1}=\sum_{k\geqslant n}(s_{k+1}-s_k)=1-s_n.\quad\quad(\clubsuit)$$ The idea is to reduce your sum to such telescoping sums. For $k=0,1,\ldots,n$ denote
$$
c_k:=(-1)^{k+1}{n\choose k}_qq^{k(k+1)/2}(q^{n+1})_{n-k},
$$
where ${n\choose k}_q=\frac{(q)_n}{(q)_k(q)_{n-k}}$ is a $q$-binomial coefficient. I claim that
$$
(-1)^{n+1}q^{(n+1)(3n+2)/2}x(qx)_n\\=c_0q^{2n+1}x+c_1q^{2n}x(q^{2n+1}x)_1+c_2q^{2n-1}x(q^{2n}x)_2+\ldots+c_{n}q^{n+1}x(q^{n+2}x)_n. \quad \quad (\heartsuit)$$ 

First of all, I show how $(\heartsuit)$ yields your formula. Applying $(\heartsuit)$ for $x=q^j$ and using $(\clubsuit)$, your sum reads as
$$
S=\sum_{i=0}^{n}c_i(1-s_{2n+1-i}).
$$
I claim that $$-\sum_{i=0}^n c_is_{2n+1-i}=s_1=\sum_{k\in \mathbb{Z}} (-1)^kq^{k(3k-1)/2}\quad\quad(\spadesuit)$$ (the last equality is Euler's Pentagonal theorem). Thus $S-s_1=\sum_{i=0}^n c_i$, which is easy to see to be a polynomial in $q$ of degree $n(3n+1)/2$. On the other hand, $S$ is divisible by $q^{(n+1)(3n+2)/2}$ (that is seen from its very definition), thus this polynomial must be equal to $-s_1 \pmod {q^{(n+1)(3n+2)/2}}=-\sum_{|k|\leqslant n} (-1)^kq^{k(3k-1)/2}$, and we get your formula for $S$.

It remains to prove $(\heartsuit)$ and $(\spadesuit)$. I start with $(\spadesuit)$. Since $s_{2n+1-i}(q^{n+1})_{n-i}=s_{n+1}$, $(\spadesuit)$ reads as 
$$
s_{n+1}\sum_{i=0}^n (-1)^{i+1}{n\choose i}_q q^{i(i+1)/2}=-s_1,
$$
or, equivalently,
$$
\sum_{i=0}^n (-1)^{i}{n\choose i}_q q^{i(i+1)/2}=(q)_{n},
$$
which is a partial case of the $q$-binomial theorem $(x)_n=\sum_{i=0}^n (-1)^{i}{n\choose i}_q q^{i(i-1)/2}x^i$ for $x=q$.

Well, it remains to prove $(\heartsuit)$. For this we rewrite $(q^{2n+2-k}x)_k$ as $(-1)^kx^kq^{k(4n+3-k)/2}(q^{-2n-1}x^{-1})_k$ and get some variant of $q$-Vandermonde convolution identity.