My attempt, although it yields an answer which isn't very nice
Take a particular set $F=\{i:1\le i\le k\}$.
Let us find the number of factors which are coprime to $\prod_{i\in F}p_i$.
For any divisor $d$ and $b_1,...,b_j \in\mathbb N_0,\,\,\, (\forall i)\,\,\,b_i\le a_i$,
$$\gcd\left(d\,\,\,,\,\,\,\prod_{i\in F}p_i\right)=1 \iff d=\prod_{i\not\in F} p_i^{b_i}$$
$\implies$ (set of all possible divisors $\neq1$ which are coprime to the product of primes indexed by $F$)
$$\left|\left\{d\neq1:(d\,\,|\,\,n)\land\gcd\left(d\,\,\,,\,\,\,\prod_{i\in F}p_i\right)=1\right\}\right|=\left(\prod_{i\not\in F} (a_i+1)\right)-1$$
Now, for $b_1,...,b_n \in\mathbb N_0$,
$$\gcd\left(\prod_{i\in J} p_i, x\right)=1 \iff \gcd\left(\prod_{i\in J} p^{b_i}, x\right)=1,$$
$\implies$ (set of all possible pairs with $f$ being composed only of primes indexed by $F$)
$$\left|\left\{(d,f):(d\,\,|\,\,n)\land(b_1,...,b_{|F|} \in\mathbb N_0,\,\,\, (\forall i)\,\,\,b_i\le a_i)\land\left(f=\prod_{i\in F}p_i^{b_i}\right)\land(\gcd\left(d,f\right)=1)\right\}\right|=\left(\left(\prod_{i\in F} (a_i+1)\right)-1\right)\left(\left(\prod_{i\not\in F} (a_i+1)\right)-1\right)$$
$$=\prod^k_{i=1}{(a_i+1)}-\prod_{i\in F}{(a_i+1)}-\prod_{i\not\in F}{(a_i+1)}+1$$
Now, we are interested in $\left|\left\{\{d,f\}:(d\,\,|\,\,n)\land(f\,\,|\,\,n)\land(\gcd(d,f)=1)\right\}\right|$. We can find this by summing the above over all possible sets $F$, but we have to take into account that everything is counted twice as we are after unordered pairs.
So the final answer, $T$, is
$$T=\left|\left\{\left\{d,f\right\}:(d\,\,|\,\,n)\land(f\,\,|\,\,n)\land(\gcd(d,f)=1)\right\}\right|=\frac{1}{2}\left|\left\{(d,f):(d\,\,|\,\,n)\land(f\,\,|\,\,n)\land(\gcd(d,f)=1)\right\}\right|$$
$$=\frac{1}{2}\sum_{F\neq\varnothing\subset\{1,...,k\}} \left(\prod^k_{i=1}{(a_i+1)}-\prod_{i\in F}{(a_i+1)}-\prod_{i\not\in F}{(a_i+1)}+1\right)$$
Let $P=\prod^k_{i=1}(a_i+1)$. Then
$$T=\frac{1}{2}\left(\left(2^k-2\right)(P+1)-\sum_{F\neq\varnothing\subset\{1,...,k\}}\left(\prod_{i\in F}{(a_i+1)}+\prod_{i\not\in F}{(a_i+1)}\right)\right)$$
$$=\left(2^{k-1}-1\right)(P+1)-\frac{1}{2}\sum_{F\neq\varnothing\subset\{1,...,k\}}\left(\prod_{i\in F}{(a_i+1)}+\prod_{i\not\in F}{(a_i+1)}\right)$$
$$=\left(2^{k-1}-1\right)(P+1)-\sum_{F\neq\varnothing\subset\{1,...,k\}}\prod_{i\in F}{(a_i+1)}$$
$$=2^{k-1}(P+1)-\sum_{F\subseteq\{1,...,k\}}\prod_{i\in F}{(a_i+1)}$$
Not sure on simplifying it further. It seems this question, Elementary symmetric polynomials and Newton's Identities are related
Example calculations:
$n=30=2^13^15^1\implies k=3, P=8, (\forall i)\,\,a_i+1=2$
So, observing that there are $2^k-2=6$ possible subsets, which all are size $1$ or $2$,
$$T=27-(3(2)+3(4))=9$$
Not sure how more complicated examples can be done without a computer for verification.
I'm not so sure this particular answer relates to $$P=\prod_{i=1}^k(2a_i+1)$$ from the link.