Tuning a PID controller to desired Overshoot, settling time and rise time

The uncompensated plant transfer function in this case is: as shown in the simulink model below:

I want to tune the PID so that the Rise time becomes 0.561 s, settling time becomes 2.6 s and overshoot narrows to 8.83%. However the PID tuner app is very inconvenient in this regard as it only varies the response time and transient behaviour, it is impossible to achieve the desired settling time and overshoot with such a PID tuner.

Answer 1

[Sam Chak]

Hi @Muhammad

I'm revisiting a control problem we encountered a year ago. Previously, we utilized the systune() method to establish PIDF gains, but the transient response exhibited significant oscillations, and the derivative term appeared unusual with an unexpectedly large time constant Tf. This time, I've discovered that a set of classical PID gains can fulfill the requirements: settling time (ts) less than 2.6, overshoot (os) of 8.83%, and a rise time (tr) of 0.561.

%% Original Plant
Gp  = tf(20, [1 4.5 64]);

%% Target Reference
Gref= tf(64.9193300346851, [1, 7.93016824066555, 31.8718356679018, 64.9193300346851]);

%% Classical PID
kp  = -0.733860787658412;
ki  =  6.597481064703570;
kd  =  0.100288660665148;
Gc  = pid(kp, ki, kd);

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

Gcl =
 
     2.006 s^2 - 14.68 s + 131.9
  ---------------------------------
  s^3 + 6.506 s^2 + 49.32 s + 131.9
 
Continuous-time transfer function.

Scl = stepinfo(Gcl)

Scl = struct with fields:
         RiseTime: 0.5610
    TransientTime: 2.0757
     SettlingTime: 2.0757
      SettlingMin: 0.9003
      SettlingMax: 1.0883
        Overshoot: 8.8301
       Undershoot: 0
             Peak: 1.0883
         PeakTime: 0.9200

stepplot(Gref), hold on
stepplot(Gcl),  grid on
stepplot(Gp),   hold off
legend('Target System', 'PID Controller', 'Original Plant', 'location', 'east')

[Emmanouil Tzorakoleftherakis]

Hi @Sam Chak,

For completeness, can you maybe share how you reached to those pid gains so that the accepted answer is complete?

Thanks

[Sam Chak]

Hi @Emmanouil Tzorakoleftherakis

Thank you for accepting my answer on behalf of the OP. The PID gains were obtained by solving the optimization problem using fmincon(). After revisiting this problem a year later, I adjusted the weights of the cost function and found a solution that yields slightly improved step performance.

For linear 2nd-order systems, achieving the desired settling time does not satisfy the rise time requirement. Thus, the rise time requirement takes precedence over the settling time. In other words, satisfying the rise time will satisfy the settling time.

For some unknown reason, systune() yields a less satisfactory solution, characterized by an oscillatory transient response, despite the TuningGoal.StepTracking being properly specified using the tau method, as demonstrated by @Paul. It is believed that systune() may be tuning the PID gains to meet internal objectives that are not known to users, aiming to achieve a balance between performance and robustness, similar to the behavior of pidtune().

%% Target Reference
zeta= 0.611315066227871;    % desired damping ratio that yields target 8.83% overshoot

wn1 = 2.29352925249462;     % natural frequency that yields target 2.6 s settling time
Gt1 = tf(wn1^2, [1, 2*zeta*wn1, wn1^2])

Gt1 =
 
          5.26
  --------------------
  s^2 + 2.804 s + 5.26
 
Continuous-time transfer function.

wn2 = 3.35700150674195;     % natural frequency that yields target 0.561 s rise time
Gt2 = tf(wn2^2, [1, 2*zeta*wn2, wn2^2]);

%% Find optimal PID gains using fmincon
[K, fval] = fmincon(@costfcn, [-0.7 6.6 0.1], [], [])

Warning: Simulation did not reach steady state. Please specify YFINAL if this system is stable and eventually settles.

Warning: Simulation did not reach steady state. Please specify YFINAL if this system is stable and eventually settles.

Warning: Simulation did not reach steady state. Please specify YFINAL if this system is stable and eventually settles.

Warning: Simulation did not reach steady state. Please specify YFINAL if this system is stable and eventually settles.

Warning: Simulation did not reach steady state. Please specify YFINAL if this system is stable and eventually settles.

Warning: Simulation did not reach steady state. Please specify YFINAL if this system is stable and eventually settles.

Warning: Simulation did not reach steady state. Please specify YFINAL if this system is stable and eventually settles.

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

K = 1×3

-0.4161 8.8316 0.1307

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fval = 4.8215e-13

%% Original Plant
Gp  = tf(20, [1 4.5 64])

Gp =
 
         20
  ----------------
  s^2 + 4.5 s + 64
 
Continuous-time transfer function.

%% Classical PID
Kp  = K(1);
Ki  = K(2);
Kd  = K(3);
Gc  = pid(Kp, Ki, Kd)

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

  with Kp = -0.416, Ki = 8.83, Kd = 0.131
 
Continuous-time PID controller in parallel form.

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

Gcl =
 
     2.614 s^2 - 8.322 s + 176.6
  ---------------------------------
  s^3 + 7.114 s^2 + 55.68 s + 176.6
 
Continuous-time transfer function.

%% Plot results
hold on
stepplot(Gt1)
stepplot(Gt2)
stepplot(Gp)
stepplot(Gcl) 
hold off
grid on
legend('Target Settling Time', 'Target Rise Time', 'Original Plant', 'PID Controller', 'location', 'east')

%% stepInfoTable
S1  = stepinfo(Gt1);
S2  = stepinfo(Gt2);
S3  = stepinfo(Gp);
S4  = stepinfo(Gcl);

stepInfoTable = struct2table([S1, S2, S3, S4]);
stepInfoTable = removevars(stepInfoTable, {'TransientTime', 'SettlingMin', 'SettlingMax', 'Undershoot', 'Peak', 'PeakTime'});
stepInfoTable.Properties.RowNames = {'Gt1', 'Gt2', 'Gp', 'Gcl'};
disp(stepInfoTable);

           RiseTime    SettlingTime    Overshoot
           ________    ____________    _________

    Gt1    0.82113           2.6          8.83  
    Gt2      0.561        1.7763          8.83  
    Gp     0.16265        1.7235        39.822  
    Gcl      0.561        1.9266          8.83  

%% Cost function
function J = costfcn(param)
    Gp  = tf(20, [1 4.5 64]);
    kp  = param(1);
    ki  = param(2);
    kd  = param(3);
    Gc  = pid(kp, ki, kd);
    Gcl = feedback(Gc*Gp, 1);
    S   = stepinfo(Gcl);
    tr  = S.RiseTime;
    os  = S.Overshoot;
    Tr  = 0.561;    % desired Rise time
    Os  = 8.83;     % desired Overshoot
    J   = 100*(tr - Tr)^2 + 50*(os - Os)^2;
end

Answer 2

[Sam Chak]

Hi @Muhammad

This is the Root Locus compensator design approach, first attempt.

s   = tf('s');
Gp  = 20/(s^2 + 4.5*s + 64)

Gp =
 
         20
  ----------------
  s^2 + 4.5 s + 64
 
Continuous-time transfer function.

controlSystemDesigner('rlocus', Gp)

The app will show root locus of .

Fig. 1: Right-click the white area and add a new design requirement.

Fig. 2: Select Percent overshoot and insert the design value.

Fig. 3: An exclusion region is indicated by the yellow shaded area, separated from the white area by the thick black lines.

Fig. 4: Add a second design requirement, select Settling time and insert the design value.

Fig. 5: A second exclusion region overlaps with the first exclusion region. A new thick black line can also be seen.

Fig. 6: Right-click the white area, click on Edit compensator, and then add an Integrator, because the Plant is a Type-0 system.

Fig. 7: Right-click the white area, click on Edit compensator, and then add a Real Pole, placing it relatively far away from the Plant's poles

Fig. 8: Right-click the white area, click on Edit compensator, and then add the Complex Zeros to cancel out the Plant's stable Complex poles. Natural frequency , and Damping ratio .

Note: This doesn't work if the Plant is unstable.

Fig. 9: The "magenta dot markers" can be seen. The Plant's stable poles are cancelled out by the Compensator's zeros. Next, the design task is to drag one of the magenta markers (on the real axis) along the root locus until a response plot stays outside of the associated exclusion regions.

Fig. 10: The "magenta dot markers" are dragged to the breakaway point (My preference).

Fig. 11: This can also be directly adjusted through entering a design value to Compensator gain.

Fig. 12: Check the Step response if the performance is satisfactory.

Fig. 13: If satisfactory, then export the designed Compensator to Workspace. With the design values, they can be entered in the Transfer function block or zpk block in Simulink.

C   = zpk(1.8*(s^2 + 4.5*s + 64)/(s*(s + 12)))

C =
 
  1.8 (s^2 + 4.5s + 64)
  ---------------------
        s (s+12)
 
Continuous-time zero/pole/gain model.

Gcl = tf(feedback(C*Gp, 1))

Gcl =
 
           36 s^2 + 162 s + 2304
  ---------------------------------------
  s^4 + 16.5 s^3 + 154 s^2 + 930 s + 2304
 
Continuous-time transfer function.

% Can check in MATLAB whether the design requirements are satisfied.
S   = stepinfo(Gcl)

S = struct with fields:
         RiseTime: 0.5599
    TransientTime: 0.9725
     SettlingTime: 0.9725
      SettlingMin: 0.9055
      SettlingMax: 1.0000
        Overshoot: 0
       Undershoot: 0
             Peak: 1.0000
         PeakTime: 2.1337

If satisfactory, find the equivalent PID gains via algebraic calculations.

% This step is unnecessary, because the Compensator works just as good!
kp   = -0.125;
ki   =  9.6;
kd   =  77/480;
Tf   =  1/12;
Cpid = pid(kp, ki, kd, Tf)

Gpid =
 
             1            s    
  Kp + Ki * --- + Kd * --------
             s          Tf*s+1 

  with Kp = -0.125, Ki = 9.6, Kd = 0.16, Tf = 0.0833
 
Continuous-time PIDF controller in parallel form.

% Proof
zpk(Cpid)

ans =
 
  1.8 (s^2 + 4.5s + 64)
  ---------------------
        s (s+12)
 
Continuous-time zero/pole/gain model.

[Sam Chak]

Hi @Muhammad

If you want to stay at the borderline design requirements, then pick these PID gain values. You should also try out and learn @Paul's tuning method. But you need to set tau , instead of .

Gp = tf(20, [1 4.5 64])

Gp =
 
         20
  ----------------
  s^2 + 4.5 s + 64
 
Continuous-time transfer function.

% PID Controller (manual via 
kp  = -1.522560;
ki  =  8.787510;
kd  =  0.508165;
N   =  4.105485;
Tf  =  1/N;
Gc  = pid(kp, ki, kd, Tf)

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

  with Kp = -1.52, Ki = 8.79, Kd = 0.508, Tf = 0.244
 
Continuous-time PIDF controller in parallel form.

Gcl = feedback(Gc*Gp, 1)

Gcl =
 
           11.27 s^2 + 50.73 s + 721.5
  ---------------------------------------------
  s^4 + 8.605 s^3 + 93.75 s^2 + 313.5 s + 721.5
 
Continuous-time transfer function.

% Desired Rise time 0.561 sec takes precedence over the desired Settling time 2.6 sec
% Desired Percent overshoot 8.83% is also met.
S   = stepinfo(Gcl)

S = struct with fields:
         RiseTime: 0.5609
    TransientTime: 1.7760
     SettlingTime: 1.7760
      SettlingMin: 0.9087
      SettlingMax: 1.0883
        Overshoot: 8.8278
       Undershoot: 0
             Peak: 1.0883
         PeakTime: 1.1890

step(Gcl, 5)

In Simulink, the following P, I, D, N values are entered into the PID controller block.

zpk(Gc)

ans =
 
  0.5637 (s^2 + 4.5s + 64)
  ------------------------
        s (s+4.105)
 
Continuous-time zero/pole/gain model.

p = pole(minreal(Gcl, 1e-4))

p = 
  -2.0527 + 2.6571i
  -2.0527 - 2.6571i

The root locus of the compensated system looks like this. As expected, because the compensator only meets the borderline design requirement, the loci of the closed-loop poles ( marked by the magenta squares) are situated at the boundary lines separating the white and yellow shaded areas.

[Paul]

Hi Sam,

I've been playing around a bit with systune and so far haven't been able to get an acceptable result for what should be a simple design. Still need to learn more about how to use that tool.

Do you think the compensator with a pole at s = 4.105 can still the be considered PID? Is there a guideline on how to pick Tf so that Gc still provides the derivative action at frequencies where it's needed but limits the gain at high freuqencies? Maybe something like 1/Tf > 10*wc, where wc is the desired open-loop crossover frequency.

[Sam Chak]

Hi @Paul,

Edit: I'm learning Tuning Goals from your example. Let's say the desired closed-loop transfer function is known.

G = tf(11.27405828, [1 4.10548454 11.27405828])

G =
 
          11.27
  ---------------------
  s^2 + 4.105 s + 11.27
 
Continuous-time transfer function.

pole(G)

ans = 
  -2.0527 + 2.6571i
  -2.0527 - 2.6571i

How should I insert location of the desired poles in this syntax?

Req = TuningGoal.Poles(location, mindecay, mindamping, maxfreq)

https://www.mathworks.com/help/control/ref/tuninggoal.poles-class.html#btjqa2o-1-location

*************************************************

Regarding the derivative filter I learned that If I design an ideal PID (without derivative filter)

,

then the filter time constant in the practical PID (with derivative filter)

,

should be selected to satisfy this condition

.

To be honest, I don't treat PID (with derivative filter) as the typical P + I + D thing, but more as a special form of compensator, which has an equivalent transfer function of

something that is similar to what is shown here:

https://www.mathworks.com/help/slcontrol/ug/when-tuning-the-pid-controller-the-d-gain-has-a-different-sign-from-the-i-gain.html

With a 2nd-order Compensator and a 2nd-order Plant, the closed-loop transfer function should be a 4th-order system without Pole–Zero cancellation. If applying the conventional design approach taught in textbooks, then we should select 4th-order system that gives the desired overshoot .

G1  = tf(1, [1 2.6784513 3.7935254 2.6784513 1])

G1 =
 
                      1
  -----------------------------------------
  s^4 + 2.678 s^3 + 3.794 s^2 + 2.678 s + 1
 
Continuous-time transfer function.

S1  = stepinfo(G1)

S1 = struct with fields:
         RiseTime: 2.6359
    TransientTime: 8.2598
     SettlingTime: 8.2598
      SettlingMin: 0.9155
      SettlingMax: 1.0883
        Overshoot: 8.8279
       Undershoot: 0
             Peak: 1.0883
         PeakTime: 6.0531

Rise Time takes precedence over the Settling Time. The desired Rise Time sec can be achieved by selecting a proper frequency .

wn  = 4.6984048;
G2  = tf(wn^4, [1 2.6784513*wn 3.7935254*wn^2 2.6784513*wn^3 wn^4]);
S2  = stepinfo(G2)

S2 = struct with fields:
         RiseTime: 0.5610
    TransientTime: 1.7580
     SettlingTime: 1.7580
      SettlingMin: 0.9155
      SettlingMax: 1.0883
        Overshoot: 8.8279
       Undershoot: 0
             Peak: 1.0883
         PeakTime: 1.2883

However, if is selected, the compensated system will have zeros on the right-half s-plane and large undershoot. Skipping the tuning process, is selected instead.

wn  = 7.5089051;
G3  = tf(wn^4, [1 2.6784513*wn 3.7935254*wn^2 2.6784513*wn^3 wn^4])

G3 =
 
                     3179
  -------------------------------------------
  s^4 + 20.11 s^3 + 213.9 s^2 + 1134 s + 3179
 
Continuous-time transfer function.

S3  = stepinfo(G3)

S3 = struct with fields:
         RiseTime: 0.3510
    TransientTime: 1.1000
     SettlingTime: 1.1000
      SettlingMin: 0.9155
      SettlingMax: 1.0883
        Overshoot: 8.8279
       Undershoot: 0
             Peak: 1.0883
         PeakTime: 0.8061

Applying the Pole Placement method (manual algebra calculation as I'm not good at setting TuningGoal and using systune) based on and converting the 2nd-order compensator to PID form:

% Practical PID controller 
kp  = -0.2203631;
ki  = 10.1815615;
kd  = 0.26916585;
Tf  = 0.06405211;
Gc  = pid(kp, ki, kd, Tf)

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

  with Kp = -0.22, Ki = 10.2, Kd = 0.269, Tf = 0.0641
 
Continuous-time PIDF controller in parallel form.

% Closed-loop system
s   = tf('s');
Gp  = 20/(s^2 + 4.5*s + 64);    % Plant
Gcl = minreal(feedback(Gc*Gp, 1))

Gcl =
 
          79.64 s^2 + 134.8 s + 3179
  -------------------------------------------
  s^4 + 20.11 s^3 + 213.9 s^2 + 1134 s + 3179
 
Continuous-time transfer function.

Designing a Prefilter to cancel out the zeros on the left-half s-plane.

% Prefilter
Gf  = tf(39.919688, [1 1.6929442 39.919688])

Gf =
 
          39.92
  ---------------------
  s^2 + 1.693 s + 39.92
 
Continuous-time transfer function.

% Step Response
Gclf = minreal(Gf*Gcl, 1e-4)    % Should behave almost similar to G3

Gclf =
 
                     3179
  -------------------------------------------
  s^4 + 20.11 s^3 + 213.9 s^2 + 1134 s + 3179
 
Continuous-time transfer function.

step(Gf*Gcl, 3)

stepinfo(Gf*Gcl)

ans = struct with fields:
         RiseTime: 0.3510
    TransientTime: 1.1000
     SettlingTime: 1.1000
      SettlingMin: 0.9153
      SettlingMax: 1.0883
        Overshoot: 8.8279
       Undershoot: 0
             Peak: 1.0883
         PeakTime: 0.8060

[Paul]

Hi Sam,

I haven't been looking at this problem too closely, still wasn't able to get systune to give me what I wanted, but I was getting closer. Anyway, there are several ways to insert an Analysis Point into a model. See the link from the doc page you linked for some examples, maybe other doc pages as well. Here's one way, using connect

G = tf(20,[1 4.5 64],'InputName','u','OutputName','y');
K = tunablePID('PID','PID');
K.InputName = 'e';
K.OutputName = 'u';
% Create the closed loop system, include an analysis point at the error
% signal
S = sumblk('e = r - y');
CL0 = connect(G,K,S,'r','y','e')

CL0 =

  Generalized continuous-time state-space model with 1 outputs, 1 inputs, 3 states, and the following blocks:
    AnalysisPoints_: Analysis point, 1 channels, 1 occurrences.
    PID: Tunable PID controller, 1 occurrences.

Type "ss(CL0)" to see the current value, "get(CL0)" to see all properties, and "CL0.Blocks" to interact with the blocks.

Then, in systune, the "location" input would 'e'

[Sam Chak]

Thank you, @Paul. Will look into it.

Answer 3

[Sulaymon Eshkabilov]

Simply double-click on PID block and its parameters window opens. Click on "Tune" button as shown in the attached screen shot - Step 1 - PID_tune_step1.jpg.

Then follow step 2 as shown in the attached file: PID_tune_step2.jpg

Alternative way is using pidtune() in MATLAB.

Good luck

Answer 4

[Paul]

The CST offers some tools to do this. Start from this doc page. Or from this doc page if you want to use Simulink.

Here is an example of the command line workflow. The result I got is not that good, so you'll have to make it work better to get the result you want

Define the plant model

G = tf(20,[1 4.5 64],'InputName','u','OutputName','y');

Define a tunable PID compensator

K = tunablePID('PID','PID');
K.InputName = 'e';
K.OutputName = 'u';

Create the closed loop system

S = sumblk('e = r - y');
CL0 = connect(G,K,S,'r','y');

Define a second order model response that meets the requirements for rise time, settling time, and overshoot.

Rtrack = TuningGoal.StepTracking('r','y',.30,8.8);
stepinfo(Rtrack.ReferenceModel)

ans = struct with fields:
         RiseTime: 0.5654
    TransientTime: 1.7893
     SettlingTime: 1.7893
      SettlingMin: 0.9077
      SettlingMax: 1.0880
        Overshoot: 8.7988
       Undershoot: 0
             Peak: 1.0880
         PeakTime: 1.1967

Use systune to design the PID gains

CL = systune(CL0,Rtrack);

Final: Soft = 3.09, Hard = -Inf, Iterations = 44

The PID controller is:

showTunable(CL)

PID =
 
             1            s    
  Kp + Ki * --- + Kd * --------
             s          Tf*s+1 

  with Kp = -0.101, Ki = 6.45, Kd = -23.9, Tf = 1.37e+04
 
Name: PID
Continuous-time PIDF controller in parallel form.

Those negative gains look peculiar!

The closed-loop step response is

stepplot(CL)

stepinfo(CL)

ans = struct with fields:
         RiseTime: 0.4233
    TransientTime: 2.0735
     SettlingTime: 2.0735
      SettlingMin: 0.8039
      SettlingMax: 1.0237
        Overshoot: 2.3675
       Undershoot: 0.0210
             Peak: 1.0237
         PeakTime: 1.5024

Presumably there are other Tuning Goals that can be applied to get the desired response.

[Hansi K]

But how do you set a tuning goal for the rise time?

TuningGoal.StepTracking('r','y',.30,8.8);

Is it about tau? 0.30 in the line above? I've already read the reference to tau but I just don't understand if there's a connection between the tau and the rise time.

[Paul]

According to the linked doc page, the reference model is a second order system with tau = 1/wn.

For a second order model, the synthesis equation for wn is

wn >= 1.8/tr (tr is the rise time requirement)

so

1/tau >= 1.8/0.561

tau <= 0.561/1.8 = 0.31

That's probably how I got the tau = 0.3 in the TuningGoal. Can't say for sure.

[Sam Chak]

Hi @Hansi K,

Given the lack of a clear description of the rise time by the OP, the definition of rise time appears to be somewhat subjective. MATLAB typically measures rise time from 10% to 90%. In contrast, many control theory textbooks use a formula defining it as the time needed for the response to rise from 0% to 100% of its final value.

%% Parameters
zeta = 0.61131525436442019; % Damping ratio that gives OS 8.83
tr   = 0.561;               % Rise time from 0 to 100%

%% Natural frequency of a linear second order system
wn   = ((1/sqrt(1 - zeta^2))*(pi - atan(sqrt(1 - zeta^2)/zeta)))/tr;
tau  = 1/wn                 % this will be the tau value I'd probably take

tau = 0.1992

%% 2nd-order transfer function that yields Rise time 0.561 s & Overshoot 8.83%
G    = tf(wn^2, [1 2*zeta*wn wn^2]) 

G =
 
          25.2
  --------------------
  s^2 + 6.137 s + 25.2
 
Continuous-time transfer function.

S    = stepinfo(G)

S = struct with fields:
         RiseTime: 0.3752
    TransientTime: 1.1880
     SettlingTime: 1.1880
      SettlingMin: 0.9087
      SettlingMax: 1.0883
        Overshoot: 8.8300
       Undershoot: 0
             Peak: 1.0883
         PeakTime: 0.7954

step(G), grid on
xline(tr, '--', sprintf('Time from 0 to 100%%: %.3f sec', tr), 'color', '#7F7F7F', 'LabelVerticalAlignment', 'bottom')

Upon analyzing the issue, it is evident that the OP aims to slow down the fast-responsive plant and reduce the overshoot, hence the inclusion of a negative sign in the proportional gain in the auto PID tuning. However, this poses an issue in the pole placement design when incorporating the Prefilter.

Gp  = tf(20, [1 4.5 64])

Gp =
 
         20
  ----------------
  s^2 + 4.5 s + 64
 
Continuous-time transfer function.

stepinfo(Gp)

ans = struct with fields:
         RiseTime: 0.1627
    TransientTime: 1.7235
     SettlingTime: 1.7235
      SettlingMin: 0.2629
      SettlingMax: 0.4369
        Overshoot: 39.8224
       Undershoot: 0
             Peak: 0.4369
         PeakTime: 0.4093

[numP, denP] = tfdata(Gp, 'v');
a1  = denP(2);
a2  = denP(3);
b   = numP(3);

%% Target Compensated System that yields Rise time 0.561 s (10-90%) & Overshoot 8.83%
G   = tf(64.9193300346851, [1, 7.93016824066555, 31.8718356679018, 64.9193300346851])

G =
 
               64.92
  --------------------------------
  s^3 + 7.93 s^2 + 31.87 s + 64.92
 
Continuous-time transfer function.

[numG, denG] = tfdata(G, 'v');
d1  = denG(2);
d2  = denG(3);
d3  = denG(4);

%% Classical PID using pole placement (but it doesn't affect zero)
kp  = -(a2 - d2)/b;
ki  = d3/b;
kd  = -(a1 - d1)/b;
Gc  = pid(kp, ki, kd)

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

  with Kp = -1.61, Ki = 3.25, Kd = 0.172
 
Continuous-time PID controller in parallel form.

%% Compensated system (with 2 zeros on the right-half plane)
Gcl = feedback(Gc*Gp, 1)

Gcl =
 
     3.43 s^2 - 32.13 s + 64.92
  --------------------------------
  s^3 + 7.93 s^2 + 31.87 s + 64.92
 
Continuous-time transfer function.

It is commonly taught that a Prefilter (Gf) can eliminate undesired zeros in the closed-loop system (Gcl). Nevertheless, some less-experienced students may unintentionally introduce unstable poles while designing the Prefilter. Mathematically, it is observed that Gcl's numerator is miraculously nullified by Gf's denominator, resulting in the Target Compensated System.

How would you design the controller instead?

%% Pre-filter
Gf  = tf(b*ki, [b*kd b*kp b*ki])

Gf =
 
            64.92
  --------------------------
  3.43 s^2 - 32.13 s + 64.92
 
Continuous-time transfer function.

Gclf= minreal(Gf*Gcl)

Gclf =
 
               64.92
  --------------------------------
  s^3 + 7.93 s^2 + 31.87 s + 64.92
 
Continuous-time transfer function.

stepinfo(Gclf)

ans = struct with fields:
         RiseTime: 0.5610
    TransientTime: 1.6437
     SettlingTime: 1.6437
      SettlingMin: 0.9152
      SettlingMax: 1.0883
        Overshoot: 8.8333
       Undershoot: 0
             Peak: 1.0883
         PeakTime: 1.2146

%% Wishful step response
step(Gclf),  grid on    

[Sam Chak]

@Hansi K, the systune() command is employed to identify the tunable parameters in the PID controller, optimizing them to best meet the requirements of the target step response that yields Rise time 0.561 and Overshoot 8.83%.

%% Original Plant
Gp  = tf(20, [1 4.5 64])

Gp =
 
         20
  ----------------
  s^2 + 4.5 s + 64
 
Continuous-time transfer function.

Gc0 = tunablePID('Gc', 'PID');

%% Initial Closed-loop system (untuned)
CL0 = feedback(Gc0*Gp, 1);
CL0.InputName  = 'r';
CL0.OutputName = 'y';

%% Reference system for target step response (Rise time 0.561 & Overshoot 8.83%)
Gref= tf(64.9193300346851, [1, 7.93016824066555, 31.8718356679018, 64.9193300346851]);

%% Tuning goal 
Req = TuningGoal.StepTracking('r', 'y', Gref);
stepinfo(Req.ReferenceModel)

ans = struct with fields:
         RiseTime: 0.5610
    TransientTime: 1.6437
     SettlingTime: 1.6437
      SettlingMin: 0.9152
      SettlingMax: 1.0883
        Overshoot: 8.8333
       Undershoot: 0
             Peak: 1.0883
         PeakTime: 1.2146

%% Tune CL0 according to tuning goals specified in Req
CL = systune(CL0, Req);

Final: Soft = 2.82, Hard = -Inf, Iterations = 48

%% Display the values of all tunable Control parameters
showTunable(CL)

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

  with Kp = -0.437, Ki = 5.23, Kd = -27.2, Tf = 2.67e+04
 
Name: Gc
Continuous-time PIDF controller in parallel form.

%% Show performances of tuned Closed-loop system
stepinfo(CL)

ans = struct with fields:
         RiseTime: 0.4516
    TransientTime: 2.2666
     SettlingTime: 2.2692
      SettlingMin: 0.8034
      SettlingMax: 1.0123
        Overshoot: 1.2300
       Undershoot: 2.4392
             Peak: 1.0123
         PeakTime: 1.6758

stepplot(Gref), hold on
stepplot(CL),   grid on
stepplot(Gp),   hold off
legend('Target System', 'Compensated', 'Original Plant', 'location', 'east')

[Hansi K]

Thx for the explanation! I think I understand it much better now. At least for this case xD

Is there a general formula for tau (or wn) when the rise time is not measured from 0% to 100%? I think my uni script uses a very strange definition of tr... Or at least another definition than MATLAB.

'Wendepunkt' means 'inflection point'

[Sam Chak]

Hi @Hansi K

I'm not entirely certain about deriving the formula for when the rise time isn't measured from 0 to 100% of its final value. However, @Paul has found the inequality quite useful. If the Rise time < requirement is loosely defined, this formula provides for the exact measured from 0 to 100%:

.

In simpler terms, it ensures that will yield the Rise time < when measured from 10% to 90% or 5% to 95%. The Wendepunkt-based rise time seems to approximate the rise time measured from 10% to 90%.

The control problem presented by the OP is somewhat unique. In my opinion, the systune() tool proves to be more versatile than the pidtune() tool, especially when the Control Designer can effectively express the Reference system for a desired step response.

The pidtune() is better suited for cases, where the 0 dB gain crossover frequency () and target phase margin (pm) are precisely specified. Nonetheless, a challenge arises as not many users are familiar with the process of converting desired rise time () or settling time () into and target phase margin. Maybe there are, but I don't remember seeing related examples in the documentation.

Here is an example illustrating how the powerful systune() tool can design the PIDF controller to compensate for a slower plant, yielding an output that precisely mirrors the Reference system for a desired step response.

%% Original Plant (2nd-order system with ts ≈ 8 and os ≈ 16%)
Gp  = tf(1, [1 1 1])

Gp =
 
       1
  -----------
  s^2 + s + 1
 
Continuous-time transfer function.

stepinfo(Gp)

ans = struct with fields:
         RiseTime: 1.6390
    TransientTime: 8.0759
     SettlingTime: 8.0759
      SettlingMin: 0.9315
      SettlingMax: 1.1629
        Overshoot: 16.2929
       Undershoot: 0
             Peak: 1.1629
         PeakTime: 3.5920

%% Set tunable parameters in the Compensator Gc
Gc0 = tunablePID('Gc', 'PID');

%% Initial Closed-loop system (untuned)
CL0 = feedback(Gc0*Gp, 1);
CL0.InputName  = 'r';
CL0.OutputName = 'y';

%% Reference system for target step response (ts < 1 & No Overshoot)
wn   = 6;   % this value gives Settling time under 1 sec
zeta = 1;   % damping ratio 1 gives zero overshoot
Gref = tf(wn^2, [1, 2*zeta*wn, wn^2])   % I want the output to look like this!

Gref =
 
        36
  ---------------
  s^2 + 12 s + 36
 
Continuous-time transfer function.

%% Tuning goal 
Req = TuningGoal.StepTracking('r', 'y', Gref);
stepinfo(Req.ReferenceModel)

ans = struct with fields:
         RiseTime: 0.5597
    TransientTime: 0.9724
     SettlingTime: 0.9724
      SettlingMin: 0.9019
      SettlingMax: 0.9999
        Overshoot: 0
       Undershoot: 0
             Peak: 0.9999
         PeakTime: 1.9879

%% Tune CL0 according to tuning goals specified in Req
CL    = systune(CL0, Req);

Final: Soft = 7.96e-13, Hard = -Inf, Iterations = 96

CLsys = getIOTransfer(CL, 'r', 'y');
Gcl   = tf(CLsys)

Gcl =
 
  From input "r" to output "y":
         36 s^2 + 36 s + 36
  ---------------------------------
  s^4 + 13 s^3 + 49 s^2 + 48 s + 36
 
Continuous-time transfer function.

%% Display the values of all tunable Control parameters
showTunable(CL)

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

  with Kp = 2.75, Ki = 3, Kd = 2.77, Tf = 0.0833
 
Name: Gc
Continuous-time PIDF controller in parallel form.

%% Show performances of tuned Closed-loop system
stepinfo(CL)

ans = struct with fields:
         RiseTime: 0.5597
    TransientTime: 0.9724
     SettlingTime: 0.9724
      SettlingMin: 0.9019
      SettlingMax: 0.9999
        Overshoot: 0
       Undershoot: 0
             Peak: 0.9999
         PeakTime: 1.9879

%% Plot result and comparison
subplot(211)
stepplot(Gref, 12), grid on
ylim([0 1.2])
legend('Target System', 'location', 'East')

subplot(212)
stepplot(CL), hold on, grid on
ylim([0 1.2])
stepplot(Gp), hold off
legend('PID Controller', 'Original Plant', 'location', 'East')

[Paul]

IIRC, the approxiation tr ~ 1.8/wn applies for "intermediate" values of zeta. In this problem, the specified overshoot was 8%, which implies a zeta of around 0.62 (which sounds like it's intermediate), which happens to correspond very closes to wn = 1.8/tr.

Here's the plot (assuming wn = 5, the plot should be agnostic to the value of wn)

zeta = 0.1:.1:0.9;  wn = 5;
for ii = 1:numel(zeta)
    S = stepinfo(tf(wn^2,[1 2*zeta(ii)*wn wn^2]));
    tr(ii) = S.RiseTime;
    Mp(ii) = S.Peak - 1;
end

figure
plot(zeta,tr*wn,zeta,Mp)
xlabel('zeta')
yline([0.08 1.8])
xline(0.62)
legend('10-90 Rise Time*wn','%OverShoot/100','Location','NorthWest')

[Senol Gulgonul]

Recently I published a preprint on PID tuning for settling time and overshoot requirements by using IAE optimization at https://arxiv.org/abs/2505.17268

By using my algorithm I reached zero overshoot with Ts=2.26s

Method Kp Ki Kd Ts PO (%) IAE

SOSTIAE 0.0000 4.9777 0.1278 2.2667 0.0000 0.6429

you can set your target Ts and PO