Since you have an input-output system, I’m not sure the package you link is the most suitable. See Home · RobustAndOptimalControl Documentation for model-order reduction of linear statespace systems with inputs, in particular, try
using ControlSystemsBase
sys = ss(matrices...)
baltrunc(sys, residual=true)
baltrunc(sys, n=desired_state_dimension, residual=true)
or
using RobustAndOptimalControl
using DSP
w1 = 1e-4 # Lower frequency limit
w2 = 1e1 # upper
fc = DSP.analogfilter(DSP.Bandpass(w1, w2), DSP.Butterworth(2))
tfc = DSP.PolynomialRatio(fc)
W = tf(DSP.coefb(tfc), DSP.coefa(tfc))
rsys, hs = frequency_weighted_reduction(sys, W, 1)
for a reduction that focuses on a particular frequency range