What is it about the process of solving OLS that requires the parameters to be linear?
Because equations which are nonlinear in their parameters can't be written as $y=X\beta$. OLS estimates $\beta$ in the equation
$$
y = X\beta +\epsilon.
$$
This is a linear relationship, so when we say that $\hat{\beta} = (X^\top X)^{-1}X^\top y$ is the optimal estimator of $\beta$, what we mean is that it's optimal in the sense that it minimizes $\|y - X\beta\|_2^2$. Minimizing $\|y - X\beta\|_2^2$ is only important if this objective is meaningful for your task; particularly, if the task isn't linear in these parameters, then the fit may be poor.
However, one reason that OLS is so flexible is that if you can find a way to represent your data in a linear way, then it is linear in the parameters, otherwise known as basis expansion.
A textbook example of a change of basis is using a polynomial basis, so you have $X_\text{polynomial} = [1, x, x^2, x^3, \dots, x^p]$. The model $X_\text{polynomial}\beta$ is linear in its parameters, but viewed as a function of $x$, it's a nonlinear polynomial.
What would happen if they were non-linear and we tried to solve using OLS?
It won't work very well!
This data's deterministic component is given by $$
y = \beta_0 + \beta_1 \sin (\beta_2 x + \beta_3) $$
which is not linear in $\beta$, the parameter vector to be estimated, because you can't write this in the form $y=X\beta$. I also add small, independent 0-mean Gaussian noise to each observation.
If we do the naive thing and assume that our output $y$ is a linear function of $x$, then we find a poor fit, in the sense that there is a large discrepancy between the estimated line (red) and the true function (blue). The model finds that the best linear approximation is a decreasing line, completely ignoring the sinusoidal behavior.
One way to try to improve the fit is to re-express $x$. Since this looks like something sinusoidal, we might try a sine function. This gives the design matrix $X_\text{sine}=[1, \sin(x)]$. This give a flatter line, but it's still not a satisfying model. Even though the model and the desired function are both sine waves, we're implicitly using $\beta_0 + \beta_1 \sin(1 \times x + 0)$ to approximate $$
y = \beta_0 + \beta_1 \sin (\beta_2 x + \beta_3).$$ This is not a good approximation, because we've fixed $\beta_2=1$ and $\beta_3=0$, so the further the true values are from these assumed values, the poorer this approximation will be.
What we really need is a way to recover all of the parameters in the function $$
y = \beta_0 + \beta_1 \sin (\beta_2 x + \beta_3), $$ but this is a nonlinear estimation task, so we need to use the appropriate tools to accommodate the nonlinearity of the $\beta$s. Nonlinear least squares is one method to achieve this, among many others.
Code
set.seed(13)
N <- 1000
x <- runif(N, -pi, pi)
f <- function(x) pi + 2 * sin(4 * x)
y <- f(x) + rnorm(N,sd=0.5)
model <- lm(y ~ x)
png("~/Desktop/nonlinear.png")
plot(x,y,col="grey")
abline(model, col="red", lwd=2, lty="dashed")
lines(sort(x), f(sort(x)), lwd=2, col="blue")
dev.off()
model2 <- lm(y ~ sin(x) )
png("~/Desktop/nonlinear2.png")
plot(x,y,col="grey")
abline(model2, col="red", lwd=2, lty="dashed")
lines(sort(x), f(sort(x)), lwd=2, col="blue")
dev.off()
