As I understand it a Thin-Plate-Spline in 2D is an interpolant function that minimizes $$ \left\{ \begin{array}{ll} \int_{\mathbb{R}^2} \left(f_{xx}^2 + 2f_{xy}^2 + f_{yy}^2\right)dxdy & \\ \text{s.t.} f(x_i,y_i) = z_i & i=1,\ldots,n \end{array} \right. $$
However according to the link above the minimizer of such problem is equivalent to solve the biharmonic equation
$$ \Delta^2f = 0 $$
Which I don't really understand. If the problem was unconstrained I agree the two are equivalent (by using the Euler Lagrange equations).
However in case of constrained interpolation and assuming R.K.H.S. I think I need to minimize the following actually
$$ U(f,\vec\lambda) = \int_{\mathbb{R}^2} \left(f_{xx}^2 + 2f_{xy}^2 + f_{yy}^2\right)dxdy - \sum_{i=1}^n \lambda_i\left(z_i - L_i f\right) $$
where $L_i$ is the evaluation functional, which exists due to the R.K.H.S. hypothesis. However I don't know how to minimize such functional (and I assume it wouldn't just yield the biharmonic equation at the end).
Can anyone either give me a reference or if it is easy enough just show the derivation of the minimizer?
Update: Based on my reading Calculus of Variations - Gelfand & Fomin I've tried to formulate correctly my problem.
I am still assuming RKHS as setting (so I can use evaluation functional), although I see (including in one of the answer given) commonly weak solutions and distribution theory is used. Suppose I want to solve:
$$ \left\{ \begin{array}{ll} \int_{\mathbb{R}^2} \left(f_{xx}^2 + 2f_{xy}^2 + f_{yy}^2\right)dxdy & \\ \text{s.t.} f(x_i,y_i) = z_i & i=1,\ldots,n \end{array} \right. $$
Observe that for $i = 1, \ldots, n$ we have can rewrite each constraint as
$$ L_i f = z_i \text{for $i = 1,\ldots,n$} $$
where $L_i$ is the evaluation functional, which exists by my hypothesis of RKHS. Using Reisz representation theorem each constraint can be rewritten in theform
$$ L_i f = \int_{\mathbb{R^2}} \omega_i(x,y) f(x,y)dxdy = z_i \text{for $i = 1,\ldots,n$}. $$ I can rewrite the constrained optimization problem as
$$ \left\{ \begin{array}{ll} \int_{\mathbb{R}^2} \left(f_{xx}^2 + 2f_{xy}^2 + f_{yy}^2\right)dxdy & \\ \text{s.t.} \int_{\mathbb{R}^2} \omega_i(x,y)f(x,y)dxdy = z_i & i=1,\ldots,n \end{array} \right. $$
I am looking at this last formulation as a special case of the isoperimetric problem (see Calculus of Variations - Gelfand & Fomin) therefore using the lagrange multipliers therefore this constrained problem should be equivalent to solve
$$ \int_{\mathbb{R}^2} \left(f_{xx}^2 + 2f_{xy}^2 + f_{yy}^2 + \sum_{i=1}^n \lambda_i\omega_i(x,y)f(x,y) \right)dxdy $$
If now I did my calculations using the Euler Equations I get the following PDE to solve
$$ \Delta^2 f + \frac{1}{2} \sum_{i=1}^n \lambda_i \omega_i = 0 $$
which make sense to me now since in the distribution theory way those $\omega_i$ would be Dirac deltas probably.
Here's a reference to a nice write-up by David Eberly. You'll see the TPS to be the Green function of the biharmonic equation with a delta function in the RHS. See eq. (15) for 1d, and (21), (28) for nD cases. The paper also talks about the Euler Lagrange equations and the bending energy which is minimized by a linear combination of the Green functions from all the centers (29): https://www.geometrictools.com/Documentation/ThinPlateSplines.pdf (the link has been stable for many years)
I too hope for a more detailed treatment, so if you find a better reference please share or kindly write your own answer.