PRÁCTICA 3 (2017)Apunte Inglés
Vista previa del texto
UNIVERSITAT POLITÈCNICA DE
Matlab Tutorial #4: State-Space
Control & Guidance
June 11, 2017
Matlab Tutorial #4
2 Suboptimal Control with Bass & Gura Method
3 Optimal Control with LQR method
4 Beyond the applied methods
List of Figures
Non-scaled Suboptimal ∆h response to ∆href = 100f t. . . . . 4
Schematic of the state-feedback system with scaling factor . . 5
Suboptimal ∆h response to ∆href = 100f t. . . . . . . . . . . . 6
Non-scaled 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 Take-Off and Landing) transport aircraft, this hands-on 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 short-period 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 short-period 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 Bass-Gura method will be implemented in order to get the required closed-loop eigenvalues shown below.
λ1,2 = −3 ± 2.2i and λ3,4 = −1 ± 1.1i Throughout the implementation of the Bass-Gura 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, closed-loop 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 closed-loop 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: Non-scaled 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 state-feedback system with scaling factor is shown in Figure 2.
4 Suboptimal Control with Bass & Gura Method Matlab Tutorial #4 Figure 2: Schematic of the state-feedback system with scaling factor When scaling the overall output, the steady-state 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 State-Space 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 Linear-Quadratic 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 non-compensated response will be depicted and second, the re-scaled 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: Non-scaled and Scaled Optimal ∆h response to ∆href = 100f t.
(left and right respectively) Fortunately, when both optimisation and re-scaling 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 (Bass-Gura 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 closed-loop 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 closed-loop 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 Bass-Gura’s Method but order of magnitude is maintained.
Mildly larger than Bass-Gura’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 Bass-Gura 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 Bass-Gura 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 ...