$$
\Delta^2 f + \frac{1}{2} \sum_{i=1}^n \lambda_i \omega_i = 0
$$$$
\Delta^2 f + \frac{1}{2} \sum_{i=1}^n \lambda_i \omega_i = 0 \hspace{10mm} (1)
$$
which make sense to me now since in the distribution theory way those $\omega_i$ would be Dirac deltas probably.
Update:
I gave it more thought so I'll add an update (more for reference for people like me who are still reading on the subject). Equation $(1)$ can be rewritten as
$$
\Delta^2 f = - \frac{1}{2} \sum_{i=1}^n \lambda_i \omega_i
$$
And I observe that such equation is defined on $\mathbb{R}^d$ (so not a boundary value problem). I've been reading through fundamental solutions and Green function, from Evans - Partial Differential Equation. In my setting I believe the use of fundamental solution is more appropriate than the green function.
To find the fundamental solution we can solve the equivalent following system
$$
\left\{
\begin{array}{l}
\Delta g = -\frac{1}{2}\sum_{i=1}^n\lambda_i \omega_i \\
\Delta f = g
\end{array}
\right.
$$
So essentially we need to solve twice the laplace equation, non homogeneous.
In my reference (Evans) how to solve the laplace equation is explained and I think a very similar observation can be done here (i.e. the fundamental solution is a function) $\Phi(x) = \Phi(\left| x \right|)$ (so radial) and solution of the non homogenous equation can be obtained as
$$
f(x) = - \frac{1}{2} \int_{\mathbb{R}^d} \Phi(x - y) \left( \sum_{i=1}^n \lambda_i \omega_i(y) \right) dy = - \frac{1}{2} \sum_{i=1}^n \lambda_i\int_{\mathbb{R}^d} \Phi(x - y) \omega_i(y) dy. \hspace{10mm} (2)
$$
Remember that each $\omega_i$ is a representation of the evaluation functional $L_i$ therefore I can write $\int_{\mathbb{R}^d} \Phi(x - y)\omega_i(y) dy = \Phi(x - x_i)$. This leads to write the solution $f$ as
$$
f(x) = -\frac{1}{2} \sum_{i=1}^n \lambda_i \Phi(x - x_i).
$$
Observe that for $j = 1 ... n$ we can write
$$
f(x_j) = -\frac{1}{2} \sum_{i=1}^n \lambda_i \Phi(x_j - x_i) = \sum_{i=1}^n \gamma \Phi(x_j - x_i)
$$
with $\gamma_i = - \frac{\lambda_i}{2}$ and this leads to a square linear system of the form
$$
\begin{pmatrix}
\Phi(0) & \Phi(x_2 - x_1) & \dots & \Phi(x_n - x_1) \\
\Phi(x_1 - x_2) & \Phi(0) & \dots & \Phi(x_n - x_2) \\
\vdots & \vdots & \ddots & \vdots \\
\Phi(x_1 - x_n) & \Phi(x_2 - x_n) & \dots & \Phi(0)
\end{pmatrix}
\begin{pmatrix}
\gamma_1 \\
\gamma_2 \\
\vdots \\
\gamma_n
\end{pmatrix} =
\begin{pmatrix}
f(x_1) \\
f(x_2) \\
\vdots \\
f(x_n)
\end{pmatrix}
$$
At this point I am not very sure how to show the system is solvable. But I know it can be solved, and I think this is overall the theory of the TPS.
