diff --git a/Manifest-v1.11.toml b/Manifest-v1.11.toml
index c46938fd2..36101ba76 100644
--- a/Manifest-v1.11.toml
+++ b/Manifest-v1.11.toml
@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.11.1"
+julia_version = "1.11.0"
 manifest_format = "2.0"
-project_hash = "f5f80681e1f41a61ae2064c46066ad6e608da783"
+project_hash = "6ea7e541a27709a6f9e32ca4b0816717a46eb307"
 
 [[deps.ARFFFiles]]
 deps = ["CategoricalArrays", "Dates", "Parsers", "Tables"]
@@ -10,6 +10,12 @@ git-tree-sha1 = "678eb18590a8bc6674363da4d5faa4ac09c40a18"
 uuid = "da404889-ca92-49ff-9e8b-0aa6b4d38dc8"
 version = "1.5.0"
 
+[[deps.ASL_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "6252039f98492252f9e47c312c8ffda0e3b9e78d"
+uuid = "ae81ac8f-d209-56e5-92de-9978fef736f9"
+version = "0.1.3+0"
+
 [[deps.AWS]]
 deps = ["Base64", "Compat", "Dates", "Downloads", "GitHub", "HTTP", "IniFile", "JSON", "MbedTLS", "Mocking", "OrderedCollections", "Random", "SHA", "Sockets", "URIs", "UUIDs", "XMLDict"]
 git-tree-sha1 = "319ade7f8fc88243369e119859a7d3a3e7e7f267"
@@ -71,9 +77,9 @@ version = "0.1.38"
 
 [[deps.Adapt]]
 deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "6a55b747d1812e699320963ffde36f1ebdda4099"
+git-tree-sha1 = "d80af0733c99ea80575f612813fa6aa71022d33a"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "4.0.4"
+version = "4.1.0"
 weakdeps = ["StaticArrays"]
 
     [deps.Adapt.extensions]
@@ -207,6 +213,12 @@ git-tree-sha1 = "aebf55e6d7795e02ca500a689d326ac979aaf89e"
 uuid = "9718e550-a3fa-408a-8086-8db961cd8217"
 version = "0.1.1"
 
+[[deps.BenchmarkTools]]
+deps = ["JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
+git-tree-sha1 = "f1dff6729bc61f4d49e140da1af55dcd1ac97b2f"
+uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+version = "1.5.0"
+
 [[deps.BitFlags]]
 git-tree-sha1 = "0691e34b3bb8be9307330f88d1a3c3f25466c24d"
 uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35"
@@ -310,6 +322,12 @@ git-tree-sha1 = "9ebb045901e9bbf58767a9f34ff89831ed711aae"
 uuid = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
 version = "0.15.7"
 
+[[deps.CodecBzip2]]
+deps = ["Bzip2_jll", "TranscodingStreams"]
+git-tree-sha1 = "e7c529cc31bb85b97631b922fa2e6baf246f5905"
+uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
+version = "0.8.4"
+
 [[deps.CodecInflate64]]
 deps = ["TranscodingStreams"]
 git-tree-sha1 = "d981a6e8656b1e363a2731716f46851a2257deb7"
@@ -445,6 +463,12 @@ git-tree-sha1 = "abe83f3a2f1b857aac70ef8b269080af17764bbe"
 uuid = "9a962f9c-6df0-11e9-0e5d-c546b8b5ee8a"
 version = "1.16.0"
 
+[[deps.DataEnvelopmentAnalysis]]
+deps = ["Distributed", "Distributions", "GLPK", "InvertedIndices", "Ipopt", "JuMP", "LinearAlgebra", "PrecompileTools", "Printf", "ProgressMeter", "Random", "SharedArrays", "SparseArrays", "Statistics", "StatsAPI", "StatsBase"]
+git-tree-sha1 = "d79c7c880808ca66f24b65d8e064b8eca43c183a"
+uuid = "a100299e-89d6-11e9-0fa0-2daf497e6a05"
+version = "0.9.1"
+
 [[deps.DataFrames]]
 deps = ["Compat", "DataAPI", "DataStructures", "Future", "InlineStrings", "InvertedIndices", "IteratorInterfaceExtensions", "LinearAlgebra", "Markdown", "Missings", "PooledArrays", "PrecompileTools", "PrettyTables", "Printf", "Random", "Reexport", "SentinelArrays", "SortingAlgorithms", "Statistics", "TableTraits", "Tables", "Unicode"]
 git-tree-sha1 = "fb61b4812c49343d7ef0b533ba982c46021938a6"
@@ -512,15 +536,15 @@ version = "0.1.10"
 
 [[deps.DiskArrays]]
 deps = ["LRUCache", "OffsetArrays"]
-git-tree-sha1 = "dea23e3d787c25fbf5f879c4501da828bc729472"
+git-tree-sha1 = "e0e89a60637a62d13aa2107f0acd169b9b9b77e7"
 uuid = "3c3547ce-8d99-4f5e-a174-61eb10b00ae3"
-version = "0.4.5"
+version = "0.4.6"
 
 [[deps.Distances]]
 deps = ["LinearAlgebra", "Statistics", "StatsAPI"]
-git-tree-sha1 = "66c4c81f259586e8f002eacebc177e1fb06363b0"
+git-tree-sha1 = "c7e3a542b999843086e2f29dac96a618c105be1d"
 uuid = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
-version = "0.10.11"
+version = "0.10.12"
 weakdeps = ["ChainRulesCore", "SparseArrays"]
 
     [deps.Distances.extensions]
@@ -742,6 +766,18 @@ git-tree-sha1 = "273bd1cd30768a2fddfa3fd63bbc746ed7249e5f"
 uuid = "38e38edf-8417-5370-95a0-9cbb8c7f171a"
 version = "1.9.0"
 
+[[deps.GLPK]]
+deps = ["GLPK_jll", "MathOptInterface"]
+git-tree-sha1 = "1d706bd23e5d2d407bfd369499ee6f96afb0c3ad"
+uuid = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
+version = "1.2.1"
+
+[[deps.GLPK_jll]]
+deps = ["Artifacts", "GMP_jll", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "fe68622f32828aa92275895fdb324a85894a5b1b"
+uuid = "e8aa6df9-e6ca-548a-97ff-1f85fc5b8b98"
+version = "5.0.1+0"
+
 [[deps.GMP_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "781609d7-10c4-51f6-84f2-b8444358ff6d"
@@ -1007,6 +1043,18 @@ git-tree-sha1 = "0dc7b50b8d436461be01300fd8cd45aa0274b038"
 uuid = "41ab1584-1d38-5bbf-9106-f11c6c58b48f"
 version = "1.3.0"
 
+[[deps.Ipopt]]
+deps = ["Ipopt_jll", "LinearAlgebra", "MathOptInterface", "OpenBLAS32_jll", "PrecompileTools"]
+git-tree-sha1 = "76315a7100e9e901d959c25b6ab6d4a37e799132"
+uuid = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
+version = "1.6.7"
+
+[[deps.Ipopt_jll]]
+deps = ["ASL_jll", "Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "MUMPS_seq_jll", "SPRAL_jll", "libblastrampoline_jll"]
+git-tree-sha1 = "a0950d209a055b3adb6d29ade5cbdf005a6bd290"
+uuid = "9cc047cb-c261-5740-88fc-0cf96f7bdcc7"
+version = "300.1400.1600+0"
+
 [[deps.IrrationalConstants]]
 git-tree-sha1 = "630b497eafcc20001bba38a4651b327dcfc491d2"
 uuid = "92d709cd-6900-40b7-9082-c6be49f344b6"
@@ -1064,6 +1112,16 @@ git-tree-sha1 = "25ee0be4d43d0269027024d75a24c24d6c6e590c"
 uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
 version = "3.0.4+0"
 
+[[deps.JuMP]]
+deps = ["LinearAlgebra", "MacroTools", "MathOptInterface", "MutableArithmetics", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays"]
+git-tree-sha1 = "c91f872c6150cf1471f9cb279f5e0dc09423bdcf"
+uuid = "4076af6c-e467-56ae-b986-b466b2749572"
+version = "1.23.3"
+weakdeps = ["DimensionalData"]
+
+    [deps.JuMP.extensions]
+    JuMPDimensionalDataExt = "DimensionalData"
+
 [[deps.JuliaVariables]]
 deps = ["MLStyle", "NameResolution"]
 git-tree-sha1 = "49fb3cb53362ddadb4415e9b73926d6b40709e70"
@@ -1083,9 +1141,9 @@ version = "1.19.3+0"
 
 [[deps.KernelAbstractions]]
 deps = ["Adapt", "Atomix", "InteractiveUtils", "MacroTools", "PrecompileTools", "Requires", "StaticArrays", "UUIDs", "UnsafeAtomics", "UnsafeAtomicsLLVM"]
-git-tree-sha1 = "04e52f596d0871fa3890170fa79cb15e481e4cd8"
+git-tree-sha1 = "e73a077abc7fe798fe940deabe30ef6c66bdde52"
 uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
-version = "0.9.28"
+version = "0.9.29"
 
     [deps.KernelAbstractions.extensions]
     EnzymeExt = "EnzymeCore"
@@ -1104,10 +1162,10 @@ uuid = "88015f11-f218-50d7-93a8-a6af411a945d"
 version = "4.0.0+0"
 
 [[deps.LLVM]]
-deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Requires", "Unicode"]
-git-tree-sha1 = "4ad43cb0a4bb5e5b1506e1d1f48646d7e0c80363"
+deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Unicode"]
+git-tree-sha1 = "d422dfd9707bec6617335dc2ea3c5172a87d5908"
 uuid = "929cbde3-209d-540e-8aea-75f648917ca0"
-version = "9.1.2"
+version = "9.1.3"
 
     [deps.LLVM.extensions]
     BFloat16sExt = "BFloat16s"
@@ -1206,9 +1264,9 @@ version = "1.11.0"
 
 [[deps.Libgcrypt_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libgpg_error_jll"]
-git-tree-sha1 = "9fd170c4bbfd8b935fdc5f8b7aa33532c991a673"
+git-tree-sha1 = "8be878062e0ffa2c3f67bb58a595375eda5de80b"
 uuid = "d4300ac3-e22c-5743-9152-c294e39db1e4"
-version = "1.8.11+0"
+version = "1.11.0+0"
 
 [[deps.Libglvnd_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg", "Xorg_libX11_jll", "Xorg_libXext_jll"]
@@ -1218,15 +1276,15 @@ version = "1.6.0+0"
 
 [[deps.Libgpg_error_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "fbb1f2bef882392312feb1ede3615ddc1e9b99ed"
+git-tree-sha1 = "c6ce1e19f3aec9b59186bdf06cdf3c4fc5f5f3e6"
 uuid = "7add5ba3-2f88-524e-9cd5-f83b8a55f7b8"
-version = "1.49.0+0"
+version = "1.50.0+0"
 
 [[deps.Libiconv_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "f9557a255370125b405568f9767d6d195822a175"
+git-tree-sha1 = "61dfdba58e585066d8bce214c5a51eaa0539f269"
 uuid = "94ce4f54-9a6c-5748-9c1c-f9c7231a4531"
-version = "1.17.0+0"
+version = "1.17.0+1"
 
 [[deps.Libtiff_jll]]
 deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"]
@@ -1285,9 +1343,9 @@ version = "1.11.0"
 
 [[deps.LoggingExtras]]
 deps = ["Dates", "Logging"]
-git-tree-sha1 = "c1dd6d7978c12545b4179fb6153b9250c96b0075"
+git-tree-sha1 = "f02b56007b064fbfddb4c9cd60161b6dd0f40df3"
 uuid = "e6f89c97-d47a-5376-807f-9c37f3926c36"
-version = "1.0.3"
+version = "1.1.0"
 
 [[deps.Lz4_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1301,6 +1359,12 @@ git-tree-sha1 = "1d2dd9b186742b0f317f2530ddcbf00eebb18e96"
 uuid = "23992714-dd62-5051-b70f-ba57cb901cac"
 version = "0.10.7"
 
+[[deps.METIS_jll]]
+deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
+git-tree-sha1 = "1fd0a97409e418b78c53fac671cf4622efdf0f21"
+uuid = "d00139f3-1899-568f-a2f0-47f597d42d70"
+version = "5.1.2+0"
+
 [[deps.MLFlowClient]]
 deps = ["Dates", "FilePathsBase", "HTTP", "JSON", "ShowCases", "URIs", "UUIDs"]
 git-tree-sha1 = "9abb12b62debc27261c008daa13627255bf79967"
@@ -1400,6 +1464,12 @@ git-tree-sha1 = "70e830dab5d0775183c99fc75e4c24c614ed7142"
 uuid = "f1f71cc9-e9ae-5b93-9b94-4fe0e1ad3748"
 version = "5.5.1+0"
 
+[[deps.MUMPS_seq_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "METIS_jll", "libblastrampoline_jll"]
+git-tree-sha1 = "85047ac569761e3387717480a38a61d2a67df45c"
+uuid = "d7ed1dd3-d0ae-5e8e-bfb4-87a502085b8d"
+version = "500.700.300+0"
+
 [[deps.MacroTools]]
 deps = ["Markdown", "Random"]
 git-tree-sha1 = "2fa9ee3e63fd3a4f7a9a4f4744a52f4856de82df"
@@ -1422,6 +1492,12 @@ deps = ["Base64"]
 uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
 version = "1.11.0"
 
+[[deps.MathOptInterface]]
+deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "DataStructures", "ForwardDiff", "JSON", "LinearAlgebra", "MutableArithmetics", "NaNMath", "OrderedCollections", "PrecompileTools", "Printf", "SparseArrays", "SpecialFunctions", "Test", "Unicode"]
+git-tree-sha1 = "c7fff4fbc415a6833988efbda7621039af2a5f36"
+uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
+version = "1.32.0"
+
 [[deps.MbedTLS]]
 deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "NetworkOptions", "Random", "Sockets"]
 git-tree-sha1 = "c067a280ddc25f196b5e7df3877c6b226d390aaf"
@@ -1483,6 +1559,12 @@ version = "0.3.4"
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
 version = "2023.12.12"
 
+[[deps.MutableArithmetics]]
+deps = ["LinearAlgebra", "SparseArrays", "Test"]
+git-tree-sha1 = "90077f1e79de8c9c7c8a90644494411111f4e07b"
+uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+version = "1.5.2"
+
 [[deps.NLSolversBase]]
 deps = ["DiffResults", "Distributed", "FiniteDiff", "ForwardDiff"]
 git-tree-sha1 = "a0b464d183da839699f4c79e7606d9d186ec172c"
@@ -1583,6 +1665,12 @@ git-tree-sha1 = "a5a5a68d079ce531b0220e99789e0c1c8c5ed215"
 uuid = "925886fa-5bf2-5e8e-b522-a9147a512338"
 version = "1.7.1"
 
+[[deps.OpenBLAS32_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"]
+git-tree-sha1 = "dd806c813429ff09878ea3eeb317818f3ca02871"
+uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2"
+version = "0.3.28+3"
+
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
@@ -1646,13 +1734,11 @@ deps = ["Compat", "FillArrays", "ForwardDiff", "LineSearches", "LinearAlgebra",
 git-tree-sha1 = "d9b79c4eed437421ac4285148fcadf42e0700e89"
 uuid = "429524aa-4258-5aef-a3af-852621145aeb"
 version = "1.9.4"
+weakdeps = ["MathOptInterface"]
 
     [deps.Optim.extensions]
     OptimMOIExt = "MathOptInterface"
 
-    [deps.Optim.weakdeps]
-    MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
-
 [[deps.OrderedCollections]]
 git-tree-sha1 = "dfdf5519f235516220579f949664f1bf44e741c5"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
@@ -1777,6 +1863,10 @@ deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 version = "1.11.0"
 
+[[deps.Profile]]
+uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
+version = "1.11.0"
+
 [[deps.ProgressMeter]]
 deps = ["Distributed", "Printf"]
 git-tree-sha1 = "8f6bc219586aef8baf0ff9a5fe16ee9c70cb65e4"
@@ -1929,6 +2019,12 @@ git-tree-sha1 = "7cf3e4a9432defe1eaa9959c1c82e85485817d07"
 uuid = "cdeec39e-fb35-4959-aadb-a1dd5dede958"
 version = "2.0.1"
 
+[[deps.SPRAL_jll]]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Hwloc_jll", "JLLWrappers", "Libdl", "METIS_jll", "libblastrampoline_jll"]
+git-tree-sha1 = "11f3da4b25efacd1cec8e263421f2a9003a5e8e0"
+uuid = "319450e9-13b8-58e8-aa9f-8fd1420848ab"
+version = "2024.5.8+0"
+
 [[deps.SQLite_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"]
 git-tree-sha1 = "004fffbe2711abdc7263a980bbb1af9620781dd9"
@@ -1954,9 +2050,9 @@ version = "1.2.1"
 
 [[deps.SentinelArrays]]
 deps = ["Dates", "Random"]
-git-tree-sha1 = "ff11acffdb082493657550959d4feb4b6149e73a"
+git-tree-sha1 = "305becf8af67eae1dbc912ee9097f00aeeabb8d5"
 uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
-version = "1.4.5"
+version = "1.4.6"
 
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -2075,9 +2171,9 @@ version = "0.1.1"
 
 [[deps.StaticArrays]]
 deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"]
-git-tree-sha1 = "eeafab08ae20c62c44c8399ccb9354a04b80db50"
+git-tree-sha1 = "777657803913ffc7e8cc20f0fd04b634f871af8f"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.9.7"
+version = "1.9.8"
 weakdeps = ["ChainRulesCore", "Statistics"]
 
     [deps.StaticArrays.extensions]
@@ -2331,9 +2427,9 @@ version = "1.6.1"
 
 [[deps.XML2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
-git-tree-sha1 = "1165b0443d0eca63ac1e32b8c0eb69ed2f4f8127"
+git-tree-sha1 = "6a451c6f33a176150f315726eba8b92fbfdb9ae7"
 uuid = "02c8fc9c-b97f-50b9-bbe4-9be30ff0a78a"
-version = "2.13.3+0"
+version = "2.13.4+0"
 
 [[deps.XMLDict]]
 deps = ["EzXML", "IterTools", "OrderedCollections"]
@@ -2349,9 +2445,9 @@ version = "1.1.41+0"
 
 [[deps.XZ_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "ac88fb95ae6447c8dda6a5503f3bafd496ae8632"
+git-tree-sha1 = "15e637a697345f6743674f1322beefbc5dcd5cfc"
 uuid = "ffd25f8a-64ca-5728-b0f7-c24cf3aae800"
-version = "5.4.6+0"
+version = "5.6.3+0"
 
 [[deps.Xorg_libX11_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libxcb_jll", "Xorg_xtrans_jll"]
@@ -2421,9 +2517,9 @@ version = "0.7.3"
 
 [[deps.YAXArrays]]
 deps = ["CFTime", "DataStructures", "Dates", "DimensionalData", "DiskArrayTools", "DiskArrays", "Distributed", "DocStringExtensions", "Glob", "Interpolations", "IntervalSets", "IterTools", "Markdown", "OffsetArrays", "OnlineStats", "Optim", "OrderedCollections", "ParallelUtilities", "ProgressMeter", "Reexport", "Requires", "Statistics", "StatsBase", "Tables", "WeightedOnlineStats", "YAXArrayBase"]
-git-tree-sha1 = "09cd01df4817c4d5f6eda845c4d75cc2855cf0b3"
+git-tree-sha1 = "02d9149f67b26256ffa1e1c3915d474fae57099e"
 uuid = "c21b50f5-aa40-41ea-b809-c0f5e47bfa5c"
-version = "0.5.11"
+version = "0.5.12"
 
 [[deps.Zarr]]
 deps = ["AWSS3", "Blosc", "CodecZlib", "DataStructures", "Dates", "DiskArrays", "HTTP", "JSON", "LRUCache", "OffsetArrays", "OpenSSL", "Pkg", "URIs", "ZipArchives"]
diff --git a/Project.toml b/Project.toml
index 69edb153c..d005a4098 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,12 +5,14 @@ version = "0.11.0"
 
 [deps]
 ArchGDAL = "c9ce4bd3-c3d5-55b8-8973-c0e20141b8c3"
+BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209"
 Bootstrap = "e28b5b4c-05e8-5b66-bc03-6f0c0a0a06e0"
 CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
 Clustering = "aaaa29a8-35af-508c-8bc3-b662a17a0fe5"
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
 CoralBlox = "f00a6e10-004a-11ef-3676-67ec635bc1a2"
 CpuId = "adafc99b-e345-5852-983c-f28acb93d879"
+DataEnvelopmentAnalysis = "a100299e-89d6-11e9-0fa0-2daf497e6a05"
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 Dates = "ade2ca70-3891-5945-98fb-dc099432e06a"
diff --git a/docs/src/assets/imgs/analysis/example_dea_fig.png b/docs/src/assets/imgs/analysis/example_dea_fig.png
new file mode 100644
index 000000000..f0ca72689
Binary files /dev/null and b/docs/src/assets/imgs/analysis/example_dea_fig.png differ
diff --git a/docs/src/usage/analysis.md b/docs/src/usage/analysis.md
index 9d47aa63e..8ae15280b 100644
--- a/docs/src/usage/analysis.md
+++ b/docs/src/usage/analysis.md
@@ -488,6 +488,42 @@ save("outcome_map.png", tf)
 
 ![Outcome mapping](../assets/imgs/analysis/outcome_map.png)
 
+### Data Envelopment Analysis
+
+Performs output-oriented (default, input-oriented can also be applied) Data Envelopment Analysis (DEA)
+given inputs X and output metrics Y. DEA is used to measure the performance of entities (scenarios),
+where inputs are converted to outputs via some process. Each scenario's "efficiency score" is calculated
+relative to an "efficiency fromtier", a region representing scenarios for which outputs cannot be further
+increased by changing inputs (scenario settings).
+
+```julia
+dom = ADRIA.load_domain("path to domain", "45")
+
+scens = ADRIA.sample(dom, 128)
+rs = ADRIA.run_scenarios(dom, scens, "45")
+
+n_scens = size(scens,1)
+
+# Get cost of deploying corals in each scenario, with user-specified function
+cost = cost_function(scens)
+
+# Get mean coral cover and shelter volume for each scenario
+s_tac = dropdims(
+    mean(ADRIA.metrics.scenario_total_cover(rs); dims=:timesteps); dims=:timesteps
+)
+s_sv =
+    dropdims(
+        mean(mean(ADRIA.metrics.absolute_shelter_volume(rs); dims=:timesteps); dims=:locations);
+        dims=(:timesteps,:locations)
+    )
+
+# Do output oriented DEA analysis seeking to maximise cover and shelter volume for minimum
+# deployment cost.
+DEA_scens = ADRIA.analysis.data_envelopment_analysis(cost, s_tac, s_sv)
+dea_fig = ADRIA.viz.data_envelopment_analysis(rs, DEA_scens)
+
+![DEA](../assets/imgs/analysis/example_dea_fig.png)
+```
 ### GUI for high-level exploration (prototype only!)
 
 ```julia
diff --git a/ext/AvizExt/viz/data_envelopment.jl b/ext/AvizExt/viz/data_envelopment.jl
new file mode 100644
index 000000000..ef752e930
--- /dev/null
+++ b/ext/AvizExt/viz/data_envelopment.jl
@@ -0,0 +1,150 @@
+using ADRIA.analysis: DEAResult
+
+"""
+    ADRIA.viz.data_envelopment_analysis(rs::ResultSet, DEA_output::DEAResult; axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(), fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(), opts::OPT_TYPE=DEFAULT_OPT_TYPE())
+    ADRIA.viz.data_envelopment_analysis(g::Union{GridLayout,GridPosition}, rs::ResultSet, DEA_output::DEAResult; axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(), opts::OPT_TYPE=DEFAULT_OPT_TYPE())
+    ADRIA.viz.data_envelopment_analysis(g::Union{GridLayout,GridPosition}, DEA_output::DEAResult; axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(), opts::OPT_TYPE=DEFAULT_OPT_TYPE())
+
+Plot results from a DEA analysis. Plots the first 2 dimensions of the effificency frontier,
+along side the technical and scale efficiencies.
+
+# Examples
+```
+# Run scenarios
+scens = ADRIA.sample(dom, 2^12)
+rs = ADRIA.run_scenarios(dom, scens, ["45"])
+
+# Compute cost an mean metrics for each scenario
+cost = cost_function(scens)
+s_tac::Vector{Float64} = Array(
+    dropdims(
+        mean(ADRIA.metrics.scenario_total_cover(rs); dims=:timesteps); dims=:timesteps
+    )
+)
+s_sv::Vector{Float64} = Array(
+    dropdims(
+        mean(mean(ADRIA.metrics.absolute_shelter_volume(rs); dims=:timesteps); dims=:sites);
+        dims=(:timesteps, :sites)
+    )
+)
+
+# Apply DEA analysis
+DEA_output = ADRIA.data_envelopment_analysis(cost, s_tac, s_sv)
+
+# Plot frontier, scale and technical efficiencies
+ADRIA.viz.data_envelopment_analysis(rs, DEA_output)
+```
+
+# Arguments
+- `rs` : ResultSet
+- `DEA_output` : output structure from a DEA analysis, carried out using `data_envelopment_analysis`
+- `axis_opts` : Additional options to pass to adjust Axis attributes
+  See: https://docs.makie.org/v0.19/api/index.html#Axis
+- `fig_opts` : Additional options to pass to adjust Figure creation
+  See: https://docs.makie.org/v0.19/api/index.html#Figure
+- `opts` : Aviz options
+    - `frontier_type`, type of frontier to plot: "CRS","VRS" or "FDH".
+    - `line_color`, color to use for best practice frontier.
+    - `data_color`, color to use for data cloud when plotting efficiency frontier.
+    - `frontier_name`, text to label efficiency frontier.
+    - `data_name`, text to label data cloud used to plot efficiency frontier.
+    - `metrics_x_lab` : Name of metric on x axis.
+    - `metrics_y_lab` : Name of metric on y axis.
+    - `scale_eff_y_lab` : Label for scale efficiency.
+    - `scale_eff_y_lab` : Label for technical efficiency.
+"""
+function ADRIA.viz.data_envelopment_analysis(
+    rs::ResultSet,
+    DEA_output::DEAResult;
+    axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
+    fig_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
+    opts::OPT_TYPE=DEFAULT_OPT_TYPE()
+)
+    f = Figure(; fig_opts...)
+    g = f[1, 1] = GridLayout()
+
+    ADRIA.viz.data_envelopment_analysis!(
+        g, rs, DEA_output; opts=opts, axis_opts=axis_opts
+    )
+
+    return f
+end
+function ADRIA.viz.data_envelopment_analysis!(
+    g::Union{GridLayout,GridPosition},
+    rs::ResultSet,
+    DEA_output::DEAResult;
+    axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
+    opts::OPT_TYPE=DEFAULT_OPT_TYPE()
+)
+    return ADRIA.viz.data_envelopment_analysis!(
+        g, DEA_output; opts=opts, axis_opts=axis_opts
+    )
+end
+function ADRIA.viz.data_envelopment_analysis!(
+    g::Union{GridLayout,GridPosition},
+    DEA_output::DEAResult;
+    axis_opts::OPT_TYPE=DEFAULT_OPT_TYPE(),
+    opts::OPT_TYPE=DEFAULT_OPT_TYPE())
+    scatter_colors = get(opts, :scatter_colors, [:red, :black])
+    legend_names = get(
+        opts, :legend_names, ["Best practice frontier", "Scenario data cloud"]
+    )
+    scale_eff_y_lab = get(opts, :scale_eff_y_lab, L"$\frac{eff_{vrs}}{eff_{crs}}$")
+    tech_eff_y_lab = get(opts, :tech_eff_y_lab, L"$\frac{1}{eff_{vrs}}$")
+    metrics_x_lab = get(opts, :metrics_x_lab, L"$metric 1$")
+    metrics_y_lab = get(opts, :metrics_y_lab, L"$metric 2$")
+
+    # Determines which returns to scale approach is used to select scenario peers
+    # (most efficient scenarios)
+    frontier_type = get(opts, :frontier_type, :vrs_peers)
+
+    Y = DEA_output.Y # Output values
+
+    # Find points on best practice frontier
+    best_practice_scens = getfield(DEA_output, frontier_type).J
+
+    scale_efficiency = DEA_output.crs_vals ./ DEA_output.vrs_vals
+
+    # Plot efficiency frontier and data cloud
+    ax_a = Axis(g[1, 1]; xlabel=metrics_x_lab, ylabel=metrics_y_lab, axis_opts...)
+    scatter!(ax_a, Y[:, 1], Y[:, 2]; color=scatter_colors[2])
+    scatter!(
+        ax_a, Y[best_practice_scens, 1], Y[best_practice_scens, 2]; color=scatter_colors[1]
+    )
+    _render_marker_legend(g, (1, 2), legend_names, scatter_colors)
+
+    # Plot the scale efficiency (ratio of efficiencies assuming CRS vs. assuming VRS)
+    ax_b = Axis(g[2, 1]; title="Scale efficiency", ylabel=scale_eff_y_lab, axis_opts...)
+    scatter!(ax_b, scale_efficiency; color=scatter_colors[2])
+    scatter!(
+        ax_b,
+        best_practice_scens,
+        scale_efficiency[best_practice_scens];
+        color=scatter_colors[1]
+    )
+
+    # Plot the technical efficiency (inverse VRS efficiencies)
+    ax_c = Axis(g[3, 1]; title="Technical efficiency", ylabel=tech_eff_y_lab, axis_opts...)
+    scatter!(ax_c, DEA_output.vrs_vals; color=scatter_colors[2])
+    scatter!(
+        ax_c,
+        best_practice_scens,
+        DEA_output.vrs_vals[best_practice_scens];
+        color=scatter_colors[1]
+    )
+
+    return g
+end
+
+function _render_marker_legend(g::Union{GridLayout,GridPosition},
+    legend_position::Tuple{Int64,Int64},
+    legend_labels::Union{Vector{Symbol},Vector{String}},
+    colors::Union{Vector{Symbol},RGBA{Float32}}
+)::Nothing
+    marker_els::Vector{MarkerElement} = [
+        MarkerElement(; color=color, marker=:circle) for color in colors
+    ]
+    Legend(g[legend_position...], marker_els, legend_labels; framevisible=false)
+
+    return nothing
+end
diff --git a/ext/AvizExt/viz/viz.jl b/ext/AvizExt/viz/viz.jl
index ea09b61bb..635a75443 100644
--- a/ext/AvizExt/viz/viz.jl
+++ b/ext/AvizExt/viz/viz.jl
@@ -95,3 +95,4 @@ include("spatial.jl")
 include("taxa_dynamics.jl")
 include("environment/dhw.jl")
 include("environment/cyclones.jl")
+include("data_envelopment.jl")
diff --git a/src/analysis/analysis.jl b/src/analysis/analysis.jl
index ca6602065..d61876196 100644
--- a/src/analysis/analysis.jl
+++ b/src/analysis/analysis.jl
@@ -114,5 +114,6 @@ include("intervention.jl")
 include("clustering.jl")
 include("rule_extraction.jl")
 include("scenario.jl")
+include("data_envelopment.jl")
 
 end
diff --git a/src/analysis/data_envelopment.jl b/src/analysis/data_envelopment.jl
new file mode 100644
index 000000000..e677ef22f
--- /dev/null
+++ b/src/analysis/data_envelopment.jl
@@ -0,0 +1,143 @@
+using DataEnvelopmentAnalysis: DataEnvelopmentAnalysis as DEA
+using ADRIA: ResultSet
+using DataFrames, YAXArrays
+
+struct DEAResult{V,V2,V3}
+    crs_vals::V # Inverse efficiencies using constant returns to scale.
+    vrs_vals::V # Inverse efficiencies using variable returns to scale.
+    fdh_vals::V # Inverse efficiencies using free disposability hull (non-convexity assumption).
+    crs_peers::V2 # Scenarios on the efficiency frontier.
+    vrs_peers::V2 # Scenarios on the efficiency frontier.
+    fdh_peers::V2 # Scenarios on the efficiency frontier.
+    X::V # Inputs
+    Y::V3 # Outputs
+end
+
+"""
+    data_envelopment_analysis(rs::ResultSet, input_function::Function, metrics...; orient::Symbol=:Output, dea_model::Function=DEA.deabigdata)::DEAResult
+    data_envelopment_analysis(rs::ResultSet, Y::AbstractArray, input_function::Function; orient::Symbol=:Output, dea_model::Function=DEA.deabigdata)::DEAResult
+    data_envelopment_analysis(X::YAXArray, metrics...; orient::Symbol=:Output, dea_model::Function=DEA.deabigdata)::DEAResult
+    data_envelopment_analysis(X::YAXArray, Y::AbstractArray; orient::Symbol=:Output, dea_model::Function=DEA.deabigdata)::DEAResult
+    data_envelopment_analysis(X::Union{Vector{Float64}, Matrix{Float64}}, Y::Matrix{Float64}; dea_model::Function=DEA.deabigdata)::DEAResult
+
+Performs output-oriented (default) Data Envelopment Analysis (DEA) given inputs X and output
+metrics Y. DEA is used to measure the performance of entities (scenarios), where inputs are
+converted to outputs via some process. Each scenario's "efficiency score" is calculated
+relative to an "efficiency fromtier", a region representing scenarios for which outputs
+cannot be further increased by changing inputs (scenario settings). Scenarios on the
+frontier serve as "benchmarks" or "peers", associated with best practice restoration
+scenarios. Scenarios with efficiencies not equal to 1 can be improved to be more efficient.
+
+# Arguments
+- `rs` : ADRIA ResultSet
+- `input_function` : function which calculates an input for each scenario (e.g. cost, effort) given the scenario
+    dataframe as input.
+- `X` : Model inputs for each scenario (usually costs).
+- `Y` : Model outputs for each scenario (metrics such as tac, rci etc.).
+- `orient` : Orientation of the analysis. Can be output oriented (`orient=:Output`), which
+    seeks to maximise outputs for a given level of input, or input oriented (`orient=:Input`),
+    which seeks to minimise an input for a given level of output.
+- `dea_model` : model to use to calculate DEA frontier (see https://javierbarbero.github.io/DataEnvelopmentAnalysis.jl/stable/)
+
+# Returns
+DEAResult, which summarizes inputs, outputs, efficiencies and peers for each scenario.
+
+# Examples
+```julia
+dom = ADRIA.load_domain("example_domain")
+scens = ADRIA.sample(dom, 128)
+rs = ADRIA.run_scenarios(dom, scens, "45")
+
+# Get cost of deploying corals in each scenario, with user-specified function
+cost = cost_function(scens)
+
+# Get mean coral cover and shelter volume for each scenario
+s_tac = dropdims(
+    mean(ADRIA.metrics.scenario_total_cover(rs); dims=:timesteps); dims=:timesteps
+)
+s_sv = dropdims(
+    mean(ADRIA.metrics.scenario_shelter_volume(rs); dims=:timesteps); dims=:timesteps
+    )
+
+# Do output oriented DEA analysis seeking to maximise cover and shelter volume for minimum
+# deployment cost.analysis.data_envelopment_analysis(cost, s_tac, s_sv)
+
+```
+# References
+1. Huguenin, J-M., 2012.
+   Data Envelopment Analysis (DEA): A pedagogical guide for decision makers in the public
+   sector. https://api.semanticscholar.org/CorpusID:188267263
+2. Pascoe, S., Cannard, T., Dowling, N.A., et. al., 2023.
+   Use of Data Envelopment Analysis (DEA) to assess management alternatives in the presence
+   of multiple objectives.
+   Marine Policy, 148, 105444.
+   https://doi.org/10.1016/j.marpol.2022.105444
+3. Pascoe, S., 2024.
+   On the use of Data Envelopment Analysis for Multi-Criteria Decision Analysis.
+   Algorithms, 17:89.
+   https://doi.org/10.3390/a17030089
+"""
+function data_envelopment_analysis(
+    rs::ResultSet,
+    input_function::Function,
+    metrics...;
+    orient::Symbol=:Output,
+    dea_model::Function=DEA.deabigdata
+)::DEAResult
+    return data_envelopment_analysis(
+        input_function(rs.inputs), metrics; orient=orient, dea_model=dea_model
+    )
+end
+function data_envelopment_analysis(
+    rs::ResultSet,
+    Y::AbstractArray,
+    input_function::Function;
+    orient::Symbol=:Output,
+    dea_model::Function=DEA.deabigdata
+)::DEAResult
+    return data_envelopment_analysis(
+        input_function(rs.inputs), Y; orient=orient, dea_model=dea_model
+    )
+end
+function data_envelopment_analysis(
+    X::YAXArray,
+    metrics...;
+    orient::Symbol=:Output,
+    dea_model::Function=DEA.deabigdata
+)::DEAResult
+    Y = Array(hcat(metrics...))
+    return data_envelopment_analysis(
+        X.data, Y; orient=orient, dea_model=dea_model
+    )
+end
+function data_envelopment_analysis(
+    X::YAXArray,
+    Y::AbstractArray;
+    orient::Symbol=:Output,
+    dea_model::Function=DEA.deabigdata
+)::DEAResult
+    return data_envelopment_analysis(
+        X.data, Array(Y); orient=orient, dea_model=dea_model
+    )
+end
+function data_envelopment_analysis(
+    X::Union{Vector{Float64},Matrix{Float64}},
+    Y::Matrix{Float64};
+    orient::Symbol=:Output,
+    dea_model::Function=DEA.deabigdata
+)::DEAResult
+    result_CRS = dea_model(X, Y; orient=orient, rts=:CRS)
+    result_VRS = dea_model(X, Y; orient=orient, rts=:VRS)
+    result_FDH = dea_model(X, Y; orient=orient, rts=:FDH)
+
+    return DEAResult(
+        1 ./ result_CRS.eff,
+        1 ./ result_VRS.eff,
+        1 ./ result_FDH.eff,
+        DEA.peers(result_CRS),
+        DEA.peers(result_VRS),
+        DEA.peers(result_FDH),
+        X,
+        Y
+    )
+end
diff --git a/src/metrics/scenario.jl b/src/metrics/scenario.jl
index 868fd4906..485290cdf 100644
--- a/src/metrics/scenario.jl
+++ b/src/metrics/scenario.jl
@@ -1,3 +1,6 @@
+using Bootstrap
+using Random
+
 """Scenario-level summaries.
 
 Note: Aggregates across the `site` dimension so trajectories over time for each scenario are
diff --git a/src/viz/viz.jl b/src/viz/viz.jl
index 96a2b657e..79ce897e6 100644
--- a/src/viz/viz.jl
+++ b/src/viz/viz.jl
@@ -53,4 +53,8 @@ function cyclone_scenario() end
 function taxonomy() end
 function taxonomy!() end
 
+# Economics
+function data_envelopment_analysis() end
+function data_envelopment_analysis!() end
+
 end  # module
diff --git a/test/analysis.jl b/test/analysis.jl
index 8e1f38e98..120a3bdb5 100644
--- a/test/analysis.jl
+++ b/test/analysis.jl
@@ -1,7 +1,7 @@
 using WGLMakie, GeoMakie, GraphMakie
 using ADRIA
 using ADRIA.metrics: total_absolute_cover
-using Statistics
+using Statistics, YAXArrays, Distributions
 
 if !@isdefined(TEST_RS)
     const TEST_DOM, TEST_N_SAMPLES, TEST_SCENS, TEST_RS = test_rs()
@@ -316,6 +316,32 @@ function test_rs_w_fig(rs::ADRIA.ResultSet, scens::ADRIA.DataFrame)
         ))
     # save("outcome_map.png", tf)
 
+    ## Data Envelopement Analysis
+    n_scens = size(scens, 1)
+
+    # Get cost of deploying corals in each scenario, with user-specified function
+    cost = YAXArray(
+        collect(range(100; stop=1000000, length=n_scens)) .+
+        rand(Uniform(1000, 2000), n_scens)
+    )
+
+    # Get mean coral cover and shelter volume for each scenario
+    s_tac = dropdims(
+        mean(ADRIA.metrics.scenario_total_cover(rs); dims=:timesteps); dims=:timesteps
+    )
+    s_sv = dropdims(
+        mean(
+            mean(ADRIA.metrics.absolute_shelter_volume(rs); dims=:timesteps);
+            dims=:locations
+        );
+        dims=(:timesteps, :locations)
+    )
+
+    # Do output oriented DEA analysis seeking to maximise cover and shelter volume for minimum
+    # deployment cost.
+    DEA_scens = ADRIA.analysis.data_envelopment_analysis(cost, s_tac, s_sv)
+    dea_fig = ADRIA.viz.data_envelopment_analysis(rs, DEA_scens)
+
     return rs
 end