PI Controller Design (Ziegler-Nichols) and Unit Step Response Graphing for 4th order system

I have a system that requires a PI controller designed using the Ziegler-Nichols Methods. I would also like to plot the unit step response in Matlab, but I am getting extremely weird results.

This is what the system looks like: $$G(s)=\frac{32}{(s+2)^4}$$

I plotted the curve in Matlab and printed the graph out to use the FOPDT approximant graphically. I know it is probably possible to do it using Matlab, but I come from a mathematics background and would rather use graphical methods.

Finding the quantities is trivial, and this is my result: $$F(s)=\frac{2}{2.519s+1}e^{-0.441s}$$ I put my graphical parameters and tuned it in Matlab. Could have done the approximation by inspection, this however should help for greater understanding. Anyway, these are the parameters: $$K_0=2$$ $$L=0.441$$ $$t=2.5188$$

Then using Ziegler-Nichol's table for tuning, I obtained the following PI controller: $$C(s)=5.1451\left(1+\frac{1}{1.4687s}\right)=\frac{7.557 s + 5.145}{1.469s}$$

I would now like to plot the closed loop step response, which should look like an unstable underdamped system that eventually stabilises. The system should oscillate drastically, such as:

This is what I get instead:

Which is totally wrong.

Here's the (relevant) matlab code:

%Unknown Transfer function

s=tf('s');

Gr=32/(s+2)^4;

% step(Gr)

% Approximation Taylor

K0=2; %Numerator gain factor

L=.4406; %Exponent power gain factor

t=2.5188; %Denominator gain factor

Go=(K0*exp(-L.*s)).*((t*s+1)^(-1));

%PI

Ti=L/.3;

Td=0;

Kp=.9*t/L; %gain

% Ki=Kp/Ti;

% Kd=Td*Kp;

Cpi=Kp*(1+1/(Ti*s));

Css=5.1451*(1+1/(1.4687*s));

%Closed Loop, Real

Gs=feedback(Gr*Css,1);

step(Gs)


(Feel free to edit this messy code list into something more tolerable)

I think my code is wrong somehow, I just don't understand how that result make sense.

The code is basically creating the 2 blocks, and then plotting for closed loop via the feedback command, unity feedback configuration.

Any tips and help is appreciated. Thanks.

Edit: Nichols-Ziegler Table and information

• Why have you gone from a 4th order TF to a Pade approximant time delay model?
– Chu
Sep 14, 2018 at 8:26
• Hi, I am using a first-order-plus-dead-time approximation, which to me just looks like the Taylor approximant, not Pade Approximant. I've tried Pade out and it worked nicely. Essentially I want to see what happens if I attempt to model it with a first order system and comment on it. Hope this clears up some things. Thanks! Sep 14, 2018 at 8:33
• Left out the code for plotting, I added it in. Sep 14, 2018 at 8:35
• Post those Ziegler-Nichols tables, I had never heard of Ziegler method dealing with dead time. Sep 14, 2018 at 9:39
• @MarkoBuršič Sure, I have updated the post with the table I am using. Sep 14, 2018 at 10:33

From what I can see the recipe table is not telling the whole truth about Kp. It is missing gain of the step response. This expression is valid as long as the step response is tending to the same value as the input step therefore their gain ratio is 1. But system G is responding to a unit step by rising to 2. So step/response ratio of the given system is 1/2.

Kp for PI controllers with Ziegler-Nichols First Method in fact is $$$$K_p=0.9⋅\frac{K_{step}}{K_0}⋅\frac{T}{L}$$$$

Because gain of a unit step is 1 this is $$$$K_p=0.9⋅\frac{1}{K_0}⋅\frac{T}{L}$$$$

Only with the step response tending to 1 so its overall gain K0 is equal to 1 as well this eventually becomes $$$$K_p=0.9⋅\frac{1}{1}⋅\frac{T}{L}=0.9\frac{T}{L}$$$$

as seen in the recipe table. But that is not the case with given sytem G. So gain ratio has to be considered as well as time ratio when obtaining Kp of G(t): $$$$K_p=0.9⋅\frac{1}{K_0}⋅\frac{T}{L}=0.9⋅\frac{1}{2}⋅\frac{2.5188}{0.441}\approx2.57$$$$

(Source: https://onlinelibrary.wiley.com/doi/pdf/10.1002/9780470029558.oth5 , page 292)

Also the penultimate line of the Matlab code

Gs=feedback(Gr*Css,1);


is looping 'Gr' which is the given system G(s) and 'Css' which is the FOPDT approximant for this very system G(s). Looping G(s) and the PI-controller 'Cpi' instead should to the trick.

Gs=feedback(Gr*Cpi,1);


Still that closed loop is unstable as it has two positive conjugated complex poles:

pole(Gs)
ans =
-4.0152 + 2.0976i
-4.0152 - 2.0976i
0.3426 + 2.0134i   % uh
0.3426 - 2.0134i   % oh
-0.6548 + 0.0000i


Using Photoshop for the graphical approximation (pixel perfect mesurement because why not) parameters of the FOPDT approximant are more like:

$$K_0=2$$ $$L=0.7$$ $$t=2.26$$

With these the step response is a stabilizing underdamped oscillation.

New Matlab code:

s=tf('s');%Time Transfer function

Gr=32/(s+2)^4;%system G(s)

% Parameters of FOPDT approximant
K0=2; %step response gain
L=.7; %step response delay time
t=2.26; %step response time constant estimate

% Ziegler-Nichols First Method
% PI controller
Ti=L/.3; %integral time
Kp=.9/K0*t/L; %proportional gain

Cpi=Kp*(1+1/(Ti*s));%transfer function of PI

%closed loop of G(s) and PI
Gs=feedback(Gr*Cpi,1);

step(Gs)