diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml
new file mode 100644
index 0000000..0f66259
--- /dev/null
+++ b/.github/workflows/CompatHelper.yml
@@ -0,0 +1,22 @@
+name: CompatHelper
+
+on:
+  schedule:
+    - cron: '00 00 * * *'
+  workflow_dispatch:
+
+jobs:
+  CompatHelper:
+    runs-on: ${{ matrix.os }}
+    strategy:
+      matrix:
+        julia-version: [1.2.0]
+        julia-arch: [x86]
+        os: [ubuntu-latest]
+    steps:
+      - name: Pkg.add("CompatHelper")
+        run: julia -e 'using Pkg; Pkg.add("CompatHelper")'
+      - name: CompatHelper.main()
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+        run: julia -e 'using CompatHelper; CompatHelper.main()'
diff --git a/.github/workflows/TagBot.yml b/.github/workflows/TagBot.yml
index d77d3a0..f389611 100644
--- a/.github/workflows/TagBot.yml
+++ b/.github/workflows/TagBot.yml
@@ -1,11 +1,20 @@
 name: TagBot
 on:
-  schedule:
-    - cron: 0 * * * *
+  issue_comment:
+    types:
+      - created
+  workflow_dispatch:
+    inputs:
+      lookback:
+        default: 3
+permissions:
+  contents: write
 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 }}
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 0000000..69c7b31
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,47 @@
+name: CI
+on:
+  push:
+    branches: [master]
+  pull_request:
+    types: [opened, synchronize, reopened]
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
+    runs-on: ${{ matrix.os }}
+    strategy:
+      fail-fast: false
+      matrix:
+        include:
+          - version: '1'
+            os: ubuntu-latest
+            arch: x64
+          - version: '1'
+            os: ubuntu-latest
+            arch: x86
+          - version: '1.6'
+            os: ubuntu-latest
+            arch: x64
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@v1
+        with:
+          version: ${{ matrix.version }}
+          arch: ${{ matrix.arch }}
+      - uses: actions/cache@v1
+        env:
+          cache-name: cache-artifacts
+        with:
+          path: ~/.julia/artifacts
+          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
+          restore-keys: |
+            ${{ runner.os }}-test-${{ env.cache-name }}-
+            ${{ runner.os }}-test-
+            ${{ runner.os }}-
+      - uses: julia-actions/julia-buildpkg@v1
+      - uses: julia-actions/julia-runtest@v1
+        with:
+          depwarn: error
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v1
+        with:
+          file: lcov.info
diff --git a/.github/workflows/documentation.yml b/.github/workflows/documentation.yml
new file mode 100644
index 0000000..8eb83cc
--- /dev/null
+++ b/.github/workflows/documentation.yml
@@ -0,0 +1,26 @@
+name: Documentation
+on:
+  push:
+    branches: [master]
+    tags: '*'
+  pull_request:
+    types: [opened, synchronize, reopened]
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    env:
+      # Fix for Plots with GR backend, see https://github.com/JuliaPolyhedra/Polyhedra.jl/pull/265
+      GKSwstype: nul
+    steps:
+      - uses: actions/checkout@v2
+      - uses: julia-actions/setup-julia@latest
+        with:
+          # Build documentation on the latest Julia 1.x
+          version: '1'
+      - name: Install dependencies
+        run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
+      - name: Build and deploy
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} # For authentication with GitHub Actions token
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} # For authentication with SSH deploy key
+        run: julia --project=docs/ docs/make.jl
diff --git a/.gitignore b/.gitignore
index 8c960ec..3f02ca7 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,3 +1,4 @@
 *.jl.cov
 *.jl.*.cov
 *.jl.mem
+Manifest.toml
diff --git a/.travis.yml b/.travis.yml
deleted file mode 100644
index d7f732e..0000000
--- a/.travis.yml
+++ /dev/null
@@ -1,28 +0,0 @@
-# Documentation: http://docs.travis-ci.com/user/languages/julia/
-language: julia
-os:
-  - linux
-  - osx
-julia:
-  - 1.0
-  - 1.1
-notifications:
-  email: false
-# Add GMP for GLPK
-addons:
-  apt_packages:
-  - libgmp-dev
-after_success:
-  # push coverage results to Coveralls
-  - julia -e 'cd(Pkg.dir("StructDualDynProg")); Pkg.add("Coverage"); using Coverage; Coveralls.submit(Coveralls.process_folder())'
-  # push coverage results to Codecov
-  - julia -e 'cd(Pkg.dir("StructDualDynProg")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())'
-jobs:
-  include:
-    - stage: "Documentation"
-      julia: 1.0
-      os: linux
-      script:
-        - julia --project=docs -e 'using Pkg; Pkg.instantiate(); Pkg.add(PackageSpec(path=pwd()))'
-        - julia --project=docs --color=yes docs/make.jl
-      after_success: skip
diff --git a/Manifest.toml b/Manifest.toml
deleted file mode 100644
index 5dec9df..0000000
--- a/Manifest.toml
+++ /dev/null
@@ -1,327 +0,0 @@
-[[ArnoldiMethod]]
-deps = ["DelimitedFiles", "LinearAlgebra", "Random", "SparseArrays", "StaticArrays", "Test"]
-git-tree-sha1 = "2b6845cea546604fb4dca4e31414a6a59d39ddcd"
-uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
-version = "0.0.4"
-
-[[Base64]]
-uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
-
-[[BenchmarkTools]]
-deps = ["JSON", "Printf", "Statistics"]
-git-tree-sha1 = "90b73db83791c5f83155016dd1cc1f684d4e1361"
-uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-version = "0.4.3"
-
-[[BinDeps]]
-deps = ["Libdl", "Pkg", "SHA", "URIParser", "Unicode"]
-git-tree-sha1 = "66158ad56b4bf6cc8413b37d0b7bc52402682764"
-uuid = "9e28174c-4ba2-5203-b857-d8d62c4213ee"
-version = "1.0.0"
-
-[[BinaryProvider]]
-deps = ["Libdl", "SHA"]
-git-tree-sha1 = "5b08ed6036d9d3f0ee6369410b830f8873d4024c"
-uuid = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
-version = "0.5.8"
-
-[[Calculus]]
-deps = ["LinearAlgebra"]
-git-tree-sha1 = "f641eb0a4f00c343bbc32346e1217b86f3ce9dad"
-uuid = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
-version = "0.5.1"
-
-[[CodecBzip2]]
-deps = ["BinaryProvider", "Libdl", "TranscodingStreams"]
-git-tree-sha1 = "5db086e510c11b4c87d05067627eadb2dc079995"
-uuid = "523fee87-0ab8-5b00-afb7-3ecf72e48cfd"
-version = "0.6.0"
-
-[[CodecZlib]]
-deps = ["BinaryProvider", "Libdl", "TranscodingStreams"]
-git-tree-sha1 = "05916673a2627dd91b4969ff8ba6941bc85a960e"
-uuid = "944b1d66-785c-5afd-91f1-9de20f533193"
-version = "0.6.0"
-
-[[CommonSubexpressions]]
-deps = ["Test"]
-git-tree-sha1 = "efdaf19ab11c7889334ca247ff4c9f7c322817b0"
-uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
-version = "0.2.0"
-
-[[Compat]]
-deps = ["Base64", "Dates", "DelimitedFiles", "Distributed", "InteractiveUtils", "LibGit2", "Libdl", "LinearAlgebra", "Markdown", "Mmap", "Pkg", "Printf", "REPL", "Random", "Serialization", "SharedArrays", "Sockets", "SparseArrays", "Statistics", "Test", "UUIDs", "Unicode"]
-git-tree-sha1 = "ed2c4abadf84c53d9e58510b5fc48912c2336fbb"
-uuid = "34da2185-b29b-5c13-b0c7-acf172513d20"
-version = "2.2.0"
-
-[[CutPruners]]
-deps = ["Compat", "MathProgBase"]
-git-tree-sha1 = "e8105c939afcb8540a8525105120b1126bd1a5aa"
-uuid = "65d46eb8-70e9-5a30-bf48-2afa3a021b8f"
-version = "0.1.0"
-
-[[DataStructures]]
-deps = ["InteractiveUtils", "OrderedCollections"]
-git-tree-sha1 = "b7720de347734f4716d1815b00ce5664ed6bbfd4"
-uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.17.9"
-
-[[Dates]]
-deps = ["Printf"]
-uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
-
-[[DelimitedFiles]]
-deps = ["Mmap"]
-uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
-
-[[DiffResults]]
-deps = ["StaticArrays"]
-git-tree-sha1 = "da24935df8e0c6cf28de340b958f6aac88eaa0cc"
-uuid = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
-version = "1.0.2"
-
-[[DiffRules]]
-deps = ["NaNMath", "Random", "SpecialFunctions"]
-git-tree-sha1 = "10dca52cf6d4a62d82528262921daf63b99704a2"
-uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
-version = "1.0.0"
-
-[[Distributed]]
-deps = ["LinearAlgebra", "Random", "Serialization", "Sockets"]
-uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
-
-[[DocStringExtensions]]
-deps = ["LibGit2", "Markdown", "Pkg", "Test"]
-git-tree-sha1 = "88bb0edb352b16608036faadcc071adda068582a"
-uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-version = "0.8.1"
-
-[[ForwardDiff]]
-deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "NaNMath", "Random", "SpecialFunctions", "StaticArrays"]
-git-tree-sha1 = "88b082d492be6b63f967b6c96b352e25ced1a34c"
-uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "0.10.9"
-
-[[GLPK]]
-deps = ["BinaryProvider", "Libdl", "MathOptInterface", "SparseArrays"]
-git-tree-sha1 = "3420033e843e140d9237238d69937a5bc7292e5a"
-uuid = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
-version = "0.12.1"
-
-[[HTTP]]
-deps = ["Base64", "Dates", "IniFile", "MbedTLS", "Sockets"]
-git-tree-sha1 = "5c49dab19938b119fe204fd7d7e8e174f4e9c68b"
-uuid = "cd3eb016-35fb-5094-929b-558a96fad6f3"
-version = "0.8.8"
-
-[[Inflate]]
-deps = ["Pkg", "Printf", "Random", "Test"]
-git-tree-sha1 = "b7ec91c153cf8bff9aff58b39497925d133ef7fd"
-uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
-version = "0.1.1"
-
-[[IniFile]]
-deps = ["Test"]
-git-tree-sha1 = "098e4d2c533924c921f9f9847274f2ad89e018b8"
-uuid = "83e8ac13-25f8-5344-8a64-a9f2b223428f"
-version = "0.5.0"
-
-[[InteractiveUtils]]
-deps = ["LinearAlgebra", "Markdown"]
-uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-
-[[JSON]]
-deps = ["Dates", "Mmap", "Parsers", "Unicode"]
-git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e"
-uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "0.21.0"
-
-[[JSONSchema]]
-deps = ["BinaryProvider", "HTTP", "JSON"]
-git-tree-sha1 = "b0a7f9328967df5213691d318a03cf70ea8c76b1"
-uuid = "7d188eb4-7ad8-530c-ae41-71a32a6d4692"
-version = "0.2.0"
-
-[[JuMP]]
-deps = ["Calculus", "DataStructures", "ForwardDiff", "LinearAlgebra", "MathOptInterface", "NaNMath", "Random", "SparseArrays", "Statistics"]
-git-tree-sha1 = "a970a86abc924f2c126cdb4978a5e8923d0e7b22"
-uuid = "4076af6c-e467-56ae-b986-b466b2749572"
-version = "0.20.0"
-
-[[LibGit2]]
-uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
-
-[[Libdl]]
-uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
-
-[[LightGraphs]]
-deps = ["ArnoldiMethod", "DataStructures", "Distributed", "Inflate", "LinearAlgebra", "Random", "SharedArrays", "SimpleTraits", "SparseArrays", "Statistics"]
-git-tree-sha1 = "a0d4bcea4b9c056da143a5ded3c2b7f7740c2d41"
-uuid = "093fc24a-ae57-5d10-9952-331d41423f4d"
-version = "1.3.0"
-
-[[LinearAlgebra]]
-deps = ["Libdl"]
-uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-
-[[Logging]]
-uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
-
-[[MacroTools]]
-deps = ["DataStructures", "Markdown", "Random"]
-git-tree-sha1 = "e2fc7a55bb2224e203bbd8b59f72b91323233458"
-uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
-version = "0.5.3"
-
-[[Markdown]]
-deps = ["Base64"]
-uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
-
-[[MathOptInterface]]
-deps = ["BenchmarkTools", "CodecBzip2", "CodecZlib", "JSON", "JSONSchema", "LinearAlgebra", "MutableArithmetics", "OrderedCollections", "SparseArrays", "Test", "Unicode"]
-git-tree-sha1 = "793416d916e8fb1d16cdd016dabe8df87c636198"
-uuid = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
-version = "0.9.10"
-
-[[MathProgBase]]
-deps = ["Compat"]
-git-tree-sha1 = "3bf2e534e635df810e5f4b4f1a8b6de9004a0d53"
-uuid = "fdba3010-5040-5b88-9595-932c9decdf73"
-version = "0.7.7"
-
-[[MbedTLS]]
-deps = ["BinaryProvider", "Dates", "Distributed", "Libdl", "Random", "Sockets", "Test"]
-git-tree-sha1 = "2d94286a9c2f52c63a16146bb86fd6cdfbf677c6"
-uuid = "739be429-bea8-5141-9913-cc70e7f3736d"
-version = "0.6.8"
-
-[[Mmap]]
-uuid = "a63ad114-7e13-5084-954f-fe012c677804"
-
-[[MutableArithmetics]]
-deps = ["LinearAlgebra", "SparseArrays", "Test"]
-git-tree-sha1 = "f1c1ceab8adc1141ac5eba3881b0d53b277949c0"
-uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
-version = "0.2.2"
-
-[[NaNMath]]
-git-tree-sha1 = "928b8ca9b2791081dc71a51c55347c27c618760f"
-uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-version = "0.3.3"
-
-[[OrderedCollections]]
-deps = ["Random", "Serialization", "Test"]
-git-tree-sha1 = "c4c13474d23c60d20a67b217f1d7f22a40edf8f1"
-uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.1.0"
-
-[[ParameterJuMP]]
-deps = ["JuMP", "MathOptInterface", "SparseArrays"]
-git-tree-sha1 = "061bbb1cfaed57f3a65c5e688285f8da9bef0178"
-uuid = "774612a8-9878-5177-865a-ca53ae2495f9"
-version = "0.1.2"
-
-[[Parsers]]
-deps = ["Dates", "Test"]
-git-tree-sha1 = "d112c19ccca00924d5d3a38b11ae2b4b268dda39"
-uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "0.3.11"
-
-[[Pkg]]
-deps = ["Dates", "LibGit2", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
-uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-
-[[Printf]]
-deps = ["Unicode"]
-uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-
-[[REPL]]
-deps = ["InteractiveUtils", "Markdown", "Sockets"]
-uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
-
-[[Random]]
-deps = ["Serialization"]
-uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-
-[[SHA]]
-uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
-
-[[Serialization]]
-uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
-
-[[SharedArrays]]
-deps = ["Distributed", "Mmap", "Random", "Serialization"]
-uuid = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
-
-[[SimpleTraits]]
-deps = ["InteractiveUtils", "MacroTools"]
-git-tree-sha1 = "2bdf3b6300a9d66fe29ee8bb51ba100c4df9ecbc"
-uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
-version = "0.9.1"
-
-[[Sockets]]
-uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
-
-[[SparseArrays]]
-deps = ["LinearAlgebra", "Random"]
-uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-
-[[SpecialFunctions]]
-deps = ["BinDeps", "BinaryProvider", "Libdl"]
-git-tree-sha1 = "3bdd374b6fd78faf0119b8c5d538788dbf910c6e"
-uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "0.8.0"
-
-[[StaticArrays]]
-deps = ["LinearAlgebra", "Random", "Statistics"]
-git-tree-sha1 = "5a3bcb6233adabde68ebc97be66e95dcb787424c"
-uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "0.12.1"
-
-[[Statistics]]
-deps = ["LinearAlgebra", "SparseArrays"]
-uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-
-[[StochOptInterface]]
-deps = ["Compat", "DocStringExtensions", "LightGraphs", "TimerOutputs"]
-git-tree-sha1 = "81b9ad624b602e8440af3fd60eabef130d7bff09"
-uuid = "ddbffdd9-5b1a-5f20-b6eb-597a58a751db"
-version = "0.0.1"
-
-[[StructJuMP]]
-deps = ["JuMP", "LinearAlgebra", "MathOptInterface", "ParameterJuMP"]
-git-tree-sha1 = "84d8135eabb9c3b63f98de5fc75fe48458cc0827"
-repo-rev = "master"
-repo-url = "https://github.com/StructJuMP/StructJuMP.jl.git"
-uuid = "34f15cae-5318-50c9-93d3-9feadd34e321"
-version = "0.1.0"
-
-[[Test]]
-deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
-uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[[TimerOutputs]]
-deps = ["Printf"]
-git-tree-sha1 = "311765af81bbb48d7bad01fb016d9c328c6ede03"
-uuid = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
-version = "0.5.3"
-
-[[TranscodingStreams]]
-deps = ["Random", "Test"]
-git-tree-sha1 = "7c53c35547de1c5b9d46a4797cf6d8253807108c"
-uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa"
-version = "0.9.5"
-
-[[URIParser]]
-deps = ["Test", "Unicode"]
-git-tree-sha1 = "6ddf8244220dfda2f17539fa8c9de20d6c575b69"
-uuid = "30578b45-9adc-5946-b283-645ec420af67"
-version = "0.4.0"
-
-[[UUIDs]]
-deps = ["Random"]
-uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
-
-[[Unicode]]
-uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
diff --git a/Project.toml b/Project.toml
index 06ff7e9..3c88ebd 100644
--- a/Project.toml
+++ b/Project.toml
@@ -2,24 +2,29 @@ name = "StructDualDynProg"
 uuid = "5d41aeb0-57df-5a94-a3c4-60ba49600646"
 
 [deps]
-Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 CutPruners = "65d46eb8-70e9-5a30-bf48-2afa3a021b8f"
-GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
+DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
-LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MathProgBase = "fdba3010-5040-5b88-9595-932c9decdf73"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StochOptInterface = "ddbffdd9-5b1a-5f20-b6eb-597a58a751db"
 StructJuMP = "34f15cae-5318-50c9-93d3-9feadd34e321"
 TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f"
 
 [compat]
-Compat = "2, 3"
-JuMP = "= 0.20"
-StructJuMP = "= 0.1"
+CutPruners = "0.3"
+Graphs = "1"
+JuMP = "1"
+StochOptInterface = "0.1"
+StructJuMP = "0.3"
+TimerOutputs = "0.5"
 julia = "1"
 
 [extras]
+GLPK = "60bf3e95-4087-53dc-ae20-288a0d20c6a6"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["GLPK", "Test"]
diff --git a/REQUIRE b/REQUIRE
deleted file mode 100644
index 8ca1215..0000000
--- a/REQUIRE
+++ /dev/null
@@ -1,8 +0,0 @@
-julia 1.0
-CutPruners
-JuMP 0.19
-StructJuMP 0.1+
-StochOptInterface
-LightGraphs
-TimerOutputs
-Compat 1
diff --git a/appveyor.yml b/appveyor.yml
deleted file mode 100644
index 28d8437..0000000
--- a/appveyor.yml
+++ /dev/null
@@ -1,35 +0,0 @@
-environment:
-  matrix:
-# - julia_version: 1.0 # Disabled to speed up CI, tested on Linux with Travis
-  - julia_version: 1.1
-
-platform:
-  - x86 # 32-bit
-  - x64 # 64-bit
-
-# This is left as an example for the next big version switch.
-#matrix:
-#  allow_failures:
-#  - julia_version: latest
-
-branches:
-  only:
-    - master
-    - /release-.*/
-
-notifications:
-  - provider: Email
-    on_build_success: false
-    on_build_failure: false
-    on_build_status_changed: false
-
-install:
-  - ps: iex ((new-object net.webclient).DownloadString("https://raw.githubusercontent.com/JuliaCI/Appveyor.jl/version-1/bin/install.ps1"))
-
-build_script:
-  - echo "%JL_BUILD_SCRIPT%"
-  - C:\julia\bin\julia -e "%JL_BUILD_SCRIPT%"
-
-test_script:
-  - echo "%JL_TEST_SCRIPT%"
-  - C:\julia\bin\julia -e "%JL_TEST_SCRIPT%"
diff --git a/docs/Project.toml b/docs/Project.toml
index 53bc6f8..3a52a5d 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -2,4 +2,4 @@
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 
 [compat]
-Documenter = "~0.19"
+Documenter = "0.27"
diff --git a/docs/make.jl b/docs/make.jl
index 980f173..5eff082 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,8 +1,9 @@
 using Documenter, StructDualDynProg
 
 makedocs(
-    format = :html,
     sitename = "StructDualDynProg",
+    # See https://github.com/JuliaDocs/Documenter.jl/issues/868
+    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
     pages = [
         "Home" => "index.md",
         "Quick Start" => "quickstart.md",
@@ -10,14 +11,10 @@ makedocs(
     ],
     # The following ensures that we only include the docstrings from
     # this module for functions define in Base that we overwrite.
-    modules = [StructDualDynProg]
+    modules = [StructDualDynProg],
 )
 
 deploydocs(
     repo   = "github.com/JuliaStochOpt/StructDualDynProg.jl.git",
-    target = "build",
-    osname = "linux",
-    julia  = "1.0",
-    deps   = nothing,
-    make   = nothing
+    push_preview = true,
 )
diff --git a/src/SDDP/SDDP.jl b/src/SDDP/SDDP.jl
index a7894da..3436897 100644
--- a/src/SDDP/SDDP.jl
+++ b/src/SDDP/SDDP.jl
@@ -1,6 +1,6 @@
 module SDDP
 
-using Compat, Compat.LinearAlgebra
+using LinearAlgebra
 
 using TimerOutputs
 
diff --git a/src/SDDP/path.jl b/src/SDDP/path.jl
index c257e66..5237ec9 100644
--- a/src/SDDP/path.jl
+++ b/src/SDDP/path.jl
@@ -53,9 +53,9 @@ mutable struct SDDPPath{TT<:SOI.AbstractTransition, SolT<:SOI.AbstractSolution}
 end
 
 function meanstdpaths(paths::Vector{<:SDDPPath}, Ktot)
-    z = Compat.reduce(append!, Vector{Float64}[x.z for x in paths], init=Float64[])
-    proba = Compat.reduce(append!, Vector{Float64}[x.proba for x in paths], init=Float64[])
-    npaths = Compat.reduce(append!, Vector{Int}[x.K for x in paths], init=Int[])
+    z = reduce(append!, Vector{Float64}[x.z for x in paths], init=Float64[])
+    proba = reduce(append!, Vector{Float64}[x.proba for x in paths], init=Float64[])
+    npaths = reduce(append!, Vector{Int}[x.K for x in paths], init=Int[])
     meanstdpaths(z, proba, npaths, Ktot)
 end
 
diff --git a/src/StructProg/StructProg.jl b/src/StructProg/StructProg.jl
index 2bf52bd..b22f717 100644
--- a/src/StructProg/StructProg.jl
+++ b/src/StructProg/StructProg.jl
@@ -1,6 +1,6 @@
 module StructProg
 
-using Compat, Compat.LinearAlgebra, Compat.SparseArrays
+using LinearAlgebra, SparseArrays
 
 using TimerOutputs
 
@@ -21,7 +21,6 @@ include("cutgen.jl")
 
 # Generic implementation of Stochastic Program
 include("cutstore.jl")
-include("solver.jl")
 include("nlds.jl")
 include("graph.jl")
 
diff --git a/src/StructProg/cutstore.jl b/src/StructProg/cutstore.jl
index e45362d..6accd33 100644
--- a/src/StructProg/cutstore.jl
+++ b/src/StructProg/cutstore.jl
@@ -56,7 +56,7 @@ end
 
 function apply!(store::CutStore{S}) where S
     if !isempty(store.bnew)
-        if store.storecuts == :Yes || (store.storecuts != :No && Compat.reduce(|, store.needstored, init=false))
+        if store.storecuts == :Yes || (store.storecuts != :No && reduce(|, store.needstored, init=false))
             store.A = [store.A; store.Anew]
             store.b = [store.b; store.bnew]
             append!(store.authors, store.authorsnew)
@@ -82,7 +82,7 @@ function needstored!(store::CutStore, nlds)
 end
 function noneedstored!(store::CutStore{S}, nlds) where S
     store.needstored[findfirst(f->f[1] === nlds, store.followers)] = false
-    if store.storecuts == :IfNeededElseDelete && !Compat.reduce(|, store.needstored, init=false)
+    if store.storecuts == :IfNeededElseDelete && !reduce(|, store.needstored, init=false)
         store.A = spzeros(S, 0, size(store.A, 2))
         store.b = S[] # See julia#22225
         store.authors = NLDS{S}[]
diff --git a/src/StructProg/graph.jl b/src/StructProg/graph.jl
index f0421d8..5b97989 100644
--- a/src/StructProg/graph.jl
+++ b/src/StructProg/graph.jl
@@ -1,6 +1,4 @@
-using LightGraphs
-
-import MathProgBase
+using Graphs
 
 mutable struct NodeData{S}
     nlds::NLDS{S}
diff --git a/src/StructProg/interface.jl b/src/StructProg/interface.jl
index 961c34a..9f8f1ce 100644
--- a/src/StructProg/interface.jl
+++ b/src/StructProg/interface.jl
@@ -2,28 +2,33 @@ struct SDDPModelData
     nodes::Vector{Union{Nothing, Int}}
 end
 
-function createnode(sp::StochasticProgram, m::Model, t, num_stages, solver, parent, pruningalgo::AbstractCutPruningAlgo, cutgen::AbstractOptimalityCutGenerator, detectlb::Bool=true, newcut::Symbol=:InvalidateSolver)
+function createnode(sp::StochasticProgram, model::StructJuMP.StructuredModel, t,
+                    num_stages, solver, parent, pruningalgo::AbstractCutPruningAlgo,
+                    cutgen::AbstractOptimalityCutGenerator,
+                    model_to_data::Dict{StructJuMP.StructuredModel, SDDPModelData},
+                    detectlb::Bool=true, newcut::Symbol=:InvalidateSolver)
     # For each model, we need to create a different node for each t.
     # We store all these node in a nodes array per model
-    if !(:SDDP in keys(m.ext))
+    if !(model in keys(model_to_data))
         nodes = Vector{Union{Nothing, Int}}(undef, num_stages)
         fill!(nodes, nothing)
-        m.ext[:SDDP] = SDDPModelData(nodes)
+        model_to_data[model] = SDDPModelData(nodes)
     end
-    nodes = m.ext[:SDDP].nodes
+    nodes = model_to_data[model].nodes
     if nodes[t] === nothing
         # The last argument contains the categories (e.g. :Cont, :Int, :Bool, ...) but it is currently unused
-        c, T, W, h, C, K, _ = StructJuMP.conicconstraintdata(m)
-        newnodedata = NodeData(NLDS{Float64}(W,h,T,K,C,c,solver,pruningalgo, newcut), parent === nothing ? 0 : SOI.get(sp, SOI.Dimension(), parent))
+        param_model = StructJuMP.ParametrizedModel(model, solver)
+        newnodedata = NodeData(NLDS{Float64}(param_model, pruningalgo, newcut), parent === nothing ? 0 : SOI.get(sp, SOI.Dimension(), parent))
         newnode = SOI.add_scenario_node!(sp, newnodedata)
         SOI.set!(sp, CutGenerator(), newnode, cutgen)
         nodes[t] = newnode
-        struc = getStructure(m)
         if t < num_stages
-            num_scen = length(struc.children)
-            for id in keys(struc.children)
-                child = createnode(sp, struc.children[id], t+1, num_stages, solver, newnode, pruningalgo, cutgen, detectlb, newcut)
-                tr = SOI.add_scenario_transition!(sp, newnode, child, struc.probability[id])
+            num_scen = length(model.children)
+            for id in keys(model.children)
+                child = createnode(sp, model.children[id], t+1, num_stages,
+                                   solver, newnode, pruningalgo, cutgen,
+                                   model_to_data, detectlb, newcut)
+                tr = SOI.add_scenario_transition!(sp, newnode, child, model.probability[id])
                 if detectlb
                     SOI.set!(sp, SOI.TransitionObjectiveValueBound(), tr, SOI.get(sp, SOI.NodeObjectiveValueBound(), child))
                 end
@@ -34,7 +39,7 @@ function createnode(sp::StochasticProgram, m::Model, t, num_stages, solver, pare
 end
 
 """
-    stochasticprogram(m::JuMP.Model, num_stages, solver,
+    stochasticprogram(model::StructJuMP.StructuredModel, num_stages, solver,
                       pruningalgo::CutPruners.AbstractCutPruningAlgo,
                       cutgen::AbstractOptimalityCutGenerator=MultiCutGenerator(),
                       detectlb::Bool=true, newcut::Symbol=:InvalidateSolver)
@@ -44,20 +49,13 @@ The master problem is assumed to have model `m` and the scenarios are considered
 The `pruningalgo` is as defined in [CutPruners](https://github.com/JuliaPolyhedra/CutPruners.jl).
 If `cutgen` is `MultiCutGenerator`, one variable `θ_i` is created for each scenario. Otherwise, if `cutgen` is `AvgCutGenerator`, only one variable `θ` is created and it represents the expected value of the objective value of the scenarios. If `cutgen` is `NoOptimalityCut` then no `θ` is created, only use this option if the objective of all models is zero except for the master model.
 """
-function SOI.stochasticprogram(m::Model, num_stages, solver,
+function SOI.stochasticprogram(model::StructJuMP.StructuredModel, num_stages, solver,
                                pruningalgo::AbstractCutPruningAlgo,
                                cutgen::AbstractOptimalityCutGenerator=MultiCutGenerator(),
                                detectlb::Bool=true, newcut::Symbol=:InvalidateSolver)
     sp = StochasticProgram{Float64}(num_stages)
-    createnode(sp, m, 1, num_stages, solver, nothing, pruningalgo, cutgen, detectlb, newcut)
-    sp
-end
-
-function clear(m::Model)
-    if :SDDP in keys(m.ext)
-        pop!(m.ext, :SDDP)
-        for (id, child) in getStructure(m).children
-            clear(child)
-        end
-    end
+    createnode(sp, model, 1, num_stages, solver, nothing, pruningalgo, cutgen,
+               Dict{StructJuMP.StructuredModel, SDDPModelData}(),
+               detectlb, newcut)
+    return sp
 end
diff --git a/src/StructProg/nlds.jl b/src/StructProg/nlds.jl
index c8454fc..6e78a4c 100644
--- a/src/StructProg/nlds.jl
+++ b/src/StructProg/nlds.jl
@@ -6,7 +6,7 @@ mutable struct Solution <: SOI.AbstractSolution
     status::Symbol
     objval
     objvalx
-    objvalxuray
+    objvalxuray # unbouded ray
     x
     xuray # unbouded ray
     θ
@@ -73,12 +73,7 @@ end
 #     σ >= 0
 #     ρ >= 0
 mutable struct NLDS{S}
-    W::AbstractMatrix{S}
-    h::AbstractVector{S}
-    T::AbstractMatrix{S}
-    K
-    C
-    c::AbstractVector{S}
+    param_model::StructJuMP.ParametrizedModel
 
     # parent solution
     x_a::AbstractVector{S}
@@ -97,8 +92,6 @@ mutable struct NLDS{S}
 
     nx::Int
     nθ::Int
-    nπ::Int
-    πs::Vector{Int}
     # Number of feasibility cuts
     nσ::Int
     # Location of each feasibility cut in the constraint matrix
@@ -108,40 +101,26 @@ mutable struct NLDS{S}
     # Location of each optimality cut in the constraint matrix
     ρs::Vector{Vector{Int}}
 
-    model
-    loaded
-    solved
     prevsol::Union{Nothing, Solution}
 
-    newcut::Symbol
     pruningalgo::AbstractCutPruningAlgo
     FCpruner::AbstractCutPruner
     OCpruners::Vector
 
-    function NLDS{S}(W::AbstractMatrix{S}, h::AbstractVector{S}, T::AbstractMatrix{S}, K, C, c::AbstractVector{S}, solver, pruningalgo::AbstractCutPruningAlgo, newcut::Symbol=:AddImmediately) where S
-        nx = size(W, 2)
+    function NLDS{S}(model::StructJuMP.ParametrizedModel, pruningalgo::AbstractCutPruningAlgo) where S
+        nx = length(model.variable_map)
         nθ = 0
-        nπ = length(h)
         localFC = CutStore{S}(nx)
         localOC = CutStore{S}(nx)
         FCpruner = CutPruner{nx, S}(pruningalgo, :≥)
         OCpruners = typeof(FCpruner)[]
-        if false
-            model = MathProgBase.ConicModel(solver)
-        else
-            model = MathProgBase.LinearQuadraticModel(solver)
-        end
-        nlds = new{S}(W, h, T, K, C, c, S[], nothing, nothing, CutStore{S}[], CutStore{S}[], localFC, localOC, Float64[], BitSet(), Float64[], nothing, AvgCutGenerator(), nx, nθ, nπ, 1:nπ, 0, Int[], Int[], Vector{Int}[], model, false, false, nothing, newcut, pruningalgo, FCpruner, OCpruners)
+        nlds = new{S}(model, S[], nothing, nothing, CutStore{S}[], CutStore{S}[], localFC, localOC, Float64[], BitSet(), Float64[], nothing, AvgCutGenerator(), nx, nθ, 0, Int[], Int[], Vector{Int}[], nothing, pruningalgo, FCpruner, OCpruners)
         addfollower(localFC, (nlds, (:Feasibility, 0)))
         addfollower(localOC, (nlds, (:Optimality, 1)))
-        nlds
+       nlds
     end
 end
 
-function NLDS{S}(W::AbstractMatrix, h::AbstractVector, T::AbstractMatrix, K, C, c::AbstractVector, solver, pruningalgo::AbstractCutPruningAlgo, newcut::Symbol=:AddImmediately) where S
-    NLDS{S}(AbstractMatrix{S}(W), AbstractVector{S}(h), AbstractMatrix{S}(T), K, C, AbstractVector{S}(c), solver, pruningalgo, newcut)
-end
-
 function add_childT!(nlds, childT)
     if childT === nothing
         @assert nlds.childT === nothing
@@ -187,16 +166,32 @@ function getobjectivebound(nlds::NLDS)
     if !isempty(nlds.θfree)
         return -Inf
     end
-    for (cone, idxs) in nlds.C
-        for i in idxs
-            if (cone == :NonNeg && nlds.c[i]  < 0) ||
-               (cone == :NonPos && nlds.c[i]  > 0) ||
-               (cone == :Free   && nlds.c[i] != 0)
-               return -Inf
-           end
+    obj = JuMP.objective_function(nlds.model.model)
+    if obj isa JuMP.VariableRef
+        obj = convert(JuMP.AffExpr, obj)
+    end
+    if obj != JuMP.AffExpr
+        # TODO quadratic objective
+        @warn "Objective bound for objective of type $(typeof(obj)) not implemented yet."
+        return -Inf
+    end
+    curlb = 0.0
+    for (coef, var) in linear_terms(obj)
+        lb, ub = MOI.Utilities.get_bounds(backend(model), index(var))
+        if coef < 0
+            bound = ub
+        elseif coef > 0
+            bound = lb
+        else
+            continue
+        end
+        if !isfinite(bound)
+            return -Inf
+        else
+            curlb += coef * bound
         end
     end
-    return Eθlb(nlds)
+    return curlb + Eθlb(nlds)
 end
 function setθbound!(nlds::NLDS, i, θlb)
     if isfinite(θlb)
@@ -264,9 +259,9 @@ function getfeasibilitycuts(nlds::NLDS)
         end
         D
     end
-    cuts_D = Compat.reduce(vcat, map(i -> f(i), 1:length(nlds.childFC)), init=nlds.localFC.A)
-    cuts_d = Compat.reduce(vcat, map(x -> x.b, nlds.childFC), init=nlds.localFC.b)
-    mycut = Compat.reduce(vcat, map(x -> veceqeqeq(x.authors, nlds), nlds.childFC), init=veceqeqeq(nlds.localFC.authors, nlds))
+    cuts_D = reduce(vcat, map(i -> f(i), 1:length(nlds.childFC)), init=nlds.localFC.A)
+    cuts_d = reduce(vcat, map(x -> x.b, nlds.childFC), init=nlds.localFC.b)
+    mycut = reduce(vcat, map(x -> veceqeqeq(x.authors, nlds), nlds.childFC), init=veceqeqeq(nlds.localFC.authors, nlds))
     (cuts_D, cuts_d, mycut)
 end
 
@@ -307,10 +302,10 @@ function getoptimalitycuts(nlds::NLDS{S}) where S
         cuts_E = spzeros(S, 0, nlds.nx + nlds.nθ)
     end
     if nlds.nθ == length(nlds.childOC)
-        cuts_E = Compat.reduce(vcat, map(f, 1:length(nlds.childOC)), init=cuts_E)
+        cuts_E = reduce(vcat, map(f, 1:length(nlds.childOC)), init=cuts_E)
     end
-    cuts_e = Compat.reduce(vcat, map(x -> x.b, nlds.childOC), init=nlds.localOC.b)
-    mycut = Compat.reduce(vcat, map(x -> veceqeqeq(x.authors, nlds), nlds.childOC), init=veceqeqeq(nlds.localOC.authors, nlds))
+    cuts_e = reduce(vcat, map(x -> x.b, nlds.childOC), init=nlds.localOC.b)
+    mycut = reduce(vcat, map(x -> veceqeqeq(x.authors, nlds), nlds.childOC), init=veceqeqeq(nlds.localOC.authors, nlds))
     (cuts_E, cuts_e, mycut)
 end
 
@@ -334,7 +329,6 @@ function notifynewcuts(nlds::NLDS{S}, A::AbstractMatrix{S}, b::AbstractVector{S}
     addstatus = addcuts!(pruner, A, b, mine)
     npushed = sum(addstatus .> ncur)
     @assert ncur + npushed == ncuts(pruner)
-    cur = nlds.nσ + sum(nlds.nρ)
     if isfc
         nlds.nσ += npushed
         is = nlds.σs
@@ -342,36 +336,24 @@ function notifynewcuts(nlds::NLDS{S}, A::AbstractMatrix{S}, b::AbstractVector{S}
         nlds.nρ[i] += npushed
         is = nlds.ρs[i]
     end
-    append!(is, zeros(Int, npushed))
+    append!(is, Vector{JuMP.ConstraintRef}(undef, npushed))
+    x = [nlds.model.variable_map[i] for i in 1:nlds.nx]
+    θ = nlds.model.parameter_map[i]
     for j in 1:nnewcuts
-        if addstatus[j] > ncur
-            cur += 1
-            is[addstatus[j]] = cur
-        end
-        if addstatus[j] > 0 && nlds.loaded
-            if 1 <= addstatus[j] <= ncur || nlds.newcut == :InvalidateSolver
-                # Need to update cut
-                nlds.loaded = false
-                nlds.solved = false
-            elseif nlds.newcut == :AddImmediately
-                idx = collect(1:nlds.nx)
-                a = A[j,:]
-                if !isfc
-                    if i == 0
-                        push!(idx, nlds.nx+1)
-                    else
-                        push!(idx, nlds.nx+i)
-                    end
-                    a = _veccat(a, one(S), true)
-                end
-                _addconstr!(nlds.model, idx, a, b[j], :NonPos)
-                nlds.solved = false
+        if addstatus[j] > 0
+            if addstatus[j] <= ncur
+                JuMP.delete(is[addstatus[j]])
+            end
+            a = A[j, :]
+            β = b[j]
+            is[addstatus[j]] = if isfc
+                @constraint(nlds.model.model, dot(a, x) ≥ β)
             else
-                error("Invalid newcut option $(nlds.newcut)")
+                @constraint(nlds.model.model, θ ≥ β - dot(a, x))
             end
+            nlds.prevsol = nothing
         end
     end
-    @assert nlds.newcut == :InvalidateSolver || cur == nlds.nσ + sum(nlds.nρ)
     checkconsistency(nlds)
 end
 
@@ -381,48 +363,16 @@ function checkconsistency(nlds)
     for i in 1:nlds.nθ
         @assert length(nlds.ρs[i]) == nlds.nρ[i]
     end
-    ρs = Compat.reduce(append!, nlds.ρs, init=Int[])
+    ρs = reduce(append!, nlds.ρs, init=Int[])
     @assert sort([nlds.πs; nlds.nπ .+ nlds.σs; nlds.nπ .+ ρs]) == collect(1:(nlds.nπ + nlds.nσ + sum(nlds.nρ)))
 end
 
-function getrhs(nlds::NLDS{S}) where S
-    checkconsistency(nlds)
-    bs = [nlds.h - nlds.T * nlds.x_a]
-    Ks = [nlds.K]
-    Kcut = []
-    cur = 0
-    b = S[]
-    if !isempty(nlds.FCpruner)
-        append!(b, nlds.FCpruner.b)
-        nlds.nσ = ncuts(nlds.FCpruner)
-        nlds.σs = cur .+ (1:nlds.nσ)
-        cur += ncuts(nlds.FCpruner)
-        push!(Kcut, (:NonPos, nlds.σs))
-    end
-    for i in 1:nlds.nθ
-        if !isempty(nlds.OCpruners[i])
-            append!(b, nlds.OCpruners[i].b)
-            nlds.nρ[i] = ncuts(nlds.OCpruners[i])
-            nlds.ρs[i] = cur .+ (1:nlds.nρ[i])
-            cur += ncuts(nlds.OCpruners[i])
-            push!(Kcut, (:NonPos, nlds.ρs[i]))
-        end
-    end
-    if !isempty(Kcut)
-        push!(bs, b)
-        push!(Ks, Kcut)
-    end
-    checkconsistency(nlds)
-    bs, Ks
-end
-
 function setparentx(nlds::NLDS, x_a::AbstractVector, xuray_a, objvalxuray_a)
     nlds.x_a = x_a
     unbounded_a = xuray_a !== nothing
     if nlds.xuray_a !== nothing || unbounded_a
-        # FIXME do better when delvars!, ... are available in MPB
-        nlds.loaded = false
-        nlds.solved = false
+        error("TODO")
+        nlds.prevsol = nothing
     end
     if unbounded_a
         nlds.xuray_a = xuray_a
@@ -435,10 +385,15 @@ function setparentx(nlds::NLDS, x_a::AbstractVector, xuray_a, objvalxuray_a)
         if unbounded_a
             nlds.loaded = false
         else
+            values = nlds.T * x_a
+            @assert length(nlds.model.parameter_map) == length(values)
+            for (index, parameter) in nlds.model.parameter_map
+                JuMP.fix(parameter, values[index])
+            end
             bs, Ks = getrhs(nlds)
             _setconstrB!(nlds.model, bs, Ks)
         end
-        nlds.solved = false
+        nlds.prevsol = nothing
     end
 end
 
@@ -459,73 +414,9 @@ function computecuts!(nlds::NLDS)
     end
 end
 
-function getcutsDE(nlds::NLDS{S}) where S
-    checkconsistency(nlds)
-    nc = nlds.nσ + sum(nlds.nρ)
-    A = spzeros(S, nc, nlds.nx + nlds.nθ)
-    if !isempty(nlds.FCpruner)
-        A[nlds.σs, 1:nlds.nx] = nlds.FCpruner.A
-    end
-    for i in 1:nlds.nθ
-        if !isempty(nlds.OCpruners[i])
-            A[nlds.ρs[i], 1:nlds.nx] = nlds.OCpruners[i].A
-            A[nlds.ρs[i], nlds.nx .+ i] .= 1
-        end
-    end
-    A
-end
-
-function load!(nlds::NLDS{S}) where S
-    if !nlds.loaded
-        bigA = nlds.W
-        if nlds.nθ > 0
-            bigA = [bigA spzeros(size(bigA, 1), nlds.nθ)]
-        end
-        if nlds.xuray_a !== nothing
-            bigA = [bigA nlds.T * nlds.xuray_a]
-        end
-
-        # Needs to be done before getcutsDE
-        bs, Ks = getrhs(nlds)
-
-        computecuts!(nlds) # FIXME what is the use of this ?
-        A = getcutsDE(nlds)
-        if nlds.xuray_a !== nothing
-            A = [A spzeros(size(A, 1), 1)]
-        end
-        bigA = [bigA; A]
-
-        bigC = nlds.C
-        bigc = nlds.c
-        if nlds.nθ > 0
-            bigC = [bigC; θC(nlds)[2]]
-            if nlds.nθ == 1
-                bigc = [bigc; 1]
-            else
-                bigc = [bigc; nlds.proba]
-            end
-        end
-        if nlds.xuray_a !== nothing
-            bigC = [bigC; (:NonNeg, [nlds.nx+nlds.nθ+1])]
-            bigc = [bigc; nlds.objvalxuray_a]
-        end
-
-        _load!(nlds.model, bigc, bigA, bs, Ks, bigC)
-        nlds.loaded = true
-    end
-end
-
-function checkstatus(status::Symbol)
-    if status == :Error
-        error("The solver reported an error")
-    elseif status == :UserLimit
-        error("The solver reached iteration limit or timed out")
-    end
-end
-
 function Eθlb(nlds::NLDS)
     if nlds.nθ == 0
-        0.
+        0.0
     else
         if nlds.nθ == length(nlds.proba)
             dot(nlds.proba, nlds.θlb)
@@ -538,92 +429,29 @@ function Eθlb(nlds::NLDS)
 end
 
 function solve!(nlds::NLDS{S}) where S
-    load!(nlds)
-    if !nlds.solved
-        MathProgBase.optimize!(nlds.model)
-        status = MathProgBase.status(nlds.model)
-        checkstatus(status)
-        sol = Solution(status, _getobjval(nlds.model) + Eθlb(nlds))
-
-        if status == :Unbounded
-            uray = _getunboundedray(nlds.model)
-            sol.xuray = uray[1:nlds.nx]
-            sol.objvalxuray = dot(nlds.c, sol.xuray)
-            sol.θuray = uray[nlds.nx .+ (1:nlds.nθ)]
-
+    if nlds.prevsol === nothing
+        JuMP.optimize!(nlds.model.model)
+        if JuMP.primal_status(nlds.model.model) == MOI.INFEASIBILITY_CERTIFICATE
+            # The dual is infeasible but this may not mean that the primal is
+            # unbounded:
             # See https://github.com/JuliaOpt/Gurobi.jl/issues/80
-            c = MathProgBase.getobj(nlds.model)
-
-            MathProgBase.setobj!(nlds.model, zero(c))
-            MathProgBase.optimize!(nlds.model)
-            newstatus = MathProgBase.status(nlds.model)
-            @assert newstatus != :Unbounded # Now the objective is 0
-            checkstatus(newstatus)
-            if newstatus == :Infeasible
-                # We discard unbounded ray, infeasibility is more important
-                status = newstatus
-                sol.xuray = nothing
-                sol.objvalxuray = nothing
-                sol.θuray = nothing
-            else
-                @assert newstatus == :Optimal
-            end
-            if newstatus == :Infeasible
-                sol.status = :Infeasible
-            end
-        else
-            newstatus = status
-        end
-
-        if newstatus == :Infeasible || status != :Unbounded
-            # if infeasible dual + λ iray is dual feasible for any λ >= 0
-            # here I just take iray to build the feasibility cut
-            dual = _getdual(nlds.model)
-            if isempty(dual)
-                error("Dual vector returned by the solver is empty while the status is $status")
-            end
-            @assert length(dual) == nlds.nπ + nlds.nσ + sum(nlds.nρ)
-
-            π = dual[nlds.πs]
-            σρ = @view dual[nlds.nπ+1:end]
-
-            σ = σρ[nlds.σs]
-            addusage!(nlds.FCpruner, σ)
-            sol.σd = dot(σ, nlds.FCpruner.b)
-
-            sol.ρe = zero(S)
-            for i in 1:nlds.nθ
-                ρ = σρ[nlds.ρs[i]]
-                addusage!(nlds.OCpruners[i], ρ)
-                sol.ρe += dot(ρ, nlds.OCpruners[i].b)
-            end
-
-            sol.πT = vec(π' * nlds.T)
-            sol.πh = dot(π, nlds.h)
-        end
-
-        # Needs to be done after since if status is unbounded I do a resolve
-        if newstatus != :Infeasible
-            primal = _getsolution(nlds.model)
-            sol.x = primal[1:nlds.nx]
-            addposition!(nlds.FCpruner, sol.x)
-            for i in 1:nlds.nθ
-                addposition!(nlds.OCpruners[i], sol.x)
-            end
-            sol.objvalx = dot(nlds.c, sol.x)
-            sol.θ = θC(nlds)[1] .+ primal[nlds.nx .+ (1:nlds.nθ)]
-        end
-
-        if status == :Unbounded
-            # It needs to be done *after* getsolution for some solver (e.g. CPLEX)
-            MathProgBase.setobj!(nlds.model, c)
+            # We set the objective to zero and reoptimize to check
+            obj = JuMP.objective_function(nlds.model.model)
+            set_objective_function(nlds.model.model, zero(AffExpr))
+            JuMP.optimize!(nlds.model.model)
+            set_objective_function(nlds.model.model, obj)
         end
-
-        nlds.prevsol = sol
+        solution = StructJuMP.optimize(nlds.model)
+        nlds.prevsol = Solution(
+            solution.feasible ? :Optimal : :Infeasible,
+            solution.objective_value, nothing, nothing,
+            [solution.variable_value[nlds.model.variable_map[i]] for i in 1:nlds.nx],
+            nothing, )
     end
+    return
 end
 
 function getsolution(nlds::NLDS)
     solve!(nlds)
-    nlds.prevsol
+    return nlds.prevsol
 end
diff --git a/src/StructProg/solver.jl b/src/StructProg/solver.jl
deleted file mode 100644
index 194680a..0000000
--- a/src/StructProg/solver.jl
+++ /dev/null
@@ -1,118 +0,0 @@
-using MathProgBase
-
-function _setconstrB!(m::MathProgBase.AbstractConicModel, b, K)
-    error("Not supported")
-end
-
-function getLPconstrbounds(bs, Ks)
-    sumlen = sum(map(length, bs))
-    for i in 1:length(bs)
-        @assert BitSet(1:length(bs[i])) == Compat.reduce(∪, map(c -> BitSet(c[2]), Ks[i]), init=BitSet())
-    end
-    lb = Vector{Float64}(undef, sumlen)
-    ub = Vector{Float64}(undef, sumlen)
-    offset = 0
-    for i in 1:length(bs)
-        b = bs[i]
-        for (cone, idx) in Ks[i]
-            #if !(cone in [:Zero, :NonPos, :NonNeg])
-            #  error("This cone is not supported")
-            #end
-            offidx = offset .+ idx
-            if cone == :Zero || cone == :NonPos
-                lb[offidx] .= b[idx]
-            else
-                lb[offidx] .= -Inf
-            end
-            if cone == :Zero || cone == :NonNeg
-                ub[offidx] .= b[idx]
-            else
-                ub[offidx] .= Inf
-            end
-        end
-        offset += length(b)
-    end
-    lb, ub
-end
-
-function _setconstrB!(m::MathProgBase.AbstractLinearQuadraticModel, bs, Ks)
-    lb, ub = getLPconstrbounds(bs, Ks)
-    MathProgBase.setconstrLB!(m, lb)
-    MathProgBase.setconstrUB!(m, ub)
-end
-
-function _addconstr!(m::MathProgBase.AbstractConicModel, idx, a, β, cone)
-    error("Not supported")
-end
-
-function _addconstr!(m::MathProgBase.AbstractLinearQuadraticModel, idx, a, β, cone)
-    lb = -Inf
-    ub = Inf
-    if cone == :Zero || cone == :NonPos
-        lb = β
-    end
-    if cone == :Zero || cone == :NonNeg
-        ub = β
-    end
-    MathProgBase.addconstr!(m, idx, Vector(a), lb, ub)
-end
-
-function _load!(model::MathProgBase.AbstractConicModel, c, A, bs, Ks, C)
-    MathProgBase.loadproblem!(model, c, A, Compat.reduce(vcat, bs, init=Float64[]), Compat.reduce(vcat, Ks, init=[]), C)
-end
-
-function _load!(model::MathProgBase.AbstractLinearQuadraticModel, c, A, bs, Ks, C)
-    lb, ub = getLPconstrbounds(bs, Ks)
-    l  = Vector{Float64}(undef, size(A, 2))
-    u  = Vector{Float64}(undef, size(A, 2))
-    @assert BitSet(1:size(A, 2)) == Compat.reduce(∪, map(c -> BitSet(c[2]), C), init=BitSet())
-    for (cone, idx) in C
-        if !(cone  in [:Free, :NonPos, :NonNeg])
-            error("This cone is not supported")
-        end
-        if cone == :Free || cone == :NonPos
-            l[idx] .= -Inf
-        else
-            l[idx] .= 0
-        end
-        if cone == :Free || cone == :NonNeg
-            u[idx] .= Inf
-        else
-            u[idx] .= 0
-        end
-    end
-    # TODO sparse objective triggers a bug in Gurobi solver, report it
-    MathProgBase.loadproblem!(model, A, l, u, Vector(c), lb, ub, :Min)
-end
-
-function _getdual(model::MathProgBase.AbstractConicModel)
-    -MathProgBase.getdual(model)
-end
-function _getdual(model::MathProgBase.AbstractLinearQuadraticModel)
-    if MathProgBase.status(model) == :Infeasible
-        MathProgBase.getinfeasibilityray(model)
-    else
-        MathProgBase.getconstrduals(model)
-    end
-end
-function _getsolution(model)
-    MathProgBase.getsolution(model)
-end
-function _getunboundedray(model::MathProgBase.AbstractConicModel)
-    error("Unbounded Ray retrieval is unsupported for conic model")
-end
-function _getunboundedray(model::MathProgBase.AbstractLinearQuadraticModel)
-    MathProgBase.getunboundedray(model)
-end
-function _getobjval(model)
-    # We never know, a solver could return something else
-    # when it is :Unbounded or :Infeasible, thinking that
-    # objval should only be called when the status is :Optimal
-    if MathProgBase.status(model) == :Unbounded
-        -Inf
-    elseif MathProgBase.status(model) == :Infeasible
-        Inf
-    else
-        MathProgBase.getobjval(model)
-    end
-end
diff --git a/src/WaitAndSee/waitandsee.jl b/src/WaitAndSee/waitandsee.jl
index ea4ca58..06593ef 100644
--- a/src/WaitAndSee/waitandsee.jl
+++ b/src/WaitAndSee/waitandsee.jl
@@ -1,6 +1,6 @@
 module WaitAndSee
 
-using Compat, Compat.SparseArrays
+using SparseArrays
 
 using StructDualDynProg
 using StochOptInterface
@@ -81,20 +81,13 @@ function SOI.optimize!(sp::SOI.AbstractStochasticProgram, algo::Algorithm,
         end
         model = MathProgBase.LinearQuadraticModel(algo.solver)
         StructDualDynProg.StructProg._load!(model, c, A, bs, Ks, C)
-        MathProgBase.optimize!(model)
-        status = MathProgBase.status(model)
-        if status == :Error
-            error("The solver reported an error")
-        elseif status == :UserLimit
-            error("The solver reached iteration limit or timed out")
-        elseif status == :Infeasible
-            error("The problem is infeasible for some realization of the uncertainty")
-        elseif status == :Unbounded
-            error("The problem is unbounded for some realization of the uncertainty")
-        else
-            @assert status == :Optimal
-            objval = MathProgBase.getobjval(model)
+        MOI.optimize!(model)
+        status = MOI.get(model, MOI.TerminationStatus())
+        if status == MOI.OPTIMAL
+            objval = MOI.get(model, MOI.ObjectiveValue())
             push!(newpaths, WaitAndSeePath(path.node, path.nlds, objval, path.proba, path.K))
+        else
+            error("The solver terminated with $status")
         end
     end
     meanstdpaths(newpaths, algo.K)
diff --git a/test/hydro_thermal_scheduling.jl b/test/hydro_thermal_scheduling.jl
index 10a6535..2a0940c 100644
--- a/test/hydro_thermal_scheduling.jl
+++ b/test/hydro_thermal_scheduling.jl
@@ -24,7 +24,7 @@
 			if s > 1
 				@constraint(model, x[s, ξ] <= x[s-1, 1] + r[ξ] - y[s, ξ])
 			else
-                @constraint(model, x[s, ξ] <= Compat.Statistics.mean(r) - y[s, ξ])
+                @constraint(model, x[s, ξ] <= Statistics.mean(r) - y[s, ξ])
 			end
 			@constraint(model, p[s, ξ] + y[s, ξ] >= d)
 			@objective(model, Min, C * p[s, ξ])
diff --git a/test/runtests.jl b/test/runtests.jl
index df5d7bf..881c6d5 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,10 +1,10 @@
-using JuMP
+using LinearAlgebra, Test
 using CutPruners
+using JuMP
 using StructJuMP
-using StructDualDynProg
-using Compat, Compat.LinearAlgebra, Test
 using StochOptInterface
 const SOI = StochOptInterface
+using StructDualDynProg
 
 # solver independent tests
 include("comp.jl")
diff --git a/test/solvers.jl b/test/solvers.jl
index 49cf741..b90ab82 100644
--- a/test/solvers.jl
+++ b/test/solvers.jl
@@ -16,12 +16,14 @@ iscpx(solver) = occursin("CplexSolver", string(typeof(solver)))
 xpr = false && try_import(:Xpress)
 clp = false && try_import(:Clp)
 isclp(solver) = occursin("ClpSolver", string(typeof(solver)))
+glp = try_import(:GLPK)
 msk = false && try_import(:Mosek)
 ismsk(solver) = occursin("MosekSolver", string(typeof(solver)))
 
 lp_solvers = Any[]
-grb && push!(lp_solvers, Gurobi.GurobiSolver(OutputFlag=0))
-cpx && push!(lp_solvers, CPLEX.CplexSolver(CPX_PARAM_SCRIND=0, CPX_PARAM_REDUCE=0))
-xpr && push!(lp_solvers, Xpress.XpressSolver(OUTPUTLOG=0))
-clp && push!(lp_solvers, Clp.ClpSolver())
-msk && push!(lp_solvers, Mosek.MosekSolver(QUIET=true))
+grb && push!(lp_solvers, optimizer_with_attributes(Gurobi.Optimizer, "OutputFlag" => 0))
+cpx && push!(lp_solvers, optimizer_with_attributes(CPLEX.Optimizer, "CPX_PARAM_SCRIND" => 0, "CPX_PARAM_REDUCE" => 0))
+xpr && push!(lp_solvers, optimizer_with_attributes(Xpress.Optimizer, "OUTPUTLOG" => 0))
+clp && push!(lp_solvers, Clp.Optimizer)
+glp && push!(lp_solvers, GLPK.Optimizer)
+msk && push!(lp_solvers, optimizer_with_attributes(Mosek.Optimizer, "QUIET" => true))