# 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

### 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

### Categories

### Community Treasure Hunt

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

Start Hunting!