Reconstructing band-limited signal's integral from discrete samples
Let me rephrase what you want to do as:


  I want to approximately calculate the integral of the product of two non-negative discrete-time signals over time. Each signal is sampled from a band-limited continuous-time signal at a sampling frequency $2N$ times its highest frequency. A method is preferred that gives about 5&nbsp;% error at most, is computationally efficient and does not require a large oversampling factor $N.$


Oversampling factor

The bandwidth of the product of two signals is equal to the sum of their bandwidths. Therefore an oversampling factor of $N = 2$ will suffice to represent the product without aliasing.

Integration

Integration is a linear time-invariant operation, basically a continuous-time filter, and applying it does not increase the bandwidth of the signal. What we want to achieve is, with the product signal as "input":

$$\xrightarrow{\text{input}}\boxed{\text{integrate}}\xrightarrow{\text{output}}$$

However, the impulse response of an integrator is of infinite bandwidth and sampling it to form the impulse response of a discrete-time integrator would cause aliasing, corrupting its frequency response. We can insert to the signal path an ideal lowpass filter that has its cutoff at the bandlimit of the input signal. It does nothing to the signal. The output will remain the same as before:

$$\xrightarrow{\text{input}}\boxed{\text{lowpass}}\xrightarrow{\text{input}}\boxed{\text{integrate}}\xrightarrow{\text{output}}$$

We can combine the lowpass filter and integrator to a single filtering operation. The impulse response of this lowpass integrator filter is the convolution of the impulse responses of the lowpass filter and the integrator filter. This shuffling of operations is enabled by the algebraic properties of convolution, namely associativity. We still get the same desired output as before:

$$\xrightarrow{\text{input}}\boxed{\text{lowpass integrate}}\xrightarrow{\text{output}}$$

The lowpass integrator impulse response:

$$\begin{align}
h[k] &amp;= \int_{-\infty}^k\frac{\sin(\pi x)}{\pi x}dx \\
\\
 &amp;= \int_0^k\frac{\sin(\pi x)}{\pi x}dx + \frac{1}{2}, \\
\end{align}$$

is sampled from the integral of the sinc function, without aliasing. Sinc is the impulse response of the ideal continuous-time low-pass filter with a cutoff frequency at $\pi$.

Figure 1. Impulse response $h[k]$ of an integrator. The impulse response continues indefinitely beyond what is shown, approaching values 0 to the left and 1 to the right.

The impulse response almost represents the process of calculating the sum of (optionally the current sample and) the past samples. As this sum becomes larger and larger, the difference to true integration becomes arbitrarily small, proportionally. The values of $h[k]$ approaching 1 for large $k$ means that for very old samples just summing them gives virtually no error in the integral estimate.

Summary

To summarize, oversample the input signals by a factor of two and calculate with a large accumulator a running sum of obtained samples of the product. This will have an arbitrarily small error compared to true integration, as enough data gets accumulated.

Analog filtering

You should also consider what the analog filtering of the current and voltage signals does to the integral of their product. This can be a big source of error. For example let's say that both signals equal $\sin^2(x).$ If this is filtered to remove all oscillation, it becomes a constant signal $\frac{1}{2}$. The average value (power) of the square of the unfiltered signal is $\frac{3}{8}$ and for the filtered signal $\frac{1}{4}$.