AM 205: lecture 20 Today: PDE optimization, constrained optimization - - PowerPoint PPT Presentation

AM 205: lecture 20

◮ Today: PDE optimization, constrained optimization example ◮ New topic: eigenvalue problems

PDE Constrained Optimization

We will now consider optimization based on a function that depends on the solution of a PDE Let us denote a parameter dependent PDE as PDE(u(p); p) = 0

◮ p ∈ Rn is a parameter vector; could encode, for example, the

flow speed and direction in a convection–diffusion problem

◮ u(p) is the PDE solution for a given p

PDE Constrained Optimization

We then consider an output functional g,1 which maps an arbitrary function v to R And we introduce a parameter dependent output, G(p) ∈ R, where G(p) ≡ g(u(p)) ∈ R, which we seek to minimize At the end of the day, this gives a standard optimization problem: min

p∈Rn G(p)

1A functional is just a map from a vector space to R

PDE Constrained Optimization

One could equivalently write this PDE-based optimization problem as min

p,u g(u)

subject to PDE(u; p) = 0 For this reason, this type of optimization problem is typically referred to as PDE constrained optimization

◮ objective function g depends on u ◮ u and p are related by the PDE constraint

Based on this formulation, we could introduce Lagrange multipliers and proceed in the usual way for constrained optimization...

PDE Constrained Optimization

Here we will focus on the form we introduced first: min

p∈Rn G(p)

Optimization methods usually need some derivative information, such as using finite differences to approximate ∇G(p)

PDE Constrained Optimization

But using finite differences can be expensive, especially if we have many parameters: ∂G(p) ∂pi ≈ G(p + hei) − G(p) h , hence we need n + 1 evaluations of G to approximate ∇G(p)! We saw from the Himmelblau example that supplying the gradient ∇G(p) cuts down on the number of function evaluations required The extra function calls due to F.D. isn’t a big deal for Himmelblau’s function, each evaluation is very cheap But in PDE constrained optimization, each p → G(p) requires a full PDE solve!

PDE Constrained Optimization

Hence for PDE constrained optimization with many parameters, it is important to be able to compute the gradient more efficiently There are two main approaches:

◮ the direct method ◮ the adjoint method

The direct method is simpler, but the adjoint method is much more efficient if we have many parameters

PDE Output Derivatives

Consider the ODE BVP −u′′(x; p) + r(p)u(x; p) = f (x), u(a) = u(b) = 0 which we will refer to as the primal equation Here p ∈ Rn is the parameter vector, and r : Rn → R We define an output functional based on an integral g(v) ≡ b

a

σ(x)u(x)dx, for some function σ; then G(p) ≡ g(u(p)) ∈ R

The Direct Method

We observe that ∂G(p) ∂pi = b

a

σ(x) ∂u ∂pi dx hence if we can compute ∂u

∂pi , i = 1, 2, . . . , n, then we can obtain

the gradient Assuming sufficient smoothness, we can “differentiate the ODE BVP” wrt pi to obtain, − ∂u ∂pi

′′

(x; p) + r(p) ∂u ∂pi (x; p) = − ∂r ∂pi u(x; p) for i = 1, 2, . . . , n

The Direct Method

Once we compute each ∂u

∂pi we can then evaluate ∇G(p) by

evaluating a sequence of n integrals However, this is not much better than using finite differences: We still need to solve n separate ODE BVPs (Though only the right-hand side changes, so could LU factorize the system matrix once and back/forward sub. for each i)

Adjoint-Based Method

However, a more efficient approach when n is large is the adjoint method We introduce the adjoint equation: −z′′(x; p) + r(p)z(x; p) = σ(x), z(a) = z(b) = 0

Adjoint-Based Method

Now, ∂G(p) ∂pi = b

a

σ(x) ∂u ∂pi dx = b

a

• −z′′(x; p) + r(p)z(x; p)

∂u ∂pi dx = b

a

z(x; p)

• − ∂u

∂pi

′′

(x; p) + r(p) ∂u ∂pi (x; p)

• dx,

where the last line follows by integrating by parts twice (boundary terms vanish because ∂u

∂pi and z are zero at a and b)

(The adjoint equation is defined based on this “integration by parts” relationship to the primal equation)

Adjoint-Based Method

Also, recalling the derivative of the primal problem with respect to pi: − ∂u ∂pi

′′

(x; p) + r(p) ∂u ∂pi (x; p) = − ∂r ∂pi u(x; p), we get ∂G(p) ∂pi = − ∂r ∂pi b

a

z(x; p)u(x; p)dx Therefore, we only need to solve two differential equations (primal and adjoint) to obtain ∇G(p)! For more complicated PDEs the adjoint formulation is more complicated but the basic ideas stay the same

Motivation: Eigenvalue Problems

A matrix A ∈ Cn×n has eigenpairs (λ1, v1), . . . , (λn, vn) ∈ C × Cn such that Avi = λvi, i = 1, 2, . . . , n (We will order the eigenvalues from smallest to largest, so that |λ1| ≤ |λ2| ≤ · · · ≤ |λn|) It is remarkable how important eigenvalues and eigenvectors are in science and engineering!

Motivation: Eigenvalue Problems

For example, eigenvalue problems are closely related to resonance

◮ Pendulums ◮ Natural vibration modes of structures ◮ Musical instruments ◮ Lasers ◮ Nuclear Magnetic Resonance (NMR) ◮ ...

Motivation: Resonance

Consider a spring connected to a mass Suppose that:

◮ the spring satisfies Hooke’s Law,2 F(t) = ky(t) ◮ a periodic forcing, r(t), is applied to the mass

2Here y(t) denotes the position of the mass at time t

Motivation: Resonance

Then Newton’s Second Law gives the ODE y′′(t) + k m

• y(t) = r(t)

where r(t) = F0 cos(ωt) Recall that the solution of this non-homogeneous ODE is the sum

• f a homogeneous solution, yh(t), and a particular solution, yp(t)

Let ω0 ≡

• k/m, then for ω = ω0 we get:3

y(t) = yh(t) + yp(t) = C cos(ω0t − δ) + F0 m(ω2

0 − ω2) cos(ωt),

3C and δ are determined by the ODE initial condition

Motivation: Resonance

The amplitude of yp(t),

F0 m(ω2

0−ω2), as a function of ω is shown

below

0.5 1 1.5 2 2.5 3 3.5 4 −20 −15 −10 −5 5 10 15 20

Clearly we observe singular behavior when the system is forced at its natural frequency, i.e. when ω = ω0

Motivation: Forced Oscillations

Solving the ODE for ω = ω0 gives yp(t) =

F0 2mω0 t sin(ω0t)

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 −5 −4 −3 −2 −1 1 2 3 4 5

t yp(t)

This is resonance!

Motivation: Resonance

In general, ω0 is the frequency at which the unforced system has a non-zero oscillatory solution For the single spring-mass system we substitute4 the oscillatory solution y(t) ≡ xeiω0t into y′′(t) + k

m

• y(t) = 0

This gives a scalar equation for ω0: kx = ω2

0mx =

⇒ ω0 =

• k/m

4Here x is the amplitude

Motivation: Resonance

Suppose now we have a spring-mass system with three masses and three springs

Motivation: Resonance

In the unforced case, this system is governed by the ODE system My′′(t) + Ky(t) = 0, where M is the mass matrix and K is the stiffness matrix M ≡   m1 m2 m3   , K ≡   k1 + k2 −k2 −k2 k2 + k3 −k3 −k3 k3   We again seek a nonzero oscillatory solution to this ODE, i.e. set y(t) = xeiωt, where now y(t) ∈ R3 This gives the algebraic equation Kx = ω2Mx

Motivation: Eigenvalue Problems

Setting A ≡ M−1K and λ = ω2, this gives the eigenvalue problem Ax = λx Here A ∈ R3×3, hence we obtain natural frequencies from the three eigenvalues λ1, λ2, λ3

Motivation: Eigenvalue Problems

The spring-mass systems we have examined so far contain discrete components But the same ideas also apply to continuum models For example, the wave equation models vibration of a string (1D)

• r a drum (2D)

∂2u(x, t) ∂t2 − c∆u(x, t) = 0 As before, we write u(x, t) = ˜ u(x)eiωt, to obtain −∆˜ u(x) = ω2 c ˜ u(x) which is a PDE eigenvalue problem

Motivation: Eigenvalue Problems

We can discretize the Laplacian operator with finite differences to

• btain an algebraic eigenvalue problem

Av = λv, where

◮ the eigenvalue λ = ω2/c gives a natural vibration frequency of

the system

◮ the eigenvector (or eigenmode) v gives the corresponding

vibration mode

Motivation: Eigenvalue Problems

We will use the Python (and Matlab) functions eig and eigs to solve eigenvalue problems

◮ eig: find all eigenvalues/eigenvectors of a dense matrix ◮ eigs: find a few eigenvalues/eigenvectors of a sparse matrix

We will examine the algorithms used by eig and eigs in this Unit