diff --git a/src/descriptor.jl b/src/descriptor.jl
index ebd4cd8..48a1af6 100644
--- a/src/descriptor.jl
+++ b/src/descriptor.jl
@@ -87,9 +87,16 @@ Compute the a balanced truncation of order `n` and the hankel singular values
 For keyword arguments, see the docstring of `DescriptorSystems.gbalmr`, reproduced below
 $(@doc(DescriptorSystems.gbalmr))
 """
-function baltrunc2(sys::LTISystem; residual=false, n=missing, kwargs...)
-    sysr, hs = DescriptorSystems.gbalmr(dss(sys); matchdc=residual, ord=n, kwargs...)
-    ss(sysr), hs
+function baltrunc2(sys::LTISystem; residual=false, n=missing, scaleY=1.0, scaleU=1.0, kwargs...)
+    # Apply scaling if needed
+    A, B, C, D = ssdata(sys)
+    # Divide by scaling factors to normalize to ~[-1,1]
+    sys_scaled = ss(A, B / scaleU, C / scaleY, scaleY \ D / scaleU, sys.timeevol)
+    sysr, hs = DescriptorSystems.gbalmr(dss(sys_scaled); matchdc=residual, ord=n, kwargs...)
+    # Multiply by scaling factors to restore original units
+    Ar, Br, Cr, Dr = ssdata(ss(sysr))
+    sys_final = ss(Ar, Br * scaleU, scaleY * Cr, scaleY * Dr * scaleU, sys.timeevol)
+    sys_final, hs
 end
 
 """
@@ -105,11 +112,15 @@ Coprime-factor reduction performs a coprime factorization of the model into \$P(
 - `kwargs`: Are passed to `DescriptorSystems.gbalmr`, the docstring of which is reproduced below:
 $(@doc(DescriptorSystems.gbalmr))
 """
-function baltrunc_coprime(sys, info=nothing; residual=false, n=missing, factorization::F = DescriptorSystems.gnlcf, kwargs...) where F
+function baltrunc_coprime(sys, info=nothing; residual=false, n=missing, factorization::F = DescriptorSystems.gnlcf, scaleY=1.0, scaleU=1.0, kwargs...) where F
+    # Apply scaling if needed
+    A, B, C, D = ssdata(sys)
+    # Divide by scaling factors to normalize to ~[-1,1]
+    sys_scaled = ss(A, B / scaleU, scaleY \ C, scaleY \ D / scaleU, sys.timeevol)
     if info !== nothing && hasproperty(info, :NM)
         @unpack N, M, NM = info
     else
-        N,M = factorization(dss(sys))
+        N,M = factorization(dss(sys_scaled))
         A,E,B,C,D = DescriptorSystems.dssdata(N)
         NM = DescriptorSystems.dss(A,E,[B M.B],C,[D M.D])
     end
@@ -128,7 +139,9 @@ function baltrunc_coprime(sys, info=nothing; residual=false, n=missing, factoriz
     Br = BN  - BM * (DMi * DN)
     Dr = (DMi * DN)
 
-    ss(Ar,Br,Cr,Dr,sys.timeevol), hs, (; NM, N, M, NMr)
+    # Multiply by scaling factors to restore original units
+    sys_final = ss(Ar, Br * scaleU, scaleY * Cr, scaleY * Dr * scaleU, sys.timeevol)
+    sys_final, hs, (; NM, N, M, NMr)
 end
 
 
@@ -139,18 +152,25 @@ Balanced truncation for unstable models. An additive decomposition of sys into `
 
 See `baltrunc2` for other keyword arguments.
 """
-function baltrunc_unstab(sys::LTISystem, info=nothing; residual=false, n=missing, kwargs...)
+function baltrunc_unstab(sys::LTISystem, info=nothing; residual=false, n=missing, scaleY=1.0, scaleU=1.0, kwargs...)
+    # Apply scaling if needed
+    A, B, C, D = ssdata(sys)
+    # Divide by scaling factors to normalize to ~[-1,1]
+    sys_scaled = ss(A, B / scaleU, scaleY \ C, scaleY \ D / scaleU, sys.timeevol)
     if info !== nothing && hasproperty(info, :stab)
         @unpack stab, unstab = info
     else
-        stab, unstab = DescriptorSystems.gsdec(dss(sys); job="stable", kwargs...)
+        stab, unstab = DescriptorSystems.gsdec(dss(sys_scaled); job="stable", kwargs...)
     end
     nx_unstab = size(unstab.A, 1)
     if n isa Integer && n < nx_unstab
         error("The model contains $(nx_unstab) poles outside the stability region, the reduced-order model must be of at least this order.")
     end
     sysr, hs = DescriptorSystems.gbalmr(stab; matchdc=residual, ord=n-nx_unstab, kwargs...)
-    ss(sysr + unstab), hs, (; stab, unstab)
+    # Multiply by scaling factors to restore original units
+    Ar, Br, Cr, Dr = ssdata(ss(sysr + unstab))
+    sys_final = ss(Ar, Br * scaleU, scaleY * Cr, scaleY * Dr * scaleU, sys.timeevol)
+    sys_final, hs, (; stab, unstab)
 end
 
 
diff --git a/test/test_reduction.jl b/test/test_reduction.jl
index 7d465af..fd24168 100644
--- a/test/test_reduction.jl
+++ b/test/test_reduction.jl
@@ -300,4 +300,37 @@ Ksr, hs, infor = baltrunc_coprime(info.Ks; n)
 controller_reduction_plot(info.Gs,info.Ks)
 controller_reduction_plot(info.Gs,info.Ks, method=:cr)
 
-# ncfmargin(P, W1*Ksr)
\ No newline at end of file
+# ncfmargin(P, W1*Ksr)
+
+
+
+# === Scaling invariance tests for model reduction ===
+using RobustAndOptimalControl: baltrunc2, baltrunc_coprime, baltrunc_unstab
+
+# Simple SISO system
+
+sys_siso = ssrand(1,1,9,proper=true)
+
+# Non-unit scaling factors
+scaleY = 5.0
+scaleU = 0.2
+
+# baltrunc2: compare frequency response with and without scaling
+sysr_noscale, _ = baltrunc2(sys_siso; n=3)
+sysr_scale, _ = baltrunc2(sys_siso; n=3, scaleY=scaleY, scaleU=scaleU)
+@test sysr_noscale.nx == sysr_scale.nx == 3
+@test hinfnorm2(minreal(sysr_noscale - sysr_scale))[1] < 1e-5
+
+
+# baltrunc_coprime: compare frequency response with and without scaling
+sysr_coprime_noscale, _, _ = baltrunc_coprime(sys_siso; n=3)
+sysr_coprime_scale, _, _ = baltrunc_coprime(sys_siso; n=3, scaleY=scaleY, scaleU=scaleU)
+@test sysr_coprime_noscale.nx == sysr_coprime_scale.nx == 3
+@test hinfnorm2(minreal(sysr_coprime_noscale - sysr_coprime_scale))[1] < 1e-5
+
+
+# baltrunc_unstab: compare frequency response with and without scaling
+sysr_unstab_noscale, _, _ = baltrunc_unstab(sys_siso; n=3)
+sysr_unstab_scale, _, _ = baltrunc_unstab(sys_siso; n=3, scaleY=scaleY, scaleU=scaleU)
+@test sysr_unstab_noscale.nx == sysr_unstab_scale.nx == 3
+@test hinfnorm2(minreal(sysr_unstab_noscale - sysr_unstab_scale))[1] < 1e-5