PID Tuning using Ziegler Nicholas - MATLAB Approach

This is an unreleased lab for undergraduate Mechatronics students to know how to practice Ziegler Nicholas method to find the PID factors using MATLAB.

1. v1.0
Misr University for Science and Technology
College of Engineering
Mechatronics Lab








PROCESS CONTROL MODULE
PID TUNING AND STABILITY
(MATLAB Simulation)




Prof. Farid A. Tolbah
Eng. Waleed A. El-Badry

2. v1.0
1. Objective
The experiment is aimed to make student acquainted with the preliminary steps to manually
tune a PID controller for a plant model by means of MATLAB.
2. Outcome
 Writing mathematical models of plants under investigation in Laplace form using MATLAB.
 Developing the mathematical Laplace representation of Ziegler-Nicholas PID controller in
MATLAB.
 Finding the critical gain (Kc) and the ultimate period (Pu) to calculate PID gains.
3. Prerequisite
Student should be familiar with the following terms:
 Closed loop system.
 System response.
 PID controller
 Ziegler-Nicholas tuning method.
Also basic understanding of MATLAB syntax is preferred.
4. The Closed loop system
The below figure represents the generic closed loop system.






Figure 1 Closed loop system
For implementation in this experiment, we are given the following plant model Wp(s):
??????
??????(??????)=
1
(10??????+1)
3
=
1
1000�
3
+300??????
2
+30??????+1

And the Ziegler-Nicholas PID is formulated as:
??????
�(??????)=??????
�(1+
1
�
????????????
+�
�??????)

Controller Plant
Feedback
FCE

Set
Point
e(t) y(t) wc
Current

Level

3. v1.0
Assuming unity feedback, redrawing the block diagram:
Figure 2 Problem block diagram
5. Developing MATLAB functions (Plant, Controller and Closed Loop)
a. Launch MATLAB software.
b. From the Home tab, select New -> Function.
c. Write down the generic plant function as shown in the following snippet:
function [ Wp ] = CreatePlant( num,den )
%CreatePlant Creates plant transfer function.
% The returned value is the system in numerator/denomerator format
%% Parameters
% num : Numerator vector (starting from highest order of
coefficients)
% den : Denomerator vector (starting from highest order of
coefficients)
% plant : Plant transfer func tion
%% EXAMPLE
% num=[1];
% den=[1 0 1];
% sys=CreatePlant(num,den)

%% Result is
% 1
% sys= ---------------
% S^2+1
%% Function implementation
syms s;
Wp=tf(num,den);
end
Snippet 1 CreatePlant function

d. Save the file.
e. Close the function file.

??????
�(1+
1
�
????????????
+�
�??????)
1
(10??????+1)
3

Set
Point
E(S) X(S)
Current

Level


Y(S)

4. v1.0

f. Repeat steps b-e for creating the following snippet for Ziegler-Nicholas generic function:
function Wc = ZieglerNicholasPID( Kc,Ti,Td )
% ZieglerNicholasPID function to generate the PID controller
transfer
%% Parameters
% Kc : Critical gain
% Ti : Reset time (minutes)
% Td : Derivative time (minutes)

%% Function implementation
s=tf('s');
Wc=Kc*(1+(1/(Ti*s))+Td*s);
end
Snippet 2 Ziger-Nicholas PID implementation

g. The final function bonds the two functions (plant and controller) to build the closed loop
system:
function sys = CLS( Wp,Wc )
%CLS Closed loop system function

%% Parameters
% Wp : Plant transfer function
% Wc : Controller transfer function
% sys : Closed Loop transfer function with assuming unity
feedback.

%% Function implementation
CLS=feedback(series(Wp,Wc),1);
end
Snippet 3 Closed loop system bonding

5. v1.0
6. Open loop system response




Figure 3 Open loop system
To plot the open loop response, perform the following steps:
a. From MATLAB command window, we will call the function CreatePlant to create the
transfer function mentioned in shown:

sys=CreatePlant(1,[1000 300 30 1]);
step(sys)

b. From the figure opened, right click on it and select characteristics -> Settling Time, Rise
Time and Steady State. Fill in the table:

















Figure 4 Characteristics of Open loop system
\
Table 1 Characteristics of open loop system
Rise Time (sec) 42.2
Settling Time (sec) 75.2
Steady State (sec) 120

1
(10??????+1)
3

Set
Point
Y(s)

6. v1.0
7. Finding the critical gain (Kc) via Nyquist plot
a. To plot the Nyquist of frequency response of the plant, write down the following code:
Wp=CreatePlant(1,[1000 300 30 1]);
nyquist(Wp);

b. Right click on the plot and select characteristics -> Minimum Stability
Margins as shown in figure



















Figure 5 Nyquist plot (Open loop)
c. Write down the gain margin Gm (in dB) and convert it to magnitude. Write down the
margin frequency Wc .










d. Calculate ??????
�=??????
??????= 8.0011 , and ??????
??????=
2??????
????????????
=
2??????
0.173
=36.32 ?????????????????? and consequently
�
??????=

7. v1.0

e. Check that Kc is the critical gain by writing down the following MATLAB code:

t=0:0.01:200;
Wp=CreatePlant(1,[1000 300 30 1]);
%Setting Kc=8, Ki=~0 and Kd=0
Wc=ZieglerNicholasPID(8,100000,0);
sys=CLS(Wp,Wc);
%plotting step response from t0=0 to tf=200 sec
step(sys,t)
Snippet 4 Plotting the system response at critical gain















Figure 6 Step response at Kc=8
8. Calculating P, PI and PID control gains
After obtaining the critical gain from the previous step, we are able to calculate the P,I and D
parameters and perform comparison of each controller type. According to Ziegler Nicholas
table:
Table 2 Ziegler Nicholas Tuning Chart
Controller Type Kp Ti (sec) Td (sec)
P 0.5*Kc = 0.5*8=4 100000 0
PI 0.45*Kc = 0.45*8=3.6 0.83*Pu=0.83*36.32=30.1 0
PID 0.59* Kc = 0.59*8=4.7 0.5*Pu=0.5*36.32=18.2 0.12*Pu =0.12*36.32=4.4

8. v1.0
Plot the step response of each controller over the plant by writing the following code:
Wp=CreatePlant(1,[1000 300 30 1]);
Wcp=ZieglerNicholasPID(4,100000,0);
Wcpi=ZieglerNicholasPID(3.6,30.1,0);
Wcpid=ZieglerNicholasPID(4.7,18.2,4.4);
t=0:0.01:500;
sys=CLS(Wp,Wcp);
step(sys,t)
hold on
sys=CLS(Wp,Wcpi);
step(sys,t)
sys=CLS(Wp,Wcpid);
step(sys,t)
legend('P','PI','PID')

Figure 7 step response of P, PI and PID controller