Estimate parameter in $x(t) = A \cdot \sin{\left(\omega \cdot t + \varphi \right)} + A_0$ through least squares method I have a bunch of data $(x_i, \, t_i)$ that is expected to distribute like $x(t) = A \cdot \sin{\left(\omega \cdot t + \varphi \right)} + A_0$.
Now I have to estimate the four parameters $(A, \, \omega, \, \varphi, \, A_0)$. So I want to minimize the formula (variances on $x_i$ are all the same, and the variances on $t$ are irrelevant): $$ f(A, \, \omega, \, \varphi, \, A_0) = \frac{1}{\sigma^2} \cdot \sum_{i=i}^n (x_i - A \cdot \sin{\left(\omega \cdot t + \varphi \right)} - A_0)^2$$
And then I have to solve the system
$$\begin{cases}
\frac{\partial f}{\partial A} = 0\\ 
\frac{\partial f}{\partial \omega} = 0 \\ 
\frac{\partial f}{\partial \varphi} = 0 \\ 
\frac{\partial f}{\partial A_0} = 0
\end{cases}
$$
And get:
$$\begin{cases}
\frac{1}{\sigma^2} \cdot \sum_{i=i}^n -2 \sin(\varphi + \omega t_i) (-A \sin(\varphi + \omega t_i) - A_0 + x_i) = 0 \\ 
\frac{1}{\sigma^2} \cdot \sum_{i=i}^n -2 A t_i \cos(\varphi + \omega t_i) (-A \sin(\varphi + \omega t_i) - A_0 + x_i) = 0 \\ 
\frac{1}{\sigma^2} \cdot \sum_{i=i}^n -2 A \cos(\varphi + \omega t_i) (-A \sin(\varphi + \omega t_i) - A_0 + x_i) = 0 \\ 
\frac{1}{\sigma^2} \cdot \sum_{i=i}^n -2 (-A \sin(\varphi + \omega t_i) - A_0 + x_i) = 0
\end{cases}
$$
But, at this point, I can't go ahead. I don't know how to solve the system because I couldn't separate parameters out of sum.
Has anybody any tips to solve the system? Or an alternative method to obtain the parameters? Any software for such minimization I can use?

An algebraically equivalent model (with the random error term $\epsilon$ explicitly added) is
$$\begin{aligned}
x(t) &= A_0 + A\left[\cos(\phi)\sin(\omega \, t) + \sin(\phi)\cos(\omega\, t)\right] + \epsilon\\
&= \beta_0 + \beta_1\sin(\omega \, t) + \beta_2 \cos(\omega\, t) + \epsilon
\end{aligned}$$
where
$$\begin{aligned}
A_0 &= \beta_0\\
A &= \sqrt{\beta_1^2 + \beta_2^2}\\
\phi &= \arg(\beta_1,\beta_2) = \operatorname{Atan2}(\beta_1,\beta_2).
\end{aligned}$$
By using the principal arctangent, the ambiguity in the original formulation is eliminated: $A$ will be non-negative and $\phi$ will be in the interval $(-\pi,\pi].$
Given any fixed value of $\omega,$ this is an ordinary least squares model for the remaining parameters $\beta_0,\beta_1,\beta_2$ and can be efficiently solved. Let the estimates associated with any such $\omega$ be $(A_0(\omega), A(\omega), \phi(\omega))$ and let the sum of squared residuals associated with those estimates be $\sigma^2(\omega).$  This reduces the problem to

Minimize $\sigma^2(\omega)$ given $\omega$ lies within a reasonable range of frequencies.

This is a one dimensional problem which needs to be solved numerically.  Confidence intervals, post hoc tests, and so on, can be conducted in the usual ways (by making distributional assumptions about $\epsilon$ and applying Maximum Likelihood, for instance).

Some readers with access to nonlinear modeling functions might think this analysis is unnecessary.  However, fitting this model correctly is tricky: $\sigma^2(\omega)$ typically has many local minima.  Even when you can supply reasonable starting values for the search it still can go awry.  Thus, it's an advantage to have an independent solution.  Having found it, you can use its estimates as starting values for your black-box solver in order to exploit the solver's added capabilities.  (The R function nls, for instance, supplies enough information for hypothesis testing.)
As an example, I generated the dataset of values plotted by the orange points below by randomly varying the values of the gray ("True") curve.  The dotted black curve was fit using the method described above (using the nlm function in R to estimate $\omega$ and its ls function for the OLS calculations).  It is a reasonable approximation to the truth.  The blue curve is the nls solution based on reasonable starting values for the parameters.  Its phase and frequency are very wrong.

To reproduce this figure, uncomment the #set.seed(17) line in the code below.  To develop some intuition for this problem, leave the line commented and play with the parameters at the beginning to vary the circumstances and study the solutions.
#
# Specify the (true) parameters.
#
A0 <- 5
A <- 2
phi <- 3*pi/4
omega <- -24
sigma <- 1  # SD of the error terms
n <- 87     # Data set size
#
# Formulate the model.
#
f <- function(t, A0, A, phi, omega) {
  A0 + A * sin(omega * t + phi)
}
#
# Specify the data.
#
X <- data.frame(t = cumsum(rexp(n, n)))
# set.seed(17)
X$x <- f(X$t, A0, A, phi, omega) + rnorm(n, 0, sigma)
#
# Plot the data and the data-generation model.
#
with(X, plot(t, x, pch=21, bg="Orange"))
curve(f(t, A0, A, phi, omega), add=TRUE, lwd=2, col="Gray", xname="t")
#
# Define the function to optimize.
#
ols <- function(omega, data) {
  fit <- lm(x ~ sin(omega * t) + cos(omega * t), data)
  list(SS=sum(residuals(fit)^2), b=coefficients(fit))
}
h <- function(omega, data) {
   ols(omega, data)$SS
}
#
# Estimate omega.
# It is crucial not to make the upper limit too large: otherwise, the 
# solution frequency will be too great.
#
# Need a reasonable starting value.  Estimate it with FFT:
#
X <- X[order(X$t), ]
y <- with(X, approx(t, x, seq(min(t), max(t), length.out=nrow(X)))$y)
omega.0 <- which.max(Mod(fft(y)[seq_len(n/2)][-1])) * 2 * pi * diff(range(X$t))

obj <- optim(omega.0, h, data=X, method="Brent",
      lower=omega.0/3, upper=omega.0*2)
omega.hat <- obj$par
#
# Compute A0, A, phi from b.
#
b <- ols(omega.hat, X)$b
A0.hat <- b["(Intercept)"]
A.hat <- sqrt(sum(b[c("sin(omega * t)", "cos(omega * t)")]^2))
phi.hat <- atan2(b["cos(omega * t)"], b["sin(omega * t)"])
#
# Plot the estimated curve.
#
curve(f(t, A0.hat, A.hat, phi.hat, omega.hat), add=TRUE, lwd=2, lty=3, xname="t")
#
# Independently apply nls as a check.
#
A0.0 <- mean(X$x)
A.0 <- diff(range(X$x))/3
phi.0 <- 0 # This is difficult to guess, so start anywhere
#
# Alternatively: Polish the nlm estimates with nls.
#
# A0.0 <- A0.hat
# A.0 <- A.hat
# phi.0 <- phi.hat
# omega.0 <- omega.hat

obj2 <- nls(x ~ A0 + A * (sin(omega * t + phi)), data=X, 
            start=list(A0=A0.0, A=A.0, phi=phi.0, omega=omega.0),
            control=list(maxiter=200))
a <- as.list(coefficients(obj2))
blue <- "#4080ff80"
curve(f(t, a$A0, a$A, a$phi, a$omega), add=TRUE, lwd=2, lty=1, col=blue, xname="t")
#
# Compare the estimates.
#
ss <-  ols(omega.hat, X)$SS
ss2 <- sum(residuals(obj2)^2)
if (ss < 0.999 * ss2) warning("nls is worse.")
#
# Plot decorations.
#
legend("top", bty="n", cex=0.9, title="Key:",
       legend=c("True", "nlm", "nls"),
       lty=c(1, 3, 1),
       col=c("Gray", "Black", blue),
       lwd=2,
       bg=NA)

