The algebraic (implicit) equation of an ellipsoid is
$ (r - C)^T Q (r - C) = 1 $
where $ r = [x, y, z]^T $ is a point on the ellipsoid surface, and $ C = [C_x, C_y, C_z]^T $ is the center of the ellipsoid. $Q$ is a $3 \times 3$ symmetric, positive definite matrix. In diagonalized form, $Q = R D R^T $ where $D$ is a diagonal matrix of the reciprocal of the square of the semi-axes of the ellipsoid, and $R$ is a rotation matrix that gives the orientation of the ellipsoid in space.
Define $X = \begin{bmatrix} r \\ 1 \end{bmatrix} $ then the equation of the ellipsoid becomes
$ X^T G X = 0 $
where
$ G = \begin{bmatrix} Q && - Q C \\ - C^T Q && C^T Q C - 1 \end{bmatrix} $
A tangent plane $ N^T X = 0 $ to the ellipsoid satisfies
$ N^T G^{-1} N = 0 $
So with the $9$ planes, solve
$0 = N_i^T G^{-1} N_i , \hspace{15pt} i = 1, 2, \dots, 9$
for the entries of $G^{-1}$ which is symmetric. This equation is linear in the entries of $G^{-1}$ and there are $10$ unknowns. Using Gauss-Jordan elimination, the entries of $G^{-1}$ can be solved as an arbitrary nonzero scalar times a constant vector. Setting the scalar to a variable $\dfrac{1}{t}$, we generate the matrix $G^{-1} = \dfrac{1}{t} G_0^{-1}$. Inverting, $G = t G_0$. Using the definition of $G$ we can solve for the scalar $t$ and thus find $Q$ and $C$.
We will have
$ G_0 = \begin{bmatrix} G_{011} && G_{012} \\ G_{012}^T && G_{022} \end{bmatrix} $
$ Q = t G_{011} $
$ Q C = - t G_{012} $
It follows that $C = - G_{011}^{-1} G_{012} $
And to determine $t$, consider the $(2, 2)$ scalar
$ t C^T G_{011} C - 1 = t G_{022} $
And this can be solve for $t$ in a straight forward manner. Now that we have the value of $t$, matrix $Q$ is known.
Finally, we diagonalize $Q$ to find the semi-axes and the orientation of the ellipsoid.
That's how I approached the problem of identification of the ellipsoid from the given $9$ planes tangent to it.
Any alternative suggestions and hints and solutions are greatly appreciated.
You are correct that the constraints from the tangent planes
$$ N_i^\top G^{-1} N_i = 0, \quad \forall\,i = 1, 2, \dots, 9 \tag{1} $$
are linear in $G^{-1}$. However, only using these constraints would also allow for the trivial solution $G^{-1} = 0$. It can be shown that
$$ G^{-1} = \begin{bmatrix} Q^{-1} - C\,C^\top & -C \\ -C^\top & -1 \end{bmatrix} \tag{2} $$
so the bottom right entry of $G^{-1}$ should be equal to minus one ($G^{-1}_{44} = -1$). Together with the six constraints from that $G$ and thus $G^{-1}$ is symmetric ($G^{-1}_{ij} = G^{-1}_{ji}, \forall\,i\neq j$) yields a total of 16 constraints, and thus uniquely defines $G^{-1}$.
Solving for $G^{-1}$ can be done in a structural way by writing all linear constraints in one system of linear equations. For this one can use the Kronecker product and vectorization, such that $(1)$ for each $i$ is equivalent to
$$ \left(N_i^\top \otimes N_i^\top\right) \text{vec}(G^{-1}) = 0. $$
Combining all these linear constraints into one set of linear equations yields
$$ \underbrace{ \begin{bmatrix} N_1^\top \otimes N_1^\top \\ \vdots \\ N_9^\top \otimes N_9^\top \\ \begin{matrix} 0 & 1 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{matrix} \end{bmatrix}}_{A\in\mathbb{R}^{16\times16}} \text{vec}(G^{-1}) = \underbrace{\begin{bmatrix} 0 \\ \vdots \\ 0 \\ -1 \end{bmatrix}}_{b}. \tag{3} $$
If each $N_i$ is chosen such that $A$ is nonsingular, then $(3)$ can be solved using
$$ \text{vec}(G^{-1}) = A^{-1} b. \tag{4} $$
If it is preferred to having to solve a smaller system of linear equations, it is also possible to use reduced coordinates $y\in\mathbb{R}^9$, such that when they are projected to $\text{vec}(G^{-1})$, for example using the following transformation
$$ \text{vec}(G^{-1}) = \underbrace{\begin{bmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}}_T y - \underbrace{\begin{bmatrix} 0 \\ \vdots \\ 0 \\ 1 \end{bmatrix}}_v, \tag{5} $$
the symmetry and $G^{-1}_{44} = -1$ constraints are always satisfied. Substituting $(5)$ into $(3)$ yields
$$ \underbrace{\begin{bmatrix} N_1^\top \otimes N_1^\top \\ \vdots \\ N_9^\top \otimes N_9^\top \end{bmatrix}\!T\,}_{\tilde{A}\in\mathbb{R}^{9\times9}}\,y = \underbrace{\begin{bmatrix} N_1^\top \otimes N_1^\top \\ \vdots \\ N_9^\top \otimes N_9^\top \end{bmatrix}\!v\,}_\tilde{b}. \tag{6} $$
And thus, similarly to $(4)$, $(6)$ can be solved for $y$ and thus $\text{vec}(G^{-1})$ using $(5)$, using
$$ \text{vec}(G^{-1}) = T\,\tilde{A}^{-1} \tilde{b} - v. \tag{7} $$
It can be noted that, besides $(7)$ being slightly less computationally expensive to solve than $(4)$, it would also be easier to expand to having more than 9 (but maybe noisy) tangent planes and obtain a least squares solution. Namely, the symmetry and $G^{-1}_{44} = -1$ constraints are hard constraints, but $(1)$ could have some uncertainty and thus doesn't have to be satisfied exactly. The least squares solution can then be obtained using
$$ \underbrace{\begin{bmatrix} N_1^\top \otimes N_1^\top \\ \vdots \\ N_M^\top \otimes N_M^\top \end{bmatrix}\!T\,}_{\hat{A}\in\mathbb{R}^{M\times9}}\,y = \underbrace{\begin{bmatrix} N_1^\top \otimes N_1^\top \\ \vdots \\ N_M^\top \otimes N_M^\top \end{bmatrix}\!v\,}_\hat{b}, \tag{8} $$
with $M>9$ the number of tangent planes, and
$$ \text{vec}(G^{-1}) = T\,\left(\hat{A}^\top \hat{A}\right)^{-1}\hat{A}^\top \hat{b} - v. \tag{9} $$