Solving Time-Dependent Optimal Control Problems in Comsol
Multiphysics by Space-Time Discretizations
Ira Neitzel1, Uwe Prüfert2, and Thomas Slawig3
1Technische Universität Berlin, DFG priority program SPP 1253
2Technische Universität Berlin, DFG research center Matheon
3Christian-Albrechts-Universität zu Kiel, DFG Cluster of Excellence The Future Ocean, and DFG
priority progam SPP 1253
Abstract
We use COMSOL Multiphysics to solve time-dependent optimal control problems for par-
tial differential equations whose optimality conditions can be formulated as a PDE. For a
class of linear-quadratic model problems we summarize known analytic results on existence
of solutions and first order optimality conditions that exhibit the typical feature of time-
dependent control problems, namely the fact that a part of the optimality system has to be
integrated backward in time. We present a strategy that is based on the treatment of the
coupled optimality system in the space-time cylinder. A brief motivation of this approach is
given by showing that the optimality system is elliptic in some sence. Numerical examples
show advantages and limits of the usage of COMSOL Multiphysics and of our approach.
Keywords
Optimal control of PDEs, finite element method, COMSOL Multiphysics.
1 Introduction
Optimal control problems subject to time-dependent partial differential equations are challenging from
the viewpoint of mathematical theory and even more so from numerical realization. Essentially, there are
two different approaches to solve such problems. The first one is the so-called "Discretize then Optimize"
strategy, where the optimal control problem is transformed into a nonlinear programming problem by
discretization. The second one is the function space based "Optimize then Discretize" strategy, that is
based on developing optimality conditions in function spaces that are discretized and solved.
For certain classes of problems it is possible to derive optimality conditions in PDE form, and the latter
strategy then involves solving systems of PDEs. It hence suggests itself to apply specialized PDE software
to solve these systems. In this paper, we aim at applying COMSOL Multiphysics for optimization, taking
advantage of the built-in routines to define, discretize and solve stationary and time-dependent PDEs via
the finite element method.
We consider simple model problems either with distributed control consisting of the objective functional
JQ(y, u) = 1
2ZQ
(y−yd)2+κ(u−ud)2dxdt, (1)
and the parabolic PDE with distributed control
yt−∆y=uin Q
∂ny+αy =gon Σ
y(t= 0) = y0in Ω
,(2)
or boundary control problems with objective functional
JΣ(y, u) = κ
2ZQ
(y−yd)2dxdt +κ
2ZΣ
(u−ud)2dsdt (3)
and the parabolic PDE with boundary control
yt−∆y=fin Q
∂ny+αy =uon Σ
y(t= 0) = y0in Ω
.(4)
Assumption 1. In this setting, Ω⊂RN,N= 1,2, is a spatial domain with sufficiently smooth boundary
∂Ω,(0, T )is a non-empty time interval, Σ := ∂Ω×(0, T ), and Q:= Ω ×(0, T ). Moreover, we consider
functions g∈L2(Σ) and y0∈L2(Ω) and controls u∈L2(Q)or u∈L2(Σ), depending on the underlying
PDE.
Short formulations of the model problems with control uand state ythen read
min JQ(y, u)subject to (2) (PQ)
and
min JΣ(y, u)subject to (4),(PΣ)
respectively.
We now briefly summarize well-established results on existence and uniqueness of optimal solutions to
Problems (PQ) and (PΣ), as well as first order necessary optimality conditions that are the basis for the
Optimize then Discretize approach. We start with the solvability of the state equation. A known result
by Wloka, [10], or Lions, [4] reads:
Theorem 1.1. For any triple (f, g, y0)∈L2(Q)×L2(Σ) ×L2(Ω) the initial-boundary value problem
yt−∆y=fin Q,
∂ny+αy =gon Σ,
y(t= 0) = y0in Ω
admits a unique solution
y∈W(0, T ) := {y∈L2(0, T ;H1(Ω))|
yt∈L2(0, T, H1(Ω)∗)}.
Here, H1(Ω)∗denotes the dual space to H1(Ω). Due to the fact that U=L2(Q)or U=L2(Σ),
respectively, are not empty, the following theorem holds by standard arguments.
Theorem 1.2. Under Assumption 1 and for JQand JΣdefined in (1), (3) and arbitrary κ > 0, the
optimal control problems defined in (PQ) and (PΣ) admit unique optimal controls u∗∈U=L2(Q)or
u∗∈U=L2(Σ), respectively.
Note that the associated optimal state y∗is uniquely determined by the optimal control u∗in either case
by Theorem 1. In the following theorem we formulate the first order necessary optimality conditions for
the control problems .
Theorem 1.3. Let u∗∈U=L2(Q)be the optimal control of Problem (PQ) and let y∗denote the
associated optimal state. Then there exists an adjoint state p∈W(0, T )as weak solution of
−pt−∆p=y∗−ydin Q
∂np+αp = 0 on Σ
p(t=T)=0 in Ω
,(5)
and the gradient equation
κ(u∗−ud) + p= 0 (6)
is fulfilled for almost all (x, t)∈Q.
Analogously, let u∗∈U=L2(Σ) be the optimal control of (PΣ) and let y∗denote the associated optimal
state. Then there exists an adjoint state p∈W(0, T )satisfying (5), and (6) is fulfilled for almost all
(x, t)∈Σ.
For more details, we refer to [4] or [9].
We point out here that optimality conditions can not easily be formulated for all types of PDE-control
problems. There is well-established theory available for linear-quadratic problems of the above types, as
2
well as for some nonlinear problems of similar structure. Also, if additional bounds on the controls are
given, well-known results can be applied. We refer for example to [4] for the treatment of these problems.
Another challenging problem are additional pointwise constraints on the state y. These types of problems
have been subject to intensive research over the last years, and the theory is far from being complete.
An attempt to handle certain state-constrained problems in COMSOL Multiphysics based on available
theory has been undertaken in [7]. We refer to this paper and the references therein for further reading.
Theorem 3 reveals a typical feature of time-dependent optimal control problems: While the state equation
(2) is an initial boundary value problem and hence a "forward-in-time" equation, the adjoint equation
(5) runs backward in time. Even though (5) can by the time transformation τ=T−tbe transformed
into an initial-boundary-value problem
pτ−∆p= ˜y∗−˜ydin Q
∂np+αp = 0 on Σ
p(τ= 0) = 0 in Ω
(7)
where ˜y∗(x, τ) = y∗(x, T −τ), the reverse time directions for state and adjoint state remain a difficulty
that needs to be taken into account when solving such problems numerically.
A somewhat classical approach to deal with this problem is to sequentially solve the state and adjoint
equation and to update the control in a gradient based optimization algorithm. We considered this
strategy in [6], where our key objective has been to realize the reverse time directions without any low
level data storage or copying.
A different solution approach and the focus of this paper is to treat the coupled optimality system in
the whole space-time cylinder by interpreting the time variable as an additional space variable, cf. also
[3]. We will explain this strategy and its implementation in COMSOL Multiphysics, conduct numerical
experiments, and comment on the applicability as well as the limits of this approach.
2 Transformation of the optimality system into a H2,1-elliptic
PDE
To justify our approach we show that the optimality system for our class of parabolic OCPs is equivalent
to a bi-harmonic H2,1-elliptic PDE. We restrict the theory here to the case of distributed control problems.
For a detailed presentation we refer to [8]. Without loss of generality, in this section we set ud≡0.
Definition 2.1. We define
¯
H2,1(Q):=H2,1(Q)∩ {u∈H2,1(Q) : ∂nu= 0 on Γ}
∩ {u∈H2,1(Q) : u(T) = u(0) = 0}.(8)
with its natural norm
kukH2,1(Q)=¡kutk2+kuk2+k∇uk2+k∆uk2¢1
/2.
This space is quite similar to the one used in [1] or [2].
For minimizing the notational effort we drop in this section the superscript ∗(which stands for optimality)
and write e.g. yinstead of y∗.
Theorem 2.2. ([8],Theorem 3.2) The optimality system given by Theorem (1.3) is equivalent to the
biharmonic PDE
−d2
dt2p+ ∆2p−2c0∆p+µc2
0+1
κ¶p=−¡d
dt yd−∆yd+c0yd¢in Q
∂n(−d
dt p−∆p+c0p) = g−∂nyd
∂np= 0 ¾on Σ(9)
−d
dt p(x, 0) −∆p(x, 0) + c0p(x, 0) = y0−yd(x, 0) on Σ0
p(x, T )=0 on ΣT.
To fit into the space ¯
H2,1(Q), we homogenize the biharmonic equation: We assume that ˆp∈H2,1(Q)is
a function that fulfills the boundary condition of the biharmonic equation (9), but not necessarily the
3
biharmonic equation itself. By v=p+ ˆp, the function vis the solution to the equation
−d2
dt2v+ ∆2v+ 2c0∆v+µc2
0+1
κ¶v=−¡d
dt yd−∆yd+c0yd¢+f(ˆp)in Q
∂n(−d
dt v−∆v+c0v)=0
∂nv= 0 ¾on Σ(10)
−d
dt v(x, 0) −∆v(x, 0) + c0v= 0 on Ω× {0}
v(x, T ) = 0 on Ω× {T}.
Lemma 2.3. ([8], Lemma 3.3) The homogenized equation (10) is equivalent to
a[v, w] = F(w)∀w∈¯
H2,1(Q)
where
a[v, w] = ZZ
Q
d
dt vd
dt w+ ∆v∆w+ 2c0∇v∇w+µc2
0+1
κ¶vw dxdt (11)
+Z
Ω
c0v(0)w(0) + ∇v(0)∇w(0) dx
is a symmetric bilinear form and F∈¡H2,1(Q)¢∗.
Lemma 2.4. The bilinear form a[v, w]is bounded and V-elliptic with respect to the space ¯
H2,1, i.e.
there are constants c1, c2>0such that
(i)
a[v, w]≤c1kvkH2,1(Q)kwkH2,1(Q)
(ii)
a[v, v]≥c2kvk2
H2,1(Q)
for all v, w ∈¯
H2,1(Q).
Proof. This lemma is a combination of Lemma 3.4 and Lemma 3.5 in [8].
Theorem 2.5. The bilinear equation
a[v, w] = F(w)
has for all F∈¡H2,1(Q)¢∗a unique solution in H2,1(Q). There is a constant c > 0such that
kvkH2,1(Q)≤ckFk(H2,1(Q))∗.
Proof. The assertion of the theorem is a direct consequence of the monotone operator theorem.
3 Treating the Reverse Time directions by Simultaneous Space-
Time Discretization
From the gradient equation (6), holding in the whole space-time domain Qor the boundary Σdepending
on the type of problem, we obtain u∗=ud−1
κp, where udand pare evaluated in the whole domain or
on the boundary, respectively. We insert this expression into the state equation (2) or (4). If the time
variable tis treated as an additional space variable we obtain boundary-value problems of the form
yt−∆y=ud−1
κp
−pt−∆p=y−yd¾in Q
∂ny+αy =g
∂np+αp = 0 ¾on Σ
y=y0in Ω× {0}
p= 0 in Ω× {T}.
4
for distributed control problems (PQ), as well as
yt−∆y=f
−pt−∆p=y−yd¾in Q
∂ny+αy =ud−1
κp
∂np+αp = 0 ¾on Σ
y=y0in Ω× {0}
p= 0 in Ω× {T}.
for boundary control problems (PΣ), i.e. we consider Qto be a purely spatial domain of dimension N+1
with boundary Σ∪Ω× {0} ∪ Ω× {T}.
The algebraic systems coming up from this discretization can be solved by the implemented solvers of
COMSOL Multiphysics without any further implementation effort, as we will see in the next section.
Naturally, a main issue of this solution approach is the dimension of the discretized system. The fact
that interpreting the time variable as additional space variable leads to an (N+ 1)-dimensional problem
and limits the applicability of this strategy to problems in at most two space dimensions. On the other
hand, this implementation strategy offers the opportunity of applying space-time adaptivity together.
Yet, one has to be aware of the fact that by the above solution strategy a parabolic system is treated
by elliptic solvers, which may induce instability issues. In a finite element discretization these may be
overcome by discontinuous ansatz functions, as they are becoming an alternative in the optimal control
community, cf. for instance [5].
4 Implementation and Numerical Examples
4.1 Distributed control
We consider first an example problem of form (PQ) with distributed control where the data is given by
Ω = (0, π),T=π,κ= 0.01,α= 0, as well as g=−sin(t),and the desired functions
yd(x, t) = sin(x) sin(t)−cos(x)(1 + π−t),
ud(x, t) = sin(x)(sin(t) + cos(t)) + 1
κcos(x)(π−t).
One can easily check that
y∗= sin(t) sin(x).
u∗= sin(t) sin(x) + cos(t) sin(x),
p∗=−κ(u∗−ud)
solves the optimality system given in Theorem 3.
We assume that the reader is familiar with all steps involved in building the fem-struc-
ture for solving a single PDE and present some details of the COMSOL Multiphysics Script code that
implements the described approach in the following. For the full code, we refer to [6].
First, we note that the time space domain Ω×(0, T ) = (0, π)×(0, π)is defined as a two dimensional
spatial domain with two "spatial" variables xand time:
fem.geom = rect2(0,pi,0,pi);
fem.sdim = {’x’ ’time’}.
Assuming that all given data is defined in the usual way we introduce a global expression for the control
u by
fem.globalexpr = {’u’ ’ud(x,time)-p/kappa’};
Moreover, we obtain for the definition of the PDE and the boundary conditions:
5
fem.equ.ga = { { {’-yx’ ’0’};
{’-px’ ’0’} } };
fem.equ.f = { {’-ytime+u’ ’ptime+y-yd(x,time)’} };
% boundaries: 1:t=0,2:x=pi,
% 3:t=pi,4:x=0
fem.bnd.ind = [1 2 3 2];
% boundary conditions:
fem.bnd.r = { {’y-y0’ 0};
{0 0};
{0 ’p’} };
fem.bnd.g = { {0 0};
{’g(time)’ ’0’};
{0 0} };
Once the definition of the problem is complete, the system can be solved by one of COMSOL Multiphysics’
implemented solution routines. For the example above, we use the linear, nonadaptive elliptic solver
femlin. We solve the problem on different grids, specified by the parameter hmax, using quadratic
finite element functions. In the following table we show the L2-errors ku∗−uhkQand ky∗−yhkQ
between the known optimal solution and the solution to the discretized problem (uh, yh)computed by
COMSOL Multiphysics on the different grids. Figures 1 show the computed optimal control uh, the
computed optimal state yh, as well as the associated adjoint state ph.
hmax ku∗−uhkQky∗−yhkQ
2−21.6417 ·10−23.2837 ·10−4
2−32.2293 ·10−33.0790 ·10−5
2−43.0615 ·10−45.0814 ·10−6
2−54.0305 ·10−54.9791 ·10−7
2−65.3730 ·10−65.9155 ·10−8
Table 1: Errors ku∗−uhkQand ky∗−yhkQ.
Figure 1: Optimal control uh, Optimal state yh, and adjoint state ph.
4.2 Boundary control
Now, we consider an example problem of form (PΣ) with boundary control, very similar to the distributed
control problem in the last section. The data is given by Ω = (0, π),T=π,κ= 0.01,α= 0, as well as
f= sin(x) cos(t) + sin(x) sin(t),
yd= sin(x) sin(t)−cos(x)(1 + π−t),
ud=−sin(t) + 1
κ(π−t).
One can easily check that
y∗= sin(t) sin(x),
u∗=−sin(t),
p∗= cos(x)(π−t)
6
solves the optimality system given in Theorem 3.
The implementation in COMSOL Multi- physics is very similar to the first example. We only need to
substitute ffor uas a source term of the state equation, as well as ufor g, since uenters the state
equation in the boundary conditions. This yields:
fem.equ.ga = { { {’-yx’ ’0’};
{’-px’ ’0’} } };
fem.equ.f = { {’-ytime+...
f(x,time)’...
’ptime+y-yd(x,time)’} };
% boundaries: 1:t=0,2:x=pi,
% 3:t=pi,4:x=0
fem.bnd.ind = [1 2 3 2];
% boundary conditions:
fem.bnd.r = { {’y-y0’ 0};
{0 0};
{0 ’p’} };
fem.bnd.g = { {0 0};
{’u’ ’0’};
{0 0} };
Similar to the distributed example, we use the linear, nonadaptive elliptic solver femlin on different
grids, specified by the parameter hmax. This time, we use linear finite element functions. The L2-errors
ku∗−uhkΣand ky∗−yhkQon the different grids are shown in Table 2. Figure 2 shows the computed
optimal control uh. The example has been constructed such that the optimal state y∗and the associated
adjoint state p∗are the same as for the distributed control problem.
hmax ku∗−uhkQky∗−yhkQ
2−22.1379 ·10−12.5347 ·10−2
2−31.0525 ·10−19.6450 ·10−3
2−44.3199 ·10−22.9205 ·10−3
2−51.7076 ·10−27.7394 ·10−4
2−65.7762 ·10−31.9415 ·10−4
Table 2: Errors for the boundary control problem.
Figure 2: Boundary control uh
4.3 An example in 2D
In this example we use the 3D capability of COMSOL Multiphysics to solve a problem in two space
dimensions. We consider the optimal control problem (PQ) where the space-time domain is defined by
Q= (0, π)2×(0, π)⊂R3and the functions yd,ud, and gare given by
yd= sin(x1) sin(x2) sin(t)−cos(x1) cos(x2)−2 cos(x1) cos(x2)(π−t)
ud= sin(x1) sin(x2) cos(t) + 2 sin(x1) sin(x2) sin(t) + 1
κcos(x1) cos(x2)(π−t)
g=−~n sin(t)(sin(x1),sin(x2))T
7
respectively. The optimal solutions are
y∗(x1, x2, t) = sin(x1) sin(x2) sin(t)
u∗(x1, x2, t) = sin(x1) sin(x2)(cos(t) + 2 sin(t))
p∗(x1, x2, t) = cos(x1) cos(x2)(π−t),
which can easily be checked by inserting them into the optimality conditions of Theorem 3.
The differences in the implementation compared to the first example are only due to the higher space
dimension of the problem. We define the 3D geometry by
% geometry and mesh:
fem.geom = block3(pi,pi,pi,...
’base’,’corner’,’pos’,[0 0 0]);
and also account for the additional dimension in the definition of the PDE and the boundary conditions:
fem.equ.ga = { { {’-yx1’ ’-yx2’ ’0’}
{’-px1’ ’-px2’ ’0’}
} };
fem.equ.f = { {’-ytime+u’ ’ptime...
+y-yd(x1,x2,time)’} };
fem.bnd.r = { {’y-y0’ 0};
{0 ’p’};
{0 0};
{0 0} };
fem.bnd.g = { {0 0};
{0 0};
{’g1(x1,time)...
-alpha*y’ ’-alpha*p’}
{’g2(x2,time)...
-alpha*y’ ’-alpha*p’}};
This time, we test the adaptive solver adaption with an initial grid determined by hauto and ngen, the
number of new grid generations, set to two. In 3D, the space-time grid consists of tetrahedrons and the
number of unknowns grows cubically when refining the grid. For that reason, we restrict our survey to
three initial grids generated using meshinit where refinement is controlled by using the parameter hauto
ranging from 7 to 5. We use again quadratic finite elements. In Figures 5–7 we present time-slice plots
of uh,yh, and ph.
Figure 3: Computed solution of Example 3.
hauto ku∗−uhkQky∗−yhkQ
73.1710 ·10−14.7920 ·10−3
61.7107 ·10−12.3017 ·10−3
55.0385 ·10−25.4455 ·10−4
Table 3: Errors to the 2D example, adaptive solver
8
5 Conclusion
We have successfully applied the finite element package COMSOL Multiphysics to simple time-dependent
optimal control problems subject to PDE constraints by utilizing an Optimize then Discretize strategy.
Figure 4: Error flow
The introduced strategy has proven to work reasonably well
for our simple example problems. We take advantage of the
fact that optimality conditions can be written in PDE form,
which allows to apply a specialized PDE solver for optimiza-
tion. The method we focused on is easily implementable and
may well serve as a first step towards optimizing a given goal
without the use of specialized optimization routines. The
proposed approach does not substitute the use of special-
ized optimization software. One reason for that is that el-
liptic solvers are used for time-dependent parabolic control
problems. This means in particular that the special role of
the time is ignored, especially the fact that the time deriva-
tive of state and adjoint state is in general in a weak space
L2(0, T, H1(Ω)∗), cf. Theorem 1. However, for cases with
higher regularity our approach is fully justified, cf. Section
2. In Figure 4, we show the error propagation through the
space-time domain for the one-dimensional distributed con-
trol problem to illustrate this behavior due to solving a singular elliptic system.
The applicability of this strategy to other problems, even if optimality conditions can be formulated in
PDE form, has to be decided on a case-to-case basis. Additional limitations are given by the size of the
problem.
References
[1] A. Borzi. Multigrid methods for parabolic distributed optimal control problems. J. Comp. Appl.
Math., 157:365-382, 2003.
[2] Guido Büttner. Ein Mehrgitterverfahren zur optimalen Steuerung parabolischer Probleme. PhD the-
sis, Technische Universität Berlin, 2004.
[3] M .Hinze ,M. Köster, and S. Turek. A Space-Time Multigrid Solver for Distributed Control of the
Time-Dependent Navier-Stokes System. SPP 1253 Report, 2008.
[4] J. L. Lions. Optimal Control of Systems Governed by Partial Differential Equations. Springer-Verlag,
Berlin, 1971.
[5] D. Meidner and B. Vexler. A priori error estimates for space-time finite element discretization of
parabolic optimal control problems Part I: Problems without control constraints. SIAM Control
Optim., 47(3):1150–1177, 2008.
[6] I. Neitzel, U. Prüfert, and T. Slawig. Strategies for time-dependent PDE control using an integrated
modeling and simulation environment. Part one: Problems without inequality constraints. Technical
Report 408, Matheon, Berlin, 2007.
[7] I. Neitzel, U. Prüfert, and T. Slawig. Strategies for time-dependent PDE control with inequality
constraints using an integrated modeling and simulation environment. Numerical Algorithms, 2008,
DOI 10.1007/s11075-008-9225-4.
[8] I. Neitzel, U. Prüfert, and T. Slawig. On Solving Parabolic Optimal Control Problems by Using
Space-Time Discretization. Technical Report 05-2009, TU Berlin, Berlin, 2009.
[9] F. Tröltzsch. Optimale Steuerung partieller Differentialgleichungen. Theorie, Verfahren und Anwen-
dungen. Vieweg, Wiesbaden, 2005.
[10] J. Wloka. Partielle Differentialgleichungen. Teubner-Verlag, Leipzig, 1982.
The complete COMSOL Multiphysics-code of all examples is available on
http://www.math.tu-berlin.de/Strategies-for-time-dependent-PDE-control
9
