A possible formulation is through gaussian processes. The unknown
smooth function $\eta(x)$ is seen as a stochastic process $Y(x)$ with
a given distribution that can be seen as a functional prior for
$\eta(x)$. Then the estimation $\widehat{\eta}(x)$ can be the
posterior mean, i.e. the expectation of $Y(x)$ conditional on the set
of inequalities $y_{i0} \le Y(x_i)$ and $Y(x_i) \le y_{i1}$. However,
it is not evident to find a software computing this efficiently.
For a more practical solution, the mid-point $y_i:=[y_{i0}+y_{i1}]/2$ can be taken as a
response at $x=x_i$ as you did. Then most least-squares based smoothing
methods can be adapted
by solving a Quadratic Programming (QP) problem instead of a
least squares problem: since QP routines
are widely available, a program will not be difficult to write.
For instance, a constrained smoothing spline
can be found by QP. Let $\mathbf{y}$ be the vector of the $n$ mid-points $y_i$,
and $\boldsymbol{\eta}$ be the vector of the unknown $\eta_i:=\eta(x_i)$.
In the usual spline smoothing, the estimate $\widehat{\boldsymbol{\eta}}$ is found by
the minimization
$$
\min_{\boldsymbol{\eta}} \: p \,\|\mathbf{y}-\boldsymbol{\eta}\|^2
+ (1-p) \,\boldsymbol{\eta}^{\mathrm{T}} \mathbf{M} \boldsymbol{\eta}
$$
where $\mathbf{M}$ is a $n\times n$ matrix with rank $n-2$ depending on the design
points $x_i$, and $0 <p < 1$
is a smoothing parameter.
For the constrained spline, we add the two sets of $n$ constraints:
$\boldsymbol{\eta} \ge \mathbf{y}_0$ and $\boldsymbol{\eta} \le \mathbf{y}_{1}$.
The unknown vector $\boldsymbol{\gamma}$
of $n-2$ "coefficients" i.e. of second order
derivatives at interior nodes, which is needed e.g. to interpolate
is related to $\boldsymbol{\eta}$ through $\boldsymbol{\eta} = \mathbf{K}\boldsymbol{\gamma}$ where the matrix $\mathbf{K}$
with dimension $n \times (n-2)$ is found (as well as $\mathbf{M}$) in the literature on smoothing splines.
At least in a first approach, the parameter $p$ can be guessed.
When $p \approx 0$, the smoothed curve will be allowed to depart from the
mid-points and to be close to either end-points to reach a greater
level of smoothing. However $p>0$ is needed to have a positive definite matrix
in the QP. The quadprog R package can be used for, say, $n \le 100$,
and larger problems can probably be decomposed in blocks.
