Resampling/Interpolating monthly rates to daily rate estimates in R I'm really not sure what to search for. If the answer to this is googleable, I'd be happy to hear what I should google.
I have a dataset of energy meter data. Readings are taken at roughly monthly intervals, though it is irregular. The data describe total energy use since the previous reading was taken. I want to resample that data to estimate usage for each month (the first to the last day). I'm looking into linear and spline interpolation, but both methods lead to high error rates. As shown below, integrating either of these methods for the known energy use of an interval results in some very large differences.
Is there a way to generate an interpolation for data like this while still being able to match the function over known ranges?
# Create Sample Meter Timeseries
ts <- data.frame(start_date = as.Date(c("2015-09-18",
                                        "2015-10-20",
                                        "2015-11-25",
                                        "2015-12-23",
                                        "2016-01-22")),
                 end_date = as.Date(c("2015-10-20",
                                      "2015-11-25",
                                      "2015-12-23",
                                      "2016-01-22",
                                      "2016-02-21")),
                 energy_use_kWh = c(3211,172,5566,8120,1344))

# Convert interval use into daily rate
ts$rate <- ts$energy_use_kWh/as.numeric(ts$end_date - ts$start_date)

# Generate midpoint of each interval
ts$mid_date <- ts$start_date + (ts$end_date - ts$start_date)/2

# Generate linear and spline interpolations between interval midpoints
approxfun_fun <- approxfun(x = ts$mid_date, y = ts$rate, rule = 2)
splinefun_fun <- splinefun(x = ts$mid_date, y = ts$rate, method = 'natural')

# Compare
ts$approx_estimate <- apply(ts[,c('start_date','end_date')], MARGIN = 1,
                            function(x) integrate(f = approxfun_fun,  
                            lower = as.Date(x[1]), upper = as.Date(x[2]))$value)

ts$spline_estimate <- apply(ts[,c('start_date','end_date')], MARGIN = 1,
                            function(x) integrate(f = splinefun_fun,  
                            lower = as.Date(x[1]), upper = as.Date(x[2]))$value)

result <- ts[c("start_date","end_date","energy_use_kWh","approx_estimate","spline_estimate")]
result$approx_error <- with(result, abs(energy_use_kWh - approx_estimate)/energy_use_kWh)
result$spline_error <- with(result, abs(energy_use_kWh - spline_estimate)/energy_use_kWh)

EDIT: I have implemented GeoMatt22's suggestion over my sample data and found some interesting results...

Taking the integral over the known billing ranges matches perfectly with the known data, so that is definitely a success. The cumulative energy use graph looks pretty reasonable, though the derivative graph definitely looks a little strange. Within a billing period, the rate change is smooth and continuous, but there are some cases of instantaneous usage rate changes, where one day it is decreasing rapidly and the next day it is increasing rapidly.
This is the best answer so far, even though the day-to-day rates look a little funky. I'd welcome improvements that address this, but I am also happy with what I've got so far.
Edit 2: I tried changing the spline method from "monoH.FC", which computes a monotone Hermite spline, to "Hyman", which computes a monotone cubic spline using Hyman filtering. The results I got are a bit more continuous, and pass my personal eye-test a little better, though it's still not nicely continuous.

Edit 3: I built a function for the C2 monotone interpolant, as suggested by GeoMatt22. It took me a long time, but I got it to work! Method is from C2 rational quadratic spline interpolation to monotonic data (1983).
If you want more info on the code, I'd be happy to share. Currently in development, but it seems to be working well. The derivative curve, which corresponds to the daily use rate, has no unsightly kinks and integrating that function over each billing interval reproduces the known values exactly.

 A: 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
