diff --git a/.gitignore b/.gitignore
index 22459730..b22046d0 100644
--- a/.gitignore
+++ b/.gitignore
@@ -45,4 +45,6 @@ PR.md
 
 # tutorial/gams
 tutorial/transport/*.gdx
-tutorial/transport/*.lst
\ No newline at end of file
+tutorial/transport/*.lst
+tests/tutorial/transport/*.gdx
+tests/tutorial/transport/*.lst
\ No newline at end of file
diff --git a/tests/tutorial/transport/dantzig_model_gams.py b/tests/tutorial/transport/dantzig_model_gams.py
new file mode 100644
index 00000000..47415c6a
--- /dev/null
+++ b/tests/tutorial/transport/dantzig_model_gams.py
@@ -0,0 +1,84 @@
+import copy
+from pathlib import Path
+
+import gams.transfer as gt
+import pandas as pd
+from gams import GamsWorkspace
+
+from ixmp4.core import Run
+
+
+def write_run_to_gams(run: Run) -> gt.Container:
+    """Writes scenario data from the Run to a GAMS container."""
+    m = gt.Container()
+    indexsets = [
+        gt.Set(
+            container=m,
+            name=indexset.name,
+            records=indexset.elements,
+            description=indexset.docs
+            if indexset.docs
+            else "",  # description is "optional", but must be str
+        )
+        for indexset in run.optimization.indexsets.list()
+    ]
+
+    for scalar in run.optimization.scalars.list():
+        gt.Parameter(
+            container=m,
+            name=scalar.name,
+            records=scalar.value,
+            description=scalar.docs if scalar.docs else "",
+        )
+
+    for parameter in run.optimization.parameters.list():
+        domains = [
+            indexset
+            for indexset in indexsets
+            if indexset.name in parameter.constrained_to_indexsets
+        ]
+        records = copy.deepcopy(parameter.data)
+        del records[
+            "units"
+        ]  # all parameters must have units, but GAMS doesn't work on them
+        gt.Parameter(
+            container=m,
+            name=parameter.name,
+            domain=domains,
+            records=records,
+            description=parameter.docs if parameter.docs else "",
+        )
+
+    return m
+
+
+# TODO I'd like to have proper Paths here, but gams can only handle str
+def solve(model_file: Path, data_file: Path, result_file: Path | None = None) -> None:
+    ws = GamsWorkspace(working_directory=Path(__file__).parent.absolute())
+    ws.add_database_from_gdx(gdx_file_name=str(data_file))
+    gams_options = ws.add_options()
+    gams_options.defines["in"] = str(data_file)
+    if result_file:
+        gams_options.defines["out"] = str(result_file)
+    job = ws.add_job_from_file(file_name=str(model_file))
+    job.run(gams_options=gams_options)
+
+
+def read_solution_to_run(run: Run, result_file: Path) -> None:
+    m = gt.Container(load_from=result_file)
+    for variable in run.optimization.variables.list():
+        # DF also includes lower, upper, scale
+        variable_data: pd.DataFrame = (
+            m.data[variable.name]
+            .records[["level", "marginal"]]
+            .rename(columns={"level": "levels", "marginal": "marginals"})
+        )
+        run.optimization.variables.get(variable.name).add(data=variable_data)
+
+    for equation in run.optimization.equations.list():
+        equation_data: pd.DataFrame = (
+            m.data[equation.name]
+            .records[["level", "marginal"]]
+            .rename(columns={"level": "levels", "marginal": "marginals"})
+        )
+        run.optimization.equations.get(equation.name).add(data=equation_data)
diff --git a/tests/tutorial/transport/test_dantzig_model_gams.py b/tests/tutorial/transport/test_dantzig_model_gams.py
new file mode 100644
index 00000000..50df0160
--- /dev/null
+++ b/tests/tutorial/transport/test_dantzig_model_gams.py
@@ -0,0 +1,131 @@
+import copy
+import shutil
+from pathlib import Path
+
+from gams.transfer import Set
+
+from ixmp4 import Platform
+
+from ...utils import all_platforms, create_dantzig_run
+from .dantzig_model_gams import read_solution_to_run, solve, write_run_to_gams
+
+
+@all_platforms
+class TestTransportTutorialLinopy:
+    # NOTE The function could be expanded by tables, equations, variables, none of which
+    #  are tested here.
+    def test_write_run_to_gams(self, test_mp, request):
+        test_mp: Platform = request.getfixturevalue(test_mp)  # type: ignore
+        run = create_dantzig_run(test_mp)
+        gams_container = write_run_to_gams(run=run)
+
+        # Should include exactly i, j, f, a, b, d
+        assert len(gams_container.data.keys()) == 6
+
+        gams_indexsets: list[Set] = []  # type: ignore
+        for indexset in run.optimization.indexsets.list():
+            gams_indexset = gams_container.data[indexset.name]
+            assert gams_indexset.name == indexset.name
+            assert gams_indexset.records["uni"].to_list() == indexset.elements
+            gams_indexsets.append(gams_indexset)
+
+        for scalar in run.optimization.scalars.list():
+            gams_scalar = gams_container.data[scalar.name]
+            assert gams_scalar.name == scalar.name
+            # Should only have one value
+            assert len(gams_scalar.records["value"]) == 1
+            assert gams_scalar.records["value"].values[0] == scalar.value
+
+        for parameter in run.optimization.parameters.list():
+            gams_parameter = gams_container.data[parameter.name]
+            assert gams_parameter.name == parameter.name
+
+            expected_domains = [
+                indexset
+                for indexset in gams_indexsets
+                if indexset.name in parameter.constrained_to_indexsets
+            ]
+            assert gams_parameter.domain == expected_domains
+
+            expected_records = copy.deepcopy(parameter.data)
+            del expected_records[
+                "units"
+            ]  # all parameters must have units, but GAMS doesn't work on them
+            assert (
+                gams_parameter.records.rename(columns={"value": "values"}).to_dict(
+                    orient="list"
+                )
+                == expected_records
+            )
+
+    def test_solve(self, test_mp, request, tmp_path):
+        test_mp: Platform = request.getfixturevalue(test_mp)  # type: ignore
+        run = create_dantzig_run(test_mp)
+        gams_container = write_run_to_gams(run=run)
+        data_file = tmp_path / "transport_data.gdx"
+        # TODO once we know where the tests land, figure out how to navigate paths
+        # Same below.
+        model_file = shutil.copy(
+            src=Path(__file__).parent.absolute() / "transport_ixmp4.gms",
+            dst=tmp_path / "transport_ixmp4.gms",
+        )
+        gams_container.write(write_to=data_file)
+
+        # Test writing to default location
+        solve(model_file=model_file, data_file=data_file)
+        default_result_file = Path(__file__).parent.absolute() / "transport_results.gdx"
+        assert default_result_file.is_file()
+
+        # Test writing to specified location
+        result_file: Path = tmp_path / "different_transport_results.gdx"  # type: ignore
+        solve(model_file=model_file, data_file=data_file, result_file=result_file)
+        assert result_file.is_file()
+
+    # TODO Maybe this test could be made more performant by receiving a run where the
+    # scenario has already been solved. However, this would make the test scenario less
+    # isolated. Also, solving the dantzig model only takes a few seconds.
+    def test_read_solution_to_run(self, test_mp, request, tmp_path):
+        test_mp: Platform = request.getfixturevalue(test_mp)  # type: ignore
+        run = create_dantzig_run(test_mp)
+        gams_container = write_run_to_gams(run=run)
+        data_file = tmp_path / "transport_data.gdx"
+        model_file = shutil.copy(
+            src=Path(__file__).parent.absolute() / "transport_ixmp4.gms",
+            dst=tmp_path / "transport_ixmp4.gms",
+        )
+        gams_container.write(write_to=data_file)
+        solve(model_file=model_file, data_file=data_file)
+        read_solution_to_run(
+            run=run,
+            result_file=Path(__file__).parent.absolute() / "transport_results.gdx",
+        )
+
+        # Test objective value
+        z = run.optimization.variables.get("z")
+        assert z.levels == [153.675]
+        assert z.marginals == [0.0]
+
+        # Test shipment quantities
+        assert run.optimization.variables.get("x").data == {
+            "levels": [50.0, 300.0, 0.0, 275.0, 0.0, 275.0],
+            "marginals": [
+                0.0,
+                0.0,
+                0.036000000000000004,
+                0.0,
+                0.009000000000000008,
+                0.0,
+            ],
+        }
+
+        # Test demand equation
+        assert run.optimization.equations.get("demand").data == {
+            "levels": [325.0, 300.0, 275.0],
+            "marginals": [0.225, 0.153, 0.126],
+        }
+
+        # Test supply equation
+        assert run.optimization.equations.get("supply").data == {
+            "levels": [350.0, 550.0],
+            "marginals": [-0.0, 0.0],
+        }
diff --git a/tests/tutorial/transport/test_dantzig_model_linopy.py b/tests/tutorial/transport/test_dantzig_model_linopy.py
index e3755bc8..09d386d8 100644
--- a/tests/tutorial/transport/test_dantzig_model_linopy.py
+++ b/tests/tutorial/transport/test_dantzig_model_linopy.py
@@ -3,75 +3,14 @@
 import xarray as xr
 
 from ixmp4 import Platform
-from ixmp4.core import Run, Unit
 
-from ...utils import all_platforms
+from ...utils import all_platforms, create_dantzig_run
 from .dantzig_model_linopy import (
     create_dantzig_model,
     read_dantzig_solution,
 )
 
 
-def create_dantzig_run(mp: Platform) -> Run:
-    """Create a Run for the transport tutorial.
-
-    Please see the tutorial file for explanation.
-    """
-    # Only needed once for each mp
-    try:
-        cases = mp.units.get("cases")
-        km = mp.units.get("km")
-        unit_cost_per_case = mp.units.get("USD/km")
-    except Unit.NotFound:
-        cases = mp.units.create("cases")
-        km = mp.units.create("km")
-        unit_cost_per_case = mp.units.create("USD/km")
-
-    # Create run and all data sets
-    run = mp.runs.create(model="transport problem", scenario="standard")
-    a_data = {
-        "i": ["seattle", "san-diego"],
-        "values": [350, 600],
-        "units": [cases.name, cases.name],
-    }
-    b_data = pd.DataFrame(
-        [
-            ["new-york", 325, cases.name],
-            ["chicago", 300, cases.name],
-            ["topeka", 275, cases.name],
-        ],
-        columns=["j", "values", "units"],
-    )
-    d_data = {
-        "i": ["seattle", "seattle", "seattle", "san-diego", "san-diego", "san-diego"],
-        "j": ["new-york", "chicago", "topeka", "new-york", "chicago", "topeka"],
-        "values": [2.5, 1.7, 1.8, 2.5, 1.8, 1.4],
-        "units": [km.name] * 6,
-    }
-
-    # Add all data to the run
-    run.optimization.indexsets.create("i").add(["seattle", "san-diego"])
-    run.optimization.indexsets.create("j").add(["new-york", "chicago", "topeka"])
-    run.optimization.parameters.create(name="a", constrained_to_indexsets=["i"]).add(
-        data=a_data
-    )
-    run.optimization.parameters.create("b", constrained_to_indexsets=["j"]).add(
-        data=b_data
-    )
-    run.optimization.parameters.create("d", constrained_to_indexsets=["i", "j"]).add(
-        data=d_data
-    )
-    run.optimization.scalars.create(name="f", value=90, unit=unit_cost_per_case)
-
-    # Create further optimization items to store solution data
-    run.optimization.variables.create("z")
-    run.optimization.variables.create("x", constrained_to_indexsets=["i", "j"])
-    run.optimization.equations.create("supply", constrained_to_indexsets=["i"])
-    run.optimization.equations.create("demand", constrained_to_indexsets=["j"])
-
-    return run
-
-
 @all_platforms
 class TestTransportTutorialLinopy:
     def test_create_dantzig_model(self, test_mp, request):
diff --git a/tests/tutorial/transport/transport_ixmp4.gms b/tests/tutorial/transport/transport_ixmp4.gms
new file mode 100644
index 00000000..861adb21
--- /dev/null
+++ b/tests/tutorial/transport/transport_ixmp4.gms
@@ -0,0 +1,70 @@
+*Basic example of transport model from GAMS model library
+
+$Title  A Transportation Problem (TRNSPORT,SEQ=1)
+$Ontext
+
+This problem finds a least cost shipping schedule that meets
+requirements at markets and supplies at factories.
+
+Dantzig, G B, Chapter 3.3. In Linear Programming and Extensions.
+Princeton University Press, Princeton, New Jersey, 1963.
+
+This formulation is described in detail in:
+Rosenthal, R E, Chapter 2: A GAMS Tutorial. In GAMS: A User's Guide.
+The Scientific Press, Redwood City, California, 1988.
+
+The line numbers will not match those in the book because of these
+comments.
+
+$Offtext
+
+Sets
+    i   canning plants
+    j   markets
+;
+
+Parameters
+    a(i)   capacity of plant i in cases
+    b(j)   demand at market j in cases
+    d(i,j) distance in thousands of miles
+    f      freight in dollars per case per thousand miles
+;
+
+* This file will read in data from a gdx file and write results to another gdx file.
+* The name of the input and output file can either be set directly from the command line,
+* e.g. `gams transport_ixmp.gms --in\=\"<name>\" --out\=\"<name>\"`.
+* If no command line parameters are given, the input and output files are set as specific below.
+
+$IF NOT set in  $SETGLOBAL in  'transport_data.gdx'
+$IF NOT set out $SETGLOBAL out 'transport_results.gdx'
+
+$GDXIN '%in%'
+$LOAD i, j, a, b, d, f
+$GDXIN
+
+Parameter c(i,j)  transport cost in thousands of dollars per case ;
+          c(i,j) = f * d(i,j) / 1000 ;
+Variables
+     x(i,j)  shipment quantities in cases
+     z       total transportation costs in thousands of dollars ;
+
+Positive Variable x ;
+
+Equations
+     cost        define objective function
+     supply(i)   observe supply limit at plant i
+     demand(j)   satisfy demand at market j ;
+
+cost ..        z  =e=  sum((i,j), c(i,j)*x(i,j)) ;
+
+supply(i) ..   sum(j, x(i,j))  =l=  a(i) ;
+
+demand(j) ..   sum(i, x(i,j))  =g=  b(j) ;
+
+Model transport /all/ ;
+
+Solve transport using lp minimizing z ;
+
+Display x.l, x.m ;
+
+Execute_unload '%out%';
\ No newline at end of file
diff --git a/tests/utils.py b/tests/utils.py
index 1dc74f8c..ecf73dd4 100644
--- a/tests/utils.py
+++ b/tests/utils.py
@@ -3,6 +3,7 @@
 import pytest
 
 from ixmp4 import DataPoint
+from ixmp4.core import Platform, Run, Unit
 
 from .conftest import SKIP_PGSQL_TESTS
 
@@ -198,3 +199,63 @@ def create_iamc_query_test_data(test_mp):
     run2.meta = {"run": 2, "test": "string", "bool": False}
 
     return [r1, r2, r3], units
+
+
+def create_dantzig_run(mp: Platform) -> Run:
+    """Create a Run for the transport tutorial.
+
+    Please see the tutorial file for explanation.
+    """
+    # Only needed once for each mp
+    try:
+        cases = mp.units.get("cases")
+        km = mp.units.get("km")
+        unit_cost_per_case = mp.units.get("USD/km")
+    except Unit.NotFound:
+        cases = mp.units.create("cases")
+        km = mp.units.create("km")
+        unit_cost_per_case = mp.units.create("USD/km")
+
+    # Create run and all data sets
+    run = mp.runs.create(model="transport problem", scenario="standard")
+    a_data = {
+        "i": ["seattle", "san-diego"],
+        "values": [350, 600],
+        "units": [cases.name, cases.name],
+    }
+    b_data = pd.DataFrame(
+        [
+            ["new-york", 325, cases.name],
+            ["chicago", 300, cases.name],
+            ["topeka", 275, cases.name],
+        ],
+        columns=["j", "values", "units"],
+    )
+    d_data = {
+        "i": ["seattle", "seattle", "seattle", "san-diego", "san-diego", "san-diego"],
+        "j": ["new-york", "chicago", "topeka", "new-york", "chicago", "topeka"],
+        "values": [2.5, 1.7, 1.8, 2.5, 1.8, 1.4],
+        "units": [km.name] * 6,
+    }
+
+    # Add all data to the run
+    run.optimization.indexsets.create("i").add(["seattle", "san-diego"])
+    run.optimization.indexsets.create("j").add(["new-york", "chicago", "topeka"])
+    run.optimization.parameters.create(name="a", constrained_to_indexsets=["i"]).add(
+        data=a_data
+    )
+    run.optimization.parameters.create("b", constrained_to_indexsets=["j"]).add(
+        data=b_data
+    )
+    run.optimization.parameters.create("d", constrained_to_indexsets=["i", "j"]).add(
+        data=d_data
+    )
+    run.optimization.scalars.create(name="f", value=90, unit=unit_cost_per_case)
+
+    # Create further optimization items to store solution data
+    run.optimization.variables.create("z")
+    run.optimization.variables.create("x", constrained_to_indexsets=["i", "j"])
+    run.optimization.equations.create("supply", constrained_to_indexsets=["i"])
+    run.optimization.equations.create("demand", constrained_to_indexsets=["j"])
+
+    return run
diff --git a/tutorial/transport/dantzig_model_gams.py b/tutorial/transport/dantzig_model_gams.py
index 00eda507..47415c6a 100644
--- a/tutorial/transport/dantzig_model_gams.py
+++ b/tutorial/transport/dantzig_model_gams.py
@@ -8,9 +8,8 @@
 from ixmp4.core import Run
 
 
-def write_run_to_gdx(run: Run, write_to: Path | str = "transport_data.gdx") -> None:
-    """Writes scenario data from the Run to a GDX file."""
-    write_to = Path(write_to).absolute()
+def write_run_to_gams(run: Run) -> gt.Container:
+    """Writes scenario data from the Run to a GAMS container."""
     m = gt.Container()
     indexsets = [
         gt.Set(
@@ -50,19 +49,23 @@ def write_run_to_gdx(run: Run, write_to: Path | str = "transport_data.gdx") -> N
             description=parameter.docs if parameter.docs else "",
         )
 
-    m.write(write_to=write_to)
+    return m
 
 
-def solve_dantzig_model(data_file: Path | str = "transport_data.gdx") -> None:
-    ws = GamsWorkspace(working_directory=Path().absolute())
-    # Data from previous step might be saved elsewhere
-    ws.add_database_from_gdx(data_file)
-    job = ws.add_job_from_file("transport_ixmp4.gms")
-    job.run()
+# TODO I'd like to have proper Paths here, but gams can only handle str
+def solve(model_file: Path, data_file: Path, result_file: Path | None = None) -> None:
+    ws = GamsWorkspace(working_directory=Path(__file__).parent.absolute())
+    ws.add_database_from_gdx(gdx_file_name=str(data_file))
+    gams_options = ws.add_options()
+    gams_options.defines["in"] = str(data_file)
+    if result_file:
+        gams_options.defines["out"] = str(result_file)
+    job = ws.add_job_from_file(file_name=str(model_file))
+    job.run(gams_options=gams_options)
 
 
-def read_dantzig_solution(run: Run) -> None:
-    m = gt.Container("transport_results.gdx")
+def read_solution_to_run(run: Run, result_file: Path) -> None:
+    m = gt.Container(load_from=result_file)
     for variable in run.optimization.variables.list():
         # DF also includes lower, upper, scale
         variable_data: pd.DataFrame = (
diff --git a/tutorial/transport/gams_transport.ipynb b/tutorial/transport/gams_transport.ipynb
index fe075a2e..1bc1fe7b 100644
--- a/tutorial/transport/gams_transport.ipynb
+++ b/tutorial/transport/gams_transport.ipynb
@@ -376,13 +376,16 @@
     "### Solve the scenario\n",
     "\n",
     "In this tutorial, we solve the scenario using GAMS. \n",
-    "The solution process requires three steps:\n",
+    "GAMS requires model data and structure to follow their proprietary formats, so we need to write out model data to this format before GAMS can handle it (\"`data_file`\").\n",
+    "For this tutorial, the structure of the model (the equations that relate the variables) has already been written to a `transport_ixmp4.gms`, so we can simply use that (\"`model_file`\").\n",
+    "After the calculation, GAMS writes the results out in the `.gdx` format again, which we need to read back into our `Run` (\"`result_file`\").\n",
+    "Thus, the solution process requires three steps:\n",
     "\n",
-    "1. Write the `run` data to a `.gdx` file that GAMS can read.\n",
+    "1. Write the `Run` data to a `.gdx` file that GAMS can read.\n",
     "2. Solve the scenario in GAMS.\n",
-    "3. Read the solution data back into our `run`.\n",
+    "3. Read the solution data back into our `Run`.\n",
     "\n",
-    "For convenience, these are pre-existing functions here. Feel free to check out their implementation in `dantzig_model_gams.py`, though, if you are interested!"
+    "For convenience, these are pre-existing functions here, so we can immediately use them. Feel free to check out their implementation in `dantzig_model_gams.py`, though, if you are interested!"
    ]
   },
   {
@@ -391,23 +394,31 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "from pathlib import Path\n",
+    "\n",
     "from tutorial.transport.dantzig_model_gams import (\n",
-    "    read_dantzig_solution,\n",
-    "    solve_dantzig_model,\n",
-    "    write_run_to_gdx,\n",
+    "    read_solution_to_run,\n",
+    "    solve,\n",
+    "    write_run_to_gams,\n",
     ")\n",
     "\n",
-    "# Per default, this will use the same directory this file is in.\n",
-    "# If you want to change that, provide a path to a `.gdx` file as `write_to`.\n",
-    "write_run_to_gdx(run=run)\n",
+    "# First, we collect all relevant Run data in a gams container:\n",
+    "gams_container = write_run_to_gams(run=run)\n",
+    "\n",
+    "# For the next steps, you can adapt the locations of the files as discussed above.\n",
+    "# If you are not sure how, please use these defaults:\n",
+    "data_file = Path(\"transport_data.gdx\").absolute()\n",
+    "model_file = Path(\"transport_ixmp4.gms\").absolute()\n",
+    "result_file = Path(\"transport_results.gdx\").absolute()\n",
+    "\n",
+    "# The gams container then allows us to write the required data format file:\n",
+    "gams_container.write(write_to=data_file)\n",
     "\n",
-    "# If you changed `write_to` from its default value before, provide the same argument\n",
-    "# here as `data_file`. The model solution will be written to the directory containing\n",
-    "# the model file in any case, though.\n",
-    "solve_dantzig_model()\n",
+    "# Now, we can actually solve the model:\n",
+    "solve(model_file=model_file, data_file=data_file, result_file=result_file)\n",
     "\n",
-    "# The model solution can thus always be read from the same location.\n",
-    "read_dantzig_solution(run=run)"
+    "# Finally, we read the solution back into our Run:\n",
+    "read_solution_to_run(run=run, result_file=result_file)"
    ]
   },
   {