You have a convex quadratic form to minimize. This can be written in the form:

$\min f(x)=\frac{1}{2} x^{T}Ax - b^{T}x + c$

where $A$ is symmetric and positive semidefinite.

The gradient of $f$ is

$\nabla f(x)=Ax-b$

Any point $x$ where $\nabla f(x)=0$ will be a minimizer. In other words, any solution to $Ax=b$ will be a minimizer. If there is no solution, then the minimization problem is unbounded. Since $A$ is positive semidefinite, $A$ may be a singular matrix.

There are lots of ways of solving such a system of equations. Since $A$ is tiny (3 by 3), using a direct factorization method is most appropriate. Since $A$ is likely singular, a good choice is to use the QR factorization with column pivoting. This gives you a permutation matrix $P$, orthogonal matrix $Q$, and an upper triangular matrix $R$ such that

$AP=QR$

If $A$ is positive semidefinite but not positive definite, then the $R$ matrix will be singular, with one or more effectively 0 entries on the diagonal. You'll have to apply some tolerance in deciding what "effectively 0" is.

To solve $Ax=b$,

$APP^{-1}x=b$

$QR(P^{-1}x)=b$

$R(P^{-1}x)=Q^{-1}b$

Since $Q$ is orthogonal $Q^{-1}=Q^{T}$

$R(P^{-1}x)=Q^{T}b$

Let $y=P^{-1}x$ or $x=Py$. Solve

$Ry=Q^{T}b$

Since $R$ is upper triangular, it's easy to solve this sytem of equations by back substitution. If this system is inconsistent, then the original minimization problem is unbounded. If the system is consistent, then take any $y$ that solves the system of equations, and let $x=Py$. In solving the system of equations you may see some entries in $y$ that are free to take on any value. A simple solution is to set these to 0.

The QR factorization with column pivoting is implemented in MATLAB, Julia, Python/Numpy, etc. At a lower level, this is available in the LAPACK library that you can call from languages like C, Fortran, etc.