Operator Methods of the Maximum Principle in Problems of Optimization of Quantum Systems

Abstract

In the class of optimal control problems for quantum systems, operator optimality conditions for control are constructed in the form of fixed-point problems in the control space. The equivalence of the obtained operator optimality conditions to the well-known Pontryagin maximum principle is shown. Based on the obtained operator forms of optimality conditions, new iterative methods for finding extreme equations satisfying the maximum principle are developed. A comparative analysis of the effectiveness of the proposed operator methods of the maximum principle with known methods is carried out on model examples of optimization of quantum systems.

1. Introduction

Mathematical formulations of topical problems related to the optimal control of quantum systems have been considered in the works of many researchers [1,2,3,4,5,6]. In the works of V.F. Krotov, V.I. Gurman, and of their followers [7,8,9], there are studied classes of controlled quantum systems described by ordinary differential controls linear in state and control with nonlinear optimality criteria. In [7], the main features of the selected class of problems are indicated. The first feature is the high dimension of the system state vector ($n\approx {10}^4-{10}^6$). The second feature is the absence of restrictions on the state, including terminal restrictions. The third feature is the use of a scalar control function characterizing the electric field. In this class, the search for an optimal solution based on standard necessary optimality conditions in the form of a boundary value problem of the maximum principle [10,11] causes significant difficulties due to the large dimension. In [7,8], the global Krotov method [12] was used as a tool for finding solutions to problems, which was compared in efficiency with the known gradient method.

In this paper, we consider and modify a new approach to finding a solution in the considered class of problems. This approach is based on the representation of optimality conditions for control in the form of operator problems about a fixed point in the space of admissible controls. This representation makes it possible to apply and modify the known methods of fixed points to find solutions to the considered problems related to the optimal control of quantum systems.

The new fixed-point approach has been used and developed for more than ten years for various classes of continuous, discrete, and discrete-continuous optimal control problems, including those involving terminal and phase constraints, mixed control functions and parameters, unfixed control process termination time, and other features.

In [13], the fixed-point approach is used to represent conditions for nonlocal improvement of control in a general class of nonlinear optimal control problems with control functions and parameters. In [14], the fixed-point approach for representing conditions for the nonlocal improvement of control is modified to the class of problems considered in [7,8,9]. The modification of the approach consists of taking into account the characteristic property of the singularity of the solutions of the problems under consideration, due to the above features.

The paper [15] describes the fixed-point approach for representing optimality conditions for control in a general class of nonlinear optimal control problems with control functions. In this paper, this approach for representing optimality conditions for control is modified and studied taking into account the characteristic property of the singularity of solutions in the considered class of optimization problems for quantum control systems.

2. Conditions for Optimality of Control

To illustrate the proposed fixed-point approach, we consider a model class of optimal control problems for quantum systems with a quadratic optimality criterion similar to the paper [14], in which new operator forms of optimality conditions have a relatively simple description:

$$\dot{x}\left(t\right)=(A+u\left(t\right)B)x\left(t\right),x\left(t_0\right)=x^0,u\left(t\right)\in U\subset R,t\in T=[t_0,t_1],$$

(1)

$$\mathsf{\Phi }\left(u\right)=\langle x\left(t_1\right),Lx\left(t_1\right)\rangle \to \underset{u\in V}{inf}.$$

(2)

The vector $x\left(t\right)=(x_1\left(t\right),...,x_n\left(t\right))$ describes the state of the system. The control $u\left(t\right)$, $t\in T$ is modeled by a piecewise continuous scalar function with values in a compact and convex set $U\subset R$. The set V denotes the corresponding set of admissible controls. The matrices A, B and L have real coefficients. The matrix L is symmetric. The initial state $x^0$ and the time interval T have fixed values.

The Pontryagin function with a conjugate variable $\psi$ in problem (1) and (2) has the form:

$$H(\psi ,x,u,t)=\langle \psi ,(A+uB)x\rangle ,\psi \in R^n.$$

The standard conjugate system is represented as:

$$\dot{\psi }\left(t\right)=-(A^T+u\left(t\right)B^T)\psi \left(t\right),t\in T,\psi \left(t_1\right)=-2Lx\left(t_1\right).$$

(3)

Let $v\in V$. Let us introduce the following notation:

-

$x(t,v)$, $t\in T$, the solution of the system (1) for $u\left(t\right)=v\left(t\right)$;

-

$\psi (t,v)$, $t\in T$, the solution of the standard conjugate system (3) for $x\left(t\right)=x(t,v)$, $u\left(t\right)=v\left(t\right)$.

Additionally, we will use the notation $P_Y$ for the projection operator on to a set $Y\subset R^k$ in the Euclidean norm:

$$P_Y\left(z\right)=arg\underset{y\in Y}{min}(\parallel y-z\parallel ),z\in R^k.$$

The projection operation is characterized by an important property that can be represented as an inequality:

$$\langle y-P_Y\left(z\right),z-P_Y\left(z\right)\rangle \le 0,y\in Y.$$

The known [10,11] necessary optimality conditions for an admissible control (maximum principle and differential maximum principle) in problems (1) and (2) are equivalent.

The condition of the maximum principle for control $u\in V$ can be represented in the form:

$$u\left(t\right)=arg\underset{w\in U}{max}\langle \psi (t,u),Bx(t,u)\rangle w,t\in T.$$

(4)

The condition of the differential maximum principle using the projection operation can be written as the following relation with the parameter $\alpha >0$:

$$u\left(t\right)=P_U(u\left(t\right)+\alpha \langle \psi (t,u),Bx(t,u)\rangle ),t\in T.$$

(5)

To fulfill the maximum principle condition (4), it suffices to check condition (5) for at least one $\alpha >0$. Conversely, condition (4) implies the fulfillment of condition (5) for all $\alpha >0$.

We define the mapping $u^{*}$ as follows:

$$u^{*}(\psi ,x)=arg\underset{w\in U}{max}\langle \psi ,Bx\rangle w,\psi \in R^n,x\in R^n.$$

We introduce a mapping $u^{\alpha }$ with a parameter $\alpha >0$ using the relation:

$$u^{\alpha }(\psi ,x,w)=P_U(w+\alpha \langle \psi ,Bx\rangle ),x\in R^n,\psi \in R^n,w\in U.$$

Using the introduced mappings, the maximum principle condition (4) can be written as:

$$u\left(t\right)=u^{*}(\psi (t,u),x(t,u)),t\in T.$$

The condition of the differential maximum principle (5) takes the following form:

$$u\left(t\right)=u^{\alpha }(\psi (t,u),x(t,u),u\left(t\right)),t\in T.$$

(6)

The well-known [10,11] approach to the search for extremal controls, i.e. satisfying the necessary optimality conditions, is the solution of the boundary value problem of the maximum principle. This problem in the considered classes, classes (1) and (2), takes the following form:

$$\dot{x}\left(t\right)=(A+u^{*}(\psi \left(t\right),x\left(t\right))B)x\left(t\right),x\left(t_0\right)=x^0,$$

(7)

$$\dot{\psi }\left(t\right)=(-A^T-u^{*}(\psi \left(t\right),x\left(t\right))B^T)\psi \left(t\right),\psi \left(t_1\right)=-2Lx\left(t_1\right).$$

(8)

Let the pair $\left(x\right(t),\psi (t\left)\right)$, $t\in T$, be a solution to the boundary value problems (7) and (8). Let us construct the output control $v\left(t\right)=u^{*}(\psi \left(t\right),x\left(t\right))$, $t\in T$. Then, by construction, we obtain the relations:

$$x\left(t\right)=x(t,v),\psi \left(t\right)=\psi (t,v),t\in T.$$

Consequently, the control $v\left(t\right)$ satisfies condition (4).

Conversely, let the control $v\in V$ be a solution to Equation (4). Let us form a pair of functions $\left(x\right(t,v),\psi (t,v\left)\right)$, $t\in T$. Then, by definition, these functions satisfy the boundary value problem (7) and (8).

Thus, the boundary value problems (7) and (8) are equivalent to the maximum principle condition (4).

Difficulties in solving the boundary value problems of the maximum principle, (7) and (8), in the general case are associated with the possible discontinuity and many meanings of the right-hand sides of the problem for the variables x, $\psi$. Even in the case of smoothness and uniqueness of the right-hand sides of the boundary value problem (7) and (8), its numerical solution by known methods (shooting method, linearization method, and finite difference method) [11] is computationally unstable due to the presence of positive real values of the eigenvalues of the corresponding Jacobian matrices.

In this paper, we consider a new approach to the search for extremal controls based on the transition from the boundary value problem of the maximum principle in the state space to equivalent operator problems on the fixed point of the maximum principle in the space of controls.

3. Operator Forms of the Maximum Principle

We define three mappings, X, $\mathsf{\Psi }$ and $V^{*}$, using the following relations:

$$X\left(v\right)=x,v\in V,x\left(t\right)=x(t,v),t\in T,$$

$$\mathsf{\Psi }\left(v\right)=\psi ,v\in V,\psi \left(t\right)=\psi (t,v),t\in T,$$

$$V^{*}(\psi ,x)=v^{*},\psi \in C\left(T\right),x\in C\left(T\right),v^{*}\left(t\right)=u^{*}(\psi \left(t\right),x\left(t\right)),t\in T,$$

where $C\left(T\right)$ is the space of functions continuous on T.

Using the above mappings, the maximum principle condition (4) can be represented as an operator equation in the form of a fixed-point problem in the control space:

$$v=V^{*}(\mathsf{\Psi }\left(v\right),X\left(v\right))=G_1^{*}\left(v\right),v\in V.$$

(9)

Construct new operator equations in the form of fixed-point problems equivalent to condition (4). Introduce the mapping $X^{*}$ as follows:

$$X^{*}\left(\psi \right)=x,\psi \in C\left(T\right),x\in C\left(T\right).$$

Here $x\left(t\right)$, $t\in T$, is the solution of the special Cauchy problem:

$$\dot{x}\left(t\right)=(A+u^{*}(\psi \left(t\right),x\left(t\right))B)x\left(t\right),x\left(t_0\right)=x^0.$$

Based on the introduced mapping $X^{*}$, consider the following operator equation:

$$v=V^{*}(\mathsf{\Psi }\left(v\right),X^{*}\left(\mathsf{\Psi }\left(v\right)\right))=G_2^{*}\left(v\right),v\in V.$$

(10)

We define the following mapping:

$${\mathsf{\Psi }}^{*}\left(x\right)=\psi ,x\in C\left(T\right),\psi \in C\left(T\right).$$

Here $\psi \left(t\right)$, $t\in T$ is the solution to the special Cauchy problem:

$$\dot{\psi }\left(t\right)=(-A^T-u^{*}(\psi \left(t\right),x\left(t\right))B^T)\psi \left(t\right),\psi \left(t_1\right)=-2Lx\left(t_1\right).$$

Consider the operator equation:

$$v=V^{*}({\mathsf{\Psi }}^{*}\left(X\left(v\right)\right),X\left(v\right))=G_3^{*}\left(v\right),v\in V.$$

(11)

Following the work [15], the operator equations, Equations (9)–(11), are equivalent to the set of admissible controls. Thus, we obtain the following statement:

Theorem 1.

Operator fixed-point problems (9)–(11) are equivalent to the boundary value problems of the maximum principle, (7) and (8).

The condition of the maximum principle in projection form (5) can also be represented in the form of equivalent operator equations on the set of admissible controls.

Introduce an additional operator $V^{\alpha }$, $\alpha >0$, by the relation:

$$V^{\alpha }(\psi ,x,v)=v^{\alpha },\psi \in C\left(T\right),x\in C\left(T\right),v\in V,$$

$$v^{\alpha }\left(t\right)=u^{\alpha }(\psi \left(t\right),x\left(t\right),v),t\in T.$$

Define the operator $X^{\alpha }$, $\alpha >0$:

$$X^{\alpha }(\psi ,v)=x^{\alpha },\psi \in C\left(T\right),v\in V,x^{\alpha }\left(t\right)=x^{\alpha }(t,\psi ,v),t\in T,$$

where $x^{\alpha }(t,\psi ,v),t\in T$, is the solution of the Cauchy problem:

$$\dot{x}\left(t\right)=(A+u^{\alpha }(\psi \left(t\right),x\left(t\right),v\left(t\right))B)x\left(t\right),x\left(t_0\right)=x^0.$$

Construct the operator ${\mathsf{\Psi }}^{\alpha },\alpha >0$:

$${\mathsf{\Psi }}^{\alpha }(x,v)={\psi }^{\alpha },x\in C\left(T\right),v\in V,{\psi }^{\alpha }\left(t\right)={\psi }^{\alpha }(t,x,v),$$

where ${\psi }^{\alpha }(t,x,v)$, $t\in T$, is the solution of the conjugate Cauchy problem:

$$\dot{\psi }\left(t\right)=(-A^T-u^{\alpha }(\psi \left(t\right),x\left(t\right),v\left(t\right))B^T)\psi \left(t\right),\psi \left(t_1\right)=-2Lx\left(t_1\right).$$

Consider three operator equations in the form of fixed-point problems:

$$v=V^{\alpha }(\mathsf{\Psi }\left(v\right),X\left(v\right),v)=G_1^{\alpha }\left(v\right),v\in V,\alpha >0,$$

(12)

$$v=V^{\alpha }(\mathsf{\Psi }\left(v\right),X^{\alpha }(\mathsf{\Psi }\left(v\right),v),v)=G_2^{\alpha }\left(v\right),v\in V,\alpha >0,$$

(13)

$$v=V^{\alpha }({\mathsf{\Psi }}^{\alpha }(X\left(v\right),v),X\left(v\right),v)=G_3^{\alpha }\left(v\right),v\in V,\alpha >0.$$

(14)

Similarly, following [15], these equations are equivalent to the set of admissible controls. Thus, the following statement holds:

Theorem 2.

Projection operator fixed-point problems (12)–(14) are equivalent to the boundary value problem of the maximum principle (7) and (8).

Let us note the following important features of the constructed projection problems on a fixed point.

1. Projection control operators, due to the properties of the projection operation, are continuous and satisfy the Lipschitz condition, in contrast to discontinuous and generally multivalued control operators based on the maximum operation in problems (9)–(11).

2. The search for extremal controls, which are solutions to the projection problems on a fixed point, (12)–(14), can be carried out for any given values of the projection parameter $\alpha >0$, including sufficiently small values.

These features of projection problems on a fixed point are essential factors for increasing the efficiency of the numerical search for extremal controls.

4. Operator Methods of the Maximum Principle

We consider the general fixed-point problem for an operator $G_E:V_E\to V_E$, acting on a set $V_E$ in a complete normalized space E with a norm ${\parallel ·\parallel }_E$,

$$v=G_E\left(v\right),v\in V_E.$$

To solve it numerically, one can apply the well-known simple iteration method with an index $k\ge 0$, which has the form:

$$v^{k+1}=G_E\left(v^k\right),v^0\in V_E.$$

The convergence of the iterative process can be analyzed using the well-known principle of compressive mappings [16].

Each operator equation from relations (9)–(14) can be considered as a fixed-point problem on the set of admissible controls in the following general form:

$$v=G\left(v\right),v\in V.$$

(15)

To solve the problem (15), an iterative process with the index $k\ge 0$ is proposed:

$$v^{k+1}=G\left(v^k\right),v\in V.$$

(16)

As an illustration of processes of the form (16), consider iterative processes for searching for extremal controls based on projection problems about a fixed point of the maximum principle (12)–(14), which, respectively, take the form with index $k\ge 0$:

$$v^{k+1}=V^{\alpha }(\mathsf{\Psi }\left(v^k\right),X\left(v^k\right),v^k)=G_1^{\alpha }\left(v^k\right),v^0\in V,\alpha >0,$$

(17)

$$v^{k+1}=V^{\alpha }(\mathsf{\Psi }\left(v^k\right),X^{\alpha }(\mathsf{\Psi }\left(v^k\right),v^k),v^k)=G_2^{\alpha }\left(v^k\right),v^0\in V,\alpha >0,$$

(18)

$$v^{k+1}=V^{\alpha }({\mathsf{\Psi }}^{\alpha }(X\left(v^k\right),v^k),X\left(v^k\right),v^k)=G_3^{\alpha }\left(v^k\right),v^0\in V,\alpha >0.$$

(19)

In the considered projection methods of the maximum principle, the projection parameter $\alpha >0$ is fixed in the iterative process of successive approximations of the control.

The complexity of implementing one iteration of each of the processes (17)–(19) is two Cauchy problems for phase and conjugate variables.

Indeed, this is obvious for the process (17). In this case, process (17) is written in the pointwise form as:

$$v^{k+1}\left(t\right)=u^{\alpha }(\psi (t,v^k),x(t,v^k),v^k\left(t\right),t),t\in T.$$

For process (18), at each k-th iteration with $k\ge 0$, we obtain the following.

After calculating the solution of the Cauchy problem $\psi (t,v^k)$, $t\in T$, we find the solution $x\left(t\right)$, $t\in T$ of the special Cauchy problem for the phase system:

$$\dot{x}\left(t\right)=(A+u^{\alpha }(\psi (t,v^k),x\left(t\right),v^k\left(t\right))B)x\left(t\right),x\left(t_0\right)=x^0.$$

Simultaneously, together with the solution of the Cauchy problem, we determine the output control according to the rule:

$$v^{k+1}\left(t\right)=u^{\alpha }(\psi (t,v^k),x\left(t\right),v^k\left(t\right)),t\in T.$$

Then, by construction, we obtain the relation:

$$x\left(t\right)=x(t,v^{k+1}),t\in T.$$

By of this equality, the iterative process (18) in pointwise form can be written in the following implicit form:

$$v^{k+1}\left(t\right)=u^{\alpha }(\psi (t,v^k),x(t,v^{k+1}),v^k\left(t\right)),t\in T.$$

(20)

Similarly, at the k-th iteration of the process (19), after calculating $x(t,v^k)$, $t\in T$, we find the solution $\psi \left(t\right)$, $t\in T$, of the special Cauchy problem for the conjugate system:

$$\dot{\psi }\left(t\right)=(-A^T-u^{\alpha }(\psi \left(t\right),x(t,v^k),v^k\left(t\right))B^T)\psi \left(t\right),\psi \left(t_1\right)=-2Lx(t_1,v^k).$$

Simultaneously the output control is constructed according to the rule:

$$v^{k+1}\left(t\right)=u^{\alpha }(\psi \left(t\right),x(t,v^k),v^k\left(t\right)),t\in T.$$

The theoretical conditions for the convergence of iterative processes (17)–(19) for sufficiently small projection parameters $\alpha >0$ can be substantiated similarly to [17] based on the formulation of requirements in problems (1) and (2), ensuring the application of the indicated principle of compressive mappings in the complete space of continuous controls or the extended complete space of measurable controls:

$$V\subset V_L=\{v\in L_{\infty }\left(T\right):v\left(t\right)\in U,t\in T\}$$

with the norm ${\parallel v\parallel }_{\infty }=ess\underset{t\in T}{sup}\parallel v\left(t\right)\parallel$, $v\in V_L$.

Unlike the standard gradient projection method, at each iteration of the proposed projection methods of the maximum principle, relaxation for the objective functional is not guaranteed.

In contrast to the global Krotov method, at each iteration of the proposed projection methods, Cauchy problems with a continuous and uniquely defined right-hand side are solved.

Let us single out other features of the proposed projection methods that are important for increasing their computational efficiency:

-

The non-locality of successive approximations of control, due to the fixed choice of the projection parameter $\alpha >0$;

-

The absence of the operation of varying control in the vicinity of the current approximation to provide improved control, which is characteristic of gradient methods;

-

The possibility of obtaining extreme controls for sufficiently small values of the projection parameter $\alpha >0$, which ensure the fundamental convergence of iterative processes.

Simple iteration methods for solving fixed-point problems (9)–(11) based on the maximization operation have a similar structure. In particular, the iterative process with the index $k\ge 0$ for searching for extremal controls based on the fixed-point problem (9) takes the form:

$$v^{k+1}=V^{*}(\mathsf{\Psi }\left(v^k\right),X\left(v^k\right))=G_1^{*}\left(v^k\right),v^0\in V,\alpha >0.$$

In pointwise form, this process is written as:

$$v^{k+1}\left(t\right)=u^{*}(\psi (t,v^k),x(t,v^k)),t\in T.$$

(21)

Note that method (21) is essentially equivalent to the well-known method of successive approximations of phase and conjugate variables [18] for solving the boundary value problems of the maximum principle, (7) and (8).

In contrast to the global Krotov method, at each iteration of the considered method (21), two simple Cauchy problems with a precomputed control are solved. In the Krotov method, at each iteration, in the general case, a special Cauchy problem with a discontinuous and multivalued right-hand side is solved.

5. Examples

The main feature of the considered class of optimal control problems for quantum systems is the property of singularity of extremal controls. This property is expressed in the existence of singular time intervals of non-zero measure for extremal controls, on which the derivative of the Pontryagin function becomes equal to zero:

$$H_u(\psi (t,u),x(t,u),u\left(t\right),t)=\langle \psi (t,u),Bx(t,u)\rangle =0.$$

On such singular time intervals, it becomes impossible to determine the values of the extremal control from the condition of the maximum principle (4).

The proposed operator methods of the maximum principle, taking into account the indicated property of singularity, are modified in specific examples under consideration and compared in terms of computational efficiency with known methods.

The computational implementation of the proposed methods of the maximum principle is characterized by the following features.

The numerical solution of phase and conjugate Cauchy problems was performed by the Runge–Kutta–Werner method of variable (5–6) order of accuracy using the DIVPRK program of the IMSL Fortran PowerStation 4.0 library [19]. The values of the controlled, phase, and conjugate variables were stored in the nodes of a fixed uniform grid $T_h$ with a sampling step $h>0$ on the interval T. In the intervals between neighboring grid nodes $T_h$, the control value was assumed to be constant and equal to the value in the left node. The numerical calculation of the fixed-point problem was carried out before the condition was fulfilled:

$$max\left\{\right|v^{k+1}\left(t\right)-v^k\left(t\right)|,t\in T_h\}\le {\epsilon }_m,$$

in which ${\epsilon }_m>0$ is the given accuracy of calculating the fixed-point problem.

Example 1.

(projection methods (17) and (18)).

The well-known model problem of control of the system of spins of quantum particles is considered [20], which can be represented in the following form:

$$\mathsf{\Phi }\left(u\right)=1-\langle x\left(t_1\right),Lx\left(t_1\right)\rangle \to inf,$$

$$L=\left(\begin{array}{cccc}{a}_1^2+b_1^2& a_1{a}_2+b_1{b}_2& 0& a_1{b}_2-b_1{a}_2\\ a_1{a}_2+b_1{b}_2& a_2^2+b_2^2& a_2{b}_1-b_2{a}_1& 0\\ 0& a_2{b}_1-b_2{a}_1& b_1^2+a_1^2& b_1{b}_2+a_1{a}_2\\ a_1{b}_2-b_1{a}_2& 0& b_1{b}_2+a_1{a}_2& b_2^2+a_2^2\end{array}\right),$$

$${\dot{x}}_1\left(t\right)=u\left(t\right)x_3\left(t\right)+x_4\left(t\right),{\dot{x}}_2\left(t\right)=x_3\left(t\right)-u\left(t\right)x_4\left(t\right),$$

$${\dot{x}}_3\left(t\right)=-u\left(t\right)x_1\left(t\right)-x_2\left(t\right),{\dot{x}}_4\left(t\right)=-x_1\left(t\right)+u\left(t\right)x_2\left(t\right),$$

$$x_1\left(0\right)=\frac{1}{\sqrt{2}},x_2\left(0\right)=\frac{1}{\sqrt{2}},x_3\left(0\right)=0,x_4\left(0\right)=0,t\in T=[0,t_1],t_1=1.5,$$

$$a_1=0.6,b_1=-0.3,a_2=0.1,b_2=\sqrt{0.54}.$$

The vector $x\left(t\right)$ describes the state of the quantum system, the function $u\left(t\right)$ characterizes the effect of an external field, $u\left(t\right)\in U=[-30,30]$, $t\in T$.

In [20], to calculate the optimal control problem under consideration, the global Krotov method was used, the efficiency of which was compared with the well-known gradient method. The control determined from physical considerations was chosen as the initial control approximation for the specified iterative methods:

$$u\left(t\right)=tg\left(2\gamma (2t-1.5)\right),t\in T,\gamma =-\frac{1}{3}arctg(-30).$$

The Pontryagin function in the problem has the form:

$$H(\psi ,x,u,t)={\psi }_1(u{x}_3+x_4)+{\psi }_2(x_3-u{x}_4)+{\psi }_3(-u{x}_1-x_2)+{\psi }_4(-x_1+u{x}_2).$$

The standard conjugate system is written as:

$${\dot{\psi }}_1\left(t\right)=u\left(t\right){\psi }_3\left(t\right)+{\psi }_4\left(t\right),{\dot{\psi }}_2\left(t\right)={\psi }_3\left(t\right)-u\left(t\right){\psi }_4\left(t\right),t\in T,$$

$${\dot{\psi }}_3\left(t\right)=-u\left(t\right){\psi }_1\left(t\right)-{\psi }_2\left(t\right),{\dot{\psi }}_4\left(t\right)=u\left(t\right){\psi }_2\left(t\right)-{\psi }_1\left(t\right),t\in T,$$

$${\psi }_1\left(t_1\right)=2(a_1^2+b_1^2)x_1\left(t_1\right)+2(a_1{a}_2+b_1{b}_2)x_2\left(t_1\right)+2(a_1{b}_2-b_1{a}_2)x_4\left(t_1\right),$$

$${\psi }_2\left(t_1\right)=2(a_1{a}_2+b_1{b}_2)x_1\left(t_1\right)+2(a_2^2+b_2^2)x_2\left(t_1\right)+2(a_2{b}_1-b_2{a}_1)x_3\left(t_1\right),$$

$${\psi }_3\left(t_1\right)=2(a_2{b}_1-b_2{a}_1)x_2\left(t_1\right)+2(b_1^2+a_1^2)x_3\left(t_1\right)+2(b_1{b}_2+a_1{a}_2)x_4\left(t_1\right),$$

$${\psi }_4\left(t_1\right)=2(a_1{b}_2-b_1{a}_2)x_1\left(t_1\right)+2(b_1{b}_2+a_1{a}_2)x_3\left(t_1\right)+2(b_2^2+a_2^2)x_4\left(t_1\right).$$

The fixed-point projection problems (12) and (13) have the same pointwise form:

$$v\left(t\right)=P_U(v\left(t\right)+\alpha ({\psi }_1(t,v)x_3(t,v)-{\psi }_2(t,v)x_4(t,v)-{\psi }_3(t,v)x_1(t,v)+{\psi }_4(t,v)x_2(t,v))).$$

The explicit iterative method of the maximum principle (17) for solving this fixed-point problem at index $k\ge 0$ has a pointwise form:

$$\begin{array}{c}{v}^{k+1}\left(t\right)=P_U(v^k\left(t\right)+\alpha ({\psi }_1(t,v^k)x_3(t,v^k)-{\psi }_2(t,v^k)x_4(t,v^k)\hfill \\ \hfill -{\psi }_3(t,v^k)x_1(t,v^k)+{\psi }_4(t,v^k)x_2(t,v^k)\left)\right).\end{array}$$

Accordingly, the implicit iterative method (20) at index $k\ge 0$ takes the form:

$$\begin{array}{c}{v}^{k+1}\left(t\right)=P_U(v^k\left(t\right)+\alpha ({\psi }_1(t,v^k)x_3(t,v^{k+1})-{\psi }_2(t,v^k)x_4(t,v^{k+1})\hfill \\ \hfill -{\psi }_3(t,v^k)x_1(t,v^{k+1})+{\psi }_4(t,v^k)x_2(t,v^{k+1})\left)\right).\end{array}$$

The calculation was carried out on a sampling grid with a step $h={10}^{-5}$ and the criterion for stopping the calculation ${\epsilon }_m={10}^{-3}$.

Table 1. The comparative results of the first four iterations.

Number s	$\mathbf{\Phi }\left({\mathit{u}}^{\mathit{s}}\right)$ GlM	$\mathbf{\Phi }\left({\mathit{u}}^{\mathit{s}}\right)$ GrM	$\mathbf{\Phi }\left({\mathit{v}}^{\mathit{s}}\right)$ PPM1	$\mathbf{\Phi }\left({\mathit{v}}^{\mathit{s}}\right)$ PPM2
0	0.7681	0.7681	0.7680	0.7680
1	0.1401	0.6911	0.6705	0.6610
2	0.0040	0.6107	0.4049	0.3881
3	0.0021	0.5421	0.2404	0.2175
4	0.0015	0.4913	0.1718	0.1191

shows the comparative results of the first four iterations of improving the objective functional with an index $s\ge 0$, starting from the starting control $v^0=u$ specified in [20]. We compare the results of the explicit projection method (PPM1) and the implicit projection method (PPM2) for the parameter $\alpha ={10}^{-2}$ with known [20] calculation data by the global method ( GlM) and the gradient method (GRM).

Figure 1 shows the final computational control $v1\left(t\right)$, $t\in T$ obtained by the PPM1 method with the number of control improvement iterations equal to 14, and the value of the functional ${\mathsf{\Phi }}^{*}\approx 0.001421$.

Figure 2 shows the final computational control $v2\left(t\right)$, $t\in T$, obtained by the PPM2 method with the number of control improvement iterations equal to 26 and the value of the functional ${\mathsf{\Phi }}^{*}\approx 0.000704$.

In [20], the final calculated value of the functional ${\mathsf{\Phi }}^{*}\approx 0.000952$ obtained by the global Krotov method at the ninth iteration of control improvement is indicated.

Thus, within the framework of the considered example, the proposed projection methods of the maximum principle allow one to achieve similar results in terms of the value of the functional with the global Krotov method. The methods differ in the number of calculated control improvement iterations. However, at each iteration of the global method, it is necessary to solve the complex Cauchy problem for phase variables with a special right-hand side based on a multivalued and discontinuous operation to the maximum for calculating control values. In the proposed fixed-point method, at each iteration, much simpler Cauchy problems with a uniquely defined and continuous right-hand side are solved, which makes the considered projection methods of the maximum principle much easier to implement than the global Krotov method. Due to the simplicity of implementation and easy adjustment of the convergence, controlled by the choice of the projection parameter $\alpha >0$, these methods can be successfully used to obtain practical initial approximations for subsequent refinement by other iterative methods for solving optimal control problems of the class under consideration.

Example 2.

(method (21)).

To illustrate the work of the maximum principle method based on the maximization operation (21), the problem from the previous example is considered.

The corresponding fixed-point problem of the maximum principle (9) has the following form:

$$v\left(t\right)=u^{*}(\psi (t,v),x(t,v)),t\in T,$$

$$u^{*}(\psi ,x)=\left\{\begin{array}{c}+30,g(\psi ,x)>0,\hfill \\ -30,g(\psi ,x)<0,\hfill \\ w\in U,g(\psi ,x)=0,\hfill \end{array}\right.$$

where $g(\psi ,x)={\psi }_1{x}_3-{\psi }_2{x}_4-{\psi }_3{x}_1+{\psi }_4{x}_2$. The iterative process (21) for solving this fixed-point problem at $k\ge 0$, respectively, takes the form:

$$v^{k+1}\left(t\right)=u^{*}(\psi (t,v^k),x(t,v^k)),t\in T,$$

with the above switching function $g(\psi ,x)$.

In the case of the existence of a time interval $[{\mathsf{\Theta }}_1,{\mathsf{\Theta }}_2]\subset T$ of a non-zero measure, where $g\left(\psi \right(t),x(t\left)\right)=0$, $t\in [{\mathsf{\Theta }}_1,{\mathsf{\Theta }}_2]$, the control $u^{*}$ is called singular on this interval. Singular controls are determined by the sequential differentiation by an argument $t\in T$ of the identity $g\left(\psi \right(t),x(t\left)\right)=0$ taking into account the phase and conjugate systems. In practical calculations, similarly to the work [20], the equality of the switching function $g(\psi ,x)$, to zero, which determines a singular mode, is understood in the sense of belonging to some small ε, neighborhood of zero, where $\epsilon >0$. Thus, we obtain the following practical calculation formula for the simple iteration method:

$$v^{k+1}\left(t\right)=\left\{\begin{array}{c}+30,g(\psi (t,v^k),x(t,v^k))>\epsilon ,\hfill \\ -30,g(\psi (t,v^k),x(t,v^k))<-\epsilon ,\hfill \\ w\in U,\left|g\right(\psi (t,v^k),x(t,v^k)\left)\right|\le \epsilon .\hfill \end{array}\right.$$

If at the time t the condition is satisfied:

$$\left|g\right(\psi (t,v^k),x(t,v^k)\left)\right|\le \epsilon ,$$

the value $w\in U$ is determined by the following rule.

The value is calculated:

$$g_k=g(\psi (t+\delta ,v^k),x(t+\delta ,v^k)),$$

for a given $\delta >0$.

If $g_k>\epsilon$, then $w=30$.

If $g_k<-\epsilon$, then $w=-30$.

If $|g_k|\le \epsilon$, then the value $w\in U$ is determined by the special control calculation rule as follows.

The value is calculated as:

$$\begin{array}{c}\hfill a_k=-{\psi }_1(t,v^k)x_4(t,v^k)-{\psi }_2(t,v^k)x_3(t,v^k)+{\psi }_3(t,v^k)x_2(t,v^k)+{\psi }_4(t,v^k)x_1(t,v^k).\end{array}$$

If $|a_k|>\epsilon$, then

1. 

The value is calculated as $c_k=\frac{b_k}{a_k}$, where:

$$\begin{array}{c}\hfill b_k=-{\psi }_1(t,v^k)x_3(t,v^k)+{\psi }_2(t,v^k)x_4(t,v^k)+{\psi }_3(t,v^k)x_1(t,v^k)-{\psi }_4(t,v^k)x_2(t,v^k).\end{array}$$

If $|c_k|\le 30$, then the value $w=c_k$.

If $|c_k|>30$, then go to step 2.

If $a_k\le \epsilon$, then:

2. 

The value is calculated as

$$\begin{array}{c}\hfill d_k={\psi }_1(t,v^k)x_2(t,v^k)-{\psi }_2(t,v^k)x_1(t,v^k)+{\psi }_3(t,v^k)x_4(t,v^k)-{\psi }_4(t,v^k)x_3(t,v^k).\end{array}$$

If $|d_k|>\epsilon$, then the value $w=0$.

If $|d_k|\le \epsilon$, then the value $w\in U$ is chosen randomly from the interval U.

In the numerical implementation of the algorithm, the value $\delta >0$ chosen equal to the grid step $h>0$.

Table 2. The comparative results of the calculation by the considered maximum principle method for the first four iterations.

Number s	$\mathbf{\Phi }\left({\mathit{u}}^{\mathit{s}}\right)$ GlM	$\mathbf{\Phi }\left({\mathit{u}}^{\mathit{s}}\right)$ GrM	$\mathbf{\Phi }\left({\mathit{v}}^{\mathit{s}}\right)$ MPM
0	0.7681	0.7681	0.7680
1	0.1401	0.6911	0.5714
2	0.0040	0.6107	0.3612
3	0.0021	0.5421	0.1904
4	0.0015	0.4913	0.1380

shows the comparative results of the calculation by the considered maximum principle method (MPM) for the first four iterations of improving the functional with an index $s\ge 0$, starting from the above starting control, with the known [20] calculation data by the global method (GlM) and the gradient method (GrM). For an adequate comparison of the methods, the $\epsilon$-neighborhood of zero was determined by the value $\epsilon =0.001$, specified in [20].

In [20], the final calculated value of the functional ${\mathsf{\Phi }}^{*}\approx 0.000952$, obtained by the global method at the ninth iteration of control improvement, is indicated. In this case, a singular section of the final control, determined according to the rules [20], is the interval $[0.0667,t_1]$.

Figure 3 shows the final computational control $v3\left(t\right)$, $t\in T$, obtained by the MPM method, with the achieved value of the functional ${\mathsf{\Phi }}^{*}\approx 0.000989$ and the number of control improvement iterations equal to 18. A singular section of the final control with the discretization grid accuracy is $[0.0693,1.4717]$.

Figure 4 shows the final computational control $v4\left(t\right)$, $t\in T$, obtained by the MPM method from the initial approximation $v^0=u1$, which was obtained by the PPM1 method in example 1. In this case, the value of the functional is ${\mathsf{\Phi }}^{*}\approx 0.000907$ with the number of control improvement iterations equal to 7. A singular section of the final control is approximately equal to $[0.0698,1.4609]$.

Figure 5 shows the final computational control $v5\left(t\right)$, $t\in T$, obtained by the MPM method from the initial approximation $v^0=u2$, which was obtained by the PPM2 method in example 1. In this case, the value of the functional is ${\mathsf{\Phi }}^{*}\approx 0.000620$, with the number of control improvement iterations equal to 3. A singular section of the final control is approximately equal to $[0.0751,1.4512]$.

The calculations performed within the framework of the model problem show a high quantitative and qualitative efficiency of the implicit projection method of the maximum principle (20), which makes it possible to accurately calculate complex singular sections of extreme controls, which are typical in optimal control problems for quantum systems of the class under consideration. The main feature of this method, which is important for increasing efficiency, is the solution at each iteration of the Cauchy problems with a special uniquely defined and continuous right-hand side, in contrast to the global Krotov method.

6. Conclusions

In the considered class of optimal control problems for quantum systems, new operator forms of the maximum principle are proposed in the form of fixed-point problems in the control space, which make it possible to effectively apply and modify the well-known apparatus of the theory and methods of fixed points for constructing iterative algorithms to find extremal controls.

The developed iterative operator methods for searching for extremal controls are characterized by the following properties:

• computational stability, in contrast to standard methods for solving the boundary value problem of the maximum principle;
• nonlocality of successive control approximations;
• the absence of a laborious procedure of needle or convex variation of the control in a small neighborhood of the considered approximation, which is typical for gradient methods;
• the numerical solution of the Cauchy problems with a continuous and uniquely defined right-hand side at each iteration of the constructed projection methods, in contrast to the well-known global Krotov method.

The indicated properties of the proposed methods for searching for extremal controls are important factors for increasing the efficiency of the numerical solution of optimal control problems for quantum systems of the class under consideration.

In quantum systems with multidimensional control, the structures of the proposed operator methods of the maximum principle and the well-known global Krotov method remain the same, but the advantage of the indicated properties of the proposed projection methods of the maximum principle increases significantly in comparison with the global method.