If we assume that each of the floating point values such as 12.28 is to be considered precise (such as if it had been written as 1228/100), then those equations are inconsistent. You can solve for three of the variables and the 4th will drop out leaving a non-zero value in a place that must be 0 for the equations to be consistent. It is not, in other words, a matter that the four variables could be rewritten as three consistent variables: the equation cannot be made to work in the form shown.
Every time you have floating point values in equations, you need to stop and question why you are using solve(). solve() is for finding exact algebraic solutions (if they exist), but floating point numbers are inherently imprecise. Is 12.28 intended to be 1228/100 or is it intended to represent some value that is in the range 1228/100 +/- 5/1000 ? And now, having that doubt about what the 12.28 really means, does the "2" of "2*r1" really mean 2 exactly or does it mean some number near 2? For quadratics the solutions are usually not all that sensitive to little fluctuations, but by the time you get to quartics the solutions are very sensitive to little fluctuations in the inputs. And if you have two quadratics together you pretty much have to assume you have the same sensitivity as if you were dealing with a quartic.
I went through your equations, varying term by term under the assumption that the term had not been given exactly, such a supposing that 12.28 might be (1228/100 + delta) for some unknown delta, or that 0*x*rr might be (0+delta)*x*rr for some unknown delta. I would alter a term, and solved to find out what the "noise" would have to be in the term in order to make the equations consistent. For example we could suggest that maybe 12.28 might result in inconsistency but maybe 12.28003582 might allow the equations to be consistent. I did that for each possibility in term, finding out how wrong the stated term would have to be in order for there to be a solution to the equations. It turns out that with those equations, in every case except perhaps one, the amount that the term would have to be wrong was fairly large compared to what the term actually is, such as -0.05*x^2 needing to be about +0.28*x^2 or -6994.94 needing to be about -15000 for there to be a solution. Most of the possibilities just stayed inherently inconsistent. Some of the possibilities led to a consistent series of equations whose solutions depended upon the magnitude of the noise (which I called "delta" when I was solving), but in every one of those cases, the consistent solutions involved expressions divided by delta, so the solutions "exploded" as the assumed noise was reduced, contradicting the assumption that there was just a little noise in the system.
The possibility that involved the least change to the parameters was that in the "i" equation, where you currently have
i=-0.05*x^2-25.82*x+r1*x+0*rr*x+271.28*r1+0*rr+B*bd-135.64*B-6994.94
if the 135.64 there was instead about 135.21, the equations could be made consistent. But this change to 135.64 only applied to that one line: if you made the same change to the 135.64 in the g and h functions, then the equations blow up again. And also when I say consistent, I mean that they could be written in terms of three equations in three variables without the 4th calculation equating two different non-zero values. In all of the cases where if you tweak a coefficient to create a consistent series of 4 equations in 4 unknowns, the solutions explode as the assumed noise is brought to zero.
It is, by the way, changes to the 135.64 on the various lines that you are most likely to be able to create a consistent series of equations.
