PID Tuning via Ziegler Nichols Method

Hi all, I am trying to tune a PID using the Ziegler NIchols method, here is the transfer function :

I started following this video : https://www.youtube.com/watch?v=n829SwSUZ_c&t=1400s&ab_channel=ChristopherLum . I got to the part where he obtains the ultimate gain Ku from simulink. He set Kd=Ki=0 and only increased Kp until he got neutral oscilations. My problem is this : my value for "Ku" is 1.745 and my response looks like this :

And if I zoom in the waveform looks like this :

Which look nothing alike the waveforms I have seen online. Now if I increase my Kp above 1.8 the response of the system looks like this :

Now my question is: what am I doing wrong ? If I use that value for Ku=1.745 and continue with the calculations the response looks like this :

Which looks pretty wrong. So what I am doing wrong ? also if i try in matlab to plot the system response i get this plot :

Below is my block scheme and step input settings :

Answer 1

[Sam Chak]

Hi @Tom

If my recollection serves me right, the Ziegler–Nichols method requires the process plant to be open-loop stable. However, in the case of your plant, , it is open-loop unstable due to a pole at the origin. To enable the application of the Ziegler–Nichols method, we aim to create a form of 'open-loop stability' from a tuning standpoint by establishing a feedback loop utilizing the original plant's output. I've denoted this modified plant as , and now the Ziegler–Nichols method can be implemented effectively.

%% Plant (open-loop unstable)
Gp  = tf(523500, [1 87.35 10470 0])

Gp =
 
           523500
  -------------------------
  s^3 + 87.35 s^2 + 10470 s
 
Continuous-time transfer function.

%% Plant (make it open-loop stable)
Gpp = feedback(Gp, 1)

Gpp =
 
                523500
  ----------------------------------
  s^3 + 87.35 s^2 + 10470 s + 523500
 
Continuous-time transfer function.

step(Gpp, 1), grid on

%% Ziegler–Nichols Gains
Ku         = 0.747;                 % Find this Ultimate Gain
G          = minreal(feedback(Ku*Gpp, 1));
step(G, 1), grid on

[y, t]     = step(G, 1);
h          = detrend(y);
[~, peaks] = findpeaks(h, 'MinPeakProminence', (max(h)-min(h))/4);
Tu         = mean(diff(t(peaks)))   % Oscillation period

Tu = 0.0614

%% Classic PID Controller synthesis from Z–N Table
Kp  = 0.6*Ku;
Ti  = Tu/2;
Td  = Tu/8;
Ki  = Kp/Ti;
Kd  = Kp*Td;
Gc  = pid(Kp, Ki, Kd)

Gc =
 
             1          
  Kp + Ki * --- + Kd * s
             s          

  with Kp = 0.448, Ki = 14.6, Kd = 0.00344
 
Continuous-time PID controller in parallel form.

%% Closed-loop system
Gcl = minreal(feedback(Gc*Gpp, 1))

Gcl =
 
             1801 s^2 + 2.346e05 s + 7.644e06
  ------------------------------------------------------
  s^4 + 87.35 s^3 + 1.227e04 s^2 + 7.581e05 s + 7.644e06
 
Continuous-time transfer function.

stepinfo(Gcl)

ans = struct with fields:
         RiseTime: 0.1395
    TransientTime: 0.3266
     SettlingTime: 0.3266
      SettlingMin: 0.8742
      SettlingMax: 1.0030
        Overshoot: 0.2976
       Undershoot: 0
             Peak: 1.0030
         PeakTime: 0.4064

step(Gcl, 1), grid on

[William Rose]

This is a very impressive answer.

[Tom]

Thanks for the repply @Sam Chak. It really helped and as @William Rose mentioned it is a very impressive answer. Thanks again, and after I changed the parameter maximum step size in Configuration Parameters the simulink model looks good again.

[Sam Chak]

You are welcome, @Tom. Also thank @William Rose.

[Siddharth Jawahar]

Hello! You could use the pidtune command to automatically tune the pid gains like so:

% Define the transfer function
s = tf('s'); % Define Laplace variable
G = 523500 / (s^3 + 87.35*s^2 + 10470*s); % Plant transfer function

% Set PID tuning options for balanced response
opts = pidtuneOptions('DesignFocus', 'balanced');

% Use pidtune to tune the PID controller with the specified options
[C, info] = pidtune(G, 'PID', opts);

% Display the tuned PID controller parameters
disp('Tuned PID Controller with Balanced Design Focus:');

Tuned PID Controller with Balanced Design Focus:

disp(C);

  pid with properties:

              Kp: 0.9552
              Ki: 7.1210
              Kd: 0.0109
              Tf: 0
        IFormula: ''
        DFormula: ''
      InputDelay: 0
     OutputDelay: 0
       InputName: {''}
       InputUnit: {''}
      InputGroup: [1x1 struct]
      OutputName: {''}
      OutputUnit: {''}
     OutputGroup: [1x1 struct]
           Notes: [0x1 string]
        UserData: []
            Name: ''
              Ts: 0
        TimeUnit: 'seconds'
    SamplingGrid: [1x1 struct]

% Display closed-loop performance metrics
disp('Tuning Info:');

Tuning Info:

disp(info);

                Stable: 1
    CrossoverFrequency: 79.5646
           PhaseMargin: 69.7955

% Analyze closed-loop response
T = feedback(C * G, 1); % Closed-loop transfer function
figure;
step(T); % Step response of the closed-loop system
grid on;
title('Closed-Loop Step Response');

Hope this helps,

Sid

Answer 2

[Sam Chak]

Hi @Tom

Regarding your original issue, you can improve the display of the waveform from the Simulink Scope by reducing the maximum step size in the Configuration Parameters dialog box. In the Simulink Editor, on the Modeling tab, click Model Settings.

By the way, I attempted to implement the Ziegler–Nichols tuning method on the initial unstable plant by determining the ultimate gain () until the output displays sustained oscillations. Surprisingly, the Ziegler–Nichols method, employing the 'some overshoot' criterion, successfully stabilizes the plant. However, it is noteworthy that the overshoot in this case amounts to 30%.

%% Plant (open-loop unstable)
Gp  = tf(523500, [1 87.35 10470 0])

Gp =
 
           523500
  -------------------------
  s^3 + 87.35 s^2 + 10470 s
 
Continuous-time transfer function.

ToF = isstable(Gp)  % 1 means stable, 0 means unstable

ToF = logical
   0

%% Ziegler–Nichols Gains
Ku         = 1.747;                 % Find this Ultimate Gain
G          = minreal(feedback(Ku*Gp, 1));
step(G, 1), grid on

[y, t]     = step(G, 1);
h          = detrend(y);
[~, peaks] = findpeaks(h, 'MinPeakProminence', (max(h)-min(h))/4);
Tu         = mean(diff(t(peaks)))   % Oscillation period

Tu = 0.0614

%% PID Controller synthesis with 'Some Overshoot' from Z–N Table
Kp  = Ku/3;
Ti  = Tu/2;
Td  = Tu/3;

Ki  = Kp/Ti;
Kd  = Kp*Td;
Gc  = pid(Kp, Ki, Kd)

Gc =
 
             1          
  Kp + Ki * --- + Kd * s
             s          

  with Kp = 0.582, Ki = 19, Kd = 0.0119
 
Continuous-time PID controller in parallel form.

%% Closed-loop system
Gcl = minreal(feedback(Gc*Gp, 1))

Gcl =
 
             6239 s^2 + 3.049e05 s + 9.931e06
  ------------------------------------------------------
  s^4 + 87.35 s^3 + 1.671e04 s^2 + 3.049e05 s + 9.931e06
 
Continuous-time transfer function.

stepinfo(Gcl)

ans = struct with fields:
         RiseTime: 0.0227
    TransientTime: 0.3922
     SettlingTime: 0.3922
      SettlingMin: 0.9047
      SettlingMax: 1.3023
        Overshoot: 30.2306
       Undershoot: 0
             Peak: 1.3023
         PeakTime: 0.0895

step(Gcl, 1), grid on

Answer 3

[종혁]

뽕빵삥ㅃ

ㅜㅇ뿡