close

Вход

Забыли?

вход по аккаунту

?

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.
Документ
Категория
Без категории
Просмотров
0
Размер файла
1 928 Кб
Теги
sdewsdweddes
1/--страниц
Пожаловаться на содержимое документа