CHAPTER 5 A COMPUTATIONAL SCHEME FOR SOLVING THE
5.4 Two-Phase Scheme for Solving TOCP
5.5.3 F-8 Fighter Aircraft
The F-8 fighter aircraft has been considered in several pioneering studies (e.g., Kaya and Noakes, 1996; Banks and Mhana, 1992; Simakov et al., 2002) and has become a standard for testing various optimal control strategies. A nonlinear dynamic model of the F-8 fighter aircraft is considered below. The model is represented in state space by the following differential equations: control input u represents the tail deflection angle. For convenience of comparison, the standard settings (Kaya and Noakes, 1996; Lee et al., 1997) are used. A control |u| ≤ 0.05236 must be found that brings the system from its initial state
[
26.7 180 , 0, 0π ]
T to the final state[
0, 0, 0]
T in minimum time.When the two-phase scheme is applied, as described in Section 5.4, the switching times computed in Phase I are 0.115, 2.067, 2.239, 4.995, and 5.282, and the terminal time is tf = 5.7417. These switching data are used to set the design variables and their corresponding bounds, and then the problem is solved by the mixed integer NLP method. Finally, the switching times for the discrete control input are 0.098, 2.027, 2.199, 4.944, and 5.265, and the terminal time tf is 5.74216. Figure 5.6 shows the comparison of the controls between Phase I and Phase II, while Figure 5.7 shows the trajectories of the states and the control of Phase I and Phase II. This example is also solved by Kaya and Noakes (1996) using the switching time computation method and by Lee et al. (1997) using the Control Parameterization Enhancing Transform (CPET) method. Table 5.1 shows the terminal time tf, switching times and the accuracy of terminal constraints computed by various methods for this problem. According to the numerical results, the two-phase scheme provides a better solution, and the accuracy of the terminal constraints is acceptable.
5.6 Summary
This chapter has proposed a novel numerical method for solving time-optimal control problems with discrete-type control inputs that include the bang-bang type most commonly encountered when the control is bounded. This two-phase computational scheme for finding a discrete optimal control for time-optimal control problems is novel because its discrete control can be more easily implemented than continuous control in practical applications. A simple example, a third-order system, was presented to demonstrate the usage of the proposed scheme. A flexible mechanism control problem and an F-8 fighter aircraft control problem were also considered and solved by application of the proposed scheme. Numerical results were obtained efficiently and accurately and provide evidence that the two-phase scheme constitutes a viable method for solving time-optimal control problems with discrete-valued controls.
Table 5.1 Results of various methods for the F-8 fight aircraft problem.
Method
t
f Switching TimesAccuracy of Terminal Constraints STC
(Kaya and Noakes, 1996)
6.3867 0.0761, 5.4672, 5.8241, 6.3867 ≤ 10-5
CPET
(Lee et al., 1997)
6.0350 2.188, 2.352, 5.233, 5.563 ≤ 10-10
Two-phase scheme 5.7422 0.098, 2.027, 2.199, 4.944, 5.265 ≤ 10-10
Table 5.2 Optimal results for the fourth-order system.
u(t0) NGP INTP NIT
J
0* Conv.Par. CPUZero 5 4.33196 1.04E-05 0.131
First 5 4.86764 5.69E-06 0.07
5
Cubic 5 4.90565 4.67E-07 0.06
Zero 15 4.30699 5.71E-09 1.382
First 7 4.28066 5.65E-06 0.35
11
Cubic 10 4.30041 1.52E-08 0.34
Zero 50 4.26239 1.38E-07 43.803
First 44 4.22087 3.50E-07 18.596
0.0
51
Cubic 40 4.22187 1.10E-08 12.659
Zero 6 4.33197 6.16E-07 0.12
First 7 4.86765 2.25E-07 0.091
5
Cubic 9 4.90560 3.72E-05 0.12
Zero 7 4.36249 3.51E-08 0.651
First 8 4.28064 1.09E-05 0.43
11
Cubic 10 4.30041 3.50E-07 0.39
Zero 49 4.26229 7.21E-08 42.872
First 42 4.22087 2.20E-06 17.315
1.0
51
Cubic 38 4.22187 2.88E-06 11.847
ū
iū
i+1J
0P
jContinuous optimum point
ū
iū
i+1New lower bound New upper bound
Original design domain
(I) (II) (III)
Figure 5.1 Conceptual layout of the branching process.
Initialization
1. Relax all discrete-valued restrictions 2. Place the resulting continuous NLP
problem on the branching tree.
3. Set the cost bound Jmax=
Is the branching tree empty?
A
Select an unexplored node from branching tree
Solve subproblem by applying the OCP Solver (AOCP)
(NLP is optimal) &&
(J0<Jmax)
Evaluate the values of criteria for selecting branch object
Decide the branch of design variable
Branching process Set new cost bound
Jmax= J0
Show results
Figure 5.2 Flow chart of the algorithm for solving discrete-valued optimal control problems.
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Mixed Discrete NLP method (tf = 0.7813) AOCP (tf = 0.7826)
(a) AOCP vs. a mixed-integer NLP method.
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0
Figure 5.3 Control trajectories for the third-order system.
0.0 1.0 2.0 3.0 4.0 5.0
Figure 5.4 State trajectories for the fourth-order system (Phase II).
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 Time(sec.)
-1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
Control input (N-m)
Phase I (tf = 4.2398) Phase II (tf=4.21794)
Figure 5.5 Control trajectories for the fourth-order system.
0.0 1.0 2.0 3.0 4.0 5.0 6.0 Time(sec.)
-0.06 -0.05 -0.04 -0.03 -0.02 -0.01 0.00 0.01 0.02 0.03 0.04 0.05 0.06
Tail deflection angle (rad.)
Phase I (tf = 5.74173) Phase II (tf = 5.74216)
Figure 5.6 Control trajectories for the F-8 fighter aircraft.
0.0 1.0 2.0 3.0 4.0 5.0 6.0
Figure 5.7 Trajectories of the states and control input for the F-8 fighter aircraft.
CHAPTER 6
ENGINEERING APPLICATIONS
6.1 Flight Level Control Problem
The flight level tracking that plays an important role in autopilot systems has received considerable attentions from many researchers (Lygeros, 2003; Lygeros et al., 1999; Cook, 1997; Tomlin et al., 1996; Etkin and Redi, 1996). A commercial aircraft‘s cruising altitude is typically assigned a flight level by air traffic control (ATC). To ensure aircraft separation, each aircraft has its own flight level separated by a few hundred feet; however, changes in flight level do happen occasionally and must be cleared by ATC. At all other times, the aircraft crew must ensure that they remain within the allowed bounds of their assigned level.
At the same time, they must also maintain limits on factors such as speed, flight path angle, and acceleration imposed by limitations of airframe and engine and passenger comfort requirements or to avoid dangerous situations such as aerodynamic stall. In this paper, the flight level tracking problem is formulated into an optimal control problem. For safety reasons, the speed of the aircraft and the flight path angle must be kept within a safe “aerodynamic envelope” (Tomlin et al., 1996) that can be translated into the dynamic constraints of the optimal control problem. A flight level tracking problem and a minimum time problem are outlined in the following sections and then solved using the proposed solver.
6.1.1 Aircraft Model
Much ATC research (e.g., Cook, 1997; Etkin and Redi, 1996) has applied a point mass model to describe aircraft motion, considering only aircraft movement in a lateral direction. In Figure 6.1, three coordinate frames are used to describe aircraft motion: Xg-Yg denotes the ground frame; Xb-Yb,the body frame; and Xw-Yw, the wind frame. In addition,
θ
,γ
, andα
denote the rotation angle between the frames; V∈ \ represents the speed of the aircraft, which is aligned with the positive Xw direction; and h is the aircraft’s altitude.The equations of the motion can be derived from the force balance relationships:
where T is the thrust exerted by the engine, D is the aerodynamic drag, and L is the aerodynamic lift. By applying basic aerodynamics, the lift (L) and drag (D) can be approximated by
where CL, CD, and c are dimension-less lift and drag coefficients, s is the wing surface area and
ρ
is the air density.According to the admissible optimal control formulation described in Section 3.4, the air model can be formulated by a three-state model with a state variable vector x(t) = [x1, x2, x3]T
= [V,
γ
, h]T and a control input vector u(t) = [u1, u2]T = [T,θ ]
T. By approximatingα
with a small angle, the equations of the motion (system equations) can be written as1 2
This model, proposed by Lygeros et al. (1999) and adopted here, extends the three dimensions of an aerodynamic envelope protection problem. Taking into the consideration of safety conditions, the aircraft speed and flight path angle are bounded in a rectangular limitation called a “safe aerodynamic envelop.” Following Tomlin et al. (1996), Lygeros (2003) proposed a simplified aerodynamic envelope that is adopted in this paper and translated into the following dynamic constraints:
min 1 max
Based on the NLP formulation described in Section 2.2, these constraints can be treated as dynamic constraints and rewritten as follows:
1 1 min 2 1 max
To illustrate the capabilities of the proposed method, the flight level tracking problem and the minimum time problem have been chosen.
Case I: Flight level tracking problem
This tracking problem is to find the controls that will maintain the system state x(t) as close as possible to the desired state r(t) in the interval [t0, tf]. The performance index for the tracking problem can be written as
0
2
0 tf
( ) ( )
( )tJ = ∫
tx t − r t
Qdt
(6.6)where Q(t) is a real symmetric n × n matrix that is positive semi-definite for all
t ∈⎣ ⎡ t
0, t
f⎤ ⎦
. The flight level tracking problem involves keeping the aircraft as near as possible to the desired level and aircraft speed. Therefore, the performance index can be represented as( ) ( ) ( )
where x1d is the desired aircraft speed, x2d is desired flight path angle and x3d is the assigned altitude.
Case II: Minimum time problem
The minimum time problem is to transfer a system from an arbitrary initial state x(t0) = x0
to a specified target set St in minimum time. The performance index for the minimum time
problem can be written as
0 0 0
tf
f t
J = − = t t ∫ d t
(6.8)where tf is the first instant of time when x(t) and St intersect. In some emergencies, the aircraft crew is asked to change their level as soon as possible.
6.1.2 Numerical examples
The following parameters, outlined here for case I, are used in both cases:
a
L = 65.3 Kg/m,a
D = 3.18 Kg/m,m = 160×10
3 Kg,g = 9.81 m/s
2,θ
min = -20°, γmin = -20,c = 6, θ
max = 25°, γmax = 25,T
min = 60×103 N,T
min = 120×103 N,V
min = 92 m/s,V
max = 170 m/s,h
min = -150 m,h
max = 150 m The initial values of the state variables arex
0 = [100, 20, -120]T (6.9)and the purpose of this problem is to find a suitable control for maintaining the flight level and keeping the aircraft altitude at the assigned level. Thus the desired states are set with following values
r(t) = [150, 0, 0]
T. (6.10)In addition to the dynamic constraints proposed in Eq. (6.5), the control inputs are also limited within the following bounds:
min 1 max
Substituting these parameters into Eqs. (6.3) and (6.7), the flight tracking problem is solved by the OCP solver. The numerical results are shown in Figure 6.2. As shown in Figure 6.2(a), all states meet the constraints, and the flight level and aircraft speed return to the
desired states. Table 6.1 shows the user subroutines for this case. Obviously, the OCP solver provides an easily usable tool for solving dynamic optimization problems.
Case II: Minimum time problem
In this problem, the aircraft crew is asked to increase their altitude in minimum time. The initial and final altitude are h0 = 0 m and hf
= 500 m, respectively. All constraints imposed on
case I are also imposed on this case. The initial state x0 = [100, 0, 0]T, the final time, tf, obtained by using the AOCP, is 73.98 seconds, and the final altitude is 499.928 m. The control histories shown in Figure 6.3(a) and (b) give the state trajectories, which, as the figure illustrates, all fall within the safe “aerodynamic envelope” (i.e., meet the dynamic constraints).6.2 Vehicle Suspension Design Problem
Many studies have treated the vehicle as a dynamic system, starting with the basic properties of vehicle suspension, the stiffness and damping coefficients (Gillespie, 1992).
Thus the design of vehicle suspension systems has received much attention in the automotive industry. Numerous researchers have examined semi-active and active vibration isolation for suspension systems. Yet, despite recent advances in active and semi-active suspension technology, vehicles with passive suspension systems still dominate current car production.
Tools must therefore be made available to vehicle designers for optimizing passive suspension systems.
The model described here is a half-car model that allows independent vertical inputs to the front and rear wheels and can thus simulate pitching and bouncing motions due to road inputs.
Two longitudinal forces, which can be positive to represent traction or negative to represent braking, are applied to the front and rear axles to simulate the effects of vehicle acceleration or deceleration. Cases of braking and accelerating while moving straight ahead are used to validate the longitudinal vehicle dynamics.
An optimal design problem in relation to vehicle suspension is considered to maximize vehicle ride performance, which may be evaluated according to passenger discomfort. The response to driver’s seat to acceleration is commonly used as the objective of suspension design. Three road profiles that excite pitch and bounce motions at a constant vehicle speed are used to calculate the optimal suspension parameters. In this optimal design problem, the objective is to minimize the extreme acceleration of the driver’s seat under a number of constraints on the dynamic response and the design parameters. The optimal design of a vehicle suspension system can be applied in diverse fields of research, including traction force control, speed control, braking system design, to name a few. In this dissertation, an emergency stop – a special case of vehicle speed control problem – is treated as a time-optimal control problem and solved by the proposed AOCP method.
6.2.1 Derivation of the Vehicle Model Half-car model
Although the quarter-car model has been commonly used in assessing vehicle ride performance, it does not fully represent the rigid body motions that a motor vehicle may exhibit. For example, the quarter-car model disregards pitching motions, which may be important, particularly when the car travels over obstructions like road bumps and potholes.
Moreover, the quarter-car model is a multi-input system that responds with both pitch motions and vertical bounce because of the longitudinal distance between the axles. These pitch and bounce motions must be understood because they provide useful information on vertical and longitudinal vibrations. As a result of these quarter-car limitations, half-car and full-car models are used in several studies on suspension. Figure 6.4 depicts a nonlinear half-car model with six degrees of freedom, modified from the model of Haug and Arora (1979). Two additional longitudinal forces, traction or braking forces, are applied to the axles, allowing the vehicle to be accelerated or decelerated. Shock absorbers are assumed to be rigidly joined to
the chassis without displacement or deflection in the longitudinal direction. Based on this assumption, the longitudinal forces only change the speed and pitch angle of vehicle. The governing equations for the vehicle can be derived from Lagrange’s equations
T T V
where T and V represent the system’s kinetic and potential energy, and Qi represent nonconservative generalized forces. In Figure 6.4, the kinetic energy of the system can be expressed as
The potential energy V of the conservative forces is
2 2
and the virtual work done by the nonconservative forces is
1( 2 3 1)( 2 3 1 the system equations describing the motion of half-vehicle model can be derived as follows
1 1 1 1 1 2 1 3 1 1 1 2 1 3 0
−
k z
2 4−k z
3 5 =0,where mi represents the masses of the seat and driver, the main body, and the wheel and axles, respectively. The parameters ki and ci represent the known stiffness and damping coefficients of the suspension system. The moment of main body inertia about its center of mass is denoted as I, while H is the vertical distance from the center of gravity (C.G.) to the ground, and L is the total length of the wheel base. The functions R1(y) and R2(y) represent displacements of the front and rear wheels, caused by undulations of the road surface on which the vehicle is traveling. Once
z
6+i =z i
i, =1,...,6 is defined, the vehicle system can be transformed into a state-space equation of the form= + +
x Ax Bu W (6.22)
where
x(t) = [z
1, z2, z3, …, z12]T represents the vector of state variables and the nonzero elements of matrices A, B and W are given as follows:( )
For safety and comfort, six dynamic constraints are imposed on the system, whose constraint equations may be written as
1( ) 1, 0 f
5 2 3 4
where θ2 to θ6 are the maximum allowable displacements and vallow is the maximum allowable speed.
Road surface displacement function
Because the dynamic response depends strongly on the vertical displacement history of the wheels on the road surface, the input road conditions are very important. Most data used in establishing the ride comfort criteria were obtained using sinusoidal inputs. Thus the road surface displacement function plotted in Figure 6.5 is defined as a sinusoidal undulation with amplitude x0 and variable half-wavelength li (Haug and Arora, 1979). The front tire displacement v( y ) at position y is thus defined as
1 1 function for the front wheel can therefore be defined as
1
where yt is the final position of the road undulation. The vertical displacement of the rear wheel has the same value as that of the front wheel but with a wheelbase lag. Therefore,
2( ) 1( )
R y =R y L− (6.32)
where R1(y) is defined in Eq. (6.31).
6.2.2 Numerical Examples
This paper uses the numerical data from Haug and Arora (1979) to validate the model specified in Section 6.2.1. The following parameters in the vehicle system equations are fixed during the calculations; m1
g = 290 lb, m
2g = 4500 lb, m
4g = m
5g = 96.6 lb, I = 41,000
lb-in-sec2, H = 20 in, L = 120 in, k4 = k5 = 1500 lb / in, vallow = 1056 in / sec (60 mph), and c4= c5 = 5 lb-sec / in. The coefficients of the suspension system are selected as design variables,
b = [k
1, k2, k3, c1, c2, c3]T. The lower and upper bounds on b are [50, 200, 200, 2, 5, 5]T and [500, 1000, 1000, 50, 80, 80]T, respectively. The maximum allowable values for the state variable constraints in Eqs. (6.23) – (6.29) are selected to be [400, 2, 5, 5, 2, 2, 1056]T. The units of z1, z2, z4, z5 and z6 are inches and those of z3 are radians.Model validation
The physical phenomenon of rigid body motion can be used to confirm the correctness of the vehicle model. Therefore, cases of braking and accelerating while traveling straight ahead are considered here to validate the longitudinal vehicle dynamics. For convenience of observation, the vehicle is assumed to travel along a straight path such that R1(y) = R2(y) = 0.
In cases of acceleration, the control problem is to determine a feasible acceleration trajectory along which a vehicle with various initial speeds can arrive at a destination in minimal time.
Hence, one additional terminal constraint is imposed:
6( )f t
z t
=y
(6.33)where yt is the destination. Similarly, one additional terminal constraint is included in cases of braking
12( ) 0f
z t
= (6.34)All the acceleration and braking test cases are transformed into time-optimal control
problems and solved by applying the proposed NLP method. Figure 6.6 and Figure 6.7 show the velocity trajectories of the vehicle with various starting speeds. As Figure 6.6 illustrates, the vehicle accelerates at the maximum allowable acceleration until the speed constraint defined in Eq. (6.29) becomes pertinent, from which point the speed is maintained. In contrast to the cases of acceleration, the vehicle decelerates with maximal allowable deceleration until it stops. Figure 6.8 shows the driver’s seat acceleration trajectory for the case of straight-ahead braking. According to these results, the vehicle motion is consistent with the motion of a rigid body, meaning that the longitudinal vehicle dynamics of the proposed model are validated.
Optimal design of the vehicle suspension system
The vertical displacement functions and system equations specified in Section 3 can be used to define an optimal suspension design problem. The driver is to be made as comfortable as possible over a range of road conditions and traveling speeds. Thus, the design objective is to minimize the maximum absolute acceleration of the driver’s seat by adjusting the vehicular suspension properties subject to the constraints that certain relative displacements do not exceed imposed limits. The objective function is therefore
0 max[0, ] 1( )
Two design cases considered by Haug and Arora (1979) are used here to examine the correctness of the proposed method. Figure 6.9 represents the road displacement profiles in the test cases. In case 1, the road surface profile includes a cavity. Case 2 involves two road displacement profiles, presented in Figure 6.9(b) and (c). The speed of the vehicle in case 1 is 450 in/sec and that in case 2 is 960 in/sec. Table 6.2 gives the optimal solutions. A comparison with the results present in the research sources (Hsieh and Arora, 1984; Haug and Arora, 1979)
shows that the results obtained by the proposed method are quite accurate.
Vehicle speed control problem
In most emergency situations, drivers must stop the vehicle quickly and safely. Changing the speed of the vehicle according to the conditions of the road and the distance from the current position to the site of accident is a vehicle speed control problem that the vehicle model and system equations derived in Section 3 can be used to solve. In this case, the initial speed of the vehicle is 880 in / sec (50 mph) and the road surface profile is as plotted in Figure 6.10. According to the definition in Section 2.2, the emergency braking problem is transformed into a time-optimal control problem that is then solved using the proposed NLP method. The minimum time, tf = 3.4 seconds, and the terminal displacement, z6 = 1585.7 inches, are obtained using the OCP solver. Figure 6.11 shows the trajectories of the vehicle speed and acceleration. Figure 6.12 plots the trajectories of the acceleration and pitch angle of the passenger seat, which are of interest to vehicle designers. The solid bold curves at the bottom of Figure 6.11 and Figure 6.12 represent the corresponding road profile. The numerical results indicate that all the constraints are satisfied and the optimal control law that solves the emergency braking problem is determined.
6.3 Summary
In this chapter, two practical applications, the flight level control problem and the vehicle suspension design problem – both highly nonlinear optimal control problems – have been formulated following the procedure suggested in this dissertation and solved by the proposed OCP solver. In the case of the flight level control problem, two common types of optimal control problem, the tracking problem and the minimum time problem, were derived to simulate practical situations. The vehicle suspension design problem provided a useful
In this chapter, two practical applications, the flight level control problem and the vehicle suspension design problem – both highly nonlinear optimal control problems – have been formulated following the procedure suggested in this dissertation and solved by the proposed OCP solver. In the case of the flight level control problem, two common types of optimal control problem, the tracking problem and the minimum time problem, were derived to simulate practical situations. The vehicle suspension design problem provided a useful