8 views (last 30 days)

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

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))

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

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))

Opportunities for recent engineering grads.

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

Start Hunting!
## 4 Comments

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/551758-why-is-my-matlab-s-bode-plot-wildly-off#comment_906556

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/551758-why-is-my-matlab-s-bode-plot-wildly-off#comment_906556

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/551758-why-is-my-matlab-s-bode-plot-wildly-off#comment_906562

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/551758-why-is-my-matlab-s-bode-plot-wildly-off#comment_906562

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/551758-why-is-my-matlab-s-bode-plot-wildly-off#comment_906583

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/551758-why-is-my-matlab-s-bode-plot-wildly-off#comment_906583

## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/551758-why-is-my-matlab-s-bode-plot-wildly-off#comment_906589

⋮## Direct link to this comment

https://in.mathworks.com/matlabcentral/answers/551758-why-is-my-matlab-s-bode-plot-wildly-off#comment_906589

Sign in to comment.