A Gradient Descent Approach to Optimal Coherent Quantum LQG Controller Design Arash Kh. Sichani, Igor G. Vladimirov, Ian R. Petersen arXiv:1502.00274v1 [quant-ph] 1 Feb 2015 Abstract This paper is concerned with the Coherent Quantum Linear Quadratic Gaussian (CQLQG) control problem of finding a stabilizing measurement-free quantum controller for a quantum plant so as to minimize an infinite-horizon mean square performance index for the fully quantum closed-loop system. In comparison with the observation-actuation structure of classical controllers, the coherent quantum feedback is less invasive to the quantum dynamics and quantum information. Both the plant and the controller are open quantum systems whose dynamic variables satisfy the canonical commutation relations (CCRs) of a quantum harmonic oscillator and are governed by linear quantum stochastic differential equations (QSDEs). In order to correspond to such oscillators, these QSDEs must satisfy physical realizability (PR) conditions, which are organised as quadratic constraints on the controller matrices and reflect the preservation of CCRs in time. The CQLQG problem is a constrained optimization problem for the steady-state quantum covariance matrix of the plant-controller system satisfying an algebraic Lyapunov equation. We propose a gradient descent algorithm equipped with adaptive stepsize selection for the numerical solution of the problem. The algorithm finds a local minimum of the LQG cost over the parameters of the Hamiltonian and coupling operators of a stabilizing PR quantum controller, thus taking the PR constraints into account. A convergence analysis of the proposed algorithm is presented. A numerical example of a locally optimal CQLQG controller design is provided to demonstrate the algorithm performance. I. INTRODUCTION Coherent quantum feedback control [11], [13] is a relatively novel quantum control paradigm which is aimed at achieving given performance specifications for quantum systems, such as internal stability and optimization of a cost functional. Such systems arise naturally in quantum physics [7] and its engineering applications (for example, nanotechnology and quantum optics [6]). The dynamic variables of quantum systems are (usually noncommuting) operators on an underlying Hilbert space which evolve according to the laws of quantum mechanics [15]. The latter make the quantum dynamics particularly sensitive to interaction with classical devices in the course of quantum measurement, as reflected in the projection postulate of quantum mechanics. In order to overcome this issue, coherent quantum control employs the idea of direct interconnection of quantum systems to be controlled (quantum plants) with other quantum systems playing the role of controllers, possibly mediated by light fields. Unlike the traditional observation-actuation control loop, this fully quantum measurement-free feedback avoids the loss of quantum information as a result of its conversion to classical signals. Quantum-optical components, such as optical cavities, beam splitters and phase shifters, make it possible to implement coherent quantum feedback governed by linear quantum stochastic differential equations (QSDEs) [17], [18], provided the latter are physically realizable (PR) as open quantum harmonic oscillators [5], [6]. The resulting PR conditions [9], [20], [21] are organized as quadratic constraints on the coefficients of the QSDEs. The PR constraints for the state-space matrices of a coherent quantum controller complicate the solution of quantum counterparts to the classical H∞ and Linear Quadratic Gaussian (LQG) control problems. The Coherent Quantum LQG (CQLQG) control problem [16] seeks for a stabilizing PR quantum controller so as to minimize a mean square performance criterion for the fully quantum closed-loop system. A numerical procedure for finding suboptimal controllers for this problem was proposed in [16], and algebraic equations for the optimal CQLQG controller were obtained in [25]. Despite the previous results, the CQLQG control problem does not lend itself to an “elegant” solution (for example, in the form of decoupled Riccati equations as in the classical case [10]) and remains a subject of research. Since the main difficulties are caused by the coupling of the equations due to the PR constraints, a conversion of the CQLQG control problem to an unconstrained problem by using Lagrange multipliers was considered in [26] for a related coherent quantum filtering problem which is a simplified feedback-free version of the CQLQG control problem. In the present paper, we develop an algorithm for the numerical solution of the CQLQG control problem by using the gradient descent method and the Hamiltonian parameterization of PR quantum controllers [25]. The latter is a different technique to handle the PR constraints by reformulating the CQLQG control problem in an unconstrained fashion. More precisely, the optimal solution is sought in the class of stabilizing PR controllers whose state-space matrices are parameterized in terms of the free Hamiltonian and coupling operators of an open quantum harmonic oscillator [5]. We obtain ordinary differential equations (ODEs) for the gradient descent in the Hilbert space of matrix-valued parameters. For this purpose, Fr´echet differentiation is used together with related algebraic techniques [3], [14], [22], [23], [25] to employ the analytic structure of the LQG cost as a composite function of the matrix-valued variables, involving Lyapunov equations. The advantage of the proposed approach is that, at intermediate steps, it produces stabilizing PR quantum controllers which can This work is supported by the Australian Research Council. The authors are with UNSW Canberra, ACT 2600, Australia. E-mail: arash kho@hotmail.com, igor.g.vladimirov@gmail.com, i.r.petersen@gmail.com. be used as gradually improving suboptimal solutions of the CQLQG control problem, and a locally optimal solution (if it exists) is achieved asymptotically by moving along negative gradient directions with a suitable choice of stepsizes. To this end, we provide an algorithm for adaptive selection of the stepsize for each iteration based on the second-order Gˆateaux derivative of the LQG cost along the gradient. However, the proposed gradient descent algorithm for the CQLQG control problem requires for its initialization a stabilizing PR quantum controller. Finding such a controller for an arbitrary given quantum plant is a nontrivial open problem. Because of the lack of a systematic solution for this quantum stabilization problem at the moment, the current version of the algorithm is initialized at a stabilizing PR quantum controller obtained by random search in the space defined by the Hamiltonian parameterization of PR controllers. Although a random search for an admissible starting point is acceptable for low-dimensional problems, the development of a more systematic solution for this issue is a subject of future research. The paper is organised as follows. Section II outlines the notation used in the paper. Sections III and IV specify the quantum plants and coherent quantum controllers being considered. Section V revisits PR conditions for linear quantum systems. Section VI formulates the CQLQG control problem. Section VII describes a gradient descent system for finding local minima in the control problem. Section VIII describes an algorithmic implementation of the gradient descent method. Section IX discusses convergence of the algorithm. Section X provides a numerical example of designing a locally optimal CQLQG controller. Section XI gives concluding remarks. Appendices A and B provide a subsidiary material on the differentiation of the LQG cost. II. NOTATION Vectors are assumed to be organized as columns unless specified otherwise, and the transpose (·)T acts on matrices with operator-valued entries as if the latter were scalars. For a vector X of operators X1 , . . . , Xr and a vector Y of operators Y1 , . . . ,Ys , the commutator matrix [X,Y T ] := XY T − (Y X T )T is an (r × s)-matrix whose ( j, k)th entry is the commutator [X j ,Yk ] := X jYk − Yk X j of the operators X j and Yk . Furthermore, (·)† := ((·)# )T denotes the transpose of the entry-wise operator adjoint (·)# . When it is applied to complex matrices, (·)† reduces to the complex conjugate transpose (·)∗ := ((·))T . T T Denoted by sym(·) := (·)+(·) and asym(·) := (·)−(·) are the symmetrizer and antisymmetrizer of matrices. Also, we denote 2 2 by Sr , Ar and Hr := Sr + iAr the subspaces of real symmetric, real antisymmetric and complex Hermitian matrices of √ order r, respectively, with i := −1 the imaginary unit. Furthermore, Ir denotes the identity matrix of order r, positive (semi-) definiteness of matrices is denoted by (<) , and ⊗ is the tensor product of spaces or operators (in particular, the Kronecker product of matrices). The adjoints and self-adjointness of linear operators acting on matrices is understood in the sense of thepFrobenius inner product hM, Ni := Tr(M ∗ N) of real or complex matrices, with the corresponding Frobenius norm kMk := hM, Mi which is the standard Euclidean norm | · | for vectors. Also, EX := Tr(ρX) denotes the quantum expectation of a quantum variable X (or a vector of such variables) over a density operator ρ which specifies the underlying quantum state. III. QUANTUM PLANT The quantum plant under consideration is an open quantum stochastic system which is coupled to another such system (playing the role of a controller), with the dynamics of both systems being affected by the environment. Both the plant and the controller are assumed to satisfy the physical realizability (PR) conditions [9], [16], [20] which will be described in Section V. The plant has n dynamic variables x1 (t), . . . , xn (t) (with n even) which are self-adjoint operators on a Hilbert space specified below. With the time arguments being omitted for brevity, the evolution of the plant state vector x := (xk )16k6n and its contribution to a p1 -dimensional output of the plant y := (yk )16k6p1 (also with self-adjoint operator-valued entries) are governed by QSDEs dx = Axdt + Bdw + Edη, dy = Cxdt + Ddw. (1) Here, A ∈ Rn×n , B ∈ Rn×m1 , C ∈ R p1 ×n , D ∈ R p1 ×m1 , E ∈ Rn×p2 are given constant matrices. Also, z := Cx (2) is a “signal part” of the plant output y, and η is a p2 -dimensional output of the controller to be described in Section IV. The external noise acting on the plant is represented by a quantum Wiener process w := (wk )16k6m1 whose entries are self-adjoint operators on a boson Fock space F1 [17] with the quantum Itˆo table dwdwT = Ω1 dt, where the matrix Ω1 ∈ Hm1 is given by Ω1 := Im1 + iJ1 < 0. Here, the matrix J1 ∈ Am1 specifies the CCRs between the entries of the plantnoise w as 0 1 [dw, dwT ] = 2iJ1 dt and (assuming that the dimension m1 is even) is given by J1 := Im1 /2 ⊗ J, with J := . −1 0 2 IV. QUANTUM CONTROLLER Consider an interconnection of the plant (1) with a coherent (that is, measurement-free) quantum controller. The latter is also an open quantum system with an n-dimensional state vector ξ := (ξk )16k6n of self-adjoint operators on another Hilbert space, which also evolve in time. The assumption that the controller has the same number of dynamic variables as the plant is adopted from the classical LQG control theory. The controller interacts with the plant and the environment according to the QSDEs dξ = aξ dt + bdω + edy, dη = cξ dt + ddω. (3) Here, a ∈ Rn×n , b ∈ Rn×m2 , c ∈ R p2 ×n , d ∈ R p2 ×m2 , e ∈ Rn×p1 are also constant matrices. Similarly to (2), the p2 -dimensional process ζ := cξ (4) is the signal part of the controller output η. The process ω in (3) is a quantum noise which effects the controller and is an m2 -dimensional quantum Wiener process (with m2 even) on another boson Fock space F2 with the quantum Ito table dωdω T = Ω2 dt, where the matrix Ω2 ∈ Hm2 is given by Ω2 := Im2 + iJ2 < 0. Here, the matrix J2 ∈ Am2 specifies the CCRs between the entries of the controller noise ω as [dω, dω T ] = 2iJ2 dt and is given by J2 := Im2 /2 ⊗ J. The plant and controller noises w and ω act on different boson Fock spaces F1 and F2 , respectively, and hence, commute with each other. Therefore, the combined quantum Wiener process w W := (5) ω of dimension m := m1 + m2 acts on the tensor product space F1 ⊗ F2 and has a block diagonal CCR matrix J: J1 0 T [dW , dW ] = 2iJdt, J := . 0 J2 (6) Furthermore, the external boson fields are assumed to be in the product vacuum state υ := υ1 ⊗υ2 , and hence, are uncorrelated. The resulting quantum Ito table of the combined Wiener process W in (5) is dW dW T = Ωdt, Ω := Im + iJ = Ω∗ < 0. (7) In the controller dynamics (3), the matrix b is the noise gain matrix, while e plays the role of the observation gain matrix, although y is not an observation signal in the classical control theoretic sense. Accordingly, the process ζ in (4) corresponds to the classical actuator signal. Now, the combined set of QSDEs (1), (3) describes the fully quantum closed-loop system shown in Fig. 1. By using a quadratic cost adopted in quantum control settings [16], [25] from classical LQG control [10], η w - quantum plant quantum controller ω y Fig. 1. The interconnected quantum plant and quantum controller form a fully quantum closed-loop system which is governed by (1), (3) and is influenced by the environment through the quantum Wiener processes w, ω. the performance of the coherent quantum controller will be described in Section VI in terms of an r-dimensional quantum process Z := Fx + Gζ . (8) Here, F ∈ Rr×n , G ∈ Rr×p2 are given weighting matrices whose entries quantify the relative importance of the state variables x1 , . . . , xn of the plant and the “actuator output” variables ζ1 , . . . , ζ p2 of the controller. The choice of F, G is specified only by control design preferences and is not subjected to physical constraints. The process Z in (8) is linearly related to the 2n-dimensional state vector x X := (9) ξ of the closed-loop system whose dynamics are governed by the QSDE dX = A X dt + BdW , 3 Z = CX (10) which is driven by the combined quantum Wiener process W from (5). The state-space matrices A ∈ R2n×2n , B ∈ R2n×m , C ∈ Rr×2n of the closed-loop system in (10) are obtained by combining (1), (3) with (5), (8), (9) and depend on the controller matrices a, b, c, e in an affine fashion: A Ec B Ed A := , B := , C := F Gc . (11) eC a eD b V. CONDITIONS FOR PHYSICAL REALIZABILITY Both the quantum plant (1) and the coherent quantum controller (3) are assumed to be physically realizable as open quantum harmonic oscillators, with initial complex separable Hilbert spaces H1 , H2 . In particular, their dynamic variables (which are self-adjoint operators on the product space H1 ⊗ H2 ⊗ F1 ⊗ F2 at any subsequent moment of time t > 0) satisfy CCRs [x, xT ] = 2iΘ1 , [ξ , ξ T ] = 2iΘ2 , [x, ξ T ] = 0, (12) where Θ1 , Θ2 ∈ An are constant nonsingular matrices. An equivalent form of the CCRs for the combined vector X from (9) is Θ 0 [X , X T ] = 2iΘ, Θ := 1 . (13) 0 Θ2 The preservation of the CCRs (12) (including the commutativity between x and ξ ) is a consequence of the unitary evolution of the isolated system formed from the plant, controller and their environment. The QSDE in (10) preserves the CCR matrix Θ in (13) in time if and only if the matrices A , B in (11) satisfy A Θ + ΘA T + BJB T = 0, (14) where J is the CCR matrix of the combined quantum Wiener process W in (5) given by (6). The relation (14) is obtained by taking the imaginary part of the algebraic Lyapunov equation (ALE) A S + SA T + BΩB T = 0 (15) (provided A is Hurwitz) for the steady-state quantum covariance matrix S := lim E(X (t)X (t)T ) = P + iΘ = S∗ < 0, t→+∞ (16) with Ω the quantum Ito matrix from (7). Here, the quantum expectation E(·) is taken over the product state ϖ ⊗ υ, where ϖ is the initial quantum state of the plant and controller on H1 ⊗ H2 , and υ is the vacuum state of the external fields on F1 ⊗ F2 . We have also used the convergence limt→+∞ EX (t) = 0 which is ensured by A being Hurwitz. The real part P := ReS (17) of the quantum covariance matrix S from (16) is the unique solution to the ALE A P + PA T + BB T = 0, (18) obtained by taking the real part of (15), and coincides with the controllability Gramian [10] of the pair (A , B). Since the left-hand side of (14) is an antisymmetric matrix of order 2n, then, by computing the diagonal (n × n)-blocks and the upper off-diagonal block of this matrix with the aid of (11), it follows that the preservation of the CCR matrix Θ in (13) is equivalent to AΘ1 + Θ1 AT + BJ1 BT + EdJ2 d T E T = 0, T T T T aΘ2 + Θ2 a + eDJ1 D e + bJ2 b = 0, T T T T (Θ1C + BJ1 D )e + E(cΘ2 + dJ2 b ) = 0; (19) (20) (21) cf. [24, Eqs. (18)–(20)]. Therefore, the fulfillment of the equalities Θ1CT + BJ1 DT = 0, T cΘ2 + dJ2 b = 0 (22) (23) is sufficient for (21). Note that (19), (22) are the conditions for physical realizability (PR) [9], [16], [20] of the quantum plant which describe the preservation of the CCR matrix Θ1 in (12) and [x, yT ] = 0. Similarly, the relations (20), (23), which describe the preservation of the CCR matrix Θ2 in (12) and [ξ , η T ] = 0, are the PR conditions for the coherent quantum controller. The PR condition (20) can be regarded as a linear equation with respect to the matrix a, and its general solution is representable as 1 (24) a = 2Θ2 R − (eDJ1 DT eT + bJ2 bT )Θ−1 2 . 2 4 Here, the matrix R ∈ Sn specifies the free Hamiltonian 21 ξ T Rξ which the PR controller would have in the absence of interaction with its surroundings; cf. [5, Eqs. (20)–(22) on pp. 8–9]. The other PR condition (23) allows the matrix c to be expressed in terms of b as c = −dJ2 bT Θ−1 (25) 2 . The coupling between the output matrix c and the noise gain matrix b makes the design of a coherent quantum controller (3) substantially different from that of the classical controllers even at the level of achieving internal stability for the closed-loop system. Indeed, if the additional quantum noise ω is effectively eliminated from the state dynamics of the quantum controller by letting b = 0, then (25) implies that c = 0, and hence, the matrix A in (11) becomes block lower triangular. In this case, the closed-loop system in (10) cannot be internally stable if A is not Hurwitz. Also note that, in the formulations of the PR conditions [9], [16], [21] for the plant and controller QSDEs (1), (3), the noise feedthrough matrices are usually given by D = I p1 0 and d = I p2 0 , with p1 6 m1 and p2 6 m2 . Such matrices D and d have full row rank and satisfy DDT = I p1 , dd T = I p2 . (26) The full row rank property of D corresponds to nondegeneracy of measurements in the classical setting, where y in (1) is an observation process. Furthermore, since det J2 6= 0 and det Θ2 6= 0, the full row rank property of d implies that the map Rn×m2 3 b 7→ c ∈ R p2 ×n , given by (25), is onto. This allows the matrix c to be assigned any value by an appropriate choice of b, which plays a part in the stabilization issue mentioned above. Although (26) simplifies the algebraic manipulations, it is the rank properties of the matrices D, d that are most important. VI. COHERENT QUANTUM LQG CONTROL PROBLEM Following [16], [25], we formulate the CQLQG control problem as that of minimizing the steady-state mean square value 1 1 (27) E := lim E(Z (t)T Z (t)) = C T C , P −→ min t→+∞ 2 2 for the criterion process Z of the closed-loop system (10) over internally stabilizing (that is, making the matrix A Hurwitz) PR quantum controllers (3) of fixed dimensions described in Sections IV, V. Here, Z T Z = ∑rk=1 Zk2 is the sum of squared entries of Z (and hence, Z T Z is a positive semi-definite self-adjoint operator) and P is the controllability Gramian of the closed-loop system given by (17), (18). The LQG cost E in (27) is a function of the triple u := (R, b, e) ∈ Sn × Rn×m2 × Rn×p1 =: U (28) which parameterizes PR quantum controllers (3) through (24), (25), with the controller noise feedthrough matrix d ∈ R p2 ×m2 being fixed and satisfying (26). Accordingly, the minimization in (27) is carried out over the set U0 := {u ∈ U : A in (11) is Hurwitz} (29) of those u which specify internally stabilizing PR quantum controllers for the quantum plant (1). For what follows, the set U on the right-hand side of (28) is endowed with the structure of a Hilbert space with the direct sum inner product h(R, b, e), (R0 , b0 , e0 )i := hR, R0 i + hb, b0 i + he, e0 i. Note that U0 in (29) is an open subset of U. VII. GRADIENT FLOW FOR THE LQG COST The gradient descent approach to the solution of the CQLQG control problem is to move with the negative gradient flow for the LQG cost function E in (27) towards a local minimum. The gradient descent can be regarded as a dynamical system governed by the ODE u(τ) ˙ = −g(u(τ)), u(0) = u0 . (30) Here, (˙) := ∂τ (·) is the derivative with respect to fictitious time τ > 0, and the gradient g(u) := ∂u E (u) = (∂R E , ∂b E , ∂e E ) (31) is the Fr´echet derivative of the LQG cost with respect to u in the Hilbert space U associated with the Hamiltonian parameterization of PR quantum controllers in (28). More precisely, the map g : U0 → U is well-defined on the open set U0 in (29). The starting point in (30) is assumed to satisfy u0 := (R0 , b0 , e0 ) ∈ U0 , (32) so that the corresponding PR controller is internally stabilizing. Unless u0 is a stationary point of E , the LQG cost is strictly • decreasing along the trajectory of the ODE (30) in view of E (u(τ)) = −kg(u(τ))k2 . The first-order Fr´echet derivative in (31) is computed in the following lemma. For its formulation, we denote by Q the observability Gramian of the pair (A , C ) which is a unique solution of the ALE A T Q + QA + C T C = 0, (33) 5 provided the matrix A in (11) is Hurwitz. Furthermore, we will use the Hankelian of the closed-loop system defined by H := QP. (34) ←n→←n→ n X X l Also, we partition (2n × 2n)-matrices X (such as P, Q, H) into (n × n)-blocks X jk as X := X11 X12 l . n 21 22 Lemma 1: For any u ∈ U0 from (29), the Fr´echet derivative (31) of the LQG cost E in (27) can be computed as ∂R E = −2sym(Θ2 H22 ), (35) ∂b E = Q21 Ed + Q22 b − ψbJ2 − χdJ2 , (36) ∂e E = H21C + Q21 BD + Q22 e − ψeDJ1 D . T T T (37) Here, ψ and χ are auxiliary (n × n)-matrices defined by ψ := asym(H22 Θ−1 2 ), T T T T χ := Θ−1 2 (H12 E + P21 F G + P22 c G G), (38) with P, Q, H the Gramians and Hankelian from (18), (33), (34). The proof of Lemma 1 is similar to that of [25, Theorem 1] and is given in Appendix A for completeness. That the trajectories of the gradient descent system in (30) will not “miss” local minima of the LQG cost is justified by the following lemma. Lemma 2: A point u∗ ∈ U0 in (29) is a stable equilibrium of the ODE (30) if and only if it is a local minimum of the LQG cost E in (27). Proof: The assertion of the lemma can be established by using [1, Theorem 3] and the analyticity [8] (rather than infinite differentiability) of the LQG cost E in a neighbourhood of any point u ∈ U0 . The analyticity follows from the representation 1 E = − vec(C T C )T (A ⊕ A )−1 vec(BB T ) (39) 2 which is obtained from (18), (27) by using the column-wise vectorization vec(·) of matrices [14], [22] and the Kronecker sum A ⊕ A := I2n ⊗ A + A ⊗ I2n of the matrix A with itself. Indeed, the representation (39) implies that E is a rational function of the entries of A , B, C in (11) which, in turn, are polynomial functions of the entries of R, b, e in view of (24), (25), and hence, E (u) is a rational function of u. Therefore, the function E (u) is analytic on the open set U0 since the matrix A ⊕ A is also Hurwitz (and hence, nonsingular) for any u ∈ U0 . In practice, the gradient descent ODE (30) is solved by using a numerical algorithm, which is the subject of the next section. VIII. GRADIENT DESCENT ALGORITHM We will now consider a numerical algorithm which implements the gradient descent method (30) for the CQLQG control problem in the form uk+1 := uk − sk g(uk ), k = 0, 1, 2, . . . . (40) This recurrence equation is initialized with matrices R0 , b0 , e0 of an internally stabilizing PR controller in (32) (see Section VIII-A). The gradient g(uk ) is computed by using Lemma 1, and the stepsize sk > 0 is chosen as described in Section VIII-B. The iterations in (40) are stopped when a termination condition is satisfied (see Section VIII-C). The ingredients of the algorithm are discussed in the subsequent sections. A. Initialization The gradient descent algorithm (40) relies on existence of an internally stabilizing PR quantum controller which can be used as an initial point. As mentioned in Introduction, the existence of such controllers (that is, nonemptiness of the set U0 in (29)) for a given quantum plant (and a systematic method of finding them) remains an open problem. In the present version of the algorithm, this quantum stabilization problem is solved by using a random search in the finite-dimensional Hilbert space U in (28). 6 B. Stepsize selection According to the conventional limited minimization rule (see, for example, [4]), the stepsize sk is chosen for each iteration of the gradient descent by solving the minimization problem sk ∈ Arg min E (uk − sg(uk )) (41) 06s6hk with a constant search horizon hk := h > 0. Here, we use the convention that E (u) := +∞ if u 6= U0 (thus discarding those controllers which are not internally stabilizing). A restricted version of the line search with a constant horizon h may suffer from the inability to adapt properly to the behaviour of the function E in its minimization over the ray {uk − sg(uk ) : s > 0} ⊂ U. In order to overcome this issue, for the stepsize selection in the gradient descent algorithm (40), we will use a modified version of the limited minimization rule with an adaptive choice of the search horizon hk in each iteration. More precisely, hk can be chosen so as to enable (41) to “capture” the minimum of E over the whole ray if E is a strictly convex quadratic function. To this end, consider the quadratically truncated Taylor series E (u − sg) = E (u) − sDg E + s2 2 D E + o(s2 ), 2 g (42) where Dv E (u) := ∂s E (u + sv)s=0 = hg(u), vi, Dv2 E (u) := ∂s2 E (u + sv) = h∂u2 E (u)(v), vi s=0 (43) (44) are the first and second-order Gˆateaux (or directional) derivatives [12] of the LQG cost at a point u ∈ U0 (specifying an internally stabilizing controller) along v ∈ U. Here, in view of (43), Dg E = kgk2 = k∂R E k2 + k∂b E k2 + k∂e E k2 > 0. (45) Also, ∂u2 E (u) := ∂u g(u) in (44) is the second-order Fr´echet derivative of E which is a self-adjoint operator on the Hilbert space U in (28) whose computation is outlined in Appendix B. Now, if Dg2 E (u) > 0, then the quadratic polynomial of s on the right-hand side of (42) (with the higher-order terms being neglected) achieves its unique minimum at a nonnegative value of s: 2 Dg E kgk2 s 2 (46) Dg E − sDg E = 2 = 2 . arg min 2 Dg E Dg E s>0 This suggests using the right-hand side of (46) as a search horizon hk in (41), provided Dg2 E (u) > 0. However, if the latter inequality does not hold, the argument, based on a quadratic approximation of the minimization problem (41), is no longer valid and needs to be amended. In this case (when Dg2 E (u) 6 0), the search horizon can be chosen so as to avoid the domination of nonlinear terms over the linear term in the quadratically truncated Taylor series for the LQG cost along the ideal gradient descent trajectory in (30): s2 + o(s2 ) 2 = E (u(τ)) − kgk2 s + Dg2 E s2 + o(s2 ), • •• E (u(τ + s)) = E (u(τ)) + (E (u)) s + (E (u)) (47) where (44) is used. For s > 0, the comparison of the absolute values kgk2 s and |Dg2 E |s2 of the linear and quadratic terms in (47) shows that the latter does not dominate the former if s6 kgk2 . |Dg2 E | (48) This inequality is closely related to the accuracy of (40) as Euler scheme for numerical integration of the ODE (30). More precisely, if the stepsizes sk > 0 in (40) are significantly smaller than the respective values of the right-hand side of (48), then uk becomes an accurate approximation of the ideal gradient descent trajectory u(τ) at fictitious time τ := s0 + . . . + sk−1 . A combination of (46) and (48) justifies the following heuristic rule for choosing the search horizon at the current point uk ∈ U0 : ! kgk2 . (49) hk := min hmax , 2 |Dg E | u=uk Here, hmax is a given positive threshold which becomes active, for example, if Dg2 E vanishes. The stepsize selection algorithm considered below, replaces the minimization problem in (41) with a different procedure which involves a finite subset of values of s from a geometric progression sk,` := hk f ` , ` = 0, 1, 2, . . . (50) 7 whose initial value hk is given by (49). The common ratio f ∈ (0, 1) is a parameter of the algorithm which affects how “densely” the progression fills the interval [0, hk ]. Now, the adaptive stepsize selection algorithm proceeds as follows: sk := sk, j , (51) where the jth element of the geometric progression (50) is chosen according to the Armijo rule [4] with a parameter σ ∈ (0, 1): (52) j := min ` > 0 : E (uk ) − E (uk − sk,` g(uk )) > σ sk,` kg(uk )k2 . Here, the subset of indices ` is nonempty since σ < 1 and lim`→+∞ sk,` = 0. The inequality in (52) is important in the convergence analysis of the gradient descent algorithm. In particular, the condition σ > 0 ensures that E (uk ) is strictly decreasing until uk achieves a stationary point of the LQG cost. Such a point is a stable equilibrium of the gradient descent only if it delivers a local minimum to the LQG cost. C. Termination condition Since the gradient descent sequence uk in (40) can converge to a local minimum of the LQG cost only asymptotically, as k → +∞, the algorithm is equipped with a termination condition (for stopping the iterations) which reflects the proximity to the limit point. More precisely, we use the following termination condition which employs the relative smallness of the gradient as specified by a dimensionless parameter ε > 0: sk kg(uk )k 6 εkuk k. (53) IX. CONVERGENCE OF THE GRADIENT DESCENT ALGORITHM As the proposed algorithm is based on the classical gradient descent approach, its convergence analysis follows a similar reasoning, which we provide below for completeness. Theorem 1: Suppose (uk )k>0 is the gradient descent sequence in (40) with the stepsize selection described by (49)–(52). Then every limit point u∗ ∈ U0 of this sequence is a stationary point of the LQG cost E , that is, g(u∗ ) = 0. Proof: Since the sequence E (uk ) > 0 is nonincreasing, it has a finite limit. Therefore, (40), (51) and the Armijo rule in (52) imply that σ sk kg(uk )k2 6 E (uk ) − E (uk+1 ) → 0 as k → +∞. Hence, in view of σ > 0, it follows that lim sk kg(uk )k2 = 0. (54) k→+∞ Now, assume that the gradient descent sequence uk has a limit point u∗ := limK 3k→+∞ uk ∈ U0 such that g(u∗ ) 6= 0, where K := {0 6 k1 < k2 < . . .} is an infinite subset of nonnegative integers which specifies the respective subsequence of uk . Then the analyticity of the LQG cost on the open set U0 implies that lim K 3k→+∞ g(uk ) = g(u∗ ) 6= 0, lim K 3k→+∞ hk = min hmax , (55) kg(u∗ )k2 |Dg2 E (u∗ )| ! > 0, (56) where use is made of (49). Note that, if Dg2 E (u∗ ) = 0, the limit in (56) is equal to hmax > 0. A combination of (55) with (54) implies that lim sk = 0. (57) K 3k→+∞ In turn, by combining (57) with (56) and recalling (51) and the condition 0 < f < 1, it follows that the indices j p := sk ln hk −ln sk log f h p = −p ln f p of the elements of the geometric progression in (50), which correspond to k p ∈ K , diverge to infinity kp as p → +∞, and hence, j p > 1 for all sufficiently large p. For all such p, the stepsize candidates sk p , j p −1 = the Armijo selection rule in (52), that is, sk sk E (uk ) − E uk − g(uk ) < σ kg(uk )k2 f f sk p f do not pass (58) for all sufficiently large k ∈ K . Upon multiplying both parts of (58) by sf and taking the limit, this inequality leads to k f sk kg(u∗ )k2 = lim E (uk ) − E uk − g(uk ) K 3k→+∞ sk f 2 6 σ lim kg(uk )k = σ kg(u∗ )k2 , (59) K 3k→+∞ where use is also made of (55) and (57). However, since σ < 1, the inequality in (59) contradicts the assumption that g(u∗ ) 6= 0. This contradiction shows that any limit point u∗ ∈ U0 of the gradient descent sequence satisfies g(u∗ ) = 0. 8 Note that the CQLQG control problem inherits a special type of symmetry from the LQG cost E which is invariant under symplectic similarity transformations of the controller variables ξ 7→ Σξ , with Σ ∈ Rn×n satisfying ΣΘ2 ΣT = Θ2 and thus preserving the CCR matrix Θ2 (see, for example, [24], [25]). Hence, the stationary points of the LQG cost are not isolated, which complicates the convergence rate analysis for the proposed gradient descent algorithm. This issue is beyond the scope of the present study and will be addressed elsewhere by using more advanced analytic tools (such as in [2] and related references). X. A NUMERICAL EXAMPLE OF OPTIMAL CQLQG CONTROLLER DESIGN The gradient descent algorithm of Section VIII was tested to find a locally optimal solution of the CQLQG control problem for a PR quantum plant in (1) with dimensions n = m2 = p1 = p2 = r = 2, m1 = 4 and randomly generated state-space matrices A, B, C, E, satisfying the PR conditions (19), (22), and the weighting matrices F, G in (8): 0.9534 −1.1165 −1.7174 −0.2189 1.9180 0.5636 A= , B= , 0.4193 1.8821 −0.6815 1.3570 0.2985 −0.3679 −1.3570 −0.2189 1 0 0 0 −0.3238 0.2779 C= , D= , E= , −0.6815 1.7174 0 1 0 0 −1.1693 −0.5966 −0.8290 −0.9665 −0.2324 −0.1608 1 0 F= , G= , d= . −1.8655 −0.0357 −0.5822 −1.0961 0 1 This plant is unstable (the eigenvalues of the matrix A are 1.4177±0.5025i). The algorithm was run with parameters hmax = 1, f = 0.5, σ = 0.9, ε = 10−6 in (49)–(53) for 10 randomly generated stabilizing PR controllers as initial points. Starting from these points, it has taken 307 to 2318 steps for the algorithm to reach the fulfillment of the termination condition, with the average number of iterations being 1075. The local minimum value of the LQG cost is Emin = 12.1026 and is achieved at the following controller parameters: −0.5611 −1.5567 1.8111 0.7201 −0.1250 4.9673 R= , b= , e= . −1.5567 1.8283 −1.4979 −3.9696 −4.4929 −1.3387 The values of the LQG cost E (uk ) for the gradient descent sequences uk are presented in Fig. 2 in the form of semilogarithmic graphs of EE(uk ) − 1. These graphs are in qualitative agreement with the relatively slow linear convergence rate, min 2 10 Magnified Region 2 10 0 Relative Deviation 10 0 10 -2 10 -2 10 -4 0 10 20 40 60 -6 10 -8 10 0 Fig. 2. The relative deviations E (uk ) Emin 500 1000 1500 Number of Steps 2000 − 1 from the minimum value of the LQG cost on a logarithmic scale versus the number of steps k. typical for gradient descent methods. However, they also show that the proposed algorithm is fairly reliable, being able to cope with poor initial approximations where the LQG cost exceeds the minimum value by an order of magnitude. XI. CONCLUSION A gradient descent algorithm has been developed for the numerical solution of the optimal CQLQG controller design problem, and its convergence has been investigated. The algorithm has been tested and it appears to be fairly reliable in a numerical example with randomly generated stabilizing PR controllers as initial points. The lack of a more systematic method for initialization and a relatively slow convergence rate are the main shortcomings of the algorithm. These issues are a subject for future research and will be tackled in subsequent publications. 9 R EFERENCES [1] P.-A.Absil, K.Kurdyka, On the stable equilibrium points of gradient systems, Sys. Contr. Lett., vol. 55, no. 7, 2006, pp. 573–577. [2] P.-A.Absil, R.Mahony, and R.Sepulchre, Optimization algorithms on matrix manifolds, Princeton University Press, 2009. [3] D.S.Bernstein, and W.M.Haddad, LQG control with an H ∞ performance bound: a Riccati equation approach, IEEE Trans. Automat. Contr., vol. 34, no. 3, 1989, pp. 293–305. [4] D.P.Bertsekas, Nonlinear Programming, Athena Scientific, Belmont, 1999. [5] S.C.Edwards, and V.P.Belavkin, Optimal quantum filtering and quantum feedback control, arXiv:quant-ph/0506018v2, August 1, 2005. [6] C.W.Gardiner, and P.Zoller, Quantum Noise, Springer, Berlin, 2004. [7] A.S.Holevo, Statistical Structure of Quantum Theory, Springer, Berlin, 2001. [8] L.H¨ormander, An Introduction to Complex Analysis in Several Variables, North-Holland, New York, 1990. [9] M.R.James, H.I.Nurdin, and I.R.Petersen, H ∞ control of linear quantum stochastic systems, IEEE Trans. Automat. Contr., vol. 53, no. 8, 2008, pp. 1787–1803. [10] H.Kwakernaak, and R.Sivan, Linear Optimal Control Systems, Wiley, New York, 1972. [11] S.Lloyd, Coherent quantum feedback. Phys. Rev. A, vol. 62, no. 2, 2000, 022108. [12] L.A.Lusternik, and V.I.Sobolev, Elements of Functional Analysis, Gordon and Breach Science Publishers, New York, 1961. [13] H.Mabuchi, Coherent-feedback quantum control with a dynamic compensator. Phys. Rev. A, vol. 78, 2008, 032323. [14] J.R.Magnus, Linear Structures, Oxford University Press, New York, 1988. [15] E.Merzbacher, Quantum Mechanics, 3rd Ed., Wiley, New York, 1998. [16] H.I.Nurdin, M.R.James, and I.R.Petersen, Coherent quantum LQG control, Automatica, vol. 45, 2009, pp. 1837–1846. [17] K.R.Parthasarathy, An Introduction to Quantum Stochastic Calculus, Birkh¨auser, Basel, 1992. [18] I.R.Petersen, Quantum linear systems theory, Proc. MTNS, Budapest, Hungary, July 5–9, 2010, pp. 2173–2184. [19] J.J.Sakurai, Modern Quantum Mechanics, Addison-Wesley, Reading, Mass., 1994. [20] A.J.Shaiju, and I.R.Petersen, On the physical realizability of general linear quantum stochastic differential equations with complex coefficients, Proc. Joint 48th IEEE Conf. Decision Control & 28th Chinese Control Conf., Shanghai, P.R. China, December 16–18, 2009, pp. 1422–1427. [21] A.J.Shaiju, and I.R.Petersen, A frequency domain condition for the physical realizability of linear quantum systems, IEEE Trans. Automat. Contr., vol. 57, no. 8, 2012, 2033–2044. [22] R.E.Skelton, T.Iwasaki, and K.M.Grigoriadis, A Unified Algebraic Approach to Linear Control Design, Taylor & Francis, London, 1998. [23] I.G.Vladimirov, and I.R.Petersen, Hardy-Schatten norms of systems, output energy cumulants and linear quadro-quartic Gaussian control, Proc. MTNS, Budapest, Hungary, July 5–9, 2010, pp. 2383–2390. [24] I.G.Vladimirov, and I.R.Petersen, A dynamic programming approach to finite-horizon coherent quantum LQG control, Proc. AUCC, Melbourne, Australia, 10–11 November, 2011, pp. 357–362. [25] I.G.Vladimirov, and I.R.Petersen, A quasi-separation principle and Newton-like scheme for coherent quantum LQG control, Sys. Contr. Lett., vol. 62, no. 7, 2013, pp. 550–559. [26] I.G.Vladimirov, and I.R.Petersen, Coherent quantum filtering for physically realizable linear quantum plants, Proc. ECC, Z¨urich, Switzerland, 17–19 July 2013, pp. 2717–2723. A PPENDIX A. Computing the Fr´echet derivative of the LQG cost From the ALEs (18), (33) and the properties of the Frobenius inner product of matrices, it follows that the LQG cost E in (27) is representable as 1 1 E (u) = hC T C , Pi = hQ, BB T i = −hH, A i, (A1) 2 2 where H is the Hankelian given by (34). In what follows, δ (·) denotes the first variation, and δX (·) is the first variation with respect to an independent matrix-valued variable X. Since the matrix R influences the LQG cost E in (A1) only through the controller matrix a in (24) which enters the matrix A in (11), then 1 1 δR E = hC T C , δR Pi = − hA T Q + QA , δR Pi 2 2 1 1 = − hQ, A δR P + (δR P)A T i = hQ, (δR A )P + PδR A T i 2 2 =hH, δR A i = hH22 , δR ai = hH22 , 2Θ2 δ Ri. (A2) Now, since R ∈ Sn , and the subspaces Sn and An are orthogonal, then the first variation of the LQG cost in (A2) takes the form δR E = 2hΘT2 H22 , δ Ri = −2hsym(Θ2 H22 ), δ Ri, which, by the definition of the Fr´echet derivative, establishes (35). By 10 a similar reasoning, the first variations of the LQG cost with respect to b and e are as follows: δb E =hH, δb A i + hQB, δ Bi + hC P, δb C i =hH22 , δb ai + hE T H12 , δb ci + h(QB)22 , δ bi + hGT FP12 + GT GcP22 , δb ci =h − asym(H22 Θ−1 2 )bJ2 + Q21 Ed + Q22 b T T T T − Θ−1 2 (H12 E + P21 F G + P22 c G G)dJ2 , δ bi, 1 δe E =hH21 , δ e Ci + hH22 , δe ai + hQ, δe (BB T )i 2 =hH21CT , δ ei + hH22 Θ−1 , asym(δ eDJ1 DT eT )i + hQB, δe Bi =hH21CT − asym(H22 Θ−1 )eDJ1 DT + (Q21 B + Q22 eD)DT , δ ei, which leads to (36) and (37) in view of (26) and (38). B. Computing the second-order Gˆateaux derivative of the LQG cost The first-order Gˆateaux derivative of the LQG cost E along the gradient g in (31) is expressed in terms of the first-order Fr´echet derivatives from (35)–(37) by using (45), provided u ∈ U0 . Hence, the second-order Gˆateaux derivative along the gradient can be computed as Dg2 E = Dg (kgk2 ) = 2hDg g, gi =2(hDg ∂R E , ∂R E i + hDg ∂b E , ∂b E i + hDg ∂e E , ∂e E i), (B1) where Dg ∂R E = − 2sym(Θ2 Dg H22 ), (B2) Dg ∂b E =Dg Q21 Ed + Dg (Q22 b) − Dg (ψb)J2 − Dg χdJ2 , Dg ∂e E =Dg H21C + Dg Q21 BD T (B3) T + Dg (Q22 e) − Dg (ψe)DJ1 DT . (B4) The second-order Gˆateaux derivative in (B1) can now be computed by using (B2)–(B4), the Leibniz product rule and the first-order Gˆateaux derivatives of P, Q. By differentiating both sides of the ALE (18) and its dual (33), it follows that the matrices Dg P and Dg Q are unique solutions of the ALEs A Dg P + Dg PA T + 2sym(Dg A P + Dg BB T ) = 0, A T Dg Q + Dg QA + 2sym(A T Dg Q + C T Dg C ) = 0. Here, in view of (11), (24), (25) and the relation Dg u = g, 0 −EdJ2 ∂b E T Θ−1 2 Dg A = , ∂e E C Dg a 0 0 , Dg B = ∂e E D ∂b E Dg C = 0 −GdJ2 ∂b E T Θ−1 2 , with Dg a = 2Θ2 ∂R E − asym(∂e E DJ1 DT eT + ∂b E J2 bT )Θ−1 2 . 11

© Copyright 2017 ExploreDoc