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 oscillators^{1} ^{2}, to stabilise power amplifiers^{3} ^{4} ^{5} and during the design of frequency dividers^{1} ^{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 twostep procedure is followed. First, an impedance $Z\!\left(j\omega\right)$ of the circuit is determined at discrete frequency points between $0$ and $f_{\mathrm{max}}$ with an AC simulation (Fig. 1). Then, the poles of $Z\!\left(j\omega\right)$ are determined in a postprocessing step.
Figure 1: $Z\!\left(j\omega\right)$ is determined by connecting a smallsignal current source $I_{ex}$ in parallel to a node in the circuit and by determining the voltage response $V\!\left(j\omega\right)$ of the circuit with an AC simulation.
In lumped circuits, $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 nonrational but still meromorphic ^{7}. This makes estimating the poles of a microwave circuit more difficult. A good fit of $Z\!\left(j\omega\right)$ can be obtained with a highorder 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 loworder 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 nonrational 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 halfplane ^{11}.
When $Z\!\left(j\omega\right)$ is square integrable, it can be split uniquely into two meromorphic parts: one part in $H_{2}$, the Hardy space of the complex right halfplane, and another part in $\overline{H}_{2}$, the orthogonal complement of $H_{2}$. The projection of $Z\!\left(j\omega\right)$ onto a basis of $\overline{H}_{2}$ is defined as its “unstable part” $Z_{\mathrm{unstab}}\!\left(j\omega\right)$ because its poles coincide with the unstable poles of $Z\!\left(j\omega\right)$.
There are only a finite amount of unstable poles, so $Z_{\mathrm{unstab}}\!\left(j\omega\right)$ is rational. To reliably estimate the unstable poles of $Z\!\left(j\omega\right)$, we propose to compute its unstable part and to approximate this unstable part with a rational function.
In^{12}, it was shown that estimates of the unstable part can be obtained by computing the Fourier coefficients of a prefiltered 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 Hankeloperatorbased 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 timedelayed Chua’s circuit (Section 3) and on simulations of a microwave oscillator (Section 4).
2. Algorithm
The response $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 \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\!\left(j\omega\right)$ from the complex plane to the unit circle:
$\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 $\alpha=\frac{2\pi f_{\mathrm{max}}}{\left(1+\sqrt{2}\right)}$. This transform maps $\omega=0$ to $\theta=\pi$ and $f_{\mathrm{max}}$ to $\theta=\frac{\pi}{4}$. The complex right halfplane 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 $\overline{H}_{2}$ is $z^{k},k\in\mathbb{N}$ ^{12} so we compute the Fourier series of $\breve Z\left(\theta\right)$ to obtain its unstable part
$\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 $f_{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 $\breve Z\left(\theta\right)$:
$\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 $\Phi_{N}=\Phi_{1..N,1..N}$. $N>50$ has yielded good results in all examples we tried. By Kronecker’s theorem for a rational symbol with $P$ poles, $\Phi_{N}$, is of rank $P$. We compute the rank using the Singular Value Decomposition of $\Phi_{N}$
$\Phi_{N}=USW^{\prime}$Ideally the number of nonzero singular values of $\Phi_{N}$ equals the number of poles of the function $Z_{\mathrm{unstab}}\!\left(j\omega\right)$. In practice due to numerical errors, the first $P$ singular values are large when compared to the others. This property is used to determine the number of unstable poles $P$ of $Z$.
Once $P$ has been determined, we can determine the circuit poles using the Principal Hankel Components (PHC) algorithm ^{14}: We split the matrices $U$, $S$ and $W$ into parts
$\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 $S_{P}$ is a diagonal matrix of size $P\times P$ which contains the singular values, $U_{P}$ is of size $N\times P$ and $W_{P}^{\prime}$ is size $P\times N$. With $U_{p}$ and $S_{P}$, we compute an estimate of the observability matrix of the approximation:
$\mathcal{O}=U_{P}S_{P}^{\frac{1}{2}}$From the observability matrix $\mathcal{O}$, we can compute the state matrix $A$ using
$\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 $\mathcal{O}_{l}$ and $\mathcal{O}_{f}$ are formed from $\mathcal{O}$ by deleting the last and first row respectively. $\cdot^{\prime}$ indicates the complex conjugate. The eigenvalues of $A$ correspond to the location of the unstable poles.
2.4. Computing the Fourier coefficients
The Fourier coefficients needed to construct $\Phi_{N}$ are given by:
$f_{k}=\frac{1}{2\pi}\int_{0}^{2\pi}\breve Z\left(\theta\right) e^{jk\theta}d\theta.\tag{4}$The function $\breve Z\left(\theta\right)$ is known on a set of $M$ sample points $\theta_{m}$. This calls for the use of an interpolation procedure for a precise evaluation of (4). When the impedance $\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 $\breve I_{m}\left(\theta\right)$ between the points $\theta_{m}$ and $\theta_{m+1}$ as:
$\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 $a_{m}$, $b_{m}$, $c_{m}$ and $d_{m}$ are determined using the constraints:
$\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_{m1}\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 $\breve I_{m}\left(\theta\right)$ are known, the integral (4) is computed with Matlab’s builtin quadrature method ^{15}.
2.5. Filtering to suppress edge effects
The simulated impedance is only known for frequencies below $f_{\mathrm{max}}$. Above that frequency, we assume that $Z\!\left(j\omega\right)=0$. This sudden jump to zero at $f_{\mathrm{max}}$ will have a strong unwanted influence on the obtained stable and unstable part ^{12}. To suppress this influence, $\breve Z\left(\theta\right)$ is multiplied by a lowpass filter with its stopband starting at $f_{\mathrm{max}}$:
$\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 $N_{H}$ Fourier coefficients of positive index:
$\breve H\left(\theta\right)=\sum_{k=0}^{N_{H}}h_{k}e^{jk\theta}$The coefficients $h_{k}$ are computed with a Remes algorithm to obtain a filter with a passband ripple of $0.02$ and a stopband rejection of $85\mathrm{dB}$ without zeroes inside the unit circle (Fig. 3). The filter has double transmission zeroes at $\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 $N_{H}=72$. With the $h_{k}$ known, the Fourier coefficients of $\breve Z_{\mathrm{filt}}$ are obtained by convolution:
$f_{k}^{\mathrm{filt}}=\sum_{l=0}^{N_{H}}f_{kl}h_{l}\tag{6}$The $f_{k}^{\mathrm{filt}}$ are then used in $\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 $f_{kN_{H}}$ to $f_{k}$ to determine $f_{k}^{\mathrm{filt}}$ exactly. This is the second advantage of the new filter.
Figure 2: Mapping from plane to disc  Figure 3: Filtering function $\breve H\left(\theta\right)$ 
2.6. Summary
Unstable poles of $Z\!\left(j\omega\right)$ are obtained with the following steps:
 Transform $Z\!\left(j\omega\right)$ to the unit circle (1)
 Compute the local rational interpolation models $M_{i}$ (5)
 Compute the Fourier coefficients of $\breve Z\left(\theta\right)$ (4)
 Use convolution to apply the filtering function (6)
 Construct $\Phi_{N}$ and compute its singular values to determine the number of unstable poles
 Use the PHC algorithm to compute the location of the unstable poles
3. Timedelayed Chua’s circuit
As a first example, we determine the unstable poles of the equilibrium solution of the timedelayed Chua’s circuit shown in Fig. 4^{16}. The impedance presented by the circuit was determined on 1500 points between 0 and $9\,\mathrm{GHz}$ (Fig. 5). A finer grid of points was placed around the strong resonances in the impedance. The singular values of $\Phi_{50}$ indicate that there are four unstable poles in the circuit (Fig. 5). Using the PHC algorithm, we estimate the following unstable poles:
$p_{1,2}=2\pi\left(0.00017\pm j2.0014\right)\cdot10^{9}$ $p_{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 $2\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
$\mathrm{error}=\max_{i}\frac{\leftp_{i,\mathrm{estimated}}p_{i,\mathrm{correct}}\right}{\leftp_{i,\mathrm{correct}}\right}=5.5\cdot10^{6}$Figure 4: Timedelayed Chua’s circuit analysed here. The impedance is determined at the location of the red dot.
Figure 5: (left) Impedance presented by the timedelayed Chua’s circuit. (middle) zoom on the impedance on a 10MHz span around 2GHz (right) the first 10 singular values of $\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\!\left(j\omega\right)$ is determined by cosimulation. Lumped models were used for the passives and transistors in the oscillator, while a Momentum simulation was used for the transmission lines. $Z\!\left(j\omega\right)$ is determined at the source of the oscillating transistor on $600$ points between $1\mathrm{MHz}$ and $100\mathrm{GHz}$ (Fig. 7). Because $Z\!\left(j\omega\right)$ is not available close to $0\mathrm{Hz}$, we first use a bandpass transformation on the data. For this circuit, $\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\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.42\mathrm{GHz}$, which gives an indication that the estimated unstable pole is correctly estimated.
Figure 6: Layout of the MMIC oscillator. The circuit impedance is determined at the location indicated by the arrow.
Figure 7: (left) Impedance of the oscillator and (right) the first five normalised singular values of $\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\!\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 timedelayed Chua’s circuit
We have the following analytical expression for the circuit impedance at the node indicated in red in Fig. 4:
$Z\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
$\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 $\tau=0$, the denominator is rational and has its root in the left halfplane. Because the circuit is realistic, its denominator will have all roots in the left halfplane for $\tau=0+\varepsilon$ ^{11}. So, to determine the amount of roots in the complex right halfplane at $\tau=0.25\mathrm{ns}$, we compute the amount of zero pairs that cross the imaginary axis in the interval $[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\left(s,z\right)$ by replacing $z=e^{s\tau}$ in the denominator quasipolynomial. The solutions of $p\left(s,z\right)=0$ where $s$ lies on the imaginary axis and $z$ on the unit circle are called critical pairs, and can be found by solving a set of nonlinear polynomial equations ^{17}. In this case, $p\left(s,z\right)$ has one critical pair. Its first three corresponding critical delays are $\tau=0.077\mathrm{ps}$, $0.16\mathrm{ns}$ and $0.33\mathrm{ns}$. Two zero pairs have crossed the imaginary axis for our delay of $\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 $2\mathrm{GHz}$, was determined using Matlab’s fsolve
,
with its starting point at $s=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 SRP19.
Footnotes

A. Suarez, Analysis and Design of Autonomous Microwave Circuits. John Wiley and Sons, 2009. ↩ ↩^{2}

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

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

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

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

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. ↩

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

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

J. Jugo, J. Portilla, A. Anakabe, A. Suarez, and J. Collantes, “Closedloop stability analysis of microwave amplifiers,” Electronics Letters, vol. 37, no. 4, pp. 226228, Feb 2001. ↩

A circuit is considered realistic when its passives are lossy and when all active devices have a finite gainbandwidth. ↩

L. Baratchart, S. Chevillard, and F. Seyfert. (2014) On transfer functions realizable with active electronic components. hal01098616. Inria Sophia Antipolis. [Research Report] RR8659. ↩ ↩^{2}

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

K. Glover, “All optimal hankelnorm approximations of linear multivariable systems and their H∞ error bounds,” Internat. J. Control, vol. 39, no. 6, pp. 11151193, 1984. ↩

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. ↩

L. Shampine, “Vectorized adaptive quadrature in matlab,” Journal of Computational and Applied Mathematics, no. 211, pp. 131140, 2008. ↩

M. Biey, F. Bonani, M. Gilli, and I. Maio, “Qualitative analysis of the dynamics of the timedelayed chua’s circuit,” IEEE Tran. on Circuits and Systems, vol. 44, no. 6, pp. 486500, Jun 1997. ↩

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