# Fast microscale acoustic streaming driven by a temperature-gradient-induced
non-dissipative acoustic body force

Wei Qiu<EMAIL_ADDRESS>Department of Biomedical Engineering, Lund
University, Ole Römers väg 3, 22363 Lund, Sweden  Jonas Helboe Joergensen
<EMAIL_ADDRESS>Department of Physics, Technical University of Denmark,
DTU Physics Building 309, DK-2800 Kongens Lyngby, Denmark  Enrico Corato
Department of Biomedical Engineering, Lund University, Ole Römers väg 3, 22363
Lund, Sweden  Henrik Bruus<EMAIL_ADDRESS>Department of Physics,
Technical University of Denmark, DTU Physics Building 309, DK-2800 Kongens
Lyngby, Denmark  Per Augustsson<EMAIL_ADDRESS>Department of
Biomedical Engineering, Lund University, Ole Römers väg 3, 22363 Lund, Sweden

(16 March 2021)

###### Abstract

We study acoustic streaming in liquids driven by a non-dissipative acoustic
body force created by light-induced temperature gradients. This thermoacoustic
streaming produces a velocity amplitude approximately 50 times higher than
boundary-driven Rayleigh streaming and 90 times higher than Rayleigh-Bénard
convection at a temperature gradient of 10 K/mm in the channel. Further,
Rayleigh streaming is altered by the acoustic body force at a temperature
gradient of only 0.5 K/mm. The thermoacoustic streaming allows for modular
flow control and enhanced heat transfer at the microscale. Our study provides
the groundwork for studying microscale acoustic streaming coupled with
temperature fields.

Acoustic streaming describes the steady time-averaged fluid motion that takes
place when acoustic waves propagate in viscous fluids. The streaming flow is
driven by a nonzero divergence in the time-averaged acoustic momentum-flux-
density tensor Lighthill (1978). Conventionally in a homogeneous fluid, this
nonzero divergence arises from two dissipation mechanisms of acoustic energy.
The first case is the so-called boundary-driven Rayleigh streaming Lord
Rayleigh (1884), in which acoustic energy is dissipated in viscous boundary
layers where the velocity of the oscillating fluid changes to match the
surface velocity of the channel walls Nyborg (1958); Hamilton  _et al._
(2003a) or of suspended objects Elder (1959); Lee and Wang (1990); Yarin  _et
al._ (1999); Tho  _et al._ (2007). This gives rise to velocity gradients that
drive the flow Schlichting (1932), and is typically observed in standing wave
fields in systems of a size comparable to the wavelength Muller  _et al._
(2013). In the second type of streaming, known as quartz wind or bulk-driven
Eckart streaming Eckart (1948), high acoustic wave attenuation yields the
gradients that drive the streaming flow, which is typically associated with
high-frequency traveling waves Eckart (1948); Nyborg (1953); Riley (2001).
Both cases have been extensively studied, and the phenomenon of acoustic
streaming continues to attract attention due its importance in processes
related to medical ultrasound Marmottant and Hilgenfeldt (2003); van der Sluis
_et al._ (2007); Wu and Nyborg (2008); Doinikov and Bouakaz (2010),
thermoacoustic engines Bailliet  _et al._ (2001); Hamilton  _et al._ (2003b),
acoustic levitation Trinh and Robey (1994); Yarin  _et al._ (1999), and
manipulation of particles and cells in microscale acoustofluidics Bruus  _et
al._ (2011); Barnkob  _et al._ (2012); Hammarström  _et al._ (2012); Collins
_et al._ (2015); Marin  _et al._ (2015); Hahn  _et al._ (2015); Guo  _et al._
(2016).

Recently, we discovered that boundary-driven streaming can be significantly
suppressed in inhomogeneous media formed by different solute molecules Karlsen
_et al._ (2018); Qiu  _et al._ (2019). This is attributed to the acoustic body
force $\bm{f}_{\mathrm{ac}}$, also originating from the nonzero divergence in
the time-averaged momentum-flux-density tensor, due to gradients in density
and compressibility in the fluid Karlsen  _et al._ (2016); Augustsson  _et
al._ (2016), which competes with the boundary-layer streaming stress. The
streaming rolls are confined to narrow regions near the channel walls, before
the inhomogeneity smears out due to diffusion and advection of the solute.
This effect enables the acoustic manipulation of submicrometer particles Qiu
_et al._ (2020) of biological relevance such as bacteria Van Assche  _et al._
(2020), as well as trapping of hot plasma in gasses Koulakis  _et al._ (2018).

In this Letter, we investigate microscale acoustic streaming in a liquid, in
which the temperature-dependent compressibility and density have been made
inhomogeneous by introducing a sustained temperature gradient. We generate
this gradient by light irradiation and absorption, and subsequently measure
the streaming driven by the temperature-gradient-induced acoustic body force
and call it thermoacoustic streaming. Using our newly developed model for
thermovisocus acoustofluidics Joergensen and Bruus (2020), the experimental
results are validated and the mechanisms responsible for the thermoacoustic
streaming are explained.

Our main findings are: (i) The thermoacoustic streaming begins to disturb the
boundary-driven Rayleigh streaming with a temperature gradient as small as 0.5
K/mm, resulting in streaming rolls with complex three-dimensional (3D)
patterns. (ii) With a temperature gradient of 10 K/mm, the thermoacoustic
streaming velocity is about 50 and 90 times higher than that of the boundary-
driven Rayleigh streaming and of the Rayleigh-Bénard convection, respectively.
(iii) In contrast to other types of acoustic streaming, the mechanism driving
the thermoacoustic streaming is non-dissipative.

The thermoacoustic streaming should be of considerable fundamental interest to
a broad community of researchers working in nonlinear acoustics,
thermoacoustics, microscale acoustofluidics, as well as heat transfer. For a
microsystem, the advective streaming flow is remarkably high compared to the
rate of thermal diffusion, and with a Péclet number Pe $\approx 1$, heat
transfer is strongly enhanced. Further, our findings pave the way for
transient or steady control of the streaming through modulations of the
temperature or the acoustic field.

Figure 1:  Sketch of the acoustofluidic silicon chip (light gray) sandwiched
between two glass layers (white) that allows optical transparency for both
heating and particle tracking. The light emitted from a blue LED below the
chip is absorbed by the water containing dye molecules and a temperature
gradient forms from low (orange) to high (yellow) temperature in the channel.
The piezoelectric transducer (dark brown) excites the resonant half-wave
pressure field $p_{1}$ (purple) at 953 kHz. A green LED shines light from
above to excite the red fluorescent light from the tracer particles, which
allows the optical recording of the tracer bead motion in a part of length
$L=1300~{}\textrm{\textmu{}m}$ of the channel of width
$W=760~{}\textrm{\textmu{}m}$, and height $H=370~{}\textrm{\textmu{}m}$.

Experimental method.— The experiments were performed using a long straight
microchannel of width $W=760\,\textrm{\textmu{}m}$ and height
$H=370\,\textrm{\textmu{}m}$ in a glass-silicon-glass chip with a
piezoelectric transducer glued underneath. The temperature gradient inside the
channel was generated by directing the focused light from a 488-nm light-
emitting-diode (LED) through water with a small amount of added dye (0.1 wt%
Orange G) that has an overlapping absorption peak with the LED, see Fig. 1,
which results in 99% of the LED light being absorbed in the liquid. The
transducer was driven at a frequency of 953 kHz with an input power of 88 mW,
which produced a standing half wave across the width with the acoustic energy
density $E_{\mathrm{ac}}=9.24$ J/m3 at room temperature. The induced streaming
was measured using general defocusing particle tracking Barnkob  _et al._
(2015); Barnkob and Rossi (2020) at 5 to 60 fps with 1.1 µm-diameter
polystyrene tracer particles (red fluorescence). After a band-pass filter
(525-550 nm), the excitation light of the tracer particles from a green LED is
barely absorbed by the dye molecules and hence does not affect the temperature
gradient in the channel. The measurements under each condition were repeated
13 to 27 times and recorded in 7800 to 40500 frames to improve the statistics.
The temperature field in the channel mid-height plane was imaged using
temperature-sensitive fluorescent dye (Rhodamine B), see the Supplemental
Material 111See Supplemental Material at [url] for details about the
temperature measurements, the comparison between measured and simulated
temperature fields, the numerical model, as well as the measured and simulated
Rayleigh-Bénard convection..

Numerical model.— The model is a pressure acoustics model including
thermoviscous boundary layers as presented in Ref. Joergensen and Bruus
(2020), which enables 3D simulations of thermoviscous acoustofluidic devices.
To simulate the long glass-silicon-glass chip, symmetry planes are exploited
to only simulate a quarter of the chip. Furthermore, the perfectly matched
layer (PML) technique Bermúdez  _et al._ (2007) is used to avoid simulating
the entire length of the chip. The solver consists of three steps: (i)
Computing the temperature field $T_{0}$ induced by the LED with an amplitude
set to match the observed temperature gradients, (ii) computing the acoustic
displacement in the solid $\bm{u}_{1}$ and the pressure $p_{1}$ in the fluid
due to an actuation on the glass, and (iii) computing the resulting acoustic
streaming field $\bm{v}_{2}$. The heating from the LED is modeled with no
absorption in the glass and an absorption parameter in the fluid selected to
absorb $99\%$ of the light passing through the chip as measured in the
experiment. The model is based on perturbation theory, but the highest
streaming velocities recorded in the experiments is found to be at the limit
of the validity of the model, because there the thermoacoustic streaming
begins to alter the temperature field $T_{0}$. For more details on the
numerical model see the Supplemental Material Note (1). Due to the inherent
difficulty in measuring energy density $E_{\mathrm{ac}}$ at high streaming
velocities when temperature gradients are present, the $E_{\mathrm{ac}}$ used
in simulations is obtained by matching the experimental streaming velocity
amplitude for each temperature gradient.

Figure 2:  (a)-(e) The measured (exp, left half) and simulated (sim, right
half) streaming patterns for $G=9.76$ K/mm ($\Delta T_{0}=3.71$ K across the
channel width $W$), observed in different planes $x$-$y$ horizontally, $x$-$z$
vertically, and $y$-$z$ vertically. The vector plot (magenta) in a given plane
is the in-plane velocity and the color plot is its magnitude from 0 (dark
blue) to 1042 $\textrm{\textmu{}m}\,\mathrm{s}^{-1}$ (yellow). The simulation
is performed at $E_{\mathrm{ac}}=23$ J/m3. Spatial bins with no experimental
data are excluded (gray). The pair of curved white lines in (b) and (d)
represent the centerlines of the two counter-rotating deformed 3D cylindrical
streaming flow rolls. The two line-plot insets show the measured (purple) and
the simulated (green) $x$ (or $y$) component $v_{2}^{x}$ (or $v_{2}^{y}$) of
the velocity along the red lines. (f) Color plot from 25.0
${}^{\circ}\\!\textrm{C}$ (black) to 30.1 ${}^{\circ}\\!\textrm{C}$ (white) of
the measured and simulated temperature $T_{0}$ in the horizontal $x$-$y$ plane
at $z=0$, see more details in the Supplemental Material Note (1). The regions
where the fluorescence intensity is affected by the channel sidewalls are
excluded (gray). (g) $T_{0}$ as in panel (f) but for the vertical $y$-$z$
plane at $x=0$, and only the simulation results are shown.

Results and discussion.— When both sound field and temperature gradients are
present, the streaming flow exhibits a complex 3D pattern. An example is shown
in Fig. 2(a)-(e), corresponding to a temperature difference $\Delta T_{0}$ =
3.71 K across the channel of width $W=760~{}\textrm{\textmu{}m}$, equivalent
to a gradient $G=2\Delta T_{0}/W=9.76$ K/mm. Here, two counter-rotating
deformed 3D cylindrical streaming flow rolls appear, whirling with a velocity
amplitude $|{\bm{v}_{2}}|=1074~{}\textrm{\textmu{}m}$/s around the pair of
curved white centerlines shown in Fig. 2(b) and (d). This velocity amplitude
is about 50 and 90 times higher than that of the boundary-driven Rayleigh
streaming and the Rayleigh-Bénard convection, respectively, under the same
driving conditions (see Supplemental Material Note (1)).

The generation of this fast streaming can be explained by the acoustic body
force $\bm{f}_{\mathrm{ac}}$ due to the temperature field induced by the blue
LED. In this experiment, the light heats the fluid from beneath, while the
silicon walls efficiently transport the heat away, thus cooling the sides of
the channel. Temperature gradients are therefore induced in all directions: In
the $x$-$y$ plane by the Gaussian profile of the light intensity and by the
cooling from the silicon sidewalls, and in the $z$-direction by light
absorption fulfilling the Beer-Lambert law. The resulting temperature field is
highest at the center of the channel bottom, as shown by the measured and
simulated temperature fields in the horizontal $x$-$y$ plane around channel
mid-height $z=0\,\textrm{\textmu{}m}$ in Fig. 2(f), and by the simulated
temperature field in the vertical $y$-$z$ cross section at
$x=0\,\textrm{\textmu{}m}$ in Fig. 2(g). The acoustic body force
$\bm{f}_{\mathrm{ac}}$ depends on the gradients in compressibility and
density, the acoustic pressure $p_{1}$, and the acoustic velocity $\bm{v}_{1}$
Karlsen  _et al._ (2016). When the inhomogeneities are created by a
temperature field, $\bm{f}_{\mathrm{ac}}$ can be expressed as,

$\displaystyle\bm{f}_{\mathrm{ac}}$
$\displaystyle=-\frac{1}{4}|p_{1}|^{2}\bm{\nabla}\kappa_{T,0}-\frac{1}{4}|\bm{v}_{1}|^{2}\bm{\nabla}\rho_{0}$
$\displaystyle=-\frac{1}{4}\Bigg{[}|p_{1}|^{2}\bigg{(}\frac{\partial\kappa_{T}}{\partial
T}\bigg{)}_{T_{0}}+|\bm{v}_{1}|^{2}\bigg{(}\frac{\partial\rho}{\partial
T}\bigg{)}_{T_{0}}\Bigg{]}\bm{\nabla}T_{0}.$ (1)

Three factors determine the action of $\bm{f}_{\mathrm{ac}}$ on the fluid. (i)
Both the compressibility and density decrease with temperature, thus
$\bm{f}_{\mathrm{ac}}$ points towards the high temperature region here the
center of the channel heated by the LED. (ii) At room temperature,
$\kappa_{T}|p_{1}|^{2}\approx\rho|\bm{v}_{1}|^{2}$ and
$\frac{1}{\kappa_{T}}|\partial_{T}\kappa_{T}|\gg\frac{1}{\rho}|\partial_{T}\rho|$,
so $\bm{f}_{\mathrm{ac}}$ is dominated by the $|p_{1}|^{2}$ compressibility
term and thus is strongest at the pressure anti-nodes at the sides of the
channel. (iii) As shown in Fig. 2(g), the temperature gradient is larger at
the bottom than at the top of the channel resulting in a stronger
$\bm{f}_{\mathrm{ac}}$ at the bottom. Consequently, in the bottom part of the
LED spot, $\bm{f}_{\mathrm{ac}}$ pulls the fluid horizontally inward to the
vertical $x$-$z$ center plane at $y=0$ and by mass conservation lets it escape
outward along the axial $x$-direction and upward along the vertical
$z$-direction. The resulting streaming flow contains the two aforementioned
deformed 3D cylindrical flow rolls, which when projected onto horizontal and
vertical planes appear as the four horizontal streaming rolls in Fig. 2(a)-(c)
strongest in the center-plane $z=0$, and as the two vertical streaming rolls
in Fig. 2(d)-(e).

Figure 3:  (a) Measured (exp, left half) and simulated (sim, right half) in-
plane streaming velocity $\bm{v}_{2}^{yz}$ (magenta vectors) and its magnitude
$|\bm{v}_{2}^{yz}|$ from 0 (dark blue) to its maximum (yellow) in the vertical
$y$-$z$ plane averaged over the vertical slab $x=0\pm
100~{}\textrm{\textmu{}m}$ for six temperature differences across the channel
ranging from $\Delta T_{0}=0$ to 3.71 K. Simulations are performed with the
energy density $E_{\mathrm{ac}}$ for which $\left|{\bm{v}_{2}^{yz}}\right|$
matches the experimental one under each $\Delta T_{0}$. Spatial bins with no
experimental data points are excluded (gray). (b) Dependence of vortex size
$\Delta$ (red) and the velocity in the $z$-direction
$\left|{\bm{v}_{2}^{z}}\right|$ (blue) on $\Delta T_{0}$. Simulated $\Delta$
is indicated by a red solid line. Light blue region shows the simulated
$\left|{\bm{v}_{2}^{z}}\right|$ under $E_{\mathrm{ac}}$ from 9.24 Pa (lower
bound) to 23 Pa (upper bound) with a reference (solid blue line) obtained at
$E_{\mathrm{ac}}=$ 18 Pa.

In contrast to the inhomogeneity created by diffusing solute molecules out of
equilibrium Karlsen  _et al._ (2018); Qiu  _et al._ (2019), the temperature
gradients here are sustained because the system is run in steady state driven
by the PZT transducer and the LED light. The resulting streaming is more than
one order of magnitude faster than the boundary-driven Rayleigh streaming, and
it is mainly due to the large non-dissipative $\bm{f}_{\mathrm{ac}}$. The
thermoacoustic driving mechanism is fundamentally different from the
dissipation mechanism of the conventional forms of acoustic streaming.

The transition from boundary-driven Rayleigh streaming to thermoacoustic
streaming is studied in the vertical $y$-$z$ cross section
($-100\,\textrm{\textmu{}m}<x<100\,\textrm{\textmu{}m}$) by gradually
increasing the output power of the LED. Following Refs. Karlsen  _et al._
(2018); Qiu  _et al._ (2019), we quantify the streaming evolution by the
vortex size $\Delta$, defined as the average of the distance from the center
of each of the two upper flow rolls to the channel ceiling at
$z=\frac{1}{2}H$. Figure 3 shows that at a zero temperature gradient $G$, the
streaming is governed by the four conventional boundary-driven Rayleigh
streaming rolls, whereas at high gradient $G\approx 3.6$ K/mm ($\Delta
T_{0}=1.36$ K), only two big temperature-gradient induced flow rolls driven by
the relatively large $\bm{f}_{\mathrm{ac}}$ is present, occupying the whole
channel cross section. In transitioning from the former to the latter
situation, the two top Rayleigh flow rolls appear to expand downwards
squeezing the bottom rolls against the channel floor at $z=-\frac{1}{2}H$.
This phenomenon can be explained by the fact that the two top (bottom)
Rayleigh flow rolls have the same (opposite) rotation direction as the two
temperature-gradient induced flow rolls. Already at $G=0.5$ K/mm ($\Delta
T_{0}=0.170$ K), $\bm{f}_{\mathrm{ac}}$ is large enough to distort the four-
flow-roll Rayleigh streaming pattern. When $\Delta T_{0}$ further increases,
the two-flow-roll thermoacoustic streaming pattern dominates, and eventually
remains unchanged, while the velocity amplitude increases almost linearly for
the investigated range, see Fig. 3(b).

The observation that thermoacoustic streaming occurs already at temperature
gradients below 0.5 K/mm calls for caution when combining acoustofluidic
devices with optical systems. For an absorbing liquid, the light in a standard
microscope is enough to induce strong velocity fields that can interfere with
the study object. While we did not record the transient build-up of the
streaming field upon activating the light, it can be noted that the
development of the temperature field, and thus the streaming field, occurs
within a few hundred milliseconds at the studied length scale which opens for
rapid modulation of local streaming fields through fast re-configurable
optical fields. Even so, the induced streaming velocity is high enough to
match the thermal diffusion time and thereby impact the heat transfer in the
structure.

Conclusion.—This letter describes a comprehensive experimental and numerical
study of the thermoacoustic streaming in liquids induced by temperature
gradients generated by light absorption in a microchannel. We have obtained a
good match between measured and simulated velocity fields in 3D. As summarized
by the main findings (i)-(iii) in the introduction, the thermoacoustic
streaming, driven by the non-dissipative acoustic body force, has a markedly
different origin than that of the conventional acoustic streaming associated
with energy dissipation. Moreover, it reaches much higher velocity amplitudes
compared to the boundary-driven Rayleigh streaming and the Rayleigh-Bénard
convection under comparable conditions. The acoustic body force relies on the
acoustic field and the gradients in compressibility and density, which is
analogous to the driving force of the classical Rayleigh-Bénard convection
relying on the gravitational field and the gradient in density.

By including the temperature dependence of both density and compressibility in
this work through our theory in Ref. Joergensen and Bruus (2020), our analysis
of thermoacoustic streaming in terms of the acoustic body force is valid for
both liquids and gases. Thus, we have extended previous related work on gases,
where compressibility effects are unimportant and therefore neglected Hamilton
_et al._ (2003b); Aktas and Ozgumus (2010); Chini  _et al._ (2014); Červenka
and Bednařík (2017); Michel and Chini (2019); Fand and Kaye (1960); Thompson
and Atchley (2005); Thompson  _et al._ (2005); Nabavi  _et al._ (2008).

This study is fundamental in scope, but also demonstrates a method with clear
potential for controlling local flows in microchannels. Further, we highlight
important implications of this phenomenon relating to heat transfer and
integration of optical fields with microscale acoustofluidic devices.

We are grateful to R. Barnkob (Technical University of Munich) and M. Rossi
(Technical University of Denmark) for providing the software GDPTlab. W.Q. was
supported by the Foreign Postdoctoral Fellowship from Wenner-Gren Foundations
and by MSCA EF Seal of Excellence IF-2018 from Vinnova, Sweden’s Innovation
Agency (Grant No. 2019-04856). This work has received funding from Independent
Research Fund Denmark, Natural Sciences (Grant No. 8021-00310B) and the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (Grant Agreement No. 852590).

## References

  * Lighthill (1978) J. Lighthill, Acoustic streaming, J Sound Vib. 61, 391 (1978).
  * Lord Rayleigh (1884) Lord Rayleigh, On the circulation of air observed in Kundt’s tubes, and on some allied acoustical problems, Philos. Trans. R. Soc. London 175, 1 (1884).
  * Nyborg (1958) W. L. Nyborg, Acoustic streaming near a boundary, J. Acoust. Soc. Am. 30, 329 (1958).
  * Hamilton  _et al._ (2003a) M. Hamilton, Y. Ilinskii, and E. Zabolotskaya, Acoustic streaming generated by standing waves in two-dimensional channels of arbitrary width, J. Acoust. Soc. Am. 113, 153 (2003a).
  * Elder (1959) S. A. Elder, Cavitation microstreaming, J. Acoust. Soc. Am. 31, 54 (1959).
  * Lee and Wang (1990) C. P. Lee and T. G. Wang, Outer acoustic streaming, J. Acoust. Soc. Am. 88, 2367 (1990).
  * Yarin  _et al._ (1999) A. Yarin, G. Brenn, O. Kastner, D. Rensink, and C. Tropea, Evaporation of acoustically levitated droplets, J. Fluid Mech. 399, 151 (1999).
  * Tho  _et al._ (2007) P. Tho, R. Manasseh, and A. Ooi, Cavitation microstreaming patterns in single and multiple bubble systems, J. Fluid Mech. 576, 191 (2007).
  * Schlichting (1932) H. Schlichting, Berechnung ebener periodischer grenzeschichtströmungen, Phys. Z. 33, 327 (1932).
  * Muller  _et al._ (2013) P. B. Muller, M. Rossi, A. G. Marin, R. Barnkob, P. Augustsson, T. Laurell, C. J. Kähler, and H. Bruus, Ultrasound-induced acoustophoretic motion of microparticles in three dimensions, Phys. Rev. E 88, 023006 (2013).
  * Eckart (1948) C. Eckart, Vortices and streams caused by sound waves, Phys. Rev. 73, 68 (1948).
  * Nyborg (1953) W. L. Nyborg, Acoustic streaming due to attenuated plane waves, J. Acoust. Soc. Am. 25, 68 (1953).
  * Riley (2001) N. Riley, Steady streaming, Annu. Rev. Fluid Mech. 33, 43 (2001).
  * Marmottant and Hilgenfeldt (2003) P. Marmottant and S. Hilgenfeldt, Controlled vesicle deformation and lysis by single oscillating bubbles, Nature 423, 153 (2003).
  * van der Sluis  _et al._ (2007) L. van der Sluis, M. Versluis, M. Wu, and P. Wesselink, Passive ultrasonic irrigation of the root canal: a review of the literature, Int. Endod. J. 40, 415 (2007).
  * Wu and Nyborg (2008) J. Wu and W. L. Nyborg, Ultrasound, cavitation bubbles and their interaction with cells, Adv. Drug Delivery Rev. 60, 1103 (2008).
  * Doinikov and Bouakaz (2010) A. A. Doinikov and A. Bouakaz, Theoretical investigation of shear stress generated by a contrast microbubble on the cell membrane as a mechanism for sonoporation, J. Acoust. Soc. Am. 128, 11 (2010).
  * Bailliet  _et al._ (2001) H. Bailliet, V. Gusev, R. Raspet, and R. A. Hiller, Acoustic streaming in closed thermoacoustic devices, J. Acoust. Soc. Am. 110, 1808 (2001).
  * Hamilton  _et al._ (2003b) M. Hamilton, Y. Ilinskii, and E. Zabolotskaya, Thermal effects on acoustic streaming in standing waves, J. Acoust. Soc. Am. 114, 3092 (2003b).
  * Trinh and Robey (1994) E. Trinh and J. Robey, Experimental study of streaming flows associated with ultrasonic levitators, Phys. Fluids 6, 3567 (1994).
  * Bruus  _et al._ (2011) H. Bruus, J. Dual, J. Hawkes, M. Hill, T. Laurell, J. Nilsson, S. Radel, S. Sadhal, and M. Wiklund, Forthcoming lab on a chip tutorial series on acoustofluidics: Acoustofluidics-exploiting ultrasonic standing wave forces and acoustic streaming in microfluidic systems for cell and particle manipulation, Lab Chip 11, 3579 (2011).
  * Barnkob  _et al._ (2012) R. Barnkob, P. Augustsson, T. Laurell, and H. Bruus, Acoustic radiation- and streaming-induced microparticle velocities determined by microparticle image velocimetry in an ultrasound symmetry plane, Phys. Rev. E 86, 056307 (2012).
  * Hammarström  _et al._ (2012) B. Hammarström, T. Laurell, and J. Nilsson, Seed particle enabled acoustic trapping of bacteria and nanoparticles in continuous flow systems, Lab Chip 12, 4296 (2012).
  * Collins  _et al._ (2015) D. J. Collins, B. Morahan, J. Garcia-Bustos, C. Doerig, M. Plebanski, and A. Neild, Two-dimensional single-cell patterning with one cell per well driven by surface acoustic waves, Nat. Commun. 6, 8686 (2015).
  * Marin  _et al._ (2015) A. Marin, M. Rossi, B. Rallabandi, C. Wang, S. Hilgenfeldt, and C. J. Kähler, Three-dimensional phenomena in microbubble acoustic streaming, Phys. Rev. Applied 3, 041001(R) (2015).
  * Hahn  _et al._ (2015) P. Hahn, I. Leibacher, T. Baasch, and J. Dual, Numerical simulation of acoustofluidic manipulation by radiation forces and acoustic streaming for complex particles, Lab Chip 15, 4302 (2015).
  * Guo  _et al._ (2016) F. Guo, Z. Mao, Y. Chen, Z. Xie, J. P. Lata, P. Li, L. Ren, J. Liu, J. Yang, M. Dao, S. Suresh, and T. J. Huang, Three-dimensional manipulation of single cells using surface acoustic waves, PNAS 113, 1522 (2016).
  * Karlsen  _et al._ (2018) J. T. Karlsen, W. Qiu, P. Augustsson, and H. Bruus, Acoustic streaming and its suppression in inhomogeneous fluids, Phys. Rev. Lett. 120, 054501 (2018).
  * Qiu  _et al._ (2019) W. Qiu, J. T. Karlsen, H. Bruus, and P. Augustsson, Experimental characterization of acoustic streaming in gradients of density and compressibility, Phys. Rev. Appl. 11, 024018 (2019).
  * Karlsen  _et al._ (2016) J. T. Karlsen, P. Augustsson, and H. Bruus, Acoustic force density acting on inhomogeneous fluids in acoustic fields, Phys. Rev. Lett. 117, 114504 (2016).
  * Augustsson  _et al._ (2016) P. Augustsson, J. T. Karlsen, H.-W. Su, H. Bruus, and J. Voldman, Iso-acoustic focusing of cells for size-insensitive acousto-mechanical phenotyping, Nat. Commun. 7, 11556 (2016).
  * Qiu  _et al._ (2020) W. Qiu, H. Bruus, and P. Augustsson, Particle-size-dependent acoustophoretic motion and depletion of micro-and nano-particles at long timescales, Phys. Rev. E 102, 013108 (2020).
  * Van Assche  _et al._ (2020) D. Van Assche, E. Reithuber, W. Qiu, T. Laurell, B. Henriques-Normark, P. Mellroth, P. Ohlsson, and P. Augustsson, Gradient acoustic focusing of sub-micron particles for separation of bacteria from blood lysate, Sci. Rep. 10, 3670 (2020).
  * Koulakis  _et al._ (2018) J. P. Koulakis, S. Pree, A. L. F. Thornton, and S. Putterman, Trapping of plasma enabled by pycnoclinic acoustic force, Phys. Rev. E 98, 043103 (2018).
  * Joergensen and Bruus (2020) J. H. Joergensen and H. Bruus, Theory of pressure acoustics with thermoviscous boundary layers and streaming in elastic cavities, arXiv.org preprint, 2012.07997 (2020), http://arxiv.org/abs/2012.07997.
  * Barnkob  _et al._ (2015) R. Barnkob, C. J. Kähler, and M. Rossi, General defocusing particle tracking, Lab Chip 15, 3556 (2015).
  * Barnkob and Rossi (2020) R. Barnkob and M. Rossi, General defocusing particle tracking: fundamentals and uncertainty assessment, Exp. Fluids 61, 110 (2020).
  * Note (1) See Supplemental Material at [url] for details about the temperature measurements, the comparison between measured and simulated temperature fields, the numerical model, as well as the measured and simulated Rayleigh-Bénard convection.
  * Bermúdez  _et al._ (2007) A. Bermúdez, L. Hervella-Nieto, A. Prieto, and R. Rodríguez, An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems, J. Comput. Phys. 223, 469 (2007).
  * Aktas and Ozgumus (2010) M. K. Aktas and T. Ozgumus, The effects of acoustic streaming on thermal convection in an enclosure with differentially heated horizontal walls, Int. J. Heat and Mass Transf. 53, 5289 (2010).
  * Chini  _et al._ (2014) G. Chini, Z. Malecha, and T. Dreeben, Large-amplitude acoustic streaming, J. Fluid Mech. 744, 329 (2014).
  * Červenka and Bednařík (2017) M. Červenka and M. Bednařík, Effect of inhomogeneous temperature fields on acoustic streaming structures in resonators, J. Acoust. Soc. Am. 141, 4418 (2017).
  * Michel and Chini (2019) G. Michel and G. P. Chini, Strong wave–mean-flow coupling in baroclinic acoustic streaming, J. Fluid Mech. 858, 536 (2019).
  * Fand and Kaye (1960) R. Fand and J. Kaye, Acoustic streaming near a heated cylinder, J. Acoust. Soc. Am. 32, 579 (1960).
  * Thompson and Atchley (2005) M. W. Thompson and A. A. Atchley, Simultaneous measurement of acoustic and streaming velocities in a standing wave using laser doppler anemometry, J. Acoust. Soc. Am. 117, 1828 (2005).
  * Thompson  _et al._ (2005) M. W. Thompson, A. A. Atchley, and M. J. Maccarone, Influences of a temperature gradient and fluid inertia on acoustic streaming in a standing wave, J. Acoust. Soc. Am. 117, 1839 (2005).
  * Nabavi  _et al._ (2008) M. Nabavi, K. Siddiqui, and J. Dargahi, Effects of transverse temperature gradient on acoustic and streaming velocity fields in a resonant cavity, Appl. Phys. Lett. 93, 051902 (2008).


#  A contract negotiation scheme for safety verification of interconnected
systems

Xiao Tan, Antonis Papachristodoulou, and Dimos V. Dimarogonas  This work was
supported in part by Swedish Research Council (VR), in part by Horizon Europe
SymAware, in part by ERC CoG LEAFHOUND, in part by EU CANOPIES Project, and in
part by Knut and Alice Wallenberg Foundation (KAW). Xiao Tan and Dimos V.
Dimarogonas are with the School of EECS, Royal Institute of Technology (KTH),
100 44 Stockholm, Sweden. Antonis Papachristodoulou is with the Department of
Engineering Science, University of Oxford, Oxford, United Kingdom. Email:
{xiaotan<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper proposes a (control) barrier function synthesis and safety
verification scheme for interconnected nonlinear systems based on assume-
guarantee contracts (AGC) and sum-of-squares (SOS) techniques. It is well-
known that the SOS approach does not scale well for barrier function synthesis
for high-dimensional systems. In this paper, we show that compositional
methods like AGC can mitigate this problem. We formulate the synthesis problem
into a set of small-size problems, which constructs local contracts for
subsystems, and propose a negotiation scheme among the subsystems at the
contract level. The proposed scheme is then implemented numerically on two
examples: vehicle platooning and room temperature regulation.

##  I Introduction

In many engineering applications, system states need to be confined to a
specific set of safe states. Designing active control to achieve this property
and verifying a given closed-loop system regarding this property are known as
safety synthesis and verification problems. Many safety-ensuring control
approaches have been proposed in the literature, including reachability
analysis [1], control barrier functions (CBF) [2], model predictive control
[3], prescribed performance control [4] among many others. In particular, when
a CBF is shown to exist, safety-ensuring feedback can be constructed, and the
safety of the system is certified [5]. Thus, there has been lots of interest
in synthesizing valid control barrier functions numerically, by, for example,
sum-of-square approaches [6, 7], learning-based approaches [8, 9], and
Hamiltonian-Jacobi reachability analysis [10]. However, most of these
approaches are limited to dynamical systems of small to moderate size, and
will become computationally intractable for large-scale systems.

Many complex, large-scale systems naturally impose an interconnected
structure. It is thus essential to exploit this structure to deal with the
numerical scalability issue. Along this line of research, the idea of
compositional reasoning has been leveraged so that one could establish
properties of the complete interconnected system by reasoning properties on
its components. As for system safety/invariance property, [11, 12] propose to
synthesize local barrier functions, establish local input-to-state safety
properties, and compose the local properties by checking a small-gain-like
condition. However, it remains unclear how to adapt local safety properties if
the condition fails. On the other hand, [13] certifies the safety property by
seeking a Lyapunov function of the interconnected system. Safety is thus
certified if a subset of the constructed Lyapunov function has no intersection
with the unsafe region. It is worth noting that the search for a Lyapunov
function is solved as a centralized semi-definite problem, and is still
computationally demanding when the size of the interconnected system becomes
larger.

In the literature of formal methods and model checking [14], the composition
of system properties is usually approached through the notion of an assume-
guarantee contract [15]. In plain words, a contract describes the behavior
that a system will exhibit (guarantees) subject to the influence of the
environment (assumptions). Originally, the main application domain of a
contract in model checking was for discrete space systems. When contracts are
applied to certify the safety of complex continuous space systems, circular
reasoning of implications might exist. This is not a trivial problem in
general, and the AGC framework is always sound only if a hierarchical
structure exists [16]. [17] introduces parameterized AGCs, laying the
foundation for finding local AGCs that can be composed of. [18] deals with
invariance properties of discrete-time linear systems. The authors show that
the composition of all local AGCs can be formulated as a linear program when
using zonotopic representation to parameterize the constraint set and input
set. In [19], the authors consider a finite transition system and propose to
determine how safe a state is by applying value iterations. The contracts are
iterated locally, yet no completeness guarantee can be asserted.

Recently, there are a few works that apply AGCs to control synthesis problems
for continuous-time systems. [20] utilizes behaviour AGC for control design
for linear systems, and [21] applies AGCs to design local feedback law under
signal temporal logic specifications.

In this work, we provide a tractable safety verification scheme for
continuous-time interconnected nonlinear systems, leveraging sum-of-square
techniques and assume-guarantee contracts. Our result is built upon [16] on
the invariance AGCs for continuous-time systems that circumvent circular
reasoning under mild assumptions. Our proposed approach consists of the
construction of local AGCs and the search for compatible AGCs. In contrast to
[13, 18], we propose to negotiate local contracts only with its neighbors, and
thus no central optimization is needed. Once a set of compatible AGCs is
returned, the safety of the interconnected system is certified. Moreover, we
show that the proposed algorithms will terminate in finite steps and find a
solution whenever one exists under relevant technical assumptions in the case
of acyclic interconnections or for homogeneous systems.

##  II Notation and Preliminaries

Notation: $\mathbb{R}^{n}$ denotes the $n$-dimensional vector space. A vector
$a=(a_{1},a_{2},\ldots,a_{n})\in\mathbb{R}^{n}$ is a column vector unless
stated otherwise. For $Z\subseteq\mathbb{R}^{n}$, we denote by $M(Z)$ the set
of continuous-time maps $z\mathrel{\mathop{\ordinarycolon}}E\to Z$, where
$E\in\\{[0,a],a\geq 0\\}\cup\\{[0,a),a>0\\}\cup\\{\mathbb{R}_{+}\\}$ is a time
interval. Given sets $X_{i}\subseteq\mathbb{R}^{n_{i}}$,
$i\in\mathcal{I}=\\{1,2,\ldots,N\\}$, the Cartesian product $X_{1}\times
X_{2}\times\cdots\times X_{N}$ is denoted by $\Pi_{i\in\mathcal{I}}X_{i}$. Let
$x\in\mathbb{R}^{n}$ be an independent variable. Denote by $\mathcal{R}[x]$
the set of polynomials in the variable $x$. We call a polynomial
$p\in\mathcal{R}[x]$ sum-of-squares if there exist polynomials
$g_{1},g_{2},\ldots,g_{N}$ in the variable $x$ such that
$p=\sum_{i=1}^{N}g_{i}^{2}.$ Denote by $\Sigma[x]$ the set of sum-of-squares
polynomials in $x$. Let
$\mathcal{R}[x_{1},x_{2},\ldots,x_{n}],\Sigma[x_{1},x_{2},\ldots,x_{n}]$
denote the sets of polynomials and SOS polynomials of independent variables
$x_{1},x_{2},\ldots,x_{n}$, respectively. Consider a directed graph
$(\mathcal{I},\mathcal{E}),\mathcal{E}\subseteq\mathcal{I}\times\mathcal{I}$.
Denote by
$N(i)=\\{j\in\mathcal{I}\mathrel{\mathop{\ordinarycolon}}(j,i)\in\mathcal{E}\\}$
the set of parent nodes of node $i$, and
$\text{Child}(i)=\\{k\in\mathcal{I}\mathrel{\mathop{\ordinarycolon}}(i,k)\in\mathcal{E}\\}$
the set of its child nodes. We call node $i$ a root node if $N(i)=\emptyset$;
node $i$ is a leaf node if $\text{Child}(i)=\emptyset$.

We first introduce the definitions of continuous-time systems, their
interconnections, and assume-guarantee contracts tailored from [16] for the
safety verification problem.

###  II-A Systems and Interconnections

In this work, we consider continuous-time systems formally defined as follows.

######  Definition 1 (Continuous-time system).

A continuous-time system $G$ is a tuple

$G=(U,W,X,Y,X^{0},\mathcal{T}),$

where the sets $U,W,X,Y,X^{0}$ represent the external input set, the internal
input set, the state set, the output set, and the initial state set,
respectively. $u\in U,w\in W,x\in X,y\in Y$ are the external input, internal
input, local state, and local output variables. $\mathcal{T}\subseteq
M(U\times W\times X\times Y)$ characterizes all the trajectories that are
described by a differential equation

$\displaystyle\dot{x}(t)$ $\displaystyle=f(x,w)+g(x,w)u$ (1)

and $o\mathrel{\mathop{\ordinarycolon}}x\mapsto y$ is the output function.

To guarantee the existence and uniqueness of the system trajectory, we
conveniently assume that the vector field and the output map are locally
Lipschitz. Now we formally define an interconnected system.

######  Definition 2.

Given $N$ subsystems $\\{G_{i}\\}_{i\in\mathcal{I}}$,
$G_{i}=(U_{i},W_{i},X_{i},Y_{i},X^{0}_{i},\mathcal{T}_{i}),\mathcal{I}=\\{1,2,\ldots,N\\},$
and a binary connectivity relation
$\mathcal{E}\subseteq\mathcal{I}\times\mathcal{I}$, we say
$\\{G_{i}\\}_{i\in\mathcal{I}}$ is _compatible for composition_ with respect
to $\mathcal{E}$ if $\Pi_{j\in N(i)}Y_{j}\subseteq W_{i}$, where
$N(i)=\\{j\mathrel{\mathop{\ordinarycolon}}(j,i)\in\mathcal{E}\\}$ is the
index set of subsystems that influence $G_{i}$. $G_{j}$ ($G_{i}$) is referred
to as a parent (child) node of $G_{i}$ ($G_{j}$).

In this definition, a set of subsystems is compatible for composition when,
for each subsystem, the output space of all parental subsystems is a subset of
its internal input space.

When the subsystems $\\{G_{i}\\}_{i\in\mathcal{I}}$ are compatible for
composition w.r.t. $\mathcal{E}$, the composed system, also referred to as the
interconnected system, is denoted by
$\langle(G_{i})_{i\in\mathcal{I}},\mathcal{E}\rangle=(U,\\{0\\},X,Y,X^{0},\mathcal{T})$,
where
$U=\Pi_{i\in\mathcal{I}}U_{i},X=\Pi_{i\in\mathcal{I}}X_{i},Y=\Pi_{i\in\mathcal{I}}Y_{i},X^{0}=\Pi_{i\in\mathcal{I}}X_{i}^{0}$.
Denote the composed state by $x$, the composed external input $u$, and the
composed output $y$. Then $(u(t),0,x(t),y(t))\in\mathcal{T}$ if and only if
for all $i\in\mathcal{I}$, there exists
$(u_{i}(t),w_{i}(t),x_{i}(t),y_{i}(t))\in\mathcal{T}_{i}$ and
$w_{i}(t)=(y_{j_{1}}(t),y_{j_{2}}(t),\ldots,y_{j_{p}}(t)),$ where
$N(i)=\\{j_{1},j_{2},\ldots,j_{p}\\}$.

###  II-B Assume-guarantee contracts for invariance

To begin with, we introduce notation that will help us define the set of all
continuous trajectories that always stay in a set. Let a nonempty set
$S\subseteq\mathbb{R}^{n}$. Define

$I_{S}^{E}=\\{z\mathrel{\mathop{\ordinarycolon}}E\to\mathbb{R}^{n}\in
M(\mathbb{R}^{n})\mathrel{\mathop{\ordinarycolon}}\forall t\in E,z(t)\in
S\\},$ (2)

where $E$ is a time interval. In the following, the superscript $E$ is
neglected as it is usually chosen as the maximal time interval of the
existence of solutions to the continuous-time system. An invariance assume-
guarantee contract (iAGC) is defined as follows:

######  Definition 3.

For a continuous-time system $G=(U,W,X,Y,X^{0},\mathcal{T})$, an _invariance
assume-guarantee contract_ (iAGC) for $G$ is a tuple
$C=(I_{\underline{W}},I_{\underline{X}},I_{\underline{Y}})$ where
$\underline{W}\subseteq W,\underline{X}\subseteq X,\underline{Y}\subseteq Y$.
We refer to $I_{\underline{W}}$ as the set of assumptions on the internal
inputs, and $I_{\underline{X}},I_{\underline{Y}}$ as the sets of guarantees on
the states and outputs. We say a system $G$ _satisfies a contract_
$C=(I_{\underline{W}},I_{\underline{X}},I_{\underline{Y}})$, denoted $G\models
C$, if there exists a feedback control
$k(\cdot,\cdot)\mathrel{\mathop{\ordinarycolon}}X\times W\to U$ such that for
all $t>0$, for all $w|_{[0,t]}\in I_{\underline{W}}$, the state and output
fulfill $x|_{[0,t]}\in I_{\underline{X}},y|_{[0,t]}\in I_{\underline{Y}}$ for
all trajectories $(u(t)=k(x,w),w(t),x(t),y(t))\in\mathcal{T}$.

A key result that establishes the compositional reasoning of the system
property is the following:

######  Lemma 1 (Compositional reasoning).

Consider an interconnected system
$\langle(G_{i})_{i\in\mathcal{I}},\mathcal{E}\rangle=(U,\\{0\\},X,Y,X^{0},\mathcal{T})$
composed of $N$ subsystems with a compatible binary connectivity relation
$\mathcal{E}$. If for each subsystem
$G_{i}=(U_{i},W_{i},X_{i},Y_{i},X_{i}^{0},\mathcal{T}_{i})$, there exists an
invariance assume-guarantee contract
$C_{i}=(I_{\underline{W}_{i}},I_{\underline{X}_{i}},I_{\underline{Y}_{i}})$
such that $G_{i}\models C_{i}$ and $\Pi_{j\in
N(i)}{I_{\underline{Y}_{j}}}\subseteq I_{\underline{W}_{i}},$ then
$\langle(G_{i})_{i\in\mathcal{I}},\mathcal{E}\rangle\models C$ with
$C=(\\{0\\},\Pi_{i\in\mathcal{I}}{I_{\underline{X}_{i}}},\Pi_{i\in\mathcal{I}}{I_{\underline{Y}_{i}}}).$

This Lemma is a special case of [16, Theorem 3] and helps circumvent possible
cyclic reasoning when deciding the forward invariance of interconnected
systems.

###  II-C System safety and barrier functions

Now we give the formal definition of safety of a continuous-time system.
Throughout this work, we refer to a _safe region_ as the collection of states
that are benign, for example, unoccupied configuration space in robotic
applications, and a _safe set_ as a subset of the safety region that is also
forward invariant.

######  Definition 4.

Given a system $G=(U,W,X,Y,X^{0},\mathcal{T})$ and a safe region
$\mathcal{Q}\subseteq X$, we say the system $G$ is _safe with respect to an
internal input set $\underline{W}$_ if and only if
$X^{0}\subseteq\mathcal{Q}$, and for all initial states $x_{0}\in X^{0}$ and
for all internal input signals $w|_{[0,t]}\in I_{\underline{W}},$ there exists
an external input signal $u|_{[0,t]}\in I_{U}$ such that $x|_{[0,t]}\in
I_{\mathcal{Q}}$ for all $t>0$.

We note that the internal input $w$ is assumed to be known via communication.
This is in contrast to the definition of a robust control invariant set where
$w$ is treated as disturbance and unknown [22], in which case the qualifier
$\exists u$ proceeds $\forall w$. One way to certify the safety of the system
is to find a (control) barrier function, also known as a barrier certificate
[5], which is given by

######  Definition 5.

A differential function $h\mathrel{\mathop{\ordinarycolon}}X\to\mathbb{R}$ for
system $G=(U,W,X,Y,X^{0},\mathcal{T})$ is called a _control barrier function_
with respect to an internal input set $\underline{W}\subseteq W$, if there
exists a class $\mathcal{K}$ function $\alpha$ such that $\forall
w\in\underline{W},\exists u\in U$ the following inequality holds for all $x$

$\nabla h(x)f(x,w)+\nabla h(x)g(x,w)u+\alpha(h(x))\geq 0.$ (3)

When the external input set is empty, i.e., $U=\emptyset$, $h$ is called a
_barrier function_ as this system has no active control. When such a (control)
barrier function is found, then the set
$\mathcal{C}=\\{x\mathrel{\mathop{\ordinarycolon}}h(x)\geq 0\\}$ is
(controlled) forward invariant, and asymptotically stable when it is
compact[2]. If $X^{0}\subseteq\mathcal{C}\subseteq\mathcal{Q}$, then system
safety is certified[2, 7]. In general, finding an invariant set $\mathcal{C}$
is computationally expensive for large nonlinear systems.

###  II-D Sum-of-squares programs

One tractable approach to deal with infinite inequalities as in (3) is via
sum-of-squares programming. A standard sum-of-squares (SOS) program takes the
following form

$\displaystyle\min_{p_{1},\ldots,p_{k}}\textstyle\sum_{j=1}^{k}a_{j}p_{j}$ (4)
$\displaystyle\text{s.t.
}b_{0}(x)+\textstyle\sum_{j=1}^{k}p_{i}b_{j}(x)\in\Sigma[x],$

where the decision variables $p_{1},\ldots,p_{k}\in\mathbb{R}$, constants
$a_{1},\ldots,a_{k}\in\mathbb{R}$ are the weights, and
$b_{0},\ldots,b_{k}\in\mathcal{R}[x]$ are given polynomials. This SOS program
is a convex optimization problem and can be equivalently transformed into a
semi-definite program (SDP). Interested readers are referred to [23, 24] for
more details.

###  II-E Problem formulation

In this work, we aim to numerically verify the safety property of
interconnected systems. The following sub-problems are considered:

  1. (P1)

For a continuous-time subsystem
$G_{i}=(U_{i},W_{i},X_{i},Y_{i},X^{0}_{i},\mathcal{T}_{i})$ and a given safe
region $\mathcal{Q}_{i}\subseteq X_{i}$, construct an invariance assume-
guarantee contract
$C_{i}=(I_{\underline{W}_{i}},I_{\underline{X}_{i}},I_{\underline{Y}_{i}})$
such that $G_{i}\models C_{i}$ and
$X_{i}^{0}\subseteq\underline{X}_{i}\subseteq\mathcal{Q}_{i}$;

  2. (P2)

For an interconnected system
$\langle(G_{i})_{i\in\mathcal{I}},\mathcal{E}\rangle=(U,\\{0\\},X,Y,X^{0},\mathcal{T})$
and a safe region $\mathcal{Q}=\Pi_{i\in\mathcal{I}}\mathcal{Q}_{i},$
$\mathcal{Q}_{i}\subseteq X_{i}$, construct an invariance contract
$C=(\\{0\\},I_{\underline{X}},I_{\underline{Y}})$ such that $G\models C$ and
$X^{0}\subseteq\underline{X}\subseteq\mathcal{Q}$.

If such a $\underline{X}$ is found, then safety of the interconnected system
is certified.

######  Assumption 1.

We assume the following:

  1. 1.

the local feedback law $u_{i}=k_{i}(x_{i},w_{i})\in U_{i}$ is known, but it
does not necessarily render the interconnected system safe;

  2. 2.

The class $\mathcal{K}$ function $\alpha(\cdot)$ in (3) is chosen to be a
linear function with constant gain $a$.

  3. 3.

The initial set $X_{i}^{0}$, safe region $\mathcal{Q}_{i}$, and the internal
input set $W_{i}$ are super-level sets of (possibly vector-valued)
differentiable functions, i.e.,
$X_{i}^{0}=\\{x_{i}\mathrel{\mathop{\ordinarycolon}}b_{i}^{0}(x_{i})\geq
0\\},\mathcal{Q}_{i}=\\{x_{i}\mathrel{\mathop{\ordinarycolon}}q_{i}(x_{i})\geq
0\\},W_{i}=\\{(y_{j_{1}},y_{j_{2}},\ldots,y_{j_{p}})\mathrel{\mathop{\ordinarycolon}}d^{i}_{j_{k}}(y_{j_{k}})\geq
0,k=1,2,\ldots,p\\},\text{where }N(i)=\\{j_{1},j_{2},\ldots,j_{p}\\}$.

  4. 4.

$b_{i}^{0},q_{i}\in\mathcal{R}[x_{i}],d^{i}_{j_{k}}(y_{j_{k}})\in\mathcal{R}[y_{j_{k}}],f_{i},g_{i},k_{i}\in\mathcal{R}[x_{i},w_{i}]$
are polynomials.

  5. 5.

The subsets of $W_{i},\mathcal{Q}_{i}$, i.e.,
$\underline{W}_{i},\underline{\mathcal{Q}}_{i}$ are chosen in the form of

$\displaystyle\underline{\mathcal{Q}}_{i}=\\{x_{i}\mathrel{\mathop{\ordinarycolon}}q_{i}(x_{i})\geq\zeta\bm{1}\text{
for some }\zeta\geq 0\\},$
$\displaystyle\underline{W}_{i}=\\{(y_{j_{1}},\ldots,y_{j_{p}})\mathrel{\mathop{\ordinarycolon}}d_{j_{k}}^{i}(y_{j_{k}})\geq\delta\bm{1}\text{
for some }\delta\geq 0\\}.$

  6. 6.

When searching for non-negative polynomials, we restrict the search to the set
of SOS polynomials up to a certain degree.

These restrictions, even though conservative, will facilitate the convergence
and completeness analysis that will be clear later. We note that the first two
assumptions are in place to avoid bilinear terms when constructing the SOS
programs. Both can be relaxed by considering iterative optimization
approaches. See [7] for more details. Assumptions 1.3 and 1.5 help to
parameterize the assumption and the guarantee sets by scalar variables
$\delta$ and $\zeta$, respectively. Assumptions 1.4 and 1.6 are standard in
the field of SOS-based system verification.

For notational simplicity, we define the set projection of an internal input
set $W_{i}$ of subsystem $G_{i}$ with respect to subsystem $G_{k}$ as
$\text{Proj}_{k}(W_{i})=\\{y_{k}\mathrel{\mathop{\ordinarycolon}}d^{i}_{k}(y_{k})\geq
0\\}$ if $k\in N(i)$, and $\text{Proj}_{k}(W_{i})=\emptyset$ otherwise. For a
mapping $o\mathrel{\mathop{\ordinarycolon}}X\to Y$, let
$o^{-1}(\underline{Y})\mathrel{\mathop{\ordinarycolon}}=\\{x\mathrel{\mathop{\ordinarycolon}}o(x)\in\underline{Y}\\}$.

##  III Proposed solutions

The proposed approach consists of 1) numerically constructing iAGCs for
subsystems by synthesizing local (control) barrier functions, and 2)
negotiating iAGCs among subsystems to certify the safety property of the
interconnected system. We also discuss the convergence properties of our
approach.

###  III-A Local barrier function and AGC construction

In this subsection, we will focus on tackling Problem (P1) for a subsystem
$G_{i}=(U_{i},W_{i},X_{i},Y_{i},X^{0}_{i},\mathcal{T}_{i})$. Under Assumption
1, the closed-loop subsystem dynamics of (1) are

$\dot{x}_{i}(t)=f_{i}(x_{i},w_{i})+g_{i}(x_{i},w_{i})k_{i}(x_{i},w_{i})\mathrel{\mathop{\ordinarycolon}}=F_{i}(x_{i},w_{i}).$
(5)

In this subsection, for the sake of notation simplicity, we will drop the
subscript $i$ when no confusion arises.

First we show the relations between a) finding a control barrier function, b)
constructing an invariance assume-guarantee contract, and c) establishing the
safety property of a subsystem.

######  Proposition 1.

Consider a continuous-time system $G=(U,W,X,Y,X^{0},\mathcal{T})$, a safe
region $\mathcal{Q}$ and an internal input set $\underline{W}\subseteq W$.
Consider the following claims:

  1. 1

there exists a CBF $h$ with respect to the internal input set $\underline{W}$.
Denote by $\mathcal{C}=\\{x\mathrel{\mathop{\ordinarycolon}}h(x)\geq 0\\}$;

  2. 2

the system $G\models C$, where
$C=(I_{\underline{W}},I_{\mathcal{C}},I_{o(\mathcal{C})})$;

  3. 3

$X^{0}\subseteq\mathcal{C}\subseteq\mathcal{Q}$;

  4. 4

the system is safe with respect to $\underline{W}$;

We have 1 $\implies$ 2; 2 and 3 $\implies$ 4.

###### Proof.

1

$\implies$ 2: 1 implies that the set $\mathcal{C}$ is forward invariant when
the signal $w(t)\in I_{\underline{W}}$, which implies 2 from Definition 3. 2
and 3 $\implies$ 4: This is straightforward according to Definition 4. ∎

Numerically, one can formulate the conditions of 1 and 3 of Proposition 1 as a
set of SOS constraints, as follows.

######  Proposition 2.

Consider a continuous-time system $G=(U,W,X,Y,X^{0},\mathcal{T})$ and a safe
region $\mathcal{Q}$. If there exist SOS polynomials
$\sigma_{init},\sigma_{safe}\in\Sigma[x]$,
$\sigma_{k}\in\Sigma[x,y_{k}],k=1,2,\ldots,p$, polynomial
$h\in\mathcal{R}(x)$, and positive $\epsilon,a,\delta$ such that

$\displaystyle h(x)-\sigma_{init}b^{0}(x)$ $\displaystyle\in\Sigma[x];$ (6a)
$\displaystyle-h(x)+\sigma_{safe}q(x)$ $\displaystyle\in\Sigma[x];$ (6b)
$\displaystyle\nabla h(x)F(x,y_{1},\ldots,y_{p})+ah(x)\hskip 8.5359pt$
$\displaystyle-\sum_{k=1}^{p}\sigma_{k}(d_{k}(y_{k})-\delta)-\epsilon$
$\displaystyle\in\Sigma[x,y_{1},\ldots,y_{p}].$ (6c)

then, letting
$\underline{W}=\\{(y_{1},\ldots,y_{k}\ldots,y_{p})\mathrel{\mathop{\ordinarycolon}}d_{k}(y_{k})\geq\delta\\}$,
1, 2, 3, 4 in Proposition 1 hold.

###### Proof.

(6c) is a SOS polynomial and thus $\nabla h(x)F(x,w)+ah(x)\geq 0,\forall
w\in\underline{W}$. This shows that claim 1 holds. Based on non-negativeness
of SOS polynomials, (6a) and (6b) imply that $\forall x,b^{0}(x)\geq 0\implies
h(x)\geq 0$, and $\forall x,h(x)\geq 0\implies q(x)\geq 0$, respectively. That
is, $X_{i}^{0}\subseteq\mathcal{C}\subseteq\mathcal{Q}$. Thus 3 holds.
Following Proposition 1, we conclude the proof. ∎

Even though (6) is only a sufficient condition for system safety, it is a
condition we can verify numerically (and efficiently when the system size is
small). For this reason, we say that $G$ is _certified to be safe_ in
$\mathcal{Q}$ w.r.t. $\underline{W}$ if condition (6) holds. We introduce the
following special sets that are useful for contract composition later. In what
follows, we take $0<\epsilon<<1$ and $a$ in (6) to be positive constants.

####  III-A1 Maximal internal input set

To quantify the largest internal input set a subsystem can tolerate while
still remaining safe, we propose the following optimization problem:

$\displaystyle\min\delta$ (7) s.t.
$\displaystyle\eqref{eq:sos_initialset},\eqref{eq:sos_safeset},\eqref{eq:sos_cbf_condiiton},\delta\geq
0,$

where the decision variables include SOS polynomials
$\sigma_{init},\sigma_{safe}\in\Sigma[x]$,
$\sigma_{k}\in\Sigma[y_{k}],k=1,2,\ldots,p$, polynomials $h\in\mathcal{R}(x)$,
and a scalar $\delta$. It should be noted that although (7) contains a
bilinear term $\sigma_{input}\delta$, this can be solved efficiently by
bisection as $\delta$ is a scalar. If (7) is feasible, denote the optimal
value by $\delta^{\star}$ and the corresponding internal input set
$\underline{W}^{\star}$. We call $\underline{W}^{\star}$ the _maximal internal
input set_ for a given subsystem and safe region $\mathcal{Q}$.

####  III-A2 Minimal safe region

Given an internal input set $\underline{W}$, to quantify the least impact of a
subsystem on its child subsystem, we propose the following optimization
problem:

$\displaystyle\max\zeta$ (8) s.t.
$\displaystyle\eqref{eq:sos_initialset},\eqref{eq:sos_cbf_condiiton},\zeta\geq
0$ $\displaystyle-h(x)+\sigma_{safe}(q(x)-\zeta)\in\Sigma[x];$

where the decision variables include SOS polynomials
$\sigma_{init},\sigma_{safe}\in\Sigma[x]$,
$\sigma_{k}\in\Sigma[y_{k}],k=1,2,\ldots,p$, polynomials $h\in\mathcal{R}(x)$,
and a scalar $\zeta$. We take $\epsilon,a$ to be positive constants. $\delta$
in (6c) is known as we assume $\underline{W}$ is given. It should be noted
that although (8) contains a bilinear term $\sigma_{safe}\zeta$, this can be
solved efficiently by bisection as $\zeta$ is a scalar. If (8) is feasible,
denote the optimal value by $\zeta^{\star}$ and the corresponding safe region
$\underline{\mathcal{Q}}^{\star}$. We call $\underline{\mathcal{Q}}^{\star}$
the _minimal safe region_ for a given $\underline{W}$.

We have the following properties about the maximal internal input set
$\underline{W}^{\star}$ and the corresponding minimal safe region
$\underline{\mathcal{Q}}^{\star}$.

######  Proposition 3.

Under Assumption 1, for a continuous-time system
$G=(U,W,X,Y,X^{0},\mathcal{T})$ and a safe region $\mathcal{Q}$, the following
results hold:

  1. 1.

If (7) is feasible for some $\delta^{\prime}\geq 0$, then (7) is also feasible
for $\delta^{\prime\prime},\delta^{\prime\prime}\geq\delta^{\prime}$. If (8)
is feasible for some $\zeta^{\prime}>0$, then (8) is also feasible for
$\zeta^{\prime\prime},0\leq\zeta^{\prime\prime}\leq\zeta^{\prime}$.

  2. 2.

Consider two safe regions
$\underline{\mathcal{Q}}^{\prime}\subseteq\underline{\mathcal{Q}}^{\prime\prime}\subseteq\mathcal{Q}$.
If (7) is feasible for the safe region $\underline{\mathcal{Q}}^{\prime}$,
then (7) is also feasible for $\underline{\mathcal{Q}}^{\prime\prime}$.
Denoting the respective optimal values by
$\delta^{\prime},\delta^{\prime\prime}$ and the corresponding internal input
sets
$\underline{\mathcal{W}}^{\prime},\underline{\mathcal{W}}^{\prime\prime}$,
then $\delta^{\prime\prime}\leq\delta^{\prime}$ and
$\underline{\mathcal{W}}^{\prime}\subseteq\underline{\mathcal{W}}^{\prime\prime}\subseteq\mathcal{W}$.

  3. 3.

Consider two internal input sets
$\underline{\mathcal{W}}^{\prime}\subseteq\underline{\mathcal{W}}^{\prime\prime}\subseteq\mathcal{W}$.
If (8) is feasible with the internal input set
$\underline{\mathcal{W}}^{\prime\prime}$, then (8) is also feasible for
$\underline{\mathcal{W}}^{\prime}$. Denoting the respective optimal values by
$\zeta^{\prime},\zeta^{\prime\prime}$ and the corresponding safe regions
$\underline{\mathcal{Q}}^{\prime},\underline{\mathcal{Q}}^{\prime\prime}$,
then $0\leq\zeta^{\prime\prime}\leq\zeta^{\prime}$ and
$\underline{\mathcal{Q}}^{\prime}\subseteq\underline{\mathcal{Q}}^{\prime\prime}\subseteq\mathcal{Q}$.

  4. 4.

If (7) is feasible, then $\underline{W}^{\star}$ is the largest internal input
set w.r.t. which $G$ is certified to be safe; if infeasible, then there exists
no $\underline{W}\subseteq W$ w.r.t. which $G$ can be certified to be safe.

  5. 5.

If (7) is feasible, letting $\underline{W}=\underline{W}^{\star}$, then (8) is
feasible and $\underline{\mathcal{Q}}^{\star}$ is the smallest safe region in
which $G$ is safe w.r.t. $\underline{W}^{\star}$.

###### Proof.

Claim 1) holds since one verifies that, when
$\delta^{\prime\prime}\geq\delta^{\prime}$, any feasible decision variables
$\sigma_{init},\sigma_{safe},\sigma_{k},k=1,2,\ldots,p,h$ to (7) for
$\delta^{\prime}$ are also feasible for the case of $\delta^{\prime\prime}$.
Similarly, when $0\leq\zeta^{\prime\prime}\leq\zeta^{\prime}$, any feasible
decision variables to (8) for $\zeta^{\prime}$ are also feasible for the case
of $\zeta^{\prime\prime}$.

For
$\underline{\mathcal{Q}}^{\prime}\subseteq\underline{\mathcal{Q}}^{\prime\prime}$
(the corresponding $\zeta^{\prime}\geq\zeta^{\prime\prime}\geq 0$), if (7) is
feasible for $\zeta^{\prime}$, then the feasible decision variables are also
feasible for the case of $\zeta^{\prime\prime}$. As (7) minimizes over
$\delta$ and the feasible set of the case $\zeta^{\prime}$ is a subset of that
of $\zeta^{\prime\prime}$, then $\delta^{\prime\prime}\leq\delta^{\prime}$,
proving Claim 2).

Similar argument applies for Claim 3). For
$\underline{\mathcal{W}}^{\prime}\subseteq\underline{\mathcal{W}}^{\prime\prime}\subseteq\mathcal{W}$
(with the corresponding $\delta^{\prime}\geq\delta^{\prime\prime}\geq 0$), if
(8) is feasible for $\delta^{\prime\prime}$, then the feasible decision
variables are also feasible for the case of $\delta^{\prime}$. As (8)
maximizes over $\zeta$ and the feasible set of the case
$\delta^{\prime\prime}$ is a subset of that of $\delta^{\prime}$, then
$\zeta^{\prime}\geq\zeta^{\prime\prime}$, proving Claim 3).

Claim 4) is true since
$\underline{W}^{\prime\prime}\subseteq\underline{W}^{\prime}$ if and only if
the corresponding $\delta^{\prime\prime}\geq\delta^{\prime}$, and the program
in (7) minimizes $\delta$. We note that the condition $\delta\geq 0$ comes
from the condition $\underline{W}\subseteq W$. If (7) is feasible, then the
feasible decision variables for $\delta=\delta^{\star}$ are also feasible for
(8) when $\zeta=0$. Further noting that
$\underline{Q}^{\prime\prime}\subseteq\underline{Q}^{\prime}$ if and only if
the corresponding $\zeta^{\prime\prime}\geq\zeta^{\prime}$ and (8) maximizes
$\zeta$, thus Claim 5) is deduced. ∎

Proposition 3’s items 2 and 3 show a monotonic relation between the internal
input sets and the safe regions. Intuitively, with a larger safe region, the
system can tolerate a larger disturbance (internal input set); with a larger
disturbance (internal input set), the most confined safe region will become
larger. Proposition 3’s items 4 and 5 further state that, for a given safe
region, $\underline{W}^{\star}$ is the largest internal input set that a
system can bear while remaining safe; for a given internal input set,
$\underline{\mathcal{Q}}^{\star}$ is the most confined influence a system has
for its child subsystems. When (7) and (8) are feasible for a subsystem
$G_{i}$, denoting the corresponding sets as
$\underline{W}_{i}^{\star},\underline{\mathcal{Q}}_{i}^{\star}$, then we
construct a local contract
$C_{i}=(I_{\underline{W}_{i}},I_{\underline{X}_{i}},I_{\underline{Y}_{i}})$,
where
$\underline{W}_{i}=\underline{W}_{i}^{\star},\underline{X}_{i}=\underline{\mathcal{Q}}_{i}^{\star}$,
and $\underline{Y}_{i}=o_{i}(\underline{X}_{i})$.

###  III-B Contract composition and negotiation

In this section, we consider the interconnected system
$G=\langle(G_{i})_{i\in\mathcal{I}},\mathcal{E}\rangle$,
$G_{i}=(U_{i},W_{i},X_{i},Y_{i},X^{0}_{i},\mathcal{T}_{i})$ with safe region
$\mathcal{Q}_{i}\subseteq X_{i}$. We have the following results on the safety
properties of the interconnected system.

######  Proposition 4.

If, for each subsystem $G_{i}$, an iAGC
$C_{i}=(I_{\underline{W}_{i}},I_{\underline{X}_{i}},I_{\underline{Y}_{i}})$
exists such that $X_{i}^{0}\subseteq\underline{X}_{i}\subseteq\mathcal{Q}_{i}$
and

$\Pi_{j\in N(i)}\underline{Y}_{j}\subseteq\underline{W}_{i},$ (9)

then the interconnected system
$\langle(G_{i})_{i\in\mathcal{I}},\mathcal{E}\rangle$ is safe.

###### Proof.

From Lemma 1, $\langle(G_{i})_{i\in\mathcal{I}},\mathcal{E}\rangle\models C$
with
$C=(\\{0\\},\Pi_{i\in\mathcal{I}}{I_{\underline{X}_{i}}},\Pi_{i\in\mathcal{I}}{I_{\underline{Y}_{i}}}).$
Note that
$\Pi_{i\in\mathcal{I}}{\underline{X}_{i}^{0}}\subseteq\Pi_{i\in\mathcal{I}}{\underline{X}_{i}}\subseteq\Pi_{i\in\mathcal{I}}{\mathcal{Q}_{i}}$,
thus $\langle(G_{i})_{i\in\mathcal{I}},\mathcal{E}\rangle$ is safe. ∎

We refer to the condition (9) as the _contract compatibility condition_ as it
indicates whether the contract of a subsystem agrees with that of its parent
subsystems. In the general case, the contracts $C_{i},i\in\mathcal{I}$ found
locally may not satisfy this condition, and we have to refine them so that (9)
holds. We call this refinement process _negotiation_. In what follows, we
consider three cases and propose several different algorithms. We note that
all algorithms are sound, but differ in finite-step termination and
completeness guarantees.

####  III-B1 Acyclic connectivity graph

In this case, we assume that there exists no cycle in the connectivity graph
$(\mathcal{I},\mathcal{E})$. In this case, the hierarchical tree structure
resembles a client-contractor relation model. For $k\in\text{Child}(i)$, we
could view $G_{k}$ as a client with an iAGC
$(I_{\underline{W}_{k}},I_{\underline{X}_{k}},I_{\underline{Y}_{k}})$, who
gives specifications on the behaviour of its parent node $G_{i}$ (viewed as
contractors) by $\underline{W}_{k}$. Based on this interpretation, we propose
Algorithm 1.

In Algorithm 1, $\mathcal{I}_{0},\mathcal{I}_{1},\mathcal{I}_{-1}$ represent
the index sets of ready-to-update, to-be-updated, and updated subsystems,
respectively. The algorithm starts with the local contract construction for
the leaf nodes. Following a bottom-up traversal along the connectivity graph,
for each subsystem $G_{i}$ in $\mathcal{I}_{0}$, Algorithm 1 first updates its
safe region $\mathcal{Q}_{i}$ such that it agrees with all its child nodes.
This is explicitly conducted in Algorithm 2, while no operation is needed for
leaf nodes. The set intersection in Algorithm 2 is again cast as a SOS
program, as follows:

$\displaystyle\hskip 56.9055pt\min_{\zeta\geq 0}\zeta$ (10) s.t.
$\displaystyle q_{i}(x_{i})-\zeta-\sigma_{k}(d^{k}_{i}\circ
o_{i}(x_{i})-\delta^{k})$ $\displaystyle\hskip
65.44133pt\in\Sigma[x_{i}],\forall k\in\text{Child}(i),$

where the decision variables include
$\sigma_{k}\in\Sigma[x_{i}],k\in\text{Child}(i),$ and a scalar $\zeta$. Recall
here $o_{i}$ is the output map of subsystem $G_{i}$,
$\text{Proj}_{i}(\underline{W}_{k})=\\{y_{i}\mathrel{\mathop{\ordinarycolon}}d^{k}_{i}(y_{i})\geq\delta^{k}\\}$.
Denoting the optimal value by $\zeta^{\prime}$ and
$\mathcal{Q}_{i}^{\prime}=\\{x\mathrel{\mathop{\ordinarycolon}}q_{i}(x_{i})\geq\zeta^{\prime}\\}$,
$\mathcal{Q}_{i}^{\prime}$ is then the largest inner-approximation of
$\bigcap_{k\in\text{Child}(i)}o_{i}^{-1}(\text{Proj}_{i}(\underline{W}_{k}))\cap\mathcal{Q}_{i}$.
Recall that the subset of $\mathcal{Q}_{i}$ is parameterized by $\zeta$ from
Assumption 1.5.

After updating the safe region, Algorithm 1 calculates the maximal internal
input set $\underline{W}_{i}^{\star}$ (Line 6), which can be seen as the least
requirement on its parent nodes as discussed in Proposition 3. Algorithm 3
then moves $G_{i}$ to $\mathcal{I}_{-1}$, and checks for every to-be-updated
subsystems whether all their child subsystems have been updated. If yes, then
that subsystem is moved to the set of ready-to-update subsystems
$\mathcal{I}_{0}$ and will be updated accordingly.

######  Proposition 5.

Consider an interconnected system with an acyclic connectivity graph
$(\mathcal{I},\mathcal{E})$. Algorithm 1 has the following properties:

  1. 1.

Algorithm 1 terminates in finite steps and returns either True or False.

  2. 2.

If Algorithm 1 returns True, then iAGCs
$C_{i}=(I_{\underline{W}_{i}},I_{\underline{X}_{i}},I_{\underline{Y}_{i}}),i\in\mathcal{I}$
satisfy the conditions in Proposition 4.

  3. 3.

If Algorithm 1 returns False, then there exist no iAGCs
$C_{i}=(I_{\underline{W}_{i}},I_{\underline{X}_{i}},I_{\underline{Y}_{i}}),i\in\mathcal{I}$
that satisfy the conditions in Proposition 4 under Assumption 1.

###### Proof.

For every iteration of Algorithm 1, it will either terminate due to
infeasibility or increase the cardinality of $\mathcal{I}_{-1}$ by $1$. Since
$|\mathcal{I}_{-1}|$ is upper bounded by $|\mathcal{I}|$, we know it has to
terminate in finite steps. This proves Claim 1).

When Algorithm 1 returns True, each subsystem has gone through Step 4 - Step
10 and computed an iAGC
$(I_{\underline{W}_{i}^{\star}},I_{\underline{X}_{i}},I_{\underline{Y}_{i}})$
in a bottom-up transverse. As Step 4 reduces the safe region of subsystem
$G_{i}$, and from (7), we have
$X_{i}^{0}\subseteq\underline{X}_{i}\subseteq\mathcal{Q}_{i}^{\prime}\subseteq\mathcal{Q}_{i}$.
Moreover, for any $k\in\text{Child}(i)$, based on Algorithm 2,
$\underline{Y}_{i}=o_{i}(\underline{X}_{i})\subseteq
o_{i}(\mathcal{Q}_{i}^{\prime})\subseteq\text{Proj}_{i}(\underline{W}_{k})$,
which proves that the contract compatibility condition (9) holds. Following
Proposition 4, we thus certify the safety of the interconnected systems. This
proves Claim 2).

We show claim 3) by contradiction. Suppose that there exist iAGCs
$C_{i}=(I_{\underline{W}_{i}},I_{\underline{X}_{i}},I_{\underline{Y}_{i}}),i\in\mathcal{I}$
that satisfy the conditions in Proposition 4. Denote the safe regions with
which such iAGCs are obtained by $\mathcal{Q}_{i}^{\prime},i\in\mathcal{I}$,
i.e., the functions defining those sets fulfill the SOS constraints in (8)
(but not necessarily minimize the size of the safe region) and the sets
fulfill the compatibility condition (9). Thanks to the tree structure, we can
start our argumentation from the leaf nodes and iteratively reason about the
nodes that are one level above, and end at the root node. From Proposition 3,
for $G_{i}$ being the leaf nodes, we have
$\underline{W}_{i}\subseteq\underline{W}_{i}^{\star}$, the set obtained from
(7). For $G_{j}$ being the nodes one level above the leaf nodes, from
Algorithm 2, we know
$\mathcal{Q}_{j}^{\prime}\subseteq\mathcal{Q}_{j}^{\prime\star}$, for that
$\mathcal{Q}_{j}^{\prime\star}$ is the largest inner-approximation of all safe
regions and that $\underline{W}_{i}\subseteq\underline{W}_{i}^{\star}$.
Following Proposition 3 item 2, we know
$\underline{W}_{j}\subseteq\underline{W}_{j}^{\star}$, where
$\underline{W}_{j}^{\star}$ is the set obtained by solving (7) with
$\mathcal{Q}_{j}^{\prime\star}$. Recursively, we thus obtain that
$(I_{\underline{W}_{i}^{\star}},I_{\underline{X}_{i}},I_{\underline{Y}_{i}}),i\in\mathcal{I}$
exists. This contradicts with the premise that Algorithm 1 returns False. Thus
Claim 3) is proven. ∎

Algorithm 1 Contract construction for acyclic graph

1:$G_{i}$, $\mathcal{Q}_{i},\forall i\in\mathcal{I}$

2:$\mathcal{I}_{0}\leftarrow$ set of leaf nodes,
$\mathcal{I}_{1}\leftarrow\mathcal{I}\setminus\mathcal{I}_{0}$,
$\mathcal{I}_{-1}\leftarrow\emptyset$.

3:while $\mathcal{I}_{0}\neq\emptyset$ do

4: for each subsystem $G_{i},i\in\mathcal{I}_{0}$ do

5: $\mathcal{Q}_{i}^{\prime}\leftarrow$ update the local safe region
$\mathcal{Q}_{i}$ by Alg. 2;

6: try

7: calculate $\delta_{i}^{\star}$ by solving (7) with safe region
$\mathcal{Q}_{i}^{\prime}$;

8: compute the corresp. iAGC
$(I_{\underline{W}_{i}^{\star}},I_{\underline{X}_{i}},I_{\underline{Y}_{i}})$

9: catch infeasible

10: return False;

11: end try

12: update $\mathcal{I}_{0},\mathcal{I}_{1},\mathcal{I}_{-1}$ by Alg. 3.

13: end for

14:end while

15:return True.

Algorithm 2 Update safe region

1:Safe region $\mathcal{Q}_{i}$ and iAGCs
$(I_{\underline{W}_{k}},I_{\underline{X}_{k}},I_{\underline{Y}_{k}})$ for all
$k\in\text{Child}(i)$.

2:$M_{i}\leftarrow\bigcap_{k\in\text{Child}(i)}o_{i}^{-1}(\text{Proj}_{i}(\underline{W}_{k}))\cap\mathcal{Q}_{i}$

3:$\mathcal{Q}_{i}^{\prime}\leftarrow$ largest inner-approximation of $M_{i}$
by (10)

4:return $\mathcal{Q}_{i}^{\prime}$.

Algorithm 3 Update $\mathcal{I}_{0},\mathcal{I}_{1}$ and $\mathcal{I}_{-1}$

1:Subsystem $G_{i}$, $\mathcal{I}_{0},\mathcal{I}_{1}$ and $\mathcal{I}_{-1}$.

2:$\mathcal{I}_{0}\leftarrow\mathcal{I}_{0}\setminus\\{i\\},\mathcal{I}_{-1}\leftarrow\mathcal{I}_{-1}\cup\\{i\\}$,

3:for each subsystem $G_{k}$, $k\in\mathcal{I}_{1}$ do

4: if $\text{Child}(k)\subseteq\mathcal{I}_{-1}$ then

5:
$\mathcal{I}_{0}\leftarrow\mathcal{I}_{0}\cup\\{k\\},\mathcal{I}_{1}\leftarrow\mathcal{I}_{1}\setminus\\{k\\}$,

6: end if

7:end for

####  III-B2 Homogeneous interconnected system

In this case, we consider the homogeneous interconnected system in the
following sense.

######  Definition 6.

An interconnected system $\langle(G_{i})_{i\in\mathcal{I}},\mathcal{E}\rangle$
is called homogeneous if $G_{i}=G_{j}$ and
$\mathcal{Q}_{i}=\mathcal{Q}_{j},\forall i,j\in\mathcal{I}$.

Algorithm 4 Contract construction for homogeneous systems

1:$G_{i},\mathcal{Q}_{i}$

2:try

3: calculate $\delta_{i}^{\star}$ by solving (7);

4: calculate $\zeta_{i}^{\star}$ by letting $\delta_{i}=\delta_{i}^{\star}$
and solving (8)

5: compute the corresp. iAGC
$C_{i}=(I_{\underline{W}_{i}^{\star}},I_{\underline{X}_{i}^{\star}},I_{\underline{Y}_{i}^{\star}})$

6:catch infeasible

7: return False;

8:end try

9:Assign all subsystem $G_{j},j\in\mathcal{I}$ with an iAGC
$C_{j}=(I_{\underline{W}_{j}^{\star}},I_{\underline{X}_{j}^{\star}},I_{\underline{Y}_{j}^{\star}})$
with
$\underline{W}_{j}^{\star}=\underline{W}_{j}^{\star},\underline{X}_{j}^{\star}=\underline{X}_{i}^{\star},\underline{Y}_{j}^{\star}=\underline{Y}_{i}^{\star}$.

10:if
$\underline{Y}_{j}^{\star}\subseteq\text{Proj}_{j}(\underline{W}_{i}^{\star})$
for all $j\in N(i)$ then

11: return True.

12:else

13: $\mathcal{Q}_{i}^{\prime}\leftarrow$ update the local safe region
$\mathcal{Q}_{i}$ by Alg. 2;

14: Goto Step 1 with updated safe region $\mathcal{Q}_{i}^{\prime}$

15:end if

Algorithm 4 starts with solving for one subsystem the maximal internal input
set and the corresponding minimal safe region. If the compatibility condition
is met, then we have verified the safety of the interconnected system;
otherwise, we will reduce the safe region by taking the set intersection in
Algorithm 2 and start the same process with the updated safe region
$\mathcal{Q}_{i}^{\prime}$.

We have the following results in this case:

######  Proposition 6.

Consider a homogeneous interconnected system as per Definition 6. Assume that
$\\{x_{i}\mathrel{\mathop{\ordinarycolon}}q_{i}(x_{i})\geq a\\}\subseteq
X_{i}^{0}$ for some $a>0$. Algorithm 4 has the following properties:

  1. 1.

Algorithm 4 terminates in finite steps and returns either True or False.

  2. 2.

If Algorithm 4 returns True, then iAGCs
$C_{i}=(I_{\underline{W}_{i}^{\star}},I_{\underline{X}_{i}^{\star}},I_{\underline{Y}_{i}^{\star}})$,
$i\in\mathcal{I}$ satisfy the conditions in Proposition 4.

  3. 3.

If Algorithm 4 returns False, then there exists no common contract
$C_{0}=(I_{\underline{W}_{0}},I_{\underline{X}_{0}},I_{\underline{Y}_{0}})$
such that $G_{i}\models C_{0},i\in\mathcal{I}$ and that the conditions in
Proposition 4 are satisfied under Assumption 1.

###### Proof.

Consider the case when $X_{i}^{0}$ is a subset of the local safe region
$\mathcal{Q}_{i}$ (otherwise, (7) yields infeasible in Step 2). There are
three possibilities: 1) (7) is infeasible, for which the algorithm terminates
with False; 2) the contract compatibility condition is satisfied, for which
the algorithm terminates with True; 3) the algorithm starts a new iteration.
For the third scenario, the local safe region $\mathcal{Q}_{i}$ is updated at
every iteration, and at each iteration, the set size gets smaller. In
particular, as the set size is parameterized by a scaler $\zeta$, which is
monotonically increasing, we know it has to converge to some number smaller
than $a$ (in which case we found compatible local contracts and the algorithm
terminates with True) or becomes larger than $a$ (in which case (7) yields
infeasible). In either case, the algorithm will terminate.

If Algorithm 4 returns True, from Step 9, the contract compatibility condition
holds for subsystem $G_{i}$. Noting that each subsystem has the same number of
parent nodes (implicit from the fact that each subsystem has the same output
set $Y_{i}=Y_{j}$ and the same internal input set $W_{i}=W_{j}$) and the same
assumption and guarantee sets (Step 8), we know this compatibility condition
also holds for all subsystems. Thus, Claim 2) is shown.

We show Claim 3) by contradiction. Suppose that a common local contract
$C_{0}=(I_{\underline{W}_{0}},I_{\underline{X}_{0}},I_{\underline{Y}_{0}})$
exists and the compatibility condition (9) holds. Denote the safe region with
which such an iAGC is obtained by $\mathcal{Q}_{0}$. That is, the functions
defining these sets fulfill the SOS constraints in (8) (but not necessarily
minimize the size of the safe region) and

$o(\mathcal{Q}_{0})\subseteq\text{Proj}_{i}(\underline{W}_{0}).$ (11)

Trivially, we know
$X_{i}^{0}\subseteq\mathcal{Q}_{0}\subseteq\mathcal{Q}_{i}$. Let the sequence
of the updated safe regions of Algorithm 4 be
$Q^{1}=\mathcal{Q}_{i},Q^{2},\ldots,Q^{M}$ (it is a finite sequence of sets
with shrinking size following Claim 1)). Without loss of generality, assume
$\mathcal{Q}^{r+1}\subsetneqq\mathcal{Q}_{0}\subseteq\mathcal{Q}^{r}$.
Following Proposition 3 items 2 and 4, we know
$\underline{W}_{0}\subseteq\left.\underline{W}^{\star}\right|_{\mathcal{Q}_{0}}\subseteq\left.\underline{W}^{\star}\right|_{\mathcal{Q}^{r}}$,
where $\left.\underline{W}^{\star}\right|_{\mathcal{Q}}$ represents the
maximal internal input set given safe region $\mathcal{Q}$. Recall $Q^{r+1}$
is obtained from Alg. 2, i.e., $Q^{r+1}$ is the largest safe region such that
$o(\mathcal{Q}^{r+1})\subseteq\text{Proj}_{i}(\left.\underline{W}^{\star}\right|_{\mathcal{Q}^{r}})$.
However, from (11), we know
$o(\mathcal{Q}_{0})\subseteq\text{Proj}_{i}(\underline{W}_{0})\subseteq\text{Proj}_{i}(\left.\underline{W}^{\star}\right|_{\mathcal{Q}^{r}})$
and $\mathcal{Q}^{r+1}\subsetneqq\mathcal{Q}_{0}$. This yields a
contradiction, which thus proves Claim 3). ∎

####  III-B3 General case

In the following, we provide a constructive approach for safety verification
of interconnected systems with general connectivity graphs and dynamics.

Algorithm 5 Contract construction for general systems

1:Subsystem $G_{i}$ and its safe region $\mathcal{Q}_{i},i\in\mathcal{I}$ and
connectivity graph $(\mathcal{I},\mathcal{E})$

2:run Algorithm 1;

3:for each $G_{i}\in\mathcal{I}_{1}$ do

4: try

5: calculate $\delta_{i}^{\star}$ by solving (7);

6: calculate $\zeta_{i}^{\star}$ by letting $\delta_{i}=\delta_{i}^{\star}$
and solving (8)

7: compute the corresp. iAGC
$(I_{\underline{W}_{i}^{\star}},I_{\underline{X}_{i}^{\star}},I_{\underline{Y}_{i}^{\star}})$

8: catch infeasible

9: return False;

10: end try

11:end for

12:if contract compatibility condition (9) does not hold then

13: update $\mathcal{Q}_{i}$ for all $i\in\mathcal{I}_{1}$ by Alg. 2

14: Goto Step 2

15:end if

16:return True.

######  Proposition 7.

For the general case, if Algorithm 5 returns True, then iAGCs
$C_{i}=(I_{\underline{W}_{i}},I_{\underline{X}_{i}},I_{\underline{Y}_{i}}),i\in\mathcal{I}$
satisfies the conditions in Proposition 4.

The proof is straightforward and omitted. We note that Algorithm 5 simplifies
to Algorithm 1 in the case of acyclic connectivity graph, and becomes
Algorithm 4 for homogeneous interconnected systems.

######  Remark 1 (Meaning of False).

It is worth highlighting that our completeness results are established under
Assumption 1, and that the algorithm returning False does not mean that the
interconnected system is unsafe. It simply means that we can not certify the
safety of the interconnected systems under the restrictions in Assumption 1.
When a False is encountered, one can group two or several subsystems together
and run the algorithms again by taking a group of subsystems as a single one.

##  IV Examples

###  IV-A Vehicular platooning: an acyclic example

Consider a vehicular platooning scenario adapted from [25] where $N+1$
autonomous vehicles are moving on a single-lane road. We assume that each
vehicle has the same length $l$, and has access to the state information of
its proceeding vehicle and the leading vehicle (hereafter referred to as the
leader). The dynamics relative to the leader is

$\displaystyle\dot{\tilde{p}}_{i}$ $\displaystyle=\tilde{v}_{i},$ (12)
$\displaystyle\dot{\tilde{v}}_{i}$
$\displaystyle=\tilde{u}_{i}-(\tilde{v}_{i}-\tilde{v}_{i-1})^{3},$

where
$\tilde{p}_{i}(t)=p_{0}(t)-p_{i}(t),\tilde{v}_{i}=v_{0}(t)-v_{i}(t),\tilde{u}_{i}(t)=u_{0}(t)-u_{i}(t)$,
$i=1,2,\ldots,N,k_{0}>0$. Here
$p_{i},v_{i},u_{i},\tilde{p}_{i},\tilde{v}_{i},\tilde{u}_{i}$ denote the
absolute position, the absolute velocity, the absolute control, the relative
position, the relative velocity, and the relative control of vehicle $i$,
respectively. We conveniently let $\tilde{v}_{0}=0$ for ease of notation.

Figure 1: Platooning scenario

(a) $X_{i}^{0}$ and $\mathcal{Q}_{i}$, $i=1,2,3$.

(b) AGC sets for vehicle $3$.

(c) AGC sets for vehicle $2$.

(d) AGC sets for vehicle $1$.

Figure 2:  Results for the platooning example.

Instead of looking at the relative dynamics in (12), we introduce a new
coordinate $x_{i}=(d_{i},\tilde{v}_{i})$ associated with vehicle $i$, where
$d_{i}=\tilde{p}_{i}-\tilde{p}_{i-1}-l$ denotes the distance between the front
of the $i$-th vehicle and the rear of its proceeding vehicle. See Figure 1 for
an illustration. The dynamics of this new state are given by

$\displaystyle\dot{d}_{i}$ $\displaystyle=\tilde{v}_{i}-\tilde{v}_{i-1}$ (13)
$\displaystyle\dot{\tilde{v}}_{i}$
$\displaystyle=-(\tilde{v}_{i}-\tilde{v}_{i-1})^{3}+\tilde{u}_{i}$

From the analysis above, we can model the vehicular platooning system as an
interconnected system. Each subsystem
$G_{i},i\in\mathcal{I}=\\{1,2,\ldots,N\\},$ is a continuous-time system

$G_{i}=(U_{i},W_{i},X_{i},Y_{i},X_{i}^{0},\mathcal{T}_{i})$

with $(d_{i},\tilde{v}_{i})$ as the state vector, and
$U_{i}=\mathbb{R},W_{i}=\left\\{\begin{smallmatrix}\emptyset,&i=1;\\\
\mathbb{R},&i\geq 2,\end{smallmatrix}\right.$,
$X_{i}=\mathbb{R}^{2},Y_{i}=\left\\{\begin{smallmatrix}\mathbb{R},&i\leq
N-1;\\\ \emptyset,&i=N,\end{smallmatrix}\right.,X_{i}^{0}\subset
X_{i},\mathcal{T}_{i}$ characterizes all the trajectories satisfying (13), and
the output map
$o_{i}\mathrel{\mathop{\ordinarycolon}}\left\\{\begin{smallmatrix}(d_{i},\tilde{v}_{i})\mapsto\tilde{v}_{i}&i\leq
N-1\\\ (d_{i},\tilde{v}_{i})\mapsto\emptyset&i=N\end{smallmatrix}\right.$. The
binary connectivity relation $\mathcal{E}$ is defined that
$(j,i)\in\mathcal{E}$ if and only if $j=i-1,i=2,3,\ldots,N$. One verifies that
$\\{G_{i}\\}_{i\in\mathcal{I}}$ is compatible for composition with respect to
$\mathcal{E}$. In the following analysis, we assume all the subsystems have
the same initial state set $X_{i}^{0}$ and safe region $\mathcal{Q}_{i}$,
which are
$X_{i}^{0}=\\{x_{i}\mathrel{\mathop{\ordinarycolon}}-x_{i}^{\top}Qx_{i}+q^{\top}x_{i}-899\geq
0\\},$ and the safe region
$\mathcal{Q}_{i}=\\{x_{i}\mathrel{\mathop{\ordinarycolon}}-x_{i}^{\top}Qx_{i}+q^{\top}x_{i}-800\geq
0\\}.$ where $Q=\begin{bmatrix}100&30\\\ 30&50\end{bmatrix}$ and
$q=(600,180)$. The initial state set and the safe region are depicted in
Fig.2(a). Each subsystem applies a local controller

$\tilde{u}_{i}=-(\tilde{v}_{i}-\tilde{v}_{i-1})-(d_{i}-3)-(d_{i}-3)^{3},i\in\mathcal{I}.$

Our task is to verify the safety of the interconnected system
$\langle(G_{i})_{i\in\mathcal{I}},\mathcal{E}\rangle$. In this scenario, we
consider $4$ vehicles ($N=3$). Algorithm 1 will be used for this example. We
adopt $d^{i}_{j_{k}}$ in $\text{Proj}_{k}(W_{i})$ in the form of
$d^{i}_{j_{k}}(x)=a^{2}-(\tilde{v}_{i-1}-b)^{2}$ to denote a bounded interval,
where $a$ and $b$ are determined accordingly.

Consider the vehicle $3$, which is the leaf node in the graph. One calculates
$d^{3}_{2}=2.439-\tilde{v}_{2}^{2}$ by projecting $\mathcal{Q}_{2}$ to the
$\tilde{v}_{2}$ coordinate. By solving (7) and (8), one obtains
$\delta^{\star}=1.704$ and $\zeta^{\star}=1.1147$. That is to say, a local
contract
$C_{3}=(I_{\underline{W}_{3}},I_{\underline{X}_{3}},I_{\underline{Y}_{3}})$
for vehicle $3$ is constructed with

$\displaystyle\underline{W}_{3}$
$\displaystyle=\\{\tilde{v}_{2}\mathrel{\mathop{\ordinarycolon}}0.735-\tilde{v}_{2}^{2}\geq
0\\},$ (14a) $\displaystyle\underline{X}_{3}$
$\displaystyle=\\{x_{3}\mathrel{\mathop{\ordinarycolon}}-x_{3}^{\top}Qx_{3}+q^{\top}x_{3}-801.115\geq
0\\}\\},$ (14b) $\displaystyle\underline{Y}_{3}$
$\displaystyle=o_{3}(\underline{X}_{3}).$ (14c)

$\underline{W}_{3}$ can be seen as the requirement from vehicle $3$ to vehicle
$2$. Following Algorithm 2 (and also by solving (10)), the updated safe region
for vehicle $2$ is
$\mathcal{Q}_{2}^{\prime}=\\{x_{2}\mathrel{\mathop{\ordinarycolon}}-x_{2}^{\top}Qx_{2}+q^{\top}x_{2}-869.85\geq
0\\}$. The initial set and the guarantee set of vehicle $3$ as well as the
updated safe region of vehicle $2$ are illustrated in Figure 2(b). We thus
follow the same procedures for vehicles $2$ and $1$, and obtain the results in
Figure 2(c) and Figure 2(d). Moreover, we calculate the assumption set
$\underline{W}_{1}=\\{\tilde{v}_{0}\mathrel{\mathop{\ordinarycolon}}0.019-\tilde{v}_{0}^{2}\geq
0\\},$ which holds true since $\tilde{v}_{0}(t)=0$ for all $t$. Thus, we have
constructed compatible local contracts for each vehicle. Following Proposition
4, we conclude that the platooning system is safe.

###  IV-B Room temperature: a homogeneous example

In the second example, we consider a room temperature regulation problem [26]
in a circular building as illustrated in Fig. 3. Each room has its temperature
$x_{i}$, which is affected by neighboring rooms, the heater, and the
environment as follows

$\displaystyle\dot{x}_{i}(t)$
$\displaystyle=\alpha(x_{i+1}+x_{i-1}-2x_{i})+\beta(t_{e}-x_{i})+\gamma(t_{h}-x_{i})u_{i},$
$\displaystyle y_{i}(t)$ $\displaystyle=x_{i},$

where $x_{i+1},x_{i-1}$ are the temperatures of room $i+1$ and $i-1$ (and we
conveniently let $x_{0}(t)=x_{N}(t),x_{N+1}(t)=x_{1}(t)$), $t_{e},t_{h}$ are
the temperatures of the environment and the heater, respectively.
$\alpha,\beta,\gamma$ are the respective conduction factors for the
neighboring room, the environment, and the heater. $u_{i}$ denotes the valve
control to the heater. Choose
$(t_{e},t_{h},\alpha,\beta,\gamma)=(-1,50,0.05,0.008,0.004)$, and

$u_{i}=0.05(x_{i+1}+x_{i-1}-2x_{i})+0.05(25-x_{i}).$

The initial set is $\mathcal{S}_{I,i}=[24,26]$ and the safe region is
$\mathcal{Q}_{i}=[20,30]$ for every room.

We can model the temperature system as an interconnected system. In
particular, each subsystem
$G_{i}=(U_{i},W_{i},X_{i},Y_{i},X_{i}^{0},\mathcal{T}_{i})$ has $x_{i}$ as the
state, $(x_{i-1},x_{i+1})$ as the internal input, $u_{i}$ as the external
input, $o_{i}(x_{i})=x_{i}$,
$U_{i},X_{i},Y_{i}=\mathbb{R},W_{i}=\mathbb{R}^{2},X_{i}^{0}=\\{x_{i}\mathrel{\mathop{\ordinarycolon}}1-(x_{i}-25)^{2}\geq
0\\}$, and
$\mathcal{Q}_{i}=\\{x_{i}\mathrel{\mathop{\ordinarycolon}}5^{2}-(x_{i}-25)^{2}\geq
0\\}$. The connectivity relation $\mathcal{E}$ is defined that
$(j,i)\in\mathcal{E}$ if and only if $j=i\pm 1,i=1,2,\ldots,N$. Per Definition
6, this is a homogeneous interconnected system and we will apply Algorithm 4
for this example.

At the first iteration, by solving (7) and (8), we obtain
$\delta^{\star}=20.575,\zeta^{\star}=0$. Thus, we have constructed a local
iAGC
$C_{i}=(I_{\underline{W}_{i}},I_{\underline{X}_{i}},I_{\underline{Y}_{i}})$
with

$\displaystyle\underline{W}_{i}=\\{(x_{i-1},x_{i+1})\mathrel{\mathop{\ordinarycolon}}-x_{j}^{2}+50x_{j}-620.575\geq
0,j=i\pm 1\\},$
$\displaystyle\underline{X}_{i}=\underline{Y}_{i}=\\{x_{i}\mathrel{\mathop{\ordinarycolon}}-x_{i}^{2}+50x_{i}-600\geq
0\\}.$

After assigning the same local contract to all subsystems, one verifies that
the contract compatibility condition (9) does not hold. According to Step 12
of Algorithm 4, we update the safe region for each room to be
$\mathcal{Q}_{i}^{\prime}=\\{x_{i}\mathrel{\mathop{\ordinarycolon}}-x_{i}^{2}+50x_{i}--620.575\geq
0\\}$ and start over. For the second iteration, we obtain local iAGC
$C_{i}=(I_{\underline{W}_{i}},I_{\underline{X}_{i}},I_{\underline{Y}_{i}})$
with

$\displaystyle\underline{W}_{i}=\\{(x_{i-1},x_{i+1})\mathrel{\mathop{\ordinarycolon}}-x_{j}^{2}+50x_{j}-622.138\geq
0,j=i\pm 1\\},$
$\displaystyle\underline{X}_{i}=\underline{Y}_{i}=\\{x_{i}\mathrel{\mathop{\ordinarycolon}}-x_{i}^{2}+50x_{i}-623.575\geq
0\\}.$

This time, one verifies that the compatibility condition (9) holds, and thus,
certifies the safety of the room temperature system. An illustration of the
assume and the guarantee sets is given in Fig. 4. We note that the computation
expense is not related to the number of rooms $N$, and only small-size SOS
optimization problems involving $3$ independent variables are to be solved.

Figure 3:  Room temperature scenario.

Figure 4:  Assume/guarantee sets for the room temperature example. Left:
iteration 1, right: iteration 2.

##  V Conclusions

In this work, we propose a safety verification scheme for interconnected
continuous-time nonlinear systems based on assume-guarantee contracts (AGCs)
and sum-of-squares (SOS) programs. The proposed scheme uses SOS optimization
to calculate local invariance AGCs by synthesizing local (control) barrier
functions, and then negotiates among neighboring subsystems at the contract
level. If the proposed algorithms find compatible local contracts, safety
property of the interconnected system is certified. We also show that the
algorithms will terminate in finite steps and will always find a solution when
one exists in the case of acyclic connectivity graphs or for homogeneous
systems. We also demonstrate the effectiveness of the proposed algorithms for
vehicle platooning and room temperature regulation examples.

## References

  * [1] S. Bansal, M. Chen, S. Herbert, and C. J. Tomlin, “Hamilton-Jacobi reachability: A brief overview and recent advances,” in _2017 IEEE 56th Annual Conference on Decision and Control (CDC)_. IEEE, 2017, pp. 2242–2253.
  * [2] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs for safety critical systems,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 8, pp. 3861–3876, 2016.
  * [3] E. F. Camacho and C. B. Alba, _Model predictive control_. Springer Science & Business Media, 2013.
  * [4] C. P. Bechlioulis and G. A. Rovithakis, “Robust adaptive control of feedback linearizable MIMO nonlinear systems with prescribed performance,” _IEEE Transactions on Automatic Control_ , vol. 53, no. 9, pp. 2090–2099, 2008.
  * [5] S. Prajna and A. Jadbabaie, “Safety verification of hybrid systems using barrier certificates,” in _International Workshop on Hybrid Systems: Computation and Control_. Springer, 2004, pp. 477–492.
  * [6] A. Clark, “Verification and synthesis of control barrier functions,” in _2021 60th IEEE Conference on Decision and Control (CDC)_ , 2021, pp. 6105–6112.
  * [7] H. Wang, K. Margellos, and A. Papachristodoulou, “Safety verification and controller synthesis for systems with input constraints,” _arXiv preprint arXiv:2204.09386_ , 2022.
  * [8] A. Robey, H. Hu, L. Lindemann, H. Zhang, D. V. Dimarogonas, S. Tu, and N. Matni, “Learning control barrier functions from expert demonstrations,” in _2020 59th IEEE Conference on Decision and Control (CDC)_. IEEE, 2020, pp. 3717–3724.
  * [9] A. Abate, D. Ahmed, A. Edwards, M. Giacobbe, and A. Peruffo, “FOSSIL: a software tool for the formal synthesis of Lyapunov functions and barrier certificates using neural networks,” in _Proceedings of the 24th International Conference on Hybrid Systems: Computation and Control_ , 2021, pp. 1–11.
  * [10] J. J. Choi, D. Lee, K. Sreenath, C. J. Tomlin, and S. L. Herbert, “Robust control barrier–value functions for safety-critical control,” in _2021 60th IEEE Conference on Decision and Control (CDC)_. IEEE, 2021, pp. 6814–6821.
  * [11] P. Jagtap, A. Swikir, and M. Zamani, “Compositional construction of control barrier functions for interconnected control systems,” in _Proceedings of the 23rd International Conference on Hybrid Systems: Computation and Control_ , 2020, pp. 1–11.
  * [12] Z. Lyu, X. Xu, and Y. Hong, “Small-gain theorem for safety verification of interconnected systems,” _Automatica_ , vol. 139, p. 110178, 2022.
  * [13] S. Coogan and M. Arcak, “A dissipativity approach to safety verification for interconnected systems,” _IEEE Transactions on Automatic Control_ , vol. 60, no. 6, pp. 1722–1727, 2014.
  * [14] P. Tabuada, _Verification and control of hybrid systems: a symbolic approach_. Springer Science & Business Media, 2009.
  * [15] A. Benveniste, B. Caillaud, D. Nickovic, R. Passerone, J.-B. Raclet, P. Reinkemeier, A. Sangiovanni-Vincentelli, W. Damm, T. A. Henzinger, K. G. Larsen _et al._ , “Contracts for system design,” _Foundations and Trends® in Electronic Design Automation_ , vol. 12, no. 2-3, pp. 124–400, 2018.
  * [16] A. Saoud, A. Girard, and L. Fribourg, “Assume-guarantee contracts for continuous-time systems,” _Automatica_ , vol. 134, p. 109910, 2021.
  * [17] E. S. Kim, M. Arcak, and S. A. Seshia, “A small gain theorem for parametric assume-guarantee contracts,” in _Proceedings of the 20th International Conference on Hybrid Systems: Computation and Control_ , 2017, pp. 207–216.
  * [18] K. Ghasemi, S. Sadraddini, and C. Belta, “Compositional synthesis via a convex parameterization of assume-guarantee contracts,” in _Proceedings of the 23rd International Conference on Hybrid Systems: Computation and Control_ , 2020, pp. 1–10.
  * [19] A. Eqtami and A. Girard, “A quantitative approach on assume-guarantee contracts for safety of interconnected systems,” in _2019 18th European Control Conference (ECC)_. IEEE, 2019, pp. 536–541.
  * [20] B. M. Shali, A. van der Schaft, and B. Besselink, “Composition of behavioural assume-guarantee contracts,” _IEEE Transactions on Automatic Control_ , 2022\.
  * [21] S. Liu, A. Saoud, P. Jagtap, D. V. Dimarogonas, and M. Zamani, “Compositional synthesis of signal temporal logic tasks via assume-guarantee contracts,” in _2022 IEEE 61st Conference on Decision and Control (CDC)_. IEEE, 2022, pp. 2184–2189.
  * [22] F. Blanchini, “Set invariance in control,” _Automatica_ , vol. 35, no. 11, pp. 1747–1767, 1999.
  * [23] S. Prajna, A. Papachristodoulou, and P. A. Parrilo, “Introducing sostools: A general purpose sum of squares programming solver,” in _Proceedings of the 41st IEEE Conference on Decision and Control, 2002._ , vol. 1. IEEE, 2002, pp. 741–746.
  * [24] M. Arcak, C. Meissen, and A. Packard, _Networks of dissipative systems: compositional certification of stability, performance, and safety_. Springer, 2016.
  * [25] S. Liu and M. Zamani, “Compositional synthesis of almost maximally permissible safety controllers,” in _2019 American Control Conference (ACC)_. IEEE, 2019, pp. 1678–1683.
  * [26] A. Girard, G. Gössler, and S. Mouelhi, “Safety controller synthesis for incrementally stable switched systems using multiscale symbolic models,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 6, pp. 1537–1549, 2015.