Polynomial methods to construct inputs for uniformly ensemble reachable linear systems

This paper is concerned with linear parameter-dependent systems and considers the notion of uniform ensemble reachability. The focus of this work is on constructive methods to compute suitable parameter-independent open-loop inputs for such systems. In contrast to necessary and suﬃcient conditions for ensemble reachability, computational methods have to distinguish between continuous-time and discrete-time systems. Based on recently derived suﬃcient conditions and techniques from polynomial approximation, we present two methods for discrete-time single-input linear systems. Moreover, we illustrate that one method can also be applied to certain continuous-time single-input systems.


Introduction
Ensemble control is an emerging field in mathematical systems and control theory referring to the task of controlling a large, potentially infinite, number of states, or systems, using a single input function or a single feedback controller. Being not a rigorously defined term it subsumes several different scenarios, where in each particular case the problems and the required techniques differ from each other. Ensemble control embraces the scenario of uncertainties in the initial data. This is modeled by a probability distribution of the initial state and the ensemble control problem leads to a transport problem of density functions and therefore to controllability and observability issues of the Liouville and the Fokker-Plank equation [5], [8], [17], [52], [54]. Another setting tackles uncertainties in the model parameters by addressing controllability problems for parameter-dependent systems, and the goal is to achieve a certain control task by using only a single or a few openloop inputs which are independent of the (usually unknown) model parameters. The notion ensemble controllability is used for this, cf. [33]. In this context we also mention the notions averaged controllability and moment controllability, cf. [55] and [51], respectively. Motivated by the controllability analysis of the Bloch equation, there is also the concept of asymptotic controllability, cf. e.g. [4], [9].
In this paper we consider families of parameter-dependent systems of the form ∂ ∂t x(t, θ) = A(θ)x(t, θ) + B(θ)u(t) where θ ∈ P is considered as a parameter and the parameter space P ⊂ C is assumed to be a Jordan arc in the complex plane, i.e. P is the image of a continuous and bijective function defined on a compact interval. Moreover, we assume that the matrix-valued functions A : P → C n×n and B : P → C n×m are continuous. Throughout the paper we will use the short notation (A, B) ∈ C n,n (P)×C n,m (P).
In addition, let C n (P) denote the space of continuous functions from P to C n . As we mainly consider discrete-time systems, we take the initial condition x(0, θ) = x 0 (θ) = 0. To express that the solutions to (1) are regarded as functions from the parameter space P to C n , we denote them by ϕ(T, u)(θ) := ϕ(T, u, 0, θ). Further, we wish to highlight the essential property that the input is independent of the parameters.
The notion of reachability we are considering in this paper is as follows: We say that a pair (A, B) ∈ C n,n (P)×C n,m (P) is uniformly ensemble reachable (from zero), if for every f ∈ C n (P) and ε > 0 there are T > 0 and u ∈ L 1 ([0, T ], C m ) or u = (u 0 , u 1 , . . . , u T −1 ), u i ∈ C m such that Note that ensemble reachability is an infinite-dimensional problem (unless P is finite) and it is equivalent to approximate reachability of the infinite dimensional linear system defined by the bounded linear matrix-multiplication operators M A : C n (P) → C n (P), M A f (θ) = A(θ)f (θ) and M B : C m → C n (P), M B u = B(θ)u.
As a consequence of the restriction that the inputs are not allowed to depend on the parameter, it is well-known that exact reachability (i.e. ε = 0) is never possible provided P is not finite, cf. [ for every T > 0, for every ε > 0 and for every pair x 0 , x 1 ∈ C n (P) there exists an input u ∈ L 1 ([0, T ], C m ) such that We note that this equivalence does not hold for discrete-time systems. Recall that a pair (A, B) ∈ C n,n (P) × C n,m (P) is called uniformly ensemble controllable (to zero), if for all x 0 ∈ C n (P) and ε > 0 there exist T ≥ 0 and an input u such that ϕ(T, u, x 0 ) ∞ < ε.
It is well-known that the notions of approximate reachability (from zero) and approximate controllability (to zero) are independent of each other, cf. [18,Lemma 4.1], and none of both implies approximate complete controllability. Also we note that this paper does not cover other notions of controllability, such as asymptotic controllabilty, cf. [45]. Moreover, recall that the application of [ where b j (θ) denotes the jth column of B(θ). In operator theory language, uniform ensemble reachability of (A, B) is equivalent to the fact that the matrix multiplication operator M A is m-multicyclic and the columns b 1 , ..., b m are cyclic vectors. For other characterizations, we refer to [12,Section 6.2] and [16, Section VI.8.a]. As these characterizations have the drawback that they are hard to check in practice, much effort has recently been spent on the derivation of pointwise testable necessary and sufficient conditions, cf. [14], [34].
The ensemble controllability problem for continuous-time parameter-dependent systems is also studied in [1], [2], [6], [32], [33], [34], [35] and [51]. Agrachev and Sarychev consider ensemble controllability for nonlinear drift-free parameterdependent systems and provide a characterization in terms of Lie-brackets. In the same direction, the work of Chen [6] also treats nonlinear systems and considers Lie extensions. We note that these approaches do not apply to the setting in this paper. In [32] a characterization for ensemble controllability for time-varying parameter-dependent linear systems is presented, which is based on the singular value decomposition of the reachability operator. Like condition (2), this condition is hard to check. The recent works [7] and [13] have shown that, if (A, B) satisfy appropriate regularity assumptions, L p -and uniform ensemble controllability, respectively, cannot hold if the parameter space P ⊂ R d , d ≥ 2 has interior points.
Problem statement The objective of this paper is as follows. We consider uniformly ensemble reachable pairs (A, B) ∈ C n,n (P) × C n,m (P). Given a desired target family f ∈ C n (P) and an ε-neighborhood of it, we treat the problem of how to find a suitable T > 0 and a suitable input u such that For continuous-time systems this problem has recently been addressed in [31], [43] and [53]. The methods in [43] and [53] are based on techniques from integral equations. Also we mention that the results in [43] require an additional assumption on the target family f so that predefined error bounds can be archived. Besides that, the work [31] uses the adjoint system to set up a continuous-time optimal control problem to derive suitable inputs. However, as the solution formula for discrete-time systems is not given in terms of an integral, the methods just mentioned cannot be applied to discrete-time systems. Moreover, we are not aware of any approaches in the literature to tackle the discrete-time case.
Novelty and main contribution The scope of the present paper is to fill this gap. More precisely we focus on single-input systems and start with treating the discrete-time case The solution at time T > 0 for inputs u = (u 0 , ..., u T −1 ) ∈ C × · · · × C is given by Noticing the restriction that the inputs u 0 , ..., u T −1 have to be independent of the parameters, the solution can be expressed in terms of the polynomial That is, we have ϕ(T, u)(θ) = p(A(θ))b(θ).
In this notation, the problem can be stated equivalently as follows. Given f ∈ C n (P) and ε > 0, how to find a suitable polynomial p such that This observation is the key to a new strategy to solve the construction problem by using tools from polynomial approximation. In this context, given the target family f ∈ C n (P) and the ε-neighborhood of it, the task of finding the amount required inputs and their values is equivalent to derive the degree and the coefficients of an approximating polynomial (in the sense of (3)). Note that this fact is significantly different from the continuous-time case, where the solution formula is given in terms of an integral operator. The contributions of this paper are made up of two parts. On the one hand, the paper provides a collection of results from (complex) polynomial approximation, for which explicit representations of the corresponding polynomials are available. More precisely, after a short recap of explicit refinements of the classical Weierstrass Approximation Theorems, the paper provides an explicit construction of the polynomial in Runge's Little Theorem, where the degree and the coefficients are determined in terms of the target function and the approximation error. To the best of the author's knowledge, such an explicit representation was not given in the literature so far. Moreover, the paper gives an explicit treatment of results which are due to Walsh. As these results are required later for the derivation of suitable inputs, the paper also provides a detailed treatment of these results (which is only implicitly contained in the literature). More specifically, we show that the construction of appropriate polynomials can be traced back to Runge's Little Theorem by using suitable conformal mappings. We note that the numerical treatment of conformal mappings is a wide area. Hence, this topic is beyond the scope of this paper and we only provide some comments to the relevant literature.
On the other hand, based on the exposition on polynomial approximation, the paper shows various methods to determine appropriate inputs. As there is no explicit characterization of uniform ensemble reachability 1 we take two sufficiency conditions, denoted by (S1) and (S2) in Section 3, derived in [14] as a starting point. In this paper we present constructive methods for each condition. Moreover, the paper shows that the methods corresponding to the sufficiency condition (S2) can also be applied to continuous-time single-input systems, while the methods for (S1) are limited to discrete-time systems.
Organization Let's quickly describe the structure of the paper. Section 2 contains the relevant results from polynomial approximation that are required for the input construction procedures. It is divided into three subsections corresponding to the results of Runge, Walsh and Weierstrass, respectively. Section 3 is devoted to the presentation of the constructive procedures to determine suitable inputs for discrete-time parameter-dependent linear single-input systems. After a short summary of relevant known results, the section has two subsections, where each subsection presents methods corresponding to one of the sufficiency conditions. Section 4 depicts that methods corresponding to sufficiency condition (S2) can also applied to continuous-time systems. Finally, Section 5 contains additional perspectives about the presented methods.
Moreover, we call C n (Ω) the set of continuous functions g : Ω → C n and we say that g ∈ C k n (Ω) if g is k-times continuously differentiable. A function g : Ω → C n is said to satisfy a Lipschitz condition if there is a L g > 0 such that for all z 1 , z 2 ∈ Ω. The set of functions satisfying a Lipschitz condition is denoted by Lip n (Ω).
If Ω ⊂ C is compact and g : Ω → C is continuous, we define For a function f : [a, b] → R, the nth Bernstein polynomial is given by For a continuous function f : ∂D → C, its kth Fourier coefficient is denoted bŷ Further, we call the nth Fejér polynomial.

Elements of Approximation Theory
In this section we discuss results from approximation theory, for which constructive proofs are available. The presented results are due to Bernstein, Runge, Weierstrass and Walsh. Note that the definite result in this context is the famous Mergelyan's Theorem, saying that a continuous function f : K → C can be uniformly approximated by polynomials if K is compact, C \ K is connected and f is analytic in the interior of K, cf. [24, Chap. III, § 2, Theorem 1]. However, its highly ingenious method of proof is not constructive. Therefore we will present some special cases which can be proved constructively and prepare the ground for the computational methods in Section 3. In Section 2.1 we recall known explicit versions of the classical Weierstrass approximation theorems. Moreover, in Section 2.2 we present an explicit treatment of the approximating polynomial in Runge's Little Theorem. In Section 2.3 we present a detailed derivation of some of Walsh's results on complex approximation. Note that these results appeared more then twenty years before Mergelyan's result.

Weierstrass Theorems
We start with the Weierstrass approximation theorems and consider the special cases where the continuous function that shall be approximated additionally satisfies a Lipschitz condition. Recall that for a function f ∈ C([a, b], R) the nth Bernstein polynomial is given by For complex-valued functions f : [a, b] → C one can apply the latter to the real and imaginary part. That is, for f (x) = g(x) + ih(x) we consider the complex Bernstein polynomial B n,f (x) = B n,g (x) + iB n,h (x). Next, we recap a result of Gzyl and Palacios that provides an explicit error bound, cf. [27].
Theorem 1 (Gzyl, Palacios (1997)) Let f : [a, b] → C satisfy a Lipschitz condition. Then, for n ≥ 3, the sequence of complex Bernstein polynomials satisfies Note that, since Gyzl and Palacios consider real-valued functions defined on the unit interval, the constants in the previous statement are adjusted. Similarly, the second or trigonometric Weierstrass approximation theorem can also be proven constructively. It is based on Fejér's Theorem and we first recall their definition. Let f : ∂D → C be a continuous function, then the the nth Fejér polynomial is given by where L f > 0 denotes the Lipschitz constant.
We note that sharper but less explicit error bounds for the second Weierstrass approximation Theorem have been derived in L. Lorch [36], S.M. Nikolski [39] and S.A. Telyakovskii [46].

Runge's Little Theorem
Another famous result from approximation theory is due to Runge. Before presenting a constructive proof, we have to fix some notation. The presentation of Runge's little Theorem is based on [26] and [42]. Let γ denote a closed (piecewise) C 1 -path. With a slight abuse of notation, we denote its trace also by γ. Then, denotes the winding number of z with respect to γ. Moreover, a closed polygon τ = [p 1 p 2 · · · p k p 1 ] composed of finitely many horizontal or vertical segments is called a grid polygon if there exists a not necessarily regular grid G ⊂ C of horizontal or vertical lines such that all vertices p 1 , ..., p k are pairwise distinct adjacent grid points of G. Moreover, let depending on the orientation of γ. We will not provide a complete proof. Instead, as we are interested in constructive methods, we will provide an explicit treatment of the degree of the approximating polynomial as well as an explicit representation of its coefficients. To best of the author's knowledge these are not available in the literature. For some analytical parts that will be left out we refer to [42].
Proof The proof is carried out in four constructive steps. We start with the construction of finitely many grid polygons in Ω \ K so that f can be represented via the Cauchy integral formula for compact sets.
Step 1: Grid polygon construction There are finitely many distinct oriented horizontal or vertical line segments τ 1 , ..., τ N of length δ < 1 Proof of Step 1: Let δ ∈ (0, 1 √ 2 d(K, ∂Ω)) and consider a grid consisting of lines that are parallel to the real and imaginary axis so that the distance between two parallel lines is δ. Since K is compact, there are finitely many boxes Q 1 , ..., Q n that intersect with K and satisfy The boundaries of the boxes Q 1 , ..., Q n consist of line segments. We now pick those boundary segments that are not common for two distinct boxes Q l and Q m with l = m. Let τ 1 , ..., τ N denote the selected segments. The construction yields that Indeed, if line segment τ k meets K, there are two boxes of the grid having τ k as a common side, a contradiction to the choice of the segments τ 1 , ..., τ N .
Let z ∈ K be a point that does not lie on any boundary of the boxes Q 1 , ..., Q n . Then there is exactly one box Q l containing z and it holds Since common sides of different boxes of the grid occur twice with different orientation, it follows that Moreover, as shown in [42, Chapter 12, § 1.1], the last equality also holds for all z ∈ K that are on the boundary of some box. Now, let ε > 0, δ > 0 and N be fixed. Moreover, let L > 0 denote the Lipschitz constant of the function ξ → f (ξ) ξ−z on τ .

Proof of
Step 3: Since The assertion then follows from observing that (7) is equivalent to Step 4: Polynomial approximation The polynomial where approximates r b uniformly on K, i.e., p satisfies Proof of Step 4: By construction, for z ∈ K, we have Since the function 1 (z−b) v+1 is holomorphic on B r (0), its series expansion around zero is given by the coefficients That is, and one has By the Cauchy estimates [26, Theorem 3.4.1], one has where in the last inequality we used that η r m (v) kl +1 < 1 since η < r. By construction, it holds δ kl > r and thus (10) is equivalent to This completes step 4.
The above construction steps yield for every z ∈ K which completes the proof.
Remark 1 (a) The set N k=1 τ k is determined only by the sets Ω and K, and the grid constructed in Step 1. The set N k=1 τ k also contains all the poles of r and is independent from the quality of approximation. (b) As ε > 0 is decreasing the number of poles of the approximating rational function is increasing. The location of the poles is not affected, i.e., the poles will still be located on N k=1 τ k . (c) It is shown in [42,Chapter 12,§ 4] that the line segments τ 1 , ..., τ N define a cycle τ consisting of finitely many grid polygons. In general, it is not possible to find a single closed path Γ such that A counterexample is given [ Theorem 4 ) Let γ be a Jordan arc. If f is continuous on γ, it can be uniformly approximated on γ by a polynomial in z.
A proof of Theorem 4 is given at the end of this section. We emphasize that the following exposition is more explicit in terms of constructing the polynomial than the ones in the literature. The significance of the subsequent exposition is that the construction of the polynomial p is traced back to Runge's result by using suitable conformal mappings. Numerical procedures for conformal mappings are an active research area and beyond the scope of this paper. At the end of this section we provide more detailed comments and references.
Before we come to the proof, we need some preparation. To do so, we have to consider sequences of Jordan curves, or more precisely we provide a notion of convergence for sequences of Jordan curves. Recall that, a Jordan curve is a homeomorphic image (within C) of the unit circle ∂D.
Following Courant [11], a sequence (Γ n ) n∈N ⊂ C of Jordan curves is said to converge to a Jordan curve Γ if (i) any accumulation point of any sequence (z n ) n∈N with z n ∈ int Γ n lies on Γ and (ii) for every ε > 0 there are n(ε) ∈ N and δ(ε) > 0 with lim ε→0 δ(ε) = 0 such that for every z ∈ Γ and for every Let Ω be a proper subdomain of C. Recall that a function f : Ω → D is called conformal, if it is holomorphic, one-to-one and onto. The following result is due to Courant, cf. [10,11].
Using a different method, Courant's Theorem was also derived by Radó [41]. We note that Radó used the Fréchet-distance to characterize convergence for sequences of Jordan curves. That these notions are in fact equivalent was shown by Markushevich. A nice and in-depth exposition for this, as well as extensions to more general domains can be found in [21].
For a given Jordan curve Γ , a sequence of Jordan curves (Γ n ) n∈N that converges to Γ from outside, i.e., Γ n ⊂ C \ int Γ , can be constructed via the grid construction in Step 1 in the proof Runge's Theorem. Indeed, as illustrated in Step 1 in the proof Runge's Theorem, it is sufficient to choose a grid consisting of horizontal and vertical lines defined by x = k 2 n y = l 2 n , k, l ∈ Z and define Γ n as the boundary of the boxes given by the vertices l + 1 2 n that intersect with int Γ . Next we recap a proof for a refinement of Runge's theorem which is also due to Walsh, cf. [49,Satz].
Theorem 6 (Walsh (1926)) Let Γ be a Jordan curve and suppose that f is holomorphic on int Γ and continuous on int Γ . Then, f can be uniformly approximated on int Γ by polynomials. That is, for every ε > 0 there is a polynomial p such that Proof Let ε > 0 and let w.l.o.g. z = 0 ∈ int Γ . Using the grid construction as in the proof of Runge's theorem there is a sequence (Γ n ) n∈N of grid polygons such that Γ n ⊂ C \ int Γ and (Γ n ) n∈N converges to Γ . Then we consider the conformal mappings Step 1: (Approximation of f by a sequence of holomorphic functions on int Γ ). We consider the holomorphic functions and shall show that for every ε > 0 there is an N ∈ N such that Indeed, for n ∈ N we consider the Jordan curve γ n := Φ n (Γ ) ⊂ D and note that for every z ∈ int Γ there exists an w ∈ int γ n so that z = Φ −1 n (w). Thus, together with the maximum modulus theorem [26,Corollary 5.3.4], it is sufficient to show that there is a K = K(ε) ∈ N such that Note that, int Γ is compact and f is continuous, and so f is uniformly continuous on int Γ . Thus there is a δ(ε) > 0 such that |f (z) − f (w)| < ε 2 for all z, w ∈ int Γ with |z − w| < δ(ε). Since the Jordan curves (Γ n ) n∈N converge to Γ we can apply Courant's Theorem 5 to the sequence (Φ −1 n ) n∈N and conclude that there is an N = N (ε) ∈ N such that which shows (11).
Step 2: (Approximation of g N via Runge's Theorem). The function is holomorphic on the compact set int Γ . Then, by Runge's Theorem 3, for ε 2 > 0 there is a polynomial p such that Consequently, for every z ∈ int Γ it follows that We will use the latter statement to prove Walsh's extension of the second Weierstrass Theorem to arbitrary Jordan curves, cf. [50, Satz I]. n − |k| nĝ where theĝ(k) denote the Fourier coefficients of g, cf. (5). Since for every z ∈ Γ there is a unique w ∈ ∂D such that z = Φ −1 (w), by Theorem 2 (a), there is a To show the claim, we shall prove that on Γ , the function Φ can be uniformly approximated on Γ by polynomials in z and the function 1 Φ can be uniformly approximated by polynomials in z and 1 z . Let δ > 0; we will specify it later, cf. (13). By Theorem 6, there is a polynomial p 1 such that Next, we show that there is a polynomial p 2 such that where a −1 is the residue of 1 Φ at zero, i.e.
To see this, recall that Φ is conformal on int Γ with Φ(0) = 0. Hence, the function 1 Φ is meromorphic on int Γ with a simple pole at 0. Moreover, the function z has a removable singularity at zero. Let Ψ : int Γ → C denote its holomorphic extension, i.e.
which is continuous on int Γ . Then, by Theorem 6, there is a polynomial p 2 such that Next, we shall show that, if we choose δ > 0 such that one has To see this, note that |Φ(z)| = 1 for all z ∈ Γ and, thus, Moreover, we recall that for all v, w ∈ C it holds |v n − w n | ≤ |v − w| n−1 k=0 |v| k |w| n−1−k .
The combination of the latter yields that for all z ∈ Γ one has Similarly, for all z ∈ Γ we get Finally, defining the polynomial we derive the following estimate This shows the assertion.
Proof of Theorem 4 Without loss of generality we assume that γ does not contain the origin 2 . Then, the Jordan arc can be extended to a Jordan curve Γ so that 0 ∈ int Γ . Also the function f can be continuously extended to Γ . This can be achieved by defining where w 1 and w 2 are the values of f at the end-points z 1 and z 2 of γ, respectively. By Caratheodory's Theorem [26, Theorem 13.2.3] there is a continuous and bijective map Φ : int Γ → D so that Φ(0) = 0 which is conformal on int Γ . Let ε > 0. By Theorem 7, for the continuous function f : Γ → C there are N = N (ε) ∈ N and a −1 ∈ C as well as polynomials p 1 , p 2 , such that the polynomial Here, the coefficients c −N +1 , ..., c N −1 are given by the Fourier coefficients corresponding to the function f • Φ −1 : ∂D → C, cf. (5), i.e.
Furthermore, on γ the function z → 1 z can be uniformly approximated by polynomials in z. To see this, let γ 0,∞ be a Jordan arc connecting the origin and ∞ such that γ 0,∞ ∩ γ = ∅. Then, Ω := C \ γ 0,∞ is simply connected, ∞ ∈ int Ω and the function z → 1 z is holomorphic on int Ω. Let ν := max Moreover, choose δ > 0 such that By Runge's Theorem 3, there is a polynomial p 3 in z such that Next, we define the polynomial Similar to the proof of Theorem 7, we get Using the triangle inequality, for every z ∈ γ, one has and |p 3 (z) + p 2 (z)| < δ + |a −1 |ν + µ as well as Thus, by (16) we have Consequently, it holds This shows the assertion.
We close this section with some comments on the construction of conformal mappings and their numerical computation. There are different approaches available in the literature. Peter Henrici presents a method to determine conformal mappings by solving a Dirichlet boundary problem. The solution of it together with its harmonic conjugate determines the conformal map, cf. [29,Theorem 16.5a]. Moreover, Henrici describes numerical techniques to solve the Dirichlet problem by means of numerical methods to solve Symm's integral equation of first kind, cf. [29, § 16.6 V]. Besides that, Tobin. A. Driscoll and Llyod N. Trefethen depict procedures, as suggested by Peter Henrici, to compute conformal maps by using the Schwarz-Christoffel mapping, cf. [15]. Another, recent approach to compute conformal maps is the Zipper-algorithm [37]. For earlier approaches to this area we refer to [22,23]. Furthermore, approximation techniques for conformal mappings based on Bergman kernels are outlined in [24] and [26]. Finally, we mention the monograph [40] which contains a profound literature review and many more references.

Discrete-time single-input systems
In this section we consider single-input pairs (A, b) ∈ C n,n (P) × C n (P), where the parameter space P is a Jordan arc. Especially, for given f ∈ C n (P) and ε > 0, we tackle the central problem of how to determine suitable inputs such that the origin is steered into the ε-neighborhood of f . The basic tools for the constructions will be the results of Bernstein, Runge, Weierstrass and Walsh from the previous section. Indeed, as pointed out in the introduction, the solution to Hence, given a target function f ∈ C n (P) and ε > 0, we shall explore the methods of Section 2 to derive a suitable polynomial p so that

Starting point: Known sufficient conditions
As a starting point, we take the necessary and sufficient conditions for uniform ensemble reachability developed in [14,Theorem 4,Corollary 3]. Since these are essential, we provide a short recap. That is, if the pair (A, b) is uniformly ensemble reachable, the following necessary conditions are satisfied: (N1) The pair (A(θ), b(θ)) is reachable for every θ ∈ P. (N2) For any pair of distinct parameters θ, θ ′ ∈ P, the spectra of A(θ) and A(θ ′ ) are disjoint: Further, (A, b) is uniformly ensemble reachable if it satisfies (N1), (N2) and one of the following sufficiency conditions: (S1) For each θ ∈ P the characteristic polynomial of A(θ) takes the form z n − (a n−1 z n−1 + · · · + a 1 z + a 0 (θ)) for some a n−1 , ..., a 1 ∈ C and a 0 ∈ C(P). (S2) A(θ) has simple eigenvalues for each θ ∈ P.
A few comments are in order. First, we note that the conditions are valid for continuous-time systems and discrete-time systems. Also, we note that a pointwise characterization for uniform ensemble reachability is currently not available. In [14,Theorem 3 (d)] it is shown that for single input pairs (A, b) a necessary condition for uniform ensemble reachability is that the eigenvalues are simple on an open dense subset of the parameter space P. In this sense condition (S2) is not far from being necessary. Furthermore, in certain cases, like the controlled harmonic oscillator, cf. [28], it happens that the sufficiency conditions (S1) and (S2) are satisfied at the same time. Besides, to the best of the author's knowledge, other general sufficient conditions are not known. Moreover, we note that, due to the condition (N2) the function a 0 in (S2) is necessarily injective. Thus, a 0 : P → a 0 (P) is one-to-one and onto and so a 0 (P) defines a Jordan arc. It will turn out that in the special case, where a 0 (P) is a subset of the real axis, based on Theorem 1, a constructive method can be obtained by using Bernstein polynomials. The general case is covered by Walsh's Theorem 4. This method turns out to be more involved and will require also the Theorems 2, 3, 6 and 7.
Before we begin with the constructions, we note that we start the constructions based on the conditions (N1), (N2) in (S1) or (S2) and discuss additional required conditions along the lines. The corresponding results are formulated subsequently to each construction. This has the advantage that the results and the corresponding polynomials can be spotted simultaneously.

Methods for condition (S1)
In this subsection, we will explore the circumstances such that the verification of the sufficient condition (S1) (cf. proof of Theorem 4 in [14]) gives rise to methods for the construction of suitable T > 0 and inputs u = (u 0 , ..., u T −1 ). Let's we assume that the pair (A, b) satisfies the conditions (N1), (N2) and (S1). Also, let f ∈ C n (P) and ε > 0 be given. By condition (N1), the matrix is continuously invertible for every θ ∈ P and one has where a 0 (θ), a 1 , ..., a n−1 denote the coefficients of the characteristic polynomial of A(θ). Writing the solution formula in terms of a polynomial p, we get where · M is a matrix norm that is submultiplicative to · and As shown in [14], for we get and, hence, Moreover, by (N2) and (S1), the function θ → a 0 (θ) is one-toone. Consequently, we get if we can find polynomials p 1 , ....p n satisfying for all k = 1, ..., n.
Depending on the structure of a 0 (P), we will present two constructive methods.

Case 1: a 0 (P) = [a, b] ⊂ R defines a compact interval
We assume that for each k = 1, ..., n the function z → (T −1 f ) k • a −1 0 (z) is Lipschitz continuous 3 and let L k >0 denote the Lipschitz constant. By Theorem 1, the condition (19) is fulfilled if we take p k as the Bernstein polynomial, cf. (4), corresponding to the function z → ( Denoting this Bernstein polynomial by B (T −1 f ) k •a −1 0 ,n k , we conclude that for Considering the monomial representation we get the following result.
Next we consider the more general problem of constructing polynomials p k satisfying (19) if a 0 (P) defines a Jordan arc. It turns out that this is more involved. We shall see how the proof of Theorem 4 gives rise to a constructive method.
Case 2: a 0 (P) defines a Jordan arc Without loss of generality 4 we assume that a 0 (P) does not contain the origin. In a first step we extent a 0 (P) to a Jordan curve Γ 0 such that the origin lies inside the Jordan curve. Then, take a continuous mapping Φ 0 : int Γ 0 → D which is conformal on int Γ 0 so that Φ 0 (0) = 0. The existence of such a mapping is ensured by Caratheodory's Theorem, cf. [26,Thm. 13.2.3].
Step 1. Extension of the target function to the Jordan curve Here, we extend the functions z → (T −1 f ) k • a −1 0 (z) to Γ 0 , where T ∈ C n,n (P) denotes the changes of coordinates introduced at the beginning of Section 3.1. This can be done as follows. For simplicity, we do not introduce a new notation. For all z ∈ Γ 0 \ a 0 (P) we define where w k,1 and w k,2 are the values of (T −1 f ) k • a −1 0 at the end-points z k,1 and z k,2 of a 0 (P), respectively. As a result, we obtain a continuous function In the following, we assume that a 0 ∈ C 2 (P) and P has a parametrization which is twice continuously differentiable. Also we suppose that the extension Γ 0 has the same properties. Then, by Kellogg's Theorem [25,Theorem 4.3], the functions h k satisfy a Lipschitz condition provided (T −1 f ) k is Lipschitz continuous 5 . We assume that this is the case and let L k > 0 denote the corresponding Lipschitz constant.
Step 2. Ansatz for the polynomial p To derive a suitable polynomial p, we take again (cf. (18)) p k z n − a n−1 z n−1 − · · · − a 1 z z k−1 .
Based on the exploration in Section 2.3 (in particular (17)), we will construct the polynomials p 1 , ..., p n as follows. For each k = 1, ..., n we consider the following Fejér-like polynomial where c l := n k − |l| n k h k (l), l = −n k + 1, ...., n k − 1 denote the Fourier coefficients of h k , cf. (5). We proceed with a simple but instructive step. By adding and subtracting appropriate terms and using the triangle inequality, one has denotes the residue of 1 Φ 0 at zero. Step 3. Derivation of n k ∈ N Applying Weierstrass Second Theorem 2 (b) to the function As Φ 0 is one-to-one and onto and taking n k ∈ N so that Hence it remains to determine the polynomials p k,1 , p k,2 , p k,3 . This will be carried out in the following three steps. To this end, we pick δ k > 0 so that Step 4. Derivation of p k,1 According to Cournat's Theorem 5, we take a Jordan curve Γ 1 such that Γ 1 ⊂ C \ int Γ 0 and we choose a continuous mapping The construction of Γ 1 and hence of Φ 1 depends on the properties of a 0 . We note that numerical methods for conformal mappings is a research area itself and, thus, it is beyond the scope of this paper. We refer to the comments provided at the end of Section 2.
To derive the polynomial p k,1 we apply Runge's Little Theorem 3 to the holomorphic function Φ 1 : int Γ 1 → D so that Then, for all z ∈ int Γ 0 , one has Further, as shown in the proof of Theorem 6 (cf. (14)), for all z ∈ Γ 0 it follows The polynomial p k,3 is then obtained by applying Runge's Little Theorem such that for all z ∈ a 0 (P) one has Then, as shown in the proof of Theorem 4, we get Thus, we have Step 7. Monomial representation of p Finally, we determine the monomial representation of the polynomial p(z) = n k=1 p k (z n + a n−1 z n−1 + · · · + a 1 z)z k−1 =: This construction yields the following Theorem 9 Let P be a Jordan arc having a twice continuously differentiable parametrization. Assume that (A, b) ∈ C 1 n,n (P) × C 1 n (P) satisfies (N1), (N2) and (S1). Let ε > 0, f ∈ C n (P) and a 0 ∈ C 2 (P) be a Jordan arc in the complex plane such that ( In particular, the result is true if (A, b) ∈ C 1 n,n (P) × C 1 n (P) and f ∈ Lip n (P).

Methods for condition (S2)
In this section we consider methods to construct suitable inputs based on the necessary conditions (N1), (N2) and the sufficiency condition (S2). Let f ∈ C n (P) and ε > 0 be given. After applying a change of coordinates T (θ), we consider the pair where λ 1 , ..., λ n denote the distinct eigenvalue Jordan arcs. Writing the solution formula in terms of a polynomial p, we get where (T −1 f ) k denotes the kth component of the vector-valued function θ → T −1 (θ)f (θ) ∈ C n . Thus, we conclude that if we can find a polynomial p so that for all k = 1, ..., n it holds Compared to the methods for (S1), we are a looking for a single polynomial working for every component of the target function. With other words, the polynomial p solves n approximation problems simultaneously, where each problem is defined on an individual domain. Inspired by [3], using the ansatz we will split the problem into n subproblems. Postponing the technicalities for a moment, on the one hand we will determine p 1 , ..., p n so that for each k ∈ {1, ..., n} the polynomial p k approximates the function (T −1 f ) k • λ −1 k over the image of λ k : P → C. Let λ k (P) denote the image of λ k . On the other hand, the polynomials q 1 , ..., q n are derived by approximating the indicator functions where λ(P) = n k=1 λ k (P).
Note that the assumptions imply that the sets λ 1 (P), ..., λ n (P) are pairwise disjoint and so h 1 , ..., h n are holomorphic functions. A precise outline for the construction of the polynomials p 1 , ..., p n and q 1 , ..., q n is presented in the following. In doing so, we will distinguish two cases. Exemplary we will discuss the construction process for polynomials p k and q k , for some k ∈ {1, ..., n}.
Case 2: λ 1 (P), ..., λ n (P) define Jordan arcs The construction process is similar to the methods for (S1). Hence, in order to avoid a lengthy presentation we will go through the analogous construction steps more quickly.
Step 1. Extension of the Jordan arcs and the target function Without loss of generality we will assume that none of the Jordan arcs contains the origin. In a first step we extent the Jordan arcs λ 1 (P), ..., λ n (P) to Jordan curves Γ 1 , ..., Γ n such that the origin lies inside of each Jordan curve. Then we take a continuous bijective mappings Φ k : int Γ k → D which are conformal on int Γ k so that Φ k (0) = 0. The existence of such mappings is ensured by Caratheodory's Theorem, cf. [26,Thm. 13.2.3].
Analogous to Step 1 in the methods for (S1), we extend the functions z → where T ∈ C n,n (P) denotes the changes of coordinates introduced at the beginning of Section 3.2. As a result, we obtain continuous functions If λ 1 , ..., λ n ∈ C 2 (P) and P has a parametrization which is twice continuously differentiable and the extensions Γ 1 , ..., Γ n have the same properties, then by Kellogg's Theorem [25,Theorem 4.3], the functions h k satisfy a Lipschitz condition provided the (T −1 f ) k are Lipschitz continuous 8 . Let L k > 0 denote the corresponding Lipschitz constant.
Step 2. Ansatz for the polynomial p As in the previous case, to derive a suitable polynomial, we set Based on the exploration in Section 2.3 (in particular (17)), for each k = 1, ..., n we make the ansatz where c l := n k − |l| n k h k (l), l = −n k + 1, ...., n k − 1 denote the Fourier coefficients corresponding to the function h k . Again we take a C \ int Γ k and we choose a continuous mapping Φ k,1 : int Γ k,1 → D which is conformal on int Γ k,1 satisfying Φ k,1 (0) = 0 and The construction of Γ k,1 and hence of Φ k,1 depends on the properties of λ k . Recall that constructive methods for conformal mappings are beyond the scope of this paper. We refer to the comments provided at the end of Section 2.
To derive the polynomial p k,1 , we apply Runge's Little Theorem 3 to the conformal mapping Φ k,1 so that Then, for all z ∈ int Γ k one has Further, as shown in the proof of Theorem 6 (cf. (14)), for all z ∈ Γ k it follows Step 5. Derivation of p k, 2 We consider the holomorphic function Ψ k : int Γ k → C, which is continuous on int Γ k . To determine p k,2 , we take another Jordan curve Γ k,2 such that Γ k,2 ⊂ C\int Γ k and we choose a continuous mapping Φ k,2 : int Γ k,2 → D which is conformal on int Γ k,2 satisfying Φ k,2 (0) = 0 and where L Ψ k > 0 denotes the Lipschitz constant of Ψ k . Then, for all z = Φ −1 k,2 (w) ∈ int Γ k we have Next, we apply Runge's Little Theorem to the holomorphic function and determine the polynomial p k,2 such that Based on the conditions (31) and (32), for every k = 1, ..., n we get Step 8. Monomial representation of p Finally, using the polynomials p 1 , ..., p n and q 1 , ..., q n , we determine the monomial representation of the polynomial This construction yields the following result.

Continuous-time single-input systems
In this section, we investigate how to get constructive methods for continuous-time single-input systems. First, we note that a direct application of the discrete-time methods above is not immediately possible. The reason is that, for u ∈ L 1 ([0, T ]) the solution to is given by which is not of the form p(A(θ))b(θ) for some polynomial p. In the following, we present an approach that enables us to trace the input construction problem for continuous-time systems back to the methods just presented in the previous section. The basic idea simply is to approximate the integral in the solution formula by a polynomial using piecewise constant inputs. The analysis will be carried out in six steps. Before we start, there are two comments in order. For continuoustime systems, the final time can be chosen arbitrarily, i.e. we can fix some T > 0 in advance. As a consequence, for continuous-time systems we are not restricted to start at the origin. Indeed, let x 0 ∈ C n (P) denote the initial states and let f ∈ C n (P) denote the final states, then usingf ∈ C n (P), defined bỹ Suppose the pair (A, b) satisfies the conditions (N1), (N2) and (S2) and let f ∈ C n (P) and ε > 0 be given. As T > 0 can chosen arbitrarily, we set T = N τ for some τ > 0 and some N ∈ N that will be specified in the construction process.
The significance of (38) is that it is independent of the number of input values u 0 , ..., u N −1 .

Remark 2
The approach in Section 4 was used in a first and direct proof of the sufficiency of (N1), (N2) and (S2) for uniform ensemble reachability for singleinput systems, cf. [28]. Indeed, since the set Ω is compact with empty interior and C \ Ω is connected, it is possible to applying Mergelyan's theorem directly to the continuous function Thus, there is a polynomial p such that |g(z) − τ p(z)| < ε ∀ z ∈ Ω.
Also we note that, the property that the interior of Ω is empty is a special case of Mergelyan's result, which was proven earlier by Lavrientev in 1934, cf. [24]. However, Lavrientev's proof is not available to the author. Maybe the proof can be used to obtain another constructive method, at least for special cases.

Comments
In this section we discuss some perspectives on the methods presented in Section 3.

Degrees of freedom
If the approximation domains are given by Jordan arcs, the methods for (S1) and (S2) have some degrees of freedom. For instance, there are plenty of ways to close a Jordan arc to a Jordan curve. This in turn frees up space for designing the conformal mappings such that the conditions (21), (23) and (28), (30) are satisfied. Thus, the computation of the conformal mappings and the design of the Jordan curves can be treated simultaneously. Besides, for pairs (A, b) where A is block diagonal and each block satisfies the sufficiency condition (S1) or (S2), the methods can be applied for every block and an overall polynomial can be constructed in the same manner the polynomials q 1 , ..., q n are determined in the methods for (S2).

Real approximation domains
In case the domains, where the approximations take place, are compact real intervals, the construction process is much easier. This can be achieved by considering a mixture of open-loop and feedback control inputs of the form K(θ)x t (θ) + u t .
Indeed, due to the necessary condition (N1), a suitable feedback term K(θ)x t (θ) can be used to shift the spectra of A(θ) to an appropriate position. That is, via K(θ) the matrices A(θ) + B(θ)K(θ) can have any spectra. Equivalently, the coefficients of the characteristic polynomials of A(θ) + B(θ)K(θ) can be arbitrarily designed. For more details on this we refer to [20,Section 6.4], [44].

Multi-input systems
The known necessary and sufficient conditions for ensemble reachability are in the multi-input case much less precise compared to the single-input systems. For special classes of systems, such as (A(θ), B(θ)) upper triangular, the problem can be traced back to the single-input case. In this context, a general but very strong sufficiency condition is that the Hermite indices of the pair (A(θ), B(θ)) are constant. Technicalities aside, these systems can be transformed into the Hermite canonical form   Ã 11 (θ) · · ·Ã 1k (θ) . . . . . .
and the construction problem can be solved by applying the single-input methods to the subpairs (Ã ii (θ),b i (θ)). We note that it is not required to use one of the methods for all subpairs. In particular, this applies to the often studied class (θA, B). For more details, we refer to [14] and the references therein.