diff --git a/docs/Project.toml b/docs/Project.toml index 3389bf931..180b42f57 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -2,6 +2,7 @@ ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" +JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" JSOSolvers = "10dff2fc-5484-5881-a0e0-c90441020f8a" NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6" NLPModelsJuMP = "792afdf1-32c1-5681-94e0-d7bf7a5df49e" diff --git a/docs/RULES_hs85_hs89.md b/docs/RULES_hs85_hs89.md new file mode 100644 index 000000000..2ed1147f9 --- /dev/null +++ b/docs/RULES_hs85_hs89.md @@ -0,0 +1,32 @@ +# OptimizationProblems.jl Rules — HS85 & HS89 Example + +This table doubles as the specification for `docs/check_rules.jl`: every row with a +`✓` in "Automated" is checked by that script (`julia --project=docs docs/check_rules.jl +...`); everything else needs a human reviewer. + +| Rule | Automated | 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 | ✓ | Source reference in the PureJuMP file header | Same | Parse header comments for a "Source" line | +| 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 (unless `:has_bounds => false`) | Same | Check for `lvar`/`uvar`/`x0` in the ADNLP file | +| Metadata | ✓ | `:name`, `:best_known_upper_bound`, and all of `get_hs85_{nvar,ncon,nlin,nnln,nequ,nineq}` in Meta file | Same | Parse Meta file for required fields and getters | +| Naming | ✓ (via File Structure) | 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 | +| JuMP Interface Consistency | ✓ | PureJuMP file uses only the modern interface (`@objective`/`@constraint`/`@expression`), never mixed with `@NLobjective`/`@NLconstraint`/`@NLexpression` | Same | Grep for `@NL*` macros in the PureJuMP file — JuMP refuses to build a model mixing both, so this crashes `MathOptNLPModel` | +| AD-Generic Signatures | ✓ | ADNLPProblems objective/constraint closures are not typed to a concrete `x::AbstractVector{T}` | Same | Grep the ADNLP file for `::AbstractVector{T}` argument annotations — these silently block ForwardDiff | +| ADNLP/PureJuMP Compatibility | ✓ | `nvar`, `x0`, `ncon`, objective, `cons()`, `lcon`/`ucon`, and `lin` agree between the two formulations at `x0` | Same | Instantiate both models and compare (this is `test/test-utils.jl`'s `test_compatibility`, run early) | +| Allocation | ✓ (informational) | `cons_nln!` allocation count reported | Same | Run `cons_nln!` twice and check `@allocated` is 0 | +| Ipopt Solve | manual | Problem solves with Ipopt (PureJuMP) | Same | Run Ipopt solver and check for solution | +| Reviewer Markdown | ✓ (presence only) | 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 | ✓ (via Header) | 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, but the automated rows and +`docs/check_rules.jl` are written to run on any problem name. +They were used as-is to catch a JuMP-model-construction crash, a ForwardDiff-breaking type annotation, and a full constraint-by-constraint ADNLP/PureJuMP mismatch in this +PR before merge. Use this table 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 000000000..cd96964ae --- /dev/null +++ b/docs/check_rules.jl @@ -0,0 +1,171 @@ +# Script to check OptimizationProblems.jl rules for HS85/HS89 +using JSON + +# All paths are resolved relative to the repository root (the parent of this file's +# directory), not the process's current working directory: this script is meant to +# be run both as `julia docs/check_rules.jl` from the repo root and `include()`d +# from docs/hs85_workflow.ipynb, whose cwd is docs/. +const REPO_ROOT = normpath(joinpath(@__DIR__, "..")) + +# Static (text-based) checks: do not require loading the package, so they run even +# if the package fails to precompile. +function check_problem_static(problem::String) + result = Dict{String, Any}() + adnlp_path = joinpath(REPO_ROOT, "src", "ADNLPProblems", "$problem.jl") + jump_path = joinpath(REPO_ROOT, "src", "PureJuMP", "$problem.jl") + meta_path = joinpath(REPO_ROOT, "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) + + adnlp_text = result["ADNLP"] ? read(adnlp_path, String) : "" + jump_text = result["PureJuMP"] ? read(jump_path, String) : "" + meta_text = result["Meta"] ? read(meta_path, String) : "" + + # Check that the header comment of the PureJuMP file (the file that consistently + # carries one across the repository) has a "Source" line in the first 20 lines, + # so the origin of the problem is traceable. + jump_lines = result["PureJuMP"] ? readlines(jump_path) : String[] + result["Header"] = + any(occursin("Source", jump_lines[i]) for i = 1:min(length(jump_lines), 20)) + + # Check that the ADNLPProblems file defines an initial point (x0), and lower/upper + # bounds (lvar/uvar) unless the Meta file explicitly says the problem has none. + has_bounds_required = !occursin(":has_bounds => false", meta_text) + result["Bounds"] = + result["ADNLP"] && + occursin("x0", adnlp_text) && + (!has_bounds_required || (occursin("lvar", adnlp_text) && occursin("uvar", adnlp_text))) + + # Check that the Meta file names the problem, records a best known upper bound, and + # defines every getter required by docs/src/contributing.md + # (get__nvar/ncon/nlin/nnln/nequ/nineq). + required_getters = ["nvar", "ncon", "nlin", "nnln", "nequ", "nineq"] + missing_getters = [g for g in required_getters if !occursin("get_$(problem)_$(g)(", meta_text)] + result["Metadata"] = + result["Meta"] && + occursin(":name => \"$problem\"", meta_text) && + occursin(":best_known_upper_bound", meta_text) && + isempty(missing_getters) + if !isempty(missing_getters) + result["MissingGetters"] = missing_getters + end + + # Check that the PureJuMP file does not mix the legacy nonlinear interface + # (@NLobjective/@NLconstraint/@NLexpression) with the modern one + # (@objective/@constraint/@expression): JuMP refuses to build a model that uses + # both, so MathOptNLPModel(model) throws instead of failing a specific test. + uses_legacy_macro = + occursin("@NLobjective", jump_text) || + occursin("@NLconstraint", jump_text) || + occursin("@NLexpression", jump_text) + result["JuMPInterfaceConsistent"] = result["PureJuMP"] && !uses_legacy_macro + + # Check that the objective/constraint closures in ADNLPProblems are not typed to a + # concrete vector element type (e.g. `x::AbstractVector{T}`), which would prevent + # ForwardDiff from calling them with Dual-typed vectors and break grad()/jacobian(). + result["ADNLPGenericSignatures"] = + result["ADNLP"] && !occursin(r"\([^)]*::AbstractVector\{T\}", adnlp_text) + + # Reviewer markdown (optional, but useful for PR traceability). + result["ReviewerMarkdown"] = + isfile(joinpath(REPO_ROOT, "docs", "review_$problem.md")) || + isfile(joinpath(REPO_ROOT, "docs", "review", "$problem.md")) + + return result +end + +# Whether OptimizationProblems and its test dependencies can be loaded, computed once. +const HAVE_DYNAMIC = try + @eval using OptimizationProblems, ADNLPModels, NLPModelsJuMP, NLPModels + true +catch err + @warn "Dynamic checks unavailable: could not load OptimizationProblems/ADNLPModels/NLPModelsJuMP/NLPModels ($err)" + false +end + +# Dynamic checks: instantiate the actual models. These catch what static text checks +# cannot, e.g. a JuMP model that throws on construction, or an ADNLPProblems model +# whose objective/constraints numerically disagree with its PureJuMP sibling. +function check_problem_dynamic(problem::String) + out = Dict{String, Any}() + if !HAVE_DYNAMIC + out["Allocation"] = "skipped (package unavailable)" + out["Compatibility"] = "skipped (package unavailable)" + return out + end + prob = Symbol(problem) + ad_mod = OptimizationProblems.ADNLPProblems + jump_mod = OptimizationProblems.PureJuMP + + if !isdefined(ad_mod, prob) + out["Allocation"] = "skipped (no ADNLPProblems.$problem)" + else + try + nlp = getfield(ad_mod, prob)(matrix_free = true) + if nlp.meta.nnln > 0 + x = nlp.meta.x0 + cx = similar(x, nlp.meta.nnln) + NLPModels.cons_nln!(nlp, x, cx) + nbytes = @allocated NLPModels.cons_nln!(nlp, x, cx) + out["Allocation"] = nbytes == 0 ? "pass (0 bytes)" : "warn ($(nbytes) bytes)" + else + out["Allocation"] = "n/a (no nonlinear constraints)" + end + catch err + out["Allocation"] = "error: $err" + end + end + + if !isdefined(ad_mod, prob) || !isdefined(jump_mod, prob) + out["Compatibility"] = "skipped (missing ADNLPProblems or PureJuMP implementation)" + else + try + nlp_ad = getfield(ad_mod, prob)(matrix_free = true) + model = getfield(jump_mod, prob)() + nlp_jump = NLPModelsJuMP.MathOptNLPModel(model; name = problem) + x0 = nlp_ad.meta.x0 + + failures = String[] + nlp_jump.meta.nvar == nlp_ad.meta.nvar || push!(failures, "nvar mismatch") + nlp_jump.meta.x0 == nlp_ad.meta.x0 || push!(failures, "x0 mismatch") + nlp_jump.meta.ncon == nlp_ad.meta.ncon || push!(failures, "ncon mismatch") + isapprox(NLPModels.obj(nlp_ad, x0), NLPModels.obj(nlp_jump, x0), rtol = 1e-6) || + push!(failures, "objective mismatch at x0") + if nlp_ad.meta.ncon > 0 && nlp_jump.meta.ncon == nlp_ad.meta.ncon + nlp_ad.meta.lcon ≈ nlp_jump.meta.lcon || push!(failures, "lcon mismatch") + nlp_ad.meta.ucon ≈ nlp_jump.meta.ucon || push!(failures, "ucon mismatch") + all(isapprox.(NLPModels.cons(nlp_ad, x0), NLPModels.cons(nlp_jump, x0), atol = 1e-6)) || + push!(failures, "cons() mismatch at x0") + nlp_ad.meta.lin == nlp_jump.meta.lin || push!(failures, "lin (linear constraint indices) mismatch") + end + out["Compatibility"] = isempty(failures) ? "pass" : "fail: " * join(failures, "; ") + catch err + out["Compatibility"] = "error: $err" + end + end + return out +end + +function check_problem(problem::String) + result = check_problem_static(problem) + merge!(result, check_problem_dynamic(problem)) + # Solving with Ipopt is not run automatically; left for manual reviewer attention. + result["IpoptSolve"] = "manual" + 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 + +if abspath(PROGRAM_FILE) == @__FILE__ + main() +end diff --git a/docs/extract_hs85.json b/docs/extract_hs85.json new file mode 100644 index 000000000..e45013a2b --- /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 000000000..53b5548ca --- /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.5, -0.5, 0.5] + }, + "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 000000000..458462304 --- /dev/null +++ b/docs/hs85_workflow.ipynb @@ -0,0 +1,462 @@ +{ + "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 the repository rule checks from `docs/check_rules.jl` (files, header, bounds,\n", + " metadata, JuMP interface consistency, AD-generic signatures).\n", + "3. Compare extracted bounds with implementation bounds.\n", + "4. Instantiate the ADNLP problem and inspect dimensions.\n", + "5. Check that ADNLPProblems and PureJuMP agree (objective, constraints, bounds,\n", + " linear/nonlinear split) and that `cons_nln!` is allocation-free.\n", + "6. Optionally solve the PureJuMP model with Ipopt (if available).\n", + "\n", + "The goal is to make the process reviewable and reproducible in one place. Steps 2\n", + "and 5 are not cosmetic: running them against the first draft of HS85/HS89 in this\n", + "PR caught a JuMP model that crashed on construction (mixed legacy/modern nonlinear\n", + "macros), a type annotation that silently broke ForwardDiff, and a full\n", + "constraint-by-constraint mismatch between the ADNLP and PureJuMP formulations —\n", + "all before a human reviewer had to notice them by hand." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "dd6abb58", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-03T22:29:28.728000Z", + "iopub.status.busy": "2026-09-03T22:29:28.728000Z", + "iopub.status.idle": "2026-09-03T22:29:49.085000Z", + "shell.execute_reply": "2026-09-03T22:29:49.085000Z" + } + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "\u001b[32m\u001b[1m Activating\u001b[22m\u001b[39m project at `C:\\Users\\kapoo\\Downloads\\OptimizationProblems.jl`\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Activated project at: C:\\Users\\kapoo\\Downloads\\OptimizationProblems.jl\\\n" + ] + } + ], + "source": [ + "using Pkg\n", + "\n", + "# Look for the repository root specifically (not just any Project.toml — docs/ has\n", + "# its own docs/Project.toml, so a plain `isfile(\"Project.toml\")` check matches there\n", + "# too when the notebook is opened with cwd == docs/).\n", + "is_repo_root(dir) = isfile(joinpath(dir, \"src\", \"OptimizationProblems.jl\"))\n", + "\n", + "repo_root = if is_repo_root(pwd())\n", + " pwd()\n", + "elseif is_repo_root(normpath(joinpath(pwd(), \"..\")))\n", + " normpath(joinpath(pwd(), \"..\"))\n", + "else\n", + " error(\"Could not locate the OptimizationProblems.jl repository root 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": 2, + "id": "cfaf1484", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-03T22:29:50.917000Z", + "iopub.status.busy": "2026-09-03T22:29:49.089000Z", + "iopub.status.idle": "2026-09-03T22:29:53.247000Z", + "shell.execute_reply": "2026-09-03T22:29:53.191000Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Loaded extraction file: C:\\Users\\kapoo\\Downloads\\OptimizationProblems.jl\\docs\\extract_hs85.json\n", + "Problem: hs85\n", + "Variables: 5\n", + "Constraints listed in JSON: 48\n", + "Uncertainties count: 0\n" + ] + } + ], + "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": 3, + "id": "7654a919", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-03T22:29:53.250000Z", + "iopub.status.busy": "2026-09-03T22:29:53.249000Z", + "iopub.status.idle": "2026-09-03T22:30:17.573000Z", + "shell.execute_reply": "2026-09-03T22:30:17.573000Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{\n", + " \"ADNLP\": true,\n", + " \"ADNLPGenericSignatures\": true,\n", + " \"Allocation\": \"pass (0 bytes)\",\n", + " \"Bounds\": true,\n", + " \"Compatibility\": \"pass\",\n", + " \"Header\": true,\n", + " \"IpoptSolve\": \"manual\",\n", + " \"JuMPInterfaceConsistent\": true,\n", + " \"Meta\": true,\n", + " \"Metadata\": true,\n", + " \"PureJuMP\": true,\n", + " \"ReviewerMarkdown\": false\n", + "}\n" + ] + } + ], + "source": [ + "include(joinpath(repo_root, \"docs\", \"check_rules.jl\"))\n", + "\n", + "checks = check_problem(\"hs85\")\n", + "println(JSON.json(checks, 2))" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "78215d1d", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-03T22:30:17.576000Z", + "iopub.status.busy": "2026-09-03T22:30:17.576000Z", + "iopub.status.idle": "2026-09-03T22:30:17.872000Z", + "shell.execute_reply": "2026-09-03T22:30:17.872000Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Bounds and x0 consistency with extraction JSON:\n", + "lvar match: true\n", + "uvar match: true\n", + "x0 match: true\n" + ] + } + ], + "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": 5, + "id": "c8224124", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-03T22:30:17.874000Z", + "iopub.status.busy": "2026-09-03T22:30:17.874000Z", + "iopub.status.idle": "2026-09-03T22:30:18.447000Z", + "shell.execute_reply": "2026-09-03T22:30:18.447000Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ADNLP model summary:\n", + "name: hs85\n", + "nvar: 5\n", + "ncon: 38 (3 linear, 35 nonlinear)\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "objective at x0: -0.9393968794311149\n", + "constraints satisfied at x0: 38 / 38\n" + ] + } + ], + "source": [ + "try\n", + " nlp = OptimizationProblems.ADNLPProblems.hs85(matrix_free = true)\n", + " x0 = nlp.meta.x0\n", + " f0 = NLPModels.obj(nlp, x0)\n", + " c0 = NLPModels.cons(nlp, x0)\n", + " # hs85 mixes linear and nonlinear constraints with different senses (>=, <=, and\n", + " # ranges), so \"feasible\" means each cons() entry falls within its own [lcon, ucon]\n", + " # bound, not a single uniform sign.\n", + " within_bounds = nlp.meta.lcon .- 1e-8 .<= c0 .<= nlp.meta.ucon .+ 1e-8\n", + "\n", + " println(\"ADNLP model summary:\")\n", + " println(\"name: \", nlp.meta.name)\n", + " println(\"nvar: \", nlp.meta.nvar)\n", + " println(\"ncon: \", nlp.meta.ncon, \" (\", nlp.meta.nlin, \" linear, \", nlp.meta.nnln, \" nonlinear)\")\n", + " println(\"objective at x0: \", f0)\n", + " println(\"constraints satisfied at x0: \", count(within_bounds), \" / \", length(within_bounds))\n", + "catch err\n", + " println(\"Skipped ADNLP instantiation step. Reason: \", err)\n", + "end" + ] + }, + { + "cell_type": "markdown", + "id": "0b76a11d", + "metadata": {}, + "source": [ + "## ADNLP / PureJuMP compatibility and allocation\n", + "\n", + "This is the check that matters most for correctness: `test/test-utils.jl` requires\n", + "the ADNLPProblems and PureJuMP formulations of a problem to agree exactly on\n", + "`nvar`, `x0`, `ncon`, `lcon`/`ucon`, `lin` (which constraint indices are linear),\n", + "and the numerical value of the objective and constraints at a given point — and it\n", + "requires `cons_nln!` to be allocation-free. `check_problem_dynamic` in\n", + "`docs/check_rules.jl` runs the same comparison so it can be caught here instead of\n", + "in CI. Running it against the first draft of `hs85`/`hs89` in this PR reported\n", + "`Compatibility => \"fail: ...\"` for every mismatched constraint, which is exactly\n", + "how the constraint reordering/sign bug fixed in this PR was tracked down." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "376b9979", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-03T22:30:18.450000Z", + "iopub.status.busy": "2026-09-03T22:30:18.450000Z", + "iopub.status.idle": "2026-09-03T22:30:21.544000Z", + "shell.execute_reply": "2026-09-03T22:30:21.544000Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "hs85" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + ": Compatibility = pass, Allocation = pass (0 bytes)\n", + "hs89" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + ": Compatibility = pass, Allocation = pass (0 bytes)\n" + ] + } + ], + "source": [ + "for problem in [\"hs85\", \"hs89\"]\n", + " dyn = check_problem_dynamic(problem)\n", + " println(problem, \": Compatibility = \", dyn[\"Compatibility\"], \", Allocation = \", dyn[\"Allocation\"])\n", + "end" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "5863d263", + "metadata": { + "execution": { + "iopub.execute_input": "2026-09-03T22:30:21.547000Z", + "iopub.status.busy": "2026-09-03T22:30:21.547000Z", + "iopub.status.idle": "2026-09-03T22:30:33.507000Z", + "shell.execute_reply": "2026-09-03T22:30:33.507000Z" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "******************************************************************************\n", + "This program contains Ipopt, a library for large-scale nonlinear optimization.\n", + " Ipopt is released as open source code under the Eclipse Public License (EPL).\n", + " For more information visit https://github.com/coin-or/Ipopt\n", + "******************************************************************************\n", + "\n", + "This is Ipopt version 3.14.19, running with linear solver MUMPS 5.8.2.\n", + "\n", + "Number of nonzeros in equality constraint Jacobian...: 0\n", + "Number of nonzeros in inequality constraint Jacobian.: 119\n", + "Number of nonzeros in Lagrangian Hessian.............: 276\n", + "\n", + "Total number of variables............................: 5\n", + " variables with only lower bounds: 0\n", + " variables with lower and upper bounds: 5\n", + " variables with only upper bounds: 0\n", + "Total number of equality constraints.................: 0\n", + "Total number of inequality constraints...............: 38\n", + " inequality constraints with only lower bounds: 20\n", + " inequality constraints with lower and upper bounds: 0\n", + " inequality constraints with only upper bounds: 18\n", + "\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 0 -9.3939688e-01 0.00e+00 8.89e-01 -1.0 0.00e+00 - 0.00e+00 0.00e+00 0\n", + " 1 -9.3945298e-01 0.00e+00 1.99e-02 -1.0 3.24e+00 - 9.90e-01 1.00e+00f 1\n", + " 2 -9.5732321e-01 0.00e+00 7.99e-04 -1.0 2.54e+02 - 9.78e-01 1.00e+00f 1\n", + " 3 -1.2162010e+00 0.00e+00 1.36e-01 -1.0 1.28e+03 - 7.51e-01 1.00e+00h 1\n", + " 4 -1.3019172e+00 0.00e+00 9.52e-03 -1.0 3.07e+03 - 9.86e-01 1.00e+00h 1\n", + " 5 -1.3363103e+00 0.00e+00 1.01e-04 -1.0 1.21e+03 - 1.00e+00 1.00e+00h 1\n", + " 6 -1.3443210e+00 0.00e+00 1.11e-06 -1.0 1.36e+02 - 1.00e+00 1.00e+00h 1\n", + " 7 -1.5109565e+00 0.00e+00 7.74e-01 -2.5 2.23e+03 - 7.51e-01 1.00e+00f 1\n", + " 8 -1.8685364e+00 4.98e-02 1.49e-01 -2.5 9.84e+03 - 6.03e-01 1.00e+00h 1\n", + " 9 -1.8858514e+00 0.00e+00 3.45e-03 -2.5 5.18e+02 - 1.00e+00 9.47e-01h 1\n", + "iter objective inf_pr inf_du lg(mu) ||d|| lg(rg) alpha_du alpha_pr ls\n", + " 10 -1.8780147e+00 0.00e+00 6.19e-05 -2.5 3.71e+02 - 1.00e+00 1.00e+00h 1\n", + " 11 -1.8799167e+00 0.00e+00 7.61e-06 -2.5 9.96e+01 - 1.00e+00 1.00e+00h 1\n", + " 12 -1.9007897e+00 0.00e+00 1.21e-01 -3.8 7.15e+02 - 9.57e-01 7.25e-01h 1\n", + " 13 -1.9045084e+00 4.95e-06 1.46e-02 -3.8 3.89e+01 - 1.00e+00 9.83e-01h 1\n", + " 14 -1.9042042e+00 0.00e+00 5.31e-08 -3.8 7.53e+00 - 1.00e+00 1.00e+00h 1\n", + " 15 -1.9051058e+00 0.00e+00 6.38e-03 -5.7 1.57e+01 - 9.74e-01 9.28e-01h 1\n", + " 16 -1.9051442e+00 0.00e+00 2.08e-09 -5.7 3.50e-01 - 1.00e+00 1.00e+00h 1\n", + " 17 -1.9051442e+00 0.00e+00 1.84e-11 -5.7 4.50e-03 - 1.00e+00 1.00e+00h 1\n", + " 18 -1.9051553e+00 0.00e+00 6.50e-09 -8.6 1.58e-01 - 1.00e+00 1.00e+00h 1\n", + "\n", + "Number of Iterations....: 18\n", + "\n", + " (scaled) (unscaled)\n", + "Objective...............: -1.9051552886872158e+00 -1.9051552886872158e+00\n", + "Dual infeasibility......: 6.5009405259171260e-09 6.5009405259171260e-09\n", + "Constraint violation....: 0.0000000000000000e+00 0.0000000000000000e+00\n", + "Variable bound violation: 6.0617809083396423e-07 6.0617809083396423e-07\n", + "Complementarity.........: 6.5213860383702599e-09 6.5213860383702599e-09\n", + "Overall NLP error.......: 6.5213860383702599e-09 6.5213860383702599e-09\n", + "\n", + "\n", + "Number of objective function evaluations = 19\n", + "Number of objective gradient evaluations = 19\n", + "Number of equality constraint evaluations = 0\n", + "Number of inequality constraint evaluations = 19\n", + "Number of equality constraint Jacobian evaluations = 0\n", + "Number of inequality constraint Jacobian evaluations = 19\n", + "Number of Lagrangian Hessian evaluations = 18\n", + "Total seconds in IPOPT = 6.413\n", + "\n", + "EXIT: Optimal Solution Found.\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Ipopt termination status: LOCALLY_SOLVED\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Objective value: -1.9051552886872158\n" + ] + } + ], + "source": [ + "try\n", + " using JuMP, Ipopt\n", + " model = OptimizationProblems.PureJuMP.hs85()\n", + " # PureJuMP constructors return a bare Model() with no optimizer attached\n", + " # (see docs/RULES_hs85_hs89.md / JuMP Interface Consistency); callers attach one.\n", + " set_optimizer(model, Ipopt.Optimizer)\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" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "julia 1.12", + "language": "julia", + "name": "julia" + }, + "language_info": { + "file_extension": ".jl", + "mimetype": "application/julia", + "name": "julia", + "version": "1.12.5" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/src/ADNLPProblems/hs85.jl b/src/ADNLPProblems/hs85.jl new file mode 100644 index 000000000..e9c666f18 --- /dev/null +++ b/src/ADNLPProblems/hs85.jl @@ -0,0 +1,205 @@ +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) + # 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, x) + # 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] + # 35 nonlinear constraints: cx[1:18] of the form c(x) >= 0, cx[19:35] of the form c(x) <= 0. + # The 3 remaining (linear) constraints of HS85 are passed separately below. + cx[1] = y2 - a[2] + cx[2] = y3 - a[3] + cx[3] = y4 - a[4] + cx[4] = y5 - a[5] + cx[5] = y6 - a[6] + cx[6] = y7 - a[7] + cx[7] = y8 - a[8] + cx[8] = y9 - a[9] + cx[9] = y10 - a[10] + cx[10] = y11 - a[11] + cx[11] = y12 - a[12] + cx[12] = y13 - a[13] + cx[13] = y14 - a[14] + cx[14] = y15 - a[15] + cx[15] = y16 - a[16] + cx[16] = y17 - a[17] + cx[17] = y4 - (T(0.28) / T(0.72)) * y5 + cx[18] = T(62212) / c17 - T(110.6) - y1 + cx[19] = y2 - b[2] + cx[20] = y3 - b[3] + cx[21] = y4 - b[4] + cx[22] = y5 - b[5] + cx[23] = y6 - b[6] + cx[24] = y7 - b[7] + cx[25] = y8 - b[8] + cx[26] = y9 - b[9] + cx[27] = y10 - b[10] + cx[28] = y11 - b[11] + cx[29] = y12 - b[12] + cx[30] = y13 - b[13] + cx[31] = y14 - b[14] + cx[32] = y15 - b[15] + cx[33] = y16 - b[16] + cx[34] = y17 - b[17] + cx[35] = T(3496) * y2 / c12 - T(21) + return cx + end + + # 3 linear constraints, given as a sparse Jacobian (rows, cols, vals): + # cx[1] = 1.5 x2 - x3 (>= 0), cx[2] = cx[3] = x2 + x3 (bounded below and above respectively) + clinrows = [1, 1, 2, 2, 3, 3] + clincols = [2, 3, 2, 3, 2, 3] + clinvals = T[1.5, -1, 1, 1, 1, 1] + + # Constraint bounds: overall cx[1:3] are the linear constraints above, + # cx[4:21] are of the form c(x) >= 0, cx[22:38] are of the form c(x) <= 0. + m = 38 + cl = zeros(T, m) + cu = fill(T(Inf), m) + cl[2], cu[3] = T(171.5), T(363.63) + cl[3] = T(-Inf) + for i = 22:38 + cl[i], cu[i] = T(-Inf), zero(T) + end + return ADNLPModels.ADNLPModel!( + f, + x0, + lvar, + uvar, + clinrows, + clincols, + clinvals, + c!, + cl, + cu; + name = "hs85", + kwargs..., + ) +end diff --git a/src/ADNLPProblems/hs89.jl b/src/ADNLPProblems/hs89.jl new file mode 100644 index 000000000..8ea0b1666 --- /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) + 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, x) + 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 000000000..ae4a2296f --- /dev/null +++ b/src/Meta/hs85.jl @@ -0,0 +1,25 @@ +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 = default_nvar, kwargs...) = 5 +get_hs85_ncon(; n::Integer = default_nvar, kwargs...) = 38 +get_hs85_nlin(; n::Integer = default_nvar, kwargs...) = 3 +get_hs85_nnln(; n::Integer = default_nvar, kwargs...) = 35 +get_hs85_nequ(; n::Integer = default_nvar, kwargs...) = 0 +get_hs85_nineq(; n::Integer = default_nvar, kwargs...) = 38 diff --git a/src/Meta/hs89.jl b/src/Meta/hs89.jl new file mode 100644 index 000000000..2dde53975 --- /dev/null +++ b/src/Meta/hs89.jl @@ -0,0 +1,25 @@ +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 = default_nvar, kwargs...) = 3 +get_hs89_ncon(; n::Integer = default_nvar, kwargs...) = 1 +get_hs89_nlin(; n::Integer = default_nvar, kwargs...) = 0 +get_hs89_nnln(; n::Integer = default_nvar, kwargs...) = 1 +get_hs89_nequ(; n::Integer = default_nvar, kwargs...) = 1 +get_hs89_nineq(; n::Integer = default_nvar, kwargs...) = 0 diff --git a/src/PureJuMP/hs85.jl b/src/PureJuMP/hs85.jl new file mode 100644 index 000000000..a03dca79d --- /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 defined with @expression + @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 + @objective( + 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 000000000..0a5d20238 --- /dev/null +++ b/src/PureJuMP/hs89.jl @@ -0,0 +1,105 @@ +# 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ⱼ ρⱼ + @objective(model, Min, sum(A[j] * ρ[j] for j = 1:30)) + + # Constraint: termA + termB = 2/15 + @constraint( + 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