I don't know if you were able to solve your problem yet, but here is a paper proposing a method for mean-preserving quadratic splines, which can take additional constraint for lower bounds :

https://www.sciencedirect.com/science/article/pii/S0038092X22007770

doi : 10.1016/j.solener.2022.10.038

I have used it for approximating daily cycles of instantaneous surface temperature from 3h-mean values, and the method is particularly robust and efficient. The authors provide the python interface to implement the method in your code.