Adam Cooman

This page contains the full text of the accepted version of the paper published at the 2018 International Microwave Symposium. The content was modified slightly to better fit the format of this website

Estimating unstable poles
in simulations of microwave circuits

Adam Cooman, Fabien Seyfert and Smain Amari

The impedance of a microwave circuit has an infinite number of poles due to the distributed elements. This complicates locating those poles with a rational approximation. In this paper, we propose an algorithm to locate the unstable poles of a circuit with distributed elements. The proposed method exploits the fact that a realistic circuit can only have a finite number of unstable poles. We first determine the unstable part whose poles coincide with the unstable poles of the circuit. A rational approximation of the unstable part is used to estimate the unstable poles.

1. Introduction

Having the ability to trace a circuit’s poles as a function of the circuit parameters is a useful design tool. Pole tracking techniques have been used for the design of oscillators1 2, to stabilise power amplifiers3 4 5 and during the design of frequency dividers1 6. The core of a pole tracking tool is a robust automatic algorithm to estimate the poles of a circuit.

To determine the poles of a circuit, a two-step procedure is followed. First, an impedance Z ⁣(jω)Z\!\left(j\omega\right) of the circuit is determined at discrete frequency points between 00 and fmaxf_{\mathrm{max}} with an AC simulation (Fig. 1). Then, the poles of Z ⁣(jω)Z\!\left(j\omega\right) are determined in a post-processing step.

simulation setup used for a closed-loop stability analysis

Figure 1: Z ⁣(jω)Z\!\left(j\omega\right) is determined by connecting a small-signal current source IexI_{ex} in parallel to a node in the circuit and by determining the voltage response V ⁣(jω)V\!\left(j\omega\right) of the circuit with an AC simulation.

In lumped circuits, Z ⁣(jω)Z\!\left(j\omega\right) is a rational function, so a rational approximation can be used to determine the circuit poles. The impedance presented by a circuit with distributed elements, like transmission lines, is non-rational but still meromorphic 7. This makes estimating the poles of a microwave circuit more difficult. A good fit of Z ⁣(jω)Z\!\left(j\omega\right) can be obtained with a high-order rational function, but the rational approximation will contain spurious poles that do not correspond to poles of the underlying function 8. To circumvent this problem, another approach proposed in 9, is to compute low-order rational approximants of the circuit’s response restricted to small frequency intervals. This local version of the rational approximation scheme, yields precise estimates of poles when these are close enough to the imaginary axis.

In this paper, we propose an algorithm that can estimate the unstable poles of a circuit without performing a rational approximation of a non-rational function. The proposed technique exploits the fact that the equilibrium solution of a realistic circuit 10 can only have a finite amount of poles in the right half-plane 11.

When Z ⁣(jω)Z\!\left(j\omega\right) is square integrable, it can be split uniquely into two meromorphic parts: one part in H2H_{2}, the Hardy space of the complex right half-plane, and another part in H2\overline{H}_{2}, the orthogonal complement of H2H_{2}. The projection of Z ⁣(jω)Z\!\left(j\omega\right) onto a basis of H2\overline{H}_{2} is defined as its “unstable part” Zunstab ⁣(jω)Z_{\mathrm{unstab}}\!\left(j\omega\right) because its poles coincide with the unstable poles of Z ⁣(jω)Z\!\left(j\omega\right).

There are only a finite amount of unstable poles, so Zunstab ⁣(jω)Z_{\mathrm{unstab}}\!\left(j\omega\right) is rational. To reliably estimate the unstable poles of Z ⁣(jω)Z\!\left(j\omega\right), we propose to compute its unstable part and to approximate this unstable part with a rational function.

In12, it was shown that estimates of the unstable part can be obtained by computing the Fourier coefficients of a pre-filtered version of the impedance function. In this paper, the method is extended so as to yield estimates of the unstable poles. In order to do so:

  • A Hankel-operator-based rational approximation technique is presented (Section 2.3).
  • Fourier coefficients are obtained via a rational interpolation scheme that reduces the errors due to sampling (Section 2.4).
  • The spectral range of the filtering function is reduced to improve the procedure’s robustness (Section 2.5).

After detailing the approach, we run it on two examples: a time-delayed Chua’s circuit (Section 3) and on simulations of a microwave oscillator (Section 4).

2. Algorithm

The response Z ⁣(jω)Z\!\left(j\omega\right) is first transformed to the unit circle (Section 2.1) to allow for a simple determination of its unstable part by means of its Fourier series (Section 2.2 & 2.4). With its Fourier expansion, the unstable poles can be recovered by means of a singular value decomposition of a truncated version of the system’s Hankel matrix (Section 2.3). These steps assume that Z(jω)Z \left(j\omega\right) is known everywhere while it is only known on a finite interval. To deal with finite data, the Fourier coefficients are convolved by the Fourier coefficients of the filtering function (Section 2.5).

2.1. Transform to the unit disc

A Mobius transform is used to transform Z ⁣(jω)Z\!\left(j\omega\right) from the complex plane to the unit circle:

Z˘(θ)=πα2ejθ1Z(jαcot(θ2))(1)\breve Z\left(\theta\right)=\sqrt{\pi\alpha}\frac{2}{e^{j\theta}-1}Z\left(j\alpha\cot\left(\frac{\theta}{2}\right)\right)\tag{1}

where α=2πfmax(1+2)\alpha=\frac{2\pi f_{\mathrm{max}}}{\left(1+\sqrt{2}\right)}. This transform maps ω=0\omega=0 to θ=π\theta=\pi and fmaxf_{\mathrm{max}} to θ=π4\theta=\frac{\pi}{4}. The complex right half-plane is mapped to the inside of the unit circle (Fig. 2). We will use ˘\breve \cdot to indicate functions on the unit disc.

2.2 Determine the unstable part

On the unit circle, a basis for H2\overline{H}_{2} is zk,kNz^{-k},k\in\mathbb{N} 12 so we compute the Fourier series of Z˘(θ)\breve Z\left(\theta\right) to obtain its unstable part

Z˘(θ)=k=+fkejkθ=k=1fkejkθZ˘unstab+k=0+fkejkθZ˘stab\breve Z\left(\theta\right)=\sum_{k=-\infty}^{+\infty}f_{k} e^{jk\theta}=\underbrace{\sum_{k=-\infty}^{-1}f_{k} e^{jk\theta}}_{\breve Z_{\mathrm{unstab}}}+\underbrace{\sum_{k=0}^{+\infty}f_{k} e^{jk\theta}}_{\breve Z_{\mathrm{stab}}}

In 12, the Fast Fourier Transform (FFT) is used to compute the fkf_{k}. This introduces an aliasing error in the method. To avoid aliasing, we compute the Fourier coefficients with the help of an interpolation function as explained in Section 2.4.

2.3. Estimating unstable poles

We use the theory introduced by Adamyan, Arov and Krein 13 to recover the location of the unstable poles on the unit disc. It uses the following infinite Hankel matrix which contains the negative Fourier coefficients of Z˘(θ)\breve Z\left(\theta\right):

Φ=[f1f2f2](2)\Phi_{\,}=\left[\begin{array}{ccc} f_{-1} & f_{-2} & \cdots\\ f_{-2} & \ddots\\ \vdots \end{array}\right]\tag{2}

In practice, we truncate the infinite Φ\Phi_{\,} to ΦN=Φ1..N,1..N\Phi_{N}=\Phi_{1..N,1..N}. N>50N>50 has yielded good results in all examples we tried. By Kronecker’s theorem for a rational symbol with PP poles, ΦN\Phi_{N}, is of rank PP. We compute the rank using the Singular Value Decomposition of ΦN\Phi_{N}

ΦN=USW\Phi_{N}=USW^{\prime}

Ideally the number of non-zero singular values of ΦN\Phi_{N} equals the number of poles of the function Zunstab ⁣(jω)Z_{\mathrm{unstab}}\!\left(j\omega\right). In practice due to numerical errors, the first PP singular values are large when compared to the others. This property is used to determine the number of unstable poles PP of ZZ.

Once PP has been determined, we can determine the circuit poles using the Principal Hankel Components (PHC) algorithm 14: We split the matrices UU, SS and WW into parts

ΦN=[UPUR][SP00SR][WPWR]\Phi_{N}=\left[\begin{array}{cc} U_{P} & U_{R}\end{array}\right]\left[\begin{array}{cc} S_{P} & 0\\ 0 & S_{R} \end{array}\right]\left[\begin{array}{c} W_{P}^{\prime}\\ W_{R}^{\prime} \end{array}\right]

such that SPS_{P} is a diagonal matrix of size P×PP\times P which contains the singular values, UPU_{P} is of size N×PN\times P and WPW_{P}^{\prime} is size P×NP\times N. With UpU_{p} and SPS_{P}, we compute an estimate of the observability matrix of the approximation:

O=UPSP12\mathcal{O}=U_{P}S_{P}^{\frac{1}{2}}

From the observability matrix O\mathcal{O}, we can compute the state matrix AA using

A=(OlOl)1OlOf\begin{align} A & =\left(\mathcal{O}_{l}^{\prime}\mathcal{O}_{l}\right)^{-1}\mathcal{O}_{l}^{\prime}\mathcal{O}_{f}\tag{3} \end{align}

where Ol\mathcal{O}_{l} and Of\mathcal{O}_{f} are formed from O\mathcal{O} by deleting the last and first row respectively. \cdot^{\prime} indicates the complex conjugate. The eigenvalues of AA correspond to the location of the unstable poles.

2.4. Computing the Fourier coefficients

The Fourier coefficients needed to construct ΦN\Phi_{N} are given by:

fk=12π02πZ˘(θ)ejkθdθ.(4)f_{k}=\frac{1}{2\pi}\int_{0}^{2\pi}\breve Z\left(\theta\right) e^{-jk\theta}d\theta.\tag{4}

The function Z˘(θ)\breve Z\left(\theta\right) is known on a set of MM sample points θm\theta_{m}. This calls for the use of an interpolation procedure for a precise evaluation of (4). When the impedance Z˘(θ)\breve Z\left(\theta\right) contains strong resonances, classical spline interpolation fails to model the response’s stiffness. We therefore introduce a rational interpolation scheme by defining I˘m(θ)\breve I_{m}\left(\theta\right) between the points θm\theta_{m} and θm+1\theta_{m+1} as:

I˘m(θ)=am(ejθ)2+bmejθ+cmejθ+dm(5)\breve I_{m}\left(\theta\right)=\frac{a_{m}\left(e^{j\theta}\right)^{2}+b_{m}e^{j\theta}+c_{m}}{e^{j\theta}+d_{m}}\tag{5}

where ama_{m}, bmb_{m}, cmc_{m} and dmd_{m} are determined using the constraints:

I˘m(θm)=Z˘(θm)I˘m(θm+1)=Z˘(θm+1)I˘m(θm+2)=Z˘(θm+2)θI˘m(θm)=θI˘m1(θm)\begin{array}{cc} \breve I_{m}\left(\theta_{m}\right)=\breve Z\left(\theta_{m}\right) & \breve I_{m}\left(\theta_{m+1}\right)=\breve Z\left(\theta_{m+1}\right)\\ \breve I_{m}\left(\theta_{m+2}\right)=\breve Z\left(\theta_{m+2}\right) & \frac{\partial}{\partial\theta}\breve I_{m}\left(\theta_{m}\right)=\frac{\partial}{\partial\theta}\breve I_{m-1}\left(\theta_{m}\right) \end{array}

The final constraint ensures that the interpolated function has a continuous first derivative at the interpolation points.

When the rational interpolation models I˘m(θ)\breve I_{m}\left(\theta\right) are known, the integral (4) is computed with Matlab’s built-in quadrature method 15.

2.5. Filtering to suppress edge effects

The simulated impedance is only known for frequencies below fmaxf_{\mathrm{max}}. Above that frequency, we assume that Z ⁣(jω)=0Z\!\left(j\omega\right)=0. This sudden jump to zero at fmaxf_{\mathrm{max}} will have a strong unwanted influence on the obtained stable and unstable part 12. To suppress this influence, Z˘(θ)\breve Z\left(\theta\right) is multiplied by a low-pass filter with its stop-band starting at fmaxf_{\mathrm{max}}:

Z˘filt=Z˘(θ)H˘(θ)\breve Z_{\mathrm{filt}}=\breve Z\left(\theta\right)\cdot \breve H\left(\theta\right)

In 12, an elliptic filter on the complex plane is used. On the unit disc, this elliptic filter has an infinite number of Fourier coefficients, which increases the aliasing error in the method. Here, we use a FIR filter with a finitely supported spectrum, and defined by its NHN_{H} Fourier coefficients of positive index:

H˘(θ)=k=0NHhkejkθ\breve H\left(\theta\right)=\sum_{k=0}^{N_{H}}h_{k}e^{jk\theta}

The coefficients hkh_{k} are computed with a Remes algorithm to obtain a filter with a pass-band ripple of 0.020.02 and a stop-band rejection of 85dB85\mathrm{dB} without zeroes inside the unit circle (Fig. 3). The filter has double transmission zeroes at θ=±π4\theta=\pm\frac{\pi}{4} to guarantee a smooth transition to zero at the edge of the interval. Good results are obtained when using a filter of order NH=72N_{H}=72. With the hkh_{k} known, the Fourier coefficients of Z˘filt\breve Z_{\mathrm{filt}} are obtained by convolution:

fkfilt=l=0NHfklhl(6)f_{k}^{\mathrm{filt}}=\sum_{l=0}^{N_{H}}f_{k-l}h_{l}\tag{6}

The fkfiltf_{k}^{\mathrm{filt}} are then used in ΦN\Phi_{N} to compute the unstable poles. With the elliptic filter used in 12, the sum in the convolution is infinite. With the new filter, this sum becomes finite. We only need Fourier coefficients fkNHf_{k-N_{H}} to fkf_{k} to determine fkfiltf_{k}^{\mathrm{filt}} exactly. This is the second advantage of the new filter.

mobius mapping used to transform from plae to the disc

Figure 2: Mapping from plane to disc

filtering function used in the analysis

Figure 3: Filtering function H˘(θ)\breve H\left(\theta\right)

2.6. Summary

Unstable poles of Z ⁣(jω)Z\!\left(j\omega\right) are obtained with the following steps:

  1. Transform Z ⁣(jω)Z\!\left(j\omega\right) to the unit circle (1)
  2. Compute the local rational interpolation models MiM_{i} (5)
  3. Compute the Fourier coefficients of Z˘(θ)\breve Z\left(\theta\right) (4)
  4. Use convolution to apply the filtering function (6)
  5. Construct ΦN\Phi_{N} and compute its singular values to determine the number of unstable poles
  6. Use the PHC algorithm to compute the location of the unstable poles

3. Time-delayed Chua’s circuit

As a first example, we determine the unstable poles of the equilibrium solution of the time-delayed Chua’s circuit shown in Fig. 416. The impedance presented by the circuit was determined on 1500 points between 0 and 9GHz9\,\mathrm{GHz} (Fig. 5). A finer grid of points was placed around the strong resonances in the impedance. The singular values of Φ50\Phi_{50} indicate that there are four unstable poles in the circuit (Fig. 5). Using the PHC algorithm, we estimate the following unstable poles:

p1,2=2π(0.00017±j2.0014)109p_{1,2}=2\pi\left(0.00017\pm j2.0014\right)\cdot10^{9} p3=2π(6.316)106p4=2π(447.22)106p_{3}=2\pi\left(6.316\right)\cdot10^{6}\qquad p_{4}=2\pi\left(447.22\right)\cdot10^{6}

The method finds the unstable pole pair around 2GHz2\mathrm{GHz}, but it also finds two unstable poles on the real axis. To verify this result, we computed the exact location of the circuit poles using its analytical expression (Appendix 1). The poles obtained with the stability analysis lie very close to the exact poles. The maximum normalised error we obtain is

error=maxipi,estimatedpi,correctpi,correct=5.5106\mathrm{error}=\max_{i}\frac{\left|p_{i,\mathrm{estimated}}-p_{i,\mathrm{correct}}\right|}{\left|p_{i,\mathrm{correct}}\right|}=5.5\cdot10^{-6}
Chua oscillator circuit

Figure 4: Time-delayed Chua’s circuit analysed here. The impedance is determined at the location of the red dot.

result of the stability analysis applied to the chua oscillator

Figure 5: (left) Impedance presented by the time-delayed Chua’s circuit. (middle) zoom on the impedance on a 10MHz span around 2GHz (right) the first 10 singular values of Φ50\Phi_{50}.

4. MMIC oscillator

As a more realistic example, we estimate the unstable poles of an MMIC oscillator from the example library of ADS Fig. 6. Z ⁣(jω)Z\!\left(j\omega\right) is determined by co-simulation. Lumped models were used for the passives and transistors in the oscillator, while a Momentum simulation was used for the transmission lines. Z ⁣(jω)Z\!\left(j\omega\right) is determined at the source of the oscillating transistor on 600600 points between 1MHz1\mathrm{MHz} and 100GHz100\mathrm{GHz} (Fig. 7). Because Z ⁣(jω)Z\!\left(j\omega\right) is not available close to 0Hz0\mathrm{Hz}, we first use a bandpass transformation on the data. For this circuit, Φ50\Phi_{50} has only one significant singular value, which indicates the presence of one unstable pole in the circuit. Its estimated location is at s=2π(0.79+j20.94)109s=2\pi\left(0.79+j20.94\right)10^{9}. Here, the exact location of the unstable pole in the equilibrium solution is unknown, so we cannot compute the error. The oscillating solution found with Harmonic Balance has a frequency of 20.42GHz20.42\mathrm{GHz}, which gives an indication that the estimated unstable pole is correctly estimated.

layout of the MMIC oscillator used in the analysis

Figure 6: Lay-out of the MMIC oscillator. The circuit impedance is determined at the location indicated by the arrow.

result of the analysis when applied to the mmic oscillator

Figure 7: (left) Impedance of the oscillator and (right) the first five normalised singular values of Φ50\Phi_{50}

5. Conclusion

In this paper, we have introduced an algorithm to estimate the unstable poles of the equilibrium solution of a microwave circuit. This is no easy task, because the impedance Z ⁣(jω)Z\!\left(j\omega\right) of a circuit with distributed elements is not rational. The impedance however only has a finite number of unstable poles: computing its unstable part via a projection step yields a rational function, the unstable poles of which are recovered by the proposed procedure. The algorithm, which computes the singular values of a Hankel matrix based on the Fourier coefficients of the filtered impedance function is shown to yield very accurate pole estimates on two examples.

The presented technique only works to estimate the unstable poles of the equilibrium solution of a circuit. The stable poles remain infinite in number, so we cannot exploit the rationality of the stable part. The same problem is encountered in periodic solutions where both the stable and unstable part have an infinite amount of poles. Extending the pole estimating technique to these cases remains a challenging task.

Appendix 1. Exact poles of the time-delayed Chua’s circuit

We have the following analytical expression for the circuit impedance at the node indicated in red in Fig. 4:

Z(s)=(n11e2sτ+n10)s+n01e2sτ+n00(d11e2sτ+d10)s+d01e2sτ+d00Z\left(s\right)=\frac{\left(n_{11}\mathrm{e}^{2s\tau}+n_{10}\right)\,s+n_{01}\mathrm{e}^{2s\tau}+n_{00}}{\left(d_{11}\mathrm{e}^{2s\tau}+d_{10}\right)\,s+d_{01}\mathrm{e}^{2s\tau}+d_{00}}

with

n11=CPRPZ0R(RS+Z0)n10=CPRPZ0R(RSZ0)n01=Z0(RP+R)(RS+Z0)n00=Z0(RP+R)(RSZ0)d11=CPR(Z0(RP+RS+Z0)+RPRS)d10=CPR(Z0(RP+RSZ0)RPRS)d01=Z0(R+RP+RS+Z0)+RS(R+RP)d00=Z0(R+RP+RSZ0)RS(R+RP)\begin{array}{rl} n_{11} & =C_{P}R_{P}Z_{0}R\left(R_{S}+Z_{0}\right)\\ n_{10} & =C_{P}R_{P}Z_{0}R\left(R_{S}-Z_{0}\right)\\ n_{01} & =Z_{0}\left(R_{P}+R\right)\left(R_{S}+Z_{0}\right)\\ n_{00} & =Z_{0}\left(R_{P}+R\right)\left(R_{S}-Z_{0}\right)\\ d_{11} & =C_{P}R\left(Z_{0}\left(R_{P}+R_{S}+Z_{0}\right)+R_{P}R_{S}\right)\\ d_{10} & =C_{P}R\left(Z_{0}\left(R_{P}+R_{S}-Z_{0}\right)-R_{P}R_{S}\right)\\ d_{01} & =Z_{0}\left(R+R_{P}+R_{S}+Z_{0}\right)+R_{S}\left(R+R_{P}\right)\\ d_{00} & =Z_{0}\left(R+R_{P}+R_{S}-Z_{0}\right)-R_{S}\left(R+R_{P}\right) \end{array}

To determine the amount of unstable poles, we use the fact that the roots of the denominator quasipolynomial vary continuously with respect to the delay τ\tau 17. For τ=0\tau=0, the denominator is rational and has its root in the left half-plane. Because the circuit is realistic, its denominator will have all roots in the left half-plane for τ=0+ε\tau=0+\varepsilon 11. So, to determine the amount of roots in the complex right half-plane at τ=0.25ns\tau=0.25\mathrm{ns}, we compute the amount of zero pairs that cross the imaginary axis in the interval [0+ε,τ][0+\varepsilon,\tau].

A value τ\tau for which a zero lies on the imaginary axis is called a critical delay. To compute the critical delays we construct the two variable polynomial p(s,z)p\left(s,z\right) by replacing z=esτz=e^{-s\tau} in the denominator quasipolynomial. The solutions of p(s,z)=0p\left(s,z\right)=0 where ss lies on the imaginary axis and zz on the unit circle are called critical pairs, and can be found by solving a set of non-linear polynomial equations 17. In this case, p(s,z)p\left(s,z\right) has one critical pair. Its first three corresponding critical delays are τ=0.077ps\tau=0.077\mathrm{ps}, 0.16ns0.16\mathrm{ns} and 0.33ns0.33\mathrm{ns}. Two zero pairs have crossed the imaginary axis for our delay of τ=0.25ns\tau=0.25\,\mathrm{ns}, so the denominator of the circuit has 4 unstable zeroes.

We observed two zeroes of the denominator quasipolynomial on the real axis, so we determined their location using dichotomy (fzero in Matlab). The location of the complex pole pair, close to 2GHz2\mathrm{GHz}, was determined using Matlab’s fsolve, with its starting point at s=j2π2109s=j2\pi 2\cdot 10^{9}.

Acknowledgements

We would like to thank department ELEC of the VUB, who provided the simulation tools for this paper through their strategic research program SRP-19.

Footnotes

  1. A. Suarez, Analysis and Design of Autonomous Microwave Circuits. John Wiley and Sons, 2009. 2

  2. F. Ramirez, M. Ponton, S. Sancho, and A. Suarez, “Stability analysis of oscillation modes in quadruple-push and rucker’s oscillators,” IEEE Tran. on Microwave Theory and Techniques, vol. 56, no. 11, Nov 2008.

  3. N. Ayllon, J. M. Collantes, A. Anakabe, I. Lizarraga, G. Soubercaze-Pun, and S. Forestier, “Systematic approach to the stabilization of multitransistor circuits,” IEEE Tran. on Microwave Theory and Techniques, vol. 59, no. 8, pp. 2073-2082, 2011.

  4. M. van Heijningen, A. P. de Hek, F. E. van Vliet, and S. Dellier, “Stability analysis and demonstration of an x-band gan power amplifier mmic,” in 2016 11th European Microwave Integrated Circuits Conference (EuMIC), Oct 2016, pp. 221-224.

  5. S. Jeon, A. Suarez, and D. B. Rutledge, “Global stability analysis and stabilization of a class-e/f amplifier with a distributed active transformer,” IEEE Tran. on Microwave Theory and Techniques, vol. 53, no. 12, pp. 3712-3722, Dec 2005.

  6. L. Pantoli, A. Suarez, G. Leuzzi, and F. D. Paolo, “Complete and systematic simulation tools for frequency divider design,” IEEE Tran. on Microwave Theory and Techniques, vol. 56, no. 11, Nov 2008.

  7. H. Carlin and P. P. Civalleri, Wideband circuit design. CRC Press, 1997.

  8. H. Stahl, “Spurious poles in padé approximation,” Journal of Computational and Applied Mathematics, no. 99, 1998.

  9. J. Jugo, J. Portilla, A. Anakabe, A. Suarez, and J. Collantes, “Closed-loop stability analysis of microwave amplifiers,” Electronics Letters, vol. 37, no. 4, pp. 226-228, Feb 2001.

  10. A circuit is considered realistic when its passives are lossy and when all active devices have a finite gain-bandwidth.

  11. L. Baratchart, S. Chevillard, and F. Seyfert. (2014) On transfer functions realizable with active electronic components. hal-01098616. Inria Sophia Antipolis. [Research Report] RR-8659. 2

  12. A. Cooman, F. Seyfert, M. Olivi, S. Chevillard, and L. Baratchart, “Modelfree closed-loop stability analysis: A linear functional approach,” IEEE Tran. on Microwave Theory and Techniques, no. 99, 2017. 2 3 4 5 6

  13. K. Glover, “All optimal hankel-norm approximations of linear multivariable systems and their H∞ error bounds,” Internat. J. Control, vol. 39, no. 6, pp. 1115-1193, 1984.

  14. S. Y. Kung and K. S. Arun, Advances in Statistical Signal Processing. JAJ Press, 1987, ch. 6. Singular value decomposition algorithms for linear system approximation and spectrum estimation.

  15. L. Shampine, “Vectorized adaptive quadrature in matlab,” Journal of Computational and Applied Mathematics, no. 211, pp. 131-140, 2008.

  16. M. Biey, F. Bonani, M. Gilli, and I. Maio, “Qualitative analysis of the dynamics of the time-delayed chua’s circuit,” IEEE Tran. on Circuits and Systems, vol. 44, no. 6, pp. 486-500, Jun 1997.

  17. X.-G. Li, S.-I. Niculescu, and A. Cela, Analytic Curve Frequency Sweeping Stability Tests for Systems with Commensurate Delays. Springer, 2015 2