# Minimax control of flexible structures using quadraticallyconstrained programming

код для вставкиСкачатьMINIMAX CONTROL OF FLEXIBLE STRUCTURES USING QUADRATICALLY CONSTRAINED PROGRAMMING by Jennifer Haggerty January 2010 A thesis submitted to the Faculty of the Graduate School of State University of New York at Buffalo in partial fulfillment of the requirements for the degree of Master of Science Department of Mechanical & Aerospace Engineering UMI Number: 1477164 All rights reserved INFORMATION TO ALL USERS The quality of this reproduction is dependent upon the quality of the copy submitted. In the unlikely event that the author did not send a complete manuscript and there are missing pages, these will be noted. Also, if material had to be removed, a note will indicate the deletion. UMI 1477164 Copyright 2010 by ProQuest LLC. All rights reserved. This edition of the work is protected against unauthorized copying under Title 17, United States Code. ProQuest LLC 789 East Eisenhower Parkway P.O. Box 1346 Ann Arbor, MI 48106-1346 T o My Family ii Acknowledgment The journey to the completion of my Master’s degree would not have been possible without the continuous support from my family. Mom and Dad, I cannot thank you enough for your love and support. You have set such great examples for me, inspiring me with all that you have accomplished and continue to do! You have taught me to always aim high and to believe in myself. With determination and dedication, anything can be done. Thank you so much, from the bottom of my heart! Natalie, you are the best sister a girl could ever have! Thank you so much for sharing in a wonderful friendship. Without your love and support, I would not be where I am today. You are a wonderful example of what I hope to be, and I feel so lucky to have your support. Ryan, I cannot thank you enough. The sacrifices you have made in supporting me throughout this crazy journey have been unbelievable. You have given me confidence when I needed it the most. You have not only improved my vocabulary, you have challenged me, offering different perspectives and supporting me with some of the most difficult decisions I have had to make. I love you so much! I am thankful every day to have you as my family. Grandma, you are my little sweetie-heart! You have helped me so much, not only with the things that you have helped me to do day-to-day, but with the iii Acknowledgment iv motivation you have given me. With all that you do, and the energy you carry with you, you have shown me how productive a person can be. You have inspired me greatly and are such a wonderful role model for me. I hope to be very much like you. I feel fortunate to have you as not only a grandmother, but as such a dear friend. I love you! Thank you so much. It is scary to say that this is my 9th year at UB. My choice to continue on with my education was made easy, due to the extraordinary people that I have come to know in my years here. My interest in research and graduate-level studies began in the DOES Lab, working under Dr. Kemper Lewis, learning about Optimization methods and DecisionBased Design. I feel very lucky to have been influenced and inspired by Dr. Lewis and all of my friends in the DOES Lab. Andy Olewnik, Scott Ferguson, Ashwin Gurnani, Tung-King See, Michael Kulok, are just some of the many amazing friends that have inspired me greatly. Most importantly, I send many thanks to my amazing mentor and advisor, Dr. Tarunraj Singh. I met Dr. Singh in the final year of my Undergraduate Education. Dr. Singh inspired me to want to pursue more than my Undergraduate degree. He gave me purpose and structure, fundamental parts of my life that were lacking when faced with the most trying times of my life. He provided me with guidance and support for which I am forever grateful. Dr. Singh, you have inspired me more that I think you know. I have a great deal of respect for you and all that you continue to accomplish. I thank God everyday for the opportunities you have given me, and all that I continue to learn from you. I feel very lucky to have had a spot in the CODE Lab working alongside some amazing students and wonderful friends. Thank you so much! Acknowledgment v I have learned a great deal from the wonderful friends that I have made in the CODE Lab. Wanseok Yang, Jayaram Gopalakrishnan, Thomas Conord, Matt Vossler, Ravi Kumar, Brandon Brown, Max Rech, Christoph Hoog Antink, Gabriel Terejanu, and Umamaheswara Reddy are just some of the friends that have supported me along the way. I especially want to thank Uma. Uma, you are brilliant and are destined to do many great things! You have been an amazing friend, offering a great deal of guidance at many points in time. I am extremely grateful to know you and wish you all the best. I feel lucky to have had such wonderful teachers, such as Dr. John Crassidis and Dr. Puneet Singla. I see in my great teachers examples of what I hope to be some day. Thank you for your inspiration! So the journey continues on towards a Ph.D. I will continue on with the motivation provided by the many inspiring people in my life. I can’t wait to see what comes next! Contents Preamble viii Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 Introduction ix 1 1.1 Prior Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 2 Review of Techniques 2.1 6 Traditional Minimax Approach . . . . . . . . . . . . . . . . . . . . . 2.1.1 7 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10 2.2 Gradient-Based Optimization . . . . . . . . . . . . . . . . . . . . . . 14 2.3 Convex Optimization Methods . . . . . . . . . . . . . . . . . . . . . . 16 2.3.1 Linear Programming Approach . . . . . . . . . . . . . . . . . 17 2.3.2 LMI Approach . . . . . . . . . . . . . . . . . . . . . . . . . . 26 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 2.4 3 QCP Minimax Formulation 35 vi Contents vii 3.1 Residual Energy Derivation . . . . . . . . . . . . . . . . . . . . . . . 35 3.2 QCP Framework . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37 3.2.1 Minimum Time Solution . . . . . . . . . . . . . . . . . . . . . 39 3.2.2 Preprocessing in Matlab and Solving with GAMs . . . . . . . 41 3.2.3 Flow of Information . . . . . . . . . . . . . . . . . . . . . . . . 42 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3 4 Numerical Examples 4.1 4.2 4.3 4.4 Harmonic Oscillator . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44 4.1.1 Uncertain Stiffness . . . . . . . . . . . . . . . . . . . . . . . . 44 4.1.2 Uncertain Stiffness and Damping . . . . . . . . . . . . . . . . 49 Floating Oscillator . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 4.2.1 Uncertain Stiffness . . . . . . . . . . . . . . . . . . . . . . . . 56 4.2.2 Uncertain Stiffness and Damping . . . . . . . . . . . . . . . . 61 Constrained Robust Control: Simple Elevator Problem . . . . . . . . 68 4.3.1 Simple Elevator Model . . . . . . . . . . . . . . . . . . . . . . 68 4.3.2 Computational Algorithm for Time-Optimal Problem . . . . . 71 4.3.3 QCP Formulation . . . . . . . . . . . . . . . . . . . . . . . . . 74 4.3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88 5 Conclusion 5.1 44 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 89 Contents viii 5.2 Contribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 5.3 Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 Bibliography 95 Abstract The focus of this thesis is on the design of a control strategy for robust control of slewing flexible structures. The minimax problem is posed in a Quadratically Constrained Programming (QCP) framework for determination of the time-optimal control profile, which minimizes the maximum magnitude of the system’s residual energy over the range of uncertain parameter values. A set of quadratic constraints are generated that merge state constraints with the quadratic form of the system’s residual energy equation and incorporate multiple system plants to represent the range of parameter uncertainty. The new formulation utilizes the most powerful large-scale optimization problem solvers, exploits a convex optimization framework, and neither requires parameterization of the control profile and corresponding timevector, nor approximation of the system’s residual energy. The technique is tested on several benchmark problems. Results approximate traditional minimax solutions and illustrate potential for use in approximating the minimax control profile for large-scale systems with parameter uncertainties. ix Chapter 1 Introduction 1.1 Prior Work Various techniques for vibration control of slewing flexible structures have been developed over the past half a century that address the input-induced vibrations (oscillations) in linear and quasi-linear systems. The development of these techniques has been motivated by the need to precisely manipulate lightweight flexible structures that experience oscillatory behavior due to the excitation of underdamped modes of the system by a reference input, such as robotic arms [2], gantry cranes [26], hard disk drives [11], spacecraft and space structures [7], high speed elevators [16], high speed tape drives [13], slosh models [31], vehicle platoon models [1] etc. While some feedback and feedforward controllers can be used to achieve rapid maneuvering in vibratory structures, problems arise with use of model-based control strategies. First, uncertainties in natural frequency and damping, and the associated underdamped modes of vibration can affect stability and performance of the controlled system. Second, unmodeled high frequency modes affect stability in some feedback 1 Chapter 1. Introduction 2 control schemes. Several methods, which subsequently incorporate robustness to modeling uncertainties, have stemmed from the concept of input shaping and have been used to address various issues in the control of vibratory systems. According to Singh [19], input shaping exploits the idea that precise position control requires minimization of vibration at the end of the maneuver. A logical approach, then, is to shape the input in order to satisfy the final value constraints, given that a designer can anticipate the response of an open-loop system to a specified input. Since the plant receiving the reference input can either represent an open-loop system or a closed-loop system, as long as the plant is stable, or marginally stable, the concept of input shaping could be applied widely. Input shaping methods first arose after the vibration control research proposed by Tallman and Smith [27] in 1957, called ”Posicast” for single-mode oscillatory systems, a method which involves summing several excited transient oscillations. The result is a steady-state output with zero transient oscillation whose magnitude matches the sum of the excitation magnitudes. Since the ”Posicast” technique relies on precise knowledge of the systems natural frequency, its application is limited to certain underdamped systems. Several years later, Singer and Seering [17] derived the same results as the ”Posicast” technique using an input shaping approach. Additionally, a technique was proposed that addressed the sensitivity of the solution to parameter modeling errors. Robustness to modeling errors was achieved with the inclusion of constraints that force the derivative of the system’s residual energy at the end of the maneuver, with respect to natural frequency or damping, to zero. Singh and Vadali [23] designed a time-delay prefilter that cancels the system poles, resulting in the same results Chapter 1. Introduction 3 as Singer and Seering. In this case, robustness is achieved by cascading multiple instances of the controller to reduce the sensitivity to parameter errors. The development of time-delay filters was later extended to include userspecified time-delays, minimization of the residual energy of multi-mode systems [22], optimal control of systems with rigid body modes and limits on the control authority (saturating controllers) [18]. The development of jerk-limited time-delay filters allowed for the determination of a reference profile which reduces the excitation of higher frequency unmodeled modes. Singh and Vadali [24] addressed robust time-optimal control with the design of a time-delay filter that cancels the poles of the system while satisfying the rigid body boundary conditions. Robustness is achieved by locating multiple zeros at the estimated locations of the poles of the systems. Liu and Wei [9] achieve robustness in time-optimal control by decoupling the equations of motion into rigid and flexible body modes and forcing the partial derivative of the decoupled states with respect to natural frequency to zero at the end of the maneuver. At this point, robustness is achieved near the nominal values of the uncertain parameters only. In 1996, the extra-insensitive input shaper was developed [25], which extended the range of parameter uncertainty addressed by constructing a desired sensitivity curve for the uncertain parameters. The zero vibration constraint at the end of the maneuver is relaxed, instead requiring that the vibration at the end of the maneuver be at a low level, yielding an improvement in robustness overall. Pao et al. [14] proposed another method that includes the probability distribution of the uncertain parameters and applies weights to the ranges of the uncertain parameters at their nominal values according to the expected modeling errors. Finally, the traditional minimax time-delay approach proposed by Singh [18] significantly improved Chapter 1. Introduction 4 robustness by minimizing the maximum magnitude of the system’s residual energy at the end of the maneuver over the range of parameter uncertainty, a crucial method, which significantly improved the level of robustness to parameter uncertainties. Until the move towards convex optimization methods for robust time-optimal vibration control, existing methods utilized gradient-based numerical methods to determine the optimal control profile. Convex optimization methods for robust vibration control began with the formulation of the minimax control problem in a linear programming framework [19]. Some of these methods define new metrics, for example, using the L1 and L∞ norms to approximate the quadratic form of the system’s residual energy at the end of the maneuver. Another method approximates the residual energy hyperspheres (the true quadratic form) using hyperplanes represented in barycentric coordinates. Vossler [32] includes a different method for linear approximation of the system’s residual energy in his paper addressing deflection limited control. Conord and Singh [5] formulated the minimax problem as a convex optimization problem using Linear Matrix Inequalities (LMI) to determine a globally optimal minimax solution, the first method neither requiring approximation of the quadratic form of the system’s residual energy, nor parameterization of the control profile. This method exploits the powerful semi-definite programming solvers available today and is the first method available to accurately approximate the robust solution given by the traditional minimax time-delay approach using a convex framework. 1.2 Overview A new approach for approximating the traditional minimax time-delay solution is presented in this thesis. Chapter 2 reviews the latest methods for determining a ro- Chapter 1. Introduction 5 bust time-optimal solution for rest-to-rest maneuvers of flexible structures, starting with the traditional minimax time-delay approach proposed by Singh, as the success of this method set a new standard for robustness to modeling uncertainties. Subsequent methods proposed for approximating the traditional minimax solution are discussed and the motivation for the new formulation presented in this thesis is developed. Chapter 3 introduces the new formulation for approximating the traditional minimax solution and is tested on several benchmark problems in Chapter 4, used in comparison with the traditional minimax time-delay approach. Finally, conclusions and future contributions are discussed in Chapter 5. Chapter 2 Review of Techniques While several tools exist to incorporate parameter uncertainty in the design of control profiles to drive vibratory systems, several issues still exist when designing controllers for large-scale systems. The various useful methods for designing robust controllers are discussed in detail in order to provide the background information needed to form the motivation of this thesis. This chapter addresses the various techniques in chronological order, as each method is developed with the goal of improving upon the prior. The move from gradient-based optimization methods to convex methods is significant, and the resulting benefits will be highlighted. It should be noted that the development of the traditional minimax approach by Singh [18] marks a turning point in significantly improving robustness in design of time-optimal control profiles for vibratory systems. Succeeding methods exist to approximate this method and further enable its application to large-scale systems. 6 Chapter 2. Review of Techniques 2.1 7 Traditional Minimax Approach Prior to Singh’s work on minimax control in 2002, in order to achieve robustness in time-delay and saturating controllers, an extra set of constraints was added to the optimization scheme to force the sensitivity states (the derivative of the cost and constraints with respect to the unknown parameters) of the system to zero at the nominal value of the system parameters. For example, as shown in [20], for a system driven by impulses, after one impulse is applied a second impulse can be applied to cancel the vibration caused by the first impulse. It is then desirable to determine the amplitudes and corresponding times at which the impulses should be applied to the system in order to achieve minimum vibration at the end of the maneuver. With a reasonable estimate of the system’s natural frequency, ω, and damping ratio, ζ, an equation can be derived to characterize the percentage of residual vibration at the end of the maneuver. This equation represents the amount of vibration caused by a sequence of impulses, relative to the vibration caused by a single, unity-magnitude impulse and is given by: V (ω, ζ) = e−ζωtn C(ω, ζ) = n X p C(ω, ζ)2 + S(ω, ζ)2 Ai eζωti cos(ω p (2.1a) 1 − ζ 2) (2.1b) 1 − ζ 2) (2.1c) i=1 S(ω, ζ) = n X Ai eζωti sin(ω p i=1 where Ai and ti represent the impulse amplitudes and time locations, respectively, and n represents the total number of impulses. With additional constraints on the control vector (e.g. the control vector must be monotonically increasing and/or the impulses must sum to one), setting Eq, (2.1a) equal to zero will allow for determination of closed-form solutions for Ai and ti in terms of ω and ζ. In order to achieve robustness Chapter 2. Review of Techniques 8 in the vicinity of the nominal parameter values, inclusion of the constraint that forces the sensitivity states of the solution to zero d V (ω, ζ) = 0 dω (2.2) reduces the sensitivity to modeling errors near the nominal parameter values. While a desirable level of robustness can be achieved, methods that include constraints on the sensitivity states of the system consider the nominal values of the system parameters only, and are therefore limited to increasing robustness in the mere vicinity of the nominal parameter values. In 2002 Singh proposed a method called minimax time-delay control [18], which minimizes the maximum magnitude of the system’s residual energy over a known range of uncertainty for the uncertain parameters. Here, with the knowledge that the uncertain parameters lie within a specified range, a controller is designed with the worst model in mind. This method allows for the inclusion of knowledge regarding the range of uncertainty for the given parameters, thus improving the effective range of insensitivity of the control profile to modeling errors. This method is discussed in the following paragraphs. An asymptotically stable mechanical system undergoing rest-to-rest maneuvers, can be represented as M ÿ + C(p)ẏ + K(p)y = Du plb ≤ p ≤ pub (2.3a) (2.3b) where K and C vary with the vector of uncertain parameters, p, with upper bounds, pub , and lower bounds, plb . In minimax time-delay control, the objective is to minimize the maximum value of the residual energy represented by the pseudo-energy function: Chapter 2. Review of Techniques 9 min max F x p 1 1 1 F (tf ) = ẏ T M ẏ + (y − yf )T K(y − yf ) + (yr − yrf )2 2 2 2 (2.4) Eq. (2.4) is evaluated at final time tf , using the corresponding desired final displacement states, yf . When y=yf , the potential energy of a hypothetical spring is zero. M is positive-definite, and K is either positive-definite or positive-semidefinite, depending on whether or not rigid body modes exist. In the case where K is positivesemidefinite, yr is the rigid body displacement and yrf is the corresponding desired final displacement. When K is positive-definite, the final term in Eq. (2.4) can be neglected. In the traditional minimax time-delay control approach [18], the control vector u is parameterized as u = A0 + n X Ai H(t − Ti ) (2.5) i=1 where H(t − Ti ) is the Heaviside function. In the case where −1 ≤ u ≤ 1, the control vector is parameterized as u=1+ n X 2(−1)i H(t − Ti ) + H(t − Tn+1 ) (2.6) i=1 for a time-optimal control problem. For a system undergoing rest-to-rest maneuvers, with initial conditions set to zero, the states of the system represented by state space model y ż = Az + Bu, z = ∈ Rn , u ∈ R1 ẏ (2.7) can be solved for using the Van Loan identity, as shown in [30]. Matrix P is constructed in order to determine the response of a linear system to a unit step input. A B P = (2.8) 0 0 Chapter 2. Review of Techniques 10 P is a Rn+1×n+1 matrix. It can be shown using the Van Loan identity that the convolution integral of the system in Eq. (2.7) subject to a unit step input can be given by the upper right hand term of matrix M , shown in Eq. (2.9). R T A(T −τ ) AT e e Bdτ 0 M = eP T = 0 I (2.9) The response of the system shown in Eq. (2.7) to the input represented by Eq. (2.5) φ = A0 eP Tn + n X Ai eP (Tn −Ti ) (2.10) i=1 is given by the first n rows of the last column of matrix M . This method of finding the system response is used to calculate the pseudo-energy function, the representation of the system’s residual energy, needed to determine the minimax solution. 2.1.1 Example The traditional minimax approach is illustrated for the harmonic oscillator shown in Figure 2.1. The equation of motion is written as mÿ + cẏ + ki y = ki u, i = 1...p (2.11) with uncertainty in stiffness, k, and the following boundary conditions. y(0) = ẏ(0) = 0, y(tf ) = 1, ẏ(tf ) = 0 (2.12) For this problem m = 1, c = 0, and k is assumed to lie between 0.7 ≤ k ≤ 1.3 with a nominal value of k = 1. The uncertain parameter k, written as ki where i = 1...p, is discretized into p = 51 different values, evenly spaced over the range of uncertainty. In the traditional minimax time-delay control approach u is parameterized as u = A0 + A1 e−sT1 + A2 e−sT2 (2.13) Chapter 2. Review of Techniques 11 Figure 2.1: Harmonic Oscillator choosing n = 2 time-delays. The initial guess supplied to the gradient-based optimizer for the unknown control vector of switch times and amplitudes is: h i h i T1 T2 A0 A1 A2 = 3 6 0.25 0.5 0.25 (2.14) The additional constraint A0 + A1 + A2 = 1 (2.15) is included to guarantee that the amplitudes of the control vector sum to one. Expressing our vector of unknown amplitudes and switch times as h i X = T1 T2 A0 A1 A2 (2.16) the Matlab f minimax solver minimizes the maximum value of the pseudo-energy function, the representation of the residual energy, shown in Eq. (2.4), a set of p functions Fi (X), expressed in terms of the unknown vector of parameters, X. This set of functions must be supplied to the optimizer and is expressed as 1 1 Fi (X) = (y(tf ) − 1)ki (y(tf ) − 1) + (ẏ(tf ))m(ẏ(tf )), i = 1...p 2 2 (2.17) Chapter 2. Review of Techniques 12 where y(tf ) and ẏ(tf ) are given by the solution of the convolution integral shown in Eq. (2.10), which is calculated in terms of the vector of unknowns, X. Eq. (2.17) is the representation of Eq. (2.4) for the given problem. The solution to the minimax 2-time-delay control problem is 0.2571 + 0.4857e−3.1591s + 0.2571e−6.3182s (2.18) and is shown in Figure 2.2 along with the system response for the nominal stiffness value. The corresponding variation of the residual energy over the vector of uncertain 1 0.8 control 0.6 position 0.4 velocity 0.2 0 0 1 2 3 Time (sec) 4 5 6 Figure 2.2: Control Input and System Response for Spring Mass System parameters, ki is shown in Figure 2.3. Figure 2.3 reveals that the residual energy is small but never exactly zero. This explains the steady-state error for the system’s Chapter 2. Review of Techniques 13 −4 4.5 x 10 4 Residual Energy, F 3.5 3 2.5 2 1.5 1 0.5 0 0.7 0.8 0.9 1 1.1 Stiffness, k 1.2 1.3 1.4 Figure 2.3: Residual Energy Variation vs. Spring Stiffness response shown in Figure 2.2. Figure 2.4 shows the variation in position along with the variation in velocity for the various stiffness values that represent the range of parameter uncertainty. The response of the nominal system is highlighted in red. Since the vector of values for the uncertain parameter, ki , is assumed to be uniformly distributed, no special emphasis is placed on driving the residual energy of the nominal system to zero. In order to eliminate the steady-state error for the nominal system, an equality constraint can easily be included in the optimization scheme to force the residual energy for the nominal system to zero. Chapter 2. Review of Techniques 14 System Response with Various Stiffness Values 1 System States 0.8 Position Variation 0.6 Velocity Variation 0.4 0.2 0 −0.2 0 1 2 3 Time (sec) 4 5 6 Figure 2.4: System Response for Spring Mass System 2.2 Gradient-Based Optimization The traditional minimax time-delay approach requires the use of gradient-based optimizers provided by MATLAB for finding optimal control profiles. In both inputshaping/time-delay control and the traditional minimax approach, the solution to the optimization problem is formulated in terms of unknown parameters. Closed-form solutions can only be found for a small set of problems. The cost and constraints are typically represented by nonlinear equations in terms of the unknown parameters being optimized for. The parameter values are then solved for using solvers such as f mincon and f minimax in Matlab, which use the following techniques to find an Chapter 2. Review of Techniques 15 optimal solution. Both the f minimax solver and the f mincon solver convert the given problem into an equivalent Nonlinear Programming problem of the form M inimize f (x) (2.19a) Subject to gi (x) = 0, i = 1...ml (2.19b) gi (x) ≤ 0, i = ml ...m (2.19c) that is solved using a sequential quadratic programming (SQP) method. The SQP method is based on quadratic approximation of the Lagrangian Function L(x, λ) = f (x) + m X λi gi (x) (2.20) i=1 which is used to form a quadratic programming (QP) subproblem. The QP subproblem is stated as 1 T d Hk d + ∇f (xk )T d 2 (2.21a) Subject to ∇gi (xk )T d + gi (xk ) = 0, i = 1...ml (2.21b) M inimize ∇gi (xk )T d + gi (xk ) = 0, i = ml ...m (2.21c) where Hk is a positive definite approximation of the Hessian matrix of the Lagrangian Function at the k th iteration. Eq. (2.21b) and Eq. (2.21c) represent a linearized version of the constraints shown in Eq. (2.19b) and Eq. (2.19c). An active set method is used to solve the QP subproblem in Eq. (2.21). This first involves calculation of a feasible point (if one exists). Then a sequence of feasible points are generated iteratively until the solution to the QP subproblem is found. A new iterate is found in the SQP scheme as xk+1 = xk + ak dk (2.22) Chapter 2. Review of Techniques 16 using step length ak , as determined by a sufficient decrease in a merit function. The approximation of the Hessian matrix, Hk , is updated by the quasi-Newton BFGS method. The f minimax and f mincon solvers are gradient-based, as seen in (2.21), and rely on an approximation of the gradients of the cost and constraints. The accuracy and speed of the optimizer can be increased with the inclusion of analytical gradients. For time-delay control, bang-bang, and bang-off-bang control profiles, closed form equations can be derived to represent the gradients of the cost and constraints. When gradients are not supplied to the optimizer, finite difference is used to approximate the gradients of the cost and constraints with respect to the parameters being optimized. Most gradient-based solvers require that an initial guess be provided to the solver in order to designate a starting point from which to move across the design space. A poor initial guess may result in local solutions only. This presents a problem in that there is no guarantee of a globally optimal solution and provided the incentive to develop the minimax formulation in a convex optimization framework for the determination of a globally optimal solution. 2.3 Convex Optimization Methods After the formulation of the traditional minimax time-delay approach, several other methods arose for approximating the minimax solution in a convex framework. These methods require the solution of linear, quadratic and semidefinite programming problems, which are convex problems, and therefore guarantee a globally optimal solution. These methods are discussed in the following sections. Chapter 2. Review of Techniques 2.3.1 17 Linear Programming Approach A Linear Programming (LP) problem is a convex optimization problem of the form M inimize cT z (2.23a) Subject to Aeq z = beq (2.23b) Az ≤ b (2.23c) where z is the parameter being optimized for and Eq. (2.23a), Eq. (2.23b), and Eq. (2.23c) represent a linear cost, linear equality constraints, and linear inequality constraints, respectively. Most solvers use the Simplex Method by Dantzig [12] to find a globally optimal solution very efficiently. In the traditional minimax approach, the residual energy, represented by the pseudo-energy function in Eq. (2.4) is a quadratic function, which does not allow for formulation as a LP problem. Several methods have been devised by Singh [19] to linearly approximate the pseudo-energy function for use in a LP framework. If a minimax problem can be formulated as min max cTi z + di , i = 1...p z i Subject to Aeq z = beq Az ≤ b (2.24a) (2.24b) (2.24c) it can be written in the form of a LP problem as M inimize f (2.25a) Subject to Aeq z = beq (2.25b) cTi z + di ≤ f, i = 1...p (2.25c) Chapter 2. Review of Techniques 18 where f is a scalar parameter to be optimized for that marks the upper-boundary of the constraints shown in Eq. (2.25c). In order to write a minimax problem in the forms shown in Eq. (2.24) and Eq. (2.25), several methods, for example, using new metrics defined by an L1 norm and an L∞ norm, are used to approximate the residual energy function. Lp Norm Metric There are several ways to represent length of a vector, the most familiar being the Euclidean Norm, or the L2 norm, which is expressed as: kxk2 = p |x1 |2 + |x2 |2 + ... + |xn |2 (2.26) The Lp norm of vector x can represent different length metrics depending on the value of p. When p is a real number, and p ≥ 1, the Lp norm of a vector, x, is defined as: 1 kxkp = (|x1 |p + |x2 |p + ... + |xn |p ) p (2.27) The L1 and L2 norms are written as kxk1 = |x1 | + |x2 | p kxk2 = |x1 |2 + |x2 |2 (2.28a) (2.28b) for a two-dimensional space. For p = ∞, the L∞ norm is given by kxkp = max{|x1 |, |x2 |, ..., |xn |} (2.29) for an n-dimensional space. The space defined by these norms is shown in Figure 2.5 for a two-dimensional vector where the norm is required to be less than the scalar quantity, f . To use the L1 norm to approximate the pseudo-energy function and formulate a LP problem of the form shown in Eq. (2.23c), the following constraints Chapter 2. Review of Techniques 19 Figure 2.5: Lp Norms in 2-Dimensional Space represent the boundary of the L1 norm space are defined, for example, for a twodimensional problem as: x1 + x2 ≤ f (2.30a) x1 − x2 ≤ f (2.30b) − x1 + x2 ≤ f (2.30c) − x1 − x2 ≤ f (2.30d) The corresponding boundary constraints for the L∞ norm are: − f ≤ x1 ≤ f (2.31a) − f ≤ x2 ≤ f (2.31b) Instead of minimizing the maximum of the quadratic representation of the system’s residual energy, a new metric is defined for minimization of the maximum of the L1 Chapter 2. Review of Techniques 20 norm and L∞ norm of the residual states over the domain of uncertainty. The residual states of the system at final time index, N , are defined as √ y y(N ) = Kx(N ) √ = ˙ ) = M ẋ(N ) ẏ y(N (2.32) which weights the position vector by the square root of the stiffness matrix and the velocity vector by the square root of the mass matrix. This results in equal contribution to the residual energy cost by the individual displacment and velocity states. For an n-dimensional problem, the minimax problem, which minimizes the L1 norm of the residual states over the domain of uncertain models can be stated in an LP form as: M inimize f #i=1,2,...,p " n X ckj yki (N ) ≤ f Subject to k=1 (2.33a) (2.33b) j=1,2,...,2n umin ≤ u ≤ umax , m = 1, 2...N (2.33c) The vertices of an n-dimensional hypercube centered at the origin, with a side length of √ 2 define the coefficients, ckj ± 1 for the given problem. The control is parameterized into N samples and p is the number of uncertain models over which the L1 norm is minimized. Now the model shown in Eq. (2.33) fits the form of the LP shown in Eq. (2.25). Chapter 2. Review of Techniques 21 The L∞ norm can be used to form a minimax problem of the form M inimize f (2.34a) Subject to − f ≤ y1i (N ) ≤ f (2.34b) − f ≤ y2i (N ) ≤ f (2.34c) .. . − f ≤ yni (N ) ≤ f − f ≤ ẏ1i (N ) ≤ f (2.34d) ∀i = 1, 2, ...p − f ≤ ẏ2i (N ) ≤ f (2.34e) (2.34f) .. . − f ≤ ẏni (N ) ≤ f umin ≤ u ≤ umax , ∀m = 1, 2...N − 1 (2.34g) (2.34h) which also fits the form shown in Eq. (2.25). Solving the problems shown in Eq. (2.33) and Eq. (2.34) require specification of the final maneuver time. Plots of the L1 norm and L∞ norm of the residual states defined in Eq. (2.32) versus the final maneuver time are monotonically decreasing, as seen in Figure 2.8 in the following section. This plot helps the designer identify the minimum maneuver time at which the norms fall below a specified value. Example Results are generated for the L1 norm LP approach and the L∞ norm LP approach for approximating the traditional minimax solution for the harmonic oscillator in Eq. (2.11) with the boundary conditions shown in Eq. (2.12). The uncertain spring stiffness is again discretized into p = 51 different values, evenly spaced over the range Chapter 2. Review of Techniques 22 of uncertainty, 0.7 ≤ k ≤ 1.3. The cost function for the L∞ norm approach is ° ° ° ° i ° y (tf ) − 1 ° ° min f = max ° ° i ° ẏ i (tf ) ° ° (2.35) ∞ which is included in the L∞ norm LP problem: M inimize f (2.36a) Subject to ẏ i (N ) ≤ f (2.36b) −ẏ i (N ) ≤ f (2.36c) √ i i k y (N ) ≤ f √ i − k y i (N ) ≤ f −1 y i (N ) y(1) N X N −1 −m = Gi + GN Hi u(m), ∀i = 1...p = 51 i ẏ i (N ) ẏ(1) m=1 (2.36d) (2.36e) (2.36f) 0 ≤ u(m) ≤ 1, ∀m = 1, 2, ..., N = 127 (2.36g) u(m) − u(m + 1) ≤ 0, ∀m = 1, 2, ...128 (2.36h) The state constraints are include in the optimization scheme with Eq. (2.36f), which is derived from the discrete time form of the system y(k + 1) = Gy(k) + Hu(K) (2.37) by writing the terminal states as a function of the prior states and the parameterized control vector, u, as yi (N ) −1 y(1) N X −1 + GN −m Hi u(m) = GN i ẏ(1) m=1 i ẏi (N ) (2.38) for the ith plant in the range of parameter uncertainty. The constraints shown in Eq. (2.36g) and Eq. (2.36h) require that the control vector is monotonically increasing and falls in the range 0 ≤ u ≤ 1. Chapter 2. Review of Techniques 23 The L1 norm LP problem for the harmonic oscillator is: M inimize f √ i i k y (n) ≤ f √ i ẏ i (N ) − k y i (N ) ≤ f √ i −ẏ i (N ) + k y i (N ) ≤ f √ i −ẏ i (N ) − k y i (N ) ≤ f −1 y i (N ) y(1) N X −1 −m = GN + GN Hi u(m), ∀i = 1...p = 51 i i ẏ(N ) ẏ(1) m=1 Subject to ẏ i (N ) + (2.39a) (2.39b) (2.39c) (2.39d) (2.39e) (2.39f) 0 ≤ u(m) ≤ 1, ∀m = 1, 2, ..., N = 128 (2.39g) u(m) − u(m + 1) ≤ 0, ∀m = 1, 2, ...127 (2.39h) with the same state constraints and control vector requirements given in the L∞ norm LP problem. The control vectors solved for using the L1 norm LP approach and the L∞ norm LP approach are plotted along with the control vector determined using the traditional minimax time-delay approach in Figure 2.6. Figure 2.8 shows how the residual state metrics vary with the final maneuver time, a plot that can assist the designer in finding the minimum maneuver time at which the metrics fall below a desirable threshold. The corresponding variations of the residual state metrics over the vector of uncertain parameters, ki are shown in Figure 2.7. As expected, both the L1 norm LP approach and the L∞ norm LP approach fail to achieve the level of robustness determined using the traditional minimax time-delay approach. This makes sense as we are approximating an L2 norm space with nonlinear constraint boundaries, using the L1 and L∞ norms, which both have linear constraint boundaries, as seen in Figure 2.5. The approximation error is significant, as Singh [19] shows that the maximum Chapter 2. Review of Techniques 24 1 0.9 Time−delay 0.8 Control Input, U L1 0.7 Linf 0.6 0.5 0.4 0.3 0.2 0 1 2 3 Time (sec) 4 5 6 Figure 2.6: Control Input for Spring Mass System magnitude of the residual state errors for the L1 minimax solution is greater than that of the L2 minimax solution. Other methods for posing the minimax problem in a LP framework have been proposed that use approximation methods to convert the quadratic residual energy constraints to linear constraints. Singh [19] uses a method that approximates the n-dimensional hypersphere, which represents the residual energy of an n-dimensional problem, using an (n − 1)-dimensional simplex (hyperplanes) written in barycentric coordinates. Another method, proposed by Singh and Vossler [32], requires a different linear approximation of the residual energy for use in a LP problem. Since large-scale Chapter 2. Review of Techniques 25 −3 3 x 10 Time−delay 2.5 Residual Energy, J L1 Linf 2 1.5 1 0.5 0 0.7 0.8 0.9 1 1.1 Stiffness, k 1.2 1.3 1.4 Figure 2.7: Residual Energy Variation vs. Spring Stiffness Value problems can be easily solved using the LP framework and the Simplex algorithm, these methods do fairly well. However, they do not adequately estimate the traditional minimax time-delay results, especially for large-scale systems. Over the past couple of years powerful software for solving large-scale semidefinite programming (SDP) problems relatively efficiently have been developed, such as the Matlab toolbox, YALMIP [10], SeDuMi [28], and SDPT3 [29]. Additionally, semi-definite programming problems, which employ the Linear Matrix Inequality (LMI), are convex optimization problems, and therefore can guarantee a globally optimal solution. This provided the motivation to pose the minimax problem in a Chapter 2. Review of Techniques 26 1.4 Residual Energy Metric, f 1.2 1 L1 Linf 0.8 0.6 0.4 0.2 0 2 3 4 5 6 7 Final Maneuver Time, tf 8 9 10 Figure 2.8: Residual Energy Metric vs. Final Maneuver Time LMI framework. The LMI approach is discussed in the following section. 2.3.2 LMI Approach The LMI approach proposed by Conord and Singh [5], was the first method able to closely approximate the minimax time-delay solution without parameterization of the control profile nor corresponding time-vector. For example, the designer is not required to specify the shape of the control vector or the corresponding number of time-delays. The control vector is given in discrete form, however, as the specified delay time, Ts , for the given problem approaches zero, the solution tends to the Chapter 2. Review of Techniques 27 traditional minimax time-delay solution, and is discussed in the following. A semi-definite programming (SDP) problem is of the form: M inimize cT y Subject to V0 + (2.40a) n X yi Vi º 0 (2.40b) i=1 where y ∈ Rn is the vector of variables being optimized and Vi are symmetric matrices. Eq. (2.40b) represents a Linear Matrix Inequality (LMI), where the symbol º represents a positive semi-definite constraint. The minimax problem is now posed in the framework shown in Eq. (2.40). For the system in Eq. (2.7), the control vector, u, is parameterized as 1 u(s) = (A0 + A1 e−sT1 + A2 e−sT2 + ... + AN e−sTN ) s (2.41) where Aj represent amplitudes and Tj represent time delays, which are assumed to be known and equal to integer multiples of the specified delay time, Ts . For an uncertain system with i = 1...p, different plants discretized over the range of uncertainty, the minimax problem is now represented as an optimization problem in terms of the unknown amplitudes Aj . The pseudo-energy function is written as min max F = ∆xTi Qi ∆xi , i = 1...p i (2.42) where ∆xi = xi (TN ) − xf represents the error between the final states and the desired final states and is written in terms of the unknown control vector amplitudes, Aj . Qi is a positive definite weighting matrix, defined by the ith stiffness and mass matrices of the system. First, it is shown how to write Eq. (2.42) in the form of an LMI shown in Eq. (2.40). Then it is shown how to write ∆xi in terms of the unknown control vector amplitudes, Aj . Chapter 2. Review of Techniques 28 The minimax problem is represented as a minimization problem using the scalar variable, γ, which represents the upper boundary of the residual energy constraint. The problem becomes M inimize F = γ (2.43a) Subject to γ − ∆xTi Qi ∆xi ≥ 0, i = 1...p (2.43b) Alb ≤ Aj ≤ Aub (2.43c) where Alb and Aub represent the lower and upper bounds, respectively, on the control vector amplitude. Using the Schur complement [33], the optimization problem in Eq. (2.43) is rewritten as M inimize F =γ T γ ∆xi Â 0, i = 1...p Subject to ∆xi Q−1 i Alb ≤ Aj ≤ Aub (2.44a) (2.44b) (2.44c) and satisfies the form of Eq. (2.40), as long as Aj are affine in ∆xi . Next, ∆xi is written in terms of Aj . In order to do so, the error states are written as: Z TN Ai TN eAi (TN −t) Bi u(t)dt − xf ∆xi = e x(0) + (2.45) 0 New variables, Ãj , are introduced, which represent the effective value of the input on each interval, such that u(τ ) = Ãj = A0 + ... + Aj (2.46) where τ ∈ [j × Ts ; (j + 1) × Ts ] for all j between 0 and N , for x(0) = 0. Now the terminal constraint on the control input can be represented as ÃN = 1. Eq. (2.45) can be rewritten as: ∆xi = N −1 X k=0 ÃZ (k+1)×Ts k×Ts ! eAi (TN −τ ) Bi Ãk dτ − xf (2.47) Chapter 2. Review of Techniques 29 Eq. (2.47) can be written as ∆xi = N −1 µZ Ts X ¶ e Ai ((N −k)×Ts −τ 0 ) Bi dτ 0 − xf (2.48) 0 k=0 if the variable τ is changed to τ 0 = τ − k × Ts , since Tj = j × Ts for all j between 0 and N . Since Ãk is a scalar, we can write ÃN −1 ! µZ X ∆xi = Ãk eAi (N −K−1)×Ts Ts ¶ Ai (Ts −τ ) e Bi dτ (2.49) 0 k=0 and it is clear that the input magnitudes appear affinely in ∆xi . The expression for ∆xi is written in matrix form as: ∆xi = Ii Ui Hi Ei − xf where, where Ii Ui Hi = Ii = [I ... I], Ã I 0 0 0 ... Ui = 0 0 0 0 ÃN −1 I (eAi Ts )N −1 .. . Hi = eAi Ts I Z Ts Ei = eAi (Ts −τ ) Bi dτ PN −1 k=0 (2.50a) (2.50b) (2.50c) (2.50d) (2.50e) 0 Ãk eAi (N −k−1)×Ts , is the summation shown in Eq. (2.49). The convolution integral, Ei , also shown in Eq. (2.49), is computed using the Van Loan Identity method discussed in Section 2.1. Matrices Ii , Ui , and Hi , are easily computed from the n-dimensional state-space model of the system, and are of dimension n × (nN ), (nN ) × (nN ), and (nN ) × n, respectively. In order to solve the SDP problem stated in Eq. (2.43), a final maneuver time must be specified. To determine the minimum time for which the residual energy Chapter 2. Review of Techniques 30 given by γ falls below an allowable threshold, the SDP is solved iteratively using the Bisection Algorithm to decrease an initial conservative guess for the final maneuver time until the allowable threshold is met. The residual energy variation with the final maneuver time is shown in Figure 2.9 for the example shown in the next section. The optimal final maneuver time determined using the traditional minimax time-delay approach is indicated by the ∗. 0.7 0.6 Residual Energy, F 0.5 0.4 0.3 0.2 0.1 0 −0.1 0 2 4 6 Final Maneuver Time, tf 8 Figure 2.9: Control Input for Spring Mass System 10 Chapter 2. Review of Techniques 31 Example The LMI formulation for approximating the traditional minimax time-delay solution is solved for the harmonic oscillator with uncertainty in stiffness shown in Eq. (2.11). Here, matrix Qi is ki 0 Qi = 0 m (2.51) where ki represents the ith manifestation of the uncertain spring stiffness. The uncertain stiffness range of 0.7 ≤ k ≤ 1.3 is evenly discretized into p = 11 values. With the control vector discretized into N = 128 different values, the following results are obtained using a final maneuver time of TN = 6.3182 to match the solution using the traditional minimax time-delay approach. Figures 2.10 and 2.11 show that the LMI formulation closely approximates the solution determined using the traditional minimax time-delay approach. LMI Solvers While much work has been done on the development of solvers for SDPs, limitations still exist when trying to find solutions to large-scale optimization problems efficiently. For example, the Matlab LMI toolbox is used to solve the minimax problem posed in the SDP framework. While there is no built-in limitation on the number of LMIs that can be specified or on the number of blocks and terms in any given LMI when using the Matlab solvers, efficiency and complexity issues arise, especially when dealing with large-scale problems where the total number of rows in the LMI system is greater than or equal to the number of decision variables. When incorporating parameter uncertainty in the LMI framework, the case when posing the minimax problem as an SDP, the number of LMI constraints quickly increases as multiple Chapter 2. Review of Techniques 32 1 Traditional Minimax LMI 0.9 Control Input, U 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0 1 2 3 4 Time (sec) 5 6 7 Figure 2.10: Control Input for Spring Mass System plants are included in the model to represent the entire specified range of parameter uncertainty, and therefore, the number of rows in the LMI system quickly exceeds the number of decision variables (defined by the level of discretization of the control profile). Matlab LMI solution methods incorporate a least-squares problem that is either solved via Cholesly factorization of the Hessian matrix or by QR factorization when numerical instabilities occur in a given problem. Memory problems quickly crop-up and efficiency decreases when using these methods to solve complex-large scale problems. These issues provide incentive for utilization of the most efficient solvers available for solving large-scale optimization problems, part of the motivation for the new formulation proposed in this thesis. Chapter 2. Review of Techniques 33 −4 x 10 Traditional Minimax LMI 4.5 4 Residual Energy, F 3.5 3 2.5 2 1.5 1 0.5 0 0.7 0.8 0.9 1 1.1 Stiffness, k 1.2 1.3 1.4 Figure 2.11: Residual Energy Variation vs. Spring Stiffness Value 2.4 Motivation Since the development of the traditional minimax approach in 2002, subsequent methods have been developed to approximate the minimax solution and further enable its application to large-scale systems. The traditional minimax approach relies on a gradient-based optimizer, so application to large-scale systems is difficult, as the method relies on a good initial guess and may result in local solutions. In addition, the traditional method requires the designer to specify the shape of the parameterization of the control profile and corresponding time vector. A move towards convex optimization methods took place with the development Chapter 2. Review of Techniques 34 of similar methods posed in the Linear Programming framework, which enabled easy application to large-scale systems. While these methods guarantee a global solution and do not require specific parameterization of the control profile, the error associated with linearly approximating the quadratic form of the residual energy is significant, yielding a solution that is less robust to parameter uncertainties than the traditional minimax solution. In 2008, the LMI approach for accurately approximating the traditional minimax solution was developed, which utilizes the powerful SDP solvers available today. The LMI approach guarantees a globally optimal solution without requiring specific parameterization of the control profile and corresponding time vector. The success of the LMI approach fueled the motivation for exploitation of the most effective solvers for large-scale optimization problems. As a result, the new approach discussed in this thesis was developed, which utilizes the powerful software available for solving optimization problems in the Quadratically Constrained Programming (QCP) framework. While powerful solvers exist for determining the minimax solution in the LMI framework, much more work has been done in developing QCP solvers commercially available for solving large-scale problems. This thesis will introduce the QCP approach for approximating the minimax solution and its potential for use with large-scale systems. Chapter 3 QCP Minimax Formulation A method is proposed to directly incorporate the quadratic form of the residual energy along with the system’s state constraints into the optimization scheme to determine a control profile that minimizes the maximum magnitude of the system’s residual energy, and is robust in the range of parameter uncertainty. This method does not require the user to define the shape of the parameterization of the control vector and will be shown to closely approximate the traditional minimax time-delay solution. The proposed approach will allow for the use of the most advanced tools capable of solving large-scale optimization problems relatively efficiently. 3.1 Residual Energy Derivation A new expression for the system’s residual energy is derived, to allow for the proposition of the minimax problem in a Quadratically Constrained Programming (QCP) framework. The new expression combines the system’s residual energy constraint with the system’s state constraints to form one constraint that is in the quadratic 35 Chapter 3. QCP Minimax Formulation 36 form. Consider a mechanical system whose second order vector differential equation is M ÿ + C ẏ + Ky = Du (3.1) where M , C and K represent the mass, damping and stiffness matrices respectively. The first order representation of this system is ẏ 0 I y 0 u = + ÿ −M −1 K −M −1 C ẏ M −1 D (3.2) which can be written in discrete time as: k x(k + 1) = G x0 + k X Gk−i Hu(i) (3.3) i=1 G and H are the discrete time state space matrices of the system. u(i) is the discrete time control vector. The states at final time, N + 1, can be written as: x(N + 1) = GN x0 + N X GN −i Hu(i) (3.4) i=1 The system maneuver time is broken up into N intervals, with sampling rate, Ts , and final maneuver time, tf = (N + 1)Ts . We can now write the system’s residual energy, J, as: K 0 {xf − x(N + 1)} 2J = {xf − x(N + 1)}T 0 M (3.5) = {xf − x(N + 1)}T Em {xf − x(N + 1)} This is a different form of the pseudo-energy function shown in Eq. (2.4). Here Em is the energy matrix and xf is the desired final value of the state vector. Combining Chapter 3. QCP Minimax Formulation 37 Eq. (3.4) and Eq. (3.5), we then obtain: 2J = uT Qu + F T u + L (3.6a) L = xTf Em xf + (GN x0 )T Em GN x0 − 2xTf Em GN x0 h i F T = 2(GN x0 − xf )T Em GN −1 H GN −2 H ... H (3.6b) where (GN −1 H)T Em GN −1 H (GN −1 H)T Em GN −2 H ... (GN −1 H)T Em H (3.6c) (GN −2 H)T Em GN −1 H (GN −2 H)T Em GN −2 H ... (GN −2 H)T Em H Q= : : : H T Em (GN −1 H) H T Em GN −2 H ... H T Em H (3.6d) An equation for the system’s residual energy has been written in terms of known quantities and the unknown control parameter, u. Eq. (3.6) now represents both the system’s residual energy and the state constraints and has combined them in a quadratic form. An optimization problem can be formulated in the quadratically constrained programming (QCP) framework. 3.2 QCP Framework The QCP problem must meet the criteria for convexity in order to guarantee a global solution and improve upon the gradient-based methods. According to [3] a convex optimization problem in standard form is represented as M inimize f0 (x) (3.7a) Subject to fi (x) ≤ 0, i = 1...m (3.7b) hi (x) = 0, j = 1...l (3.7c) where f0 ...fm are convex functions. In order to have a convex optimization problem, a convex objective function must be minimized over a convex set. Thus, the inequality Chapter 3. QCP Minimax Formulation 38 constraint functions fi (x), i = 1...m comprise a convex set. In addition, the equality constraint functions hi (x) = aTi x − bi must be affine. The feasible set of a convex optimization problem must be convex, since both the feasible set is the intersection of convex sets that cover the domain of the problem, and the intersection of convex sets is also convex. The domain of the problem is a convex set with m + 1 convex functions consisting of the objective function and the inequality constraint functions and l hyperplanes aTi x − bi = 0, i = 1...l. As stated in [6], a convex QCP problem is defined as M inimize cT x (3.8a) Subject to xT Qi x + 2qiT x + γi ≤ 0, i = 1...b (3.8b) Aj x ≤ bj , j = 1...w (3.8c) Ck x = dk , k = 1...z (3.8d) xlb ≤ x ≤ xub (3.8e) where the decision vector x ∈ Rn , and the parameters Qi º 0 (i.e. Qi is positive semidefinite), Qi ∈ Rnxn , c ∈ Rn , γi ∈ R, and qi ∈ Rn , for all i = 1...p. Eq. (3.8b) is a set of b convex quadratic constraints, Eq. (3.8c) and Eq. (3.8e) compose a convex polyhedron, and Eq. (3.8d), z hyperplanes. The expression shown in Eq. (3.8) fits the form of a convex problem stated in Eq. (3.7). The residual energy expression shown in Eq. (3.6) is now posed in the convex QCP framework, Eq. (3.8), to determine an approximation to the minimax time-delay control solution. In the case of an asymptotically stable mechanical system undergoing restto-rest maneuvers, which varies over a range of uncertain parameters, as shown in Eq. (2.3), the corresponding representation of the pseudo-energy function shown in Chapter 3. QCP Minimax Formulation 39 Eq. (3.6) can be written as J = uT Q(p)u + F T (p)u + L(p) (3.9) since Q, F , and L now vary with the vector of uncertain parameters, p. When we discretize the range of each uncertain parameter to cover the full-range of uncertainty and determine Q, F , and L for each and every combination of uncertain parameters, this results in m different system models. We now have a representation of the pseudoenergy function for a range of system models, which covers the full range of parameter uncertainty for a given problem. The minimax time-delay control solution for u, which best satisfies the full range of parameter uncertainty can now be determined by the following optimization problem: M inimize z (3.10a) Subject to z ≥ uT Qj u + FjT u + Lj , ∀j = 1..m (3.10b) where z is a parameter to be optimized along with vector u and represents the upperbound of the quadratically constrained region, the way in which the minimax problem is posed as a minimization problem. Since Q, F , and L are derived from the stiffness and mass matrices, K and M for the system model, along with the discrete time state matrices, G and H, it can be shown that Qj º 0 and the constraints shown in Eq. (3.10b) satisfy the form of the QCP constraints shown in Eq. (3.8b). Eq. (3.10) is a convex QCP program that represents an approximation to the minimax time-delay control solution, when solved iteratively for the optimal final time tf . 3.2.1 Minimum Time Solution The solution to the QCP program shown in Eq. (3.10) represents the minimax solution for the set of system models generated for a given final maneuver time tf . In order Chapter 3. QCP Minimax Formulation 40 Figure 3.1: Reduction of the Final Maneuver Time using the Bisection Algorithm to determine the time-optimal solution that approximates the minimax time-delay approach, the Bisection Algorithm is used to determine the minimum final maneuver time, which satisfies the specified allowable residual energy criteria. Residual energy decreases monotonically, as shown in Figure 3.1. With a conservative initial estimate, a zero-order optimization scheme, such as Golden Section or Bisection Algorithm, can be used to decrease the final maneuver time accordingly until the threshold of allowable residual energy is crossed. Figure 3.1 shows how the initial conservative guess for the final maneuver time, shown at tf 1 evolves according to the Bisection algorithm. The optimal solution is shown with the ?, which marks the minimum time at which the residual energy, J, Chapter 3. QCP Minimax Formulation 41 falls below the allowable threshold of the residual energy, Jallow , at the end of the maneuver. The Bisection algorithm simply involves evaluation of the residual energy function, J, at the beginning, middle and end of each time interval. The initial time interval, [lb1 , ub1 ], is specified by the designer. The middle of the current interval is found, and new lower and upper bounds are chosen as two of either the current upper bound, current lower bound, and the middle of the interval, based on evaluation of the residual energy at the current upper bound, current lower bound, and the middle of the interval. New interval bounds are assigned, [lb2 , ub2 ], such that the residual energy at the lower bound falls above the allowable residual energy, Jallow , and the residual energy at the upper bound falls below the allowable residual energy, Jallow . This guarantees that the optimal point remains inside the current time interval at all times. The designer specifies a tolerance on the size of the time interval that brackets the optimal point to signal adequate convergence to the optimal final maneuver time. 3.2.2 Preprocessing in Matlab and Solving with GAMs Derivation of Q(j), F (j), and L(j) for each system studied is easily done in Matlab. The optimization problem shown in Eq. (3.10) is formulated and solved using the CPLEX solver in GAMs. Data is passed back and forth between Matlab and GAMs using the ”Matlab and GAMS: Interfacing Optimization and Visualization Software” by Michael C. Ferris (University of Wisconsin). The QCP can be solved using the GAMS/CPLEX solver, which is based on the Cplex Callable Library from the ILOG CPLEX Division. For solving QCPs, CPLEX implements the Cplex Barrier method, which transforms the QCPs into Second Order Cone Programs (SOCPs) and outputs a primal solution. Chapter 3. QCP Minimax Formulation 3.2.3 42 Flow of Information Figure 3.2 shows how information is passed between Matlab and GAMs when solving the QCP minimax problem. Figure 3.2: QCP Minimax Approach Chapter 3. QCP Minimax Formulation 3.3 43 Conclusion A set of quadratic constraints has been generated that incorporate the quadratic form of the system’s residual energy, along with the system’s state constraints. These constraints are used in a convex QCP framework to determine a global solution to the minimax problem. This method utilizes the powerful large-scale optimization solvers available today, and accurately approximates the minimax time-delay solution. Numerical results are shown in the next section for comparison with the traditional minimax approach. Chapter 4 Numerical Examples In this section the QCP minimax approach is compared to the traditional minimax time-delay approach discussed in Section 2.1. The techniques are first tested on two benchmark problems: a harmonic oscillator and a floating oscillator. Finally, the QCP minimax approach is applied to a simple elevator model for determination of efficacy in application to large-scale systems. 4.1 Harmonic Oscillator The first problem in this section addresses uncertainty in a single parameter, in this case, the stiffness in a harmonic oscillator model. Next, the QCP approach is tested on a harmonic oscillator with uncertainty in both stiffness and damping. 4.1.1 Uncertain Stiffness The first problem addressed is that of uncertainty in stiffness in the harmonic oscillator shown in Eq. (2.11), with the boundary conditions shown in Eq. (2.12). Again, 44 Chapter 4. Numerical Examples 45 we choose m = 1, c = 0, and k is assumed to be the only uncertain parameter. With a nominal value of k = 1, k is assumed to lie between 0.7 ≤ k ≤ 1.3. The uncertain parameter k is discretized into p = 51 different values, evenly spaced over the range of uncertainty. In the traditional minimax time-delay control approach u is parameterized as u = A0 + A1 e−sT1 + A2 e−sT2 (4.1) The initial guess provided to the gradient-based optimizer for the unknown control vector switch times and amplitudes is: h i h i T1 T2 A0 A1 A2 = 3 6 0.25 0.5 0.25 (4.2) The solution is determined to be 0.2571 + 0.4857e−3.1591s + 0.2571e−6.3182s (4.3) for the minimax 2-time-delay control problem. In the QCP approach, matrices Qj , Fj , and Lj , for j = 1..51, are derived to cover the full range of parameter uncertainty. The solution to the minimax QCP formulation for the harmonic oscillator requires specification of the final maneuver time. In this case, the optimal final maneuver time is taken from the result for the 2time-delay minimax controller in Eq. (4.3), to allow for accurate comparison without the use of the Bisection Algorithm. Upper and lower bounds on the control vector, ul and uu respectively, are specified as 0 = ul ≤ u ≤ uu = 1 (4.4) uj+1 − uj ≤ 0, ∀j = 1..(N − 1) (4.5) The additional constraint Chapter 4. Numerical Examples 46 is added to the QCP problem to guarantee that u is monotonically increasing, a requirement necessary to guarantee that the QCP solution meets the same shape requirements given by the parameterization of the 2-time-delay minimax control profile. Results for both the 2-time-delay minimax controller and the QCP minimax controller are shown in Figures 4.1-4.4. The slight differences in position and velocity System Response with Nominal Stiffness Value 1 QCP Position 0.8 QCP Velocity System States Time−Delay Position Time−Delay Velocity 0.6 0.4 0.2 0 −0.2 0 1 2 3 Time (sec) 4 5 6 Figure 4.1: System Response for comparison of the 2-Time-Delay Minimax solution and the QCP Minimax Solution of the nominal system, when comparing approaches in Figure 4.1, can be explained by the minor differences in the optimal control profiles determined using the two methods. Chapter 4. Numerical Examples 47 1 QCP 0.9 TD Control Input, U 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0 1 2 3 4 Time (sec) 5 6 7 Figure 4.2: Control Input vs. Time for comparison of the 2-Time-Delay Minimax solution and the QCP Minimax Solution The control profiles differ, as shown in Figure 4.2, due to the discretization error in the QCP formulation. Since the traditional minimax time-delay approach determines the exact instances of the optimal time-delays, the discrete representation of the control profile used in the QCP formulation may not allow for a change in input amplitude at the exact optimal time. As the level of discretization increases, the QCP solution will tend towards the traditional minimax time-delay solution. Chapter 4. Numerical Examples 48 −4 4.5 x 10 4 Residual Energy, J 3.5 3 2.5 2 1.5 QCP Time−Delay 1 0.5 0 0.7 0.8 0.9 1 1.1 Stiffness, k 1.2 1.3 1.4 Figure 4.3: Residual Energy vs. Stiffness for comparison of the 2-Time-Delay Minimax solution and the QCP Minimax Solution Figure 4.3 shows how closely the QCP minimax approach approximates the solution to the traditional minimax time-delay approach. With the control vector discretized into 128 different evenly spaced values, excellent results are achieved. Chapter 4. Numerical Examples 49 0.7 0.6 Residual Energy, J 0.5 0.4 0.3 0.2 0.1 0 −0.1 0 2 4 6 Final Maneuver Time, tf 8 10 Figure 4.4: Residual Energy vs. Final Maneuver Time for the QCP Minimax Approach The variation in the residual energy cost with the final maneuver time is shown to be monotonically decreasing in Figure 4.4. This allows for the use of the Bisection Algorithm to determine the optimal final maneuver time. 4.1.2 Uncertain Stiffness and Damping The same approach is carried out for the harmonic oscillator shown in 2.1, with uncertainty in both stiffness and damping. The range of uncertainty for stiffness remains the same, 0.7 ≤ k ≤ 1.3 with a nominal value of k = 1. The value of Chapter 4. Numerical Examples 50 damping now falls in the range of uncertainty, 0.1 ≤ c ≤ 0.3, with a nominal value c = 0.2. k and c are discretized into 15 different values each, yielding m = 225 different system models included in the optimization scheme. The initial guess provided to the gradient-based optimizer used in the 2-time-delay minimax control problem is h i h i T1 T2 A0 A1 A2 = 3 6 0.3 0.4 0.2 (4.6) for the parameterization shown in Eq. (4.1). Constraints on the control vector shown in Eq. (4.4) and Eq. (4.5) are included in the formulation. The solution for the minimax 2-time-delay control profile with uncertainty in both stiffness and damping is u = 0.3360 + 0.4739e−3.1604s + 0.1901e−6.3296s Results are shown in Figures 4.5-4.8. (4.7) Chapter 4. Numerical Examples 51 System Response with Nominal Stiffness and Damping Values 1.2 1 QCP Position QCP Velocity Time−Delay Position Time−Delay Velocity System States 0.8 0.6 0.4 0.2 0 −0.2 0 1 2 3 Time (sec) 4 5 6 Figure 4.5: System Response for comparison of the 2-Time-Delay Minimax solution and the QCP Minimax Solution In Figure 4.5 the QCP approach is shown to very closely approximate the results obtained using the traditional minimax time-delay approach. The steady-state error shown is expected, as the system response is generated for the nominal parameter values. A constraint was not included to force the residual energy to zero at the nominal parameter values, although an additional constraint could be easily included to do so. Chapter 4. Numerical Examples 52 1 Time−Delay QCP 0.9 Control Input, U 0.8 0.7 0.6 0.5 0.4 0 1 2 3 4 Time (sec) 5 6 7 Figure 4.6: Control Input vs. Time for comparison of the 2-Time-Delay Minimax solution and the QCP Minimax Solution A very small error exists between the control vector generated using the QCP approach and the time-delay approach, as seen in Figure 4.6. This is again due to the discretization of the control vector used in the QCP approach, which was chosen as N = 128. The error can be further reduced with a finer discretization. Square Root of Residual Energy Chapter 4. Numerical Examples 53 0.025 0.02 0.015 0.01 0.005 0 0.4 1.4 0.3 1.2 1 0.2 Damping, c 0.1 0.8 Stiffness, k Figure 4.7: Square Root of Residual Energy vs. Damping vs. Stiffness for the QCP Minimax Solution Figure 4.7 shows the 3-dimensional variation of the square root of the system’s residual energy with the values of damping and stiffness. Chapter 4. Numerical Examples 54 0.6 Residual Energy, J 0.5 0.4 0.3 0.2 0.1 0 −0.1 0 1 2 3 4 5 Final Maneuver Time, tf 6 7 8 Figure 4.8: Residual Energy vs. Final Maneuver Time for the QCP Minimax Solution The variation of the residual energy with final maneuver time, shown in Figure 4.8, is monotonically decreasing as the final maneuver time is increased. This confirms that the Bisection Algorithm can be used to find the minimum time at which the allowable residual energy threshold is met for the case of uncertainty in multiple parameters. 4.2 Floating Oscillator In this section, the QCP minimax approach and the traditional minimax time-delay approach are applied to the benchmark problem given by the floating oscillator shown in Figure 4.9 below. The approaches are first applied to the system with uncertainty Chapter 4. Numerical Examples 55 in the stiffness value. Second, the methods are tested on the system with uncertainty in both stiffness and damping. The equation of motion for the floating oscillator is: Figure 4.9: Floating Oscillator m1 0 y 1 y˙ k −k y¨ c −c 0 1 = u 1 + 1 + 0 −k k y2 −c c y˙2 m2 y¨2 (4.8) The control vector is bounded, resulting in the constraint: −1 ≤ u ≤ 1 (4.9) The boundary conditions considered are: ˙ = 0, y1 (0) = y2 (0) = y˙1 (0) = y2 (0) ˙ f ) = y2 (t ˙ f) = 0 y1 (tf ) = y2 (tf ) = 1, y1 (t (4.10) A rigid body mode exists, making Qj positive semi-definite. A virtual spring is added to the first mass, so that the potential energy of the first mass with respect to the final position is zero when the desired final position is satisfied, and Qj is now positive Chapter 4. Numerical Examples 56 definite. The new state space representation of the virtual system is m1 0 y¨1 c −c y˙1 (k + kp ) −k y 1 + + 1 = u 0 m2 y¨2 −c c y˙2 −k k y2 0 (4.11) where kp = 1 is the stiffness of the virtual spring added to the system. 4.2.1 Uncertain Stiffness For the first problem, m = 1, c = 0, and k is assumed to be the only uncertain parameter. With a nominal value of k = 1, k is assumed to lie between 0.7 ≤ k ≤ 1.3. The uncertain parameter k is discretized into 51 different values, evenly spaced over the range of uncertainty. Matrices Qj , Fj , and Lj are derived for j = 1..51. Since the control vector is constrained such that abs(u) ≤ 1 (4.12) the control vector is parameterized as in Eq. (2.6), with 5 switches or n = 5. The initial guess supplied to the gradient-based optimizer for the unknown control vector of switch times and amplitudes is h i h i T1 T2 T3 T4 T5 T6 = 1 2 3 4 5 6 (4.13) which yields 1 − 2e−0.7256s + 2e−1.6913s − 2e−2.9593s + 2e−4.2247s − 2e−5.1892s + e−5.9093s (4.14) as the solution to the minimax 5-switch control problem. The solution to the minimax QCP formulation for the floating oscillator requires specification of the final maneuver time. The final maneuver time is taken from the result for the 5-switch minimax controller, to allow for accurate comparison without use of the Bisection Algorithm. Chapter 4. Numerical Examples 57 Results for the floating oscillator with uncertainty in stiffness are shown in Figures 4.10-4.14 Mass 1 Position 1.5 QCP Time−Delay 1 0.5 0 0 1 2 3 Time (sec) 4 5 6 0 1 2 3 Time (sec) 4 5 6 1 Velocity 0.5 0 −0.5 −1 Figure 4.10: Response of Mass 1 for comparison of the 2-Time-Delay Minimax solution and the QCP Minimax Solution With the control vector discretized into 256 different values, Figure 4.10 and Figure 4.11 show that the system response for the QCP solution is very close to the traditional minimax time-delay solution. The differences in movement between the first and second masses are seen in Figure 4.10 and Figure 4.11. The position and velocity of the second mass vary in Chapter 4. Numerical Examples Mass 2 1.5 Position 58 1 QCP Time−Delay 0.5 0 0 1 2 3 Time (sec) 4 5 6 0 1 2 3 Time (sec) 4 5 6 0.4 Velocity 0.3 0.2 0.1 0 Figure 4.11: Response of Mass 2 for comparison of the 2-Time-Delay Minimax solution and the QCP Minimax Solution a relatively smooth manner due to the damper placed between the masses. In the derivation of the pseudo-energy function in the QCP framework, a hypothetical spring was attached to the first mass to deal with the rigid body displacement. This pseudospring was only used to create a virtual potential energy for use in the pseudo-energy function and limit the vibration due to the rigid body displacement of the system when determining the optimal control profile. The virtual spring was not included in the simulation of the system response. Chapter 4. Numerical Examples 59 1.5 Time−Delay QCP 1 Control Input, U 0.5 0 −0.5 −1 −1.5 0 1 2 3 Time (sec) 4 5 6 Figure 4.12: Control Input vs. Time for comparison of the 2-Time-Delay Minimax solution and the QCP Minimax Solution Figure 4.12 reveals that the QCP approach closely estimates the switch-times for the bang-bang control profile, yet does not represent the perfect bang-bang shape. This error is due to the disretization of the control profile and is reduced with a finer level of discretization. This illustrates potential for use of the QCP approach to assist in estimating the shape of the parameterization of the control profile used in the traditional minimax approach. Chapter 4. Numerical Examples 60 −3 1.6 x 10 1.4 Residual Energy, J 1.2 1 0.8 0.6 0.4 Time−Delay QCP 0.2 0 0.7 0.8 0.9 1 1.1 Stiffness, k 1.2 1.3 1.4 Figure 4.13: Residual Energy vs. Stiffness for the QCP Minimax Solution The variation of the residual energy over the range of uncertainty is shown in Figure 4.13. The QCP approach very closely approximates the residual energy variation shown using the traditional minimax time-delay approach. Chapter 4. Numerical Examples 61 0.4 Residual Energy, J 0.3 0.2 0.1 0 −0.1 0 1 2 3 4 5 Final Maneuver Time, tf 6 7 8 Figure 4.14: Residual Energy vs. Final Maneuver Time for the QCP Minimax Solution Figure 4.14 verifies that the Bisection Algorithm can be used to determine the minimum final maneuver time for which the residual energy falls below an allowable threshold. 4.2.2 Uncertain Stiffness and Damping The same approach is carried out for the floating oscillator shown in 4.9, with uncertainty in both stiffness and damping. The range of uncertainty for stiffness remains the same, 0.7 ≤ k ≤ 1.3, with a nominal value of k = 1. The value of damping now Chapter 4. Numerical Examples 62 falls in the uncertain range, 0.1 ≤ c ≤ 0.3, with a nominal value of c = 0.2. k is discretized into 15 values and c into 15 different values, to yield m = 225 system models included in the optimization scheme. The initial guess supplied to the gradient-based optimizer is h i h i T1 T2 T3 T4 T5 T6 = 1 3 4 5 6 7 (4.15) for the 5-switch minimax controller, parameterized as shown in Eq. (2.6). Upper and lower bounds are set on the control vector as shown in Eq. (4.9). The solution for the minimax 5-switch control profile with uncertainty in both stiffness and damping is: u = 1 − 2e−0.8656s + 2e−1.9229s − 2e−3.1176s + 2e−4.4006s − 2e−5.2825s + e−5.8754s (4.16) Results for the floating oscillator with uncertainty in both stiffness and damping are shown in Figures 4.15-4.19. Chapter 4. Numerical Examples Mass 1 1.5 Position 63 Time−Delay QCP 1 0.5 0 0 1 2 3 Time (sec) 4 5 6 0 1 2 3 Time (sec) 4 5 6 1 Velocity 0.5 0 −0.5 −1 Figure 4.15: Response of Mass 1 for comparison of the 2-Time-Delay Minimax solution and the QCP Minimax Solution Figures 4.15 and 4.16 show that the QCP minimax approach fairly accurately determines a control profile that approximates the response of the system to the traditional minimax time-delay solution. Slight differences exist due to the discretization of the control profile used in the QCP minimax approach that does not allow for the perfect bang-bang structure given by the traditional minimax time-delay solution. The error shown can be reduced with a finer discretization of the control profile. Results were generated using a discretization of N = 128. Chapter 4. Numerical Examples Mass 2 1.5 Position 64 1 Time−Delay QCP 0.5 0 0 1 2 3 Time (sec) 4 5 6 0 1 2 3 Time (sec) 4 5 6 0.5 Velocity 0.4 0.3 0.2 0.1 0 Figure 4.16: Response of Mass 2 for comparison of the 2-Time-Delay Minimax solution and the QCP Minimax Solution Chapter 4. Numerical Examples 65 1.5 Time−Delay QCP 1 Control Input 0.5 0 −0.5 −1 −1.5 0 1 2 3 Time 4 5 6 Figure 4.17: Control Input vs. Time for comparison of the 2-Time-Delay Minimax solution and the QCP Minimax Solution Figure 4.17 shows that the QCP approach fairly accurately estimates the switching structure determined using the traditional minimax approach. The QCP results could, therefore, be used to estimate the parameterization of the control profile used in the traditional minimax time-delay approach. Square Root of Residual Energy Chapter 4. Numerical Examples 66 0.05 0.04 0.03 0.02 0.01 0 0.4 1.4 0.3 1.2 1 0.2 Damping, c 0.1 0.8 Stiffness, k Figure 4.18: Square Root of Residual Energy vs. Damping vs. Stiffness for the QCP Minimax Solution Chapter 4. Numerical Examples 67 Residual Energy, J 0.4 0.3 0.2 0.1 0 −0.1 0 1 2 3 4 5 Final Maneuver Time, tf 6 7 8 Figure 4.19: Residual Energy vs. Final Maneuver Time for the QCP Minimax Solution The variation of the residual energy with the stiffness and damping values are shown in Figure 4.18. Figure 4.19 verifies that the Bisection algorithm can be used to determine minimum final maneuver time for which the residual energy falls below an allowable threshold. Chapter 4. Numerical Examples 4.3 68 Constrained Robust Control: Simple Elevator Problem A computational algorithm for time-optimal planning of elevators subject to a dynamic model and constraints is proposed in [16]. It is verified that the QCP approach, when applied using the nominal stiffness value alone, determines the time-optimal solution. The QCP minimax approach is then applied to the simple elevator model with linear state constraints for determination of a robust control profile. While a small scale elevator model is used, it should be noted that the following can be applied to a larger scale model, such as the state-space model with 54 states and one input that was fitted to the frequency response data of an elevator and involves a lumped parameter model of the ropes [15]. 4.3.1 Simple Elevator Model A free body diagram of the elevator model is shown in Figure 4.20. The input u is the drive sheave torque. The drive sheave rotational position is given by θ, with radius r. m1 and m2 represent the masses of the counterweight and elevator car, respectively. The drive sheave inertia is represented by I. Spring stiffnesses are given by k1 , k2 , and k. Several assumptions are made about this model in [16]: 1) a single hoist motor drives the elevator and is coupled to the drive sheave; 2) segments of the hoist and comp ropes that wrap around the drive and comp sheaves are inextensible during vertical movements of the car/counterweight; 3) segments of the hoist and comp ropes that lie between the pulleys and move in the vertical direction are made up of an extensible and rigid element; 4) the comp and drive sheaves are Chapter 4. Numerical Examples 69 Figure 4.20: Free Body Diagram of Simple Elevator Model massless and frictionless. Under these assumptions, the system is composed of six states: position and velocity of the elevator car, position and velocity of the counterweight, and position and velocity of the drive sheave taken at a point on the hoist cable just before making contact with the drive sheave. The equations of motion are Chapter 4. Numerical Examples 70 given by m1 x¨1 = Fs1 − Fs5 − m1 g (4.17a) m2 x¨2 = Fs2 − Fs4 − m2 g (4.17b) I θ̈ = u + r[Fs1 − Fs2 ] − µθ̇ 0 = Fs5 − Fs4 (4.17c) (4.17d) µ is the drive sheave damping constant. The spring forces are given by: Fs1 = k1 (x6 − x1 ) (4.18a) Fs2 = k2 (x3 − x2 ) (4.18b) Fs4 = k(x2 − x4 ) (4.18c) Fs5 = k(x1 − x5 ) (4.18d) Using the kinematic compatibility equations, rθ = x3 = −x6 and x4 = −x5 , the following representation of the system is derived. x¨1 0 0 0 x˙1 m1 0 0 0 m2 0 x¨2 + 0 0 0 x˙2 µ I x¨3 x˙3 0 0 r 0 0 r k k1 (k1 + k2 ) x1 2 k k x2 (k + + k ) −k + 2 p 2 2 2 rk1 −rk2 r(k1 + k2 ) x3 −m1 0 h i h i = 0 u + −m2 g 0 1 (4.19) A pseudo spring with stiffness kp = 1 is added to m2 , the mass that represents the elevator car, to eliminate the rigid body mode of vibration and transforms the positive semi-definite stiffness matrix into a positive definite stiffness matrix. The nominal Chapter 4. Numerical Examples 71 parameter values are taken from [16] as: k1 = k2 = k = 10000 N m m1 = 100 kg m2 = 110 kg I = 100 kgm2 (4.20) r = 1m g = 10 m s2 µ=0 A number of constraints are typically applied to the elevator model. These constraints vary, depending on the input actuator characteristics and the comfort features desired by the passengers. For this problem, linear constraints are placed on the trajectory of the elevator car and the input is bounded. Typically, nonlinear constraints, such as box constraints on power, are also implemented. In order to include nonlinear constraints, Schlemmer [16] linearizes the nonlinear constraints for implementation as linear constraints. Purely linear constraints are considered in this thesis. 4.3.2 Computational Algorithm for Time-Optimal Problem The approach used in [16] requires that the system be driven by a single input (torque at the hoist motor), and that the elevator hoist dynamics be adequately modeled by a linear system (shown in [15]). A special optimality condition is then derived and used to determine the time-optimal trajectory. The special optimality condition determines a bang-bang control solution for systems with general mixed constraints on the inputs and states, a generalization of the well known classical results for time-optimal control Chapter 4. Numerical Examples 72 for linear systems under box constraints on the inputs([4],[8]). In the derivation in [16], the dynamic system model is converted into a canonical form and the time-optimal problem is mapped into the space of the canonical variables. Using Pontryagin’s minimum principle, it is shown that everywhere along the optimal path, at least one inequality constraint must be active. For a single input system, such as the high-rise elevator problem, a single active constraint completely prescribes the solution. The computational algorithm is shown now, which determines the switching structure of the constraints and the corresponding switching times. Consider a linear time-invariant system of the form ẋ = Ax + Bu + c (4.21) with constraints given by C(x, u) ≤ 0 (4.22a) D(x) ≤ 0 (4.22b) ψ(x(t0 ), x(tf )) = 0 (4.22c) where x(t) ∈ Rn is the state vector, u(t) ∈ Rm is the control vector, and the matrices A ∈ Rn×n , B ∈ Rn×m , and c ∈ Rn are constants. Eq. (4.22a) represents mixed state and control constraints, while Eq. (4.22b) represents state constraints, both of which are given by sufficiently smooth vector functions. For the system given by Eq. (4.21), the time-optimal planning problem is Chapter 4. Numerical Examples 73 written as: M inimize tf (4.23a) Subject to C(x, u) ≤ 0 (4.23b) D(x) ≤ 0 (4.23c) ψ(x(t0 ), x(tf )) = 0 (4.23d) The optimal control problem (OP) in Eq. (4.23) is reformulated as a Sequential Optimization Problem (SOP), where the final time tf is fixed at each iteration, until the minimum optimal final time is determined. The scalar variable S(t), is defined as the distance away from the constraints. The SOP is represented as: M inimize S(tf ) (4.24a) ˙ =0 Subject to S(t) (4.24b) Cj (x(t), u(t)) ≤ S(t), ∀j (4.24c) Dk (x(t)) ≤ S(t), ∀k (4.24d) The distance away from the constraint at the final time is minimized. Eq. (4.24b) guarantees that one constraint is active at any given time when the cost is equal to zero. Eq. (4.24c) represents mixed state and control constraints, while Eq. (4.24d) represents state constraints. The cost given in Eq. (4.24a), determines whether or not the optimal solution of the SOP is feasible and/or optimal for the optimal control (OP) problem. Three cases exist: 1) if the cost is negative, it can be concluded that the optimal solution to the SOP is feasible, but not optimal for the OP; 2) if the cost is positive, it can be concluded that the optimal solution to the SOP is not feasible for the OP, and that no feasible solution to the OP exists with the specified final time Chapter 4. Numerical Examples 74 tf ; 3) if the cost equals zero, at least one constraint is active at any given time and the optimality condition for the OP has been satisfied. The SOP can be solved by adapting tf such that the appropriate change in sign of S(tf ) is achieved. When S(tf ) = 0 and tf can no longer be decreased without causing S(tf ) to become positive, the optimal solution to the overall optimization problem has been found. 4.3.3 QCP Formulation A QCP approach is used to approximate the time-optimal solution for the simple elevator problem and is shown below. The time-optimal solution determined using the SOP approach discussed in the previous section is compared to the time-optimal solution determined using a QCP formulation. This first involves derivation of the residual energy expression used in the QCP approach for the simple elevator problem. Once results are generated for the time-optimal problem, robustness is considered with the application of the QCP minimax approach for determining a control profile that is robust to parameter uncertainty. As seen in Eq. (4.19), the simple elevator model has the following form: M ÿ + C ẏ + Ky = D1 u + D2 g The first order representation of this system is y ẏ 0 I = ÿ −M −1 K −M −1 C ẏ 0 0 + −1 u + −1 g M D1 M D2 (4.25) (4.26) Chapter 4. Numerical Examples 75 which can be written in discrete time as: k x(k + 1) = G x0 + k X k−i G H1 u(i) + k X i=1 Gk−i H2 g (4.27) i=1 G, H1 and H2 are the discrete time state space matrices of the system. u(i) is the discrete time control vector. The states at final time, N + 1, can be written as: N x(N + 1) = G x0 + N X N −i G H1 u(i) + i=1 N X GN −i H2 g (4.28) i=1 Substituting Eq. (4.28) into Eq. (3.5), the residual energy can be represented as 2J = uT Qu + F T u + L (4.29a) where L = xTf Em xf + (GN x0 )T Em GN x0 − 2xTf Em GN x0 (4.29b) + f2 Em f2T − 2xTf Em f2 + 2(GN x0 )T Em f2 F T = 2(GN x0 − xf )T Em f1 + 2f1T Em f2 (4.29c) (GN −1 H1 )T Em GN −1 H1 (GN −1 H1 )T Em GN −2 H1 ... (GN −1 H1 )T Em H1 (GN −2 H1 )T Em GN −1 H1 (GN −2 H1 )T Em GN −2 H1 ... (GN −2 H1 )T Em H1 Q= : : : T N −1 T N −2 T H1 Em (G H1 ) H1 Em G H1 ... H1 Em H1 (4.29d) h i f1 = GN −1 H1 GN −2 H1 ... H1 h i f2 = GN −1 H2 g GN −2 H2 g ... H2 g (4.29e) (4.29f) This residual energy derivation can now be used in the QCP formulation shown in Chapter 3 for approximating the time-optimal control profile, as well as the minimax time-delay control profile. Chapter 4. Numerical Examples 4.3.4 76 Results Time-Optimal Control Time-optimal control of the simple elevator model is considered first. The SOP approach used in [16] is compared to the QCP approach for determining a timeoptimal control profile for the simple elevator model shown in Eq. (4.19) with the nominal parameter values given in Eq. (4.20). The two constraints considered are taken from [16]. A constraint is placed on the maximum/minimum allowable drive sheave torque and is given as: −500 N N ≤ u ≤ 500 m m (4.30) A constraint is placed on the elevator car velocity and is given as: −5 m m ≤ x2 ≤ 5 s s (4.31) The SOP problem for determining the time-optimal solution for the system in Eq. (4.19) with the constraints in Eq. (4.30) and Eq. (4.31) is written as: M inimize S(tf ) (4.32a) ˙ =0 Subject to S(t) (4.32b) S(t) ≥ x2 (t) − 5 (4.32c) S(t) ≤ x2 (t) + 5 (4.32d) x(t0 ) = x(1) (4.32e) x(tf ) = xf (4.32f) −500 ≤ u ≤ 500 (4.32g) Using the expression for the states at final time, N + 1, in Eq. (4.28), x(t) is written Chapter 4. Numerical Examples 77 in terms of u(t). Eq. (4.32) can now be written as: M inimize S(tf ) (4.33a) ˙ =0 Subject to S(t) (4.33b) N S(t) ≥ [G x( 1) + S(t) ≤ [GN x( 1) + N X i=1 N X N −i G H1 u(i) + GN −i H1 u(i) + i=1 xf = GN x( 1) + N X i=1 GN −i H1 u(i) + N X i=1 N X GN −i H2 g]row 5 − 5 (4.33c) GN −i H2 g]row 5 + 5 (4.33d) i=1 N X GN −i H2 g (4.33e) i=1 −500 ≤ u ≤ 500 (4.33f) where row 5 is considered in the expression for the states at final time, N + 1, used in the constraints in Eq. (4.33c) and Eq. (4.33d), since row 5 determines the position of the elevator car velocity given in the state-space model for the system. The constraints on the initial and final positions given in Eq. (4.32e) and Eq. (4.32f) are combined to form the constraint in Eq. (4.33e). In order to determine the time-optimal solution using the QCP approach, the residual energy expression in Eq. (4.29a) was used to solve the QCP problem of the form Eq. (3.10), while considering the nominal parameter values only. This means that only one quadratic constraint representing the residual energy at the end of the maneuver for the nominal system was included in the QCP problem. The initial and final positions considered are x(1) = [x1 (0) x2 (0) x3 (0) ẋ1 (0) ẋ2 (0) ẋ3 (0)] = [0 0 0 0 0 0] (4.34) xf = [x1 (tf ) x2 (tf ) x3 (tf ) ẋ1 (tf ) ẋ2 (tf ) ẋ3 (tf )] = [−30 30 30 0 0 0] for the time-optimal problem. The time-optimal final maneuver time was determined to be 10.0838 seconds with the control vector discretized into N = 256 different val- Chapter 4. Numerical Examples 78 ues. Results for both the QCP approach and the SOP approach are shown in Figures 4.21-4.24. These results verify that the QCP approach accurately determines the time-optimal solution for the system with linear state constraints. Counterweight 0 Position −10 SOP QCP −20 −30 −40 0 2 4 6 Time (sec) 8 10 0 2 4 6 Time (sec) 8 10 2 Velocity 0 −2 −4 −6 Figure 4.21: Time-Optimal Response of the Counterweight for Comparison of the SOP approach and QCP approach Chapter 4. Numerical Examples 79 Elevator Car 40 Position 20 SOP QCP 0 −20 0 2 4 6 Time (sec) 8 10 0 2 4 6 Time (sec) 8 10 6 Velocity 4 2 0 −2 Figure 4.22: Time-Optimal Response of the Elevator Car for Comparison of the SOP approach and QCP approach Figure 4.22 shows that the elevator car velocity remains below the required threshold of, 5 ms . Chapter 4. Numerical Examples 80 Hoist Rope 40 Position 30 20 SOP QCP 10 0 0 2 4 6 Time (sec) 8 10 0 2 4 6 Time (sec) 8 10 5 Velocity 4 3 2 1 0 Figure 4.23: Time-Optimal Response of the Hoist Rope for Comparison of the SOP approach and QCP approach Chapter 4. Numerical Examples 81 600 Control Input, U 400 200 QCP SOP 0 −200 −400 −600 0 2 4 6 8 10 Time (sec) Figure 4.24: Time-Optimal Control Profile for Comparison of the SOP approach and QCP approach Comparing Figures 4.22-4.24, it is clear that the fast switching of the control profile occurs during the span of the maneuver when the elevator car velocity is maintained at the upper threshold of the maximum allowable velocity. Without constraints on the elevator car velocity, fewer switches would be required to complete the desired time-optimal maneuver. Chapter 4. Numerical Examples 82 Robust Time-Optimal Control The QCP minimax approach is applied to the simple elevator model for determination of a control profile that is robust to uncertainty in spring stiffness and satisfies linear constraints on the elevator car velocity. The springs represented by stiffnesses, k1 , k2 , and k are assumed to fall in the range of uncertainty, 7000 ≤ ki ≤ 13000, each with a nominal value of 10000. The range of uncertainty covers a thirty percent deviation in the positive and negative directions for the spring stiffness values. The range of uncertainty is discretized into 51 different values, evenly spaced over the range of uncertainty. Results are shown in Figures 4.25-4.30. Counterweight Position 20 0 −20 −40 0 2 4 6 8 10 Time (sec) 12 14 16 0 2 4 6 8 10 Time (sec) 12 14 16 2 Velocity 0 −2 −4 −6 Figure 4.25: Response of the Counterweight using the QCP Minimax Approach Chapter 4. Numerical Examples 83 Figures 4.25-4.27 show the system response with the nominal stiffness value. A constraint was not included to force the residual energy to zero for the nominal value, although a constraint could be easily included. The counterweight, elevator car, and hoist rope move smoothly to the desired position, with a slight steady-state error that is expected with the minimax solution. Elevator Car 30 Position 20 10 0 −10 0 2 4 6 8 10 Time (sec) 12 14 16 0 2 4 6 8 10 Time (sec) 12 14 16 6 Velocity 4 2 0 −2 Figure 4.26: Response of the Elevator Car using the QCP Minimax Approach Figure 4.26 shows that the elevator car velocity remains below the required threshold of, 5 ms . Chapter 4. Numerical Examples 84 Hoist Rope 40 Position 20 0 −20 0 2 4 6 8 10 Time (sec) 12 14 16 0 2 4 6 8 10 Time (sec) 12 14 16 6 Velocity 4 2 0 −2 Figure 4.27: Response of the Hoist Rope using the QCP Minimax Approach Chapter 4. Numerical Examples 85 600 400 Control Input 200 0 −200 −400 −600 0 2 4 6 8 Time 10 12 14 16 Figure 4.28: Control Profile for the Simple Elevator Model and the QCP Minimax Approach The control profile reveals a bang-bang solution for the minimax problem, an expected result for a minimum-time solution. The inclusion of a constraint on the allowable elevator car velocity greatly increases the number of constraints, since a constraint is included for each time step. The large number of additional constraints, combined with the many constraints required to guarantee a robust solution, leads to a large number of switches in the control profile shown in Figure 4.29. Chapter 4. Numerical Examples 86 180 160 Residual Energy, J 140 120 100 80 60 40 20 0 0.7 0.8 0.9 1 Stiffness, k 1.1 1.2 1.3 4 x 10 Figure 4.29: Residual Energy vs. Spring Stiffness for the Simple Elevator Model and the QCP Minimax Approach The variation of the residual energy with spring stiffness value is shown in Figure 4.29. Results were generated for 51 different values evenly spaced throughout the range of uncertainty. Chapter 4. Numerical Examples 87 650 600 550 Residual Energy, J 500 450 400 350 300 250 200 150 5 10 15 20 25 Final Maneuver Time, tf 30 35 40 Figure 4.30: Residual Energy vs. Final Maneuver Time for the Simple Elevator Model and the QCP Minimax Approach Figure 4.30 determined a minimum final maneuver time of 16 seconds, for which the residual energy falls below an allowable threshold of 170. A more accurate final maneuver time can be determined with a finer discretization of the unknown control vector and plot of residual energy vs. final maneuver time. Chapter 4. Numerical Examples 4.4 88 Conclusion Results for the harmonic oscillator and floating oscillator show that the QCP minimax approach provides a fairly accurate approximation to the results obtained using the traditional minimax approach for small-scale systems. The results generated for the simple elevator problem reveal potential for application of the QCP minimax approach to large-scale systems with linear state and control constraints for determining a control profile that is robust to parameter uncertainties. Chapter 5 Conclusion 5.1 Summary Various methods have been developed for determining control profiles that are robust to parameter uncertainties. These concepts stem from the idea of input shaping, an idea that exploits the fact that the designer can anticipate the response of an open-loop system to a specified input. It makes sense to shape the input in order to satisfy final value constraints which minimize the vibration at the end of maneuver, a condition necessary to achieve precise position control. The design of input shapers/time-delay filters requires knowledge of the system model parameters, such as stiffness and/or damping, so insensitivity to modeling errors is important. In order to achieve robustness to parameter uncertainty, input shaping methods enhanced robustness by including constraints on the sensitivity states of the system, achieved in several different ways, for example, by placing multiple zeros at the locations of the poles of the system in the case of robust time-optimal control, or forcing the partial derivatives of the systems states with respect to natural frequency 89 Chapter 5. Conclusion 90 to zero at the end of the maneuver. These methods improve robustness in the vicinity of the nominal parameter values of the system. With the development of time-delay filters came the concept of minimization of the residual energy at the end of the maneuver, and the definition of the pseudo-energy function, a metric especially useful for achieving robustness to parameter uncertainty in vibration minimization of multi-mode systems. In 2002, Singh developed the traditional time-delay approach for minimization of the maximum of the pseudo-energy function at the end of the maneuver, a crucial method for greatly improving robustness to modeling uncertainties, which includes a known range of parameter uncertainty and allows for the attainment of robustness beyond the mere vicinity of the nominal parameter values. Excellent results are achieved with the traditional minimax time-delay approach, however, several issues hinder the application to large-scale systems. The traditional minimax time-delay approach requires knowledge of the shape of the parameterization of the control vector amplitudes and corresponding vector of time-delays. Developing the parameterization of the control vector for large-scale systems can be difficult. Even if a reasonable parameterization can be assumed, the traditional approach relies on the use of gradient-based solvers for solution to the optimal control problem. Without a good initial guess, the solver may result in a local solution only. This provided the motivation to develop methods to approximate the traditional minimax time-delay results in a convex framework, in order to guarantee a global solution and remove the reliance on a good initial guess. Several methods then arose to approximate the traditional minimax time-delay solution in a convex framework, without the reliance on knowledge of the shape of the parameterization of the control profile. These methods involve the formulation of the Chapter 5. Conclusion 91 minimax problem as a LP optimization problem. Since the pseudo-energy function is in a quadratic form, new metrics, for example, using an L1 norm and an L∞ norm were implemented to define pseudo-energy constraints at the end of maneuver in a linear form. While the LP problems could be solved relatively efficiently without specific control vector parameterization or a good initial guess, the solutions do not adequately approximate results achieved when using the true quadratic form of the system’s residual energy defined by the pseudo-energy function in the traditional minimax approach. Additional methods were developed to approximate the quadratic pseudo-energy function in an LP framework, but these methods fail to achieve a reasonable approximation of the quadratic form, for large-scale systems when the size of the optimization problem grows along with both the size of the system itself and the level of approximation required to minimize approximation error. In 2008, Conord and Singh formulated the minimax problem in a SDP framework using LMI’s to form quadratic constraints on the residual energy at the end of the maneuver, a method which reaps all of the benefits of the LP approaches, and also allows for inclusion of the true quadratic form of the pseudo-energy function. The LMI approach to the minimax problem gave excellent results, which adequately approximate the results given by the traditional minimax time-delay approach. Since powerful solvers exist for solving SDP’s, this extended the application of the traditional minimax approach to large-scale systems. The success of the LMI minimax approach provided motivation for the development of the QCP minimax approach for application to large-scale systems. While much work has been done to develop solvers for large-scale SDP’s, efficiency is still limited with the application to large-scale systems. In the LMI minimax approach, the number of LMI constraints quickly exceeds the size of the unknown con- Chapter 5. Conclusion 92 trol vector, especially as the level of parameter uncertainty increases along with the size of the system. SDP solvers can encounter efficiency and complexity issues, thus, limiting the efficacy in large-scale systems. The QCP minimax approach was developed in this thesis to exploit the most powerful software available for solving large scale optimization problems and allow for accurate approximation of the results given by the traditional minimax time-delay approach. 5.2 Contribution Like the LMI minimax approach, the QCP minimax approach utilizes a convex optimization framework to guarantee a globally optimal solution and neither requires knowledge of the shape of the control vector parameterization nor approximation of the true quadratic form of the system’s residual energy. Since much work has been done to develop solvers for large-scale optimization problems in the QCP framework, this method provides a useful alternative for approximating the excellent results obtained using the traditional minimax time-delay approach for large-scale systems, and could allow for application to large-scale systems where SDP solvers encounter efficiency and complexity problems. 5.3 Future Work Several issues remain open for further investigation. These issues are discussed in the following: Chapter 5. Conclusion 93 1) It could be beneficial to study the effects of the change in the desired system maneuver, as different initial/final state constraints seem to have an effect on the shape of the optimal control profile in systems with flexible modes of vibrations only. 2) The quadratic pseudo-energy function provides a convex cost that quantifies the error between the actual and desired final states at the end of the maneuver, in order to guarantee a global solution. The stiffness and mass matrices are used to weigh the contributions of the square of the position and velocity errors at the end of the maneuver. It could be useful to determine a rigorous approach for choosing the form of the pseudo-energy function that best fits a given problem and further investigate the effect of different weights on the quadratic position and velocity terms used in the definition of the pseudo-energy function. Kumar and Singh recently proposed a method for assigning modal weights to the residual energy cost, to scale the potential energy contribution of the modes of vibration that are considered to be the most significant for a given problem. These modal weighted minimax input shapers improved the performance when compared to the traditional minimax time-delay approach, which equally weighs all of the modes. 3) For systems with rigid-body modes, the effects of the value of the pseudo-springs added to the system needs to be better understood. A rigorous study as to the relationship between pseudo-stiffness and the optimal solution could be performed, with the hope of determining a formal method for choosing the pseudo-stiffness values for a given problem. It could also be useful to determine an additional method for dealing with the rigid-body displacement that fits into the QCP or LMI framework. The QCP solvers can handle semi-definite quadratic constraints, which might be useful in dealing with systems with rigid-body modes. However, a new method for constraining Chapter 5. Conclusion 94 the rigid-body displacement would need to be developed. Since the SDP formulation of the minimax problem requires the use of the Schur Complement for posing the LMI constraints, the LMI formulation cannot handle semi-definite constraints. 4) For nonlinear systems, the QCP framework could be utilized to extend the Sequential Linear Programming (SLP) approach in [21] to include second-order terms in the taylor series approximation of the nonlinear system used in the determination of a time-optimal control profile. Bibliography [1] H. S. Bae and J. C. Gerdes, Command Modification using Input Shaping for Automated Highway Systems with Heavy Trucks, Proceedings of the 2003 American Control Conference, vol. 1, IEEE, pp. 54–9. [2] W. L. Ballhaus, S. M. Rock, and A. E. Bryson Jr, Optimal Control of a TwoLink Flexible Robotic Manipulator using Time-Varying Controller Gains. Initial Experiments, Advances in the Astronautical Sciences, vol. 78, Univelt Inc, pp. 335–351. [3] S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004. [4] A. E. Bryson, Applied optimal Control, (1969). [5] T. Conord and T. Singh, Robust Input Shaper Design using Linear Matrix Inequalities, Proceedings of the IEEE International Conference on Control Applications (2008), 1470–1475. [6] D. Goldfarb and G. Iyengar, Robust Convex Quadratically Constrained Programs, Mathematical Programming, Series B 97, no. 3, 495–515. 95 Bibliography 96 [7] J. L. Junkins, Z. H. Rahman, and H. Bang, Near-Minimum-Time Maneuvers of Flexible Vehicles. a Liapunov Control Law Design Method, AIAA, Washington, DC, USA, pp. 299–310. [8] D. E. Kirk, Optimal Control Theory: An Introduction, ser. Electrical Engineering (1970). [9] Q. Liu and B. Wie, Robust Time-Optimal Control of Uncertain Flexible Spacecraft, Journal of Guidance, Control, and Dynamics 15, no. 3, 597–604. [10] J. Lofberg, Yalmip, http://control.ee.ethz.ch/joloef/yalmip.php, 2006. [11] Denny K. Miu and Sudarshan P. Bhat, Minimum Power and Minimum Jerk Position Control and its Applications in Computer Disk Drives, IEEE Transactions on Magnetics 27, no. 6 pt 1, 4471–4475. [12] J. C. Nash, The (Dantzig) Simplex Method for Linear Programming, Computing in Science and Engineering 2 (2000), no. 1, 29–31. [13] S. P. Panda and Lu Yan, Tutorial on Control Systems Design in Tape Drives, Proceedings of the 2003 American Control Conference, vol. 1, IEEE, pp. 1–17. [14] L. Y. Pao, T. N. Chang, and E. Hou, Input Shaper Designs for Minimizing the Expected Level of Residual Vibration in Flexible Structures, Proceedings of the 1997 American Control Conference, vol. 6, American Autom. Control Council, pp. 3542–6. [15] R. Roberts, Control of High-Rise/High-Speed Elevators, Proceedings of the 1998 American Control Conference, vol. 6, American Autom. Control Council, pp. 3440–4. Bibliography 97 [16] M. Schlemmer and S. K. Agrawal, A Computational Approach for Time-Optimal Planning of High-Rise Elevators, IEEE Transactions on Control Systems Technology 10 (2002), no. 1, 105–11. [17] N. C. Singer and W. P. Seering, Preshaping Command Inputs to Reduce System Vibration, Transactions of the ASME. Journal of Dynamic Systems, Measurement and Control 112, no. 1, 76–82. [18] T. Singh, Minimax Design of Robust Controllers for Flexible Systems, Journal of Guidance, Control, and Dynamics 25 (2002), no. 5, 868–75. , Optimal Reference Shaping for Dynamical Systems, CRC Press, 2010. [19] [20] T. Singh and W. Singhose, Tutorial on Input Shaping/Time-Delay Control of Maneuvering Flexible Structures, Proceedings of the American Control Conference, vol. 3, Institute of Electrical and Electronics Engineers Inc., 2002, pp. 1717– 1731. [21] T. Singh and P. Singla, Sequential Linear Programming for Design of TimeOptimal Controllers, Proceedings of the 46th IEEE Conference on Decision and Control, IEEE, 2008, pp. 4755–60. [22] T. Singh and S. R. Vadali, Robust Time-Delay Control of Multimode Systems, International Journal of Control 62, no. 6, 1319–39. [23] , Robust Time-Delay Control, Journal of Dynamic Systems, Measurement and Control, Transactions of the ASME 115 (1993), no. 2 A, 303–306. [24] , Robust Time-Optimal Control: Frequency Domain Approach, Journal of Guidance, Control, and Dynamics 17 (1994), no. 2, 346–53. Bibliography 98 [25] W. Singhose, S. Derezinski, and N. Singer, Extra-Insensitive Input Shapers for Controlling Flexible Spacecraft, Journal of Guidance, Control, and Dynamics 19, no. 2, 385–91. [26] W. E. Singhose, L. J. Porter, and W. P. Seering, Input Shaped Control of a Planar Gantry Crane with Hoisting, Proceedings of the 1997 American Control Conference, vol. 1, American Autom. Control Council, pp. 97–100. [27] O. J. M. Smith, Posicast Control of Damped Oscillatory Systems, Institute of Radio Engineers – Proceedings 45, no. 9, 1249–1255. [28] Jos F. Strum, Sedumi, http://sedumi.mcmaster.ca/, 2006. [29] K. C. Toh, R. H. Tutuncu, and M. J. Todd, Sdpt3, http://www.math.cmu.edu/reha/sdpt3.html. [30] C. F. Van Loan, Computing Integrals involving the Matrix Exponential, vol. AC23, pp. 395–404. [31] R. Venugopal and D. S. Bernstein, State Space Modeling and Active Control of Slosh, Proceedings of the 1996 IEEE International Conference on Control Applications., IEEE, pp. 1072–7. [32] M. Vossler and T. Singh, Minimax Defection/Energy-Limited Controllers for Flexible Structures: Linear Programming Approach, IEEE International Conference on Control Applications, 2008, pp. 1109–1114. [33] V. A. Yakubovich, Linear Matrix Inequalities in System and Control Theory (S. Boyd, L. E. Ghaoui, E. Feron, and V. Balakrishan), SIAM Review 37 (1995), no. 3, 479–479, Compilation and indexing terms, Copyright 2009 Elsevier Inc.

1/--страниц