How to approximately solve system of linear equations that has no solutions

When I prepare data for Bayesian network meta analysis using model published in Woods B.S, Hawkins N et al - Network meta-analysis on the log-hazard scale, combining count and hazard ratio statistics accounting for multi-arm trials; A tutorial(2010), I need to solve a system of (usually trivial) simultaneous equations that might be over-specified, like the one below:

$\left\{ \begin{array}{ll} \text{hr}_1 = 0 \\ \text{hr}_2 - \text{hr}_1 & = \text{HR}_{2,\,1}\\ \text{hr}_3 - \text{hr}_1 & = \text{HR}_{3,\,1}\\ \text{hr}_3 - \text{hr}_2 & = \text{HR}_{3,\,2} \end{array} \right.$

where $\text{HR}_{i,\,j}$ is a published hazard rate, and $\text{hr}_k$ is latent cumulative hazard that I want to estimate.

What people usually do, is ignoring random equation(s) from the above set and solve the rest algebraically. But I think that a better solution would be to systematically (not randomly) produce all distinct and algebraically solvable truncated equation sets, and than average all the results. It is doable, of course, and I can program it in e.g. R. (I'd need to take all possible permutations of truncation of this equation, decide which is solvable, solve them, and average the results).

In the abundance of the packages already existing in R, maybe someone has already published functions that I might use to this end? The problem is, that I don't know how to find them. Maybe at least you can help my creativity and help me with the search phrase I could put to Google?

You can also think of your problem as follows. You want to find the $x$ that comes closest to solving $$Ax = b,$$ while noting that this system doesn't actually have any solutions. Now, one way to formalize this goal is to look for the $x$ that minimizes $$||Ax - b||_2^2.$$ Writing it like this, you just have a least squares problem. The solution for $x$ is given by $$x = (A^TA)^{-1}A^Tb.$$