Add hs85 and hs89 problems

Author: arnavk23

Project.toml

@@ -10,13 +10,15 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 DataFrames = "1"
 JLD2 = "0.5, 0.6"
 JuMP = "^1.15"
 Requires = "1"
 SpecialFunctions = "2"
+StaticArrays = "1.9.18"
 julia = "1.6"
 
 [extras]

src/ADNLPProblems/hs85.jl

@@ -0,0 +1,164 @@
+export hs85
+
+function hs85(; type::Type{T}=Float64, kwargs...) where {T}
+    a = T[0, 17.505, 11.275, 214.228, 7.458, 0.961, 1.612, 0.146, 107.99,
+          922.693, 926.832, 18.766, 1072.163, 8961.448, 0.063, 71084.33, 2802713]
+    b = T[0, 1053.6667, 35.03, 665.585, 584.463, 265.916, 7.046, 0.222,
+          273.366, 1286.105, 1444.046, 537.141, 3247.039, 26844.086, 0.386,
+          140000, 12146108]
+    c10 = T(12.3) / T(752.3)
+
+    # Variable bounds and starting point (exact from the standard model)
+    lvar = T[704.4148, 68.6, 0.0, 193.0, 25.0]
+    uvar = T[906.3855, 288.88, 134.75, 287.0966, 84.1988]
+    x0   = T[900.0, 80.0, 115.0, 267.0, 27.0]
+
+    # Best known value ≈ -1.90513375
+    function f(x::AbstractVector{T})
+        # All intermediates (identical to those used in constraints)
+        y1  = x[2] + x[3] + T(41.6)
+        c1  = T(0.024) * x[4] - T(4.62)
+        y2  = T(12.5) / c1 + T(12)
+        c2  = T(0.0003535) * x[1]^2 + T(0.5311) * x[1] + T(0.08705) * y2 * x[1]
+        c3  = T(0.052) * x[1] + T(78) + T(0.002377) * y2 * x[1]
+        y3  = c2 / c3
+        y4  = T(19) * y3
+        c4  = T(0.04782) * (x[1] - y3) + T(0.1956) * (x[1] - y3)^2 / x[2] +
+              T(0.6376) * y4 + T(1.594) * y3
+        c5  = T(100) * x[2]
+        c6  = x[1] - y3 - y4
+        c7  = T(0.95) - c4 / c5
+        y5  = c6 * c7
+        y6  = x[1] - y5 - y4 - y3
+        c8  = (y5 + y4) * T(0.995)
+        y7  = c8 / y1
+        y8  = c8 / T(3798)
+        c9  = y7 - T(0.0663) * y7 / y8 - T(0.3153)
+        y9  = T(96.82) / c9 + T(0.321) * y1
+        y10 = T(1.29) * y5 + T(1.258) * y4 + T(2.29) * y3 + T(1.71) * y6
+        y11 = T(1.71) * x[1] - T(0.452) * y4 + T(0.58) * y3
+        c11 = T(1.75) * y2 * T(0.995) * x[1]
+        c12 = T(0.995) * y10 + T(1998)
+        y12 = c10 * x[1] + c11 / c12
+        y13 = c12 - T(1.75) * y2
+        y14 = T(3623) + T(64.4) * x[2] + T(58.4) * x[3] + T(146312) / (y9 + x[5])
+        c13 = T(0.995) * y10 + T(60.8) * x[2] + T(48) * x[4] -
+              T(0.1121) * y14 - T(5095)
+        y15 = y13 / c13
+        y16 = T(148000) - T(331000) * y15 + T(40) * y13 - T(61) * y15 * y13
+        c14 = T(2324) * y10 - T(28740000) * y2
+        y17 = T(14130000) - T(1328) * y10 - T(531) * y11 + c14 / c12
+        c15 = y13 / y15 - y13 / T(0.52)
+        c16 = T(1.104) - T(0.72) * y15
+        # c17 not needed for objective
+
+        return -T(5.843e-7) * y17 +
+               T(1.17e-4) * y14 +
+               T(2.358e-5) * y13 +
+               T(1.502e-6) * y16 +
+               T(0.0321) * y12 +
+               T(0.004324) * y5 +
+               T(1e-4) * c15 / c16 +
+               T(37.48) * y2 / c12 +
+               T(0.1365)
+    end
+
+      # Constraint function (48 inequalities, all of the form c(x) >= 0)
+      function c!(cx::AbstractVector{T}, x::AbstractVector{T})
+            # All intermediates (identical to those used in objective)
+            y1  = x[2] + x[3] + T(41.6)
+            c1  = T(0.024) * x[4] - T(4.62)
+            y2  = T(12.5) / c1 + T(12)
+            c2  = T(0.0003535) * x[1]^2 + T(0.5311) * x[1] + T(0.08705) * y2 * x[1]
+            c3  = T(0.052) * x[1] + T(78) + T(0.002377) * y2 * x[1]
+            y3  = c2 / c3
+            y4  = T(19) * y3
+            c4  = T(0.04782) * (x[1] - y3) + T(0.1956) * (x[1] - y3)^2 / x[2] + T(0.6376) * y4 + T(1.594) * y3
+            c5  = T(100) * x[2]
+            c6  = x[1] - y3 - y4
+            c7  = T(0.95) - c4 / c5
+            y5  = c6 * c7
+            y6  = x[1] - y5 - y4 - y3
+            c8  = (y5 + y4) * T(0.995)
+            y7  = c8 / y1
+            y8  = c8 / T(3798)
+            c9  = y7 - T(0.0663) * y7 / y8 - T(0.3153)
+            y9  = T(96.82) / c9 + T(0.321) * y1
+            y10 = T(1.29) * y5 + T(1.258) * y4 + T(2.29) * y3 + T(1.71) * y6
+            y11 = T(1.71) * x[1] - T(0.452) * y4 + T(0.58) * y3
+            c11 = T(1.75) * y2 * T(0.995) * x[1]
+            c12 = T(0.995) * y10 + T(1998)
+            y12 = c10 * x[1] + c11 / c12
+            y13 = c12 - T(1.75) * y2
+            y14 = T(3623) + T(64.4) * x[2] + T(58.4) * x[3] + T(146312) / (y9 + x[5])
+            c13 = T(0.995) * y10 + T(60.8) * x[2] + T(48) * x[4] - T(0.1121) * y14 - T(5095)
+            y15 = y13 / c13
+            y16 = T(148000) - T(331000) * y15 + T(40) * y13 - T(61) * y15 * y13
+            c14 = T(2324) * y10 - T(28740000) * y2
+            y17 = T(14130000) - T(1328) * y10 - T(531) * y11 + c14 / c12
+            c15 = y13 / y15 - y13 / T(0.52)
+            c16 = T(1.104) - T(0.72) * y15
+            c17 = y9 + x[5]
+            # Constraints
+            cx[1]  = T(1.5) * x[2] - x[3]
+            cx[2]  = y1 - T(213.1)
+            cx[3]  = T(405.23) - y1
+            cx[4]  = y2  - a[2]
+            cx[5]  = y3  - a[3]
+            cx[6]  = y4  - a[4]
+            cx[7]  = y5  - a[5]
+            cx[8]  = y6  - a[6]
+            cx[9]  = y7  - a[7]
+            cx[10] = y8  - a[8]
+            cx[11] = y9  - a[9]
+            cx[12] = y10 - a[10]
+            cx[13] = y11 - a[11]
+            cx[14] = y12 - a[12]
+            cx[15] = y13 - a[13]
+            cx[16] = y14 - a[14]
+            cx[17] = y15 - a[15]
+            cx[18] = y16 - a[16]
+            cx[19] = y17 - a[17]
+            cx[20] = b[2]  - y2
+            cx[21] = b[3]  - y3
+            cx[22] = b[4]  - y4
+            cx[23] = b[5]  - y5
+            cx[24] = b[6]  - y6
+            cx[25] = b[7]  - y7
+            cx[26] = b[8]  - y8
+            cx[27] = b[9]  - y9
+            cx[28] = b[10] - y10
+            cx[29] = b[11] - y11
+            cx[30] = b[12] - y12
+            cx[31] = b[13] - y13
+            cx[32] = b[14] - y14
+            cx[33] = b[15] - y15
+            cx[34] = b[16] - y16
+            cx[35] = b[17] - y17
+            cx[36] = y4 - (T(0.28) / T(0.72)) * y5
+            cx[37] = T(21) - T(3496) * y2 / c12
+            cx[38] = T(62212) / c17 - T(110.6) - y1
+            # 10 box constraints (5 lower, 5 upper)
+            cx[39] = x[1] - lvar[1]
+            cx[40] = x[2] - lvar[2]
+            cx[41] = x[3] - lvar[3]
+            cx[42] = x[4] - lvar[4]
+            cx[43] = x[5] - lvar[5]
+            cx[44] = uvar[1] - x[1]
+            cx[45] = uvar[2] - x[2]
+            cx[46] = uvar[3] - x[3]
+            cx[47] = uvar[4] - x[4]
+            cx[48] = uvar[5] - x[5]
+            return cx
+      end
+      # Constraint bounds: all inequalities of the form c(x) >= 0
+      m = 48
+      cl = zeros(T, m)
+      cu = fill(T(Inf), m)
+      return ADNLPModels.ADNLPModel!(
+            f, x0, lvar, uvar,
+            c!, cl, cu;
+            name = "hs85",
+            kwargs...
+      )
+end

src/ADNLPProblems/hs89.jl

@@ -0,0 +1,86 @@
+export hs89
+
+function hs89(; type::Type{T} = Float64, kwargs...) where T
+    # First 30 positive roots of tan(μ) = μ
+    mu = T[
+        0.8603335890193798,  3.425618459481728,  6.437298179171945,  9.529334405361963,
+       12.645287223856588, 15.771284874815820, 18.902409956860000, 22.036496727938500,
+       25.172446326646600, 28.309642854452000, 31.447714637546200, 34.586424215288900,
+       37.725612827776500, 40.865170330488000, 44.005017920830800, 47.145097736761000,
+       50.285366337773600, 53.425790477394600, 56.566344279821500, 59.707007305335400,
+       62.847763194454400, 65.988598698490300, 69.129502973895200, 72.270467060308900,
+       75.411483488848100, 78.552545984242900, 81.693649235601600, 84.834788718042200,
+       87.975960552493200, 91.117161394464700
+    ]
+
+    # Precomputed coefficients A_j = 2 sin(μ_j) / (μ_j + sin(μ_j) cos(μ_j))
+    A = [2 * sin(mu[j]) / (mu[j] + sin(mu[j]) * cos(mu[j])) for j in 1:30]
+
+    # Objective: φ(x) = ∑_{j=1}^{30} A_j ρ_j(x)
+    function f(x::AbstractVector{T})
+        s = zero(T)
+        r = x[1]^2 + x[2]^2 + x[3]^2
+        for j = 1:30
+            μ² = mu[j]^2
+            exp_r = exp(-μ² * r)
+            ρ = - (exp_r + 2*exp(-μ²*(x[2]^2 + x[3]^2)) + 2*exp(-μ²*x[3]^2) + 1) / μ²
+            s += A[j] * ρ
+        end
+        return s
+    end
+
+    # Equality constraint c(x) = 0
+    # Full expression with cross terms (double sum over i < j)
+    function c!(cx::AbstractVector{T}, x::AbstractVector{T})
+        r = x[1]^2 + x[2]^2 + x[3]^2
+        termA = zero(T)
+        termB = zero(T)
+        # Compute termA and termB directly, no heap allocation
+        for j = 1:30
+            μ = mu[j]
+            μ² = μ^2
+            exp_r  = exp(-μ² * r)
+            exp_r23 = exp(-μ² * (x[2]^2 + x[3]^2))
+            exp_r3  = exp(-μ² * x[3]^2)
+            ρj = - (exp_r + 2*exp_r23 + 2*exp_r3 + 1) / μ²
+            termA += A[j]^2 * ρj^2 * (sin(2*μ)/(2*μ) + one(T)) / 2
+        end
+        for i = 1:29
+            μi = mu[i]
+            μi² = μi^2
+            exp_ri  = exp(-μi² * r)
+            exp_ri23 = exp(-μi² * (x[2]^2 + x[3]^2))
+            exp_ri3  = exp(-μi² * x[3]^2)
+            ρi = - (exp_ri + 2*exp_ri23 + 2*exp_ri3 + 1) / μi²
+            for j = i+1:30
+                μj = mu[j]
+                μj² = μj^2
+                exp_rj  = exp(-μj² * r)
+                exp_rj23 = exp(-μj² * (x[2]^2 + x[3]^2))
+                exp_rj3  = exp(-μj² * x[3]^2)
+                ρj = - (exp_rj + 2*exp_rj23 + 2*exp_rj3 + 1) / μj²
+                denom_plus  = μi + μj
+                denom_minus = μi - μj
+                sin_plus  = iszero(denom_plus)  ? one(T) : sin(denom_plus)/denom_plus
+                sin_minus = iszero(denom_minus) ? one(T) : sin(denom_minus)/denom_minus
+                termB += A[i] * A[j] * ρi * ρj * (sin_plus + sin_minus)
+            end
+        end
+        cx[1] = termA + termB - T(2)/15
+        return cx
+    end
+
+    # Starting point (common in literature / CUTE)
+    x0 = T[0.5, -0.5, 0.5]
+
+    # One equality constraint c(x) = 0
+    lcon = ucon = T[0]
+
+    return ADNLPModels.ADNLPModel!(
+        f, x0,
+        c!, lcon, ucon;
+        name = "hs89",
+        lin = Int[],   # no linear constraints
+        kwargs...
+    )
+end

src/Meta/hs85.jl

@@ -0,0 +1,23 @@
+hs85_meta = Dict(
+  :nvar => 5,
+  :variable_nvar => false,
+  :ncon => 48,
+  :variable_ncon => false,
+  :minimize => true,
+  :name => "hs85",
+  :has_equalities_only => false,
+  :has_inequalities_only => true,
+  :has_bounds => true,
+  :has_fixed_variables => false,
+  :objtype => :other,
+  :contype => :general,
+  :best_known_lower_bound => -Inf,
+  :best_known_upper_bound => -1.90513375,
+  :is_feasible => true,
+  :defined_everywhere => missing,
+  :origin => :unknown,
+)
+get_hs85_nvar(; n::Integer = 5, kwargs...) = 5
+get_hs85_ncon(; n::Integer = 48, kwargs...) = 48
+get_hs85_nlin(; n::Integer = 48, kwargs...) = 0
+get_hs85_nnln(; n::Integer = 48, kwargs...) = 48

src/Meta/hs89.jl

@@ -0,0 +1,23 @@
+hs89_meta = Dict(
+  :nvar => 3,
+  :variable_nvar => false,
+  :ncon => 1,
+  :variable_ncon => false,
+  :minimize => true,
+  :name => "hs89",
+  :has_equalities_only => false,
+  :has_inequalities_only => true,
+  :has_bounds => false,
+  :has_fixed_variables => false,
+  :objtype => :other,
+  :contype => :general,
+  :best_known_lower_bound => -Inf,
+  :best_known_upper_bound => 1.36265681,
+  :is_feasible => true,
+  :defined_everywhere => missing,
+  :origin => :unknown,
+)
+get_hs89_nvar(; n::Integer = 3, kwargs...) = 3
+get_hs89_ncon(; n::Integer = 1, kwargs...) = 1
+get_hs89_nlin(; n::Integer = 1, kwargs...) = 0
+get_hs89_nnln(; n::Integer = 1, kwargs...) = 1