Should we always minimize squared deviations if we want to find the dependency of mean on features? It is known that for any distribution a minimization of squared deviations from a value $c$ gives mean of the distribution. In other words, if we generate many values $y_i$ using our distribution, and then try to find $c$ that minimizes the following expression:
$
\lim_{N \to \infty}  \frac{1}{N}\sum_{i}^N (y_i - c)^2
$
then we find that the minimum is achieved when $c = \mu$, where $\mu$ is the mean of the given distribution. And for any $N$ the above minimization gives the sample mean.

Now, let's assume tat we have a distribution which depends on a single real-valued feature $x$, and we have a data-set generated by this distribution: ($x_1$, $y_1$), ($x_2$, $y_2$), $\dots$, ($x_n$, $y_n$). Then we try to find a function that gives dependency of the mean of this distribution on the feature:
$
\mu = \mu(x, p_1, p_2, \dots, p_k)
$
where $p_i$ are model parameters that we want to find.
Obviously, that we can get the desired function (values of the parameters) if we choose them by minimization of squared deviations:
$
\frac{1}{N} \sum_{i}^{N} [y_i - \mu(x_i, p_1, p_2, \dots, p_k)]^2
$
In other words, we can find the correct answer (correct values of the parameters $p_i$) under condition that $N$ goes to infinity and that the real function, that we are searching, is covered by the space of functions that we consider.
To summarize, the minimization of squared deviations can give not only mean of a given constant distribution, it can also gives correct dependency of the mean on a feature.

In reality $N$ is never infinity and can be rather small. As a result minimization of squared deviation will not give exactly correct values of the model parameters $p_i$. For another sample of the same size we will give other estimations of the parameters $p_i$. In general case, we can say that our estimations of the parameters are somehow distributed.
If we, instead of squared deviations minimize some other measure, we will get other distributions of our model parameters $p_i$. So, different measure of accuracy give different distributions in the parameters space for a given fixed sample size.
These distributions are characterized by their "width" and possible a "shift" relative to the correct point (true values of the model parameters).
Can it be the case that a use of something different from squared deviation gives a better estimate of model parameters (for example more accurate (smaller width) and with smaller or no systematic from the correct model parameters?

Your estimator is the OLS estimator in nonlinear regression
Your problem is essentially just the OLS estimation problem in nonlinear regression.  To see this, suppose you have nonlinear regression model of the form:
$$Y_i = \mu(x_i, \mathbf{p}) + \varepsilon_i
\quad \quad \quad \quad \quad 
\varepsilon_i \sim \text{IID Dist}(\text{Mean} = 0),$$
with error terms that are independent of the explanatory variables.$^\dagger$  Taking conditional expectations of both sides yields the equation:
$$\mu(x_i, \mathbf{p}) = \mathbb{E}(Y_i|X_i=x_i, \mathbf{p}).$$
This shows that the function $\mu$ is the true regression function which represents the conditional expected value of the response variable conditional on the explanatory variable.  Now, if you have $n$ observations from this model, then the OLS estimator for $\mathbf{p}$ is defined by the optimisation requirement:
$$\mu(x_i, \hat{\mathbf{p}}_n) = \min_\mathbf{p} \sum_{i=1}^n (y_i-\mu(x_i, \mathbf{p}))^2,$$
which is the optimisation requirement in your question.  So, your problem is about OLS estimation of the parameter $\mathbf{p}$, and you want to know the consistency properties of the OLS estimator $\hat{\mathbf{p}}_n$.
Consisteny properties in regression models is something that has been examined extensively in the statistical literature.  Under broad conditions on the sequence of explanatory variables, the estimator $\hat{\mathbf{p}}_n$ is a consistent estimator for $\mathbf{p}$.  For OLS estimation in nonlinear regression, these conditions are extensions of the "Grenander conditions" for consistency in linear regression (see e.g., Richardson and Bhattacharyya 1990).  The exact conditions are quite technical, but heuristically, they require the sequence of explanatory variables $x_1,x_2,x_3,...$ to be such that the "influence" of any finite set of data points tends to zero as $n \rightarrow \infty$.
If you have the mathematical background to do so, I recommend you read the linked paper to get an understanding of the consistency conditions in nonlinear regression.  For a simpler place to get started you can have a look at some answers on this site that look at the Grenander conditions for OLS consistency in linear regression (see e.g., here).

$^\dagger$ In this formulation of the nonlinear regression model I allow any error distribution with zero mean.  In the special case where the error distribution is normal the OLS estimator will correspond to the MLE.  The general case is used here because you have not specified an error distribution in your problem.
