How to find the value of PID using ziggler 2 method

Hi, i have a tf below

G1 =

350 s + 20000

------------------------------------------------------

70000 s^4 + 192500 s^3 + 3.7e07 s^2 + 5.25e07 s + 3e09

how do i find the Kp, Ki, Kd using the ziggler 2 method, i've tried by PID but i couldn't find anything.

[Raka Adhwa]

and also the graphic

[Mathieu NOE]

see

A MATLAB Approach to study different types of Ziegler-Nichols P-I-D Controller Tuning Algorithm - File Exchange - MATLAB Central (mathworks.com)

Answer 1

[Sam Chak]

Hi @Raka Adhwa

Unsure if you are still interested in the problem. But if the model is available, then the analysis can be performed to learn and design the PID compensator in a more meaningful way.

s  = tf('s');
Gp = (350*s + 20000)/(70000*s^4 + 192500*s^3 + 3.7e07*s^2 + 5.25e07*s + 3e09);
Gp = minreal(Gp)

Gp =
 
                0.005 s + 0.2857
  ---------------------------------------------
  s^4 + 2.75 s^3 + 528.6 s^2 + 750 s + 4.286e04
 
Continuous-time transfer function.

steady_state_Gp = 20000/(3e09)

steady_state_Gp = 6.6667e-06

step(Gp)

From the Denominator of the Plant

it is clear that is a Type-0 system. Thus, introducing an Integral action

will make

becoming a Type-1 system and from the closed-loop system transfer function

the steady-state error will be 0 because .

The closed-loop system characteristic polynomial is

.

We want to select an appropriate value for the Integral gain so that the closed-loop system behaves like the desired 5th-order system with the following Characteristic Polynomial:

% Desired 5th-order Characteristic Polynomial
cp5 = conv([1 2 1], [1 3 3 1])

cp5 = 1×6

1 5 10 10 5 1

Comparing with 4th term from the desired 5th-order Characteristic Polynomial and solve for :

fun = @Kigain1;
x0  = [5000, 5];
x   = fsolve(fun, x0)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

x = 1×2

1.0e+03 * 4.6687 0.0042

Gc  = x(1)/s;
Gcl = minreal(feedback(Gc*Gp, 1))

Gcl =
 
                       23.34 s + 1334
  --------------------------------------------------------
  s^5 + 2.75 s^4 + 528.6 s^3 + 750 s^2 + 4.288e04 s + 1334
 
Continuous-time transfer function.

step(Gcl)

Zero steady-state error for a unit step input is expected, but the settling time at 300 s is probably too long. So, we use the 3rd term from the desired 5th-order Characteristic Polynomial and solve for :

fun = @Kigain2;
x0  = [70000, 7];
x   = fsolve(fun, x0)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

x = 1×2

1.0e+04 * 7.1106 0.0007

Gc  = x(1)/s;
Gcl = minreal(feedback(Gc*Gp, 1))

Gcl =
 
                       355.5 s + 2.032e04
  ------------------------------------------------------------
  s^5 + 2.75 s^4 + 528.6 s^3 + 750 s^2 + 4.321e04 s + 2.032e04
 
Continuous-time transfer function.

step(Gcl)

The settling time is significantly shorter, and we know that high value of in the range of should be selected. We can test multiple values for and use the for loops to run repetitive computations:

Ki  = 6e4:1.5e4:12e4                    % generate evenly-spaced multiple values for Ki

Ki = 1×5

60000 75000 90000 105000 120000

t   = 0:0.01:20;

for j   = 1:5
    Gc  = Ki(j)/s;                      % integral compensator
    Gcl = minreal(feedback(Gc*Gp, 1));  % closed-loop system, Gcl
    [y(:, j), t] = step(Gcl, t);        % generate data for step response
    info(j)      = stepinfo(Gcl);       % performance of the step response
end

plot(t, y), grid on, xlabel('t, [sec]'), ylabel('Response')
title('Plot of Step Response Curves')
legend(strcat('K_i = ', num2str(Ki')), 'location', 'Best')

One can evaluate the performance in terms of Rise time, Settling time, and Percentage overshoot to decide on the range and selection. For example, the performance of the closed-loop system #5 (green curve) can be obtained:

info(5)

ans = struct with fields:
         RiseTime: 2.0618
    TransientTime: 8.2600
     SettlingTime: 8.2600
      SettlingMin: 0.8150
      SettlingMax: 1.0217
        Overshoot: 2.1661
       Undershoot: 0
             Peak: 1.0217
         PeakTime: 6.0432

--------------------------------------------

Functions that are used to initially determine the value for :

function F = Kigain1(x)

F(1) = 750 - 10*x(2)^3;
F(2) = 0.2857*x(1) - x(2)^5;

end

function F = Kigain2(x)

F(1) = 528.6 - 10*x(2)^2;
F(2) = 0.2857*x(1) - x(2)^5;

end