## 2 Problem Formulation  
In this paper, MIMO discrete time linear time-invariant systems are considered. Usually, the underlying phenomena that are examined are inherently continuous systems, but they can be sufficiently approximated with discrete time systems, and all the discussion is carried out with the assumption that the real system is modeled with a discrete system using appropriate sampling. Formally, let an element of the transfer function matrix be a generic transfer function [22]:  
\[H(z,\zeta)=\frac{\sum\limits_{k=0}^{m}b_{k}(\zeta)z^{-k}}{1+\sum\limits_{k=1}^{ n}a_{k}(\zeta)z^{-k}} \tag{1}\]  
where  
* \(\zeta\in\mathbb{R}\) is the parameter;
* \(H(z,\zeta)\) is the transfer function at a certain \(\zeta\) value;* \(z\) is the complex variable of the function;
* \(\{b_{k}(\zeta)\}\) are the numerator coefficients, where \(b_{k}(\zeta)\ \in\mathbb{R}\);
* \(\{a_{k}(\zeta)\}\) are the denominator coefficients, where \(a_{k}(\zeta)\in\mathbb{R}\);
* \(m\) is the order of the numerator;
* \(n\) is the order of the denominator. The following assumptions are made on \(H(z,\zeta)\):  
* The transfer function remains stable for all \(\zeta\) ;
* It does not have any repeated poles;
* \(m<n\), ensuring a proper transfer function;
* \(\{b_{k}(\zeta)\}\) and \(\{a_{k}(\zeta)\}\) are continuous function of \(\zeta\);
* The sampling frequency is chosen to be high enough to ensure that all poles are mapped to the right-hand side of the complex plane while still remaining within the unit circle. It is further assumed that for an element of the transfer function matrix that is \(H(z,\zeta)\), there is a partial fractional expansion,  
\[H(z,\zeta)=\frac{r_{1}(\zeta)}{z-p_{1}(\zeta)}+\frac{r_{2}(\zeta)}{z-p_{2}( \zeta)}+\ldots+\frac{r_{n}(\zeta)}{z-p_{n}(\zeta)} \tag{2}\]  
where \(p_{i}(\zeta)\) and \(r_{i}(\zeta)\), for \(i=1,\ldots,n\), represent the poles and residues. The model \(H(z,\zeta)\) is assumed to be difficult or impractical to evaluate for arbitrary values of \(\zeta\). Instead, it is only available at \(k\) specific points, given by \(G_{i}=H(z,\zeta_{i})\) for \(i=1,\ldots,k\). Consequently, the behavior of \(H(z,\zeta)\) remains unknown for any \(\zeta\) other than these given points. Given that \(\zeta\) is a scalar parameter, it is important to note that the models have a natural ordering. For instance, one can identify adjacent known models of \(H(z,\zeta)\ \zeta_{i}\leq\zeta\leq\zeta_{i+1}\) as \(G_{i}\) and \(G_{i+1}\). This ordering will be utilized in the following discussion. This paper presents an interpolation algorithm that preserves stability and ensures that the deviation of the interpolated model from the known adjacent models remains within a guaranteed upper bound in terms of the \(H_{\infty}\) norm. Formally, there exists a function \(f\) associated with the proposed interpolation method such that the following holds:  
\[\left\|G(\zeta)-G_{i}\right\|_{\infty}\leq f(G_{i+1},G_{i})\] \[\left\|G(\zeta)-G_{i+1}\right\|_{\infty}\leq f(G_{i+1},G_{i}), \quad\zeta_{i}\leq\zeta\leq\zeta_{i+1} \tag{3}\]  
where \(G(\zeta)\) is the interpolated model. The derivation of \(f\) is presented in Section 4. The proposed method is beneficial when it is known a priori that the given system remains stable for all parameter values and that the unknown function \(H(z,\zeta)\), within the range \(\zeta_{i}\leq\zeta\leq\zeta_{i+1}\), exhibits a behavior \"similar\" to \(G_{i}\) and \(G_{i+1}\). Such conditions frequently arise in practical applications. For instance, while the dynamics of an airplane vary with airspeed, it is known that the system remains stable within a certain speed range, and these dynamics change gradually. Similarly, geometrical modifications in a mechanical component often (though not always) result in predictable changes in its dynamic behavior. If the dynamics of the system are known for specific configurations, physical considerations suggest that other configurations will not exhibit arbitrary deviations.