Fixed-Endpoint Optimal Control of Bilinear Ensemble SystemsThis work was supported in part by the National Science Foundation under the awards CMMI-1301148, CMMI-1462796, and ECCS-1509342.

# Fixed-Endpoint Optimal Control of Bilinear Ensemble Systems††thanks: This work was supported in part by the National Science Foundation under the awards CMMI-1301148, CMMI-1462796, and ECCS-1509342.

Shuo Wang Department of Electrical and Systems Engineering, Washington University, St. Louis, Missouri, 63130, USA (swang35@wustl.edu).    Jr-Shin Li Department of Electrical and Systems Engineering, Washington University, St. Louis, Missouri, 63130, USA (jsli@wustl.edu). Questions, comments, or corrections to this document may be directed to this email address.
###### Abstract

Optimal control of bilinear systems has been a well-studied subject in the area of mathematical control. However, techniques for solving emerging optimal control problems involving an ensemble of structurally identical bilinear systems are underdeveloped. In this work, we develop an iterative method to effectively and systematically solve these challenging optimal ensemble control problems, in which the bilinear ensemble system is represented as a time-varying linear ensemble system at each iteration and the optimal ensemble control law is then obtained by the singular value expansion of the input-to-state operator that describes the dynamics of the linear ensemble system. We examine the convergence of the developed iterative procedure and pose optimality conditions for the convergent solution. We also provide examples of practical control designs in magnetic resonance to demonstrate the applicability and robustness of the developed iterative method.

Key words. Ensemble control, iterative methods, sweep method, fixed-endpoint problems, bilinear systems, optimality conditions, magnetic resonance.

AMS subject classifications.

## 1 Introduction

Newly emerging fields in science and engineering, such as systems neuroscience, synchronization engineering, and quantum science and technology, give rise to new classes of optimal control problems that involve underactuated manipulation of individual and collective behavior of dynamic units in a large ensemble. Representative examples include neural stimulation for alleviating the symptoms of neurological disorders such as Parkinson’s disease, where a population of neurons in the brain is affected by a small number of electrodes [Ching2013a]; pulse designs for exciting and transporting quantum systems between desired states, where an ensemble of quantum systems is driven by a single or multiple pulses in a pulse sequence [Brent06, Li_PNAS11]; and the engineering of dynamical structures for complex oscillator networks, where sequential patterns of a network of nonlinear rhythmic elements are created and altered by a mild global waveform [Kiss:2007:1886-1889]. Solving these nontraditional and large-scale underactuated control problems requires the development of systematic and computationally tractable and effective methods.

Among these emerging control problems, in this paper, we will study fixed-endpoint optimal control problems involving bilinear ensemble systems, which arise from the domain of quantum control [Chen14] and appear in a variety of other different fields, such as cancer chemotherapy [Schaettler2002] and robotics [Becker12]. The control of bilinear systems has been a well-studied subject in the area of mathematical control. From Pontryagin’s maximum principle to spectral collocation methods, a wide variety of theoretical and computational methods have been developed to solve optimal control problems of bilinear systems [Aganovic1994, Mohler2000]. In particular, the numerical methods are in principle categorized into direct, e.g., pseudospectral methods [Gong06, Ross04], and indirect approaches, e.g., indirect transcription method [Rao09] and shooting methods [Bertolazzi2005]. Implementing these existing numerical methods to solve optimal control problems involving an ensemble, i.e., a large number (finitely or infinitely many) or a parameterized family, of bilinear systems may encounter low efficiency, slow convergence, and instability issues, because most of these methods rely on suitable discretization of the continuous-time dynamics into a large-scale nonlinear program (LSNLP). In addition, the global constraint for such an optimal ensemble control problem, in which each individual system receives the same control input, makes the discretized LSNLP very restrictive and intractable to solve or even to find a feasible solution [Li_TAC12_QCP].

On the other hand, optimal control problems involving a linear system, or a linear ensemble system, are often computationally tractable and analytically solvable for many special cases, such as the linear quadratic regulator (LQR) [Brockett70] and the minimum-energy control of harmonic oscillator ensembles [Li_TAC11]. This suggests a bypass to solve optimal control problems of bilinear ensemble systems through solving that of linear ensemble systems and motivates the development of the iterative method in this work. The central idea is to represent the bilinear ensemble system as a linear ensemble system at each iteration, and then feasibly calculate the optimal control and trajectory for each iteration until a convergent solution is found. Iterative methods have been introduced and adopted to deal with diverse control design problems, including the free-endpoint quadratic optimal control of bilinear systems [Hofer1988] and optimal state tracking for nonlinear systems [Cimen2004], while the fixed-endpoint problems along with the emerging problems that involve controlling a bilinear ensemble system remain unexplored.

In this paper, we combine the idea of the aforementioned iterative method with our previous work on optimal control of linear ensemble systems to construct an iterative algorithm for solving optimal control problems involving a time-invariant bilinear ensemble system of the form,

 ddtX(t,β)=A(β)X(t,β)+B(β)u(t)+(m∑i=1ui(t)Bi(β))X(t,β),

where denotes the state, with compact and a positive integer, is the control, and the matrices , , and , , for .

This paper is structured as follows. In the next section, we present the developed iterative method for fixed-endpoint optimal control of a time-invariant bilinear system, where we introduce a sweep method that accounts for the terminal condition based on the notion of flow mapping from the optimal control theory. In Section LABEL:sec:convergence, we examine the convergence of the iterative method using the fixed-point theorem. In Section LABEL:sec:optimality, we propose the conditions for global optimality of the convergent solution. Then, in Section LABEL:sec:ensemble, we extend the developed iterative method to solve optimal control problems involving bilinear ensemble systems and show the convergence of the method. Finally, examples and simulations of practical control design problems are illustrated in Section LABEL:sec:examples to demonstrate the applicability and robustness of the developed iterative procedure.

## 2 Iterative method for optimal control of bilinear systems

We start with considering a fixed-endpoint, finite-time, quadratic optimal control problem involving a time-invariant bilinear system of the form

 min J=12∫tf0[xT(t)Qx(t)+uT(t)Ru(t)]dt, s.t. ˙x=Ax+Bu+[m∑i=1uiBi]x, \hb@xt@.01(P1) x(0)=x0,x(tf)=xf,

where is the state and is the control; , , and are constant matrices; is positive definite and is positive semi-definite; and are the initial and the desired terminal state, respectively. We first represent the time-invariant bilinear system in (LABEL:eq:oc1) as a time-varying linear system,

 ˙x(t)=Ax+Bu+[n∑j=1xj(t)Nj]u, \hb@xt@.01(2.1)

in which we write the bilinear term with the element of , for , and . Then, we solve this optimal control problem by Pontryagin’s maximum principle. The Hamiltonian of this problem is

 H(x,u,λ)=12(xTQx+uTRu)+λT{Ax+[B+(n∑j=1xjNj)]u}, \hb@xt@.01(2.2)

where is the co-state vector. The optimal control is then obtained by the necessary condition, , given by

 u∗=−R−1(B+n∑j=1xjNj)Tλ, \hb@xt@.01(2.3)

and the optimal trajectory of the state and the co-state satisfy, for ,

 ˙xi =[Ax]i−[(B+n∑j=1xjNj)R−1(B+n∑j=1xjNj)Tλ]i, \hb@xt@.01(2.4) ˙λi =−[Qx]i−[ATλ]i+λT{NiR−1(B+n∑j=1xjNj)T+(B+n∑j=1xjNj)R−1NTi}λ \hb@xt@.01(2.5)

with the boundary conditions and , where , and , , are the component of the associated vectors. By the following change of variables,

 ~Aij=Aij−[(NjR−1(B+n∑j=1xjNj)T+(B+n∑j=1xjNj)R−1NTj)λ]i, \hb@xt@.01(2.6) ~BR−1~BT=BR−1BT−(n∑j=1xjNj)R−1(n∑j=1xjNj)T, \hb@xt@.01(2.7) ~Q=Q, \hb@xt@.01(2.8)

we can rewrite (LABEL:eq:state) and (LABEL:eq:costate) into the form

 ˙x =~Ax−~BR−1~BTλ,x(0)=x0,x(tf)=xf, \hb@xt@.01(2.9) ˙λ =−~Qx−~ATλ, \hb@xt@.01(2.10)

which coincides with the canonical form of the state and co-state equations characterizing the optimal trajectories for the analogous optimal control problem involving the time-invariant linear system [Schaettler13]. In this way, the optimal state and co-state trajectories for the optimal control problem (LABEL:eq:oc1) involving a time-invariant bilinear system are now expressed in terms of the equations related to a time-varying linear system as in (LABEL:eq:state1) and (LABEL:eq:costate1).

Using this “linear-system representation” together with the Sweep method [Schaettler13, Bryson75], we will solve the optimal control problem (LABEL:eq:oc1) in an iterative manner. Specifically, we will consider at each iteration the fixed-endpoint linear quadratic optimal control problem,

 min J=12∫tf0[(x(k+1))T(t)Qx(k+1)(t)+(u(k+1))T(t)Ru(k+1)(t)]dt, s.t. ˙x(k+1)(t)=~A(k)x(k+1)+~B(k)u(k+1) \hb@xt@.01(P2) x(k+1)(0)=x0,x(k+1)(tf)=xf,

by treating the previous trajectory as a known quantity, where denotes the iteration. In the following sections, we will introduce the Sweep method and present the iterative procedure.

### 2.1 Sweep method for fixed-endpoint problems

Observe that in (LABEL:eq:state1) and (LABEL:eq:costate1) there are two boundary conditions for the state while none for the co-state . It requires implementing specialized computational methods, such as shooting methods, to solve such a two-point boundary value problem, which in general involve intensive numerical optimizations. Here, we adopt the idea of the Sweep method by letting

 λ(t)=K(t)x(t)+S(t)ν, \hb@xt@.01(2.11)

with , where and is the multiplier, a constant associated with the terminal constraint , which in this case is . From the transversality condition in Pontryagin’s maximum principle, we know that because there is no terminal cost and . Moreover, if is chosen to satisfy the Riccati equation

 ˙K(t)=−Q−~ATK(t)−K(t)~A+K(t)~BR−1~BTK(t), \hb@xt@.01(2.12)

with the terminal condition , then satisfies the matrix differential equation

 ˙S(t)=−(~AT−K(t)~BR−1~BT)S(t), \hb@xt@.01(2.13)

with the terminal condition , by taking the time derivative of (LABEL:eq:lam) and using (LABEL:eq:state1), (LABEL:eq:costate1) and (LABEL:eq:kRiccati). In addition, in order to fulfill the terminal condition at time , the multiplier associated with must satisfy

 xf=ST(t)x(t)+P(t)ν \hb@xt@.01(2.14)

for all , where obeys the matrix differential equation

 ˙P(t)−ST(t)~BR−1~BTS(t)=0, \hb@xt@.01(2.15)

with the terminal condition . It follows from (LABEL:eq:endpoint) using that

 ν=[P(0)]−1[xf−ST(0)x0], \hb@xt@.01(2.16)

provided is invertible for . More details about the Sweep method based on the notion of flow mapping are provided in Appendix LABEL:appd:sweep.

### 2.2 Iteration procedure

The optimal solution of the problem (LABEL:eq:oc1) is characterized by the homogeneous time-varying linear system described in (LABEL:eq:state1) and (LABEL:eq:costate1), and we will solve for and via an iterative procedure, which is based on analytical expressions and requires no numerical optimizations. To proceed this, we write (LABEL:eq:state1) and (LABEL:eq:costate1) as the iteration equations,

 ˙x(k+1)=~A(k)x(k+1)−~B(k)R−1(~B(k))Tλ(k+1), \hb@xt@.01(2.17) ˙λ(k+1)=−~Q(k)x(k+1)−(~A(k))Tλ(k+1), \hb@xt@.01(2.18)

with identical boundary conditions and for all , where , , and are defined according to (LABEL:eq:A), (LABEL:eq:B), and (LABEL:eq:Q), by

 ~A(k)ij=Aij−[(NjR−1(B+n∑j=1x(k)jNj)T+(B+n∑j=1x(k)jNj)R−1NTj)λ(k)]i, \hb@xt@.01(2.19) ~B(k)R−1(~B(k))T=BR−1BT−(n∑j=1x(k)jNj)R−1(n∑j=1x(k)jNj)T, \hb@xt@.01(2.20) ~Q(k)=Q. \hb@xt@.01(2.21)

Applying the Sweep method introduced in Section LABEL:sec:sweep, we let for , where satisfies the Riccati equation

 ˙K(k+1)=−Q(k)−K(k+1)~A(k)−(~A(k))TK(k+1)+K(k+1)~B(k)R−1(~B(k))TK(k+1), \hb@xt@.01(2.22)

with the boundary condition , and follows

 ˙S(k+1)=−[(~A(k))T−K(k+1)~B(k)R−1(~B(k))T]S(k+1),S(k+1)(tf)=I. \hb@xt@.01(2.23)

Moreover, the multiplier satisfies

 ν(k+1)=[P(k+1)(0)]−1[xf−(S(k+1))T(0)x0], \hb@xt@.01(2.24)

where is invertible (see Lemma LABEL:lem:P_inverse in Section LABEL:sec:convergence) and satisfies the dynamic equation

 ˙P(k+1)=(S(k+1))T~B(k)R−1(~B(k))TS(k+1) \hb@xt@.01(2.25)

with the terminal condition . Then, the optimal control (LABEL:eq:u) for the original Problem (LABEL:eq:oc1) can be expressed as

 u∗(t)=−R−1[B+n∑j=1x∗j(t)Nj]T[K∗(t)x∗(t)+S∗(t)ν∗], \hb@xt@.01(2.26)

if this iterative procedure is convergent, where , , and .

###### Remark 1

The iterative method can be initialized by conveniently using the optimal control of the system involving only the linear part of the bilinear system in (LABEL:eq:oc1), i.e., the LQR control. That is, the solution to the homogeneous system

 ˙x(0)=Ax(0)−BR−1BTλ(0),x(0)(0)=x0,x(0)(tf)=xf, ˙λ(0)=−ATλ(0).

However, the linear system may be uncontrollable so that the desired transfer between and is impossible and the LQR solution does not exist. In such a case, any state trajectory with the endpoints and can be a feasible initial trajectory of the iterative procedure.

### 2.3 A special case: minimum-energy control of bilinear systems

Before analyzing the convergence of the iterative method, we illustrate the procedure using the example of minimum-energy control of bilinear systems, which is a special case of Problem (LABEL:eq:oc1) with . Consider the following fixed-endpoint optimal control problem,

 min J=12∫tf0uT(t)Ru(t)dt, s.t. ˙x(t)=Ax+Bu+[n∑j=1xj(t)Nj]u, \hb@xt@.01(P3) x(0)=x0,x(tf)=xf.

The Hamiltonian of this problem is , where is the co-state vector. The optimal control is of the form as in (LABEL:eq:u) and the optimal state and co-state trajectories satisfy (LABEL:eq:state1) and (LABEL:eq:costate1), respectively, with . The respective iteration equations follow (LABEL:eq:x_k) and (LABEL:eq:lambda_k) with for all .

Following the iterative method presented in Section LABEL:sec:iterative, we represent the costate , , and the matrix satisfies the Riccati equation,

 ˙K(k+1)= −K(k+1)~A(k)−(~A(k))TK(k+1)+K(k+1)~B(k)R−1(~B(k))TK(k+1), \hb@xt@.01(2.27)

with the terminal condition , which has the trivial solution, , , and for . This gives

 λ(k+1)(t)=S(k+1)(t)ν(k+1), \hb@xt@.01(2.28)

and satisfies

 ˙S(k+1)=−(~A(k))TS(k+1),S(k+1)(tf)=I. \hb@xt@.01(2.29)

In addition, the multiplier associated with the terminal constraint is expressed as in (LABEL:eq:nuk). Combining (LABEL:eq:lambda_min_energy) with (LABEL:eq:u) gives the minimum-energy control at the iteration,

 (u∗)(k+1)(t)=−R−1[B+n∑j=1x(k+1)jNj]TS(k+1)ν(k+1). \hb@xt@.01(2.30)

Note that the auxiliary variable at each iteration satisfies (LABEL:eq:Pk), and thus

 P(k+1)(0)=−Φ~A(k)(tf,0)W(k+1)ΦT~A(k)(tf,0),

where is the transition matrix for the homogeneous equation (LABEL:eq:S(k+1)) and

 W(k+1)=∫tf0Φ~A(k)(0,σ)~B(k)R−1(~B(k))TΦT~A(k)(0,σ)dσ

is the controllability Gramian for the time-varying linear system as in Problem (LABEL:eq:ocme), or, equivalently, as in (LABEL:eq:state1) and (LABEL:eq:costate1) with . Moreover, the closed-loop expression in (LABEL:eq:ocume) is consistent with the open-loop expression of the minimum-energy control in terms of the controllability Gramian, that is,

 (u∗)(k+1)=R−1[B+n∑j=1x(k+1)jNj]TΦT~A(k)(0,t)(W(k+1))−1ξ(k), \hb@xt@.01(2.31)

where .

## 3 Convergence of the Iterative Method

Following the iterative algorithm described in Section LABEL:sec:iterative, we expect to find the optimal control for Problem (LABEL:eq:oc1), provided the iterations are convergent. In this section, we show that the convergence of this algorithm is pertinent to the controllability of the linear system considered at each iteration and depends on the choice of the weight matrix . In Section LABEL:sec:ensemble, we will extend this iterative method to solve optimal control problems involving bilinear ensemble systems.

To facilitate the proof, we introduce the following mathematical tools. Considering the Banach spaces, , , and with the norms

 ∥x∥α =supt∈[0,tf][∥x(t)∥exp(−αt)], for  x∈X, \hb@xt@.01(3.1) ∥y∥α =supt∈[0,tf][∥y(t)∥exp(−α(tf−t))], for  y∈Y, \hb@xt@.01(3.2) ∥z∥α =supt∈[0,tf][∥z(t)∥exp(−α(tf−t))], for  z∈Z, \hb@xt@.01(3.3)

in which for and for , and the parameter serves as an additional degree of freedom to control the rate of convergence [Hofer1988], we define the operators , , and that characterize the dynamics of , , and as described in Section LABEL:sec:iterative, given by

 ddtT1[x,K,Sν](t) =~A(x(t),K(t),S(t)ν)T1[x,K,Sν](t)−~B(x(t))R−1~BT(x(t))T3[x,K,Sν](t) −~B(x(t))R−1~BT(x(t))T2[x,K,Sν](t)T1[x,K,Sν](t), \hb@xt@.01(3.4) T1[x,K,Sν](0) =x0 ddtT2[x,K,Sν](t) =−Q+T2[x,K,Sν](t)~B(x(t))R−1~BT(x(t))T2[x,K,Sν](t) −T2[x,K,Sν](t)~A(x(t),K(t),S(t)ν)−~AT(x(t),K(t),S(t)ν)T2[x,K,Sν](t) \hb@xt@.01(3.5) T2[x,K,Sν](tf) =0, ddtT3[x,K,Sν](t) =−[~AT(x(t),K(t),S(t)ν)−T2[x,K,Sν](t)~B(x(t))R−1~BT(x(t))]⋅ T3[x,K,Sν](t), \hb@xt@.01(3.6) T3[x,K,Sν](tf) =ν(T1[x,K,Sν],T2[x,K,Sν],T3[x,K,Sν])

where is the multiplier satisfying (LABEL:eq:nuk). With these definitions and the following lemma, the convergence of the iterative method can be developed using the fixed-point theorem.

###### Lemma 3.1

The matrix as in (LABEL:eq:Pk) is nonsingular over at each iteration if and only if the time-varying linear system in Problem (LABEL:eq:oc2) is controllable over [Schaettler13].

Proof: See Appendix LABEL:appd:Pnonsingular.

###### Theorem 3.2

Consider the iterative method with the iterations evolving according to

 x(k+1)(t) =T1[x(k),K(k),S(k)ν(k)](t), \hb@xt@.01(3.7) K(k+1)(t) =T2[x(k),K(k),S(k)ν(k)](t), \hb@xt@.01(3.8) S(k+1)(t)ν(k+1) =T3[x(k),K(k),S(k)ν(k)](t), \hb@xt@.01(3.9)

where the operators , , and are defined in (LABEL:eq:T1), (LABEL:eq:T2), and (LABEL:eq:T3), respectively. If at each iteration the linear system as in (LABEL:eq:oc2) is controllable, then , , and are contractive. Furthermore, starting with a triple of feasible trajectories , the iteration procedure is convergent, and the sequences , and converge to the unique fixed points, , , and , respectively.

Proof: Because the linear system in (LABEL:eq:oc2) is controllable at each iteration , by Lemma LABEL:lem:P_inverse the matrix defined in (LABEL:eq:Pk) is invertible and hence the multiplier expressed in (LABEL:eq:nuk) is well-defined. Then, we have, at time , , since .

From (LABEL:eq:Ak) and (LABEL:eq:Bk), for each fixed , we obtain the bounds

 ∥~A(k+1)−~A(k)∥ ≤[n∑i=1∥Gi∥2]1/2∥λ(k+1)−λ(k)∥ +∥[n∑i,j=1∥Hij∥2]1/2{∥λ(k+1)∥∥x(k+1)−x(k)∥+∥x(k)∥∥λ(k+1)−λ(k)∥}, ∥~B(k+1)R−1(~B (k+1))T−~B(k)R−1(~B(k))T∥≤∥[n∑i,j=1∥Hij∥2]1/2∥(x(k+1))2−(x(k))2∥,

where and , and from (LABEL:eq:Qk), we have for all . Substituting (LABEL:eq:lambda_k) into the above inequalities, we can write these bounds in terms of , and , given by

 ∥~A(k+1)−~A(k) ∥≤{[n∑i=1∥Gi∥2]1/2+∥x(k)∥[n∑i,j=1∥Hij∥2]1/2}⋅ {∥K(k+1)∥∥x(k+1)−x(k)∥+∥K(k+1)−K(k)∥∥x(k)∥+∥S(k+1)ν(k+1)−S(k)ν(k)∥} \hb@xt@.01(3.10) +[n∑i,j=1∥Hij∥2]1/2∥K(k+1)∥∥x(k+1)∥+∥S(k+1)ν(k+1)∥∥x(k+1)−x(k)∥, ∥~B(k+1)R−1(~B (k+1))T−~B(k)R−1(~B(k))T∥≤[n∑i,j=1∥Hij∥2]1/2⋅ {∥x(k+1)∥+∥x(k)∥}∥x(k+1)−x(k)∥. \hb@xt@.01(3.11)

In addition, the solution to (LABEL:eq:Sk) is given by

 S(k+1)(t) =Φ−[(~A(k))T−K(k+1)~B(k)R−1(~B(k))T](t,tf)S(k+1)(tf) =ΦT[~A(k)−~B(k)R−1(~B(k))TK(k+1)](tf,t), \hb@xt@.01(3.12)

where denotes the transition matrix associated with the homogeneous system (LABEL:eq:Sk) and . Then, we have

 ∥S(k+1)(t)−S(k)(t)∥ ≤∫Tt∥[(S(k+1))T(t)]−1∥∥S(k+1))T(σ)∥[∥~A(k)(σ)−~A(k−1)(σ)∥ +∥~B(k)R−1(~B(k))T(σ)−~B(k−1)R−1(~B(k−1))T(σ)∥∥K(k+1)(