Firstly, note that hybridisation theory as applied to transition metal complexes is an incorrect, flawed theory. It is an attempt to rationalise experimental observation (e.g. geometries), but in the process it invokes models of bonding which are highly unrealistic and demonstrably untrue. See: Why is it wrong to use the concept of hybridization for transition metal complexes?
Now, with that disclaimer out of the way, let's examine how that flawed theory works. In traditional hybridisation theory, tetrahedral and square planar 4-coordinate complexes are considered to be $\mathrm{sp^3}$ and $\mathrm{dsp^2}$ respectively. Since $\ce{[Ni(CO)4]}$ is tetrahedral and $\ce{[Ni(CN)4]^2-}$ is square planar, one can assign the hybridisations accordingly. If you have a textbook that still mentions this outdated theory, it should have a table telling you what hybridisation each geometry corresponds to. Here is an example from Housecroft's Inorganic Chemistry 4th ed. (p 667). To their credit, there is a subsequent section debunking this theory.
These are generally based on post hoc justifications. Essentially, somebody looked at the geometry and then tried to piece together which orbitals could be combined to form such a geometry. So, I do not think there is a good way to predict them a priori. This formula you have dug up does not work very well, even for main group elements. For transition metals the factors deciding geometries are much more complex, and simply cannot be captured in one line of mathematics.
The real problem with the approach above is that you need to know the geometry before you can assign a hybridisation. (Don't even try to do it the other way round - it simply does not work because hybridisation is an incomplete theory. It is like trying to predict wave-particle duality using classical mechanics.)
If you cannot predict the geometry of the complex, then you cannot predict the hybridisation. You could argue that cyanide is a strong field ligand and in conjunction with a $\mathrm{d^8}$ metal, a square planar geometry is preferred. On the other hand, in $\ce{[Ni(CO)4]}$, the central atom is $\mathrm{d^{10}}$ and hence even with a strong-field ligand like carbon monoxide, a tetrahedral geometry is favoured, simply based on steric grounds.
However, this does nothing but introduce more questions. The term "strong-field" has no meaning outside of molecular orbital theory. Hybridisation theory has no inkling of what strong- or weak-field ligands are, and consequently cannot explain why $\mathrm{d^8}$ Ni(II) only prefers to go square planar with strong-field ligands. On top of that, it is not at all obvious from hybridisation theory how the electronic configuration of the central metal should favour a particular geometry. It's simply not possible for your formula to explain any of this.
So, to answer your questions directly:
Is there a systematic and logical way of knowing hybridisation of coordination compounds?
Yes. You use a more advanced theory (MOT) to predict the geometry, and then you use the table above to find the hybridisation. That's a step backwards though, you already used MOT and there's no reason why you should go back to hybridisation.
Why this formula does not work for coordination compounds?
Because hybridisation theory is an incomplete theory and cannot capture the subtle factors influencing geometry of coordination compounds.
Main group compounds
I guess I will elaborate a bit more on main group compounds. The formula you quoted seems to work for main group compounds, but let's get two things out of the way.
It gives the correct hybridisations as espoused by old theory. $\ce{SF6}$ is not $\mathrm{sp^3d^2}$ hybridised because d-orbital involvement in bonding is negligible. So, this formula does not actually work.
It does not predict the hybridisation without predicting the geometry.
Point 1 is important, but let's not go into that now. Point 2 may seem confusing to you. Clearly it is predicting the hybridisation, simply from the number of valence electrons, etc. etc.?
The truth is that this formula only "works" because these quantities (number of valence electrons, number of bond pairs, ...) allow one to predict the geometry, using VSEPR theory. For example, the structure of $\ce{CH4}$ is tetrahedral, and VSEPR theory is sufficient to predict that. Once the geometry has been established as tetrahedral, the hybridisation can then be correctly assigned as $\mathrm{sp^3}$.
So, this formula is skipping a step in between. In fact, the correct sequence is to predict the geometry using VSEPR theory, and then obtain the hybridisation. Now, with transition metal complexes, VSEPR does not work. Thus, it is no surprise that that formula also does not work.
