Nowadays, modeling dynamical systems often results in high-order systems (i.e. with 10,000 or more equations). In order to ensure a numerical simulation within a reasonable time, the given dynamical system of equations is reduced to a system of the same form, which allows a solution with a greatly reduced computational time. Often it is required that the reduced system has the same properties as the unreduced model; important properties in this context are especially stability and passivity. Furthermore, the approximation error should be as small as possible. Above all, the model reduction methods should be numerically stable and efficient and ideally end automatically with a given error tolerance.