# Why is my MATLAB's bode plot wildly off?

8 views (last 30 days)
Niket Parikh on 20 Jun 2020
Commented: Niket Parikh on 22 Jun 2020
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 Show 1 older comment
Niket Parikh on 20 Jun 2020
@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 on 20 Jun 2020
Share the ss equation in LaTeX form and the Mathematica input just to be sure you have the same equations used.
Niket Parikh on 20 Jun 2020
Alright.
The ss eqn.: The ss eqn. in Mathematica: Star Strider on 20 Jun 2020
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))
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.).
.

Star Strider on 21 Jun 2020
My pleasure!
Paul on 21 Jun 2020
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 on 22 Jun 2020
This was very helpful, Paul. I will keep these in mind. Thank you!

Paul on 20 Jun 2020
Edited: Paul on 20 Jun 2020
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))

#### 1 Comment

Niket Parikh on 20 Jun 2020
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?

R2018a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!