It looks like your set of endpoints and control points can be
any set of points in the plane.
This means that the $order$ of the points
is critical,
so that the generated curve goes through the points
in a specified order.
This is much different than the ordinary interpolation problem,
where the points
of the form $(x_i, y_i)$
are ordered so that
$x_i < x_{i+1}$.
As I read your desire,
if you gave a set of points on a circle
ordered by the angle of the line
from the center to each point,
you would want the result to be
a curve close to the circle.
There are a number of ways this could be done.
I will assume that
you have $n+1$ points
and your points are $(x_i, y_i)_{i=0}^n$.
The first way I would do this
is to separately parameterize the curve
by arc length,
with $d_i = \sqrt{(x_i-x_{i-1})^2+(y_i-y_{i-1})^2}$
for $i=1$ to $n$,
so $d_i$ is the distance from
the $i-1$-th point to the
$i$-th point.
For a linear fit,
for each $i$ from $1$ to $n$,
let $t$ go from
$0$ to $d_i$
and construct separate curves
$X_i(t)$ and $Y_i(t)$
such that
$X_i(0) = x_{i-1}$,
$X_i(d_i) = x_i$,
and
$Y_i(0) = y_{i-1}$,
$Y_i(d_i) = y_i$.
Then piece these together.
For a smoother fit,
do a spline curve
through each of
$(T_i, x_i)$
and
$(T_i, y_i)$
for $i=0$ to $n$,
where
$T_0 = 0$
and
$T_i = T_{i-1}+d_i$.
To get a point for any $t$ from
$0$ to $T_n$,
find the $i$ such that
$T_{i-1} \le t \le T_i$
and then,
using the spline fits
for $x$ and $y$
(instead of the linear fit),
get the $x$ and $y$ values from their fits.
Note that
$T_i$ is the cumulative length
from $(x_0, y_0)$
to $(x_i, y_i)$,
and $T_n$ is the total length of the line segments
joining the consecutive points.
To keep the curves from
not getting too wild,
you might look up "splines under tension".
Until you get more precise,
this is as far as I can go.
Your thinking is correct -- you need to use a parametric curve (either a spline or a single polynomial) that gives $x$ and $y$ as functions of some parameter $t$. To compute the curve, you must first assign a parameter value $t_i$ to each point $P_i$. As you say, these parameter values are entirely fabricated -- they are not an intrinsic part of the original interpolation problem. Furthermore, different choices for these parameter values will produce significantly different curves. The usual approach is to choose them so that $t_{i+1}-t_i$ is equal to the chord-length between the points $P_i$ and $P_{i+1}$.
Using arclength seems like the natural approach, as you say, but this doesn't work because you don't know the arclengths until you've constructed the curve -- a chicken-and-egg problem.
There's a bit more info in this answer.
A good reference is this web site.
Best Answer
There was a paper written by Deddi et al (2000) on interpolation with curvature constraints using quadratic Bezier curves. When doing curvature constraints, you also need to include some form of curve length minimization to avoid smooth but completely unusable cases.
There is also the curvature minimizing Clough-Tocher scheme in 2D, though it is non parametric.
It is a tricky problem, and although the goal sound appealing, the fact that it is not widely used could suggest that it is not very reliable/practical in real applications. If you can accept that in order to smooth your curve, the interpolation may not pass through all the points (which is acceptable for noisy data), the splprep routine adapted from FITPACK might still be the simplest choice.