Application of polynomial chaos in stability and ... - Semantic Scholar
Automatica 42 (2006) 789 – 795 www.elsevier.com/locate/automatica
Brief paper
Application of polynomial chaos in stability and control夡 Franz S. Hover ∗ , Michael S. Triantafyllou Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Received 29 October 2004; received in revised form 13 December 2005; accepted 3 January 2006
Abstract The polynomial chaos of Wiener provides a framework for the statistical analysis of dynamical systems, with computational cost far superior to Monte Carlo simulations. It is a useful tool for control systems analysis because it allows probabilistic description of the effects of uncertainty, especially in systems having nonlinearities and where other techniques, such as Lyapunov’s method, may fail. We show that stability of a system can be inferred from the evolution of modal amplitudes, covering nearly the full support of the uncertain parameters with a finite series. By casting uncertain parameters as unknown gains, we show that the separation of stochastic from deterministic elements in the response points to fast iterative design methods for nonlinear control. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Polynomial chaos; Control applications
1. Introduction Wiener’s polynomial chaos (Wiener, 1938) is fundamentally a framework for separating stochastic components of a system response from deterministic components. It derives from the Cameron–Martin theorem (Cameron & Martin, 1947), which establishes that a random process with finite second-order moments can be decomposed into an infinite, convergent series of polynomials in a random variable. The Hermite orthogonal polynomials were the first such series considered, and are known to be the optimum basis for Gaussian distributions, because the weighting function of the Hermite polynomials is the same as the probability density function (PDF). Modern analysis began with Ghanem and Spanos (1991), who studied the mechanics of structures having uncertain physical parameters, but with well-defined distributions. Polynomial chaos extended to the time domain enables the study of dynamical systems, and this approach has been applied successfully to large-scale problems in computational fluid mechanics 夡 This paper was not presented at any IFAC meeting. This paper was recommended for publication in the revised form by Associate Editor Henk Nijmeijer under the direction of Editor Hassan Khalil. ∗ Corresponding author. Tel.: +1 617 253 6762; fax: +1 617 253 8125. E-mail addresses: [email protected] (F.S. Hover), [email protected] (M.S. Triantafyllou).
0005-1098/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.automatica.2006.01.010
(Xiu, Lucor, Su, & Karniadakis, 2002). Treatment of other important distributions (e.g., uniform, Poisson) is achieved via the Askey scheme, discussed by Xiu and Karniadakis (2002). Each distribution type relates to a specific family of orthogonal polynomials, so that in many practical cases the optimum convergence properties can be realized. Functionally, we find that the polynomial chaos allows the deterministic parts of a response to be efficiently pre-computed, and the stochastic parts to be injected subsequently, at low computational cost. Hence, the polynomial chaos is an efficient alternative to Monte Carlo simulations for complex systems. The potential of such a tool applied to control engineering is significant, especially considering nonlinear systems, for which limited methods are available for design and analysis of even seemingly simple problems. Indeed, nonlinear systems like the one considered in this paper are stabilized today using feedback linearization where possible, Lyapunov functions where possible, and trial and error tuning with Monte Carlo testing. The first approach may place unrealistic measurement demands on the designer, and the second approach sometimes simply fails. The last approach, of course the most onerous from a computational point of view, often remains as the only viable path to progress in the design. In the case where a controller design has been made through some justifiable method or procedure, it is also desirable to
790
F.S. Hover, M.S. Triantafyllou / Automatica 42 (2006) 789 – 795
understand the distribution of the response in terms of the uncertainties that exist. The probabilistic point of view in both open- and closed-loop response has an intuitive appeal. For example, the manufacturing accuracy of physical components is typically expressed as a standard deviation of a (often normal) distribution. Given that stability is known, how does a specific distribution in the plant parameters map to the range of responses? The same might be said for initial conditions: the initial conditions for which a system must be stabilized sometimes occupy a distribution, with known statistics. The polynomial chaos enables one to directly describe the statistics of the system response. These concepts differ with the established procedure in control theory of defining uncertainty with a uniform probability distribution, i.e., having equally likely values between a hard minimum and a hard maximum. The polynomial chaos expansion appears to be a completely new technique for these problems, and has not been reported in the literature. As we will show, the method is noteworthy for the fact that random variables with large variance can be accommodated. We focus here on two main topics, suitable for general nonlinear systems: 1. Prediction of short-term statistics and stability, given parameters or initial conditions with known distribution. 2. Controller design, enabled by treating the unknown gains as uncertain parameters. The paper is focussed on four examples, highlighting the convergence properties and the computational tradeoffs that occur. We focus on the formulation employing Hermite polynomials with Gaussian distribution, which is not always suitable for stability analysis, since the normal distribution has infinite support. Nonetheless, if the approach is used in a pragmatic way, effective indicators of stability can be constructed. For tighter stability analysis, the uniform distribution and its associated Legendre polynomials are more directly relevant within the classical description of parametric uncertainty. We also note that the case where an uncertain parameter or variable is driven, perhaps through a dynamical system, by a random process, the Karhunen–Loeve spectral decomposition can be incorporated into the polynomial chaos (Lucor, Su, & Karniadakis, 2004; Masri, Smyth, & Traina, 1998). To our knowledge, this is the first time that polynomial chaos methods are applied to control problems. Because feedback control exists largely to limit the effects of uncertainty and disturbances in dynamical systems, however, it would seem to be a natural fit.
of Abramowitz & Stegun, 1972): G=
The Cameron–Martin theorem states that any second-order (i.e., finite variance), one-dimensional process G can be written as a weighted sum of Hermite polynomials in the Gaussian random variable (taken to have zero mean and unity variance), and that this series converges in the L2 sense (using the notation
gi Hei ().
(1)
i=0
Intuitively, the distribution G is constructed from modes which are Hermite polynomials: He0 () = 1, He1 () = , He2 () = 2 − 1, . . . . The Hermite polynomials possess an inner product defined with the Gaussian kernel: ∞ 2 e− /2 Hei ()Hej () d = eij ij , (2) −∞
where ij is the Kronecker delta function. Note that orthogonality holds only for the double product, and that the nonzero inner products are positive integers. Higher products are similarly defined, e.g.,
∞ −∞
e−
2 /2
Hei ()Hej ()Hek () d = eij k .
(3)
As a specific example of how the polynomial chaos can be used, let us consider the following system, which is the basis of the numerical results presented in the next section: x˙ = x + ( + c)y + xy, y˙ = x,
(4)
where , , are known constants and c is a random variable. When = 0, this system if stable is either overdamped or is a second-order oscillator, with natural frequency n = − − c, and damping ratio = −/2n . The regions of stability for the linear uncertain system are of course well understood—they correspond with those cases where + c < 0 and < 0. When = 0, however, this becomes a much more difficult problem, and is typical of a class of relatively low-order polynomial vector fields. This system resists easy analysis by the first or the second methods of Lyapunov, including Krasovski’s theorem. For example, with the Lyapunov function candidate V ( x ) = x 2 /2 + (−c − )y 2 /2, we ob2 ˙ tain V ( x ) = ( + y)x , a suggestive but inconclusive result. Monte Carlo simulation would appear to answer the question of performance and stability, but only indirectly since not all possible initial conditions and uncertainties c can be simulated. In this problem, the random variable c is expressed as a finite summation of modes, up to order P, along with the states: c=
2. The method
∞
P
ci Hei (),
i=0
y(t) =
P
yi (t)Hei (),
i=0
x(t) =
P i=0
xi (t)Hei ().
(5)
F.S. Hover, M.S. Triantafyllou / Automatica 42 (2006) 789 – 795
Monte Carlo 1000
α:
-1
y0(0):
0
800
β:
-4
y1(0):
1
0 0
N:
104
600
γ: c0:
400
c1:
0
200 0
-1
-0.5
0
0.5
mean: σ:
0.001568 0.1724
P:
1
1
Polynomial Chaos 1000 800
791
with the properties of time-invariant linear systems. Fig. 1 shows that the polynomial chaos is successful in generating the correct PDF (shown as a histogram) for the case of y(0) = (that is, y0 (0) = 0, y1 (0) = 1, y2 (0) = 0, . . .), x(0) = 0, no parametric uncertainty, and no nonlinear effects. In this case, only the second mode (He1 ) is nonzero, and so P = 1 is sufficient. We consider here the case of one random dimension only, but the extension to multiple random variables is straightforward, albeit computationally more expensive. Making this restriction, we can have either one random parameter or one random initial condition, but not both. We choose x(t =0)=0 in all calculations below, and consider the cases of random y(t = 0) and random c separately. 3. Numerical examples
600
3.1. Parametric uncertainty in the linear case
400 200 0
-1
-0.5
0
0.5
1
mean:
0
σ:
0.1733
y(t=3) Fig. 1. The case of a linear dynamics, with random initial position y; Mode 1 is the only one excited. Polynomial chaos accurately reconstructs a Gaussian distribution of y(t = 3), as shown in comparison with Monte Carlo simulation having 10,000 trials.
Substituting into the dynamic equation, we obtain P
x˙i (t)Hei () =
i=0
P
xi (t)Hei ()
i=0
+ +
P
ci Hei ()
P
xi (t)Hei ()
i=0 P
y˙i (t)Hei () =
i=0
P
yj (t)Hej ()
j =0
i=0
+
P
P
yj (t)Hej (),
j =0
xi (t)Hei ().
(6)
i=0
Next, employing the inner product, we multiply by the kth mode and perform a Galerkin projection: P P 1 ci yj (t)eij k x˙k (t) = xk (t) + yk (t) + ekk i=0 j =0
P P + xi (t)yj (t)eij k , ekk i=0 j =0
y˙k (t) = xk (t).
(7)
This construction is very similar to that given by Li and Ghanem (1998), for the Duffing oscillator. In the case of no parametric uncertainty (ci = 0, i > 0) and without the nonlinearity ( = 0), each mode evolves independently of the others, in accordance
We first turn to the situation of a linear dynamics, with Gaussian uncertainty in the parameter c. Coupling between the modal states yi (t) and the modes of the unknown parameter c occurs, and hence a higher number of modes has to be considered than in the first case above. As shown in Fig. 2, with P = 12, the polynomial chaos accurately reflects the asymmetric distribution of y(t = 3). The short-term time response is dominated by the first (mean) mode decaying from the initial value; the second mode essentially gives the root variance of the PDF. Near t = 1.7, the variance is small because y(t) is near its peak value. The short-term log plot further shows that the higher modes computed grow progressively smaller by about an order of magnitude per mode. This ordering is typical in the short term. We note, however, that the Hermite polynomials grow rapidly as Hei () = 2Hei−1 () − 2(i − 1)Hei−2 ()
(8)
and that this level of separation in the successive amplitudes is not observed in the present case. As a result, the higher modes will dominate for certain realizations. The capability of the polynomial chaos to accurately produce short-term statistics also applies to unstable responses; this condition was illustrated by Xiu and Karniadakis (2002) for a first-order system. Considering the long-term evolution, decay of P computed modes—to an arbitrarily low-amplitude level—is strongly suggestive of stability because any realization is a summation of modes, and hence decays as well. The time integration in Fig. 2 is performed using adaptive step-size Runge–Kutta routine with relative tolerance 10−6 , and rescaling to maintain accuracy. This example brings up the important issue of finite series truncation error. While the modes for the illustrated case with P = 12 approach zero exponentially in the long term, there do exist values of the random parameter c for which instability occurs, because c is a Gaussian random variable, and therefore has infinite tails. The “spring” term in the governing equation here is + c, or −4 + 0.5 ( = N (0, 1)), and hence the minimum value of which destabilizes the system is eight. The probability of exceeding eight standard deviations above the mean is about 6 × 10−16 ; for most practical purposes, the system is
792
F.S. Hover, M.S. Triantafyllou / Automatica 42 (2006) 789 – 795 Monte Carlo
2000
yi(t)
1
0
0 -0.5
y(t)
2000
0
1000
0
0.5
-4
c1 :
0.5
P:
0
0.5
1
1.5
2
2.5
0 -0.5
3
γ: 0
y1(0): 0
Polynomial Chaos
-0.5 0
seconds
0.5
12
mean: 0.1567 σ: 0.0626
y(t=3)
0
0
0 4
-5
8
-10
log10(abs(yi(t)))
log10(abs(yi(t)))
β:
0
N: 10000 mean: 0.1556 σ: 0.0635
1 0.5
-1
c0:
y0(0): 1
1000
0.5
α:
-5
-10 0 4
-15
-15
8 12
-20 0
0.5
1
1.5
2
2.5
3
seconds
0
10
20
30
40
50
60
seconds
Fig. 2. The case of a random parameter in the linear system, with standard deviation c = 0.5. Top left: 100 Monte Carlo runs, compared with the modal trajectories; only the first few modes are visible. Top right: histograms of the Monte Carlo simulation and from the polynomial chaos. Bottom left: log modal amplitudes for the short-term integration. Bottom right: log modal amplitudes for long-term integration. Mode numbers are listed on right-hand side of the last two graphs.
Table 1 Referencing Fig. 2, the maximum standard deviation of the zero-mean Gaussian variable c that still achieves decay of modes in the long-term, decreases as the number of modes P increases P
max
1 − F (−/)
4 8 12 16
1.31 0.83 0.65 0.55
1.1 × 10−3 7.2 × 10−7 3.8 × 10−10 1.8 × 10−13
The probability of this variable exceeding the critical value (c = 4) for system stability is given in the last column.
stable. Table 1 shows that the maximum value of for which long-term decay of the modes is achieved depends strongly on P. The concern is clearest for the case P = 4, for which decay of the modes holds not only for the nominal case = 0.5, but up to 1.31. Only 4/1.31 = 3.05 standard deviations of this variable will be tolerated by the actual system; this critical level is exceeded in about 0.1 percent of cases, and hence the decay of modes fails as a convincing stability indicator for P = 4. When P is larger, however, we see that accuracy improves, reaching a reasonable instability likelihood of 1.8 × 10−13 with 16 modes. These considerations raise the question of how the polynomial chaos can be practically applied in the general case. First, we note that the normal distribution and the Hermite chaos expansion are specific to this paper, and that different distribution
types with their corresponding orthogonal polynomials (Xiu & Karniadakis, 2002) may be both more accurate and more computationally efficient for suitable problems. At the same time, choosing the best chaos expansion can be difficult, because the uncertainty may have finite support, while the response may not, and vice versa. In these cases especially a high mode number P is required; low-order calculations can be expected to capture major transient effects, but do not represent the distributions well enough to accurately portray long-term trajectories. The situation is not unlike that in perturbation theory, wherein a slowly converging or even divergent series solution may result, yet still provide a useful solution by optimally truncating the series (Bender & Orszag, 1978).
3.2. Random initial condition in the nonlinear case Next we examine the case of a random initial condition and the nonlinear term in the governing equation. The parametric uncertainty in c is replaced with uncertainty in the initial condition y, that is, we set y1 (t = 0) = 1 and yj (t = 0) = 0, for j = 1. Recalling that with V ( x ) = x 2 /2 + (−c − )y 2 /2, V˙ ( x )=(+y)x 2 . V is not a true Lyapunov function, but it does suggest stability when y(t) < − / when > 0, or y(t) < 1 if we set =−1 and =1. Fig. 3 shows that the polynomial chaos accurately predicts the short-term distribution of y(t). Additionally, the system is indicated stable for all initial conditions
F.S. Hover, M.S. Triantafyllou / Automatica 42 (2006) 789 – 795
793
Monte Carlo
2000
yi(t)
1
0 -0.5
y(t)
2
0
-4
c1 :
0
0.5
γ: 1
y1(0): 1
1
Polynomial Chaos
2000
P:
-1
12
1000
-2 0 -0.5
-3 0
0.5
1
1.5
2
2.5
3
0
seconds
0.5
1
mean: 0.03339 σ: 0.149
y(t=3)
0
0
0 4
-5 8
-10 12
-15
log10(abs(yi(t)))
log10(abs(yi(t)))
β:
0
N: 10000 mean: 0.0359 σ: 0.15
1 0
-1
c0:
y0(0): 0
1000
0
α:
-5
-10 0
-15
4 8
-20 0
0.5
1
1.5
2
2.5
3
0
10
20
seconds
30
40
50
60
seconds
Fig. 3. Nonlinear case, with no parameter uncertainty, but uncertainty in the initial condition. Description is the same as for Fig. 2.
in y, subject to x(t =0)=0. This conclusion on stability follows from the same modal decay result described for the linear case. 3.3. Random parameter in the nonlinear case We now consider the nonlinear system with random variable c; for this case, we increased to a value of 5, so as to bring out more of the nonlinear effects, especially those which change the shape of the PDF’s. Fig. 4 shows a similar performance as in the prior cases. The polynomial chaos accurately generates the short-term PDF, and indicates stability of the system response to the initial conditions. The examples shown so far confirm that nonsmooth PDF’s are well described by a Hermite chaos (i.e., using Gaussian bases) of reasonable order, but that the accuracy of the result can rapidly deteriorate if the number of modes is too low. 3.4. Controller design Consider a design situation for nonlinear systems without successful linearization or Lyapunov-based controller design: Monte Carlo techniques are used in an iterative approach, using a scalar metric likely based on several properties of the transient response, as well as stability. The method proceeds via classical function minimization techniques, such as the Nelder–Meade simplex, which admits complex error functions without known derivatives. In this way, the space of potential gain sets is navigated so as to optimize the performance metric. In this context, the polynomial chaos examples above point to control design, and not just analysis, if random parameters
are replaced with available gains. In other words, we could have interpreted c in the example case above as a (proportional) control gain. As we show next, polynomial chaos allows the time integration of the deterministic parts of the response to be completely computed beforehand, so that the metric evaluation requires no integration. Fig. 5 illustrates this for the system under discussion. The function shown is T J= y(t)2 (1 + c2 ) dt, (9) 0
an LQR-like cost where the control penalty is , and the underlining of c denotes a specific realization of the random variable. The evaluation of this cost is made trivial by the polynomial chaos because, using the expansion,
T
J= 0
P
(yi (t)Hei (c))2 (1 + c2 ) dt
i=0
= (1 + c ) 2
P P
Hei (c)Hej (c)
P P
yi (t)yj (t) dt
0
i=0 j =0
= (1 + c2 )
T
Lij (c)Mij .
(10)
i=0 j =0
The definite integrals Mij can be computed independently of c. Then, for a given value of c taken from the distribution, Lij ( c ) is computed from the Hermite polynomials and hence J is given directly. The complete mapping from c to J is made solving differential equations only once—to obtain the set of
794
F.S. Hover, M.S. Triantafyllou / Automatica 42 (2006) 789 – 795 Monte Carlo
1000 1
yi(t)
0
0 -0.05
y(t)
1 0.5
-1
β:
-4
c0:
0
c1:
0.5
y0(0): 1
500
0.5
α:
0
γ: 5
y1(0): 0
N: 10000 mean: 0.005302 σ: 0.02585
0.05
Polynomial Chaos
1000
P:
0
12
500 -0.5 0
1
2
3
4
5
0 -0.05
6
0
seconds
y(t=6)
0
0
0 4
-5 8 12
-10
log10(abs(yi(t)))
log10(abs(yi(t)))
0.05
mean: 0.005422 σ: 0.02579
-5
-10 0 4
-15
8 -15
-20 0
1
2
3
4
5
6
0
10
seconds
20
30
40
50
60
seconds
Fig. 4. Nonlinear case, with parametric uncertainty and fixed initial condition. Description is the same as for Fig. 2. 3 nonlinear, p.c. nonlinear, M.C. linear, p.c. linear, M.C.
2.5
J
2
1.5
1
0.5
-4
-3
-2
-1
0
1
c
Fig. 5. Linear quadratic cost J as a function of gain c using polynomial chaos for the example system with P = 9, T = 4.5, = −1, = −1, = −5, = 0.1. Results are shown for the nonlinear and linear ( = 0) optimizations; results from polynomial chaos and Monte Carlo simulations are overlaid. We use c0 = 0, c1 = 1 and y0 (0) = 1, y1 (0) = 0 in solving the deterministic equations.
yi (t). As in previous cases, we have considered only one random variable, the control gain, and so the initial condition is fixed in this analysis. This is an inconsequential simplification, for if multidimensional polynomial chaos is employed, then the distribution of initial conditions in x and y sets the initial conditions of xi and yi . We see in Fig. 5 that the results considering nonlinear and linear systems are quite different. The nonlinear system has a minimum cost of 1.22 at c = −2.0, while for the linear system, we have a minimum cost of 0.82 at c = −1.0. If this optimum
c from the linear case is employed in the nonlinear system, the cost is 1.34. The mismatch is not unexpected, and motivates design procedures that can take into account the nonlinearity in a quantitative way. The figure shows both the norm computed through the polynomial chaos and the norm via simulation; they are generally in agreement to the fourth significant digit. In this example, we used a low-order expansion (P = 9) to accurately study the short-term system behavior, but the modes are divergent in the long term. This unlikely but still useful condition results from the fundamental tradeoff between stability and accuracy mentioned previously for distributions with infinite support. Accurate description of the system response for random values far from the mean can only occur with high mode number P. Higher P, however, implies also that the tails are more accurately reflected in the response, i.e., the modes are unstable. In the present example, P was fixed and the distribution parameters (specifically c1 ) were chosen so that the range of gains appropriate to the control problem were within several standard deviations, i.e., in a numerically well-behaved neighborhood, of the mean value. The accuracy was checked by making parallel time-domain simulations, and the comparison indicated in the figure is poor more than several deviations away from c’s mean value. 4. Summary We demonstrated that the polynomial chaos method applies to two major areas in nonlinear systems analysis: assessing transient response and long-term stability. The former is illustrated by calculation of explicit probability density functions at
F.S. Hover, M.S. Triantafyllou / Automatica 42 (2006) 789 – 795
a given time, and the latter by decay of the modal amplitudes, for an arbitrarily large number of modes P. In both cases, the uncertainty can occur in either or both of the model parameters and the initial conditions. For random parameters with infinite support, stability can only be investigated in a pragmatic way due to truncation errors, and other chaos expansions may be used when specific stability bounds are the main concern. The construction of statistics requires random samples, but stability is a one-time calculation, namely the temporal evolution of the modes: if they are converging to zero, then stability is indicated. Recasting a parametric uncertainty problem as a gain selection problem allows polynomial chaos to be applied in control system design. Not only can statistics of the response be quickly calculated at any point in time, but so can statistics of functions of the state, including norms, peak value, time to arrive, control expense, and so on. With the particular expansion we chose for this paper, i.e. Hermite/Gaussian, we are left still with the question of long-term behavior for the chosen gain or gains, because design to minimize a short-term cost with a truncated series does not prove stability. A straightforward procedure to subsequently demonstrate stability, however, would be to implement the chosen gains as fixed values, and perform the polynomial chaos analysis again to assess the long-term trajectories.
Acknowledgements Work was supported by the Office of Naval Research’s National Naval Responsibility (NNR) program, grant number N00014-02-1-0623, monitored by Dr. Sharon Curtin-Beerman. We also acknowledge ONR’s Autonomous Operations Future Naval Capability (AOFNC) program, grant N00014-02-C0202, monitored by J. Valentine. References Abramowitz, M., & Stegun, I. A. (1972). Handbook of mathematical tables. New York: Dover Publications. Bender, C. M., & Orszag, S. A. (1978). Advanced mathematical methods for scientists and engineers. New York: McGraw-Hill.
795
Cameron, R. H., & Martin, W. T. (1947). The orthogonal development of non-linear functionals in series of Fourier–Hermite functionals. Annals of Mathematics, 48, 385–392. Ghanem, R., & Spanos, P. (1991). Stochastic finite elements: A spectral approach. New York: Springer. Li, R., & Ghanem, R. (1998). Adaptive polynomial chaos expansions applied to statistics of extremes in nonlinear random vibration. Probability Engineering Mechanics, 13, 125–136. Lucor, D., Su, C.-H., & Karniadakis, G. E. (2004). Generalized polynomial chaos and random oscillators. International Journal for Numerical Methods in Engineering, 60, 571–596. Masri, S. F., Smyth, A. W., & Traina, M. I. (1998). Probabilistic representation and transmission of nonstationary processes in multi-degree-of-freedom systems. Journal of Applied Mechanics, 65, 398–409. Wiener, N. (1938). The homogeneous chaos. American Journal of Mathematics, 60, 897–936. Xiu, D., & Karniadakis, G. Em. (2002). The Wiener–Askey polynomial chaos for stochastic differential equations. SIAM Journal on Scientific Computing, 24, 619–644. Xiu, D., Lucor, D., Su, C.-H., & Karniadakis, G. Em. (2002). Stochastic modelling of flow-structure interactions using generalized polynomial chaos. Journal of Fluids Engineering, 124, 51–59. Franz Hover received his M.S. and Sc.D. degrees from the Woods Hole Oceanographic Institution/Massachusetts Institute of Technology Joint Program in Applied Ocean Physics and Engineering. He is currently a Principal Research Engineer with the Department of Mechanical Engineering at the Massachusetts Institute of Technology, with research interests in the area of applied control and robotics, fluid mechanics, and offshore structures. Dr. Hover is a Member of the IEEE and the ASME. Michael Triantafyllou is a Professor of Mechanical and Ocean Engineering at the Massachusetts Institute of Technology, Director of the Center for Ocean Engineering, and Director of the Testing Tank and Marine Hydrodynamics Laboratories. His research interests are in the dynamics and control of marine vehicles, and biomimetics. He is a Member of the American Physical Society, the Society of Naval Architects and Marine Engineers, and the International Society of Offshore and Polar Engineers.
Brief paper
Application of polynomial chaos in stability and control夡 Franz S. Hover ∗ , Michael S. Triantafyllou Department of Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, MA 02139, USA Received 29 October 2004; received in revised form 13 December 2005; accepted 3 January 2006
Abstract The polynomial chaos of Wiener provides a framework for the statistical analysis of dynamical systems, with computational cost far superior to Monte Carlo simulations. It is a useful tool for control systems analysis because it allows probabilistic description of the effects of uncertainty, especially in systems having nonlinearities and where other techniques, such as Lyapunov’s method, may fail. We show that stability of a system can be inferred from the evolution of modal amplitudes, covering nearly the full support of the uncertain parameters with a finite series. By casting uncertain parameters as unknown gains, we show that the separation of stochastic from deterministic elements in the response points to fast iterative design methods for nonlinear control. 䉷 2006 Elsevier Ltd. All rights reserved. Keywords: Polynomial chaos; Control applications
1. Introduction Wiener’s polynomial chaos (Wiener, 1938) is fundamentally a framework for separating stochastic components of a system response from deterministic components. It derives from the Cameron–Martin theorem (Cameron & Martin, 1947), which establishes that a random process with finite second-order moments can be decomposed into an infinite, convergent series of polynomials in a random variable. The Hermite orthogonal polynomials were the first such series considered, and are known to be the optimum basis for Gaussian distributions, because the weighting function of the Hermite polynomials is the same as the probability density function (PDF). Modern analysis began with Ghanem and Spanos (1991), who studied the mechanics of structures having uncertain physical parameters, but with well-defined distributions. Polynomial chaos extended to the time domain enables the study of dynamical systems, and this approach has been applied successfully to large-scale problems in computational fluid mechanics 夡 This paper was not presented at any IFAC meeting. This paper was recommended for publication in the revised form by Associate Editor Henk Nijmeijer under the direction of Editor Hassan Khalil. ∗ Corresponding author. Tel.: +1 617 253 6762; fax: +1 617 253 8125. E-mail addresses: [email protected] (F.S. Hover), [email protected] (M.S. Triantafyllou).
0005-1098/$ - see front matter 䉷 2006 Elsevier Ltd. All rights reserved. doi:10.1016/j.automatica.2006.01.010
(Xiu, Lucor, Su, & Karniadakis, 2002). Treatment of other important distributions (e.g., uniform, Poisson) is achieved via the Askey scheme, discussed by Xiu and Karniadakis (2002). Each distribution type relates to a specific family of orthogonal polynomials, so that in many practical cases the optimum convergence properties can be realized. Functionally, we find that the polynomial chaos allows the deterministic parts of a response to be efficiently pre-computed, and the stochastic parts to be injected subsequently, at low computational cost. Hence, the polynomial chaos is an efficient alternative to Monte Carlo simulations for complex systems. The potential of such a tool applied to control engineering is significant, especially considering nonlinear systems, for which limited methods are available for design and analysis of even seemingly simple problems. Indeed, nonlinear systems like the one considered in this paper are stabilized today using feedback linearization where possible, Lyapunov functions where possible, and trial and error tuning with Monte Carlo testing. The first approach may place unrealistic measurement demands on the designer, and the second approach sometimes simply fails. The last approach, of course the most onerous from a computational point of view, often remains as the only viable path to progress in the design. In the case where a controller design has been made through some justifiable method or procedure, it is also desirable to
790
F.S. Hover, M.S. Triantafyllou / Automatica 42 (2006) 789 – 795
understand the distribution of the response in terms of the uncertainties that exist. The probabilistic point of view in both open- and closed-loop response has an intuitive appeal. For example, the manufacturing accuracy of physical components is typically expressed as a standard deviation of a (often normal) distribution. Given that stability is known, how does a specific distribution in the plant parameters map to the range of responses? The same might be said for initial conditions: the initial conditions for which a system must be stabilized sometimes occupy a distribution, with known statistics. The polynomial chaos enables one to directly describe the statistics of the system response. These concepts differ with the established procedure in control theory of defining uncertainty with a uniform probability distribution, i.e., having equally likely values between a hard minimum and a hard maximum. The polynomial chaos expansion appears to be a completely new technique for these problems, and has not been reported in the literature. As we will show, the method is noteworthy for the fact that random variables with large variance can be accommodated. We focus here on two main topics, suitable for general nonlinear systems: 1. Prediction of short-term statistics and stability, given parameters or initial conditions with known distribution. 2. Controller design, enabled by treating the unknown gains as uncertain parameters. The paper is focussed on four examples, highlighting the convergence properties and the computational tradeoffs that occur. We focus on the formulation employing Hermite polynomials with Gaussian distribution, which is not always suitable for stability analysis, since the normal distribution has infinite support. Nonetheless, if the approach is used in a pragmatic way, effective indicators of stability can be constructed. For tighter stability analysis, the uniform distribution and its associated Legendre polynomials are more directly relevant within the classical description of parametric uncertainty. We also note that the case where an uncertain parameter or variable is driven, perhaps through a dynamical system, by a random process, the Karhunen–Loeve spectral decomposition can be incorporated into the polynomial chaos (Lucor, Su, & Karniadakis, 2004; Masri, Smyth, & Traina, 1998). To our knowledge, this is the first time that polynomial chaos methods are applied to control problems. Because feedback control exists largely to limit the effects of uncertainty and disturbances in dynamical systems, however, it would seem to be a natural fit.
of Abramowitz & Stegun, 1972): G=
The Cameron–Martin theorem states that any second-order (i.e., finite variance), one-dimensional process G can be written as a weighted sum of Hermite polynomials in the Gaussian random variable (taken to have zero mean and unity variance), and that this series converges in the L2 sense (using the notation
gi Hei ().
(1)
i=0
Intuitively, the distribution G is constructed from modes which are Hermite polynomials: He0 () = 1, He1 () = , He2 () = 2 − 1, . . . . The Hermite polynomials possess an inner product defined with the Gaussian kernel: ∞ 2 e− /2 Hei ()Hej () d = eij ij , (2) −∞
where ij is the Kronecker delta function. Note that orthogonality holds only for the double product, and that the nonzero inner products are positive integers. Higher products are similarly defined, e.g.,
∞ −∞
e−
2 /2
Hei ()Hej ()Hek () d = eij k .
(3)
As a specific example of how the polynomial chaos can be used, let us consider the following system, which is the basis of the numerical results presented in the next section: x˙ = x + ( + c)y + xy, y˙ = x,
(4)
where , , are known constants and c is a random variable. When = 0, this system if stable is either overdamped or is a second-order oscillator, with natural frequency n = − − c, and damping ratio = −/2n . The regions of stability for the linear uncertain system are of course well understood—they correspond with those cases where + c < 0 and < 0. When = 0, however, this becomes a much more difficult problem, and is typical of a class of relatively low-order polynomial vector fields. This system resists easy analysis by the first or the second methods of Lyapunov, including Krasovski’s theorem. For example, with the Lyapunov function candidate V ( x ) = x 2 /2 + (−c − )y 2 /2, we ob2 ˙ tain V ( x ) = ( + y)x , a suggestive but inconclusive result. Monte Carlo simulation would appear to answer the question of performance and stability, but only indirectly since not all possible initial conditions and uncertainties c can be simulated. In this problem, the random variable c is expressed as a finite summation of modes, up to order P, along with the states: c=
2. The method
∞
P
ci Hei (),
i=0
y(t) =
P
yi (t)Hei (),
i=0
x(t) =
P i=0
xi (t)Hei ().
(5)
F.S. Hover, M.S. Triantafyllou / Automatica 42 (2006) 789 – 795
Monte Carlo 1000
α:
-1
y0(0):
0
800
β:
-4
y1(0):
1
0 0
N:
104
600
γ: c0:
400
c1:
0
200 0
-1
-0.5
0
0.5
mean: σ:
0.001568 0.1724
P:
1
1
Polynomial Chaos 1000 800
791
with the properties of time-invariant linear systems. Fig. 1 shows that the polynomial chaos is successful in generating the correct PDF (shown as a histogram) for the case of y(0) = (that is, y0 (0) = 0, y1 (0) = 1, y2 (0) = 0, . . .), x(0) = 0, no parametric uncertainty, and no nonlinear effects. In this case, only the second mode (He1 ) is nonzero, and so P = 1 is sufficient. We consider here the case of one random dimension only, but the extension to multiple random variables is straightforward, albeit computationally more expensive. Making this restriction, we can have either one random parameter or one random initial condition, but not both. We choose x(t =0)=0 in all calculations below, and consider the cases of random y(t = 0) and random c separately. 3. Numerical examples
600
3.1. Parametric uncertainty in the linear case
400 200 0
-1
-0.5
0
0.5
1
mean:
0
σ:
0.1733
y(t=3) Fig. 1. The case of a linear dynamics, with random initial position y; Mode 1 is the only one excited. Polynomial chaos accurately reconstructs a Gaussian distribution of y(t = 3), as shown in comparison with Monte Carlo simulation having 10,000 trials.
Substituting into the dynamic equation, we obtain P
x˙i (t)Hei () =
i=0
P
xi (t)Hei ()
i=0
+ +
P
ci Hei ()
P
xi (t)Hei ()
i=0 P
y˙i (t)Hei () =
i=0
P
yj (t)Hej ()
j =0
i=0
+
P
P
yj (t)Hej (),
j =0
xi (t)Hei ().
(6)
i=0
Next, employing the inner product, we multiply by the kth mode and perform a Galerkin projection: P P 1 ci yj (t)eij k x˙k (t) = xk (t) + yk (t) + ekk i=0 j =0
P P + xi (t)yj (t)eij k , ekk i=0 j =0
y˙k (t) = xk (t).
(7)
This construction is very similar to that given by Li and Ghanem (1998), for the Duffing oscillator. In the case of no parametric uncertainty (ci = 0, i > 0) and without the nonlinearity ( = 0), each mode evolves independently of the others, in accordance
We first turn to the situation of a linear dynamics, with Gaussian uncertainty in the parameter c. Coupling between the modal states yi (t) and the modes of the unknown parameter c occurs, and hence a higher number of modes has to be considered than in the first case above. As shown in Fig. 2, with P = 12, the polynomial chaos accurately reflects the asymmetric distribution of y(t = 3). The short-term time response is dominated by the first (mean) mode decaying from the initial value; the second mode essentially gives the root variance of the PDF. Near t = 1.7, the variance is small because y(t) is near its peak value. The short-term log plot further shows that the higher modes computed grow progressively smaller by about an order of magnitude per mode. This ordering is typical in the short term. We note, however, that the Hermite polynomials grow rapidly as Hei () = 2Hei−1 () − 2(i − 1)Hei−2 ()
(8)
and that this level of separation in the successive amplitudes is not observed in the present case. As a result, the higher modes will dominate for certain realizations. The capability of the polynomial chaos to accurately produce short-term statistics also applies to unstable responses; this condition was illustrated by Xiu and Karniadakis (2002) for a first-order system. Considering the long-term evolution, decay of P computed modes—to an arbitrarily low-amplitude level—is strongly suggestive of stability because any realization is a summation of modes, and hence decays as well. The time integration in Fig. 2 is performed using adaptive step-size Runge–Kutta routine with relative tolerance 10−6 , and rescaling to maintain accuracy. This example brings up the important issue of finite series truncation error. While the modes for the illustrated case with P = 12 approach zero exponentially in the long term, there do exist values of the random parameter c for which instability occurs, because c is a Gaussian random variable, and therefore has infinite tails. The “spring” term in the governing equation here is + c, or −4 + 0.5 ( = N (0, 1)), and hence the minimum value of which destabilizes the system is eight. The probability of exceeding eight standard deviations above the mean is about 6 × 10−16 ; for most practical purposes, the system is
792
F.S. Hover, M.S. Triantafyllou / Automatica 42 (2006) 789 – 795 Monte Carlo
2000
yi(t)
1
0
0 -0.5
y(t)
2000
0
1000
0
0.5
-4
c1 :
0.5
P:
0
0.5
1
1.5
2
2.5
0 -0.5
3
γ: 0
y1(0): 0
Polynomial Chaos
-0.5 0
seconds
0.5
12
mean: 0.1567 σ: 0.0626
y(t=3)
0
0
0 4
-5
8
-10
log10(abs(yi(t)))
log10(abs(yi(t)))
β:
0
N: 10000 mean: 0.1556 σ: 0.0635
1 0.5
-1
c0:
y0(0): 1
1000
0.5
α:
-5
-10 0 4
-15
-15
8 12
-20 0
0.5
1
1.5
2
2.5
3
seconds
0
10
20
30
40
50
60
seconds
Fig. 2. The case of a random parameter in the linear system, with standard deviation c = 0.5. Top left: 100 Monte Carlo runs, compared with the modal trajectories; only the first few modes are visible. Top right: histograms of the Monte Carlo simulation and from the polynomial chaos. Bottom left: log modal amplitudes for the short-term integration. Bottom right: log modal amplitudes for long-term integration. Mode numbers are listed on right-hand side of the last two graphs.
Table 1 Referencing Fig. 2, the maximum standard deviation of the zero-mean Gaussian variable c that still achieves decay of modes in the long-term, decreases as the number of modes P increases P
max
1 − F (−/)
4 8 12 16
1.31 0.83 0.65 0.55
1.1 × 10−3 7.2 × 10−7 3.8 × 10−10 1.8 × 10−13
The probability of this variable exceeding the critical value (c = 4) for system stability is given in the last column.
stable. Table 1 shows that the maximum value of for which long-term decay of the modes is achieved depends strongly on P. The concern is clearest for the case P = 4, for which decay of the modes holds not only for the nominal case = 0.5, but up to 1.31. Only 4/1.31 = 3.05 standard deviations of this variable will be tolerated by the actual system; this critical level is exceeded in about 0.1 percent of cases, and hence the decay of modes fails as a convincing stability indicator for P = 4. When P is larger, however, we see that accuracy improves, reaching a reasonable instability likelihood of 1.8 × 10−13 with 16 modes. These considerations raise the question of how the polynomial chaos can be practically applied in the general case. First, we note that the normal distribution and the Hermite chaos expansion are specific to this paper, and that different distribution
types with their corresponding orthogonal polynomials (Xiu & Karniadakis, 2002) may be both more accurate and more computationally efficient for suitable problems. At the same time, choosing the best chaos expansion can be difficult, because the uncertainty may have finite support, while the response may not, and vice versa. In these cases especially a high mode number P is required; low-order calculations can be expected to capture major transient effects, but do not represent the distributions well enough to accurately portray long-term trajectories. The situation is not unlike that in perturbation theory, wherein a slowly converging or even divergent series solution may result, yet still provide a useful solution by optimally truncating the series (Bender & Orszag, 1978).
3.2. Random initial condition in the nonlinear case Next we examine the case of a random initial condition and the nonlinear term in the governing equation. The parametric uncertainty in c is replaced with uncertainty in the initial condition y, that is, we set y1 (t = 0) = 1 and yj (t = 0) = 0, for j = 1. Recalling that with V ( x ) = x 2 /2 + (−c − )y 2 /2, V˙ ( x )=(+y)x 2 . V is not a true Lyapunov function, but it does suggest stability when y(t) < − / when > 0, or y(t) < 1 if we set =−1 and =1. Fig. 3 shows that the polynomial chaos accurately predicts the short-term distribution of y(t). Additionally, the system is indicated stable for all initial conditions
F.S. Hover, M.S. Triantafyllou / Automatica 42 (2006) 789 – 795
793
Monte Carlo
2000
yi(t)
1
0 -0.5
y(t)
2
0
-4
c1 :
0
0.5
γ: 1
y1(0): 1
1
Polynomial Chaos
2000
P:
-1
12
1000
-2 0 -0.5
-3 0
0.5
1
1.5
2
2.5
3
0
seconds
0.5
1
mean: 0.03339 σ: 0.149
y(t=3)
0
0
0 4
-5 8
-10 12
-15
log10(abs(yi(t)))
log10(abs(yi(t)))
β:
0
N: 10000 mean: 0.0359 σ: 0.15
1 0
-1
c0:
y0(0): 0
1000
0
α:
-5
-10 0
-15
4 8
-20 0
0.5
1
1.5
2
2.5
3
0
10
20
seconds
30
40
50
60
seconds
Fig. 3. Nonlinear case, with no parameter uncertainty, but uncertainty in the initial condition. Description is the same as for Fig. 2.
in y, subject to x(t =0)=0. This conclusion on stability follows from the same modal decay result described for the linear case. 3.3. Random parameter in the nonlinear case We now consider the nonlinear system with random variable c; for this case, we increased to a value of 5, so as to bring out more of the nonlinear effects, especially those which change the shape of the PDF’s. Fig. 4 shows a similar performance as in the prior cases. The polynomial chaos accurately generates the short-term PDF, and indicates stability of the system response to the initial conditions. The examples shown so far confirm that nonsmooth PDF’s are well described by a Hermite chaos (i.e., using Gaussian bases) of reasonable order, but that the accuracy of the result can rapidly deteriorate if the number of modes is too low. 3.4. Controller design Consider a design situation for nonlinear systems without successful linearization or Lyapunov-based controller design: Monte Carlo techniques are used in an iterative approach, using a scalar metric likely based on several properties of the transient response, as well as stability. The method proceeds via classical function minimization techniques, such as the Nelder–Meade simplex, which admits complex error functions without known derivatives. In this way, the space of potential gain sets is navigated so as to optimize the performance metric. In this context, the polynomial chaos examples above point to control design, and not just analysis, if random parameters
are replaced with available gains. In other words, we could have interpreted c in the example case above as a (proportional) control gain. As we show next, polynomial chaos allows the time integration of the deterministic parts of the response to be completely computed beforehand, so that the metric evaluation requires no integration. Fig. 5 illustrates this for the system under discussion. The function shown is T J= y(t)2 (1 + c2 ) dt, (9) 0
an LQR-like cost where the control penalty is , and the underlining of c denotes a specific realization of the random variable. The evaluation of this cost is made trivial by the polynomial chaos because, using the expansion,
T
J= 0
P
(yi (t)Hei (c))2 (1 + c2 ) dt
i=0
= (1 + c ) 2
P P
Hei (c)Hej (c)
P P
yi (t)yj (t) dt
0
i=0 j =0
= (1 + c2 )
T
Lij (c)Mij .
(10)
i=0 j =0
The definite integrals Mij can be computed independently of c. Then, for a given value of c taken from the distribution, Lij ( c ) is computed from the Hermite polynomials and hence J is given directly. The complete mapping from c to J is made solving differential equations only once—to obtain the set of
794
F.S. Hover, M.S. Triantafyllou / Automatica 42 (2006) 789 – 795 Monte Carlo
1000 1
yi(t)
0
0 -0.05
y(t)
1 0.5
-1
β:
-4
c0:
0
c1:
0.5
y0(0): 1
500
0.5
α:
0
γ: 5
y1(0): 0
N: 10000 mean: 0.005302 σ: 0.02585
0.05
Polynomial Chaos
1000
P:
0
12
500 -0.5 0
1
2
3
4
5
0 -0.05
6
0
seconds
y(t=6)
0
0
0 4
-5 8 12
-10
log10(abs(yi(t)))
log10(abs(yi(t)))
0.05
mean: 0.005422 σ: 0.02579
-5
-10 0 4
-15
8 -15
-20 0
1
2
3
4
5
6
0
10
seconds
20
30
40
50
60
seconds
Fig. 4. Nonlinear case, with parametric uncertainty and fixed initial condition. Description is the same as for Fig. 2. 3 nonlinear, p.c. nonlinear, M.C. linear, p.c. linear, M.C.
2.5
J
2
1.5
1
0.5
-4
-3
-2
-1
0
1
c
Fig. 5. Linear quadratic cost J as a function of gain c using polynomial chaos for the example system with P = 9, T = 4.5, = −1, = −1, = −5, = 0.1. Results are shown for the nonlinear and linear ( = 0) optimizations; results from polynomial chaos and Monte Carlo simulations are overlaid. We use c0 = 0, c1 = 1 and y0 (0) = 1, y1 (0) = 0 in solving the deterministic equations.
yi (t). As in previous cases, we have considered only one random variable, the control gain, and so the initial condition is fixed in this analysis. This is an inconsequential simplification, for if multidimensional polynomial chaos is employed, then the distribution of initial conditions in x and y sets the initial conditions of xi and yi . We see in Fig. 5 that the results considering nonlinear and linear systems are quite different. The nonlinear system has a minimum cost of 1.22 at c = −2.0, while for the linear system, we have a minimum cost of 0.82 at c = −1.0. If this optimum
c from the linear case is employed in the nonlinear system, the cost is 1.34. The mismatch is not unexpected, and motivates design procedures that can take into account the nonlinearity in a quantitative way. The figure shows both the norm computed through the polynomial chaos and the norm via simulation; they are generally in agreement to the fourth significant digit. In this example, we used a low-order expansion (P = 9) to accurately study the short-term system behavior, but the modes are divergent in the long term. This unlikely but still useful condition results from the fundamental tradeoff between stability and accuracy mentioned previously for distributions with infinite support. Accurate description of the system response for random values far from the mean can only occur with high mode number P. Higher P, however, implies also that the tails are more accurately reflected in the response, i.e., the modes are unstable. In the present example, P was fixed and the distribution parameters (specifically c1 ) were chosen so that the range of gains appropriate to the control problem were within several standard deviations, i.e., in a numerically well-behaved neighborhood, of the mean value. The accuracy was checked by making parallel time-domain simulations, and the comparison indicated in the figure is poor more than several deviations away from c’s mean value. 4. Summary We demonstrated that the polynomial chaos method applies to two major areas in nonlinear systems analysis: assessing transient response and long-term stability. The former is illustrated by calculation of explicit probability density functions at
F.S. Hover, M.S. Triantafyllou / Automatica 42 (2006) 789 – 795
a given time, and the latter by decay of the modal amplitudes, for an arbitrarily large number of modes P. In both cases, the uncertainty can occur in either or both of the model parameters and the initial conditions. For random parameters with infinite support, stability can only be investigated in a pragmatic way due to truncation errors, and other chaos expansions may be used when specific stability bounds are the main concern. The construction of statistics requires random samples, but stability is a one-time calculation, namely the temporal evolution of the modes: if they are converging to zero, then stability is indicated. Recasting a parametric uncertainty problem as a gain selection problem allows polynomial chaos to be applied in control system design. Not only can statistics of the response be quickly calculated at any point in time, but so can statistics of functions of the state, including norms, peak value, time to arrive, control expense, and so on. With the particular expansion we chose for this paper, i.e. Hermite/Gaussian, we are left still with the question of long-term behavior for the chosen gain or gains, because design to minimize a short-term cost with a truncated series does not prove stability. A straightforward procedure to subsequently demonstrate stability, however, would be to implement the chosen gains as fixed values, and perform the polynomial chaos analysis again to assess the long-term trajectories.
Acknowledgements Work was supported by the Office of Naval Research’s National Naval Responsibility (NNR) program, grant number N00014-02-1-0623, monitored by Dr. Sharon Curtin-Beerman. We also acknowledge ONR’s Autonomous Operations Future Naval Capability (AOFNC) program, grant N00014-02-C0202, monitored by J. Valentine. References Abramowitz, M., & Stegun, I. A. (1972). Handbook of mathematical tables. New York: Dover Publications. Bender, C. M., & Orszag, S. A. (1978). Advanced mathematical methods for scientists and engineers. New York: McGraw-Hill.
795
Cameron, R. H., & Martin, W. T. (1947). The orthogonal development of non-linear functionals in series of Fourier–Hermite functionals. Annals of Mathematics, 48, 385–392. Ghanem, R., & Spanos, P. (1991). Stochastic finite elements: A spectral approach. New York: Springer. Li, R., & Ghanem, R. (1998). Adaptive polynomial chaos expansions applied to statistics of extremes in nonlinear random vibration. Probability Engineering Mechanics, 13, 125–136. Lucor, D., Su, C.-H., & Karniadakis, G. E. (2004). Generalized polynomial chaos and random oscillators. International Journal for Numerical Methods in Engineering, 60, 571–596. Masri, S. F., Smyth, A. W., & Traina, M. I. (1998). Probabilistic representation and transmission of nonstationary processes in multi-degree-of-freedom systems. Journal of Applied Mechanics, 65, 398–409. Wiener, N. (1938). The homogeneous chaos. American Journal of Mathematics, 60, 897–936. Xiu, D., & Karniadakis, G. Em. (2002). The Wiener–Askey polynomial chaos for stochastic differential equations. SIAM Journal on Scientific Computing, 24, 619–644. Xiu, D., Lucor, D., Su, C.-H., & Karniadakis, G. Em. (2002). Stochastic modelling of flow-structure interactions using generalized polynomial chaos. Journal of Fluids Engineering, 124, 51–59. Franz Hover received his M.S. and Sc.D. degrees from the Woods Hole Oceanographic Institution/Massachusetts Institute of Technology Joint Program in Applied Ocean Physics and Engineering. He is currently a Principal Research Engineer with the Department of Mechanical Engineering at the Massachusetts Institute of Technology, with research interests in the area of applied control and robotics, fluid mechanics, and offshore structures. Dr. Hover is a Member of the IEEE and the ASME. Michael Triantafyllou is a Professor of Mechanical and Ocean Engineering at the Massachusetts Institute of Technology, Director of the Center for Ocean Engineering, and Director of the Testing Tank and Marine Hydrodynamics Laboratories. His research interests are in the dynamics and control of marine vehicles, and biomimetics. He is a Member of the American Physical Society, the Society of Naval Architects and Marine Engineers, and the International Society of Offshore and Polar Engineers.
