Fix a bug where nonlinear objective terms are not added correctly

lib/nleval.cpp

@@ -690,7 +690,6 @@ void NonlinearFunctionEvaluator::analyze_hessian_structure(
 void NonlinearFunctionEvaluator::eval_objective(const double *x, double *y)
 {
 	double obj = 0.0;
-	double temp = 0.0;
 
 	double *p = this->p.data();
 
@@ -708,18 +707,17 @@ void NonlinearFunctionEvaluator::eval_objective(const double *x, double *y)
 			{
 				auto &x_indices = inst.xs;
 				auto &p_indices = inst.ps;
-				evaluator.f_eval.p(x, p, &temp, x_indices.data(), p_indices.data());
+				evaluator.f_eval.p(x, p, &obj, x_indices.data(), p_indices.data());
 			}
 		}
 		else
 		{
 			for (const auto &inst : inst_vec)
 			{
 				auto &x_indices = inst.xs;
-				evaluator.f_eval.nop(x, &temp, x_indices.data());
+				evaluator.f_eval.nop(x, &obj, x_indices.data());
 			}
 		}
-		obj += temp;
 	}
 
 	y[0] += obj;

src/pyoptinterface/_src/ipopt.py

@@ -98,6 +98,7 @@ def compile_functions_c(model: "Model", jit_compiler: TCCJITCompiler):
             np=np,
             indirect_x=True,
             indirect_p=True,
+            add_y=True,
         )
         if autodiff_structure.has_jacobian:
             jacobian_name = name + "_jacobian"
@@ -193,6 +194,7 @@ def compile_functions_llvm(model: "Model", jit_compiler: LLJITCompiler):
             np=np,
             indirect_x=True,
             indirect_p=True,
+            add_y=True,
         )
         export_functions.append(f_name)
         if autodiff_structure.has_jacobian:

tests/test_nlp_clnlbeam.py

@@ -0,0 +1,47 @@
+import pyoptinterface as poi
+from pyoptinterface import nlfunc
+
+import pytest
+
+
+def test_clnlbeam(ipopt_model_ctor):
+    model = ipopt_model_ctor()
+
+    N = 1000
+    h = 1 / N
+    alpha = 350
+
+    t = model.add_variables(range(N + 1), lb=-1.0, ub=1.0)
+    x = model.add_variables(range(N + 1), lb=-0.05, ub=0.05)
+    u = model.add_variables(range(N + 1))
+
+    def obj(vars):
+        return 0.5 * h * (vars.u2**2 + vars.u1**2) + 0.5 * alpha * h * (
+            nlfunc.cos(vars.t2) + nlfunc.cos(vars.t1)
+        )
+
+    obj_f = model.register_function(obj)
+    for i in range(N):
+        model.add_nl_objective(
+            obj_f, nlfunc.Vars(t1=t[i], t2=t[i + 1], u1=u[i], u2=u[i + 1])
+        )
+
+    def con(vars):
+        return vars.x2 - vars.x1 - 0.5 * h * (nlfunc.sin(vars.t2) + nlfunc.sin(vars.t1))
+
+    con_f = model.register_function(con)
+    for i in range(N):
+        model.add_nl_constraint(
+            con_f, nlfunc.Vars(t1=t[i], t2=t[i + 1], x1=x[i], x2=x[i + 1]), eq=0.0
+        )
+
+    for i in range(N):
+        model.add_linear_constraint(
+            t[i + 1] - t[i] - 0.5 * h * u[i + 1] - 0.5 * h * u[i], poi.Eq, 0.0
+        )
+
+    model.optimize()
+
+    objective_value = model.get_model_attribute(poi.ModelAttribute.ObjectiveValue)
+
+    assert objective_value == pytest.approx(350.0, abs=1e-8)

tests/test_nlp_rocket.py

@@ -73,7 +73,7 @@ def dynamics_eq(vars):
                 T2=T[i + 1],
                 step=step,
             ),
-            eq=[0.0, 0.0, 0.0],
+            eq=0.0,
         )
 
     # Boundary conditions