I am given $N$ plane equations in the $ax+by+cz+d=0$ form, meaning we have their normals $a, b, c: ||(a, b, c)|| = 1$ and distances from the origin ($d$). In general, the planes are similar, so the coefficients do not differ very much. I do not have any points on the given planes, but I guess they can be sampled.
I need to find an equation of the plane which is closest to them all, i.e represents the cluster of the planes in the best way. To put it differently, I need to merge all the planes into one. I am not sure what metric to use to measure plane fitness.
I've tried taking average of all plane coefficients, which produces a good result of the normal, but the last coefficient ($d$) seems wrong, the resulting plane has a small offset from the cluster.
It loos like a minimization problem subjected to constraints. Generate $n$ points from each of the $N$ planes, now you'll have a dataset of $\{(z_i, x_i, y_i)\}_{i=1}^{Nn}$ points in $\mathbb{R}^3$, normalizing the constraint w.r.t to $c\neq0 $, (assume, for now that $c$ is known or estimable apriori) $$||(a,b,c)||=1/c\sqrt{a^2/c^2 + b^2/c^2 +1} =1/c\sqrt{\beta_1^2 + \beta_2^2 +1} =1 \to g(\beta_1, \beta_2)=c.$$
So using Lagrange multiplier, your problem is $$ \min_{\beta \in \mathbb{R}^3} \frac{1}{nN}\sum_{i=1}^{nN}(z_i - \beta_0 - \beta_1x_i - \beta_2y_i)^2 - \lambda (g(\beta_1, \beta_2) - c). $$ It looks like a variation of LASSO regression. So take a look at LASSO and Ridge regressions. It may help.