See discussions, stats, and author profiles for this publication at: [Link]

net/publication/387265702

Numerical aspect of uniform ensemble controllability for linear systems

Preprint · December 2024

DOI: 10.13140/RG.2.2.30147.52004

CITATIONS READS
0 107

2 authors:

Martin Lazar Jérôme Lohéac

University of Dubrovnik University of Lorraine
51 PUBLICATIONS 397 CITATIONS 48 PUBLICATIONS 434 CITATIONS

SEE PROFILE SEE PROFILE

All content following this page was uploaded by Martin Lazar on 20 December 2024.

The user has requested enhancement of the downloaded file.

 Numerical aspect of uniform ensemble controllability for
linear systems
Martin Lazar1 Jérôme Lohéac2
December 20, 2024

Abstract
We consider the ensemble controllability problem for a linear time-invariant system
ẋ(t, θ) = A(θ)x(t, θ) + B(θ)u(t), where A and B are continuous matrices with respect to
the parameter θ, which belongs to some compact set Θ ⊂ R. Given any continuous initial
state datum θ 7→ x0 (θ) and any continuous target state θ 7→ x1 (θ), we investigate the numer-
ical computation of a θ-independent open loop control u such that x0 is steered, in a given
time T > 0, at a distance ε > 0 of x1 in the uniform norm (with respect to the parameter).
We approach the problem both theoretically and numerically. Using the Fenchel-
Rockafellar duality, we first prove the existence and uniqueness of the ensemble control of
a minimal L2 norm. The numerical recovery of the optimal control is obtained by solving
the dual problem, which consists in the unconstrained minimization of a non-differentiable
functional in the space of Radon measures.

Keywords: parameter dependent systems, ensemble controllability, linear control system, Fenchel-
Rockafellar duality, numerical algorithms.

Contents
1 Introduction 2
2 Preliminaries 3
3 Minimal L2 -norm controls 5
3.1 Optimality conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
3.2 Continuity properties of the minimizers . . . . . . . . . . . . . . . . . . . . . . . . 10
4 Systems with a finite number of parameter values 14
5 Discretized problem and convergence 18
6 Numerical strategies 22
6.1 Greedy coordinate descent . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
6.2 Reduced minimization problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
6.3 Proposed numerical strategies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
7 Numerical examples 27
A Proof of Proposition 6.1 32
B Some unreachable states 35

1 Faculty of Electrical Engineering and Applied Computing, University of Dubrovnik, Croatia

2 Université de Lorraine, CNRS, CRAN UMR 7039 Nancy, France

1
 1 1 Introduction
2 The concept of ensemble controllability [1, 8, 16, 21, 23] needs no more publicity nowadays with
3 regard to its practical importance. Its aim is to control systems whose natural dynamic depends
4 on some unknown or uncertain parameters in a robust and computationally efficient manner. In
5 particular, the goal is to find a parameter invariant control that performs well for each particular
6 realization of the parameters. This kind of problem arises for instance in the transport of quantum
7 particles, where the goal is to develop external excitation that can simultaneously steer the ensemble
8 of systems with variation in their internal parameters from an initial state to a desired final
9 state [20]. One can also see [3] where it is used to study the controllability of the Bloch equation,
10 for an ensemble of non-interacting half-spins, in a static magnetic field, with dispersion in the
11 Larmor frequency. Other real life practical example can be found, as for instance in non-holonomic
12 systems theory where ensemble controllability is used to derive an approximate steering algorithm
13 for a non-holonomic unicycle in the presence of model perturbation [4]. Nowadays, however,
14 the most widely application is related to machine learning, in particular to training processes in
15 supervised learning. By using the dynamic perspective of neural networks, the problem of their
16 training can be recast as a control problem in which the network coefficients (weights and biases)
17 stand for data(parameter)-independent controls [13].
18 From the generic notion of ensemble controllability, many other notions emerged as the one of
19 uniform ensemble controllability [16, 20, 21], uniform ensemble reachability, Lq -ensemble reachabil-
20 ity [11]. Even if this topic is studied actively, it is far from being well understood. The necessary and
21 sufficient condition of ensemble controllability can be obtained as a special form of the generaliza-
22 tion of the Kalman rank theorem to infinite dimensional Banach spaces (e.g. [28, Theorem 3.1.1]).
23 However, such a characterization is primarily theoretical and not suitable for checking the control-
24 lability properties of a particular system. In the above-mentioned non-exhaustive references, the
25 authors have derived necessary and/or sufficient conditions for linear time invariant systems to be
26 uniformly ensemble controllable in some special cases. From these references, especially [8, 10, 11],
27 it could be seen that the topology of Θ plays a crucial role in the controllability issue. Roughly
28 speaking, for regular enough parameter dependent systems, uniform ensemble controllability can-
29 not hold if the compact set Θ is homeomorphic to some compact subset of Rd (with d ⩾ 2) with
30 non-empty interior.
31 Most papers on the topic analyze conditions under which a system is ensemble controllable, but
32 do not consider the associated optimization problem of finding a feasible control of minimal norm,
33 nor do they deal with numerical recovery of the solution. An approximate construction of a feasible
34 control is provided only in some special cases (e.g. [3, 4]). In [15] the authors consider an associated
35 penalization problem in which the constraint of reaching the prescribed target with some a priori
36 given precision is replaced by an additional term in the cost functional that penalizes deviation of
37 the final state from the target. Consequently, such a problem consists in the unconstrained mini-
38 mization of a smooth functional and is numerically easier to handle. In [25], the uniform ensemble
39 control problem is replaced by a simultaneous one (consisting of hitting a reachable target). The
40 author discusses the strategies by which the latter can be approximated by its discretized version
41 (which is equivalent to an exact control problem for a large, finite-dimensional system). Let us also
42 mention [26], where an ensemble control is designed using polynomial approximations. However,
43 the results of [26] seem only applicable for single input systems.
44 In [18] we have characterized the optimal ensemble control in the discrete case (i.e., for a finite
45 number of parameters). Then, assuming Lipschitz regularity, we showed how to obtain a sub-
46 optimal control for the continuous case. In the present article, we directly tackle the case of a
47 continuous set Θ, and derive optimality conditions. Based on them we propose several numerical
48 strategies for recovering the minimal norm control. More precisely, we prove the existence and
49 uniqueness of a minimal norm ensemble control, and express it through the solution of the corre-

2
50 sponding dual problem, which consists in the unconstrained minimization of a non-differentiable
51 functional in the space of Radon measures. The numerical strategies that we propose in this paper
52 are obtained using the optimality condition on the adjoint problem.
53 The paper is organized as follows. In the next section, we recall the notion and some basic
54 results on uniform ensemble controllability. In Section 3, we introduce the optimization problem of
55 finding a uniform ensemble control of a minimal L2 norm and derive the associated dual problem.
56 We demonstrate the existence of a unique solution for both problems and list their properties. In
57 particular, the continuity of solutions with respect to the problem entries is obtained. Section 4
58 deals with a special but important case of discretized parameters (obtaining only a finite number of
59 values), where we also discuss the sparsity properties of the optimal measure. The convergence of
60 solutions when passing from discrete to continuous parameter setting is demonstrated in Section 5.
61 Numerical iterative algorithms are presented in Section 6, while their performances are checked
62 and discussed on several particular examples in Section 7. Appendices A and B contain some
63 technical and related results.
64 Notation. The following notations are used in the paper.
65 N = {0, 1, 2, . . .} stands for the set of natural numbers including zero, while N∗ denotes the set of
66 positive integers. R+ is the set of nonnegative real numbers, and we set R∗+ = R+ \ {0}. For a
67 matrix M ∈ Rn×m , M ∗ ∈ Rm×n is the transpose matrix of M . |·| stands for the Euclidean norm
68 in Rn , and ⟨·, ·⟩ is the scalar product on Rn . The open (respectively closed) ball of Rn centered on
69 some x and of radius ε > 0 is denoted by B(x, ε) (respectively B(x, ε)).
70 For X a Banach space and f : X → R ∪ {∞}, we define dom f = {x ∈ X | f (x) < ∞}, and cont f
71 is the set of points x ∈ dom f where f is continuous. For S ⊂ X, int S is the interior of S.
n
72 Given n ∈ N∗ , by C(Θ) we denote the space of continuous complex vector functions on Θ ⊆ R
n
73 equipped with L norm, i.e., C(Θ) = C(Θ; Rn ) and ∥f ∥C(Θ)n = supθ |f (θ)|. B∞ (y, ε) is the closed
∞
n
74 ball in the space C(Θ) centered in y and of radius ε.
75 M(Θ) stands for the set of Radon measures supported on Θ, while by M(Θ)n we
n
76 denote the dual space to C(Θ) with the dual product given by ⟨µ, f ⟩M(Θ)n ,C(Θ)n =
Pn ∗ n×m m
i=1 Θ fi (θ) dµi (θ) = Θ f (θ) dµ(θ). In addition, for F ∈ C(Θ)
, we define F ∗ ♭ µ ∈ M(Θ)
R R
77
∗ ∗
by ⟨F ∗ ♭ µ, f ⟩M(Θ)m ,C(Θ)m = ⟨µ, F f ⟩M(Θ)n ,C(Θ)n = Θ f (θ) F (θ) dµ(θ).
R
78

79 2 Preliminaries
80 Let Θ be a compact subset of R and consider for every constant θ ∈ Θ the system described by

ẋ(t, θ) = A(θ)x(t, θ) + B(θ)u(t) (t > 0, θ ∈ Θ). (2.1)

n×n n×m
81 In (2.1), we assume that A ∈ C(Θ) and B ∈ C(Θ) . The input u is independent of the
82 parameter θ, and ẋ is the derivative of x with respect to the time variable t. Given an initial state
m
83 x0 ∈ C(Θ)n and an input u ∈ L2 (0, T ) , the solution of (2.1) at time t, is denoted by x(t, θ; x0 , u).
84 More precisely, we have,
Z t
x(t, θ; x0 , u) = etA(θ) x0 (θ) + e(t−τ )A(θ) B(θ)u(t) dt (t ⩾ 0, θ ∈ Θ). (2.2)
0

85 Definition 2.1 (Uniform Ensemble Controllability). The system (2.1) is said to be uniformly
86 ensemble controllable if for any x0 , x1 ∈ C(Θ)n and any ε > 0, there exist a time T > 0 and a
m
87 control u ∈ L2 (0, T ) , such that

sup x(T, θ; x0 , u) − x1 (θ) ⩽ ε.

θ∈Θ

3
88 The terminology uniform in the above definition is used to emphasize the fact that the set of
89 output functions is the continuous one endowed with the uniform norm. This definition actually
90 makes sense since for any continuous initial state datum x0 and any θ-independent input u ∈
m
91 L2 (0, T ) , the output trajectory is continuous with respect to the parameter θ.

92 Remark 2.2. The systems considered in the present paper focus on the case where A(θ) and B(θ)
93 are real matrices. However, all the obtained results can be applied to complex matrix, having in
94 mind that if A(θ) ∈ Cn×n and B(θ) ∈ Cn×m  are complex  matrices, thenthe complex system can
ℜA −ℑA ℜB −ℑB
95 be identified to a real system with matrices ∈ R2n×2n and ∈ R2n×2m ,
ℑA ℜA ℑB ℜB
96 where ℜ and ℑ denote the real and imaginary part respectively.
97 Let us also mention that it is classical that the set of parameter Θ is assumed to be a subset of C.
98 In this paper we only consider the case where Θ ⊂ R. This assumption is made, firstly, because
99 it avoids technical difficulties and, secondly, because, as said in the introduction, from [8, 10, 11]
100 typical sets Θ ∈ C for which uniform ensemble controllability can be expected are one dimensional
101 objects.

102 The following facts can be found in [11].

103 Proposition 2.3. If (2.1) is uniformly ensemble controllable, then it is uniformly ensemble con-
n
104 trollable in any time T > 0, i.e., for any x0 , x1 ∈ C(Θ) , any T > 0 and any ε > 0, there exists a
m
105 control u ∈ L (0, T ) , such that supθ∈Θ x(T, θ; x , u) − x1 (θ) ⩽ ε.
2 0

106 Proposition 2.4. If (2.1) is uniformly ensemble controllable, then for every θ ∈ Θ, the pair
107 (A(θ), B(θ)) is controllable.
108 Remark 2.5. The direct consequence of the last result is that the control operator can not vanish
109 in any point, i.e. B(θ) ̸= 0 for every parameter θ.
110 More generally, we have the following exact controllability property.

111 Proposition 2.6. Assume (2.1) is uniformly ensemble controllable. For every K ∈ N∗ and every
112 θ1 , . . . , θK ∈ Θ, two by two distinct, the pair (A, B) is controllable, where
   
A(θ1 ) B(θ1 )
A= .. ∈R
Kn×Kn
, and B =  ...  ∈ RKn×m .
   
.
A(θK ) B(θK )

113 Let us also mention the unique continuation property for the adjoint system (cf. [10]). To this
114 end, we denote by M(Θ) the set of Radon measures supported in Θ, and we define the control to
m n
115 final state operator ΦT ∈ L(L2 (0, T ) , C(Θ) ) by
Z T
m
(ΦT u)(θ) = e(T −t)A(θ) B(θ)u(t) dt (θ ∈ Θ, u ∈ L2 (0, T ) ), (2.3)
0

116 i.e., (ΦT u)(θ) = x(T, θ, 0, u) with x given by (2.2). When there is no fear of ambiguity, we shall
117 use abbreviation xT = ΦT u.
n m
118 Note that the adjoint of ΦT is Φ∗T ∈ L(M(Θ) , L2 (0, T ) ) defined by
Z D E
∗ ∗ ∗
(Φ∗T µ)(t) = B(θ) e(T −t)A(θ) dµ(θ) = e(T −t)A ♭ µ, B
Θ M(Θ)n ,C(Θ)n×m
n
(t ∈ (0, T ), µ ∈ M(Θ) ). (2.4)

4
119 Proposition 2.7. The system (2.1) is uniformly ensemble controllable in time T > 0 if and only if
 Z 
n ∗ (T −t)A(θ)∗
µ ∈ M(Θ) | ∀t ∈ [0, T ], B(θ) e dµ(θ) = 0 = {0}. (2.5)
Θ

120 By the last proposition the uniform ensemble controllability implies the unique continuation
121 property of the adjoint operator, i.e.,

Φ∗T µ = 0 =⇒ µ = 0. (2.6)
n n
122 In particular, the Gramian operator defined as ΛT := ΦT Φ∗T ∈ L(M(Θ) ; C(Θ) ) is injective.

123 3 Minimal L2 -norm controls

124 3.1 Optimality conditions
125 The aim of this section is to characterize the minimal L2 -norm control steering the system (2.1)
126 from x0 to a distance ε > 0 of x1 in a given time T > 0. More precisely, given T > 0, ε > 0 and
n
127 x0 , x1 ∈ C(Θ) , we aim to find a minimizer of
1 2
min 2 ∥u∥L2 (0,T )m ,
T
(P )
Z
m
u ∈ L (0, T ) 2
and sup y(θ) − e(T −t)A(θ) B(θ)u(t) dt ⩽ ε.
θ∈Θ 0

n
128 where y ∈ C(Θ) is the final target adjusted by the solution of the homogeneous equation, i.e.,

y(θ) = x1 (θ) − eT A(θ) x0 (θ) (θ ∈ Θ).

129 When there is no fear of ambiguity, we shall also refer to y as (adjusted) target.
m n
130 We rewrite the problem (P ) by introducing functions f ∈ C(L2 (0, T ) ) and g : C(Θ) →
131 R+ ∪ {∞} defined as,
2 m
f (u) = 1
∥u∥L2 (0,T )m
2 (u ∈ L2 (0, T ) ),
(
0 if ∥ξ − y∥L∞ (Θ)n ⩽ ε, n
g(ξ) = (ξ ∈ C(Θ) ).
∞ otherwise
n
132 In particular, g is the indicator function of the ball B∞ (y, ε) ⊂ C(Θ) .
133 Hence, (P ) can be expressed as

min f (u) + g(ΦT u).

u∈L2 (0,T )m

134 Note that f and g are convex functions, and our aim is to approach the optimization problem (P )
135 by the Fenchel-Rockafellar duality. To this end, we consider the Fenchel conjugates of f and g
136 defined as
  m
2
f ∗ (v) = sup ⟨v, u⟩L2 (0,T )m − f (u) = 21 ∥v∥L2 (0,T )m (v ∈ L2 (0, T ) ),
u∈L2 (0,T )m
 
g ∗ (µ) = sup ⟨µ, ξ⟩M(Θ)n ,C(Θ)n − g(ξ) = sup ⟨µ, ξ⟩M(Θ)n ,C(Θ)n
ξ∈C(Θ)n ξ∈C(Θ)n
∥ξ−y∥L∞ (Θ)n ⩽ε
n
= ⟨µ, y⟩M(Θ)n ,C(Θ)n + ε ∥µ∥M(Θ)n (µ ∈ M(θ) ).

5
 n
137 Here, and in the rest of this paper, for µ = (µ1 , . . . , µm ) ∈ M(Θ) , we set
n
X
∥µ∥M(Θ)n = sup ⟨µ, ξ⟩M(Θ)n ,C(Θ)n = sup ⟨µi , ξi ⟩M(Θ),C(Θ) .
ξ∈C(Θ)n ξ1 ,...,ξn ∈C(Θ) i=1
∥ξ∥L∞ (Θ)n ⩽1 supΘ n 2
P
i=1 |ξi | ⩽1

138 Specially, if the components of µ are written

Din the form
 µi E = fi♭ ν for some i-independent scalar
and positive measure ν, then ∥µ∥M(Θ)n = = Θ |f (θ)| dν(θ). In the
pP R
2
139
i fi ♭ ν, 1
M(Θ),C(Θ)
special case of µ Pbeing a combination of disjoint Dirac masses µ = k mk δθk for mk ∈ Rn , its
P
140

141 norm is given by k |mk |.

142 Theorem 3.1. Assume (2.1) is uniformly ensemble controllable. Then for every y ∈ C(Θ)n and
m
143 every ε > 0, the minimization problem (P ) admits a unique minimizer uopt ∈ L2 (0, T ) given by
144 the expression
uopt = Φ∗T µopt ,
n
145 where µopt ∈ M(Θ) minimizes the dual functional
2 n
J(µ) = 1
2 ∥Φ∗T µ∥L2 (0,T )m − ⟨µ, y⟩M(Θ)n ,C(Θ)n + ε ∥µ∥M(Θ)n (µ ∈ M(Θ) ). (3.1)

146 Furthermore, the minimizer µopt is unique, and satisfies

2
∥uopt ∥L2 (0,T )m = ⟨µopt , y⟩M(Θ)n ,C(Θ)n − ε ∥µopt ∥M(Θ)n , (3.2)
147
supp µopt ⊂ {θ ∈ Θ | |(ΦT uopt )(θ) − y(θ)| = ε} (3.3)
n
148 and, for every test function υ ∈ C(Θ) such that ΦT uopt + υ ∈ B∞ (y, ε),

⟨µopt , υ⟩M(Θ)n ,C(Θ)n ⩾ 0. (3.4)

149 Finally, µopt is the zero measure if and only if y ⩽ ε.

150 Remark 3.2. As |(ΦT uopt )(θ) − y(θ)| is the distance of the optimal final state from the target x1 ,
151 the characterization (3.3) restricts the support of the optimal measure to the parameters θ ∈ Θ
152 for which the optimal final state lies exactly on the boundary of the target ball B∞ (x1 , ε).
153 Proof. We first observe that f is strictly convex, and dom g is not empty. Furthermore, we observe
m
154 that u ∈ L2 (0, T ) | ΦT u ∈ dom g is a convex set and g is lower semi-continuous. These facts
155 ensure the existence and uniqueness of a minimizer of (P ) (see e.g. [24, Theorem 2.19]). Note
m
156 also that cont g = int dom g and dom f = L2 (0, T ) , hence, the uniform ensemble controllability
157 assumption guaranties that ΦT dom f ∩ cont g ̸= ∅. Thus, according to Fenchel duality (see e.g. [5,
158 Theorem 4.4.3], [24, Theorem 3.51] or [12]), we have

min f (u) + g(ΦT u) = max −f ∗ (Φ∗T µ) − g ∗ (−µ) = − min f ∗ (Φ∗T µ) + g ∗ (−µ)

u∈L2 (0,T )m µ∈M(Θ)n µ∈M(Θ)n
2
=− min 1 ∥Φ∗T µ∥L2 (0,T )m − ⟨µ, y⟩M(Θ)n ,C(Θ)n + ε ∥µ∥M(Θ)n .
µ∈M(Θ)n 2

n
159 In addition (3.1) admits a minimizer µopt ∈ M(Θ) , and we have

f (uopt ) + g(ΦT uopt ) = −f ∗ (Φ∗T µopt ) − g ∗ (−µopt ) (3.5)

160 and
Φ∗T µopt ∈ ∂f (uopt ) and − µopt ∈ ∂g(ΦT uopt ),

6
161 where ∂ stands for the sub-differential. As f is a smooth function, its sub-differential is a singleton,
162 and from the first of the above inclusions we get that uopt = Φ∗T µopt .
163 Noticing that f ∗ = f , g(ΦT uopt ) = 0, we get from (3.5), 2f (uopt ) = −g ∗ (−µopt ), which leads
164 to (3.2).
165 We also have that xTopt := ΦT uopt ∈ C(Θ)n is such that xTopt − y L∞ (Θ)n ⩽ ε, i.e.,
166 g(xTopt ) = 0. Then, −µopt ∈ ∂g(xTopt ) implies that for every v ∈ C(Θ)n , g(xTopt + v) =
167 g(xTopt + v) − g(xTopt ) ⩾ ⟨−µopt , v⟩M(Θ)n ,C(Θ)n . In particular, if xTopt + v − y L∞ (Θ)n ⩽ ε, we
168 have 0 ⩽ ⟨µopt , v⟩M(Θ)n ,C(Θ)n , thus proving (3.4).
169 From here we obtain the characterization (3.3) of the support of µopt . Indeed, for ev-
ery v ∈ C(Θ)n such that supp v ⊂ θ ∈ Θ | xTopt (θ) − y(θ) < ε , there exists α > 0 such

170

171 that xTopt + αv − y L∞ (Θ)n < ε and xTopt − αv − y L∞ (Θ)n < ε. Hence, we shall have 0 ⩽
172 ⟨µopt , v⟩M(Θ)n ,C(Θ)n and 0 ⩽ − ⟨µopt , v⟩M(Θ)n ,C(Θ)n , leading to 0 = ⟨µopt , v⟩M(Θ)n ,C(Θ)n .
173 Dealing with the uniqueness of µopt , assume by contradiction that µ1 and µ2 are minimizer
174 of J. By the uniqueness of uopt , we have uopt = Φ∗T µ1 = Φ∗T µ2 . According to Proposition 2.7, we
175 conclude that µ1 = µ2 .
176 Finally, if ∥y∥C(Θ)n ⩽ ε (meaning that the free dynamics brings the state within the target
177 ball) by the definition (3.1) of the functional J and the Cauchy-Schwarz inequality, we have that
n
178 for every measure µ ∈ M(Θ)

J(µ) ⩾ − ∥µ∥M(Θ)n ∥y∥C(Θ)n + ε ∥µ∥M(Θ)n ⩾ 0 = J(0),

179 i.e., the functional J is nonnegative and µopt = 0 is its unique minimizer.
180 The reverse implication holds trivially, as if ∥y∥C(Θ)n > ε then a nonzero control is required in
181 order to satisfy the constraint (of reaching the target ball). Using classical optimization techniques
182 arguments, it can be even shown that in that case ∥µopt ∥M(Θ)n is bounded from below by a positive
183 constant depending linearly on the difference ∥y∥C(Θ)n − ε (e.g. [17, Proposition 2.3]).
184 Proposition 3.3. Using the assumptions and notations introduced in Theorem 3.1, the optimal
185 measure µopt is of the form

dµopt (θ) = α(θ) y(θ) − xTopt (θ) d |µopt | (θ), (3.6)

Pn
186 where |µopt | = i=1 |µopt,i | ∈ M(Θ) is the sum of total variations of components of the vector
187 measure µopt , while α : Θ → R+ is an integrable (with respect to |µopt |) nonnegative scalar function.
188 Proof. Step 1. We first show that

⟨µopt , υ⟩M(Θ)n ,C(Θ)n ⩾ 0, (3.7)

189 for every test function υ ∈ C(Θ)n such that υ(θ), y(θ) − xTopt (θ) ⩾ 0 whenever (xTopt )(θ) − y(θ) =
190 ε, θ ∈ Θ.
191 Suppose first that v satisfies the last inequality strictly, i.e., that v(θ), y(θ) − xTopt (θ) > 0
192 for every θ ∈ Θ satisfying |xopt (θ) − y(θ)| = ε. Then we show there exists a θ independent
193 scalar α > 0 small enough such that ∥xopt + αv − y∥L∞ (Θ)n ⩽ ε. Indeed, let us set Θopt =
194 {θ ∈ Θ | |xopt (θ) − y(θ)| = ε}, since θ 7→ |xopt (θ) − y(θ)| is continuous, we have that Θopt is
195 compact. On this set, we have ⟨v(θ), y(θ) − xopt (θ)⟩ > 0, using the uniform continuity on Θ of this
196 function, we get the existence of δ > 0 such that ⟨v(θ), y(θ) − xopt (θ)⟩ > 2δ for every θ ∈ Θopt ,
197 and there also exists η > 0 such that ⟨v(θ), y(θ) − xopt (θ)⟩ > δ for every θ ∈ Θopt + B(0, η). Using
198 the fact that Θ \ (Θopt + B(0, η)) is compact, we also obtain the existence of εe < ε such that
199 |xopt (θ) − y(θ)| ⩽ εe for every θ ∈ Θ \ (Θopt + B(0, η)). All in all, we have for every α > 0,

7
200 • for θ ∈ Θopt + B(0, η),

2 2 2
|xopt (θ) + αv(θ) − y(θ)| = |xopt (θ) − y(θ)| − 2α ⟨v(θ), y(θ) − xopt (θ)⟩ + α2 |v(θ)|
2
⩽ ε2 − 2αδ + α2 ∥v∥L∞ (Θ)n .

201 • for θ ∈ Θ \ (Θopt + B(0, η)),

2 2
xTopt (θ) + αv(θ) − y(θ) ⩽ εe2 + 2αe
ε ∥v∥L∞ (Θ)n + α2 ∥v∥L∞ (Θ)n .

202 Hence, taking α > 0 small enough, we have xTopt (θ) + αv(θ) − y(θ) ⩽ ε for every θ ∈ Θ. By (3.4)
203 this ensures that ⟨µopt , v⟩M(θ)n ,C(Θ)n ⩾ 0 for every v ∈ C(Θ)n such that v(θ), y(θ) − xTopt (θ) > 0
204 for every θ ∈ Θopt .
205 It remains to prove that (3.7) holds true, even if v(θ), y(θ) − xTopt (θ) ⩾ 0 for every θ ∈ Θopt . In
fact, if v satisfies the last condition, then v(θ) + α y(θ) − xTopt (θ) , y(θ) − xTopt (θ) > 0 for every


206

θ ∈ Θopt and every α > 0. This ensures that 0 ⩽ µopt , v(θ) + α y(θ) − xTopt (θ) M(Θ)n ,C(Θ)n , and


207

208 hence, taking the limit α → 0, we obtain the desired result.

209 Step 2. Now we prove the relation (3.6)
210 First, let us note that µopt is absolutely continuous with respect to the measure |µopt |. By the
211 Lebesgue-Radon-Nikodym theorem (see e.g. [14]) there exists an integrable (with respect to |µopt |)
212 vector function φ such that dµopt (θ) = φ(θ) d|µopt |(θ).
213 It remains to prove that φ = α(y − xTopt ) for some α : Θ → R+ which is integrable (with respect
214 to |µopt |). To this end, without loosing generality, let us write φ as

φ(θ) = α(θ)(y(θ) − xTopt (θ)) + w(θ), (3.8)

⊥
where w(θ) ∈ y(θ) − xTopt (θ) , while α is an arbitrary (not necessarily non-negative) scalar

215

216 function. Without loss of generality, we can assume that supp φ ⊂ supp µopt , i.e., we can assume
217 that α and w vanish on Θ\supp µopt . Taking into account that y −xTopt is continuous on Θ and that
218 y − xTopt = ε on supp µopt , we get that α(θ) = y(θ) − xTopt (θ), φ(θ) /ε2 for every θ ∈ supp µopt .
This ensures that α is |µopt |-measurable, and then w = φ − α y − xTopt is also |µopt |-measurable.


219

Define υ = ψ y − xTopt , where ψ ∈ C(Θ) is an arbitrary non-negative, scalar test function. υ is

220

221 therefore an eligible test function in (3.7), implying

Z
α(θ)(y(θ) − xTopt (θ)) + w(θ) , ψ(θ) y(θ) − xTopt (θ) d |µopt | (θ)



0 ⩽ ⟨µopt , υ⟩M(Θ)n ,C(Θ)n =
Θ
Z D E
2 2
= ψ(θ) y(θ) − xTopt (θ) α(θ) d |µopt | (θ) = α♭ |µopt |, ψ y − xTopt .
Θ M(Θ),C(Θ)

222 As |µopt | is a non-negative measure, by varying ψ and taking again into account that y − xTopt = ε
223 on supp µopt , it follows α ⩾ 0 for |µopt |-a.e. θ. Taking ψ = 1, we also get that α is |µopt |-integrable.
2 2
224 Furthermore, (3.2) together with ∥uopt ∥L2 (0,T )m = ∥Φ∗T µopt ∥L2 (0,T )m = µopt , xTopt M(Θ)n ,C(Θ)n
225 leads to
ε ∥µopt ∥M(Θ)n = µopt , y − xTopt M(Θ)n ,C(Θ)n .

At this stage we exploit the form (3.8) and put µopt = α(y − xTopt ) + w ♭ |µopt | in the last


226

227 expression. By taking into account that w(θ) is orthogonal to y(θ) − xTopt (θ), and recalling that

8
228 y − xTopt = ε on supp µopt , we get
Z Z
2
ε α(θ)(y(θ) − xTopt (θ)) + w(θ) d |µopt | (θ) = α(θ) (y(θ) − xTopt (θ)) d |µopt | (θ)
Θ Θ
Z
=ε α(θ) (y(θ) − xTopt (θ)) d |µopt | (θ) (3.9)
Θ

229 As α is a non-negative function and α(y − xTopt ) + w ⩾ α y − xTopt , the equality (3.9) implies
230 α(y − xTopt ) + w = α y − xTopt for |µopt |-a.e. θ, that is to say that w = 0 for |µopt | (and conse-
231 quently for µopt ) a.e. θ.
232 Remark 3.4. The last result covers the well-known behavior of the optimal control in the deter-
233 ministic case, i.e., when Θ is a singleton. In that case, unless the optimal control is trivial, the
234 solution of the dual problem is a vector having the same direction as y − xTopt (see e.g. [19]).
235 In general, when a functional of the form (3.1) is defined on a Hilbert space, its minimizer hopt is
236 either trivial or a solution to the Euler-Lagrange equation which formally reads as:
hopt
y − xTopt = ε .
∥hopt ∥

237 In that case it follows directly that the minimizer has the same direction as y − xTopt . This kind of
238 result is stronger than the one obtained in the last proposition, and it corresponds to a constant
239 function α in (3.6). The reason why such approach is not possible in our setting is that the norm
240 term ∥µ∥M(Θ)n entering the functional J is not differentiable.
241 As a consequence of the last proposition, we obtain the following characterization of the optimal
242 measure in the scalar case.
243 Corollary 3.5. Using the assumptions and notations introduced in Theorem 3.1, let us suppose
244 n = m = 1.
245 Then the control operator does not change sign, i.e., there exists ϵ ∈ {−1, 1} and β ∈ C(Θ) such
246 that β(θ) > 0 and B(θ) = ϵβ(θ) for every θ ∈ Θ. 
247 Set Θ+ = θ ∈ Θ | y(θ) − xTopt (θ) = ε and Θ− = θ ∈ Θ | y(θ) − xTopt (θ) = −ε . Then
248 supp µopt ⊂ Θ− ∪ Θ+ and µopt is nonnegative (respectively, nonpositive) on Θ+ (respectively Θ− ).
249 In addition, if Θ− = ∅ (respectively Θ+ = ∅), then either uopt ≡ 0 or sign(uopt (t)) = ϵ
250 (respectively sign(uopt (t)) = −ϵ) for every t ∈ [0, T ].
251 Proof. The fact that B ∈ C(Θ) does not vanish directly follows from Proposition 2.4. Being a
252 continuous function, its sign is constant.
253 The support and sign conditions of µopt directly follow from (3.3) and Proposition 3.3, respectively.
254 For the sign of uopt , let us assume that Θ− = ∅ (the proof for other case is identical). We have for
255 every t ∈ [0, T ],
Z Z
(T −t)A∗ (Θ) ∗
uopt (t) = e B(θ) dµopt (θ) = ϵ e(T −t)A (Θ) β(θ) dµopt (θ),
Θ Θ+

∗
256 We conclude by noticing that, for every θ ∈ Θ, e(T −t)A (Θ) β(θ) > 0 and µopt is nonnegative. In
257 particular, uopt ≡ 0 if and only if µopt = 0 and, according to Theorem 3.1, this holds if and only if
258 ∥y∥C(Θ)n ⩽ ε. Otherwise, the optimal control has a constant sign equal to ϵ for every t ∈ [0, T ].
259 The rest of this subsection provides results by which one can construct a family of optimization
260 problem (P ) for which the minimizers of (3.1) can be explicitly calculated. This will allow us to

9
261 compare the results obtained by the numerical procedure to be presented in the Section 6 with the
262 exact ones and to verify efficiency and correctness of the proposed algorithms (see Section 7).
263 The idea is to invert the problem: we specify a positive measure µ e and calculate x e,
eT = ΦT Φ∗T µ
264 the final state of the system run by the control u ∗
e = ΦT µ e. The adjusted target y is then carefully
265 chosen as a center of a square whose one vertex is xeT (note that the geometrical representation of
n
266 a ball in C(Θ) is actually a square). In particular, in accordance to Proposition 3.3, the difference
267 y−x eT has the same direction as µ e for every parameter value θ.
268 Proposition 3.6. Using the assumptions and notations introduced in Theorem 3.1 take ε > 0,
269 η ∈ M(Θ) a nonnegative scalar measure, and f ∈ C(Θ)n such that |f | does not vanish. Define
270 µ eT = ΦT Φ∗T µ
e = f ♭ η, x e and y = x e is the minimizer of J defined by (3.1).
eT + ε |ff | . Then µ
271 Proof. Let us first observe that for every µ ∈ M(Θ)n , we have
2
J(µ) = 1
2 ∥Φ∗T µ∥L2 (0,T )m − ⟨µ, y⟩M(Θ)n ,C(Θ)n + ε ∥µ∥M(Θ)n
 D E 
2
= 21 ∥Φ∗T µ∥L2 (0,T )m − ⟨Φ∗T µ, Φ∗T µ
e⟩L2 (0,T )m + ε ∥µ∥M(Θ)n − µ, |ff |
M(Θ)n ,C(Θ)n

272 and Z 
⟨f, f ⟩
Z
2 2
µ) = − 12 ∥Φ∗T µ
J(e e∥L2 (0,T )m + ε dη = − 12 ∥Φ∗T µ
|f | dη − e∥L2 (0,T )m .
Θ Θ |f |
273 Let µopt ∈ M(Θ)n be the minimizer of J, we have J(µopt ) ⩽ J(e µ), that is to say that
 D E 
1 ∗ 2 f
2 ∥ΦT (µopt − µ
e)∥L2 (0,T )m ⩽ −ε ∥µopt ∥M(Θ)n − µopt , |f |
n n
.
M(Θ) ,C(Θ)
D E
274 Obviously, we have ∥µopt ∥M(Θ)n ⩾ µopt , |ff | , from which we can conclude that we
M(Θ)n ,C(Θ)n
275 necessarily have e) = 0. We conclude that µopt = µ
Φ∗T (µopt
−µ e using the unique continuation
276 property (2.6) of the adjoint.
277 As a direct corollary, we have the following result.
278 Corollary 3.7. Using the assumptions and notations introduced in Theorem 3.1. Suppose n = 1,
279 let ε > 0, ϵ ∈ {−1, 1}, η ∈ M(Θ) a nonnegative measure, y = ϵ(ΦT Φ∗T η + ε). The minimizer
280 of (3.1) is ϵη.
281 Proof. This is Proposition 3.6 with f = ϵ.
282 From Proposition 3.6 and Corollary 3.7, we see that there is no reason that the optimal mea-
283 sure µopt is sparse, i.e., the Lebesgue measure of its support is in general not zero. This might look
284 surprisingly at the first glance, as µopt is a minimizer of the functional containing the ∥·∥M(Θ)n
285 term. This non-sparsity pitfall will be discussed in Section 4 in the case where the parameter set
286 is of a finite cardinality.

287 3.2 Continuity properties of the minimizers

288 The characterization of the solution through the adjoint problem, that is obtained in the previous
289 subsection, allows us also to demonstrate the continuity of solutions to optimization problem (P )
290 with respect to the tolerance ε, the time T and the adjusted target y. More precisely, we consider
291 a sequence of functionals
2
Jk (µ) = 1
2 Φ∗Tk µ L2 (0,T )m
− ⟨µ, yk ⟩M(Θ)n ,C(Θ)n + εk ∥µ∥M(Θ)n (µ ∈ M(Θ)n ). (3.10)

292 where:

10
293 (i) yk → y strongly in C(Θ)n ,
294 (ii) Tk → T in R∗+
295 (iii) εk → ε in R∗+ .
296 By Theorem 3.1, for each k ∈ N, the functional Jk admits the unique minimizer µk satisfying
297 the relation
2
Φ∗Tk µk L2 (0,T )m = ⟨µk , yk ⟩M(Θ)n ,C(Θ)n − εk ∥µk ∥M(Θ)n (3.11)
k

298 and the control uk = Φ∗Tk µk is the solution of the optimal ensemble control problem

1 2
min ∥u∥L2 (0,T )m ,
2
m (Pk )
u ∈ L2 (0, T ) and ∥yk − ΦTk u∥L∞ (Θ)n ⩽ εk .

299 We are going to prove that (µk )k and (uk )k converge respectively vaguely and strongly (in the
300 appropriate spaces) to the limit points µ and u, respectively, where µ is the minimizer of the limit
301 functional (3.1) and u is the solution of the ensemble optimal control problem (P ).
302 Note that the convergence of time horizons Tk implies strong convergence of the operators ΦTk
303 and Φ∗Tk . The following technical result characterizes strong convergence of arbitrary operators on
304 Banach spaces.
305 Lemma 3.8. Let (Lk )k be a sequence of linear operators in L(X, Y ), with X and Y being Banach
306 spaces. The following statements are equivalent.
307 (a) Lk −→ L strongly in L(X, Y ), i.e.,

(∀ u ∈ X) Lk u −→ Lu in Y.

308 (b) For every weakly ∗ converging sequence (vk )k in Y ′ it holds

∗ ∗
vk −⇀ v =⇒ L∗k vk −⇀ L∗ v in X ′ .

309 (c) For every strongly converging sequence (uk )k in X it holds

uk −→ u =⇒ Lk uk −→ Lu in Y.

310 Proof. (a) ⇒ (b): Let (vk )k converge weakly ∗ to v ∈ Y ′ . By employing the strong convergence of
311 the operators Lk , for every u ∈ X we get

⟨L∗k vk , u⟩X ′ ,X = ⟨vk , Lk u⟩Y ′ ,Y −→ ⟨v, Lu⟩Y ′ ,Y = ⟨L∗ v, u⟩X ′ ,X ,

∗
312 i.e., L∗k vk −⇀ L∗ v in X ′ .
313

314 (b) ⇒ (c): We first prove the weak convergence of the sequence (Lk uk )k . To this end, note that
315 for an arbitrary v ∈ Y ′ we have

⟨v, Lk uk ⟩Y ′ ,Y = ⟨L∗k v, uk ⟩X ′ ,X −→ ⟨L∗ v, u⟩X ′ ,X = ⟨v, Lu⟩Y ′ ,Y .

316 Here we have employed the strong convergence of (uk ) and the assumption (b). In particular, we
317 have
∥Lu∥Y ⩽ lim inf ∥Lk uk ∥Y . (3.12)
k→∞

318 On the other hand, for each k there exists vk ∈ Y ′ such that
2
∥Lk uk ∥Y = ⟨vk , Lk uk ⟩Y ′ ,Y , (3.13)

11
319 and ∥vk ∥Y ′ = ∥Lk uk ∥Y (see e.g. [6, Corollary 1.3.]).
320 Due to the weak convergence of (Lk uk )k , the sequence of norms ∥Lk uk ∥Y is bounded, in sequel
321 implying boundedness of (vk )k in Y ′ . In particular, up to a subsequence, (vk )k converge weakly ∗
322 to some v ∈ Y ′ , and it holds

∥v∥Y ′ ⩽ lim inf ∥vk ∥Y ′ = lim inf ∥Lk uk ∥Y . (3.14)

k→∞ k→∞

323 From here it follows

2
∥Lk uk ∥Y = ⟨vk , Lk uk ⟩Y ′ ,Y = ⟨L∗k vk , uk ⟩X ′ ,X −→ ⟨L∗ v, u⟩X ′ ,X = ⟨v, Lu⟩Y ′ ,Y , (3.15)

324 where we passed to the limit by using the assumption (b) again.
325 In particular, by employing (3.15) and (3.14), we obtain
2
lim ∥Lk uk ∥Y = Y ′⟨ v, Lu ⟩Y ⩽ ∥v∥Y ′ ∥Lu∥Y ⩽ lim ∥Lk uk ∥Y ∥Lu∥Y ,
k→∞ k→∞

326 implying limk→∞ ∥Lk uk ∥Y ⩽ ∥Lu∥Y . Together with (3.12) this finishes the proof.
327

328 (c) ⇒ (a): Trivial (take a constant sequence uk = u).

329 Before proving the convergence, we first demonstrate the boundedness properties of the mini-
330 mizers.
331 Lemma 3.9. The sequence (µk )k of minimizers of functionals Jk given by (3.10) is bounded
332 in M(Θ)n .
333 Proof. We argue by contradiction and assume that

∥µk ∥M(Θ)n −→ ∞ (3.16)

334 (up to a non-relabelled subsequence).

335 Dividing the relation (3.11) by ∥µk ∥M(Θ)n and denoting νk = µk / ∥µk ∥M(Θ)n we get

2
∥µk ∥M(Θ)n Φ∗Tk νk L2 (0,T )m
= ⟨νk , yk ⟩M(Θ)n ,C(Θ)n − εk . (3.17)

336 As all the terms on the right-hand side are bounded (with respect to k), the contradictory assump-
337 tion (3.16) implies
m
Φ∗Tk νk → 0 strongly in L2 (0, T ) . (3.18)
338 On the other hand, denoting by ν the vague limit of (νk ) (up to a subsequence), the strong
m
339 convergence of operators ΦTk and Lemma 3.8 imply Φ∗Tk νk −⇀ Φ∗T ν weakly in L2 (0, T ) . Together
340 with (3.18), we get ΦT ν = 0, while the ensemble controllability assumption (2.5) implies ν = 0.
∗

341 Going back to (3.17) and passing to the limit we get

0 > εk − ⟨νk , yk ⟩M(Θ)n ,C(Θ)n −→ ε.

342 Here we use that ⟨νk , yk ⟩M(Θ)n ,C(Θ)n is the dual product of a strongly convergent functions with
343 measures converging vaguely to zero, thus vanishing on the limit.
344 The last relation provides the desired contradiction which, in turn, implies the result.
345 Based on the last two lemmas, we are going to obtain the main result of this paragraph.

12
346 Proposition 3.10. The sequence of minimizers µk of the functionals Jk given by (3.10) converges
n
347 vaguely in M(Θ) to µopt (the minimizer of the limit functional (3.1)), and ∥µk ∥M(Θ)n converges
348 to ∥µopt ∥M(Θ)n . In addition, the sequence of optimal ensemble controls uk = Φ∗Tk µk converges
m
349 strongly in L2 (0, T ) to the solution of the limit problem (P ).
350 Proof. By Lemma 3.9, there exists a (non-relabelled) subsequence of the minimizers µk converging
n
351 vaguely to some measure µL ∈ M(Θ) .
352 In the first part of the proof, we show that the limit µL is the unique minimizer of the limit
353 functional (3.1). Thus, it suffices to show
n
J(µL ) ⩽ J(µ), µ ∈ M(Θ) . (3.19)
354 To this end, we explore the convergence properties of the sequence of minimal functional values
 
2
lim inf Jk (µk ) = lim inf 12 Φ∗Tk µk L2 (0,T )m − ⟨µk , yk ⟩M(Θ)n ,C(Θ)n + εk ∥µk ∥M(Θ)n . (3.20)
k→∞ k→∞

355 Due to the weak lower semi-continuity of norms in Banach spaces, the last term is bounded from
356 below by ε ∥µL ∥M(Θ)n .
357 The second term on the right-hand side of (3.20) is the product of a strongly and vaguely
358 convergent sequence, thus ⟨µk , yk ⟩M(Θ)n ,C(Θ)n → ⟨µL , y⟩M(Θ)n ,C(Θ)n .
359 As for the first term entering functional Jk by Lemma 3.8 and the strong convergence of the
360 operators ΦTk , we have Φ∗Tk µk −⇀ Φ∗T µL weakly in L2 (0, T )n , implying
2 2
lim inf Φ∗Tk µk L2 (0,T )m
⩾ ∥Φ∗T µL ∥L2 (0,T )m . (3.21)
k→∞

361 Consequently, from (3.20) we get

2
lim inf Jk (µk ) ⩾ 1
2 ∥Φ∗T µL ∥L2 (0,T )m − ⟨µL , y⟩ + ε ∥µL ∥M(Θ)n = J(µL )
k→∞

362 As µk minimizes Jk , for an arbitrary µ ∈ M(Θ)n we have

lim inf Jk (µ) ⩾ J(µL ). (3.22)
k→∞

363 On the other hand, employing the assumed convergences, (i), (ii) and (iii), for each fixed µ we have
2
Jk (µ) = 1
2 Φ∗Tk µ
L2 (0,T )m
− ⟨µ, yk ⟩M(Θ)n ,C(Θ)n + εk ∥µ∥M(Θ)n
(3.23)
1 ∗ 2
→ 2 ∥ΦT µ∥L2 (0,T )m − ⟨µ, y⟩M(Θ)n ,C(Θ)n + ε ∥µ∥M(Θ)n = J(µ).
∗
364 Together with (3.22) this implies (3.19), and in particular µk −⇀ µopt .
365 In the second part of the proof, we will show that the sequence ∥µk ∥M(Θ)n converges
366 to ∥µopt ∥M(Θ)n .
367 Taking µ = µopt in (3.23) we have,
J(µopt ) = lim Jk (µopt ) ⩾ lim Jk (µk ).
k→∞ k→∞

368 Combining the obtained inequality with (3.22) (remember that µopt = µL ) we get the convergence
369 of the optimal values Jk (µk ) → J(µopt ). Together with (3.21) this convergence implies,
lim inf ∥µk ∥M(Θ)n ⩽ ∥µopt ∥M(Θ)n .
k→∞

370 which together with the vague convergence of the sequence (µk ) implies the desired result.
371 The strong convergence of optimal ensemble controls uk = Φ∗Tk µk follows now directly from the
372 strong convergence of the operators ΦTk and the vague convergence of the measures µk (cf. Lem-
373 ma 3.8 (c)).

13
374 4 Systems with a finite number of parameter values
375 In this section, we restrict the analysis to the case of a finite number of parameters, i.e., card Θ < ∞.
376 Such situation deserves a special attention, not only because in some applications the parameter
377 assumes only discrete and bounded values, but also because the most general situations (i.e., those
378 in which the parameter ranges over a continuous set of values) are numerically treated and analyzed
379 as the limit case of discretized problems as the number of elements in the parameter set goes to
380 infinity.
381 The aim is to provide a geometrical description of the problem and its solution, which will
382 allow us to discuss the conditions under which the optimal measure (i.e., the minimizer of the
383 functional J given by (3.1)) is sparse.
384 To this end, we assume that the parameter dependent system (2.1) is uniformly ensemble
385 controllable, and that we deal with a parameter set Θ of finite cardinality. As we consider finite
386 dimensional systems only, these assumptions imply that the Gramian operator ΛT = ΦT Φ∗T ∈
n n
387 L(M(Θ) ; C(Θ) ) is not just injective (due to Proposition 2.7), but also of full range. Indeed,
n n
388 due to the assumed finite cardinality of the parameter set, both spaces M(Θ) and C(Θ) are
389 finite dimensional and are isomorphic to RnK (where K = card Θ ∈ N∗ denotes the number of
390 parameters), thus the statement follows from the rank-nullity theorem.
391 For c ⩾ 0, we introduce the subset Ec ⊂ C(Θ)n consisting of all the final states xT that can be
392 reached with a control of norm less or equal to c, i.e.,
n m
o
Ec := xT ∈ C(Θ)n | ∃u ∈ L2 (0, T ) such that ∥u∥L2 (0,T )m ⩽ c and xT = ΦT u .

393 As the minimal norm control reaching any state xT is of the HUM form Φ∗T µ for some µ ∈ M(Θ)
394 (see e.g. [29]), the introduced set can be characterized as follows
n o
Ec = xT ∈ C(Θ)n | ∃µ ∈ M(Θ)n such that xT = ΛT µ and ∥Φ∗T µ∥L2 (0,T )m ⩽ c ,

395 where ΛT stands for the Gramian ΦT Φ∗T . As the ensemble controllability assumption ensures the
396 injectivity of the Gramian, the optimal control norm can be rewritten as
2
∥Φ∗T µ∥L2 (0,T )m = ⟨µ, ΦT Φ∗T µ⟩M(Θ)n ,C(Θ)n = Λ−1 T
T x ,x
T
M(Θ)n ,C(Θ)n
,

397 leading to the following characterization of the set Ec :

n o
Ec = xT ∈ C(Θ)n | Λ−1 T
T x ,x
T
n
M(Θ) ,C(Θ) n ⩽ c2 . (4.1)

398 In other words, Ec is a sublevel set of a smooth function F (xT ) = Λ−1 T

T x ,x
T
. In particular,
399 {Ec , c ⩾ 0} is a nested family of closed, convex sets with smooth boundary ∂Ec = {xT | F (xT ) =
400 c2 }. The normal functional at the point xT ∈ ∂Ec is nEc (xT ) = ∇F (xT ) = 2Λ−1
T x .
T

401 Theorem 4.1. Let xT , x1 ∈ C(Θ)n be such that x1 (θ) − xT (θ) ⩽ ε for every θ ∈ Θ.
402 Then xT is the optimal final state of the problem (P ) with the target x1 and initial datum x0 = 0
403 if and only if the normal functional nEc (xT ) is of the form

nEc xT = α(x1 − xT ) ♭ ν, (4.2)

404 where α ∈ C(Θ) and ν ∈ M(Θ) are respectively a non-negative scalar function and a non-negative
405 measure, and α(θ) = 0 whenever x1 (θ) − xT (θ) < ε.
406 Here Ec is the sublevel set defined by (4.1) with the constant c = Φ∗T Λ−1
T x
T
L2 (0,T )m
being the
407 minimal norm required to reach xT .

14
408 Remark 4.2. By the choice of the level c, it follows that xT ∈ ∂Ec . If additionally the relation (4.2)
409 holds, i.e., if xT is the optimal final state of the problem (P ), then c = Φ∗T Λ−1
T x
T
L2 (0,T )m
is the
410 smallest level such that the intersection of the sublevel set Ec with the target box is non-empty
(see Figure 1).

nEc (xT )
1

0.5

x1
0
−2 −1.5 −1 −0.5 0 0.5 1 1.5
xT
−0.5

Ec

Figure 1: Set Ec for n = 1 and K = 2, together with the target ball, the optimal final state xT
and the normal functional nEc (xT ). Each x ∈ C(Θ) is presented by a point (x(θ1 ), x(θ2 )) ∈ R2 .
For this plot, we have set A(θ1 ) = 1, A(θ2 ) = 5, B(θ1 ) = 1, B(θ2 ) = 1/50, T = 1, ε = 1/4, and
the plotted set Ec is for c = 1.
411

412 Proof of Theorem 4.1. =⇒ Let us assume that xT is the optimal final state of the problem (P ).
413 In that case 21 nEc (xT ) = Λ−1
T x
T
is the optimal measure and the result follows directly from
414 Proposition 3.3 (ensuring that it is of the form α(x1 − xT )|µopt |, with α and |µopt | being a non-
415 negative scalar function and measure, respectively) and the description of the support of the
416 optimal measure (3.3).
417 ⇐= Assume that (4.2) holds. In order to prove the statement, it is enough to show that xT is the
418 only element in the intersection of the sublevel set Ec and the ε ball around the target x1 ∈ C(Θ)n .
419 We assume the contrary, i.e., that there exists some x eT ̸= xT lying in the above intersection.
420 We split the parameter set into two disjoint subsets, defined as

Θb = θ ∈ Θ | x1 (θ) − xT (θ) = ε and Θi = θ ∈ Θ | x1 (θ) − xT (θ) < ε .

 

421 By the assumption made on α, Θi is a zero measure set for nEc (xT ). On the other hand, for any
422 θ ∈ Θb , the vector xT (θ) lies on the boundary of the ε ball (in Rn ) around the target x1 (θ), which,
423 together with the fact that x eT (θ) belongs to the same ball, implies that

x1 − xT (θ), x eT − xT (θ) ⩾ 0 (4.3)

(θ ∈ Θb ).

424 Consequently, we get

n Ec x T , x

T
e − xT = nEc xT , x

T
M(Θ)n ,C(Θ)n
e − xT M(Θ1 )n ,C(Θ1 )n
D E
= α x1 − xT eT − xT



ν, x (4.4)
♭ M(Θ1 )n ,C(Θ1 )n
1 T T T



= α♭ ν, x − x · x e −x M(Θ1 )n ,C(Θ1 )n
⩾ 0,

425 where the final inequality follows by (4.3) and the assumed non-negativity of the measure α♭ ν.

15
426 As nEc (xT ) is the normal functional, the only way for x
eT ∈ Ec ∩ B∞ (x1 , ε) to satisfy the
e − x M(Θ)n ,C(Θ)n = 0, i.e., that x
inequality (4.4) is that nEc x , x eT belongs to the tangential
T

T T
427

428 space TEc (xT ). As the convex set Ec intersects the tangential space TEc (xT ) only in its boundary
429 points, it follows that x
eT ∈ ∂Ec . In such way we obtain two distinct points x eT and xT lying in
430 ∂Ec ∩ B∞ (x , ε). However, that would imply existence of two different feasible controls of the same
1

431 norm, c, which contradicts the uniqueness of the solution to problem (P ).

432 Remark 4.3. Note that the relation (4.2) implies that the optimal measure

µopt = Λ−1 T 1 T
(4.5)


T x = 2 n Ec x

433 has, in every point θ ∈ Θ, the same direction as (x1 − xT )(θ), which is in accordance with the
434 results of Proposition 3.3.

435 The last result allows an alternative, more direct proof of Proposition 3.6, in the case where
436 the parameter set is of finite cardinal, by which one can construct a series of problems for which
437 the exact solution is known.
438 Alternative proof of Proposition 3.6. With the notations and assumptions of Proposition 3.6, let
439 c = Φ∗T Λ−1
T xeT L2 (0,T )m be the minimal norm required to reach the state x
eT . Then for the normal
440 functional nEc (e
xT ) we have
 
|f |
eT = 2Λ−1 eT = 2f ♭ η = 2 eT



n Ec x T x y−x ♭ η,
ε

441 where η ∈ M(Θ) and f ∈ C(Θ)n are the nonnegative scalar measure and the continuous vector
442 function from the statement of Proposition 3.6 such that Λ−1 eT = f ♭ η.
T x
443 From here we see that the normal nEc (e
xT ) is of the form (4.2) (note that y = x1 for the zero
444 initial datum). The statement now follows from Theorem 4.1.
445 The introduced sets Ec allow us to describe the set of targets for which the associated optimal
446 control vanishes for some parameter θ ∈ Θ. This allows us to construct examples for which the
447 exact solution is known. This is important in numerical implementations, where the detection of
448 active parameters plays a crucial role in constructing the optimal measure.
449 For simplicity, the discussion is provided for scalar systems, i.e., we assume n = 1.
450 Given some optimal final state xT , we know that the corresponding optimal measure is of the
451 form (4.5). In particular, its i-th component (with i ∈ {1, . . . , K} and K = card Θ) is going to be
452 zero measure if and only if xT = ΛT µopt lies in the hyperplane Hi = [Λ1 , . . . , Λi−1 , Λi+1 . . . ΛK ]
453 spanned by all the columns of the Gramian matrix except the i-th one. This happens exactly for
454 the states xT for which the normal derivative nEc (xT ) = 2µopt is orthogonal to the i-th coordinate
455 axes (Figure 2). √
456 As any target x1 lies at most K 2 ε/2 from the optimal final state, this implies that the set of
457 targets Ni for which the i-th component
√ of the optimal measure equals zero is contained within the
458 hyperstrip Hiε := Hi + B(0, K 2 ε/2). More precise characterization of this set can be obtained in
459 the case K = 2, i.e., when the parameter obtains only two values. In that case, the set of targets
460 for which the i-th component of the optimal measure is zero is given by

nEc (xT )
   
Ni = x T + ε + βe i | x T
∈ Hi , |β| ⩽ 1 ∪ B∞ (0, ε),
|nEc (xT )|

461 and depicted on Figure 3.

16
 6

0
−4 −2 0 2 4

H2 −2

−4
H1

−6

Figure 2: The hyperplanes H1 and H2 for n = 1 and K = 2. Arrows represent nEc (xT ) for
different values of c and well-chosen values of xT . The parameters used for this plot are the ones
of Figure 1.

0
−4 −2 0 2 4
−1
N2
−2

−3

N1

Figure 3: Sets N1 and N2 for n = 1 and K = 2. The parameters used for this plot are the ones
of Figure 1.

17
462 For numerical implementations, it would be very useful to know whether the given target
463 belongs to some set Ni . This is because revealing the support of the optimal measure is often the
464 most challenging part of the numerical process, while, as we shall see in Section 7 (see in particular
465 Examples 7.1 and 7.2), tuning the weights of active Dirac distributions turns out to be easier.
466 However, the description of the sets Ni requires computation of the Gramian operator, which is a
467 computationally very demanding task (therefore construction of the Gramian is generally avoided
468 and instead various iterative procedures are used for solving control optimization problems). For
469 this reason, the above description of targets for which the optimal measure vanishes at some points
470 cannot directly accelerate numerical procedures. However, it contributes to a better understanding
471 of the problem and its solution, and allows us to describe some characteristics of the optimal
472 measures, in particular the sparsity.
473 In general, we can say that for the targets that are not contained in the union of hyperstrips
474 ∪i Hiε the optimal measure is going to be fully supported. Of course, as hyperstrips Hiε , i ∈
475 {1, . . . , K} cover just a smaller portion of the whole output space RK , this will be the case for
476 the most of the targets. This might look surprisingly, having in mind that the term ε ∥µ∥M(Θ)n
477 entering the dual functional (3.1) is supposed to imply the sparsity in the support of the optimal
478 measure. However, in the discrete case, the norm ∥·∥M(Θ)n is equivalent to the ℓ1 norm (see
479 comments preceding Theorem 3.1), while the size of the measure support (i.e., number of active
480 Dirac distributions) is given by the ℓ0 norm. Although the ℓ1 norm might be used to recover sparse
481 solutions under suitable assumptions (see e.g. [7]), in general the ℓ0 term is the one that induces
482 sparsity.

483 5 Discretized problem and convergence

484 In this section, we show that the solution to the problem (P ) can be arbitrarily well approximated
485 by a solution to a discretized problem taking into account only a finite number of parameter
486 values. This is of crucial importance from the numerical point of view, as the minimization of the
n
487 dual functional in the latter case (originally defined on M(Θ) ) reduces to the minimization of a
488 function on R Kn
for some large enough K ∈ N . ∗

489 To this end, let h > 0 and consider Θh a compact subset of Θ such that dist(θ, Θh ) ⩽ h for
490 every θ ∈ Θ. The discretized version of (P ) is
1 2
min 2 ∥u∥L2 (0,T )m ,
Z T
m
u ∈ L2 (0, T ) and sup x1 (θ) − eT A(θ) x0 (θ) − e(T −t)A(θ) B(θ)u(t) dt ⩽ ε.
θ∈Θh 0
(Ph )
491 Following [18], we have the following results.
n
492 Theorem 5.1. Assume that (2.1) is uniformly ensemble controllable. Given y ∈ C(Θ) and
m
493 ε > 0, we set uopt ∈ L2 (0, T ) (respectively uh ∈ L2 (0, T )) the unique minimizer of (P ) (respec-
n n
494 tively (Ph )). Let us also set µopt ∈ M(Θ) and µh ∈ M(Θh ) the unique minimizers of J defined
n n
495 by (3.1) over M(Θ) and M(Θh ) respectively.
m
496 We have ∥uh ∥L2 (0,T )m ⩽ ∥uopt ∥L2 (0,T )m , for every h > 0, and in addition, uh → uopt in L2 (0, T )
497 and µh vaguely converges to µopt as h → 0.
498 Before giving the proof of this result, let us first recall that according to Proposition 2.4,
499 if (2.1) is uniformly ensemble controllable with parameter set Θ, then (2.1) is uniformly ensemble
500 controllable with parameter set Θh . Hence, the existence and uniqueness of uh and µh (and of uopt
501 and µopt ) directly follow from Theorem 3.1. Let us also mention that the strong convergence

18
502 of (uh )h has already been given in [18, §2.2], under the additional assumption that A, B and y are
503 Lipschitz function with respect to θ. It seems that this regularity assumption is not required, and
504 we give this proof below.
505 Proof of Theorem 5.1. We identify µh to a measure in M(Θ) supported on Θh . Since Θh ⊂ Θ, it
506 is obvious that ∥uh ∥L2 (0,T )m ⩽ ∥uopt ∥L2 (0,T )m .
Step 1. Let us show that for every sequence (hk )k∈N ∈ (R∗+ ) converging to 0, ∥µhk ∥M(Θ)n is
N
507

508 uniformly bounded with respect to k.

Assume by contradiction that there exist a sequence (hk )k∈N ∈ (R∗+ ) converging to 0 such that
N
509

510 ∥µhk ∥M(Θ)n → ∞ as k → ∞. We then set λk = ∥µhk ∥M(Θ)n , and we assume without loss of
511 generality that λk > 0 for every k ∈ N, and we set ςk = µhk /λk . By vague compactness (see
512 e.g. [2, Chapter IV]), there exist a subsequence (still denoted by (ςk )k∈N ) such that (ςk )k is vaguely
n
513 convergent to ς0 ∈ M(Θ) . According to eq. (3.2), we have
2
λk ∥Φ∗T ςk ∥L2 (0,T ) = ⟨ςk , y⟩M(Θ)n ,C(Θ)n − ε ∥ςk ∥M(Θ)n = ⟨ςk , y⟩M(Θ)n ,C(Θ)n − ε.

514 This ensures that ∥Φ∗T ςk ∥L2 (0,T ) → 0 as k → ∞, i.e., Φ∗T ς0 = 0, and hence, ς0 = 0 by Propo-
2
515 sition 2.7. Hence, for k large enough, we have, λk ∥Φ∗T ςk ∥L2 (0,T ) = ⟨ςk , y⟩M(Θ)n ,C(Θ)n − ε < 0,
2
516 which leads to a contradiction with λk ∥Φ∗T ςk ∥L2 (0,T ) ⩾ 0.
Step 2. Let us consider a sequence (hk )k∈N ∈ (R∗+ ) such that hk → 0 as k → ∞, and assume
N
517
n
518 that (µhk )k∈N converges vaguely in M(Θ) . We aim to prove that (µhk )k∈N converges vaguely
519 to µopt , and that the corresponding sequence of controls (uhk )k∈N converges strongly to uopt in
m
520 L2 (0, T ) .
m
521 Let us denote by µ0 the vague limit of (µhk )k and we set u0 = Φ∗T µ0 ∈ L2 (0, T ) . For every
m
522 φ ∈ L (0, T ) , we have as k → ∞,
2

⟨uhk , φ⟩L2 (0,T )m = ⟨µhk , ΦT φ⟩M(Θ)n ,C(Θ)n −→ ⟨µ0 , ΦT φ⟩M(Θ)n ,C(Θ)n = ⟨u0 , φ⟩L2 (0,T )m .

523 This ensures the weak convergence of (uhk )k∈N to u0 , in particular, we have ∥u0 ∥L2 (0,T )m ⩽
524 lim inf k→∞ ∥uhk ∥L2 (0,T )m ⩽ ∥uopt ∥L2 (0,T )m .
525 Since A and B are continuous on the compact set Θ, there are uniformly continuous and
526 bounded. In addition, it is easy to see that there exists a constant C > 0, depending only on T
527 and on the bounds of A and B, such that

e(T −t)A(θ1 ) B(θ1 ) − e(T −t)A(θ0 ) B(θ0 ) ⩽ C (|A(θ1 ) − A(θ0 )| + |B(θ1 ) − B(θ0 )|)
(t ∈ [0, T ], θ0 , θ1 ∈ Θ).
m
528 This ensures that for every v ∈ L2 (0, T ) and every θ0 , θ1 ∈ Θ we have
√
|(ΦT v) (θ1 ) − (ΦT v) (θ0 )| ⩽ T C ∥v∥L2 (0,T )m (|A(θ1 ) − A(θ0 )| + |B(θ1 ) − B(θ0 )|) .

529 The last inequality clearly implies the uniform continuity of ΦT uhk , more precisely, for every δ > 0,
530 there exists η(δ) > 0 (with η(δ) → 0 as δ → 0) such that for every k ∈ N and every θ0 , θ1 ∈ Θ
531 satisfying |θ1 − θ0 | ⩽ δ, we have

|(ΦT uhk ) (θ1 ) − (ΦT uhk ) (θ0 )| ⩽ η(δ).

532 By taking larger η(δ) if needed, we can also assume that

|y(θ1 ) − y(θ0 )| ⩽ η(δ).

19
533 Indeed, y is also uniformly continuous on Θ, since y is continuous on the compact set Θ.
534 For every θ ∈ Θ and every k ∈ N, let us set θhk ∈ Θhk such that |θ − θhk | ⩽ hk . We recall that
535 since θhk ∈ Θhk , we have |(ΦT uhk )(θhk ) − y(θhk )| ⩽ ε. We then have, for every θ ∈ Θ and every
536 k ∈ N,
|(ΦT u0 ) (θ) − y(θ)| ⩽ |(ΦT u0 ) (θ) − (ΦT uhk ) (θ)| + |(ΦT uhk ) (θ) − (ΦT uhk ) (θhk )|
+ |(ΦT uhk ) (θhk ) − y(θhk )| + |y(θhk ) − y(θ)|
⩽ |(ΦT u0 ) (θ) − (ΦT uhk ) (θ)| + 2η(hk ) + ε.
537 Taking the limit k → ∞, we conclude that |(ΦT u0 ) (θ) − y(θ)| ⩽ ε. This ensures that u0 solves the
538 uniform ensemble controllability problem with parameter set Θ. Thus, we have ∥u0 ∥L2 (0,T )m ⩾
539 ∥uopt ∥L2 (0,T )m . But, we already know that ∥u0 ∥L2 (0,T )m ⩽ ∥uopt ∥L2 (0,T )m from the weak conver-
540 gence of the sequence (uhk )k . This ensures that ∥u0 ∥L2 (0,T )m = ∥uopt ∥L2 (0,T )m , and by unique-
541 ness of the minimal L2 -norm control, u0 = uopt . Finally, from ∥uopt ∥L2 (0,T )m ⩾ ∥uhk ∥L2 (0,T )m
542 and ∥u0 ∥L2 (0,T )m ⩽ lim inf k→∞ ∥uhk ∥L2 (0,T )m , we conclude that limk→∞ ∥uhk ∥L2 (0,T )m =
m
543 ∥uopt ∥L2 (0,T )m . This ensures that uhk → uopt strongly in L2 (0, T ) as k → ∞.
544 The fact that u0 = uopt also ensures that Φ∗T µ0 = Φ∗T µopt and hence, µ0 = µopt by Propo-
545 sition 2.7.
546 Conclusion. Assume by contradiction that µh is not vaguely convergent to µopt as
547 h → 0. Then, there exists f ∈ C(Θ)n such that ⟨µh , f ⟩M(Θ)n ,C(Θ)n does not converge to
⟨µopt , f ⟩M(Θ)n ,C(Θ)n as h → 0. In particular, there exists r > 0 and a sequence (hk )k∈N ∈ (R∗+ )
N
548

549 converging to 0 such that ⟨µhk − µopt , f ⟩M(Θ)n ,C(Θ)n > r. But according to step 1, (µhk )k∈N
550 is uniformly bonded. Using vague compactness, up to the extraction of a subsequence (µhk )k is
551 vaguely convergent. Finally, the step 2 ensures that the only possible vague limit is µopt . This
552 leads to a contradiction, and concludes the proof.
553 Remark 5.2. Let us mention that the arguments used in the first step of the proof of Theorem 5.1
554 are similar to the ones used in Lemma 3.9. This type of similarity is expected since in Theorem 5.1,
555 we prove continuity of the minimizers of J with respect to the parameter set Θ, while Lemma 3.9
556 is an ingredient to prove continuity with respect to T , ε and y (see Proposition 3.10).
557 Let us now consider a discretized version of the problem. To this end, we consider K ∈ N∗ ,
558 θ1 , . . . , θK ∈ Θ be two by two distinct, ΘhK = {θ1 , . . . , θK }, with hK = maxθ∈Θ dist(θ, ΘhK ). Let
PK
559 µhK be the minimizer of J defined by (3.1) over ΘhK , then µhK = k=1 wk δθk , where δθk is the
K
560 atomic mass located at θk , and w = (wk )k∈{1,...,K} ∈ (Rn ) is the unique minimizer of
K 2 K K
Z T
1 X ∗ (T −t)A(θk )∗
X X
JΘhK (w) = B(θk ) e wk dt − ⟨wk , y(θk )⟩ + ε |wk | (5.1)
2 0 k=1 k=1 k=1

561 and we have,

K
X ∗ ∗
uhK (t) = Φ∗T µhK = B(θk ) e(T −t)A(θk ) wk (t ∈ [0, T ]).
k=1
562 In order to prepare § 6.1, let us give the following results.
K
563 Proposition 5.3. w = (wk )k∈{1,...,K} ∈ (Rn ) is a minimizer of JΘhK if and only if for every
564 l ∈ {1, . . . , K}, wl ∈ Rn is a minimizer of
* Z T +
1
b ∈ Rn 7→ w
Jl : w b∗ Λθl w
b − w,
b y(θl ) − e(T −t)A(θl ) B(θl )ul (t) dt + ε |w|
b ,
2 0

20
565 with
K Z T
X ∗ ∗ ∗ ∗
ul (t) = B(θk ) e(T −t)A(θk ) wk and Λθl = e(T −t)A(θl ) B(θl )B(θl ) e(T −t)A(θl ) dt.
k=1 0
k̸=l

566 Proof. According to Theorem A.1, w = (wk )k∈{1,...,K} is a minimizer of JΘhK if and only if for
(
wk , if k ̸= l,
567 every l ∈ {1, . . . , K}, wl is a minimizer of w n
e with w
b ∈ R 7→ JΘhK (w) ek = .
w,
b otherwise
568 But, for every l ∈ {1, . . . , K}, we have,

K 2 K K
1 T X
Z
∗ (T −t)A(θk )∗
X X
JΘhK (w) = B(θk ) e wk dt − ⟨wk , y(θk )⟩ + ε |wk |
2 0
k=1 k=1 k=1

1 T
Z 2
∗ ∗
= B(θl ) e(T −t)A(θl ) wl dt
2 0
Z T* K
+
∗ (T −t)A(θl )∗ ∗ (T −t)A(θk )∗
X
+ B(θl ) e wl , B(θk ) e wk dt − ⟨wl , y(θl )⟩ + ε |wl |
0 k=1
k̸=l
2
Z T K K K
1 X ∗ ∗ X X
+ B(θk ) e(T −t)A(θk ) wk dt − ⟨wk , y(θk )⟩ + ε |wk | .
2 0 k=1 k=1 k=1
k̸=l k̸=l k̸=l

569 It is then clear that minimizers of w e are minimizers of Jl .

b 7→ JΘhK (w)

570 Let us now give condition on the minimizers of Jl defined in Proposition 5.3.
571 Proposition 5.4. Let Λ ∈ Rn×n , be a positive definite matrix, z ∈ Rn and ε > 0, the minimizer
572 popt ∈ Rn of J : p ∈ Rn 7→ 12 p∗ Λp − ⟨p, z⟩ + ε |p| over Rn is unique and satisfies popt = 0 if |z| < ε,
573 and otherwise popt = ρq where,
574 • |q| = 1;
575 • shrink(|z| , ε)/CΛ ⩽ ρ ⩽ shrink(|z| , ε)/cΛ ;
576 • (εI + ρΛ)q = z,
577 where 0 < cΛ ⩽ CΛ are such that
2 2
cΛ |p| ⩽ p∗ Λp ⩽ CΛ |p| (p ∈ Rn )

578 and 
s + µ
 if s < −µ,
shrink(s, µ) = 0 if − µ ⩽ s ⩽ µ, (s ∈ R, µ ∈ R+ ).

s−µ if µ < s


579 Remark 5.5. When n = 1, we have popt = shrink(z, ε)/Λ.

580 Proof. First, J is strictly convex and coercive, hence admits a unique minimizer popt on Rn .
581 Furthermore, popt is characterized by 0 ∈ ∂J (popt ). But, we have
(
p
Λp − z + ε |p| if p ̸= 0,
∂J (p) =
B(0, ε) − z if p = 0,

21
582 recall that B(0, ε) ⊂ Rn is the closed ball of Rn centered on 0 and of radius ε. In particular, we
583 have 0 ∈ ∂J (0) if |z| ⩽ ε. Let us then assume that |z| > ε, we then have popt ̸= 0, and we set
584 ρ = |popt | and q = popt /ρ. We thus have (εI + ρΛ)q = z, and hence ρ > 0 shall be such that
2
2 −1
585 1 = |q| = (εI + ρΛ) z . This leads to the bounds on ρ.

586 Before entering the algorithmic computation of a minimizer of JΘhK , let us give the sub-
587 differential of JΘhK . For every k ∈ {1, . . . , K}, we have,
 Z T
wk
y(θ ) − e(T −t)A(θk ) B(θk )u(t) dt + ε , if wk ̸= 0,

 k

∂wk JΘhK (w) = Z0 T |wk|

y(θk ) − e(T −t)A(θk ) B(θk )u(t) dt + B (0, ε) , if wk = 0,



0
PK ∗ (T −t)A(θk )∗
588 with u(t) = k=1 B(θk ) e wk . In particular, if 0 ∈ ∂wk JΘhK (w), we have
Z T !
(T −t)A(θk )
shrink y(θk ) − e B(θk )u(t) dt , ε = 0.
0

589 In addition, when wk = 0, we have

Z T !
(T −t)A(θk )
shrink y(θk ) − e B(θk )u(t) dt , ε = min {|z| | z ∈ ∂wk J(w)} .
0

590 Alternatively to the algorithm proposed in § 6.1, one can use a direct optimization based on
591 the optimality conditions given in Proposition 3.3. We have the following direct corollary of Pro-
592 position 3.3.
K
593 Corollary 5.6. w = (wk )k∈{1,...,K} ∈ (Rn ) is a minimizer of JΘhK if there exists (αk )k∈{1,...,K} ∈
594 RK such that for every k ∈ {1, . . . , K},
595 • αk ⩾ 0, and αk = 0 if |y(θk ) − (ΦT Φ∗T µhK ) (θk )| < ε;
596 • wk = αk (y(θk ) − (ΦT Φ∗T µhK ) (θk )),
PK
597 where we have set µkK = k=1 wk δθk .
598 These condition will lead to the minimization problem given in § 6.2.

599 6 Numerical strategies

600 In this section, we propose three ways to numerically compute a minimizer of JΘh .
601 The first one, S.0 is basically using existing solvers without taking in consideration differentia-
602 bility issues that emerge from non-smoothness of the dual functional. The second one, S.1 is based
603 on the optimality conditions given in Corollary 5.6, see § 6.2. The last ones, S.2.0 and S.2.1, use
604 the greedy coordinate descent method proposed in [22], see § 6.1.
605 Finally, the strategies S.2.0 and S.2.1 can also be combined with the strategies S.0 and S.1. All
606 the possibility are given in § 6.3, and compared on some numerical examples in Section 7.

607 6.1 Greedy coordinate descent

608 Inspired from [22], we propose the following greedy coordinate descent algorithm (Algorithm 1).
609 This algorithm is a coordinate descent algorithm, but instead of considering each coordinate se-
610 quentially, we select at each iteration the coordinate for which the decrease of the cost function is
611 the most important.
612 The convergence of Algorithm 1 is ensured by the following result.

22
 K
613 Proposition 6.1. Let f1 , . . . , fK ∈ C(Rn ) and F ∈ C 1 ((Rn ) ) be all convex. We consider the cost
614 function
K
X
J(w1 , . . . , wK ) = F (w1 , . . . , wK ) + fk (wk ) (w1 , . . . , wK ∈ Rn ),
k=1

615 and we assume that J is coercive and strictly convex.

K
616 Given w0 = (w10 , . . . , wK
0
) ∈ (Rn ) , we define the sequences (ki )i∈N ∈ {1, · · · , K}N and (wi )i∈N ∈
 N
K
617 (Rn ) by:
 
ki ∈ argmin minn J(w1i , . . . , wk−1
i
, w, i
b wk+1 i
, . . . , wK ) , (6.1)
k∈{1,...,K} w∈R
b

618 and

wki+1
i
∈ argmin J(w1i , . . . , wki i −1 , w,
b wki i +1 , . . . , wK
i
),
w∈R
b n

wki+1 = i
wk (∀k ∈ {1, . . . , K} \ {ki }),
K
619 Then limi→∞ wi = wopt ∈ (Rn ) , where wopt is the minimizer of J.
620 The proof of this result follows the arguments given in [22, Section 2.4] and is given in Appen-
621 dix A for the sake of completeness. Note that in [22] the proof is given for n = 1.
622 Remark 6.2. We expected that the convergence result would also hold by selecting

ki ∈ argmax dist 0, ∂wk J(wi ) . (6.2)

k∈{1,...,K}

623 But, we did not succeed to prove this result without additional regularity assumptions on the
624 functions fk . However, in practice, it seems that using the selection rule (6.2), give similar results
625 compared to the selection rule (6.1), but is much faster (see examples in Section 7).

Algorithm 1 Greedy coordinate descent for uniform ensemble control.

Require: T > 0, ε > 0, (A, B, x0 , x1 ) ∈ C(Θ)n×n × C(Θ)n×m × C(Θ)n × C(Θ)n , with n, m ∈ N∗
Require: Θh ⊂ Θ with K := card(Θh ) ∈ N∗ and Θh = {θ1 , . . . , θK }
Require: η > 0 ▷ A tolerance parameter.
Z T
m
Ensure: u ∈ L2 (0, T ) is such that sup x1 (θ) − eT A(θ) x0 (θ) − e(T −t)A(θ) B(θ)u(t) dt ⩽ ε+δ
θ∈Θh 0

23
 Algorithm 1 (continued)
1: for all k ∈ {1, . . . , K} do
2: wk ← 0
3: yk ← x1 (θk ) − eT A(θk ) x0 (θk )
4: Gk ← shrink(|yk | , ε) ▷ Gk = dist {0, ∂wk JΘh (w)}.
5: end for
K
X ∗ ∗
6: u(t) ← 0 ▷ u(t) = B(θk ) e(T −t)A(θk ) wk .
k=1
7: c←0 ▷ c = JΘh (w).
8: while sup |yk | > ε + η do
k∈{1,...,K}
9: for all l ∈ {1, . . . , K} do ▷ Look for the best descent coordinate.
10: if Gl ̸= 0 then
11: for all k ∈ {1, . . . , K} do w bk ← wk
12: end for ∗
∗
13: ub(t) ← B(θl ) e(T −t)A(θl ) wl
Z T
14: ybl ← yl + e(T −t)A(θl ) B(θl )b
u(t) dt
0
Z T K
X ∗ ∗
▷ We have ybl = x1 (θl ) − eT A(θl ) x0 (θl ) − e(T −t)A(θl ) B(θl ) B(θk ) e(T −t)A(θk ) wk dt.
0 k=1
k̸=l
15: wbl ← Approximate Control(ε, T, A(θl ), B(θl ), ybl ) ▷ See Algorithm 2.
16: cl ← JΘh (w)
b b
17: cl ← c
else b
18: end if
19: end for
20: l ← an element of argmax b
ck ▷ l is the best coordinate to update.
k∈{1,...,K}
∗ ∗
21: ub(t) ← B(θl ) e(T −t)A(θl ) wl
Z T
22: ybl ← yl + e(T −t)A(θl ) B(θl )b
u(t) dt
0
23: wl ← Approximate Control(ε, T, A(θl ), B(θl ), ybl )
K
∗ (T −t)A(θl )∗
X ∗ ∗
24: u(t) ← u(t) − u
b(t) + B(θl ) e wl ▷ We have u(t) = B(θk ) e(T −t)A(θk ) wk .
k=1
25: c ← JΘh (w)
26: for all k ∈ {1, . . . , K} do
Z T
27:
1
yk ← x (θk ) − e T A(θk ) 0
x (θk ) − e(T −t)A(θk ) B(θk )u(t) dt
0
28: if wk = 0 then Gk ← shrink(|yk | , ε)
wk
29: else Gk ← yk − ε
|wk |
30: end if
31: end for
32: end while
33: return u(t)

626 Remark 6.3. Here are some comments on Algorithm 1.

627 • The loop, lines 9 to 19, can be done in parallel.

24
628 • What is computed in the loop, lines 9 to 19, can be stored and reused lines 21 to 23 and 25.
629 • Algorithm 2 will be called many times, and hence has to be fast.
630 • This algorithm is adapted from [22], where the aim is to find sparse solutions. Thus, as we
631 will see on some examples in Section 7, will lead to sparse minimizers. Consequently, it might
632 not be optimal when we are in the situation described in Proposition 3.6 and Corollary 3.7,
633 i.e., when the optimal measure is fully supported. However, we shall see that the controls
634 obtained in this way still provide a fairy good approximation of the optimal one.

Algorithm 2 Computation of an approximate control for the system ẋ = Ax + Bu.

Require: T > 0, ε > 0 and (A, B, y) ∈ Rn×n × Rn×m × Rn with n, m ∈ N∗
Z T
∗
Ensure: w ∈ Rn is such that y − e(T −t)A Bu(t) dt ⩽ ε with u(t) = B ∗ e(T −t)A w
0
1: function Approximate Control(ε, T, A, B, y)
2: if |y| ⩽ ε then w ← 0 ▷ See Remark 5.5.
3: else ▷ In that case, we have w ̸= 0.
4: if n = 1 Zthen
T
∗
5: Λ← e(T −t)A BB ∗ e(T −t)A dt
0
6: w ← shrink(y, ε)/Λ ▷ See Remark 5.5.
7: else if the control Gramian is “computable” then ▷ e.g. for n “small”.
Z T
∗
8: Λ← e(T −t)A BB ∗ e(T −t)A dt
0
9: (r, q) ← solution of r > 0, q ∈ Rn , |q| = 1 and (εI + rΛ)q(r) = y
10: w ← rq(r) ▷ See Proposition 5.4.
11: else !
1 T ∗ (T −t)A∗ 2
Z
12: w ← argmin B e ν dt − ⟨ν, y⟩ + ε |ν|
ν∈Rn 2 0
13: end if
14: end if
15: return w
16: end function

635 Remark 6.4. • In Algorithm 2 (line 9), finding a solution of r > 0, q ∈ Rn , |q| = 1 and
636 (εI + rΛ)q(r) = y can be done using a dichotomy search.
637 • To find a minimizer of J (line 12 of Algorithm 2), one can use a gradient descent method.
638 In fact, note that J is only not differentiable at 0, and note that 0 cannot be the minimizer
639 in the situation of line 12 this is because |y| > ε. A way to avoid the evaluation of J (ν) =
1 T ∗ (T −t)A∗ 2
ν dt − ⟨ν, y⟩ + ε |ν| at 0 is to initialize the gradient descent method with
R
640
2 0 B e
641 ν ∈ Rn such that J (ν) < J (0) = 0. This is for instance possible by setting ν = αy, with
2 RT ∗ 2
642 α = shrink(|y| , ε |y|)/C, with C = 0 B ∗ e(T −t)A y dt. We thus have ν ̸= 0 and J (ν) < 0,
643 see Remark 5.5.
644 • Computing the Gramian matrix might be resource consuming. Hence, instead of computing
645 it in Algorithm 2, it can be computed and stored in Algorithm 1.

25
646 6.2 Reduced minimization problem
647 For this strategy, we will use the optimality condition given in Corollary 5.6. More precisely, we
648 consider the following minimization problem (recall that we have set y(θ) = x1 (θ) − eT A(θ) x0 (θ)),

K 2
1 X ∗
min B(θk )∗ e(T −·)A(θk ) wk (6.3a)
2
k=1 L2 (0,T )m

649 with respect to w1 , . . . , wK ∈ Rn , subject to the following constraints


T 2 2
 y(θk ) − y (θk ) ⩽ ε ,
 
2
∀k ∈ {1, ..., K}, ∃αk ⩾ 0 such that: αk y(θk ) − y T (θk ) − ε2 ⩾ 0, (6.3b)


wk = αk y(θk ) − y T (θk ) ,


650 where we have set

K
!
Z T
∗ (T −t)A(θl )∗
X
T
y (θk ) = e (T −t)A(θk )
B(θk ) B(θl ) e wl dt. (6.3c)
0 l=1

651 Observe that the last three constraints

PKcoincide with the optimality condition given in Corollary 5.6,
∗
652 and the control is given by u(t) = k=1 B(θk )∗ e(T −t)A(θk ) wk . Using Corollary 5.6, it is obvious
653 that minimizers of JΘhK coincide with minimizers of (6.3).
654 Remark 6.5. In practice the minimization problem (6.3) can be expressed as a minimization only
655 on (αk )k∈{1,...,K} , however, from a numerical point of view, it is easier to consider the minimiza-
656 tion problem (6.3) as a minimization with respect to (wk , αk , y(θk ))k∈{1,...,K} since gradient and
657 Hessians can be easily computed.

658 6.3 Proposed numerical strategies

659 Based on the previous consideration, we can propose the following numerical strategies.
660 S.0. Use existing nonlinear solvers, and try to directly minimize JΘhK , without considering po-
661 tential issues related to the non-differentiability of the cost functional. In the examples of
662 Section 7, we use the fminunc function of Matlab with the quasi-Newton algorithm.
663 S.1. Use existing nonlinear solvers to find a minimizer of (6.3). In the examples of Section 7, we
664 use the fmincon function of Matlab with the interior-point algorithm.
665 In addition, one can use algorithms based on Algorithm 1. This leads to two strategies:
666 S.2.0. Use directly Algorithm 1.
667 S.2.1. Use Algorithm 1 but instead of using the selection rule given by (6.1), we use the one given
668 by (6.2) (see Remark 6.2).
669 As we will see in Section 7, with strategies S.2.0 and S.2.1, it appears that many steps are between
670 non-zeros values. Hence, we propose a refinement of these two algorithms by running S.0 or S.1 after
671 some component of (wk )k switches from null to a non-null value, or vice versa. More precisely,
672 using notations of Algorithm 1, if wl was 0 and is updated to a non-null value or if wl was
673 non-null and is updated to 0, we run S.0 or S.1 for the same minimization problem, but with
674
e = {θk ∈ Θh , ∀k ∈ {1, . . . , K} s.t. wk ̸= 0} in place of Θh . In other words, we try to update all
Θ
675 the non-zero values of (wk )k in one shoot. This leads to the four additional strategies, S.2.i–S.j
676 (with i, j ∈ {0, 1}).
677 Let us recall that we have analytical guaranties of convergence only for the strategy S.2.0.

26
678 7 Numerical examples
679 Unless explicitly said otherwise, for all the following examples
680 • the final time is T = 1;
681 • Θ = [1, 2], and Θh is a uniform discretization of Θ with step h = 10−2 , leading to K =
682 card Θh = 101;
683 • ε = 5 × 10−2 ;
684 • the time integrals are computed using the trapezoidal rule with a uniform time step of 10−4 .
685 In particular, to compute ΦT , we use the trapezoidal rule, and the matrix exponentials are
686 computed numerically. Another alternative could have been to solve the ordinary differential
687 equation using a numerical scheme;
688 • the stopping criterion on the gradient (i.e., δ in Algorithm 1) is set to 10−12 , we also set a
689 step tolerance of 10−20 (i.e., step sizes shall be greater than this number) and a maximal
690 number of iterations of 105 .
691 Examples 7.1 and 7.2 correspond to situation where the optimal control and the adjoint measure
692 are known (see Proposition 3.6 and Corollary 3.7). Example 7.3 corresponds to a situation where
693 the target cannot be reached (see Appendix B). Finally, Example 7.4 tackled the uniform ensemble
694 controllability for a discretized heat equation.
695 In the tables used for the comparison of the different numerical strategies, G = (G1 , . . . , GK ) ∈
696 (R+ )K is given by lines 27 to 29 of Algorithm 1 that is Gk = dist (0, ∂wk JΘh (w)).
697 The codes are run using Matlab version R2020b on a 12 core processor of 64 bits and 3.20 GHz.

698 Example 7.1. We consider a first basic example, where n = m = 1, A(θ) = −θ and B(θ) = 1.
699 It is well-known that this system is ensemble controllable (see [21, 27]). We consider x0 (θ) = 0,
700 dµopt (θ) = χω (θ) dθ, where ω is an open subset of Θ. According to Corollary 3.7, we know that
701 uopt = Φ∗T µopt is the minimizer of (P ) when the target is set to x1 (θ) = (ΦT Φ∗T µopt ) (θ) + ε.
702 In practice, Psince Θh (introduced in Section 5) is chosen as a uniform discretization of Θ, we
703 set µopt,h = h θ∈Θh ∩ω δθ , and it is clear that µopt,h is vaguely convergent to µopt . Also, the
704 target x1h is set accordingly, i.e., x1h (θ) = (ΦT Φ∗T µopt,h ) (θ) + ε, and the discretized optimal control
705 is uopt,h = Φ∗T µopt,h .
706 The numerical simulations are run with ω = (4/3, 7/4) (we recall that the other parameters
707 are given in the preamble of Section 7). A comparison of the different strategies listed in § 6.3 is
708 given in Table 1.

CPU running nb. ∥uh − uopt,h ∥L2 ∥µh − µopt,h ∥M

Strategy max G
time (s) time (s) iter. ∥uopt,h ∥L2 ∥µopt,h ∥M
S.0 6.07 0.63 7 2.553×10− 5 1.561×10−3 0.959
S.1 4.59 2.82 83 1.293×10−13 4.384×10−8 1.271
S.2.0 1409.58 201.93 100000 2.387×10− 8 2.467×10−5 1.952
S.2.0–S.0 39.29 8.49 11 8.161×10−10 4.432×10−6 1.952
S.2.0–S.1 11.19 2.47 14 7.452×10−15 6.833×10−9 1.902
S.2.1 198.18 20.01 100000 3.389×10− 7 1.140×10−3 1.800
S.2.1–S.0 1.06 0.19 14 3.043×10− 9 1.617×10−5 1.905
S.2.1–S.1 4.28 1.19 12 3.757×10−14 2.951×10−8 1.857
(a) Errors, time of computation and number of iterations.

Table 1: Comparison of the strategies listed in § 6.3 for Example 7.1.

27
 Strategy Stopping reason
S.0 Cannot decrease the objective function along the current search direction.
S.1 Local minimum possible. Constraints satisfied.
S.2.0 Maximal number of iteration reached.
S.2.0–S.0 Cost is not decreasing.
S.2.0–S.1 Cost is not decreasing.
S.2.1 Maximal number of iteration reached.
S.2.1–S.0 Cost is not decreasing.
S.2.1–S.1 Cost is not decreasing.
(b) Stopping reasons.

Table 1: (continued)

709 In Table 1b, for strategies S.2.i–S.j (i, j ∈ {0, 1}) the last cost modifications is of order 10−17 ,
710 which is close to the numerical 0 (indeed, by denoting copt the minimal cost, we have numerically
711 copt + 10−18 = copt ). Let us also mention that for all simulations, we have x1 (θ) − x(T, θ) = ε
712 up to an error of order 10−8 . We also see that the strategy S.2.0–S.1 gives the best results and
713 strategy S.2.1–S.0 is the fastest one. As claimed in § 6.3, for strategy S.2.0 and S.2.1, most of the
714 time is spent in updating non-zero values in fact for S.2.0, for all 105 steps the selected index is
715 in {1, 55, 101}, and for S.2.1, after some iterations the selected index loops between 1, 42 and 101
716 and then between 42, 69 and 101. In Table 1a, we see that even if the control is closed to the
717 optimal ones, the adjoint is rather far from the optimal one, this is also illustrated on Figure 4. As
718 we can see from Figure 4 and Table 1a, wide variation on the adjoint does not necessarily produce
719 large variations on the control. One can see from Figure 4, that the strategies S.0 and S.1 lead to
720 adjoints that are fully supported on Θh , while the other strategies lead to adjoints for which the
721 support consist in only few points of Θh . These facts will be discussed later on.

0.14
0.01 µopt µopt
µ µ
0.12
0.008
0.1

0.006 0.08

0.06
0.004

0.04
0.002
0.02

0 0
1 1.2 1.4 1.6 1.8 2 1 1.2 1.4 1.6 1.8 2
θ θ
(a) Strategy S.0. (b) Strategy S.1.
0.4 µopt 0.4 µopt µopt
µ µ 0.2 µ
0.35 0.35

0.3 0.3
0.15
0.25 0.25

0.2 0.2 0.1

0.15 0.15

0.1 0.1 0.05

0.05 0.05

0 0 0
1 1.2 1.4 1.6 1.8 2 1 1.2 1.4 1.6 1.8 2 1 1.2 1.4 1.6 1.8 2
θ θ θ
(c) Strategy S.2.0. (d) Strategy S.2.0–S.0. (e) Strategy S.2.0–S.1.

Figure 4: Adjoints obtained with strategies listed in § 6.3 for Example 7.1.

28
 0.2
µopt µopt µopt
µ µ 0.2 µ
0.25

0.15 0.2
0.15

0.15
0.1
0.1

0.1

0.05 0.05
0.05

0 0 0
1 1.2 1.4 1.6 1.8 2 1 1.2 1.4 1.6 1.8 2 1 1.2 1.4 1.6 1.8 2
θ θ θ
(f) Strategy S.2.1. (g) Strategy S.2.1–S.0. (h) Strategy S.2.1–S.1.

Figure 4: (continued)

722 Example 7.2. We consider the system and data given in Example n 7.1, except that
o we set µopt =
723 δθe1 + 2δθe2 . In that case, for i ∈ {1, 2}, we set θei,h ∈ argmin θei − θ , θ ∈ Θh and µopt,h =
724 δθe1,h + 2δθe2,h .
725 The numerical simulation is run with θe1 = 7/6 and θe2 = 17/10 (we recall that the other
726 parameters are given in the preamble of Section 7). A comparison of the different strategies listed
in § 6.3 is given in Table 2.

CPU running nb. ∥uh − uopt,h ∥L2 ∥µh − µopt,h ∥M

Strategy max G
time (s) time (s) iter. ∥uopt,h ∥L2 ∥µopt,h ∥M
S.0 4.11 0.44 19 2.431×10− 8 2.421×10−5 1.960
S.1 2.92 1.78 76 2.058×10−13 2.675×10−8 1.981
S.2.0 1478.55 213.11 100000 7.575×10− 7 1.065×10−4 2.000
S.2.0–S.0 0.55 0.10 10 4.783×10− 9 1.999×10−5 2.000
S.2.0–S.1 5.60 1.38 11 2.553×10−12 3.665×10−7 2.000
S.2.1 194.10 19.49 100000 2.495×10− 6 3.447×10−4 2.000
S.2.1–S.0 0.28 0.04 5 2.188×10− 8 4.409×10−5 2.000
S.2.1–S.1 3.01 0.78 8 9.378×10−12 5.331×10−7 2.000
(a) Errors, time of computation and number of iterations

Strategy Stopping reason

S.0 Cannot decrease the objective function along the current search direction.
S.1 Local minimum possible. Constraints satisfied.
S.2.0 Maximal number of iteration reached.
S.2.0–S.0 Step size is less than the step tolerance.
S.2.0–S.1 Step size is less than the step tolerance.
S.2.1 Maximal number of iteration reached.
S.2.1–S.0 Cost is not decreasing.
S.2.1–S.1 Cost is not decreasing.
(b) Stopping reasons.

Table 2: Comparison of the strategies listed in § 6.3 for Example 7.2.

727

728 We see form Table 2 that almost all the conclusion made for Example 7.1 also holds for Exam-
729 ple 7.1, except that here the strategy S.1 give slightly better results than the strategy S.2.0–S.1. As
730 for Example 7.1, we display on Figure 5 the adjoints numerically computed. The comments made

29
731 in Example 7.1 also applies here, wide variation on the adjoint does not necessarily produce large
732 variations on the control, and the strategies S.0 and S.1 lead to adjoints that are fully supported
733 on Θh , while the other strategies lead to adjoints for which the support consist in only few points
of Θh .

2 µopt 2 µopt
µ µ

1.5 1.5

1 1

0.5 0.5

0 0
1 1.2 1.4 1.6 1.8 2 1 1.2 1.4 1.6 1.8 2
θ θ
(a) Strategy S.0. (b) Strategy S.1.
µopt µopt 2 µopt
µ µ µ
2 2

1.5
1.5 1.5

1
1 1

0.5 0.5
0.5

0 0 0
1 1.2 1.4 1.6 1.8 2 1 1.2 1.4 1.6 1.8 2 1 1.2 1.4 1.6 1.8 2
θ θ θ
(c) Strategy S.2.0. (d) Strategy S.2.0–S.0. (e) Strategy S.2.0–S.1.
2 µopt µopt 2 µopt
µ µ µ
2

1.5 1.5
1.5

1 1
1

0.5 0.5 0.5

0 0 0
1 1.2 1.4 1.6 1.8 2 1 1.2 1.4 1.6 1.8 2 1 1.2 1.4 1.6 1.8 2
θ θ θ
(f) Strategy S.2.1. (g) Strategy S.2.1–S.0. (h) Strategy S.2.1–S.1.

Figure 5: Adjoints obtained with strategies listed in § 6.3 for Example 7.1.
734

735 Based on the results given in Examples 7.1 and 7.2, we can make the following comments.
736 • All strategy lead to some acceptable control.
737 • Strategies S.0 and S.1 lead to adjoint measure, which are everywhere non-zero. In fact, they
738 are not adapted for sparsity, contrary to strategies S.2.0 and S.2.1 and their derivatives. In
739 addition, for strategy S.1 since an interior-point algorithm is used, the fact that the obtained
740 solution is non-zero everywhere is expected.
741 • Even if the strategy S.2.1–S.1 did not give the best results, it is relatively fast and provide
742 good results for the control. We will then use this strategy for the following examples.

743 Example 7.3. We consider the system and data given in Example 7.1, but we set x1 (θ) = 1. In
744 that case, it is known from Proposition B.1 that the target state in not reachable. Hence, it is
745 expected that the norm of the minimal L2 -norm control goes to ∞ as ε goes to 0.
746 This fact is illustrated on Figure 6. On Figure 7, we also display the obtained results for
747 ε = 2 × 10−4 (we recall that the other parameters are given in the preamble of Section 7). In

30
748 particular, Figures 7b and 7d confirm the claim of Corollary 3.5, i.e., the sign of the optimal measure
749 coincide with the sign of x1 − xT and the measure is only supported on the set of parameters θ such
750 that x1 (θ) − xT (θ) = ε. The results for this example have been obtained using strategy S.2.1–S.1.

5
7 × 105
4
6 × 105
3
5 × 105
4 × 105
2 3 × 105
2 × 105
1
105 2
1.8
0 0 1.6
10−6 10−5 10−4 10−3 10−2 10−1 100 10−6 10−5 1.4 θ
ε 10−4 10−3 1.2
ε 10−2 10−1
100 1
(a) Norm of the control (black) and number of Dirac
masses in the adjoint state (gray), with respect to ε. (b) |µh | with respect to ε.

Figure 6: Norm of the control, number of Dirac masses in the adjoint state and absolute value of
the adjoint state, with respect to ε, for Example 7.3.

14 u ε x1 (·) − x(T, ·)

12 0.00015

10 0.0001

8 5 × 10−5

6 0

4 −5 × 10−5

2 −0.0001
0 −0.00015
−2 −ε
0 0.2 0.4 0.6 0.8 1 1 1.2 1.4 1.6 1.8 2
t θ

(a) Obtained control (∥u∥L2 (0,T ) ≃ 3.713092). (b) Error to the target.
1.4

1.2 µ
600
1
400
0.8
200

0.6 0

0.4 −200

0.2 −400

0 −600
0 0.2 0.4 0.6 0.8 1
−800
t
1 1.2 1.4 1.6 1.8 2
1 1.2 1.4 1.6 1.8 2 θ
θ

(c) Time dependent solutions for some values of θ. (d) Obtained adjoint measure.
On this graph, we plot x1 (θ) − x(·, θ).

Figure 7: Results obtained with strategy S.2.1–S.1 for the system given in Example 7.3, with
ε = 2 × 10−4 .

31
751 Example 7.4. For this example, we consider a discretized version of the family indexed by θ ∈ Θ
752 of systems of 1D heat equation:

ẏ(t, θ, x) = ∂x2 y(t, θ, x) − (θ − 1)y(t, τ, x) + u(t, x) (t > 0, x ∈ (0, 1)), (7.1a)

y(t, θ, 0) = y(t, θ, 1) = 0 (t > 0), (7.1b)
0
y(0, θ, x) = y (θ, x) (x ∈ (0, 1)). (7.1c)

753 One can refer to [9] for some ensemble controllability results related to parabolc systems.
754 Given n ∈ N∗ , using centered finite differences, the discretized version of (7.1) is

Ẏ (t, θ) = A(θ)Y (t, θ) + B(θ)U (t), Y (0, θ) = Y 0 (θ), (7.2)

755 with
−2
 
1 0 ... 0
.. .. .. .. 
. . . . 

 1
1 
.. .. ..

A(θ) = . . . − (θ − 1)In ∈ Rn×n and B(θ) = In ∈ Rn×n
 
(n + 1)2 
 0 0
.. .. .. ..

. . .
 
 . 1
0 ... 0 1 −2

756 and where Yi (t, θ) ≃ y(t, θ, n+2

i i
), Ui (t) ≃ u(t, n+2 ) and Yi0 (θ) = y 0 (θ, n+2
i
).
q
757 Let us mention that using the trapezoidal rule, ∥y(t, θ, ·)∥L2 (0,1) is approximated by 1
n+1 |Y (t, θ)|
758 hence, the value of ε will be adapted accordingly. For the control problem, we set y (θ, x) = 0

759 sin(x + θ) and null target, i.e., y 1 (θ) = 0. The goal is to find a control u ∈ L2 (0, T ; L2 (0, 1))
760 such that supθ∈Θ ∥y(T, θ, ·)∥L2 (0,1) ⩽ ε, ε, T, Θ and discretization parameters are defined in the
761 preamble of Section 7. For the numerical illustration, we use the strategy S.2.1–S.1, with n = 19.
762 The obtained results are displayed on Figure 8.

763 A Proof of Proposition 6.1

764 For this proof, we follow [22, Section 2.4]. We consider a cost function of the form
K
X
J(w1 , . . . , wK ) = F (w1 , . . . , wK ) + fk (wk ) (w1 , . . . , wK ∈ Rn ),
k=1

765 where f1 , . . . , fK are continuous convex function and F is a convex function of class C 1 . For every
K−1
766 k ∈ {1, . . . , K}, we define w qk = (w1 , . . . , wk−1 , wk+1 , . . . , wK ) ∈ (Rn ) , and Jk,wqk (wk ) = J(w)
767 and Fk,wqk (wk ) = F (w).
768 Theorem A.1. Using the notations and assumptions introduced above, w = (w1 , . . . , wk ) is a
769 minimizer of J if and only if for every k ∈ {1, . . . , K}, wk is a minimizer of Jk,wqk .
770 Proof. If w is a minimizer of J, then for every δ ∈ Rn and every k, we have

J(w1 , . . . , wk−1 , wk + δ, wk+1 , . . . , wk ) ⩾ J(w1 , . . . , wk−1 , wk , wk+1 , . . . , wK ),

771 That is to say that Jk,wqk (wk + δ) ⩾ Jk,wqk (wk ), i.e., wk minimizes Jk,wqk .
772 Reciprocally, if wk minimizes Jk,wqk , for every δ ∈ Rn and every t ∈ [0, 1], we have,

0 ⩽ Jk,wqk (wk + tδ) − Jk,wqk (wk ) = Fk,wqk (wk + tδ) − Fk,wqk (wk ) + fk (wk + tδ) − fk (wk ).

32
 8 0.06
6 0.04
10 4
8 2 0.06 0.02
6 0 0
−2 0.04
4 −0.02
−4 0.02 −0.04
2 −6
0 0 −0.06
−2 −0.02
−4 −0.04
1 2
−6 −0.06
0.8 1.8
−8
0 0.6 1.6
0.2 t 0
0.4 0.4 0.2 1.4 θ
0.6 0.2 0.4
x 0.8 x 0.6 1.2
1 0 0.8
1 1

(a) Obtained control (∥u∥L2 ((0,T )×(0,1)) ≃ 3.4857). (b) State at final time.
1.8

1.6
16000
∥µθ ∥L2 (0,1)
1.4
14000
1.2
12000
1

0.8 10000

0.6 8000

0.4 6000

0.2
4000
ε
0
0 0.2 0.4 0.6 0.8 1 2000
t
0
1 1.2 1.4 1.6 1.8 2
1 1.2 1.4 1.6 1.8 2 θ
θ

(c) L2 -norm if the time dependent solutions for some (d) L2 -norm of the obtained adjoint measure.
values of θ. On this graph, we plot ∥y(t, θ, ·)∥L2 (0,1) .

Figure 8: Results obtained with strategy S.2.1–S.1 for the system given in Example 7.4.

773 But, since fk is convex, we have fk (wk +tδ) = fk ((1 − t)wk + t(wk + δ)) ⩽ (1−t)fk (wk )+tfk (wk +
774 δ), and we get

fk (wk + tδ) − fk (wk ) ⩽ t (fk (wk + δ) − fk (wk )) (t ∈ [0, 1]).

775 This ensures

1
0⩽ (Fk,wqk (wk + tδ) − Fk,wqk (wk )) + fk (wk + δ) − fk (wk ) (t ∈ (0, 1)).
t
776 Since F is of class C 1 , Fk,wqk is also of class C 1 and taking the limit t → 0 in the above equation,
777 we get
0 ⩽ ⟨∇Fk,wqk (wk ), δ⟩ + fk (wk + δ) − fk (wk ).
778 Now, assume that w = (w1 , . . . , wK ) is such that wk minimizes Fk,wqk for every k ∈ {1, . . . , K}.

33
 K
779 For every δ = (δ1 , . . . , δK ) ∈ (Rn ) , we have, using the convexity (and C 1 regularity) of F ,
K
X
J(w + δ) − J(w) = F (w + δ) − F (w) + (fk (wk + δk ) − fk (wk ))
k=1
K
X
⩾ ⟨∇F (w), δ⟩ + (fk (wk + δk ) − fk (wk ))
k=1
K
X
= (⟨∇wk F (w), δk ⟩ + fk (wk + δk ) − fk (wk ))
k=1
K
X
= (⟨∇Fk,wqk (wk ), δk ⟩ + fk (wk + δk ) − fk (wk ))
k=1
⩾ 0.

780 That is to say that w is a minimizer of J.

781 Lemma A.2. Assume J is strictly convex and coercive. Then, for every k ∈ {1, . . . , K}, the
K−1
782 function gk : (Rn ) → Rn defined by gk (w)
q = argmin Jk,wq is continuous.

783 Proof. Firstly, since we assumed that J is strictly convex and coercive, we have that Jk,wq is also
K−1
784 strictly convex and coercive, this holds for every wq ∈ (Rn ) and every k ∈ {1, . . . , K}. This
785 ensures the existence and uniqueness of minimizers of Jk,wq , i.e., the function gk is well-defined.
 N
K−1 K−1
786 Let wq ∈ (Rn ) and (wqi )i∈N ∈ (Rn ) be a sequence converging to w, q for every
787 i ∈ N, we define wki = gk (w qi ) and wk = gk (w). q We then have, Jk,wq i (wki ) ⩽ Jk,wq i (wk ), and
788 limi→∞ Jk,wq i (wk ) = Jk,wq (wk ). This ensures that the sequence (Jk,wq i (wki ))i is bounded. But
789 since J is coercive (and recalling that Jk,wq i (wki ) = J(wi )), this implies that the sequence (wki )i is
790 bounded, and hence admits a convergent subsequence (still denote by (wki )i ), and we denote by
791 wk∞ its limit. We then have Jk,wq (wk∞ ) = limi→∞ Jk,wq i (wki ) ⩽ limi→∞ Jk,wq i (wk ) = Jk,wq (wk ),
792 this implies wk = wk∞ , by uniqueness of the minimizer of Jk,wq . We have thus proved that
793 limi→∞ gk (wqi ) = gk (w).
q

794 We are now in position to prove the convergence of the algorithm.

795 Proof of Proposition 6.1. It is trivial that

min J ⩽ J(wi+1 ) ⩽ J(wi ) (i ∈ N). (A.1)

(Rn )K

796 This ensures that (J(wi ))i is bounded and convergent. Due to the coercivity of J, we also have, by
797 compactness, the existence of ϕ : N → N increasing and w∞ ∈ (Rn )K such that limi→∞ wϕ(i) =
798 w∞ . The continuity of J ensures that J(w∞ ) = limi→∞ J(wϕ(i) ) = limi→∞ J(wi ) (the last equality
799 follows from the uniqueness of limits).
K
800 For every k ∈ {1, . . . , K} and every w ∈ (Rn ) , we define ∆Jk (w) = J(w1 , . . . , wK ) −
801 qk ), wk+1 , . . . , wK ) ⩾ 0, where gk is defined in Lemma A.2. From Lemma A.2,
J(w1 , . . . , wk−1 , gk (w
802 gk is continuous, hence ∆Jk is also continuous. We then have ∆Jk (w∞ ) = limi→∞ ∆Jk (wϕ(i) ).
803 We also have, by definition of the sequence,

0 ⩽ ∆Jk (wϕ(i) ) ⩽ ∆Jkϕ(i) (wϕ(i) ) = J(wϕ(i) ) − J(wϕ(i)+1 ) (i ∈ N, k ∈ {1, · · · , K}).

34
804 which goes to 0 as i goes to ∞. Hence, taking the limit i → ∞, we conclude that

0 = ∆Jk (w∞ ) (k ∈ {1, · · · , K}).

805 This ensures, using the definition of ∆Jk and Theorem A.1, that w∞ is the minimizer of J. Note
806 that the existence and uniqueness of the minimizer of J is ensured by the coercivity and strict
807 convexity of J.
808 Finally, the uniqueness of minimizers, the coercivity of J, and inequality (A.1) lead to
809 limi→∞ wi = w∞ .

810 B Some unreachable states

811 In this paragraph, we will show, in the case m = 1, that for every k ∈ N, θ ∈ Θ 7→ A(θ)k B(θ) ∈ Rn
812 is unreachable for the system (2.1). To get this result, we will assume that the system (2.1) is
813 uniformly ensemble controllable, and that Θ has a nonempty interior. Let us first define the set of
814 reachable points in time T > 0,
( )
Z T
RA,B (T ) = θ ∈ Θ 7→ e(T −t)A(θ) B(θ)u(t) dt | u ∈ L2 (0, T )m .
0

815 Before entering the proof

 of this result, let us mention that (2.1) is uniformly ensemble control-
816 lable if and only if Span θ ∈ Θ 7→ A(θ)k B(θ)v | k ∈ N, v ∈ Rm is dense in C(Θ). This fact can
817 be found in [28] (see also [11]).
818 Let us first prove our result in the case n = 1.
819 Proposition B.1. Let Θ ⊂ R be a compact set and assume that there exists θ0 ∈ Θ such that
820 θ0 is an accumulation point of Θ. Let a, b ∈ C(Θ) and assume that (a, b) is uniformly ensemble
821 controllable. Then for every T > 0, and every k ∈ N, ak b ̸∈ Ra,b (T ).

822 Proof. Let us set Θ e = a(Θ), we have Θ e ⊂ R is a compact set, θe0 = a(θ0 ) ∈ Θ e is an accumulation
823 point of Θ,
e and we set ı e the identity map on Θ.
Θ
e Since the pair (a, b) is ensemble controllable, we
824 necessarily have that a is injective (see [11, Proposition 4]) and b(θ) ̸= 0 for every θ ∈
Θ (see Re-
825 mark 2.5). This ensures that if y ∈ Ra,b (T ) then y ◦ a−1 ∈ RıΘe ,b◦a−1 (T ) = b ◦ a−1 RıΘe ,1 (T ).
826 Consequently, it is enough to prove this result for a = ıΘ and b = 1.
827 Let us assume by contradiction that there exists k ∈ N and u ∈ L2 (0, T ) such that
R T (T −t)θ
828
0
e u(t) dt = θk for every θ ∈ Θ. Since Θ admits an accumulation point, by analytic
829 continuation, we have
Z T
e(T −t)θ u(t) dt = θk (θ ∈ R).
0
830 For every ℓ ∈ N, we also have by derivation with respect to θ,
Z T
(T − t)k+1+ℓ e(T −t)θ u(t) dt = 0 (θ ∈ R).
0

831 Finally, using Stone-Weierstrass Theorem, we conclude that

(T − t)k+1 e(T −t)θ u(t) = 0 (θ ∈ R, t ∈ [0, T ] a.e.)

832 and using the fact that (T − t)k+1 e(T −t)θ does not vanish on (0, T ), we conclude that u = 0. This
RT
833 leads to a contradiction with 0 e(T −t)θ u(t) dt = θk .

35
834 Even if the above result is given in the real context, it can be easily extended to the complex
835 case (i.e., with a(θ), b(θ) ∈ C, complex controls and Θ ⊂ C). Its complex version will be used in
836 the next result.
837 We are now in position to give the general result.

838 Proposition B.2. Let Θ ⊂ R be a compact set and assume that int Θ ̸= ∅. Let A ∈ C(Θ)n×n and
839 B ∈ C(Θ)n , and assume that (A, B) is uniformly ensemble controllable. Then for every T > 0 and
840 every k ∈ N, we have
Ak B ̸∈ RA,B (T ).

841 Proof. Following [10, Lemma A.2], there exists θ0 ∈ int Θ, δ > 0, λ ∈ C([θ0 − δ, θ0 + δ], C) and
842 v ∈ C([θ0 − δ, θ0 + δ], C)n such that v(θ) ̸= 0 and A(θ)∗ v(θ) = λ(θ)v(θ) for every θ ∈ [θ0 − δ, θ0 + δ].
843 We now notice that if the system (2.1) is ensemble controllable on Θ, then the system given by the
844 pair (λ, v ∗ B) is also uniformly ensemble controllable on [θ0 −δ, θ0 +δ]. According to Proposition B.1,
845 we have that λk v ∗ B ̸∈ Rλ,v∗ B (T ) for every T > 0. This clearly ensures that Ak B ̸∈ RA,B (T ) for
846 every T > 0.

847 Remark B.3. Let us mention that in the case m > 1, the situation is not so trivial.
848 For instance, let T > 0, v ∈ C ∞ ([0, T ]) such that v (k) (0) = v (k) (T ) = 0 for every k ∈ N, and such
RT R −1
T
849 that 0 e(T −t)θ v(t) dt ̸= 0, and define β(θ) = 0 e(T −t)θ v(t) dt . For every k ∈ N, we have
T
θk ∈ Rθ,(1 β) (T ). In fact, we clearly have 1 = 0 e(T −t)θ β(θ)v(t) dt, from this relation and by
R
850
RT
851 integration by parts, one can deduce that θk = 0 e(T −t)θ β(θ)v (k) (t) dt for every k ∈ N.

852 References
853 [1] A. Agrachev, Y. Baryshnikov, and A. Sarychev. Ensemble controllability by Lie algebraic
854 methods. ESAIM, Control Optim. Calc. Var., 22(4):921–938, 2016.
855 [2] H. Bauer. Measure and integration theory. Transl. from the German by Robert B. Burckel,
856 volume 26 of De Gruyter Stud. Math. Berlin: de Gruyter, 2001.
857 [3] K. Beauchard, J.-M. Coron, and P. Rouchon. Controllability issues for continuous-spectrum
858 systems and ensemble controllability of Bloch equations. Commun. Math. Phys., 296(2):525–
859 557, 2010.
860 [4] A. Becker and T. Bretl. Approximate steering of a unicycle under bounded model perturbation
861 using ensemble control. IEEE Transactions on Robotics, 28(3):580–591, 2012.

862 [5] J. M. Borwein and Q. J. Zhu. Techniques of variational analysis, volume 20 of CMS Books
863 Math./Ouvrages Math. SMC. New York, NY: Springer, 2005.
864 [6] H. Brezis. Functional analysis, Sobolev spaces and partial differential equations. Universitext.
865 Springer, New York, 2011.
866 [7] E. J. Candès, J. Romberg, and T. Tao. Robust uncertainty principles: exact signal reconstruc-
867 tion from highly incomplete frequency information. IEEE Trans. Inform. Theory, 52(2):489–
868 509, 2006.
869 [8] X. Chen. Controllability issues of linear ensemble systems over multi-dimensional parameter-
870 ization spaces, 2020. preprint.

36
871 [9] B. Danhane and J. Lohéac. Ensemble controllability of parabolic type equations. Syst. Control
872 Lett., 183:10, 2024. Id/No 105683.
873 [10] B. Danhane, J. Lohéac, and M. Jungers. Conditions for uniform ensemble output control-
874 lability, and obstruction to uniform ensemble controllability. Math. Control Relat. Fields,
875 14(3):1128–1175, 2024.
876 [11] G. Dirr and M. Schönlein. Uniform and Lq -ensemble reachability of parameter-dependent
877 linear systems. J. Differ. Equations, 283:216–262, 2021.
878 [12] I. Ekeland and R. Témam. Convex analysis and variational problems., volume 28 of Classics
879 Appl. Math. Philadelphia, PA: Society for Industrial and Applied Mathematics, unabridged,
880 corrected republication of the 1976 English original edition, 1999.
881 [13] C. Esteve, B. Geshkovski, D. Pighin, and E. Zuazua. Large-time asymptotics in deep learning,
882 2021.

883 [14] G. B. Folland. Real analysis. Pure and Applied Mathematics (New York). John Wiley &
884 Sons, Inc., New York, second edition, 1999. Modern techniques and their applications, A
885 Wiley-Interscience Publication.
886 [15] P. A. Guth, K. Kunisch, and S. S. Rodrigues. Tracking optimal feedback control under
887 uncertain parameters. Phys. D, 467:Paper No. 134245, 13, 2024.

888 [16] U. Helmke and M. Schönlein. Uniform ensemble controllability for one-parameter families of
889 time-invariant linear systems. Syst. Control Lett., 71:69–77, 2014.
890 [17] E. Kreindler and P. Sarachik. Constrained linear-quadratic optimization problems with
891 parameter-dependent entries. Journal of Optimization Theory and Applicationsl, 198(2):781–
892 804, 2023.

893 [18] M. Lazar and J. Lohéac. Chapter 8 - control of parameter dependent systems. In E. Trélat
894 and E. Zuazua, editors, Numerical Control: Part A, volume 23 of Handbook of Numerical
895 Analysis, pages 265–306. Elsevier, 2022.
896 [19] M. Lazar and C. Molinari. Optimal distributed control of linear parabolic equations by spectral
897 decomposition. Optimal Control Appl. Methods, 42(4):891–926, 2021.
898 [20] J.-S. Li and N. Khaneja. Control of inhomogeneous quantum ensembles. Phys. Rev. A,
899 73:030302, Mar 2006.
900 [21] J.-S. Li and J. Qi. Ensemble control of time-invariant linear systems with linear parameter
901 variation. IEEE Trans. Autom. Control, 61(10):2808–2820, 2016.

902 [22] Y. Li and S. Osher. Coordinate descent optimization for ℓ1 minimization with application to
903 compressed sensing; a greedy algorithm. Inverse Probl. Imaging, 3(3):487–503, 2009.
904 [23] W. Miao, G. Cheng, and J.-S. Li. On uniform ensemble controllability of diagonalizable linear
905 ensemble systems. preprint, 2021.

906 [24] J. Peypouquet. Convex optimization in normed spaces. Theory, methods and examples.
907 SpringerBriefs Optim. Cham: Springer, 2015.
908 [25] M. Schönlein. Computation of open-loop inputs for uniformly ensemble controllable systems.
909 Math. Control Relat. Fields, 12(3):813–829, 2022.

37
 910 [26] M. Schönlein. Polynomial methods to construct inputs for uniformly ensemble reachable linear
911 systems. Math. Control Signals Syst., 36(2):251–296, 2024.
912 [27] M. Schönlein and U. Helmke. Controllability of ensembles of linear dynamical systems. Math.
913 Comput. Simul., 125:3–14, 2016.
914 [28] R. Triggiani. Controllability and observability in Banach space with bounded operators. SIAM
915 J. Control 13, 462–491 (1975)., 1974.
916 [29] E. Zuazua. Controllability and observability of partial differential equations: some results
917 and open problems. In Handbook of differential equations: evolutionary equations. Vol. III,
918 Handb. Differ. Equ., pages 527–621. Elsevier/North-Holland, Amsterdam, 2007.

38

View publication stats