Latest revision as of 00:54, 24 December 2022

Abstract

In this paper we consider a system of three nonlinear ordinary differential equations that model a three Josephson Junctions Array (JJA). Our goal is to stabilize the system around an unstable equilibrium employing an optimal control approach. We first define the cost functional and calculate its differential by using perturbation analysis and variational calculus. For the computational solution of the optimality system we consider a conjugate gradient algorithm for quadratic functionals in Hilbert spaces, combined with a finite difference discretization of the involucrated differential equations.

Key word: Josephson Junction, Cryogenic Memory, Conjugate Gradient Algorithms, Stabilization, Optimal Control.

MSC[2010] 65K10, 49M25.

1 Introduction

In this article we discuss the numerical simulation of a control approach to stabilize a Josephson Junction Array (JJA) around an unstable steady state. Our methodology relies significantly on a conjugate gradient algorithm operating in a well-chosen control space. It has been shown recently that such an array can carry out the functions of a memory cell operations ([1,2,3,4,5]), consequently this array is frequently referred as a Josephson junction array memory (JJAM).

A Josephson junction is a quantum interference device that consists of two super-conductors coupled by a weak link that may be, for example, an insulator or a ferromagnetic material. A Josephson junction can carry current without resistance (this current is called supercurrent) and such a device is an example of a macroscopic quantum phenomenon. In 1962, B.D. Josephson proposed the equations governing the dynamics of the Josephson junction effect, namely:

where ${\textstyle \hbar =h/2\pi }$, ${\textstyle h}$ being the Planck constant, ${\textstyle -e}$ is the electric charge of the electron, ${\textstyle V}$ and ${\textstyle I}$ are the voltage and current across the junction, respectively, ${\textstyle \phi }$ is the phase difference across the junction, and, finally, the constant ${\textstyle I_{c}}$ is the critical current across the junction. The Josephson junction technology offers numerous applications and could potentially provide sound alternatives for computer memory devices. In a 2005 US National Security Agency (NSA) report it has been alluded that, as transistors were rapidly approaching their physical limits, their most likely successors would be cryogenic devices based on Josephson junctions [6]. Indeed, (cf. [7]), ``Single flux quantum (SFQ) circuits produce small current pulses that travel at about ${\textstyle 1/3}$ the speed of light. Superconducting passive transmission lines (PTL) are also able to transmit the pulses with extremely low losses". These are currently the fastest switching digital circuits, since (cf. [8]) ``Josephson junctions, the superconducting switching devices, switch quickly (${\textstyle \sim 1}$ ps), dissipate very little energy per switch (${\textstyle <10^{-19}J}$), and communicate information via current pulses that propagate over superconducting transmission lines nearly without loss". Recently it was suggested that a small array consisting of inductively coupled Josephson junctions possesses all the basic functions (WRITE, READ, RESET) of a memory cell ([1], [3] and the references therein). In the two above articles, stable zero-voltage states were identified and basic memory cell operations based on the transitions between the equilibrium states were identified.

The equations modeling the dynamics of the inductively coupled Josephson junction array can be written in the following dimensionless form (${\textstyle t=1}$ corresponds to ${\textstyle 4.15ps}$)

where ${\textstyle I_{j}=i_{j}+ad_{j}}$, ${\textstyle i_{j}}$ being a direct current and ${\textstyle ad_{j}}$ an additional energy as a current pulse. Plausible values for the various quantities in the model are:

For those regimes where the second order derivatives can be neglected (meaning the junctions are (relatively) highly dissipative), (2) reduces to

The Read/Write operations can be performed by applying appropriate Gaussian pulses ([1], [4], [5] ) in order to drive the system from an equilibrium configuration to another one. In [9], we went beyond Gaussian pulses and investigated an optimal controllability approach in ${\textstyle L^{2}}$ to perform these transitions between equilibrium configurations. These approach is closely related to the methodology discussed in [10] for systems modelled by partial differential equations.

For the driving currents and coupling parameters provided in (3) and (4), the ordinary differential equation systems (5) (and (2) with ${\textstyle {\dfrac {d{\boldsymbol {\phi }}}{dt}}(0)={\boldsymbol {0}}}$) have several steady state solutions, consisting of phase triplets. In many cases, the steady state phases are near the values ${\textstyle 2\pi n_{k}}$, for certain integers ${\textstyle n_{k}}$, and, according to Braiman, et. al., ([1]), equilibrium junction phases can be defined by their offsets ${\textstyle \sigma _{k}}$ from the negative cosine function's minima, as shown in the following equation

Following Braiman et al., ([1]), Table 1 includes all the possible triplets of stable steady state offsets, where ${\textstyle 2\geq n_{1}\geq n_{2}\geq n_{3}}$. For the steady states, we assumed, without loss of generality, that ${\textstyle n_{3}=0}$, since all the phases can be shifted by the same multiple of ${\textstyle 2\pi }$. Also we include in Table 1 three unstable steady states, namely: ${\textstyle \{n_{1},n_{2},n_{3}\}=\{1,0,0(u)\},}$ ${\textstyle \{2,1,0(u)\},}$ ${\textstyle \{2,2,0(u)\}}$. Here the description of some unstable equilibria is artificial (in the sense that the terms ${\textstyle \sigma _{i}/2\pi }$ are not small and its only purpose is to be consistent with the description of stable states).

For the set of dc currents and coupling parameters defined by (3) and (4), no stable steady state exists when the first junction phase is about ${\textstyle 2\pi }$ (respectively ${\textstyle 4\pi }$) greater than both of the other two junction phases (see row where ${\textstyle \{n_{1},n_{2},n_{3}\}=\{1,0,0\}}$ (respectively ${\textstyle \{n_{1},n_{2},n_{3}\}=\{2,0,0\}}$)). Following Braiman, et. al., [1], for memory cell demonstration it is sufficient to manipulate the system within the particular sets of states, namely states ${\textstyle \{0,0,0\}}$ and ${\textstyle \{1,1,0\}}$. The circuit of three coupled Rapid Single Joint (RSJ) unit circuits associated with model (5) can operate as a memory cell if a set of operators can transition the system to clearly defined states and can output a signal that discriminates memory states. The value ${\textstyle n_{k}}$ as presented in equation (6), will describe the location of the ${\textstyle k^{th}}$ junction phase. When all three junction phases are in the same sinusoidal potential well, the system will be considered in the `0' state, ${\textstyle \{n_{1},n_{2},n_{3}\}=\{0,0,0\}.}$ When the phases of the first and second junctions are shifted to the next potential well (about ${\textstyle 2\pi }$ greater), the cell will be in the `1' state, ${\textstyle \{n_{1},n_{2},n_{3}\}=\{1,1,0\}}$. These two states correspond to the first and third rows of Table 1. Figure 1 shows the phases of the junctions when the systems starts from zero initial conditions. The junctions settle into the steady state `0' where all phases are close to the same multiple of ${\textstyle 2\pi }$. For convenience, in several of the following figures the phases are normalized by ${\textstyle 2\pi }$ and in all the graphs that include different colors, color blue, red and green is associated with junction 1, 2 and 3, respectively.

Figure 1: Time series of the phases scaled by ${\displaystyle 2\pi }$, with zero initial conditions.

In this work we investigate a controllability approach in order to stabilize the phase junctions of (5) around an unstable equilibrium. Figures 2 and 3 show the behavior of the system when the initial conditions are approximations of the unstable equilibrium ${\textstyle \{n_{1},n_{2},n_{3}\}}$=${\textstyle \{1,0,0(u)\}}$, ${\textstyle \{2,1,0(u)\}}$, ${\textstyle \{2,2,0(u)\}}$.

Figure 2: Evolution to state ${\displaystyle \{1,1,0(s)\}}$(left), ${\displaystyle \{2,1,0(s)\}}$(right) from an approximation to the unstable equilibrium ${\displaystyle \{1,0,0(u)\}}$(left), ${\displaystyle \{2,1,0(u)\}}$(right).

Figure 3: Evolution to ${\displaystyle \{2,2,0(s)\}}$ from an approximation to the unstable equilibrium ${\displaystyle \{2,2,0(u)\}}$.

The organization of this paper is as follows. In Section 2 we focus on a methodology to stabilize the JJAM system around an unstable equilibrium configuration, via a controllability approach. Section 3 is concerning the practical aspects of this methodology. In Section 4 we show the numerical results and in Section 5 we give the conclusions.

2 Formulation of the optimal control problem for the stabilization of the JJAM and the conjugate gradient algorithm

2.1 The approach

As we explained in the Introduction, Figures 2 and 3 show the behavior of the system (5) when the initial conditions are approximations of the unstable equilibrium configurations ${\textstyle \{n_{1},n_{2},n_{3}\}}$= ${\textstyle \{1,0,0(u)\}}$, ${\textstyle \{2,1,0(u)\}}$, ${\textstyle \{2,2,0(u)\}}$. Our goal in this subsection is to stabilize the system (5), around an unstable equilibrium, controlling it via one, two or three junctions i.e., we want to maintain the phase values in Figures 2 and 3 almost constant along time. This constant value must be the initial value on each case. The approach taken here is the following classical one, namely,

(a) Linearize the model (5) in the neighborhood of an (unstable) equilibrium of the system.

(b) Compute an optimal control for the linearized model.

(c) Apply the above control to the nonlinear system.

Let us consider an unstable state ${\textstyle {\boldsymbol {\theta }}=\{\theta _{1},\theta _{2},\theta _{3}\}}$ of (5) and a small initial variation ${\textstyle \delta {\boldsymbol {\theta }}}$ of ${\textstyle {\boldsymbol {\theta }}}$. The perturbation ${\textstyle \delta {\boldsymbol {\phi }}}$ of the steady-state ${\textstyle {\boldsymbol {\theta }}}$ in ${\textstyle (0,T)}$ satisfies approximately the following linear model

At least one of the eigenvalues of the jacobian of system (7) is positive. This means that the system can develop blow up phenomena (in infinite time). Figures 4, 5 and 6 show the solution of (7) for the three unstable equilibriums given by ${\textstyle \{n_{1},n_{2},n_{3}\}=\{1,0,0(u)\},\{2,1,0(u)\},\{2,2,0(u)\}}$ in Table 1, respectively. The used value for the initial perturbation was ${\textstyle \delta {\boldsymbol {\theta }}=[1e-2,1e-2,1e-2]}$.

Figure 4: Solution ${\displaystyle \mathbf {y} =\delta {\boldsymbol {\phi }}}$ of (7) for the unstable equilibrium given by ${\displaystyle \{n_{1},n_{2},n_{3}\}=\{1,0,0(u)\}}$ and ${\displaystyle \delta {\boldsymbol {\theta }}=[1e-2,1e-2,1e-2]}$.

Figure 5: Solution ${\displaystyle \mathbf {y} =\delta {\boldsymbol {\phi }}}$ of (7) for the unstable equilibrium given by ${\displaystyle \{n_{1},n_{2},n_{3}\}=\{2,1,0(u)\}}$ and ${\displaystyle \delta {\boldsymbol {\theta }}=[1e-2,1e-2,1e-2]}$.

Figure 6: Solution ${\displaystyle \mathbf {y} =\delta {\boldsymbol {\phi }}}$ of (7) for the unstable equilibrium given by ${\displaystyle \{n_{1},n_{2},n_{3}\}=\{2,2,0(u)\}}$ and ${\displaystyle \delta {\boldsymbol {\theta }}=[1e-2,1e-2,1e-2]}$.

It is clear that the model (7) is no longer valid if ${\textstyle \delta {\boldsymbol {\theta }}}$ becomes too large but the idea behind considering the linearized model (7) is to use it to compute a control action preventing ${\textstyle \delta {\boldsymbol {\phi }}}$ from becoming too large (and possibly driving ${\textstyle \delta {\boldsymbol {\phi }}}$ to zero) and hope that the computed control will also stabilize the original nonlinear system (5) with initial condition

2.2 The Linear Control Problem

Using the notation ${\textstyle \mathbf {y} }$ ${\textstyle =}$ ${\textstyle \delta {\boldsymbol {\phi }}}$, we shall (try to) stabilize the controlled variant of system (7) (the three junctions are used to control):

via the following control formulation:

where

and ${\textstyle ||\mathbf {y} ||^{2}=\left\vert y_{1}\right\vert ^{2}+\left\vert y_{2}\right\vert ^{2}+\left\vert y_{3}\right\vert ^{2}}$.

2.3 Optimality Conditions and Conjugate Gradient Solution for Problem (10)

2.3.1 Generalities and Synopsis

Let us denote by ${\textstyle DJ(\mathbf {v} )}$ the differential of ${\textstyle J}$ at ${\textstyle \mathbf {v} \in {\mathcal {U}}=L^{2}(0,T;3)=(L^{2}(0,T))^{3}}$. Since ${\textstyle {\mathcal {U}}}$ is a Hilbert Space for the scalar product defined by

the action ${\textstyle \left\langle DJ(\mathbf {v} ),\mathbf {w} \right\rangle }$ of ${\textstyle DJ(\mathbf {v} )}$ on ${\textstyle \mathbf {w} \in {\mathcal {U}}}$ can also be written as

with ${\textstyle DJ(\mathbf {v} )\mathbf {(} t\mathbf {)} \in {\mathcal {U}}}$. If ${\textstyle \mathbf {u} }$ is the solution of problem (10), it is characterized [from convexity arguments (see, e.g., [11])] by

Since the cost function ${\textstyle J}$ is quadratic and since the state model (7) is linear, ${\textstyle DJ(\mathbf {v} )}$ is in fact an affine function of ${\textstyle \mathbf {v} }$, implying in turn from (12) that ${\textstyle \mathbf {u} }$ is the solution of a linear equation in the control space ${\textstyle {\mathcal {U}}}$. In abstract form, equation (12) can be written as

and from the properties (need to be shown) of the operator ${\textstyle \mathbf {v} \rightarrow DJ(\mathbf {v} )-DJ(\mathbf {0} )}$, problem ${\textstyle DJ(\mathbf {u} )=\mathbf {0} }$ could be solved by a (quadratic case)-conjugate gradient algorithm operating in the space ${\textstyle {\mathcal {U}}}$. The practical implementation of the above algorithm would requires the explicit knowledge of ${\textstyle DJ(\mathbf {v} )}$.

2.3.2 Computing DJ(v)

We compute the differential ${\textstyle DJ(\mathbf {v} )}$ of the cost function ${\textstyle J}$ at ${\textstyle \mathbf {v} }$ assumming that ${\textstyle {\mathcal {U}}=L^{2}(0,T;3)}$. To achieve that goal we will use a perturbation analysis. Let us consider thus a perturbation ${\textstyle \delta \mathbf {v} }$ of the control variable ${\textstyle \mathbf {v} }$. We have then,

where in (13):

(i) We denote by a${\textstyle \mathbf {\cdot } }$b the dot product of two vectors a and b.

(ii) The function ${\textstyle \delta \mathbf {y} }$ ${\textstyle =\{\delta y_{1},\delta y_{2},\delta y_{3}\}}$ is the solution of the following initial value problem, obtained by perturbation of (9):

where ${\textstyle C_{1}=\cos \theta _{1}}$, ${\textstyle C_{2}=\cos \theta _{2},}$ ${\textstyle C_{3}=\cos \theta _{3}.}$ In matrix-vector form (14) reads as follows:

where ${\textstyle C}$ and ${\textstyle \Gamma }$ are diagonal matrices with coefficients ${\textstyle C_{i}}$ and ${\textstyle \gamma _{i}}$, ${\textstyle i=1,2,3}$, respectively, and

We introduce now a vector-valued function ${\textstyle \mathbf {p} =\{p_{1},p_{2},p_{3}\}}$ that we assume differentiable over ${\textstyle (0,T)}$; multiplying by ${\textstyle \mathbf {p} }$ both sides of the differential equation in (15), and integrating over ${\textstyle (0,T)}$ we obtain, after integrating by parts, and taking into account the symmetry of the matrices ${\textstyle \Gamma }$ and ${\textstyle K}$, the following equation

Now suppose that the vector-valued function ${\textstyle \mathbf {p} }$ is solution of the following (adjoint system):

It follows from (16), (17) that

Combining (13) and (18) we obtain that

which implies in turn that

Remark 1. So far, we have been assuming that a control ${\textstyle v_{i}}$ is acting on the ${\textstyle i}$-th junction. The method we used to compute ${\textstyle DJ}$ would still apply if one considers controlling the transition from ${\textstyle \mathbf {y} (0)=\delta {\boldsymbol {\phi }}(0)=\delta {\boldsymbol {\theta }}}$ to ${\textstyle \mathbf {y} (T)=\delta {\boldsymbol {\phi }}(T)=\mathbf {0} }$ using only one or only two controls. The same apply for the results in the next sections

2.3.3 Optimality conditions for problem (10)

Let ${\textstyle \mathbf {u} }$ be the solution of problem (10) and let us denote by ${\textstyle \mathbf {y} }$ (respectively, ${\textstyle \mathbf {p} }$) the corresponding solution of the state system (9) (respectively, of the adjoint system (17)). It follows from Subsubsection 3.3.2 that ${\textstyle DJ(\mathbf {u} )=\mathbf {0} }$ is equivalent to the following (optimality) system:

Conversely, it can be shown (see, e.g., [11]) that the system (21)–(25) characterizes ${\textstyle \mathbf {u} }$ as the solution (necessarily unique here) of the control problem (10). The optimality conditions (21)–(25) will play a crucial role concerning the iterative solution of the control problem (10).

2.3.4 Functional equations satisfied by the optimal control

Our goal here is to show that the optimality condition ${\textstyle DJ(\mathbf {u} )=\mathbf {0} }$ can also be written as

where the linear operator ${\textstyle \mathbf {A} }$ is a strongly elliptic and symmetric isomorphism from ${\textstyle {\mathcal {U}}}$ into itself (an automorphism of ${\textstyle {\mathcal {U)}}}$ and where ${\textstyle {\boldsymbol {\beta }}\in {\mathcal {U}}}$. A candidate for ${\textstyle \mathbf {A} }$ is the linear operator from ${\textstyle {\mathcal {U}}}$ into itself defined by

where ${\textstyle \mathbf {p(v)} }$ is obtained from ${\textstyle \mathbf {v} }$ via the successive solution of the following two problems:

which is forward in time, and the adjoint system (17) in ${\textstyle (0,T)}$ which is backward in time. To see that operator ${\textstyle \mathbf {A} }$ is a symmetric and strongly elliptic isomorphism from ${\textstyle {\mathcal {U}}}$ into itself consider ${\textstyle \mathbf {v} ^{1}}$, ${\textstyle \mathbf {v} ^{2}}$ belonging to ${\textstyle {\mathcal {U}}}$ and define ${\textstyle \mathbf {y} ^{i}}$, ${\textstyle \mathbf {p} ^{i}}$ by