From f428189ff5f0ca5e3bb4a2f3cbb4abdfa8b8a55f Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Tue, 17 Mar 2026 06:33:26 +0530 Subject: [PATCH 1/9] Add hs85 and hs89 problems --- Project.toml | 2 + src/ADNLPProblems/hs85.jl | 164 ++++++++++++++++++++++++++++++++++++++ src/ADNLPProblems/hs89.jl | 86 ++++++++++++++++++++ src/Meta/hs85.jl | 23 ++++++ src/Meta/hs89.jl | 23 ++++++ src/PureJuMP/hs85.jl | 109 +++++++++++++++++++++++++ src/PureJuMP/hs89.jl | 95 ++++++++++++++++++++++ 7 files changed, 502 insertions(+) create mode 100644 src/ADNLPProblems/hs85.jl create mode 100644 src/ADNLPProblems/hs89.jl create mode 100644 src/Meta/hs85.jl create mode 100644 src/Meta/hs89.jl create mode 100644 src/PureJuMP/hs85.jl create mode 100644 src/PureJuMP/hs89.jl diff --git a/Project.toml b/Project.toml index f69a86a5..c2a178e8 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ 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" @@ -17,6 +18,7 @@ JLD2 = "0.5, 0.6" JuMP = "^1.15" Requires = "1" SpecialFunctions = "2" +StaticArrays = "1.9.18" julia = "1.6" [extras] diff --git a/src/ADNLPProblems/hs85.jl b/src/ADNLPProblems/hs85.jl new file mode 100644 index 00000000..01644ead --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/ADNLPProblems/hs89.jl b/src/ADNLPProblems/hs89.jl new file mode 100644 index 00000000..d6ec668b --- /dev/null +++ b/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 \ No newline at end of file diff --git a/src/Meta/hs85.jl b/src/Meta/hs85.jl new file mode 100644 index 00000000..eba8e319 --- /dev/null +++ b/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 diff --git a/src/Meta/hs89.jl b/src/Meta/hs89.jl new file mode 100644 index 00000000..13ac47fc --- /dev/null +++ b/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 diff --git a/src/PureJuMP/hs85.jl b/src/PureJuMP/hs85.jl new file mode 100644 index 00000000..a07f2df6 --- /dev/null +++ b/src/PureJuMP/hs85.jl @@ -0,0 +1,109 @@ +## Hock and Schittkowski problem number 85 +# +# Source: +# Problem 85 in +# W. Hock and K. Schittkowski, +# Test examples for nonlinear programming codes, +# Lectures Notes in Economics and Mathematical Systems 187, +# Springer Verlag, Heidelberg, 1981. +# +# classification QGR-P1-(1,...,6) +# +# Implementation: AI/JSO, 03/2026 + +export hs85 + +"HS85 model" +function hs85_jump() + m = Model(Ipopt.Optimizer) + set_attribute(m, "tol", 1e-10) + set_attribute(m, "print_level", 5) + + # Decision variables + @variable(m, 704.4148 ≤ x1 ≤ 906.3855) + @variable(m, 68.6 ≤ x2 ≤ 288.88) + @variable(m, 0.0 ≤ x3 ≤ 134.75) + @variable(m, 193.0 ≤ x4 ≤ 287.0966) + @variable(m, 25.0 ≤ x5 ≤ 84.1988) + + # Intermediates (as @expression or @variable + @NLconstraint) + @expression(m, y1, x2 + x3 + 41.6) + @expression(m, c1, 0.024 * x4 - 4.62) + @expression(m, y2, 12.5 / c1 + 12) + @expression(m, c2, 0.0003535 * x1^2 + 0.5311 * x1 + 0.08705 * y2 * x1) + @expression(m, c3, 0.052 * x1 + 78 + 0.002377 * y2 * x1) + @expression(m, y3, c2 / c3) + @expression(m, y4, 19 * y3) + @expression(m, c4, 0.04782 * (x1 - y3) + 0.1956 * (x1 - y3)^2 / x2 + + 0.6376 * y4 + 1.594 * y3) + @expression(m, c5, 100 * x2) + @expression(m, c6, x1 - y3 - y4) + @expression(m, c7, 0.95 - c4 / c5) + @expression(m, y5, c6 * c7) + @expression(m, y6, x1 - y5 - y4 - y3) + @expression(m, c8, (y5 + y4) * 0.995) + @expression(m, y7, c8 / y1) + @expression(m, y8, c8 / 3798) + @expression(m, c9, y7 - 0.0663 * y7 / y8 - 0.3153) + @expression(m, y9, 96.82 / c9 + 0.321 * y1) + @expression(m, y10, 1.29 * y5 + 1.258 * y4 + 2.29 * y3 + 1.71 * y6) + @expression(m, y11, 1.71 * x1 - 0.452 * y4 + 0.58 * y3) + @expression(m, c11, 1.75 * y2 * 0.995 * x1) + @expression(m, c12, 0.995 * y10 + 1998) + @expression(m, y12, (12.3/752.3) * x1 + c11 / c12) + @expression(m, y13, c12 - 1.75 * y2) + @expression(m, y14, 3623 + 64.4 * x2 + 58.4 * x3 + 146312 / (y9 + x5)) + @expression(m, c13, 0.995 * y10 + 60.8 * x2 + 48 * x4 - + 0.1121 * y14 - 5095) + @expression(m, y15, y13 / c13) + @expression(m, y16, 148000 - 331000 * y15 + 40 * y13 - 61 * y15 * y13) + @expression(m, c14, 2324 * y10 - 28740000 * y2) + @expression(m, y17, 14130000 - 1328 * y10 - 531 * y11 + c14 / c12) + @expression(m, c15, y13 / y15 - y13 / 0.52) + @expression(m, c16, 1.104 - 0.72 * y15) + @expression(m, c17, y9 + x5) + + # Bounds on y_i + a = [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 = [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] + + for i in 2:17 + yi = Symbol("y$i") + @constraint(m, getfield(Main, yi) >= a[i]) + @constraint(m, getfield(Main, yi) <= b[i]) + end + + # Other inequalities + @constraint(m, 1.5 * x2 - x3 >= 0) + @constraint(m, y1 >= 213.1) + @constraint(m, y1 <= 405.23) + @constraint(m, y4 >= (0.28 / 0.72) * y5) + @constraint(m, 3496 * y2 / c12 <= 21) + @constraint(m, 62212 / c17 - 110.6 >= y1) + + # Objective + @NLobjective(m, Min, + -5.843e-7 * y17 + + 1.17e-4 * y14 + + 2.358e-5 * y13 + + 1.502e-6 * y16 + + 0.0321 * y12 + + 0.004324 * y5 + + 1e-4 * (c15 / c16) + + 37.48 * y2 / c12 + + 0.1365 + ) + + # Good starting point helps a lot + set_start_value(x1, 900.0) + set_start_value(x2, 80.0) + set_start_value(x3, 115.0) + set_start_value(x4, 267.0) + set_start_value(x5, 27.0) + + return m +end \ No newline at end of file diff --git a/src/PureJuMP/hs89.jl b/src/PureJuMP/hs89.jl new file mode 100644 index 00000000..a36876fa --- /dev/null +++ b/src/PureJuMP/hs89.jl @@ -0,0 +1,95 @@ +# Hock and Schittkowski problem number 89 +# +# Source: +# Problem 89 in +# W. Hock and K. Schittkowski, +# Test examples for nonlinear programming codes, +# Lectures Notes in Economics and Mathematical Systems 187, +# Springer Verlag, Heidelberg, 1981. +# +# classification QGR-P1-(1,...,6) +# +# Implementation: AI/JSO, 03/2026 + +export hs89 + +"HS89 model" +function hs89_jump( + optimizer = Ipopt.Optimizer, + optimizer_attributes = Dict( + "tol" => 1e-10, + "print_level" => 5, + "max_iter" => 10000, + "acceptable_tol" => 1e-8 + ) +) + model = Model(optimizer) + + # Apply solver-specific options + for (k, v) in optimizer_attributes + set_optimizer_attribute(model, k, v) + end + + # Variables (no simple bounds in standard HS89, but loose ones help numerics) + @variable(model, x1) + @variable(model, x2) + @variable(model, x3) + + # First 30 positive roots of tan(μ) = μ + mu = [ + 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 + ] + + # Coefficients Aⱼ + A = [2 * sin(mu[j]) / (mu[j] + sin(mu[j]) * cos(mu[j])) for j in 1:30] + + # Precompute nothing — JuMP + AD will handle it + + # Helper expressions for ρⱼ (makes model more readable) + @expression(model, ρ[j=1:30], + let μ² = mu[j]^2, + r = x1^2 + x2^2 + x3^2, + r23 = x2^2 + x3^2, + r3 = x3^2; + -( exp(-μ² * r) + 2*exp(-μ² * r23) + 2*exp(-μ² * r3) + 1 ) / μ² + end + ) + + # Objective: ∑ Aⱼ ρⱼ + @NLobjective(model, Min, + sum(A[j] * ρ[j] for j in 1:30) + ) + + # Constraint: termA + termB = 2/15 + @NLconstraint(model, eq, + sum( + A[j]^2 * ρ[j]^2 * (sin(2*mu[j])/(2*mu[j]) + 1)/2 + for j in 1:30 + ) + + + sum( + sum( + A[i] * A[j] * ρ[i] * ρ[j] * ( + sin(mu[i] + mu[j]) / (mu[i] + mu[j]) + + sin(mu[i] - mu[j]) / (mu[i] - mu[j]) + ) + for j in i+1:30 + ) + for i in 1:29 + ) + == 2/15 + ) + + # Good starting point (from CUTE / literature) + set_start_value(x1, 0.5) + set_start_value(x2, -0.5) + set_start_value(x3, 0.5) + return model +end \ No newline at end of file From 763e55bb0ebaf9eb7253307315b02bfb11a939f0 Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Tue, 17 Mar 2026 06:46:38 +0530 Subject: [PATCH 2/9] doc files --- docs/RULES_hs85_hs89.md | 22 ++++++++++++++ docs/check_rules.jl | 49 ++++++++++++++++++++++++++++++ docs/extract_hs85.json | 66 +++++++++++++++++++++++++++++++++++++++++ docs/extract_hs89.json | 21 +++++++++++++ 4 files changed, 158 insertions(+) create mode 100644 docs/RULES_hs85_hs89.md create mode 100644 docs/check_rules.jl create mode 100644 docs/extract_hs85.json create mode 100644 docs/extract_hs89.json diff --git a/docs/RULES_hs85_hs89.md b/docs/RULES_hs85_hs89.md new file mode 100644 index 00000000..c5da5011 --- /dev/null +++ b/docs/RULES_hs85_hs89.md @@ -0,0 +1,22 @@ +# OptimizationProblems.jl Rules — HS85 & HS89 Example + +| Rule | HS85 | HS89 | How to Check | +|------|------|------|--------------| +| File Structure | src/ADNLPProblems/hs85.jl, src/PureJuMP/hs85.jl, src/Meta/hs85.md | src/ADNLPProblems/hs89.jl, src/PureJuMP/hs89.jl, src/Meta/hs89.md | Check for all three files per problem | +| Header | HS-style header: source, classification, implementation date | Same | Parse header comments for required fields | +| Mathematical Expressions | All intermediates, constraints, objective explicit, match paper | Same | Compare expressions to paper and extracted JSON | +| Variable Bounds | Explicit bounds and x0 in ADNLP and PureJuMP | Same | Check for bounds and x0 in both files | +| Metadata | References, classification, origin, scalable, etc. in Meta file | Same | Parse Meta file for required metadata | +| Naming | hs85.jl, hs85_jump.jl, hs85.md | hs89.jl, hs89_jump.jl, hs89.md | Check file/function names | +| Allocation | Constraint function minimizes allocations (relaxed for complex) | Same | Run allocation test or allow allocation | +| Ipopt Solve | Problem solves with Ipopt (PureJuMP) | Same | Run Ipopt solver and check for solution | +| Reviewer Markdown | Summary, PDF screenshot, extraction uncertainties, test results | Same | Generate reviewer markdown file | +| Duplication | No duplicate (by name, structure, metadata) | Same | Check for similar names/metadata | +| Traceability | Origin clear and referenced in header/Meta | Same | Check Meta/header | +| Scalability | Marked if scalable (not for hs85/hs89) | Same | Check Meta/implementation | +| Multiple Problems | One PDF = one problem (for hs85/hs89) | Same | Allow multiple if needed | +| Uncertainty | Warn if extraction unclear (not for hs85/hs89) | Same | Parse extraction JSON for uncertainties | + +--- + +This table is tailored to HS85 and HS89. Use as a checklist for similar problems and PRs. diff --git a/docs/check_rules.jl b/docs/check_rules.jl new file mode 100644 index 00000000..70716b04 --- /dev/null +++ b/docs/check_rules.jl @@ -0,0 +1,49 @@ +# Script to check OptimizationProblems.jl rules for HS85/HS89 +using JSON + +function check_problem(problem::String) + result = Dict() + # File structure + result["ADNLP"] = isfile("src/ADNLPProblems/$problem.jl") + result["PureJuMP"] = isfile("src/PureJuMP/$problem.jl") + result["Meta"] = isfile("src/Meta/$problem.md") + # Header + header_ok = false + if result["ADNLP"] + header = readlines("src/ADNLPProblems/$problem.jl")[1:10] + header_ok = any(occursin("Hock and Schittkowski", h) for h in header) + end + result["Header"] = header_ok + # Bounds/x0 + bounds_ok = false + if result["ADNLP"] + file = String(read("src/ADNLPProblems/$problem.jl")) + bounds_ok = occursin("lvar", file) && occursin("uvar", file) && occursin("x0", file) + end + result["Bounds"] = bounds_ok + # Metadata + meta_ok = false + if result["Meta"] + meta = String(read("src/Meta/$problem.md")) + meta_ok = occursin("Source", meta) && occursin("classification", meta) + end + result["Metadata"] = meta_ok + # Allocation (relaxed) + result["Allocation"] = false + # Ipopt solve (not checked here) + result["IpoptSolve"] = "manual" + # Reviewer markdown + result["ReviewerMarkdown"] = isfile("docs/review_$problem.md") + return result +end + +function main() + problems = ["hs85", "hs89"] + results = Dict() + for p in problems + results[p] = check_problem(p) + end + println(JSON.json(results)) +end + +main() diff --git a/docs/extract_hs85.json b/docs/extract_hs85.json new file mode 100644 index 00000000..e45013a2 --- /dev/null +++ b/docs/extract_hs85.json @@ -0,0 +1,66 @@ +{ + "problem": "hs85", + "source": "Hock and Schittkowski, Problem 85", + "variables": ["x1", "x2", "x3", "x4", "x5"], + "bounds": { + "lvar": [704.4148, 68.6, 0.0, 193.0, 25.0], + "uvar": [906.3855, 288.88, 134.75, 287.0966, 84.1988], + "x0": [900.0, 80.0, 115.0, 267.0, 27.0] + }, + "objective": "-5.843e-7 * y17 + 1.17e-4 * y14 + 2.358e-5 * y13 + 1.502e-6 * y16 + 0.0321 * y12 + 0.004324 * y5 + 1e-4 * c15 / c16 + 37.48 * y2 / c12 + 0.1365", + "constraints": [ + "1.5 * x2 - x3 >= 0", + "y1 - 213.1 >= 0", + "405.23 - y1 >= 0", + "y2 - a2 >= 0", + "y3 - a3 >= 0", + "y4 - a4 >= 0", + "y5 - a5 >= 0", + "y6 - a6 >= 0", + "y7 - a7 >= 0", + "y8 - a8 >= 0", + "y9 - a9 >= 0", + "y10 - a10 >= 0", + "y11 - a11 >= 0", + "y12 - a12 >= 0", + "y13 - a13 >= 0", + "y14 - a14 >= 0", + "y15 - a15 >= 0", + "y16 - a16 >= 0", + "y17 - a17 >= 0", + "b2 - y2 >= 0", + "b3 - y3 >= 0", + "b4 - y4 >= 0", + "b5 - y5 >= 0", + "b6 - y6 >= 0", + "b7 - y7 >= 0", + "b8 - y8 >= 0", + "b9 - y9 >= 0", + "b10 - y10 >= 0", + "b11 - y11 >= 0", + "b12 - y12 >= 0", + "b13 - y13 >= 0", + "b14 - y14 >= 0", + "b15 - y15 >= 0", + "b16 - y16 >= 0", + "b17 - y17 >= 0", + "y4 - (0.28 / 0.72) * y5 >= 0", + "21 - 3496 * y2 / c12 >= 0", + "62212 / c17 - 110.6 - y1 >= 0", + "x1 - lvar1 >= 0", + "x2 - lvar2 >= 0", + "x3 - lvar3 >= 0", + "x4 - lvar4 >= 0", + "x5 - lvar5 >= 0", + "uvar1 - x1 >= 0", + "uvar2 - x2 >= 0", + "uvar3 - x3 >= 0", + "uvar4 - x4 >= 0", + "uvar5 - x5 >= 0" + ], + "metadata": { + "classification": "QGR-P1-(1,...,6)", + "implementation": "AI/JSO, 03/2026" + }, + "uncertainties": [] +} diff --git a/docs/extract_hs89.json b/docs/extract_hs89.json new file mode 100644 index 00000000..172f3476 --- /dev/null +++ b/docs/extract_hs89.json @@ -0,0 +1,21 @@ +{ + "problem": "hs89", + "source": "Hock and Schittkowski, Problem 89", + "variables": ["x1", "x2", "x3"], + "bounds": { + "lvar": "none", + "uvar": "none", + "x0": "none" + }, + "objective": "sum_{j=1}^{30} A_j * rho_j(x)\nwhere rho_j(x) = - (exp(-mu_j^2 * r) + 2*exp(-mu_j^2*(x2^2 + x3^2)) + 2*exp(-mu_j^2*x3^2) + 1) / mu_j^2\nr = x1^2 + x2^2 + x3^2\nA_j = 2*sin(mu_j)/(mu_j + sin(mu_j)*cos(mu_j))\nmu_j: first 30 positive roots of tan(mu) = mu", + "constraints": [ + "c(x) = termA + termB - 2/15 = 0", + "termA = sum_{j=1}^{30} A_j^2 * rho_j(x)^2 * (sin(2*mu_j)/(2*mu_j) + 1)/2", + "termB = sum_{i=1}^{29} sum_{j=i+1}^{30} A_i * A_j * rho_i(x) * rho_j(x) * (sin(mu_i+mu_j)/(mu_i+mu_j) + sin(mu_i-mu_j)/(mu_i-mu_j))" + ], + "metadata": { + "classification": "QGR-P1-(1,...,6)", + "implementation": "AI/JSO, 03/2026" + }, + "uncertainties": [] +} From 726ca3298227acc4254022a271d2d9057b858005 Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Tue, 31 Mar 2026 03:01:01 +0530 Subject: [PATCH 3/9] changes --- docs/check_rules.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/docs/check_rules.jl b/docs/check_rules.jl index 70716b04..5c56927e 100644 --- a/docs/check_rules.jl +++ b/docs/check_rules.jl @@ -3,36 +3,36 @@ using JSON function check_problem(problem::String) result = Dict() - # File structure + # Check that the required files exist for the problem: ADNLPProblems, PureJuMP, and Meta. result["ADNLP"] = isfile("src/ADNLPProblems/$problem.jl") result["PureJuMP"] = isfile("src/PureJuMP/$problem.jl") result["Meta"] = isfile("src/Meta/$problem.md") - # Header + # Check that the header of the ADNLPProblems file for this problem contains a reference to "Hock and Schittkowski" in the first 10 lines, to ensure proper attribution/documentation. header_ok = false if result["ADNLP"] header = readlines("src/ADNLPProblems/$problem.jl")[1:10] header_ok = any(occursin("Hock and Schittkowski", h) for h in header) end result["Header"] = header_ok - # Bounds/x0 - bounds_ok = false + # Check that the ADNLPProblems file for problem defines lower bounds (lvar), upper bounds (uvar), and an initial point (x0).bounds_ok = false if result["ADNLP"] file = String(read("src/ADNLPProblems/$problem.jl")) bounds_ok = occursin("lvar", file) && occursin("uvar", file) && occursin("x0", file) end result["Bounds"] = bounds_ok - # Metadata + # Check that the Meta file for this problem contains both a "Source" and a "classification" entry. meta_ok = false if result["Meta"] meta = String(read("src/Meta/$problem.md")) meta_ok = occursin("Source", meta) && occursin("classification", meta) end result["Metadata"] = meta_ok - # Allocation (relaxed) + # Allocation check is currently not implemented (set to false as a placeholder). result["Allocation"] = false - # Ipopt solve (not checked here) + # Ipopt solve check is not performed by this script. + # Marked as "manual" for reviewer attention. result["IpoptSolve"] = "manual" - # Reviewer markdown + # Check that a reviewer markdown file exists for this problem in docs/review_$problem.md. result["ReviewerMarkdown"] = isfile("docs/review_$problem.md") return result end From f10898a8c4b9a538b07b05906f7e4161f286da4f Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Tue, 31 Mar 2026 03:22:48 +0530 Subject: [PATCH 4/9] run JuliaFormatter --- src/ADNLPProblems/hs85.jl | 344 +++++++++++++++++++++----------------- src/ADNLPProblems/hs89.jl | 177 +++++++++++--------- src/PureJuMP/hs85.jl | 199 +++++++++++++--------- src/PureJuMP/hs89.jl | 151 +++++++++-------- 4 files changed, 486 insertions(+), 385 deletions(-) diff --git a/src/ADNLPProblems/hs85.jl b/src/ADNLPProblems/hs85.jl index 01644ead..dc64284a 100644 --- a/src/ADNLPProblems/hs85.jl +++ b/src/ADNLPProblems/hs85.jl @@ -1,164 +1,192 @@ 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) +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] + # 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 + # 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 + 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 \ No newline at end of file + # 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 diff --git a/src/ADNLPProblems/hs89.jl b/src/ADNLPProblems/hs89.jl index d6ec668b..3c360ee0 100644 --- a/src/ADNLPProblems/hs89.jl +++ b/src/ADNLPProblems/hs89.jl @@ -1,86 +1,111 @@ 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 - ] +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] + # 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 = 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 + # 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 + # 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] + # Starting point (common in literature / CUTE) + x0 = T[0.5, -0.5, 0.5] - # One equality constraint c(x) = 0 - lcon = ucon = T[0] + # 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 \ No newline at end of file + return ADNLPModels.ADNLPModel!( + f, + x0, + c!, + lcon, + ucon; + name = "hs89", + lin = Int[], # no linear constraints + kwargs..., + ) +end diff --git a/src/PureJuMP/hs85.jl b/src/PureJuMP/hs85.jl index a07f2df6..dd55570d 100644 --- a/src/PureJuMP/hs85.jl +++ b/src/PureJuMP/hs85.jl @@ -15,95 +15,128 @@ export hs85 "HS85 model" function hs85_jump() - m = Model(Ipopt.Optimizer) - set_attribute(m, "tol", 1e-10) - set_attribute(m, "print_level", 5) + m = Model(Ipopt.Optimizer) + set_attribute(m, "tol", 1e-10) + set_attribute(m, "print_level", 5) - # Decision variables - @variable(m, 704.4148 ≤ x1 ≤ 906.3855) - @variable(m, 68.6 ≤ x2 ≤ 288.88) - @variable(m, 0.0 ≤ x3 ≤ 134.75) - @variable(m, 193.0 ≤ x4 ≤ 287.0966) - @variable(m, 25.0 ≤ x5 ≤ 84.1988) + # Decision variables + @variable(m, 704.4148 ≤ x1 ≤ 906.3855) + @variable(m, 68.6 ≤ x2 ≤ 288.88) + @variable(m, 0.0 ≤ x3 ≤ 134.75) + @variable(m, 193.0 ≤ x4 ≤ 287.0966) + @variable(m, 25.0 ≤ x5 ≤ 84.1988) - # Intermediates (as @expression or @variable + @NLconstraint) - @expression(m, y1, x2 + x3 + 41.6) - @expression(m, c1, 0.024 * x4 - 4.62) - @expression(m, y2, 12.5 / c1 + 12) - @expression(m, c2, 0.0003535 * x1^2 + 0.5311 * x1 + 0.08705 * y2 * x1) - @expression(m, c3, 0.052 * x1 + 78 + 0.002377 * y2 * x1) - @expression(m, y3, c2 / c3) - @expression(m, y4, 19 * y3) - @expression(m, c4, 0.04782 * (x1 - y3) + 0.1956 * (x1 - y3)^2 / x2 + - 0.6376 * y4 + 1.594 * y3) - @expression(m, c5, 100 * x2) - @expression(m, c6, x1 - y3 - y4) - @expression(m, c7, 0.95 - c4 / c5) - @expression(m, y5, c6 * c7) - @expression(m, y6, x1 - y5 - y4 - y3) - @expression(m, c8, (y5 + y4) * 0.995) - @expression(m, y7, c8 / y1) - @expression(m, y8, c8 / 3798) - @expression(m, c9, y7 - 0.0663 * y7 / y8 - 0.3153) - @expression(m, y9, 96.82 / c9 + 0.321 * y1) - @expression(m, y10, 1.29 * y5 + 1.258 * y4 + 2.29 * y3 + 1.71 * y6) - @expression(m, y11, 1.71 * x1 - 0.452 * y4 + 0.58 * y3) - @expression(m, c11, 1.75 * y2 * 0.995 * x1) - @expression(m, c12, 0.995 * y10 + 1998) - @expression(m, y12, (12.3/752.3) * x1 + c11 / c12) - @expression(m, y13, c12 - 1.75 * y2) - @expression(m, y14, 3623 + 64.4 * x2 + 58.4 * x3 + 146312 / (y9 + x5)) - @expression(m, c13, 0.995 * y10 + 60.8 * x2 + 48 * x4 - - 0.1121 * y14 - 5095) - @expression(m, y15, y13 / c13) - @expression(m, y16, 148000 - 331000 * y15 + 40 * y13 - 61 * y15 * y13) - @expression(m, c14, 2324 * y10 - 28740000 * y2) - @expression(m, y17, 14130000 - 1328 * y10 - 531 * y11 + c14 / c12) - @expression(m, c15, y13 / y15 - y13 / 0.52) - @expression(m, c16, 1.104 - 0.72 * y15) - @expression(m, c17, y9 + x5) + # Intermediates (as @expression or @variable + @NLconstraint) + @expression(m, y1, x2 + x3 + 41.6) + @expression(m, c1, 0.024 * x4 - 4.62) + @expression(m, y2, 12.5 / c1 + 12) + @expression(m, c2, 0.0003535 * x1^2 + 0.5311 * x1 + 0.08705 * y2 * x1) + @expression(m, c3, 0.052 * x1 + 78 + 0.002377 * y2 * x1) + @expression(m, y3, c2 / c3) + @expression(m, y4, 19 * y3) + @expression(m, c4, 0.04782 * (x1 - y3) + 0.1956 * (x1 - y3)^2 / x2 + 0.6376 * y4 + 1.594 * y3) + @expression(m, c5, 100 * x2) + @expression(m, c6, x1 - y3 - y4) + @expression(m, c7, 0.95 - c4 / c5) + @expression(m, y5, c6 * c7) + @expression(m, y6, x1 - y5 - y4 - y3) + @expression(m, c8, (y5 + y4) * 0.995) + @expression(m, y7, c8 / y1) + @expression(m, y8, c8 / 3798) + @expression(m, c9, y7 - 0.0663 * y7 / y8 - 0.3153) + @expression(m, y9, 96.82 / c9 + 0.321 * y1) + @expression(m, y10, 1.29 * y5 + 1.258 * y4 + 2.29 * y3 + 1.71 * y6) + @expression(m, y11, 1.71 * x1 - 0.452 * y4 + 0.58 * y3) + @expression(m, c11, 1.75 * y2 * 0.995 * x1) + @expression(m, c12, 0.995 * y10 + 1998) + @expression(m, y12, (12.3/752.3) * x1 + c11 / c12) + @expression(m, y13, c12 - 1.75 * y2) + @expression(m, y14, 3623 + 64.4 * x2 + 58.4 * x3 + 146312 / (y9 + x5)) + @expression(m, c13, 0.995 * y10 + 60.8 * x2 + 48 * x4 - 0.1121 * y14 - 5095) + @expression(m, y15, y13 / c13) + @expression(m, y16, 148000 - 331000 * y15 + 40 * y13 - 61 * y15 * y13) + @expression(m, c14, 2324 * y10 - 28740000 * y2) + @expression(m, y17, 14130000 - 1328 * y10 - 531 * y11 + c14 / c12) + @expression(m, c15, y13 / y15 - y13 / 0.52) + @expression(m, c16, 1.104 - 0.72 * y15) + @expression(m, c17, y9 + x5) - # Bounds on y_i - a = [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] + # Bounds on y_i + a = [ + 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 = [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] + b = [ + 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, + ] - for i in 2:17 - yi = Symbol("y$i") - @constraint(m, getfield(Main, yi) >= a[i]) - @constraint(m, getfield(Main, yi) <= b[i]) - end + for i = 2:17 + yi = Symbol("y$i") + @constraint(m, getfield(Main, yi) >= a[i]) + @constraint(m, getfield(Main, yi) <= b[i]) + end - # Other inequalities - @constraint(m, 1.5 * x2 - x3 >= 0) - @constraint(m, y1 >= 213.1) - @constraint(m, y1 <= 405.23) - @constraint(m, y4 >= (0.28 / 0.72) * y5) - @constraint(m, 3496 * y2 / c12 <= 21) - @constraint(m, 62212 / c17 - 110.6 >= y1) + # Other inequalities + @constraint(m, 1.5 * x2 - x3 >= 0) + @constraint(m, y1 >= 213.1) + @constraint(m, y1 <= 405.23) + @constraint(m, y4 >= (0.28 / 0.72) * y5) + @constraint(m, 3496 * y2 / c12 <= 21) + @constraint(m, 62212 / c17 - 110.6 >= y1) - # Objective - @NLobjective(m, Min, - -5.843e-7 * y17 + - 1.17e-4 * y14 + - 2.358e-5 * y13 + - 1.502e-6 * y16 + - 0.0321 * y12 + - 0.004324 * y5 + - 1e-4 * (c15 / c16) + - 37.48 * y2 / c12 + - 0.1365 - ) + # Objective + @NLobjective( + m, + Min, + -5.843e-7 * y17 + + 1.17e-4 * y14 + + 2.358e-5 * y13 + + 1.502e-6 * y16 + + 0.0321 * y12 + + 0.004324 * y5 + + 1e-4 * (c15 / c16) + + 37.48 * y2 / c12 + + 0.1365 + ) - # Good starting point helps a lot - set_start_value(x1, 900.0) - set_start_value(x2, 80.0) - set_start_value(x3, 115.0) - set_start_value(x4, 267.0) - set_start_value(x5, 27.0) + # Good starting point helps a lot + set_start_value(x1, 900.0) + set_start_value(x2, 80.0) + set_start_value(x3, 115.0) + set_start_value(x4, 267.0) + set_start_value(x5, 27.0) - return m -end \ No newline at end of file + return m +end diff --git a/src/PureJuMP/hs89.jl b/src/PureJuMP/hs89.jl index a36876fa..15897292 100644 --- a/src/PureJuMP/hs89.jl +++ b/src/PureJuMP/hs89.jl @@ -15,81 +15,96 @@ export hs89 "HS89 model" function hs89_jump( - optimizer = Ipopt.Optimizer, - optimizer_attributes = Dict( - "tol" => 1e-10, - "print_level" => 5, - "max_iter" => 10000, - "acceptable_tol" => 1e-8 - ) + optimizer = Ipopt.Optimizer, + optimizer_attributes = Dict( + "tol" => 1e-10, + "print_level" => 5, + "max_iter" => 10000, + "acceptable_tol" => 1e-8, + ), ) - model = Model(optimizer) + model = Model(optimizer) - # Apply solver-specific options - for (k, v) in optimizer_attributes - set_optimizer_attribute(model, k, v) - end + # Apply solver-specific options + for (k, v) in optimizer_attributes + set_optimizer_attribute(model, k, v) + end - # Variables (no simple bounds in standard HS89, but loose ones help numerics) - @variable(model, x1) - @variable(model, x2) - @variable(model, x3) + # Variables (no simple bounds in standard HS89, but loose ones help numerics) + @variable(model, x1) + @variable(model, x2) + @variable(model, x3) - # First 30 positive roots of tan(μ) = μ - mu = [ - 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 - ] + # First 30 positive roots of tan(μ) = μ + mu = [ + 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, + ] - # Coefficients Aⱼ - A = [2 * sin(mu[j]) / (mu[j] + sin(mu[j]) * cos(mu[j])) for j in 1:30] + # Coefficients Aⱼ + A = [2 * sin(mu[j]) / (mu[j] + sin(mu[j]) * cos(mu[j])) for j = 1:30] - # Precompute nothing — JuMP + AD will handle it + # Precompute nothing — JuMP + AD will handle it - # Helper expressions for ρⱼ (makes model more readable) - @expression(model, ρ[j=1:30], - let μ² = mu[j]^2, - r = x1^2 + x2^2 + x3^2, - r23 = x2^2 + x3^2, - r3 = x3^2; - -( exp(-μ² * r) + 2*exp(-μ² * r23) + 2*exp(-μ² * r3) + 1 ) / μ² - end - ) + # Helper expressions for ρⱼ (makes model more readable) + @expression( + model, + ρ[j = 1:30], + let μ² = mu[j]^2, r = x1^2 + x2^2 + x3^2, r23 = x2^2 + x3^2, r3 = x3^2; + -(exp(-μ² * r) + 2*exp(-μ² * r23) + 2*exp(-μ² * r3) + 1) / μ² + end + ) - # Objective: ∑ Aⱼ ρⱼ - @NLobjective(model, Min, - sum(A[j] * ρ[j] for j in 1:30) - ) + # Objective: ∑ Aⱼ ρⱼ + @NLobjective(model, Min, sum(A[j] * ρ[j] for j = 1:30)) - # Constraint: termA + termB = 2/15 - @NLconstraint(model, eq, - sum( - A[j]^2 * ρ[j]^2 * (sin(2*mu[j])/(2*mu[j]) + 1)/2 - for j in 1:30 - ) - + - sum( - sum( - A[i] * A[j] * ρ[i] * ρ[j] * ( - sin(mu[i] + mu[j]) / (mu[i] + mu[j]) + - sin(mu[i] - mu[j]) / (mu[i] - mu[j]) - ) - for j in i+1:30 - ) - for i in 1:29 - ) - == 2/15 - ) + # Constraint: termA + termB = 2/15 + @NLconstraint( + model, + eq, + sum(A[j]^2 * ρ[j]^2 * (sin(2*mu[j])/(2*mu[j]) + 1)/2 for j = 1:30) + sum( + sum( + A[i] * + A[j] * + ρ[i] * + ρ[j] * + (sin(mu[i] + mu[j]) / (mu[i] + mu[j]) + sin(mu[i] - mu[j]) / (mu[i] - mu[j])) for + j = (i + 1):30 + ) for i = 1:29 + ) == 2/15 + ) - # Good starting point (from CUTE / literature) - set_start_value(x1, 0.5) - set_start_value(x2, -0.5) - set_start_value(x3, 0.5) - return model -end \ No newline at end of file + # Good starting point (from CUTE / literature) + set_start_value(x1, 0.5) + set_start_value(x2, -0.5) + set_start_value(x3, 0.5) + return model +end From f30253121215a1dd043955bbcfc30604274728d8 Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Mon, 6 Apr 2026 04:49:48 +0530 Subject: [PATCH 5/9] add notebook --- docs/check_rules.jl | 54 +++++++++++++++++++++++++++++----------- docs/hs85_workflow.ipynb | 0 2 files changed, 39 insertions(+), 15 deletions(-) create mode 100644 docs/hs85_workflow.ipynb diff --git a/docs/check_rules.jl b/docs/check_rules.jl index 5c56927e..59465e04 100644 --- a/docs/check_rules.jl +++ b/docs/check_rules.jl @@ -3,28 +3,52 @@ using JSON function check_problem(problem::String) result = Dict() + adnlp_path = "src/ADNLPProblems/$problem.jl" + jump_path = "src/PureJuMP/$problem.jl" + meta_path = "src/Meta/$problem.jl" # Check that the required files exist for the problem: ADNLPProblems, PureJuMP, and Meta. - result["ADNLP"] = isfile("src/ADNLPProblems/$problem.jl") - result["PureJuMP"] = isfile("src/PureJuMP/$problem.jl") - result["Meta"] = isfile("src/Meta/$problem.md") + result["ADNLP"] = isfile(adnlp_path) + result["PureJuMP"] = isfile(jump_path) + result["Meta"] = isfile(meta_path) # Check that the header of the ADNLPProblems file for this problem contains a reference to "Hock and Schittkowski" in the first 10 lines, to ensure proper attribution/documentation. header_ok = false - if result["ADNLP"] - header = readlines("src/ADNLPProblems/$problem.jl")[1:10] - header_ok = any(occursin("Hock and Schittkowski", h) for h in header) + if result["ADNLP"] || result["PureJuMP"] + adnlp_header_ok = false + jump_header_ok = false + if result["ADNLP"] + adnlp_lines = readlines(adnlp_path) + k = min(length(adnlp_lines), 20) + adnlp_header_ok = any(occursin("Hock and Schittkowski", adnlp_lines[i]) for i in 1:k) + end + if result["PureJuMP"] + jump_lines = readlines(jump_path) + k = min(length(jump_lines), 20) + jump_header_ok = any(occursin("Hock and Schittkowski", jump_lines[i]) for i in 1:k) + end + header_ok = adnlp_header_ok || jump_header_ok end result["Header"] = header_ok - # Check that the ADNLPProblems file for problem defines lower bounds (lvar), upper bounds (uvar), and an initial point (x0).bounds_ok = false + # Check that the ADNLPProblems file for problem defines lower bounds (lvar), upper bounds (uvar), and an initial point (x0). + bounds_ok = false if result["ADNLP"] - file = String(read("src/ADNLPProblems/$problem.jl")) - bounds_ok = occursin("lvar", file) && occursin("uvar", file) && occursin("x0", file) + adnlp_file = String(read(adnlp_path)) + has_x0 = occursin("x0", adnlp_file) + has_bounds_in_code = occursin("lvar", adnlp_file) && occursin("uvar", adnlp_file) + has_bounds_required = true + if result["Meta"] + meta_file = String(read(meta_path)) + if occursin(":has_bounds => false", meta_file) + has_bounds_required = false + end + end + bounds_ok = has_x0 && (has_bounds_in_code || !has_bounds_required) end result["Bounds"] = bounds_ok # Check that the Meta file for this problem contains both a "Source" and a "classification" entry. meta_ok = false if result["Meta"] - meta = String(read("src/Meta/$problem.md")) - meta_ok = occursin("Source", meta) && occursin("classification", meta) + meta = String(read(meta_path)) + meta_ok = occursin(":name => \"$problem\"", meta) && occursin(":best_known_upper_bound", meta) end result["Metadata"] = meta_ok # Allocation check is currently not implemented (set to false as a placeholder). @@ -32,18 +56,18 @@ function check_problem(problem::String) # Ipopt solve check is not performed by this script. # Marked as "manual" for reviewer attention. result["IpoptSolve"] = "manual" - # Check that a reviewer markdown file exists for this problem in docs/review_$problem.md. - result["ReviewerMarkdown"] = isfile("docs/review_$problem.md") + # Check whether reviewer markdown exists (optional, but useful for PR traceability). + result["ReviewerMarkdown"] = isfile("docs/review_$problem.md") || isfile("docs/review/$problem.md") return result end function main() - problems = ["hs85", "hs89"] + problems = isempty(ARGS) ? ["hs85", "hs89"] : ARGS results = Dict() for p in problems results[p] = check_problem(p) end - println(JSON.json(results)) + println(JSON.json(results, 2)) end main() diff --git a/docs/hs85_workflow.ipynb b/docs/hs85_workflow.ipynb new file mode 100644 index 00000000..e69de29b From 355f0686e236319dd1840df65a7b2fe3e661bbe8 Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Mon, 6 Apr 2026 04:50:19 +0530 Subject: [PATCH 6/9] update --- docs/hs85_workflow.ipynb | 247 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 247 insertions(+) diff --git a/docs/hs85_workflow.ipynb b/docs/hs85_workflow.ipynb index e69de29b..0f4096ff 100644 --- a/docs/hs85_workflow.ipynb +++ b/docs/hs85_workflow.ipynb @@ -0,0 +1,247 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "ff939ce8", + "metadata": {}, + "source": [ + "# HS85 End-to-End Workflow (PR #412 style)\n", + "\n", + "This notebook demonstrates one complete flow for adding/verifying a problem in `OptimizationProblems.jl` using **HS85** as a focused test case.\n", + "\n", + "Flow covered:\n", + "1. Load extraction artifact (`docs/extract_hs85.json`).\n", + "2. Run repository rule checks (files, header, bounds, metadata).\n", + "3. Compare extracted bounds with implementation bounds.\n", + "4. Instantiate the ADNLP problem and inspect dimensions.\n", + "5. Optionally solve the PureJuMP model with Ipopt (if available).\n", + "\n", + "The goal is to make the process reviewable and reproducible in one place." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dd6abb58", + "metadata": { + "vscode": { + "languageId": "julia" + } + }, + "outputs": [], + "source": [ + "using Pkg\n", + "\n", + "repo_root = if isfile(joinpath(pwd(), \"Project.toml\"))\n", + " pwd()\n", + "elseif isfile(joinpath(pwd(), \"..\", \"Project.toml\"))\n", + " normpath(joinpath(pwd(), \"..\"))\n", + "else\n", + " error(\"Could not locate Project.toml from current working directory: $(pwd())\")\n", + "end\n", + "\n", + "Pkg.activate(repo_root)\n", + "\n", + "using JSON\n", + "using OptimizationProblems\n", + "using ADNLPModels\n", + "using NLPModels\n", + "\n", + "println(\"Activated project at: \", repo_root)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cfaf1484", + "metadata": { + "vscode": { + "languageId": "julia" + } + }, + "outputs": [], + "source": [ + "extract_path = joinpath(repo_root, \"docs\", \"extract_hs85.json\")\n", + "extract = JSON.parsefile(extract_path)\n", + "\n", + "println(\"Loaded extraction file: \", extract_path)\n", + "println(\"Problem: \", extract[\"problem\"])\n", + "println(\"Variables: \", length(extract[\"variables\"]))\n", + "println(\"Constraints listed in JSON: \", length(extract[\"constraints\"]))\n", + "println(\"Uncertainties count: \", length(extract[\"uncertainties\"]))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7654a919", + "metadata": { + "vscode": { + "languageId": "julia" + } + }, + "outputs": [], + "source": [ + "function first_lines_have_pattern(path::String, pattern::String; n::Int = 20)\n", + " if !isfile(path)\n", + " return false\n", + " end\n", + " lines = readlines(path)\n", + " k = min(length(lines), n)\n", + " return any(occursin(pattern, lines[i]) for i in 1:k)\n", + "end\n", + "\n", + "function check_problem_files_and_rules(repo_root::String, problem::String)\n", + " adnlp_path = joinpath(repo_root, \"src\", \"ADNLPProblems\", \"$(problem).jl\")\n", + " jump_path = joinpath(repo_root, \"src\", \"PureJuMP\", \"$(problem).jl\")\n", + " meta_path = joinpath(repo_root, \"src\", \"Meta\", \"$(problem).jl\")\n", + "\n", + " file_ok = Dict(\n", + " \"ADNLP\" => isfile(adnlp_path),\n", + " \"PureJuMP\" => isfile(jump_path),\n", + " \"Meta\" => isfile(meta_path),\n", + " )\n", + "\n", + " header_ok = first_lines_have_pattern(adnlp_path, \"Hock and Schittkowski\")\n", + "\n", + " adnlp_text = file_ok[\"ADNLP\"] ? read(adnlp_path, String) : \"\"\n", + " bounds_ok = occursin(\"lvar\", adnlp_text) && occursin(\"uvar\", adnlp_text) && occursin(\"x0\", adnlp_text)\n", + "\n", + " meta_text = file_ok[\"Meta\"] ? read(meta_path, String) : \"\"\n", + " metadata_ok = occursin(\":name => \\\"$(problem)\\\"\", meta_text) &&\n", + " occursin(\":best_known_upper_bound\", meta_text)\n", + "\n", + " return Dict(\n", + " \"problem\" => problem,\n", + " \"files\" => file_ok,\n", + " \"header\" => header_ok,\n", + " \"bounds\" => bounds_ok,\n", + " \"metadata\" => metadata_ok,\n", + " \"all_core_checks_pass\" => all(values(file_ok)) && header_ok && bounds_ok && metadata_ok,\n", + " )\n", + "end\n", + "\n", + "checks = check_problem_files_and_rules(repo_root, \"hs85\")\n", + "println(JSON.json(checks, 2))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "78215d1d", + "metadata": { + "vscode": { + "languageId": "julia" + } + }, + "outputs": [], + "source": [ + "function parse_numeric_vector_literal(text::String, name::String)\n", + " m = match(Regex(name * raw\"\\s*=\\s*T\\[(.*?)\\]\", \"s\"), text)\n", + " m === nothing && error(\"Could not locate vector literal for $(name)\")\n", + " body = m.captures[1]\n", + " tokens = split(replace(body, '\\n' => ' '), ',')\n", + " vals = Float64[]\n", + " for t in tokens\n", + " s = strip(t)\n", + " isempty(s) && continue\n", + " push!(vals, parse(Float64, s))\n", + " end\n", + " return vals\n", + "end\n", + "\n", + "adnlp_text = read(joinpath(repo_root, \"src\", \"ADNLPProblems\", \"hs85.jl\"), String)\n", + "code_lvar = parse_numeric_vector_literal(adnlp_text, \"lvar\")\n", + "code_uvar = parse_numeric_vector_literal(adnlp_text, \"uvar\")\n", + "code_x0 = parse_numeric_vector_literal(adnlp_text, \"x0\")\n", + "\n", + "json_lvar = Float64.(extract[\"bounds\"][\"lvar\"])\n", + "json_uvar = Float64.(extract[\"bounds\"][\"uvar\"])\n", + "json_x0 = Float64.(extract[\"bounds\"][\"x0\"])\n", + "\n", + "println(\"Bounds and x0 consistency with extraction JSON:\")\n", + "println(\"lvar match: \", code_lvar == json_lvar)\n", + "println(\"uvar match: \", code_uvar == json_uvar)\n", + "println(\"x0 match: \", code_x0 == json_x0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c8224124", + "metadata": { + "vscode": { + "languageId": "julia" + } + }, + "outputs": [], + "source": [ + "try\n", + " nlp = OptimizationProblems.ADNLPProblems.hs85()\n", + " x0 = nlp.meta.x0\n", + " f0 = NLPModels.obj(nlp, x0)\n", + " c0 = NLPModels.cons(nlp, x0)\n", + "\n", + " println(\"ADNLP model summary:\")\n", + " println(\"name: \", nlp.meta.name)\n", + " println(\"nvar: \", nlp.meta.nvar)\n", + " println(\"ncon: \", nlp.meta.ncon)\n", + " println(\"objective at x0: \", f0)\n", + " println(\"min constraint value at x0 (for c(x) >= 0 style constraints): \", minimum(c0))\n", + "catch err\n", + " println(\"Skipped ADNLP instantiation step. Reason: \", err)\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5863d263", + "metadata": { + "vscode": { + "languageId": "julia" + } + }, + "outputs": [], + "source": [ + "try\n", + " using JuMP, Ipopt\n", + " model = OptimizationProblems.PureJuMP.hs85_jump()\n", + " optimize!(model)\n", + " println(\"Ipopt termination status: \", termination_status(model))\n", + " println(\"Objective value: \", objective_value(model))\n", + "catch err\n", + " println(\"Skipped solve step. Reason: \", err)\n", + " println(\"This is expected if Ipopt/JuMP are not available in the current notebook environment.\")\n", + "end" + ] + }, + { + "cell_type": "markdown", + "id": "46e5809a", + "metadata": {}, + "source": [ + "## Manual vs AI Intervention Notes\n", + "\n", + "What this notebook automates well:\n", + "- Structural checks (file existence, headers, basic metadata patterns).\n", + "- Consistency checks between extraction JSON and implementation bounds/x0.\n", + "- Basic model instantiation and optional solve sanity-check.\n", + "\n", + "What still needs manual review:\n", + "- Ambiguous PDF interpretation and formula disambiguation.\n", + "- Naming choices and duplicate detection decisions.\n", + "- Semantic fidelity of objective/constraints against original publication.\n", + "\n", + "This balance matches PR #412 workflow: AI accelerates drafting/checking, while humans validate mathematical intent and project policy decisions." + ] + } + ], + "metadata": { + "language_info": { + "name": "python" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} From dbd831160165d73bd87aab6348123433cd8df543 Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Mon, 6 Apr 2026 05:24:53 +0530 Subject: [PATCH 7/9] Apply suggestions from code review Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- docs/RULES_hs85_hs89.md | 4 ++-- docs/extract_hs89.json | 2 +- src/Meta/hs89.jl | 4 ++-- src/PureJuMP/hs85.jl | 6 ++---- src/PureJuMP/hs89.jl | 5 ++++- 5 files changed, 11 insertions(+), 10 deletions(-) diff --git a/docs/RULES_hs85_hs89.md b/docs/RULES_hs85_hs89.md index c5da5011..6b8555a9 100644 --- a/docs/RULES_hs85_hs89.md +++ b/docs/RULES_hs85_hs89.md @@ -2,12 +2,12 @@ | Rule | HS85 | HS89 | How to Check | |------|------|------|--------------| -| File Structure | src/ADNLPProblems/hs85.jl, src/PureJuMP/hs85.jl, src/Meta/hs85.md | src/ADNLPProblems/hs89.jl, src/PureJuMP/hs89.jl, src/Meta/hs89.md | Check for all three files per problem | +| File Structure | src/ADNLPProblems/hs85.jl, src/PureJuMP/hs85.jl, src/Meta/hs85.jl | src/ADNLPProblems/hs89.jl, src/PureJuMP/hs89.jl, src/Meta/hs89.jl | Check for all three files per problem | | Header | HS-style header: source, classification, implementation date | Same | Parse header comments for required fields | | Mathematical Expressions | All intermediates, constraints, objective explicit, match paper | Same | Compare expressions to paper and extracted JSON | | Variable Bounds | Explicit bounds and x0 in ADNLP and PureJuMP | Same | Check for bounds and x0 in both files | | Metadata | References, classification, origin, scalable, etc. in Meta file | Same | Parse Meta file for required metadata | -| Naming | hs85.jl, hs85_jump.jl, hs85.md | hs89.jl, hs89_jump.jl, hs89.md | Check file/function names | +| Naming | hs85.jl in ADNLPProblems, hs85.jl in PureJuMP, hs85.jl in Meta | hs89.jl in ADNLPProblems, hs89.jl in PureJuMP, hs89.jl in Meta | Check file/function names against repository conventions | | Allocation | Constraint function minimizes allocations (relaxed for complex) | Same | Run allocation test or allow allocation | | Ipopt Solve | Problem solves with Ipopt (PureJuMP) | Same | Run Ipopt solver and check for solution | | Reviewer Markdown | Summary, PDF screenshot, extraction uncertainties, test results | Same | Generate reviewer markdown file | diff --git a/docs/extract_hs89.json b/docs/extract_hs89.json index 172f3476..12668aa5 100644 --- a/docs/extract_hs89.json +++ b/docs/extract_hs89.json @@ -5,7 +5,7 @@ "bounds": { "lvar": "none", "uvar": "none", - "x0": "none" + "x0": [0.0, 0.0, 0.0] }, "objective": "sum_{j=1}^{30} A_j * rho_j(x)\nwhere rho_j(x) = - (exp(-mu_j^2 * r) + 2*exp(-mu_j^2*(x2^2 + x3^2)) + 2*exp(-mu_j^2*x3^2) + 1) / mu_j^2\nr = x1^2 + x2^2 + x3^2\nA_j = 2*sin(mu_j)/(mu_j + sin(mu_j)*cos(mu_j))\nmu_j: first 30 positive roots of tan(mu) = mu", "constraints": [ diff --git a/src/Meta/hs89.jl b/src/Meta/hs89.jl index 13ac47fc..5163b9bb 100644 --- a/src/Meta/hs89.jl +++ b/src/Meta/hs89.jl @@ -5,8 +5,8 @@ hs89_meta = Dict( :variable_ncon => false, :minimize => true, :name => "hs89", - :has_equalities_only => false, - :has_inequalities_only => true, + :has_equalities_only => true, + :has_inequalities_only => false, :has_bounds => false, :has_fixed_variables => false, :objtype => :other, diff --git a/src/PureJuMP/hs85.jl b/src/PureJuMP/hs85.jl index dd55570d..2b364f26 100644 --- a/src/PureJuMP/hs85.jl +++ b/src/PureJuMP/hs85.jl @@ -14,10 +14,8 @@ export hs85 "HS85 model" -function hs85_jump() - m = Model(Ipopt.Optimizer) - set_attribute(m, "tol", 1e-10) - set_attribute(m, "print_level", 5) +function hs85() + m = Model() # Decision variables @variable(m, 704.4148 ≤ x1 ≤ 906.3855) diff --git a/src/PureJuMP/hs89.jl b/src/PureJuMP/hs89.jl index 15897292..c939819d 100644 --- a/src/PureJuMP/hs89.jl +++ b/src/PureJuMP/hs89.jl @@ -21,8 +21,11 @@ function hs89_jump( "print_level" => 5, "max_iter" => 10000, "acceptable_tol" => 1e-8, - ), + ); + n::Int = 3, + kwargs..., ) + _ = n, kwargs model = Model(optimizer) # Apply solver-specific options From a4c4bdc58525eedc428306e1eaf4c95c2929a2b1 Mon Sep 17 00:00:00 2001 From: Arnav Kapoor Date: Mon, 6 Apr 2026 05:32:14 +0530 Subject: [PATCH 8/9] Update src/PureJuMP/hs89.jl Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- src/PureJuMP/hs89.jl | 19 ++++++++----------- 1 file changed, 8 insertions(+), 11 deletions(-) diff --git a/src/PureJuMP/hs89.jl b/src/PureJuMP/hs89.jl index c939819d..74988492 100644 --- a/src/PureJuMP/hs89.jl +++ b/src/PureJuMP/hs89.jl @@ -18,19 +18,16 @@ function hs89_jump( optimizer = Ipopt.Optimizer, optimizer_attributes = Dict( "tol" => 1e-10, - "print_level" => 5, - "max_iter" => 10000, - "acceptable_tol" => 1e-8, - ); - n::Int = 3, - kwargs..., + optimizer = nothing, + optimizer_attributes = nothing, ) - _ = n, kwargs - model = Model(optimizer) + model = optimizer === nothing ? Model() : Model(optimizer) - # Apply solver-specific options - for (k, v) in optimizer_attributes - set_optimizer_attribute(model, k, v) + # Apply solver-specific options only when an optimizer and attributes are provided + if optimizer_attributes !== nothing + for (k, v) in optimizer_attributes + set_optimizer_attribute(model, k, v) + end end # Variables (no simple bounds in standard HS89, but loose ones help numerics) From d30b774f52339e881b14725f317772be59848a29 Mon Sep 17 00:00:00 2001 From: arnavk23 Date: Mon, 6 Apr 2026 05:49:23 +0530 Subject: [PATCH 9/9] adding copilot suggestion --- docs/hs85_workflow.ipynb | 47 +++++++++++++-------------------------- src/ADNLPProblems/hs85.jl | 15 ++----------- src/Meta/hs85.jl | 8 +++---- src/PureJuMP/hs85.jl | 8 +++---- src/PureJuMP/hs89.jl | 14 ++++++------ 5 files changed, 32 insertions(+), 60 deletions(-) diff --git a/docs/hs85_workflow.ipynb b/docs/hs85_workflow.ipynb index 0f4096ff..2c4671c4 100644 --- a/docs/hs85_workflow.ipynb +++ b/docs/hs85_workflow.ipynb @@ -23,11 +23,7 @@ "cell_type": "code", "execution_count": null, "id": "dd6abb58", - "metadata": { - "vscode": { - "languageId": "julia" - } - }, + "metadata": {}, "outputs": [], "source": [ "using Pkg\n", @@ -54,11 +50,7 @@ "cell_type": "code", "execution_count": null, "id": "cfaf1484", - "metadata": { - "vscode": { - "languageId": "julia" - } - }, + "metadata": {}, "outputs": [], "source": [ "extract_path = joinpath(repo_root, \"docs\", \"extract_hs85.json\")\n", @@ -75,11 +67,7 @@ "cell_type": "code", "execution_count": null, "id": "7654a919", - "metadata": { - "vscode": { - "languageId": "julia" - } - }, + "metadata": {}, "outputs": [], "source": [ "function first_lines_have_pattern(path::String, pattern::String; n::Int = 20)\n", @@ -129,11 +117,7 @@ "cell_type": "code", "execution_count": null, "id": "78215d1d", - "metadata": { - "vscode": { - "languageId": "julia" - } - }, + "metadata": {}, "outputs": [], "source": [ "function parse_numeric_vector_literal(text::String, name::String)\n", @@ -169,11 +153,7 @@ "cell_type": "code", "execution_count": null, "id": "c8224124", - "metadata": { - "vscode": { - "languageId": "julia" - } - }, + "metadata": {}, "outputs": [], "source": [ "try\n", @@ -197,16 +177,12 @@ "cell_type": "code", "execution_count": null, "id": "5863d263", - "metadata": { - "vscode": { - "languageId": "julia" - } - }, + "metadata": {}, "outputs": [], "source": [ "try\n", " using JuMP, Ipopt\n", - " model = OptimizationProblems.PureJuMP.hs85_jump()\n", + " model = OptimizationProblems.PureJuMP.hs85()\n", " optimize!(model)\n", " println(\"Ipopt termination status: \", termination_status(model))\n", " println(\"Objective value: \", objective_value(model))\n", @@ -238,10 +214,17 @@ } ], "metadata": { + "kernelspec": { + "display_name": "julia 1.11", + "language": "julia", + "name": "julia" + }, "language_info": { - "name": "python" + "name": "julia", + "version": "3.11.9" } }, "nbformat": 4, "nbformat_minor": 5 } + diff --git a/src/ADNLPProblems/hs85.jl b/src/ADNLPProblems/hs85.jl index dc64284a..13562068 100644 --- a/src/ADNLPProblems/hs85.jl +++ b/src/ADNLPProblems/hs85.jl @@ -95,7 +95,7 @@ function hs85(; type::Type{T} = Float64, kwargs...) where {T} T(0.1365) end - # Constraint function (48 inequalities, all of the form c(x) >= 0) + # Constraint function (38 nonlinear 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) @@ -171,21 +171,10 @@ function hs85(; type::Type{T} = Float64, kwargs...) where {T} 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 + m = 38 cl = zeros(T, m) cu = fill(T(Inf), m) return ADNLPModels.ADNLPModel!(f, x0, lvar, uvar, c!, cl, cu; name = "hs85", kwargs...) diff --git a/src/Meta/hs85.jl b/src/Meta/hs85.jl index eba8e319..0cbd12a7 100644 --- a/src/Meta/hs85.jl +++ b/src/Meta/hs85.jl @@ -1,7 +1,7 @@ hs85_meta = Dict( :nvar => 5, :variable_nvar => false, - :ncon => 48, + :ncon => 38, :variable_ncon => false, :minimize => true, :name => "hs85", @@ -18,6 +18,6 @@ hs85_meta = Dict( :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 +get_hs85_ncon(; n::Integer = 38, kwargs...) = 38 +get_hs85_nlin(; n::Integer = 38, kwargs...) = 0 +get_hs85_nnln(; n::Integer = 38, kwargs...) = 38 diff --git a/src/PureJuMP/hs85.jl b/src/PureJuMP/hs85.jl index 2b364f26..f9db0bfc 100644 --- a/src/PureJuMP/hs85.jl +++ b/src/PureJuMP/hs85.jl @@ -14,7 +14,7 @@ export hs85 "HS85 model" -function hs85() +function hs85(args...; kwargs...) m = Model() # Decision variables @@ -100,10 +100,10 @@ function hs85() 12146108, ] + y = Any[y1, y2, y3, y4, y5, y6, y7, y8, y9, y10, y11, y12, y13, y14, y15, y16, y17] for i = 2:17 - yi = Symbol("y$i") - @constraint(m, getfield(Main, yi) >= a[i]) - @constraint(m, getfield(Main, yi) <= b[i]) + @constraint(m, y[i] >= a[i]) + @constraint(m, y[i] <= b[i]) end # Other inequalities diff --git a/src/PureJuMP/hs89.jl b/src/PureJuMP/hs89.jl index 74988492..b03773f3 100644 --- a/src/PureJuMP/hs89.jl +++ b/src/PureJuMP/hs89.jl @@ -14,15 +14,18 @@ export hs89 "HS89 model" -function hs89_jump( - optimizer = Ipopt.Optimizer, - optimizer_attributes = Dict( - "tol" => 1e-10, +function hs89( + args...; optimizer = nothing, optimizer_attributes = nothing, + kwargs..., ) model = optimizer === nothing ? Model() : Model(optimizer) + if optimizer !== nothing && optimizer_attributes === nothing + optimizer_attributes = Dict("tol" => 1e-10) + end + # Apply solver-specific options only when an optimizer and attributes are provided if optimizer_attributes !== nothing for (k, v) in optimizer_attributes @@ -72,9 +75,6 @@ function hs89_jump( # Coefficients Aⱼ A = [2 * sin(mu[j]) / (mu[j] + sin(mu[j]) * cos(mu[j])) for j = 1:30] - # Precompute nothing — JuMP + AD will handle it - - # Helper expressions for ρⱼ (makes model more readable) @expression( model, ρ[j = 1:30],