Suppose the 2D version of the Whittaker smoothing spline
https://eigenvector.com/wp-content/uploads/2020/01/WhittakerSmoother.pdf
minimizing
$$ \phi=(X-X_{s})^{T}(X-X_{s})^{T}+(Y-Y_{s})^{T}(Y-Y_{s})+X_{s}^{T}D_{1}^{T}D_{1}X_{s}+Y_{s}^{T}D_{1}^{T}D_{1}Y_{s}, $$
where
$$ D_{1}=\left[\begin{array}{cccccc} -1 & 1 & 0 & \cdots & 0 & 0\\ 0 & -1 & 1 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ 0 & 0 & 0 & \cdots & -1 & 1 \end{array}\right], $$
is the difference matrix, and $X,Y$ are coordinates of control points $p_i$.
Using the condition $\frac{\partial\phi}{\partial X_{s}}=0,$ $\frac{\partial\phi}{\partial Y_{s}}=0$, we get the following spline equations
$$ X_{s}=(I+D_{1}^{T}D_{1})^{-1}X,\qquad Y_{s}=(I+D_{1}^{T}D_{1})^{-1}Y, $$
where $I$ is the unit matrix. How to ensure $C^{0}$ continuity at joints for two splines given by the control points $(p_{1},...,p_{k})$ and $(q_{1}=p_{k},...,q_{l})$; see the Figure, please.
Neither the same end points $q_{1}=p_{k}$ nor the same derivatives $q_{1}^{\prime}=p_{k}^{\prime}$ are sufficient.
Thanks for your help.
Update:
Currently, I am able to solve the problem approximately, using the weight matrices $W_{1}$, $W_{2}$. The refined version of function $$ \phi=W_{1}\left[(X-X_{s})^{T}(X-X_{s})^{T}+(Y-Y_{s})^{T}(Y-Y_{s})\right]+W_{2}\left[X_{s}^{T}D_{1}^{T}D_{1}X_{s}+Y_{s}^{T}D_{1}^{T}D_{1}Y_{s}\right], $$ leads to the spline equations $$ X_{s}=(W_{1}+W_{2}D_{2}^{T}D_{2})^{-1}W_{1}X,\qquad Y_{s}=(W_{1}+W_{2}D_{2}^{T}D_{2})^{-1}W_{1}Y. $$ Initialize $W_{1}=W_{2}=I$ and set $W_{1}(4,4)=\infty$, $W_{2}(4,4)=0$ for the first spline, and $W_{1}(1,1)=\infty$, $W_{2}(1,1)=0$ for the second spline. For the results, see the next Figure.
Unfortunately, this is only an indirect solution of the problem...