The number of conjugacy classes of pairs of commuting elements in a finite group G is the cardinality of the set {c(a,b) | a,b in G and ab=ba} where c(a,b) = {(gag^(-1),gbg^(-1)) | g in G}.
It is equal to the number of conjugacy classes within the centralizers of class representatives of G.
This reformulation was employed in the sequence-generating program.
It is also equal to the rank of the modular fusion category Z(Rep(G)), the Drinfeld center of Rep(G).
These reformulations are explained in the linked MathOverflow posts.
REFERENCES
A. Davydov, Bogomolov multiplier, double class-preserving automorphisms, and modular invariants for orbifolds. J. Math. Phys. 55 (2014), no. 9, 092305, 13 pp.