Second Order Methods for Optimal Control of Time-Dependent Fluid Flow.
Second order methods for open loop optimal control problems governed by the two-dimensional instationary Navier--Stokes equations are investigated. Optimality systems based on a Lagrangian formulation and adjoint equations are derived. The Newton and quasi-Newton methods as well as various variants of SQP methods are developed for applications to optimal flow control, and their complexity in terms of system solves is discussed. Local convergence and rate of convergence are proved. A numerical example illustrates the feasibility of solving optimal control problems for two-dimensional instationary Navier--Stokes equations by second order numerical methods in a standard workstation environment.
Introduction
This research is devoted to the analysis of second methods for solving optimal
control problems involving the time dependent Navier Stokes equations. Thus we
consider
min J(y; u) over (y; u)
subject to> > <
@y
@
in
Here
is a bounded domain in R 2 , with su-ciently smooth boundary @
The
nal time T > 0 and the initial condition y 0 are xed. The vector valued variable y
and the scalar valued variable p represent the velocity and the pressure of the
uid.
Further u denotes the control variable and B the control operator. The precise
functional analytic setting of problem (1.1), (1.2) will be given in Section 2. For
the moment it su-ces to say that typical cost functionals include tracking type
and functionals involving the vorticity of the
uid
jcurl y(t; )j 2
where  > 0 and z are given. For the following discussion it will be convenient to
formally represent all the equality constraints involved in (1.2) by ^ e(y;
that (1.1), (1.2) can be expressed in the form
min J(y; u) over (y; u)
subject to
In this form solving (1.1), (1.2) appears at rst to be a standard task, see [AM,
SECOND ORDER METHODS IN FLOW CONTROL 3
and the references given there. However, the formidable
size of (1.1), (1.2) and the goal of analyzing second order methods necessitate
an independent analysis.
For second order methods applied to optimal control problems two classes can
be distinguished depending on whether (y; p) in (1.1), (1.2) are considered as independent
variables or as functions of the control variable u. In the former case
represents an explicit constraint for the optimization problem whereas
in the latter case serves the purpose of describing the evaluation
of (y; p) as a function of u. In fact (P ) can be expressed as the reduced problem
where y(u) is implicitly dened via
To obtain a second order method in the case when (y; p) are considered as independent
variables one can derive the optimality system for (P ) and apply the
Newton algorithm to the optimality system. This is referred to as the sequential
quadratic programming (SQP){method. Alternatively, if (y; p) are considered as
functions of u, then Newton's method can be applied to
directly. The relative
merits of these two approaches will be discussed in Section 4. To anticipate some of
this discussion let us point out that the dierence in numerical eort between these
two methods is rather small. In fact, after proper rearrangements the dierence
in computational cost per iteration of the SQP{method for (P ) and the Newton
method for
solving either the linearized equation (1.2) or the full
nonlinear equation itself. In view of the time{dependence of either of these two
equations an iterative procedure is used for their solution so that the dierence
between solving the linearized and nonlinear equation per sweep is not so signi-
cant. A second consideration that may in
uence the choice between SQP{method
or Newton{method applied to
and u 0 for (y; p; u) can clearly be used independently of each other in the SQP{
method, where the states are decoupled from the controls. It is sometimes hinted at
that this decoupling is not only important for the initialization but also during the
iteration and that as a consequence the SQP{method may require fewer iterations
than Newton's method for
As we shall see below, the variables y and p
can be initialized independently from u 0 also in the Newton method. Specically,
are available it is not necessary to abandon (y
As for the choice of the initial guess (y
is to rely on one of the suboptimal strategies that were developed in the
recent past to obtain approximate solutions to (1.1), (1.2). We mention reduced
order techniques [IR], POD{based methods [HK, KV, LT] and the instantaneous
control method [CTMK, BMT, CHK]. As another possibility one can carry out
some gradient steps before one switches to the Newton iteration.
Let us brie
y comment on some related contributions. In [AT] optimality systems
are derived for problems of the type (1.1), (1.2). A gradient technique is proposed
in [GM] for the solution of (1.1), (1.2). Similarly in [B] gradient techniques are
analyzed for a boundary control problem related to (1.1), (1.2). In [FGH] the
authors analyze optimality systems for exterior boundary control problems. One of
the few contributions focusing on second order methods for optimal control of
uids
are given in [GB, H]. These works are restricted to stationary problems, however.
This paper, on the other hand focuses on second order methods for time dependent
problems. We show that despite the di-culties due to the size of (1.1),
(1.2) and the fact that the optimality systems contains a two point boundary value
problem, forward in time for the primal- and backwards in time for the adjoint vari-
ables, second order methods are computationally feasible. We establish that the
initial approximation to the reduced Hessian is only a compact perturbation of the
Hessian at the minimizer. In addition we give conditions for second order su-cient
optimality conditions of tracking type problems. These results imply superlinear
convergence of quasi{Newton as well as SQP{methods. While the present paper
focuses on distributed control problems in a future paper we plan to address the
case of velocity control along the boundary.
The paper is organized as follows. Section 2 contains the necessary analytic
prerequisites. First and second order derivatives of the cost functional with respect
to the control are computed in Section 3. The fourth section contains a comparison
of second order methods to solve (1.1), (1.2). In Section 5 convergence of quasi{
Newton method and SQP{methods applied to
P ) is analyzed. Numerical results
for the Newton{method and comparisons to a gradient method are contained in
Section 6.
2. The optimal control problem
In this section we consider the optimal control problem (1.1), (1.2) in the abstract
subject to e(y;
To dene the spaces and operators arising in (2.1) we
assume
to be a bounded
domain in R 2 with Lipschitz boundary and introduce the solenoidal spaces
with the superscripts denoting closures in the respective norms. Further we dene
W endowed with the norm
H
equipped with the norm
H
denoting the dual space of V . Here
is an abbreviation for L 2 (0;
that up to a set of measure zero in (0; T ) elements can be identied
with elements in C([0; T can be identied with
elements in C([0; T ]; V ). In (2.1) further U denotes the Hilbert space of controls
R is the cost functional which is assumed to be bounded
from below, weakly lower semi-continuous, twice Frechet dierentiable with locally
Lipschitzean second derivative, and radially unbounded in u, i.e.
. Furthermore, the control space U is identied with
SECOND ORDER METHODS IN FLOW CONTROL 5
its dual U  . To simplify the notation for the second derivative we also assume that
the functional J can be decomposed as
The nonlinear mapping
is dened by
Comparing (1.1), (1.2) to (2.1) we note that
the conservation of mass, as well as the boundary condition are realized in the choice
of the space W while the dynamics are described by the condition e(y; In
variational form the constraints in (2.1) can be equivalently expressed as:
given u 2 U nd y 2 W such that
The following existence result for the Navier{Stokes equations in dimension two is
well known ([CF, L, T], Chapter III).
Proposition 2.1. There exists a constant C such that for every u 2 U there exists
a unique element
and
U
From Proposition (2.1) we conclude that with respect to existence (2.1) is equivalent
to
subject to u 2 U;
Theorem 2.1. Problem (2.1) admits a solution (y
Proof. With the above formalism the proof is quite standard and we only give a
brief outline. Since J is bounded from below there exists a minimizing sequence
un g in W  U . Due to the radial unboundedness property
of J in u and Proposition 2.1 the sequence f(y n ; un )g is bounded in W  U and
hence there exists a subsequence with a weak limit (y  ; u  ) 2 W  U . Weak lower
semi-continuity of (y; u) ! J(y; u) implies that
and it remains to show that y This can be achieved by passing to the
limit in (2.3) with (y; u) replaced by (y(u n ); un ).
We shall also require the following result concerning strong solutions to the Navier-Stokes
equation, ([T], Theorem III. 3.10).
Proposition 2.2. If y then for every u 2 U the
solution Moreover, for every
bounded set U in U
is bounded in H 2;1 (Q):
6 MICHAEL HINZE AND KARL KUNISCH
We shall frequently refer to the linearized Navier-Stokes system and the adjoint
equations given next:
in
a.e. on (0; T ];
and
in
a.e. on [0; T
Proposition 2.3. Let y
3 ]. Then (2.5) admits a unique variational solution v 2 W and (2.6) has a
unique variational solution w and the
rst equation in (2.6) holding in L  (V  ) \ W  . Moreover, the following estimates
hold.


iii. jwj L 2

If in addition y
iv. jwj L 2
For
@
solutions v of (2.5) and w of (2.6) are elements of H 2;1 (Q) and satisfy the a-priori
estimates
v.

and
vi.
Proof. The proof is sketched in the Appendix.
3. Derivatives
In this section representations for the rst and second derivatives of ^
J appropriate
for the treatment of (2.4) by the Newton and quasi{Newton method are derived.
We shall utilize the notation
Proposition 3.1. The operator continuously differentiable
with Lipschitz continuous second derivative. The action of the rst two
derivatives of e 1 are given by
he 1
where (v; r) 2 X .
SECOND ORDER METHODS IN FLOW CONTROL 7
Proof. Since e 2 is linear we restrict our attention to e 1 . Let
dened by
and recall that, due to the assumption
that
for all (u; v; To argue local Lipschitz continuity of e,
let
pZ Tjy ~
Here and below C denotes a constant independent of x; ~
x and . Due to the
continuous embedding of W into L 1 (H) we have
jx ~
R T
Using Holder's inequality this further implies the estimate
jx ~
and consequently, after redening C one last time
This estimate establishes the local Lipschitz continuity of e. To verify that the
formula for e x given above represents the Frechet - derivative of e we estimate
R T
sup
R T
Cjy ~
R T
and Frechet - dierentiability of e follows. To show Lipschitz continuity of the rst
derivative let x; ~
x and (v; r) be in X and estimate
R T
R T
Cjy ~
The expression for the second derivative can be veried by an estimate analogous
to the one for the rst derivative. The second derivative is independent of the point
at which it is taken and thus it is necessarily Lipschitz continuous.
From (3.2) it follows that for  2 L 2 (V ) and w 2 W the mapping
is an element of W  . In Section 4 we shall use the fact that  can also be identied
with an element of L 4
Lemma 3.1. For  2 L 2 (V ) and w 2 W the functional  can be identied with
an element in W  \ L 4=3 (V  ).
Proof. To argue that  2 L 4=3 using (3.2)
where k is the embedding constant of V into H . This gives the claim .
Proposition 3.2. Let is a homeo-
morphism. Moreover, if the inverse of its adjoint e
is applied to an
element setting (w; w
we have w and w is the variational solution to (2.6).
Proof. Due to Proposition 3.1, e y (x) is a bounded linear operator. By the closed
range theorem the claim follows once it is argued that (2.5) has a unique solution
. This is a direct consequence of Proposition 2.3, i.
and ii. The assertion concerning the adjoint follows from the same proposition, iii.
and its proof.
As a consequence of Propositions 3.1 and 3.2 and the implicit function theorem
the rst derivative of the mapping u ! y(u) at u in direction -u is given by
y (x)e u (x)-u;
u). By the chain rule we thus obtain
Introducing the variable
we obtain utilizing Proposition 2.3 iii. with
representation for the rst derivative of
Here
is the variational solution
of
where the rst equation holds in L 4=3 (V  ) \ W  .
The computation of the second derivative of ^
J is more involved.
Let (-u; -v) 2 U  U and note that the second derivative of u ! y(u) from U to
W can be expressed as
SECOND ORDER METHODS IN FLOW CONTROL 9
By the chain rule, and since W
We introduce the Lagrangian
and the matrix operator
y (x)e u (x)
We observe that the second derivative of L with respect to x can be expressed as
The above computation for ^
J 00 (u) together with (3.4) imply that
4. Second order methods
This section contains a description and a comparison of second order methods
to solve (2.1). Throughout u  denotes a (local) solution to (2.1).
4.1. Newton{and quasi{Newton algorithm. For the sake of reference let us
specify the Newton algorithm.
Algorithm 4.1. (Newton Algorithm).
1. Choose u
2. Do until convergence
ii) update u
Let us consider the linear system in 2. i). Its dimension is that of the control space
U . From the characterization of the Hessian ^
we conclude that its evaluation
requires as many solutions to the linearized Navier{Stokes equation (3.4) with
appropriate right hand sides as is the dimension of U . If U is innite dimensional
then an appropriate discretization must be carried out. Let us assume now that the
dimension of U is large so that direct evaluation of ^
not feasible. In this
case 2. i) must be solved iteratively, e. g. by a conjugate gradient technique. We
shall then refer to 2. i) as the "inner " loop as opposed to the do{loop in 2. which
is the "outer" loop of the Newton algorithm. The inner loop at iteration level k of
the outer loop requires to
with
(ii) iteratively evaluate the action of ^
j , the j{th iterate of the inner
loop on the k{th level of the outer loop.
The iterate
j can be evaluated by successively applying the steps
a) solve in L 2
b) evaluate J yy (x)v
c) solve in W  for w
d) and nally set q := J uu -u +B  w.
We recall that  1 2 L 2 (V ) and that for s 2 W
he 1
Z TZ
Moreover, by Lemma 3.1 the functional appearing in b) is an element of W  \
Hence by Proposition 2.3 the adjoint equation in c) can equivalently be
rewritten as
where the rst equation holds in W  \ L 4=3 (V  ). Summarizing, for the outer
iteration of the Newton method one Navier{Stokes solve for y(u k ) and one linearized
Navier{Stokes solve for are required. For the inner loop one forward ({in time)
as well as one backwards linearized Navier{Stokes solve per iteration is necessary.
Concerning initialization we observe that if initial guesses (y are
available (with y 0 not necessarily y(u 0 )) then alternatively to the initialization in
Algorithm 4.1 this information can be used advantageously to compute the adjoint
variable  1 required for the initial guess for the right hand side of the linear system
as well as to carry out steps a) - c) for the evaluation of the Hessian. There is no
necessity to recompute y(u 0 ) from u 0 .
To avoid the di-culties of evaluating the action of the exact Hessian in Algorithm
4.1 one can resort to quasi{Newton algorithms. Here we specify one of the most
prominent candidates, the BFGS{method. For w and z in U we dene the rank{one
operator
z 2 L(U ), the action of which is given by
(w
In the BFGS{method the Hessian ^
J 00 at u  is approximated by a sequence of operators
Algorithm 4.2. (BFGS{Algorithm)
1. Choose u
SECOND ORDER METHODS IN FLOW CONTROL 11
2. Do until convergence
ii) update u
Note that the BFGS{algorithm requires no more system solves than the gradient
algorithm applied to (2.1), which is one forward solution of the nonlinear equation
to obtain y(u k ) and one backward solve of the linearized equation (3.7) obtain the
adjoint variable (u k ).
In order to compare Newton's method to the SQP method derived in the next
section we rewrite the update step 2. i) in Algorithm 4.1. To begin with we observe
that the right hand side in the update step can be written with the help of the
adjoint variable  from and the operator T (x) dened in (3.10) as
J
J
where we dropped the iteration indices. As a consequence, with
(3.3) the update can be written as
-y
J
so that
-y
J
holds. Since e x
closed and we have the sequence of identities
x
Thus there exists - 2 Z such that
e
-y
J
Using this equation together with the denition of -y, Newton's update may be
rewritten as
2-y
4.2. Basic SQP{method. Here we regard (2.1) as a minimization problem of the
functional J over the space X subject to the explicit constraint
SQP-algorithm consists in applying Newton's method to the rst order optimality
system
where the Lagrangian L is dened in (3.9).
With x  denoting a solution to problem (P), e x surjective by Proposition
3.2, and hence there exists a Lagrange multiplier   2 Z which is even unique such
that (4.4) holds. The SQP{method will be well dened and locally second order
convergent, if in addition to the surjectivity of e x the following second order
optimality condition holds.
There exists  > 0 such that
If (H1) holds then, due to the regularity properties of e there exists a neighborhood
uniformly positive denite on ker(e x (x)) for every
Algorithm 4.3. (SQP{algorithm)
1. Choose
2. Do until convergence
solve
ii) update
Just as for Newton's method step 2. i.) is the di-cult one. While in contrast
to Newton's method neither the Navier{Stokes equation nor its linearization needs
to be solved, the dimension of the system matrix which is twice the dimension of
the state plus the dimension of the control space is formidable for applications in
uid mechanics. In addition from experience with Algorithm 4.3 for other optimal
control problems, see [KA, V] for example, it is well known that preconditioning
techniques must be applied to solve (4.5) e-ciently. As a preconditioner one might
consider the (action of the) operator
H is the inverse to the (discretized) instationary Stokes operator or
the (discretized) linearization of the Navier{Stokes equation at the state y k , either
one with homogenous boundary conditions.
One iteration of the preconditioned version of Algorithm 4.3 therefore requires
two linear parabolic solves, one forward and one backwards in time. As a con-
sequence, even with the application of preconditioning techniques, the numerical
expense counted in number of parabolic system solves is less for the SQP{method
than for Newton's method. However, the number of iterations of iterative methods
applied to solve the system equations in Algorithms 4.1 and 4.3 strongly depends on
the system dimension, which is much larger for Algorithm 4.3 than for Algorithm
4.1.
To further compare the structure of the Newton and the SQP{methods let us
assume for an instance that x k is feasible for the primal equation, i. e. e(x k
and feasible for the adjoint equation (3.5), i. e. e
SECOND ORDER METHODS IN FLOW CONTROL 13
Then the right hand side of (4.5) has the form@J u
A
and comparing to the computation at the end of section 4.1 we observe that the
linear systems describing the Newton and the SQP{methods coincide. In general
the nonlinear primal and the linearized adjoint equation will not be satised by the
iterates of the SQP{method and we therefore refer to the SQP{method as an outer
or unfeasible method, while the Newton method is a feasible one.
4.3. Reduced SQP{method. The idea of the reduced SQP{method is to replace
(4.5) with an equation in ker e x (x), so that the reduced system is of smaller dimension
than the original one. To develop the reduced system we follow the lines of
[KS]. Recall the denition of T
Note that A is a right{inverse to e x (x). In fact, we have
y (x)e u (x)v
By Proposition 3.2 and due to B 2 L(U; L 2 (V  )) the operator T (x) is an isomorphism
from U to ker e x (x) and hence the second equality in (4.5) given by
can be expressed as
Using this in the rst equality of (4.5) we nd
x
Applying T  (x) to this last equation and ii) from above implies that if -u is a
solution coordinate of (4.5) then it also satises
Once -u is computed from (4.8) then -y and - can be obtained from (4.7) (which
requires one forward linear parabolic solve) and the rst equation in (4.5) (another
backwards linear parabolic solve).
Let us note that if x is feasible then the rst term on the right hand side of (4.8)
is zero and (4.8) is identical to step 2. i) in Newton's Algorithm 4.1.
This again re
ects the fact that Newton's method can be viewed as an SQP{
method that obeys the feasibility constraint It also points at the fact
that the amount of work (measured in equation solves) for the inner loop coincides
for both the Newton and the reduced SQP{methods. The signicant dierence
between the two methods lies in the outer iteration. To make this evident we next
specify the reduced SQP{algorithm.
14 MICHAEL HINZE AND KARL KUNISCH
Algorithm 4.4. (Reduced SQP{algorithm).
1. Choose x
2. Do until convergence
i) Lagrange multiplier update: solve
e
ii) Solve
iii) update
Note that in the algorithm that we specied we did not follow the procedure outlined
above for the update of the Lagrange variable. In fact for reduced SQP{methods
there is no "optimal" update strategy for . The two choices described above are
natural and frequently used. To implement Algorithm 4.4 two linear parabolic
systems have to be solved in steps 2. i) and 2. ii) ) and, in addition two linear
parabolic systems are necessary to evaluate the term involving the operator A on
the right hand side of 2. ii) ). In applications this term is often neglected since it
vanishes at x  .
The reduced SQP{method and Newton's method turn out to be very similar.
Let us discuss the points in which they dier:
Most signicantly the velocity eld is updated by means of the nonlinear
equation in Newton's method and via the linearized equation in the reduced
SQP{method.
ii) The right hand sides of the linear systems dier due to the appearance of the
term involving the operator A. As mentioned above this term is frequently
not implemented.
iii) Formally there is a dierence in the initialization procedure in that y 0 is
chosen independently from u 0 in the reduced SQP{method and y
Newton's method. However, as explained in section 4.1 above, if a good initial
guess y 0 independent from y(u 0 ) is available, it can be utilized in Newton's
method as well.
5. Convergence analysis
We present local convergence results for the algorithms introduced in Section 4
for cost functionals of separable type (2.2). For this purpose it will be essential
to derive conditions that ensure positive deniteness of ^
(H1). The
key to these conditions are the a-priori estimates of Proposition 2.3. We shall also
prove that the dierence ^
compact. This property is required
for the rate of convergence analysis of quasi-Newton methods. In our rst result we
assert positive deniteness of the Hessian provided that J y (x) is su-ciently small,
a condition which is applicable to tracking-type problems.
SECOND ORDER METHODS IN FLOW CONTROL 15
Lemma 5.1. (Positive deniteness of Hessian)
Let u 2 U and assume that J yy positive semi-denite
and J uu (x) 2 L(U) be positive denite, where Then, the Hessian
(u) is positive denite provided that jJ y (x)j L 2 (V  ) is su-ciently small.
Proof: We recall from (3.11) that
is the solution to (3.7). It follows that
e
Here we note that for -u 2 U the functional
is an element of W  . Since J yy (x) is assumed to be positive denite and J uu (x) is
positive denite the result will follow provided the operator norm of
R := e
can be bounded by jJ y (x)j L 2 (V  ) . Straightforward estimation gives
From Proposition 2.3 we conclude that
To estimate
yy (x)(; );  1 (x)ik L(W;W  ) we recall that for
he 1
ZZ
Using (3.2) and the continuity of the embedding W ,! L 1 (H) we may estimate
with a constant C independent of g and h. Therefore,
where we applied iii. in Proposition 2.3 to (3.7).
Lemma 5.2. Let x 2 X and denote by the function dened in (3.5).
Then, under the assumptions of Lemma 5.1 on J condition (H1) is satised with
replaced by (x; ).
Proof. Let (v; u) 2 N (e x (x)). Then v solves (2.5) with
to Proposition 2.3, v 2 W and satises
be chosen such that J uu (x)(u; u)  -juj 2
U for all u 2 U . We nd
pT
Here and below C denotes a generic constant independent of (v; u) and
Due to (3.5) and Proposition 2.3
These estimates imply
and combined with (5.4) the claim follows.
Lemma 5.3. If B 2 L(U; L 2 (H)), then the dierence
is compact for every u 2 U .
Proof. Utilizing (5.2) we may rewrite
It will be shown that both summands dene compact operators
on U . For this purpose let U be a bounded subset of U . Utilizing
Proposition 2.3 we conclude that
y (x)e
is a bounded subset of W and hence of L 2 (V ). Since by assumption J is twice
continuously Frechet dierentiable with respect to y from L 2 (V ) to R it follows
that J yy (S) is a bounded subset of L 2 (V  ). Proposition 2.3, iii. implies that
consequently e
y (J yy (S)) is bounded in W 2
)g. Since W 2
4=3 is compactly embedded in L 2 (H) [CF] and B 2 L(U; L 2 (H))
it follows from the fact that e
that
is pre-compact in U .
Let us turn to the second addend in (5.5). Due to Lemma 3.1 and its proof the
set
is a bounded subset of W  \ L 4=3 (V  ). It follows, utilizing Proposition 2.3 that
is a bounded subset of W 2
4=3 H . As above the assumption that B 2 L(U; L 2 (H))
implies that
SECOND ORDER METHODS IN FLOW CONTROL 17
is precompact in U and the lemma is veried.
The following lemma concerning the operators T (x) and A(x) dened in (3.10) and
(4.6) will be required for the analysis of the reduced SQP-method.
Lemma 5.4. The mappings Let x 7! A(x) from X to L(Z  ; X) and x 7! T (x)
from X to L(U; X) are Frechet dierentiable with Lipschitz continuous derivatives.
Proof. An immediate consequence of i., ii. in Proposition 2.3 and the identities ii)
and iii) in Section (4.3) together with the dierentiability properties of the mapping
x 7! e x (x).
We are now in the position to prove local convergence for the algorithms discussed
in Section 4. Throughout we assume that (y  ; u  ) is a local solution to (2.1)
and set y In addition to the general conditions on J , B
and e we require
positive semi-denite, J uu
positive denite, and jJ y su-ciently small.
With (H2) holding (H1) is satised due to Lemma 5.1. In particular a second
order su-cient optimality condition holds and (y  ; u  ) is a strict local solution to
(2.1). The following theorem follows from well known results on Newton's algorithm

Theorem 5.1. If (H2) holds then there exist a neighbourhood U(u  ) such that
for every u 0 2 U(u  ) the iterates fu n gn2N of Newton's Algorithm 4.1 converge
quadratically to u  .
Theorem 5.2. If (H2) holds then there exist a neighbourhood U(u  ) and  > 0
such that for all u positive denite operators H 0 2 L(U) with
the BFGS method of Algorithm 4.2 converges linearly to u  . If in addition B 2
then the convergence is super-linear.
Proof: For the rst part of the theorem we refer to [GR, Section4], for example.
For the second claim we observe that the dierence ^
by Lemma 5.3, so that the claim follows from [GR, Theorem 5.1], see also [KS1].
Theorem 5.3. Assume that (H2) holds and let   be the Lagrange multiplier
associated to x  . Then there exist a neighbourhood U(x  ;   )  X  Z such that
for all 4.3 is well dened and the iterates
converge quadratically to
Proof: Since J and e are twice continuously dierentiable with Lipschitz continuous
second derivative, e x surjective by Proposition 3.2 and (H1) is satised,
second order convergence of the SQP-method follows from standard results, see for
instance [IK].
We now turn to the reduced SQP-method.
Theorem 5.4. Assume that (H1) holds and let   denote the Lagrange multiplier
associated to x  . Then there exist a neighbourhood U(x  )  X such that for all
reduced SQP-algorithm 4.4 is well dened and its iterates fx n gn2N
converge two-step quadratically to x  , i.e.
for some positive constant C independent of k 2 N.
Proof: First note that (H1) implies positive deniteness of T
in a neighbourhood ~
U(x  ) of x  . By Lemma 5.4 the mappings x 7! T (x) and
x 7! A(x) are Frechet dierentiable with Lipschitz continuous derivatives. Fur-
thermore, utilizing Proposition 2.3, iii. and Lemma A.2 it can be shown that the
mapping x 7! (x) is locally Lipschitz continuous, where  is dened through (3.5).
This, in particular, implies for the Lagrange multiplier updates  k the estimate
where the constant C is positive and depends on x  and on supfjJ yy (x)j L(L 2 (V );L 2 (V
U(x  )g. Altogether, the assumptions for Corollary 3.6 in [K] are met and there
exists a neighbourhood ^
claim follows.
6. Numerical results
Here we present numerical examples that should rst of all demonstrate the
feasibility of utilizing Newton's method for optimal control of the two-dimensional
instationary Navier-Stokes equations in a workstation environment despite the formidable
size of the optimization problem. The total number of unknowns (primal-,
adjoint-, and control variables) in Example 1 below, for instance, is of order 2.210 6 .
The time horizon could still be increased or the mesh size decreased by utilizing
reduced storage techniques at the expense of additional cpu-time, but we shall not
pursue this aspect here. The control problem is given by (1.1), (1.2) with cost
function J dened by
Qo
Qc
with
c
and
subsets
of
denoting the control and observation volumes, respectively. In our examples
Re =400 and B is the indicator function of Q c . The results for
Newton's method will be compared to those of the gradient algorithm, which we
recall here for the sake of convenience.
Algorithm 6.1. Gradient Algorithm
1. choose u 0 ,
2. Set d
d)
3. Set
4. 2.
SECOND ORDER METHODS IN FLOW CONTROL 19
Given a control u the evaluation of the gradient of J at a point u amounts to
solving (1.2) for the state y and (3.7) for the adjoint variable . Implementing
a stepsize rule to determine an approximation of   is numerically expensive as
every evaluation of the functional J at a control u requires solving the instationary
Navier-Stokes equations with right hand side Bu.
We compare two possibilities for computing approximations to the optimal step
size   . For this purpose let us consider for search direction d 2 U the
solutions of the systems
and
is the associated adjoint variable.
1. For a given search direction d 2 U interpolate the function I() by a quadratic
polynomial using the values I(0); I 0 (0) and I 00 (0), i.e.
and use the unique zero
of the equation I 0 as approximation of   , with w given by (6.3).
2. Use the linearization of the mapping  7! y(u + d) at
in the cost functional J . This results in the quadratic approximation
I 2 () := J(y(u) d)
of the functional I(). Now use the unique root
of the equation I 0
as approximation of   , with v given in (6.2).
The denominator of
(u)d; di U . From (5.1) with u  replaced
by u it follows as in the proof of Theorem 5.1 that it is positive, provided that the
state y(u) is su-ciently close to z in L 2 (H).
Let us note that the computation of
1 requires the solution of linearized Navier-Stokes
equations forward and backward in time, whereas that of
only requires
one solve of the linearized Navier-Stokes equations. In addition, a numerical comparison
shows that the step-size guess
performs better than
both with respect
to the number of iterations in the gradient method and with respect to computational
time. For the numerical results presented below we therefore use the step
size proposal
. Thus, every iteration of the gradient algorithm amounts to solving
the nonlinear Navier-Stokes equations forward in time and the associated adjoint
equations backward in time for the computation of the gradient, and to solving
linearized Navier-Stokes equations forward in time for the step size proposal.
The inner iteration of Newton's method is performed by the conjugate gradient
method, the choice of which is justied in a neighbourhood of a local solution u
of the optimal control problem by the positive deniteness of ^
desired state z is su-ciently close to the optimal state y(u  ).
For the numerical tests the target
ow z is given by the Stokes
ow with boundary
condition z tangential direction, see Fig. 1. The termination criterion for
the j-th iterate u k
j in the conjugate gradient method is chosen as
min<

The initialization for Newton's method was u 0 := 0.

Figure

1. Control target, Stokes
ow in the cavity
The discretization of the Navier-Stokes equations, its linearization and adjoint
was carried out by using parts of the code developed by Bansch in [BA], which is
based on Taylor-Hood nite elements for spatial discretization. As time step size
we took which resulted in 160 grid points for the time grid, and 545
pressure and 2113 velocity nodes for the spatial discretization. All computations
were performed on a DEC-ALPHA TM station 500.
Iteration CG-steps
6 19 4.686819e-6 0.032 1.480534e-3

Table

1. Performance of Newton's method for Example 1
Example 1 We rst present the results
for
c

Table

conrms super-linear convergence of the in-exact Newton method. To
SECOND ORDER METHODS IN FLOW CONTROL 21
achieve the the same accuracy as Newton's method the gradient algorithm requires
iterations. The computing time with Newton's method is approximately
minutes whereas the gradient method requires 110 minutes. This demonstrates
the superiority of Newton's method over the gradient algorithm for this example.
For larger values of  and coarser time and space grids the dierence in computing
time is less drastic. In fact this dierence increases with decreasing  and increasing
mesh renement. As expected a signicant amount of computing time is spent for
read-write actions of the variables to the hard-disc in the sub-problems.
In

Figures

2, 3, 4 the evolution of the cost functional, the dierence to the Stokes
ow and the control as a function of time are documented. It can be observed that
Newton's method tends to over-estimate the control in the rst iteration step,
whereas the gradient algorithm approximates the optimal control from below, see

Figure

4. Graphically there is no signicant change after the second iteration for
Newton's method. These comments hold for quite a wide range of values for .
In Fig. 5 the uncontrolled
ow together with the controlled
ow and the control
action at the end of the time interval are presented.2.8e-02

Figure

2. Newton's method (6 Iterations) (top) versus Gradient
algorithm (96 Iterations), Re=400, Evolution of cost
functional for relative
In the previous example the observation
volume
and the control
volume
c
each cover the whole spatial domain. From the practical point of view this is not
feasible. However, from the numerical standpoint this is a complicated situation,
since the inhomogeneities in the primal and adjoint equations are large.
22 MICHAEL HINZE AND KARL KUNISCH1.4e-02

Figure

3. Newton's method (6 Iterations) (top) versus Gradient
algorithm (96 Iterations), Re=400, Evolution of dier-
ence to Stokes
ow for relative
We next present two numerical examples with dierent observation and control
volumes. This results in smaller control and observation volumes than in Example
1, and thus the primal and adjoint equations are numerically simpler to solve.
Example 2
Here
and
1:). The
spatial and temporal discretizations as well as the parameter  are the same as
in Example 1. Newton's method takes 15 minutes cpu-time and its convergence
statistics are presented in Tab. 2. The gradient algorithm needs 25 iterations and
26 minutes cpu to reduce the value of the cost functional from ^
to
Iteration CG-steps

Table

2. Performance of Newton's method for Example 2
Example 3
Here
and
0:7). Again, the
discretization of the spatial and the time domain as well as the parameter  are the
SECOND ORDER METHODS IN FLOW CONTROL 231.8e+00

Figure

4. Newton's method (6 Iterations) (top) versus Gradient
algorithm (96 Iterations), Re=400, Evolution of control
for relative
same as in Example 2. The gradient algorithm needs 38 iterations to reduce the
value of the cost functional from ^
It
takes about 80 minutes cpu-time. We also implemented the Polak-Ribiere variant
of the conjugate gradient algorithm. It converges after 37 iterations and yields
a slightly better reduction of the residual. The amount of of cpu-time needed is
nearly equal to that taken by the gradient algorithm. Newton's method is faster.
It converges within 7 iterations to the approximate solution and needs
cpu-time. The average cpu-time for the inner iteration loop is 7.5 minutes. As in
the previous examples the average cost of a conjugate gradient iteration in the inner
loop decreases with decreasing residual of the outer-iteration loop. The results are
depicted in Tab. 3.


Appendix

A. Proof of Proposition 2.3
In the proof of Proposition 2.3 we make frequent use of the following Lemma.
Lemma A.1. For
Z

Figure

5. Results for from top to bottom: uncontrolled
controlled
ow at control force at
Iteration CG-steps

Table

3. Performance of Newton's method for Example 3
there exists a positive constant C such that
with a positive constant C.
Proof. In [T1].
Lemma A.2. There exists a positive constant C such that
Proof. Since the proof is identical to that of Lemma 3.1.
Note that the power 4=3 in the previous estimate cannot be improved by requiring
Su-cient conditions for (ru) t v +(u  r)v 2 L 2 (V  ) are given by requiring
in addition that u or v 2 L 1 (V ).
Proof of Proposition 2.3. Existence and uniqueness of a solution to (2.5) can be
shown following the lines of the existence and uniqueness proof for the instationary
two-dimensional Navier-Stokes equations in [T, Chap.III]. In the following we sketch
the derivation of the necessary a-priori estimates.
i. Test with respect to time, use b(u; v;
estimate using Young's inequality and the rst estimate in (A.1). This results
in

After integration from 0 to t Gronwall's inequality gives

Using (A.3) in (A.2), the Cauchy-Schwarz inequality yields
26 MICHAEL HINZE AND KARL KUNISCH
Combining (A.3) and (A.4) yields the rst claim.
ii. Test (2.5) with  2 V pointwise in time and estimate using the Cauchy-Schwarz
inequality and the rst estimate in (A.1). This gives
Z
which implies
This, together with y 2 W  L 1 (H), and the estimates (A.3), (A.4) gives ii.
Combining i. and ii. implies

iii. For y 2 W we introduce the bounded linear operator A(y) 2 L(W;Z  ) by
Note that A(y) coincides with e y (x) of Section 3. Due to i. this operator
admits a continuous inverse A(y) 1 2 L(Z  ; W ). For the adjoint A(y)  2
there exists for every
solution (w; to
From i. and ii. together with the fact that
we have
By Lemma A.2 and the assumption that g 2 L 4=3 (V  ) the mapping
is an element of L  (V  ), with  2 [1; 4=3]. From (A.8) we therefore conclude
that w t 2 L  (V  ). Together with w
[DL5]. From (A.8) we deduce that the rst equation in (2.6) is well dened
in L  (V  ). Referring to (A.8) a third time and utilizing the fact that w(T )
is well dened in H it follows that w(T
there exists a constant C such that
Combining this estimate with (A.9) implies the estimate in iii.
iv. If y utilizing
(A.8) we nd that w t 2 L 2 (V  ). Moreover, by (A.1) we have

Together with (A.9) this gives the desired estimate in iv.
v. Test (2.5) with v pointwise in time and utilize Young's inequality and the
last estimate in (A.1) to obtain

Integration from 0 to t together with (A.4) results in
so that Gronwall's inequality gives
Using this in (A.11) yields

To estimate jv t j L 2 (H) test (2.5) with  2 V and use the last estimate in (A.1).
This gives
Z
so that y 2
together with (A.13) and (A.14) implies

Therefore,

which is v. The estimation for jwj L 1
2 ) is similar to that for
In order to cope with b(; in the estimation of
utilizes the third estimate in (A.1) to obtain the estimate vi.
28 MICHAEL HINZE AND KARL KUNISCH



--R

The Lagrange-Newton method for state constrained optimal control problems
On some control problems in uid mechanics

Numerical solution of a ow-control problem: vorticity reduction by dynamic boundary action


Instantaneous control of backward-facing-step ows Preprint No
Feedback control for unsteady ow and its application to the stochastic Burgers equation
Mathematical Analysis and Numerical Methods for Science and Technology

Boundary value problems and optimal boundary control for the Navier-Stokes system: the two-dimensional case
Numerical Methods for Nonlinear Variational Problems
Optimal control of two-and three-dimensional incompressible Navier-Stokes Flows
The local convergence of Broyden-like methods in Lipschitzean problems in Hilbert spaces
The velocity tracking problem for Navier-Stokes ows with bounded distributed controls

Formulation and analysis of a sequential quadratic programming method for the optimal Dirichlet boundary control of Navier-Stokes ow

Control strategies for uid ows - optimal versus suboptimal control
Augmented Lagrangian-SQP-methods for nonlinear optimal control problems of tracking type
Optimal control of thermally convected uid ow

Reduced SQP methods for parameter identi


Mesh independence of the gradient projection method for optimal control problems
Control of Burgers equation by a reduced order approach using Proper Orthogonal Decomposition Bericht Nr.
Mathematical Topics in Fluid Mechanics I
Modelling and control of physical processes using proper orthogonal decomposition




--TR

--CTR
Kerstin Brandes , Roland Griesse, Quantitative stability analysis of optimal solutions in PDE-constrained optimization, Journal of Computational and Applied Mathematics, v.206 n.2, p.908-926, September, 2007
optimal control Navier-Stokes equations;newton method;second order sufficient optimality;SQP method