Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am using Mathematica to go through the examples and exercises on the book Modern Control Systems, 6th edition by Dorf. On page 605, there is a table (Table 8.5) with the Bode plot for several transfer functions. In what follows there is a piece of code that attempts to build the very same table.

Here is the code:

With[{ τ1 = 20, τ2 = 2, τ3 = 0.4, τ4 = 0.05, τa = 10, τb = 1, k = 10}, 
     Grid[
       Partition[
         Table[ BodePlot[ sys, PlotLabel->sys, GridLines -> Automatic], { sys, 
                { k/(s τ1 + 1), (k(s τa + 1))/(s(s τ1 + 1)(s τ2 + 1)), 
                  k/((s τ1 + 1)(s τ2 + 1)), k/s^2, k/((s τ1 + 1)(s τ2 + 1)(s τ3 + 1)),
                  k/(s^2 (s τ1 + 1)), k/s, (k(s τa + 1))/(s^2 (s τ1 + 1)), 
                  k/(s(s τ1 + 1)), k/s^3, k/(s(s τ1 + 1)(s τ2 + 1)), 
                 (k (s τa + 1))/s^3, (k (s τa + 1)(s τb + 1))/s^3, 
                 (k (s τa + 1))/(s^2 (s τ1 + 1)(s τ2 + 1)),
                 (k (s τa + 1)(s τb + 1))/(s(s τ1 + 1)(s τ2 + 1)(s τ3 + 1)(s τ4 + 1)) } 
                                                                         }
              ],  2], Frame->All, Spacings->6] ]

a bunch of Bode plots

All the transfer functions with 1/s^n ( n > 1 ) give the wrong result as far as the phase plot is concerned. Is there a simple way to fix this? Wolfram does not have a time line to go through the problem and sort it out.

share|improve this question
1  
Could you explain why "the result" is wrong ??? I edited your question, you can see how it should be done clicking edit under your question. –  Artes Jan 4 '13 at 13:13
    
Dear Artes - Many thanks. –  Ed Mendes Jan 4 '13 at 16:20
    
Dear Artes - Many thanks. Wrong means that it does not follow the standard convention when plotting a Bode Diagram. Check all the most adopted textbooks on Control Systems and see the Bode Diagram for the transfer function 1/s^2. I bet that in all of them the phase is -180 and not +180 as Mathematica BodePlot shows you. Why don't they follow the convention? 1 pole at s=0 gives -90, two poles at s=0 gives -180, 3 poles at s=0 gives -270 ... –  Ed Mendes Jan 4 '13 at 16:27
    
@EdMendes If you want to write a comment to a selected user write @name (e.g. to me @Artes). Have you read Details and Options in the documentation pages of BodePlot ? It says : option specifications include: opt->val use val for both the magnitude and the phase plot. If you find a correct solution to the problem you can answer your own question as well. This is a welcome practice. –  Artes Jan 5 '13 at 21:27
add comment

2 Answers

up vote 4 down vote accepted

To address the above comment by Ed

However if we change the example. Matlab - bode(tf(10*[10 1],[1 0 0 0])) - phase is negative (-270 to -180). Mathematica - See plots above - phase from +90 to +180. If instead of 10(10s+1)/sˆ3, one uses 10(10s+1)/(s+0.0001)ˆ3, the phase is negative

This below is a direct implementation of the phase plot part of Bode, using ArcTan. I used the following 2 transfer functions to compare with Mathematica BodePlot: 10(10s+1)/(s+0.0001)ˆ3 and 10(10s+1)/(s)ˆ3 and the result does show that there is a sudden phase change shift by 180 which does not show when using straight calculations using ArcTan to find the phase. Conclusion: There seems to be some convention used that causes this change as I would have expected it to match the ArcTan direct method.

Mathematica BodePlot phase diagram

Clear[s];
expr1 = (100 s + 10)/(s)^3;
expr2 = (100 s + 10)/(s + 0.0001)^3;

Grid[{{expr1, expr2},
  BodePlot[TransferFunctionModel[#, s], GridLines -> Automatic, 
      ImageSize -> 300, PlotLayout -> "List", 
      FrameLabel -> {{{"magnitude (db)", None}, {None,"Bode plot"}}, 
                    {{"phase(deg)",None}, {"Frequency (rad/sec)", None}}}, 
      ScalingFunctions -> {{"Log10", "dB"}, {"Log10", "Degree"}},
      PlotRange -> {{{0.001, 10}, All}, {{0.001, 10}, All}}][[2]] & /@ {expr1, expr2}}]

Mathematica graphics

Direct implementation of the phase plot using ArcTan

Clear[s];
ticks[min_, max_] := Table[{i, Superscript[10, i]}, {i, Ceiling[min], Floor[max], 1}];

makePhasePlot[expr_, s_Symbol] := Module[{ex, w, re, im, data},

   ex = expr /. s -> (w I);
   re = ComplexExpand[Re[ex]];
   im = ComplexExpand[Im[ex]];

   data =Table[{Log[10, w], 180/Pi ArcTan[re, im]}, {w, 0.001, 10, 0.001}];

   ListPlot[data,
    Joined -> True,
    PlotRange -> All,
    FrameTicks -> {{Automatic, Automatic}, {ticks, Automatic}},
    AxesOrigin -> {0, 0},
    Frame -> True,
    ImageSize -> 300,
    FrameLabel -> {{"angle(deg)", None}, {"Frequency (rad/sec)", None}},
    Axes -> False]
   ];

expr1 = (100 s + 10)/s^3;
expr2 = (100 s + 10)/(s + 0.0001)^3;
Grid[{{expr1, expr2},makePhasePlot[#, s] & /@ {expr1, expr2}}, Spacings -> {3, 0}]

Mathematica graphics

share|improve this answer
    
Thanks Nasser. That is exactly the point what I was trying to get across. There is a conflict within BodePlot that other control toolboxes do not seem to have (at least not in this case). Thanks to Mr. Wizard and Mr. James Cunnane for showing that I should not have used the word decent. –  Ed Mendes Jan 5 '13 at 12:30
add comment

No error seen in those phase plots.

For a transfer function of the form K/s^n, i.e. n poles at the origin, we expect constant phase of (-90 * n) degrees, plus or minus some integer multiple of 360 degrees - which is exactly what your Mathematica plots show.

share|improve this answer
    
I disagree - although -180 e +180 are the same thing - Mathematica BodePlot does not follow the usual standard convention when plotting a Bode Diagram - If you have a pole at s=0 the phase is -90 (270 sounds really weird). If you have two poles at s=0, one would expect -90 + -90 = -180 and not +180 (this means that the phase margin is 0 and not 360). Any decent control toolbox returns the angle in the standard convention. –  Ed Mendes Jan 4 '13 at 16:17
    
@EdMendes Can you provide a reference for your assertion? –  Mr.Wizard Jan 4 '13 at 16:37
    
@EdMendes Contrariwise - see MATLAB Bode plot documentation; the first and last Bode plots all show phases greater than 180 degrees. Unless your comment is intended to suggest that MATLAB is not included in "any decent control toolbox". –  James Cunnane Jan 4 '13 at 18:05
    
Many thanks for pointing out the Matlab documentation page. Please note that I am not talking about "phases greater than 180". Let us use the Matlab example - H = tf([1 0.1 7.5],[1 0.12 9 0 0]); bode(H) - In my mac Pro it shows negative phases (Actually the phase starts from -180 (due to the poles at s=0) goes near to -45 and comes back to -180. –  Ed Mendes Jan 4 '13 at 22:04
    
Let us do the same thing with Mathematica. Thanks to you I can get the same result if I use BodePlot[TransferFunctionModel[(s^2 + 0.1*s + 7.5)/(s^4 + 0.12*s^3 + 9*s^2), s]]. That is great! Many thanks. However if we change the example. Matlab - bode(tf(10*[10 1],[1 0 0 0])) - phase is negative (-270 to -180). Mathematica - See plots above - phase from +90 to +180. If instead of 10(10s+1)/sˆ3, one uses 10(10s+1)/(s+0.0001)ˆ3, the phase is negative. –  Ed Mendes Jan 4 '13 at 22:23
show 1 more comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.