# Overcoming Distribution Shifts in Plug-and-Play Methods with Test-Time
Training

Edward P. Chandler, Shirin Shoushtari, Jiaming Liu,
M. Salman Asif, and Ulugbek S. Kamilov
Washington University in St. Louis, MO 63130, USA
University of California Riverside, CA 92521, USA
{e.p.chandler, s.shirin, jiaming.liu<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Plug-and-Play Priors (PnP) is a well-known class of methods for solving
inverse problems in computational imaging. PnP methods combine physical
forward models with learned prior models specified as image denoisers. A
common issue with the learned models is that of a performance drop when there
is a distribution shift between the training and testing data. Test-time
training (TTT) was recently proposed as a general strategy for improving the
performance of learned models when training and testing data come from
different distributions. In this paper, we propose PnP-TTT as a new method for
overcoming distribution shifts in PnP. PnP-TTT uses deep equilibrium learning
(DEQ) for optimizing a self-supervised loss at the fixed points of PnP
iterations. PnP-TTT can be directly applied on a single test sample to improve
the generalization of PnP. We show through simulations that given a sufficient
number of measurements, PnP-TTT enables the use of image priors trained on
natural images for image reconstruction in magnetic resonance imaging (MRI).

##  1 Introduction

Many computational imaging problems can be formulated as _inverse problems_ ,
where the goal is to recover an unknown image from a set of noisy
measurements. It is common to solve inverse problems by integrating the
measurement model characterizing the response of the imaging instrument with a
regularizer infusing prior knowledge on the unknown image. There has been
considerable recent interest in using deep learning (DL) for designing data-
driven image priors [1, 2, 3]. DL methods eliminate the need for explicit
prior modeling by learning a mapping from measurements to target images using
convolutional neural networks (CNN).

Model-based DL (MBDL) is an extension to traditional DL that integrates the
image prior defined through a CNN with the knowledge of the measurement
models. For example, plug-and-play priors (PnP) is a well-known MBDL approach
that uses pre-trained image denoiser as priors [4, 5, 3]. Other MBDL widely-
used MBDL approaches include deep unfolding (DU) and deep equilibrium (DEQ)
learning, both of which rely on the integration of the measurement model
during the training of the image prior [6, 7, 8, 9, 10]. While both DU and DEQ
interpret iterations of image reconstruction as neural network layers, the
memory complexity of DEQ is independent of the number of unfolded iterations.

Much of the existing research on MBDL has focused on the scenarios where the
statistical distribution of the training data matches that of the testing
data. While this strategy has led to significant theoretical and algorithmic
innovations, it does not address the issue of the performance gap due to data
distribution shifts. For example, image priors trained with a specific
distribution in PnP, performs poorly on samples from different distributions
[11]. Thus, distribution shifts limit the applicability of priors pre-trained
for one class to another one.

Domain adaptation refers to a class of DL techniques for improving the
performance of a learned model on a target task containing insufficient
annotated data by using the knowledge learned by the model from another
related task with adequate labeled data [12, 13]. _Test-time training (TTT)_
was recently proposed as a domain adaptation strategy based on self-supervised
optimization of the learned model utilizing only test-time measurements [14].
The TTT strategy was also recently used in the context of imaging inverse
problems to address domain shifts in end-to-end image reconstruction with DL
for accelerated magnetic resonance imaging (MRI) [15].

In this paper, we investigate TTT in the context of PnP methods. We propose
_PnP-TTT_ as a method for overcoming the performance gap in PnP due to data
distribution shifts. PnP-TTT uses DEQ to update the weights of the CNN prior
in PnP at test-time. The DEQ learning in PnP-TTT is used to minimize a self-
supervised loss at the fixed points of PnP iterations for one test sample. We
also present numerical results showing that DEQ training in PnP-TTT can
significantly boost the performance of the shifted priors. We evaluate the
proposed method on image reconstruction for compressed sensing MRI (CS-MRI),
where we recover MRI images from subsampled Fourier measurements. Our results
show that given enough measurements, PnP-TTT can close the gap due to
distribution shift between test and training data. It is worth mentioning that
our method can also be applied to other tasks and different variants of PnP,
highlighting its broader applicability for inverse problems in computational
imaging.

Algorithm 1 Test-Time Training for Plug-and-Play Methods

0: forward model $\mathbf{A}$, PnP initialization $\boldsymbol{x}_{0}$,
measurement $\boldsymbol{y}$, denoiser ${\mathsf{D}}_{\theta}$, and
numIter$\geq 0$

$\boldsymbol{x}^{*}_{0}=\texttt{PnP}(\boldsymbol{x}_{0},\boldsymbol{y};{\mathsf{D}}_{\theta})$

for $i=1$ to numIter do

$l=\texttt{LOSS}(\mathbf{A}\boldsymbol{x}^{*}_{i-1},\boldsymbol{y})$

$\texttt{DEQ\\_GRAD}(l,\theta)$ {Update parameters $\theta$ using DEQ}

$\boldsymbol{x}^{*}_{i}=\texttt{PnP}(\boldsymbol{x}_{0},\boldsymbol{y};{\mathsf{D}}_{\theta})$

end for

return $\boldsymbol{x}^{*}_{i}$

##  2 Background

###  2.1 Inverse Problems

We consider the problem of recovering an image ${\bm{x}}\in\mathbb{C}^{n}$
from its noisy measurement ${\bm{y}}={\bm{A}}{\bm{x}}+{\bm{e}}$, where
${\bm{A}}\in\mathbb{C}^{m\times n}$ is the _measurement operator_ and
${\bm{e}}\in\mathbb{C}^{m}$ is additive white Gaussian noise (AWGN). We can
formulate the problem as a regularized optimization problem

${\bm{x}^{\ast}}=\mathop{\mathsf{arg\,min}}_{{\bm{x}}\in\mathbb{R}^{n}}f({\bm{x}})\quad\text{with}\quad
f({\bm{x}})=g({\bm{x}})+h({\bm{x}}),$ (1)

where $g$ is the _data-fidelity_ term used to ensure the consistency of the
solution with the measurement and $h$ is the _regularization_ term that
infuses prior knowledge. For example, the least-squares loss is a widely-used
data-fidelity term
$g({\bm{x}})=\frac{1}{2}\|{\bm{y}}-{\bm{A}}{\bm{x}}\|_{2}^{2}$ and total
variation (TV) is commonly used as the regularizer [16].

###  2.2 Plug-and-Play Priors

PnP framework includes a family of methods that incorporate the measurement
model with CNN denoisers to solve inverse problems [3]. PnP methods can be
interpreted as a fixed-point iteration of some high-dimensional operator where
the CNN takes the role of the prior. For example, the _proximal gradient
method (PGM)_ variant of PnP can be expressed as

$\displaystyle{\bm{x}}^{k}={\mathsf{T}}_{{\bm{\theta}}}({\bm{x}}^{k-1})\quad\text{with}\quad{\mathsf{T}}_{{\bm{\theta}}}\,\coloneqq\,{\mathsf{D}}_{{\bm{\theta}}}({\mathsf{I}}-\gamma\nabla
g),$ (2)

where ${\mathsf{D}}_{\bm{\theta}}$ is the denoiser, $g$ is the data-fidelity
term, $\nabla g$ is the gradient of $g$, ${\mathsf{I}}$ is the identity
mapping, and $\gamma>0$ is the step-size. The PnP method in (2) is commonly
refered to as PnP-PGM.

Figure 1: Evaluation of PnP-TTT for different sampling ratios in accelerated MRI. The left-most chart displays the best PSNR performance achieved by PnP-TTT at different sampling ratios. The remaining five charts show PSNR over all examples at each TTT iteration. Note that the best performance is above the lower baseline for all the sampling ratios; however, TTT eventually overfits to the test-time measurement, reducing performance. Additionally, note that at larger sampling ratios, the performance of PnP-TTT prior can surpass that of the matched prior due to the DEQ training. Table 1: PSNR (dB) values for accelerated MRI with matched, mismatched, and PnP-TTT priors. Radial CS Ratio | $10\%$ | $20\%$ | $30\%$ | $40\%$ | $50\%$
---|---|---|---|---|---
| PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM | PSNR | SSIM
Natural prior | $30.81$ | $0.9351$ | $33.87$ | $0.9648$ | $35.29$ | $0.9721$ | $36.37$ | $0.9759$ | $37.25$ | $0.9782$
MRI Prior | $31.67$ | $0.9468$ | $34.9$ | $0.9707$ | $36.51$ | $0.9773$ | $37.67$ | $0.9807$ | $38.57$ | $0.9828$
PnP-TTT (Ours) | $30.97$ | $0.938$ | $34.97$ | $0.9718$ | $37.03$ | $0.9796$ | $38.58$ | $0.984$ | $39.96$ | $0.9873$
PnP-TTT $-$ Natural | $0.16$ | $0.0029$ | $1.1$ | $0.007$ | $1.74$ | $0.0075$ | $2.21$ | $0.0081$ | $2.71$ | $0.0091$

###  2.3 Deep Equilibrium Models

DEQ is a recent approach for training MBDL architectures in a memory-efficient
way [9]. DEQ uses implicit differentiation for training possibly infinite-
depth networks by backpropagating through the fixed points of an operator. For
the operator defined in eq. (2), the output of the DEQ is implicitly expressed
as

${\overline{\bm{x}}}={\mathsf{T}}_{\bm{\theta}}({\overline{\bm{x}}}),$ (3)

where ${\mathsf{T}}_{\bm{\theta}}$ is the operator parameterized by
${\bm{\theta}}$, and ${\overline{\bm{x}}}$ is the fixed point acquired using
fixed point iterations in the forward pass of DEQ. The connection of DEQ and
PnP has inspired end-to-end training of CNN denoisers as model dependant
priors in many imaging problems such as MRI [9] and computed tomography (CT)
[10].

The prior ${\mathsf{D}}_{\bm{\theta}}$ in DEQ is trained by minimizing the
loss between the fixed points from eq. (3) and the ground truth image
${\bm{x}^{\ast}}$

$\ell({\bm{\theta}})=\frac{1}{2}\|{\mathsf{T}}_{\theta}({\overline{\bm{x}}})-{\bm{x}^{\ast}}\|_{2}^{2}.$
(4)

Implicit differentiation of the fixed points yields the gradient of the loss
with respect to ${\bm{\theta}}$ in the backward pass as the following

$\nabla\ell({\bm{\theta}})=(\nabla_{\bm{\theta}}{\mathsf{T}}_{\bm{\theta}}({\overline{\bm{x}}}))^{\mathsf{T}}\left({\mathsf{I}}-\nabla_{\bm{x}}{\mathsf{T}}_{{\bm{\theta}}}({\overline{\bm{x}}})\right)^{-{\mathsf{T}}}({\overline{\bm{x}}}-{\bm{x}^{\ast}}),$
(5)

where ${\mathsf{I}}$ is the identity mapping and $\ell$ is the loss in eq.
(5).

###  2.4 Test-Time Training

Current PnP methods are built on the premise that the prior represents the
same distribution as that of the desired solution. However, it is common to
observe distribution shifts between training and testing data. In some
scenarios, there are insufficient training samples to train a DL network as
the prior, and as a result, alternative priors trained on a shifted
distribution are used. The distribution shift results in suboptimal
reconstruction performance. TTT has been proposed to reduce the performance
gap due to distribution shift in various tasks [14]. The key idea of TTT is to
update the shifted model’s weight at test-time by minimizing a self-supervised
loss

${\bm{\theta}}^{\ast}=\mathop{\mathsf{arg\,min}}_{{\bm{\theta}}}\ell_{\textsf{\tiny
sup}}\left({\mathsf{D}}_{{\bm{\theta}}}({\bm{y}}),{\bm{y}}\right),$ (6)

where ${\mathsf{D}}_{\bm{\theta}}$ is the neural network and ${\bm{y}}$ is a
test sample. Depending on the selection of $\ell_{\textsf{\tiny sup}}$, TTT
has shown improved performance in many imaging tasks. For example, it can be
used to improve the MRI reconstruction using DL models trained in an end-to-
end matter on shifted distributions [15]. In this scenario, the self-
supervised loss proposed is the following

$\ell_{\textsf{\tiny
sup}}({\bm{\theta}})=\frac{\left\|{\bm{A}}{\mathsf{D}}_{{\bm{\theta}}}\left({\bm{A}}^{\dagger}{\bm{y}}\right)-{\bm{y}}\right\|_{1}}{\|{\bm{y}}\|_{1}},$
(7)

where ${\bm{A}}$ is the measurement model, ${\bm{A}}^{\dagger}$ is the
Hermitian transpose, and ${\bm{y}}$ is the test-time measurement. Note that as
opposed to (4), TTT in (7) does not need ground truth reconstruction to
compute $\ell_{\textsf{\tiny sup}}$. Losses other than the normalized
$\ell_{1}$ norm can work as well [15].

###  2.5 Our contribution

We propose PnP-TTT as a novel approach for enhancing the performance of image
reconstruction for PnP methods under distribution shifts. Our approach
involves domain adaptation for a shifted pre-trained image prior through TTT
using DEQ to close the distribution gap, which only requires a test single
measurement. Our results show that PnP-TTT can improve the performance
significantly given sufficient measurement for shifted priors with minimal
computational cost.

##  3 Method

We now present our method for domain adaption of image priors in PnP. We
consider the PnP-PGM algorithm in eq. (2) and run it until its convergence. In
practice, we find that about $100$ iterations of PnP-PGM are sufficient in our
configuration. We can update the weights of the image prior on a test
measurement by minimizing the following self-supervised loss

$\ell_{\textsf{\tiny
sup}}({\bm{\theta}})=\left\|{\bm{A}}{\mathsf{T}}_{\bm{\theta}}({\overline{\bm{x}}})-{\bm{y}}\right\|_{2}^{2},$
(8)

where ${\overline{\bm{x}}}$ is the fixed-point of PnP-PGM defined in eq. (3)
and ${\mathsf{T}}_{\theta}$ is the operator defined in (2). We use the DEQ to
compute the gradient of $\ell_{\textsf{\tiny sup}}$ at test-time using
implicit differentiation.

We follow a method similar to [17, 18] to train the image priors using the
DnCNN architecture, with batch normalization layers replaced with spectral
normalization to control the Lipschitz constant of the denoisers. DnCNN is
trained as a denoiser for AWGN level $\sigma=5$. During the training stage we
do not use DU or DEQ so that the learned prior model is purely an image
denoiser. We use 400 CBSD to train natural prior on grayscale images of size
$180\times 180$ [19]. We train MRI priors on MRI brain images of size
$256\times 256$ [20].

For test-time training, we initialize PnP-PGM with ${\bm{x}}^{0}=\bm{0}$ and
$100$ iterations in the forward pass of DEQ, using the denoising prior
acquired as described above. We use Nesterov acceleration [21], and set
stepsize $\gamma=1$. We use $100$ iterations and Anderson acceleration in the
backward pass of DEQ [22]. We allow TTT to run for 50 iterations, using SGD to
update the parameters $\theta$ with a step size of $1\times 10^{-5}$. At
inference, using the adapted prior, we again run for $100$ iterations with
Nesterov acceleration in PnP-PGM with step size $\gamma=1$. Note that once all
$50$ TTT iterations are performed for a particular measurement, $\theta$ is
reset to the non-domain-adapted weights. Since the goal of PnP-TTT is to
overcome the performance gap from a distribution shift between train- and
test-time, performing TTT on as many measurements are available at test-time
may be beneficial. Future experiments could examine if there is any
performance improvement when using multiple measurements instead of only one
for PnP-TTT.

The measurement model for a single-coil, accelerated MRI with radial Fourier
sampling can be modeled as ${\bm{A}}={\bm{M}}{\bm{F}}$, where ${\bm{M}}$ is
the diagonal sampling matrix and ${\bm{F}}$ is the Fourier transform. We
investigate five different sampling ratios $(m/n)$ in the experiments. For the
experiments reported here, we consider noiseless scenario; however, we expect
similar performance of PnP-TTT under moderate amoungs of noise.

Figure 2: Visual evaluation of various priors at different CS ratios for CS-
MRI problem with corresponding PSNR (upper left) and SSIM (upper right). Note
that natural priors (top row) performs suboptimally compared to matched MRI
prior (bottom row). Additionally, note the improvement due to the usage of
PnP-TTT framework for all CS ratios (middle row).

##  4 Results

We test our proposed method by reconstructing ten brain MRI images selected
from the test dataset of [20] with mismatched DnCNN prior trained on natural
images. Due to distribution shift, natural priors demonstrate suboptimal
performance for the MRI task. To establish a performance baseline, we compare
the result of the proposed method with those obtained by mismatched natural
prior and matched MRI prior. Specifically, we consider the performance
achieved by natural prior as the lower baseline, and that achieved by the MRI
prior as the upper baseline. Our proposed PnP-TTT seeks to enhance the
performance of a mismatched natural prior so as to approach that of a matched
MRI prior.

Table 1 reports the best results achieved for five CS ratios: $10,\ 20,\ 30,\
40,\ \text{and}\ 50$. It can be seen that PnP-TTT can close the performance
gap for CS ratios of $20$ and more, while for a CS ratio of $10$, it can make
an improvement compared to the lower baseline (mismatched natural prior). The
reconstruction quality is quantified using peak signal-to-noise ratio (PSNR)
in dB and the structural similarity index measure (SSIM).

Figure 1 illustrates two empirical results: (_a_) The empirical performance of
PnP-TTT at testing for different CS ratios (_left figure_), (_b_) The
empirical performance during test-time training (_remaining five figures_).
Note that during test-time training, the prior can overfit to the measurement.
Thus, in practice it is necessary to hold out some measurements to use early
stopping during TTT [15], although we have not included such results in this
paper. The visual results can be found in Figure 2. It can be seen both
empirically and visually that PnP-TTT can shorten the gap due to distribution
shift, close it completely, or go beyond closing it given the CS ratios.

##  5 Conclusion

We present PnP-TTT as a novel framework for closing the performance gap that
arises due to mismatched priors in imaging inverse problems solved by PnP
framework. PnP-TTT achieves this by adapting the mismatched priors during the
testing phase by using DEQ training to update the weights of the mismatched
priors. One of the main advantage of PnP-TTT is that one can use mismatched
priors on a shifted distribution without the need to do additional training on
other samples from the shifted distribution. Instead, the prior can simply be
adapted to the test-time measurements. Our results show that PnP-TTT can
significantly enhance the performance, achieving performance comparable to
that of using a matched prior during inference. Furthermore, this work
demonstrate that priors from different tasks can be used interchangeably in
scenarios with shifted distribution without the loss of performance.

## Acknowledgement

This material is partially based upon work supported by the NSF CAREER award
under grant CCF-2043134 and by the Gordon and Betty Moore Foundation grant
11396.

## References

  * [1] M. T. McCann, K. H. Jin, and M. Unser,  “Convolutional neural networks for inverse problems in imaging: A review,”  IEEE Signal Process. Mag., vol. 34, no. 6, pp. 85–95, 2017.
  * [2] A. Lucas, M. Iliadis, R. Molina, and A. K. Katsaggelos,  “Using deep neural networks for inverse problems in imaging: Beyond analytical methods,”  IEEE Signal Process. Mag., vol. 35, no. 1, pp. 20–36, Jan. 2018.
  * [3] U. S. Kamilov, C. A. Bouman, G. T. Buzzard, and B. Wohlberg,  “Plug-and-play methods for integrating physical and learned models in computational imaging,”  IEEE Signal Process. Mag., vol. 40, no. 1, pp. 85–97, Jan. 2023.
  * [4] S. V. Venkatakrishnan, C. A. Bouman, and B. Wohlberg,  “Plug-and-play priors for model based reconstruction,”  IEEE Glob. Conf. Signal Inf. Process., pp. 945–948, 2013.
  * [5] S. Sreehari, S. V. Venkatakrishnan, B. Wohlberg, G. T. Buzzard, L. F. Drummy, J. P. Simmons, and C. A. Bouman,  “Plug-and-play priors for bright field electron tomography and sparse interpolation,”  IEEE Trans. Comput. Imag., vol. 2, no. 4, pp. 408–423, December 2016.
  * [6] S. A. Hosseini, B. Yaman, S. Moeller, M. Hong, and M. Akcakaya,  “Dense recurrent neural networks for accelerated MRI: History-cognizant unrolling of optimization algorithms,”  IEEE J. Sel. Topics Signal Process., vol. 14, no. 6, pp. 1280–1291, Oct. 2020.
  * [7] G. Ongie, A. Jalal, C. A. Metzler, R. G. Baraniuk, A. G. Dimakis, and R. Willett,  “Deep learning techniques for inverse problems in imaging,”  IEEE J. Sel. Areas Inf. Theory, vol. 1, no. 1, pp. 39–56, May 2020.
  * [8] V. Monga, Y. Li, and Y. C. Eldar,  “Algorithm unrolling: Interpretable, efficient deep learning for signal and image processing,”  IEEE Signal Process. Mag., vol. 38, no. 2, pp. 18–44, Mar. 2021.
  * [9] D. Gilton, G. Ongie, and R. Willett,  “Deep equilibrium architectures for inverse problems in imaging,”  IEEE Trans. Comput. Imag., vol. 7, pp. 1123–1133, 2021.
  * [10] J. Liu, X. Xu, W. Gan, S. Shoushtari, and U. S. Kamilov,  “Online deep equilibrium learning for regularization by denoising,”  in Proc. Advances Neural Inf. Process. Syst., New Orleans, LA, 2022.
  * [11] Shirin Shoushtari, Jiaming Liu, Yuyang Hu, and Ulugbek S Kamilov,  “Deep model-based architectures for inverse problems under mismatched priors,”  IEEE J. Sel. Areas Inf. Theory, 2022.
  * [12] H. Venkateswara, S. Chakraborty, and S. Panchanathan,  “Deep-learning systems for domain adaptation in computer vision: Learning transferable feature representations,”  IEEE Signal Process. Mag., vol. 34, no. 6, pp. 117–129, Nov. 2017.
  * [13] Abolfazl Farahani, Sahar Voghoei, Khaled Rasheed, and Hamid R. Arabnia,  “A brief review of domain adaptation,”  in Advances in Data Science and Inf. Eng., Cham, 2021, pp. 877–894, Springer International Publishing.
  * [14] Yu Sun, Xiaolong Wang, Zhuang Liu, John Miller, Alexei Efros, and Moritz Hardt,  “Test-time training with self-supervision for generalization under distribution shifts,”  in Proc. 37th Int. Conf. on Mach. Learn., 13–18 Jul 2020, vol. 119, pp. 9229–9248.
  * [15] Mohammad Zalbagi Darestani, Jiayu Liu, and Reinhard Heckel,  “Test-time training can close the natural distribution shift performance gap in deep learning based compressed sensing,”  in Proc. 39th Int. Conf. on Mach. Learn., 17–23 Jul 2022, vol. 162, pp. 4754–4776.
  * [16] A. Beck and M. Teboulle,  “Fast gradient-based algorithm for constrained total variation image denoising and deblurring problems,”  IEEE Trans. Image Process., vol. 18, no. 11, pp. 2419–2434, November 2009.
  * [17] Kai Zhang, Wangmeng Zuo, Yunjin Chen, Deyu Meng, and Lei Zhang,  “Beyond a Gaussian denoiser: Residual learning of deep CNN for image denoising,”  IEEE Trans. Image Process., vol. 26, no. 7, pp. 3142–3155, 2017.
  * [18] Takeru Miyato, Toshiki Kataoka, Masanori Koyama, and Yuichi Yoshida,  “Spectral normalization for generative adversarial networks,”  in Int. Conf. on Learn. Represent., 2018.
  * [19] D. Martin, C. Fowlkes, D. Tal, and J. Malik,  “A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics,”  in Proc. 8th Int. Conf. Comput. Vis., July 2001, vol. 2, pp. 416–423.
  * [20] J. Zhang and B. Ghanem,  “Ista-net: Interpretable optimization-inspired deep network for image compressive sensing,”  IEEE Conf. Comput. Vis. Pattern Recognit., pp. 1828–1837, 2018.
  * [21] Y. Nesterov,  Introductory lectures on convex optimization: A basic course, vol. 87,  Springer Science & Business Media, 2003.
  * [22] D. G. Anderson,  “Iterative procedures for nonlinear integral equations,”  J. ACM, vol. 12, no. 4, pp. 547–560, 1965.


# Optimal motion of a scallop: some case studies

Rosario Maggistro1 and Marta Zoppello2 1Department of Management, Università
Ca’ Foscari Venezia, Fondamenta S. Giobbe, 873, 30121 Cannaregio, Venezia,
Italy<EMAIL_ADDRESS>of Computer Science, Università di
Verona, Strada le Grazie, 15, 37134, Verona, Italy<EMAIL_ADDRESS>
###### Abstract

In this paper we focus on a two-link swimmer called scallop which moves
changing dynamics between two fluids regimes. We address and solve explicitly
two optimal control problems, the minimum time one and the minimum quadratic
cost needed to move the swimmer between two fixed positions using a periodic
control. Considering only one switching in the dynamics and exploiting the
structure of the equation of motion we are able to split the problem into
simpler ones. We solve explicitly each sub-problem obtaining a discontinuous
global solution. Then we approximate it through a suitable sequence of
continuous functions. Finally, we show numerical simulations suggesting that
to switch less times is the best strategy for both costs.

###### Index Terms:

Scallop swimmer, switching dynamics, optimal control, Pontryagin’s Maximum
Principle.

##  I INTRODUCTION

The study of locomotion strategies in fluids and its link with the
construction of artificial devices that can self-propel in fluids have
generated, in the recent years, considerable interest by different research
communities. Theories of swimming generally utilize either low Reynolds number
approximation [1, 2] and Resistive Force Theory (RTF) [10], or the assumption
of inviscid ideal fluid dynamics (high Reynolds number) [6, 20]. These two
different regimes are also distinct in terms of the mechanism of locomotion
[8]–[9]. In this paper we focus on swimmers immersed in these two kind of
fluids which produce a linear dynamics. In particular we consider the system
describing the motion of a scallop for which several strategies were proposed
in order to overcome the so called Scallop theorem. According to it a swimmer
that moves opening and closing periodically its valves, cannot achieve any net
motion both in a viscous and ideal irrotational fluid, because of the time
reversibility of the equations [14, 2]. Some authors in order to solve the no
net motion problem, added, for example, a degree of freedom introducing the
Purcell swimmer [14], or the three sphere swimmer [13]. Others, instead,
supposed the scallop immersed in a non Newtonian fluid, in which the viscosity
is not constant, ending up with a non reversible dynamics [7]–[15]. Here,
inspired by an idea of A. Bressan, we use the approach proposed in [16] where
a change in the fluid’s regime between the opening and closing of the valves
is supposed. Moreover an hysteresis operator is introduced through a
thermostat, see Fig. 2 (see [17] for mathematical models for hysteresis), to
model a delay in this change of fluid’s regime. This operator has been used in
other contexts to model the switching dynamics phenomenon [4]. In our
framework the use of this operator on the one hand permits both to overcome
the Scallop theorem and make the system completely controllable. On the other
hand it requires the use of a continuous periodic control (angular velocity)
whose explicit expression is difficult to derive. In this paper we address the
optimal time and optimal quadratic cost problem. This has been done for other
kind of swimmers, for example the Purcell three link one [11, 12].

Our main result is to provide an explicit periodic solution of two optimal
control problems, i.e. the minimum time and the linear quadratic optimal
control one, so that the scallop can move between two fixed positions. The
existence of such a solution is guaranteed thanks to the controllability of
the scallop system with switching dynamics and only relaxing the hypothesis on
the regularity of the control. Moreover, due to the structure of the equations
of motion, we reduce to the case of only one switching and overcome the
difficulty of managing the changing dynamics splitting the optimal control
problem into simpler sub-problems which, that can be solved in the relaxed
framework. Thus using the Pontryagin Maximum Principle (PMP) [3, 19], we solve
explicitly the relaxed sub-problems.Then we recover a continuous periodic
solution building a suitable sequence of continuous controls which approximate
in $L^{1}$ the (possible) discontinuous optimal one. In the quadratic case
where only the square of the $L^{2}$-norm of the control is take into account
as a cost function, we prove that the control that minimize that cost is the
same one that solve the minimum time problem.
Moreover, through numerical simulations we give some hints on what happens in
the case of $n$ switchings in the dynamics. More precisely we find an upper
bound for the cost with $n$ switchings and show numerical simulations that
suggest that performing less switchings is the best strategy.

##  II Setting of the problem

###  II-A The model

We recall the scallop swimmer introduced in [16] and present its equation of
motion as a system of ordinary differential equations. The swimmer is modelled
as an articulated rigid body immersed in a fluid that changes its
configuration. It is composed by two rigid valves of elliptical shape of mass
$m$, major semiaxes $a$ and minor $b$ with $b<<a$. The valves are joined in
order that they can be opened and closed. Moreover this body is constrained to
move along the $\vec{e}_{x}$-axis and is symmetric with respect to it. Finally
we will neglect the interaction between the two valves. In order to determine
completely its state we need the position of the juncture point $x$ and the
orientation of each valve $\theta$ with respect to the $\vec{e}_{x}$-axis.

Figure 1: The scallop configuration

The temporal evolution of these coordinates is obtained solving the Newton’s
equations coupled with the Navier-Stokes equations relative to the surrounding
fluid. We will face this problem considering the body as immersed in two kinds
of different fluids: one viscous at low Reynolds number in which we neglect
the effects of inertia, and another one ideal inviscid and irrotational, in
which we neglect the viscous forces in the Navier-Stokes equations. In the
viscous case we consider the two valves of the scallop as one dimensional
links and use the local drag approximation of the Resistive Force Theory [10]
to compute the total drag force exerted by the fluid on the swimmer. According
to it the density of the force $\vec{f}_{i}$ acting on the $i$-th segment at
the point of arc-length $s$, is assumed to depend linearly on the velocity of
that point. It is defined by

$\vec{f}_{i}(s):=-\xi\left(\vec{v}_{i}(s)\cdot\vec{e}_{i}\right)\vec{e}_{i}-\eta\left(\vec{v}_{i}(s)\cdot\vec{e}_{i}^{\bot}\right)\vec{e}_{i}^{\bot},$
(1)

where $\xi$ and $\eta$ are respectively the drag coefficients in the
directions parallel, $\vec{e}_{i}$, and perpendicular, $\vec{e}_{i}^{\bot}$,
to each link. Integrating the density of force and neglecting inertia, (sum of
forces equal zero), we obtain the following dynamics (see [16] for details)

$\dot{x}=V_{1}(\theta)\dot{\theta}=\frac{a\eta\sin(\theta)}{\xi\cos^{2}(\theta)+\eta\sin^{2}(\theta)}\dot{\theta}$
(2)

In the ideal case, instead, it has been proven in [18], [20] that the system
body $+$ fluid is geodetic with Lagrangian given by the sum of the kinetic
energy of the body ($T^{b}$) and the one of the fluid ($T^{f}$).

$T^{tot}=T^{b}+T^{f}$

Moreover since inertial forces dominates over the viscous ones, in order to
derive the kinetic energy of the fluid we will make use of the concept of
added mass [5]. To this end we consider each valve as a thin ellipse.
Following a procedure introduced by Alberto Bressan in [18], in order to end
up with a control system we perform a partial Legendre transformation on the
kinetic energy defining

$p=\frac{\partial T^{tot}}{\partial\dot{x}}$

from which we derive the equation of motion for $x$

$\dot{x}=V_{2}(\theta)\dot{\theta}=\frac{2a\dot{\theta}\sin\theta(m+m_{22})}{2(m+m_{11}\cos^{2}\theta+m_{22}\sin^{2}\theta)}$
(3)

where $m_{ii}$, $i=1\cdots 3$, are the added mass coefficients of the ellipse
[5].
Finally, to set notation, from now on we name $F_{i}$ the primitives of the
dynamics $V_{i}$ $i=1,2$.

###  II-B Thermostatic-like case and controllability

Our idea is now to use the angular velocity of opening and closing of the
valves as a control, and according to [16], it is possible to overcome the
famous Scallop theorem, swithching between the two fluid regimes. More
precisely if the valves are opening ($\dot{\theta}>0$) we suppose that the
scallop is immersed in an ideal fluid; instead when the valves are closing
($\dot{\theta}<0$) we assume the scallop immersed in a viscous fluid. This
idea is inspired by [7] where the fluid has a pseudoelastic nature that helps
the valve opening but resist the valve closing. Furthermore in [16] it has
been proven that the system is completely controllable using a periodic
control $\dot{\theta}$ only introducing a delayed thermostatic rule. We
suppose that the dynamics $V$ depends on the angle
$\theta\in]0,\frac{\pi}{2}[$, and also depends on a discrete variable
$w\in\left\\{1,2\right\\}$, which characterizes the fluid regime ($2$ ideal or
$1$ viscous) whose evolution is governed by a delayed thermostatic rule,
$h_{\varepsilon}[\cdot]$, subject to the evolution of the control $u$. In Fig.
2 the behavior of such a rule is explained, correspondingly to the choice of a
fixed threshold parameter $\varepsilon>0$ (see [17] for a complete
explanation). This rule models a situation in which the change between the two
fluid regimes is not abrupt from opening to closure of the valves, but the
system remains in the fluid regime of origin as long as the angular velocity
has reached the threshold $\varepsilon$, then switches.

Figure 2:  The thermostatic approximation

The controlled evolution is then given by

$\begin{cases}\dot{x}(t)=V_{w(t)}(\theta(t))u(t),\\\ \dot{\theta}(t)=u(t)\\\
w(t)=h_{\varepsilon}[u](t)\\\ x(0)=x_{0},\ \ \theta(0)=\theta_{0}\\\
u(0)=u_{0},\ \ w(0)=w_{0}\end{cases}$ (4)

##  III Minimum time optimal control problem

In this section we study a minimum time problem for the scallop, we provide
the conditions for which there exists an optimal solution and compute it
explicitly.
We assume that the swimmer starts at the initial configuration
$(x(0),\theta(0),u(0))$ and we want to find a swimming strategy that minimizes
the time to reach the final configuration $(x(t_{f}),\theta(t_{f}),u(t_{f}))$.
Note that in [16] we consider three possible cases: i)
$u_{0}\in[-\varepsilon,\varepsilon]\,,w_{0}=2$; ii)
$u_{0}<-\varepsilon\,,w_{0}=~{}1$; iii) $u_{0}>\varepsilon\,,w_{0}=2$. Here,
we will consider only the case i) observing that the other two cases can be
treated in analogous way. Then we have

$\begin{cases}\displaystyle\inf_{u}\int_{0}^{t_{f}}ds\\\
\dot{x}(t)=V_{w}(\theta(t))u(t),\qquad\forall t\in[0,t_{f}],\\\
\dot{\theta}(t)=u(t),\qquad\forall t\in[0,t_{f}],\\\
w(t)=h_{\varepsilon}[u](t)\\\ u\in C^{0}[0,t_{f}],\qquad
u(0)=u(t_{f})=u_{0},\\\ u(t)\in U:=[-\varepsilon,\varepsilon],\qquad\forall
t\in[0,t_{f}],\\\ x(0)=x_{0},\qquad x(t_{f})=x_{f},\\\
\theta(0)=\theta(t_{f})=\theta_{0}\par\end{cases}$ (5)

It has been already proved in [16] that the system (5) is controllable by
using a control $u$ in $C^{0}[0,t_{f}]$. However, since $u\in C^{0}[0,t_{f}]$,
the theorem of Filippov-Cesary [19] can not be applied. Indeed it is based on
approximating the optimal solution with a minimizing sequence that should
converge exactly to the optimal solution. In the hypotheses of Filippov-Cesary
theorem the space in which the controls live is $L^{\infty}(\cdot)$, and the
dynamics evaluated in the sequence of minimizing controls converges weakly-
star in $L^{\infty}(\cdot)$ to the dynamics evaluated in the optimal control.
This is due to the fact that the space of admissible controls is complete with
respect to the $L^{1}(\cdot)$ topology. In our case instead we have that the
space of admissible controls is $C^{0}[0,t_{f}]$ which is not complete with
respect to the $L^{1}(0,t_{f})$ topology. Therefore, the theorem can not be
used to prove the existence of the optimal solution of (5). Actually we have
the following

###### Proposition 1

There is no solution of the minimum time optimal control problem (5).

###### Proof:

Suppose for absurd that there exists an optimal solution
$(x^{*},\theta^{*},u^{*})\in C^{0}[0,t_{f}]\times C^{1}[0,t_{f}]\times
C^{0}[0,t_{f}]$ of the problem (5) with $x^{*}=x(u^{*})$ and
$\theta^{*}=\theta(u^{*})$. This means that there exists a minimizing
trajectories sequence $(x_{n},\theta_{n})_{n}\in C^{0}[0,t_{f}]\times
C^{1}[0,t_{f}]$ (with $x_{n}=x(u_{n})$ and $\theta_{n}=\theta(u_{n})$) such
that

$\lim_{n\to\infty}\|(x_{n},\theta_{n})-(x^{*},\theta^{*})\|_{\infty}=0.$ (6)

Condition (6) implies that $u_{n}\to u^{*}$ in $L^{1}(0,t_{f})$ from which
follows that $u^{*}$ not necessarily belongs to $C^{0}[0,t_{f}]$, since
$C^{0}[0,t_{f}]$ is not complete with respect to the $L^{1}(0,t_{f})$
topology. This fact contradicts the hypothesis, hence the thesis holds. ∎

As consequnce of the Proposition 1, to obtain an optimal solution for (5) our
strategy consists in relaxing the hypothesis on the regularity of the control,
to prove the existence of a solution and then try to approximate it with
sequences of continuous functions which converge to the relaxed solution. The
relaxed formulation of (5) is

$\begin{cases}\displaystyle\inf_{u}\int_{0}^{t_{f}}ds\\\
\dot{x}(t)=V_{w}(\theta(t))u(t),\qquad\forall t\in[0,t_{f}],\\\
\dot{\theta}(t)=u(t),\qquad\forall t\in[0,t_{f}],\\\ u\in
L^{\infty}[0,t_{f}],\ u(t)\in U:=[-\varepsilon,\varepsilon],\ \forall
t\in[0,t_{f}],\\\ w(t)=h_{\varepsilon}[u](t)\\\ x(0)=x_{0},\qquad
x(t_{f})=x_{f},\\\ \theta(0)=\theta(t_{f})=\theta_{0}\\\ \end{cases}$ (7)

Note that in (7), since the control $u$ is in $L^{\infty}$, we do not consider
the initial and final condition on the control $u$ because now we are
interesting on its value on positive measure intervals instead of on some
points. Moreover, even if the control is no more continuous the system remains
controllable (as it is shown in the final numerical simulations in [16]) since
$C^{0}[0,t_{f}]\subset L^{\infty}[0,t_{f}]$. Now by applying Filippov-Cesary
Theorem ([19]), there exists a minimal time such that the constraints are
satisfied i.e., the infimum in (7) is a minimum.
Now, the aim is to find an explicit solution of (7). At first we suppose that
in order to arrive in $x_{f}$ starting from $x_{0}$ the dynamics switches from
$2$ to $1$ only one time as in case i) in [16]. Therefore we can split (7) in
two simpler sub-problems. Indeed, fixed $\Delta
x=F_{2}(\theta(t_{1}))-F_{1}(\theta(t_{1}))$ according to [16] there exists
$\theta_{1}$ such that $\theta(t_{1})=\theta_{1}$, with $t_{1}$ the first
switching instant, such that the corresponding $\Delta x$ is exactly
$x_{f}-x_{0}$. Then we have

$(P1)\begin{cases}\displaystyle\min_{u}t_{1}\\\
\dot{\theta}(t)=u(t),\qquad\forall t\in[0,t_{1}],\\\
u(t)\in[-\varepsilon,\varepsilon],\qquad\forall t\in[0,t_{1}],\\\
\theta(0)=\theta_{0}\quad\theta(t_{1})=\theta_{1}\end{cases}$
$(P2)\begin{cases}\displaystyle\min_{u}(t_{f}-t_{1})\\\
\dot{\theta}(t)=u(t),\qquad\forall t\in[0,t_{f}-t_{1}],\\\
u(t)\in[-\varepsilon,\varepsilon],\quad\forall t\in[0,t_{f}-t_{1}],\\\
\theta(0)=\theta_{1}\quad\theta(t_{f}-t_{1})=\theta_{0}.\\\ \end{cases}$

We start studying the problem $(P1)$. We apply the Pontryagin Maximum
Principle (PMP) to the associated Hamiltonian

$H(\theta,p)=1+pu,$ (8)

and we get

$\dot{p}=-\frac{\partial H}{\partial\theta}=0\Longrightarrow p(t)=p_{0}.$ (9)

By the stationarity condition $H(\theta(t_{1}),p(t_{1}))=1+p_{0}u(t_{1})$ with
$\theta(t_{1})=\theta_{1}$ we get
Case 1: if $\theta_{1}<\theta_{0}$ then $u(t_{1})=-\varepsilon$ and
$p_{0}=1/\varepsilon$;
Case 2: if $\theta_{1}>\theta_{0}$ then $u(t_{1})=\varepsilon$ and
$p_{0}=-1/\varepsilon$.
Accordingly we have

$\small u=\begin{cases}\varepsilon&\text{if}\quad p_{0}<0\\\
-\varepsilon&\text{if}\quad p_{0}>0,\end{cases}\,,\
\theta(t)=\begin{cases}\varepsilon t+\theta_{0}&\text{if}\quad p_{0}<0,\\\
-\varepsilon t+\theta_{0}&\text{if}\quad p_{0}>0\end{cases}$ (10)

and

$t_{1}=\begin{cases}(\theta_{0}-\theta_{1})/\varepsilon&\text{if}\quad\theta_{1}<\theta_{0},\\\
(\theta_{1}-\theta_{0})/\varepsilon&\text{if}\quad\theta_{1}>\theta_{0}.\end{cases}$
(11)

Now we consider the problem $(P2)$ and let $t_{2}=t_{f}-t_{1}$. Proceeding as
in (8)-(9) and by the stationarity condition
$H(\theta(t_{2}),p(t_{2}))=1+p_{0}u(t_{2})$ with $\theta(t_{2})=\theta_{0}$
follows that
Case 1: if $\theta_{1}<\theta_{0}$ then $u(t_{2})=\varepsilon$ and
$p_{0}=-1/\varepsilon$;
Case 2: if $\theta_{1}>\theta_{0}$ then $u(t_{2})=-\varepsilon$ and
$p_{0}=1/\varepsilon$.
Hence

$u=\begin{cases}\varepsilon&\text{if}\quad p_{0}<0\\\
-\varepsilon&\text{if}\quad p_{0}>0,\end{cases}\,,\
\theta(t)=\begin{cases}\varepsilon t+\theta_{1}&\text{if}\quad p_{0}<0,\\\
-\varepsilon t+\theta_{1}&\text{if}\quad p_{0}>0\end{cases}$ (12)

and

$t_{2}=\begin{cases}(\theta_{0}-\theta_{1})/\varepsilon&\text{if
}\theta_{1}<\theta_{0},\\\ (\theta_{1}-\theta_{0})/\varepsilon&\text{if
}\theta_{1}>\theta_{0}.\end{cases}$

Therefore we get in both cases

$t_{f}=t_{1}+t_{2}=\frac{2|\theta_{0}-\theta_{1}|}{\varepsilon}$ (13)

We have obtained the bang bang optimal control (10)-(12) supposing that we do
not need a continuous control, indeed the existence of a continuous optimal
control is not ensured. Nevertheless we can use a continuous approximation of
that bang bang control using piecewise linear functions as shown in Fig. 3 and
computed below.

Figure 3:  Continuous approximation of the control in (10)-(12)

In the situation described by Case 1 for the problem $(P1)$ let

$u_{k}(t)=\begin{cases}-kt(\varepsilon+u_{0})+u_{0}&0\leq t\leq 1/k,\\\
-\varepsilon&1/k\leq t\leq\tilde{t}_{1_{k}},\end{cases}$ (14)
$\theta_{k}(t)=\begin{cases}\displaystyle\frac{-kt^{2}}{2}(\varepsilon+u_{0})+u_{0}t+\theta_{0}&0\leq
t\leq 1/k,\\\
\displaystyle-\varepsilon(t-1/k)+\frac{u_{0}-\varepsilon}{2k}+\theta_{0}&1/k\leq
t\leq\tilde{t}_{1_{k}}.\end{cases}$ (15)

be the continuous approximation of the control and the corresponding
trajectory in (10) . Note that now we have to consider again the fact that in
(5) $u(0)=u_{0}=u(t_{f})$. Therefore we can compute the new switching time
$\tilde{t}_{1_{k}}$ imposing the following

$\displaystyle\theta_{k}(\tilde{t}_{1_{k}})=\theta_{1}\Longleftrightarrow\tilde{t}_{1_{k}}=\frac{u_{0}+\varepsilon+2k(\theta_{0}-\theta_{1})}{2k\varepsilon}.$

It is clearly different from the optimal switching time $t_{1}$ (11) but
converges to it as $k\to\infty$.
In Case 2 for $(P1)$ and Case 1 and 2 for $(P2)$, similar approximations to
ones in (14)–(15) hold.
Therefore we compute the time $\tilde{t}_{2_{k}}$ imposing the final condition
for $\theta_{k}$

$\displaystyle\theta_{k}(\tilde{t}_{2_{k}})=\theta_{0}\Longleftrightarrow\tilde{t}_{2_{k}}=\frac{\varepsilon+k(\theta_{0}-\theta_{1})}{k\varepsilon},$

Finally, putting together the values of $\tilde{t}_{1_{k}}$ and
$\tilde{t}_{2_{k}}$ we get the final time

$\displaystyle
t_{f_{k}}=\tilde{t}_{1_{k}}+\tilde{t}_{2_{k}}=\frac{u_{0}+3\varepsilon+4k(\theta_{0}-\theta_{1})}{2k\varepsilon}$
(16)

In both cases above, the approximation through piecewise linear function
implies that:

$\displaystyle\lim_{k\to+\infty}u_{k}=u\ \text{in}\ L^{1}([0,t_{f}])$
$\displaystyle\lim_{k\to+\infty}t_{f_{k}}=t_{f}.$

Moreover we have the following proposition.

###### Proposition 2

The trajectory $\theta_{k}$ converges uniformly to $\theta$, i.e.

$\lim_{k\to+\infty}\theta_{k}=\theta\ \text{in}\ L^{\infty}([0,t_{f}])$

###### Proof:

We have just shown that $\lim_{k\to+\infty}u_{k}=~{}u$ in $L^{1}([0,t_{f}])$
and $u_{k}$ are equibounded in $L^{\infty}([0,t_{f}])$. Thus we have the
following estimate

$||\theta_{k}-\theta||_{\infty}\leq\int_{0}^{t_{f_{k}}}|u_{k}(s)-u(s)|\,ds$

where we assume to continuously extend the control $u$ in the interval
$[t_{f},t_{f_{k}}]$. Therefore we conclude for the $L^{1}$ convergence of
$u_{k}$ to $u$. ∎

##  IV Linear Quadratic optimal control problem

In this section we minimize a linear quadratic cost instead of the minimum
time as in the previous section. In particular, we consider the following
problem

$\begin{cases}\displaystyle\inf_{u\in[-\varepsilon,\varepsilon]}\int_{0}^{t_{f}}Au^{2}(s)+B\theta(s)^{2}\,ds,\quad
A,B>0\\\ \dot{x}(t)=V_{w}(\theta(t))u(t)\qquad\forall t\in[0,t_{f}],\\\
\dot{\theta}(t)=u(t),\qquad\forall t\in[0,t_{f}],\\\
w(t)=h_{\varepsilon}[u](t),\\\ u\in C^{0}[0,t_{f}],\qquad
u(0)=u(t_{f})=u_{0},\\\ u(t)\in U:=[-\varepsilon,\varepsilon],\qquad\forall
t\in[0,t_{f}],\\\ x(0)=x_{0},\qquad x(t_{f})=x_{f},\\\
\theta(0)=\theta_{0},\qquad\theta(t_{f})=\theta_{0}.\end{cases}$ (17)

that is a generalization of the linear quadratic cost associated to the
linearized dynamics, thus in some sense can be considered as a generalized
energy. In order to compute the optimal solution we proceed as before: first
we relax the problem considering controls $u\in L^{2}(0,t_{f})$. Then
exploiting the properties of $\Delta x$ as function of $\theta(t_{1})$, since
there is a unique switching dynamics and the running cost is nonnegative, we
divide the integral in (17) in two parts. The first, in the time interval
$[0,t_{1}]$, the second one in $[t_{1},t_{f}]$. Hence consider the first
problem

$\begin{cases}\displaystyle\min_{u\in[-\varepsilon,\varepsilon]}J[u]=\min_{u\in[-\varepsilon,\varepsilon]}\int_{0}^{t_{1}}Au^{2}(s)+B\theta(s)^{2}\,ds\\\
\dot{\theta}(t)=u(t),\qquad\forall t\in[0,t_{1}],\\\ u\in
U=[-\varepsilon,\varepsilon],\quad u\in L^{2}[0,t_{1}],\\\
\theta(0)=\theta_{0},\qquad\theta(t_{1})=\theta_{1}.\end{cases}$ (18)

where we don’t consider the initial and final condition on the control $u$
because now it is $L^{2}(0,t_{1})$. Let us now focus on the case
$\theta_{1}>\theta_{0}$, (the reverse inequality can be treated similarly) and
apply the PMP considering the associated Hamiltonian

$H(\theta,u,p)=Au^{2}+B\theta^{2}+pu$

Then, through the necessary condition for optimality, we get $\forall
t\in[0,t_{1}]$

$\begin{cases}p(t)=-2\theta_{0}\sqrt{AB}\exp\left(\frac{\sqrt{B}}{\sqrt{A}}t\right),\\\
u(t)=\frac{\sqrt{B}}{\sqrt{A}}\theta_{0}\exp\left(\frac{\sqrt{B}}{\sqrt{A}}t\right),\end{cases}\quad\text{if}\quad\theta_{1}>\theta_{0}>0$

Now we have to consider the constraint $|u(t)|\leq\varepsilon$. Since the
control $u$ is positive we get the following condition on $t$

$t\leq\sqrt{\frac{A}{B}}\log\Bigl{(}\sqrt{\frac{A}{B}}\frac{\varepsilon}{\theta_{0}}\Bigr{)}.$
(19)

In (19) if $\varepsilon$ is small compared to $A,B,\theta_{0}$ one gets that
$\log\Bigl{(}\sqrt{\frac{A}{B}}\frac{\varepsilon}{\theta_{0}}\Bigr{)}<0$,
therefore for $t>0$ we have

$u(t)=\varepsilon.$ (20)

Thus considering the condition $\theta(t_{1})=\theta_{1}$ we have

$t_{1}=(\theta_{1}-\theta_{0})/{\varepsilon}$

The second problem is completely analogous to the first one (18) but in the
time interval $[t_{1},t_{f}]$ and with reversed initial and final conditions
on $\theta$. Therefore the solution can be computed in a similar way

$\begin{cases}p(t)=2\theta_{1}\sqrt{AB}\exp\left(-\frac{\sqrt{B}}{\sqrt{A}}t\right),\\\
u(t)=-\frac{\sqrt{B}}{\sqrt{A}}\theta_{1}\exp\left(-\frac{\sqrt{B}}{\sqrt{A}}t\right),\end{cases}\quad\text{if}\quad\theta_{1}>\theta_{0}>0$

Imposing again the constraint $|u(t)|\leq\varepsilon$ and following the
considerations made for the first problem we obtain a condition on $t$ which
leads to the following control

$u(t)=\begin{cases}-\varepsilon\quad
t<-\sqrt{\frac{A}{B}}\log\Bigl{(}\sqrt{\frac{A}{B}}\frac{\varepsilon}{\theta_{1}}\Bigr{)},\\\
-\frac{\sqrt{B}}{\sqrt{A}}\theta_{1}\exp\left(-\frac{\sqrt{B}}{\sqrt{A}}t\right)\quad
t\geq-\sqrt{\frac{A}{B}}\log\Bigl{(}\sqrt{\frac{A}{B}}\frac{\varepsilon}{\theta_{1}}\Bigr{)}\end{cases}$

Consequently we compute the corresponding final time

$t_{2}=t_{f}-t_{1}=\begin{cases}\frac{\sqrt{A}}{\sqrt{B}}\log\left(\frac{\theta_{1}}{\theta_{0}}\right)&\text{if}\quad\theta_{1}>\varepsilon\sqrt{\frac{A}{B}}>\theta_{0},\\\
(\theta_{1}-\theta_{0})/{\varepsilon}&\text{if}\quad\theta_{1}>\theta_{0}>\varepsilon\sqrt{\frac{A}{B}},\end{cases}$

which gives us

$t_{f}=\begin{cases}\frac{\theta_{1}-\theta_{0}}{\varepsilon}+\sqrt{\frac{A}{B}}\log\left(\frac{\theta_{1}}{\theta_{0}}\right)&\text{if}\quad\theta_{1}>\sqrt{\frac{A}{B}}\varepsilon>\theta_{0},\\\
2(\theta_{1}-\theta_{0})/{\varepsilon}&\text{if}\quad\theta_{1}>\theta_{0}>\sqrt{\frac{A}{B}}\varepsilon.\\\
\end{cases}$

Note that the solution of (17) $(\theta(t),u(t))$ for $t\in[0,t_{f}]$ is not
continuous in $t_{1}\in(0,t_{f})$, hence we can consider a continuous
approximation of the control $u$ through piecewise linear functions as done
for the minimum time problem. For that we have to reconsider the initial and
final condition $u(0)=u_{0}=u(t_{f})$. In this way we get $u_{k},\theta_{k}\in
C^{0}[0,t_{f_{k}}]$ and $t_{f_{k}}$ such that

$\lim_{k\to+\infty}u_{k}=u\ \text{in}\
L^{2}([0,t_{f}]),\qquad\lim_{k\to+\infty}t_{f_{k}}=t_{f}.$

Note that the above convergence of $u_{k}$ in $L^{2}$ is also valid in $L^{1}$
due to the boundedness of $[0,t_{f}]$. Moreover the following hold

###### Proposition 3

The trajectory $\theta_{k}$ and the cost $J_{k}=J[u_{k}]$ converge uniformly
to $\theta$ and $J$ respectively, i.e.

$\displaystyle\lim_{k\to+\infty}\theta_{k}=\theta\ \text{in}\
L^{\infty}([0,t_{f}]),$ (21) $\displaystyle\lim_{k\to+\infty}J_{k}=J\
\text{in}\ L^{\infty}([0,t_{f}]).$

###### Proof:

The proof of the uniform converge of $\theta_{k}$ to $\theta$ is analogous to
one of Proposition 2. While for the convergence of $J_{k}$

$\begin{split}\|J_{k}-J\|_{\infty}&\leq\int_{0}^{t_{f_{k}}}A|u_{k}^{2}(s)-u_{(}s)|\,ds\\\
&+\int_{0}^{t_{f_{k}}}B|\theta_{k}^{2}(s)-\theta(s)^{2}|\,ds\\\ &\leq
A\|u_{k}+u\|_{\infty}\int_{0}^{t_{f_{k}}}|u_{k}(s)-u_{(}s)|\,ds\\\
&+B\|\theta_{k}+\theta\|_{\infty}\int_{0}^{t_{f_{k}}}|\theta_{k}-\theta|\,ds.\end{split}$

Therefore we conclude for the $L^{1}$ convergence of $u_{k}$ to $u$ and the
$L^{\infty}$ convergence of $\theta_{k}$ to $\theta$ ∎

###  IV-A Case $B=0$

Let us now consider a particular case of the optimal control problem (17) i.e
we fix $B=0$. Also in this case we can split the problem into two subproblems.
Considering the first time interval the cost is

$J[u]=\int_{0}^{t_{1}}u(s)^{2}\,ds.$

Notice that since $\dot{\theta}=u$ and using the Jensen inequality we have the
following

$J[u]=\int_{0}^{t_{1}}\dot{\theta}(s)^{2}\,ds\geq
t_{1}\Bigl{(}\frac{1}{t_{1}}\int_{0}^{t_{1}}\dot{\theta}(s)\,ds\Bigr{)}^{2}=\frac{(\theta_{1}-\theta_{0})^{2}}{t_{1}}.$
(22)

Moreover, since $u\in[-\varepsilon,\varepsilon]$, the equality holds only if
the control $u$ is constant, either $-\varepsilon$ or $\varepsilon$. Therefore
we conclude that in this case the problem is equivalent to the minimal time
one with

$u(t)=\begin{cases}\varepsilon&\text{if}\quad\theta_{1}>\theta_{0}\\\
-\varepsilon&\text{if}\quad\theta_{0}>\theta_{1}\end{cases}$ (23)

and consequently $t_{1}=(\theta_{1}-\theta_{0})/\varepsilon.$ The same
considerations are valid for the second time interval, therefore the entire
problem is equivalent to the minimal time one treated in Section III. This
means that if we constrain the control in the compact set
$[-\varepsilon,\varepsilon]$ the optimal strategy for minimizing the energy is
equivalent to the one which minimize the time to go from the initial to the
final configuration.

###  IV-B Case on $n$ switchings

Recall that for the problem (7) and (18) we analyse the simple case with only
a switching. Let now suppose to perform $n$ switches to go from $x_{0}$ to
$x_{f}$. Therefore we divide $\Delta x$ in $n$ sub-intervals assuming that
only one switching is required to move from the start point to the end point
of each one. Thus

$J^{n}(u,t)=\sum_{i=1}^{n}\int_{0}^{t_{i}}\ell_{0}(\theta,u,s)\,ds+n$

with $\ell_{0}(\theta,u,s)=1$ in the minimum time case and
$\ell_{0}(\theta,u,s)=u(s)^{2}$ in the minimum energy one, is the total cost
to minimize using $n$ switchings with $t_{i}$ the minimum time used to cross
the subinterval $i$. Then the following hold

###### Proposition 4

$\min_{n,u}J^{n}(u,t)\leq\min_{u}\left\\{J^{1}(u,t)+1\right\\},$

where $J^{1}(u,t)$ is the cost to minimize using only one switching.

###### Proof:

The proof follows by the definition of minimum of a set. ∎

In computing the explicit solution both of the minimum time problem and the
minimum quadratic cost one we have supposed that to reach $x_{f}$ from $x_{0}$
only one switching in the dynamics is required. Now, we will show numerically
that our assumption on the number of switching is also the optimal strategy
for both the studied problems. We fix $\Delta x=x_{f}-x_{0}=10$ and use the
following parameters

  * •

$a=10,\ b=0.1$, $\xi=1,\ \eta=2$;

  * •

$m=1,\ m_{11}=a\pi,\ m_{22}=b\pi$;

  * •

$\theta_{0}=\pi/6,\ A=1,\ B=0,\ n=30$.

We divide $\Delta x$ into $n$ subintervals in each of them we compute the
optimal time (and the control $L^{2}$ norm) and then we sum them up obtaining
the total time (and the total quadratic cost) needed to cross the interval
$[x_{0},x_{f}]$. Finally we add to the total cost the number of switchings
obtaining $J^{n}(u,t)$. Plotting each cost in function of the number $n$ of
switchings that we perform, we get that the greater the number of switchings,
the greater the time to reach $x_{f}$ is. Same conclusion also holds for the
energy. See Fig. 4.

(a) (b)

Figure 4:  Plot of the time (a) and of the cost (b) to realize the
displacement $\Delta_{x}$ as function of the number $n$ of switchings.

##  V CONCLUSIONS

In this paper we explicitly derive the solution of two different optimal
control problems (the minimum time one and the linear quadratic one) related
to the motion of the scallop swimmer with one switching in the dynamics. We
approximate the possible discontinuous solution with a suitable sequence of
periodic continuous functions. Comparing the two problem we find that the
optimal strategy for minimizing the $L^{2}$ norm of the control is equivalent
to the one which minimize the time to reach the final configuration. The more
general case with $n$ switching in the dynamics is then treated and numerical
simulations showing the relations between $n$ and a generalized cost function
$J^{n}$ are provided. In future research, inspired by the numerical results,
we plan to give an analytic proof about the above relations and study other
optimal control problems.

## ACKNOWLEDGMENT

The authors thank Fabio Bagagiolo and Antonio Marigonda for many enlighting
and helpful discussions.

## References

  * [1] F. Alouges, A. DeSimone, L. Giraldi, and M. Zoppello, “Self-propulsion of slender micro-swimmers by curvature control: N-link swimmers”, International Journal of Non-Linear Mechanics, vol. 56 pp.132-141, 2013.
  * [2] F. Alouges, A. DeSimone, A. Lefebre, “Biological Fluid Dynamics:Swimming at low Reynolds number”, _Springer of Complexity and System Science_ , 2009
  * [3] A. A. Agrachev, _Non linear and optimal control theory_ , Springer Verlag, 2008.
  * [4] F. Bagagiolo, R. Maggistro, “Hybrid Thermostatic Approximation of Junctions for some Optimal Control Problems on Networks”, arxiv.org/abs/1710.07160
  * [5] F. Bessel _Investigations on the length of the seconds pendulum_ Königlichen Akademie der Wissenschaften, Berlin, (1828).
  * [6] T. Chambrion and A. Munnier, “Locomotion and control of a self-propelled shape-changing body in a fluid”, _J. Nonlin. Sci._ , vol. 21, pp. 325-385, 2011.
  * [7] J. Y. Cheng , E. DeMont, “Hydrodynamics of scallop locomotion: unsteady fluid forces on clapping shells”, _J Fluid Mech_ , vol. 317, pp. 73-90, 1996.
  * [8] S. Childress, _Mechanics of swimming and flying_ , Cambridge University Press, Cambridge 1981.
  * [9] M. J. Lighthill, _Mathematical biofluiddynamics_ , Society of Industrial and Applied Mathematics, Philadelphia, 1975.
  * [10] J. Gray and J. Hancock, “The propulsion of sea-urchin spermatozoa”, Journal of Experimental Biology, vol. 32, pp. 802-814, 1955.
  * [11] L. Giraldi, P. Martinon, and M. Zoppello, “Controllability and optimal strokes for n-link microswimmer”, _IEEE 52th Conference on Decision and Control_ , 2013.
  * [12] L. Giraldi, P. Martinon, and M. Zoppello, “Optimal design of the three-link purcell swimmer”, Phys. Rev. E, vol. 91, 2015.
  * [13] R. Golestanian, A. Ajdari, “Analytic results for the three-sphere swimmer at low Reynolds number”, _Phys Rev E_ , vol. 7, 036308, 2008.
  * [14] E. M. Purcell, “Life at low Reynolds number”, _Am J Phys_ , vol. 45, pp. 3-11, 1977.
  * [15] T. Qiu et al, “Swimming by reciprocal motion at low Reynolds number”, _Nat Commun_ vol. 5, 5119. doi:10.1038/ncomms6119, 2014.
  * [16] F. Bagagiolo, R. Maggistro, M. Zoppello, “Swimming by switching”, _Meccanica_ , vol. 52 , no. 14, pp. 3499-3511, 2017.
  * [17] A. Visintin, _Differential models of hysteresis_ , Springer, Heidelberg, 1994.
  * [18] A. Bressan, “Impulsive control of Lagrangian systems and locomotion in fluids”, _Discrete Contin Dyn Syst_ , vol. 20, no. 1, pp. 1-35, 2008.
  * [19] E. Trélat, _Contrôle optimal, Théorie et applications_ , Vuibert, Paris, 2005.
  * [20] F. Cardin, M. Zoppello, “Swim-like motion of bodies immersed in an ideal fluid”, _ESAIM COCV_ , doi.org/10.1051/cocv/2017028, 2017.