diff --git a/src/Algorithms/ASB07/ASB07.jl b/src/Algorithms/ASB07/ASB07.jl
index ad90e736dd..2970f51ec1 100644
--- a/src/Algorithms/ASB07/ASB07.jl
+++ b/src/Algorithms/ASB07/ASB07.jl
@@ -1,6 +1,6 @@
 @with_kw struct ASB07{N, AM} <: AbstractContinuousPost
     δ::N
-    approx_model::AM=CorrectionHullApproximation(order=10, exp_method="base")
+    approx_model::AM=CorrectionHull(order=10, exp_method=:base)
     max_order::Int=10
 end
 
diff --git a/src/Algorithms/BFFPSV18/BFFPSV18.jl b/src/Algorithms/BFFPSV18/BFFPSV18.jl
index e3266016a7..d64c145d41 100644
--- a/src/Algorithms/BFFPSV18/BFFPSV18.jl
+++ b/src/Algorithms/BFFPSV18/BFFPSV18.jl
@@ -17,7 +17,7 @@ TODO
 @with_kw struct BFFPSV18{N, ST, AM, IDX, BLK, RBLK, CBLK, PT} <: AbstractContinuousPost
     δ::N
     setrep::ST  # remove ?
-    approx_model::AM=ForwardApproximation(sih_method="concrete", exp_method="base", phi2_method="base")
+    approx_model::AM=ForwardApproximation(sih_method=:concrete, exp_method=:base, phi2_method=:base)
     vars::IDX
     block_indices::BLK
     row_blocks::RBLK
diff --git a/src/Algorithms/GLGM06/GLGM06.jl b/src/Algorithms/GLGM06/GLGM06.jl
index 45dcb6cbb3..4d6d9ead03 100644
--- a/src/Algorithms/GLGM06/GLGM06.jl
+++ b/src/Algorithms/GLGM06/GLGM06.jl
@@ -1,5 +1,5 @@
 """
-    GLGM06{N, AM<:AbstractApproximationModel} <: AbstractContinuousPost
+    GLGM06{N, AM} <: AbstractContinuousPost
 
 Implementation of Girard - Le Guernic - Maler algorithm for reachability of
 uncertain linear systems using zonotopes.
@@ -22,10 +22,11 @@ The type fields are:
 
 See [xxx] and [yyy]
 """
-@with_kw struct GLGM06{N, AM<:AbstractApproximationModel} <: AbstractContinuousPost
+@with_kw struct GLGM06{N, AM} <: AbstractContinuousPost
     δ::N
     # nota: la opcion set_operations="zonotope" es ignorada (?)
-    approx_model::AM=ForwardApproximation(sih_method="concrete", exp_method="base", phi2_method="base", set_operations="zonotope")
+    approx_model::AM=Forward(sih_method=:concrete, exp_method=:base,
+                             phi2_method=:base, set_operations=:zonotope)
     max_order::Int=10
 end
 
diff --git a/src/Continuous/discretization.jl b/src/Continuous/discretization.jl
index f3ecf7da03..4770d0ff0d 100644
--- a/src/Continuous/discretization.jl
+++ b/src/Continuous/discretization.jl
@@ -1,41 +1,40 @@
+"""
+    abstract type AbstractApproximationModel end
+
+Abstract supertype for all approximation models.
+"""
 abstract type AbstractApproximationModel end
-const AAModel = AbstractApproximationModel
-
-export ForwardApproximation,
-       BackwardApproximation,
-       DiscreteApproximation,
-       CorrectionHullApproximation
-
-@with_kw struct ForwardApproximation <: AbstractApproximationModel
-    exp_method::String="base"
-    set_operations::String="lazy"
-    sih_method::String="concrete"
-    phi2_method::String="base"
+
+@with_kw struct Forward <: AbstractApproximationModel
+    exp_method::Symbol=:base
+    set_operations::Symbol=:lazy
+    sih_method::Symbol=:concrete
+    phi2_method::Symbol=:base
 end
 
-@with_kw struct BackwardApproximation <: AbstractApproximationModel
-    exp_method::String="base"
-    set_operations::String="lazy"
-    sih_method::String="concrete"
-    phi2_method::String="base"
+@with_kw struct Backward <: AbstractApproximationModel
+    exp_method::Symbol=:base
+    set_operations::Symbol=:lazy
+    sih_method::Symbol=:concrete
+    phi2_method::Symbol=:base
 end
 
 # no bloating
-struct DiscreteApproximation <: AbstractApproximationModel
+struct Discrete <: AbstractApproximationModel
 #
 end
 
-@with_kw struct CorrectionHullApproximation <: AbstractApproximationModel
+@with_kw struct CorrectionHull <: AbstractApproximationModel
    order::Int=10
-   exp_method::String="base"
+   exp_method::Symbol=:base
 end
 
 function _default_approximation_model(ivp::IVP{<:AbstractContinuousSystem})
-    return ForwardApproximation()
+    return Forward()
 end
 
 # homogeneous case
-function discretize(ivp::IVP{<:CLCS, <:LazySet}, δ::Float64, alg::ForwardApproximation)
+function discretize(ivp::IVP{<:CLCS, <:LazySet}, δ::Float64, alg::Forward)
     A = state_matrix(ivp)
     X0 = initial_state(ivp)
     ϕ = _exp(A, δ, alg.exp_method)
@@ -43,7 +42,7 @@ function discretize(ivp::IVP{<:CLCS, <:LazySet}, δ::Float64, alg::ForwardApprox
     Phi2A_abs = Φ₂(A_abs, δ, alg.phi2_method)
 
     # "forward" algorithm, uses E⁺
-    @assert alg.sih_method == "concrete"
+    @assert alg.sih_method == :concrete
     # TODO : specialize, add option to compute the concrete linear map
     Einit = symmetric_interval_hull(Phi2A_abs * symmetric_interval_hull((A * A) * X0))
 
@@ -54,7 +53,7 @@ function discretize(ivp::IVP{<:CLCS, <:LazySet}, δ::Float64, alg::ForwardApprox
 end
 
 # inhomogeneous case
-function discretize(ivp::IVP{<:CLCCS, <:LazySet}, δ::Float64, alg::ForwardApproximation)
+function discretize(ivp::IVP{<:CLCCS, <:LazySet}, δ::Float64, alg::Forward)
     A = state_matrix(ivp)
     X0 = initial_state(ivp)
     X = stateset(ivp)
@@ -63,7 +62,7 @@ function discretize(ivp::IVP{<:CLCCS, <:LazySet}, δ::Float64, alg::ForwardAppro
     A_abs = _elementwise_abs(A)
     Phi2A_abs = Φ₂(A_abs, δ, alg.phi2_method)
 
-    @assert alg.sih_method == "concrete"
+    @assert alg.sih_method == :concrete
     # TODO : specialize, add option to compute the concrete linear map
     Einit = symmetric_interval_hull(Phi2A_abs * symmetric_interval_hull((A * A) * X0))
 
@@ -91,7 +90,7 @@ _correction_hull(A::AbstractMatrix, t, p) = correction_hull(IntervalMatrix(A), t
 # implements: Ω0 = CH(X0, exp(A*δ) * X0) ⊕ F*X0
 # where F is the correction (interval) matrix
 # if A is an interval matix, the exponential is overapproximated
-function discretize(ivp::IVP{<:CLCS, <:LazySet}, δ::Float64, alg::CorrectionHullApproximation)
+function discretize(ivp::IVP{<:CLCS, <:LazySet}, δ::Float64, alg::CorrectionHull)
     A = state_matrix(ivp)
     X0 = initial_state(ivp)
     X = stateset(ivp)
diff --git a/src/Continuous/exponentiation.jl b/src/Continuous/exponentiation.jl
index c12976f208..c90822eab2 100644
--- a/src/Continuous/exponentiation.jl
+++ b/src/Continuous/exponentiation.jl
@@ -20,7 +20,7 @@ using LinearAlgebra: checksquare
 @inline _exp(A::IdentityMultiple) = IdentityMultiple(exp(A.M.λ), size(A, 1))
 
 """
-    _exp(A::AbstractMatrix, δ::Float64, method::String)
+    _exp(A::AbstractMatrix, δ::Float64, method::Symbol)
 
 Compute the matrix exponential ``e^{Aδ}``.
 
@@ -28,15 +28,15 @@ Compute the matrix exponential ``e^{Aδ}``.
 
 - `A`       -- matrix
 - `δ`       -- step size
-- `method`  -- the method used to take the matrix exponential of `A`;
+- `method`  -- symbol with the method used to take the matrix exponential of `A`;
                possible options are:
 
-    - `"base"` -- use the scaling and squaring method implemented in Julia standard
-                  library; see `?exp` for details
-    - `"lazy"` -- return a lazy wrapper type around the matrix exponential using
-                  the implementation `LazySets.SparseMatrixExp`
-    - `"pade"` -- apply Pade approximant method to compute the matrix exponential
-                  of a sparse matrix (requires `Expokit`)
+    - `:base` -- use the scaling and squaring method implemented in Julia standard
+                 library; see `?exp` for details
+    - `:lazy` -- return a lazy wrapper type around the matrix exponential using
+                 the implementation `LazySets.SparseMatrixExp`
+    - `:pade` -- apply Pade approximant method to compute the matrix exponential
+                 of a sparse matrix (requires `Expokit`)
 
 ### Output
 
@@ -48,15 +48,15 @@ If the algorithm `"lazy"` is used, evaluations of the action of the matrix
 exponential are done with the `expmv` implementation from `Expokit`
 (but see `LazySets#1312` for the planned generalization to other backends).
 """
-function _exp(A::AbstractMatrix, δ::Float64, method::String)
+function _exp(A::AbstractMatrix, δ::Float64, method::Symbol)
     n = checksquare(A)
-    if method == "base"
+    if method == :base
         return _exp(A * δ) # TODO use dots ? (requires MathematicalSystems#189 for IdentityMultiple)
 
-    elseif method == "lazy"
+    elseif method == :lazy
         return _exp_lazy(A * δ)
 
-    elseif method == "pade"
+    elseif method == :pade
         @requires Expokit
         return _exp_pade(A * δ)
 
@@ -122,15 +122,15 @@ function Φ₁(A::AbstractMatrix, δ::Float64, method::String)
     n = checksquare(A)
     B = _Aδ_3n(A, δ, n)
 
-    if method == "base"
+    if method == :base
         P = _exp(B)
         return P[1:n, (n+1):2*n]
 
-    elseif method == "lazy"
+    elseif method == :lazy
         P = _exp_lazy(B)
         return sparse(get_columns(P, (n+1):2*n)[1:n, :])
 
-    elseif method == "pade"
+    elseif method == :pade
         P = _exp_pade(B)
         return P[1:n, (n+1):2*n]
 
@@ -185,23 +185,23 @@ It can be shown that `Φ₂(A, δ) = P[1:n, (2*n+1):3*n]`.
 International Conference on Computer Aided Verification. Springer, Berlin,
 Heidelberg, 2011.
 """
-function Φ₂(A::AbstractMatrix, δ::Float64, method::String)
-    if method == "inverse"
+function Φ₂(A::AbstractMatrix, δ::Float64, method::Symbol)
+    if method == :inverse
         return _Φ₂_inverse(A, δ)
     end
 
     n = checksquare(A)
     B = _Aδ_3n(A, δ, n)
 
-    if method == "base"
+    if method == :base
         P = _exp(B)
         return P[1:n, (2*n+1):3*n]
 
-    elseif method == "lazy"
+    elseif method == :lazy
         P = _exp_lazy(B)
         return sparse(get_columns(P, (2*n+1):3*n)[1:n, :])
 
-    elseif method == "pade"
+    elseif method == :pade
         P = _exp_pade(B)
         return P[1:n, (2*n+1):3*n]
 
diff --git a/src/Initialization/exports.jl b/src/Initialization/exports.jl
index 5ffb0378e9..3ac3c28d32 100644
--- a/src/Initialization/exports.jl
+++ b/src/Initialization/exports.jl
@@ -8,6 +8,12 @@ export
     TMJets,
     ASB07,
 
+# Approximation models
+    Forward,
+    Backward,
+    Discrete,
+    CorrectionHull,
+
 # Flowpipes
     flowpipe,
     Flowpipe,