Differential Evolution for Quantum Robust
Control: Algorithm, Applications and Experiments
Daoyi Dong, Xi Xing, Hailan Ma, Chunlin Chen, Zhixin Liu, Herschel Rabitz
Abstract—Robust control design for quantum systems has
been recognized as a key task in quantum information technology, molecular chemistry and atomic physics. In this paper,
an improved differential evolution algorithm of msMS DE is
proposed to search robust fields for various quantum control
problems. In msMS DE, multiple samples are used for fitness
evaluation and a mixed strategy is employed for mutation
operation. In particular, the msMS DE algorithm is applied to
the control problem of open inhomogeneous quantum ensembles and the consensus problem of a quantum network with
uncertainties. Numerical results are presented to demonstrate
the excellent performance of the improved DE algorithm for
these two classes of quantum robust control problems. Furthermore, msMS DE is experimentally implemented on femtosecond laser control systems to generate good signals of twophoton absorbtion and control fragmentation of halomethane
molecules CH2 BrI. Experimental results demonstrate excellent
performance of msMS DE in searching effective femtosecond
laser pulses for various tasks.

I. I NTRODUCTION
Controlling quantum systems has become a fundamental
task in promising quantum technology [1], and it is relevant
to many emerging areas such as atomic physics, molecular chemistry, and quantum information science [2], [3],
[4]. Some control methods, such as optimal control theory
[1], learning control algorithms [2] and Lyapunov control
approaches [5], have been developed for quantum systems.
Among these methods, learning control is a powerful approach for many complex quantum control tasks [6] and has
achieved great success in laser control of molecules since
the method was presented in the seminal paper [7]. Many
quantum learning control problems could be formulated as an
optimization problem and a learning algorithm is employed
to search for an optimal control field for achieving desired
performance. Gradient algorithms have been demonstrated
to be a good candidate for numerically finding an optimal
This work was supported by the National Natural Science Foundation
of China (No. 61374092), the Australian Research Council’s Discovery
Projects funding scheme under Project DP130101658, NSF (CHE-1464569)
and ARO (W911NF-16-1-0014).
D. Dong is with the School of Engineering and Information Technology,
University of New South Wales, Canberra, ACT 2600, Australia, and the
Department of Chemistry, Princeton University, Princeton, NJ 08544, USA.
(email: daoyidong@gmail.com).
X. Xing is with the Department of Chemistry, Princeton University,
Princeton, NJ 08544, USA. (email: xxing@princeton.edu).
H. Ma and C. Chen is with the Department of Control and Systems
Engineering, School of Management and Engineering, Nanjing University,
Nanjing 210093, China.(e-mail: clchen@nju.edu.cn).
Z.X. Liu is with the Key Laboratory of Systems and Control, Academy of
Mathematics and Systems Science, Chinese Academy of Sciences, Beijing
100190, China (e-mail: lzx@amss.ac.cn).
H. Rabitz is with the Department of Chemistry, Princeton University,
Princeton, NJ 08544, USA. (email: hrabitz@princeton.edu).

field due to their high efficiency [8]. In many practical
applications, the gradient information may not be easy to
obtain and some complex quantum control problems may
have local optimum. For these situations, it is necessary to
employ stochastic searching algorithms to find a good control
field. Genetic algorithm (GA) has been widely used in the
area of quantum control and has achieved success in closedloop learning control of molecular systems in the laboratory
[2]. In this paper, we focus on robust control problems of
quantum systems and explore the use of differential evolution
(DE) [9] for searching robust control fields.
DE is a competitive form of evolutionary computation
and has shown great performance for many complex optimization problems [10]. Recently, it has also been used for
solving quantum control problems [11], [12]. For example,
Zahedinejad et al. [11], [13] proposed a subspace-selective
self-adaptive differential evolution (SUSSADE) algorithm to
achieve a high-fidelity single-shot Toffoli gate and singleshot three-qubit gates. DE methods have been adopted to
fulfil a desired state transition by designing optimal control
fields for an open quantum ensemble [14]. Zahedinegad et al.
[15] investigated several promising evolution algorithms and
found DE outperformed GA and particle swarm optimization
(PSO) for hard quantum control. In this paper, we employ DE
algorithms to solve several classes of quantum robust control
problems. The robustness of quantum systems has been
recognized as an essential requirement for the development
of practical quantum technology, and robust control methods
could provide enhanced robustness performance for quantum
systems [16], [17]. In particular, we propose an msMS DE
(multiple-samples and Mixed-Strategy DE) algorithm where
a mixed strategy and an average performance with multiple
samples are employed. The msMS DE is used for three
classes of quantum robust control problems: control of inhomogeneous open quantum ensembles, consensus of quantum
networks with uncertainties and experimental fragmentation
control using femtosecond laser pulses [18].
A quantum ensemble consists of a large number of single
quantum systems (e.g., spin systems or molecules) and
every individual system is referred to as a member of the
ensemble [19]. In practical applications, the members of
a quantum ensemble could have variations in the parameters characterizing the system dynamics. Such a quantum
ensemble is called an inhomogeneous quantum ensemble
[19], [20]. Inhomogeneous quantum ensembles have wide
applications in fields ranging form long-distance quantum
communication to magnetic-resonance imaging [21], [22].
Hence, it is highly desirable to design control laws that can

steer an inhomogeneous ensemble from an initial state to a
desired state when variations exist in the system parameters.
Some results have been presented for controllability analysis
and control design of inhomogeneous quantum ensembles.
For example, Li and co-workers [20], [21] presented a
series of results to analyze controllability and design optimal control laws for inhomogeneous spin ensembles. A
Lyapunov control design approach has been proposed to
asymptotically stabilize a spin ensemble around a uniform
state of spin +1/2 or −1/2 [23]. Chen et al. [19] presented
a sampling-based learning control method to achieved high
fidelity control of inhomogeneous quantum ensembles. In
these results, decoherence and dissipation were usually not
considered. The existence of decoherence and dissipation
may irreversibly lead a quantum ensemble to becoming an
open system [24], and the manipulation of inhomogeneous
open quantum ensembles becomes more challenging than
that without considering decoherence. For an inhomogeneous
quantum ensemble, we cannot employ different control fields
to control individual members. A practical solution is to find
robust control fields that can drive all of the members in
the ensemble into a given target state. In this paper, we
employ an msMS DE algorithm to search such robust control
fields for inhomogeneous open quantum ensembles aiming at
achieving enhanced control performance.
Another problem under consideration is to drive a quantum
network into a consensus state even in the presence of uncertainties. Achieving consensus is one of primary objectives
in distributed control and coordination of classical (nonquantum) networked systems [25]. Consensus usually means
that all of the nodes in a network hold the same state. Recent
development in quantum technology has made it significant
and feasible to analyze quantum networks where each node
(agent) represents a quantum system such as a photon, an
electron, a spin system or a superconducting quantum bit
(qubit). Consensus of quantum networks may have potential
applications in promising quantum communication networks,
distributed quantum computation and one-way quantum computation [26], [27]. Since the nodes in a quantum network
are described by quantum mechanics, some unique characteristics such as quantum entanglement and measurement
backaction should be carefully considered, and the analysis
and control of quantum networks raise new challenges.
Some results have been presented for the consensus problem
of quantum networks. For example, Sepulchre et al. [28]
generalized consensus algorithms to noncommutative spaces
and analyzed the asymptotic convergence to the consensus
state of a fully mixed state. Ticozzi and co-workers [26], [29],
[30] presented a series of results on consensus of quantum
networks including several different definitions for quantum
consensus, quantum gossip algorithms, and quantum consensus results within a group-theoretic framework. Shi et
al. [27] presented a systematic investigation on consensus
of quantum networks with continuous-time dynamics within
the framework of graph theory. In this paper, we consider
the basic problem of finding a robust control law to steer a

quantum network to a reduced state consensus (as defined
in [26]) and would not consider the distributed solutions
to achieving quantum consensus. In particular, we employ
the proposed msMS DE for driving a superconducting qubit
network with uncertainties into a reduced state consensus.
Femtosecond (fs) (1fs = 10−15 second) laser [18] has
found wide applications in controlling the molecular dynamics because of its short pulse duration, which is comparable
to the time scales of electronic and nuclear motions of a
molecule. The temporal structures of a femtosecond pulse
could be manipulated by pulse shaping techniques [18],
which is typically achieved by modulating the phase and/or
amplitude of the laser frequency components, spatially separated in the Fourier plane, with a computer programmable
spatial light modulator (SLM) before recombination into a
“shaped pulse”. The advances in both fs lasers and pulse
shaping techniques have spurred recent interest in laserselective chemistry and quantum control using shaped fs
laser pulses. In quantum control experiments, a practical
approach is to use closed-loop learning control [2] to find
an optimal field that can steer the quantum system towards
the desired outcome, which has achieved a lot of successes.
An evolutionary algorithm (e.g., GA) is often employed to
assist the search of an optimal pulse. To the best of our
knowledge, there have not been any experimental results
reported where DE algorithms are used on femtosecond laser
control systems. There are few quantum control experiments
using femtosecond laser pulses that have investigated the
robustness subject to variations in the control. In this work,
we employ the msMS DE algorithm to an experimental
quantum control problem, where the goal is to identify a
robust solution (shaped fs laser pulse) that can maximize
the CH2 Br+ /CH2 I+ product ratio from the fragmentation
of CH2 BrI molecule in a time-of-flight mass spectrometry
(TOF-MS). The msMS DE algorithm is also used to identify
the transform limited (TL) pulse via optimizing the two
photon absorption (TPA) signal, which is carried out prior to
the fragmentation experiment.
The main contributions of this paper are summarized as
follows:
•

•

•

Motivated by solving three classes of quantum robust
control problems, an improved DE algorithm of msMS DE is proposed where an average fitness function of
multiple samples is utilized and a mixed strategy of
mutation is employed.
The control problem of inhomogeneous open quantum
ensembles is investigated and the msMS DE algorithm
is used to learn robust control fields for inhomogeneous
open quantum ensembles. Numerical results show that
the control fields learned by msMS DE usually have
improved robustness performance compared with those
learned by basic DE and GA.
The task of driving a quantum network with uncertainties to a consensus state is investigated and robust control fields can be found by employing the msMS DE algorithm. Numerical results demonstrate that msMS DE

•

has excellent performance in searching robust control
fields for achieving consensus of quantum networks.
Several experiments are implemented on femtosecond
laser control systems where the msMS DE algorithm
is employed to find effective femtosecond laser pulses
for generating excellent TPA signal and achieving good
fragmentation control of molecules. These experiments
present the first result of testing DE for femtosecond
laser quantum control systems as well as realize the
sampling-based learning control method [19], [31], [32]
in the area of quantum control.

The paper is organized as follows. Section II formulates
three classes of quantum robust control problems under
consideration. Section III briefly introduces the basic DE and
then presents a systematic description of the msMS DE algorithm. Numerical results on quantum ensemble control and
quantum network consensus are provided in Section IV. In
Section V, we present experimental results on femtosecond
laser control systems. Conclusions and discussions are given
in Section VI.
II. T HREE CLASSES

OF QUANTUM CONTROL PROBLEMS

For an inhomogeneous open quantum ensemble, the
Hamiltonian can be described in the form of
M

HΘ (t) = g0 (θ0 )H0 + ∑ g j (θ j )u j (t)H j .

Let Θ = (θ0 , θ1 , ..., θM ) and the functions g j (θ j ) ( j =
0, 1, . . . , M) characterize possible inhomogeneities. For example, g0 (θ0 ) corresponds to inhomogeneity in the free
Hamiltonian (e.g., due to chemical shift in NMR). g j (θ j )
( j = 1, ..., M) can characterize possible multiplicative noises
in the control fields or imprecise parameters in the dipole
approximation. We assume that g j (θ j ) ( j = 0, 1, ..., M) are
continuous functions of θ j and the parameters θ j could be
time-dependent and θ j ∈ [1 − E j , 1 + E j ]. For simplicity, we
assume that g j (θ j ) = θ j , the nominal value of θ j is 1 and
E0 = ... = E j = ... = EM = E.
For an open quantum system in (1), we may define
a coherent vector as y := (tr(Ul ρ ), tr(U2 ρ ), ..., tr(Um ρ ))⊤ ,
where iU1 , iU2 , ... iUm (m = n2 − 1) are orthogonal generators
of the special unitary group SU(n) with degree n . Its density
operator can be written as:

ρ=

A. Control of inhomogeneous open quantum ensembles
We consider an inhomogeneous quantum ensemble where
the members of the ensemble could have variations in
the parameters that characterize the system dynamics. For
example, the spins of an ensemble in nuclear magneticresonance (NMR) experiments may encounter large dispersion in the strength of the applied radio frequency field (rf
inhomogeneity) and there also exist variations in the natural
frequencies of these spins (Larmor dispersion) [20]. Several
methods have been presented to design control laws for
inhomogeneous quantum ensembles when dissipation and
decoherence were not considered [19], [21].
For a practical quantum ensemble, it may be unavoidable
to interact with its environment and each member in the
quantum ensemble should be dealt as an open quantum
system. For an open quantum system, its state can be described by a positive Hermitian density operator ρ satisfying
tr(ρ ) = 1. Under the assumption of a short environmental
correlation time permitting the neglect of memory effects, a
Markovian master equation for ρ (t) can be used to describe
the dynamics of an open quantum system interacting with
its environment [24]. Markovian master equations in the
Lindblad form are described as [3], [33]
i
ρ̇ (t) = − [H(t), ρ (t)] + ∑ γk D[Lk ]ρ (t),
h̄
k

I 1 m
+ ∑ yl Ul .
n 2 l=1

(3)

Substituting (3) into (1), we can obtain the evolution of the
coherent vector y as:
M

ẏ = (LH0 + LD )y + ∑ u j LH j y + l0 ,

(4)

j=1

where the superoperators LH0 , LD , LH j ( j = 1, 2, ..., M) and
the term l0 are explained in detail in [34], [35]. We choose
the objective function to be maximized as follows [34]:
J(u) = 1 −

n
k y f − y(T ) k2 ,
8(n − 1)

(5)

where k x k2 = xT x is a vector norm and it is clear that J(u) ∈
[0, 1]. And y f and y(T ) are the target state and the final
state of the quantum system in terms of coherent vector,
respectively.
The control problem of an inhomogeneous open quantum
ensemble can be formulated as:
max J(u) := max E[JΘ (u)]
u

u

s.t.
(1)

√
where i = −1, H(t) is the system Hamiltonian, h̄ is reduced
Planck constant (we set h̄ = 1 in this paper), the non-negative
coefficient γk specify the relevant relaxation rates, Lk are
appropriate Lindblad operators and
1
1
D[Lk ]ρ = (Lk ρ L†k − L†k Lk ρ − ρ L†k Lk ).
2
2

(2)

j=1






ẏΘ (t) = (θ0 LH0 + LD + θ j

M

∑ u j (t)LH j )yΘ (t) + l0

j=1


yΘ (0) = y0 , t ∈ [0, T ]




θ j ∈ [1 − E, 1 + E], j = 0, 1, ...M

(6)

where E[JΘ (u)] denotes the average performance function
regarding parameter inhomogeneities Θ .

B. Consensus in quantum networks

The problem can be formulated as follows:

Achieving quantum consensus is a primary objective in the
investigation of quantum networks. Existing results presented
some distributed solutions to quantum consensus problems. For example, Mazzarella et al. [26] proposed a quantum gossip iteration algorithm where discrete-time quantum
swapping operations between arbitrary two nodes are used to
make a quantum network achieving consensus. A graphical
method has been developed in [27] to build the connection
between quantum consensus and its classical counterpart, and
asymptotical convergence results on achieving a consensus
state have been presented for a class of quantum networks
with continuous-time Markovian dynamics. Here, we would
not intend to develop a distributed algorithm for quantum
consensus. In contrast, we consider how to design a robust
control field to drive a quantum network from an initial state
into a consensus state with high fidelity when uncertainties or
inaccuracies may exist in the system dynamics. We consider
the type of Reduced State Consensus that was defined in
[26]. We denote H a Hilbert space, A ⊗ B returns the tensor
product of A and B, and |ai with Dirac representation is
a vector in H , i.e., |ai ∈ H [4]. In order to present the
definition of reduced state consensus, we need to use the
concept of partial trace defined as follows.
Definition 1 (Partial trace): [4] Let HA and HB be the
state spaces of two quantum systems A and B, respectively.
Their composite system is described as a density operator
ρ AB. The partial trace over system B, denoted as TrHB is
given in the following form
TrHB (|a1 iha2 | ⊗ |b1ihb2 |) = |a1 iha2 |Tr(|b1 ihb2 |),

(7)

where the vectors |a1 i, |a2 i ∈ HA , and the vectors |b1 i, |b2 i
∈ HB . When the composite system is in the state ρ AB , the
reduced density operator for system A is defined as ρA =
TrHB (ρ AB ) and the reduced density operator for system B is
defined as ρB = TrHA (ρ AB ).
Reduced state consensus for a quantum network can be
defined as follows.
Definition 2 (Reduced State Consensus): [26] A quantum network consisting of m nodes with the state ρ̄ is in
Reduced State Consensus (RSC) if

ρ̄1 = ρ̄2 = ... = ρ̄m ,
where ρ̄ j = Tr⊗k6= j Hk (ρ̄ ) ( j = 1, 2, ..., m) are defined as the
reduced density operator for node j and can be calculated
according to Definition 1.
We aim to steer a quantum network into a consensus
state defined in Definition 2. In practical applications, the
existence of noise (including extrinsic and intrinsic), inaccuracies (e.g., inaccurate operation in the coupling between
nodes) and fluctuations (e.g., fluctuations in control fields)
in quantum networks is unavoidable. We assume that the
Hamiltonian with uncertainties can be written into
M

HΘ (t) = θ0 H0 + ∑ θ j u j (t)H j .
j=1

(8)

max J(u) := max E[JΘ (u, ρ̄ )]
u

u

s.t.
"
#

M




ρ̇ (t) = −i θ0 H0 + ∑ θ j u j (t)H j , ρ (t)

(9)

j=1


ρ (0) = ρ 0 , t ∈ [0, T ]




θ j ∈ [1 − E, 1 + E], j = 0, 1, 2, ...M

where E[JΘ (u, ρ̄ ] denotes the average performance function
with respect to the parameter fluctuations Θ and the target
consensus state ρ̄ , and E ∈ [0, 1] are the bounds of the
parameter uncertainties.
C. Femtosecond laser quantum control
We consider the experimental control of molecular fragmentation using shaped femtosecond laser pulses. Here,
CH2 BrI is chosen as the target molecule. As a family member
of Halomethane molecules, whose dissociative product plays
a central role in ozone depletion, CH2 BrI has attracted
wide attention because of its importance in environmental
chemistry. In addition, it is one of the simplest prototype
molecules containing different bonds, a stronger C-Br bond
and a weaker C-I bond, which is ideal for the study of controlling the selective bond-breaking. Under strong fields with
femtosecond laser pulses, CH2 BrI molecules will undergo
ionization and dissociation, and their charged products can
be separated and detected with a TOF-MS. In particular, we
choose to optimize the photoproduct ratio of CH2 Br+ /CH2 I+
as our control objective, which corresponds to breaking the
strong bond versus the weak bond. We apply closed-loop
learning control, using our proposed msms DE algorithm, to
search for a robust ultrafast laser pulse that maximizes this
ratio.
In closed-loop learning control, the learning process can
be conceptually demonstrated as follows. First, one applies
trial input pulses to the molecules to be controlled and
observes the results. Second, a learning algorithm suggests
better control inputs based on the prior experiments. Third,
one applies “better” control inputs to new molecules [2]. This
approach has been employed to explore the quantum control
landscape [36] to find the optimal control strategy where
the control performance function J(u)reaches its maximum.
In order to achieve good performance, we first need to
identify a reference phase mask on the SLM that give shortest
transform limited (TL) pulse, which can be obtained from
optimizing the signal of two-photon absorption (TPA). We
first use the proposed msms DE algorithm to search for a
good control to obtain highest TPA signal. Then we apply
the same algorithm to search for a good control for the
fragment ratio of CH2 Br+ /CH2 I+ . The consideration of
robustness with multiple samples (MS) in DE would also
ensure good transferability of the experimental results or
photonic reagents [37]. That is, an optimal pulse identified
from one laser system would also perform well (if not optimal) when transferred to another system despite the minor

differences or uncertainties lying in the control parameters
(i.e., the spectral phases on the SLM).
III. D IFFERENTIAL

EVOLUTION AND

msMS DE

ALGORITHM

DE was initially proposed in 1990s by Storn and Price
[9], [39] and has derived many variants. DE algorithms have
witnessed wide applications in diverse domains of science
and engineering [10], [38]. In this section, we first briefly
introduce the conventional DE algorithm, analyze DE with
variant strategies and control parameters, and then propose
an improved DE algorithm of msMS DE for these quantum
robust control problems presented in Section II.
A. DE Algorithm
In DE, the individual trial solutions (which constitute a
population) are termed as parameter vectors or genomes,
usually represented in a vector X = [x1 , x2 , · · · , xD ]T where
each parameter xi is a real number. Solving an optimization
problem using DE is basically a search for a parameter vector
to minimize or maximize a fitness function (or objective
function) f (X). The following operations are used to evolve a
population of D-dimensional parameter vectors until a “best”
individual is generated and found.
(a) Initialization. DE searches for a global optimum
point in a D-dimensional real parameter space ℜD . Here,
we denote the population at the current generation as
j
j
j
Xi,G = (x1i,G , · · · , xD
i,G ), i = 1, ..., NP and let xi,G ∈ [xmin , xmax ],
( j = 1, 2, ..., D) since these parameters are usually related to
physical variables with relevant bounds. We usually initialize
the population (at G = 0) as follows [43]:
j

j

j

j
xi,0 = xmin + rand(0, 1) · (xmax
− xmin ),

j = 1, 2, ..., D, (10)

where rand(0, 1) is a uniformly distributed random number.
(b) Mutation. In DE, the key to “mutation” is to generate
a difference vector by choosing three other distinct parameter
vectors from the current generation (say, Xr1 , Xr2 , Xr3 ). The
indices r1 , r2 , r3 ∈ {1, ..., NP} are mutually exclusive integers
randomly generated within the range [1, NP] and r1 , r2 , r3 6= i.
The donor vector Vi,G+1 are generated by
Vi,G = Xr1 ,G + F · (Xr2 ,G − Xr3 ,G ),

(11)

where the scaling factor F is a positive control parameter
typically in the interval [0.4, 1].
(c) Crossover. To enhance the potential diversity of the
population, a crossover operation comes into play after
the mutation. The DE family has two types of crossover
operations (i.e., exponential and binomial). The binomial
(uniform) crossover is outlined as
( j
vi,G , if rand( j) ≤ CR or j = rand(1, D),
j
ui,G =
(12)
j
xi,G , if rand( j) > CR and j 6= rand(1, D),
where j = 1, 2, ..., D and rand( j) ∈ [0, 1] is a uniform random
number. CR is a user-specified constant within the range
[0, 1), rand(1, D) ∈ {1, 2, ..., D} is a randomly chosen index.

(d) Selection. To keep the population size constant over
subsequent generations, DE uses the following selection
operation to determine whether the target vector or the trial
vector survives to the next generation:
(
Ui,G , if f (Ui,G ) ≥ f (Xi,G ),
Xi,G+1 =
(13)
Xi,G , otherwise.
If the new trial vector yields an equal or lower value of the
objective function, it replaces the corresponding target vector
in the next generation; otherwise the target vector survives.
Usually, it is the mutation operation that demarcates
one DE scheme from another. The DE strategy with the
mutation in (11) is referred to as “DE/rand/1” using the
notation “DE/x/y”, where x represents a string denoting the
base vector to be perturbed, y is the number of difference
vector considered for perturbation of x. Furthermore, when
crossover is also considered, the notation “DE/x/y/z” is used,
where z stands for the type of crossover (bin: binomial, exp: exponential). DE variants with different mutation
strategies usually have different performance for solving
optimization problems. “DE/rand/1/bin” is a commonly used
strategy and it usually shows slow convergence speed and
strong exploration capacity. The strategies relying on the
best solution found so far such as “DE/rand-to-best/1/bin”,
“DE/best/2/bin”, usually have a fast convergence speed and
perform well when solving unimodal problems. However,
they are more likely to get stuck at a local optimum and
lead to premature convergence for multimodal problems.
Two-difference-vectors-based strategies may result in better
perturbations than one-difference-vector-based strategies. For
control parameters of DE, Storn and Price [39] have proposed
that a good initial choice of F was 0.5 and the range of
F is usually set [0.4, 1]. The crossover rate CR may be
CR ∈ [0, 1]. Several results also proposed the techniques of
self-adaptation to automatically find an optimal set of control
parameters [10], [40] to provide improved performance.
B. msMS DE algorithm
When implementing DE, users need to determine the
appropriate mutation-trial strategies and parameter settings
to ensure the success of the algorithm [40]. It is a highcost practice to perform a trial-and-error search for the most
appropriate trial vector generation strategy and fine-tune its
associated control parameter values, i.e., the values of CR,
F for a given problem. Moreover, during different stages of
evolution, different trial vector generation strategies coupled
with specific control parameter values can be more effective
and single strategy may result in premature convergence
thus leading to a failure in complex problems such as
nonseparable and multimodal functions [10], [41], [42]. Also,
several variants of DE utilizing the idea of mixed strategies
such as SaDE [40] and EPDE [43] have been proposed and
exhibited good performance. Our numerical results show that
DE with a single strategy might be enough for easy problems
while DE variants with mixed strategies might emerge as a

promising candidate for quantum control problems with multimodal landscapes. Existing results of sampling-based learning control [19], [31], [44] have shown that the employment
of an average objective function with multiple samples can
provide improved performance for quantum robust control
problems. Inspired by these observations, we adopt a mixed
strategy and an average performance of multiple samples to
present an improved DE algorithm (i.e., msMS DE) for these
quantum robust control problems outlined in Section II.
We first choose one mutation scheme from a pool of
strategy candidates where several mutation schemes with
effective yet diverse characteristics are equally distributed.
Subsequently, a binomial crossover operation is performed
on the corresponding mutant vector to generate the trial
vector. Note that we assign various values of F and CR for
each individual during the current generation to increase the
diversity of the population. To construct the candidate pool,
we investigate several commonly used mutation strategies
[42] and select four strategies with distinct capabilities at
different stages of evolution as follows:
DE/rand/1:
Vi = Xr1 + F · (Xr2 − Xr3 ).

(14)

DE/rand to best/2:
Vi = Xi + F · (Xbest − Xi ) + F · (Xr1 − Xr2 ) + F · (Xr3 − Xr4 ).
(15)
DE/rand/2:
Vi = Xr1 + F · (Xr2 − Xr3 ) + F · (Xr4 − Xr5 ).

(16)

DE/current-to-rand/1:
Vi = Xi + K · (Xr1 − Xi ) + F · (Xr2 − Xr3 ).

(17)

The indices r1 , r2 , r3 , r4 and r5 are mutually exclusive integers
randomly chosen from the range [1, NP] and all of them
are different from the index i. Xbest is the best individual
vector with the best fitness (i.e., the lowest objective function
value for a minimization problem) in the population. The
control parameter K in the strategy DE/current-to-rand/1 is
set as K = 0.5 to eliminate one additional parameter. As
for the crossover operation, the first three mutation schemes
are combined with a binomial crossover operation, while
the fourth scheme directly generates trial vectors without
crossover.
In the proposed msMS DE algorithm, the parameter F is
approximated by a normal distribution with mean value 0.5
and standard deviation 0.3, denoted by N(0.5, 0.3). It is easy
to verify that values of F fall into he range [−0.4, 1.4] with
probability of 0.997 which helps maintain both exploitation
(with small F values) and exploration (with larger F values).
Similarly, we assume CR obeys a normal distribution denoted
by N(0.5,0.1) and the small standard deviation 0.1 is enough
to guarantee that most values of CR lies in [0, 1] [40]. Consequently, a set of F and CR values are randomly sampled
from normal distribution (denoted by Normrnd) and applied
to each target vector in the current population. We may obtain
some extraordinary values far from [-0.4,1.4] for scale factor

F and we usually accept them to increase diversity. While
the crossover rate has probabilistic meaning for the chance
of survival, we should abandon those fall outside [0, 1], and
generate another valid parameter by CR = N(0.5, 0.1) to
guarantee the practical meaning of crossover.
The msMS DE method is proposed for three classes of
quantum control tasks. In order to design appropriate control
laws to achieve good robustness performance, we integrate
the idea of sampling-based learning control [19] into the
msMS DE algorithm. To begin with, we prepare N samples
Θk = (θ0 , θ1 , ..., θM ) (k = 1, 2, . . . , N) with different values of
the uncertainty parameters. We compute the fitness values of
these sample vectors . Then, we evaluate the average fitness
value f¯ for these samples, and f¯ is defined as follows
1 N
f¯(Ui,G ) = ∑ f (Ui,G , Θk ),
N k=1
1 N
f¯(Xi,G ) = ∑ f (Xi,G , Θk ).
N k=1
The algorithmic description of the msMS DE is presented in
Algorithm 1.
Remark 1: In simulation, after we obtain the nominal
value of an individual, we may generate the other samples by
perturbing the nominal value and then calculate the average
fitness function. In experiment, we need to run an experiment
for each sample to measure the fitness of each sample, and
then calculate the average fitness. Usually, a larger number
of samples N may lead to a better robustness performance
[19]. However, the computational or experimental time will
significantly increase with the increase of N. In this paper, we
use three samples for each uncertainty parameter for saving
computational and experimental time.
Remark 2: After performing the mutation operation, we
obtain new donor vectors, and some of them might survive
into the next generation and serve as parents. Therefore, we
add a procedure where each vector is evaluated in view of
boundary constraints. If any parameter of the vector falls
beyond the pre-defined lower or upper bounds, we will
replace it with a random value within the allowed range.
Remark 3: In our proposed msMS DE, we preset a maximum generation Gmax as the termination criterion. During
the implementation of the algorithm, the population evolves
until the learning process reaches G = Gmax . In numerical
examples, we let Gmax = 50000. In experimental examples,
we choose Gmax = 150.
IV. N UMERICAL

RESULTS FOR ENSEMBLE CONTROL AND

QUANTUM NETWORK CONSENSUS

A. Control of open inhomogeneous two-level quantum ensembles
We consider a specific inhomogeneous open two-level
ensemble with inhomogeneous parameter bound E = 0.2.
Members of the ensemble are governed by the following
Hamiltonian:
1
H(t) = θ0 σz + θ1 u(t)(cos ϕσx + sin ϕσy ),
(18)
2

Algorithm 1.
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:

Algorithmic description of msMS DE

Set the generation number G = 0
for i = 1 to NP do
for j = 1 to D do
j
j
j
j
xi,G = xmin + rand(0, 1) · (xmax − xmin)
end for
end for
Initialize fitness f (Xi,G ) and evaluate vector by f
Mark the best vector with maximum f as Xbest,G
repeat (for each generation G = 1, 2, . . . , Gmax )
repeat (for each vector Xi , i = 1, 2, . . . , NP)
Set parameter Fi,G = Normrnd(0.5, 0.3)
Randomly generate a real number pp∈ [0, 1],
if pp > 0 and pp ≤ 0.25 then flag=1
Vi,G = Xr1 ,G + Fi,G · (Xr2 ,G − Xr3 ,G )

end if
if pp > 0.25 and pp ≤ 0.5 then flag=2

Vi,G = Xi,G + Fi,G · (Xbest,G − Xi,G ) + Fi,G · (Xr1 ,G − Xr2,G )

18: +Fi,G · (Xr3 ,G − Xr4 ,G )
19:
end if
20:
21:

if pp > 0.5 and pp ≤ 0.75 then flag=3

22:

end if
if pp > 0.75 and pp ≤ 1 then flag=4

23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:
45:
46:
47:
48:
49:
50:
51:
52:
53:

Vi,G = Xr1 ,G + Fi,G · (Xr2 ,G − Xr3 ,G ) + Fi,G · (Xr4 ,G − Xr5 ,G )

Vi,G = Xi,G + K · (Xr1 ,G − Xi,G ) + Fi,G · (Xr2 ,G − Xr3,G )

end if
for i = 1 to D do
j
j
j
j
if vi,G
> xmax
or vi,G
< xmin
then
j
j
j
j
vi,G = xmin + rand(0, 1) · (xmax
− xmin
)
end if
end for
Set parameter CRi,G = Normrnd(0.5, 0.1)
while CRi,G < 0 or CRi,G > 1 do
CRi,G = Normrnd(0.5, 0.1)
end while
if flag=1 or flag=2 or flag=3 then
for j = 1 to D do
j
j
ui,G
= vi,G
, if (rand[0, 1] ≤ CRi,G or j = jrand )
j

j

ui,G = xi,G , otherwise

end for
end if
if flag=4 then Ui,G = Vi,G
end if
for each sample Θk , (k = 1, 2, ..., N) do
evaluate the fitness function f (Ui,G , Θk )
end for
Compute f¯(Ui,G ) = N1 ∑Nk=1 f (Ui,G , Θk )
if f¯(Ui,G ) ≥ f¯(Xi,G ) then
Xi,G+1 ← Ui,G , f¯(Xi,G+1 ) ← f¯(Ui,G ).
end if
Renew the best vector Xbest,G and i ← i + 1
until i = NP
G ← G+1
until G = Gmax

where ϕ ∈ [0, 2π ], u(t) ∈ [−10, 10], and the Pauli operators
are defined as:






0 1
0 −i
1 0
, σy =
, σz =
. (19)
σx =
1 0
i 0
0 −1
Let γk = 1 and the Lindblad operators are given by [45]






0 0
0 0.2
0.2 0
L1 =
, L2 =
, L3 =
.
0.1 0
0 0
0 0
(20)
For a two-level quantum ensemble, U1 ,U2 and U3 in (3)
can be chosen as σx , σy and σz . The coherent vector for the
density matrix is the Bloch vector
r = (x, y, z)T = (tr(σx ρ ), tr(σy ρ ), tr(σz ρ )).
The dynamical equation for r can be written as




−0.045
−θ0
0
0
 r(t) +  0 
θ0
−0.045
0
ṙ(t) = 
0
0
−0.05
0.03
0
0
−2 sin ϕ
0
0
2 cos ϕ  r(t).
+θ1 u(t) 
2 sin ϕ −2 cos ϕ
0
(21)
The average fitness function is given as
E[JΘ (u)] =

1
1
∑[1 − 4 k r f − rθ0 ,θ1 (T ) k2],
N∑
θ θ
0

(22)

1

where N is the total number of the chosen samples. An upper
bound of the fitness function is 1 although we do not know
the maximum that can be achieved. In msMS DE algorithm,
we choose three samples for each parameter, and here we
have N = 9. During learning control of the inhomogeneous
quantum ensemble, we employ DE algorithms to seek for
the optimal control u∗ (t). Then, we apply the optimal control
field to additional samples with inhomogeneous parameters
(θ0 , θ1 ) following uniform distributions within [0.8, 1.2] to
test its performances. We assume that the initial state and
the target state are, respectively,




1 0
0 0
ρ0 =
,
ρf =
.
(23)
0 0
0 1
The target time T = 10 and the time interval [0, T ] is
equally divided into D = 200 time steps, and ∆t = 0.05. The
population size is set as NP = 50 for all the algorithms in
this example. The simulation is implemented on MATLAB
platform (version 8.3.0.532). The hardware environment for
simulation is Intel(R)-Core(TM) i7-6700K CPU, dominant
frequency @4.00GHz, and 16G(ARM).
To demonstrate the performance of the proposed msMS DE algorithm for the control problem of inhomogeneous
quantum ensembles, we make performance comparison between it and ms DE (DE with multiple samples, i.e., using
the average fitness function of multiple samples) with various
parameters. To begin with, we present the results for the
traditional DE (i.e., “DE/rand/1/bin”) using multiple samples
with three typical sets of control parameters. Three cases
with different control parameters are labeled as “ms DE1”

B. Consensus in superconducting qubit networks

J(u)

0.96
0.95
ms_DE1
ms_DE2
ms_DE3

0.94

(a)
0.93
0

10000

20000

30000

40000

50000

Iterations

J(u)

0.98
0.96
ms_DE1
msMS_DE
GA

0.94
0.92

(b)

0.9
0

10000

20000

30000

40000

50000

Iterations

Fig. 1.
(a) The training performance of the two-level open quantum
ensemble via ms DE1 (F = 0.9,CR = 0.1), ms DE2 (F = 0.9,CR = 0.9)
and ms DE3 (F = 0.5,CR = 0.3). (b) The training performance of ms DE1,
GA and msMS DE.

(F = 0.9,CR = 0.1), “ms DE2” (F = 0.9,CR = 0.9), and
“ms DE3” (F = 0.5, CR = 0.3), and the training performance
is presented in Fig. 1(a). It is clear that ms DE1 and ms DE3
have better performance than ms DE2 for the quantum
control problem. ms DE1 can achieve the highest fitness
0.9566 among these three cases. We then compare the training performance of ms DE1, GA and msMS DE, and the
results are illustrated in Fig. 1(b). The msMS DE algorithm
achieve the highest fitness Jmax = 0.9798, while ms DE1 and
GA converge to a maximum value of 0.9566 and 0.9667,
respectively. A comparison of testing performance for 2000
additional samples and training time between DE1 (using one
sample), ms DE1, ms DE2, ms DE3, GA (with crossover
probability Pc = 0.8 and mutation probability Pm = 0.05)
and msMS DE in Table 1 shows that msMS DE is superior
to ms DE and DE1. More numerical results also show
that msMS DE usually can find the control field with the
best robustness among these algorithms because msMS DE
employs mixed mutation strategies as well as average performance using multiple samples. ms DE1, ms DE2, ms DE3,
GA and msMS DE also take similar time to find an optimal
solution for the ensemble control problem. For example,
msMS DE takes 9 hours 20 minutes and GA takes 10 hours
18 minutes and 14 seconds.
Algorithm
DE1
ms DE1
ms DE2
ms DE3
GA
msMS DE

parameters
F = 0.9,CR = 0.1, N = 1
F = 0.9,CR = 0.1, N = 9
F = 0.9,CR = 0.9, N = 9
F = 0.5,CR = 0.3, N = 9
Pc = 0.8, Pm = 0.05, N = 9
F = N(0.5, 0.3),
CR = N(0.5, 0.1), N = 9

training time
1h10m47s
9h27m47s
9h20m15s
9h40m5s
10h18m14s
9h20m0s

J̄(u)
0.9408
0.9610
0.9537
0.9601
0.9691
0.9803

TABLE I
P ERFORMANCE COMPARISON OF DIFFERENT ALGORITHMS

The nodes in a quantum network could be photons, electrons, or other quantum systems. In this section, we consider
a quantum network that consists of superconducting qubits
as its nodes. Superconducting quantum circuits based on
Josephson junctions are one promising candidate for building
hardwares of quantum computers. These macroscopic circuits can behave quantum mechanically like artificial atoms that
can also be used to test the laws of quantum mechanics on
macroscopic systems [46]. Superconducting qubits have been
widely investigated theoretically and implemented experimentally since they could be easily embedded in nanometerscale electronic devices and scaled up to provide a large number of qubits for quantum computation [47]. Manipulation of
superconducting qubits can be achieved by adjusting external
parameters such as currents and voltages or by tuning the
coupling between two superconducting qubits [48].
In superconducting circuit, the charging energy EC and
the Josephson coupling energy EJ have significant effect on
the quantum mechanical behavior of a Josephson-junction
circuit. Different kinds of superconducting qubits can be
realized according to the regimes of EJ /EC . For example,
a charge qubit can form when EC ≫ EJ . In practical applications, the Josephson junction in the charge qubit is usually replaced by a dc superconducting quantum interference
device (SQUID) with low inductance and a magnetic flux
[49]. The equivalent Hamiltonian of a charge qubit can be
described as [50], [51]
H = Fz (Vg )σz − Fx (Φ )σx ,

(24)

where Fz (Vg ) can be adjusted through external voltage Vg ,
and Fx (Φ ) corresponds to a tunable effective coupling with
the external magnetic flux Φ in the SQUID. Hence, Fz (Vg )
and Fx (Φ ) are related to external control fields.
Now, consider a quantum network consisting of three
superconducting qubits with control fields acting on all
(12)
(23)
qubits. We denote σx = σx ⊗ σx ⊗ I, σx = I ⊗ σx ⊗ σx ,
(13)
σx = σx ⊗ I ⊗ σx . Its free Hamiltonian can be described as
(12)

H0 = ω12 σx
(1)

(23)

+ ω23σx

(13)

+ ω13 σx

(2)

.

(25)

(3)

Denote σx = σx ⊗ I ⊗ I, σx = I ⊗ σx ⊗ I, σx = I ⊗ I ⊗ σx ,
(1)
(2)
(3)
and σz = σz ⊗ I ⊗ I, σz = I ⊗ σz ⊗ I, σz = I ⊗ I ⊗ σz . We
have the control Hamiltonian in the following form
(1)

(1)

(2)

(2)

(3)

(3)

Hu (t) = ux1 σx + uz1 σz + ux2 σx + uz2 σz + ux3 σx + uz3 σz .
(26)
Our task is to drive the quantum network from an arbitrary
initial state (usually three qubits having different reduced
states) to a consensus state. Furthermore, if we withdraw the
external control fields, the quantum network will remain in
the consensus state under the free Hamiltonian. Denote 1n as
an n-dimensional matrix with all of its elements being 1. Let
the target state be ρ̄ = 81 18 . We have the following result.
Proposition 1: The state ρ̄ = 18 18 is a consensus state
for the three qubit network. Also, ρ̄ is invariant under

1

J(u)

0.9
msMS_DE
DE1

0.8

0.7

0

10000

20000

30000

40000

50000

Iterations
Fig. 2. The initial and target (reduced) states of the three qubits on the
Bloch sphere.

(12)

(23)

the action of free Hamiltonian H0 = ω12 σx + ω23 σx +
(13)
ω13 σx .
Proof: For ρ̄ = 18 18 , we can calculate the reduced states
for three notes as follows:
1
1
1
ρ̄1 = 12 , ρ̄2 = 12 , ρ̄3 = 12 .
2
2
2
It is clear that ρ̄1 = ρ̄2 = ρ̄3 , that is, the state ρ̄ is a consensus
state for the three qubit network according to Definition 2.
A direct calculation shows that [H0 , ρ̄ ] = 0. Hence,

ρ̄˙ = [H0 , ρ̄ ] = 0.
That is, ρ̄ is invariant under the action of free Hamiltonian
H0 .
The chosen target state ρ̄ is a symmetric state from which
the superconducting qubit network will keep invariant with
only free Hamiltonian H0 . The initial state is set as ρ̄ 0 =
ρ10 ⊗ ρ20 ⊗ ρ30 where


 1



1 0
0 0
− 12
0
2
ρ10 =
, ρ20 =
,
ρ
=
.
1
3
0 0
0 1
− 12
2
The initial and target states of qubit network are illustrated
in Fig. 2.
In practical applications, there may exist fluctuations in
magnetic fields and electric fields in superconducting qubits.
The practical control Hamiltonian is assumed to be
(1)
(1)
(2)
(2)
Hu (t) =θx ux1 σx + θz uz1 σz + θx ux2 σx + θz uz2 σz
(3)
(3)
+ θx ux3 σx + θz uz3 σz .

(27)

We apply the msMS DE algorithm to search for the robust control fields to reach a consensus state in the above
quantum network. Stimulation parameters are set as: the
population size is set as NP = 100, the time internal [0, 20]
ns is equally divided into 100 smaller time intervals (i.e.,
D = 100), the control terms ux1 , uz1 , ux2 , uz2 , ux3 , uz3 ∈ [0, 1]
GHz. Let ω12 = ω23 = ω13 = 0.1 GHz. We assume that
θx ∈ [0.98, 1.02] and θz ∈ [0.98, 1.02] (i.e., E = 0.02). For
each uncertainty parameter, we choose three samples and
we have N = 9 samples for training. We employ DE1
(one sample) for comparison. The training performance of
driving qubit network is illustrated in Fig. 3. As we can
see, msMS DE achieves a rather high fitness Jmax = 0.9988
(where 1 is an upper bound), while DE1 achieves the fitness

Fig. 3. The training performance of superconducting qubit network via
DE1 and msMS DE.

of Jmax = 0.9561. For the case θx = θz = 1, Fig. 4 shows
the reduced states of three qubits from different trajectories
asymptotically converging to the same trajectory using the
control learned from msMS DE. Based on the control fields
from DE1 and msMS ED, we test 2000 additional samples
regarding the trace distance defined as (for i = 1, 2, 3)
r
1
1
1
1
1 2
||ρi − 12 ||Tr = Tr (ρi − 12 )† (ρi − 12 ) = ∑ |λ j |
2
2
2
2
2 j=1
where λ j are eigenvalues of ρi − 12 12 . Since the maximum
trace distance between two quantum states may be 1, we
define the relative error between two quantum states ρi and
ρ j as ||ρi − ρ j ||Tr × 100%. Fig. 5 shows that the relative
error between each qubit and its target state always keeps
below 1.2% for the case using msMS DE while the relative
error between each qubit and its target state may exceed
10.0% for the case using DE1. We further show the trace
distances between the quantum states of different qubits after
the control Hamiltonian is withdrawn in Fig. 6. It is clear
that the relative errors between the reduced states are always
below 2.0% for the case of msMS DE while the relative
errors may exceed 12.0% for the case of DE1. The results
demonstrate that the approximate consensus state achieved
using msMS DE has much better stability than that obtained
using DE1.
V. E XPERIMENTAL RESULTS ON

FEMTOSECOND LASER
CONTROL SYSTEMS

The following quantum control experiments were carried
out in Department of Chemistry at Princeton University.
A. Experimental setup
The experimental setup contains three major components:
1) a fs laser system, 2) a pulse shaper and 3) a time-of-flight
mass spectrometry (TOF-MS). Briefly, the fs laser system
(KMlab, Dragon) consists of a Ti:sapphire oscillator and
a amplifier, which produces 1 mJ, 25 fs pulses centered
at 790 nm. The laser pulses are introduced into a pulse
shaper with a programmable dual-mask liquid crystal spatial
light modulator (SLM). The SLM has the capability of
independent phase and amplitude modulation and has 640
pixels with 0.2 nm/pixel resolution [52], [53]. Typically,
every 8 adjacent pixels are bundled together to form an array
of 80 “grouped pixels”, which are the control variables. Each

1

trace distance

r3=(0,0,1)

qubit1
qubit2
qubit3

o

qubit1−qubit2
qubit2−qubit3
qubit1−qubit3

0.15

0.1

0.05
20

25

30

35

40

t

o

z

(a)

0

o

−1
1

trace distance

0.02

rtarget=(1,0,0)

r2=(−1,0,0)

qubit1−qubit2
qubit2−qubit3
qubit1−qubit3

0.01

(b)

o
20

25

30

1

r1=(0,0,−1)

0

−1

y

−1

35

40

t

0

x

Fig. 4.
The reduced states of three qubits from different trajectory
asymptotically converge to the same trajectory for θx = θz = 1 using the
control learned by msMS DE.

Fig. 6. The trace distances between different qubits in the qubit network
with free Hamiltonian. (a) Average evolution curves of trace distances
between three qubits for 2000 samples (DE1). (b) Average evolution curves
of trace distances between three qubits for 2000 samples (msMS DE).

trace distance

B. Optimization of TPA signal
1
0.8
qubit1
qubit2
qubit3

0.6
0.4

(a)

0.2
0

0

10

20

30

40

trace distance

t
1
0.8
qubit1
qubit2
qubit3

0.6
0.4

(b)

0.2
0

0

10

20

30

40

t
Fig. 5. The testing performance of the qubit network. (a) Average evolution
curves of trace distances of three qubits for 2000 samples using the control
field learned by DE1. (b) Average evolution curves of trace distances of
three qubits for 2000 samples using the control field learned by msMS DE.

control variable can have a phase value between 0 and 2π ,
and an amplitude value between 0 and 1. In this experiment,
we do phase-only control, with all the amplitude values
fixed at 1. The shaped laser pulses out of the shaper are
focused into a vacuum chamber, where photoionization and
photofragmentation occurs. The fragment ions are separated
with a set of ion lens and passing through a TOF tube
before being collected with an MCP detector. The MS signals
are recorded with a fast oscilloscope, which accumulates 1
second with 3000 laser shots each time before sending to
a personal computer for further analysis. A small fraction
of the beam (< 5%) is separated from the main beam and
focused into a GaP photodiode (Thorlab, DET25K), which
collects signals arising from two photon absorption.

A preliminary task is to optimize the two-photon absorption (TPA) signal, which is a convenient way to identify a
shortest pulse that removes the residual highorder dispersion
in the amplifier output. The parameter setting is as follows:
D = 80, NP = 30 and N = 3. The control variable are the
phases bounded in between 0 and 2π . The fitness function
corresponds to the TPA signal and an average fitness of
three samples was used. Three samples for each individual
were selected as follows: The first sample came from the
T
current individual, denoted as Xi1 = [x1i , x2i , . . . , x80
i ] , the
second sample was generated by adding a random fluctuation between 0 and 0.1π to each component of the
current individual, i.e., Xi2 = [x1i + 0.05rand(0, 1) × 2π , x2i +
T
0.05rand(0, 1) × 2π , . . . , x80
i + 0.05rand(0, 1) × 2π ] , and the
3
1
third sample was selected as Xi = [xi − 0.05rand(0, 1) ×
T
2π , x2i − 0.05rand(0, 1) × 2π , . . ., x80
i − 0.05rand(0, 1) × 2π ] .
This means that each control variable is permitted to have
up to 5% (of the maximum phase) additive noise. The
interaction between the algorithm and the modulation of
SLM is accomplished by LabVIEW system design software.
An experimentally reasonable termination condition of 150
generations (iterations) is used. For 150 iterations, it takes
around five and a half hours to run the experiment. For each
generation, a total of 90,000 signal measurements were made.
The experimental result is shown in Figure 7 where TPA
signal during each iteration is presented in Fig. 7(a) and the
optimized phases of 80 control variables for the final optimal
result is given in Fig. 7(b). After 150 generations, the best
average TPA signal for three samples can reach 1.35.
C. Fragmentation control
We consider the fragmentation control of molecules
CH2 BrI, where the fitness is defined as the photofragment
ratio of CH2 Br+ /CH2 I+ , while the control variables are the
phases. The parameter setting is the same as that in Section

Fig. 7. Experimental result for optimizing the TPA signal using msMS DE.
(a) TPA signal vs iterations, where ‘Best’ represents the maximum fitness
and ‘Average’ represents the average fitness of all individuals during each
iteration. (b) Optimized phases of 80 control variables for the final optimal
result corresponding to the maximum fitness.

V.B: D = 80 and NP = 30. DE algorithms were employed to
optimize the phases of 80 control variables. Here, we apply
both DE1 and msMS DE for comparison.
Figure 8 shows the experimental results using DE1, where
the ratio CH2 Br+ /CH2 I+ as the fitness function is presented
in Fig. 8(a) and the optimized phases of 80 control variables
for the final optimal result is given in Fig. 8(b). In Fig.
8(b), ‘Best’ represents the maximum fitness and ‘Average’
represents the average fitness of all individuals during each
iteration. With 150 iterations, DE1 can find an optimized
pulse to make CH2 Br+ /CH2 I+ achieve 2.41.
Figure 9 shows the results from the msMS DE algorithm,
in which three samples are measured in each experiment.
Three samples for each individual were selected using the
same method as that in the experiments of optimizing TPA
signals. With 150 iterations, msMS DE can find an optimal
pulse for making the average CH2 Br+ /CH2 I+ of three
samples to achieve 2.67. The experimental results are shown
in Fig. 9, where the average ratio CH2 Br+ /CH2 I+ of three
samples as the fitness function is presented in Fig. 9(a) and
the optimized phases of 80 control variables for the final
optimal result is given in Fig. 9(b).
After we obtained the optimal femtosecond control pulses
using DE1 and msMS DE, we can test the performance
of the optimal pulses. The testing results are shown in
Fig. 10, where Fig. 10(a) is the average TOF signal of
100 testing results with random noises between −7.5% and
+7.5% (with respect to the maximum phase 2π ) for the
femtosecond pulse optimized by DE1. In other words, if
T
we denote the best individual as Xbest = [x1b , x2b , . . . , x80
b ] ,
k
1
these 100 testing samples can be written as Xs = [xb +
0.075(2rand(0, 1) − 1) × 2π , x2b + 0.075(2rand(0, 1) − 1) ×
T
2π , . . ., x80
b + 0.075(2rand(0, 1) − 1) × 2π ] . Fig. 10(b) shows
the average TOF signal of 100 testing results for the
femtosecond pulse optimized by msMS DE. The average
CH2 Br+ /CH2 I+ of the 100 testing samples can achieve

Fig. 8. Experimental result for optimizing the ratio between the product
CH2 Br+ and the product CH2 I+ using DE1. (a) Ratio CH2 Br+ /CH2 I+
vs iterations, where ‘Best’ represents the maximum fitness and ‘Average’
represents the average fitness of all individuals during each iteration. (b)
Optimized phases of 80 control variables for the final optimal result
corresponding to the maximum fitness.

Fig. 9.
Experimental result for optimizing the ratio between the
product CH2 Br+ and the product CH2 I+ using msMS DE. (a) Ratio
+
CH2 Br /CH2 I+ vs iterations, where ‘Best’ represents the maximum fitness
and ‘Average’ represents the average fitness of all individuals during each
iteration. (b) Optimized phases of 80 control variables for the final optimal
result corresponding to the maximum fitness.

2.61 for the pulse from msMS DE while the average
CH2 Br+ /CH2 I+ is only 2.12 for the pulse from DE1. It is
clearly shown that msMS DE outperforms DE1 in terms of
reaching a better objective fitness value in this experiment.
VI. CONCLUSION
In order to solve three classes of quantum robust control
problems, we have proposed improved msMS DE using
multiple samples for fitness evaluation and a mixed strategy
for mutation. The msMS DE algorithm shows excellent
performance for the control problem of open inhomogeneous
quantum ensembles and the consensus problem of quantum
networks with uncertainties. We have experimentally implemented msMS DE on femtosecond laser control systems
in the laboratory to generate good TPA signal and control

Fig. 10. Experimental results for TOF signal where m represents mass
and q represents charge. (a) The average TOF signal of 100 testing results
with random noises between −7.5% and +7.5% for the femtosecond pulse
optimized by DE1. (b) The average TOF signal of 100 testing results
with random noises between −7.5% and +7.5% for the femtosecond pulse
optimized by msMS DE.

fragmentation of chemical molecules. In future research,
there is plenty of room for exploring the use of DE for emerging quantum control engineering. For example, it is worth
adapting more efficient DE algorithms for high-dimensional
quantum control problems [54]. In the laboratory, many
quantum control problems involve multiple even many objectives that need to be optimized, and it is also worth adapting
multi-objective and many-objective evolution optimization
algorithms [55], [56] for these challenging quantum control
problems.
R EFERENCES
[1] D. Dong, and I. R. Petersen, “Quantum control theory and applications:
a survey,” IET Control Theory & Applications, vol. 4, no. 12, pp. 26512671, 2010.
[2] H. Rabitz, R. De Vivie-Riedle, M. Motzkus, and K. Kompa, “Whither
the future of controlling quantum phenomena?” Science, vol. 288, no.
5467, pp. 824-828, 2000.
[3] H. M. Wiseman, and G. J. Milburn, “Quantum Measurement and
Control,” Cambridge, England: Cambridge University Press, 2010.
[4] M. A. Nielsen, and I. L. Chuang, “Quantum Computation and Quantum Information,” Cambridge, U.K.: Cambridge University. Press,
2000.
[5] S. Kuang, D. Dong, and I. R. Petersen, “Rapid Lyapunov control
of finite-dimensional quantum systems” Automatica, in press, 2017,
online: quant-ph, arXiv, https://arxiv.org/abs/1512.02545.
[6] C. Chen, D. Dong, H.-X. Li, J. Chu, and T.-J. Tarn, “Fidelity
based probabilistic Q-learning for control of quantum systems,” IEEE
Transactions on Neural Networks and Learning Systems, vol. 25, no.
5, pp. 920-933, 2014.
[7] R. S. Judson, and H. Rabitz, “Teaching lasers to control molecules,”
Physical Review Letters, vol. 68, pp. 1500-1503, 1992.
[8] N. Khaneja, T. Reiss, C Kehlet, T. Schulte-Herbrüggen, and S. J.
Glaser, “Optimal control of coupled spin dynamics: design of NMR
pulse sequences by gradient ascent algorithms,” Journal of Magnetic
Resonance, vol. 172, no. 2, pp. 296-305, 2005.
[9] R. Storn, and K. Price, “Differential evolution-a simple and efficient
heuristic for global optimization over continuous spaces,” Journal of
Global Optimization, vol. 11, no. 4, pp. 341-359, 1997.
[10] S. Das, and P. N. Suganthan, “Differential evolution: A survey of the
state-of-the-art,” IEEE Transactions on Evolutionary Computation, vol.
15, no. 1, pp. 4-31, 2011.

[11] E. Zahedinejad, J. Ghosh, and B. C. Sanders, “High-fidelity singleshot Toffoli gate via quantum control,” Physics Review Letters, vol.
114, no. 20, p. 200502, 2015.
[12] H. Ma, C. Chen, and D. Dong, “Differential evolution with equallymixed strategies for robust control of open quantum systems,” IEEE
International Conferernce on Systems, Man and Cybernetics, pp. 20552060, Hong Kong, October 9-12, 2015.
[13] E. Zahedinejad, J. Ghosh, and B. C. Sanders, “Desgining high-fidelity
single-shot three-qubit gates: A machine-learning approach,” Physics
Review Applied, vol. 6, p. 054005, 2016.
[14] Y. Sun, H. Ma, C. Wu, C. Chen, and D. Dong, “Ensemble control of
open quantum systems using differential evolution,” 2015 10th Asian
Control Conference (ASCC), pp. 1-6, 2015.
[15] E. Zahedinejad, S. Schirmer, and B. C. Sanders, “Evolutionary algorithms for hard quantum control,” Physics Review A, vol. 90, no. 3, p.
032310, 2014.
[16] M. R. James, H. I. Nurdin, and I. R. Petersen, “H ∞ control of
linear quantum stochastic systems,” IEEE Transactions on Automatic
Control, vol. 53, pp. 1787-1803, 2008.
[17] D. Dong and I. R. Petersen, “Sliding mode control of two-level
quantum systems,” Automatica, 48(5): 725-735, 2012.
[18] X. Xing, R. Rey-de-Castro and H. Rabitz, “Assessment of optimal
control mechanism complexity by experimental landscape Hessian
analysis: fragmentation of CH2 BrI,” New Journal of Physics, vol. 16,
p. 125004, 2014.
[19] C. Chen, D. Dong, R Long, and H. A. Rabitz, “Sampling-based learning control of inhomogeneous quantum ensembles,” Physics Review
A, vol. 89, no. 2, p. 023402, 2014.
[20] J. S. Li, and N. Khaneja, “Control of inhomogeneous quantum
ensembles,” Physics Review A, vol. 73, no. 3, p. 030302, 2006.
[21] J. S. Li, J. Ruths, T. Y. Yu, H. Arthanari, and G. Wagner, “Optimal
pulse design in quantum control: A unified computational method,”
Proceedings of the National Academy of Sciences of the USA, vol.
108, no. 5, pp. 1879-1884, 2011.
[22] L. M. Duan, M. D. Lukin, J. I. Cirac, and P. Zoller, “Long-distance
quantum communication with atomic ensembles and linear optics,”
Nature, vol. 414, no. 6862, pp. 413-418, 2001.
[23] K. Beauchard, P.S.P. da Silva, and P. Rouchon, “Stabilization for an
ensemble of half-spin systems,” Automatica, vol. 48, pp. 68-76, 2012.
[24] H. P. Breuer, and F. Petruccione, “The Theory of Open Quantum
Systems,” Oxford University Press, 2002.
[25] M. Mesbahi, and M. Egerstedt, “Graph theoretic methods in multiagent
networks,” Princeton University Press, 2010.
[26] L. Mazzarella, A. Sarlette, and F. Ticozzi, “Consensus for quantum
networks: from symmetry to gossip iterations,” IEEE Transactions on
Automatic Control, vol. 60, no. 1, pp. 158-172, 2015.
[27] G. Shi, D. Dong, I. R. Petersen, and K. H. Johansson, “Reaching
a quantum consensus: master equations that generate symmetrization
and synchronization,” IEEE Transactions on Automatic Control, vol.
61, no. 2, pp. 374-387, 2016.
[28] R. Sepulchre, A.Sarlette, and P. Rouchon, “Consensus in noncommutative spaces,” in Proc. 49th IEEE Conf. Decision Control,
Atlanta, GA, USA, Dec. 15-17, 2010, pp. 6596-6601.
[29] L. Mazzarella, F. Ticozzi, and A. Sarlette, “Extending robustness and
randomization from consensus to symmetrization algorithms,” SIAM
Journal on Control and Optimization, vol. 53, no. 4, pp. 2076-2099,
2015.
[30] F. Ticozzi, “Symmetrizing quantum dynamics beyond gossip-type
algorithms,” Automatica, vol. 74, pp. 38-46, 2016.
[31] D. Dong, M. A. Mabrok, I. R. Petersen, B. Qi, C. Chen, and H.
Rabitz, “Sampling-based learning control for quantum systems with
uncertainties,” IEEE Transactions on Control Systems Technology, vol.
23, pp. 2155-2166, 2015.
[32] C. Wu, B. Qi, C. Chen, and D. Dong, “Robust learning
control
design
for
quantum
unitary
transformations,”
IEEE Transactions on Cybernetics, in press, 2017, online
http://ieeexplore.ieee.org/document/7579636/.
[33] G. Lindblad, “On the generators of quantum dynamical semigroups,”
Communications in Mathmatical Physics, vol. 48, no. 2, pp. 119-130,
1976.
[34] F. Yang, S. Cong, R. Long, T. S. Ho, R. Wu, and H. Rabitz, “Exploring
the transition-probability-control landscape of open quantum systems:
Application to a two-level case,” Physics Review A, vol. 88, no. 3, p.
033420, 2013.

[35] G. Kimura, “The Bloch vector for N-level systems,” Physics Letters
A, vol. 314, pp. 339-349, 2003.
[36] H. Rabitz, M. M. Hsieh, and C. M. Rosenthal, “Quantum optimally
controlled transition landscapes,” Science, vol. 303, no. 5666, pp.
1998-2001, 2004.
[37] K. M. Tibbetts, X. Xing and H. Rabitz, “Laboratory transferability
of optimally shaped laser pulses for quantum control,” Journal of
Chemical Physics, vol. 140, p. 074302, 2014.
[38] N. M. Hamza, D. L. Essam, and R. A. Sarker, “Constraint consensus
mutation-based differential evolution for constrained optimization,”
IEEE Transactions on Evolutionary Computation, vol. 20, no. 3, pp.
447-459, 2016.
[39] R. Storn, and K. Price, “Differential evolution-a simple and efficient
adaptive scheme for global optimization over continuous spaces,”
ICSI, USA, Techonology Rep. TR-95-012, 1995 [Online]. Avaliable:
http://icsi.berkeley.edu/ storn/litera.html
[40] A. K. Qin, V. L. Huang, and P. N. Suganthan, “Differential evolution
algorithm with strategy adaptation for global numerical optimization,”
IEEE Transactions on Evolutionary Computation, vol. 13, no. 2, pp.
398-417, 2009.
[41] F. Neri, and V. Tirronen, “Recent advances in differential evolution: A
survey and experimental analysis,” Artificial Intelligence Review, vol.
33, no. 1-2, pp. 61-106, 2010.
[42] R. L. Becerra, and C. A. Coello Coello, “Cultured differential evolution for constrained optimization,” Computing Methods in Applied
Mechanics and Engineering, vol. 195, no. 33-36, pp. 4303-4322, 2006.
[43] R. Mallipeddi, P. N. Suganthan, Q. K. Pan, and M. F. Tasgetirenc,
“Differential evolution algorithm with ensemble of parameters and
mutation strategies,” Applied Soft Computing, vol. 11, no. 2, pp. 16791696, 2011.
[44] C. Chen, D. Dong, B. Qi, I. R. Petersen and H.
Rabitz,
“Quantum
ensemble
classification:
A samplingbased learning control approach,” IEEE Transactions on
Neural Networks and Learning Systems, in press, online
http://ieeexplore.ieee.org/stamp/stamp.jsp?arnumber=7439835.
[45] H. Jirari, and W. Pötz, “Optimal coherent control of dissipative N-level
systems,” Physics Review A, vol. 72, no. 1, p. 013409, 2005.
[46] J. Q. You, and F. Nori, “Superconducting circuits and quantum
information,” Physics Today, vol. 58, no. 11, pp. 42-47, 2005.
[47] Y. Makhlin, G. Scöhn, and A. Shnirman, “Josephson-junction qubits
with controlled couplings,” Nature, vol. 398, no. 6725, pp. 305-307,
1999.
[48] R. C. Bialczak, M. Ansmann, M. Hofheinz, M. Lenander, E. Lucero,
M. Neeley, A. D. O’Connell, D. Sank, H. Wang, M. Weides, J. Wenner,
T. Yamamoto, A. N. Cleland, and J. M. Martinis, “Fast tunable coupler
for superconducting qubits,” Physics Review Letters, vol. 106, no. 6,
p. 060501, 2011.
[49] J. Q. You, J. S. Tsai, and F. Nori, “Controllable manipulation and entanglement of macroscopic quantum states in coupled charge qubits,”
Physics Review B, vol. 68, no. 2, p. 024510, 2003.
[50] D. Dong, C. Chen, B. Qi, I. R. Petersen and F. Nori, “Robust manipulation of superconducting qubits in the presence of fluctuations,”
Scientific Reports, 5: 7873, 2015.
[51] D. Dong, C. Wu, C. Chen, B. Qi, I. R. Petersen and F. Nori,
“Learning robust pulses for generating universal quantum gates,”
Scientific Reports, 6: 36090, 2016.
[52] K. M. Tibbetts, X. Xing and H. Rabitz, “Systematic trends in photonic
reagent induced reactions in a homologous chemical family,” Journal
of Physical Chemistry A, vol. 117, pp. 8025-8215, 2013.
[53] K. M. Tibbetts, X. Xing and H. Rabitz, “Optimal control of molecular
fragmentation with homologous families of photonic reagents and
chemical substrates,” Physical Chemistry Chemical Physics, vol. 15,
pp. 18012-18022, 2013.
[54] P. Palittapongarnpim, P. Wittek, E. Zahedinejad, S. Vedaie, and B.
C. Sanders, “Learning in quantum control: High-dimensional global
optimization for noisy quantum dynamics”, learning, arxiv: https://arxiv.org/abs/1607.03428.
[55] H. Zhang, A. Zhou, S. Song, Q. Zhang, X. Z. Gao and J. Zhang, “A
self-organizing multiobjective evolutionary algorithm,” IEEE Transactions on Evolutionary Computation, vol. 20, no. 5, pp. 792-806, 2016.
[56] Y. Yuan, H. Xu, B. Wang, and X. Yao, “A new dominance relationbased evolutionary algorithm for many-objective optimization,” IEEE
Transactions on Evolutionary Computation, vol. 20, no. 1, pp. 16-37,
2016.

