diff --git a/lib/ModelingToolkitTearing/src/reassemble.jl b/lib/ModelingToolkitTearing/src/reassemble.jl
index c3ddb66..812ccaf 100644
--- a/lib/ModelingToolkitTearing/src/reassemble.jl
+++ b/lib/ModelingToolkitTearing/src/reassemble.jl
@@ -503,7 +503,7 @@ function generate_system_equations!(state::TearingState, neweqs::Vector{Equation
     offset = 1
     findnextfn = let diff_vars_set = diff_vars_set, solved_vars_set = solved_vars_set,
         diff_to_var = diff_to_var, ispresent = ispresent
-        j -> !(j in diff_vars_set || j in solved_vars_set) && diff_to_var[j] === nothing &&
+        j -> !(j in diff_vars_set || j in solved_vars_set || j in extra_vars) && diff_to_var[j] === nothing &&
             ispresent(j)
     end
     for (i, v) in enumerate(var_ordering)
@@ -1028,6 +1028,15 @@ function update_simplified_system!(
     unknown_idxs = filter(filterer, eachindex(state.fullvars))
     unknowns = state.fullvars[unknown_idxs]
     unknowns = [unknowns; extra_unknowns]
+    # `generate_system_equations!` will include the extra unknowns in `var_ordering` due to
+    # the `setdiff` at the end. Excluding them leads to issues when running `singularity_check`
+    # on the simplified `state` in `InitializationProblem`. However, `ispresent` can filter
+    # some of them out of `unknown_idxs` if they don't occur in the selected equations. As
+    # a result, we still add them here. `unique!` serves to remove duplicates. Otherwise,
+    # `tearing_hacks` may count an array variable occuring `length + 1` times in
+    # unknowns+observables, and will thus not add the `x ~ [x[1], x[2]...]` array hack
+    # observed equation.
+    unique!(unknowns)
     if StateSelection.is_only_discrete(structure)
         # Algebraic variables are shifted forward by one, so we backshift them.
         _unknowns = SymbolicT[]
diff --git a/lib/ModelingToolkitTearing/test/runtests.jl b/lib/ModelingToolkitTearing/test/runtests.jl
index 5e3116e..5488f7e 100644
--- a/lib/ModelingToolkitTearing/test/runtests.jl
+++ b/lib/ModelingToolkitTearing/test/runtests.jl
@@ -220,3 +220,12 @@ end
     @test length(equations(ts)) == 1
     @test length(ts.fullvars) == 2
 end
+
+@testset "Unknowns in underdetermined systems aren't counted twice" begin
+    @variables x[1:3]
+    @mtkcompile sys = System([x[1] ~ 2x[2]+3x[3], x[2]^3 + x[3]^3 - 9 ~ 0], [x[1], x[2], x[3]], []) fully_determined=false
+    # Without the fix in https://github.com/JuliaComputing/StateSelection.jl/pull/80
+    # `unknowns` prior to running `tearing_hacks` will contain `x[3]` twice, so
+    # the observed equation `x ~ [x[1], x[2], x[3]]` won't be added.
+    @test isequal(observables(sys)[2], x)
+end