Using a local optimizer is not sufficient for global optimization.
Naïvely using a local optimizer (such as Levenberg–Marquardt) to find the parameters of a non-identified model is not a good idea because you have no assurance that the solution achieves a global minimum. Some models (like neural networks) simply ignore this fact, but they do so at the cost of never being able to interpret the coefficients. It's clear from the question and comment thread that you really care about this specific parameterization of the model, so you want to estimate $a,b,c,\lambda$ and not some alternative parameterization. So we should use a method that will do a good job of estimating them.
The ordinary logistic regression model is not sufficient because the model is not identified.
Let's define $X$ as the polynomial design matrix:
$$
X = \begin{bmatrix}
1 & x_1^2 & x_1 x_2 & x_2^2 & x_1 & x_2
\end{bmatrix}
$$
where $x_i$ is a column of the original data.
The standard logistic regression estimates $\hat y = \mathbb{P}(y=1)$ as a linear function on the logit scale, and estimates a single, distinct scalar coefficient $\beta_i$ to each column of $X$.
$$
\text{logit}(\mathbb{P}(y=1)) = X\beta
$$
Maximizing the likelihood (equiv. minimizing the cross entropy) has a unique optimal solution as long as $X$ is full rank and the data do not exhibit perfect separation. Moreover, the minimization of the loss is convex (equiv. maximization of the likelihood is concave).
So the relationship between the linear-in-parameters model and the nonlinear model is
$$
\begin{bmatrix}
\beta_0 \\
\beta_1 \\
\beta_2 \\
\beta_3 \\
\beta_4 \\
\beta_5
\end{bmatrix} = \begin{bmatrix}
c \\
a \lambda^2 \\
2 a \lambda (1 - \lambda) \\
a (1 - \lambda)^2 \\
b \lambda \\
b (1 - \lambda)
\end{bmatrix}.
$$
If we stop here, and just estimate the linear-in-parameters model, we may not recover $a,b$, which are the desired quantities.
Enforcing some constraints such as $a > 0$ and $b > 0$ and $ 0 < \lambda < 1$ will guarantee that we can write down a system of equations that is linear in $\log \beta_i$:
$$
\begin{bmatrix}
\log \beta_1 \\
\log \beta_2 \\
\log \beta_3 \\
\log \beta_4 \\
\log \beta_5 \\
\end{bmatrix}
=
\begin{bmatrix}
0 & 1 & 2 & 0 & 0 \\
1 & 1 & 1 & 1 & 0 \\
0 & 1 & 0 & 2 & 0 \\
0 & 0 & 1 & 0 & 1 \\
0 & 0 & 0 & 1 & 1 \\
\end{bmatrix}
\begin{bmatrix}
\log 2 \\
\log a \\
\log \lambda \\
\log 1 - \lambda \\
\log b
\end{bmatrix}
$$
However, this has its own problems:
- (i) we would need specialized software to enforce the constraints and
- (ii) there is no guarantee that the system is consistent.
We could fix (ii) by finding a least-error estimate of $a,b$, but would still require the constrained optimization in (i). In any event, the constraints do not appear in OP's question, so the the constrained model may not even correspond to a valid solution to OP's problem.
For fixed $\lambda$, the model is identified.
Instead, we can rewrite the proposed model to avoid these complications, at the cost of a slightly more involved estimation procedure:
$$
X
\begin{bmatrix}
c \\
a \lambda^2 \\
2 a \lambda (1 - \lambda) \\
a (1 - \lambda)^2 \\
b \lambda \\
b (1 - \lambda)
\end{bmatrix} = X \Lambda \begin{bmatrix}c \\ a \\ a \\ a \\ b \\ b \end{bmatrix}
$$
where $\Lambda$ is a diagonal matrix $$\Lambda = \begin{bmatrix}
1 & & & & & \\
& \lambda^2 & & & & \\
& &2 \lambda (1 - \lambda) & & & \\
& & &(1 - \lambda)^2 & & \\
& & & &\lambda & \\
& & & & &(1 - \lambda)
\end{bmatrix}.$$
The proposed model is not linear in its parameters, so it cannot be solved using standard logistic regression software. (The parameters are not affine.)
However, this factorization gives us a significant hint: conditional on fixing $\Lambda$, the model is linear in its parameters. So we can take $X\Lambda$ as the design matrix, and then either
- use specialized optimization software that enforces $a=\beta_1 = \beta_2 = \beta_3$ and $b=\beta_4 = \beta_5$; or
- use a standard logistic regression software to estimate $\beta$ and then work backward to the estimates of $a,b$ by minimizing error in the linear system.
It's up to you to choose the approach that you prefer.
All of this rests on choosing $\Lambda$. The matrix $\Lambda$ has only 1 parameter. A straightforward approach is to wrap your chosen logistic regression procedure inside of a one-dimensional optimizer (e.g. bisection search), and choosing the result that achieves the best fit according to whatever criterion you desire. So you'll be testing a sequence of $\Lambda^{(k)}$ matrices, and estimating coefficients $\beta^{(k)}$ (or $a^{(k)}, b^{(k)}$) for each one. The final model is the $\Lambda^{(k)}$ and $\beta^{(k)}$ (or $a^{(k)}, b^{(k)}$) that has the best fit.
This approach is broadly similar to Do statisticians assume one can't over-water a plant, or am I just using the wrong search terms for curvilinear regression?
