PI controller does not have the effect I would expect from the calculations

I am designing a PI controller for the following plant

K = 11041666;
B = 100;
M = 1019;
num = [K];
den = [M B K];
sys = tf(num,den);

I need to design for 10% overshoot and a settling time of 5 seconds

The closed loop transfer function is as follows:

syms k m b s kp ki
G = k/(m*s^2+b*s+k);
Gc = (kp*s+ki)/s;
T = (kp*s*k+ki*k)/(m*s^3+b*s^2+k*s*(1+kp)+ki*k)

These are the calculations for Kp and Ki

overshoot = 10/100;
Zeta = (-log(overshoot))/sqrt(pi^2+(log(overshoot)^2));
Wn = 5/(4*Zeta);

Ki = (Wn^3*M)/K;
Kp = ((2.15*Wn^2*M)/K)-1;

However with these values this is the response I get

The block diagram

The controller internals

The step input

Yellow if the target and blue is the system response

I am making a mistake somewhere but I have been unable to find it. Any assistance would be greatly appreciated

[Aquatris]

Any reason why you are not using a PID? Your system is lightly damped. I doubt you can achieve both stability and performance with a PI controller. You would need a damping term in your controller.

K = 11041666;
B = 100;
M = 1019;
num = [K];
den = [M B K];
sys = tf(num,den);
damp(sys)

                                                                       
         Pole              Damping       Frequency      Time Constant  
                                       (rad/seconds)      (seconds)    
                                                                       
 -4.91e-02 + 1.04e+02i     4.71e-04       1.04e+02         2.04e+01    
 -4.91e-02 - 1.04e+02i     4.71e-04       1.04e+02         2.04e+01    

The formulas for the overshoot etc are valid if the system is stable already. They do not tell if a system is stable or not. Also I am not sure about the formulas you use to get Kp and Ki from Wn,M and K values.

[Wynand]

I was recommened to use the PI controller but have no real reason to. I will redo the calculations for a PID controller

Answer 1

[Sam Chak]

Hi @Wynand

Please see my explanations below.

%% Plant
K           = 11041666;
B           = 100;
M           = 1019;
num         = [K];
den         = [M B K];
Plant       = tf(num, den)

Plant =
 
           1.104e07
  ---------------------------
  1019 s^2 + 100 s + 1.104e07
 
Continuous-time transfer function.

Step Response of the Uncontrolled Plant:

The response is highly oscillatory, and the settling time takes roughly 100 seconds (which is way too long). No real elevator behaves in this manner. This suggests that the modeling is incorrect, an issue that should be addressed by the person who posted the unverified solution in the original thread or by posing a new question to modeling experts.

figure(1)
step(Plant), grid on, grid minor

Controller Design:

The damping ratio has been calculated correctly; however, the natural frequency of the desired closed-loop system has not been. It seems that you may have attempted to use an approximation formula but forgot to verify its correct application. The precise value (not an approximation) of ω is 1.18516590360739 rad/s. Additionally, a 1st-order PI controller is insufficient to provide adequate compensation for a generic 2nd-order plant. Therefore, a 2nd-order PID controller should be employed.

%% Kontroler (Polish)
overshoot   = 10/100;
Zeta        = (-log(overshoot))/sqrt(pi^2+(log(overshoot)^2))

Zeta = 0.5912

Wn          = 3.9/(5*Zeta)

Wn = 1.3195

Wn          = 5/(4*Zeta)    % super precise value is 1.18516590360739 rad/s

Wn = 2.1145

Ki          = (Wn^3*M)/K

Ki = 8.7250e-04

Kp          = ((2.15*Wn^2*M)/K)-1

Kp = -0.9991

Kontroler   = pid(Kp, Ki)

Kontroler =
 
             1 
  Kp + Ki * ---
             s 

  with Kp = -0.999, Ki = 0.000873
 
Continuous-time PI controller in parallel form.

Analysis of the Closed-Loop System:

After designing the controller, the closed-loop poles should be checked to ensure that all real parts are negative. Clearly, at least one pole has a positive real part, as indicated by the response, which explodes, as shown in the figure below.

%% Closed-loop system
ClosedLoop  = minreal(feedback(Kontroler*Plant, 1))

ClosedLoop =
 
          -1.083e04 s + 9.454
  -----------------------------------
  s^3 + 0.09814 s^2 + 9.613 s + 9.454
 
Continuous-time transfer function.

pole(ClosedLoop)    % check for positive real part

ans = 
   0.4074 + 3.1923i
   0.4074 - 3.1923i
  -0.9129 + 0.0000i

figure(2)
step(ClosedLoop), grid on, grid minor

Suggested Control Solution:

If you wish to continue pursuing a control solution for this questionable elevator model, consider using the following PID control configuration. No overshoot is always preferable for the passengers' experience.

%% PID controller
Ts       = 5;                               % Desired Settling Time
[Gc, Gh] = chakpid(Plant, Ts)               % Gc in forward path, Gh in feedback path

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

  with Kp = 0.000113, Ki = 4.42e-05, Kd = 7.22e-05
 
Continuous-time PID controller in parallel form.

Gh =
 
  2.875 s^2 - 1.385e04 s + 0.6122
  -------------------------------
      s^2 + 1.565 s + 0.6122
 
Continuous-time transfer function.

%% Closed-loop System
Gcl      = minreal(feedback(Gc*Plant, Gh))

Gcl =
 
    0.7824 s^4 + 2.449 s^3 + 2.874 s^2 + 1.499 s + 0.2932
  ---------------------------------------------------------
  s^5 + 3.912 s^4 + 6.122 s^3 + 4.79 s^2 + 1.874 s + 0.2932
 
Continuous-time transfer function.

S        = stepinfo(Gcl)                    % Performances

S = struct with fields:
         RiseTime: 2.8083
    TransientTime: 5.0000
     SettlingTime: 5.0000
      SettlingMin: 0.9005
      SettlingMax: 0.9986
        Overshoot: 0
       Undershoot: 0
             Peak: 0.9986
         PeakTime: 8.3900

figure(3)
step(Gcl, round(3*Ts, 0)), grid on, grid minor

%% Double compensation control scheme
function [C, H] = chakpid(P, Ts)
    %  Gp is a 2nd-order Plant
    %  Ts is the desired settling time

    [numP, denP] = tfdata(minreal(P), 'v');
    a1  = denP(2);
    a2  = denP(3);
    b   = numP(3);
    Tc  = -log(0.02)/Ts;                        % Desired time constant

    %% The formulas
    k1  = (3*(b^2)*((Tc/b)^2) - a2)/b;      
    k2  = (b^2)*((Tc/b)^3);
    k3  = (3*b*(Tc/b) - a1)/b;       
    k4  = 2*b*((Tc/b)^2);                       % P-gain of PID controller
    k5  = k2;                                   % I-gain of PID controller
    k6  = Tc/b;                                 % D-gain of PID controller

    C   = pid(k4, k5, k6);                      % Classical PID controller
    H   = minreal(tf([k3 k1 k2], [k6 k4 k5]));  % Compensator at the Sensor path
end

[Wynand]

Thank you for your time.

I have followed another approach to the design but the respone time is way too fast.

From first principles the values for the parameters were found to be:

I determined the closed loop transfer function to be:

syms Kp Ki Kd s m k b Wn
T = (Kp +Ki/s+s*Kd)/(m*s^2+b*s+k+Kp+Ki/s+s*Kd)

Then using the ideal 3rd order transfer function equation I calculated the PID parameters:

ITF = s^3+1.75*Wn*s^2+2.15*Wn^2*s+Wn^3;

K = 12421874;
M = 1019;
B = 100;
Wn = sqrt(K/M)

Wn = 110.4095

Kd = ((1.75)*(Wn*M)-B)

Kd = 1.9679e+05

Ki = (M*(Wn^3))

Ki = 1.3715e+09

Kp = ((2.15*(Wn^2)*M)-K)

Kp = 1.4285e+07

Implementing that into my controller results in the following graph

By manually tuning a gain (0.01) placed before the controller, I can get to the results I want, but I prefer to not have to tune it by hand at all. I also see small oscillations that are still present in the position response that I would like to remove completely. Any tips?

[Sam Chak]

Hi @Wynand

Your design methodology is somewhat correct from a stabilization perspective, as you have only placed the closed-loop poles according to the Hurwitz polynomial of your ideal transfer function. Therefore, stabilization is guaranteed. However, the transient response is not guaranteed due to the non-ideal nature of the numerator in the actual closed-loop transfer function.

If you compare the characteristic polynomial of the actual closed-loop transfer function with that of the ideal transfer function, it appears that your computed PID gain has not successfully placed the closed-loop poles at the desired locations, possibly due to incorrect formulas.

Nonetheless, I cannot replicate your result with a gain factor of 0.01. A gain factor of appears to provide a satisfactory response (under 5 seconds), albeit rather chattering.

It is important to understand that Control Theory indicates that a high value of ω leads to a very short settling time, and vice versa, provided that the system is stable.

K   = 12421874;
M   = 1019;
B   = 100;
num = [K];
den = [M B K];
Gp  = tf(num, den);
Wn  = sqrt(K/M)     % No control theory books will tell you to use this formula

Wn = 110.4095

%% Ideal Transfer Function
s   = tf('s');
ITF = Wn^3/(s^3 + 1.75*Wn*s^2 + 2.15*Wn^2*s + Wn^3)

ITF =
 
                 1.346e06
  ---------------------------------------
  s^3 + 193.2 s^2 + 2.621e04 s + 1.346e06
 
Continuous-time transfer function.

pole(ITF)   % locations of desired closed-loop poles

ans = 
   1.0e+02 *

  -0.5752 + 1.1793i
  -0.5752 - 1.1793i
  -0.7818 + 0.0000i

figure(1)
step(ITF), grid on, grid minor

Kd  = ((1.75)*(Wn*M)-B);
Ki  = (M*(Wn^3));
Kp  = ((2.15*(Wn^2)*M)-K);
Gc  = pid(Kp, Ki, Kd)

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

  with Kp = 1.43e+07, Ki = 1.37e+09, Kd = 1.97e+05
 
Continuous-time PID controller in parallel form.

%% Closed-loop System (without scaling factor)
Gcl = minreal(feedback(Gc*Gp, 1))

Gcl =
 
     2.399e09 s^2 + 1.741e11 s + 1.672e13
  ------------------------------------------
  s^3 + 2.399e09 s^2 + 1.741e11 s + 1.672e13
 
Continuous-time transfer function.

pole(Gcl)   % locations of actual closed-loop poles

ans = 
   1.0e+09 *

  -2.3989 + 0.0000i
  -0.0000 + 0.0000i
  -0.0000 - 0.0000i

%% Closed-loop System (with scaling factor)
Gcl = minreal(feedback(1e-9*Gc*Gp, 1))

Gcl =
 
      2.399 s^2 + 174.1 s + 1.672e04
  ---------------------------------------
  s^3 + 2.497 s^2 + 1.236e04 s + 1.672e04
 
Continuous-time transfer function.

S   = stepinfo(Gcl)

S = struct with fields:
         RiseTime: 1.6022
    TransientTime: 2.9365
     SettlingTime: 2.9365
      SettlingMin: 0.9004
      SettlingMax: 1.0000
        Overshoot: 0.0036
       Undershoot: 0
             Peak: 1.0000
         PeakTime: 5.2782

figure(2)
step(Gcl), grid on, grid minor

[Wynand]

Thank you for the response. The transfer function has a K as the numerator but I have been instructed to use a 1 as the numerator. In which case the step response is an attenuating oscillation around -1

[Wynand]

Good morning @Sam Chak

Thank you for your insights.

I have re-evaluated over and over but I am missing a few steps it seems

syms m b k kp ki kd s

G = k/(m*s^2+b*s+k)
Gc = (kp*s + ki*s^2 + kd)/s
T = (k*(kp*s + ki*s^2 + kd))/((m*s^3+b*s^2+k*s+kp*s + ki*s^2 + kd))

From this and the Ideal ITAE 3rd order transfer function I get the following expressions for Kp, Ki and Kd

Zeta = 1; 
Ts = 5;
Wn = 4/(Zeta*Ts);
M = 1019;
B = 100;
K = 12421875;

Kp = (2.15)*(Wn^2)*(M/K)-1

Kp = -0.9999

Ki = (Wn^3)*(M/K)

Ki = 4.2001e-05

Kd = ((1.75)*(M)*(Wn)-B)*(1/K)

Kd = 1.0680e-04

Am I making a mistake somewhere or did I misinterpret your previous answers?