< Day Day Up > |
7.4 Summary
In this chapter, I presented two MIMO control design examples. The first (Example 7.2) was an aircraft control problem in which the plant had two inputs (aircraft angle of attack α and engine thrust T) and three outputs (the measured aircraft speed V, the measured flight path angle γ, and the measured distance of the aircraft above or below the glideslope herr). This plant has significant cross-coupling between the inputs and outputs, so breaking the model into separate SISO systems would not be a satisfactory approach. I used the pole placement design method to develop an observer-controller for this system that provides satisfactory performance.
Example 7.2: Aircraft glideslope control.
Consider a large passenger jet aircraft on approach to landing. By sensing radio signals transmitted from the ground-based instrument landing system (ILS), the aircraft is able to determine how far it is above or below the desired flight path, which is called the glideslope. The design task is to develop an autopilot control system to fly the aircraft in the vertical plane along the glideslope while maintaining its speed at a desired constant setting. I will ignore the side-to-side motion and roll orientation of the aircraft and assume those flight aspects are handled adequately by other control systems. I will also ignore the effects of wind and the actions required to initially position the aircraft on the glideslope and to land the aircraft when it reaches the runway.
The system configuration is shown in Figure 7.1, where x is the forward direction in aircraft body coordinates. The engine thrust T (in units of newtons) acts in the x direction. The angle of attack α (radians) is the angle between the aircraft velocity vector V (m/s) and the body x-axis. The flight path angle γ (radians) is the angle between the horizontal and the velocity vector V. As drawn in Figure 7.1, γ has a negative value. The drag force D (newtons) acts in the negative direction of the velocity vector. The lift force L (newtons) acts in an upward direction perpendicular to the velocity vector.
Figure 7.1: Aircraft and glideslope.
The glideslope is a fixed line in space radiating upward from the runway at an angle of γgs (radians) relative to the horizontal. The vertical distance from the aircraft to the glideslope is the height error, herr (meters). The goal of the autopilot is to drive herr to zero and simultaneously maintain the commanded velocity Vc (m/s).
This plant has three outputs: the measured aircraft speed V, the measured flight path angle γ, and the measured distance of the aircraft above or below the glideslope herr. The plant inputs are the throttle setting and the elevator command. The elevator is the moveable surface along the back edge of the horizontal tail that controls the aircraft angle of attack. These two plant inputs enable control of the aircraft speed and vertical motion.
However, there is significant cross-coupling between the two inputs and the three outputs. An increase in the throttle setting causes the speed to increase but also causes the aircraft to climb in altitude. An upward deflection of the elevator causes the aircraft to climb but also results in a decrease in speed. The use of multiple SISO systems in the design process would not account for this cross-coupling, so a MIMO design approach is preferable.
The response time of the engines to throttle commands and the airframe to elevator commands is fast compared with changes in the aircraft speed and altitude rate. To keep the model simple, those effects will be ignored. Instead, I will change the assumed plant inputs to be the engine thrust T and the angle of attack a.
This simplified system is described by the nonlinear equations shown in Eq. 7.1 [1]. The first row of Eq. 7.1 is Newton's law applied along the direction of the velocity vector. The second row uses the forces perpendicular to the velocity vector to determine the rate of change of the flight path angle. The third row describes the rate of change of the altitude error on the basis of the aircraft velocity, flight path angle, and glideslope angle.
| (7.1) |
The lift and drag force terms (L(α, V) and D(α, V), in newtons) of Eq. 7.1 are shown in Eq. 7.2, in which ρ is the air density and S is the aircraft reference area.
| (7.2) |
The parameters appearing in Eqs. 7.1 and 7.2 are defined in Table 7.1. These values represent a nominal flight condition during approach for a landing.
|
|
|
|
|---|---|---|---|
| Aircraft mass | 190,000 | kg |
| Gravitational acceleration | 9.81 | m/s2 |
| Glideslope angle | 3 0.052359 | degrees radians |
Zero-angle of attack lift coefficient | 0.8212 | Dimensionless | |
Lift coefficient due to angle of attack | 5.105 | 1/rad | |
| Zero-angle of attack drag coefficient | 0.025455 | Dimensionless |
K | Drag coefficient due to lift | 0.04831 | Dimensionless |
| Air density | 1.225 | kg/m3 |
| Aircraft reference area | 427.8 | m2 |
I begin the control system design procedure by creating a linear model of the system. Because the plant model of Eqs. 7.1 and 7.2 is nonlinear, I must linearize the model about an appropriate equilibrium point. Recall that the equilibrium condition occurs when all the derivatives in Eq. 7.1 are zero.
The equilibrium of interest occurs when the aircraft is trimmed on the glideslope (γ = -3° (-0.05236 radians), herr = 0) and moving at the commanded velocity V = Vc = 81.8 m/s. The values of a and T at the equilibrium condition can be found by setting the left-hand side of Eq. 7.1 to zero and performing an iterative search to solve the first two rows for the two unknowns α and T. This equilibrium solution occurs when α = 2.6857° (0.04688 radians) and T = 42,387 newtons.
For this example, the equilibrium solution was obtained from a Simulink model of the nonlinear system. The trim() command was used at the MATLAB command line to determine the equilibrium conditions. See Glideslope.mdl in the Examples\Ch07 directory of the accompanying CD-ROM for the nonlinear Simulink model and dotrim.m for the trim procedure.
With the equilibrium conditions determined, the next step is to develop a linear plant model. A standard engineering approach is to use calculus techniques to linearize the nonlinear plant model about the equilibrium point. An alternative is to execute linmod() at the MATLAB command prompt to derive the linearized model from the Simulink diagram. The latter approach was used in this case. In addition to determining the equilibrium solution, the dotrim.m script extracts the linear plant model from the diagram. The resulting model appears in Eq. 7.3, where x = [V γ herr]T and u = [α T]T.
| (7.3) |
I should mention a couple of points of numerical interest regarding Eq. 7.3. First, the lower left element of the A matrix is extremely small and probably should be zero. It is not exactly zero because finite-precision floating-point mathematics is subject to small errors such as this. Second, the two upper elements in the right column of the B matrix are quite small as well but should not be assumed to be zero. These elements are much smaller than the others because the magnitude of variations in T is about six orders of magnitude larger than the variations in a. It's usually best to use the models produced by linmod() directly and not try to tweak them by setting apparently insignificant matrix elements to zero.
I will use the pole placement design technique to develop a MIMO observer-controller for the linear plant model of Eq. 7.3. The control system will use the configuration shown in Figure 5.1, except that the reference inputs are zero. With no reference input, the controller functions as a regulator and attempts to maintain all the states at zero deviation from the equilibrium values.
| Note | It is important to remember that the inputs and states for the model of Eq. 7.3 are defined as deviations from the equilibrium values. For instance, if the instantaneous thrust is 42,400 N, the T input to Eq. 7.3 would be (42,400 - 42,387) = 13 N. Similarly, if the state V = 0.2, the actual velocity represented by that value is (81.8 + 0.2) = 82.0 m/s. |
The design specifications for the controller are given as a settling time of 2 seconds and a damping ratio of 0.8. The output of the design process will be the controller gain matrix K and the observer gain matrix L, which will be used in the controller structure of Figure 5.1.
The first step is to create a state-space plant model.
>> a = [-0.018001 -9.7966 0; ...
0.0029251 -0.0062765 0; ...
0 81.912 0];
>> b = [-4.8374 5.2573e-6; ...
0.57862 3.0149e-9; ...
0 0];
>> c = eye(3);
>> d = 0;
>> ssplant = ss(a, b, c, d);
The plant's controllability and observability must be checked as described in Chapter 5. Those tests will not be repeated here. This plant model passes both of the tests.
This is a third-order plant, so three closed-loop pole locations must be selected. The select_poles() function performs this task.
>> n = 3;
>> t_settle = 2;
>> damp_ratio = 0.8;
>> p = select_poles(n, t_settle, damp_ratio);
The controller gain is computed with a call to place().
>> K = place(ssplant.a, ssplant.b, p)
Next, the observer poles must be selected. I will multiply the closed-loop poles by 3 to locate the observer poles. A second call to place() determines the observer gain matrix L.
>> obs_pole_mult = 3;
>> q = obs_pole_mult * p;
>> L = place(ssplant.a', ssplant.c', q)'
The design of the controller and observer gains is now complete. The next step is to form the state-space observer-controller.
>> ssobsctrl = ss(ssplant.a-L*ssplant.c,
[L ssplant.b-L*ssplant.d], -K, 0);
Form an augmented plant model that passes the plant inputs as additional outputs. This modified model is needed because the observer-controller requires the plant inputs in addition to the plant outputs. In the following statements, n is the number of plant states (3 in this case)and r is the number of inputs (2 in this case).
>> r = size(b, 2); % Number of inputs
>> n = size(a, 1); % Number of states
>> ssplant_aug = ss(ssplant.a, ssplant.b, [ssplant.c; zeros(r, n)], ...
[ssplant.d; eye(r)]);
Form a closed-loop system consisting of the augmented plant model and the observer-controller.
>> sscl = feedback(ssplant_aug, ssobsctrl, +1);
Plot the pole locations along with the design specifications. The resulting plot appears in Figure 7.2.
Figure 7.2: Closed-loop pole locations.
>> plot_poles(sscl, t_settle, damp_ratio);
The second example (Example 7.3) was an inverted pendulum on a cart. This system has one input (the cart driving force F) and two outputs (the cart position x and the pendulum angle θ). I used optimal design techniques to develop an observer-controller for this system.
Example 7.3: Inverted pendulum on a cart.
This example involves a plant with one input and two outputs. I use optimal control design techniques (Chapter 6) to develop the observer-controller.
An inverted pendulum on a cart is shown in Figure 7.3. The input to this plant is the cart driving force F (newtons). An electric motor provides the force. The motor is assumed to have a fast response compared with the cart and pendulum, so I do not model it.
Figure 7.3: Inverted pendulum on a cart.
This is a fourth-order plant. Two of the system states are the cart's horizontal position x (meters) and the angle of the pendulum θ (radians) with respect to the vertical. The two remaining states are the cart velocity (m/s) and the pendulum angular velocity (rad/s).
The pendulum is a thin cylindrical rod with length l = 1 meter and mass mp = 2.5 kilograms. It pivots freely about the pin attaching it to the cart. The cart mass mc = 5 kilograms. The dynamics of the wheels will be ignored.
The control system must keep the pendulum inverted by balancing it by moving the cart and simultaneously driving the cart to a commanded position. These objectives conflict to some degree, so this presents an interesting control design problem. The controller design goal is to provide a cart position settling time of 1.5 seconds with satisfactory damping characteristics.
The first step is to develop a plant model for this system. One approach is to use the equations of motion of rigid bodies to analyze the system and produce a set of nonlinear system equations. An alternative is to use software tools to perform this analysis for you. I take the latter approach.
SimMechanics is an add-on product for Simulink from The MathWorks, Inc., that supports the modeling of mechanical systems such as the inverted pendulum and cart. The user constructs block diagrams of a mechanical system, in which the blocks define the properties of rigid bodies and the joints that connect them.
The inverted pendulum system contains two bodies: the cart and the pendulum. For modeling purposes, we assume the cart is connected to the ground by a sliding joint. The pendulum is connected to the cart by a pivot joint. The Simulink/SimMechanics block diagram of the inverted pendulum appears in Figure 7.4.
Figure 7.4: Simulink/SimMechanics model of the inverted pendulum on a cart.
The diagram of Figure 7.4 is contained in the file Pendulum.mdl in the Examples\Ch07 directory of the accompanying CD-ROM. The input to this model is the cart force F and the outputs are y = [θ x ]T. This model is nonlinear, so to design a controller, I must select an equilibrium point and linearize the system about it. The equilibrium point of interest is when all of the states are zero, which occurs when the cart is stopped at the origin of the x-axis and the pendulum is straight up and stationary.
The linmod() command extracts the linear model from the Simulink diagram. The following commands extract the model matrices and create a state-space plant model.
>> [a,b,c,d] = linmod('Pendulum');
>> ssplant = ss(a,b,c,d);
The resulting linear model is shown in Eq. 7.4. The states are x = [θ x ]T, the outputs are y = [θ x ]T, and the input is u = F.
| (7.4) |
The controller reference input is the desired cart position (meters), and the output is the commanded cart drive force F (newtons). As mentioned previously, the dynamics of the power amplifier and drive motor are assumed to be fast relative to the plant dynamics and will be ignored.
I will use optimal design techniques to develop controller and observer gains for this system. The control system structure is shown in Figure 5.1. I will develop the controller gain K by the LQR method, and the observer will function as a steady-state Kalman filter.
I begin with the controller gain design. This is a fourth-order system with one input, so the Q weighting matrix is 4 × 4 and Ris a scalar. Because only the relative values of Q and R affect the results, I fix R = 1 and vary the diagonal elements of Q during the design process.
As a first attempt, create a linear plant model, set Q to be the identity matrix, and design a controller gain with lqr().
>> Q = eye(4);
>> R = 1;
>> K = lqr(ssplant, Q, R)
This sequence of commands produces the controller gain K = [-165.9 -1 -38.2 -4.5].
Form a closed-loop system using the controller gain.
>> cl_sys = feedback(ssplant, -K*inv(ssplant.c), +1);
| Note | In the previous statement, the inverse of the plant C matrix was used to recover the states from the plant outputs. This is possible in this example because the output equation y = Cx can be rewritten as x = C-1y since C is a square and invertible matrix. This simplification is not generally available because C does not usually possess these properties. |
Plot the closed-loop pole locations.
>> t_settle = 1.5;
>> damp_ratio = 0.7;
>> plot_poles(cl_sys, t_settle, damp_ratio);
The resulting plot appears in Figure 7.5. As the figure shows, two of the poles are much slower than the design specifications.
Figure 7.5: Closed-loop pole locations for initial controller design.
Try increasing the element of Q corresponding to the x state by a factor of 10.
>> Q = diag([1 10 1 1]);
Repeat the steps to compute the controller gain, form a closed-loop system, and plot the closed-loop pole locations. As the new plot will show, the slowest pole locations have moved to the left, but not far enough. After a few more iterations, the following Q matrix is found to produce a response satisfying the design requirement.
>> Q = diag([1 1e6 1 1]);
This system has the gain vector K = [-2012 -1000 -520 -629]. The closed-loop poles, along with 1.5-second settling time and 0.7 damping ratio constraints, appear in Figure 7.6.
Figure 7.6: Closed-loop pole locations with settling time and damping ratio constraints.
Observe in Figure 7.6 that two of the poles are much faster (with a real part of -9.55) than the other two poles (which have a real part of -3.84). The faster poles correspond to motion of the pendulum relative to the cart. The slower poles relate to moving the cart to the commanded position.
The next step is to design the observer gain. Assume that the cart position x and pendulum angle θ are measured with sensors that provide quantized samples. It will be necessary to estimate the cart velocity and pendulum angular rate from those measurements. This linear system model is shown in Eq. 7.5, where
x = [θ x ]T, y = [θ x]T, and u = F.
| (7.5) |
Assume the system contains nonlinearities, such as friction and random disturbances applied to the cart and pendulum. A simplified model of these combined effects is a random force applied to the cart in addition to the force generated by the drive motor. In this model, the G matrix equals the B matrix shown in Eq. 7.5, and the H matrix is zero. Let the standard deviation of this random force be 0.1 newtons. The process noise covariance is then QN = 0.12.
Assume the quantization least significant bit (LSB) for the cart position is 0.01 meters and the quantization of θ is 0.0044 radians (0.25°). This results in measurement noise covariance
With this information, the steady-state Kalman observer gain is computed with the following MATLAB commands.
>> a = ssplant.a; b = ssplant.b;
>> c = [1 0 0 0; 0 1 0 0]; d = [0; 0];
>> g = b;
>> h = d;
>> obs_plant = ss(a, [b g], c, [d h]);
>> QN = 0.1^2;
>> RN = [0.01^2/12 0; 0 0.0044^2/12];
>> [kest, L] = kalman(obs_plant, QN, RN);
The resulting gain matrix is
Form a state-space observer-controller.
>> ssobsctrl = ss(a-L*c, [L b-L*d], -K, 0);
Augment the plant model to pass its inputs as additional outputs.
>> r = size(b, 2); % Number of inputs
>> n = size(a, 1); % Number of states
>> ssplant_aug = ss(a, b, [c; zeros(r, n)], [d; eye(r)]);
The reference input provides a commanded value for the x state. The next step is to determine the feedforward gain N that results in zero steady-state error to a step input. To do this, form the state-space model of Eq. 7.5 and extract from it a SISO system that describes the transfer function between the plant input and the x output. Note that the x state is the second element of the y vector.
>> ssplant2 = ss(a, b, c, d);
>> sys_r_to_x = ssplant2(2,1);
Compute the feedforward gain with the technique described in Section 5.8.
>> Nxu = inv([sys_r_to_x.a sys_r_to_x.b; ...
sys_r_to_x.c sys_r_to_x.d]) * [zeros(n,1); 1];
>> Nx = Nxu(1:n);
>> Nu = Nxu(end);
>> N = Nu + K*Nx;
The value of N resulting from this computation is -1000.
Form a closed-loop system with positive feedback.
>> sscl = N*feedback(ssplant_aug, ssobsctrl, +1);
Plot the pole locations of the closed-loop system (Figure 7.7).
Figure 7.7: Closed-loop pole locations including the observer-controller.
>> plot_poles(sscl, t_settle, damp_ratio);
Notice that two of the observer poles in Figure 7.7 are slower than the settling time requirement, though all of them satisfy the damping ratio requirement.
Finally, display the step response of the linear closed-loop system, which is the response to a step of 1 meter in the x cart position command. The first command shown below specifies names for the sscl model inputs and outputs. The second command plots the step response, including input and output names. The resulting plot appears in Figure 7.8.
Figure 7.8: Closed-loop step response.
>> set(sscl,'InputName','R (m)', 'OutputName', {'Theta (rad)', ...
'X (m)', 'Force (N)'});
>> step(sscl)
As you can see from Figure 7.8, a step command of 1 meter in r results in an initial cart force of 1,000 N in the -x direction. This motion deflects the pendulum in the +θ direction before the cart begins moving in the +x direction. The cart then moves to the x = 1 position in a manner that brings the pendulum back into vertical balance just as the cart arrives at its destination.
Observe that the θ plot (the top trace in Figure 7.8) shows a peak amplitude of 37° (0.65 radians) at 0.16 seconds after the step() command. This large angular deflection could violate the assumptions of the linear model to such a degree that the control system would fail to function properly with the nonlinear plant. The linear model assumes the angle θ is small enough that sinθ ≈ θ and cos θ ≈ 1, which might not be reasonable for this large of a deflection. In Chapter 9 (control system implementation and testing), I explore testing strategies for determining whether a controller will behave appropriately under all anticipated operating conditions.
These examples demonstrate the application of the state-space design techniques introduced in Chapters 5 (pole placement) and 6 (optimal design) to MIMO plants. As demonstrated by these examples, the basic design approach does not change very much in moving from a SISO plant to a MIMO plant. The state-space design approaches are ideal for use with complex MIMO plants as well as with simpler SISO systems.
No comments:
Post a Comment