How to fix BodePlot that comes with Mathematica?

Question

I am using Mathematica to go through the examples and exercises on the book Modern Control Systems by Dorf, 6h edition. 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] ]

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.

Answer 1

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}}]

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}]

Answer 2

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.