Why is my MATLAB's bode plot wildly off?

I will give out all the details in case it is relevant. I have a MIMO state space system. I find its bode plot using MATLAB and separately using Mathematica. The plot from MATLAB is wildly off compared to the correct bode plot. However, the Mathematica plot is quite close to the correct one. Here's the interesting thing. The transfer function whose bode plot is being taken in both the softwares is calculated to be the same. Have I made some mistake or is it numerical errors contributing to the problem? Here is relevant part of my code:

A=[0,0,0,0,1,0,0,0;
    0,0,0,0,0,1,0,0;
    0,0,0,0,0,0,1,0;
    0,0,0,0,0,0,0,1;
    -297.3,163.5,0,0,0,0,0,0;
    162.9,-267.2,104.2,0,0,0,0,0;
    0,57.8,-74.2,16.4,0,0,0,0;
    0,0,16.4,-16.4,0,0,0,0];
B=[0,0,0,0,0;
    0,0,0,0,0;
    0,0,0,0,0;
    0,0,0,0,0;
    131.4,0.046,0,0,0;
    0,0,0.045,0,0;
    0,0,0,0.025,0;
    0,0,0,0,0.25];
S=[1,0,0,0,0,0,0,0;
    0,0,0,1,0,0,0,0];
D=zeros(2,5);
sys=ss(A,B,S,D)
systf=tf(sys);
s1a1=systf(1,2);
bode(s1a1);

I have attached the plots.

Upper one is MATLAB and the other Mathematica

[madhan ravi]

Wildly off, what does mean? You told everything but forgot to share the SS equation xD.

[Niket Parikh]

@madhanravi The magnitude plot shows peak values at almost +50 dB while the actual plot never crosses zero dB. Also, the relative size of the resonant peaks is also not correct.

I have edited my question to include the ss equation.

[madhan ravi]

Share the ss equation in LaTeX form and the Mathematica input just to be sure you have the same equations used.

[Niket Parikh]

Alright.

The ss eqn.:

The ss eqn. in Mathematica:

Answer 1

[Star Strider]

This is a MIMO (5x2) system. Be careful that you plot the same input-outpair in MATLAB as Mathematica.

The Mathematica plot is Input 1 to Output 1 (or so it appears). So forget about the conversion to the transfer function and just do this:

[mag,phs,w] = bode(sys);
mag11 = squeeze(mag(1,1,:));                % Choose Input(1) To Output(1)

figure
semilogx(w, 20*log10(mag11))
xlabel('Frequency (rad/s)')
ylabel('Magnitude (dB)')
title('Input_1 To Output_1')

That plot looks to me to be the same as the Mathematica plot, other than the Mathematica plot being in terms of Hz (so Mathematica knows something MATLAB doesn’t.).

.

[Paul]

Can you post this plot? It doesn't look at like the Mathematica plot to me (even after noting the Mathematica plot is in Hz). The peaks in your plot are at +- 130 dB or so, The Mathematica doesn't get over -20 dB.

[Star Strider]

The plot obviously depends on choosing the same input-output pair, as I mentioned. That is not specifically stated with respect to the Mathematica plot and the MATLAB plot. I cannot determine from the description what are plotted in Mathematica.

The magnitudes of the pole peaks may depend on the frequency resolution. Higher frequency resolution will result in higher peak amplitudes, since the exact peaks are more likely to be plotted. I was considering the morphology of the plot, and the locations of the peaks.

I also mentioned that there is no need to convert to the transfer function first.

.

[Niket Parikh]

Thank you for the answer. So, higher frequency resolution should give better plots? Because in my case for very high resolution (i.e., no. of points) the plot deviates significantly from the correct one. Any idea why this happens?

[Star Strider]

My pleasure!

Higher frequency resolution will give more detailed plots. If you provide a ‘w’ (radian frequency) input vector, you can control where the ‘sys’ object is evaluated, and the frequency resolution. (The linspace or related logspace function is best for this.) Higher frequency resolution is more likely to result in plotting the outputs closer to the eigenfrequencies of the ‘A’ matrix, and so the amplitudes of the peaks will be higher.

I am not certain what you consider to be the ‘correct’ one. If you plot it as:

[mag,phs,w] = bode(sys);
mag11 = squeeze(mag(1,2,:));                % Choose Input(1) To Output(2)

figure
semilogx(w/(2*pi), 20*log10(mag11))
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
title('Input_1 To Output_2')
grid
xlim([0 70])

they look the same to me, except for the peak amplitudes, which are higher in the MATLAB plot for the reason I described. I have no idea how Mathematica chooses its frequencies, or if it uses a short-window moving-average filter before it executes the plot.

That code produces this plot:

The frequency axis is slightly truncated (to match the Mathematica plot), and the frequency axis converted from rad/s to Hz.

[Niket Parikh]

But that's the problem. I am reading a paper wherein the bode plot is given and I am assuming that to be correct and in that the magnitude always remains negative dB and doesn't cross 0 dB. However, as I increase my frequency resolution, my plot crosses 0 dB. Isn't something going wrong?

[Star Strider]

No, nothing is wrong. You are seeing the correct result.

The frequency resolution governs the output, so (unless you are also using a moving-average filter on the output before you plot it), the frequencies at which the system is evaluated will produce the result at those frequencies. You might want to specifically avoid the ‘A’ matrix eigenfrequencies to avoid the possibility of very high (perhaps infinite) results. Experiment with different ‘w’ vectors to see different results.

MATLAB optimises the frequency vector depending on the system. I consider this to be one of its strengths. You can override it with your own frequency vector.

Try this:

wv = logspace(-2, +2, 500)*2*pi;            % Radian Frequency Vector
bode(sys)
[mag,phs,w] = bode(sys, wv);
mag11 = squeeze(mag(1,2,:));                % Choose Input(1) To Output(2)

figure
semilogx(w/(2*pi), 20*log10(mag11))
xlabel('Frequency (Hz)')
ylabel('Magnitude (dB)')
title('Input_1 To Output_2')
grid
xlim([0 70])

Then experiment by making the frequency vector longer by significantly increasing the last argument to logspace. If you increase it past 5E+4, the pole peaks increase in amplitude to be significantly greater than 0 dB. This the same system, simply evaluated at different frequencies.

.

[Paul]

Star,

I misunderstood what your were concluding from your code snippet. You did show no need to convert to tf first; I was just showing that, as an alternative, one can subscript on the ss object directly if desired, instead of computing the full 2x5 frequency response and then pulling out the (1,1) element (or the (1,2) element) from that result. Apologies if you thought I was trampling on your response.

Niket,

In general higher frequency resolution will give you a better plot. You said that with more resolution you get deviation from the "correct one." Again, what are you using as the "correct one" for comparison? Try this:

[m,p,w]=bode(sys(1,2));m=squeeze(m);
ww=logspace(0,2,1000000);
[mm,pp]=bode(sys(1,2),ww);mm=squeeze(mm);
semilogx(w,db(m),ww,db(mm)),grid
set(gca,'xlim',[12 13])

Unless you show the code you're running and the results you're getting, it will be difficult to help much more.

Also, keep in mind that for your problem the the magnitude plot is a bit misleading. From the structure of your system, you know that the poles of the transfer function are the sqrts of the eigenvuales of the lower left block of the A matrix:

>> [cplxpair(eig(sys(1,2))) cplxpair([sqrt(eig(sys(1,2).a(5:8,1:4))); -sqrt(eig(sys(1,2).a(5:8,1:4)))])]

ans =

   7.9298e-17 - 2.1297e+01i   0.0000e+00 - 2.1297e+01i
   7.9298e-17 + 2.1297e+01i   0.0000e+00 + 2.1297e+01i
   5.8486e-16 - 1.2467e+01i   0.0000e+00 - 1.2467e+01i
   5.8486e-16 + 1.2467e+01i   0.0000e+00 + 1.2467e+01i
   1.1932e-15 - 6.2044e+00i   0.0000e+00 - 6.2044e+00i
   1.1932e-15 + 6.2044e+00i   0.0000e+00 + 6.2044e+00i
  -1.4658e-15 - 2.7612e+00i   0.0000e+00 - 2.7612e+00i
  -1.4658e-15 + 2.7612e+00i   0.0000e+00 + 2.7612e+00i

In reality, the poles are purely imaginary and the magnitude response would have a hole at those imaginary frequencies. But there is some round off error happening how CST computes the frequency response, so you get things like this:

  >> [m,p]=bode(sys(1,2),imag(ans(2,2)));squeeze(m)

ans =

   5.1822e+10
   

instead of Inf. So we know that Matlab won't give you the theoretically correct answer at theoretical singular frequency. You might ask how far away do you have to get from that frequency to get an numerically reliable result. I suspect it's not too far, but I don't know how much.

[Paul]

Star,

Isn't mag(1,2,:) the magnitude of the transfer function from input 2 to output 1?

Also, I think we were both typing at the same time, which is why our responses are so similar.

[Star Strider]

The output is the first dimension, the input is the second dimension, and the frequency is the third dimension. The important aspect here however is the morphology of the plot and what ‘correct’ implies in that context.

A magnitude of 5.1822E+10 is +214 dB, so considerably greater than 0 dB. That is my point. The frequencies at which the system is evaluated determine how the system outputs are plotted.

[Paul]

So mag(1,2,:) is the response from input 2 to output 1. But your code comment and the title of your plot say the opposite, which could cause confusion to someone not following closely.

My example with the 5.1822E+10 was directed at the OP's comment about resolution, not at yours (I'd link back to those comments if I knew how). I was basically saying that the OP's problem has some interesting aspects that need to be understood in order to decide what is "correct," and I think we are in agreement on that point, and pretty much everything else.

[Star Strider]

I am not concerned about the title of my plot. Changing my subscript references and my plot title do not change the point of my Answer.

My point exactly.

[Niket Parikh]

Star and Paul, I think I get my mistake now. Could you please tell me if this is correct? The transfer function will become unbounded at the frequencies corresponding to poles since this has poles on the imaginary axis. Hence in this case, increasing number of points gave me ever increasing peaks which is as expected. I was wrong in sticking to one of the plots as 'correct' because the correct one would have inf values at pole frequencies.

[Star Strider]

That is correct as I understand it.

Evaluating a system with poles on the imaginary axis will produce an infinite result at the eigenfrequencies of the ‘A’ matrix. (If they are to the left of the imaginary axis, they will be less than infinite, the amplitude depending on the magnitude of the negative real part.) Evaluating the system at an increasing number of points will increase the probability that they will be evaluated at or very close to the eigenfrequencies of the ‘A’ matrix.

With respect to it being ‘correct’, The system has Inf values at the pole frequencies. The Bode plot result depends on the frequencies where the system is evaluated.

[Niket Parikh]

@Star Things are clear now. Also, didn’t know that eigenfrequencies of matrix A are poles of TF so that was nice to know!

[Star Strider]

Great!

I am happy that you got all this sorted.

Have fun as you explore control systems with MATLAB. The Control Systems Toolbox (and related Toolboxes) make that much easier.

[Niket Parikh]

Sure! Thank you once again.

[Star Strider]

My pleasure!

[Paul]

Niket,

As you now know, the eigenvalues of the A-matrix are the poles of the TF, or the TF-matrix in the MIMO case such as yours. However, the poles and zeros that are used by the bode function to compute the response may not be exactly what they are supposed to be based on the physics of the problem that you're starting from. After all, the eigenvalues of your A-matrix should be exactly on the imaginary axis. But they are computed as having small real parts, which have an effect on the results. I'm pretty sure the same is true for the zeros of sys(1,2). Also, bode does some other scaling to try to improve the accuracy of your result. The bottom line is that the numerical result computed by bode is going to be different than what the physics says it should be, and it's up to you to understand your system well enough to know if such differences are worth worrying about. This is true regardless of how much frequency resolution you use to make your bode plot. For example, I don't think bode will ever return Inf magnitude for your system for any frequency input, even though we know that it should based on the structure of your A-matrix.

Another issue that you should be aware of is as follows. Suppose I plot the bode plot for some system hhh.

[m,p,w]=bode(hhh); semilogx(w,db(squeeze(m))); set(gca,'XLim',[0.8 1.2]), grid

The resulting plot is:

Something funny is happeing at w = 1 rad/sec. Let's look at the frequency vector that bode used around that critical frequency:

>> w(76:86)

ans =

     9.976358055610334e-01
     9.999999507491011e-01
     9.999999568092839e-01
     9.999999745683338e-01
     9.999999805285186e-01
     1.000000000000000e+00
     1.001642434954680e+00
     1.007627067240600e+00
     1.007627073347005e+00
     1.007627091241505e+00
     1.007627097247149e+00

So bode is trying to get very high resolution around w = 1 rad/sec. I could increase that resolution even further manually if I wanted, and the plot would likely change around w = 1. But would that be getting me closer to the "correct" answer? Let's look at the poles and zeros:

>> [z,p]=zpkdata(hhh);sort([z{:} p{:}])

ans =

     -4.890917532085481e-12 - 9.999999911204747e-01i     -7.088774012231625e-14 - 9.999999850602916e-01i
     -4.890917532085481e-12 + 9.999999911204747e-01i     -7.088774012231625e-14 + 9.999999850602916e-01i
      4.891111821114791e-12 - 1.000000008879525e+00i      7.098835408392290e-14 - 1.000000014839710e+00i
      4.891111821114791e-12 + 1.000000008879525e+00i      7.098835408392290e-14 + 1.000000014839710e+00i
      

We see that hhh has double poles very close +-1i and double zeros very close to +-1i. But they don't quite exactly cancel and so bode delivers the magnitude response as shown, with very rapid changes around w = 1. So the question is then, does hhh represent the actual system that I care about? Or should there actually have been a perfect pole/zero cancellation at +-1i, but the computations that led to hhh followed by the computations in bode are causing some peculiarities in the result. That's a judgement call that I'd have to make based on my fundamental understanding of the underlying system. If I thought that the poles and zeros should have cancelled, I can start to take appropriate action to improve things. See doc minreal for example.

In summary, listen to everying Star told you about using higher frequency resolution. In addition, you still need to use solid engineering judgement to decide if the result you're getting is what it should really be compared to the physical system you're modeling, regardless of the resolution you use.

In many cases, this type of issue won't really be important. But it can be an issue in fields where very lightly damped or undamped systems are common. Sometimes your system can even appear to be unstable when it really isn't.

[Niket Parikh]

This was very helpful, Paul. I will keep these in mind. Thank you!

Answer 2

[Paul]

How do you know what the correct Bode plot is that your taking as your reference for comparison?

Basically repeating from Star's comment for completeness: In Matlab, you are geting the Bode plot from the second input to the first output. Is that what you want? From the name s1a1 and from the title on the Mathematica plot, it looks like you may have wanted the Bode plot from the first input to first output. And be careful with the Hz vs. rad/s thing.

The bode function in the CST takes great care in selecting its frequency vector if not specified by the user (as in your case) in an attempt to make sure it captures important dyanmics. In your case, all of the poles and zeros of s1a1 are basically on the imaginary axis. If you don't catch exactly the right frequency to evaluate you'll miss a peak. Try this and see what you get:

bode(s1a1,logspace(-1,2,300))

What frequency vector did Mathematica use?

Finally, once you have the state space model, there's no reason convert to tf first. In fact, I'm quite certain that documentation for older versions of the CST specifically recommended NOT doing this. Don't know if that's the case now, but I'm pretty sure that if you start with a ss representation, you should just use it and not bother converting to the other forms. If all you wanted is the Bode plot from the second intput to the first output you get that directly:

bode(sys(1,2))

[Niket Parikh]

Thank you, this worked! Also, thanks for pointing out that I shouldn't convert to tf at all, didn't know that plots can be obtained directly from SS representation.

Also, I just observed that just like taking less points in logspace, taking too many also leads to inaccurate plots. Are you aware of why such a thing happens?