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/docs/RULES_hs85_hs89.md b/docs/RULES_hs85_hs89.md new file mode 100644 index 00000000..6b8555a9 --- /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.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 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 | +| 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..59465e04 --- /dev/null +++ b/docs/check_rules.jl @@ -0,0 +1,73 @@ +# Script to check OptimizationProblems.jl rules for HS85/HS89 +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(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"] || 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 + if result["ADNLP"] + 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(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). + result["Allocation"] = false + # Ipopt solve check is not performed by this script. + # Marked as "manual" for reviewer attention. + result["IpoptSolve"] = "manual" + # 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 = isempty(ARGS) ? ["hs85", "hs89"] : ARGS + results = Dict() + for p in problems + results[p] = check_problem(p) + end + println(JSON.json(results, 2)) +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..12668aa5 --- /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": [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": [ + "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": [] +} diff --git a/docs/hs85_workflow.ipynb b/docs/hs85_workflow.ipynb new file mode 100644 index 00000000..2c4671c4 --- /dev/null +++ b/docs/hs85_workflow.ipynb @@ -0,0 +1,230 @@ +{ + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "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": {}, + "outputs": [], + "source": [ + "try\n", + " using JuMP, Ipopt\n", + " model = OptimizationProblems.PureJuMP.hs85()\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": { + "kernelspec": { + "display_name": "julia 1.11", + "language": "julia", + "name": "julia" + }, + "language_info": { + "name": "julia", + "version": "3.11.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} + diff --git a/src/ADNLPProblems/hs85.jl b/src/ADNLPProblems/hs85.jl new file mode 100644 index 00000000..13562068 --- /dev/null +++ b/src/ADNLPProblems/hs85.jl @@ -0,0 +1,181 @@ +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 (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) + 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 + return cx + end + # Constraint bounds: all inequalities of the form c(x) >= 0 + m = 38 + 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 new file mode 100644 index 00000000..3c360ee0 --- /dev/null +++ b/src/ADNLPProblems/hs89.jl @@ -0,0 +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, + ] + + # 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 + 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 diff --git a/src/Meta/hs85.jl b/src/Meta/hs85.jl new file mode 100644 index 00000000..0cbd12a7 --- /dev/null +++ b/src/Meta/hs85.jl @@ -0,0 +1,23 @@ +hs85_meta = Dict( + :nvar => 5, + :variable_nvar => false, + :ncon => 38, + :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 = 38, kwargs...) = 38 +get_hs85_nlin(; n::Integer = 38, kwargs...) = 0 +get_hs85_nnln(; n::Integer = 38, kwargs...) = 38 diff --git a/src/Meta/hs89.jl b/src/Meta/hs89.jl new file mode 100644 index 00000000..5163b9bb --- /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 => true, + :has_inequalities_only => false, + :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..f9db0bfc --- /dev/null +++ b/src/PureJuMP/hs85.jl @@ -0,0 +1,140 @@ +## 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(args...; kwargs...) + m = Model() + + # 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, + ] + + y = Any[y1, y2, y3, y4, y5, y6, y7, y8, y9, y10, y11, y12, y13, y14, y15, y16, y17] + for i = 2:17 + @constraint(m, y[i] >= a[i]) + @constraint(m, y[i] <= 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 diff --git a/src/PureJuMP/hs89.jl b/src/PureJuMP/hs89.jl new file mode 100644 index 00000000..b03773f3 --- /dev/null +++ b/src/PureJuMP/hs89.jl @@ -0,0 +1,110 @@ +# 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( + 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 + set_optimizer_attribute(model, k, v) + end + 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 = 1:30] + + @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 = 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 = 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