Lyapunov Control of High-Dimensional Closed Quantum Systems Based on Particle Swarm Optimization

For high-dimensional closed quantum systems, this paper proposes a novel quantum Lyapunov control scheme based on the particle swarm optimization algorithm and achieves a high-probability population transfer of the system to a non-isolated target eigenstate in the LaSalle invariant set under usual smooth Lyapunov control laws. Via a quadratic Lyapunov function with unknown parameters, a control law with the unknown parameters is designed; based on the LaSalle invariance principle and the energy-level connectivity graph, the stability of the system is analyzed; by using the particle swarm optimization algorithm, a set of optimal parameters is obtained to achieve the control goal. In particular, we propose a path planning method based on the energy-level connectivity graph to determine the initial values of the unknown parameters, which is such that the optimization algorithm can efficiently and conveniently find a set of optimal solutions of the unknown parameters. Numerical simulation experiments on a five-level quantum system and a three-qubit system demonstrate that the proposed Lyapunov control scheme based on the particle swarm optimization algorithm has a good control effect.


I. INTRODUCTION
Since Tarn, Rabitz and Ong et al. started studying quantum control theory in 1980s [1]- [3], research on quantum control theory has made great progress and breakthroughs and quantum control has been successfully applied in quantum optics, atomic physics, quantum chemistry, quantum computing, and quantum communication [4]- [8], e.g., the preparation and manipulation of pure states that act as information carriers in the process of quantum computing, the manipulation of single atoms in the field of atomic physics, the preparation of compressed coherent states in the field of quantum optics, and the preparation and manipulation of entangled states in quantum communication. It should be addressed that such success of quantum control partly attributes to the extension of classical control theory to microscopic quantum systems, e.g., quantum optimal control [9], [10], quantum Lyapunov The associate editor coordinating the review of this manuscript and approving it for publication was Sotirios Goudos .
In the Lyapunov control methods, an appropriate Lyapunov function is first selected, and then a control law is designed by guaranteeing the monotone decrement of the Lyapunov function along the system trajectory and is exerted on the controlled quantum system. In particular, for closed quantum systems, the corresponding Lyapunov control laws are usually designed in a closed-loop way and implemented in an open-loop way. Dependent on different Lyapunov functions, two classes of Lyapunov control methods have appeared for closed systems, i.e., Lyapunov control methods based on deterministic Lyapunov functions [11], [24]- [27] and those based on Lyapunov functions with unknown VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ parameters [12], [28], [31]. The control laws obtained in these methods are usually smooth control laws and their effectiveness depends on the isolation of the target state in the LaSalle invariant set. When the target state is isolated in the LaSalle invariant set, it is easy to achieve the convergence of the system to the target state. When the target state is not isolated, it is difficult to achieve convergence to the target state since the designed control law cannot distinguish the target state from other states in the invariant set.
To solve the problem of non-isolated target states, different approaches have been proposed. For example, [11] proposed an asymptotic approximation method by using quantum adiabatic theorem and tracking a reference trajectory. Reference [24] proposed an approximation method based on an implicit Lyapunov function. Reference [29] proposed a switching Lyapunov control law determined by a switching sequence and realized approximate convergence to the target state. Via a quadratic Lyapunov function with unknown parameters, a control law with the unknown parameters was designed in [12], where the unknown parameters was roughly constructed by analyzing the largest invariant set via the LaSalle principle and therefore an approximate transfer to the non-isolated target eigenstate was achieved. However, in order to realize a high-probability population transfer to the target state, it is necessary to further determine specific values of the unknown parameters by numerical simulations and repetitious trials. For high-dimensional quantum systems, [30] proposed a multilayer Lyapunov control scheme based on the energy-level connectivity graph by analyzing the frequency selectivity characteristic of the Lyapunov control law and achieved a layer-by-layer high-probability population transfer to the target eigenstate. However, the amplitude of the control law is required to be sufficiently small to induce the frequency selectivity of population transfer. It should be noticed that small amplitude control fields and the multilayer scheme itself greatly lengthen the period of state transfer, which will reduce the robustness of quantum systems to uncertainties such as decoherence and perturbations caused by the environment. On the other hand, [13] first discussed the problem of rapid Lyapunov control of quantum systems, and proposed a standard bang-bang Lyapunov control law (the control function is a symbolic function) under the constraint of strength type. However, the discontinuous control law cannot guarantee a high-probability population transfer to the target state. Reference [28] proposed two rapidly convergent switching Lyapunov control laws for two-level quantum systems and two approximate bang-bang Lyapunov control laws for general quantum systems. Unfortunately, they still cannot induce high-probability population transfer to a non-isolated target eigenstate for high-dimensional quantum systems.
This paper further considers the problem of highprobability population transfer of high-dimensional closed quantum systems to non-isolated target eigenstates. Different from the multilayer transfer in [30], we hope to achieve a one-time high-probability population transfer to the target state, which can speed up the control process. To this end, we choose the Lyapunov function with unknown parameters to design a control law, and use an appropriate optimization algorithm to determine those parameters. Since the particle swarm optimization (PSO) algorithm has the advantages of fast search speed, being easily implemented, fewer adjustable parameters, and possessing memorability in the search process, we select the PSO algorithm to optimize the unknown parameters in the control law. The PSO algorithm was first proposed by Kenney and Eberhart in 1995 [35] and is an iteration-based optimization tool. At present, this algorithm has been improved in the aspects such as parameter setting [36], particle diversity [37], and population structure [38]. In this paper, we will utilize a standard model for the PSO algorithm to search for a set of optimal parameters and achieve the transfer of the quantum system to the target state. To improve the efficiency of the algorithm, this paper also will propose a path planning method based on the energylevel connectivity graph to obtain the good initial values of the unknown parameters to be optimized.
The main contributions of this paper are as follows. 1) We propose a novel quantum Lyapunov control scheme based on particle swarm optimization for high-dimensional closed quantum systems by combining Lyapunov control and optimization calculation and efficiently solve the problem of high-probability population transfer of the system to a non-isolated target eigenstate in the LaSalle invariant set.
2) Based on the energy-level connectivity graph for the quantum system, we obtain the inequality relationship satisfied by the unknown parameters in the control law and therefore generate good initial values of the parameters (i.e., initial values of particles) conveniently, which greatly improves the efficiency of the PSO algorithm's finding a set of optimal solutions of the unknown parameters. It should be pointed out that the method of determining appropriate initial values of the unknown parameters via the energy-level connectivity graph has not yet been reported in the existing references related to quantum optimal control. 3) Under an inequality constraint, we present an initialization method for the particle swarm and a search algorithm for the optimal solutions of the unknown parameters. 4) The proposed control scheme is verified on a five-level quantum system and a three-qubit (8-dimensional) system and the desired eigenstates are prepared successfully.
The rest part of this paper is organized as follows. In Section II, we present the quantum system model under consideration, describe the control task of this paper, and design a control law with the unknown parameters. Also, we will summarize the stability results based on the energylevel connectivity graph and the framework of the basic PSO algorithm. In Section III, we determine the unknown parameters via the PSO algorithm to achieve high-probability population transfer to the target state, where a path planning method based on the energy-level connectivity graph will be proposed to obtain the good initial values of the unknown parameters. We perform numerical simulation experiments on a five-level quantum system and a three-qubit system to verify the control effect of the proposed control scheme in Section IV. In Section V, we present the conclusion of this paper.

A. CONTROL PROBLEM
In this paper, we consider N -dimensional closed quantum systems. Assume that the quantum system under consideration is controllable and its state evolution is described by the following Schrödinger equation: where |ψ(t) is a wave function that describes the state of the N -dimensional quantum system and can be represented by a complex vector in C N ; ι = √ −1 is the imaginary unit;h is the reduced Planck constant and set to be 1 in this paper as usual; u κ (κ = 1, 2, . . . , r) are external real-valued control fields; H 0 is the internal Hamiltonian of the system; H k is the control Hamiltonian that describes the interaction between the system and the control field u κ . In the energy representation, H 0 can be expressed by a diagonal matrix, i.e., H 0 = diag[λ 1 , λ 2 , . . . , λ N ]. The transition frequency between two energy levels is defined as ω ab λ a −λ b (a, b ∈ {1, 2, . . . , N }). For convenience, we denote the jth eigenstate of H 0 associated with the eigenvalue λ j as |λ j [0, . . . , 0, 1, 0, · · · , 0] T with the jth element being 1 and the other elements being 0.
The control task of this paper is to steer the N -dimensional quantum system in (1) to an eigenstate of H 0 , |λ f , f ∈ {1, 2, . . . , N }. We assume that the system (1) satisfies the following conditions: Equations (2) and (3) mean that H 0 is non-degenerate (i.e., the diagonal elements of H 0 are mutually distinct) and the transition frequencies between different eigenstates are mutually distinct. For convenience, we only consider the case of one control Hamiltonian hereafter, i.e., κ = 1. The generalization to multiple control Hamiltonians is straightforward.

B. CONTROLLER DESIGN
Let us consider the following Lyapunov function: where P is a Hermitian operator to be constructed. The first-order time derivative of (4) along the trajectory of the system (1) can be calculated aṡ According to the Lyapunov stability theory, to guaranteė V ≤ 0 in (5), we let [H 0 , P] = 0 (6) and design the control law as where K 1 is a positive constant. Equation (6) implies that P is a diagonal matrix and the matrices H 0 and P have the same normalized eigenvectors. Thus, the matrix P can be written as where P m is the mth eigenvalue of P, and |λ m is the normalized eigenvector of H 0 corresponding to λ m . The first equal sign of (8) gives a corresponding relationship between the eigenvalues P m (m = 1, . . . , N ) of P and the eigenvectors |λ m (m = 1, . . . , N ) of H 0 . In each eigenstate |λ m , the value of the Lyapunov function is equal to P m , i.e., V (|λ m ) = P m . Here, the diagonal elements P 1 , P 2 , . . . , P N of P are the unknown parameters that need to be constructed. As for the relationship between the extrema of the Lyapunov function (4) and the matrix P, [31] proved the following lemma: Lemma 1: The set of critical points of the Lyapunov function (4) defined on the unit complex sphere ψ|ψ = 1 is given by the normalized eigenvectors of P. The eigenvectors associated with the largest eigenvalue of P are the maximum points of V , the eigenvectors associated with the smallest eigenvalue are the minimum points of V and all other eigenvectors are saddle points.
According to Lemma 1, we design the diagonal element P f of P associated with the target eigenstate |λ f as a minimum eigenvalue so that the system can be steered to the target state as the Lyapunov function decreases to its minimum value P f .

C. STABILITY RESULTS BASED ON ENERGY-LEVEL CONNECTIVITY GRAPHS
The quantum system (1) under the action of the control law (7) is a nonlinear autonomous system. We can use the LaSalle invariance principle to analyze its stability. This principle ensures that the system (1) with the control law (7) necessarily converges to the largest invariant set E contained in M {|ψ :V (|ψ ) = 0}. By calculating the largest invariant set, [30] gave the following stability theorem.
The energy-level connectivity graph can be used to assist in expressing the largest invariant set that the system converges to. Let us denote the energy-level connectivity graph associated with the system (1) as G = (V, E), where the set of vertices V is defined as V {|λ 1 , |λ 2 , . . . , |λ N }, and the set of edges E as E {(|λ j , |λ l ) : j = l, (H 1 ) jl = 0}. Evidently, the energy-level connectivity graph is determined by the control Hamiltonian H 1 . The condition (H 1 ) jl = 0 means that the jth and lth eigenstates in the energy-level connectivity graph are directly coupled. We call the nonzero element (H 1 ) jl (j < l) the weight of the edge (|λ j , |λ l ) or the coupling coefficient of the eigenstates |λ j and |λ l . For simplicity in notation, we also use the number ''j'' to stand for the jth vertex (i.e., the eigenstate |λ j ) in the connectivity graph. Thus, via the energy-level connectivity graph, the largest invariant set in Theorem 1 can be further simplified in the following theorem.
Theorem 2 [30]: For the system (1) satisfying conditions (2) and (3) and with the control law (7), the largest invariant set E contained in M can be expressed as the union of all single eigenstates of the system (1) and all single invariant subsets generated by those eigenstates that are not directly connected in their energy-level connectivity graph or that are associated with the same degenerate eigenvalue of P.
According to Theorem 2, for any initial state, the system converges to a certain invariant subset of the largest invariant set E. Further, when the diagonal elements of P are mutually distinct, all invariant subsets of E are fully determined by the control Hamiltonian H 1 . If the target eigenstate is directly coupled with all other eigenstates in the connectivity graph, then it is easily known that the target state is isolated in the invariant set E. Otherwise, the target state is non-isolated in (a certain invariant subset of) E. When the target state |ψ f is isolated in E, we can design the diagonal elements of P to satisfy: where V (|ψ(0) ) represents the initial value of the Lyapunov function V . It is easily verified that (9) can always guarantees the exact convergence of the system to the target state. When the target state is non-isolated in (a certain invariant subset of) E, it is difficult to construct the relation similar to (9) in order to achieve convergence to the target state. In this case, it is necessary to introduce an appropriate optimization algorithm to obtain a set of optimal solutions of the diagonal elements of P and therefore achieve a high-probability population transfer of the system to the target state.

D. PARTICLE SWARM OPTIMIZATION
In this paper, we use the PSO algorithm to find the diagonal elements of P to achieve a high-probability population transfer to a non-isolated target eigenstate in E. In the PSO algorithm, each particle can be regarded as a flight individual in the search space and has a position vector and a velocity vector. The current position of each particle is a candidate solution of the optimization problem. We may calculate the fitness value of the current position via a fitness function to evaluate the quality of the current position. Generally speaking, the fitness function can be easily obtained from the objective function to be optimized. The flight process of a particle is the search process of this individual. The flight velocity of each particle represents the direction and increment of the particle's movement in the next iteration and can be dynamically adjusted according to the current optimal position of this particle and the current optimal position of the whole particle swarm. Here, the optimal position that a particle has achieved so far is called the current individual best position (or the current individual optimal solution), and the optimal individual best position in the whole particle swarm so far is called the current global best position (or the current global optimal solution). In what follows, we describe the update formulas of the position and velocity of each particle in the basic PSO algorithm [36], [39].
Suppose that there are m particles in a D-dimensional search space. We represent the position and velocity vectors of the ith particle as and respectively. Denote the current individual best position of the ith particle as and the current global best position of the whole particle swarm as p(g) = (p 1 (g), p 2 (g), . . . , p D (g)).
In the iteration process, the velocity and position of the ith particle are updated according to the following formulas: x k+1 where i = 1, . . . , m; d = 1, . . . , D; ω is an inertia coefficient; r 1 and r 2 are random numbers in [0, 1]; c 1 and c 2 are learning factors. The first term on the right-hand side of (14) is the inertia term, which indicates that the particle has a tendency to maintain its original velocity. The second term is local heuristic information, which indicates that the particle is approaching its own historical best position. The third term is global heuristic information and indicates that the particle is approaching the historical best position of the whole population group. The basic procedure of the PSO algorithm can be divided into the following steps.
Step 1. Randomly initialize the particle swarm, including the number of particles m and the initial values of position x(i) and velocity v(i) of each particle.
Step 2. Select an appropriate fitness function and calculate the fitness value f (x(i)) of each particle.
Step 3. For each particle x(i), compare the fitness values f (x(i)) and f (p(i)). If f (x(i)) > f (p(i)), then update the current f (p(i)) with f (x(i)) and update p(i) with x(i).
Step 4. Similarly, for each particle x(i), compare the fitness values f (x(i)) and f (p(g)). If f (x(i)) > f (p(g)), then update the current f (p(g)) with f (x(i)) and update p(g) with x(i).
Step 5. Update the velocity and position of each particle according to (14) and (15). Step 6. Judge whether the termination condition is satisfied.
If satisfied, then the iteration ends and the optimal solution is output; otherwise, jump to Step 2. Here, the termination condition is usually a preset maximum number of iterations or a desired fitness value.

III. DETERMINATION OF UNKNOWN PARAMETERS
As stated in section II-C, for an isolated target state in the largest invariant set E, we may construct the unknown parameters in P according to (9) to achieve convergence to the target state. But for a non-isolated target state in E, it is difficult to find the diagonal elements of P. In this section, we solve this problem via the PSO algorithm in section II-D to achieve a high-probability population transfer to the target state. For the PSO algorithm to efficiently find the diagonal elements of P, we first propose a path planning method based on the energy-level connectivity graph to determine the initial values of the parameters P 1 , P 2 , . . . , P N , and then use the PSO algorithm to search for and find the optimal solutions of these parameters.

A. CHOICE OF INITIAL VALUES
The energy-level connectivity graph represents the connection relationship between energy eigenstates. According to the energy-level connectivity graph, we have the following observations: Observation 1: The population transfer of the system to the target state can be regarded as a process where the populations in all energy eigenstates transfer to the target state along various paths in the energy-level connectivity graph.
Observation 2: A path whose each edge has a larger weight modulus is easier to induce the population transfer of the system to the target state.
Based on observations 1 and 2, we may plan one path/some paths that covers/cover all eigenstates such that the system population transfers to the target state along the planned path/paths. To this end, we present several concepts related to paths. In the energy-level connectivity graph associated with a quantum system, the number of eigenstates covered by a path is called the covering number of this path. If a path (or a combination of several branch paths) to an eigenstate covers all eigenstates and has the least covering number, then this path (or this combination of the several branch paths) is called a complete path to this eigenstate. The product of weights of all edges of a path is called the weight of this path. The modulus of the weight of a path is called the strength of this path, denoted by S. The smallest value of strengths of all branch paths of a complete path is called the strength of this complete path, denoted by S c .
Thus, based on the connectivity graph, we can easily choose a complete high-strength path as a superior path. Then, according to the correspondence between all eigenstates and the diagonal elements of P, we set the diagonal elements of P in a decreasing order along the superior path from far to near and therefore obtain a set of initial values of the diagonal elements of P. In this paper, this method of obtaining the initial values of the unknown parameters from the energy-level connectivity graph is called the path planning method.
Example 1: Assume that a quantum system is controlled by only one external control field and its internal and control Hamiltonians are given as and respectively. Further, the matrix P is chosen as a nondegenerate matrix. The target eigenstate is given as |λ 3 . According to Theorem 2, we can easily write the largest invariant set E as E = span{|λ 1 , |λ 3 } ∪ span{|λ 2 , |λ 3 } ∪ span{|λ 3 , |λ 4 }∪span{|λ 4 , |λ 5 }, where span{|λ j , |λ l } represents the subspace spanned by |λ j and |λ l . Its energy-level connectivity graph is drawn in FIGURE 1, where the number beside each edge is the weight of the corresponding edge. It is clear that the target state is non-isolated in the largest invariant set. It can be seen from FIGURE 1 that there are six complete paths to |λ 3 , i.e., Their strengths (i.e., moduli of weights of the six paths) are S c 1 = 0.5, S c 2 = 0.2, S c 3 = 0.05, S c 4 = 0.01, S c 5 = 0.1, and S c 6 = 0.1, respectively. Therefore, we choose Path 1 as a superior path and set the diagonal elements of P to satisfy P 1 > P 4 > P 2 > P 5 > P 3 . According to such an inequality, the initial values of P 1 , P 2 , . . . , P N can be obtained.
Usually, we are inclined to use a complete path with only one branch since it provides a more definite inequality constraint for the unknown parameters.

B. PARAMETER OPTIMIZATION BASED ON THE PSO ALGORITHM
According to the path planning method in section III-A, we can only obtain the initial values of the unknown parameters (the diagonal elements of P). In order to achieve a highprobability population transfer of the system to the target state, we further use the PSO algorithm to optimize these unknown parameters in this section.
In the PSO algorithm, we need to choose an appropriate fitness function according to the optimization objective. Since our control objective is to achieve a high-probability population transfer of the high-dimensional quantum system to an eigenstate of its internal Hamiltonian, it is natural to select the transition probability from the system state to the target state at the end of control as our fitness function, i.e., where |ψ(T ) denotes the system state at a prespecified terminate moment T and |λ f is the desired target state. We hope that the transition probability in (16) can reach or approach the maximum value 1.
In order to use the PSO algorithm to search for a set of optimal values of the parameters P 1 , P 2 , . . . , P N , we construct the position vector of an N -dimensional particle as (P 1 , P 2 , . . . , P N ) (x 1 (i), x 2 (i), . . . , x N (i)) = x(i), (17) where x(i) (i = 1, 2, . . . , m) represents the position vector of the ith particle and m is the number of particles; P j is the jth diagonal element of P.
Suppose that the inequality relationship among P 1 , P 2 , . . ., P N obtained from the path planning method is where v j ∈ {1, 2, . . . , N }, and v 1 = v 2 = · · · = v N . Then, all components of each particle's position vector should satisfy the inequality (18) in the initialization of the particle swarm and the iteration process of the PSO algorithm. Letting with j = 1, . . . , N and i = 1, . . . , m, we can construct an N -dimensional search space. Each particle moves in this space and its each position component takes values in [a, b].
In the following, we present a particle swarm initialization method under the constraint (18). For the ith (i = 1, 2, . . . , m) particle, we select the range of x ν 1 (i) as [c, b], where c ∈ (a, b) is used to restrict the lower bound of x ν 1 (i). Then, under the constraint condition (18), we can initialize every component of the ith particle as follows: Step 1. According to x ν 1 (i) = c + ν 1 rand(1), randomly generate a number in [c, b] as an initial value of x ν 1 (i), and denoted as z ν 1 (i), where rand(1) represents a random number in [0, 1] and ν 1 is a positive number not greater than b − c.
After completing the random initialization of the particle swarm and the selection of the fitness function, we can search for the optimal values of the unknown parameters via the PSO algorithm in section II-D. Since the inequality in (18) is a constraint condition in the iteration process of the algorithm, we need to verify whether (18) is satisfied for every particle in each generation population. If a particle in a certain iteration does not meet the constraint condition in (18), then we set the current position of the particle as its historical optimal value. Thus, we can write a complete algorithm procedure of this paper as follows.
Step 1. Initialize the particle swarm according to the initialization method in this section.
Step 3. For each particle x(i), compare the fitness value f (x(i)) and the previous individual extremum f (p(i)).
, then update f (p(i)) with f (x(i)) and update p(i) with x(i).
Step 4. For each particle x(i), compare the fitness value f (x(i)) and the previous global extremum f (p(g)).
If f (x(i)) > f (p(g)), then update f (p(g)) with f (x(i)) and update p(g) with x(i).
Step 5. Update the velocity and position of each particle according to (14) and (15). Step 6. Move each particle to the next position.
Step 7. Check if the current position of each particle is still inside the search space and the constraint condition (18) is satisfied. If the current position of a particle is outside the space or the condition (18) is not satisfied, then set the current position of this particle to its best value within the space.
Step 8. Judge whether the termination condition is satisfied. If satisfied, then the iteration ends and the optimal solution is output; otherwise, jump to Step 2.

IV. NUMERICAL SIMULATIONS
In this section, we perform simulation experiments on a five-level quantum system and a three-qubit system to verify the effectiveness of the Lyapunov control scheme with particle swarm optimization proposed in this paper, and at the same time prepare desired target states.

A. A FIVE-LEVEL QUANTUM SYSTEM
Consider a five-level quantum system, whose internal and control Hamiltonians are respectively. Assume that the initial state of the system is given as |ψ(0) = 1 |λ 4 and the target eigenstate is |λ 3 . Based on the control Hamiltonian H 1 , we can draw the energy-level connectivity graph as shown in FIGURE 2. The largest invariant set of this system and the complete paths to the target state are the same as that in Example 1. With the weights in FIGURE 2, it is easily known that the strength of each complete path is 1. Therefore, we can choose an arbitrary complete path to the target state as a superior path. Here, we choose the path 4 → 2 → 1 → 5 → 3. Thus, according to the path planning method, the diagonal elements of P satisfy We first verify the effect of the path planning method. Under the inequality constraint in (20), we choose the matrix P as P = diag [1.4, 1.8, 0.6, 2.0, 1.3] and meanwhile take the coefficient K 1 = 0.07 in the control field (7). The simulation results are shown in FIGURE 3 and FIGURE 4. It can be seen from FIGURE 3 that the system trajectory finally converges to span{|λ 1 , |λ 3 } and the population in the target eigenstate (i.e., the transition probability from the system state to the target state) reaches around 0.73.
Next, we use the PSO algorithm to determine the matrix P. The number of particles is chosen as m = 10. By letting every position component of each particle takes values in [0.1, 2.5], a search space for the particle swarm is established. According to the initialization method in section III-B, we initialize the particle swarm by selecting c = 2 and letting x 4 (i) = 2 + 4 rand(1), (1) (16) is set as T = 700 a.u.. In order to compare with the path planning method, we choose K 1 = 0.15 in the control law (7) so that the maximum strength of the control law (7) can reach the maximum strength of the control field in FIGURE 4. Through the execution of the PSO algorithm in section III-B, an optimal matrix P is found, i.e., P = diag[1.7894, 1.8308, 1.5229, 2.0388, 1.6616]. With this P, VOLUME 8, 2020  The simulation shows that the population in the target state |λ 3 reaches 0.9999. It can be seen from FIGURE 3 and FIGURE 5 that the PSO algorithm further improves the control effect. Thus, the control scheme based on the PSO algorithm achieves a satisfactory high-probability population transfer of the system to the target eigenstate.

B. A THREE-QUBIT QUANTUM SYSTEM
We further consider a quantum system of three qubits with the internal and control Hamiltonians as We assume that the target state is |λ 1 and the initial state is given as |ψ (0)  With the control Hamiltonian H 1 , we can obtain the energy-level connectivity graph of this three-qubit quantum  system, which is shown in FIGURE 7. From FIGURE 7 and Theorem 2, it follows that the target state |λ 1 is nonisolated in the largest invariant set of the system. We choose the following complete path to the target state |λ 1 as our superior path: 3 → 7 → 5 → 6 → 8 → 4 → 2 → 1. According to the path planning method, we can obtain the inequality relationship between the diagonal elements of P as In order to use the PSO algorithm to determine the matrix P, we choose the number of particles as m = 15, the learning factors as c 1 = 2 and c 2 = 2, and the inertia coefficient as ω = 0.5. The search space of the particle swarm is defined by guaranteeing that every component of each particle's position vector takes values in [0.1, 2.5]. We initialize the particle swarm by following the initialization method in section III-B, where c = 2, x 3 (i) = 2 + 3 rand(1), FIGURE 9. Under the matrix P determined by the PSO algorithm, the evolution curve of the control field for the three-qubit quantum system. (1), . . ., x 1 (i) = x 2 (i) − 1 rand(1) (i = 1, . . . , 15), with 3 = 0.1, 7 = 0.5, 5 = 0.3, and Further, we set the control gain in (7) as K 1 = 0.15, the terminate moment in the fitness function (16) as T = 200 a.u., and the termination condition for the PSO algorithm as the maximum number of iterations 50. By running the PSO algorithm in section III-B, we find an optimal matrix that achieves a high-probability population transfer of the system to the target state, i.e., P = diag[0.7071, 1.2599, 2.0866, 1.3369, 1.9342, 1.7239, 1.9971, 1.6638]. With this optimal matrix P, the simulation results are shown in FIGURE 8 and FIGURE 9. According to FIGURE 8 and simulation data, the population in the target state (i.e., the transition probability from the system state to the target state) reaches 0.9985 at the terminate moment T = 200 a.u..

V. CONCLUSION
This paper has solved the problem of high-probability population transfer of high-dimensional closed quantum systems under Lyapunov control to a non-isolated target eigenstate in the largest invariant set by using the particle swarm optimization algorithm. We have designed a control law via a Lyapunov function with unknown parameters and proposed a path planning method based on the energy-level connectivity graph to determine the inequality relationship satisfied by those unknown parameters. According to the obtained inequality relation, we have generated the initial values of particles in the PSO algorithm, which enables a relatively good control effect. Accordingly, we have proposed a particle swarm initialization method and a PSO algorithm under the inequality constraint. These make the PSO algorithm more efficient to find a set of optimal solutions of the unknown parameters and therefore easier to achieve a satisfactory highprobability population transfer of the system to the target state. It should be addressed that the optimal solutions of the unknown parameters found by the PSO algorithm are not unique since the PSO algorithm is a random search algorithm in nature. We have done numerical simulation experiments on a five-level system and a three-qubit system. Simulation results show that the proposed path planning method does provide good initial values for the PSO algorithm and that the proposed Lyapunov control scheme based on particle swarm optimization has a good control performance. Therefore, such a new control scheme can be used to satisfactorily solve the equilibrium manifold problem caused by a non-isolated target state in the largest invariant set under usual Lyapunov control and achieve a satisfactory transfer of the system to the target state.
Finally, we would like to point out that this paper only use the standard PSO algorithm to search for a set of optimal solutions of the unknown parameters. Some improved PSO algorithms can be used to further improve the efficiency and performance of the algorithm. In the PSO algorithm, we select the transfer probability from the system state to the target state at a specified terminal moment as a fitness function to evaluate the quality of each particle. In fact, the terminal moment and the control energy can be optimized. In addition, one also can consider introducing rapid Lyapunov control into the control scheme in this paper to further improve the speed of state transfer.