My strategy is to define every piece $f_k(x)$ of $f(x)$ so that it starts exactly at $(x_{k-1}, d_{k-1})$ and ends at $(x_k,d_k)$, where it meets the next piece ensuring continuity. If this can be done, it means that the complete function not only goes close to the target points, but it touches them all, making the required sum null (thus minimum). Moreover it wouldn't be necessary to force the condition $\sum d_i = 0$.
So, for the $k^{th}$ piece I am imposing the conditions:
$$
\begin{cases}
A_k \sin(\pi x_{k-1} + \theta_k) = d_{k-1}\\
A_k \sin(\pi x_k + \theta_k) = d_k
\end{cases}
$$
This is not a linear problem, but thinking naively we have two conditions to be satisfied tuning two parameters ($A_k$ and $\theta_k$) so that it is worth trying. For $k=1$ and $k=N+1$ only the innermost side is constrained so it is trivial to find compliant pieces $f_1(x)$ and $f_{N+1}(x)$.
I'll extract $\theta_k$ from the first condition. Due to the nature of the $\sin$ function and its inverse, I have two possibilities:
$$
\theta_k = \arcsin \left(\frac{d_{k-1}}{A_k}\right)-\pi x_{k-1}\\
\theta'_k = \pi - \arcsin \left(\frac{d_{k-1}}{A_k}\right)-\pi x_{k-1}
$$
I've found that the second one ends up giving the same function as the first, so I'll look only at $\theta_k$. Putting it in the second condition I have
$$\require{cancel}
A_k \sin \left( \underbrace{\pi x_k - \pi x_{k-1}}_{:= \ \beta} + \arcsin \left(\frac{d_{k-1}}{A_k}\right) \right)=d_k\\
A_k \cos(\beta) \sin \left(\arcsin\left(\frac{d_{k-1}}{A_k}\right)\right)+A_k\sin(\beta)\cos\left(\arcsin\left(\frac{d_{k-1}}{A_k}\right)\right)=d_k\\
\cancel{A_k}\cos(\beta)\frac{d_{k-1}}{\cancel{A_k}}\pm A_k\sin(\beta)\sqrt{1-\frac{{d_{k-1}}^2}{A_k^2}}=d_k\\
\cos(\beta)d_{k-1}\pm \sin(\beta)\sqrt{A_k^2-d^2_{k-1}}=d_k\\
\pm\sqrt{A_k^2-d^2_{k-1}}=\frac{d_k-\cos(\beta)d_{k-1}}{\sin(\beta)}\\
A_k^2=\left(\frac{d_k-\cos(\beta)d_{k-1}}{\sin(\beta)}\right)^2+d^2_{k-1}\\
A_k^2=\frac{1}{\sin^2(\beta)}\left(d_k^2-2\cos(\beta)d_k d_{k-1}+\cos^2(\beta)d^2_{k-1}+\sin^2(\beta)d^2_{k-1}\right)\\
A_k^\pm=\pm\frac{1}{\sin(\beta)}\sqrt{d_k^2+d^2_{k-1}-2\cos(\beta)d_k d_{k-1}}
$$
Now I can go back to plug these parameters in the first condition to get
$$
\theta_k^\pm=\pm\arcsin \left( \frac{d_{k-1}\sin(\beta)}{\sqrt{d_k^2+d^2_{k-1}-2\cos(\beta)d_k d_{k-1}}} \right)-\pi x_{k-1}
$$
Finally for each $k=2 \ldots N$ we have:
$$
f^\pm_k(x) = A^\pm_k \sin(\pi x + \theta^\pm_k)
$$
where it turns out that only one of the two satisfies the rightmost condition too. Gluing together the correct ones I build the complete $f(x)$.
I've tested my result generating some random sets $\{ x_i \}$ (uniform distribution in $(0,1)$ ) and $\{ d_i \}$ (normal distribution around $0$) for different $N$'s. Here below the plots (dashed lines when $f^-_k$ was chosen). Matlab's code is available if needed.
