# How to fix BodePlot that comes with Mathematica?

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


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.

-
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

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


-
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