This is an interesting question and as @user100927 has correctly commented, a necessary condition for this to be possible is:
$$\int_{x_i}^{x_{i+1}}f(x)dx ≤ y_{i+1}−y_i$$ for all $i$.
To prove this condition let us denote by $s(x)$ the derivative of the interpolating function $S(x)$, and by $F(x)+c$ the anti-derivative (or indefinite integral) of $f(x)$.
We define the function $h(x) = s(x)-f(x)$, which by our restriction is positive. Then $H(x) = S(x)-F(x) + c$, the anti-derivative of $h(x)$, is a monotone increasing function (since it's derivative is positive).
Thus, since $x_i < x_{i+1}$, we get $H(x_{i+1}) - H(x_i) \geq 0$.
However,
$$H(x_{i+1}) - H(x_i) = (S(x_{i+1}) - F(x_{i+1})) - (S(x_i)-F(x_i)) =\\ (S(x_{i+1}) - S(x_{i})) - (F(x_{i+1}) - F(x_i))$$
and from the interpolation condition and the fundamental theorem of calculus
we get: $y_{i+1}-y_i - \int_{x_i}^{x_{i+1}}f(x)dx \geq 0$, which is our necessary condition.
So the answer to your question "Is this possible?" is "not generally".
However, assuming that the condition is met, we can think of a method to construct an interpolating function that meets the requirements.
For the rest of the answer I assume $f(x)$ is a function that has a known anti-derivative $F(x)$ (e.g., a polynomial or a spline). If not, we bound it from above by such a function and consider that to be our $f$.
I denote by $C_{i+1}$ the positive value from the above proof, i.e., $C_{i+1} = y_{i+1}-y_i - \int_{x_i}^{x_{i+1}}f(x)dx \geq 0$.
We can then compute $S(x)$ in the interval $[x_i, x_{i+1}]$ to be:
$$S(x)=y_i + F(x)-F(x_i) +C(x)$$
where $C(x)$ is a function with the following properties:
- $C(x_i) = 0$
- $C(x_{i+1}) = C_{i+1}$
- $C'(x) \geq 0$ for all $x \in [x_i, x_{i+1}]$.
From (1) we get the interpolation at $x_i$, from (2) we get the interpolation at $x_{i+1}$, and from (3) we get $S'(x) \geq f(x)$.
These three conditions can easily be met, for example, with the function:
$C(x) = \frac{C_{i+1}}{x_{i+1}-x_i}(x-x_i)$.
This will result in a continuous interpolating function, but the function will not be smooth.
We can get a smooth $C^1$ function if we set $C(x)$ to be a cubic Hermite spline, where condition (3) can be met if we set the end derivatives (at $x_i$ and $x_{i+1}$) to zero. Assigning $t=\frac{x-x_i}{x_{i+1}-x_i}$, this leads to the function (see here): $C(x)=C_{i+1}(-2t^3 + 3t^2)$
Interestingly, we cannot in general achieve a smooth function using quadratic polynomial segments. For a counter example think of the data: $y_0=0, y_1=0, y_2=1, y_3=1$ for the corresponding $x$ values $x_i = i$, and $f(x) \equiv 0$. It
can be shown (e.g., using Rolle's theorem) that a parabola interpolating identical $y$-values with a non-negative derivative must be a constant. Therefore, $S(x)$ over $[0,1]$ will have to be identically zero, and similarly over $[2,3]$ will be identically one. But this means that for $S(x)$ to be smooth it must pass a parabola starting at $(1,0)$ and ending at $(2,1)$ with zero derivative at each end point and this is impossible.
To summarize: Given a function $f(x)$ that we can compute its anti-derivative $F(x)$, and assuming the interpolation points meet the necessary condition (if not we will discover this along the way), we can compute the required smooth interpolation function (assuming $F(x)$ is smooth) by setting $S(x)=y_i + F(x)-F(x_i) +C(x)$ for every interval $[x_i, x_{i+1}]$. Where $C(x)$ is a cubic polynomial segment - a Hermite spline over $[x_i, x_{i+1}]$ that satisfies $C(x_i)=0$, $C(x_{i+1})=C_{i+1}$, and $C'(x_i)=C'(x_{i+1})=0$.
This result cannot, in general, be achieved by a regular spline (i.e., a quadratic $C^1$ spline). If we wish to achieve $C^2$ continuity, we can do it with degree-5 polynomial segments (I haven't shown this, but this can be done by constraining the second derivatives at the end points to also be zero). I suspect, although I haven't worked out an example so this is just an educated guess, that similar to the quadratic case, $C^2$ continuity cannot be achieved with less than degree-5.
