If I understand the question correctly, the idea is to resample the energy-usage rates conservatively. To ensure conservative resampling, you should resample the extensive quantity ("mass" = cumulative energy used) rather than the intensive quantity ("density" = usage rate). This is very similar to how resampling a probability density is tricky, but resampling a cumulative distribution is straightforward (i.e. no coordinate-change adjustment is required).
In the current case, we have a time series of cumulative energy usage $(E_k,t_k)$. The original question is phrased in terms of the average energy-usage rate (power), which is the ratio of first differences, i.e. $\bar{r}_k=\Delta E_k/\Delta t_k$, and notes difficulty in conservative resampling of the $\bar{r}_k$ time series.
However, if we consider the energy series itself, then we have
$$\Delta E_k = E_{k+1}-E_k = \int_{t_k}^{t_{k+1}}r[t]dt = \bar{r}_k\Delta t_k$$
The energy series itself $E(t)$ can be interpolated using any monotone scheme (e.g. linear interpolation, or monotone cubic). The monotone requirement ensures that it will always be non-decreasing through time: $E_k\leq E_{k+\phi}\leq E_{k+1}$ for $\phi\in[0,1]$.
Once this is done, the new energy series can be differenced to get average usage rates over the new time intervals.
Summary: the average usage rate is by definition $\bar{r}=\Delta E/\Delta t$, so if you interpolate the cumulative energy $E(t)$ (monotonically) then you will automatically have conservative results.
I do not use R, but skimming the help, it looks like you can do something like:
nt <- length(ts$start_date)
t <- c(ts$start_date,ts$end_date[nt])
E <- cumsum(c(0,ts$energy_use))
Espline <- splinefun(x = t, y = E, method = 'monoH.FC')
dEdt_spline <- function(t) Espline( t , deriv = 1 )
Then you can evaluate the average power consumption as $\langle r\rangle_{t\in[t_1,t_2]} = \frac{E(t_2)-E(t_1)}{t_2-t_1}$, and the "instantaneous" power consumption with $E'(t)$.
(Note: I quickly tried this on R-fiddle and it seemed to work, but your integrate test still did not work. I strongly believe this must be due to some code error, either on my part or in the R libraries. That is, by design $E(t)$ has been fit with a monotone interpolating spline, which has an analytic derivative, that itself has an analytic integral, as they are piecewise polynomial functions. Most likely the inconsistency is due to my lack of R knowledge, or it could be due to numerical approximations used in the spline calculus functions.)
Update: As expected, the above was due to my lack of R knowledge, as shown in updated question. (I had literally never written any R before Googling to do the above, so not too shabby!) Note also that as seen there, the monotone cubic spline functions will have a discontinuous second derivative (seen as kinks in the $E'(t)$ plots). This could be avoided with a monotone C2 interpolant (e.g. this), though I do not know what R package this might be in.
Note on monotone interpolation: This simply means that the interpolation does not introduce any new local maxima/minima not already present in the data. For example the following picture from Wikipedia demonstrates how the standard cubic spline is not monotone
Public Domain, https://en.wikipedia.org/w/index.php?curid=9051137
Best Answer
Here is one approach using cut() to create the appropriate hourly factors and ddply() from the plyr library for calculating the means.