Separate tests for inequality and equality constraints

Author: Vaibhavdixit02

src/systems/optimization/optimizationsystem.jl

@@ -177,7 +177,12 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
                           expression = Val{false})
 
     obj_expr = toexpr(equations(sys))
-    pairs_arr = p isa SciMLBase.NullParameters ? [Symbol(_s) => Expr(:ref, :x, i) for (i, _s) in enumerate(dvs)] : [[Symbol(_s) => Expr(:ref, :x, i) for (i, _s) in enumerate(dvs)]..., [Symbol(_p) => p[i] for (i, _p) in enumerate(ps)]...]
+    pairs_arr = p isa SciMLBase.NullParameters ?
+                [Symbol(_s) => Expr(:ref, :x, i) for (i, _s) in enumerate(dvs)] :
+                [
+        [Symbol(_s) => Expr(:ref, :x, i) for (i, _s) in enumerate(dvs)]...,
+        [Symbol(_p) => p[i] for (i, _p) in enumerate(ps)]...,
+    ]
     rep_pars_vals!(obj_expr, pairs_arr)
     if grad
         grad_oop, grad_iip = generate_gradient(sys, checkbounds = checkbounds,
@@ -208,10 +213,11 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
 
     if length(sys.constraints) > 0
         @named cons_sys = NonlinearSystem(sys.constraints, dvs, ps)
-        cons = generate_function(cons_sys, checkbounds = checkbounds, linenumbers = linenumbers,
-                                                expression = Val{false})[1]
-        cons_j = generate_jacobian(cons_sys; expression=Val{false}, sparse=sparse)[2]
-        cons_h = generate_hessian(cons_sys; expression = Val{false}, sparse =sparse)[2]
+        cons = generate_function(cons_sys, checkbounds = checkbounds,
+                                 linenumbers = linenumbers,
+                                 expression = Val{false})[1]
+        cons_j = generate_jacobian(cons_sys; expression = Val{false}, sparse = sparse)[2]
+        cons_h = generate_hessian(cons_sys; expression = Val{false}, sparse = sparse)[2]
 
         cons_expr = toexpr(equations(cons_sys))
         rep_pars_vals!.(cons_expr, Ref(pairs_arr))
@@ -225,24 +231,24 @@ function DiffEqBase.OptimizationProblem{iip}(sys::OptimizationSystem, u0map,
         end
 
         _f = DiffEqBase.OptimizationFunction{iip}(f,
-                                              SciMLBase.NoAD();
-                                              grad = _grad,
-                                              hess = _hess,
-                                              hess_prototype = hess_prototype,
-                                              cons = cons,
-                                              cons_j = cons_j,
-                                              cons_h = cons_h,
-                                              cons_jac_prototype = cons_jac_prototype,
-                                              cons_hess_prototype = cons_hess_prototype,
-                                              expr = obj_expr,
-                                              cons_expr = cons_expr)
+                                                  SciMLBase.NoAD();
+                                                  grad = _grad,
+                                                  hess = _hess,
+                                                  hess_prototype = hess_prototype,
+                                                  cons = cons,
+                                                  cons_j = cons_j,
+                                                  cons_h = cons_h,
+                                                  cons_jac_prototype = cons_jac_prototype,
+                                                  cons_hess_prototype = cons_hess_prototype,
+                                                  expr = obj_expr,
+                                                  cons_expr = cons_expr)
     else
         _f = DiffEqBase.OptimizationFunction{iip}(f,
-                                              SciMLBase.NoAD();
-                                              grad = _grad,
-                                              hess = _hess,
-                                              hess_prototype = hess_prototype,
-                                              expr = obj_expr)
+                                                  SciMLBase.NoAD();
+                                                  grad = _grad,
+                                                  hess = _hess,
+                                                  hess_prototype = hess_prototype,
+                                                  expr = obj_expr)
     end
 
     OptimizationProblem{iip}(_f, u0, p; lb = lb, ub = ub, kwargs...)
@@ -315,16 +321,22 @@ function OptimizationProblemExpr{iip}(sys::OptimizationSystem, u0,
     ub = varmap_to_vars(ub, dvs; check = false, tofloat = false, use_union)
 
     obj_expr = toexpr(equations(sys))
-    pairs_arr = p isa SciMLBase.NullParameters ? [Symbol(_s) => Expr(:ref, :x, i) for (i, _s) in enumerate(dvs)] : [[Symbol(_s) => Expr(:ref, :x, i) for (i, _s) in enumerate(dvs)]..., [Symbol(_p) => p[i] for (i, _p) in enumerate(ps)]...]
+    pairs_arr = p isa SciMLBase.NullParameters ?
+                [Symbol(_s) => Expr(:ref, :x, i) for (i, _s) in enumerate(dvs)] :
+                [
+        [Symbol(_s) => Expr(:ref, :x, i) for (i, _s) in enumerate(dvs)]...,
+        [Symbol(_p) => p[i] for (i, _p) in enumerate(ps)]...,
+    ]
     rep_pars_vals!(obj_expr, pairs_arr)
 
     if length(sys.constraints) > 0
         @named cons_sys = NonlinearSystem(sys.constraints, dvs, ps)
-        cons = generate_function(cons_sys, checkbounds = checkbounds, linenumbers = linenumbers,
-                                                expression = Val{false})[1]
-        cons_j = generate_jacobian(cons_sys; expression=Val{false}, sparse=sparse)[2]
+        cons = generate_function(cons_sys, checkbounds = checkbounds,
+                                 linenumbers = linenumbers,
+                                 expression = Val{false})[1]
+        cons_j = generate_jacobian(cons_sys; expression = Val{false}, sparse = sparse)[2]
 
-        cons_h = generate_hessian(cons_sys; expression = Val{false}, sparse =sparse)[2]
+        cons_h = generate_hessian(cons_sys; expression = Val{false}, sparse = sparse)[2]
 
         cons_expr = toexpr(equations(cons_sys))
         rep_pars_vals!.(cons_expr, Ref(pairs_arr))

test/optimizationsystem.jl

@@ -1,4 +1,5 @@
-using ModelingToolkit, SparseArrays, Test, Optimization, OptimizationOptimJL, OptimizationMOI, Ipopt, AmplNLWriter, Ipopt_jll
+using ModelingToolkit, SparseArrays, Test, Optimization, OptimizationOptimJL,
+      OptimizationMOI, Ipopt, AmplNLWriter, Ipopt_jll
 
 @variables x y
 @parameters a b
@@ -49,14 +50,32 @@ sol = solve(prob, BFGS(initial_stepnorm = 0.0001), allow_f_increases = true)
 sol = solve(prob2, BFGS(initial_stepnorm = 0.0001), allow_f_increases = true)
 @test sol.minimum < -1e9
 
-prob = OptimizationProblem(sys2, [x => 0.0, y => 0.0], [a => 1.0, b => 100.0], lcons = [-1.0, -1.0], ucons = [500.0, 500.0], grad = true, hess = true)
+#inequality constraint
+prob = OptimizationProblem(sys2, [x => 0.0, y => 0.0], [a => 1.0, b => 100.0],
+                           lcons = [-1.0, -1.0], ucons = [500.0, 500.0], grad = true,
+                           hess = true)
 sol = solve(prob, IPNewton(), allow_f_increases = true)
 @test sol.minimum < 1.0
 sol = solve(prob, Ipopt.Optimizer())
 @test sol.minimum < 1.0
 sol = solve(prob, AmplNLWriter.Optimizer(Ipopt_jll.amplexe))
 @test sol.minimum < 1.0
 
+#equality constraint
+cons2 = [0.0 ~ x^2 + y^2]
+sys2 = OptimizationSystem(loss, [x, y], [a, b], name = :sys2, constraints = cons2)
+prob = OptimizationProblem(sys2, [x => 0.0, y => 0.0], [a => 1.0, b => 1.0], lcons = [1.0],
+                           ucons = [1.0], grad = true, hess = true)
+sol = solve(prob, IPNewton())
+@test sol.minimum < 1.0
+@test prob.f.cons(sol.minimizer, [1.0, 1.0]) ≈ [1.0]
+sol = solve(prob, Ipopt.Optimizer())
+@test sol.minimum < 1.0
+@test prob.f.cons(sol.minimizer, [1.0, 1.0]) ≈ [1.0]
+sol = solve(prob, AmplNLWriter.Optimizer(Ipopt_jll.amplexe))
+@test sol.minimum < 1.0
+@test prob.f.cons(sol.minimizer, [1.0, 1.0]) ≈ [1.0]
+
 rosenbrock(x, p) = (p[1] - x[1])^2 + p[2] * (x[2] - x[1]^2)^2
 x0 = zeros(2)
 _p = [1.0, 100.0]