Technische Universit¨at Berlin
Institut f¨ur Mathematik
Numerical solution of optimal control problems with
convex control constraints
Preprint 31-2005
Preprint-Reihe des Instituts f¨ur Mathematik
Technische Universit¨at Berlin
Report 31-2005 November 2005
Numerical solution of optimal control problems
with convex control constraints
Daniel Wachsmuth
Abstract. We study optimal control problems with vector-valued controls.
As model problem serves the optimal distributed control of the instationary
Navier-Stokes equations. In the article, we propose a solution strategy to
solve optimal control problems with pointwise convex control constraints. It
involves a SQP-like step with an imbedded active-set algorithm. The efficiency
of that method is demonstrated in numerical examples and compared to the
primal-dual active-set strategy for box-constraints.
Keywords. Optimal control, convex control constraints, set-valued mappings,
active-set strategy, Navier-Stokes equations.
AMS subject classifications. Primary 49M05, Secondary 26E25, 49K20.
1. Introduction
In this article, we want to investigate solution strategies to solve optimal control
problems with pointwise convex control constraints. In flow control, any control
regardless of distributed or boundary control is a vector-valued function. That
means, there are many possibilities to formulate control constraints. Here, we will
study a general concept of such control constraints. As an example of an optimal
control problem with vector-valued controls we chose the following problem of
optimal distributed control of the instationary Navier-Stokes equations in two
dimensions. To be more specific, we want to minimize the following quadratic
objective functional:
(1) J(y, u) = αT
2ZΩ
|y(x, T)−yT(x)|2dx+αQ
2ZQ
|y(x, t)−yQ(x, t)|2dxdt
+αR
2ZQ
|curl y(x, t)|2dxdt+γ
2ZQ
|u(x, t)|2dxdt
subject to the instationary Navier-Stokes equations
(2)
yt−ν∆y+ (y· ∇)y+∇p=uin Q,
div y= 0 in Q,
y(0) = y0in Ω,
2 D. Wachsmuth
and to the control constraints u∈Uad with set of admissible controls defined by
(3) Uad ={u∈L2(Q)2:u(x, t)∈U(x, t) a.e. on Q}.
Here, Ω is a bounded domain in R2,Qdenotes the time-space cylinder Q:=
Ω×(0, T). The set of admissible controls is generated by the set-valued mapping
U,U:Q7→ 2R2. The conditions imposed on the various ingredients of the optimal
control problem are specified in Sections 2.1 and 2.2, see assumptions (A) and
(AU).
Now, let us discuss several choices of the control constraint U. The so-called
box constraints are often used. There, Uis a rectangle defined by
(4a) U(x, t) = [ua,1(x, t), ub,1(x, t)] ×[ua,2(x, t), ub,2(x, t)]
with given functions ua, ub. For the analysis of optimal control problems for the
non-stationary Navier-Stokes equations using this particular type of control con-
straints, we refer to Hinze and Hinterm¨uller [10], Tr¨oltzsch and Wachsmuth [19],
Ulbrich [20], and Wachsmuth [23].
However, these box constraints are not the only choice for vector-valued con-
trols. For instance, if one wants to bound the R2-norm of the control, one gets a
nonlinear constraint
(4b) |u(x, t)|=pu1(x, t)2+u2(x, t)2≤ρ(x, t),
which gives U(x, t) = {v∈R2:|v| ≤ ρ(x, t)}.
What happens if the control is not allowed to act in all possible directions
but only in directions of a segment with an opening angle less than 180o? Using
polar coordinates ur(x, t) and uφ(x, t) for the control vector u(x, t), this can be
formulated as
(4c) 0 ≤ur(x, t)≤ψ(uφ(x, t), x, t),
where the function ψmodels the shape of the set of allowed control actions.
Here, we will use the constraint in the general form
(4d) u(x, t)∈U(x, t).
All other concepts mentioned above are included in this kind of constraint. One
only has to choose the right U: a box, a circle, or a segment of the plane. We will
impose no further assumptions on Uexcept of the natural ones: non-emptyness,
closedness, convexity, and measurability. We do not need a representation of Uby
inqualities like in (4b) and (4c).
Optimal control problems with such control constraints are rarely investi-
gated in literature. Second-order necessary conditions for problems with the control
constraint u(ξ)∈U(ξ) were proven by P´ales and Zeidan [15] involving second-order
admissible variations. Second-order necessary as well as sufficient conditions were
established in Bonnans [4], Bonnans and Shapiro [5], and Dunn [7]. However, the
set of admissible controls has to be polygonal and independent of ξ, i.e. U(ξ)≡U.
This results were extended by Bonnans and Zidani [6] to the case of finitely many
convex contraints gi(u(ξ)) = 0, i= 1,...,l. A second-order sufficient optimality
Convex control constraints 3
condition for the general case (4d) was proven by the author in [22]. There also
strongly active sets were employed to reduce the subspace on which a coercivity
assumption has to hold. State constraints of the form y(x, t)∈Care considered
in the recent research paper by Griesse and Reyes [8].
The plan of the article is as follows. In Section 2 we collect results concerning
the solvability of the state equation and the optimal control problem. Necessary
optimality conditions are stated in Section 3. Section 4 is devoted to the derivation
of a solution strategy. There, a new active-set algorithm is described. Numerical
experiments that confirm the efficiency of the proposed algorithm are presented
in Section 5.
2. Statement of the optimal control problem
At first, we introduce some notations and results that we will need later on. To
begin with, we define the spaces of solenoidal or divergence-free functions H:=
{v∈L2(Ω)2: div v= 0},V:= {v∈H1
0(Ω)2: div v= 0}.These spaces are
Hilbert spaces with scalar products (·,·)Hand (·,·)Vrespectively. The dual of V
with respect to the scalar product of Hwe denote by V0.
We will work with the standard spaces of abstract functions from [0, T] to
a real Banach space X,Lp(0, T ;X), endowed with their natural norms. In the
sequel, we will identify the spaces Lp(0, T;Lp(Ω)2) and Lp(Q)2for 1 <p<∞,
and denote their norm by kukp:= |u|Lp(Q)2.The usual L2(Q)2-scalar product we
denote by (·,·)Qto avoid ambiguity.
To deal with the time derivative in (2), we introduce the common spaces of
functions ywhose time derivatives ytexist as abstract functions,
Wα(0, T;V) := {y∈L2(0, T;V) : yt∈Lα(0, T;V0)}, W(0, T) := W2(0, T;V),
where 1 ≤α≤2. Endowed with the norm
kykWα(0,T ;V):= kykL2(V)+kytkLα(V0),
these spaces are Banach spaces. Every function of W(0, T) is, up to changes on
sets of zero measure, equivalent to a function of C([0, T ], H), and the imbedding
W(0, T),→C([0, T], H) is continuous, cf. [13].
2.1. The state equation
Before we start with the discussion of the state equation, we specify the require-
ments for the various ingredients describing the optimal control problem. In the
sequel, we assume that the following conditions are satisfied:
(A)
1. Ω ⊂R2is domain with Lipschitz boundary Γ := ∂Ω, such
that Ωis locally on one side of Γ,
2. y0, yT∈H,yQ∈L2(Q)2,
3. αT, αQ, αR≥0,
4. γ, ν > 0.
4 D. Wachsmuth
The assumptions on the set-valued mapping Uare given in the next section. Now,
we will briefly summarize known facts about the solvability of the instationary
Navier-Stokes equations (2). To specify the problem setting, we introduce a linear
operator A:L2(0, T;V)7→ L2(0, T;V0) by
ZT
0
h(Ay)(t), v(t)iV0,V dt:= ZT
0
(y(t), v(t))Vdt,
and a nonlinear operator Bby
ZT
0B(y)(t), v(t)V0,V dt:= ZT
0ZΩ
2
X
i,j=1
yi(t)∂yj(t)
∂xi
vj(t) dxdt.
For instance, the operator Bis continuous and twice Frech´et-differentiable as op-
erator from W(0, T ) to L2(0, T;V0).
Now, we introduce the notation of weak solutions for the instationary Navier-
Stokes equations (2) in the Hilbert space setting.
Definition 2.1 (Weak solution).Let f∈L2(0, T;V0)and y0∈Hbe given. A
function y∈L2(0, T;V)with yt∈L2(0, T;V0)is called weak solution of (2) if
(5) yt+νAy +B(y) = f,
y(0) = y0.
Results concerning the solvability of (5) are standard, cf. [17] for proofs and
further details.
Theorem 2.2 (Existence and uniqueness of solutions).For every source term f∈
L2(0, T;V0)and initial value y0∈H, the equation (5) has a unique solution
y∈W(0, T). Moreover, the mapping (y0, f)7→ yis locally Lipschitz continuous
from H×L2(0, T ;V0)to W(0, T).
2.2. Set-valued functions
Before we begin with the formulation of the optimal control problem with inclusion
constraints, we will specify the notation and assumptions for the admissible set
U(·). It is itself a mapping from the control domain Qto the set of subsets of R2,
it is a so-called set-valued mapping. We will use the notation U:Q R2.
The controls are taken from the space L2(Q)2, so it is natural to require
the fulfillment of (4d) for (only) almost all (x, t)∈Q. And we have to impose at
least some measurability conditions on the mapping U. In the sequel, we will work
with measurable set-valued mappings. For an excellent — and for our purposes
complete — introduction we refer to the textbook by Aubin and Frankowska [2].
Convex control constraints 5
Once and for all, we specify the requirements for the function U, which defines
the control constraints.
(AU)
The set-valued function U:Q R2satisfies:
1. Uis a measurable set-valued function.
2. The images of Uare non-empty, closed, and convex a.e.
on Q. That is, the sets U(x, t)are non-empty, closed and
convex for almost all (x, t)∈Q.
3. There exists a function fU∈L2(Q)2with fU(x, t)∈
U(x, t)a.e. on Q.
Please note, we did not impose any conditions on the sets U(x, t) that are beyond
convexity such as boundedness or regularity of the boundaries ∂U(x, t). Assump-
tions (i) and (ii) guarantee that there exists a measurable selection of U, i.e. a
measurable single-valued function fMwith fM(x, t)∈U(x, t) a.e. on Q. The ex-
istence of a square integrable, admissible function is then ensured by the third
assumption. This implies that the set of admissible control is non-empty. Further-
more, it is also convex and closed in L2(Q)2.
By conditon (iii) it follows also that the pointwise projection on Uad of a
L2-function is itself a L2-function, see e.g. [22].
Corollary 2.3. Let be given a function u∈L2(Q)2. Then the function vdefined
pointwise a.e. by
v(x, t) = ProjU(x,t)(u(x, t))
is also in L2(Q)2. Further, if for some p≥2the functions uand fUare in Lp(Q)2,
then the projection vis in Lp(Q)2as well.
The assumption (AU) is as general as possible. In the case that the set-
valued function Uis a constant function, i.e. U(x, t)≡U0, we can give a simpler
characterization: it suffices that U0is non-empty, closed, and convex.
2.3. Existence of optimal controls
Before we can think about existence of solution, we have to specify which problem
we want to solve. We will assume that conditions (A) of Section 2.1 are satified.
Moreover, we assume that U(·) fulfills the pre-requisite (AU). So we end up with
the following optimization problem
(6a) min J(y, u)
subject to the state equation
(6b) yt+νAy +B(y) = uin L2(0, T;V0),
y(0) = y0,
and the control constraint
(6c) u∈Uad,
where Uad is given by (3).
6 D. Wachsmuth
Under the assumptions above, the optimal control problem (6) is solvable.
We recall that in Section 2.1 the regularization parameter γis supposed to be
greater than zero. One can prove existence even with γ= 0 under the additional
condition of boundedness of Uad in L2.
Theorem 2.4. The optimal control problem admits a - global optimal - solution
¯u∈Uad with associated state ¯y∈W(0, T ).
3. First-order necessary conditions
The necessary optimality conditions for the optimal control problem discussed in
the present article differ less from the conditions that can be found in the literature.
However, we will repeat the exact statement for convenience of the reader.
Theorem 3.1 (Necessary condition).Let ¯ube locally optimal in L2(Q)2with as-
sociated state ¯y=y(¯u). Then there exists a unique Lagrange multiplier ¯
λ∈
W4/3(0, T;V), which is the weak solution of the adjoint equation
(7) −¯
λt+νA¯
λ+B0(¯y)∗¯
λ=αQ(¯y−yQ) + αRcurl∗curl ¯y
¯
λ(T) = αT(¯y(T)−yT).
Moreover, the variational inequality
(8) (γ¯u+¯
λ, u −¯u)Q≥0∀u∈Uad
is satisfied.
Similar as in the box-constrained case, we can reformulate the variational
inequality (8). The projection representation of the optimal control is now realized
using the admissible sets U(·)
(9) ¯u(x, t) = ProjU(x,t)−1
γ¯
λ(x, t)a.e. on Q.
Secondly, another formulation of the variational inequality uses the normal cone
NUad (¯u). The normal cone admits a pointwise representation as Uad itself, see
e.g. [2, 16]:
NUad (u) = v∈L2(Q)2:v(x, t)∈ NU(x,t)(u(x, t)) a.e. on Q.
Then inequality (8) can be written equivalently as the inclusion
(10) ν¯u+¯
λ+NUad (¯u)30.
It will allow us to write the optimality system as a generalized equation.
4. Solution strategy
Now, we are going to describe our strategy to solve optimal control problems with
pointwise convex control constraints. The starting point of our considerations is
Newtons method applied to generalized equations.
Convex control constraints 7
4.1. Generalized Newton’s method
In the sequel, we want to apply Newton’s method to the optimality system. This
system can be written as a generalized equation. Let us define a function Fby
(11) F(y, u, λ) =
yt+νAy +B(y)
y(0)
−λt+νAλ +B0(y)∗λ
λ(T)
γu +λ
−
u
y0
αQ(y−yQ) + αR~
curl curl y
αT(y(T)−yT)
0
.
Then a triple (¯y, ¯u, ¯
λ) fulfills the necessary optimality conditions consisting of the
state equation (6b), the adjoint equation (7), and the variational inequality (8) if
and only if it fulfills the generalized equation
(12) F(¯y, ¯u, ¯
λ) + 0,0,0,0,NUad (¯u)T30.
Now, we will apply the generalized Newton method to equation (12). If the control
constraints are box constraints then this method is equivalent to the SQP-method,
see e.g. [18]. Given iterates (yn, un, λn) we have to solve the linearized generalized
equation
(13) F(yn, un, λn) + F0(yn, un, λn)[(y−yn, u −un, λ −λn)]
+0,0,0,0,NUad (¯u)T30
in every step. It turns out that this equation is the optimality system of the linear-
quadratic optimal control problem.
(14a) min Jn(y, u) = αT
2|y(T)−yd|2
H+αQ
2ky−yQk2
2+αR
2kcurl yk2
2+γ
2kuk2
2
−bQ(y−yn, y −yn, λn)
subject to the linearized state equation
(14b) yt+νAy +B0(yn)y=u+B(yn)
y(0) = y0
and the control constraint
(14c) u∈Uad.
Let us emphasize the following observation. In the generalized equation (12) only
N(x) represents the control constraint. And it is the only term that was not lin-
earized in (13). Consequently, the subproblem (14) is subject to the same control
constraint u(x, t)∈U(x, t) as the original non-linear problem.
That means, the control constraint is not linearized, even if it is written as an
inequality like u1(x, t)2+u2(x, t)2≤ρ(x, t)2. This is a difference to the standard
SQP-method for optimization problems with nonlinear inequalities. An inequality
g(x)≤0 would be linearized to g(xn) + g0(xn)(x−xn)≤0.
8 D. Wachsmuth
Now, it remains to explain how to solve the optimization problem (14). The
main difficulty here, is the convex control constraint (14c).
4.2. Active-set strategy
As for the box-constrained case we will use an active set algorithm. It is very
similar to the primal-dual active-set strategy introduced by [3]. The algorithm we
propose here tries to solve the projection representation of optimal controls
¯u(x, t) = ProjU(x,t)−1
γ¯
λ(x, t)a.e. on Q.
The algorithm to solve the subproblem (14) works as follows. We denote the control
iterates of step kby uk. The state ykis the solution of (14b) with right-hand side
uk, and λkis the solution of the adjoint equation associated to (14b).
Algorithm 4.1. Take a starting guess u0with associated state y0and adjoint λ0.
Set k= 0.
1. Given uk, yk, λk. Determine the active set Ak+1 =Ak+1(λk)by
Ak+1 := (x, t) : −1
γλk(x, t)6∈ U(x, t).
2. Minimize the functional Jngiven by (14a) subject to the linearized state
equation (14b) and to the control constraints
u|Ak+1 = ProjU(x,t)−1
γλk(x, t)
u|Q\Ak+1 free.)
Denote the solution by (˜u, ˜y, ˜
λ).
3. Project uk+1 := ProjU(x,t)−1
γ˜
λ(x, t), compute yk+1 and λk+1.
4. If
uk+1 −ProjUad −1
γλk+1
2< ε →ready,
else set k:= k+ 1 and go back to 1.
In the first step of the algorithm, the active set is determined. The control
constraint at a particular point (x, t) is considered active if −1/γ λk(x, t) does not
belong to U(x, t). In the second step, a linear-quadratic optimization problem is
solved. It involves no inequality constraints, since the control is fixed on the active
set and the control is free on the inactive set. After that, the optimal adjoint state
of that problem is projected on the admissible set to get the new control. The
algorithm stops if the residual in the projection representation is small enough.
The algorithm is very similar to the primal-dual method for box-constrained
problems. But there are some fine differences. Our algorithm uses only one ac-
tive set A. The primal-dual method exploits the fact that the box-constraints are
formed by independent inequalities and therefore it uses active sets Ai,i= 1,2
for each inequality.
Let us explain the behaviour of both algorithms in the detection of the ac-
tive sets for the simple box constraint |ui| ≤ 1, i= 1,2. Let us assume that
Convex control constraints 9
−1/γ λk(x0, t0) = (2,0), i.e. −1/γ λk(x0, t0) is not admissible. Then in our method
the point (x0, t0) will belong to the active set Ak+1 in the next step. And at this
point the value (1,0) is prescribed for uk+1:uk+1(x0, t0) = (1,0). In the primal-
dual method, the point (x0, t0) will be added to the active set Ak+1
1associated to
the inequality |u1| ≤ 1. It will not belong to the active set Ak+1
2for the second
inequality |u2| ≤ 1 since this inequality was not violated. And for the inner prob-
lem we will get the constraints uk+1
1(x0, t0) = 1 and uk+1
2(x0, t0) = free, that is
the control is allowed to vary in tangential directions on the right side of the box!
Furthermore, in the primal-dual active set method the new control iterate
uk+1 was taken as the solution of the inner problem without projecting it. Then
the lower constraint was taken to be active, if uk(x, t) + 1
cmk(x, t)< ua(x, t) with
the multiplier of the control constraint mk(x, t) = −(γuk(x, t) + λk(x, t)). This
works well, since in the box-constrained case there are only two possibilities of
an active control constraint: upper and lower. Now, for convex constraints this is
completely different. If we know that the constraint is active, what value should
we prescribe on the active set? Here, we cannot take simply uaor ub. We have
to compute a point on the boundary of the active set U(x, t). To this end, the
projection of −1/γλ is taken as the value of uon the active set.
The primal-dual method for linear-quadratic optimal control problems with
box constraints is known to be equivalent the semi-smooth Newton method [11]. If
we apply that method to the non-smooth equation u= ProjUad −1
γλ, we would
end up with a different algorithm. As for the primal-dual active-set strategy, we
have to allow tangential variations of the control on the active sets.
5. Numerical results
Now, let us report about the performance of the proposed algorithm. We present
results for an optimal control problem, where the solution is known, to study the
convergence speed and approximation order.
5.1. Problem setting
The computational domain was chosen to be the unit square Ω = (0,1)2with final
time T= 1. We want to minimize the functional
1
2Z1
0ZΩ
|y(x, t)−yd(x, t)|2dxdt+γ
2Z1
0ZΩ
|u(x, t)|2dxdt
subject to the instationary Navier-Stokes equations on Ω ×(0,1) with distributed
control
yt−ν∆y+ (y· ∇)y+∇p=u+fin Q,
div y= 0 in Q,
y= 0 on Σ,
y(0) = y0in Ω.
10 D. Wachsmuth
and subject to some control constraints to be specified later on. We will impose box
as well as convex constraints to compare the behaviour of numerical algorithms.
Let us construct a triple of state, control and adjoint, that satisfies the first-
order optimality system. We chose as state and adjoint state
¯y(x, t) = e−νt sin2(πx1) sin(πx2) cos(πx2)
−sin2(πx2) sin(πx1) cos(πx1)
and
¯
λ(x, t) = e−νt −e−νsin2(πx1) sin(πx2) cos(πx2)
−sin2(πx2) sin(πx1) cos(πx1).
Regardless of the choice of Uad, the control is computed using the projection
formula as
¯u= ProjUad −1
γ¯
λ(x, t).
All other quantities are now chosen in such a way that ¯yand ¯
λare the solutions
of the state and adjoint equations respectively:
f= ¯yt−ν∆¯y+ (¯y· ∇)¯y−¯u,
y0= ¯y(0),
yd= ¯y−−¯
λt−ν∆¯
λ−(¯y· ∇)λ+ (∇¯y)Tλ.
5.2. Discretization
The continuous problem was discretized using Taylor-Hood finite elements with
different mesh sizes. The coarsest grid consists of 256 triangles with 545 velocity
and 145 pressure nodes yielding a mesh size h0= 0.125. Further, we use the semi-
implicit Euler scheme for time integration with a equidistant time discretization
with different step lengths, where in the coarsest discretization we set τ0= 0.01.
We computed solutions of the optimal control problem for different spatial
and time discretizations. The coarse grid was refined uniformly, i.e. each triangle
was divided into four congruent triangles. It leads to a reduction of the mesh size
from h0to h1:= h0/2, and h2:= h0/4. The time step was shortened from τ0= 0.01
to τ1:= τ0/4 = 0.0025, and τ2:= τ0/16 = 0.000625 to get a uniform reduction of
the approximation errors connected with the spatial and time discretization. See
also Table 1, where all these values are summarized.
triangles velocity pressure Mesh size Time step
nodes nodes h τ
Coarse 256 545 145 0.125 0.01
1024 2113 545 0.0625 0.0025
Fine 4096 8321 2113 0.03125 0.000625
Table 1. Discretization parameters
Convex control constraints 11
The control is approximated by piecewise continuous functions in time. For
the spatial discretization, we used piecewise constant functions, too. Here, one can
expect an approximation of the optimal control with order hin the L2as well as
L∞-norm, see Arada, Casas and Tr¨oltzsch [1]. We also investigated the effects of
a post-processing step due to Meyer and R¨osch [14]: the obtained adjoint state ¯
λh
is used in a projection step to get a better approximation of the optimal control
by wh= ProjUad (−1/γ¯
λh).
We used the SQP-method without any globalization to solve the problem with
box constraints. The constrained SQP-subproblems were solved by the primal-dual
active-set method using the method of conjugate gradients (CG) for the inner loop.
Since those subproblems are linear-quadratic optimization problems, this active-set
strategy can be interpreted as a semi-smooth Newton method, see Hinterm¨uller,
Ito, and Kunisch [11], to solve the non-smooth equation u= ProjUad −1
γλ(u),
compare (9). Here, λ(u) denotes the adjoint state for a given control uof the sub-
problem (14). This method is known to converge locally with super-linear conver-
gence rate [11, 20, 21] if the quadratic form L00 is coercive i.e. a sufficient optimality
condition holds. Under some strong assumptions it converges even globally [12].
To solve the optimal control problem with convex constraints, we employ the
solution strategy proposed in the previous section: generalized Newton method
with active-set strategy as inner loop. The quadratic subproblems of the active-set
method were again solved by the CG algorithm.
In all examples, the stopping criteria of the nested methods are balanced in
the following way. The outer method was stopped if the norm of the difference of
two successive iterates was much smaller than the norm of the difference of the
actual iterate to the solution, i.e. if
kyn−yn−1k∞+kun−un−1k∞+kλn−λn−1k∞
≤εEX kyn−¯yk∞+kun−¯uk∞+kλn−¯
λk∞
was satisfied. We used εEX = 10−1in the computations.
The primal-dual active set method was stopped if either the active sets of
two successive control iterates coincide or the error in the variational inequality
given by
φ(u) =
u−ProjUad −1
γλ
2
is reduced by a factor of 0.1. The proposed active-set strategy was stopped if
the latter reduction criterion was fulfilled. In this method, it is not appropriate to
compare active sets, since the values of the control on the active set are not uniquely
determined as in the box-constrained case. The innermost iteration procedure —
the CG method — was stopped if the norm of the residual was reduced by a factor
of 0.01. The initial guesses u0and y0for control and state were set to zero in all
computations.
12 D. Wachsmuth
5.3. Results
Now, let us report about the results for convex control constraints compared to
the box-constraints considered above. The parameters of the example are set to
γ= 0.01 and ν= 0.1. We computed solutions for the box constraint
|ui(x, t)| ≤ 1.0
and for the convex constraint in polar coordinates
0≤ur(x, t)≤ψ(uφ(x, t)).
The function ψ(φ) is given as the spline interpolation with cubic splines and peri-
odic boundary conditons of the function values in Table 2. The admissible control
set is drawn on the right of Table 2. It is chosen such that the admissible set
is comparable to the box constrained case. The projection on that set was com-
puted by Newtons method. That means the computation of the projection is more
costly than in the box-constraint case. Hovever, this additional computation time
is negligible compared with one forward solve of the partial differential equation.
φ0π
4
π
2
3
4π π 5
4π3
2π7
4π
ψ(φ) 1.0 1.0 0.64 0.44 0.4 0.44 0.64 1.0
Table 2. Convex constraint
The effect of the different shapes of the admissible set can be seen in Figure 1.
There is a large part of the control constraints active. The convex constraint (right)
gives a smoother control although the constraint is active on a large part of the
domain. This is an advantage of using convex and smoother admissible sets than
the box.
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Figure 1. Control snapshots at t= 0.05
Figure 2 shows the errors k¯uh−¯uk2(4) and k¯uh−¯uk∞(5) for the the same
optimization problem with box constraints and convex constraints respectively.
Convex control constraints 13
Here, the abscissae show the number of the discretizations, where ’1’ corresponds
to the coarse and ’3’ to the finest discretization. The axis of ordinates corresponds
to the value of the error norms. Auxiliary lines are plotted that visualizes the rates
h+τand h+τ2. As one can see, the approximation order is not affected by the
type of the constraint. The same holds for the post-processed control, see Figure 3.
1 1.5 2 2.5 3 3.5 4
10−4
10−3
10−2
10−1
100
τ+h
τ+h2
||uh−u||2
||uh−u||∞
1 1.5 2 2.5 3 3.5 4
10−4
10−3
10−2
10−1
100
τ+h
τ+h2
||uh−u||2
||uh−u||∞
Figure 2. Difference to solution: box constraints, convex constraints
1 1.5 2 2.5 3 3.5 4
10−4
10−3
10−2
10−1
100
τ+h
τ+h2
||wh−u||2
||wh−u||∞
1 1.5 2 2.5 3 3.5 4
10−4
10−3
10−2
10−1
100
τ+h
τ+h2
||wh−u||2
||wh−u||∞
Figure 3. Postprocessed control: box constraints, convex constraints
Let us compare briefly the convergence of the active set algorithms: the
primal-dual active set method for the box constrained problem and the active set
method proposed above to solve the convex constrained problem. In all our com-
putations both methods showed a similar behaviour, which can be seen in Table 3.
Although they solved optimal control problems with different control constraints,
they needed almost the same number of outer (SQP/generalized Newton) and in-
ner (active-set) iterations. Also the residual u−Proj(−1/γ λ) depicted by ’Res’
decreased in the same way. So, we can say that our algorithm is as efficient as the
well-known primal dual active set method. There is also hope to give convergence
proofs for our method analogous to the proofs in [9, 12].
References
[1] N. Arada, E. Casas, and F. Tr¨oltzsch. Error estimates for a semilinear elliptic control
problem. Computational Optimization and Application, 23:201–229, 2002.
14 D. Wachsmuth
Box constraints Convex constraints
SQP PD-AS gNewton AS
It kun−un−1k∞Res It kun−un−1k∞Res
1 0.1564 1 0.9912 ·10−1
0.4823 ·10−10.2327 ·10−1
0.1128 ·10−10.3475 ·10−2
2 0.1128 ·10−12 0.3498 ·10−2
0.1221 0.4971 ·10−30.1311 0.1976 ·10−3
3 0.6967 ·10−33 0.2530 ·10−3
0.1148 ·10−10.1561 ·10−50.3741 ·10−4
0.1186 ·10−10.7047 ·10−5
4 0.1891 ·10−54 0.7035 ·10−5
0.4105 ·10−60.4643 ·10−60.1291 ·10−5
0.2804 ·10−30.2308 ·10−6
Table 3. Comparison of the algorithms
[2] J.-P. Aubin and H. Frankowska. Set-valued analysis. Birkh¨auser, Boston, 1990.
[3] M. Bergounioux, K. Ito, and K. Kunisch. Primal-dual strategy for constrained opti-
mal control problems. SIAM J. Control Optim., 37:1176–1194, 1999.
[4] J. F. Bonnans. Second-order analysis for control constrained optimal control prob-
lems of semilinear elliptic equations. Appl. Math. Optim., 38:303–325, 1998.
[5] J. F. Bonnans and A. Shapiro. Perturbation analysis of optimization problems.
Springer, New York, 2000.
[6] J. F. Bonnans and H. Zidani. Optimal control problems with partially polyhedric
constraints. SIAM J. Control Optim., 37:1726–1741, 1999.
[7] J. C. Dunn. Second-order optimality conditions in sets of L∞functions with range
in a polyhedron. SIAM J. Control Optim., 33(5):1603–1635, 1995.
[8] R. Griesse and J.C. de los Reyes. State-constrained optimal control of the stationary
Navier-Stokes equations. submitted, 2005.
[9] M. Hinterm¨uller. A primal-dual active set algorithm for bilaterally control con-
strained optimal control problems. Quarterly of Applied Mathematics, LXI 1:131–
161, 2003.
[10] M. Hinterm¨uller and M. Hinze. A SQP-semi-smooth Newton-type algorithm applied
to control of the instationary Navier-Stokes system subject to control constraints.
SIAM J. Optim. To appear.
[11] M. Hinterm¨uller, K. Ito, and K. Kunisch. The primal-dual active set strategy as a
semismooth Newton method. SIAM J. Optim., 13:865–888, 2003.
[12] K. Kunisch and A. R¨osch. Primal-dual active set strategy for a general class of
constrained optimal control problems. SIAM J. Optim., 13:321–334, 2002.
Convex control constraints 15
[13] J. L. Lions and E. Magenes. Non-homogeneous boundary value problems and appli-
cations, volume I. Springer, Berlin, 1972.
[14] C. Meyer and A. R¨osch. Superconvergence properties of optimal control problems.
SIAM J. Control Optim., 43(3):970–985, 2004.
[15] Zs. P´ales and V. Zeidan. Optimum problems with measurable set-valued constraints.
SIAM J. Optim., 11:426–443, 2000.
[16] R. T. Rockafellar. Conjugate duality and optimization. SIAM, Philadelphia, 1974.
[17] R. Temam. Navier-Stokes equations. North Holland, Amsterdam, 1979.
[18] F. Tr¨oltzsch. On the Lagrange-Newton-SQP method for the optimal control of semi-
linear parabolic equations. SIAM J. Control Optim., 38:294–312, 1999.
[19] F. Tr¨oltzsch and D. Wachsmuth. Second-order sufficient optimality conditions for
the optimal control of Navier-Stokes equations. ESAIM: COCV, 2005. To appear.
[20] M. Ulbrich. Constrained optimal control of Navier-Stokes flow by semismooth New-
ton methods. Systems & Control Letters, 48:297–311, 2003.
[21] M. Ulbrich. Semismooth Newton methods for operator equations in function spaces.
SIAM J. Optim., 13:805–842, 2003.
[22] D. Wachsmuth. Sufficient second-order optimality conditions for convex control con-
straints. Preprint 30-2004, Institut f¨ur Mathematik, TU Berlin, submitted, 2004.
[23] D. Wachsmuth. Regularity and stability of optimal controls of instationary Navier-
Stokes equations. Control and Cybernetics, 34:387–410, 2005.
Daniel Wachsmuth,
Institut f¨ur Mathematik,
Technische Universit¨at Berlin,
Str. des 17. Juni 136,
D-1062 Berlin, Germany.
