Unable to tune a pitch-damper autopilot gains with looptune

Hello everyone,

I am working with a non-linear 3DoF model of a buisness jet aircraft. I have implemented the open-loop dynamics correctly in Simulink, generating the operating point for a given horizontal steady flight. You can see a glimpse of the main block, with its inputs and outputs.

Now, I have to implement an autopilot. In particular, I am trying to implement the TECS (Total Energy Control System). This control architecture has an outerloop indipendent from the aircraft dynamics and an innerloop which instead depends on the aircraft dynamics. The outputs of the outerloop will be a commanded pitch angle and a commanded thrust, so I have to implement the innerloops giving these inputs.

Following this guide Tuning of Gain-Scheduled Three-Loop Autopilot - MATLAB & Simulink - MathWorks Italia, in order to implement the autopilot, I have prepared ad second model in which the state propagator recives the linearized state around the same trim condition from the open-loop model. The code is the following

% Linearizing the open loop around trim

% Specify the model name
model = '3DoF_open_loop';
% Linearize 
sys = linearize(model,'3DoF_open_loop/State Propagator',op);

% Swapping blocks
BlockSubs = struct('Name',...
   'Closed_loop_linearized_AC3DoF_model/Linearized State Propagator','Value', sys);

% Loading the trimmed values of the inputs from the open-loop to closed
% loop
in_cl = Simulink.SimulationInput('Closed_loop_linearized_AC3DoF_model');
in_cl = in_cl.setExternalInput(inputs);

in_cl.applyToModel;

and this is a glimpse of the model. The linearized state propagator state was initialized manually with the trim values of theta and alpha previously computed.

The thrust controller was preatty easy to implement and tune with a PI controller block

Now it is the turn of the pitch-damper autopilot

that I have implemented in this way

For this moment, I am negleting the actuator dynamics. I have tryed to tune the two gains in this way with looptune

ST0 = slTuner('Closed_loop_linearized_AC3DoF_model', {'Kp', 'Kd'}, BlockSubs);

addPoint(ST0,{'theta_cmd','delta_e_cmd','theta_out', 'q_out'});

% Design requirements
wc = [3,12];   % bandwidth
TrackReq = TuningGoal.Tracking('theta_cmd','theta_out',0.300);  % tracking

Controls = 'delta_e_cmd';
Measurements = {'theta_out', 'q_out'};

[ST,gam,Info] = looptune(ST0,Controls,Measurements,wc,TrackReq);

figure();
Tr1 = getIOTransfer(ST,'theta_cmd','theta_out');
step(Tr1,5)
grid

writeBlockValue(ST)

The problem is that the optimization does not converge

Final: Peak gain = 1e+03, Iterations = 1

I have tryed also with different requirements, but I get always the same results. In your opinion, are there any errors because of the implementation of the autopilots? How can I improve the results? Thank you.

[Sam Chak]

Hi @Leonardo Molino,

Could you please provide the mathematical model for the open-loop linearized Aircraft 3-DoF? The linear model should be a SISO (Single Input Single Output) system, which can be represented as either a transfer function or a state-space model. Tuning it in MATLAB would likely be more convenient.

[Leonardo Molino]

Hello @Sam Chak, thank you for your answer.

Yes, I can provide you the LTF from the elevator deflection to the pitch angle. I have determined it through the Model Linearizer app in the open loop Simulink model

Size:  1 inputs, 1 outputs, 5 states

Linearization Result:

  From input "u1" to output "y1":
  
           -11.64 s^3 - 17.38 s^2 - 0.2759 s - 0.0001568
  ---------------------------------------------------------------
  s^5 + 2.167 s^4 + 36.36 s^3 + 0.3578 s^2 + 0.3471 s + 0.0005279
 
Name: Linearization around the operating point "op"
Continuous-time transfer function.

And this is the Impulse response with the phugoid

and the Short period

Answer 1

[Sam Chak]

Hi @Leonardo Molino,

Ogata's formula for the design of the Prefilter is effective, but it only applies when the Compensator is a PID controller in standard form and the Plant has no zeros. It's worth noting that for the Ideal PID controller form in Simulink, the value for need to be calculated algebraically.

... Simulink's PID controller block in Ideal form

Example:

Edit: The code has been updated with a Plant that produces over 20% overshoot in the step response. To eliminate this undesired overshoot while maintaining the desired characteristics of the original plant, the Prefilter and PID controller in standard form can be employed. However, it is important to note that Ogata's formulas are effective only when the tuned values of and are both positive real numbers. Otherwise, the Prefilter will contain unstable poles.

%% Plant, Gp
Gp  = tf(5, [1, 2, 5])          % produces over 20% overshoot

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

%% PID controller in standard form, Gc
Kp  = 1.42976215428325;         % Proportional gain
Ti  = 0.877233144300333;        % Integral constant
Td  = 0.564721629217525;        % Derivative constant
Gc  = pidstd(Kp, Ti, Td)

Gc =
 
             1      1          
  Kp * (1 + ---- * --- + Td * s)
             Ti     s          

  with Kp = 1.43, Ti = 0.877, Td = 0.565
 
Continuous-time PID controller in standard form

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

Gcl =
 
  4.0371 (s^2 + 1.771s + 2.019)
  -----------------------------
           (s+2.012)^3
 
Continuous-time zero/pole/gain model.

%% Prefilter, Gf
Gf  = tf(1, [Ti*Td, Ti, 1])     % Ogata's formula (see A-8-8)

Gf =
 
              1
  -------------------------
  0.4954 s^2 + 0.8772 s + 1
 
Continuous-time transfer function.

%% Command compensated system, Gcc
Gcc = series(Gf, Gcl)

Gcc =
 
    8.1493 (s^2 + 1.771s + 2.019)
  ----------------------------------
  (s+2.012)^3 (s^2 + 1.771s + 2.019)
 
Continuous-time zero/pole/gain model.

Gcc = minreal(Gcc)

Gcc =
 
    8.1493
  -----------
  (s+2.012)^3
 
Continuous-time zero/pole/gain model.

%% Plot results and comparison
stepplot(Gp, 12), hold on
stepplot(Gcl),    grid on
stepplot(Gcc),    hold off
legend('Gp', 'Gcl', 'Gcc', 'location', 'East', 'fontsize', 18)

S1  = stepinfo(Gp);
S2  = stepinfo(Gcl);
S3  = stepinfo(Gcc);
stepInfoTable = struct2table([S1, S2, S3]);
stepInfoTable = removevars(stepInfoTable, {'SettlingMin', 'SettlingMax', 'Undershoot', 'Peak', 'PeakTime'});
stepInfoTable.Properties.RowNames = {'Gp', 'Gcl', 'Gcc'};
stepInfoTable

stepInfoTable = 3×4 table
           RiseTime    TransientTime    SettlingTime    Overshoot
           ________    _____________    ____________    _________

    Gp     0.69034        3.7352           3.7352        20.787  
    Gcl     2.2137        3.5542           3.5542             0  
    Gcc     2.0972        3.7352           3.7352             0  

When the dynamics of the actuator is taken into account...

In the case of reference tracking problems, the dynamics of the actuator can be combined with the plant's dynamics , resulting in . However, if there is an external disturbance at the input of the plant, tuning the Compensator alone won't be effective for disturbance rejection. In such cases, the actuator's dynamics should be grouped with the Compensator Gc, forming . Logically, it's important to keep in mind that the parameters in should be fixed as non-tunable.

[Leonardo Molino]

Dear @Sam Chak I really appreciate you for you help and support.

Since this is my very first time I try to implement an AP (also in theory, not only on MatLab/Simulink), I can tolerate the procedure that you have listed in the previous message.

In particular, I have exploited your 'hint' to tune the AP with a 'target model'. I have also included some actuator dynamics to the whole thing, so I will share to you my guess.

First of all, since I have also to perform some gain-scheduling, I have preferred to work directly with my model. In the following picture, there is the Pitch damper

As you can see, I have introduced the actuator dynamics (the time constant comes from Nelsen's book for buisness jet aircraft) and a saturation block. The code that will look for the gains of the PID is the following

PD0 = slTuner('Closed_loop_linearized_AC3DoF_model','Pitch damper');

% Tuning goal 
Gtr = tf(5, [1, 2, 5])

The tuning goal comes from one of you examples, as you will see the algorithm will not struggle to track this dynamics. The other goal you have mentioned (Gtr = tf(0.04, [1, 0.4, 0.04]) % <-- the settling time of this system is about 30 sec), since is very slow, will produce some oscillatory behaviour in the transient. I don't know why, maybe the actuator dynamics?

S1 = stepinfo(Gtr);

Req = TuningGoal.StepTracking('theta_cmd', 'theta_out', Gtr);
stepinfo(Req.ReferenceModel);

PD  = systune(PD0, Req);

Final: Soft = 2.98, Hard = -Inf, Iterations = 20

Pitch_damper_values =

1 s

Kp + Ki * --- + Kd * --------

s Tf*s+1

with Kp = -4.18, Ki = -0.831, Kd = 4.26, Tf = 1.09

Name: Pitch_damper

Continuous-time PIDF controller in parallel form.

The same procedure also applies for the Throttle controller

TC0 = slTuner('Closed_loop_linearized_AC3DoF_model','Throttle controller');

% The tuning goal is the Pitch damper, so...
Req1 = TuningGoal.StepTracking('TH_cmd', 'Thrust', Gtr);
stepinfo(Req1.ReferenceModel);

% Let it tune, let it tune, I am one with the Plant and Controller!
TC  = systune(TC0, Req1);

Final: Soft = 0.0383, Hard = -Inf, Iterations = 20

Throttle_controller_values =

1 s

Kp + Ki * --- + Kd * --------

s Tf*s+1

with Kp = 0.00043, Ki = 1.23e-07, Kd = -0.000215, Tf = 0.498

Name: Throttle_controller

Continuous-time PIDF controller in parallel form.

So, despite a 20% overshooting, I am pretty satisfied with the rise time and the settling time, as you also can see from the simulation

This procedure makes sense to you? The 'target following' tuning was pretty straightforward!

[Sam Chak]

Hi @Leonardo Molino

I'm glad to hear that you learned from my systune example and followed the design procedure correctly. If you found the systune solution helpful, please consider clicking 'Accept' ✔️ on the answer. Additionally, you can show your appreciation by voting 👍 for my other answers as a token of support for my knowledge sharing. Your support is greatly appreciated!

Regarding the 20% overshoot issue, I have proposed two models below. Both models exhibit good settling times of approximately 4 seconds.

%% Target reference model 1
Gtr1 = tf(5, [1, 2*sqrt(5), 5])             % 2nd-order system

Gtr1 =
 
          5
  -----------------
  s^2 + 4.472 s + 5
 
Continuous-time transfer function.

%% Target reference model 2 (control effort is less demanding than Model 1)
wn   = 2.012;
Gtr2 = tf(wn^3, [1, 3*wn, 3*wn^2, wn^3])    % 3rd-order system

Gtr2 =
 
                8.145
  ---------------------------------
  s^3 + 6.036 s^2 + 12.14 s + 8.145
 
Continuous-time transfer function.

step(Gtr1), hold on
step(Gtr2), grid on
legend('Model 1', 'Model 2', 'fontsize', 16)

[Leonardo Molino]

Hi @Sam Chak, thank you again for you help!

Do you have any other advice/suggestion in order to tune also the outer loop?

It looks like this

The paper from which I have taken this concept design suggests to place KEP = KTP = 1 and KTI = KEI. I guess that the unkown must be determined imposing a constraint on the bandwidth, taking into account the one of the innerloops. Probably I need again looptune?

[Sam Chak]

Don't mention it @Leonardo Molino. The general rule is to ensure that the inner control loop responds faster than the outer control loop. Unfortunately, I'm not familiar with the system in the block diagram. It appears to be a cross-coupled system, where the dynamics of the states affect each other. I can't determine whether you need looptune() or not at the moment. However, I can provide an example for tuning a Multi-loop Control System using systune(). You can find it here:

https://www.mathworks.com/help/control/ug/tuning-multi-loop-control-systems.html

[Leonardo Molino]

@Sam Chak so, is it wrong to first design the innerloops as I did in the previous answer and then the outerloops? 😔

[Sam Chak]

Hi @Leonardo Molino

Based on the image, it is not clear which loop corresponds to the inner loop and which one corresponds to the outer loop. If you are utilizing an auto-tuning approach, it is possible to tune both loops simultaneously. However, translating specific performance requirements into tuning goals can often be challenging. Personally, I prefer to mathematically tune controllers for SISO systems, as I find it to be a more efficient approach.

For instance, the optimal LQR (Linear Quadratic Regulator) algorithm does not provide inputs where I can directly specify the desired overshoot and settling time for the closed-loop system. Conventional designers still have to painstakingly tune the Q and R weights.

By the way, were you able to achieve the desired performance requirements with the two proposed models?

[Leonardo Molino]

Hi @Sam Chak, sorry for the late answer but I was trying to figure out how to solve this problem.

So, yes, the Pitch damper and the Throttle controller do have the desired performances whitout any prefiltering, just via applying the "tuning goal" requirement. These two controllers is what I call (and the paper about TECS also) the "inner loops". Now, the TECS is a MIMO controller (it is what I call the outer loop), as you have said in the previous answer it recives multiple errors (the energy rate error and the energy rate distribution error) to generate a thrust (normalized) command and a pitch command. These outputs must be delivered to the inner loops, the one that I have tuned with your approach.

Now, before to design this MIMO controller, I have attempted the classic SISO approach, with an Altitude loop and a Speed loop as outer loops. I have implemented everything in this way

At first, I tune the innermost loops, "unplugging" the outer ones with the same code that I have posted previously. Then, I "plug" the outerloops and I tune them in the same way of the inner ones, but now I ask for a "slower" dynamics

Tune_inner = 0;

AL0     = slTuner('Closed_loop_linearized_AC3DoF_model','Altitude loop');
Gtr_out = tf(0.04, [1, 0.4, 0.04]); %Overshoot-free, 30secs dynamics
S4      = stepinfo(Gtr_out);

Req2 = TuningGoal.StepTracking('Alt_cmd', 'Alt_out', Gtr_out);
AL   = systune(AL0, Req2);

Final: Soft = 1.16, Hard = -Inf, Iterations = 23

T3 = getIOTransfer(AL, 'Alt_cmd', 'Alt_out');

figure();
step(Gtr_out, T3)
title('Step Response, Altitude loop')

% Update the block
writeBlockValue(AL)

SL0 = slTuner('Closed_loop_linearized_AC3DoF_model', 'Speed loop');

Req3 = TuningGoal.StepTracking('V_cmd', 'V_out', Gtr_out);
SL   = systune(SL0, Req3);

Final: Soft = 1.05, Hard = -Inf, Iterations = 20

T4 = getIOTransfer(SL, 'V_cmd', 'V_out');

figure();
step(Gtr_out, T4)
title('Step Response, Speed loop')

% Update the block
writeBlockValue(SL)

figure();
bode(T3)
grid on
title('Altitude loop')

omega_pitch = bandwidth(T3);
disp(['Altitude loop bandwidth is ', num2str(omega_pitch), ' rad/s'])

Altitude loop bandwidth is 0.1026 rad/s

figure();
bode(T4)
grid on
title('Speed loop')

omega_throttle = bandwidth(T4);
disp(['Speed loop bandwidth is ', num2str(omega_throttle), ' rad/s'])

Speed loop bandwidth is 0.10347 rad/s

Now, this approach works quite well, but I need the MIMO TECS design. Now, I think that this example Tuning of a Two-Loop Autopilot - MATLAB & Simulink - MathWorks Italia is showing something very close to what I need. I can exploit two MIMO controllers that uses energy rate error and the energy distribution error to generate the commanded pitch and the commanded thrust. I can tune them with the same second order dynamics, as required by the paper. What do you think about that?

19870017485.pdf (nasa.gov)

[Sam Chak]

Hi @Leonardo Molino

If the inner control loops have been successfully designed, then I believe you can proceed with constructing the TECS MIMO control architecture and apply the systune design methodology according to the desired tuning goals. Since this is a MIMO system, I cannot predict the outcome. It will be necessary to evaluate the performance of the systune approach to determine if it can achieve the desired outer control objectives.

[Leonardo Molino]

Hi @Sam Chak, thanks again for you help. Unfortunately, I cannot figure out how to implement this MIMO control logic, I will ask if I can keep the classic SISO autopilot autothrottle logic.

I have really appreciated your help with the classic tuning.

[Sam Chak]

Hi @Leonardo Molino, I apologize for not being able to provide further assistance regarding the MIMO case. Nevertheless, I encourage you to construct the TECS model in Simulink according to your understanding of the correct approach. Then, proceed to tune the controllers using the tools in MATLAB/Simulink you are comfortable with. Let's see if any error messages arise in the process.

[Leonardo Molino]

Hi @Sam Chak, don't worry it is not your fault!

The biggest problem with this TECS MIMO is that I don't know how to tell MatLab/Simulink to tune these gains

The Pitch controller and the Throttle controller were designed successfuly thanks to your help.

The issue right now is to find KTI = KEI, given KEP = KTP = 1. The idea is that both the energy rate error , which eventually generates the thrust command, and the distribution rate error, which eventually generates the pitch command, go to zero simultaneously with more or less the same dynamics. In general, we need

Moreover, the flightpath angle command and the norm. acceleration command are genereted in this way

in which I have other two gains (Kv and Kh) which must be the same.

My very first attempt was just to build the two errors and feeding them to two PI(s) controllers. The two controllers are then tuned with systune via giving a StepTracking command, but MatLab says that it was unable to enforce the CL-stability, so I guess this is the wrong approach...

[Leonardo Molino]

Hi @Sam Chak! After several days of struggling, I succeded in implementig and manually tuning the MIMO controller. The problem is that now the whole system is unstable at 10^4 seconds... I have some poles and zeros which slightly fall on the right side of the imaginary axis. Is there any work around to improve the stability? Maybe with a lead-phase compensator?

Thank you

[Sam Chak]

Hi @Leonardo Molino

MIMO systems often exhibit high-order characteristics due to the coupling effects from other states. If you use a low-order PID controller, the autotuner may struggle to find stabilizing gains. One approach you can try is using a slightly higher-order compensator. Ideally, if the compensator has the same order as the MIMO system, it should be capable of stabilizing it. However, in some cases, the transient response may appear peculiar, but this can typically be addressed using a pre-filter.

I also recommend creating a new question and describe what you’re trying to tune on the MIMO system. By doing so, 'fresh' unanswered questions with concise subject lines can often attract the attention, especially from experts who specialize in designing lead-lag compensators.

[Leonardo Molino]

Hi @Sam Chak, sorry if I bother you again after weeks, but I have some issues in tuning the SISO innerloops if I change the operating point. In particular, if I switch altitude/speed and I leave unchanged the tracking goal requirement (specified as transfer function), then I get this error

Final: Failed to enforce closed-loop stability (max Re(s) = 2.7e-05)

Failed to stabilize closed-loop system. Try different initial values or increase the number of random starts using "systuneOptions".

If I do a manual adjustment of the PIDs, I can find the gains correctly. I've manually tried several starting points and what I get as "optimal" is always this answer more or less

Altitude = 0m, Speed = 150 m/s

Altitude = 1000m, Speed = 150 m/s

Is there any script/tool that I can exploit to tell via code that I want always this response from theta_cmd to theta_out?

[Sam Chak]

Hi @Leonardo Molino

To be honest, I haven't fully explored the tuning capability of 'systune()'. I usually prefer to apply my own model-based design approach, using algebraic formulas to directly determine the controller gains for low-order stable SISO linear systems.

My experience with the auto-tuning tool for high-order systems is comparable to using 'fmincon()'. The success rate of finding the optimal values depends on the initial guess values and the feasibility of achieving the desired closed-loop response specified in 'TuningGoal.StepTracking()'.

Since your manually-tuned PID gains are capable of stabilizing the system, you can use these gains as the initial guess values for the auto-tuner. In the worst-case scenario, 'systune()' may indicate that the initial guess point appears to be a local minimum or solution because the gradient of your objective function is close to zero.

If you continue to encounter difficulties, I recommend deriving and obtaining the SISO model, whether in the form of a state-space representation or a transfer function. Working with accurate mathematical models tends to provide designers with a sense of control security, as the states are under their complete control and can be predicted deterministically.

[Leonardo Molino]

@Sam Chak Thank you! I will try

Answer 2

[Sam Chak]

Hi @Leonardo Molino

The transfer function of your aircraft appears to be a 5th-order system. Given its complexity, it's not surprising that a relatively low-order PD controller struggles to stabilize the system effectively. However, I have some preliminary design values for a 5th-order compensator and a 7th-order prefilter that aim to bring the pitch response to settle within 10 seconds, which is typical for a business jet aircraft. If I successfully tune the compensator using the 'looptune' command (which I'm not very familiar with), I will update this answer with the results.

Additionally, please note that the coefficients of your plant are displayed with only four significant digits by default, which is used for designing this compensator. If possible, you may consider saving the transfer function data from the Workspace into a '.mat' file and attaching it here for further analysis.

Edit: The code has been updated to incorporate the actual linearized system, and a new set of design values for both the Compensator and Prefilter have been provided. It is important to note that the old Compensator destabilizes the new Plant, and similarly, the new Compensator cannot stabilize the old Plant either.

%% Plant, Gp
s   = tf('s');
Gp0 = (-11.64*s^3 - 17.38*s^2 - 0.2759*s - 0.0001568)/(s^5 + 2.167*s^4 + 36.36*s^3 + 0.3578*s^2 + 0.3471*s + 0.0005279)

Gp0 =
 
           -11.64 s^3 - 17.38 s^2 - 0.2759 s - 0.0001568
  ---------------------------------------------------------------
  s^5 + 2.167 s^4 + 36.36 s^3 + 0.3578 s^2 + 0.3471 s + 0.0005279
 
Continuous-time transfer function.

load('delta_e_to_theta_lin.mat')
sys = ss(linsys1.A, linsys1.B, linsys1.C, linsys1.D);   % actual linearized system
Gp  = tf(sys)                                           % convert to transfer function

Gp =
 
           -11.64 s^3 - 17.38 s^2 - 0.2759 s - 0.0001568
  ---------------------------------------------------------------
  s^5 + 2.167 s^4 + 36.36 s^3 + 0.3578 s^2 + 0.3471 s + 0.0005279
 
Continuous-time transfer function.

step(Gp), grid on

%% Compensator, Gc
% old values
% cz  = [-1.07356000305645    + 5.93188324681864i;
%        -1.07356000305645    - 5.93188324681864i;
%        -0.00392143279578472 + 0.0976523949655871i;
%        -0.00392143279578472 - 0.0976523949655871i];
% cp  = [-46.8270041939968;
%         22.6535491689391 + 39.7360641466516i;
%         22.6535491689391 - 39.7360641466516i;
%        -1.47708667716432;
%        -0.0146074667170639];
% ck  = 8414.3275211565;

% new values
cz  = [-1.07379571191125    + 5.93147247179768i;
       -1.07379571191125    - 5.93147247179768i;
       -0.00392224270006151 + 0.0976581789345803i;
       -0.00392224270006151 - 0.0976581789345803i];
cp  = [-46.8206642188242    +                0i;
        22.650593013026     + 39.7307101862808i;
        22.650593013026     - 39.7307101862808i;
        -1.47703886330065   +                0i;
        -0.0146093282841061 +                0i];
ck  = 8410.96876788618;

Gc  = tf(zpk(cz, cp, ck))

Gc =
 
      8411 s^4 + 1.813e04 s^3 + 3.058e05 s^2 + 2570 s + 2919
  --------------------------------------------------------------
  s^5 + 3.011 s^4 - 27.16 s^3 + 9.789e04 s^2 + 1.461e05 s + 2113
 
Continuous-time transfer function.

%% Closed-loop transfer function, Gcl (with zeros)
% Gcl = feedback(Gc*Gp0, 1);    % unstable if the old Plant is used
Gcl = feedback(Gc*Gp, 1)

Gcl =
 
       -9.79e04 s^7 - 3.572e05 s^6 - 3.877e06 s^5 - 5.35e06 s^4 - 1.63e05 s^3 - 5.149e04 s^2 - 805.9 s - 0.4577
  ------------------------------------------------------------------------------------------------------------------
  s^10 + 5.179 s^9 + 15.72 s^8 + 32.86 s^7 + 50.72 s^6 + 59.44 s^5 + 52.8 s^4 + 34.7 s^3 + 16 s^2 + 4.665 s + 0.6579
 
Continuous-time transfer function.

%% Pre-filter, Gf
% old values
% fz  = [];
% fp  = [-1.07356000305645    + 5.93188324681866i;
%        -1.07356000305645    - 5.93188324681866i;
%        -1.47708635966585;
%        -0.00392143279578489 + 0.0976523949655872i;
%        -0.00392143279578489 - 0.0976523949655872i;
%        -0.0154505273816399;
%        -0.00059026071883032];
% fk  =  -6.71759286360342e-06;

% new values
fz  = [];
fp  = [-1.07379571191125     +   5.93147247179768i;
       -1.07379571191125     -   5.93147247179768i;
       -1.47703854570304     +                  0i;
       -0.00392224270006147  + 0.0976581789345804i;
       -0.00392224270006147  - 0.0976581789345804i;
       -0.0154526245974611   +                  0i;
       -0.000590080807089891 +                  0i];
fk  =  -6.72030185739573e-06;

Gf  = tf(zpk(fz, fp, fk))

Gf =
 
                                         -6.72e-06
  ---------------------------------------------------------------------------------------
  s^7 + 3.649 s^6 + 39.6 s^5 + 54.65 s^4 + 1.665 s^3 + 0.526 s^2 + 0.008232 s + 4.675e-06
 
Continuous-time transfer function.

%% Command Compensated System, Gcc (zeros are eliminated)
Gcc = minreal(series(Gf, Gcl))

Gcc =
 
                                                        0.6579
  ------------------------------------------------------------------------------------------------------------------
  s^10 + 5.179 s^9 + 15.72 s^8 + 32.86 s^7 + 50.72 s^6 + 59.44 s^5 + 52.8 s^4 + 34.7 s^3 + 16 s^2 + 4.665 s + 0.6579
 
Continuous-time transfer function.

S   = stepinfo(Gcc)

S = struct with fields:
         RiseTime: 3.9434
    TransientTime: 9.9997
     SettlingTime: 9.9997
      SettlingMin: 0.9016
      SettlingMax: 1.0149
        Overshoot: 1.4865
       Undershoot: 0
             Peak: 1.0149
         PeakTime: 13.5954

Ts  = round(S.SettlingTime);
step(Gcc, 3*Ts), grid on

[Leonardo Molino]

Hello @Sam Chak, thank you for you answer.

First of all, I would like to ask you where I can read more about the "compensator-prefilter" technique that you have mentioned in your answer.

In addition to this, I will attach to this message the linear analysis .mat file from the Simulink Model Linearizer app.

Thank you

[Sam Chak]

Hi @Leonardo Molino

I have revised the code in my answer. Please review it.

The design concept of the Compensator–Prefilter is elaborated upon in the chapter titled "The Design of Feedback Control Systems" in the book "Modern Control Systems" by Dorf and Bishop. Ogata's book "Modern Control Engineering" also addresses this design technique in the chapter titled "PID Controllers and Modified PID Controllers", specifically in Examples A-8-8 to A-8-10. However, Ogata refers to it as the "Input Filter".

I frequently draw an analogy between the Compensator and Shampoo, and the Prefilter and Conditioner. The Compensator serves a primary role akin to cleansing the system, facilitating the removal of unwanted or unstable poles. Consequently, this action often introduces zeros in the closed-loop transfer function, some of which may yield undesirable effects. Conversely, akin to a Conditioner applied after shampooing, the Prefilter is employed to enhance the system's transient response by mitigating the impact of negative zeros through the manipulation of the filter's poles.

Example:

Let's consider a Double Integrator system, .

To achieve a critically-damped step response (without overshoot), the control signal u can be designed using a Proportional-Derivative (PD) controller . This configuration results in the desired mass-damper-spring-like structure .

In MATLAB and Simulink, it's common for students to directly employ a PID controller, inadvertently introducing a zero into the closed-loop transfer function. However, the effect of this zero can be readily mitigated by implementing a Prefilter at the command input port.

s   = tf('s');

%% Plant, Gp
Gp  = 1/s^2

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

%% PD controller, Gc
Kp  = 1;
Kd  = 2;
Gc  = pid(Kp, 0, Kd)

Gc =
 
             
  Kp + Kd * s
             

  with Kp = 1, Kd = 2
 
Continuous-time PD controller in parallel form.

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

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

%% Prefilter (a.k.a Input Filter), Gf
Gf  = 1/(2*s + 1)

Gf =
 
     1
  -------
  2 s + 1
 
Continuous-time transfer function.

%% Command compensated system (as coined by Ogata), Gcc
Gcc = zpk(series(Gf, Gcl))

Gcc =
 
      (s+0.5)
  ---------------
  (s+1)^2 (s+0.5)
 
Continuous-time zero/pole/gain model.

Gcc = tf(minreal(Gcc))

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

step(Gcl), hold on
step(Gcc), grid on
legend('Gcl (without prefilter)', 'Gcc (with prefilter)', 'location', 'east')

[Leonardo Molino]

Dear @Sam Chak, thank you for your answer, it is very clear and well exposed.

I have read about the "Prefilter-Compensator" technique in the Ogata's book. If you don't mind, I have some other questions that I would like to ask you.

If I want to generate the prefilter TF for my own, shoud I consider the formula of this example

from the Ogata's book (Example A-8-8, as you suggested)?

If the answer to the latter question is yes, the gains are just the ones that came from the PID mask block after tuning on Simulink?

If I want to indroduce also the actuator dynamics, can I consider the elevator TF within the block?

Answer 3

[Sam Chak]

Hi @Leonardo Molino

If you can tolerate a slightly slower response, the unconventional PID controller can still achieve a relatively satisfactory performance. I must admit that I haven't quite figured out how to effectively use the 'looptune()' command, as it tends to generate frustrating messages regarding the connection of port names. The lack of sufficient examples in the documentation makes it challenging to learn 😤. Therefore, I opted to use the more user-friendly 'systune()' command instead 😄.

format long g

%% Plant
load('delta_e_to_theta_lin.mat')
sys = ss(linsys1.A, linsys1.B, linsys1.C, linsys1.D);   % actual linearized system
Gp  = tf(sys)

Gp =
 
           -11.64 s^3 - 17.38 s^2 - 0.2759 s - 0.0001568
  ---------------------------------------------------------------
  s^5 + 2.167 s^4 + 36.36 s^3 + 0.3578 s^2 + 0.3471 s + 0.0005279
 
Continuous-time transfer function.

S1  = stepinfo(Gp);

%% Compensator Gc with tunable PID parameters
Gc0 = tunablePID('Gc', 'PID');

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

%% Target Response
Gtr = tf(0.04, [1, 0.4, 0.04])  % <-- the settling time of this system is about 30 sec

Gtr =
 
         0.04
  ------------------
  s^2 + 0.4 s + 0.04
 
Continuous-time transfer function.

S2  = stepinfo(Gtr);

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

%% Let it tune, let it tune, I am one with the Plant and Controller!
CL  = systune(CL0, Req);

Final: Soft = 1.7, Hard = -Inf, Iterations = 426

%% Display tuned PID Controller
Kp  = CL.Blocks.Gc.Kp.Value;
Ki  = CL.Blocks.Gc.Ki.Value;
Kd  = CL.Blocks.Gc.Kd.Value;
Tf  = CL.Blocks.Gc.Tf.Value;
Gc  = pid(Kp, Ki, Kd, Tf)

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

  with Kp = 15.9, Ki = -0.172, Kd = -1.5e+03, Tf = 93.2
 
Continuous-time PIDF controller in parallel form.

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

Gcl =
 
          2.197 s^5 + 3.299 s^4 + 0.1029 s^3 + 0.03262 s^2 + 0.0005104 s + 2.898e-07
  -------------------------------------------------------------------------------------------
  s^7 + 2.178 s^6 + 38.58 s^5 + 4.047 s^4 + 0.4539 s^3 + 0.03688 s^2 + 0.000516 s + 2.898e-07
 
Continuous-time transfer function.

ssO = dcgain(Gcl)   % Steady-state output response

ssO = 
     1

S3  = stepinfo(Gcl);

%% Plot results and comparison
stepInfoTable = struct2table([S1, S2, S3]);
stepInfoTable = removevars(stepInfoTable, {'SettlingMin', 'SettlingMax', 'Undershoot', 'Peak', 'PeakTime'});
stepInfoTable.Properties.RowNames = {'Plant', 'Target System', 'Closed-loop TF'};
stepInfoTable

stepInfoTable = 3×4 table
                          RiseTime         TransientTime        SettlingTime          Overshoot    
                      ________________    ________________    ________________    _________________

    Plant             0.42708291413902    1244.01348852939                 NaN     1677.10726435862
    Target System      16.791687222308    29.1705588453553    29.1705588453553                    0
    Closed-loop TF     22.488320591896    32.9827752129389    32.9827752129389    0.595147170485277

stepplot(Gtr,  90), hold on
stepplot(Gcl), grid on
legend('Target System', 'Closed-loop TF', 'location', 'East')