diff --git a/README.md b/README.md
index 549aa4ae..b30eea87 100644
--- a/README.md
+++ b/README.md
@@ -263,6 +263,14 @@ REISE.run_scenario(;
     interval=24, n_interval=3, start_index=1, outputfolder="output",
     inputfolder=pwd(), num_segments=3)
 ```
+If another solver is desired, it can be passed via the `optimizer_factory` argument, e.g.:
+```julia
+import GLPK
+REISE.run_scenario(;
+    interval=24, n_interval=3, start_index=1, outputfolder="output",
+    inputfolder=pwd(), optimizer_factory=GLPK.Optimizer)
+```
+Be sure to pass the factory itself (e.g. `GLPK.Optimizer`) rather than an instance (e.g. `GLPK.Optimizer()`). See the [JuMP.Model documentation] for more information.
 
 ## Usage (Python)
 
@@ -665,3 +673,4 @@ Penalty for ending the interval with less stored energy than the start, or rewar
 
 [Gurobi.jl]: https://github.com/JuliaOpt/Gurobi.jl#installation
 [Julia Package Manager]: https://julialang.github.io/Pkg.jl/v1/managing-packages/
+[JuMP.Model documentation]: https://jump.dev/JuMP.jl/stable/solvers/#JuMP.Model-Tuple{Any}
diff --git a/src/REISE.jl b/src/REISE.jl
index af931949..ac20dd42 100644
--- a/src/REISE.jl
+++ b/src/REISE.jl
@@ -29,15 +29,18 @@ Run a scenario consisting of several intervals.
 'n_interval' specifies the number of intervals in a scenario.
 'start_index' specifies the starting hour of the first interval, to determine
     which time-series data should be loaded into each intervals.
-'outputfolder' specifies where to store the results. This folder will be
-    created if it does not exist at runtime.
 'inputfolder' specifies where to load the relevant data from. Required files
     are 'case.mat', 'demand.csv', 'hydro.csv', 'solar.csv', and 'wind.csv'.
+'outputfolder' specifies where to store the results. Defaults to an `output`
+    subdirectory of inputfolder. This folder will be created if it does not exist at
+    runtime.
+'optimizer_factory' is the solver used for optimization. If not specified, Gurobi is
+    used by default.
 """
 function run_scenario(;
         num_segments::Int=1, interval::Int, n_interval::Int, start_index::Int,
         inputfolder::String, outputfolder::Union{String, Nothing}=nothing,
-        threads::Union{Int, Nothing}=nothing)
+        threads::Union{Int, Nothing}=nothing, optimizer_factory=nothing)
     # Setup things that build once
     # If outputfolder not given, by default assign it inside inputfolder
     isnothing(outputfolder) && (outputfolder = joinpath(inputfolder, "output"))
@@ -45,7 +48,12 @@ function run_scenario(;
     isdir(outputfolder) || mkdir(outputfolder)
     stdout_filepath = joinpath(outputfolder, "stdout.log")
     stderr_filepath = joinpath(outputfolder, "stderr.err")
-    env = Gurobi.Env()
+    if isnothing(optimizer_factory)
+        optimizer_factory = Gurobi.Env()
+        solver_kwargs = Dict("Method" => 2, "Crossover" => 0)
+    else
+        solver_kwargs = Dict()
+    end
     case = read_case(inputfolder)
     storage = read_storage(inputfolder)
     println("All scenario files loaded!")
@@ -56,7 +64,6 @@ function run_scenario(;
         "storage" => storage,
         "interval_length" => interval,
         )
-    solver_kwargs = Dict("Method" => 2, "Crossover" => 0)
     # If a number of threads is specified, add to solver settings dict
     isnothing(threads) || (solver_kwargs["Threads"] = threads)
     println("All preparation complete!")
@@ -64,11 +71,13 @@ function run_scenario(;
     println("Redirecting outputs, see stdout.log & stderr.err in outputfolder")
     redirect_stdout_stderr(stdout_filepath, stderr_filepath) do
         # Loop through intervals
-        interval_loop(env, model_kwargs, solver_kwargs, interval, n_interval,
-                      start_index, inputfolder, outputfolder)
+        interval_loop(optimizer_factory, model_kwargs, solver_kwargs, interval,
+                      n_interval, start_index, inputfolder, outputfolder)
         GC.gc()
-        Gurobi.finalize(env)
-        println("Connection closed successfully!")
+        if isa(optimizer_factory, Gurobi.Env)
+            Gurobi.finalize(optimizer_factory)
+            println("Connection closed successfully!")
+        end
     end
 end
 
diff --git a/src/loop.jl b/src/loop.jl
index 5e839430..6da27e8d 100644
--- a/src/loop.jl
+++ b/src/loop.jl
@@ -1,9 +1,25 @@
+"""Convert a dict with string keys to a NamedTuple, for python-eqsue kwargs splatting"""
+function symbolize(d::Dict{String,Any})::NamedTuple
+    return (; (Symbol(k) => v for (k,v) in d)...)
+end
+
+
+function new_model(factory_like)::Union{JuMP.Model, JuMP.MOI.AbstractOptimizer}
+    if isa(factory_like, Gurobi.Env)
+        return JuMP.direct_model(Gurobi.Optimizer(factory_like))
+    else
+        return JuMP.Model(factory_like)
+    end
+end
+
+
 """
-    interval_loop(env, model_kwargs, solver_kwargs, interval, n_interval,
+    interval_loop(factory_like, model_kwargs, solver_kwargs, interval, n_interval,
                   start_index, inputfolder, outputfolder)
 
 Given:
-- a Gurobi environment `env`
+- optimizer instantiation object `factory_like`:
+    either a Gurobi environment or something that can be passed to JuMP.Model
 - a dictionary of model keyword arguments `model_kwargs`
 - a dictionary of solver keyword arguments `solver_kwargs`
 - an interval length `interval` (hours)
@@ -15,7 +31,7 @@ Given:
 Build a model, and run through the intervals, re-building the model and/or
 re-setting constraint right-hand-side values as necessary.
 """
-function interval_loop(env::Gurobi.Env, model_kwargs::Dict,
+function interval_loop(factory_like, model_kwargs::Dict,
                        solver_kwargs::Dict, interval::Int,
                        n_interval::Int, start_index::Int,
                        inputfolder::String, outputfolder::String)
@@ -45,10 +61,9 @@ function interval_loop(env::Gurobi.Env, model_kwargs::Dict,
             if storage_enabled
                 model_kwargs["storage_e0"] = storage.sd_table.InitialStorage
             end
-            m_kwargs = (; (Symbol(k) => v for (k,v) in model_kwargs)...)
-            m = JuMP.direct_model(Gurobi.Optimizer(env))
+            m = new_model(factory_like)
             JuMP.set_optimizer_attributes(m, pairs(solver_kwargs)...)
-            m, voi = _build_model(m; m_kwargs...)
+            m, voi = _build_model(m; symbolize(model_kwargs)...)
         elseif i == 2
             # Build a model with an initial ramp constraint
             model_kwargs["initial_ramp_enabled"] = true
@@ -56,10 +71,9 @@ function interval_loop(env::Gurobi.Env, model_kwargs::Dict,
             if storage_enabled
                 model_kwargs["storage_e0"] = storage_e0
             end
-            m_kwargs = (; (Symbol(k) => v for (k,v) in model_kwargs)...)
-            m = JuMP.direct_model(Gurobi.Optimizer(env))
+            m = new_model(factory_like)
             JuMP.set_optimizer_attributes(m, pairs(solver_kwargs)...)
-            m, voi = _build_model(m; m_kwargs...)
+            m, voi = _build_model(m; symbolize(model_kwargs)...)
         else
             # Reassign right-hand-side of constraints to match profiles
             bus_demand = _make_bus_demand(case, interval_start, interval_end)
@@ -125,8 +139,9 @@ function interval_loop(env::Gurobi.Env, model_kwargs::Dict,
                 results = get_results(f, voi, model_kwargs["case"])
                 break
             elseif ((status in numeric_statuses)
+                    & isa(factory_like, Gurobi.Env)
                     & !("BarHomogeneous" in keys(solver_kwargs)))
-                # if BarHomogeneous is not enabled, enable it and re-build
+                # if Gurobi, and BarHomogeneous is not enabled, enable it and re-solve
                 solver_kwargs["BarHomogeneous"] = 1
                 println("enable BarHomogeneous")
                 JuMP.set_optimizer_attribute(m, "BarHomogeneous", 1)
@@ -134,24 +149,23 @@ function interval_loop(env::Gurobi.Env, model_kwargs::Dict,
                     & !("load_shed_enabled" in keys(model_kwargs)))
                 # if load shed not enabled, enable it and re-build the model
                 model_kwargs["load_shed_enabled"] = true
-                m_kwargs = (; (Symbol(k) => v for (k,v) in model_kwargs)...)
                 println("rebuild with load shed")
-                m = JuMP.direct_model(Gurobi.Optimizer(env))
+                m = new_model(factory_like)
                 JuMP.set_optimizer_attributes(m, pairs(solver_kwargs)...)
-                m, voi = _build_model(m; m_kwargs...)
+                m, voi = _build_model(m; symbolize(model_kwargs)...)
                 intervals_without_loadshed = 0
-            elseif !("BarHomogeneous" in keys(solver_kwargs))
-                # if BarHomogeneous is not enabled, enable it and re-build
+            elseif (isa(factory_like, Gurobi.Env)
+                    & !("BarHomogeneous" in keys(solver_kwargs)))
+                # if Gurobi, and BarHomogeneous is not enabled, enable it and re-solve
                 solver_kwargs["BarHomogeneous"] = 1
                 println("enable BarHomogeneous")
                 JuMP.set_optimizer_attribute(m, "BarHomogeneous", 1)
             elseif !("load_shed_enabled" in keys(model_kwargs))
                 model_kwargs["load_shed_enabled"] = true
-                m_kwargs = (; (Symbol(k) => v for (k,v) in model_kwargs)...)
                 println("rebuild with load shed")
-                m = JuMP.direct_model(Gurobi.Optimizer(env))
+                m = new_model(factory_like)
                 JuMP.set_optimizer_attributes(m, pairs(solver_kwargs)...)
-                m, voi = _build_model(m; m_kwargs...)
+                m, voi = _build_model(m; symbolize(model_kwargs)...)
                 intervals_without_loadshed = 0
             else
                 # Something has gone very wrong
@@ -199,10 +213,9 @@ function interval_loop(env::Gurobi.Env, model_kwargs::Dict,
                 # delete! will work here even if the key is not present
                 delete!(solver_kwargs, "BarHomogeneous")
                 delete!(model_kwargs, "load_shed_enabled")
-                m_kwargs = (; (Symbol(k) => v for (k,v) in model_kwargs)...)
-                m = JuMP.direct_model(Gurobi.Optimizer(env))
+                m = new_model(factory_like)
                 JuMP.set_optimizer_attributes(m, pairs(solver_kwargs)...)
-                m, voi = _build_model(m; m_kwargs...)
+                m, voi = _build_model(m; symbolize(model_kwargs)...)
             end
         end
     end