Ajuste de dos curvas con regresión lineal / no lineal

Necesito ajustar dos curvas (que ambas deben pertenecer a funciones cúbicas) en un conjunto de puntos con JuMP.

He terminado de ajustar una curva, pero estoy luchando por ajustar 2 curvas en el mismo conjunto de datos.

Pensé que si podía distribuir puntos a las curvas, por lo que si cada punto solo se podía usar una vez, podría hacerlo como a continuación, pero no funcionó. (Sé que puedo usar cosas mucho más complicadas, quiero que sea simple).

Esta es una parte de mi código actual:

# cubicFunc is a two dimensional array which accepts cubicFunc[x,degree]

@variable(m, mult1[1:4]) // 0:3 because it's cubic
@variable(m, mult2[1:4]) // 0:3 because it's cubic

@variable(m, 0 <= includeIn1[1:numOfPoints] <= 1, Int)
@variable(m, 0 <= includeIn2[1:numOfPoints] <= 1, Int)

# some kind of hack to force one of them to 0 and other one to 1
@constraint(m, loop[i in 1:numOfPoints], includeIn1[i] + includeIn2[i] == 1)

@objective(m, Min, sum( (yPoints - cubicFunc*mult1).*includeIn1 .^2 ) + sum( (yPoints - cubicFunc*mult2).*includeIn2 .^2 ))

Pero da varios errores dependiendo de lo que estoy intentando; *includeIn1y, .*includeIn1no funciona, he intentado hacerlo @NLobjectivepero me dio ~ 50 líneas de errores, etc.

¿Es realista mi idea? ¿Puedo incluirlo en el código?

Cualquier ayuda será muy apreciada. Muchas gracias.

Puede escribir el problema, por ejemplo, así:

using JuMP, Ipopt

m = Model(with_optimizer(Ipopt.Optimizer))

@variable(m, mult1[1:4])
@variable(m, mult2[1:4])
@variable(m, 0 <= includeIn1[1:numOfPoints] <= 1)
@variable(m, 0 <= includeIn2[1:numOfPoints] <= 1)

@NLconstraint(m, loop[i in 1:numOfPoints], includeIn1[i] + includeIn2[i] == 1)

@NLobjective(m, Min, sum(includeIn1[i] * (yPoints[i] - sum(cubicFunc[i,j]*mult1[j] for j in 1:4)) ^2 for i in 1:numOfPoints) +
                     sum(includeIn2[i] * (yPoints[i] - sum(cubicFunc[i,j]*mult2[j] for j in 1:4)) ^2 for i in 1:numOfPoints))

optimize!(m)

Dadas las restricciones includeIn1y includeIn2será 1o 0en óptimo (si no lo son, esto significa que no importa a qué grupo asigne el punto), por lo que no tenemos que restringirlas para que sean binarias. También uso un solucionador no lineal ya que el problema no parece ser posible reformular como una tarea de optimización lineal o cuadrática.

Sin embargo, doy el código anterior solo como un ejemplo de cómo puedes escribirlo. La tarea que ha formulado no tiene un mínimo local único (que es global entonces), sino varios mínimos locales. Por lo tanto, el uso de solucionadores convexos no lineales estándar que admite JuMP solo encontrará un óptimo local (no necesariamente uno global). Para buscar óptimos globales, debe cambiar a solucionadores globales como, por ejemplo, https://github.com/robertfeldt/BlackBoxOptim.jl .