Calculating total power for the signal
How can I calculate the total power? It is &quot;the variance plus the mean squared&quot;, so I should just calculate $(\sum_{i}^N x[i]) / N$ for the signal of length $N$?

$(\sum_{i}^N x[i]\color{red}{^2}) / N$ is the average power, not the total energy. Just the sum $\sum_{i}^N x[i]\color{red}{^2}$ is what you're looking for.

Or just sum the PSD (like in this link)?

and

Should I perform this calculation in the time domain (just time series like above), or should I calculate FFT and PSD (Power Spectral Density) first and calculate on PSD (which is also just an array of numbers, if I understand correctly)?

The sum of the power in time domain (as above) and in frequency domain are the same, so that's just another way to compute the same (probably just more labor-intense).

I have seen quite a few formulas for the total power. For finite number of samples and discrete measurements I have seen (∑Nix[i]) divided by N, N+1, 2N or 2N+1. Which one should I use? Does it have something to do with bias/unbiased sampling, i.e. in statistics we add +1 to the divisor to get unbiased calculations?

Since you're really not looking for the power but the energy, neither. You don't have to divide at all.