PRÁCTICA 3 (2017)
Apunte InglésUniversidad  Universidad Politécnica de Cataluña (UPC) 
Grado  Ingeniería de Aeronavegación  3º curso 
Asignatura  Control y Guiaje 
Año del apunte  2017 
Páginas  18 
Fecha de subida  25/06/2017 
Descargas  8 
Subido por  areig 
Vista previa del texto
UNIVERSITAT POLITÈCNICA DE
CATALUNYA
Enginyeria d’Aeronavegació
Matlab Tutorial #4: StateSpace
Autopilot Design
Control & Guidance
June 11, 2017
Spring Term
Albertí, Joan
Nogués, Sergi
Reig, Anna
6GX3
CONTENTS
Matlab Tutorial #4
Contents
1 Introduction
2
2 Suboptimal Control with Bass & Gura Method
3
3 Optimal Control with LQR method
7
4 Beyond the applied methods
10
5 Annex
14
List of Figures
1
2
3
4
5
6
7
Nonscaled Suboptimal ∆h response to ∆href = 100f t. . . . . 4
Schematic of the statefeedback system with scaling factor . . 5
Suboptimal ∆h response to ∆href = 100f t. . . . . . . . . . . . 6
Nonscaled and Scaled Optimal ∆h response to ∆href = 100f t.
(left and right respectively) . . . . . . . . . . . . . . . . . . . 9
Aircraft modes . . . . . . . . . . . . . . . . . . . . . . . . . . 11
Natural frequency vs damping ratio . . . . . . . . . . . . . . . 11
Updated response using new eigenvalues. . . . . . . . . . . . . 12
1
Introduction
1
Matlab Tutorial #4
Introduction
Based on the study of a STOL (Short TakeOff and Landing) transport aircraft, this handson activity is intended to use state feedback to design an
autopilot able to maintaining a constant altitude. In this way, an utterly
control of the aircraft will be acquired.
However, for the sake of simplicity, some assumptions have been assumed
such as a fixed forward speed as well as the exclusion of the actuator dynamics. Hence, the state equation can be set out.
For our case study, the state equation is represented by the shortperiod
approximation as well as the kinematic equation representing the change in
vertical height presented as follows:
∆h˙ = uo (∆θ − ∆α)
(1)
Therefore, once the prior assumptions are applied and adding the vertical velocity equation to the shortperiod equations, we obtain the following fourthorder system:
∆α(t)
˙
−1.397
1
0 0
∆α(t)
−0.124
∆q(t))
−5.47 −3.27 0 0 ∆q(t) −13.2
˙
·
+
· ∆δe (t)
˙ = 0
∆θ(t)
1
0 0 ∆θ(t) 0
˙
−400
0
400 0
∆h(t)
0
∆h(t)
(2)
Where the units are American Units (US), thus feet, radians and radians per
second are used. Parameters are the following:
• α: angle of attack (radians).
• q: rotation speed (radians/second).
• θ: pitch angle (radians).
• h: altitude (feet).
2
Suboptimal Control with Bass & Gura Method
2
Matlab Tutorial #4
Suboptimal Control with Bass & Gura Method
First, the BassGura method will be implemented in order to get the required
closedloop eigenvalues shown below.
λ1,2 = −3 ± 2.2i and λ3,4 = −1 ± 1.1i
Throughout the implementation of the BassGura method with Matlab using
the characteristic polynomial of the state space matrix in open loop, the matrix of controllability associated to it as well as the characteristic polynomial
from the eigenvalues, the algorithm leads to obtain the needed gain values.
k1 = 1.2801
k2 = −0.2645
k3 = −2.3442
k4 = −0.0043
At first glance, it could be said that the obtained gain values are pretty
reasonable since they are not really large. At a later stage, this fact will
avoid taking big control action for the given state matrix. In other words,
closedloop poles will not be very further away from open loop poles. Once
know the gain values, the corresponding control law can be computed as a
function of the rate of the already defined variables.
(3)
u = −k · xT
∆δe (t) = 1.2801 · ∆α(t) − 0.2645 · ∆q(t) − 2.3442 · ∆θ(t) − 0.0043 · ∆h(t) (4)
Furthermore, reasonable values corresponding to a cruise phase have to be
selected and tested for ∆θ(t), ∆q(t), ∆α(t) and ∆h(t) = 100f t leading to a
certain elevator deflection ∆δe (t) value. Therefore, the selected inputs and
elevator deflection are the following:
∆θ(t) = 1°
∆q(t) = 1
deg
s
∆α(t) = 2°
∆h(t) = 100f t
∆δe (t) = 24.69°
Undoubtedly, the obtained value for the elevator deflection is not a coherent value because normally elevators cannot exceed ±15°. Nevertheless, the
sign is correct because the aircraft is to experiment a positive change of altitude, thus its nose must point upwards (i.e deflection must be positive so
downwards to create the proper moment). This fact allows us to see that
3
Suboptimal Control with Bass & Gura Method
Matlab Tutorial #4
even having chosen relative small input values, the elevator deflection is outside the typical range. Hence, neither the computed gain values are not the
appropriated ones and consequently nor the used method for the studied case.
In addition, we would like to get the ∆h closedloop response to a reference
value of ∆href = 100f t. In order to do that, we will have to apply a feedback
using the previously obtained gain values as it follows:
A∗ = A − B · k T
(5)
Figure 1: Nonscaled Suboptimal ∆h response to ∆href = 100f t.
As we can see in Figure 1, the response of the aircraft is unacceptable because
the steady state error is too large having a magnitude of 104 . Indeed, this
is normal because we have not already scaled the reference input so that
the output equals the reference in steady state. Thus, the basic schematic
corresponding to a statefeedback system with scaling factor is shown in
Figure 2.
4
Suboptimal Control with Bass & Gura Method
Matlab Tutorial #4
Figure 2: Schematic of the statefeedback system with scaling factor
When scaling the overall output, the steadystate error should be reduced.
Hence, once the prior mentioned precompensator scale factor is found by a
given function (rscale.m) the response is depicted in Figure 3.
Now, the response is completely reasonable as it reaches the wanted href
at a settling time value equal to 4.52s. Besides the large steady state error
has been removed. However, this response is not perfectly perfect as we
can appreciate, it has some acceptable overshoot (5.70%) and consequently a
peak not really far from the href . At a later stage, it will be eliminate using
a more sophisticated method and accurate limits of our variables. The most
interesting values of the response are detailed below.
• Overshoot=5.70%
• Peak =105.7 ft
• Peak Time=3.35 s
• Settling Time=4.52 s
5
Suboptimal Control with Bass & Gura Method
Matlab Tutorial #4
Figure 3: Suboptimal ∆h response to ∆href = 100f t.
6
Optimal Control with LQR method
3
Matlab Tutorial #4
Optimal Control with LQR method
So far, no constraints have been imposed to limit angles. However, if we
want to optimise the response of the aircraft, some restrictions need to be
imposed. So, to obtain an optimal StateSpace control system design, the
Linear Quadratic Regular (LQR) has to be applied and thus, the next index
performance has to be minimised.
tf
(6)
(xt Qx + ut Ru)dt → minimised
J=
to
Where,
• Q is the weighting matrix relating the different x states
• R is the matrix relating the vector u which is a control vector able to
drive the state x0 (t0 ) to the desired final value xf (tf ).
That is for our case:
∞
J=
0
∆α
∆αmax
2
+
∆θ
∆θmax
2
+
∆h
∆hmax
2
+
∆δ
∆δmax
2
dt
(7)
Where it is imposed that,
∆θmax = 10°
∆αmax = 8°
∆hmax = 120f t
∆δemax = 15°
Once the parameters have been imposed, we have to define and fulfil matrix Q
and R belonging to the performance index so to implement the minimisation
throughout Matlab. Therefore, Q and R matrixes can be defined as below
and applying the LQR, the associated gain values can be obtained.
1
0
0
0
∆α2max
0
0
0
0
1
Q=
R=
(8)
1
0 ∆θ2
0
∆δe 2max
0
max
1
0
0
0
∆h2
max
k1 = 0.4611
k2 = −0.3729
k3 = −2.5966
7
k4 = −0.0022
Optimal Control with LQR method
Matlab Tutorial #4
Because of the imposed constraints and the LinearQuadratic Regulator method
used, as seen above, we get slightly different gain values in comparison to the
ones obtained in suboptimal control (i.e in terms of numerical values). However, from one method to another, their order of magnitude is maintained.
Notice that we are asked to use the same values we have used in the first
section to calculate the optimal control law and hence, the order of magnitude
of the elevator deflection. These gain values yield to the following optimal
control law:
∆δopt = 0.4611 · ∆α(t) − 0.3729 · ∆q(t) − 2.5966 · ∆θ(t) − 0.0022 · ∆h(t) (9)
Therefore, testing the aforementioned values of ∆θ(t), ∆q(t), ∆α(t) and
∆h(t), we get an optimal elevator deflection ∆δopt = 14.55°. The optimal
elevator deflection value is considerably smaller than before taking into account that the angle has decreased by approximately 10°being now between
the typical ranges (i.e ± 15°) The positive sign is the appropriate again since
we want the aircraft’s nose pointing upwards requiring a downwards elevator’s deflection.
Then, we are asked to get the ∆h response to a reference value of ∆href =
100f t. First, the noncompensated response will be depicted and second, the
rescaled response will be represented too.
As it can be appreciated in the next figures, for the step resulting from the
application of the LQR method without scaling, a large steady state error is
obtained again compared to the step with scaling. Nevertheless, when comparing with the previous method, although the steady state error has the
same order of magnitude (i.e 104 ), it is greater when using LQR whose final
value is obtained between 40000 and 50000, while, for the suboptimal case,
it is obtained between 20000 and 25000.
8
Optimal Control with LQR method
Matlab Tutorial #4
Figure 4: Nonscaled and Scaled Optimal ∆h response to ∆href = 100f t.
(left and right respectively)
Fortunately, when both optimisation and rescaling are applied as seen in
Figure 4 (right), the response is smoother with an overshoot practically negligible (i.e 0.075% against 5.70% for the suboptimal method). However, it is
slower having a greater settling time (i.e 7.5s against 4.52s for the suboptimal
method). These adjectives qualifying the response are more important than
it can be thought. Having this response will mean that the passengers of the
flight will barely notice changes on their body when stabilising the aircraft
as the aircraft has a mild and progressive behaviour.
Finally, everything explained so far has been outlined in the following table.
Bass & Gura Method
(with scaling factor)
Steady State Error
0
Overshoot
5.70%
Settling Time
4.52s
Peak
105ft at 3.35s
LQR Method
(with scaling factor)
0
0.075%
7.5s
Negligible
Table 1: Comparison between Bass & Gura Method and LQR Method with
scaling factor applied
9
Beyond the applied methods
4
Matlab Tutorial #4
Beyond the applied methods
So far, results for Bass & Gura and Linear Quadratic Regulator Method have
been obtained, analysed and then compared. Now we are asked to check the
influence of the desired eigenvalues (BassGura method) on the response of
the system.
To know what are the desired eigenvalues, first we must deduce some information from the initial ones. As seen below, is possible to extract the
damping factor and the natural frequencies associated to each pair of initial eigenvalues. In this way, these characteristics terms will tell us how to
improve them.
λ1,2 = −3 ± 2.2i
→
λ2 + 6λ + 13.8
λ3,4 = −1 ± 1.1i
→
λ2 + 2λ + 2.2
⇒ ξ = 0.81 and ωn = 3.72 rad/s
⇒ ξ = 0.67 and ωn = 1.5 rad/s
Therefore, since we have two couples of complex poles, the overall response
will be dominated by the slowest and the most damped component (i.e ωn =
1.5 rad/s and ξ = 0.81).
Before finding the values that improve the response of the system, it should
be considered the following statements:
• Larger gains lead to bigger control action for a given state vector, so
we are looking for closedloop poles to be further away from open loop
poles. Recall:
k = (V W T )−1 (α − α
¯)
(10)
But control input can become very large if gains are too large, and
exceeds servo actuator capability to respond (physical limits): actuator
then saturated system may not perform as expected (more powerful
actuator needed).
• It should be avoided the excitation of unmodeled structural modes by
increasing too much the closedloop frequency response.
10
Beyond the applied methods
Matlab Tutorial #4
Figure 5: Aircraft modes
• The following figure shows how eigenvalues should be placed:
Figure 6: Natural frequency vs damping ratio
Considering these statements, we can corroborate that the previous eigenvalues lead to an acceptable response, however, since the natural frequency
is below 2, the response is a bit slow. For satisfactory eigenvalues, it is
needed a natural frequency between 2 rad/s and 5 rad/s, if possible between 3 rad/s and 4 rad/s to improve their placement and subsequently
its behaviour in front of a step. Concerning the damping factor, it should be
close to 1 so to have the minimal oscillations. Consequently, the proposed
11
Beyond the applied methods
Matlab Tutorial #4
eigenvalues should be a little bit more at the left side of the imaginary plane
(i.e to increase the natural frequency) and a little bit closer to the real axis
(i.e to increase the damping factor). According to all said before, varying
the location of the poles as exposed, the step response obtained should be
better.
Once known the influence of these desired eigenvalues, let’s calculate them.
λ1,2 = −3 ± 1.5i
→
λ2 + 6λ + 11.25
⇒ ξ = 0.9 and ωn = 3.35 rad/s
λ3,4 = −4 ± 1.1i
→
λ2 + 8λ + 17.2
⇒ ξ = 0.96 and ωn = 4.14 rad/s
Now, we have achieved a better behaviour within the satisfactory zone by
means of the relocation of the poles. Indeed, we have reached an improved
response in which the climb is smoother and slightly faster. A depiction of
the result is shown next.
Figure 7: Updated response using new eigenvalues.
12
Beyond the applied methods
Matlab Tutorial #4
To sum up all said before and to highlight the advantages and disadvantages
of each method, the following table has been set up.
Bass & Gura
Feedback/gain
Elevator
deflection
Step
response
Step response
(after scaling)
Design
complexity
Reasonable
and not really large.
Coherent,
within limits
Unacceptable.
Too large steady
state error
(magnitude of 10000)
Reasonable:
Overshoot=5.70%
Peak=105.7 ft
Peak Time=3.35 s
Settling Time=4.52 s
Not much
LQR
Slightly different gain values
than BassGura’s Method
but order of magnitude
is maintained.
Mildly larger than BassGura’s.
An increase of approximately 3°
Unnacceptable.
Steady state error has same
order of magnitude
but twice bigger.
Smoother with an overshoot
practically negligible
and slower having
a greater settling time.
More complex than BassGura
Table 2: Advantages and disadvantages of both methods used
Depending on what the system is going to be used for, we will be more
interested in BassGura or in LQR but if we are looking for optimal control,
LQR is better as its step response is more adequate for an airplane’s flight
level change.
13
Annex
5
Matlab Tutorial #4
Annex
Main code of the practice
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
%% PART 1.1
%% Matrix A and B
A=[−1.397 1 0 0; −5.47 −3.27 0 0; 0 1 0 0; −400 0 400 0];
B=[−0.124; −13.2; 0; 0];
C=[0 0 0 1];
%% Characteristic polynomial of A in openloop
charac_pol_A=poly(A);
eig_A=eig(A);
a=[charac_pol_A(2) charac_pol_A(3) charac_pol_A(4) charac_pol_A
(5)];
%% Controllability matrix
controllability=[B A*B A^2*B A^3*B];
%% Characteristic polynomial with desired coefficients
roots=[−3+2.2i −3−2.2i −1+1.1i −1−1.1i];
desired_coef=poly(roots);
a_desired=[desired_coef(2) desired_coef(3) desired_coef(4)
desired_coef(5)];
%% State space: Bass&Gura
W=[1 charac_pol_A(2) charac_pol_A(3) charac_pol_A(4)
0 1 charac_pol_A(2) charac_pol_A(3)
0 0 1 charac_pol_A(2)
0 0 0 1];
Z=((controllability*W)')^(−1);
a_dif=(a_desired−a);
k=Z*a_dif'; %%%%%%%%%%%%%%%
%% PART 1.2
%% Control law
%u=−k^t*x
14
Annex
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
Matlab Tutorial #4
alpha=2*pi/180;
q=1*pi/180; %angular velocity
theta=1*pi/180; %pitch
x=[alpha;q;theta;100];
control_law=−k'*x;
delta_e=control_law*180/pi; %%%%%%%%%%
%magnitud: 11.53 degrees
%signe: positive
%% PART 1.3
%% Plot delta_h
h_ref=100;
sys = ss(A,B,h_ref*C,0);
A_feedback=A−B.*k';
sys_feedback=ss(A_feedback,B,h_ref*C,0);
%opt.StepAmplitude=100;
figure(1); %%%%%%%%%%%%%%%%%
step(sys_feedback);
%% PART 1.4
%% Plot delta_h rescale
Nbar=rscale(A,B,C,0,k);
sys_feedback_rscale=ss(A_feedback,B*Nbar,h_ref*C,0);
figure(2); %%%%%%%%%%%%%%%%%%%%
step(sys_feedback_rscale);
%% PART 2.5
%% LQR−Optimal control law−same values as PART 1.2
alpha_max=8*pi/180;
q_max=0;
theta_max=10*pi/180;
delta_h_max=120;
delta_e_max=15*pi/180;
Q=[1/(alpha_max)^2 0 0 0; 0 0 0 0; 0 0 1/(theta_max)^2 0; 0 0 0
1/(delta_h_max)^2];
67 R=1/(delta_e_max)^2;
68 [optimal_k]=lqr(A,B,Q,R);
69
15
Annex
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
Matlab Tutorial #4
%% Control law
%u=−k^t*x
alpha_opt=2*pi/180;
q_opt=pi/180; %angular velocity
theta_opt=pi/180; %pitch
x_opt=[alpha_opt;q_opt;theta_opt;100];
control_law_opt=−optimal_k*x_opt;
delta_e_opt=control_law_opt*180/pi; %%%%%%%%%%
%magnitud: 14.5472 degrees
%signe: positive
%% PART 2.6
%% Plot delta_h
sys = ss(A,B,h_ref*C,0);
A_feedback_opt=A−B.*optimal_k;
sys_feedback_opt=ss(A_feedback_opt,B,h_ref*C,0);
figure(3); %%%%%%%%%%%%%%%%%%%%%%
step(sys_feedback_opt);
%% Plot delta_h rescale
Nbar_opt=rscale(A,B,C,0,optimal_k');
sys_feedback_opt_rscale=ss(A_feedback_opt,B*Nbar_opt,h_ref*C,0)
;
92 figure(4); %%%%%%%%%%%%%%%%%%%%
93 step(sys_feedback_opt_rscale);
Rescale function
1 function[Nbar]=rscale(a,b,c,d,k)
2
% Given the single−input linear system:
3
%
.
4
%
x = Ax + Bu
5
%
y = Cx + Du
6
% and the feedback matrix K,
7
%
8
% the function rscale(sys,K) or rscale(A,B,C,D,K)
16
Annex
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42 end
Matlab Tutorial #4
% finds the scale factor N which will
% eliminate the steady−state error to a step reference
% for a continuous−time, single−input system
% with full−state feedback using the schematic below:
%
%
/−−−−−−−−−\
%
R
+
u  .

%
−−−> N −−−>() −−−−> X=Ax+Bu −−> y=Cx −−−> y
%
−
\−−−−−−−−−/
%


%
<−−−− K <−−−−
%
%8/21/96 Yanjie Sun of the University of Michigan
%
under the supervision of Prof. D. Tilbury
%6/12/98 John Yook, Dawn Tilbury revised
%11/22/16 Revised by Adeline de Montlaur
narginchk(2,5); % Validate number of input arguments
% −−− Determine which syntax is being used −−−
% nagrin: Number of function input arguments
if (nargin==2) % System form
[A,B,C,D] = ssdata(a); % Quick access to state−space
data
K=b;
elseif (nargin==5) % A,B,C,D matrices
A=a; B=b; C=c; D=d; K=k;
else error('Input must be of the form (sys,K) or (A,B,
C,D,K)')
end;
% compute Nbar
s = size(A,1);
Z = [zeros([1,s]) 1];
N = inv([A,B;C,D])*Z';
Nx = N(1:s);
Nu = N(1+s);
Nbar=Nu + K'*Nx;
17
...