The AAA algorithm provides a natural way to identify LTI (linear time-invariant) system parameters such as poles, zeros and DC gain from Bode plots. For example, consider the 4th order system $$ G(s) = 2\frac{(1+105s)(1+\frac{1.4}{0.05}s+\frac{1}{0.05^2}s^2)}{(1+100s)(1+\frac{1.4}{0.04}s+\frac{1}{0.04^2}s^2)(1+s)} $$
Nc = [8.4e4 6.68e3 2.66e2 2]; % numerator coefficients Dc = [6.25e4 6.6625e4 4.26e3 1.36e2 1]; % denominator coefficients N = @(s) Nc*s.^[3:-1:0]'; D = @(s) Dc*s.^[4:-1:0]'; G = @(s) N(s(:))./D(s(:)); % system transfer function
for which we have the following:
pol = roots(Dc), zer = roots(Nc), DCgain = abs(G(0))
pol =
-1.000000000000001 + 0.000000000000000i
-0.028000000000000 + 0.028565713714171i
-0.028000000000000 - 0.028565713714171i
-0.010000000000000 + 0.000000000000000i
zer =
-0.035000000000000 + 0.035707142142714i
-0.035000000000000 - 0.035707142142714i
-0.009523809523810 + 0.000000000000000i
DCgain =
2
Let's sample some values in the frequency range $10^{-4}\leq \omega \leq 10^2$, and draw the Bode plots of magnitude and phase:
w = logspace(-4,2,3000);
mag = abs(G(i*w)); ph = -angle(G(i*w));
LW = 'linewidth'; LO = 'location'; SW = 'southwest';
subplot(211), semilogx(w,20*log10(mag),'b-'), grid on
title('Magnitude (dB)')
subplot(212), semilogx(w,ph*180/pi,'b-'), grid on
title('Phase (degrees)')
Given magnitude and phase, an approximation $H(s)$ for $G(s)$ is readily obtained by AAA approximation of the complex signal. Samples are mirrored in order to enforce symmetry.
wA = [-fliplr(w) w]; magA = [fliplr(mag) mag]; phA = [-fliplr(ph) ph];
GA = magA.*exp(i*phA); % complex signal
[H,polA,resA,zerA] = aaa(GA,i*wA);
polA, zerA, DCgainA = abs(H(0)) % relevant parameters
subplot(211), hold on, semilogx(w,20*log10(abs(H(i*w))),'k--')
legend('G(s)','AAA',LO,SW)
subplot(212), hold on, semilogx(w,angle(H(i*w))*180/pi,'k--')
legend('G(s)','AAA',LO,SW)
polA = -0.999999999999993 + 0.000000000000045i -0.027999999999904 + 0.028565713714183i -0.009999999998222 - 0.000000000000566i -0.028000000001871 - 0.028565713713665i zerA = -0.035000000000000 + 0.035707142142714i -0.035000000000000 - 0.035707142142714i -0.009523809523810 - 0.000000000000000i DCgainA = 1.999999999999999
Also, $H(s)$ shows negligible errors in initial data:
err_mag = norm(mag-abs(H(i*w)),inf) err_ph = norm(ph-angle(H(i*w)),inf)
err_mag =
4.884981308350689e-15
err_ph =
2.569992829659640e-15
The following means of recomputing poles will play a key role in what follows. In the case when they aren't computed in complex conjugate pairs, which is a necessary condition for $H(s)$ to represent a physical model, we discard the imaginary parts of the corresponding polynomials and extract the new roots (poles). Here, poles and zeros do match:
[NcA,DcA] = residue(resA,polA,[]); polA = roots(real(DcA)), zerA = roots(real(NcA))
polA = -0.999999999999993 + 0.000000000000000i -0.028000000000887 + 0.028565713713924i -0.028000000000887 - 0.028565713713924i -0.009999999998222 + 0.000000000000000i zerA = -0.035000000000529 + 0.035707142142959i -0.035000000000529 - 0.035707142142959i -0.009523809523871 + 0.000000000000000i
Now let's complicate things a little bit. A reduced order approximation, useful in applications to simplify analysis and control design, is obtained by calling the AAA algorithm with a low degree. 20 Lawson iterations under the hood place our scarce resource (poles) at best, though not necessarily in complex conjugate pairs, hence we force a recomputation, and eventually solve a least-squares problem. This idea is a variant of the AAA-LS method introduced in [1].
Going back to our approximation problem, with a 2nd order reduction we expect two real distinct poles for the reduced $H_r(s)$. Note that in this case a straightforward zero-pole (over)simplification, leaving with $1/(1+s)$, wouldn't really be acceptable, and AAA-LS finds a decent compromise:
[~,polAr] = aaa(GA,i*wA,'degree',2);
polAr = roots(real(poly(polAr))); % pole recomputation
d = min(abs(i*wA(:)-polAr.'),[],1);
Q = d./(i*wA(:)-polAr.');
c = Q\GA.'; % solve LS problem, new residues c
Hr = @(s) [d./(s(:)-polAr.')]*c;
[NAr] = residue(c,polAr,[]);
zerAr = roots(real(NAr)), polAr, DCgainAr = abs(Hr(0))
subplot(211), semilogx(w,20*log10(abs(Hr(i*w))),'c-'), hold off
legend('G(s)','AAA','reduced order AAA',LO,SW)
subplot(212), semilogx(w,angle(Hr(i*w))*180/pi,'c-'), hold off
legend('G(s)','AAA','reduced order AAA',LO,SW)
zerAr = -0.294963170363907 polAr = -0.798508392918573 -0.035811189966051 DCgainAr = 2.036869044252427
To see how good AAA-LS actually is, consider the scalar example with noise found in [2], i.e. $$ f(s) = (s-1)/(s^2+s+2). $$ The function is sampled at 500 logarithmically spaced points in the interval [0.1,10], and then normally distributed noise with a standard deviation of $10^{-2}$ is added:
Nc = [1 -1]; Dc = [1 1 2];
N = @(s) Nc*s.^[1 0]';
D = @(s) Dc*s.^[2 1 0]';
f = @(s) N(s(:))./D(s(:));
w = logspace(-1,1,500); mag = abs(f(i*w)); ph = -angle(f(i*w));
mag = mag+0.01*randn(1,length(mag)); ph = ph+0.01*randn(1,length(ph));
subplot(211), semilogx(w,20*log10(mag),'r-'), grid on
title('Magnitude (dB)')
subplot(212), semilogx(w,ph*180/pi,'r-'), grid on
title('Phase (degrees)')
We compute a rational approximant of degree only 2 using the above method. The AAA-LS approximant shows no significant deviations from the measurements, at least in the eyeball norm. This time the number of Lawson iterations is increased, to enhance their filtering effect:
wn = [-fliplr(w) w]; magn = [fliplr(mag) mag]; phn = [-fliplr(ph) ph];
fn = magn.*exp(i*phn);
[~,poln] = aaa(fn,i*wn,'degree',2,'lawson',30);
poln(find(real(poln)>0)) ~= -1; % force system stability
poln = roots(real(poly(poln)));
dn = min(abs(i*wn(:)-poln.'),[],1);
Qn = dn./(i*wn(:)-poln.');
cn = Qn\fn.';
Hn = @(s) [dn./(s(:)-poln.')]*cn;
subplot(211), hold on, semilogx(w,20*log10(abs(Hn(i*w))),'b-'), hold off
legend('Noisy data','AAA approximant',LO,SW)
subplot(212), hold on, semilogx(w,angle(Hn(i*w))*180/pi,'b-'), hold off
legend('Noisy data','AAA approximant',LO,SW)
The poles are decently approximated, even in these perturbed conditions, as shown by the coefficients of the denominator:
Dcn = poly(poln)
Dcn = 1.000000000000000 1.001449519203945 1.998358534005181
The AAA-LS approximant effectively estimates the additive noise rather accurately, both in magnitude and in phase:
subplot(211), loglog(w,abs(mag(:)-abs(Hn(i*w))),'r-',LW,.5), grid on
title('Estimated noise in magnitude'), axis([min(w) max(w) 1e-5 1e-1]);
subplot(212), loglog(w,abs(ph(:)-angle(Hn(i*w))),'r-',LW,.5), grid on
title('Estimated noise in phase'), axis([min(w) max(w) 1e-5 1e-1]);
[1] S. Costa and L. N. Trefethen, AAA-least squares rational approximation and solution of Laplace problems, Proceedings of the 8ECM, 2021.
[2] I. V. Gosea and S. Güttel, Algorithms for the rational approximation of matrix-valued functions, arXiv:2003.06410v2, 2021.