Clear Filters
Clear Filters

Unable to tune a pitch-damper autopilot gains with looptune

19 views (last 30 days)
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.
  2 Comments
Sam Chak
Sam Chak on 12 Mar 2024
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
Leonardo Molino on 12 Mar 2024
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

Sign in to comment.

Accepted Answer

Sam Chak
Sam Chak on 13 Mar 2024
Edited: Sam Chak on 14 Mar 2024
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.
  16 Comments
Sam Chak
Sam Chak on 4 Apr 2024
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.

Sign in to comment.

More Answers (2)

Sam Chak
Sam Chak on 12 Mar 2024
Edited: Sam Chak on 13 Mar 2024
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
  3 Comments
Sam Chak
Sam Chak on 13 Mar 2024
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
Leonardo Molino on 13 Mar 2024
Edited: Leonardo Molino on 13 Mar 2024
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?

Sign in to comment.


Sam Chak
Sam Chak on 13 Mar 2024
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')

Products


Release

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!