diff --git a/Project.toml b/Project.toml
index 245837b3..04f05c77 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 MultivariatePolynomials = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
+Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
@@ -18,5 +19,6 @@ LinearAlgebra = "<0.0.1, 1.4"
 MacroTools = "0.5"
 MultivariatePolynomials = "0.3, 0.4, 0.5"
 RecipesBase = "0.6, 0.7, 0.8, 1.0, 1.1, 1.2"
+Requires = "1.3.1"
 SparseArrays = "<0.0.1, 1.4"
 julia = "1.5"
diff --git a/docs/Project.toml b/docs/Project.toml
index c25c14b6..5330a05c 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -6,4 +6,4 @@ LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043"
 [compat]
 Documenter = "1"
 DocumenterCitations = "1.3"
-LazySets = "2, 3, 4"
+LazySets = "5"
diff --git a/docs/src/lib/methods.md b/docs/src/lib/methods.md
index 109d87b2..3c4c9a96 100644
--- a/docs/src/lib/methods.md
+++ b/docs/src/lib/methods.md
@@ -55,6 +55,7 @@ isblackbox(::AbstractSystem)
 isnoisy(::AbstractSystem)
 iscontrolled(::AbstractSystem)
 isconstrained(::AbstractSystem)
+isparametric(::AbstractSystem)
 state_matrix(::AbstractSystem)
 input_matrix(::AbstractSystem)
 noise_matrix(::AbstractSystem)
diff --git a/docs/src/lib/types.md b/docs/src/lib/types.md
index e7784b16..8ae0c445 100644
--- a/docs/src/lib/types.md
+++ b/docs/src/lib/types.md
@@ -57,6 +57,8 @@ SecondOrderConstrainedAffineControlContinuousSystem
 SecondOrderConstrainedLinearControlContinuousSystem
 SecondOrderContinuousSystem
 SecondOrderConstrainedContinuousSystem
+LinearParametricContinuousSystem
+LinearControlParametricContinuousSystem
 ```
 
 ## Discrete Systems
@@ -94,6 +96,8 @@ SecondOrderConstrainedAffineControlDiscreteSystem
 SecondOrderConstrainedLinearControlDiscreteSystem
 SecondOrderDiscreteSystem
 SecondOrderConstrainedDiscreteSystem
+LinearParametricDiscreteSystem
+LinearControlParametricDiscreteSystem
 ```
 
 #### Discretization Algorithms
diff --git a/docs/src/man/systems.md b/docs/src/man/systems.md
index 9519dcbf..852e5866 100644
--- a/docs/src/man/systems.md
+++ b/docs/src/man/systems.md
@@ -99,6 +99,8 @@ However in this table we only included continuous system types for brevity.
 |SOACS|[`SecondOrderAffineContinuousSystem`](@ref)|
 |SOCACCS|[`SecondOrderConstrainedAffineControlContinuousSystem`](@ref)|
 |SOCLCCS|[`SecondOrderConstrainedLinearControlContinuousSystem`](@ref)|
+|LPCS|[`LinearParametricContinuousSystem`](@ref)|
+|LCPCS|[`LinearControlParametricContinuousSystem`](@ref)|
 
 The following table summarizes the equation represented by each system type
 (the names are given in abbreviated form). Again, discrete systems are not included. The column *Input constraints* is `yes` if the structure can model input or noise constraints (or both).
@@ -137,3 +139,5 @@ The following table summarizes the equation represented by each system type
 |``Mx'' + Cx' + Kx = Bu, x ∈ X, u ∈ U``|yes|yes|SOCLCCS|
 |``Mx'' + Cx' + f_i(x) = f_e``|no|no|SOCS|
 |``Mx'' + Cx' + f_i(x) = f_e``, x ∈ X, u ∈ U``|yes|yes|SOCCS|
+|``x' = A(θ)x, θ ∈ Θ``| no|no|LPCS|
+|``x' = A(θ)x + B(θ)u, θ ∈ Θ``| no|no|LPCS|
diff --git a/src/MathematicalSystems.jl b/src/MathematicalSystems.jl
index f620044f..cee8b36e 100644
--- a/src/MathematicalSystems.jl
+++ b/src/MathematicalSystems.jl
@@ -1,5 +1,6 @@
 module MathematicalSystems
 
+using Requires
 using LinearAlgebra, SparseArrays
 using LinearAlgebra: checksquare
 import RecipesBase
@@ -55,7 +56,8 @@ export islinear,
        isblackbox,
        isnoisy,
        iscontrolled,
-       isconstrained
+       isconstrained,
+       isparametric
 
 #====================================
 Concrete Types for Continuous Systems
diff --git a/src/abstract.jl b/src/abstract.jl
index a562dba1..99c04fd0 100644
--- a/src/abstract.jl
+++ b/src/abstract.jl
@@ -213,6 +213,18 @@ can also be applied to `typeof(s)`.
 """
 function isconstrained(::AbstractSystem) end
 
+"""
+    isparametric(s::AbstractSystem)
+
+Determines if the system `s` contains one or more symbolic parameters in its
+dynamics (for example, matrix entries that can take values from a
+specified set).
+
+The result of this function only depends on the
+system type, not the value, and can also be applied to `typeof(s)`.
+"""
+function isparametric(::AbstractSystem) end
+
 """
     AbstractMap
 
diff --git a/src/macros.jl b/src/macros.jl
index 4be98940..25bf5c5c 100644
--- a/src/macros.jl
+++ b/src/macros.jl
@@ -127,7 +127,7 @@ Return the system type whose field names match those in `fields`.
 ### Input
 
 - `AT`     -- abstract system type, which can be either `AbstractContinuousSystem`
-              or `AbstractDiscreSystem`
+              or `AbstractDiscreteSystem`
 - `fields` -- tuple of field names
 
 ### Output
diff --git a/src/systems.jl b/src/systems.jl
index 6618894f..60457c1a 100644
--- a/src/systems.jl
+++ b/src/systems.jl
@@ -40,6 +40,9 @@ for (Z, AZ) in ((:ContinuousIdentitySystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -119,6 +122,9 @@ for (Z, AZ) in ((:ConstrainedContinuousIdentitySystem, :AbstractContinuousSystem
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -205,6 +211,9 @@ for (Z, AZ) in ((:LinearContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -295,6 +304,9 @@ for (Z, AZ) in ((:AffineContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -337,7 +349,8 @@ for (Z, AZ) in ((:LinearControlContinuousSystem, :AbstractContinuousSystem),
         struct $(Z){T,MTA<:AbstractMatrix{T},MTB<:AbstractMatrix{T}} <: $(AZ)
             A::MTA
             B::MTB
-            function $(Z)(A::MTA, B::MTB) where {T,MTA<:AbstractMatrix{T},MTB<:AbstractMatrix{T}}
+            function $(Z)(A::MTA,
+                          B::MTB) where {T,MTA<:AbstractMatrix{T},MTB<:AbstractMatrix{T}}
                 if checksquare(A) != size(B, 1)
                     throw(DimensionMismatch("incompatible dimensions"))
                 end
@@ -387,6 +400,9 @@ for (Z, AZ) in ((:LinearControlContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -477,6 +493,9 @@ for (Z, AZ) in ((:ConstrainedLinearContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -574,6 +593,9 @@ for (Z, AZ) in ((:ConstrainedAffineContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -615,7 +637,8 @@ ConstrainedAffineDiscreteSystem
 for (Z, AZ) in ((:AffineControlContinuousSystem, :AbstractContinuousSystem),
                 (:AffineControlDiscreteSystem, :AbstractDiscreteSystem))
     @eval begin
-        struct $(Z){T,MTA<:AbstractMatrix{T},MTB<:AbstractMatrix{T},VT<:AbstractVector{T}} <: $(AZ)
+        struct $(Z){T,MTA<:AbstractMatrix{T},MTB<:AbstractMatrix{T},VT<:AbstractVector{T}} <:
+               $(AZ)
             A::MTA
             B::MTB
             c::VT
@@ -674,6 +697,9 @@ for (Z, AZ) in ((:AffineControlContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -715,7 +741,8 @@ AffineControlDiscreteSystem
 for (Z, AZ) in ((:ConstrainedAffineControlContinuousSystem, :AbstractContinuousSystem),
                 (:ConstrainedAffineControlDiscreteSystem, :AbstractDiscreteSystem))
     @eval begin
-        struct $(Z){T,MTA<:AbstractMatrix{T},MTB<:AbstractMatrix{T},VT<:AbstractVector{T},ST,UT} <:
+        struct $(Z){T,MTA<:AbstractMatrix{T},MTB<:AbstractMatrix{T},VT<:AbstractVector{T},ST,
+                    UT} <:
                $(AZ)
             A::MTA
             B::MTB
@@ -783,6 +810,9 @@ for (Z, AZ) in ((:ConstrainedAffineControlContinuousSystem, :AbstractContinuousS
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -890,6 +920,9 @@ for (Z, AZ) in ((:ConstrainedLinearControlContinuousSystem, :AbstractContinuousS
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -935,7 +968,8 @@ for (Z, AZ) in ((:LinearDescriptorContinuousSystem, :AbstractContinuousSystem),
         struct $(Z){T,MTA<:AbstractMatrix{T},MTE<:AbstractMatrix{T}} <: $(AZ)
             A::MTA
             E::MTE
-            function $(Z)(A::MTA, E::MTE) where {T,MTA<:AbstractMatrix{T},MTE<:AbstractMatrix{T}}
+            function $(Z)(A::MTA,
+                          E::MTE) where {T,MTA<:AbstractMatrix{T},MTE<:AbstractMatrix{T}}
                 if size(A) != size(E)
                     throw(DimensionMismatch("incompatible dimensions"))
                 end
@@ -982,6 +1016,9 @@ for (Z, AZ) in ((:LinearDescriptorContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -1076,6 +1113,9 @@ for (Z, AZ) in ((:ConstrainedLinearDescriptorContinuousSystem, :AbstractContinuo
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -1173,6 +1213,9 @@ for (Z, AZ) in ((:PolynomialContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -1217,7 +1260,8 @@ for (Z, AZ) in ((:ConstrainedPolynomialContinuousSystem, :AbstractContinuousSyst
             statedim::Int
             X::ST
             function $(Z)(p::VPT, statedim::Int,
-                          X::ST) where {T,PT<:AbstractPolynomialLike{T},VPT<:AbstractVector{PT},ST}
+                          X::ST) where {T,PT<:AbstractPolynomialLike{T},VPT<:AbstractVector{PT},
+                                        ST}
                 if statedim != MultivariatePolynomials.nvariables(p)
                     throw(DimensionMismatch("the state dimension $(statedim) does not match the number of state variables"))
                 end
@@ -1271,6 +1315,9 @@ for (Z, AZ) in ((:ConstrainedPolynomialContinuousSystem, :AbstractContinuousSyst
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -1352,6 +1399,9 @@ for (Z, AZ) in ((:BlackBoxContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -1435,6 +1485,9 @@ for (Z, AZ) in ((:ConstrainedBlackBoxContinuousSystem, :AbstractContinuousSystem
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -1519,6 +1572,9 @@ for (Z, AZ) in ((:BlackBoxControlContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -1609,6 +1665,9 @@ for (Z, AZ) in ((:ConstrainedBlackBoxControlContinuousSystem, :AbstractContinuou
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -1663,7 +1722,8 @@ for (Z, AZ) in ((:NoisyLinearContinuousSystem, :AbstractContinuousSystem),
         struct $(Z){T,MTA<:AbstractMatrix{T},MTD<:AbstractMatrix{T}} <: $(AZ)
             A::MTA
             D::MTD
-            function $(Z)(A::MTA, D::MTD) where {T,MTA<:AbstractMatrix{T},MTD<:AbstractMatrix{T}}
+            function $(Z)(A::MTA,
+                          D::MTD) where {T,MTA<:AbstractMatrix{T},MTD<:AbstractMatrix{T}}
                 if checksquare(A) != size(D, 1)
                     throw(DimensionMismatch("incompatible dimensions"))
                 end
@@ -1713,6 +1773,9 @@ for (Z, AZ) in ((:NoisyLinearContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -1814,6 +1877,9 @@ for (Z, AZ) in ((:NoisyConstrainedLinearContinuousSystem, :AbstractContinuousSys
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -1857,7 +1923,8 @@ NoisyConstrainedLinearDiscreteSystem
 for (Z, AZ) in ((:NoisyLinearControlContinuousSystem, :AbstractContinuousSystem),
                 (:NoisyLinearControlDiscreteSystem, :AbstractDiscreteSystem))
     @eval begin
-        struct $(Z){T,MTA<:AbstractMatrix{T},MTB<:AbstractMatrix{T},MTD<:AbstractMatrix{T}} <: $(AZ)
+        struct $(Z){T,MTA<:AbstractMatrix{T},MTB<:AbstractMatrix{T},MTD<:AbstractMatrix{T}} <:
+               $(AZ)
             A::MTA
             B::MTB
             D::MTD
@@ -1916,6 +1983,9 @@ for (Z, AZ) in ((:NoisyLinearControlContinuousSystem, :AbstractContinuousSystem)
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -1959,7 +2029,8 @@ NoisyLinearControlDiscreteSystem
 for (Z, AZ) in ((:NoisyConstrainedLinearControlContinuousSystem, :AbstractContinuousSystem),
                 (:NoisyConstrainedLinearControlDiscreteSystem, :AbstractDiscreteSystem))
     @eval begin
-        struct $(Z){T,MTA<:AbstractMatrix{T},MTB<:AbstractMatrix{T},MTD<:AbstractMatrix{T},ST,UT,
+        struct $(Z){T,MTA<:AbstractMatrix{T},MTB<:AbstractMatrix{T},MTD<:AbstractMatrix{T},ST,
+                    UT,
                     WT} <: $(AZ)
             A::MTA
             B::MTB
@@ -2031,6 +2102,9 @@ for (Z, AZ) in ((:NoisyConstrainedLinearControlContinuousSystem, :AbstractContin
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -2146,6 +2220,9 @@ for (Z, AZ) in ((:NoisyAffineControlContinuousSystem, :AbstractContinuousSystem)
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -2265,6 +2342,9 @@ for (Z, AZ) in ((:NoisyConstrainedAffineControlContinuousSystem, :AbstractContin
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -2358,6 +2438,9 @@ for (Z, AZ) in ((:NoisyBlackBoxControlContinuousSystem, :AbstractContinuousSyste
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -2457,6 +2540,9 @@ for (Z, AZ) in ((:NoisyConstrainedBlackBoxControlContinuousSystem, :AbstractCont
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -2604,6 +2690,9 @@ for (Z, AZ) in ((:SecondOrderLinearContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -2710,6 +2799,9 @@ for (Z, AZ) in ((:SecondOrderAffineContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -2754,8 +2846,9 @@ Discrete-time second order constrained linear control system of the form:
 """
 SecondOrderConstrainedLinearControlDiscreteSystem
 
-for (Z, AZ) in ((:SecondOrderConstrainedLinearControlContinuousSystem, :AbstractContinuousSystem),
-                (:SecondOrderConstrainedLinearControlDiscreteSystem, :AbstractDiscreteSystem))
+for (Z, AZ) in
+    ((:SecondOrderConstrainedLinearControlContinuousSystem, :AbstractContinuousSystem),
+     (:SecondOrderConstrainedLinearControlDiscreteSystem, :AbstractDiscreteSystem))
     @eval begin
         struct $(Z){T,MTM<:AbstractMatrix{T},
                     MTC<:AbstractMatrix{T},
@@ -2837,6 +2930,9 @@ for (Z, AZ) in ((:SecondOrderConstrainedLinearControlContinuousSystem, :Abstract
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -2883,8 +2979,9 @@ Discrete-time second order constrained affine control system of the form:
 """
 SecondOrderConstrainedAffineControlDiscreteSystem
 
-for (Z, AZ) in ((:SecondOrderConstrainedAffineControlContinuousSystem, :AbstractContinuousSystem),
-                (:SecondOrderConstrainedAffineControlDiscreteSystem, :AbstractDiscreteSystem))
+for (Z, AZ) in
+    ((:SecondOrderConstrainedAffineControlContinuousSystem, :AbstractContinuousSystem),
+     (:SecondOrderConstrainedAffineControlDiscreteSystem, :AbstractDiscreteSystem))
     @eval begin
         struct $(Z){T,MTM<:AbstractMatrix{T},
                     MTC<:AbstractMatrix{T},
@@ -2908,7 +3005,8 @@ for (Z, AZ) in ((:SecondOrderConstrainedAffineControlContinuousSystem, :Abstract
                                         VT<:AbstractVector{T},
                                         ST,
                                         UT}
-                if !(checksquare(M) == checksquare(C) == checksquare(K) == size(B, 1) == length(d))
+                if !(checksquare(M) == checksquare(C) == checksquare(K) == size(B, 1) ==
+                     length(d))
                     throw(DimensionMismatch("incompatible dimensions"))
                 end
                 return new{T,MTM,MTC,MTK,MTB,VT,ST,UT}(M, C, K, B, d, X, U)
@@ -2969,6 +3067,9 @@ for (Z, AZ) in ((:SecondOrderConstrainedAffineControlContinuousSystem, :Abstract
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -3075,6 +3176,9 @@ for (Z, AZ) in ((:SecondOrderContinuousSystem, :AbstractContinuousSystem),
             function isconstrained(::$T)
                 return false
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
@@ -3136,7 +3240,8 @@ for (Z, AZ) in ((:SecondOrderConstrainedContinuousSystem, :AbstractContinuousSys
             U::UT
 
             function $(Z)(M::MTM, C::MTC, fi::FI, fe::FE, X::ST,
-                          U::UT) where {T,MTM<:AbstractMatrix{T},MTC<:AbstractMatrix{T},FI,FE,ST,UT}
+                          U::UT) where {T,MTM<:AbstractMatrix{T},MTC<:AbstractMatrix{T},FI,FE,
+                                        ST,UT}
                 if checksquare(M) != checksquare(C)
                     throw(DimensionMismatch("incompatible dimensions"))
                 end
@@ -3186,6 +3291,194 @@ for (Z, AZ) in ((:SecondOrderConstrainedContinuousSystem, :AbstractContinuousSys
             function isconstrained(::$T)
                 return true
             end
+            function isparametric(::$T)
+                return false
+            end
         end
     end
 end
+
+function __init__()
+    @require LazySets = "b4f0291d-fe17-52bc-9479-3d1a343d9043" begin
+        using .LazySets: MatrixZonotope
+
+        export LinearParametricContinuousSystem,
+               LinearParametricDiscreteSystem,
+               LinearControlParametricContinuousSystem,
+               LinearControlParametricDiscreteSystem
+
+        @doc """
+            LinearParametricContinuousSystem
+
+        Continuous-time linear  parametric system of the form:
+
+        ```math
+            x(t)' = A(\\theta) x(t), \\theta ∈ \\Theta \\; \\forall t
+        ```
+        where ``A(θ)`` belongs to a continuous set of matrices, e.g., an interval
+        matrix, matrix zonotope, or other convex matrix sets.
+
+        ### Fields
+        - `A`      --  parametric state matrix
+        """
+        LinearParametricContinuousSystem
+
+        @doc """
+            LinearParametricDiscreteSystem
+
+        Discrete-time linear parametric system of the form:
+        ```math
+            x_{k+1} = A(\\theta) x_k, \\theta ∈ \\Theta \\; \\forall k
+        ```
+        where ``A(θ)`` belongs to a continuous set of matrices, e.g., an interval
+        matrix, matrix zonotope, or other convex matrix sets.
+
+        ### Fields
+        - `A`      --  parametric state matrix
+        """
+        LinearParametricDiscreteSystem
+
+        for (Z, AZ) in ((:LinearParametricContinuousSystem, :AbstractContinuousSystem),
+                        (:LinearParametricDiscreteSystem, :AbstractDiscreteSystem))
+            BaseName = Symbol(replace(string(Z), "Parametric" => ""))
+
+            @eval begin
+                struct $(Z){T,MA<:MatrixZonotope{T}} <: $(AZ)
+                    A::MA
+                end
+
+                statedim(s::$Z) = size(s.A, 1)
+                inputdim(::$Z) = 0
+                noisedim(::$Z) = 0
+                state_matrix(s::$Z) = s.A
+
+                function $(BaseName)(M::MatrixZonotope)
+                    return $Z(M)
+                end
+            end
+
+            for T in [Z, Type{<:eval(Z)}]
+                @eval begin
+                    function islinear(::$T)
+                        return true
+                    end
+                    function isaffine(::$T)
+                        return true
+                    end
+                    function ispolynomial(::$T)
+                        return false
+                    end
+                    function isnoisy(::$T)
+                        return false
+                    end
+                    function iscontrolled(::$T)
+                        return false
+                    end
+                    function isconstrained(::$T)
+                        return false
+                    end
+                    function isparametric(::$T)
+                        return true
+                    end
+                end
+            end
+        end
+
+        @doc """
+            LinearControlParametricContinuousSystem
+
+        Continuous-time linear parametric control system of the form:
+
+        ```math
+            x(t)' = A(θ) x(t) + B(θ) u(t), \\theta ∈ \\Theta \\; \\forall t
+        ```
+
+        where ``A(θ)`` and ``B(θ)`` belong to continuous sets of matrices, e.g., an interval
+        matrix, matrix zonotope, or other convex matrix sets.
+
+        ### Fields
+        - `A`      --  parametric state matrix
+        - `B`      --  parametric input matrix
+        """
+        LinearControlParametricContinuousSystem
+
+        @doc """
+            LinearControlParametricDiscreteSystem 
+
+        Discrete-time linear parametric control system of the form:
+
+        ```math
+            x_{k+1} = A(θ) x_k + B(θ) u_k, \\theta ∈ \\Theta \\; \\forall k
+        ```
+
+        where ``A(θ)`` and ``B(θ)`` belong to continuous sets of matrices, e.g., an interval
+        matrix, matrix zonotope, or other convex matrix sets.
+
+        ### Fields
+        - `A`      --  parametric state matrix
+        - `B`      --  parametric input matrix
+        """
+        LinearControlParametricDiscreteSystem
+
+        for (Z, AZ) in
+            ((:LinearControlParametricContinuousSystem, :AbstractContinuousSystem),
+             (:LinearControlParametricDiscreteSystem, :AbstractDiscreteSystem))
+            BaseName = Symbol(replace(string(Z), "Parametric" => ""))
+
+            @eval begin
+                struct $(Z){T,MTA<:MatrixZonotope{T},MTB<:MatrixZonotope{T}} <: $(AZ)
+                    A::MTA
+                    B::MTB
+                    function $(Z)(A::MTA,
+                                  B::MTB) where {T,MTA<:MatrixZonotope{T},
+                                                 MTB<:MatrixZonotope{T}}
+                        if checksquare(A) != size(B, 1)
+                            throw(DimensionMismatch("incompatible dimensions"))
+                        end
+                        return new{T,MTA,MTB}(A, B)
+                    end
+                end
+
+                function $(Z)(A::Number, B::Number)
+                    return $(Z)(hcat(A), hcat(B))
+                end
+
+                statedim(s::$Z) = size(s.A, 1)
+                inputdim(s::$Z) = size(s.B, 2)
+                noisedim(::$Z) = 0
+                state_matrix(s::$Z) = s.A
+                input_matrix(s::$Z) = s.B
+
+                function $(BaseName)(MA::MatrixZonotope, MB::MatrixZonotope)
+                    return $Z(MA, MB)
+                end
+            end
+
+            for T in [Z, Type{<:eval(Z)}]
+                @eval begin
+                    function islinear(::$T)
+                        return true
+                    end
+                    function isaffine(::$T)
+                        return true
+                    end
+                    function ispolynomial(::$T)
+                        return false
+                    end
+                    function isnoisy(::$T)
+                        return false
+                    end
+                    function iscontrolled(::$T)
+                        return true
+                    end
+                    function isconstrained(::$T)
+                        return false
+                    end
+                    function isparametric(::$T)
+                        return true
+                    end
+                end
+            end
+        end
+    end
+end
\ No newline at end of file
diff --git a/test/@system.jl b/test/@system.jl
index 4225ea6d..7ed991b3 100644
--- a/test/@system.jl
+++ b/test/@system.jl
@@ -1,7 +1,9 @@
 @testset "_corresponding_type" begin
     for AS in [AbstractContinuousSystem, AbstractDiscreteSystem]
         for S in subtypes(AS)
-            @test MathematicalSystems._corresponding_type(AS, fieldnames.(S)) == S
+            if S ∉ vcat(PARAMETRIC_CTYPES, PARAMETRIC_DTYPES) #TODO refactor _corresponding_type
+                @test MathematicalSystems._corresponding_type(AS, fieldnames.(S)) == S
+            end
         end
     end
 end
diff --git a/test/Project.toml b/test/Project.toml
index d162affa..b2981b71 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -10,6 +10,6 @@ TypedPolynomials = "afbbf031-7a57-5f58-a1b9-b774a0fad08d"
 
 [compat]
 Aqua = "0.8.9"
-LazySets = "1, 2, 3, 4"
+LazySets = "5"
 RecipesBase = "1"
 TypedPolynomials = "0.3, 0.4"
diff --git a/test/continuous.jl b/test/continuous.jl
index cb6fa398..fc9e9c4d 100644
--- a/test/continuous.jl
+++ b/test/continuous.jl
@@ -55,7 +55,7 @@ Es = [e][:, :];
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -73,7 +73,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -91,7 +91,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = LinearContinuousSystem(a)
@@ -112,7 +112,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test !islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = AffineContinuousSystem(a, c)
@@ -133,7 +133,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = LinearControlContinuousSystem(a, b)
@@ -154,7 +154,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = ConstrainedLinearContinuousSystem(a, Xs)
@@ -175,7 +175,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test !islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = ConstrainedAffineContinuousSystem(a, c, Xs)
@@ -196,7 +196,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test !islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = ConstrainedAffineControlContinuousSystem(a, b, c, Xs, Us)
@@ -216,7 +216,7 @@ end
     @test affine_term(s) == c
     for s in [s, typeof(s)]
         @test !islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     a = 1.0
@@ -243,7 +243,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
 
     @testset "initial value problem composite type" begin
@@ -300,7 +300,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = ConstrainedLinearDescriptorContinuousSystem(a, e, Xs)
@@ -330,7 +330,7 @@ end
     @test mapping(s) == [p]
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
 
     @test TypedPolynomials.nvariables(s) == 2
@@ -357,7 +357,7 @@ end
     @test mapping(s) == [p]
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
 
     @test TypedPolynomials.nvariables(s) == 2
@@ -390,7 +390,7 @@ end
 
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -415,7 +415,7 @@ end
     @test stateset(s) == H
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -441,7 +441,7 @@ end
     @test mapping(s) == vanderpol_controlled!
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
-        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -467,7 +467,7 @@ end
     @test dx ≈ [0.0, -1.0 + u[1]]
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
-        @test !isnoisy(s) && iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -489,7 +489,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
+        @test isnoisy(s) && !iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = NoisyLinearContinuousSystem(a, d)
@@ -510,7 +510,7 @@ end
     @test noiseset(s) == W
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test isnoisy(s) && !iscontrolled(s) && isconstrained(s)
+        @test isnoisy(s) && !iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = NoisyConstrainedLinearContinuousSystem(a, d, Xs, Ws)
@@ -531,7 +531,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s)
+        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = NoisyLinearControlContinuousSystem(a, b, d)
@@ -552,7 +552,7 @@ end
     @test noiseset(s) == W
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test isnoisy(s) && iscontrolled(s) && isconstrained(s)
+        @test isnoisy(s) && iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = NoisyConstrainedLinearControlContinuousSystem(a, b, d, Xs, Us, Ws)
@@ -573,7 +573,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test isaffine(s) && !islinear(s) && !ispolynomial(s) && !isblackbox(s)
-        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s)
+        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = NoisyAffineControlContinuousSystem(a, b, c, d)
@@ -594,7 +594,7 @@ end
     @test noiseset(s) == W
     for s in [s, typeof(s)]
         @test isaffine(s) && !islinear(s) && !ispolynomial(s) && !isblackbox(s)
-        @test isnoisy(s) && iscontrolled(s) && isconstrained(s)
+        @test isnoisy(s) && iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = NoisyConstrainedAffineControlContinuousSystem(a, b, c, d, Xs, Us, Ws)
@@ -617,7 +617,7 @@ end
     @test mapping(s) == add_one
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
-        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s)
+        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -636,7 +636,7 @@ end
     @test noiseset(s) == W
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
-        @test isnoisy(s) && iscontrolled(s) && isconstrained(s)
+        @test isnoisy(s) && iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -664,6 +664,7 @@ end
     @test !isnoisy(s)
     @test !iscontrolled(s)
     @test !isconstrained(s)
+    @test !isparametric(s)
 
     s = SecondOrderLinearContinuousSystem(2, 3, 4)
     @test mass_matrix(s) == hcat(2)
@@ -682,6 +683,7 @@ end
     @test !isnoisy(s)
     @test !iscontrolled(s)
     @test !isconstrained(s)
+    @test !isparametric(s)
 
     s = SecondOrderAffineContinuousSystem(2, 3, 4, 6)
     @test mass_matrix(s) == hcat(2)
@@ -702,6 +704,7 @@ end
     @test !isnoisy(s)
     @test iscontrolled(s)
     @test isconstrained(s)
+    @test !isparametric(s)
 
     s = SecondOrderConstrainedLinearControlContinuousSystem(2, 3, 4, 5, X1, U1)
     @test mass_matrix(s) == hcat(2)
@@ -724,6 +727,7 @@ end
     @test !isnoisy(s)
     @test iscontrolled(s)
     @test isconstrained(s)
+    @test !isparametric(s)
 
     s = SecondOrderConstrainedAffineControlContinuousSystem(2, 3, 4, 5, 6, X1, U1)
     @test mass_matrix(s) == hcat(2)
@@ -751,6 +755,7 @@ end
     @test !isnoisy(s)
     @test !iscontrolled(s)
     @test !isconstrained(s)
+    @test !isparametric(s)
 
     s = SecondOrderContinuousSystem(2, 3, fi, fe)
     @test mass_matrix(s) == hcat(2)
@@ -773,6 +778,7 @@ end
     @test !isnoisy(s)
     @test !iscontrolled(s)
     @test isconstrained(s)
+    @test !isparametric(s)
 
     s = SecondOrderConstrainedContinuousSystem(2, 3, fi, fe, X, U)
     @test mass_matrix(s) == hcat(2)
@@ -780,3 +786,59 @@ end
     @test stateset(s) === X && inputset(s) === U
     @test s.fi == fi && s.fe == fe
 end
+
+# ==================
+# Parametric systems
+# ==================
+
+@testset "Linear parametric continuous systems" begin
+    @static if isdefined(@__MODULE__, :LazySets)
+        using LazySets: MatrixZonotope
+        
+        Ac = [1.0 0.0; 0.0 1.0]
+        A1 = [0.1 0.0; 0.0 0.0]
+        A2 = [0.0 0.0; 0.0 0.2]
+        A = MatrixZonotope(Ac, [A1, A2])
+
+        s = LinearParametricContinuousSystem(A)
+        @test isa(s, LinearParametricContinuousSystem)
+
+        # shortcut constructor
+        s = LinearContinuousSystem(A)
+        @test isa(s, LinearParametricContinuousSystem)
+
+        @test statedim(s) == 2
+        @test inputdim(s) == 0
+        @test noisedim(s) == 0
+        @test state_matrix(s) === A
+        @test islinear(s)
+        @test isparametric(s)
+        @test isaffine(s)
+        @test !isnoisy(s)
+        @test !iscontrolled(s)
+        @test !isconstrained(s)
+        
+        # control case
+        Bc = hcat([1.0; 0.5])
+        B1 = hcat([0.05; 0.0])
+        B = MatrixZonotope(Bc, [B1])
+
+        sc = LinearControlParametricContinuousSystem(A, B)
+        @test isa(sc, LinearControlParametricContinuousSystem)
+        # shortcut constructor for control
+        sc = LinearControlContinuousSystem(A, B)
+        @test isa(sc, LinearControlParametricContinuousSystem)
+
+        @test statedim(sc) == 2
+        @test inputdim(sc) == 1
+        @test noisedim(sc) == 0
+        @test state_matrix(sc) === A
+        @test input_matrix(sc) === B
+        @test islinear(sc)
+        @test isparametric(sc)
+        @test isaffine(sc)
+        @test iscontrolled(sc)
+        @test !isnoisy(sc)
+        @test !isconstrained(sc)
+    end
+end
diff --git a/test/discrete.jl b/test/discrete.jl
index 9bbc5dc8..27ccdb40 100644
--- a/test/discrete.jl
+++ b/test/discrete.jl
@@ -53,7 +53,7 @@ Es = [e][:, :];
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -71,7 +71,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -89,7 +89,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = LinearDiscreteSystem(a)
@@ -110,7 +110,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test !islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = AffineDiscreteSystem(a, c)
@@ -131,7 +131,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = LinearControlDiscreteSystem(a, b)
@@ -152,7 +152,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = ConstrainedLinearDiscreteSystem(a, Xs)
@@ -173,7 +173,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test !islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = ConstrainedAffineDiscreteSystem(a, c, Xs)
@@ -194,7 +194,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = ConstrainedLinearControlDiscreteSystem(a, b, Xs, Us)
@@ -216,7 +216,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = LinearDescriptorDiscreteSystem(a, e)
@@ -238,7 +238,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = ConstrainedLinearDescriptorDiscreteSystem(a, e, Xs)
@@ -260,7 +260,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -279,7 +279,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && ispolynomial(s) && !isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -298,7 +298,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
-        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && !iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -317,7 +317,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
-        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s)
+        @test !isnoisy(s) && iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -339,7 +339,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test isnoisy(s) && !iscontrolled(s) && !isconstrained(s)
+        @test isnoisy(s) && !iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = NoisyLinearDiscreteSystem(a, d)
@@ -360,7 +360,7 @@ end
     @test noiseset(s) == W
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test isnoisy(s) && !iscontrolled(s) && isconstrained(s)
+        @test isnoisy(s) && !iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = NoisyConstrainedLinearDiscreteSystem(a, d, Xs, Ws)
@@ -381,7 +381,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s)
+        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = NoisyLinearControlDiscreteSystem(a, b, d)
@@ -402,7 +402,7 @@ end
     @test noiseset(s) == W
     for s in [s, typeof(s)]
         @test islinear(s) && isaffine(s) && !ispolynomial(s) && !isblackbox(s)
-        @test isnoisy(s) && iscontrolled(s) && isconstrained(s)
+        @test isnoisy(s) && iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = NoisyConstrainedLinearControlDiscreteSystem(a, b, d, Xs, Us, Ws)
@@ -423,7 +423,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test isaffine(s) && !islinear(s) && !ispolynomial(s) && !isblackbox(s)
-        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s)
+        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = NoisyAffineControlDiscreteSystem(a, b, c, d)
@@ -444,7 +444,7 @@ end
     @test noiseset(s) == W
     for s in [s, typeof(s)]
         @test isaffine(s) && !islinear(s) && !ispolynomial(s) && !isblackbox(s)
-        @test isnoisy(s) && iscontrolled(s) && isconstrained(s)
+        @test isnoisy(s) && iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
     # Scalar System
     scalar_sys = NoisyConstrainedAffineControlDiscreteSystem(a, b, c, d, Xs, Us, Ws)
@@ -466,7 +466,7 @@ end
     @test isnothing(noiseset(s))
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
-        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s)
+        @test isnoisy(s) && iscontrolled(s) && !isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -485,7 +485,7 @@ end
     @test noiseset(s) == W
     for s in [s, typeof(s)]
         @test !islinear(s) && !isaffine(s) && !ispolynomial(s) && isblackbox(s)
-        @test isnoisy(s) && iscontrolled(s) && isconstrained(s)
+        @test isnoisy(s) && iscontrolled(s) && isconstrained(s) && !isparametric(s)
     end
 end
 
@@ -509,3 +509,57 @@ end
         @test ui isa Hyperrectangle
     end
 end
+
+# ==================
+# Parametric systems
+# ==================
+
+@testset "Linear parametric discrete systems" begin
+    @static if isdefined(@__MODULE__, :LazySets)
+        using LazySets: MatrixZonotope
+        
+        Ac = [1.0 0.0; 0.0 1.0]
+        A1 = [0.1 0.0; 0.0 0.0]
+        A2 = [0.0 0.0; 0.0 0.2]
+        A = MatrixZonotope(Ac, [A1, A2])
+
+        s = LinearParametricDiscreteSystem(A)
+        @test isa(s, LinearParametricDiscreteSystem)
+        # shortcut constructor
+        s = LinearDiscreteSystem(A)
+        @test isa(s, LinearParametricDiscreteSystem)
+
+        @test statedim(s) == 2
+        @test inputdim(s) == 0
+        @test noisedim(s) == 0
+        @test state_matrix(s) === A
+        @test islinear(s)
+        @test isparametric(s)
+        @test isaffine(s)
+        @test !isnoisy(s)
+        @test !iscontrolled(s)
+        @test !isconstrained(s)
+        
+        Bc = hcat([1.0; 0.5])
+        B1 = hcat([0.05; 0.0])
+        B = MatrixZonotope(Bc, [B1])
+
+        sc = LinearControlParametricDiscreteSystem(A, B)
+        @test isa(sc, LinearControlParametricDiscreteSystem)
+        # shortcut constructor
+        sc = LinearControlDiscreteSystem(A, B)
+        @test isa(sc, LinearControlParametricDiscreteSystem)
+
+        @test statedim(sc) == 2
+        @test inputdim(sc) == 1
+        @test noisedim(sc) == 0
+        @test state_matrix(sc) === A
+        @test input_matrix(sc) === B
+        @test islinear(sc)
+        @test isparametric(sc)
+        @test isaffine(sc)
+        @test iscontrolled(sc)
+        @test !isnoisy(sc)
+        @test !isconstrained(sc)
+    end
+end
diff --git a/test/discretize.jl b/test/discretize.jl
index e2a80ba3..6973be0c 100644
--- a/test/discretize.jl
+++ b/test/discretize.jl
@@ -1,15 +1,19 @@
 
 @testset "Convert continuous to discrete type" begin
     for dtype in subtypes(AbstractDiscreteSystem)
-        ctype = eval.(Meta.parse.(replace(string(dtype), "Discrete" => "Continuous")))
-        @test MathematicalSystems._complementary_type(dtype) == ctype
-        @inferred MathematicalSystems._complementary_type(dtype)
+        if dtype ∉ PARAMETRIC_DTYPES
+            ctype = eval.(Meta.parse.(replace(string(dtype), "Discrete" => "Continuous")))
+            @test MathematicalSystems._complementary_type(dtype) == ctype
+            @inferred MathematicalSystems._complementary_type(dtype)
+        end
     end
 
     for ctype in subtypes(AbstractContinuousSystem)
-        dtype = eval.(Meta.parse.(replace(string(ctype), "Continuous" => "Discrete")))
-        @test MathematicalSystems._complementary_type(ctype) == dtype
-        @inferred MathematicalSystems._complementary_type(ctype)
+        if ctype ∉ PARAMETRIC_CTYPES
+            dtype = eval.(Meta.parse.(replace(string(ctype), "Continuous" => "Discrete")))
+            @test MathematicalSystems._complementary_type(ctype) == dtype
+            @inferred MathematicalSystems._complementary_type(ctype)
+        end
     end
 end
 
@@ -37,8 +41,8 @@ end
     CTYPES = filter(x -> (occursin("Linear", string(x)) || occursin("Affine", string(x))) &&
                         !occursin("Descriptor", string(x)), subtypes(AbstractContinuousSystem))
 
-    # this test doesn't apply for second order systems
-    filter!(x -> x ∉ SECOND_ORDER_CTYPES, CTYPES)
+    # this test doesn't apply for second order systems and parametric systems
+    filter!(x -> x ∉ SECOND_ORDER_CTYPES && x ∉ PARAMETRIC_CTYPES, CTYPES)
 
     DTYPES = MathematicalSystems._complementary_type.(CTYPES)
     n_types = length(CTYPES)
@@ -78,8 +82,8 @@ end
     CTYPES = filter(x -> (occursin("Linear", string(x)) || occursin("Affine", string(x))) &&
                         !occursin("Descriptor", string(x)), subtypes(AbstractContinuousSystem))
 
-    # this test doesn't apply for second order systems
-    filter!(x -> x ∉ SECOND_ORDER_CTYPES, CTYPES)
+    # this test doesn't apply for second order systems and parametric systems
+    filter!(x -> x ∉ SECOND_ORDER_CTYPES && x ∉ PARAMETRIC_CTYPES, CTYPES)
 
     DTYPES = MathematicalSystems._complementary_type.(CTYPES)
     n_types = length(CTYPES)
diff --git a/test/runtests.jl b/test/runtests.jl
index 8f5ac25b..968ff7c1 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -11,6 +11,13 @@ SECOND_ORDER_CTYPES = [SecondOrderAffineContinuousSystem,
                        SecondOrderConstrainedLinearControlContinuousSystem,
                        SecondOrderLinearContinuousSystem]
 
+PARAMETRIC_CTYPES = [LinearParametricContinuousSystem,
+                     LinearControlParametricContinuousSystem]
+
+
+PARAMETRIC_DTYPES = [LinearParametricDiscreteSystem,
+                     LinearControlParametricDiscreteSystem]
+
 @testset "Utilities" begin
     include("utilities.jl")
 end