As a first tool of analysis, analyze the multicollinearity in your design matrix , i.e the X matrix of your linear model y=Xb+e, compute the Condition number of X (condition number =sigmamax/sigmamin), sigma corresponds to the singular values of X matrix, in general condition number of 1000 or more is a case of strong multicollinearity , in which case you have to choose an alternate experimental design for data collection.

secondly, do a canonical analysis, plot the response surface and analyze the nature of the stationary point, the coefficients of the quadratic terms appear negative, so the Hessian matrix is most probably negative definite/Negative semi-definite or indefinite, solve for the eigenvalues of the associated Hessian matrix.

thirdly, over the chosen experimental range of your data collection, do experiments at one or more factor settings not originally a part of the initial data collection design, compare these experimental values with the value predicted by your regression model, and observe whether the experimental value lies within the confidence limits(preferable 95%) of the predicted mean yield at the corresponding factor setting(s).

You can always look at my work in this regard. However, if you want more information please do not hesitate to write

When all the independent variables equal zero yield will be -51.9. Since this result is a little silly, but remember that the equation is only good within the range of your data. If you lowest rate of nitrogen is 50 kilos/hectare, then estimating the effect at zero kilos/ha might not work so well.

Try getting a 95% confidence interval about each coefficient.

Make sure that the model assumptions are satisfied.

Do you have X1 through X32 plus all pairwise interaction terms? That is a large number of variables to work with. Alternatively, there is a theoretical basis for these select variables and the others that could be present are not relevant to the problem.