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

<!-- metadata commented in wiki content
==Stabilizing a Josephson Junction Array  Memory  around an unstable equilibrium:   A control approach using a first order state model ==

'''Jorge López López 

Universidad Juárez Autónoma de Tabasco, México.

Correspondig author: jorge.lopez@hotmail.com



Héctor Juárez Valencia.

Universidad Autónoma Metropolitana, Iztapalapa

hect@xanum.uam.mx.



Roland Glowinski.'''
-->

==Abstract==

'''In this paper we consider a system of three nonlinear ordinary diffe-rential equations that model a three Josephson Junctions Array (JJA). Our goal is to stabilize the system around an unstable equili-brium 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 Algo-rithms, Stabilization, Optimal Control.

MSC[2010] 65K10, 49M25.

==2 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 (<span id='citeF-1'></span><span id='citeF-2'></span><span id='citeF-3'></span><span id='citeF-4'></span><span id='citeF-5'></span>[[#cite-1|[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:

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\{ \begin{array}{c}V = \frac{\displaystyle \hbar }{\displaystyle 2e}\frac{\displaystyle d\phi }{\displaystyle dt}, \\ I = I_{c}\sin \phi , \end{array} \right. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (1)
|}

where <math display="inline">\hbar = h / 2\pi </math>, <math display="inline">h</math> being the Planck constant, <math display="inline">-e</math> is the electric charge of the electron, <math display="inline">V</math> and <math display="inline">I</math> are the voltage and current across the junction, respectively, <math display="inline">\phi </math> is the phase difference across the junction, and, finally, the constant <math display="inline">I_{c}</math> 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 <span id='citeF-6'></span>[[#cite-6|[6]]]. Indeed, (cf. <span id='citeF-7'></span>[[#cite-7|[7]]]), ``Single flux quantum (SFQ) circuits produce small current pulses that travel at about <math display="inline">1/3</math> 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. <span id='citeF-8'></span>[[#cite-8|[8]]]) ``Josephson junctions, the superconducting switching devices, switch quickly (<math display="inline">\sim 1</math> ps), dissipate very little energy per switch (<math display="inline"><10^{-19}J</math>), 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 (<span id='citeF-1'></span>[[#cite-1|[1]]], <span id='citeF-3'></span>[[#cite-3|[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 (<math display="inline">t=1</math> corresponds to <math display="inline">4.15ps</math>)

<span id="eq-2"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
|}

where <math display="inline">I_{j}=i_{j}+ad_{j}</math>, <math display="inline">i_{j}</math> being a direct current and <math display="inline">ad_{j}</math> an additional energy as a current pulse. Plausible values for the various quantities in the model are:

<span id="eq-3"></span>
<span id="eq-4"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\gamma _{1} =0.7,\gamma _{2}=1.1,\gamma _{3}=0.7,i_{1}=1,i_{2}=0.8,i_{3}=-1,  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
|-
| style="text-align: center;" | <math> \kappa _{1} =\kappa _{2}=0.1.  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (4)
|}
|}

For those regimes where the second order derivatives can be neglected (meaning the junctions are (relatively) highly dissipative), ([[#eq-2|2]]) reduces to

<span id="eq-5"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
|}

The ''Read/Write'' operations can be performed by applying appropriate Gaussian pulses (<span id='citeF-1'></span>[[#cite-1|[1]]], <span id='citeF-4'></span>[[#cite-4|[4]]], <span id='citeF-5'></span>[[#cite-5|[5]]] ) in order to drive the system from an equilibrium configuration to another one.  In  <span id='citeF-9'></span>[[#cite-9|[9]]], we went beyond Gaussian pulses and investigated an optimal  controllability approach in <math display="inline">L^2</math> to perform these transitions between  equilibrium confi-gurations. These approach is closely related to the methodology discussed in <span id='citeF-10'></span>[[#cite-10|[10]]] for systems modelled by partial differential equations.

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

<span id="eq-6"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\theta _{k}=2\pi n_{k}+\sigma _{k}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
|}

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


{|  class="floating_tableSCP wikitable" style="text-align: center; margin: 1em auto;min-width:50%;"
|+ style="font-size: 75%;" |<span id='table-1'></span>Table. 1 Stable and unstable equilibrium configurations for the JJAM system defined by (2) and (3). The symbols <math>(u)</math> and <math>(s)</math> qualify the unstable and stable equilibria, respectively.
|- style="border-top: 2px solid;"
| colspan='1' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>n_{1}</math>
| colspan='1' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>n_{2}</math>
| colspan='1' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>n_{3}</math>
| colspan='1' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>\sigma _{1}/2\pi </math>
| colspan='1' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>\sigma _{2}/2\pi </math>
| colspan='1' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>\sigma _{3}/2\pi </math>
| colspan='1' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>\theta _{1}</math>
| colspan='1' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>\theta _{2}</math>
| colspan='1' style="border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>\theta _{3}</math>
|- style="border-top: 2px solid;"
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0 (s)</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0.1992</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0.1187</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>-0.1552</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>1.2517</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.7458</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.9752</math> 
|- style="border-top: 2px solid;"
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>1</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0 (s)</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | 
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | 
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | no
| style="border-left: 2px solid;border-right: 2px solid;" | stable 
| style="border-left: 2px solid;border-right: 2px solid;" | equil 
| style="border-left: 2px solid;border-right: 2px solid;" |  
|- style="border-top: 2px solid;"
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>1</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>1</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0 (s)</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0.1810</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0.0338</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>-0.0515</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>7.4207</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>6.4958</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.3236</math> 
|- style="border-top: 2px solid;"
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>2</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0 (s)</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | 
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | 
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | no
| style="border-left: 2px solid;border-right: 2px solid;" | stable   
| style="border-left: 2px solid;border-right: 2px solid;" |  equil  
| style="border-left: 2px solid;border-right: 2px solid;" |    
|- style="border-top: 2px solid;"
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>2</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>1</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0 (s)</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0.0661</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0.1164</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>-0.0437</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>12.9821</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>7.0148</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>-0.2746</math> 
|- style="border-top: 2px solid;"
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>2</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>2</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0 (s)</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0.1671</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>-0.0443</math>
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | <math>0.0333</math>
| style="border-left: 2px solid;border-right: 2px solid;" | <math>13.6161</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>12.2884</math> 
| style="border-left: 2px solid;border-right: 2px solid;" | <math>0.2092</math> 
|- style="border-top: 2px solid;"
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | -
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | -
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | -
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | -
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | -
| colspan='1' style="text-align: left;border-left: 2px solid;border-right: 2px solid;border-left: 2px solid;border-right: 2px solid;" | -
| style="border-left: 2px solid;border-right: 2px solid;" | - 
| style="border-left: 2px solid;border-right: 2px solid;" | - 
| style="border-left: 2px solid;border-right: 2px solid;" | - 

|}
For the set of dc currents and coupling parameters defined by ([[#eq-3|3]]) and ([[#eq-4|4]]), no stable  steady state exists when the first junction phase is about <math display="inline">2\pi </math> (respectively <math display="inline">4\pi </math>) greater than both of the other two junction phases (see row where <math display="inline">\{ n_{1},n_{2},n_{3}\} =\{ 1,0,0\} </math> (respectively <math display="inline">\{ n_{1},n_{2},n_{3}\} =\{ 2,0,0\} </math>)).  Following Braiman,  et. al., <span id='citeF-1'></span>[[#cite-1|[1]]], for memory cell demonstration it is sufficient to manipulate the system within the particular sets of states, namely states <math display="inline">\{ 0,0,0\} </math> and <math display="inline">\{ 1,1,0\} </math>. The circuit of three coupled Rapid Single Joint (RSJ) unit circuits associated with model ([[#eq-5|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 <math display="inline">n_{k}</math> as presented in equation ([[#eq-6|6]]), will describe the location of the <math display="inline">k^{th}</math> junction phase. When all three junction phases are in the same sinusoidal potential well, the system will be considered in the `0' state, <math display="inline">\{ n_{1},n_{2},n_{3}\} =\{ 0,0,0\} .</math> When the phases of the first and second junctions are shifted to the next potential well (about <math display="inline">2\pi </math> greater), the cell will be in the `1' state, <math display="inline">\{ n_{1},n_{2},n_{3}\} =\{ 1,1,0\} </math>. These two states correspond to the first and third rows of Table [[#table-1|1]].  Figure [[#img-1|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 <math display="inline">2\pi </math>. For convenience, in several of the following figures the phases are normalized by <math display="inline">2\pi </math> and in all the graphs that include different colors, color blue, red and green is associated with junction 1, 2 and 3, respectively.

<div id='img-1'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_LOPEZ_262416069-cfigure2.png|307px|Time series of the phases scaled by 2π, with zero initial conditions.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 1:''' Time series of the phases scaled by <math>2\pi </math>, with zero initial conditions.
|}

In this work we investigate a controllability approach  in order to stabilize the phase junctions of ([[#eq-5|5]]) around an unstable equilibrium. Figures [[#img-2|2]] and [[#img-3|3]] show the behavior of the system when the initial conditions are approximations of the unstable equilibrium <math display="inline">\{ n_{1},n_{2},n_{3}\} </math>=<math display="inline">\{ 1,0,0(u)\} </math>, <math display="inline">\{ 2,1,0(u)\} </math>, <math display="inline">\{ 2,2,0(u)\} </math>.

<div id='img-2'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_LOPEZ_262416069-cfigure3.png|300px|]]
|[[Image:Draft_LOPEZ_262416069-cfigure4.png|300px|Evolution to state \{ 1,1,0(s)\} (top), \{ 2,1,0(s)\} (bottom) from an approximation to the unstable equilibrium \{ 1,0,0(u)\} (top), \{ 2,1,0(u)\} (bottom).]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 2:''' Evolution to state <math>\{ 1,1,0(s)\} </math>(top), <math>\{ 2,1,0(s)\} </math>(bottom) from an approximation to the unstable equilibrium <math>\{ 1,0,0(u)\} </math>(top), <math>\{ 2,1,0(u)\} </math>(bottom).
|}

<div id='img-3'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_LOPEZ_262416069-cfigure5.png|307px|Evolution to \{ 2,2,0(s)\}  from an approximation to the unstable equilibrium \{ 2,2,0(u)\} .]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 3:''' Evolution to <math>\{ 2,2,0(s)\} </math> from an approximation to the unstable equilibrium <math>\{ 2,2,0(u)\} </math>.
|}

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.

==3 Formulation of the optimal control problem for the stabiliza-tion of the JJAM and the conjugate gradient algorithm==

===3.1 The approach===

As we explained in the Introduction, Figures [[#img-2|2]] and [[#img-3|3]] show  the behavior of the system ([[#eq-5|5]]) when the initial conditions are approximations of the unstable  equilibrium configurations <math display="inline">\{ n_{1},n_{2},n_{3}\} </math>= <math display="inline">\{ 1,0,0(u)\} </math>, <math display="inline">\{ 2,1,0(u)\} </math>, <math display="inline">\{ 2,2,0(u)\} </math>. Our goal in this  subsection is to stabilize the system ([[#eq-5|5]]), around an unstable equilibrium, controlling it via one, two or three junctions i.e., we want to maintain the phase values in Figures [[#img-2|2]] and [[#img-3|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 ([[#eq-5|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 <math display="inline">\boldsymbol{\theta } = \{ \theta _1,\theta _2,\theta _3  \} </math> of ([[#eq-5|5]]) and a small initial variation <math display="inline">\delta \boldsymbol{\theta } </math> of <math display="inline">\boldsymbol{\theta } </math>. The perturbation <math display="inline">\delta \boldsymbol{\phi } </math> of the steady-state <math display="inline">\theta </math> in <math display="inline">(0,T)</math> satisfies approximately the following linear model

<span id="eq-7"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
|}

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

<div id='img-4'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_LOPEZ_262416069-Ufigure1.png|280px|Solution y=δϕ  of ([[#eq-7|7]]) for the unstable equilibrium given by \{ n₁,n₂,n₃\} =\{ 1,0,0(u)\}  and δθ=[ 1e-2, 1e-2, 1e-2].]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 4:''' Solution <math>\mathbf{y}=\delta \boldsymbol{\phi } </math> of ([[#eq-7|7]]) for the unstable equilibrium given by <math>\{ n_{1},n_{2},n_{3}\} =\{ 1,0,0(u)\} </math> and <math>\delta \boldsymbol{\theta }=[ 1e-2, 1e-2, 1e-2]</math>.
|}

<div id='img-5'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_LOPEZ_262416069-Ufigure2.png|280px|Solution y=δϕ  of ([[#eq-7|7]]) for the unstable equilibrium given by \{ n₁,n₂,n₃\} =\{ 2,1,0(u)\}  and δθ=[ 1e-2, 1e-2, 1e-2].]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 5:''' Solution <math>\mathbf{y}=\delta \boldsymbol{\phi } </math> of ([[#eq-7|7]]) for the unstable equilibrium given by <math>\{ n_{1},n_{2},n_{3}\} =\{ 2,1,0(u)\} </math> and <math>\delta \boldsymbol{\theta }=[ 1e-2, 1e-2, 1e-2]</math>.
|}

<div id='img-6'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_LOPEZ_262416069-Ufigure3.png|280px|Solution y=δϕ  of ([[#eq-7|7]]) for the unstable equilibrium given by \{ n₁,n₂,n₃\} =\{ 2,2,0(u)\}  and δθ=[ 1e-2, 1e-2, 1e-2].]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 6:''' Solution <math>\mathbf{y}=\delta \boldsymbol{\phi } </math> of ([[#eq-7|7]]) for the unstable equilibrium given by <math>\{ n_{1},n_{2},n_{3}\} =\{ 2,2,0(u)\} </math> and <math>\delta \boldsymbol{\theta }=[ 1e-2, 1e-2, 1e-2]</math>.
|}

It is clear that the model ([[#eq-7|7]]) is no longer valid if <math display="inline">\delta \boldsymbol{\theta } </math> becomes too large but the idea behind considering the linearized model ([[#eq-7|7]]) is to use it to compute a control action preventing <math display="inline">\delta  \boldsymbol{\phi } </math> from becoming too large (and possibly driving <math display="inline">\delta \boldsymbol{\phi } </math> to zero) and hope that the computed control will also stabilize the original nonlinear system ([[#eq-5|5]]) with initial condition

<span id="eq-8"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\boldsymbol{\phi }(0) = \boldsymbol{\theta } + \delta \boldsymbol{\theta }.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
|}

===3.2 The Linear Control Problem===

Using the notation <math display="inline">\mathbf{y}</math> <math display="inline">=</math> <math display="inline">\delta \boldsymbol{\phi } </math>, we shall (try to) stabilize the controlled variant of system ([[#eq-7|7]]) (the three junctions are used to control):

<span id="eq-9"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
|}

via the following control formulation:

<span id="eq-10"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
|}

where

<span id="eq-11"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>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},  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
|}

and <math display="inline">|| \mathbf{y}|| ^{2}=\left\vert y_{1}\right\vert  ^{2}+\left\vert y_{2}\right\vert ^{2}+\left\vert y_{3}\right\vert ^{2}</math>.

===3.3 Optimality Conditions and Conjugate Gradient Solution for Problem ([[#eq-10|10]])===

====3.3.1 Generalities and Synopsis====

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

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> (\mathbf{v,w})_{\mathcal{U}}=\int \limits _{0}^{T}\mathbf{v(}t\mathbf{)\cdot  w(}t\mathbf{)}\hbox{ }dt,  </math>
|}
|}

the action <math display="inline">\left\langle DJ(\mathbf{v}),\mathbf{w}\right\rangle </math> of <math display="inline">DJ(\mathbf{v})</math> on <math display="inline">\mathbf{w}\in \mathcal{U}</math> can also be written as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \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{,}  </math>
|}
|}

with <math display="inline">DJ(\mathbf{v})\mathbf{(}t\mathbf{)}\in \mathcal{U}</math>. If <math display="inline">\mathbf{u}</math> is the solution of problem ([[#eq-10|10]]), it is characterized [from convexity arguments (see, e.g., <span id='citeF-11'></span>[[#cite-11|[11]]])] by

<span id="eq-12"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>DJ(\mathbf{u})=\mathbf{0.}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
|}

Since the cost function <math display="inline">J</math> is quadratic and since the state model ([[#eq-7|7]]) is linear, <math display="inline">DJ(\mathbf{v})</math> is in fact an affine function of <math display="inline">\mathbf{v}</math>, implying in turn from ([[#eq-12|12]]) that <math display="inline">\mathbf{u}</math> is the solution of a linear equation in the control space <math display="inline">\mathcal{U}</math>. In abstract form, equation ([[#eq-12|12]]) can be written as

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> DJ(\mathbf{u})-DJ(\mathbf{0})=-DJ(\mathbf{0}),  </math>
|}
|}

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

====3.3.2 <span id='lb-3.3.2'></span>Computing DJ(v)====

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

<span id="eq-13"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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} </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
|}

where in ([[#eq-13|13]]):

(i) We denote by '''a'''<math display="inline">\mathbf{\cdot }</math>'''b''' the dot product of two vectors '''a''' and '''b'''.

(ii) The function <math display="inline">\delta \mathbf{y}</math> <math display="inline">=\{ \delta y_{1},\delta y_{2},\delta  y_{3}\} </math> is the solution of the following initial value problem, obtained by perturbation of ([[#eq-9|9]]):

<span id="eq-14"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
|}

where <math display="inline">C_{1}=\cos \theta _{1}</math>, <math display="inline">C_{2}=\cos \theta _{2},</math> <math display="inline">C_{3}=\cos \theta  _{3}.</math> In matrix-vector form ([[#eq-14|14]]) reads as follows:

<span id="eq-15"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
|}

where <math display="inline">C</math> and <math display="inline">\Gamma </math> are diagonal matrices with coefficients <math display="inline">C_i</math> and <math display="inline">\gamma _i</math>, <math display="inline">i =1,2,3</math>, respectively, and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> 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).  </math>
|}
|}

We introduce now a vector-valued function <math display="inline">\mathbf{p}=\{ p_{1},p_{2},p_{3}\} </math> that we assume differentia-ble over <math display="inline">(0,T)</math>; multiplying by <math display="inline">\mathbf{p}</math> both sides of the differential equation in ([[#eq-15|15]]), and integrating over <math display="inline">(0,T)</math> we obtain, after integrating by parts, and taking into account the symmetry of the matrices <math display="inline">\Gamma </math> and <math display="inline">K</math>, the following equation

<span id="eq-16"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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).   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
|}

Now suppose that the vector-valued function <math display="inline">\mathbf{p}</math> is solution of the following ('''adjoint system'''):

<span id="eq-17"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
|}

It follows from ([[#eq-16|16]]), ([[#eq-17|17]]) that

<span id="eq-18"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>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.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
|}

Combining ([[#eq-13|13]]) and ([[#eq-18|18]]) we obtain that

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
|}

which implies in turn that

<span id="eq-20"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>DJ(\mathbf{v})=\mathbf{v}+\mathbf{p}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
|}

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

====3.3.3 Optimality conditions for problem ([[#eq-10|10]])====

Let <math display="inline">\mathbf{u}</math> be the solution of problem ([[#eq-10|10]]) and let us denote by <math display="inline">\mathbf{y}</math> (respectively, <math display="inline">\mathbf{p}</math>) the corresponding solution of the state system ([[#eq-9|9]]) (respectively, of the adjoint system ([[#eq-17|17]])). It follows from Subsubsection 2.3.2 that <math display="inline">DJ(\mathbf{u})=\mathbf{0}</math> is equivalent to the following (optimality) system:

<span id="eq-21"></span>
<span id="eq-22"></span>
<span id="eq-23"></span>
<span id="eq-24"></span>
<span id="eq-25"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{u} =-\mathbf{p},   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
|-
| style="text-align: center;" | <math> \Gamma \frac{d\mathbf{y}}{dt}+(K+C)\mathbf{y} =\mathbf{u}, \hbox{ in }(0,T),   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
|-
| style="text-align: center;" | <math> \mathbf{y}(0) =\delta \boldsymbol{\theta } ,   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
|-
| style="text-align: center;" | <math> -\mathbf{\Gamma }\frac{d\mathbf{p}}{dt}\hbox{ }+(K+C)\mathbf{p}\hbox{ } \mathbf{=}k_{1}\mathbf{y,}\hbox{ in }(0,T),   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
|-
| style="text-align: center;" | <math> \mathbf{\Gamma p(}T\mathbf{)} \mathbf{=}k_{2}\mathbf{y}(T).  </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
|}
|}

Conversely, it can be shown (see, e.g., <span id='citeF-11'></span>[[#cite-11|[11]]]) that the system ([[#eq-21|21]])&#8211;([[#eq-25|25]]) characterizes <math display="inline">\mathbf{u}</math> as the solution (necessarily unique here) of the control problem ([[#eq-10|10]]). The optimality conditions ([[#eq-21|21]])&#8211;([[#eq-25|25]]) will play a crucial role concerning the iterative solution of the control problem ([[#eq-10|10]]).

====3.3.4 Functional equations satisfied by the optimal control====

Our goal here is to show that the optimality condition <math display="inline">DJ(\mathbf{u})=\mathbf{0}</math> can also be written as

<span id="eq-26"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{Au}=\mathbf{\beta },   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
|}

where the linear operator <math display="inline">\mathbf{A}</math> is a strongly elliptic and symmetric isomorphism from <math display="inline">\mathcal{U}</math> into itself (an automorphism of <math display="inline">\mathcal{U)}</math> and where <math display="inline">\mathbf{\beta }\in \mathcal{U}</math>. A candidate for <math display="inline">\mathbf{A}</math> is the linear operator from <math display="inline">\mathcal{U}</math> into itself defined by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{Av}=\mathbf{v+p(v),}  </math>
|}
|}

where <math display="inline">\mathbf{p(v)}</math> is obtained from <math display="inline">\mathbf{v}</math> via the successive solution of the following two problems:

<span id="eq-27"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
|}

which is forward in time, and the adjoint system ([[#eq-17|17]]) in <math display="inline">(0,T)</math> which is backward in time. To see that operator <math display="inline">\mathbf{A}</math> is a symmetric and strongly elliptic isomorphism from <math display="inline">\mathcal{U}</math> into itself consider <math display="inline">\mathbf{v}^{1}</math>, <math display="inline">\mathbf{v}^{2}</math> belonging to <math display="inline">\mathcal{U}</math> and define <math display="inline">\mathbf{y}^{i}</math>, <math display="inline">\mathbf{p}^{i}</math> by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \mathbf{y}^{i}=\mathbf{y}(\mathbf{v}^{i}),\hbox{ }\mathbf{p}^{i}=\mathbf{p}(\mathbf{v}^{i}),\hbox{ }i=1,2\hbox{;}  </math>
|}
|}

we have then

<span id="eq-28"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
|}

On the other hand, we have, starting with the differential equation in ([[#eq-27|27]]) and using integration by parts, that

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>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 </math>
|-
| style="text-align: center;" | <math> =\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 </math>
|-
| style="text-align: center;" | <math> =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. </math>
|}
|}

This implies that

<span id="eq-29"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
|}

Combining ([[#eq-28|28]]) with ([[#eq-29|29]]) we obtain

<span id="eq-30"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
|}

Relation ([[#eq-30|30]]) imply the symmetry of <math display="inline">\mathbf{A}</math>; we also have

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \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{,}  </math>
|}
|}

which implies the strong ellipticity of <math display="inline">\mathbf{A}</math> over <math display="inline">\mathcal{U}</math>. The linear operator <math display="inline">\mathbf{A}</math>, being continuous and strongly elliptic over <math display="inline">\mathcal{U}</math>, is an automorphism of <math display="inline">\mathcal{U}</math>. To identify the right hand side <math display="inline">\mathbf{\beta }</math> of equation ([[#eq-26|26]]), we introduce <math display="inline">Y_{0}</math> and <math display="inline">P_{0}</math> defined as the solutions of

<span id="eq-31"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
|}

and

<span id="eq-32"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
|}

Suppose now that <math display="inline">\mathbf{y}</math> and <math display="inline">\mathbf{p}</math> satisfies the optimality conditions. Define

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\overline{\mathbf{y}} =\mathbf{y-Y}_{0}, </math>
|-
| style="text-align: center;" | <math> \overline{\mathbf{p}} =\mathbf{p-P}_{0}. </math>
|}
|}

Subtracting ([[#eq-31|31]]) and ([[#eq-32|32]]) from ([[#eq-27|27]]) and ([[#eq-17|17]]) we obtain

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
|}

and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
|}

From the definition of operator <math display="inline">\mathbf{A}</math> it follows that

<span id="eq-35"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{Au=u+}\overline{\mathbf{p}}=\mathbf{-p+}\overline{\mathbf{p}}=-\mathbf{P}_{0};   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
|}

the right hand side of ([[#eq-35|35]]) is the vector <math display="inline">\mathbf{\beta }</math> that we were looking for.

From the properties of <math display="inline">\mathbf{A}</math>, problem ([[#eq-35|35]]) can be solved by a conjugate gradient algorithm operating in the Hilbert space <math display="inline">\mathcal{U} </math>. This algorithm will be described in the following subsubsection.

====3.3.5 Conjugate gradient solution of the control problem ([[#eq-10|10]])====

Problem ([[#eq-35|35]]) can be written in variational form as

<span id="eq-36"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{u} \in \mathcal{U}\hbox{ }(=L^{2}(0,T;3)),  </math>
|-
| style="text-align: center;" | <math> \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{ .}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
|}

From the symmetry and <math display="inline">\mathcal{U}</math>-ellipticity of the bilinear form

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \left\{\mathbf{v},\mathbf{w}\right\}\rightarrow \int \limits _{0}^{T}(\mathbf{Av})\cdot \mathbf{w}\hbox{ }dt,  </math>
|}
|}

the variational problem ([[#eq-36|36]]) is a particular case of the following class of linear variational problems:

<span id="eq-37"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\left\{ \begin{array}{c}u\in \mathcal{V}\hbox{,} \\  a(u,v)=L(v),\forall v\in \mathcal{V}\hbox{,}\end{array}\right.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
|}

where

(i) <math display="inline">  \mathcal{V}</math> is a real Hilbert space for the scalar product <math display="inline">(\cdot  ,\cdot )</math> and the correspondent norm <math display="inline">\left\Vert \cdot \right\Vert ;</math>

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

(iii) <math display="inline">L</math> : <math display="inline">\mathcal{V\rightarrow  \mathbb{R} }</math> is linear and continuous.

If the above properties hold, then problem ([[#eq-37|37]]) has a unique solution (see <span id='citeF-10'></span>[[#cite-10|[10]]]); in fact, the symmetry of <math display="inline">a(\cdot ,\cdot )</math> is not necessary in order to have a unique solution to problem ([[#eq-37|37]]); however, the symmetry of <math display="inline">a(\cdot ,\cdot )</math>, combined with the other properties, allows problem ([[#eq-37|37]]) to be solved by the following conjugate gradient algorithm:

* '''Step 1'''. <math display="inline">u^{0}\in \mathcal{V}</math> is given.

* '''Step 2'''. Solve

<span id="eq-38"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
|}

and set

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>w^{0}=g^{0}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
|}

* '''Step 3'''. For <math display="inline">q\geq 0</math>, <math display="inline">u^{q}</math>, <math display="inline">g^{q}</math> and <math display="inline">w^{q}</math> being known, compute <math display="inline">u^{q+1}</math> , <math display="inline">g^{q+1}</math> and <math display="inline">w^{q+1}</math> as follows: Compute

<span id="eq-40"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\rho ^{q}=\left\Vert g^{q}\right\Vert ^{2}/a(w^{q},w^{q}),   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
|}

and take

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

u^{q+1}=u^{q}-\rho ^{q}w^{q}. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
|}

Solve

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
|}

and compute

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\gamma ^{q}=\left\Vert g^{q+1}\right\Vert ^{2}/\left\Vert g^{q}\right\Vert  ^{2}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
|}

<span id="eq-44"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

w^{q+1}=g^{q+1}+\gamma ^{q}w^{q}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
|}

* '''Step 4'''. Do n=n+1 and go to ([[#eq-40|40]]).

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

* Suppose

<span id="eq-45"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{u}^{0}\hbox{ is given in }L^{2}(0,T;3)\hbox{ }(\mathbf{u}^{0}=\mathbf{0}\hbox{ for example}).   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
|}

* Solve

<span id="eq-46"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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.,   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
|}

and then

<span id="eq-47"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
|}

* Set

<span id="eq-48"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{g}^{0}=\mathbf{u}^{0}+\mathbf{p}^{0}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
|}

* If <math display="inline">\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],</math> take <math display="inline">\mathbf{u}=\mathbf{u}^{0}</math>; otherwise, set

<span id="eq-49"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{w}^{0}=\mathbf{g}^{0}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
|}

For <math display="inline">q\geq 0</math>, <math display="inline">\mathbf{u}^{q}</math>, <math display="inline">\mathbf{g}^{q}</math> and <math display="inline">\mathbf{w}^{q}</math> being known, we compute <math display="inline">\mathbf{u}^{q+1},</math> <math display="inline">\mathbf{g}^{q+1}</math> and if necessary, <math display="inline">\mathbf{w}^{q+1}</math> as follows:

* Solve

<span id="eq-50"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
|}

and then

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
|}

Set

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\overline{\mathbf{g}}^{q}=\mathbf{w}^{q}+\overline{\mathbf{p}}^{q}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (52)
|}

and

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\rho _{q}=\left\Vert \mathbf{g}^{q}\right\Vert ^{2}/\int \limits _{0}^{T}\overline{\mathbf{g}}^{q}\mathbf{\cdot w}^{q}dt. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
|}

* Set

<span id="eq-54"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{u}^{q+1}=\mathbf{u}^{q}-\rho _{q}\mathbf{w}^{q}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
|}

* Set

<span id="eq-55"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{g}^{q+1}=\mathbf{g}^{q}-\rho _{q}\overline{\mathbf{g}}^{q}.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
|}

* If <math display="inline">\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},</math> take <math display="inline">\mathbf{u}=\mathbf{u}^{q+1}</math>; otherwise, compute

<span id="eq-56"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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},  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
|}

and set

<span id="eq-57"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{w}^{q+1}=\mathbf{g}^{q+1}+\gamma _{q}\mathbf{w}^{q}.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
|}

* Do <math display="inline">q+1\longrightarrow q</math> and return to ([[#eq-50|50]]) .

* End of algorithm.

The state variable <math display="inline">\mathbf{y}^{q}</math> can be actualized simultaneously with the control <math display="inline">\mathbf{u}^{q}</math> since for all <math display="inline">q</math> we have

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
|}

so,

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
|}

which, by definition of <math display="inline">\overline{\mathbf{y}}^{q},</math>  implies that

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \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}.  </math>
|}
|}

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

==4 Discrete formulation of the optimal control problem==

===4.1 <span id='lb-4.1'></span>Finite difference approximation of problem ([[#eq-10|10]]) when \mathcalU=L²(0,T;3)===

We approximate ([[#eq-10|10]]) when <math display="inline">\mathcal{U}=(L^{2}(0,T))^{3}</math> by

<span id="eq-60"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
|}

where:

* <math display="inline">\Delta t=T/N</math> with <math display="inline">N</math> a "large" positive integer.

* <math display="inline">\mathcal{U}^{\Delta t}=(\mathbb{R} ^{3})^{N}.</math>

The cost functional <math display="inline">J^{\Delta t}</math> is defined by

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> 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},  </math>
|}
|}

with <math display="inline">\mathbf{v}=\{ \mathbf{v}^{n}\} _{n=1}^{N}</math> and <math display="inline">\{ \mathbf{y}^{n}\} _{n=1}^{N}</math> obtained from <math display="inline">\mathbf{v}</math> and <math display="inline">\delta \boldsymbol{\theta }</math> via the following discrete variant of ([[#eq-9|9]]):

<span id="eq-61"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{y}^{0}=\delta \boldsymbol{\theta } ,   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
|}

and for <math display="inline">n=1,....,N</math>

<span id="eq-62"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Gamma \frac{\mathbf{y}^{n}-\mathbf{y}^{n-1}}{\Delta t}+(K+C)\mathbf{y}^{n}=\mathbf{v}^{n}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
|}

To compute <math display="inline">\mathbf{y}^{n}</math> we have thus to solve a linear system of the following type:

<span id="eq-63"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>(\Gamma +\Delta t(K+C))\mathbf{y}^{n}=RHS^{n}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (63)
|}

The matrix <math display="inline">\Gamma +\Delta t(K+C)</math>  being a <math display="inline">3\times 3</math> matrix symmetric and positive definite, solving ([[#eq-63|63]]) is easy.

===4.2 Optimality Conditions and conjugate gradient solution of ([[#eq-60|60]])===

====4.2.1 <span id='lb-4.2.1'></span>Computing DJ<sup>∆t</sup>(v)====

Assuming that one wants to use the conjugate gradient algorithm ([[#eq-38|38]])-([[#eq-44|44]]) to solve the discrete problem ([[#eq-60|60]]), we compute first <math display="inline">DJ^{\Delta t}(\mathbf{v})</math>. On <math display="inline">\mathcal{U}^{\Delta t}=(\mathbb{R} ^{3})^{N}</math> we will use the following inner-product

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> (\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}.  </math>
|}
|}

We have

<span id="eq-64"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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},  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (64)
|}

with <math display="inline">\{ \delta \mathbf{y}^{n}\} _{n=0}^{N}</math> obtained by perturbation of ([[#eq-61|61]]), ([[#eq-62|62]]), that is

<span id="eq-65"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\delta \mathbf{y}^{0}=\mathbf{0},   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
|}

and for <math display="inline">n=1,....,N</math>

<span id="eq-66"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Gamma \frac{\delta \mathbf{y}^{n}-\delta \mathbf{y}^{n-1}}{\Delta t}+(K+C)\delta \mathbf{y}^{n} =\mathbf{\delta v}^{n} .\hbox{ }   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
|}

Let us introduce <math display="inline">\{ \mathbf{p}^{n}\} _{n=1}^{N}\in (\mathbb{R} ^{3})^{N}.</math> Taking the dot product of <math display="inline">\mathbf{p}^{n}</math> with each side of equation ([[#eq-66|66]]), we obtain after summation and multiplication by <math display="inline">\Delta t:</math>

<span id="eq-67"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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{ }  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
|}

Applying a discrete integration by parts to relation ([[#eq-67|67]]), and considering that <math display="inline">\delta \mathbf{y}^{0}=\mathbf{0,}</math> we obtain:

<span id="eq-68"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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{ }  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
|}

or equivalently

<span id="eq-69"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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{ }  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
|}

Suppose that <math display="inline">\{ \mathbf{p}^{n}\} _{n=1}^{N+1}</math> verifies the following discrete adjoint system:

<span id="eq-70"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Gamma \mathbf{p}^{N+1}=k_{2}\mathbf{y}^{N},   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (70)
|}

and for <math display="inline">n=N,...,1</math>

<span id="eq-71"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\Gamma \frac{\mathbf{p}^{n}-\mathbf{p}^{n+1}}{\Delta t}+(K+C)\mathbf{p}^{n}=k_{1}\mathbf{y}^{n}\mathbf{.}   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (71)
|}

It follows from ([[#eq-64|64]]), ([[#eq-69|69]])-([[#eq-71|71]]) that

<span id="eq-72"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\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},  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (72)
|}

that is

<span id="eq-73"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>DJ^{\Delta t}(\mathbf{v})=\{ \mathbf{v}^{n}+\mathbf{p}^{n}\} _{n=1}^{N}.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (73)
|}

====4.2.2 Optimality conditions for ([[#eq-60|60]])====

The optimality conditions for the discrete problem ([[#eq-60|60]]) are

<span id="eq-74"></span>
<span id="eq-75"></span>
<span id="eq-76"></span>
<span id="eq-77"></span>
<span id="eq-78"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{u}^{n} =-\mathbf{p}^{n}\hbox{ , }n=1,...,N,   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (74)
|-
| style="text-align: center;" | <math> \mathbf{y}^{0} =\delta \boldsymbol{\theta } ,   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (75)
|-
| style="text-align: center;" | <math> \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,   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (76)
|-
| style="text-align: center;" | <math> \Gamma \mathbf{p}^{N+1} =k_{2}\mathbf{y}^{N},   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (77)
|-
| style="text-align: center;" | <math> \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.   </math>
| style="width: 5px;text-align: right;white-space: nowrap;" | (78)
|}
|}

====4.2.3 Functional equation for the discrete control solution of ([[#eq-60|60]])====

Following the sketch for the continuous case we can show that the discrete version <math display="inline">\mathbf{A}^{\Delta t}</math> of operator <math display="inline">\mathbf{A}</math> and the discrete version <math display="inline">\mathbf{\beta }^{\Delta t}</math> of <math display="inline">\mathbf{\beta }</math> satisfies the equation

<span id="eq-79"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>\mathbf{A}^{\Delta t}\mathbf{u}^{\Delta t}=\mathbf{\beta }^{\Delta t},  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (79)
|}

where <math display="inline">\mathbf{u}^{\Delta t}</math> is the discrete control satisfying the optimality condition ([[#eq-74|74]]). Operator <math display="inline">\mathbf{A}^{\Delta t}</math> enjoys the same properties than the continuous version: symmetric, strongly elliptic and continuous, allowing us to use a conjugate gradient like ([[#eq-38|38]])-([[#eq-44|44]]) to solve ([[#eq-79|79]]).

====4.2.4 Conjugate gradient solution of the discrete control problem ([[#eq-60|60]])====

Using <math display="inline">\mathbf{y}_{q}^{n}=\{ \mathbf{y}_{iq}^{n}\} _{i=1}^{3}</math> to denote the discrete value of the vector-valued function <math display="inline">\mathbf{y}</math> at time <math display="inline">n\Delta t</math> and iteration <math display="inline">q</math>; similarly, <math display="inline">\mathbf{u}_{q}^{n}</math> will denote the discrete value of the control <math display="inline">\mathbf{u}</math> at time <math display="inline">n\Delta t</math> and iteration <math display="inline">q</math>, the conjugate gradient algorithm ([[#eq-38|38]])-([[#eq-44|44]]) to solve the finite dimensional problem ([[#eq-60|60]]) reads as follow:

* Suppose

<span id="eq-80"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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)}.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (80)
|}

* Compute <math display="inline">\{ \mathbf{y}_{0}^{n}\} _{n=0}^{N}=\{ \{ \mathbf{y}_{i0}^{n}\} _{i=1}^{3}\} _{n=0}^{N}</math> and <math display="inline">\{ \mathbf{p}_{0}^{n}\} _{n=1}^{N+1}=\{ \{ \mathbf{p}_{i0}^{n}\} _{i=1}^{3}\} _{n=1}^{N+1}</math> via the solution of

<span id="eq-81"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (81)
|}

and then

<span id="eq-82"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (82)
|}

* Set

<span id="eq-83"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{g}_{0}=\{ \mathbf{g}_{0}^{n}\} _{n=1}^{N}=\{ \mathbf{u}_{0}^{n}+\mathbf{p}_{0}^{n}\} _{n=1}^{N}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (83)
|}

* If

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \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{,}  </math>
|}
|}

take <math display="inline">\mathbf{u}^{\Delta t}=\mathbf{u}_{0}</math>; otherwise, set

<span id="eq-84"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{w}_{0}=\mathbf{g}_{0}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (84)
|}

For <math display="inline">q\geq 0</math>, <math display="inline">\mathbf{u}_{q}</math>, <math display="inline">\mathbf{g}_{q}</math> and <math display="inline">\mathbf{w}_{q}</math> being known, the last two different from <math display="inline">\mathbf{0}</math>, we compute <math display="inline">\mathbf{u}_{q+1} </math>, <math display="inline">\mathbf{g}_{q+1}</math> and, if necessary <math display="inline">\mathbf{w}_{q+1}</math> as follows:

* Solve

<span id="eq-85"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (85)
|}

and then

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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. </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (86)
|}

Set

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\overline{\mathbf{g}}_{q}=\mathbf{w}_{q}+\overline{\mathbf{p}}_{q}, </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (87)
|}

and

<span id="eq-88"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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}).   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (88)
|}

* Compute

<span id="eq-89"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{u}_{q+1}=\mathbf{u}_{q}-\rho _{q}\mathbf{w}_{q},  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (89)
|}

and

<span id="eq-90"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{g}_{q+1}=\{ \mathbf{g}_{q+1}^{n}\} _{n=1}^{N}=\mathbf{g}_{q}-\rho _{q}\overline{\mathbf{g}}_{q}.   </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (90)
|}

* If

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \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},  </math>
|}
|}

take <math display="inline">\mathbf{u}^{\Delta t}=\mathbf{u}_{q+1}</math>; otherwise, compute

<span id="eq-91"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\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}},  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (91)
|}

and set

<span id="eq-92"></span>
{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math>

\mathbf{w}_{q+1}=\mathbf{g}_{q+1}+\gamma _{q}\mathbf{w}_{q}.  </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | (92)
|}

* Do <math display="inline">q+1\longrightarrow q</math> and return to ([[#eq-85|85]]) .

* End of algorithm

Similar to the continuous case, we can deduce that if <math display="inline">\mathbf{y}_{q}=\mathbf{y}(\mathbf{u}_{q})</math> then

{| class="formulaSCP" style="width: 100%; text-align: left;" 
|-
| 
{| style="text-align: left; margin:auto;width: 100%;" 
|-
| style="text-align: center;" | <math> \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}.  </math>
|}
|}

==5 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 <math display="inline">tol^2 = 10^{-16}</math> for the stopping criteria in conjugate gradient algorithm, and <math display="inline">\Delta t = 10^{-2}</math>  for solving the differential systems. In the next two subsections we use as <math display="inline">\theta </math> the equilibrium given by <math display="inline">\{ n_{1},n_{2},n_{3}\} =\{ 1,0,0(u)\} </math> and <math display="inline">\delta \boldsymbol{\theta }=[1e-2,1e-2,1e-2]</math>; the time interval  under consideration being <math display="inline">[0,2]</math>. Finally we apply iteratively this control process to stabilize in the  longer interval <math display="inline">[0,T]</math> using 10 subintervals of length 2.

===5.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 <math display="inline">k_1</math> and <math display="inline">k_2</math>). In Figures [[#img-7|7]] and [[#img-8|8]] are shown the respective results.

<div id='img-7'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_LOPEZ_262416069-Ufigure8J23L16.png|276px|]]
|[[Image:Draft_LOPEZ_262416069-Ufigure7J23L16.png|276px|u<sup>∆t</sup> (top) and \| y<sup>∆t</sup>(⋅)\|  (bottom) for several values of k₁ and k₂. The unstable equilibrium θ is given by \{ n₁,n₂,n₃\} =\{ 1,0,0(u)\} . The junctions used to control are 2 and 3.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 7:''' <math>\mathbf{u}^{\Delta t}</math> (top) and <math>\| \mathbf{y}^{\Delta t}(\cdot )\| </math> (bottom) for several values of <math>k_1</math> and <math>k_2</math>. The unstable equilibrium <math>\boldsymbol{\theta }</math> is given by <math>\{ n_{1},n_{2},n_{3}\} =\{ 1,0,0(u)\} </math>. The junctions used to control are 2 and 3.
|}

<div id='img-8'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_LOPEZ_262416069-Ufigure9J23L16.png|384px|ln\| g<sub>q</sub><sup>∆t</sup>\|  for several values of k₁ and k₂. The unstable equilibrium θ is given by \{ n₁,n₂,n₃\} =\{ 1,0,0(u)\} . The junctions used to control are 2 and 3.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 8:''' <math>ln\| \mathbf{g}_q^{\Delta t}\| </math> for several values of <math>k_1</math> and <math>k_2</math>. The unstable equilibrium <math>\boldsymbol{\theta }</math> is given by <math>\{ n_{1},n_{2},n_{3}\} =\{ 1,0,0(u)\} </math>. The junctions used to control are 2 and 3.
|}

===5.2 Controlling via three junctions===

Figures [[#img-9|9]] and [[#img-10|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.

<div id='img-9'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_LOPEZ_262416069-Ufigure11J123L16.png|360px|]]
|[[Image:Draft_LOPEZ_262416069-Ufigure10J123L16.png|360px|u<sup>∆t</sup> (top) and \| y<sup>∆t</sup>(⋅)\|  (bottom) for several values of k₁ and k₂. The unstable equilibrium θ is given by \{ n₁,n₂,n₃\} =\{ 1,0,0(u)\} . The junctions used to control are 1, 2 and 3.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 9:''' <math>\mathbf{u}^{\Delta t}</math> (top) and <math>\| \mathbf{y}^{\Delta t}(\cdot )\| </math> (bottom) for several values of <math>k_1</math> and <math>k_2</math>. The unstable equilibrium <math>\boldsymbol{\theta }</math> is given by <math>\{ n_{1},n_{2},n_{3}\} =\{ 1,0,0(u)\} </math>. The junctions used to control are 1, 2 and 3.
|}

<div id='img-10'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_LOPEZ_262416069-Ufigure12J123L16.png|384px|ln\| g<sub>q</sub><sup>∆t</sup>\|  for several values of k₁ and k₂. The unstable equilibrium θ is given by \{ n₁,n₂,n₃\} =\{ 1,0,0(u)\} . The junctions used to control are 1, 2 and 3.]]
|- style="text-align: center; font-size: 75%;"
| colspan="1" | '''Figure 10:''' <math>ln\| \mathbf{g}_q^{\Delta t}\| </math> for several values of <math>k_1</math> and <math>k_2</math>. The unstable equilibrium <math>\boldsymbol{\theta }</math> is given by <math>\{ n_{1},n_{2},n_{3}\} =\{ 1,0,0(u)\} </math>. The junctions used to control are 1, 2 and 3.
|}

For the particular values <math display="inline">k_1=1E+3</math> and <math display="inline">k_2=1E+3</math> of the penalty parameters  we show in Figure [[#img-11|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 [[#img-9|9]])

<div id='img-11'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_LOPEZ_262416069-Ufigure15t12L16.png|240px|]]
|[[Image:Draft_LOPEZ_262416069-Ufigure17t12L16.png|240px|Optimal controls for the linear system (top), and Euclidean norm of the controlled solution ϕ<sup>∆t</sup> of the nonlinear system  (bottom).]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 11:''' Optimal controls for the linear system (top), and Euclidean norm of the controlled solution <math>\boldsymbol{\phi }^{\Delta t}</math> of the nonlinear system  (bottom).
|}

In Figure [[#img-12|12]] we show the solution for the linear and  nonlinear model, using controls in Figure [[#img-11|11]] (we continue from <math display="inline">t=2</math> to <math display="inline">t=20</math> with no control).

<div id='img-12'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_LOPEZ_262416069-Ufigure13J123L16.png|360px|]]
|[[Image:Draft_LOPEZ_262416069-Ufigure14J123L16.png|360px|Extended controlled solution y<sub>i</sub><sup>∆t</sup> of the linear perturbation model (top), and extended controlled solution ϕ<sub>i</sub><sup>∆t</sup> of the nonlinear system (bottom). After t=2 the controls are extended as cero.]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 12:''' Extended controlled solution <math>y_i^{\Delta t}</math> of the linear perturbation model (top), and extended controlled solution <math>\phi _i^{\Delta t}</math> of the nonlinear system (bottom). After <math>t=2</math> the controls are extended as cero.
|}

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

To control during the whole time interval <math display="inline">(0,20)</math> we have divided the time interval into subintervals of smaller length <math display="inline">\Delta T = T/Q=2</math>, and we denote <math display="inline">q\Delta T</math> by <math display="inline">T_q</math> for <math display="inline">q = 1, . . .,Q</math>; we proceed then as follows:

*   For <math display="inline">q = 0</math>, we denote by <math display="inline">\mathbf{y}_{0}</math> the difference <math display="inline">\boldsymbol{\phi }_0-\boldsymbol{\theta }</math>, and we solve the  associated linear control problem ([[#eq-9|9]])-([[#eq-11|11]]) in <math display="inline">[0,T_1]</math>; let us denote by <math display="inline">\mathbf{u}_{1}</math> the corresponding control. This control is injected in ([[#eq-5|5]]) with initial condition ([[#eq-8|8]]) and <math display="inline">\delta \boldsymbol{\theta }=\mathbf{y}_{0}</math>, and we denote by <math display="inline">\mathbf{y}_{1}</math> the difference <math display="inline">\boldsymbol{\phi }(T_1)-\boldsymbol{\theta }</math>.
*  For <math display="inline">q > 0</math>, we denote by <math display="inline">\mathbf{y}_{q}</math> the difference <math display="inline">\boldsymbol{\phi }(T_q)-\boldsymbol{\theta }</math>; we solve the  associated linear control problem ([[#eq-9|9]])-([[#eq-11|11]]) in <math display="inline">[T_q, T_{q+1}]</math>, with <math display="inline">\mathbf{y}_{0}</math> replaced by <math display="inline">\mathbf{y}_{q}</math>, and we denote by <math display="inline">\mathbf{u}_{q+1}</math> the corresponding optimal control. The control <math display="inline">\mathbf{u}_{q+1}</math> is injected in ([[#eq-5|5]]) with initial condition ([[#eq-8|8]]) and <math display="inline">\delta \boldsymbol{\theta }=\mathbf{y}_{q}</math>,  and we denote by <math display="inline">\mathbf{y}_{q+1}</math> the difference <math display="inline">\boldsymbol{\phi }(T_{q+1})-\boldsymbol{\theta }</math>.
*  We do <math display="inline">q = q + 1</math> and we repeat the process.

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

<div id='img-13'></div>
{| class="floating_imageSCP" style="text-align: center; border: 1px solid #BBB; margin: 1em auto; width: 100%;max-width: 100%;"
|-
|[[Image:Draft_LOPEZ_262416069-Ufigure18t28L16.png|300px|]]
|[[Image:Draft_LOPEZ_262416069-Ufigure19t28L16.png|300px|]]
|-
| colspan="2"|[[Image:Draft_LOPEZ_262416069-Ufigure20t28L16.png|300px|Controls calculated each two seconds (top); Euclidean norm of the controlled solution y<sup>∆t</sup> using the controls in the top (middle); Euclidean norm of the controlled solution ϕ<sup>∆t</sup> using the controls in the top (bottom).]]
|- style="text-align: center; font-size: 75%;"
| colspan="2" | '''Figure 13:''' Controls calculated each two seconds (top); Euclidean norm of the controlled solution <math>\mathbf{y}^{\Delta t}</math> using the controls in the top (middle); Euclidean norm of the controlled solution <math>\boldsymbol{\phi }^{\Delta t}</math> using the controls in the top (bottom).
|}

==6 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 (<math display="inline">t\in [0, 30]</math>) 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===

<div id="cite-1"></div>
'''[[#citeF-1|[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.

<div id="cite-2"></div>
'''[[#citeF-2|[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.

<div id="cite-3"></div>
'''[[#citeF-3|[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.

<div id="cite-4"></div>
'''[[#citeF-4|[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.

<div id="cite-5"></div>
'''[[#citeF-5|[5]]]'''  N. Nair and Y. Braiman. (2018) A Ternary Memory Cell Using Small Josephson Junction Arrays, Superconductor Science and Technology (31).

<div id="cite-6"></div>
'''[[#citeF-6|[6]]]'''  F.&nbsp;Bedard and N.&nbsp;Welker and G.&nbsp;R. Cotter and M.&nbsp;A. Escavage and J.&nbsp;T. Pinkston. (2005) Superconducting Technology Assessment, National Security Agency, Office of Corporate Assessment, url = https://www.nitrd.gov/pubs/nsa/sta.pdf.

<div id="cite-7"></div>
'''[[#citeF-7|[7]]]'''  F.&nbsp;A. Holmes and L.&nbsp;Ripple and M.&nbsp;A. Manheimer. (2013) Energy-Efficient Superconducting Computing-Power Budgets and Requirements, IEEE Transactions on Applied Superconductivity , 23 (3).

<div id="cite-8"></div>
'''[[#citeF-8|[8]]]'''  IARPA. (2013) ''Broad Agency Announcement'': IARPA-BAA-13-05, Cryogenic Computing Complexity (C3) Program, IARPA.

<div id="cite-9"></div>
'''[[#citeF-9|[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.

<div id="cite-10"></div>
'''[[#citeF-10|[10]]]'''  R.&nbsp;Glowinski and J.L.&nbsp;Lions and J.&nbsp;W. He. (2008) Exact and Approximate Controllability for Distributed Parameter Systems, Cambridge University Press, Cambridge, UK.

<div id="cite-11"></div>
'''[[#citeF-11|[11]]]'''  J.&nbsp;L. Lions. (1971) Optimal Control of Systems Governed by Partial Differential Equations, Springer Verlag, New York.