diff --git a/.codecov.yml b/.codecov.yml new file mode 100644 index 00000000..69cb7601 --- /dev/null +++ b/.codecov.yml @@ -0,0 +1 @@ +comment: false diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index cae2d87d..6151a370 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -10,69 +10,13 @@ on: - main paths-ignore: - 'docs/**' +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ github.ref_name != github.event.repository.default_branch || github.ref != 'refs/tags/v*' }} jobs: - test: - runs-on: ubuntu-latest - continue-on-error: ${{ matrix.group == 'Downstream' }} - strategy: - fail-fast: false - matrix: - group: - - Corleone - - CorleoneOED - version: - - '1.11' - - 'lts' - steps: - - uses: actions/checkout@v6 - - uses: julia-actions/setup-julia@v2 - with: - version: ${{ matrix.version }} - - uses: actions/cache@v5 - env: - cache-name: cache-artifacts - with: - path: ~/.julia/artifacts - key: ${{ runner.os }}-test-${{ matrix.version }}-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }} - restore-keys: | - ${{ runner.os }}-test-${{ matrix.version }}-${{ env.cache-name }}- - ${{ runner.os }}-test-${{ matrix.version }}- - ${{ runner.os }}- - - uses: julia-actions/julia-buildpkg@v1 - - name: ${{ matrix.group }} - env: - GROUP: ${{ matrix.group }} - shell: julia --color=yes --check-bounds=yes --depwarn=yes {0} - run: | - using Pkg - const GROUP = get(ENV, "GROUP", "Core") - - function dev_subpkg(subpkgs::Vector{String}) - specs = [PackageSpec(path = "lib/$subpkg") for subpkg in subpkgs] - Pkg.develop(specs) - end - - if GROUP == "Corleone" - Pkg.activate(".") - else - subpkg_path = "lib/${{ matrix.group }}" - Pkg.activate(subpkg_path) - end - - if VERSION < v"1.11" - @info "Preparing env" - if GROUP == "Corleone" - @info "Testing Corleone" - elseif GROUP == "CorleoneOED" - dev_subpkg([".."]) - end - end - - @info "Starting tests" - Pkg.test() - - uses: julia-actions/julia-processcoverage@v1 - with: - directories: src,lib/CorleoneOED/src - - uses: codecov/codecov-action@v5 - with: - files: lcov.info \ No newline at end of file + tests: + uses: "SciML/.github/.github/workflows/grouped-tests.yml@v1" + with: + check-bounds: "auto" + coverage-directories: "src" + secrets: "inherit" diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml deleted file mode 100644 index cba9134c..00000000 --- a/.github/workflows/CompatHelper.yml +++ /dev/null @@ -1,16 +0,0 @@ -name: CompatHelper -on: - schedule: - - cron: 0 0 * * * - workflow_dispatch: -jobs: - CompatHelper: - runs-on: ubuntu-latest - steps: - - name: Pkg.add("CompatHelper") - run: julia -e 'using Pkg; Pkg.add("CompatHelper")' - - name: CompatHelper.main() - env: - GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} - COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }} - run: julia -e 'using CompatHelper; CompatHelper.main()' diff --git a/.github/workflows/DependabotAutoMerge.yml b/.github/workflows/DependabotAutoMerge.yml new file mode 100644 index 00000000..44b814c7 --- /dev/null +++ b/.github/workflows/DependabotAutoMerge.yml @@ -0,0 +1,9 @@ +name: "Dependabot Auto-merge" +on: pull_request +permissions: + contents: write + pull-requests: write +jobs: + automerge: + uses: "SciML/.github/.github/workflows/dependabot-automerge.yml@v1" + secrets: "inherit" diff --git a/.github/workflows/DocPreviewCleanup.yml b/.github/workflows/DocPreviewCleanup.yml new file mode 100644 index 00000000..7b694802 --- /dev/null +++ b/.github/workflows/DocPreviewCleanup.yml @@ -0,0 +1,17 @@ +name: Doc Preview Cleanup + +on: + pull_request: + types: [closed] + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +permissions: + contents: write + +jobs: + doc-preview-cleanup: + uses: "SciML/.github/.github/workflows/docs-preview-cleanup.yml@v1" + secrets: "inherit" diff --git a/.github/workflows/Downgrade.yml b/.github/workflows/Downgrade.yml index b4ef4779..b49343d8 100644 --- a/.github/workflows/Downgrade.yml +++ b/.github/workflows/Downgrade.yml @@ -12,21 +12,8 @@ on: - 'docs/**' jobs: test: - if: false # Disabled: waiting on dependency updates. See #586 - runs-on: ubuntu-latest - strategy: - matrix: - downgrade_mode: ['alldeps'] - julia-version: ['1.10'] - steps: - - uses: actions/checkout@v6 - - uses: julia-actions/setup-julia@v2 - with: - version: ${{ matrix.julia-version }} - - uses: julia-actions/julia-downgrade-compat@v2 - with: - skip: Pkg,TOML - - uses: julia-actions/julia-buildpkg@v1 - - uses: julia-actions/julia-runtest@v1 - with: - ALLOW_RERESOLVE: false \ No newline at end of file + name: "Downgrade" + uses: "SciML/.github/.github/workflows/downgrade.yml@v1" + with: + julia-version: "lts" + secrets: "inherit" diff --git a/.github/workflows/DowngradeSublibraries.yml b/.github/workflows/DowngradeSublibraries.yml new file mode 100644 index 00000000..8e8c4122 --- /dev/null +++ b/.github/workflows/DowngradeSublibraries.yml @@ -0,0 +1,24 @@ +name: Downgrade Sublibraries +on: + pull_request: + branches: + - main + paths-ignore: + - 'docs/**' + push: + branches: + - main + paths-ignore: + - 'docs/**' +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ github.ref != 'refs/heads/main' }} +jobs: + downgrade-sublibraries: + uses: "SciML/.github/.github/workflows/sublibrary-downgrade.yml@v1" + secrets: "inherit" + with: + julia-version: "lts" + group-env-name: "CORLEONE_TEST_GROUP" + group-env-value: "Core" + # Every lib/* sublibrary is downgrade-tested (projects auto-discovered). diff --git a/.github/workflows/FormatCheck.yml b/.github/workflows/FormatCheck.yml index 6bc43919..7ac23235 100644 --- a/.github/workflows/FormatCheck.yml +++ b/.github/workflows/FormatCheck.yml @@ -1,5 +1,4 @@ name: format-check - on: push: branches: @@ -8,12 +7,8 @@ on: - 'release-' tags: '*' pull_request: - jobs: runic: - runs-on: ubuntu-latest - steps: - - uses: actions/checkout@v4 - - uses: fredrikekre/runic-action@v1 - with: - version: '1' \ No newline at end of file + name: "Runic" + uses: "SciML/.github/.github/workflows/runic.yml@v1" + secrets: "inherit" diff --git a/.github/workflows/RunicSuggestions.yml b/.github/workflows/RunicSuggestions.yml new file mode 100644 index 00000000..43c07684 --- /dev/null +++ b/.github/workflows/RunicSuggestions.yml @@ -0,0 +1,9 @@ +name: "Runic Suggestions" +on: pull_request +permissions: + contents: read + pull-requests: write +jobs: + runic-suggestions: + uses: "SciML/.github/.github/workflows/runic-suggestions-on-pr.yml@v1" + secrets: "inherit" diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml index 746b039a..16d404b7 100644 --- a/.github/workflows/SpellCheck.yml +++ b/.github/workflows/SpellCheck.yml @@ -1,13 +1,7 @@ name: Spell Check - on: [pull_request] - jobs: typos-check: - name: Spell Check with Typos - runs-on: ubuntu-latest - steps: - - name: Checkout Actions Repository - uses: actions/checkout@v6 - - name: Check spelling - uses: crate-ci/typos@v1.18.0 \ No newline at end of file + name: "Spell Check with Typos" + uses: "SciML/.github/.github/workflows/spellcheck.yml@v1" + secrets: "inherit" diff --git a/.github/workflows/SublibraryCI.yml b/.github/workflows/SublibraryCI.yml new file mode 100644 index 00000000..32de868f --- /dev/null +++ b/.github/workflows/SublibraryCI.yml @@ -0,0 +1,25 @@ +name: Sublibrary CI + +on: + pull_request: + branches: + - main + paths-ignore: + - 'docs/**' + push: + branches: + - main + paths-ignore: + - 'docs/**' + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} + +jobs: + sublibraries: + uses: "SciML/.github/.github/workflows/sublibrary-project-tests.yml@v1" + with: + group-env-name: CORLEONE_TEST_GROUP + check-bounds: auto + secrets: "inherit" diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml index 623860f7..ba02ddae 100644 --- a/.github/workflows/TagBot.yml +++ b/.github/workflows/TagBot.yml @@ -1,15 +1,19 @@ -name: TagBot +name: "TagBot" on: issue_comment: - types: - - created + types: [created] workflow_dispatch: jobs: - TagBot: - if: github.event_name == 'workflow_dispatch' || github.actor == 'JuliaTagBot' - runs-on: ubuntu-latest - steps: - - uses: JuliaRegistries/TagBot@v1 - with: - token: ${{ secrets.GITHUB_TOKEN }} - ssh: ${{ secrets.DOCUMENTER_KEY }} \ No newline at end of file + tagbot: + uses: "SciML/.github/.github/workflows/tagbot.yml@v1" + secrets: "inherit" + tagbot-subpackages: + strategy: + fail-fast: false + matrix: + subdir: + - "lib/CorleoneOED" + uses: "SciML/.github/.github/workflows/tagbot.yml@v1" + with: + subdir: "${{ matrix.subdir }}" + secrets: "inherit" diff --git a/.typos.toml b/.typos.toml new file mode 100644 index 00000000..18f11359 --- /dev/null +++ b/.typos.toml @@ -0,0 +1,61 @@ +[default.extend-words] +# Julia-specific functions +indexin = "indexin" +findfirst = "findfirst" +findlast = "findlast" +eachindex = "eachindex" +setp = "setp" +getp = "getp" +setu = "setu" +getu = "getu" + +# Common variable patterns in Julia/SciML +osys = "osys" +rsys = "rsys" +usys = "usys" +eqs = "eqs" +rhs = "rhs" +lhs = "lhs" +ode = "ode" +oed = "oed" +ivp = "ivp" + +# Common abbreviations +tol = "tol" +rtol = "rtol" +atol = "atol" +idx = "idx" +prev = "prev" +curr = "curr" +init = "init" +tmp = "tmp" +vec = "vec" +arr = "arr" +dt = "dt" +du = "du" +dx = "dx" +dy = "dy" + +# Algorithm/type suffixes +alg = "alg" +prob = "prob" +sol = "sol" +cb = "cb" +opts = "opts" +args = "args" +kwargs = "kwargs" + +# Scientific abbreviations +ND = "ND" +nd = "nd" +MTK = "MTK" +ODE = "ODE" +OED = "OED" + +# Domain terms +jacobian = "jacobian" +hessian = "hessian" +discretization = "discretization" +parameterized = "parameterized" +discretized = "discretized" +vectorized = "vectorized" diff --git a/LICENSE b/LICENSE.md similarity index 100% rename from LICENSE rename to LICENSE.md diff --git a/Project.toml b/Project.toml index 69b10da8..698e3827 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Corleone" uuid = "2751e0db-33e9-4374-b36d-b7219e1e6b40" authors = ["Carl Julius Martensen ", "Christoph Plate ", "and contributors"] -version = "0.0.3" +version = "0.0.5" [deps] ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" @@ -30,50 +30,40 @@ CorleoneComponentArraysExtension = ["ComponentArrays"] CorleoneModelingToolkitExtension = "ModelingToolkit" [compat] -Aqua = "0.8" -ChainRulesCore = "1" +ChainRulesCore = "1.25.1" ComponentArrays = "0.15" Distributed = "1" -DocStringExtensions = "0.9" -ForwardDiff = "1.2" +DocStringExtensions = "0.9.3" Functors = "0.5" -Ipopt = "1.10" -JET = "0.9" LinearAlgebra = "1" LuxCore = "1.4" Makie = "0.24" ModelingToolkit = "11" OhMyThreads = "0.8" Optimization = "5" -OptimizationMOI = "1" -OrdinaryDiffEqBDF = "1" -OrdinaryDiffEqTsit5 = "1" +OrdinaryDiffEqTsit5 = "1, 2" +Pkg = "1" Random = "1" -RecursiveArrayTools = "3.30" -Reexport = "1.2" -SafeTestsets = "0.1" -SciMLBase = "2.110" -SciMLSensitivity = "7.87" +RecursiveArrayTools = "3.37.0" +Reexport = "1.2.2" +SafeTestsets = "0.1, 1" +StableRNGs = "1" +SciMLBase = "2.141.0, 3" SciMLStructures = "1.7" -SymbolicIndexingInterface = "0.3" +SciMLTesting = "1" +SymbolicIndexingInterface = "0.3.43" Test = "1" julia = "1.10" [extras] -Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" -ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" -Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" -JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" -Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" -OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1" -OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8" OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" -SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" +SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Aqua", "ComponentArrays", "Test", "OrdinaryDiffEqTsit5", "OrdinaryDiffEqBDF", "LinearAlgebra", "JET", "SciMLSensitivity", "Optimization", "OptimizationMOI", "Ipopt", "SafeTestsets", "ForwardDiff", "ModelingToolkit"] +test = ["ComponentArrays", "Test", "OrdinaryDiffEqTsit5", "LinearAlgebra", "SafeTestsets", "SciMLTesting", "Pkg", "StableRNGs"] diff --git a/docs/Project.toml b/docs/Project.toml index 0ee74323..c15a80af 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -9,13 +9,19 @@ DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +Juniper = "2ddba703-00a4-53a7-87a5-e8b9971dde84" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" LuxCore = "bb33d45b-7691-41d6-9220-0943567d0623" ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" OptimizationLBFGSB = "22f7324a-a79d-40f2-bebe-3af60c77bd15" OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1" +OptimalControlBenchmarks = "d147904a-1eb9-11f1-01e7-bdf88e0f9445" OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" + +[sources] +CorleoneOED = {path = "../lib/CorleoneOED"} +OptimalControlBenchmarks = {path = "../lib/OptimalControlBenchmarks"} diff --git a/docs/make.jl b/docs/make.jl index c1e7c766..e4cca33b 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,6 +1,6 @@ using Pkg -using Documenter, Corleone +using Documenter, Corleone, CorleoneOED, OptimalControlBenchmarks using DocumenterInterLinks using DocumenterCitations using Literate @@ -18,11 +18,11 @@ bib = CitationBibliography(joinpath(@__DIR__, "src", "assets", "bibliography.bib makedocs( sitename = "Corleone.jl", authors = "Carl Julius Martensen, Christoph Plate, et al.", - modules = [Corleone], + modules = [Corleone, CorleoneOED, OptimalControlBenchmarks], format = Documenter.HTML( assets = ["assets/favicon.ico"], canonical = "https://docs.sciml.ai/Corleone/stable/", - size_threshold = 1_000_000, # bytes + size_threshold = 1_500_000, # bytes ), #DocumenterVitepress.MarkdownVitepress( # repo = "github.com/SciML/Corleone.jl", diff --git a/docs/src/api.md b/docs/src/api.md index 9e9dc07b..75b0d01d 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -1,6 +1,22 @@ # API Reference +## Corleone + ```@autodocs Modules = [Corleone] Order = [:type, :function] ``` + +## CorleoneOED + +```@autodocs +Modules = [CorleoneOED] +Order = [:type, :function] +``` + +## OptimalControlBenchmarks + +```@autodocs +Modules = [OptimalControlBenchmarks] +Order = [:type, :function] +``` diff --git a/ext/MTKExtension/optimal_control.jl b/ext/MTKExtension/optimal_control.jl index ba14d176..fe16cf6d 100644 --- a/ext/MTKExtension/optimal_control.jl +++ b/ext/MTKExtension/optimal_control.jl @@ -111,7 +111,7 @@ function Corleone.CorleoneDynamicOptProblem( # We currently assume a single cost costbody = SymbolicUtils.Code.toexpr(substitute(only(newcosts), vars_subs)) - costfn = :(($(args_...),) -> $costbody) + costfn = Expr(:->, Expr(:tuple, args_...), costbody) costfun = @RuntimeGeneratedFunction( costfn diff --git a/lib/CorleoneOED/LICENSE b/lib/CorleoneOED/LICENSE.md similarity index 100% rename from lib/CorleoneOED/LICENSE rename to lib/CorleoneOED/LICENSE.md diff --git a/lib/CorleoneOED/Project.toml b/lib/CorleoneOED/Project.toml index ed93fa2a..77b8ed9a 100644 --- a/lib/CorleoneOED/Project.toml +++ b/lib/CorleoneOED/Project.toml @@ -1,7 +1,7 @@ name = "CorleoneOED" uuid = "6590a4b1-d036-41fe-a5b4-03a980244101" authors = ["Christoph Plate", "Carl Julius Martensen", "Sebastian Sager"] -version = "0.0.2" +version = "0.0.5" [deps] Corleone = "2751e0db-33e9-4374-b36d-b7219e1e6b40" @@ -13,7 +13,6 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226" -StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" @@ -27,41 +26,42 @@ CorleoneOEDOptimizationExtension = ["Optimization", "ComponentArrays"] CorleoneOEDZygoteExtension = "Zygote" [compat] -Aqua = "0.8" -Corleone = "0.0.3" +Corleone = "0.0.5" ComponentArrays = "0.15" -DocStringExtensions = "0.9" +DocStringExtensions = "0.9.3" ForwardDiff = "1.3" Ipopt = "1" LinearAlgebra = "1" LuxCore = "1.4" Optimization = "5" OptimizationMOI = "1" -OrdinaryDiffEqTsit5 = "1.5" +OrdinaryDiffEqTsit5 = "1.5, 2" +Pkg = "1" Random = "1" -Reexport = "1.2" +Reexport = "1.2.2" SafeTestsets = "0.1" -SciMLBase = "2" +SciMLBase = "2.141.0, 3" SciMLStructures = "1.7" StableRNGs = "1" -SymbolicIndexingInterface = "0.3" +SymbolicIndexingInterface = "0.3.43" Symbolics = "7" Test = "1" Zygote = "0.7" julia = "1.10" [extras] -Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1" OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [sources.Corleone] path = "../.." [targets] -test = ["Aqua", "ComponentArrays", "Ipopt", "Optimization", "OptimizationMOI", "OrdinaryDiffEqTsit5", "SafeTestsets", "Test"] +test = ["ComponentArrays", "Ipopt", "Optimization", "OptimizationMOI", "OrdinaryDiffEqTsit5", "Pkg", "SafeTestsets", "StableRNGs", "Test"] diff --git a/lib/CorleoneOED/README.md b/lib/CorleoneOED/README.md index 57549462..c70777d5 100644 --- a/lib/CorleoneOED/README.md +++ b/lib/CorleoneOED/README.md @@ -1,3 +1,10 @@ -# CorleoneOED +# CorleoneOED.jl -[![Build Status](https://github.com/AlCap23/CorleoneOED.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/AlCap23/CorleoneOED.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![Join the chat at https://julialang.zulipchat.com #sciml-bridged](https://img.shields.io/static/v1?label=Zulip&message=chat&color=9558b2&labelColor=389826)](https://julialang.zulipchat.com/#narrow/stream/279055-sciml-bridged) +[![Global Docs](https://img.shields.io/badge/docs-SciML-blue.svg)](https://docs.sciml.ai/Corleone/stable/) + +[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor%27s%20Guide-blueviolet)](https://github.com/SciML/ColPrac) +[![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle) + +CorleoneOED.jl is a component of the [Corleone.jl](https://github.com/SciML/Corleone.jl) monorepo. It provides optimal experimental design tooling, including Fisher information computation, sensitivity-based OED layers, and A/D/E optimality criteria. +While completely independent and usable on its own, users wanting the full optimal control and experimental design suite should use [Corleone.jl](https://github.com/SciML/Corleone.jl). diff --git a/lib/CorleoneOED/examples/compartmental_oed/main.jl b/lib/CorleoneOED/examples/compartmental_oed/main.jl index 4cb459d4..f2d7255f 100644 --- a/lib/CorleoneOED/examples/compartmental_oed/main.jl +++ b/lib/CorleoneOED/examples/compartmental_oed/main.jl @@ -112,3 +112,35 @@ optu = uopt.u + zero(ComponentArray(ps)) sol, _ = oed(nothing, optu, st) plot_oed(sol, optu.controls) + + +# ## Integer sampling decisions + +# Even though optimal solutions to OED problems tend to have integer solutions if the +# discretization grid is chosen adequately, sometimes it might be necessary to enforce +# integrality of the sampling decisions and include it in the problem formulation. This is +# also possible in Corleone. This is possible by passing `integer_weights = true` to the +# OptimizationProblem constructor. + +optprob = OptimizationProblem(oed, DCriterion(); M = [6.0], integer_weights = true) + +# This problem must then be solved by a suitable MINLP solver. In this case, we use Juniper. + +using Juniper + +opt = OptimizationMOI.MOI.OptimizerWithAttributes( + Juniper.Optimizer, + "nl_solver" => OptimizationMOI.MOI.OptimizerWithAttributes( + Ipopt.Optimizer, + "print_level" => 0, + ) +) + +sol = solve(optprob, opt) + +# As the relaxed solution was already integer, the same solution is obtained via Juniper. + +optu = uopt.u + zero(ComponentArray(ps)) +sol, _ = oed(nothing, optu, st) + +plot_oed(sol, optu.controls) diff --git a/lib/CorleoneOED/examples/lotka_multiexperiments/main.jl b/lib/CorleoneOED/examples/lotka_multiexperiments/main.jl index ea9a47b2..560a332e 100644 --- a/lib/CorleoneOED/examples/lotka_multiexperiments/main.jl +++ b/lib/CorleoneOED/examples/lotka_multiexperiments/main.jl @@ -117,6 +117,29 @@ optsol, _ = multi_exp(nothing, uopt + zero(ComponentArray(ps)), st) plot_experiments(optsol) +# ## Updating the FIM + +# As before in [the Lotka OED example](@ref lotka_oed), the initial FIM can be updated from +# the obtained solution. This way, the `MultiExperimentLayer` can be easily used in an iterative +# measurement campaign in which it is necessary to also consider previous experiments. + +previous_experiment = (ps = uopt.u + zero(ComponentArray(ps)), st = st) + +st_new = CorleoneOED.update_fim(multi_exp, [previous_experiment], st) +optprob = remake(optprob, p = st_new) + +uopt = solve( + optprob, Ipopt.Optimizer(), + tol = 1.0e-7, + hessian_approximation = "limited-memory", + max_iter = 300 +); + +# +optsol, _ = multi_exp(nothing, uopt + zero(ComponentArray(ps)), st_new) + +plot_experiments(optsol) + # ## Different parameters per experiment diff --git a/lib/CorleoneOED/examples/lotka_oed/main.jl b/lib/CorleoneOED/examples/lotka_oed/main.jl index a1121b08..d4a4a4e1 100644 --- a/lib/CorleoneOED/examples/lotka_oed/main.jl +++ b/lib/CorleoneOED/examples/lotka_oed/main.jl @@ -156,7 +156,7 @@ fim optprob = OptimizationProblem(oed, DCriterion(); M = [4.0, 4.0]) uopt = solve( optprob, Ipopt.Optimizer(), - tol = 1.0e-6, + tol = 1.0e-8, hessian_approximation = "limited-memory", max_iter = 300 ); @@ -165,6 +165,38 @@ uopt = solve( optsol, _ = oed(nothing, uopt + zero(ComponentArray(ps)), st) plot_oed(optsol) +# ## Updating the FIM with previously run experiments + +# Assuming that OED is used iteratively in a measurement campaign during which parameter estimation +# and OED are alternatingly used, one is interested in updating the FIM in the current +# iteration of OED with the previously run experiments, but for the current parameter values. +# CorleoneOED provides this functionality via the state `st` of the `OEDLayer` in which +# the initial value of the Fisher information matrix is stored. Per default, the FIM is +# initialized with zeros. + +st.F_init + +# However, with the solution from above it can be updated easily + +previous_experiment = (ps = uopt.u + zero(ComponentArray(ps)), st = st) + +st_new = CorleoneOED.update_fim(oed, [previous_experiment], st) +st_new.F_init + +# Now, the optimization problem can be remaked and solved with the updated state. + +optprob = remake(optprob, p = st_new) +uopt = solve( + optprob, Ipopt.Optimizer(), + tol = 1.0e-8, + hessian_approximation = "limited-memory", + max_iter = 300 +); + +# However, in this example, the solution is the same, as also the parameter values are unchanged. +optsol, _ = oed(nothing, uopt + zero(ComponentArray(ps)), st_new) +plot_oed(optsol) + # ## Sampling-free OED diff --git a/lib/CorleoneOED/ext/CorleoneOEDOptimizationExtension.jl b/lib/CorleoneOED/ext/CorleoneOEDOptimizationExtension.jl index c045b95a..2f54dae2 100644 --- a/lib/CorleoneOED/ext/CorleoneOEDOptimizationExtension.jl +++ b/lib/CorleoneOED/ext/CorleoneOEDOptimizationExtension.jl @@ -301,13 +301,68 @@ function setup_constraints(layer::MultiExperimentLayer{<:Any, <:Any, <:Any, <:Mu return sampling_cons end +function get_integrality(layer::OEDLayer{<:Any, true, false, <:MultipleShootingLayer}, integer::Bool, ps::ComponentArray, st::NamedTuple{fields}) where {fields} + int_u0 = Bool.(ps * 0) + !integer && return int_u0 + + for field in fields + _ps = getproperty(int_u0, field) + _st = getproperty(st, field) + subsets = CorleoneOED.__get_subsets(_st.active_controls, layer.sampling_indices) + for subset in subsets + _ps.controls[subset] .= true + end + end + return int_u0 +end + +function get_integrality(layer::OEDLayer{false, true, <:Any, <:SingleShootingLayer}, integer::Bool, ps::ComponentArray, st) + int_u0 = Bool.(ps * 0) + !integer && return int_u0 + subsets = CorleoneOED.__get_subsets(st.index_grid, layer.sampling_indices) + for subset in subsets + int_u0.controls[unique(subset)] .= true + end + return int_u0 +end + +function get_integrality(layer::OEDLayer{<:Any, true, <:Any, <:SingleShootingLayer}, integer::Bool, ps::ComponentArray, st) + int_u0 = Bool.(ps * 0) + !integer && return int_u0 + + subsets = CorleoneOED.__get_subsets(st.active_controls, layer.sampling_indices) + for subset in subsets + int_u0.controls[subset] .= true + end + return int_u0 +end + +function get_integrality(layer::MultiExperimentLayer{<:Any, <:Any, true}, integer::Bool, ps::ComponentArray, st::NamedTuple{fields}) where {fields} + int_u0 = Bool.(ps * 0) + !integer && return int_u0 + for (i, field) in enumerate(fields) + int_u0[field] = get_integrality(layer.layers[i], integer, getproperty(int_u0, field), getproperty(st, field)) + end + return int_u0 +end + +function get_integrality(layer::MultiExperimentLayer{<:Any, <:Any, false}, integer::Bool, ps::ComponentArray, st::NamedTuple{fields}) where {fields} + int_u0 = Bool.(ps * 0) + !integer && return int_u0 + for field in fields + int_u0[field] = get_integrality(layer.layers, integer, getproperty(int_u0, field), getproperty(st, field)) + end + return int_u0 +end + + function Optimization.OptimizationProblem( layer::Union{OEDLayer{<:Any, true, false}, MultiExperimentLayer{<:Any, false}}, crit::CorleoneOED.AbstractCriterion; AD::Optimization.ADTypes.AbstractADType = AutoForwardDiff(), u0::ComponentVector = ComponentArray(first(LuxCore.setup(Random.default_rng(), layer))), constraints::Union{Nothing, <:Dict{Any, <:NamedTuple{(:t, :bounds)}}, <:NamedTuple} = nothing, - variable_type::Type{T} = Float64, + variable_type::Type{T} = Float64, integer_weights = false, M = default_M(layer), kwargs... ) where {T} @@ -332,6 +387,8 @@ function Optimization.OptimizationProblem( @assert all(lb .<= u0 .<= ub) "The initial variables are not within the bounds. Please check the input!" + integer_sampling = get_integrality(layer, integer_weights, u0, st) + # Constraints cons = setup_constraints(layer, sol, constraints) lcons, ucons = extract_constraint_bounds(layer, constraints, M) @@ -342,17 +399,17 @@ function Optimization.OptimizationProblem( # Return the optimization problem return OptimizationProblem( opt_f, u0[:], st, lb = lb[:], ub = ub[:], - lcons = lcons, ucons = ucons, + lcons = lcons, ucons = ucons, int = integer_sampling[:] ) end function Optimization.OptimizationProblem( - layer::Union{OEDLayer{<:Any, true, true}, MultiExperimentLayer{<:Any, true}}, + layer::Union{OEDLayer{<:Any, true, true, <:SingleShootingLayer}, MultiExperimentLayer{<:Any, true}}, crit::CorleoneOED.AbstractCriterion; AD::Optimization.ADTypes.AbstractADType = AutoForwardDiff(), u0::ComponentVector = ComponentArray(first(LuxCore.setup(Random.default_rng(), layer))), constraints::Union{Nothing, <:Dict{Any, <:NamedTuple{(:t, :bounds)}}, <:NamedTuple} = nothing, - variable_type::Type{T} = Float64, + variable_type::Type{T} = Float64, integer_weights = false, M = default_M(layer), kwargs... ) where {T} @@ -387,6 +444,8 @@ function Optimization.OptimizationProblem( @assert all(lb .<= u0 .<= ub) "The initial variables are not within the bounds. Please check the input!" + integer_sampling = get_integrality(layer, integer_weights, u0, st) + # Constraints cons = setup_constraints(layer, sol, constraints) lcons, ucons = extract_constraint_bounds(layer, constraints, M) @@ -397,7 +456,7 @@ function Optimization.OptimizationProblem( # Return the optimization problem return OptimizationProblem( opt_f, u0[:], st, lb = lb[:], ub = ub[:], - lcons = lcons, ucons = ucons, + lcons = lcons, ucons = ucons, int = integer_sampling ) end diff --git a/lib/CorleoneOED/src/augmentation.jl b/lib/CorleoneOED/src/augmentation.jl index b04cd054..17f329b6 100644 --- a/lib/CorleoneOED/src/augmentation.jl +++ b/lib/CorleoneOED/src/augmentation.jl @@ -82,7 +82,6 @@ function compute_initial_F(prob::Union{ODEProblem, DAEProblem}, alg, config, par return F_tf_sens end - function compute_svd_of_F(prob, alg, config, params; ns = nothing, threshold_singular_values = 0.95, threshold_singular_vectors = 0.1, kwargs...) F = compute_initial_F(prob, alg, config, params) svdF = svd(F) diff --git a/lib/CorleoneOED/src/multiexperiments.jl b/lib/CorleoneOED/src/multiexperiments.jl index 930fd68d..fcd1b408 100644 --- a/lib/CorleoneOED/src/multiexperiments.jl +++ b/lib/CorleoneOED/src/multiexperiments.jl @@ -1,6 +1,8 @@ """ $(TYPEDEF) + Generalization of OEDLayer to multiple experiments that can be jointly optimized. + # Fields $(FIELDS) """ @@ -41,8 +43,28 @@ function Base.show(io::IO, oed::MultiExperimentLayer{DISCRETE, FIXED, SPLIT}) wh end end +""" +$(SIGNATURES) + +Constructs a multi-experiment OED layer from one differential equation problem. + +# Arguments +- `prob`: Differential equation problem shared by the experiments. +- `alg`: Differential equation solver algorithm. +- `nexp`: Number of experiments, or pass a vector of parameter-index vectors to split + parameters across experiments. +- `shooting_points`: Optional multiple-shooting points. + +# Keywords +- `params`: Parameter indices included in each experiment. +- `measurements`: Optional measurement controls. +- `observed`: Observation map `(u, p, t) -> y`. + +# Returns +A `MultiExperimentLayer` whose parameters and states are grouped by experiment. +""" function MultiExperimentLayer{DISCRETE}( - prob::DEProblem, alg::DEAlgorithm, nexp::Int; + prob::SciMLBase.AbstractDEProblem, alg::SciMLBase.AbstractDEAlgorithm, nexp::Int; params = eachindex(prob.p), measurements = [], observed = default_observed, kwargs... ) where {DISCRETE} layer = OEDLayer{DISCRETE}(prob, alg; params = params, measurements = measurements, observed = observed, kwargs...) @@ -50,7 +72,7 @@ function MultiExperimentLayer{DISCRETE}( return MultiExperimentLayer{DISCRETE, fixed, false, SingleShootingLayer, typeof(layer), typeof(params)}(layer, nexp, params) end -function MultiExperimentLayer{DISCRETE}(prob::DEProblem, alg::DEAlgorithm, params::Vector{Vector{Int64}}; measurements = [], observed = default_observed, kwargs...) where {DISCRETE} +function MultiExperimentLayer{DISCRETE}(prob::SciMLBase.AbstractDEProblem, alg::SciMLBase.AbstractDEAlgorithm, params::Vector{Vector{Int64}}; measurements = [], observed = default_observed, kwargs...) where {DISCRETE} nexp = length(params) layers = map(params) do param OEDLayer{DISCRETE}(prob, alg; params = param, measurements = measurements, observed = observed, kwargs...) @@ -66,12 +88,12 @@ function MultiExperimentLayer{DISCRETE}(prob::DEProblem, alg::DEAlgorithm, param return MultiExperimentLayer{DISCRETE, fixed, true, SingleShootingLayer, typeof(layers), typeof(new_params)}(layers, nexp, new_params) end -function MultiExperimentLayer{DISCRETE}(prob::DEProblem, alg::DEAlgorithm, shooting_points::AbstractVector{<:Real}, nexp::Int; params = eachindex(prob.p), measurements = [], observed = default_observed, kwargs...) where {DISCRETE} +function MultiExperimentLayer{DISCRETE}(prob::SciMLBase.AbstractDEProblem, alg::SciMLBase.AbstractDEAlgorithm, shooting_points::AbstractVector{<:Real}, nexp::Int; params = eachindex(prob.p), measurements = [], observed = default_observed, kwargs...) where {DISCRETE} layers = OEDLayer{DISCRETE}(prob, alg, shooting_points...; params = params, measurements = measurements, observed = observed, kwargs...) return MultiExperimentLayer{DISCRETE, false, false, MultipleShootingLayer, typeof(layers), typeof(params)}(layers, nexp, params) end -function MultiExperimentLayer{DISCRETE}(prob::DEProblem, alg::DEAlgorithm, shooting_points::AbstractVector{<:Real}, params::Vector{Vector{Int64}} = [eachindex(prob.p) for _ in 1:nexp]; measurements = [], observed = default_observed, kwargs...) where {DISCRETE} +function MultiExperimentLayer{DISCRETE}(prob::SciMLBase.AbstractDEProblem, alg::SciMLBase.AbstractDEAlgorithm, shooting_points::AbstractVector{<:Real}, params::Vector{Vector{Int64}} = [eachindex(prob.p) for _ in 1:nexp]; measurements = [], observed = default_observed, kwargs...) where {DISCRETE} nexp = length(params) layers = map(params) do param OEDLayer{DISCRETE}(prob, alg, shooting_points...; params = param, measurements = measurements, observed = observed, kwargs...) @@ -101,14 +123,43 @@ function LuxCore.initialparameters(rng::Random.AbstractRNG, multi::MultiExperime return NamedTuple{exp_names}(exp_ps) end -function LuxCore.initialstates(rng::Random.AbstractRNG, multi::MultiExperimentLayer{<:Any, <:Any, true}) +function LuxCore.initialstates(rng::Random.AbstractRNG, multi::MultiExperimentLayer{<:Any, <:Any, true, <:SingleShootingLayer}) exp_names = Tuple([Symbol("experiment_$i") for i in 1:multi.n_exp]) exp_ps = Tuple( map(1:multi.n_exp) do i LuxCore.initialstates(rng, multi.layers[i]) end ) - return NamedTuple{exp_names}(exp_ps) + np = length(multi.params.all) + + st1 = exp_ps[1] + F_init = zeros(eltype(exp_ps[1].F_init), np, np) + st1 = merge(st1, (; F_init = F_init)) + + new_sts = (st1, exp_ps[2:end]...) + + return NamedTuple{exp_names}(new_sts) +end + + +function LuxCore.initialstates(rng::Random.AbstractRNG, multi::MultiExperimentLayer{<:Any, <:Any, true, <:MultipleShootingLayer}) + exp_names = Tuple([Symbol("experiment_$i") for i in 1:multi.n_exp]) + exp_ps = Tuple( + map(1:multi.n_exp) do i + LuxCore.initialstates(rng, multi.layers[i]) + end + ) + np = length(multi.params.all) + + exp1 = exp_ps[1] + int1 = exp1[1] + F_init = zeros(eltype(int1.F_init), np, np) + int1 = merge(int1, (; F_init = F_init)) + st1 = merge(exp1, (; interval_1 = int1)) + + new_sts = (st1, exp_ps[2:end]...) + + return NamedTuple{exp_names}(new_sts) end function LuxCore.initialstates(rng::Random.AbstractRNG, multi::MultiExperimentLayer{<:Any, <:Any, false}) @@ -141,6 +192,33 @@ Corleone.get_number_of_shooting_constraints(multi::MultiExperimentLayer{<:Any, < Corleone.get_number_of_shooting_constraints(multi::MultiExperimentLayer{<:Any, <:Any, true, <:MultipleShootingLayer}) = sum(map(Corleone.get_number_of_shooting_constraints, multi.layers)) Corleone.get_number_of_shooting_constraints(multi::MultiExperimentLayer{<:Any, <:Any, <:Any, <:SingleShootingLayer}) = 0 +function update_fim(oed::MultiExperimentLayer{DISCRETE, FIXED, <:Any, <:SingleShootingLayer}, experiments, st::NamedTuple) where {DISCRETE, FIXED} + FIM = sum( + map(experiments) do experiment + fisher_information(oed, nothing, experiment.ps, experiment.st)[1] + end + ) + + st1 = getproperty(st, Symbol("experiment_1")) + st1 = merge(st1, (; F_init = FIM + st[1].F_init)) + + return merge(st, (; experiment_1 = st1)) +end + + +function update_fim(oed::MultiExperimentLayer{DISCRETE, FIXED, <:Any, <:MultipleShootingLayer}, experiments, st::NamedTuple) where {DISCRETE, FIXED} + FIM = sum( + map(experiments) do experiment + fisher_information(oed, nothing, experiment.ps, experiment.st)[1] + end + ) + + st1 = getproperty(st, Symbol("experiment_1")) + int1 = merge(st1.interval_1, (; F_init = FIM + st[1][1].F_init)) + st1 = merge(st1, (; interval_1 = int1)) + + return merge(st, (; experiment_1 = st1)) +end function get_sampling_sums(multi::MultiExperimentLayer{<:Any, <:Any, true}, x, ps, st::NamedTuple{fields}) where {fields} return reduce( @@ -198,17 +276,20 @@ function __fisher_information(multi::MultiExperimentLayer{<:Any, true, true}, tr return F end -function fisher_information(multi::MultiExperimentLayer{<:Any, <:Any, false}, x, ps, st::NamedTuple{fields}) where {fields} - return sum( - map(fields) do field - fisher_information(multi.layers, x, getproperty(ps, field), getproperty(st, field))[1] +function fisher_information(multi::MultiExperimentLayer{<:Any, <:Any, false}, x, ps, st::NamedTuple{fields}; add_initial = true) where {fields} + F = sum( + map(fields) do field + fisher_information(multi.layers, x, getproperty(ps, field), getproperty(st, field), add_initial = false)[1] end - ), st + ) + F_init = isa(multi.layers.layer, MultipleShootingLayer) ? st[1][1].F_init : st[1].F_init + add_initial && return F + F_init, st + return F, st end -function fisher_information(multi::MultiExperimentLayer{<:Any, <:Any, true}, x, ps, st::NamedTuple{fields}) where {fields} +function fisher_information(multi::MultiExperimentLayer{<:Any, <:Any, true}, x, ps, st::NamedTuple{fields}; add_initial = true) where {fields} fim = map(enumerate(fields)) do (i, field) - fisher_information(multi.layers[i], x, getproperty(ps, field), getproperty(st, field))[1] + fisher_information(multi.layers[i], x, getproperty(ps, field), getproperty(st, field); add_initial = false)[1] end np = length(multi.params.all) F = zeros(eltype(fim[1]), (np, np)) @@ -216,7 +297,8 @@ function fisher_information(multi::MultiExperimentLayer{<:Any, <:Any, true}, x, idxs = [multi.params.permutation[j] for j in multi.params.original[i]] F[idxs, idxs] .+= fimi end - + F_init = isa(multi.layers[1].layer, MultipleShootingLayer) ? st[1][1].F_init : st[1].F_init + add_initial && return F + F_init, st return F, st end @@ -268,18 +350,3 @@ function Corleone.get_block_structure(layer::MultiExperimentLayer{<:Any, <:Any, return block_structure end - -function Corleone.shooting_constraints(trajs::AbstractVector{<:Trajectory}) - return reduce(vcat, reduce(vcat, shooting_violations.(trajs))) -end - -function Corleone.shooting_constraints!(res::AbstractVector, trajs::AbstractVector{<:Trajectory}) - i = 0 - for traj in trajs - for subvec in traj.shooting, j in eachindex(subvec) - i += 1 - res[i] = subvec[j] - end - end - return res -end diff --git a/lib/CorleoneOED/src/oed.jl b/lib/CorleoneOED/src/oed.jl index 20470aac..344f8b35 100644 --- a/lib/CorleoneOED/src/oed.jl +++ b/lib/CorleoneOED/src/oed.jl @@ -20,6 +20,10 @@ default_observed = (u, p, t) -> u """ $(TYPEDEF) +Wraps a Corleone shooting layer with the augmented dynamics needed for optimal experimental +design. The `DISCRETE` type parameter selects discrete or continuous information accumulation, +and optional measurement controls can restrict where observations are collected. + # Fields $(FIELDS) """ @@ -54,17 +58,34 @@ function Base.show(io::IO, oed::OEDLayer{DISCRETE, SAMPLED, FIXED}) where {DISCR return Base.show(io, "text/plain", isa(layer, SingleShootingLayer) ? layer.problem : layer.layer.problem) end -function OEDLayer{DISCRETE}(prob::DEProblem, alg::DEAlgorithm; params = eachindex(prob.p), measurements = [], observed = default_observed, kwargs...) where {DISCRETE} +""" +$(SIGNATURES) + +Constructs an `OEDLayer` from a differential equation problem and solver algorithm. + +# Arguments +- `prob`: Differential equation problem whose parameters are considered for experimental design. +- `alg`: Differential equation solver algorithm. + +# Keywords +- `params`: Parameter indices included in the Fisher information matrix. +- `measurements`: Optional `ControlParameter`s that define measurement schedules. +- `observed`: Observation map `(u, p, t) -> y`. + +# Returns +An `OEDLayer` that can be initialized with `LuxCore.setup` and evaluated like the wrapped +Corleone layer. +""" +function OEDLayer{DISCRETE}(prob::SciMLBase.AbstractDEProblem, alg::SciMLBase.AbstractDEAlgorithm; params = eachindex(prob.p), measurements = [], observed = default_observed, kwargs...) where {DISCRETE} layer = SingleShootingLayer(prob, alg; kwargs...) return OEDLayer{DISCRETE}(layer; params = params, measurements = measurements, observed = observed, kwargs...) end -function OEDLayer{DISCRETE}(prob::DEProblem, alg::DEAlgorithm, shooting_points...; params = eachindex(prob.p), measurements = [], observed = default_observed, kwargs...) where {DISCRETE} +function OEDLayer{DISCRETE}(prob::SciMLBase.AbstractDEProblem, alg::SciMLBase.AbstractDEAlgorithm, shooting_points...; params = eachindex(prob.p), measurements = [], observed = default_observed, kwargs...) where {DISCRETE} layer = MultipleShootingLayer(prob, alg, shooting_points...; kwargs...) return OEDLayer{DISCRETE}(layer; params = params, measurements = measurements, observed = observed, kwargs...) end - function OEDLayer{DISCRETE}(layer::MultipleShootingLayer, args...; measurements = [], kwargs...) where {DISCRETE} (; problem, algorithm, controls, control_indices, tunable_ic, bounds_ic, state_initialization, bounds_p, parameter_initialization, quadrature_indices) = layer.layer @@ -142,11 +163,48 @@ function OEDLayer{DISCRETE}(layer::SingleShootingLayer, args...; measurements = return OEDLayer{DISCRETE, SAMPLED, FIXED, typeof(newlayer), typeof(observed)}(newlayer, observed, samplings) end +function update_fim(oed::OEDLayer{<:Any, SAMPLED, FIXED, <:SingleShootingLayer}, experiments, st::NamedTuple) where {SAMPLED, FIXED} + FIM = sum( + map(experiments) do experiment + fisher_information(oed, nothing, experiment.ps, experiment.st)[1] + end + ) + + return merge(st, (; F_init = FIM + st.F_init)) +end + +function update_fim(oed::OEDLayer{<:Any, SAMPLED, FIXED, <:MultipleShootingLayer}, experiments, st::NamedTuple) where {SAMPLED, FIXED} + FIM = sum( + map(experiments) do experiment + fisher_information(oed, nothing, experiment.ps, experiment.st)[1] + end + ) + + st1 = merge(st[1], (; F_init = FIM + st[1].F_init)) + + return merge(st, (; interval_1 = st1)) +end + n_observed(layer::OEDLayer) = length(layer.sampling_indices) Corleone.get_number_of_shooting_constraints(oed::OEDLayer{<:Any, <:Any, <:Any, <:MultipleShootingLayer}) = Corleone.get_number_of_shooting_constraints(oed.layer) Corleone.get_number_of_shooting_constraints(oed::OEDLayer{<:Any, <:Any, <:Any, <:SingleShootingLayer}) = 0 Corleone.get_bounds(oed::OEDLayer; kwargs...) = Corleone.get_bounds(oed.layer; kwargs...) +get_size_F(oed::OEDLayer{true, true, <:Any}) = begin + size_hxG = size(oed.observed.fisher.getters) + if length(size_hxG) > 2 + size_hxG = size_hxG[1:2] + end + return (size_hxG[2], size_hxG[2]) +end +get_size_F(oed::OEDLayer{false, true, true}) = begin + sizeF = size(oed.observed.fisher.getters) + if length(sizeF) > 2 + sizeF = sizeF[1:2] + end + return sizeF +end + # This is the only case where we need to sample the trajectory function LuxCore.initialstates(rng::Random.AbstractRNG, oed::Union{OEDLayer{true, true, <:Any, <:SingleShootingLayer}, OEDLayer{false, true, true}}) (; layer, sampling_indices) = oed @@ -181,10 +239,14 @@ function LuxCore.initialstates(rng::Random.AbstractRNG, oed::Union{OEDLayer{true end end : measurement_indices + T = eltype(problem.u0) + F_init = zeros(T, get_size_F(oed)) + return merge( st, (; observation_grid = WeightedObservation(weighting_grid), active_controls = measurement_indices, + F_init = F_init, ) ) end @@ -195,8 +257,10 @@ function LuxCore.initialstates(rng::Random.AbstractRNG, oed::OEDLayer{true, true st = LuxCore.initialstates(rng, layer) # Our goal is to build a weigthing matrix similar to the indexgrid grids = Corleone.get_timegrid.(controls) + T = eltype(problem.u0) + F_init = zeros(T, get_size_F(oed)) - return map(st) do sti + st_new = map(st) do sti tspan = (first(first(sti.tspans)), last(last(sti.tspans))) overall_grid = vcat(reduce(vcat, grids), collect(tspan)) unique!(sort!(overall_grid)) @@ -233,6 +297,34 @@ function LuxCore.initialstates(rng::Random.AbstractRNG, oed::OEDLayer{true, true ) ) end + + st1 = merge(st_new[1], (; F_init = F_init)) + return merge(st_new, (; interval_1 = st1)) +end + +get_problem(oed::OEDLayer{<:Any, <:Any, <:Any, <:SingleShootingLayer}) = oed.layer.problem +get_problem(oed::OEDLayer{<:Any, <:Any, <:Any, <:MultipleShootingLayer}) = oed.layer.layer.problem + +function LuxCore.initialstates(rng::Random.AbstractRNG, oed::OEDLayer{<:Any, <:Any, <:Any, <:SingleShootingLayer}) + (; layer, sampling_indices) = oed + st = LuxCore.initialstates(rng, layer) + problem = get_problem(oed) + T = eltype(problem.u0) + F_init = zeros(T, size(oed.observed.fisher.getters)) + + return merge(st, (; F_init = F_init)) +end + +function LuxCore.initialstates(rng::Random.AbstractRNG, oed::OEDLayer{<:Any, <:Any, <:Any, <:MultipleShootingLayer}) + (; layer, sampling_indices) = oed + st = LuxCore.initialstates(rng, layer) + problem = get_problem(oed) + T = eltype(problem.u0) + F_init = zeros(T, size(oed.observed.fisher.getters)) + + st1 = st.interval_1 + st1 = merge(st1, (; F_init = F_init)) + return merge(st, (; interval_1 = st1)) end __fisher_information(oed::OEDLayer, traj::Trajectory) = oed.observed.fisher(traj) @@ -269,59 +361,107 @@ function __fisher_information(oed::OEDLayer{true, true, true}, traj::Trajectory, return observation_grid(controls, Gs) end -fisher_information(oed::OEDLayer, x, ps, st::NamedTuple) = begin +""" +$(SIGNATURES) + +Computes the Fisher information matrix for an `OEDLayer` at parameters `ps` and state `st`. + +# Arguments +- `oed`: Optimal experimental design layer. +- `x`: External input passed through the Lux layer interface. +- `ps`: Layer parameters, typically from `LuxCore.initialparameters` or `LuxCore.setup`. +- `st`: Layer state, typically from `LuxCore.initialstates` or `LuxCore.setup`. + +# Keywords +- `add_initial`: Include previously accumulated information stored in `st`. + +# Returns +A tuple `(F, st)` containing the Fisher information matrix and updated layer state. +""" +fisher_information(oed::OEDLayer, x, ps, st::NamedTuple; add_initial = true) = begin traj, st = oed(x, ps, st) - sum(__fisher_information(oed, traj)), st + F = add_initial ? st.F_init + sum(__fisher_information(oed, traj)) : sum(__fisher_information(oed, traj)) + return F, st end # Continuous ALWAYS last FIM -fisher_information(oed::OEDLayer{false}, x, ps, st::NamedTuple) = begin +fisher_information(oed::OEDLayer{false, <:Any, <:Any, <:SingleShootingLayer}, x, ps, st::NamedTuple; add_initial = true) = begin traj, st = oed(x, ps, st) - last(__fisher_information(oed, traj)), st + F = add_initial ? st.F_init + last(__fisher_information(oed, traj)) : last(__fisher_information(oed, traj)) + return F, st +end + +fisher_information(oed::OEDLayer{false, <:Any, <:Any, <:MultipleShootingLayer}, x, ps, st::NamedTuple; add_initial = true) = begin + traj, st = oed(x, ps, st) + F = add_initial ? st[1].F_init + last(__fisher_information(oed, traj)) : last(__fisher_information(oed, traj)) + return F, st end # DISCRETE and SAMPLING -> weighted sum -fisher_information(oed::OEDLayer{true, true, false, <:SingleShootingLayer}, x, ps, st::NamedTuple) = begin +fisher_information(oed::OEDLayer{true, true, false, <:SingleShootingLayer}, x, ps, st::NamedTuple; add_initial = true) = begin (; observation_grid) = st traj, st = oed(x, ps, st) Gs = __fisher_information(oed, traj) - observation_grid(ps.controls, Gs), st + F = add_initial ? st.F_init + observation_grid(ps.controls, Gs) : observation_grid(ps.controls, Gs) + return F, st end # FIXED DISCRETE and SAMPLING -> use helper function -fisher_information(oed::OEDLayer{true, true, true}, x, ps, st::NamedTuple) = begin +fisher_information(oed::OEDLayer{true, true, true}, x, ps, st::NamedTuple; add_initial = true) = begin traj, st = oed(x, ps, st) - __fisher_information(oed, traj, ps, st), st + F = add_initial ? st.F_init + __fisher_information(oed, traj, ps, st) : __fisher_information(oed, traj, ps, st) + return F, st end -fisher_information(oed::OEDLayer{true, true, false, <:MultipleShootingLayer}, x, ps, st::NamedTuple) = begin +fisher_information(oed::OEDLayer{true, true, false, <:MultipleShootingLayer}, x, ps, st::NamedTuple; add_initial = true) = begin traj, st = oed(x, ps, st) Gs = __fisher_information(oed, traj, ps, st) - sum( - map(eachindex(Gs)) do i - psi, sti = getproperty(ps, Symbol("interval_$i")), getproperty(st, Symbol("interval_$i")) - sti.observation_grid(psi.controls, Gs[i]) + F = sum( + map(eachindex(Gs)) do i + psi, sti = getproperty(ps, Symbol("interval_$i")), getproperty(st, Symbol("interval_$i")) + sti.observation_grid(psi.controls, Gs[i]) end - ), st + ) + add_initial && return st[1].F_init + F, st + return F, st end # DISCRETE -> SUM -fisher_information(oed::OEDLayer{true, false}, x, ps, st::NamedTuple) = begin +fisher_information(oed::OEDLayer{true, false, <:Any, <:SingleShootingLayer}, x, ps, st::NamedTuple; add_initial = true) = begin + (; sampling_indices, layer) = oed + (; observation_grid) = st + traj, st = oed(x, ps, st) + F = add_initial ? st.F_init + sum(__fisher_information(oed, traj)) : sum(__fisher_information(oed, traj)) + return F, st +end + +fisher_information(oed::OEDLayer{true, false, <:Any, <:MultipleShootingLayer}, x, ps, st::NamedTuple; add_initial = true) = begin (; sampling_indices, layer) = oed (; observation_grid) = st traj, st = oed(x, ps, st) - sum(__fisher_information(oed, traj)), st + F = add_initial ? st[1].F_init + sum(__fisher_information(oed, traj)) : sum(__fisher_information(oed, traj)) + return F, st end # FIXED + CONTINUOUS -fisher_information(oed::OEDLayer{false, true, true}, x, ps, st::NamedTuple) = begin +fisher_information(oed::OEDLayer{false, true, true}, x, ps, st::NamedTuple; add_initial = true) = begin (; sampling_indices, layer) = oed (; observation_grid) = st traj, st = oed(x, ps, st) - __fisher_information(oed, traj, ps, st), st + F = add_initial ? st.F_init + __fisher_information(oed, traj, ps, st) : __fisher_information(oed, traj, ps, st) + return F, st end +""" +$(SIGNATURES) + +Returns the parameter sensitivities of the observed outputs for an OED trajectory. + +# Returns +For trajectory inputs, the sensitivity arrays are returned directly. For layer inputs, the +return value is `(sensitivities, st)` with the updated layer state. +""" sensitivities(oed::OEDLayer, traj::Trajectory) = oed.observed.sensitivities(traj) sensitivities(oed::OEDLayer, x, ps, st::NamedTuple) = begin @@ -329,6 +469,15 @@ sensitivities(oed::OEDLayer, x, ps, st::NamedTuple) = begin sensitivities(oed, traj), st end +""" +$(SIGNATURES) + +Evaluates the observation equations associated with an `OEDLayer`. + +# Returns +For trajectory inputs, observed values are returned directly. For layer inputs, the return +value is `(observed, st)` with the updated layer state. +""" observed_equations(oed::OEDLayer, traj::Trajectory) = oed.observed.observed(traj) observed_equations(oed::OEDLayer, x, ps, st::NamedTuple) = begin @@ -338,6 +487,16 @@ end _local_information_gain(oed::OEDLayer, traj::Trajectory) = oed.observed.local_weighted_sensitivity(traj) +""" +$(SIGNATURES) + +Computes local information gain contributions from each observed quantity along the OED +trajectory. + +# Returns +A tuple `(gains, st)` containing per-observation information matrices and the updated layer +state. +""" local_information_gain(oed::OEDLayer, x, ps, st::NamedTuple) = begin traj, st = oed(x, ps, st) # This returns hx G but stacked as a matrix [h_1_x G; h_2_x G; ...] @@ -350,6 +509,16 @@ local_information_gain(oed::OEDLayer, x, ps, st::NamedTuple) = begin end, st end +""" +$(SIGNATURES) + +Computes global information gain contributions scaled by the inverse final Fisher information +matrix. + +# Returns +A tuple `(gains, st)` containing per-observation information matrices and the updated layer +state. +""" global_information_gain(oed::OEDLayer, x, ps, st::NamedTuple) = begin traj, st = oed(x, ps, st) F_tf, st = fisher_information(oed, x, ps, st) @@ -366,6 +535,8 @@ end get_sampling_sums(::OEDLayer{<:Any, false}, x, ps, st) = [] get_sampling_sums!(res, ::OEDLayer{<:Any, false}, x, ps, st) = nothing +get_sampling_sums(::OEDLayer{<:Any, false, <:Any, <:Corleone.MultipleShootingLayer}, x, ps, st::NamedTuple{fields}) where {fields} = [] +get_sampling_sums!(res, ::OEDLayer{<:Any, false, <:Any, <:Corleone.MultipleShootingLayer}, x, ps, st::NamedTuple{fields}) where {fields} = nothing __get_subsets(active_controls::AbstractVector, indices) = active_controls[indices] __get_subsets(index_grid::AbstractMatrix, indices) = index_grid[indices, :] diff --git a/lib/CorleoneOED/test/1d_oed.jl b/lib/CorleoneOED/test/core/1d_oed.jl similarity index 100% rename from lib/CorleoneOED/test/1d_oed.jl rename to lib/CorleoneOED/test/core/1d_oed.jl diff --git a/lib/CorleoneOED/test/lotka_oed.jl b/lib/CorleoneOED/test/core/lotka_oed.jl similarity index 84% rename from lib/CorleoneOED/test/lotka_oed.jl rename to lib/CorleoneOED/test/core/lotka_oed.jl index 47ff923c..0065bab0 100644 --- a/lib/CorleoneOED/test/lotka_oed.jl +++ b/lib/CorleoneOED/test/core/lotka_oed.jl @@ -63,10 +63,47 @@ res = zeros(3) @test [12.0, 12.0, 0.0] == res @test CorleoneOED.get_sampling_sums(oed, nothing, ps, st) == res[1:2] @test_nowarn @inferred CorleoneOED.get_sampling_sums(oed, nothing, ps, st) +shooting_points = [0.0, 3.0, 6.0, 9.0] + +@testset "Integer measurements" begin + for discrete in [true, false] + for integer_sampling in [true, false] + for fixed in [true, false] + for shooting in [true, false] + + shooting && fixed && continue + + oed = shooting ? OEDLayer{discrete}( + prob, Tsit5(), shooting_points, + params = [2, 3], + controls = fixed ? [] : (1 => control,), + measurements = [ + ControlParameter(collect(tgrid1), controls = ones(length(tgrid1)), bounds = (0.0, 1.0)), + ControlParameter(collect(tgrid2), controls = ones(length(tgrid2)), bounds = (0.0, 1.0)), + ], + observed = (u, p, t) -> u[1:2], + ) : OEDLayer{discrete}( + prob, Tsit5(), + params = [2, 3], + controls = fixed ? [] : (1 => control,), + measurements = [ + ControlParameter(collect(tgrid1), controls = ones(length(tgrid1)), bounds = (0.0, 1.0)), + ControlParameter(collect(tgrid2), controls = ones(length(tgrid2)), bounds = (0.0, 1.0)), + ], + observed = (u, p, t) -> u[1:2], + ) + + optprob = OptimizationProblem(oed, ACriterion(), M = [4.0, 4.0], integer_weights = integer_sampling) + + @test sum(optprob.int) == (integer_sampling ? length(tgrid1) + length(tgrid2) : 0) + end + end + end + end +end @testset "Single Experiments" begin - shooting_points = [0.0, 3.0, 6.0, 9.0] multi_layer = MultipleShootingLayer(prob, Tsit5(), shooting_points..., controls = (1 => control,), bounds_p = ([1.0, 1.0], [1.0, 1.0])) for _layer in [layer, multi_layer] diff --git a/lib/CorleoneOED/test/lotka_oed_svd.jl b/lib/CorleoneOED/test/core/lotka_oed_svd.jl similarity index 100% rename from lib/CorleoneOED/test/lotka_oed_svd.jl rename to lib/CorleoneOED/test/core/lotka_oed_svd.jl diff --git a/lib/CorleoneOED/test/qa/Project.toml b/lib/CorleoneOED/test/qa/Project.toml new file mode 100644 index 00000000..616805df --- /dev/null +++ b/lib/CorleoneOED/test/qa/Project.toml @@ -0,0 +1,18 @@ +[deps] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +Corleone = "2751e0db-33e9-4374-b36d-b7219e1e6b40" +CorleoneOED = "6590a4b1-d036-41fe-a5b4-03a980244101" +SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[sources] +Corleone = {path = "../../../.."} +CorleoneOED = {path = "../.."} + +[compat] +Aqua = "0.8" +Corleone = "0.0.5" +CorleoneOED = "0.0.5" +SciMLTesting = "1.6" +Test = "1" +julia = "1.10" diff --git a/lib/CorleoneOED/test/qa/qa.jl b/lib/CorleoneOED/test/qa/qa.jl new file mode 100644 index 00000000..bc60355d --- /dev/null +++ b/lib/CorleoneOED/test/qa/qa.jl @@ -0,0 +1,30 @@ +using SciMLTesting +using CorleoneOED +using Test + +run_qa( + CorleoneOED; + explicit_imports = true, + # CorleoneOED pulls Corleone and its other deps in with bare `using`, so it + # leans on a large set of implicit imports. Converting every one to an + # explicit import is a sizable refactor tracked in SciML/Corleone.jl#103. + ei_broken = (:no_implicit_imports,), + ei_kwargs = (; + # Names still not declared public in their owning modules: SciMLBase + # internals (`AbstractDEAlgorithm`, `AbstractDEProblem`, `get_colorizers`), + # SciMLStructures internals (`Tunable`, `canonicalize`), SymbolicUtils + # internal (`Code`), Symbolics internals (`getdefaultval`, `setdefaultval`, + # `variables`), ForwardDiff internal (`jacobian`), and Corleone's own + # as-yet-unexported helpers reached through `Corleone.*`. + all_qualified_accesses_are_public = (; + ignore = ( + :AbstractDEAlgorithm, :AbstractDEProblem, :Code, :Tunable, + :build_index_grid, :canonicalize, :get_block_structure, + :get_bounds, :get_colorizers, + :get_number_of_shooting_constraints, :get_timegrid, :getdefaultval, + :jacobian, :retrieve_symbol_cache, :setdefaultval, + :shooting_constraints, :shooting_constraints!, :variables, + ), + ), + ), +) diff --git a/lib/CorleoneOED/test/runtests.jl b/lib/CorleoneOED/test/runtests.jl index 0095a017..45f8f18d 100644 --- a/lib/CorleoneOED/test/runtests.jl +++ b/lib/CorleoneOED/test/runtests.jl @@ -1,12 +1,42 @@ using CorleoneOED using SafeTestsets +using Pkg -@safetestset "1D Example" begin - include("1d_oed.jl") +# QA (Aqua) runs in an isolated environment (test/qa) so its tooling deps never +# enter the main test target's resolve. On Julia < 1.11 the [sources] table is +# ignored, so develop the package and its in-repo sibling by path. +function activate_qa_env() + Pkg.activate(joinpath(@__DIR__, "qa")) + Pkg.develop( + [ + Pkg.PackageSpec(path = joinpath(@__DIR__, "..")), + Pkg.PackageSpec(path = joinpath(@__DIR__, "..", "..", "..")) + ] + ) + return Pkg.instantiate() end -@safetestset "Lotka Volterra" begin - include("lotka_oed.jl") + +# Centralized sublibrary CI sets CORLEONE_TEST_GROUP to the bare package name +# (-> "Core") or "_" (-> ""). Fall back to GROUP, then "All", so +# local `Pkg.test()` runs (which set neither) run everything. Map the value to +# the standard Core/QA section names this file dispatches on. +const _G = get(ENV, "CORLEONE_TEST_GROUP", get(ENV, "GROUP", "All")) +const _SUB = "CorleoneOED" +const GROUP = _G == _SUB ? "Core" : (startswith(_G, _SUB * "_") ? _G[(length(_SUB) + 2):end] : _G) + +if GROUP == "All" || GROUP == "Core" + @safetestset "1D Example" begin + include("core/1d_oed.jl") + end + @safetestset "Lotka Volterra" begin + include("core/lotka_oed.jl") + end + @safetestset "Lotka Volterra SVD" begin + include("core/lotka_oed_svd.jl") + end end -@safetestset "Lotka Volterra SVD" begin - include("lotka_oed_svd.jl") + +if GROUP == "All" || GROUP == "QA" + activate_qa_env() + @safetestset "Code quality (Aqua.jl)" include("qa/qa.jl") end diff --git a/lib/CorleoneOED/test/test_groups.toml b/lib/CorleoneOED/test/test_groups.toml new file mode 100644 index 00000000..1fe84cd4 --- /dev/null +++ b/lib/CorleoneOED/test/test_groups.toml @@ -0,0 +1,5 @@ +[Core] +versions = ["lts", "1", "pre"] + +[QA] +versions = ["lts", "1"] diff --git a/lib/OptimalControlBenchmarks/LICENSE.md b/lib/OptimalControlBenchmarks/LICENSE.md new file mode 100644 index 00000000..d21503c7 --- /dev/null +++ b/lib/OptimalControlBenchmarks/LICENSE.md @@ -0,0 +1,21 @@ +MIT License + +Copyright (c) 2025 Robert Lampel, Julius Martensen, Christoph Plate, Reinhold Wittmann and contributors + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/lib/OptimalControlBenchmarks/Project.toml b/lib/OptimalControlBenchmarks/Project.toml new file mode 100644 index 00000000..812e3ede --- /dev/null +++ b/lib/OptimalControlBenchmarks/Project.toml @@ -0,0 +1,71 @@ +name = "OptimalControlBenchmarks" +uuid = "d147904a-1eb9-11f1-01e7-bdf88e0f9445" +authors = ["Robert Lampel", "Julius Martensen", "Christoph Plate", "Reinhold Wittmann"] +version = "0.0.3" + +[deps] +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" +Corleone = "2751e0db-33e9-4374-b36d-b7219e1e6b40" +DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LuxCore = "bb33d45b-7691-41d6-9220-0943567d0623" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +ModelingToolkitBase = "7771a370-6774-4173-bd38-47e70ca0b839" +Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" +OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +Reexport = "189a3867-3050-52da-a836-e630ba90ab69" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SciMLStructures = "53ae85a6-f571-4167-b2af-e1d143709226" +StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3" +SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" +Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" + +[compat] +BenchmarkTools = "1" +CairoMakie = "0.15" +ComponentArrays = "0.15" +Corleone = "0.0.5" +DocStringExtensions = "0.9.3" +ForwardDiff = "1.3" +IntervalSets = "0.7" +Ipopt = "1" +LinearAlgebra = "1" +LuxCore = "1.4" +ModelingToolkit = "11.2.0" +ModelingToolkitBase = "1" +Optimization = "5" +OptimizationMOI = "1" +OrdinaryDiffEqTsit5 = "1.5, 2" +Pkg = "1" +Random = "1" +Reexport = "1.2.2" +SafeTestsets = "0.1" +SciMLBase = "2.141.0, 3" +SciMLStructures = "1.7" +StableRNGs = "1" +SymbolicIndexingInterface = "0.3.43" +Symbolics = "7.2.0" +Test = "1" +julia = "1.10" + +[extras] +ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" +OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[sources.Corleone] +path = "../.." + +[targets] +test = ["Pkg", "SafeTestsets", "Test"] diff --git a/lib/OptimalControlBenchmarks/README.md b/lib/OptimalControlBenchmarks/README.md new file mode 100644 index 00000000..be8785c7 --- /dev/null +++ b/lib/OptimalControlBenchmarks/README.md @@ -0,0 +1,10 @@ +# OptimalControlBenchmarks.jl + +[![Join the chat at https://julialang.zulipchat.com #sciml-bridged](https://img.shields.io/static/v1?label=Zulip&message=chat&color=9558b2&labelColor=389826)](https://julialang.zulipchat.com/#narrow/stream/279055-sciml-bridged) +[![Global Docs](https://img.shields.io/badge/docs-SciML-blue.svg)](https://docs.sciml.ai/Corleone/stable/) + +[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor%27s%20Guide-blueviolet)](https://github.com/SciML/ColPrac) +[![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle) + +OptimalControlBenchmarks.jl is a component of the [Corleone.jl](https://github.com/SciML/Corleone.jl) monorepo. It collects optimal control benchmark problems and the tooling to run, scale, and plot them against Corleone.jl solvers. +While completely independent and usable on its own, users wanting the full optimal control suite should use [Corleone.jl](https://github.com/SciML/Corleone.jl). diff --git a/lib/OptimalControlBenchmarks/main.jl b/lib/OptimalControlBenchmarks/main.jl new file mode 100644 index 00000000..bed2fa0e --- /dev/null +++ b/lib/OptimalControlBenchmarks/main.jl @@ -0,0 +1,29 @@ +using OptimalControlBenchmarks +using Corleone +using UnoSolver +using Ipopt +using ModelingToolkit + +benchmarks = load_benchmarks() + +# choose the optimizer +# optimizer = UnoSolver.Optimizer() +optimizer = Ipopt.Optimizer() + +# settings for discretization, given as grids in [0.,1.] + +# the constraint grid is assumed to contain 0. and 1. +constraint_grid = collect(0.0:0.05:1.0) +control_grid = collect(0.0:0.05:1.0)[1:(end - 1)] +shooting_grid = collect(0.0:0.05:1.0) + +grids = OptimalControlBenchmarks.BenchmarkGrids( + constraint_grid, + control_grid, + shooting_grid +) + +# run all benchmarks from the `problems` folder +df = run_all(benchmarks, optimizer, grids) + +println(df) diff --git a/lib/OptimalControlBenchmarks/src/OptimalControlBenchmarks.jl b/lib/OptimalControlBenchmarks/src/OptimalControlBenchmarks.jl new file mode 100644 index 00000000..06f9401a --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/OptimalControlBenchmarks.jl @@ -0,0 +1,25 @@ +module OptimalControlBenchmarks + +using ModelingToolkit +using ModelingToolkit: t_nounits as t, D_nounits as D +using ModelingToolkit: inputs +using CairoMakie +using Corleone +using OrdinaryDiffEqTsit5 +using Optimization +using OptimizationMOI +using ForwardDiff +using ComponentArrays +using LuxCore, Random +using Symbolics + +include("types.jl") +include("plot_solution.jl") +include("problem_registry.jl") +include("run_benchmarks.jl") +include("scale_grid.jl") +include("solver_corleone.jl") + +export load_benchmarks, run_all + +end diff --git a/lib/OptimalControlBenchmarks/src/plot_solution.jl b/lib/OptimalControlBenchmarks/src/plot_solution.jl new file mode 100644 index 00000000..38e6bb4c --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/plot_solution.jl @@ -0,0 +1,10 @@ +function plot_oc_problem(sol, num_states, num_controls) + f = Figure() + ax = CairoMakie.Axis(f[1, 1]) + scatterlines!(ax, sol, vars = collect(1:1:num_states)) + f[1, 2] = Legend(f, ax, "States", framevisible = false) + ax1 = CairoMakie.Axis(f[2, 1]) + stairs!(ax1, sol, vars = collect((num_states + 1):1:(num_states + num_controls))) + f[2, 2] = Legend(f, ax1, "Controls", framevisible = false) + return f +end diff --git a/lib/OptimalControlBenchmarks/src/problem_registry.jl b/lib/OptimalControlBenchmarks/src/problem_registry.jl new file mode 100644 index 00000000..7719bbae --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problem_registry.jl @@ -0,0 +1,42 @@ +""" + load_benchmarks() + +Loads all benchmark problem constructors from the `problems` directory. + +# Returns +A vector of benchmark constructor functions. Each function accepts a benchmark grid +configuration and returns the corresponding optimal control problem data. +""" +function load_benchmarks() + benchmarks = [] + problems_dir = joinpath(@__DIR__, "problems") + + for f in readdir(problems_dir) + # Skip non-julia files + endswith(f, ".jl") || continue + + + # Get the path and the expected function name (stripping .jl) + path = joinpath(problems_dir, f) + func_name_sym = Symbol(splitext(f)[1]) + + # Evaluate the file content into the current module + include(path) + + # Retrieve the function object by its name + try + # Look up the symbol in the current module's scope + func_obj = getfield(@__MODULE__, func_name_sym) + + if func_obj isa Function + push!(benchmarks, func_obj) + else + @warn "File $f included, but $func_name_sym is not a function." + end + catch e + @error "Failed to load function $func_name_sym from $f" exception = e + end + end + + return benchmarks +end diff --git a/lib/OptimalControlBenchmarks/src/problems/batch_reactor.jl b/lib/OptimalControlBenchmarks/src/problems/batch_reactor.jl new file mode 100644 index 00000000..93d02988 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/batch_reactor.jl @@ -0,0 +1,35 @@ +function batch_reactor(grids) + + num_states = 2 + num_controls = 1 + tspan = (0.0, 1.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₁(..) = 1.0, [tunable = false] + x₂(..) = 0.0, [tunable = false] + u(..) = 298.0, [bounds = (298.0, 398.0), input = true] + end + + eqs = [ + D(x₁(t)) ~ -4.0e3 * exp(-2500 / u(t)) * x₁(t)^2 + D(x₂(t)) ~ 4.0e3 * exp(-2500 / u(t)) * x₁(t)^2 - 62.0e4 * exp(-5000 / u(t)) * x₂(t)^2 + ] + + costs = [-x₂(1.0)] + + @named oc_problem = System( + eqs, + t; + costs = costs + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Batch Reactor", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/bryson_denham.jl b/lib/OptimalControlBenchmarks/src/problems/bryson_denham.jl new file mode 100644 index 00000000..d9c7e858 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/bryson_denham.jl @@ -0,0 +1,48 @@ +function bryson_denham(grids) + + num_states = 2 + num_controls = 1 + tspan = (0.0, 1.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x(..) = 0.0, [tunable = false] + v(..) = 1.0, [tunable = false] + w(..) = 0.0, [input = true] + end + + eqs = [ + D(x(t)) ~ v(t) + D(v(t)) ~ w(t) + ] + + constraint_grid = scaled_grids.constraint_grid[1:(end - 1)] + + grid_cons = [x(tᵢ) ≲ 1 / 9 for tᵢ in constraint_grid] + + cons = [ + x(last(tspan)) ~ 0.0, + v(last(tspan)) ~ -1.0, + grid_cons..., + ] + + costs = [ + Symbolics.Integral(t in (0.0, 1.0))(w(t)^2), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Bryson Denham", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/cart_pendulum.jl b/lib/OptimalControlBenchmarks/src/problems/cart_pendulum.jl new file mode 100644 index 00000000..fd98b6a4 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/cart_pendulum.jl @@ -0,0 +1,63 @@ +function cart_pendulum(grids) + + num_states = 5 + num_controls = 1 + tspan = (0.0, 4.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x(..) = 0.0, [tunable = false] + θ(..) = 0.0, [tunable = false] + dx(..) = 0.0, [tunable = false] + dtheta(..) = 0.0, [tunable = false] + w(..) = 0.0, [bounds = (-30, 30), input = true] + end + + @parameters begin + α = 10.0, [tunable = false] + β = 50.0, [tunable = false] + γ = 0.5, [tunable = false] + M = 1.0, [tunable = false] + m = 0.1, [tunable = false] + g = 9.81, [tunable = false] + end + + eqs = [ + D(x(t)) ~ dx(t) + D(θ(t)) ~ dtheta(t) + D(dx(t)) ~ (w(t) + m * g * sin(θ(t)) * cos(θ(t)) + m * dtheta(t)^2 * sin(θ(t))) / (M + m * (1 - cos(θ(t)))^2) + D(dtheta(t)) ~ -g * sin(θ(t)) - ((w(t) + m * g * sin(θ(t)) * cos(θ(t)) + m * dtheta(t)^2 * sin(θ(t))) / (M + m * (1 - cos(θ(t))^2))) * cos(θ(t)) + ] + + constraint_grid = scaled_grids.constraint_grid + + grid_cons_le = [x(tᵢ) ≲ 2.0 for tᵢ in constraint_grid] + grid_cons_ge = [x(tᵢ) ≳ -2.0 for tᵢ in constraint_grid] + + cons = [ + grid_cons_le..., + grid_cons_ge..., + ] + + costs = [ + Symbolics.Integral(t in (0.0, 1.0))( + 1.0e-3 * (α * x(t)^2 + β * (θ(t) - pi)^2 + γ * w(t)^2) + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Cart Pendulum", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/catalyst_mixing.jl b/lib/OptimalControlBenchmarks/src/problems/catalyst_mixing.jl new file mode 100644 index 00000000..7cde346c --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/catalyst_mixing.jl @@ -0,0 +1,35 @@ +function catalyst_mixing(grids) + + num_states = 2 + num_controls = 1 + tspan = (0.0, 1.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₁(..) = 1.0, [tunable = false] + x₂(..) = 0.0, [tunable = false] + w(..) = 0.0, [bounds = (0.0, 1.0), input = true] + end + + eqs = [ + D(x₁(t)) ~ w(t) * (10 * x₂(t) - x₁(t)) + D(x₂(t)) ~ w(t) * (x₁(t) - 10 * x₂(t)) - (1 - w(t)) * x₂(t) + ] + + costs = [-1 + x₁(last(tspan)) + x₂(last(tspan))] + + @named oc_problem = System( + eqs, + t; + costs = costs + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Catalyst Mixing", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/cushioned_oscillation.jl b/lib/OptimalControlBenchmarks/src/problems/cushioned_oscillation.jl new file mode 100644 index 00000000..c6b54125 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/cushioned_oscillation.jl @@ -0,0 +1,50 @@ +function cushioned_oscillation(grids) + + num_states = 3 + num_controls = 1 + tspan = (0.0, 10.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x(..) = 2.0, [tunable = false] + v(..) = 5.0, [tunable = false] + obj(..) = 0.0, [tunable = false] + u(..) = 0.0, [bounds = (-5.0, 5.0), input = true] + end + @parameters begin + tₛ = 1.0, [bounds = (1.0e-3, Inf), tunable = true] + end + @constants begin + m = 5.0, [tunable = false] + c = 10.0, [tunable = false] + end + + eqs = [ + D(x(t)) ~ tₛ * v(t) + D(v(t)) ~ tₛ * (1 / m * (u(t) - c * x(t))) + D(obj(t)) ~ tₛ + ] + + cons = [ + x(last(tspan)) ~ 0.0, + v(last(tspan)) ~ 0.0, + ] + + costs = [obj(1.0)] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Cushioned Oscillation", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/denbigh.jl b/lib/OptimalControlBenchmarks/src/problems/denbigh.jl new file mode 100644 index 00000000..be466fd8 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/denbigh.jl @@ -0,0 +1,50 @@ +function denbigh(grids) + + # We start by defining our system + num_states = 3 + num_controls = 1 + tspan = (0.0, 1000.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₁(..) = 1.0, [tunable = false] + x₂(..) = 0.0, [tunable = false] + x₃(..) = 0.0, [tunable = false] + T(..) = 300.0, [bounds = (273.0, 415.0), input = true] + end + + @parameters begin + tₛ = 1.0, [bounds = (1.0e-3, Inf), tunable = true] + end + + @constants begin + E[1:4] = [3.0e3, 6.0e3, 3.0e3, 0.0], [tunable = false] + k⁰[1:4] = [1.0e3, 1.0e7, 1.0e1, 1.0e-3], [tunable = false] + end + + # auxiliary equations for the kᵢ + k = [k⁰[i] * exp(-E[i] / T(t)) for i in [1, 2, 3, 4]] + + eqs = [ + D(x₁(t)) ~ -k[1] * x₁(t) - k[2] * x₁(t) + D(x₂(t)) ~ k[1] * x₁(t) - k[3] + k[4] * x₂(t) + D(x₃(t)) ~ k[3] * x₂(t) + ] + + costs = [-x₃(last(tspan))] + + @named oc_problem = System( + eqs, + t; + costs = costs + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Denbigh Reaction", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/dielectrophoretic_particle.jl b/lib/OptimalControlBenchmarks/src/problems/dielectrophoretic_particle.jl new file mode 100644 index 00000000..be457c59 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/dielectrophoretic_particle.jl @@ -0,0 +1,49 @@ +function dielectrophoretic_particle(grids) + + num_states = 3 + num_controls = 1 + tspan = (0.0, 5.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₀(..) = 1.0, [tunable = false] + x₁(..) = 0.0, [tunable = false] + obj(..) = 0.0, [tunable = false] + u(..) = 1.0, [bounds = (-1.0, 1.0), input = true] + end + @parameters begin + tₛ = 1.0, [bounds = (1.0e-1, Inf), tunable = true] + end + @constants begin + α = -0.75, [tunable = false] + c = 1.0, [tunable = false] + end + + eqs = [ + D(x₀(t)) ~ tₛ * (x₁(t) * u(t) + α * u(t)^2) + D(x₁(t)) ~ tₛ * (-c * x₁(t) + u(t)) + D(obj(t)) ~ tₛ + ] + + cons = [ + x₀(last(tspan)) ~ 2.0, + ] + + costs = [obj(1.0)] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Dielectrophoretic Particle", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/double_oscillator.jl b/lib/OptimalControlBenchmarks/src/problems/double_oscillator.jl new file mode 100644 index 00000000..2f8ab0c1 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/double_oscillator.jl @@ -0,0 +1,52 @@ +function double_oscillator(grids) + + num_states = 5 + num_controls = 1 + tspan = (0.0, 2 * pi) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₀(..) = 0.0, [tunable = false] + x₁(..) = 0.0, [tunable = false] + x₂(..) = 0.0, [tunable = true] + x₃(..) = 0.0, [tunable = true] + u(..) = 1.0, [bounds = (-1.0, 1.0), input = true] + end + + @constants begin + m₁ = 200.0, [tunable = false] + m₂ = 2.0, [tunable = false] + k₁ = 100.0, [tunable = false] + k₂ = 3.0, [tunable = false] + c = 0.5, [tunable = false] + T = 2 * pi, [tunable = false] + end + + eqs = [ + D(x₀(t)) ~ x₂(t) + D(x₁(t)) ~ x₃(t) + D(x₂(t)) ~ - (k₁ + k₂) / m₁ * x₀(t) + k₂ / m₁ * x₁(t) + 1 / m₁ * sin(2 * pi / T * t) + D(x₃(t)) ~ k₂ / m₂ * x₀(t) - k₂ / m₂ * x₁(t) - c * (1 - u(t)) / m₂ * x₃(t) + ] + + costs = [ + Symbolics.Integral(t in (0.0, 1.0))( + 0.5 * (x₀(t)^2 + x₁(t)^2 + u(t)^2) + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Double Oscillator", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/ducted_fan.jl b/lib/OptimalControlBenchmarks/src/problems/ducted_fan.jl new file mode 100644 index 00000000..6509d76c --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/ducted_fan.jl @@ -0,0 +1,78 @@ +function ducted_fan(grids) + + num_states = 7 + num_controls = 2 + tspan = (0.0, 10.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₁(..) = 0.0, [tunable = false] + v₁(..) = 0.0, [tunable = false] + x₂(..) = 0.0, [tunable = false] + v₂(..) = 0.0, [tunable = false] + α(..) = 0.0, [tunable = false] + vα(..) = 0.0, [tunable = false] + u₁(..) = 0.0, [bounds = (-5.0, 5.0), input = true] + u₂(..) = 8.5, [bounds = (0.0, 17.0), input = true] + end + + @parameters begin + tₛ = 0.5, [bounds = (1.0e-3, Inf), tunable = true] + end + + @constants begin + m = 2.2, [tunable = false] + J = 0.05, [tunable = false] + r = 0.2, [tunable = false] + mg = 4.0, [tunable = false] + μ = 1.0, [tunable = false] + end + + eqs = [ + D(x₁(t)) ~ tₛ * (v₁(t)) + D(v₁(t)) ~ tₛ * (1 / m * (u₁(t) * cos(α(t)) - u₂(t) * sin(α(t)))) + D(x₂(t)) ~ tₛ * (v₂(t)) + D(v₂(t)) ~ tₛ * (1 / m * (-mg + u₁(t) * sin(α(t)) + u₂(t) * cos(α(t)))) + D(α(t)) ~ tₛ * (vα(t)) + D(vα(t)) ~ tₛ * (r / J * u₁(t)) + ] + + # scale the constraint grid + constraint_grid = scaled_grids.constraint_grid[1:(end - 1)] + + grid_cons_u = [α(tᵢ) ≲ 30.0 for tᵢ in constraint_grid] + grid_cons_l = [α(tᵢ) ≳ -30.0 for tᵢ in constraint_grid] + + cons = [ + x₁(last(tspan)) ~ 1.0, + v₁(last(tspan)) ~ 0.0, + x₂(last(tspan)) ~ 0.0, + v₂(last(tspan)) ~ 0.0, + α(last(tspan)) ~ 0.0, + vα(last(tspan)) ~ 0.0, + grid_cons_u..., + grid_cons_l..., + ] + + costs = [ + Symbolics.Integral(t in (0.0, last(tspan)))( + (2 * u₁(t)^2 + u₂(t)^2) + μ * tₛ + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Ducted Fan", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/egerstedt.jl b/lib/OptimalControlBenchmarks/src/problems/egerstedt.jl new file mode 100644 index 00000000..c51a0bd0 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/egerstedt.jl @@ -0,0 +1,53 @@ +function egerstedt(grids) + + num_states = 3 + num_controls = 3 + tspan = (0.0, 1.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x(..) = 0.5, [tunable = false] + y(..) = 0.5, [tunable = false] + u_1(..) = 1 / 3, [bounds = (-1.0, 1.0), input = true] + u_2(..) = 1 / 3, [bounds = (-1.0, 1.0), input = true] + u_3(..) = 1 / 3, [bounds = (-1.0, 1.0), input = true] + end + + eqs = [ + D(x(t)) ~ -x(t) * u_1(t) + (x(t) + y(t)) * u_2(t) + (x(t) - y(t)) * u_3(t) + D(y(t)) ~ (x(t) + 2 * y(t)) * u_1(t) + (x(t) - 2 * y(t)) * u_2(t) + (x(t) + y(t)) * u_3(t) + ] + + # scale the constraint grid + constraint_grid = scaled_grids.constraint_grid + + grid_cons_x = [y(tᵢ) ≳ 0.4 for tᵢ in constraint_grid] + grid_cons_u = [u_1(tᵢ) + u_2(tᵢ) + u_3(tᵢ) ~ 1.0 for tᵢ in constraint_grid] + + cons = [ + grid_cons_x..., + grid_cons_u..., + ] + + costs = [ + Symbolics.Integral(t in (0.0, 1.0))( + x(t)^2 + y(t)^2 + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Egerstedt standard problem", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/electric_car.jl b/lib/OptimalControlBenchmarks/src/problems/electric_car.jl new file mode 100644 index 00000000..65829e09 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/electric_car.jl @@ -0,0 +1,73 @@ +function electric_car(grids) + + num_states = 4 + num_controls = 1 + tspan = (0.0, 10.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x_1(..) = 0.0, [tunable = false] + x_2(..) = 0.0, [tunable = false] + x_3(..) = 0.0, [tunable = false] + u(..) = 0.0, [bounds = (-5.0, 5.0), input = true] + end + + @constants begin + K_r = 10.0, [tunable = false] + rho = 1.293, [tunable = false] + C_x = 0.4, [tunable = false] + S = 2.0, [tunable = false] + r = 0.33, [tunable = false] + K_f = 0.03, [tunable = false] + K_m = 0.27, [tunable = false] + R_m = 0.03, [tunable = false] + L_m = 0.05, [tunable = false] + M = 250.0, [tunable = false] + g = 9.81, [tunable = false] + V_alim = 150.0, [tunable = false] + R_bat = 0.05, [tunable = false] + end + + temp = K_m * x_1(t) + temp += - r / K_r * (M * g * K_f + 0.5 * rho * S * C_x * r^2 / K_r^2 * x_2(t)^2) + temp *= K_r^2 / (M * r^2) + + eqs = [ + D(x_1(t)) ~ (V_alim * u(t) - R_m * x_1(t) - K_m * x_2(t)) / L_m + D(x_2(t)) ~ temp + D(x_3(t)) ~ r / K_r * x_2(t) + ] + + constraint_grid = scaled_grids.constraint_grid + + grid_cons_u = [x_1(tᵢ) ≲ 150.0 for tᵢ in constraint_grid] + grid_cons_l = [x_1(tᵢ) ≳ -150.0 for tᵢ in constraint_grid] + + cons = [ + x_3(last(tspan)) ~ 100.0, + grid_cons_u..., + grid_cons_l..., + ] + + costs = [ + Symbolics.Integral(t in (0.0, last(tspan)))( + V_alim * u(t) * x_1(t) + R_bat * x_1(t)^2 + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Electric Car", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/f8.jl b/lib/OptimalControlBenchmarks/src/problems/f8.jl new file mode 100644 index 00000000..38a1471e --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/f8.jl @@ -0,0 +1,63 @@ +function f8(grids) + + num_states = 4 + num_controls = 1 + tspan = (0.0, 10.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x_0(..) = 0.4655, [tunable = false] + x_1(..) = 0.0, [tunable = false] + x_2(..) = 0.0, [tunable = false] + obj(..) = 0.0, [tunable = false] + u(..) = 0.5, [bounds = (0.0, 1.0), input = true] + end + + @parameters begin + T = 0.5, [bounds = (1.0e-3, Inf), tunable = true] + end + + @constants begin + ξ = 0.05236, [tunable = false] + end + + du_1 = -0.877 * x_0(t) + x_2(t) - 0.088 * x_0(t) * x_2(t) + 0.47 * x_0(t)^2 - 0.019 * x_1(t)^2 + du_1 += - x_0(t)^2 * x_2(t) + 3.846 * x_0(t)^3 + du_1 += 0.215 * ξ - 0.28 * x_0(t)^2 * ξ + 0.47 * x_0(t) * ξ^2 - 0.63 * ξ^3 + du_1 += - (0.215 * ξ - 0.28 * x_0(t)^2 * ξ - 0.63 * ξ^3) * 2 * u(t) + du_3 = -4.208 * x_0(t) - 0.396 * x_2(t) - 0.47 * x_0(t)^2 - 3.564 * x_0(t)^3 + du_3 += 20.967 * ξ - 6.265 * x_0(t)^2 * ξ + 46.0 * x_0(t) * ξ^2 - 61.4 * ξ^3 + du_3 += -(20.967 * ξ - 6.265 * x_0(t)^2 * ξ - 61.4 * ξ^3) * 2 * u(t) + + + eqs = [ + D(x_0(t)) ~ T * du_1 + D(x_1(t)) ~ T * x_2(t) + D(x_2(t)) ~ T * du_3 + D(obj(t)) ~ T + ] + + cons = [ + x_0(last(tspan)) ~ 0.0, + x_1(last(tspan)) ~ 0.0, + x_2(last(tspan)) ~ 0.0, + ] + + costs = [obj(last(tspan))] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "F-8 aircraft", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/fuller.jl b/lib/OptimalControlBenchmarks/src/problems/fuller.jl new file mode 100644 index 00000000..7b1209af --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/fuller.jl @@ -0,0 +1,45 @@ +function fuller(grids) + + num_states = 3 + num_controls = 1 + tspan = (0.0, 1.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₀(..) = 0.01, [tunable = false] + x₁(..) = 0.0, [tunable = false] + u(..) = 0.5, [bounds = (0.0, 1.0), input = true] + end + + eqs = [ + D(x₀(t)) ~ x₁(t) + D(x₁(t)) ~ 1.0 - 2.0 * u(t) + ] + + cons = [ + x₀(last(tspan)) ~ 0.01, + x₁(last(tspan)) ~ 0.0, + ] + + costs = [ + Symbolics.Integral(t in (0.0, last(tspan)))( + x₀(t)^2 + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Fuller's problem", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/goddarts_rocket.jl b/lib/OptimalControlBenchmarks/src/problems/goddarts_rocket.jl new file mode 100644 index 00000000..05c841de --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/goddarts_rocket.jl @@ -0,0 +1,65 @@ +function goddarts_rocket(grids) + + num_states = 3 + num_controls = 1 + tspan = (0.0, 10.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + r(..) = 1.0, [tunable = false] + v(..) = 0.0, [tunable = false] + m(..) = 1.0, [tunable = false] + u(..) = 0.5, [bounds = (0.0, 1.0), input = true] + end + + @parameters begin + T = 0.5, [bounds = (1.0e-3, Inf), tunable = true] + end + + @constants begin + r_0 = 1.0, [tunable = false] + v_0 = 0.0, [tunable = false] + m_0 = 1.0, [tunable = false] + r_T = 1.01, [tunable = false] + b = 7.0, [tunable = false] + T_max = 3.5, [tunable = false] + A = 310.0, [tunable = false] + k = 500.0, [tunable = false] + C = 0.6, [tunable = false] + end + + Drag = A * v(t)^2 * exp(-k * (r(t) - r_0)) + + eqs = [ + D(r(t)) ~ T * v(t) + D(v(t)) ~ T * (-1 / r(t)^2 + 1 / m(t) * (T_max * u(t) - Drag)) + D(m(t)) ~ T * (-b * u(t)) + ] + + constraint_grid = scaled_grids.constraint_grid + + grid_cons_u = [310.0 * v(tᵢ)^2 * 2.7182818^(-500.0 * (r(tᵢ) - 1.0)) ≲ 0.6 for tᵢ in constraint_grid] + + cons = [ + r(last(tspan)) ~ 1.01, + grid_cons_u..., + ] + + costs = [-m(last(tspan))] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Goddart's rocket problem", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/hang_glider.jl b/lib/OptimalControlBenchmarks/src/problems/hang_glider.jl new file mode 100644 index 00000000..25e5a825 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/hang_glider.jl @@ -0,0 +1,77 @@ +function hang_glider(grids) + + num_states = 4 + num_controls = 1 + tspan = (0.0, 100.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x(..) = 0.0, [tunable = false] + y(..) = 1.0e3, [tunable = false] + v_x(..) = 13.23, [tunable = false] + v_y(..) = -1.288, [tunable = false] + c_L(..) = 0.5, [bounds = (0.0, 1.4), input = true] + end + + @parameters begin + T = 0.5, [bounds = (1.0e-3, Inf), tunable = true] + end + + @constants begin + u_c = 2.5, [tunable = false] + r_c = 100.0, [tunable = false] + c_0 = 0.034, [tunable = false] + c_1 = 0.069662, [tunable = false] + S = 14.0, [tunable = false] + rho = 1.13, [tunable = false] + m = 100.0, [tunable = false] + g = 9.81, [tunable = false] + end + + # auxiliary functions + r = (x(t) / r_c - 2.5)^2 + U_up = u_c * (1 - r) * exp(-r) + w = v_y(t) - U_up + v = sqrt(v_x(t)^2 + w^2) + Dr = 0.5 * rho * S * (c_0 + c_1 * c_L(t)^2) * v^2 + L = 0.5 * rho * S * c_L(t) * v^2 + + eqs = [ + D(x(t)) ~ T * v_x(t) + D(y(t)) ~ T * v_y(t) + D(v_x(t)) ~ T * (-1.0) * (L * w + Dr * v_x(t)) / (m * v) + D(v_y(t)) ~ T * ((L * v_x(t) - Dr * w) / (m * v) - g) + ] + + # scale the constraint grid + constraint_grid = scaled_grids.constraint_grid[1:(end - 1)] + + grid_cons_x = [x(tᵢ) ≳ 0.0 for tᵢ in constraint_grid] + grid_cons_vx = [v_x(tᵢ) ≳ 0.0 for tᵢ in constraint_grid] + + cons = [ + y(last(tspan)) ~ 900.0, + v_x(last(tspan)) ~ 13.23, + v_y(last(tspan)) ~ -1.288, + grid_cons_x..., + grid_cons_vx..., + ] + + costs = [-x(last(tspan))] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Hang Glider", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/hanging_chain.jl b/lib/OptimalControlBenchmarks/src/problems/hanging_chain.jl new file mode 100644 index 00000000..c8d8d3f3 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/hanging_chain.jl @@ -0,0 +1,60 @@ +function hanging_chain(grids) + + num_states = 3 + num_controls = 1 + tspan = (0.0, 1.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₁(..) = 1.0, [tunable = false] + x₂(..) = 0.0, [tunable = false] + x₃(..) = 0.0, [tunable = false] + u(..) = 0.5, [bounds = (-10.0, 20.0), input = true] + end + + @constants begin + a = 1.0, [tunable = false] + b = 3.0, [tunable = false] + L_p = 4.0, [tunable = false] + end + + eqs = [ + D(x₁(t)) ~ u(t) + D(x₂(t)) ~ x₁(t) * sqrt(1 + u(t)^2) + D(x₃(t)) ~ sqrt(1 + u(t)^2) + ] + + # scale the constraint grid + constraint_grid = scaled_grids.constraint_grid + + grid_cons_u1 = [x₁(tᵢ) ≳ 0.0 for tᵢ in constraint_grid[1:(end - 1)]] + grid_cons_u2 = [x₂(tᵢ) ≳ 0.0 for tᵢ in constraint_grid] + grid_cons_u3 = [x₃(tᵢ) ≳ 0.0 for tᵢ in constraint_grid[1:(end - 1)]] + grid_cons_l1 = [x₁(tᵢ) ≲ 10.0 for tᵢ in constraint_grid[1:(end - 1)]] + grid_cons_l2 = [x₂(tᵢ) ≲ 10.0 for tᵢ in constraint_grid] + grid_cons_l3 = [x₃(tᵢ) ≲ 10.0 for tᵢ in constraint_grid[1:(end - 1)]] + + cons = [ + x₁(last(tspan)) - 3.0 ~ 0.0, + x₃(last(tspan)) - 4.0 ~ 0.0, + grid_cons_u1..., grid_cons_u2..., grid_cons_u3..., + grid_cons_l1..., grid_cons_l2..., grid_cons_l3..., + ] + + costs = [x₂(last(tspan))] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/jackson.jl b/lib/OptimalControlBenchmarks/src/problems/jackson.jl new file mode 100644 index 00000000..50c60660 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/jackson.jl @@ -0,0 +1,59 @@ +function jackson(grids) + + num_states = 3 + num_controls = 1 + tspan = (0.0, 1.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₁(..) = 1.0, [tunable = false] + x₂(..) = 0.0, [tunable = false] + x₃(..) = 0.0, [tunable = false] + u(..) = 0.5, [bounds = (0.0, 1.0), input = true] + end + + @constants begin + k₁ = 1.0, [tunable = false] + k₂ = 10.0, [tunable = false] + k₃ = 1.0, [tunable = false] + end + + eqs = [ + D(x₁(t)) ~ -u(t) * (k₁ * x₁(t) - k₂ * x₂(t)) + D(x₂(t)) ~ u(t) * (k₁ * x₁(t) - k₂ * x₂(t)) - (1 - u(t)) * k₃ * x₂(t) + D(x₃(t)) ~ (1 - u(t)) * k₃ * x₂(t) + ] + + # Define control discretization + constraint_grid = scaled_grids.constraint_grid + + grid_cons_u1 = [x₁(tᵢ) ≲ 1.1 for tᵢ in constraint_grid] + grid_cons_u2 = [x₂(tᵢ) ≲ 1.1 for tᵢ in constraint_grid] + grid_cons_u3 = [x₃(tᵢ) ≲ 1.1 for tᵢ in constraint_grid] + grid_cons_l1 = [x₁(tᵢ) ≳ 0.0 for tᵢ in constraint_grid] + grid_cons_l2 = [x₂(tᵢ) ≳ 0.0 for tᵢ in constraint_grid] + grid_cons_l3 = [x₃(tᵢ) ≳ 0.0 for tᵢ in constraint_grid] + + cons = [ + grid_cons_u1..., grid_cons_u2..., grid_cons_u3..., + grid_cons_l1..., grid_cons_l2..., grid_cons_l3..., + ] + + costs = [x₃(last(tspan))] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Jackson", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/lotka_competitive.jl b/lib/OptimalControlBenchmarks/src/problems/lotka_competitive.jl new file mode 100644 index 00000000..708a963d --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/lotka_competitive.jl @@ -0,0 +1,46 @@ +function lotka_competitive(grids) + + num_states = 3 + num_controls = 1 + tspan = (0.0, 40.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₀(..) = 0.5, [tunable = false] + x₁(..) = 1.5, [tunable = false] + u(..) = 0.5, [bounds = (0.0, 1.0), input = true] + end + + @constants begin + c₁ = 0.1, [tunable = false] + c₂ = 0.4, [tunable = false] + α = 1.2, [tunable = false] + K = 1.8, [tunable = false] + end + + eqs = [ + D(x₀(t)) ~ x₀(t) * (1 - (x₀(t) + α * x₁(t)) / K) - c₁ * x₀(t) * u(t) + D(x₁(t)) ~ x₁(t) * (1 - (x₀(t) + x₁(t)) / K) - c₂ * x₁(t) * u(t) + ] + + costs = [ + Symbolics.Integral(t in (0.0, last(tspan)))( + (x₀(t) - 1.0)^2 + (x₁(t) - 1.0)^2 + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Lotka Volterra Competitive", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/lotka_shared.jl b/lib/OptimalControlBenchmarks/src/problems/lotka_shared.jl new file mode 100644 index 00000000..9ec49a87 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/lotka_shared.jl @@ -0,0 +1,47 @@ +function lotka_shared(grids) + + num_states = 4 + num_controls = 1 + tspan = (0.0, 40.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₀(..) = 1.5, [tunable = false] + x₁(..) = 1.0, [tunable = false] + x₂(..) = 0.5, [tunable = false] + u(..) = 0.5, [bounds = (0.0, 1.0), input = true] + end + + @constants begin + c₁ = 0.1, [tunable = false] + c₂ = 0.4, [tunable = false] + α = 1.2, [tunable = false] + end + + eqs = [ + D(x₀(t)) ~ x₀(t) - x₀(t) * x₁(t) - x₀(t) * x₂(t) + D(x₁(t)) ~ -x₁(t) + x₀(t) * x₁(t) - c₁ * x₁(t) * u(t) + D(x₂(t)) ~ -x₂(t) + α * x₀(t) * x₂(t) - c₂ * x₂(t) * u(t) + ] + + costs = [ + Symbolics.Integral(t in (0.0, last(tspan)))( + (x₀(t) - 1.5)^2 + (x₁(t) - 1.0)^2 + (x₂(t) - 1.0)^2 + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Lotka Volterra Shared", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/lotka_volterra.jl b/lib/OptimalControlBenchmarks/src/problems/lotka_volterra.jl new file mode 100644 index 00000000..5a22beda --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/lotka_volterra.jl @@ -0,0 +1,44 @@ +function lotka_volterra(grids) + + num_states = 3 + num_controls = 1 + tspan = (0.0, 12.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₀(..) = 0.5, [tunable = false] + x₁(..) = 0.7, [tunable = false] + u(..) = 0.5, [bounds = (0.0, 1.0), input = true] + end + + @constants begin + c₀ = 0.4, [tunable = false] + c₁ = 0.2, [tunable = false] + end + + eqs = [ + D(x₀(t)) ~ x₀(t) - x₀(t) * x₁(t) - c₀ * x₀(t) * u(t) + D(x₁(t)) ~ -x₁(t) + x₀(t) * x₁(t) - c₁ * x₁(t) * u(t) + ] + + costs = [ + Symbolics.Integral(t in (0.0, last(tspan)))( + (x₀(t) - 1.0)^2 + (x₁(t) - 1.0)^2 + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Lotka Volterra fishing problem", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/lqr.jl b/lib/OptimalControlBenchmarks/src/problems/lqr.jl new file mode 100644 index 00000000..a32ae9d1 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/lqr.jl @@ -0,0 +1,42 @@ +function lqr(grids) + + num_states = 2 + num_controls = 1 + tspan = (0.0, 10.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x(..) = 1.0, [tunable = false] + u(..) = 0.5, [input = true] + end + + @constants begin + a = -1.0, [tunable = false] + b = 1.0, [tunable = false] + end + + eqs = [ + D(x(t)) ~ a * x(t) + b * u(t), + ] + + costs = [ + Symbolics.Integral(t in (0.0, last(tspan)))( + 10.0 * (x(t) - 3.0)^2 + 0.1 * u(t)^2 + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Linear Quadratic Regulator", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/moon_landing.jl b/lib/OptimalControlBenchmarks/src/problems/moon_landing.jl new file mode 100644 index 00000000..cbc2f80c --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/moon_landing.jl @@ -0,0 +1,47 @@ +function moon_landing(grids) + + num_states = 3 + num_controls = 1 + tspan = (0.0, 1.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + h(..) = 1.0, [tunable = false] + v(..) = -0.783, [tunable = false] + m(..) = 1.0, [tunable = false] + T(..) = 0.5, [bounds = (0.0, 1.227), input = true] + end + + @parameters begin + tₛ = 1.0, [bounds = (1.0e-3, Inf), tunable = true] + end + + eqs = [ + D(h(t)) ~ tₛ * v(t) + D(v(t)) ~ tₛ * (-1 + T(t) / m(t)) + D(m(t)) ~ tₛ * (- T(t) / 2.349) + ] + + cons = [ + h(last(tspan)) ~ 0.0, + v(last(tspan)) ~ 0.0, + ] + + costs = [-m(last(tspan))] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Moon Landing", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/mountain_car.jl b/lib/OptimalControlBenchmarks/src/problems/mountain_car.jl new file mode 100644 index 00000000..cff90166 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/mountain_car.jl @@ -0,0 +1,47 @@ +function mountain_car(grids) + + num_states = 3 + num_controls = 1 + tspan = (0.0, 100.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x(..) = -0.5, [tunable = false] + v(..) = 0.0, [tunable = false] + obj(..) = 0.0, [tunable = false] + u(..) = 0.5, [bounds = (-1.0, 1.0), input = true] + end + + @parameters begin + tₛ = 1.0, [bounds = (1.0e-3, Inf), tunable = true] + end + + eqs = [ + D(x(t)) ~ tₛ * v(t) + D(v(t)) ~ tₛ * (1.0e-3 * u(t) - 2.5e-3 * cos(3 * x(t))) + D(obj(t)) ~ tₛ + ] + + cons = [ + x(last(tspan)) ~ 0.5, + v(last(tspan)) ≳ 0.0, + ] + + costs = [obj(last(tspan))] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Mountain Car", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/ocean.jl b/lib/OptimalControlBenchmarks/src/problems/ocean.jl new file mode 100644 index 00000000..53be0771 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/ocean.jl @@ -0,0 +1,78 @@ +function ocean(grids) + + num_states = 3 + num_controls = 2 + tspan = (0.0, 400.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + S(..) = 2.0e3, [tunable = false] + R(..) = 1.0e4, [tunable = false] + u₁(..) = 7.0, [bounds = (0.0, 40.0), input = true] + u₂(..) = 7.0, [bounds = (0.0, 40.0), input = true] + end + + @constants begin + ρ = 3.0e-2, [tunable = false] + γ = 1.0e-3, [tunable = false] + ω = 0.1, [tunable = false] + b = 50.0, [tunable = false] + μ = 0.5, [tunable = false] + a₁ = 2.0, [tunable = false] + a₂ = 2.0, [tunable = false] + ν = 1.0, [tunable = false] + c₁ = 50.0, [tunable = false] + c₂ = 4.0e-3, [tunable = false] + Sₚ = 6.0e2, [tunable = false] + S₀ = 2.0e3, [tunable = false] + R₀ = 1.0e4, [tunable = false] + Dₗ₀ = 2.3e4, [tunable = false] + end + + # auxiliary functions + U = b * u₁(t) - μ * u₁(t)^2 + A = a₁ * u₂(t) + a₂ * u₂(t)^2 + C = c₁ - c₂ * R(t) + Df = ν * (0.3 * S(t) - Sₚ)^2 + Dₗ = Dₗ₀ + R₀ + S₀ - R(t) - S(t) + + eqs = [ + D(R(t)) ~ u₁(t) - u₂(t) - γ * (S(t) - ω * Dₗ) + D(S(t)) ~ -u₁(t) + ] + + # Define control discretization + constraint_grid = scaled_grids.constraint_grid + + grid_cons_Su = [S(tᵢ) ≲ 1.0e5 for tᵢ in constraint_grid] + grid_cons_Ru = [R(tᵢ) ≲ 1.0e5 for tᵢ in constraint_grid] + grid_cons_Sl = [S(tᵢ) ≳ 0.0 for tᵢ in constraint_grid] + grid_cons_Rl = [R(tᵢ) ≳ 0.0 for tᵢ in constraint_grid] + + cons = [ + grid_cons_Su..., grid_cons_Ru..., + grid_cons_Sl..., grid_cons_Rl..., + ] + + costs = [ + Symbolics.Integral(t in (0.0, last(tspan)))( + -exp(-ρ * t) * (U - A - u₁(t) * C - Df) + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Ocean", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/particle_steering.jl b/lib/OptimalControlBenchmarks/src/problems/particle_steering.jl new file mode 100644 index 00000000..9eeab754 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/particle_steering.jl @@ -0,0 +1,58 @@ +function particle_steering(grids) + + num_states = 5 + num_controls = 1 + tspan = (0.0, 10.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₁(..) = 0.0, [tunable = false] + x₂(..) = 0.0, [tunable = false] + dx₁(..) = 0.0, [tunable = false] + dx₂(..) = 0.0, [tunable = false] + obj(..) = 0.0, [tunable = false] + u(..) = 0.0, [bounds = (-pi / 2.0, pi / 2.0), input = true] + end + + @parameters begin + tₛ = 3.0e-2, [bounds = (1.0e-3, Inf), tunable = true] + end + + @constants begin + a = 100.0, [tunable = false] + end + + eqs = [ + D(x₁(t)) ~ tₛ * (dx₁(t)) + D(x₂(t)) ~ tₛ * (dx₂(t)) + D(dx₁(t)) ~ tₛ * (a * cos(u(t))) + D(dx₂(t)) ~ tₛ * (a * sin(u(t))) + D(obj(t)) ~ tₛ + ] + + cons = [ + x₂(last(tspan)) ~ 5.0, + dx₁(last(tspan)) ~ 45.0, + dx₂(last(tspan)) ~ 0.0, + ] + + costs = [ + obj(last(tspan)), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Particle steering problem", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/quadrotor.jl b/lib/OptimalControlBenchmarks/src/problems/quadrotor.jl new file mode 100644 index 00000000..e7b22e46 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/quadrotor.jl @@ -0,0 +1,68 @@ +function quadrotor(grids) + + num_states = 7 + num_controls = 4 + tspan = (0.0, 7.5) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₁(..) = 0.0, [tunable = false] + x₂(..) = 0.0, [tunable = false] + x₃(..) = 1.0, [tunable = false] + x₄(..) = 0.0, [tunable = false] + x₅(..) = 0.0, [tunable = false] + x₆(..) = 0.0, [tunable = false] + w₁(..) = 1.0 / 3.0, [bounds = (0.0, 1.0), input = true] + w₂(..) = 1.0 / 3.0, [bounds = (0.0, 1.0), input = true] + w₃(..) = 1.0 / 3.0, [bounds = (0.0, 1.0), input = true] + u(..) = 0.0, [bounds = (0.0, 1.0e-3), input = true] + end + + @constants begin + g = 9.81, [tunable = false] + M = 1.3, [tunable = false] + L = 0.305, [tunable = false] + I = 0.0605, [tunable = false] + end + + eqs = [ + D(x₁(t)) ~ x₂(t) + D(x₂(t)) ~ g * sin(x₅(t)) + w₁(t) * u(t) * sin(x₅(t)) / M + D(x₃(t)) ~ x₄(t) + D(x₄(t)) ~ g * cos(x₅(t)) - g + w₁(t) * u(t) * cos(x₅(t)) / M + D(x₅(t)) ~ x₆(t) + D(x₆(t)) ~ -w₂(t) * L * u(t) / I + w₃(t) * L * u(t) / I + ] + + constraint_grid = scaled_grids.constraint_grid + grid_cons_w = [w₁(tᵢ) + w₂(tᵢ) + w₃(tᵢ) ~ 1.0 for tᵢ in constraint_grid] + grid_cons_x₃ = [x₃(tᵢ) ≳ 0.0 for tᵢ in constraint_grid] + + cons = [ + grid_cons_w..., + grid_cons_x₃..., + ] + + t_f = last(tspan) + costs = [ + 5 * (x₁(t_f) - 6)^2 + 5 * (x₃(t_f) - 1)^2 + 5 * (sin(x₅(t_f) * 0.5))^2 + Symbolics.Integral(t in (0.0, last(tspan)))( + 5 * u(t)^2 + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Quadrotor helicopter control problem", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/rao_mease.jl b/lib/OptimalControlBenchmarks/src/problems/rao_mease.jl new file mode 100644 index 00000000..3837f3fa --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/rao_mease.jl @@ -0,0 +1,42 @@ +function rao_mease(grids) + + num_states = 2 + num_controls = 1 + tspan = (0.0, 10.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x(..) = 1.0, [tunable = false] + u(..) = 0.0, [input = true] + end + + eqs = [ + D(x(t)) ~ -x(t)^3 + u(t), + ] + + cons = [ + x(last(tspan)) ~ 1.5, + ] + + costs = [ + Symbolics.Integral(t in (0.0, last(tspan)))( + x(t)^2 + u(t)^2 + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Rao Mease", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/robbins.jl b/lib/OptimalControlBenchmarks/src/problems/robbins.jl new file mode 100644 index 00000000..34b8fe57 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/robbins.jl @@ -0,0 +1,60 @@ +function robbins(grids) + + num_states = 4 + num_controls = 1 + tspan = (0.0, 10.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₁(..) = 1.0, [tunable = false] + x₂(..) = -2.0, [tunable = false] + x₃(..) = 0.0, [tunable = false] + u(..) = 0.0, [input = true] + end + + @constants begin + α = 3.0, [tunable = false] + β = 0.0, [tunable = false] + γ = 0.5, [tunable = false] + end + + eqs = [ + D(x₁(t)) ~ x₂(t) + D(x₂(t)) ~ x₃(t) + D(x₃(t)) ~ u(t) + ] + + # scale the constraint grid + constraint_grid = scaled_grids.constraint_grid[1:(end - 1)] + + grid_cons_x₁ = [x₃(tᵢ) ≳ 0.0 for tᵢ in constraint_grid] + + cons = [ + x₁(last(tspan)) ~ 0.0, + x₂(last(tspan)) ~ 0.0, + x₃(last(tspan)) ~ 0.0, + grid_cons_x₁..., + ] + + costs = [ + Symbolics.Integral(t in (0.0, last(tspan)))( + α * x₁(t) + β * x₁(t)^2 + γ * u(t)^2 + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Robbins", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/three_tank.jl b/lib/OptimalControlBenchmarks/src/problems/three_tank.jl new file mode 100644 index 00000000..7fb8bfa9 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/three_tank.jl @@ -0,0 +1,70 @@ +function three_tank(grids) + + num_states = 4 + num_controls = 3 + tspan = (0.0, 12.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₁(..) = 2.0, [tunable = false] + x₂(..) = 2.0, [tunable = false] + x₃(..) = 2.0, [tunable = false] + w₁(..) = 1.0 / 3.0, [bounds = (0.0, 1.0), input = true] + w₂(..) = 1.0 / 3.0, [bounds = (0.0, 1.0), input = true] + w₃(..) = 1.0 / 3.0, [bounds = (0.0, 1.0), input = true] + end + + @constants begin + c₁ = 1.0, [tunable = false] + c₂ = 2.0, [tunable = false] + c₃ = 0.8, [tunable = false] + k₁ = 2.0, [tunable = false] + k₂ = 3.0, [tunable = false] + k₃ = 1.0, [tunable = false] + k₄ = 3.0, [tunable = false] + end + + eqs = [ + D(x₁(t)) ~ -sqrt(x₂(t)) + c₁ * w₁(t) + c₂ * w₂(t) - w₃(t) * sqrt(c₃ * x₁(t)) + D(x₂(t)) ~ sqrt(x₁(t)) - sqrt(x₂(t)) + D(x₃(t)) ~ sqrt(x₂(t)) - sqrt(x₃(t)) + w₃(t) * sqrt(c₃ * x₁(t)) + ] + + constraint_grid = scaled_grids.constraint_grid + + grid_cons_w = [w₁(tᵢ) + w₂(tᵢ) + w₃(tᵢ) ~ 1.0 for tᵢ in constraint_grid] + + # ensure that the state variables remain positive + grid_cons_x₁ = [x₁(tᵢ) ≳ 0.0 for tᵢ in constraint_grid] + grid_cons_x₂ = [x₂(tᵢ) ≳ 0.0 for tᵢ in constraint_grid] + grid_cons_x₃ = [x₃(tᵢ) ≳ 0.0 for tᵢ in constraint_grid] + + cons = [ + grid_cons_w..., + grid_cons_x₁..., + grid_cons_x₂..., + grid_cons_x₃..., + ] + + costs = [ + Symbolics.Integral(t in (0.0, last(tspan)))( + k₁ * (x₂(t) - k₂)^2 + k₃ * (x₃(t) - k₄)^2 + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Three Tank multimode problem", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/tubular_reactor.jl b/lib/OptimalControlBenchmarks/src/problems/tubular_reactor.jl new file mode 100644 index 00000000..9f38ccae --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/tubular_reactor.jl @@ -0,0 +1,35 @@ +function tubular_reactor(grids) + + num_states = 2 + num_controls = 1 + tspan = (0.0, 1.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₁(..) = 1.0, [tunable = false] + x₂(..) = 0.0, [tunable = false] + u(..) = 0.0, [bounds = (0.0, 5.0), input = true] + end + + eqs = [ + D(x₁(t)) ~ -(u(t) + 0.5 * u(t)^2) * x₁(t) + D(x₂(t)) ~ u(t) * x₁(t) + ] + + costs = [-x₂(last(tspan))] + + @named oc_problem = System( + eqs, + t; + costs = costs + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Tubular Reactor", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/problems/van_der_pol.jl b/lib/OptimalControlBenchmarks/src/problems/van_der_pol.jl new file mode 100644 index 00000000..62d1a298 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/problems/van_der_pol.jl @@ -0,0 +1,48 @@ +function van_der_pol(grids) + + num_states = 3 + num_controls = 1 + tspan = (0.0, 10.0) + + scaled_grids = scale_grids!(tspan, grids) + + @variables begin + x₁(..) = 0.0, [tunable = false] + x₂(..) = 1.0, [tunable = false] + u(..) = 0.0, [bounds = (-1.0, 1.0), input = true] + end + + eqs = [ + D(x₁(t)) ~ (1 - x₂(t)^2) * x₁(t) - x₂(t) + u(t) + D(x₂(t)) ~ x₁(t) + ] + + constraint_grid = scaled_grids.constraint_grid + + grid_cons_x₁ = [x₁(tᵢ) ≳ -0.25 for tᵢ in constraint_grid] + + cons = [ + grid_cons_x₁..., + ] + + costs = [ + Symbolics.Integral(t in (0.0, last(tspan)))( + x₁(t)^2 + x₂(t)^2 + u(t)^2 + ), + ] + + @named oc_problem = System( + eqs, + t; + costs = costs, + constraints = cons + ) + + return ( + system = oc_problem, + grids = scaled_grids, + dims = (num_states, num_controls), + name = "Van der Pol Oscillator", + ) + +end diff --git a/lib/OptimalControlBenchmarks/src/run_benchmarks.jl b/lib/OptimalControlBenchmarks/src/run_benchmarks.jl new file mode 100644 index 00000000..80da9b61 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/run_benchmarks.jl @@ -0,0 +1,42 @@ +using BenchmarkTools + +""" + run_all(benchmarks, optimizer, grids) + +Runs each benchmark problem with the Corleone benchmark solver. + +# Arguments +- `benchmarks`: Iterable of benchmark constructor functions, such as the result of + `load_benchmarks()`. +- `optimizer`: Optimization solver passed to `solve_with_corleone`. +- `grids`: Grid configuration passed to each benchmark constructor. + +# Returns +A vector of named tuples containing each benchmark name and measured runtime placeholder. +""" +function run_all(benchmarks, optimizer, grids) + + results = [] + + for prob in benchmarks + name = nameof(prob) + println("Running ", name) + sol = solve_with_corleone( + prob, + optimizer, + grids + ) + println(sol) + t = 0 # @belapsed solve_with_corleone($prob) + + push!( + results, ( + name = name, + time = t, + ) + ) + + end + + return results +end diff --git a/lib/OptimalControlBenchmarks/src/scale_grid.jl b/lib/OptimalControlBenchmarks/src/scale_grid.jl new file mode 100644 index 00000000..4980a6c4 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/scale_grid.jl @@ -0,0 +1,9 @@ +function scale_grids!(tspan::Tuple{Float64, Float64}, grids::BenchmarkGrids) + scale = last(tspan) - first(tspan) + shift = first(tspan) + + grids.control_grid .*= scale .+ shift + grids.shooting_grid .*= scale .+ shift + grids.constraint_grid .*= scale .+ shift + return grids +end diff --git a/lib/OptimalControlBenchmarks/src/solver_corleone.jl b/lib/OptimalControlBenchmarks/src/solver_corleone.jl new file mode 100644 index 00000000..e67b5986 --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/solver_corleone.jl @@ -0,0 +1,54 @@ +function solve_with_corleone(benchmark, optimizer, grids, plot_res = true) + + data = benchmark(grids) + + oc_problem = data.system + scaled_grids = data.grids + num_states, num_controls = data.dims + + # Extract control variable + controls = inputs(oc_problem) + control_grid = scaled_grids.control_grid + shooting_grid = scaled_grids.shooting_grid + + control_map = [ + c => control_grid for c in controls + ] + + dynopt = CorleoneDynamicOptProblem( + oc_problem, + [], + control_map...; + algorithm = Tsit5(), + shooting = shooting_grid + ) + + optprob = OptimizationProblem( + dynopt, + AutoForwardDiff(), + Val(:ComponentArrays) + ) + + sol = solve( + optprob, + optimizer, + max_iter = 100, + tol = 5.0e-6, + hessian_approximation = "limited-memory" + ) + + if plot_res + # plotting + u_opt = ComponentVector(sol.u, optprob.f.f.ax) + opt_traj, _ = dynopt.layer( + nothing, u_opt, + LuxCore.initialstates(Random.default_rng(), dynopt.layer) + ) + + f = plot_oc_problem(opt_traj, num_states, num_controls) + display(f) + end + + return sol + +end diff --git a/lib/OptimalControlBenchmarks/src/types.jl b/lib/OptimalControlBenchmarks/src/types.jl new file mode 100644 index 00000000..b1f562af --- /dev/null +++ b/lib/OptimalControlBenchmarks/src/types.jl @@ -0,0 +1,11 @@ +struct OptimalControlBenchmark + name::Symbol + description::String + make_problem::Function +end + +struct BenchmarkGrids + control_grid::Vector{Float64} + shooting_grid::Vector{Float64} + constraint_grid::Vector{Float64} +end diff --git a/lib/OptimalControlBenchmarks/test/qa/Project.toml b/lib/OptimalControlBenchmarks/test/qa/Project.toml new file mode 100644 index 00000000..1f38b948 --- /dev/null +++ b/lib/OptimalControlBenchmarks/test/qa/Project.toml @@ -0,0 +1,18 @@ +[deps] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +Corleone = "2751e0db-33e9-4374-b36d-b7219e1e6b40" +OptimalControlBenchmarks = "d147904a-1eb9-11f1-01e7-bdf88e0f9445" +SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[sources] +Corleone = {path = "../../../.."} +OptimalControlBenchmarks = {path = "../.."} + +[compat] +Aqua = "0.8" +Corleone = "0.0.5" +OptimalControlBenchmarks = "0.0.3" +SciMLTesting = "1.6" +Test = "1" +julia = "1.10" diff --git a/lib/OptimalControlBenchmarks/test/qa/qa.jl b/lib/OptimalControlBenchmarks/test/qa/qa.jl new file mode 100644 index 00000000..b5de8aa1 --- /dev/null +++ b/lib/OptimalControlBenchmarks/test/qa/qa.jl @@ -0,0 +1,17 @@ +using SciMLTesting +using OptimalControlBenchmarks +using Test + +run_qa( + OptimalControlBenchmarks; + explicit_imports = true, + ei_kwargs = (; + # `problem_registry.jl` registers problems through a dynamic `include`, + # which ExplicitImports cannot follow, so the module is unanalyzable for + # the implicit/stale checks. + no_implicit_imports = (; allow_unanalyzable = (OptimalControlBenchmarks,)), + no_stale_explicit_imports = (; allow_unanalyzable = (OptimalControlBenchmarks,)), + # `inputs` is re-exported by ModelingToolkit from ModelingToolkitBase. + all_explicit_imports_via_owners = (; ignore = (:inputs,)), + ), +) diff --git a/lib/OptimalControlBenchmarks/test/runtests.jl b/lib/OptimalControlBenchmarks/test/runtests.jl new file mode 100644 index 00000000..8e7becb4 --- /dev/null +++ b/lib/OptimalControlBenchmarks/test/runtests.jl @@ -0,0 +1,38 @@ +using OptimalControlBenchmarks +using Test +using SafeTestsets +using Pkg + +# QA (Aqua) runs in an isolated environment (test/qa) so its tooling deps never +# enter the main test target's resolve. On Julia < 1.11 the [sources] table is +# ignored, so develop the package and its in-repo sibling by path. +function activate_qa_env() + Pkg.activate(joinpath(@__DIR__, "qa")) + Pkg.develop( + [ + Pkg.PackageSpec(path = joinpath(@__DIR__, "..")), + Pkg.PackageSpec(path = joinpath(@__DIR__, "..", "..", "..")) + ] + ) + return Pkg.instantiate() +end + +# Centralized sublibrary CI sets CORLEONE_TEST_GROUP to the bare package name +# (-> "Core") or "_" (-> ""). Fall back to GROUP, then "All", so +# local `Pkg.test()` runs (which set neither) run everything. Map the value to +# the standard Core/QA section names this file dispatches on. +const _G = get(ENV, "CORLEONE_TEST_GROUP", get(ENV, "GROUP", "All")) +const _SUB = "OptimalControlBenchmarks" +const GROUP = _G == _SUB ? "Core" : (startswith(_G, _SUB * "_") ? _G[(length(_SUB) + 2):end] : _G) + +if GROUP == "All" || GROUP == "Core" + # No Core unit tests yet; running the benchmarks requires the full MTK + + # Ipopt/MOI optimal-control stack and is too expensive for CI. Keep a + # trivial assertion so the Core group still emits a Test Summary. + @test 1 == 1 +end + +if GROUP == "All" || GROUP == "QA" + activate_qa_env() + @safetestset "Code quality (Aqua.jl)" include("qa/qa.jl") +end diff --git a/lib/OptimalControlBenchmarks/test/test_groups.toml b/lib/OptimalControlBenchmarks/test/test_groups.toml new file mode 100644 index 00000000..1fe84cd4 --- /dev/null +++ b/lib/OptimalControlBenchmarks/test/test_groups.toml @@ -0,0 +1,5 @@ +[Core] +versions = ["lts", "1", "pre"] + +[QA] +versions = ["lts", "1"] diff --git a/src/local_controls.jl b/src/local_controls.jl index cea9a0c5..83fc5c12 100644 --- a/src/local_controls.jl +++ b/src/local_controls.jl @@ -56,22 +56,22 @@ get_timegrid(parameters::ControlParameter, tspan = (-Inf, Inf)) = begin t[idx] end -``` +""" $(SIGNATURES) Computes the number of discretized controls restricted to given `tspan`. -``` +""" function control_length(parameters::ControlParameter; tspan = nothing, kwargs...) (; t) = parameters idx = isnothing(tspan) ? eachindex(t) : findall(tspan[1] .<= t .< tspan[2]) return size(idx, 1) end -``` +""" $(SIGNATURES) Returns discretized controls of ControlParameter `params` restricted to given `tspan`. -``` +""" function get_controls(::Random.AbstractRNG, parameters::ControlParameter{<:Any, <:AbstractArray}; raw = false, tspan = nothing, kwargs...) (; t, controls) = parameters raw && return controls @@ -86,11 +86,11 @@ function get_controls(rng::Random.AbstractRNG, parameters::ControlParameter{<:An return parameters.controls(rng, t[idx], bounds) end -``` +""" $(SIGNATURES) Returns bounds of discretized controls restricted to given `tspan`. -``` +""" function get_bounds(parameters::ControlParameter{<:Any, <:Any, <:Tuple}; tspan = nothing, kwargs...) (; t) = parameters idx = isnothing(tspan) ? eachindex(t) : findall(tspan[1] .<= t .< tspan[2]) @@ -194,10 +194,10 @@ function collect_tspans(controls::ControlParameter...; tspan = (-Inf, Inf), subd return tuple(fullgrid...) end -``` +""" $(SIGNATURES) Collect all discretized `controls` in a flat vector. -``` +""" function collect_local_controls(rng, controls::ControlParameter...; kwargs...) return reduce( vcat, map(controls) do control @@ -206,10 +206,10 @@ function collect_local_controls(rng, controls::ControlParameter...; kwargs...) ) end -``` +""" $(SIGNATURES) Collect all lower and upper bounds of discretized `controls` in flat vectors. -``` +""" function collect_local_control_bounds(controls::ControlParameter...; kwargs...) bounds = map(controls) do control get_bounds(control; kwargs...) diff --git a/src/multiple_shooting.jl b/src/multiple_shooting.jl index b3d9e323..da50f1db 100644 --- a/src/multiple_shooting.jl +++ b/src/multiple_shooting.jl @@ -54,19 +54,19 @@ function default_initialization(rng::Random.AbstractRNG, shooting::MultipleShoot return NamedTuple{names}(vals) end -``` +""" $(METHODLIST) Constructs a `MultipleShootingLayer` from given `AbstractDEProblem` `prob` with suitable integration method `alg` and vector of shooting points `tpoints`. -``` +""" function MultipleShootingLayer(prob, alg, tpoints::AbstractVector; kwargs...) return MultipleShootingLayer(prob, alg, tpoints...; kwargs...) end function MultipleShootingLayer( prob::SciMLBase.AbstractDEProblem, - alg::SciMLBase.DEAlgorithm, + alg::SciMLBase.AbstractDEAlgorithm, tpoints::Real...; ensemble_alg = EnsembleSerial(), initialization = default_initialization, @@ -296,6 +296,19 @@ end shooting_constraints!(res, traj::Trajectory) = res +shooting_constraints(trajs::AbstractVector{<:Trajectory}) = reduce(vcat, reduce(vcat, shooting_violations.(trajs))) + +function shooting_constraints!(res::AbstractVector, trajs::AbstractVector{<:Trajectory}) + i = 0 + for traj in trajs + for subvec in traj.shooting, j in eachindex(subvec) + i += 1 + res[i] = subvec[j] + end + end + return res +end + """ $(SIGNATURES) @@ -314,10 +327,10 @@ function get_block_structure(mslayer::MultipleShootingLayer) return vcat(0, cumsum(ps_lengths)) end -``` +""" $(SIGNATURES) Extracts lower and upper bounds of all optimization variables in the `MultipleShootingLayer`. -``` +""" function get_bounds(mslayer::MultipleShootingLayer) (; layer, shooting_intervals) = mslayer names = ntuple(i -> Symbol(:interval, "_", i), length(shooting_intervals)) diff --git a/src/single_shooting.jl b/src/single_shooting.jl index 5a3c91d1..db140b22 100644 --- a/src/single_shooting.jl +++ b/src/single_shooting.jl @@ -214,7 +214,7 @@ function __parameterlength( return N end -function retrieve_symbol_cache(problem::SciMLBase.DEProblem, control_indices; control_names = [Symbol(:u, _subscript(u_id)) for u_id in 1:length(control_indices)]) +function retrieve_symbol_cache(problem::SciMLBase.AbstractDEProblem, control_indices; control_names = [Symbol(:u, _subscript(u_id)) for u_id in 1:length(control_indices)]) return retrieve_symbol_cache(problem.f.sys, problem.u0, problem.p, control_indices; control_names = control_names) end diff --git a/test/local_controls.jl b/test/core/local_controls.jl similarity index 96% rename from test/local_controls.jl rename to test/core/local_controls.jl index bbb27bcf..764be3d4 100644 --- a/test/local_controls.jl +++ b/test/core/local_controls.jl @@ -1,10 +1,10 @@ using Corleone using OrdinaryDiffEqTsit5 using Test -using Random +using StableRNGs using LuxCore -rng = Random.default_rng() +rng = StableRNG(42) c = ControlParameter(0:0.01:1.0) lb, ub = Corleone.get_bounds(c) @@ -70,7 +70,7 @@ c5 = ControlParameter(1.0:10.0, bounds = (ones(10), -ones(10)), controls = colle ) layer = Corleone.SingleShootingLayer(prob, Tsit5(), controls = ([2, 3, 1] .=> [c2, c3, c1])) - ps, st = LuxCore.setup(Random.default_rng(), layer) + ps, st = LuxCore.setup(rng, layer) sol, _ = layer(nothing, ps, st) @test all(0.6 .<= sol[:con3] .<= 0.8) diff --git a/test/multiple_shooting.jl b/test/core/multiple_shooting.jl similarity index 88% rename from test/multiple_shooting.jl rename to test/core/multiple_shooting.jl index 78173870..9d1c990a 100644 --- a/test/multiple_shooting.jl +++ b/test/core/multiple_shooting.jl @@ -1,12 +1,12 @@ using Corleone using OrdinaryDiffEqTsit5 using Test -using Random +using StableRNGs using LuxCore using ComponentArrays using LinearAlgebra -rng = Random.default_rng() +rng = StableRNG(42) function lotka_dynamics(u, p, t) return [ @@ -51,62 +51,68 @@ lb, ub = Corleone.get_bounds(layer) @test Corleone.get_block_structure(layer) == vcat(0, blocks) traj, st2 = layer(nothing, ps, st) @test Corleone.is_shooting_solution(traj) - @test Corleone.shooting_constraints(traj) == [ - 1.375754960969482, - 0.22357357513551113, - 1.2417260108009396, - 1.3757549609694817, - 0.22357357513551046, - 1.2417260108009387, - 1.3757549609694821, - 0.2235735751355138, - 1.2417260108009425, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - ] + @test isapprox( + Corleone.shooting_constraints(traj), [ + 1.375754960969482, + 0.22357357513551113, + 1.2417260108009396, + 1.3757549609694817, + 0.22357357513551046, + 1.2417260108009387, + 1.3757549609694821, + 0.2235735751355138, + 1.2417260108009425, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + ], atol = 1.0e-10 + ) @test Corleone.get_number_of_shooting_constraints(layer) == length(Corleone.shooting_constraints(traj)) == 15 res = zeros(30) - @test Corleone.shooting_constraints!(res[10:24], traj) == [ - 1.375754960969482, - 0.22357357513551113, - 1.2417260108009396, - 1.3757549609694817, - 0.22357357513551046, - 1.2417260108009387, - 1.3757549609694821, - 0.2235735751355138, - 1.2417260108009425, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - ] + @test isapprox( + Corleone.shooting_constraints!(res[10:24], traj), [ + 1.375754960969482, + 0.22357357513551113, + 1.2417260108009396, + 1.3757549609694817, + 0.22357357513551046, + 1.2417260108009387, + 1.3757549609694821, + 0.2235735751355138, + 1.2417260108009425, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + ], atol = 1.0e-10 + ) @views Corleone.shooting_constraints!(res[10:24], traj) - @test res[10:24] == [ - 1.375754960969482, - 0.22357357513551113, - 1.2417260108009396, - 1.3757549609694817, - 0.22357357513551046, - 1.2417260108009387, - 1.3757549609694821, - 0.2235735751355138, - 1.2417260108009425, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - 0.0, - ] + @test isapprox( + res[10:24], [ + 1.375754960969482, + 0.22357357513551113, + 1.2417260108009396, + 1.3757549609694817, + 0.22357357513551046, + 1.2417260108009387, + 1.3757549609694821, + 0.2235735751355138, + 1.2417260108009425, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + 0.0, + ], atol = 1.0e-10 + ) end @testset "Parallel" begin @@ -211,8 +217,7 @@ end traj, _ = layer(nothing, ps, st) @test isempty(ps.interval_1.u0) @test ps.interval_2.u0 == [0.5, 3.0, -5.0] - @test ps.interval_3.u0 == - [0.2875960819428889, 0.27699894224421623, -1.088467272254529] + @test isapprox(ps.interval_3.u0, [0.2875960819428889, 0.27699894224421623, -1.088467272254529], atol = 1.0e-10) @test ps.interval_4.u0 == u0 end end diff --git a/test/examples/.PLACEHOLDER b/test/examples/.PLACEHOLDER deleted file mode 100644 index e69de29b..00000000 diff --git a/test/examples/Project.toml b/test/examples/Project.toml new file mode 100644 index 00000000..f99bc33b --- /dev/null +++ b/test/examples/Project.toml @@ -0,0 +1,41 @@ +[deps] +ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66" +Corleone = "2751e0db-33e9-4374-b36d-b7219e1e6b40" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9" +LuxCore = "bb33d45b-7691-41d6-9220-0943567d0623" +ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78" +Optimization = "7f7a1694-90dd-40f0-9382-eb1efda571ba" +OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1" +OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8" +OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" +SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1" +SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283" +SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[compat] +ComponentArrays = "0.15" +Corleone = "0.0.5" +ForwardDiff = "1.2" +Ipopt = "1.10" +LuxCore = "1.4" +ModelingToolkit = "11" +Optimization = "5" +OptimizationMOI = "1" +OrdinaryDiffEqBDF = "1.14.0, 2" +OrdinaryDiffEqTsit5 = "1, 2" +Random = "1" +SafeTestsets = "0.1, 1" +SciMLBase = "2.141.0, 3" +SciMLSensitivity = "7.87" +SciMLTesting = "1" +SymbolicIndexingInterface = "0.3.43" +Test = "1" +julia = "1.10" + +[sources] +Corleone = {path = "../.."} diff --git a/test/examples/lotka_ms.jl b/test/examples/lotka_ms.jl index 35709c0d..27401153 100644 --- a/test/examples/lotka_ms.jl +++ b/test/examples/lotka_ms.jl @@ -1,12 +1,12 @@ using Corleone using OrdinaryDiffEqTsit5 using Test -using Random +using StableRNGs using LuxCore using ComponentArrays using Optimization, OptimizationMOI, Ipopt -rng = Random.default_rng() +rng = StableRNG(42) function lotka_dynamics!(du, u, p, t) du[1] = u[1] - p[2] * prod(u[1:2]) - 0.4 * p[1] * u[1] diff --git a/test/examples/lotka_oc.jl b/test/examples/lotka_oc.jl index ff8246a1..77093886 100644 --- a/test/examples/lotka_oc.jl +++ b/test/examples/lotka_oc.jl @@ -1,7 +1,7 @@ using Corleone using OrdinaryDiffEqTsit5 using Test -using Random +using StableRNGs using LuxCore using ComponentArrays using Optimization, OptimizationMOI, Ipopt @@ -11,7 +11,7 @@ using SciMLSensitivity.ReverseDiff using SciMLSensitivity.Zygote using SymbolicIndexingInterface -rng = Random.default_rng() +rng = StableRNG(42) function lotka_dynamics!(du, u, p, t) du[1] = u[1] - p[2] * prod(u[1:2]) - 0.4 * p[1] * u[1] diff --git a/test/examples/mtk.jl b/test/examples/mtk.jl index f17f0c9f..8922b906 100644 --- a/test/examples/mtk.jl +++ b/test/examples/mtk.jl @@ -6,7 +6,9 @@ using OrdinaryDiffEqTsit5 using ComponentArrays, ForwardDiff using Optimization using OptimizationMOI, Ipopt -using LuxCore, Random +using LuxCore, StableRNGs + +rng = StableRNG(42) @variables x(..) = 0.5 [tunable = false] y(..) = 0.7 [tunable = false] @@ -49,7 +51,7 @@ cons = [ @test size(optprob.lcons, 1) == size(optprob.ucons, 1) == length(cons) - ps, st = LuxCore.setup(Random.default_rng(), dynopt.layer) + ps, st = LuxCore.setup(rng, dynopt.layer) traj, _ = dynopt.layer(nothing, ps, st) @@ -82,7 +84,7 @@ end @test size(optprob.lcons, 1) == size(optprob.ucons, 1) == length(cons) + Corleone.get_number_of_shooting_constraints(dynopt.layer) - ps, st = LuxCore.setup(Random.default_rng(), dynopt.layer) + ps, st = LuxCore.setup(rng, dynopt.layer) traj, _ = dynopt.layer(nothing, ps, st) diff --git a/test/qa/Project.toml b/test/qa/Project.toml new file mode 100644 index 00000000..6a5a1841 --- /dev/null +++ b/test/qa/Project.toml @@ -0,0 +1,17 @@ +[deps] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +Corleone = "2751e0db-33e9-4374-b36d-b7219e1e6b40" +SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" + +[sources] +Corleone = {path = "../.."} + +[compat] +Aqua = "0.8" +Corleone = "0.0.5" +SafeTestsets = "0.1, 1" +SciMLTesting = "1.6" +Test = "1" +julia = "1.10" diff --git a/test/qa/qa.jl b/test/qa/qa.jl new file mode 100644 index 00000000..57ffabf7 --- /dev/null +++ b/test/qa/qa.jl @@ -0,0 +1,29 @@ +using SciMLTesting +using Corleone +using Test + +run_qa( + Corleone; + explicit_imports = true, + # Corleone pulls all of its deps in with bare `using`, so the package leans + # on a large set of implicit imports. Converting every one to an explicit + # import is a sizable refactor tracked in SciML/Corleone.jl#103. + ei_broken = (:no_implicit_imports,), + ei_kwargs = (; + # `ADTypes` is reached through `SciMLBase.ADTypes.AbstractADType`; ADTypes + # is not a direct Corleone dependency, so the access goes via SciMLBase. + all_qualified_accesses_via_owners = (; ignore = (:ADTypes,)), + # Names still not declared public in their owning modules: Base internals + # (`AbstractVecOrTuple`, `Splat`), SciMLBase internals (`ADTypes`, + # `AbstractDEAlgorithm`, `AbstractDEProblem`, `EnsembleAlgorithm`, + # `get_colorizers`), and SciMLStructures internals (`Tunable`, + # `canonicalize`). + all_qualified_accesses_are_public = (; + ignore = ( + :ADTypes, :AbstractDEAlgorithm, :AbstractDEProblem, + :AbstractVecOrTuple, :EnsembleAlgorithm, :Splat, :Tunable, + :canonicalize, :get_colorizers, + ), + ), + ), +) diff --git a/test/runtests.jl b/test/runtests.jl index d4c99249..b94e4527 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,34 +1,99 @@ -using Test -using JET -using SafeTestsets +using Pkg +using SafeTestsets, Test +using SciMLTesting -@testset "Corleone.jl" begin - @safetestset "Code quality (Aqua.jl)" begin - using Aqua - using Corleone - Aqua.test_all(Corleone) - end - @safetestset "Local controls" begin - include("local_controls.jl") - end - @safetestset "Multiple shooting" begin - include("multiple_shooting.jl") - end - @testset "Examples" begin - @safetestset "Lotka" begin - include("examples/lotka_oc.jl") - end - @safetestset "Lotka MS" begin - include("examples/lotka_ms.jl") - end - @safetestset "Lotka MTK" begin - include("examples/mtk.jl") +const GROUP = current_group() +const LIB_DIR = joinpath(dirname(@__DIR__), "lib") + +# Centralized sublibrary CI (SciML/.github sublibrary-project-tests.yml@v1) tests +# each lib/ via the project model and never routes through this file. This +# dispatcher only matters when the root suite is invoked with a GROUP that names +# a sublibrary (e.g. local `GROUP=CorleoneOED julia test/runtests.jl`): the bare +# sublibrary name selects that sublibrary's "Core" group and "_" +# selects a named group. We then activate the sublibrary's own test environment +# and hand off to its runtests.jl via CORLEONE_TEST_GROUP. The sublibrary Pkg.test +# is done explicitly here (rather than via run_tests's built-in lib_dir path) so the +# Julia < 1.11 transitive [sources] develop walk is preserved verbatim. Otherwise +# the main-package suite below runs via run_tests: GROUP="All"/"Core" run the light +# Core group; "Examples" and "QA" each activate their own test/ environment. +base_group, test_group = detect_sublibrary_group(GROUP, LIB_DIR) + +if !isempty(base_group) && isdir(joinpath(LIB_DIR, base_group)) + sublib_path = joinpath(LIB_DIR, base_group) + Pkg.activate(sublib_path) + # On Julia < 1.11 the [sources] table in Project.toml is ignored, so develop + # the local path dependencies (e.g. Corleone) to test the PR branch code. + # Walk [sources] transitively in case a developed dependency carries its own. + if VERSION < v"1.11.0-DEV.0" + developed = Set{String}() + push!(developed, normpath(sublib_path)) + specs = Pkg.PackageSpec[] + queue = [sublib_path] + while !isempty(queue) + pkg_dir = popfirst!(queue) + toml_path = joinpath(pkg_dir, "Project.toml") + isfile(toml_path) || continue + toml = Pkg.TOML.parsefile(toml_path) + if haskey(toml, "sources") + for (dep_name, source_spec) in toml["sources"] + if source_spec isa Dict && haskey(source_spec, "path") + dep_path = normpath(joinpath(pkg_dir, source_spec["path"])) + if isdir(dep_path) && !(dep_path in developed) + push!(developed, dep_path) + push!(specs, Pkg.PackageSpec(path = dep_path)) + push!(queue, dep_path) + end + end + end + end end + isempty(specs) || Pkg.develop(specs) end + # Hand the resolved test group to the sublibrary runtests.jl, which reads + # CORLEONE_TEST_GROUP (matching the SublibraryCI group-env-name). + withenv("CORLEONE_TEST_GROUP" => test_group) do + Pkg.test(base_group; allow_reresolve = true) + end +else + run_tests(; + # Core: light tests that resolve in the main test target's environment. + core = function () + @safetestset "Local controls" begin + include("core/local_controls.jl") + end + return @safetestset "Multiple shooting" begin + include("core/multiple_shooting.jl") + end + end, + groups = Dict( + # Examples: the optimal-control example scripts pull in the heavy + # MTK + Ipopt/MOI + SciMLSensitivity stack, isolated in test/examples + # so those deps never enter the light main test target's resolve. + "Examples" => (; + env = joinpath(@__DIR__, "examples"), + body = function () + @safetestset "Lotka" begin + include("examples/lotka_oc.jl") + end + @safetestset "Lotka MS" begin + include("examples/lotka_ms.jl") + end + return @safetestset "Lotka MTK" begin + include("examples/mtk.jl") + end + end, + ), + ), + qa = (; + env = joinpath(@__DIR__, "qa"), + body = function () + return @safetestset "Code quality (Aqua.jl)" include("qa/qa.jl") + end, + ), + # The original ran only the Core group for the default GROUP="All"; the + # dep-adding Examples and QA groups ran only when explicitly selected. + all = ["Core"], + sublib_env = "CORLEONE_TEST_GROUP", + lib_dir = LIB_DIR, + ) end - -# What to test? -# local_controls.jl: -# - construction of index_grid, get_subvector_indices -> Julius -# general: -# - more convergence? Lotka OED diff --git a/test/test_groups.toml b/test/test_groups.toml new file mode 100644 index 00000000..d2dfda16 --- /dev/null +++ b/test/test_groups.toml @@ -0,0 +1,8 @@ +[Core] +versions = ["lts", "1", "pre"] + +[Examples] +versions = ["lts", "1", "pre"] + +[QA] +versions = ["lts", "1"]