I am not even sure I have a question anymore (I will just give up)... in the past month or so I have been researching cubic Bézier curves. The idea was to find a fit through data points, using piecewise Bézier curves and the Euler-Lagrange equation to minimize an "energy" proportional with the square of the length of the (total) curve (elastic energy used to stretch) plus the integral of the square of the curvature $\kappa$ of the curve (energy required to bend the curve).
Edit: it took me three whole days to type this "question". Please don't get mad for wasting your time if you choose to read it.
Here is what I have done so far with respect to the length of the curve.
Let $$\begin{equation}\begin{aligned}\vec{P}(t) &= (1-t)^3 \vec{P_0} + 3(1-t)^2t \vec{C_0} + 3(1-t)t^2 \vec{C_1} + t^3 \vec{P_1} \\ &= \vec{P_0} + 3(\vec{C_0} - \vec{P_0})t + 3(\vec{C_1} - 2\vec{C_0} + \vec{P_0})t² + (\vec{P_1} - 3\vec{C_1} + 3\vec{C_0} - \vec{P_0})t³\end{aligned} \end{equation}\tag{1}$$ be a cubic Bézier curve defined by its begin and end points $\vec{P_0}$ and $\vec{P_1}$, and control points $\vec{C_0}$ and $\vec{C_1}$.
Instead of $\{\vec{P_0}, \vec{C_0}, \vec{C_1}, \vec{P_1}\}$ to define the curve, lets use $\{\vec{B}, \vec{V}, \vec{A}, \vec{J}\}$ where $$\begin{align}\vec{B}&=\vec{P_0} \\ \vec{V}&=3(\vec{C_0} - \vec{P_0}) \\ \vec{A}&=6((\vec{C_1} - \vec{P_0}) - 2(\vec{C_0} - \vec{P_0})) \\ \vec{J}&=6((\vec{P_1} - \vec{P_0}) - 3(\vec{C_1} - \vec{C_0}))\end{align} \tag{2}$$
This is not random, you will thank me later; the idea behind it is that the coefficients other than the constant term represent the velocity, acceleration and jerk (if $t$ was "time") with respect to $t$, which only depends on the shape of the curve, not its position. Therefore those coefficients should be expressed in terms of differences between any two of $\{\vec{P_0}, \vec{C_0}, \vec{C_1}, \vec{P_1}\}$.
Substituting this into $(1)$ gives, $$\vec{P}(t) = \vec{B} + \vec{V}t + \frac{1}{2}\vec{A}t^2 + \frac{1}{6}\vec{J}t^3 \tag{3}$$
Then the chord length of the curve between $\vec{P_0}$ and $\vec{P_1}$ is given by
$$chord\_length = \int_0^1 \lVert \vec{P'}(t) \rVert \,dt \tag{4}$$
where $\vec{P'}(t)$ is the derivative of $\vec{P}(t)$ with respect to $t$, also called the velocity:
$$\vec{P'}(t) = \vec{V} + \vec{A}t + \frac{1}{2}\vec{J}t^2 \tag{5}$$
In order to calculate the norm $\lVert \vec{P'}(t) \rVert$ we can take the square root of the dot product with itself: $$\lVert \vec{P'}(t) \rVert = \sqrt{\vec{P'}(t) \cdot \vec{P'}(t)} \tag{6}$$
Working out this dot product we get a fourth degree polynomial in $t$: $$\vec{P'}(t) \cdot \vec{P'}(t) = c_4 t^4 + c_3 t^3 + c_2 t^2 + c_1 t + c_0 \tag{7}$$ where $$\begin{aligned}c_0 &= \lVert \vec{V} \rVert^2 \\ c_1 &= 2 (\vec{V} \cdot \vec{A}) \\ c_2 &= \vec{V} \cdot \vec{J} + \lVert \vec{A} \rVert^2 \\ c_3 &= \vec{A} \cdot \vec{J} \\ c_4 &= \frac{1}{4}\lVert \vec{J} \rVert^2\end{aligned} \tag{8}$$
Note how all of this holds for vectors of dimension two or more. We're not talking about $x$ and $y$ or whatever.
We might also be interested in a polynomial in $u$ after a substitution $t = u + s$, where we'll integrate from $u=-s$ to $u=-s+1$. After substitution we end up with: $$\begin{align}chord\_length &= \int_{-s}^{-s+1} \sqrt{c_4' u^4 + c_3' u^3 + c_2' u^2 + c_1' u + c_0'} \space \,du \\ &= \sqrt{c_0'} \int_{-s}^{-s+1} \sqrt{\frac{c_4'}{c_0'}u^4 + \frac{c_3'}{c_0'}u^3 + \frac{c_2'}{c_0'} u^2 + \frac{c_1'}{c_0'} + 1}\space\space \,du \end{align} \tag{9}$$
where $$\begin{aligned}c_0' &= c_0 + c_1 s + c_2 s^2 + c_3 s^3 + c_4 s^4 \\ c_1' &= c_1 + 2 c_2 s + 3 c_3 s^2 + 4 c_4 s^3 \\ c_2' &= c_2 + 3 c_3 s + 6 c_4 s^2 \\ c_3' &= c_3 + 4 c_4 s \\ c_4' &= c_4\end{aligned} \tag{10}$$
Let $$\begin{aligned}d_1 &= \frac{c_1'}{c_0'} \\ d_2 &= \frac{c_2'}{c_0'} \\ d_3 &= \frac{c_3'}{c_0'} \\ d_4 &= \frac{c_4'}{c_0'}\end{aligned} \tag{11}$$
So that the problem of finding the chord length becomes integration of $\sqrt{1 + d_1 u + d_2 u^2 + d_3 u^3 + d_4 u^4}$.
Since there is no known algebraic expression for that, we can proceed by finding the Taylor series of this square root around $u=0$ ($t=s$). Lets write $$\sqrt{1 + d_1 u + d_2 u^2 + d_3 u^3 + d_4 u^4} = 1 + \sum_{n=1}^{\infty} \frac{1}{n}g_n u^n \tag{12}$$
then the following recursive relationship for the $g_n$ can be found (with three days of hard work): $$\begin{aligned} g_1 &= +\frac{1}{2}d_1 \\ g_2 &= -\frac{1}{2}d_1 g_1 + \frac{1}{1}d_2 \\ g_3 &= -\frac{3}{4}d_1 g_2 - \frac{0}{1}d_2 g_1 + \frac{3}{2}d_3 \\ g_4 &= -\frac{5}{6}d_1 g_3 - \frac{1}{2}d_2 g_2 + \frac{1}{2}d_3 g_1 + \frac{2}{1}d_4 \\ g_5 &= -\frac{7}{8}d_1 g_4 - \frac{2}{3}d_2 g_3 - \frac{1}{4}d_3 g_2 + \frac{1}{1}d_4 g_1 \\ &\space\space\vdots \notag \\ g_n &= -\frac{2n-3}{2n-2}d_1 g_{n-1} - \frac{n-3}{n-2}d_2 g_{n-2} - \frac{2n-9}{2n-6}d_3 g_{n-3} - \frac{n-6}{n-4}d_4 g_{n-4}\end{aligned} \tag{13}$$
This Taylor series should converge if $u$ is close enough to zero. In fact, it turns out that for tame curves where, say the deviation from a straight line is less than a third of the distance between $\vec{P_0}$ and $\vec{P_1}$, one can calculate the chord length with a single integral taking $s=1/2$ and integrating from $-1/2$ to $+1/2$. Otherwise we might have to do it piecewise and add up a few chord length pieces. For example one could pick six pieces where $s$ is one of $\{\frac{1}{12}, \frac{3}{12}, \frac{5}{12}, \frac{7}{12}, \frac{9}{12}, \frac{11}{12} \}$ and then integrate (six times) from $-\frac{1}{12}$ to $+\frac{1}{12}$. Thus, in general we are interested in obtaining the integral between $-\delta$ and $+\delta$ where $\delta$ is some fixed, positive real less than or equal $1/2$.
Combining equation $(9)$ through $(13)$ we can write, $$chord\_length\_piece = \sqrt{c_0'} \int_{-\delta}^{\delta}\left(1+\sum_{n=1}^{\infty}\frac{1}{n}g_n u^n\right)\,du \tag{14}$$
Assuming $\delta$ is small enough that the sum converges, we can swap the sum and integration and get $$chord\_length\_piece = 2\delta\sqrt{c_0'} \left( 1 + \sum_{even\ n\geq 2}^{\infty} \frac{\delta^n}{n(n+1)}g_n \right) \tag{15}$$
Note that only the $g_n$ for even $n$ are left because the odd ones cancel. For example $\int_{-\delta}^{\delta} g_1 u\,du = \left[\frac{1}{2}g_1u^2\right]_{-\delta}^{\delta} = \frac{1}{2}g_1\delta^2 - \frac{1}{2}g_1(-\delta)^2 = 0$. The other terms are doubled, hence the leading 2. I also moved one factor of $\delta$ to the front, making all exponents of $\delta$ equal to the index of $g_n$ again. For a straight line, the sum is zero and $\vec{c_0'} = \lVert \vec{V_s} \rVert^2$ the square of the velocity at $t=s$.
The total chord length then can be found by summing $(15)$ for all values of $s$ involved. Let $\delta = \frac{1}{2N}$ for some positive integer $N$, the number of pieces we will cut the curve into. Then we need $N$ values of $s$: $$s = (2k + 1)\delta \qquad\text{ for } k = \mathbin{{0}{...}{N}{-}{1}}$$
and $$\begin{aligned}chord\_length &= 2 \delta \sum_{k=0}^{N-1} \sqrt{c_0'^{(s)}} \left( 1 + \sum_{even\ n\geq 2}^{\infty} \frac{\delta^n}{n(n+1)}g_n^{(s)} \right) \\ &= 2 \delta \sum_{k=0}^{N-1} \sqrt{c_0'^{(s)}} + 2 \delta \sum_{even\ n\geq 2}^{\infty} \frac{\delta^n}{n(n+1)} \sum_{k=0}^{N-1} g_n^{(s)} \end{aligned} \tag{16}$$
Recall that the $g_n$ are a function of the $d_i$ which in turn are a function of the $c_i'$ which are polynomials in $s$ (hence the notation $g_n^{(s)})$ with coefficients $c_i$, where $s = (2k + 1)\delta$. Since $\delta$ is a known, fixed constant, the second term of equation $(16)$ is a multivariate polynomial in the $c_i$.
The first term contains the $\sqrt{c_0'^{(s)}}$, which will carry over to a lot of other terms after squaring the length; but is still straightforward to differentiate. Note that $c_0'$ is the square of the velocity at $t=s$, so that if we make $N$ arbitrary large (causing $\delta$ to go to zero) the first sum turns into the original integral (equation $(4)$).
At this point I feel like giving up though: assuming I can come up with a differentiable equation for the bending energy part, I am starting to think this is going to be a dead-end. In order to calculate the ideal curve I'd have to minimize the energy of the whole curve, all Bézier curves combined, where one curve influences the next only due to the demand that the curve must have $G_1$ continuity.
Am I right to give up, or does anyone have an idea that can be done in a reasonable amount of time?
Not really an answer, but far too long for a comment. A few remarks …
Short answer … no, I don’t think you should give up. This seems like a solvable problem, to me.
If you decide that writing your own code is not worth the effort, there are plenty of spline fitting packages available. Here is one. I have never used it, but I have great respect for its author.