Unable to tune a pitch-damper autopilot gains with looptune
24 views (last 30 days)
Show older comments
Leonardo Molino
on 12 Mar 2024
Commented: Leonardo Molino
on 5 Apr 2024
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
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.
Accepted Answer
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
%% 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)
%% Closed-loop system, Gcl
Gcl = zpk(feedback(Gc*Gp, 1))
%% Prefilter, Gf
Gf = tf(1, [Ti*Td, Ti, 1]) % Ogata's formula (see A-8-8)
%% Command compensated system, Gcc
Gcc = series(Gf, Gcl)
Gcc = minreal(Gcc)
%% 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
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
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.
More Answers (2)
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)
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
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))
%% Closed-loop transfer function, Gcl (with zeros)
% Gcl = feedback(Gc*Gp0, 1); % unstable if the old Plant is used
Gcl = feedback(Gc*Gp, 1)
%% 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))
%% Command Compensated System, Gcc (zeros are eliminated)
Gcc = minreal(series(Gf, Gcl))
S = stepinfo(Gcc)
Ts = round(S.SettlingTime);
step(Gcc, 3*Ts), grid on
3 Comments
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
%% PD controller, Gc
Kp = 1;
Kd = 2;
Gc = pid(Kp, 0, Kd)
%% Closed-loop system, Gcl
Gcl = feedback(Gc*Gp, 1)
%% Prefilter (a.k.a Input Filter), Gf
Gf = 1/(2*s + 1)
%% Command compensated system (as coined by Ogata), Gcc
Gcc = zpk(series(Gf, Gcl))
Gcc = tf(minreal(Gcc))
step(Gcl), hold on
step(Gcc), grid on
legend('Gcl (without prefilter)', 'Gcc (with prefilter)', 'location', 'east')
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)
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
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);
%% 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)
%% Closed-loop system
Gcl = feedback(Gc*Gp, 1)
ssO = dcgain(Gcl) % Steady-state output response
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
stepplot(Gtr, 90), hold on
stepplot(Gcl), grid on
legend('Target System', 'Closed-loop TF', 'location', 'East')
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!