Stabilizing a Josephson Junction Array Memory around an unstable equilibrium: A control approach using a first order state model

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:

$\left\{{\begin{array}{c}V={\frac {\displaystyle \hbar }{\displaystyle 2e}}{\frac {\displaystyle d\phi }{\displaystyle dt}},\\I=I_{c}\sin \phi ,\end{array}}\right.$

(1)

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}$)

$\left\{{\begin{array}{c}{\frac {\displaystyle d^{2}\phi _{1}}{\displaystyle dt^{2}}}+\gamma _{1}{\frac {\displaystyle d\phi _{1}}{\displaystyle dt}}+\kappa _{1}(\phi _{1}-\phi _{2})+\sin \phi _{1}=I_{1},\\{\frac {\displaystyle d^{2}\phi _{2}}{\displaystyle dt^{2}}}+\gamma _{2}{\frac {\displaystyle d\phi _{2}}{\displaystyle dt}}+\kappa _{1}(\phi _{2}-\phi _{1})+\kappa _{2}(\phi _{2}-\phi _{3})+\sin \phi _{2}=I_{2},\\{\frac {\displaystyle d^{2}\phi _{3}}{\displaystyle dt^{2}}}+\gamma _{3}{\frac {\displaystyle d\phi _{3}}{\displaystyle dt}}+\kappa _{2}(\phi _{3}-\phi _{2})+\sin \phi _{3}=I_{3},\end{array}}\right.$

(2)

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:

$\gamma _{1}=0.7,\gamma _{2}=1.1,\gamma _{3}=0.7,i_{1}=1,i_{2}=0.8,i_{3}=-1,$ (3) $\kappa _{1}=\kappa _{2}=0.1.$ (4)

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

$\left\{{\begin{array}{c}\gamma _{1}{\frac {\displaystyle d\phi _{1}}{\displaystyle dt}}+\kappa _{1}(\phi _{1}-\phi _{2})+\sin \phi _{1}=I_{1},\\\gamma _{2}{\frac {\displaystyle d\phi _{2}}{\displaystyle dt}}+\kappa _{1}(\phi _{2}-\phi _{1})+\kappa _{2}(\phi _{2}-\phi _{3})+\sin \phi _{2}=I_{2},\\\gamma _{3}{\frac {\displaystyle d\phi _{3}}{\displaystyle dt}}+\kappa _{2}(\phi _{3}-\phi _{2})+\sin \phi _{3}=I_{3}.\end{array}}\right.$

(5)

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

$\theta _{k}=2\pi n_{k}+\sigma _{k}.$

(6)

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).

Table. 1 Stable and unstable equilibrium configurations for the JJAM system defined by (2) and (3). The symbols $(u)$ and $(s)$ qualify the unstable and stable equilibria, respectively.

$n_{1}$	$n_{2}$	$n_{3}$	$\sigma _{1}/2\pi$	$\sigma _{2}/2\pi$	$\sigma _{3}/2\pi$	$\theta _{1}$	$\theta _{2}$	$\theta _{3}$
$0$	$0$	$0(s)$	$0.1992$	$0.1187$	$-0.1552$	$1.2517$	$0.7458$	$-0.9752$
$1$	$0$	$0(s)$			no	stable	equil
$1$	$1$	$0(s)$	$0.1810$	$0.0338$	$-0.0515$	$7.4207$	$6.4958$	$-0.3236$
$2$	$0$	$0(s)$			no	stable	equil
$2$	$1$	$0(s)$	$0.0661$	$0.1164$	$-0.0437$	$12.9821$	$7.0148$	$-0.2746$
$2$	$2$	$0(s)$	$0.1671$	$-0.0443$	$0.0333$	$13.6161$	$12.2884$	$0.2092$
$1$	$0$	$0(u)$	$0.0908$	$0.3591$	$-0.4140$	$6.8539$	$2.2566$	$-2.6015$
$2$	$1$	$0(u)$	$0.1011$	$0.4539$	$-0.0125$	$13.2016$	$9.1355$	$-0.0786$
$2$	$2$	$0(u)$	$0.1576$	$-0.1029$	$-0.5938$	$13.5568$	$11.9196$	$-3.7436$

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 $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 $\{1,1,0(s)\}$(left), $\{2,1,0(s)\}$(right) from an approximation to the unstable equilibrium $\{1,0,0(u)\}$(left), $\{2,1,0(u)\}$(right).

Figure 3: Evolution to $\{2,2,0(s)\}$ from an approximation to the unstable equilibrium $\{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

$\left\{{\begin{array}{c}\gamma _{1}{\frac {\displaystyle d\delta \phi _{1}}{\displaystyle dt}}+\kappa _{1}(\delta \phi _{1}-\delta \phi _{2})+\delta \phi _{1}\cos \theta _{1}=0,\\\gamma _{2}{\frac {\displaystyle d\delta \phi _{2}}{\displaystyle dt}}+\kappa _{1}(\delta \phi _{2}-\delta \phi _{1})+\kappa _{2}(\delta \phi _{2}-\delta \phi _{3})+\delta \phi _{2}\cos \theta _{2}=0,\\\gamma _{3}{\frac {\displaystyle d\delta \phi _{3}}{\displaystyle dt}}+\kappa _{2}(\delta \phi _{3}-\delta \phi _{2})+\delta \phi _{3}\cos \theta _{3}=0,\\\delta {\boldsymbol {\phi }}(0)=\delta {\boldsymbol {\theta }}.\end{array}}\right.$

(7)

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 $\mathbf {y} =\delta {\boldsymbol {\phi }}$ of (7) for the unstable equilibrium given by $\{n_{1},n_{2},n_{3}\}=\{1,0,0(u)\}$ and $\delta {\boldsymbol {\theta }}=[1e-2,1e-2,1e-2]$.

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

Figure 6: Solution $\mathbf {y} =\delta {\boldsymbol {\phi }}$ of (7) for the unstable equilibrium given by $\{n_{1},n_{2},n_{3}\}=\{2,2,0(u)\}$ and $\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

${\boldsymbol {\phi }}(0)={\boldsymbol {\theta }}+\delta {\boldsymbol {\theta }}.$

(8)

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):

$\left\{{\begin{array}{c}\gamma _{1}{\frac {\displaystyle dy_{1}}{\displaystyle dt}}+\kappa _{1}(y_{1}-y_{2})+y_{1}\cos \theta _{1}=v_{1},\\\gamma _{2}{\frac {\displaystyle dy_{2}}{\displaystyle dt}}+\kappa _{1}(y_{2}-y_{1})+\kappa _{2}(y_{2}-y_{3})+y_{2}\cos \theta _{2}=v_{2},\\\gamma _{3}{\frac {\displaystyle dy_{3}}{\displaystyle dt}}+\kappa _{2}(y_{3}-y_{2})+y_{3}\cos \theta _{3}=v_{3},\\\mathbf {y} (0)=\delta {\boldsymbol {\theta }}.\end{array}}\right.$

(9)

via the following control formulation:

$\left\{{\begin{array}{c}\mathbf {u} \in {\mathcal {U}},\\J(\mathbf {u} )\leq J(\mathbf {v} ),\forall \mathbf {v} \in {\mathcal {U}},\end{array}}\right.$

(10)

where

$J(\mathbf {v} )={\frac {1}{2}}\int \limits _{0}^{T}||\mathbf {v} ||^{2}dt+{\frac {k_{1}}{2}}\int \limits _{0}^{T}||\mathbf {y} ||^{2}dt+{\frac {k_{2}}{2}}||\mathbf {y} (T)||^{2},$

(11)

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

$(\mathbf {v,w} )_{\mathcal {U}}=\int \limits _{0}^{T}\mathbf {v(} t\mathbf {)\cdot w(} t\mathbf {)} {\hbox{ }}dt,$

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

$\left\langle DJ(\mathbf {v} ),\mathbf {w} \right\rangle =\int \limits _{0}^{T}DJ(\mathbf {v} )\mathbf {(} t\mathbf {)\cdot w(} t\mathbf {)} {\hbox{ }}dt,{\hbox{ }}\forall \mathbf {v,w} \in {\mathcal {U}}{\hbox{,}}$

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

$DJ(\mathbf {u} )=\mathbf {0} .$

(12)

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

$DJ(\mathbf {u} )-DJ(\mathbf {0} )=-DJ(\mathbf {0} ),$

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,

${\begin{array}{c}\delta J(\mathbf {v} )=\int \limits _{0}^{T}DJ(\mathbf {v} )\cdot \delta \mathbf {v} dt=\int \limits _{0}^{T}\mathbf {v\cdot } \delta \mathbf {v} dt+k_{1}\int \limits _{0}^{T}\mathbf {y} \cdot \delta \mathbf {y} dt+k_{2}\mathbf {y} (T)\cdot \delta \mathbf {y} (T),\end{array}}$

(13)

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):

$\left\{{\begin{array}{c}\gamma _{1}{\frac {\displaystyle d\delta y_{1}}{\displaystyle dt}}+\kappa _{1}(\delta y_{1}-\delta y_{2})+C_{1}\delta y_{1}=\delta v_{1}{\hbox{, in }}(0,T),\\\gamma _{2}{\frac {\displaystyle d\delta y_{2}}{\displaystyle dt}}+\kappa _{1}(\delta y_{2}-\delta y_{1})+\kappa _{2}(\delta y_{2}-\delta y_{3})+C_{2}\delta y_{2}=\delta v_{2}{\hbox{, in }}(0,T),\\\gamma _{3}{\frac {\displaystyle d\delta y_{3}}{\displaystyle dt}}+\kappa _{\hbox{2}}(\delta y_{3}-\delta y_{2})+C_{3}\delta y_{3}=\delta v_{3}{\hbox{, in }}(0,T),\\\mathbf {\delta y} (0)=\mathbf {0} ,\end{array}}\right.$

(14)

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:

$\left\{{\begin{array}{c}\Gamma {\frac {\displaystyle d}{\displaystyle dt}}\delta \mathbf {y} +(K+C)\delta \mathbf {y} =\delta \mathbf {v} {\hbox{, in }}(0,T),\\\mathbf {\delta y} (0)=\mathbf {0} .\end{array}}\right.$

(15)

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

$K=\left({\begin{array}{ccc}\kappa _{1}&-\kappa _{1}&0\\-\kappa _{1}&\kappa _{1}+\kappa _{2}&-\kappa _{2}\\0&-\kappa _{2}&\kappa _{2}\end{array}}\right).$

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

$\Gamma \mathbf {p} (T)\cdot \delta \mathbf {y} (T)+\int \limits _{0}^{T}\left\{-\Gamma {\frac {d}{dt}}\mathbf {p} {\hbox{ }}+(K+C)\mathbf {p} {\hbox{ }}\right\}\mathbf {\cdot } \delta \mathbf {y} {\hbox{ }}dt=\int \limits _{0}^{T}\mathbf {p\cdot } \delta \mathbf {v} {\hbox{ }}dt,{\hbox{ in }}(0,T).$

(16)

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

$\left\{{\begin{array}{c}-\mathbf {\Gamma } {\frac {\displaystyle d}{\displaystyle dt}}\mathbf {p} {\hbox{ }}+(K+C)\mathbf {p} {\hbox{ }}=k_{1}\mathbf {y,} \\\mathbf {\Gamma p(} T\mathbf {)=} k_{2}\mathbf {y} (T).\end{array}}\right.$

(17)

It follows from (16), (17) that

$k_{2}\mathbf {y} (T)\cdot \delta \mathbf {y} (T)+k_{1}\int \limits _{0}^{T}\mathbf {y\cdot } \delta \mathbf {y} {\hbox{ }}dt\mathbf {=} \int \limits _{0}^{T}\mathbf {p\cdot } \delta \mathbf {v} {\hbox{ }}dt.$

(18)

Combining (13) and (18) we obtain that

$\int \limits _{0}^{T}DJ(\mathbf {v} )\cdot \delta \mathbf {v} dt=\int \limits _{0}^{T}\mathbf {v\cdot } \delta \mathbf {v} dt+k_{1}\int \limits _{0}^{T}\mathbf {y} \cdot \delta \mathbf {y} dt+k_{2}\mathbf {y} (T)\cdot \delta \mathbf {y} (T)=\int \limits _{0}^{T}(\mathbf {v} +\mathbf {p} )\cdot \delta \mathbf {v} dt,$

(19)

which implies in turn that

$DJ(\mathbf {v} )=\mathbf {v} +\mathbf {p} .$

(20)

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:

$\mathbf {u} =-\mathbf {p} ,$ (21) $\Gamma {\frac {d\mathbf {y} }{dt}}+(K+C)\mathbf {y} =\mathbf {u} ,{\hbox{ in }}(0,T),$ (22) $\mathbf {y} (0)=\delta {\boldsymbol {\theta }},$ (23) $-\mathbf {\Gamma } {\frac {d\mathbf {p} }{dt}}{\hbox{ }}+(K+C)\mathbf {p} {\hbox{ }}\mathbf {=} k_{1}\mathbf {y,} {\hbox{ in }}(0,T),$ (24) $\mathbf {\Gamma p(} T\mathbf {)} \mathbf {=} k_{2}\mathbf {y} (T).$ (25)

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

$\mathbf {Au} ={\boldsymbol {\beta }},$

(26)

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

$\mathbf {Av} =\mathbf {v+p(v),}$

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

$\left\{{\begin{array}{c}\Gamma {\frac {\displaystyle d}{\displaystyle dt}}\mathbf {y} +(K+C)\mathbf {y} =\mathbf {v} ,{\hbox{ in }}(0,T),\\\mathbf {y} (0)=\mathbf {0} ,\end{array}}\right.$

(27)

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

$\mathbf {y} ^{i}=\mathbf {y} (\mathbf {v} ^{i}),{\hbox{ }}\mathbf {p} ^{i}=\mathbf {p} (\mathbf {v} ^{i}),{\hbox{ }}i=1,2{\hbox{;}}$

we have then

$\int \limits _{0}^{T}(A\mathbf {v} ^{1})\cdot \mathbf {v} ^{2}{\hbox{ }}dt=\int \limits _{0}^{T}(\mathbf {v} ^{1}+\mathbf {p} ^{1})\cdot \mathbf {v} ^{2}{\hbox{ }}dt=\int \limits _{0}^{T}\mathbf {v} ^{1}\cdot \mathbf {v} ^{2}{\hbox{ }}dt+\int \limits _{0}^{T}\mathbf {p} ^{1}\cdot \mathbf {v} ^{2}{\hbox{ }}dt.$

(28)

On the other hand, we have, starting with the differential equation in (27) and using integration by parts, that

$0=\int \limits _{0}^{T}\left(\Gamma {\frac {d}{dt}}\mathbf {y} ^{2}+(K+C)\mathbf {y} ^{2}-\mathbf {v} ^{2}\right)\cdot \mathbf {p} ^{1}{\hbox{ }}dt$ $=\mathbf {\Gamma p} ^{1}\mathbf {(} T\mathbf {)\cdot y} ^{2}(T)+\int \limits _{0}^{T}\left(-\Gamma {\frac {d}{dt}}\mathbf {p} ^{1}+(K+C)\mathbf {p} ^{1}\right)\cdot \mathbf {y} ^{2}dt-\int \limits _{0}^{T}\mathbf {v} ^{2}\cdot \mathbf {p} ^{1}dt$ $=k_{2}\mathbf {y} ^{1}(T)\mathbf {\cdot y} ^{2}(T)+\int \limits _{0}^{T}k_{1}\mathbf {y} ^{1}\cdot \mathbf {y} ^{2}dt-\int \limits _{0}^{T}\mathbf {v} ^{2}\cdot \mathbf {p} ^{1}dt.$

This implies that

$\int \limits _{0}^{T}\mathbf {v} ^{2}\cdot \mathbf {p} ^{1}dt=k_{2}\mathbf {y} ^{1}(T)\mathbf {\cdot y} ^{2}(T)+\int \limits _{0}^{T}k_{1}\mathbf {y} ^{1}\cdot \mathbf {y} ^{2}dt.$

(29)

Combining (28) with (29) we obtain

$\int \limits _{0}^{T}(\mathbf {Av} ^{1})\cdot \mathbf {v} ^{2}{\hbox{ }}dt=\int \limits _{0}^{T}\mathbf {v} ^{1}\cdot \mathbf {v} ^{2}{\hbox{ }}dt+k_{2}\mathbf {y} ^{1}(T)\mathbf {\cdot y} ^{2}(T)+\int \limits _{0}^{T}k_{1}\mathbf {y} ^{1}\cdot \mathbf {y} ^{2}dt.$

(30)

Relation (30) imply the symmetry of ${\textstyle \mathbf {A} }$; we also have

$\int \limits _{0}^{T}(\mathbf {Av} )\cdot \mathbf {v} {\hbox{ }}dt\geq \int \limits _{0}^{T}\left\vert \mathbf {v} \right\vert ^{2}{\hbox{ }}dt,{\hbox{ }}\forall \mathbf {v} \in {\mathcal {U}}{\hbox{,}}$

which implies the strong ellipticity of ${\textstyle \mathbf {A} }$ over ${\textstyle {\mathcal {U}}}$. The linear operator ${\textstyle \mathbf {A} }$, being continuous and strongly elliptic over ${\textstyle {\mathcal {U}}}$, is an automorphism of ${\textstyle {\mathcal {U}}}$. To identify the right hand side ${\textstyle {\boldsymbol {\beta }}}$ of equation (26), we introduce ${\textstyle \mathbf {Y} _{0}}$ and ${\textstyle \mathbf {P} _{0}}$ defined as the solutions of

$\left\{{\begin{array}{c}\Gamma {\frac {\displaystyle d}{\displaystyle dt}}\mathbf {Y} _{0}+(K+C)\mathbf {Y} _{0}=\mathbf {0} {\hbox{, in }}(0,T),\\\mathbf {Y} _{0}(0)=\delta {\boldsymbol {\theta }},\end{array}}\right.$

(31)

and

$\left\{{\begin{array}{c}-\mathbf {\Gamma } {\frac {\displaystyle d}{\displaystyle dt}}\mathbf {P} _{0}{\hbox{ }}+(K+C)\mathbf {P} _{0}\mathbf {=} k_{1}\mathbf {Y} _{0}\mathbf {,} \\\mathbf {\Gamma \mathbf {P} _{0}(} T\mathbf {)=} k_{2}\mathbf {Y} _{0}(T).\end{array}}\right.$

(32)

Suppose now that ${\textstyle \mathbf {y} }$ and ${\textstyle \mathbf {p} }$ satisfies the optimality conditions. Define

${\overline {\mathbf {y} }}=\mathbf {y-Y} _{0},$ ${\overline {\mathbf {p} }}=\mathbf {p-P} _{0}.$

Subtracting (31) and (32) from (27) and (17) we obtain

$\left\{{\begin{array}{c}\Gamma {\frac {\displaystyle d}{\displaystyle dt}}{\overline {\mathbf {y} }}+(K+C){\overline {\mathbf {y} }}=\mathbf {u} {\hbox{, in }}(0,T),\\{\overline {\mathbf {y} }}(0)=\mathbf {0} ,\end{array}}\right.$

(33)

and

$\left\{{\begin{array}{c}-\mathbf {\Gamma } {\frac {\displaystyle d}{\displaystyle dt}}{\overline {\mathbf {p} }}{\hbox{ }}+(K+C){\overline {\mathbf {p} }}\mathbf {=} k_{1}{\overline {\mathbf {y} }}\mathbf {,} \\\mathbf {\Gamma {\overline {\mathbf {p} }}(} T\mathbf {)=} k_{2}{\overline {\mathbf {y} }}(T).\end{array}}\right.$

(34)

From the definition of operator ${\textstyle \mathbf {A} }$ it follows that

$\mathbf {Au=u+} {\overline {\mathbf {p} }}=\mathbf {-p+} {\overline {\mathbf {p} }}=-\mathbf {P} _{0};$

(35)

the right hand side of (35) is the vector ${\textstyle {\boldsymbol {\beta }}}$ that we were looking for.

From the properties of ${\textstyle \mathbf {A} }$, problem (35) can be solved by a conjugate gradient algorithm operating in the Hilbert space ${\textstyle {\mathcal {U}}}$. This algorithm will be described in the following subsubsection.

2.3.5 Conjugate gradient solution of the control problem (10)

Problem (35) can be written in variational form as

$\mathbf {u} \in {\mathcal {U}}{\hbox{ }}(=L^{2}(0,T;3)),$ $\int \limits _{0}^{T}(\mathbf {Au} )\cdot \mathbf {v} {\hbox{ }}dt=-\int \limits _{0}^{T}\mathbf {P} _{0}\cdot \mathbf {v} {\hbox{ }}dt,\forall \mathbf {v} \in {\mathcal {U}}{\hbox{ .}}$

(36)

From the symmetry and ${\textstyle {\mathcal {U}}}$-ellipticity of the bilinear form

$\left\{\mathbf {v} ,\mathbf {w} \right\}\rightarrow \int \limits _{0}^{T}(\mathbf {Av} )\cdot \mathbf {w} {\hbox{ }}dt,$

the variational problem (36) is a particular case of the following class of linear variational problems:

$\left\{{\begin{array}{c}u\in {\mathcal {V}}{\hbox{,}}\\a(u,v)=L(v),\forall v\in {\mathcal {V}}{\hbox{,}}\end{array}}\right.$

(37)

where

(i) ${\textstyle {\mathcal {V}}}$ is a real Hilbert space for the scalar product ${\textstyle (\cdot ,\cdot )}$ and the correspondent norm ${\textstyle \left\Vert \cdot \right\Vert ;}$

(ii) ${\textstyle a}$ : ${\textstyle {\mathcal {V\times V\rightarrow \mathbb {R} }}}$ is bilinear, continuous, symmetric and ${\textstyle {\mathcal {V}}}$-elliptic, i.e., ${\textstyle \exists \alpha >0}$ such that ${\textstyle a(v,v)\geq \alpha \left\Vert v\right\Vert ^{2},\forall v\in {\mathcal {V}}}$;

(iii) ${\textstyle L}$ : ${\textstyle {\mathcal {V\rightarrow \mathbb {R} }}}$ is linear and continuous.

If the above properties hold, then problem (37) has a unique solution (see [10]); in fact, the symmetry of ${\textstyle a(\cdot ,\cdot )}$ is not necessary in order to have a unique solution to problem (37); however, the symmetry of ${\textstyle a(\cdot ,\cdot )}$, combined with the other properties, allows problem (37) to be solved by the following conjugate gradient algorithm:

• Step 1. ${\textstyle u^{0}\in {\mathcal {V}}}$ is given.

• Step 2. Solve

$\left\{{\begin{array}{c}g^{0}\in {\mathcal {V}}{\hbox{,}}\\(g^{0},v)=a(u^{0},v)-L(v),\forall v\in {\mathcal {V}}{\hbox{,}}\end{array}}\right.$

(38)

and set

$w^{0}=g^{0}.$

(39)

• Step 3. For ${\textstyle q\geq 0}$, ${\textstyle u^{q}}$, ${\textstyle g^{q}}$ and ${\textstyle w^{q}}$ being known, compute ${\textstyle u^{q+1}}$ , ${\textstyle g^{q+1}}$ and ${\textstyle w^{q+1}}$ as follows:

$\rho ^{q}=\left\Vert g^{q}\right\Vert ^{2}/a(w^{q},w^{q}),$

(40)

and take

$u^{q+1}=u^{q}-\rho ^{q}w^{q}.$

(41)

Solve

$\left\{{\begin{array}{c}g^{q+1}\in {\mathcal {V}}{\hbox{,}}\\(g^{q+1},v)=(g^{q},v)-\rho ^{q}a(w^{q},v),\forall v\in {\mathcal {V}}{\hbox{,}}\end{array}}\right.$

(42)

and compute

$\gamma ^{q}=\left\Vert g^{q+1}\right\Vert ^{2}/\left\Vert g^{q}\right\Vert ^{2},$

(43)

$w^{q+1}=g^{q+1}+\gamma ^{q}w^{q}.$

(44)

• Step 4. Do n=n+1 and go to (40).

Let us recall that ${\textstyle {\mathcal {U}}=L^{2}(0,T;3)}$ is a Hilbert space for the inner-product ${\textstyle \{\mathbf {v} ,\mathbf {w} \}}$ ${\textstyle \longrightarrow }$ ${\textstyle \int \limits _{0}^{T}\mathbf {v\cdot w} dt}$ and the associated norm ${\textstyle \mathbf {v} }$ ${\textstyle \longrightarrow }$ ${\textstyle {\sqrt {\int \limits _{0}^{T}\left\vert \mathbf {v} \right\vert ^{2}dt}}}$, implying that problem (10)-(36) can be solved applying the conjugate gradient algorithm (38)-(44). The above algorithm takes the following form:

• Suppose

$\mathbf {u} ^{0}{\hbox{ is given in }}L^{2}(0,T;3){\hbox{ }}(\mathbf {u} ^{0}=\mathbf {0} {\hbox{ for example}}).$

(45)

• Solve

$\left\{{\begin{array}{c}\Gamma {\frac {d}{dt}}\mathbf {y} ^{0}+(K+C)\mathbf {y} ^{0}=\mathbf {u} ^{0}{\hbox{, in }}(0,T),\\\mathbf {y} ^{0}(0)=\delta {\boldsymbol {\theta }}.\end{array}}\right.,$

(46)

and then

$\left\{{\begin{array}{c}-\mathbf {\Gamma } {\frac {d}{dt}}\mathbf {p} ^{0}{\hbox{ }}+(K+C)\mathbf {p} ^{0}{\hbox{ }}\mathbf {=} k_{1}\mathbf {\mathbf {y} } ^{0}\mathbf {,} {\hbox{ in }}(0,T),\\\mathbf {\Gamma p} ^{0}\mathbf {(} T\mathbf {)=} k_{2}\mathbf {y} ^{0}(T).\end{array}}\right.$

(47)

• Set

$\mathbf {g} ^{0}=\mathbf {u} ^{0}+\mathbf {p} ^{0}.$

(48)

• If ${\textstyle \int \limits _{0}^{T}\left\vert \mathbf {g} ^{0}\right\vert ^{2}dt\leq tol^{2}\max[1,\int \limits _{0}^{T}\left\vert \mathbf {u} ^{0}\right\vert ^{2}dt],}$ take ${\textstyle \mathbf {u} =\mathbf {u} ^{0}}$; otherwise, set

$\mathbf {w} ^{0}=\mathbf {g} ^{0}.$

(49)

For ${\textstyle q\geq 0}$, ${\textstyle \mathbf {u} ^{q}}$, ${\textstyle \mathbf {g} ^{q}}$ and ${\textstyle \mathbf {w} ^{q}}$ being known, we compute ${\textstyle \mathbf {u} ^{q+1},}$ ${\textstyle \mathbf {g} ^{q+1}}$ and if necessary, ${\textstyle \mathbf {w} ^{q+1}}$ as follows:

• Solve

$\left\{{\begin{array}{c}\Gamma {\frac {d}{dt}}{\overline {\mathbf {y} }}^{q}+(K+C){\overline {\mathbf {y} }}^{q}=\mathbf {w} ^{q}{\hbox{ in }}(0,T),\\\mathbf {y} ^{0}(0)=\mathbf {0} .\end{array}}\right.$

(50)

and then

$\left\{{\begin{array}{c}-\mathbf {\Gamma } {\frac {d}{dt}}{\overline {\mathbf {p} }}^{q}{\hbox{ }}+(K+C){\overline {\mathbf {p} }}^{q}{\hbox{ }}\mathbf {=} k_{1}{\overline {\mathbf {y} }}^{q}\mathbf {,} {\hbox{ in }}(0,T),\\\mathbf {\Gamma } {\overline {\mathbf {p} }}^{q}\mathbf {(} T\mathbf {)=} k_{2}{\overline {\mathbf {y} }}^{q}(T).\end{array}}\right.$

(51)

Set

${\overline {\mathbf {g} }}^{q}=\mathbf {w} ^{q}+{\overline {\mathbf {p} }}^{q},$

(52)

and

$\rho _{q}=\left\Vert \mathbf {g} ^{q}\right\Vert ^{2}/\int \limits _{0}^{T}{\overline {\mathbf {g} }}^{q}\mathbf {\cdot w} ^{q}dt.$

(53)

• Set

$\mathbf {u} ^{q+1}=\mathbf {u} ^{q}-\rho _{q}\mathbf {w} ^{q}.$

(54)

• Set

$\mathbf {g} ^{q+1}=\mathbf {g} ^{q}-\rho _{q}{\overline {\mathbf {g} }}^{q}.$

(55)

• If ${\textstyle {\frac {\int \limits _{0}^{T}\left\vert \mathbf {g} ^{q+1}\right\vert ^{2}dt}{\max[\int \limits _{0}^{T}\left\vert \mathbf {g} ^{0}\right\vert ^{2}dt,\int \limits _{0}^{T}\left\vert \mathbf {u} ^{q+1}\right\vert ^{2}dt]}}\leq tol^{2},}$ take ${\textstyle \mathbf {u} =\mathbf {u} ^{q+1}}$; otherwise, compute

$\gamma _{q}={\frac {\int \limits _{0}^{T}\left\vert \mathbf {g} ^{q+1}\right\vert ^{2}dt}{\int \limits _{0}^{T}\left\vert \mathbf {g} ^{q}\right\vert ^{2}dt}},$

(56)

and set

$\mathbf {w} ^{q+1}=\mathbf {g} ^{q+1}+\gamma _{q}\mathbf {w} ^{q}.$

(57)

• Do ${\textstyle q+1\longrightarrow q}$ and return to (50) .

• End of algorithm.

The state variable ${\textstyle \mathbf {y} ^{q}}$ can be actualized simultaneously with the control ${\textstyle \mathbf {u} ^{q}}$ since for all ${\textstyle q}$ we have

$\left\{{\begin{array}{c}\Gamma {\frac {\displaystyle d}{\displaystyle dt}}\mathbf {y} ^{q}+(K+C)\mathbf {y} ^{q}=\mathbf {u} ^{q}{\hbox{, in }}(0,T),\\\mathbf {y} ^{q}(0)=\delta {\boldsymbol {\theta }}.\end{array}}\right.$

(58)

so,

$\left\{{\begin{array}{c}\Gamma {\frac {\displaystyle d}{\displaystyle dt}}(\mathbf {y} ^{q+1}-\mathbf {y} ^{q})+(K+C)(\mathbf {y} ^{q+1}-\mathbf {y} ^{q})=\mathbf {u} ^{q+1}-\mathbf {u} ^{q}=-\rho _{q}\mathbf {w} ^{q}{\hbox{, in }}(0,T),\\(\mathbf {y} ^{q+1}-\mathbf {y} ^{q})(0)=\mathbf {0} ,\end{array}}\right.$

(59)

which, by definition of ${\textstyle {\overline {\mathbf {y} }}^{q},}$ implies that

$\mathbf {y} ^{q+1}-\mathbf {y} ^{q}=-\rho _{q}{\overline {\mathbf {y} }}^{q}\Longrightarrow \mathbf {y} ^{q+1}=\mathbf {y} ^{q}-\rho _{q}{\overline {\mathbf {y} }}^{q}.$

The practical implementation of algorithm (45)-(57), via a finite difference discretization of problem (10), will be discussed in the following section.

3 Discrete formulation of the optimal control problem

3.1 Finite difference approximation of problem (10)

We approximate (10) when ${\textstyle {\mathcal {U}}=(L^{2}(0,T))^{3}}$ by

$\left\{{\begin{array}{c}\mathbf {u} ^{\Delta t}\in {\mathcal {U}}^{\Delta t},\\J^{\Delta t}(\mathbf {u} ^{\Delta t})\leq J^{\Delta t}(\mathbf {v} ),\forall \mathbf {v} \in {\mathcal {U}}^{\Delta t},\end{array}}\right.$

(60)

where:

• ${\textstyle \Delta t=T/N}$ with ${\textstyle N}$ a "large" positive integer.

• ${\textstyle {\mathcal {U}}^{\Delta t}=(\mathbb {R} ^{3})^{N}.}$

The cost functional ${\textstyle J^{\Delta t}}$ is defined by

$J^{\Delta t}(\mathbf {v} )={\frac {\Delta t}{2}}\sum _{n=1}^{N}||\mathbf {v} ^{n}||^{2}+{\frac {k_{1}\Delta t}{2}}\sum _{n=1}^{N}||\mathbf {y} ^{n}||^{2}+{\frac {k_{2}}{2}}||\mathbf {y} ^{N}||^{2},$

with ${\textstyle \mathbf {v} =\{\mathbf {v} ^{n}\}_{n=1}^{N}}$ and ${\textstyle \{\mathbf {y} ^{n}\}_{n=1}^{N}}$ obtained from ${\textstyle \mathbf {v} }$ and ${\textstyle \delta {\boldsymbol {\theta }}}$ via the following discrete variant of (9):

$\mathbf {y} ^{0}=\delta {\boldsymbol {\theta }},$

(61)

and for ${\textstyle n=1,....,N}$

$\Gamma {\frac {\mathbf {y} ^{n}-\mathbf {y} ^{n-1}}{\Delta t}}+(K+C)\mathbf {y} ^{n}=\mathbf {v} ^{n}.$

(62)

To compute ${\textstyle \mathbf {y} ^{n}}$ we have thus to solve a linear system of the following type:

$(\Gamma +\Delta t(K+C))\mathbf {y} ^{n}=RHS^{n}.$

(63)

The matrix ${\textstyle \Gamma +\Delta t(K+C)}$ being a ${\textstyle 3\times 3}$ matrix symmetric and positive definite, solving (63) is easy.

3.2 Optimality Conditions and conjugate gradient solution of (60)

3.2.1 Computing DJ^∆t(v)

Assuming that one wants to use the conjugate gradient algorithm (38)-(44) to solve the discrete problem (60), we compute first ${\textstyle DJ^{\Delta t}(\mathbf {v} )}$. On ${\textstyle {\mathcal {U}}^{\Delta t}=(\mathbb {R} ^{3})^{N}}$ we will use the following inner-product

$(\mathbf {v} ,\mathbf {w} )_{\Delta t}=\Delta t\sum _{n=1}^{N}\mathbf {v} ^{n}\cdot \mathbf {w} ^{n},{\hbox{ }}\forall \mathbf {v=} \{\mathbf {v} ^{n}\}_{n=1}^{N}{\hbox{, }}\mathbf {w=} \{\mathbf {w} ^{n}\}_{n=1}^{N}\in (\mathbb {R} ^{3})^{N}.$

We have

$\delta J^{\Delta t}(\mathbf {v} )=\Delta t\sum _{n=1}^{N}\mathbf {v} ^{n}\cdot \delta \mathbf {v} ^{n}+k_{1}\Delta t\sum _{n=1}^{N}\mathbf {y} ^{n}\cdot \delta \mathbf {y} ^{n}+k_{2}\mathbf {y} ^{N}\cdot \delta \mathbf {y} ^{N},$

(64)

with ${\textstyle \{\delta \mathbf {y} ^{n}\}_{n=0}^{N}}$ obtained by perturbation of (61), (62), that is

$\delta \mathbf {y} ^{0}=\mathbf {0} ,$

(65)

and for ${\textstyle n=1,....,N}$

$\Gamma {\frac {\delta \mathbf {y} ^{n}-\delta \mathbf {y} ^{n-1}}{\Delta t}}+(K+C)\delta \mathbf {y} ^{n}=\mathbf {\delta v} ^{n}.{\hbox{ }}$

(66)

Let us introduce ${\textstyle \{\mathbf {p} ^{n}\}_{n=1}^{N}\in (\mathbb {R} ^{3})^{N}.}$ Taking the dot product of ${\textstyle \mathbf {p} ^{n}}$ with each side of equation (66), we obtain after summation and multiplication by ${\textstyle \Delta t:}$

$\Delta t\sum _{n=1}^{N}\Gamma {\frac {\delta \mathbf {y} ^{n}-\delta \mathbf {y} ^{n-1}}{\Delta t}}\cdot \mathbf {p} ^{n}+\Delta t\sum _{n=1}^{N}(K+C)\delta \mathbf {y} ^{n}\cdot \mathbf {p} ^{n}=\Delta t\sum _{n=1}^{N}\delta \mathbf {v} ^{n}\cdot \mathbf {p} ^{n}.{\hbox{ }}$

(67)

Applying a discrete integration by parts to relation (67), and considering that ${\textstyle \delta \mathbf {y} ^{0}=\mathbf {0,} }$ we obtain:

$\Gamma \mathbf {p} ^{N+1}\cdot \delta \mathbf {y} ^{N}+\Delta t\sum _{n=1}^{N}\Gamma {\frac {\mathbf {p} ^{n}-\mathbf {p} ^{n+1}}{\Delta t}}\cdot \delta \mathbf {y} ^{n}+\Delta t\sum _{n=1}^{N}(K+C)\mathbf {p} ^{n}\cdot \delta \mathbf {y} ^{n}=\Delta t\sum _{n=1}^{N}\delta \mathbf {v} ^{n}\cdot \mathbf {p} ^{n},{\hbox{ }}$

(68)

or equivalently

$\Gamma \mathbf {p} ^{N+1}\cdot \delta \mathbf {y} ^{N}+\Delta t\sum _{n=1}^{N}\left\{\Gamma {\frac {\mathbf {p} ^{n}-\mathbf {p} ^{n+1}}{\Delta t}}+(K+C)\mathbf {p} ^{n}\right\}\cdot \delta \mathbf {y} ^{n}=\Delta t\sum _{n=1}^{N}\delta \mathbf {v} ^{n}\cdot \mathbf {p} ^{n}.{\hbox{ }}$

(69)

Suppose that ${\textstyle \{\mathbf {p} ^{n}\}_{n=1}^{N+1}}$ verifies the following discrete adjoint system:

$\Gamma \mathbf {p} ^{N+1}=k_{2}\mathbf {y} ^{N},$

(70)

and for ${\textstyle n=N,...,1}$

$\Gamma {\frac {\mathbf {p} ^{n}-\mathbf {p} ^{n+1}}{\Delta t}}+(K+C)\mathbf {p} ^{n}=k_{1}\mathbf {y} ^{n}\mathbf {.}$

(71)

It follows from (64), (69)-(71) that

$\delta J^{\Delta t}(\mathbf {v} )=\Delta t\sum _{n=1}^{N}\mathbf {v} ^{n}\cdot \delta \mathbf {v} ^{n}+k_{1}\Delta t\sum _{n=1}^{N}\mathbf {y} ^{n}\cdot \delta \mathbf {y} ^{n}+k_{2}\mathbf {y} ^{N}\cdot \delta \mathbf {y} ^{N}=\Delta t\sum _{n=1}^{N}(\mathbf {v} ^{n}+\mathbf {p} ^{n})\cdot \delta \mathbf {v} ^{n},$

(72)

that is

$DJ^{\Delta t}(\mathbf {v} )=\{\mathbf {v} ^{n}+\mathbf {p} ^{n}\}_{n=1}^{N}.$

(73)

3.2.2 Optimality conditions for (60)

The optimality conditions for the discrete problem (60) are

$\mathbf {u} ^{n}=-\mathbf {p} ^{n}{\hbox{ , }}n=1,...,N,$ (74) $\mathbf {y} ^{0}=\delta {\boldsymbol {\theta }},$ (75) $\Gamma {\frac {\mathbf {y} ^{n}-\mathbf {y} ^{n-1}}{\Delta t}}+(K+C)\mathbf {y} ^{n}=\mathbf {u} ^{n}{\hbox{ in }}(0,T),{\hbox{ }}n=1,...,N,$ (76) $\Gamma \mathbf {p} ^{N+1}=k_{2}\mathbf {y} ^{N},$ (77) $\Gamma {\frac {\mathbf {p} ^{n}-\mathbf {p} ^{n+1}}{\Delta t}}+(K+C)\mathbf {p} ^{n}=k_{1}\mathbf {y} ^{n}{\hbox{ , }}n=N,...,1.$ (78)

3.2.3 Functional equation for the discrete control solution of (60)

Following the sketch for the continuous case we can show that the discrete version ${\textstyle \mathbf {A} ^{\Delta t}}$ of operator ${\textstyle \mathbf {A} }$ and the discrete version ${\textstyle {\boldsymbol {\beta }}^{\Delta t}}$ of ${\textstyle {\boldsymbol {\beta }}}$ satisfies the equation

$\mathbf {A} ^{\Delta t}\mathbf {u} ^{\Delta t}={\boldsymbol {\beta }}^{\Delta t},$

(79)

where ${\textstyle \mathbf {u} ^{\Delta t}}$ is the discrete control satisfying the optimality condition (74). Operator ${\textstyle \mathbf {A} ^{\Delta t}}$ enjoys the same properties than the continuous version: symmetric, strongly elliptic and continuous, allowing us to use a conjugate gradient like (38)-(44) to solve (79).

3.2.4 Conjugate gradient solution of the discrete control problem (60)

Using ${\textstyle \mathbf {y} _{q}^{n}=\{y_{iq}^{n}\}_{i=1}^{3}}$ to denote the discrete value of the vector-valued function ${\textstyle \mathbf {y} }$ at time ${\textstyle n\Delta t}$ and iteration ${\textstyle q}$; similarly, ${\textstyle \mathbf {u} _{q}^{n}}$ will denote the discrete value of the control ${\textstyle \mathbf {u} }$ at time ${\textstyle n\Delta t}$ and iteration ${\textstyle q}$, the conjugate gradient algorithm (38)-(44) to solve the finite dimensional problem (60) reads as follow:

• Suppose

$\mathbf {u} _{0}=\{\{u_{i0}^{n}\}_{i=1}^{3}\}_{n=1}^{N}{\hbox{ is given in }}{\mathcal {U}}_{ad}^{\Delta t}=(\mathbb {R} ^{3})^{N}{\hbox{ (}}\mathbf {u} _{0}=\mathbf {0} {\hbox{ for example)}}.$

(80)

• Compute ${\textstyle \{\mathbf {y} _{0}^{n}\}_{n=0}^{N}=\{\{y_{i0}^{n}\}_{i=1}^{3}\}_{n=0}^{N}}$ and ${\textstyle \{\mathbf {p} _{0}^{n}\}_{n=1}^{N+1}=\{\{p_{i0}^{n}\}_{i=1}^{3}\}_{n=1}^{N+1}}$ via the solution of

$\left\{{\begin{array}{c}\mathbf {y} _{0}^{0}=\delta {\boldsymbol {\theta }},\\{\hbox{for }}n=1,....,N{\hbox{ solve}}\\\Gamma {\frac {\displaystyle \mathbf {y} _{0}^{n}-\mathbf {y} _{0}^{n-1}}{\displaystyle \Delta t}}+(K+C)\mathbf {y} _{0}^{n}=\mathbf {u} _{0}^{n},\end{array}}\right.$

(81)

and then

$\left\{{\begin{array}{c}\Gamma \mathbf {p} ^{N+1}=k_{2}\mathbf {y} ^{N},\\{\hbox{for }}n=N,....,1{\hbox{ solve}}\\\Gamma {\frac {\displaystyle \mathbf {p} _{0}^{n}-\mathbf {p} _{0}^{n+1}}{\displaystyle \Delta t}}+(K+C)\mathbf {p} _{0}^{n}=k_{1}\mathbf {y} _{0}^{n}\mathbf {.} \end{array}}\right.$

(82)

• Set

$\mathbf {g} _{0}=\{\mathbf {g} _{0}^{n}\}_{n=1}^{N}=\{\mathbf {u} _{0}^{n}+\mathbf {p} _{0}^{n}\}_{n=1}^{N}.$

(83)

• If

${\frac {\Delta t\sum _{n=1}^{N}\left\vert \mathbf {g} _{0}^{n}\right\vert ^{2}}{\max[1,\Delta t\sum _{n=1}^{N}\left\vert \mathbf {u} _{0}^{n}\right\vert ^{2}]}}\leq tol^{2},{\hbox{ with }}\left\vert \mathbf {g} _{0}^{n}\right\vert ^{2}=\left\vert g_{10}^{n}\right\vert ^{2}+\left\vert g_{20}^{n}\right\vert ^{2}+\left\vert g_{30}^{n}\right\vert ^{2}{\hbox{,}}$

take ${\textstyle \mathbf {u} ^{\Delta t}=\mathbf {u} _{0}}$; otherwise, set

$\mathbf {w} _{0}=\mathbf {g} _{0}.$

(84)

For ${\textstyle q\geq 0}$, ${\textstyle \mathbf {u} _{q}}$, ${\textstyle \mathbf {g} _{q}}$ and ${\textstyle \mathbf {w} _{q}}$ being known, the last two different from ${\textstyle \mathbf {0} }$, we compute ${\textstyle \mathbf {u} _{q+1}}$, ${\textstyle \mathbf {g} _{q+1}}$ and, if necessary ${\textstyle \mathbf {w} _{q+1}}$ as follows:

• Solve

$\left\{{\begin{array}{c}{\overline {\mathbf {y} }}_{q}^{0}=\mathbf {0} ,\\{\hbox{for }}n=1,....,N{\hbox{ solve}}\\\Gamma {\frac {\displaystyle {\overline {\mathbf {y} }}_{q}^{n}-{\overline {\mathbf {y} }}_{q}^{n-1}}{\displaystyle \Delta t}}+(K+C){\overline {\mathbf {y} }}_{q}^{n}=\mathbf {w} _{q}^{n},\end{array}}\right.$

(85)

and then

$\left\{{\begin{array}{c}\Gamma {\overline {\mathbf {p} }}_{q}^{N+1}=k_{2}{\overline {\mathbf {y} }}_{q}^{N},\\{\hbox{for }}n=N,....,1{\hbox{ solve}}\\\Gamma {\frac {\displaystyle {\overline {\mathbf {p} }}_{q}^{n}-{\overline {\mathbf {p} }}_{q}^{n+1}}{\displaystyle \Delta t}}+(K+C){\overline {\mathbf {p} }}_{q}^{n}=k_{1}{\overline {\mathbf {y} }}_{q}^{n}\mathbf {.} \end{array}}\right.$

(86)

Set

${\overline {\mathbf {g} }}_{q}=\mathbf {w} _{q}+{\overline {\mathbf {p} }}_{q},$

(87)

and

$\rho _{q}=\Delta t\sum _{n=1}^{N}\left\vert \mathbf {g} _{q}^{n}\right\vert ^{2}/(\Delta t\sum _{n=1}^{N}{\overline {\mathbf {g} }}_{q}^{n}\mathbf {\cdot w} _{q}^{n}).$

(88)

• Compute

$\mathbf {u} _{q+1}=\mathbf {u} _{q}-\rho _{q}\mathbf {w} _{q},$

(89)

and

$\mathbf {g} _{q+1}=\{\mathbf {g} _{q+1}^{n}\}_{n=1}^{N}=\mathbf {g} _{q}-\rho _{q}{\overline {\mathbf {g} }}_{q}.$

(90)

• If

${\frac {\Delta t\sum _{n=1}^{N}\left\vert \mathbf {g} _{q+1}^{n}\right\vert ^{2}}{\max[\Delta t\sum _{n=1}^{N}\left\vert \mathbf {g} _{0}^{n}\right\vert ^{2},\Delta t\sum _{n=1}^{N}\left\vert \mathbf {u} _{q+1}^{n}\right\vert ^{2}]}}\leq tol^{2},$

take ${\textstyle \mathbf {u} ^{\Delta t}=\mathbf {u} _{q+1}}$; otherwise, compute

$\gamma _{q}={\frac {\sum _{n=1}^{N}\left\vert \mathbf {g} _{q+1}^{n}\right\vert ^{2}}{\sum _{n=1}^{N}\left\vert \mathbf {g} _{q}^{n}\right\vert ^{2}}},$

(91)

and set

$\mathbf {w} _{q+1}=\mathbf {g} _{q+1}+\gamma _{q}\mathbf {w} _{q}.$

(92)

• Do ${\textstyle q+1\longrightarrow q}$ and return to (85) .

• End of algorithm

Similar to the continuous case, we can deduce that if ${\textstyle \mathbf {y} _{q}=\mathbf {y} (\mathbf {u} _{q})}$ then

$\mathbf {y} _{q+1}-\mathbf {y} _{q}=-\rho _{q}{\overline {\mathbf {y} }}_{q}\Longrightarrow \mathbf {y} _{q+1}=\mathbf {y} _{q}-\rho _{q}{\overline {\mathbf {y} }}_{q}.$

4 Numerical Results

In previous sections we have described a methodology and the respective practical algorithms to use a control on each joint in order to stabilize the linear JJAM perturbation system around an unstable equilibrium. It is easy to simplify the procedure and algorithm to the case when we want to control via only (any combina-tion of) two junctions or via only one junction. However, since (according to the experiments) it is necessary to control via at least two junctions in order to stabilize the system around an unstable equilibrium, we show only the results when two and three junctions are used to control the system model. For the calculations we used ${\textstyle tol^{2}=10^{-16}}$ for the stopping criteria in conjugate gradient algorithm, and ${\textstyle \Delta t=10^{-2}}$ for solving the differential systems. In the next two subsections we use as ${\textstyle \theta }$ the equilibrium given by ${\textstyle \{n_{1},n_{2},n_{3}\}=\{1,0,0(u)\}}$ and ${\textstyle \delta {\boldsymbol {\theta }}=[1e-2,1e-2,1e-2]}$; the time interval under consideration being ${\textstyle [0,2]}$. Finally we apply iteratively this control process to stabilize in the longer interval ${\textstyle [0,T]}$ using 10 subintervals of length 2.

4.1 Controlling via two junctions

When controlling via two junctions only the case when using junctions 2 and 3 is a successful one (for all values of ${\textstyle k_{1}}$ and ${\textstyle k_{2}}$). In Figures 7 and 8 are shown the respective results.

Figure 7: $\mathbf {u} ^{\Delta t}$ (left) and $\|\mathbf {y} ^{\Delta t}(\cdot )\|$ (right) for several values of $k_{1}$ and $k_{2}$. The unstable equilibrium ${\boldsymbol {\theta }}$ is given by $\{n_{1},n_{2},n_{3}\}=\{1,0,0(u)\}$. The junctions used to control are 2 and 3.

Figure 8: $Ln\|\mathbf {g} _{q}^{\Delta t}\|$ for several values of $k_{1}$ and $k_{2}$. The unstable equilibrium ${\boldsymbol {\theta }}$ is given by $\{n_{1},n_{2},n_{3}\}=\{1,0,0(u)\}$. The junctions used to control are 2 and 3.

4.2 Controlling via three junctions

Figures 9 and 10 show the results when the three junctions are used to control. As we can see, for all values of the penalty parameters, the linear system is controlled.

Figure 9: $\mathbf {u} ^{\Delta t}$ (left) and $\|\mathbf {y} ^{\Delta t}(\cdot )\|$ (right) for several values of $k_{1}$ and $k_{2}$. The unstable equilibrium ${\boldsymbol {\theta }}$ is given by $\{n_{1},n_{2},n_{3}\}=\{1,0,0(u)\}$. The junctions used to control are 1, 2 and 3.

Figure 10: $Ln\|\mathbf {g} _{q}^{\Delta t}\|$ for several values of $k_{1}$ and $k_{2}$. The unstable equilibrium ${\boldsymbol {\theta }}$ is given by $\{n_{1},n_{2},n_{3}\}=\{1,0,0(u)\}$. The junctions used to control are 1, 2 and 3.

For the particular values ${\textstyle k_{1}=1E+3}$ and ${\textstyle k_{2}=1E+3}$ of the penalty parameters we show in Figure 11 the controls and the norm of the solution of the nonlinear system (the behavior of the solution of the linear system is shown in Figure 9)

Figure 11: Optimal controls (left), and Euclidean norm of the controlled solution ${\boldsymbol {\phi }}^{\Delta t}$ of the nonlinear system (right).

In Figure 12 we show the solution for the linear and nonlinear model, using controls in Figure 11 (we continue from ${\textstyle t=2}$ to ${\textstyle t=20}$ with no control).

Figure 12: Extended controlled solution $y_{i}^{\Delta t}$ of the linear perturbation model (left), and extended controlled solution $\phi _{i}^{\Delta t}$ of the nonlinear system (right). After $t=2$ the controls are extended as cero.

4.3 Controlling via three junctions in the interval [0,20]

To control during the whole time interval ${\textstyle (0,20)}$ we have divided the time interval into subintervals of smaller length ${\textstyle \Delta T=T/Q=2}$, and we denote ${\textstyle q\Delta T}$ by ${\textstyle T_{q}}$ for ${\textstyle q=1,...,Q}$; we proceed then as follows:

• For ${\textstyle q=0}$, we denote by ${\textstyle \mathbf {y} _{0}}$ the difference ${\textstyle {\boldsymbol {\phi }}_{0}-{\boldsymbol {\theta }}}$, and we solve the associated linear control problem (9)-(11) in ${\textstyle [0,T_{1}]}$; let us denote by ${\textstyle \mathbf {u} _{1}}$ the corresponding control. This control is injected in (5) with initial condition (8) and ${\textstyle \delta {\boldsymbol {\theta }}=\mathbf {y} _{0}}$, and we denote by ${\textstyle \mathbf {y} _{1}}$ the difference ${\textstyle {\boldsymbol {\phi }}(T_{1})-{\boldsymbol {\theta }}}$.
• For ${\textstyle q>0}$, we denote by ${\textstyle \mathbf {y} _{q}}$ the difference ${\textstyle {\boldsymbol {\phi }}(T_{q})-{\boldsymbol {\theta }}}$; we solve the associated linear control problem (9)-(11) in ${\textstyle [T_{q},T_{q+1}]}$, with ${\textstyle \mathbf {y} _{0}}$ replaced by ${\textstyle \mathbf {y} _{q}}$, and we denote by ${\textstyle \mathbf {u} _{q+1}}$ the corresponding optimal control. The control ${\textstyle \mathbf {u} _{q+1}}$ is injected in (5) with initial condition (8) and ${\textstyle \delta {\boldsymbol {\theta }}=\mathbf {y} _{q}}$, and we denote by ${\textstyle \mathbf {y} _{q+1}}$ the difference ${\textstyle {\boldsymbol {\phi }}(T_{q+1})-{\boldsymbol {\theta }}}$.
• We do ${\textstyle q=q+1}$ and we repeat the process.

The above time partitioning method has been applied with ${\textstyle {\boldsymbol {\phi }}_{0}={\boldsymbol {\theta }}+\delta {\boldsymbol {\theta }}}$ and ${\textstyle \delta {\boldsymbol {\theta }}=[1e-2,1e-2,1e-2]}$, the time interval under consideration being ${\textstyle [0,20]}$; we have used ${\textstyle \Delta T=2.0}$. After ${\textstyle t=20}$, we have taken ${\textstyle \mathbf {v} =\mathbf {0} }$ in (9) and in (5) to observe the evolution of the suddenly uncontrolled linear and nonlinear systems. The results are reported in Figure 13. We observe that the system is practically stabilized for ${\textstyle 1\leq t\leq 20}$, but if one stops controlling, the small residual perturbations of the system at ${\textstyle t=20}$, are sufficient to destabilize the linear and nonlinear systems and induces the nonlinear one to transition to a stable equilibrium in finite time.

Figure 13: Controls calculated each two seconds (top-left); Euclidean norm of the controlled solution $\mathbf {y} ^{\Delta t}$ using the controls in the top (top-right); Euclidean norm of the controlled solution ${\boldsymbol {\phi }}^{\Delta t}$ using the controls in the top-left (bottom).

5 Conclusions

We have stabilized the phases of a JJAM model, around an unstable equilibrium by using the classical approach: linearize the state model around the unstable equilibrium; control the linear model in order to stabilize it around the unstable equilibrium; apply the linear control to the nonlinear model and hope this control will also stabilize it. Since the time interval is large (${\textstyle t\in [0,30]}$) in certain applications, and could be that the linear control do not stabilize the nonlinear model, we subdivi-ded the original interval into subintervals and calculate iteratively the linear control on each subinterval, obtaining a piecewise control that stabilize not only the linear model but also the nonlinear one. For an efficient calculation of the control we formulated an operational linear equation satisfied by the control. The associated operator is self-adjoint and elliptic, so a conjugate gradient algorithm for quadratic functional was used.

BIBLIOGRAPHY

[1] Y. Braiman, B. Neschke, N. Nair, N. Ima, and R. Glowinski. (2016) Memory States in Small Arrays of Josephson Junctions, PHYSICAL REVIEW E (94), 052223: 1-13.

[2] Y. Braiman and N. Nair and J. Rezac and N. Imam. (2016) Memory Cell Operation Based on Small Josephson Junction Arrays, Superconductor Science and Technology (129), 124003: 1-15.

[3] J.D. Rezac and N. Imam and Y. Braiman. (2017) Parameter optimization for transitions between memory states in small arrays of Josephson junctions, PHYSICA A (474), 267-281.

[4] Harvey, Roland and Qu, Zhihua. (2018) Control of Cryogenic Memory State Transitions in a Josephson Junction Array, 2018 Annual American Control Conference (ACC), 5671-5676.

[5] N. Nair and Y. Braiman. (2018) A Ternary Memory Cell Using Small Josephson Junction Arrays, Superconductor Science and Technology (31).

[6] F. Bedard and N. Welker and G. R. Cotter and M. A. Escavage and J. T. Pinkston. (2005) Superconducting Technology Assessment, National Security Agency, Office of Corporate Assessment, url = https://www.nitrd.gov/pubs/nsa/sta.pdf.

[7] F. A. Holmes and L. Ripple and M. A. Manheimer. (2013) Energy-Efficient Superconducting Computing-Power Budgets and Requirements, IEEE Transactions on Applied Superconductivity , 23 (3).

[8] IARPA. (2013) Broad Agency Announcement: IARPA-BAA-13-05, Cryogenic Computing Complexity (C3) Program, IARPA.

[9] Roland Glowinski and Jorge López and Héctor Juárez and Yehuda Braiman. (2020) On the controllability of transitions between equilibrium states in small inductively coupled arrays of Josephson junctions: A computational approach, Journal of Computational Physics (403), 109023, doi = https://doi.org/10.1016/j.jcp.2019.109023, url = http://www.sciencedirect.com/science/article/pii/S0021999119307296.

[10] R. Glowinski and J.L. Lions and J. W. He. (2008) Exact and Approximate Controllability for Distributed Parameter Systems, Cambridge University Press, Cambridge, UK.

[11] J. L. Lions. (1971) Optimal Control of Systems Governed by Partial Differential Equations, Springer Verlag, New York.