&=-{\lnoh \tN u, \gd v \rnohs}+\overline{\lnoh \tN v, \gd u \rnohs}.
\end{align}

[1]
H. Abels, G. Grubb, and I. G. Wood, Extension theory and
Kreĭn-type resolvent formulas for nonsmooth boundary value problems,
J. Funct. Anal., 266 (2014), pp. 4037–4100.

[2]
N. I. Akhiezer and I. M. Glazman, Theory of linear operators in
Hilbert space, Dover Publications, Inc., New York, 1993.
Translated from the Russian and with a preface by Merlynd Nestell,
Reprint of the 1961 and 1963 translations, Two volumes bound as one.

[3]
S. Albeverio, F. Gesztesy, R. Høegh-Krohn, and H. Holden, Solvable models in quantum mechanics, AMS Chelsea Publishing, Providence,
RI, second ed., 2005.
With an appendix by Pavel Exner.

[4]
S. Albeverio and K. Pankrashkin, A remark on Krein's resolvent
formula and boundary conditions, J. Phys. A, 38 (2005), pp. 4859–4864.

[5]
A. Alonso and B. Simon, The Birman-Kreĭn-Vishik theory of
selfadjoint extensions of semibounded operators, J. Operator Theory, 4
(1980), pp. 251–270.

[6]
W. O. Amrein and D. B. Pearson, $M$ operators: a generalization of
Weyl-Titchmarsh theory, J. Comput. Appl. Math., 171 (2004), pp. 1–26.

[7]
Y. Arlinskiĭ, Boundary triplets and maximal accretive
extensions of sectorial operators, in Operator methods for boundary value
problems, vol. 404 of London Math. Soc. Lecture Note Ser., Cambridge Univ.
Press, Cambridge, 2012, pp. 35–72.

[8]
V. I. Arnold, On a characteristic class entering into conditions of
quantization, Funkcional. Anal. i Priložen., 1 (1967), pp. 1–14.

[9]
height 2pt depth -1.6pt width 23pt, Sturm theorems and
symplectic geometry, Funktsional. Anal. i Prilozhen., 19 (1985), pp. 1–10,

[10]
A. Avila, On the spectrum and Lyapunov exponent of limit periodic
Schrödinger operators, Comm. Math. Phys., 288 (2009), pp. 907–918.

[11]
C. Bandle, Isoperimetric inequalities and applications, vol. 7 of
Monographs and Studies in Mathematics, Pitman (Advanced Publishing Program),
Boston, Mass.-London, 1980.

[12]
J. Behrndt, S. Hassi, and H. de Snoo, Boundary Value Problems, Weyl
Functions, and Differential Operators, Springer International Publishing,

[13]
J. Behrndt and M. Langer, Boundary value problems for elliptic
partial differential operators on bounded domains, J. Funct. Anal., 243
(2007), pp. 536–565.

[14]
J. Behrndt and T. Micheler, Elliptic differential operators on
Lipschitz domains and abstract boundary value problems, J. Funct. Anal.,
267 (2014), pp. 3657–3709.

[15]
J. Behrndt, J. Rohleder, and S. Stadler, Eigenvalue inequalities for
Schrödinger operators on unbounded Lipschitz domains, J. Spectr.
Theory, 8 (2018), pp. 493–508.

[16]
G. Berkolaiko and P. Kuchment, Introduction to quantum graphs,
vol. 186 of Mathematical Surveys and Monographs, American Mathematical
Society, Providence, RI, 2013.

[17]
G. Berkolaiko, P. Kuchment, and U. Smilansky, Critical partitions
and nodal deficiency of billiard eigenfunctions, Geom. Funct. Anal., 22
(2012), pp. 1517–1540.

[18]
M. S. Birman, On the theory of self-adjoint extensions of positive
definite operators, Mat. Sb. N.S., 38(80) (1956), pp. 431–450.

[19]
B. Booss-Bavnbek and K. Furutani, The Maslov index: a functional
analytical definition and the spectral flow formula, Tokyo J. Math., 21
(1998), pp. 1–34.

[20]
B. Booss-Bavnbek and C. Zhu, General spectral flow formula for fixed
maximal domain, Cent. Eur. J. Math., 3 (2005), pp. 558–577.

[21]
B. Booß Bavnbek and C. Zhu, The Maslov index in weak symplectic
functional analysis, Ann. Global Anal. Geom., 44 (2013), pp. 283–318.

[22]
height 2pt depth -1.6pt width 23pt, The Maslov index
in symplectic Banach spaces, Mem. Amer. Math. Soc., 252 (2018), pp. x+118.

[23]
B. M. Brown, G. Grubb, and I. G. Wood, $M$-functions for closed
extensions of adjoint pairs of operators with applications to elliptic
boundary problems, Math. Nachr., 282 (2009), pp. 314–347.

[24]
M. Brown, J. Hinchcliffe, M. Marletta, S. Naboko, and I. Wood, The
abstract Titchmarsh-Weyl $M$-function for adjoint operator pairs and
its relation to the spectrum, Integral Equations Operator Theory, 63 (2009),
pp. 297–320.

[25]
M. Brown, M. Marletta, S. Naboko, and I. Wood, Boundary triplets and
$M$-functions for non-selfadjoint operators, with applications to elliptic
PDEs and block operator matrices, J. Lond. Math. Soc. (2), 77 (2008),
pp. 700–718.

[26]
S. E. Cappell, R. Lee, and E. Y. Miller, On the Maslov index,
Comm. Pure Appl. Math., 47 (1994), pp. 121–186.

[27]
F. Chardard, F. Dias, and T. J. Bridges, Fast computation of the
Maslov index for hyperbolic linear systems with periodic coefficients, J.
Phys. A, 39 (2006), pp. 14545–14557.

[28]
F. Chardard, F. Dias, and T. J. Bridges, Computing the Maslov
index of solitary waves. I. Hamiltonian systems on a four-dimensional
phase space, Phys. D, 238 (2009), pp. 1841–1867.

[29]
height 2pt depth -1.6pt width 23pt, On the Maslov
index of multi-pulse homoclinic orbits, Proc. R. Soc. Lond. Ser. A Math.
Phys. Eng. Sci., 465 (2009), pp. 2897–2910.

[30]
height 2pt depth -1.6pt width 23pt, Computing the
Maslov index of solitary waves, Part 2: Phase space with dimension
greater than four, Phys. D, 240 (2011), pp. 1334–1344.

[31]
S. Clark, F. Gesztesy, and M. Mitrea, Boundary data maps for
Schrödinger operators on a compact interval, Math. Model. Nat.
Phenom., 5 (2010), pp. 73–121.

[32]
S. Clark, F. Gesztesy, R. Nichols, and M. Zinchenko, Boundary data
maps and Krein's resolvent formula for Sturm-Liouville operators on a
finite interval, Oper. Matrices, 8 (2014), pp. 1–71.

[33]
G. M. Coclite, A. Favini, C. G. Gal, G. R. Goldstein, J. A. Goldstein,
E. Obrecht, and S. Romanelli, The role of Wentzell boundary
conditions in linear and nonlinear analysis, in Advances in nonlinear
analysis: theory methods and applications, vol. 3 of Math. Probl. Eng.
Aerosp. Sci., Camb. Sci. Publ., Cambridge, 2009, pp. 277–289.

[34]
G. M. Coclite, A. Favini, G. R. Goldstein, J. A. Goldstein, and
S. Romanelli, Continuous dependence on the boundary conditions for the
Wentzell Laplacian, Semigroup Forum, 77 (2008), pp. 101–108.

[35]
G. Cox, C. K. R. T. Jones, Y. Latushkin, and A. Sukhtayev, The
Morse and Maslov indices for multidimensional Schrödinger operators
with matrix-valued potentials, Trans. Amer. Math. Soc., 368 (2016),
pp. 8145–8207.

[36]
G. Cox, C. K. R. T. Jones, and J. L. Marzuola, A Morse index
theorem for elliptic operators on bounded domains, Comm. Partial
Differential Equations, 40 (2015), pp. 1467–1497.

[37]
G. Cox, C. K. R. T. Jones, and J. L. Marzuola, Manifold
decompositions and indices of Schrödinger operators, Indiana Univ.
Math. J., 66 (2017), pp. 1573–1602.

[38]
J. L. Dalecki and M. G. Krein, Stability of solutions of
differential equations in Banach space, American Mathematical Society,
Providence, R.I., 1974.
Translated from the Russian by S. Smith, Translations of Mathematical
Monographs, Vol. 43.

[39]
D. Damanik, J. Fillman, and M. Lukic, Limit-periodic continuum
Schrödinger operators with zero measure Cantor spectrum, J. Spectr.
Theory, 7 (2017), pp. 1101–1118.

[40]
J. Deng and C. Jones, Multi-dimensional Morse index theorems and a
symplectic view of elliptic boundary value problems, Trans. Amer. Math.
Soc., 363 (2011), pp. 1487–1508.

[41]
V. Derkach and M. Malamud, Extension theory of symmetric operators
and boundary value problems (In Russian), vol. 104 of Transactions of
Math. Inst. of Ukrainian Natl. Academy, Math. Inst. of Ukrainian Natl.
Academy, 2017.

[42]
V. A. Derkach and M. M. Malamud, Generalized resolvents and the
boundary value problems for Hermitian operators with gaps, J. Funct.
Anal., 95 (1991), pp. 1–95.

[43]
K.-J. Engel and R. Nagel, One-parameter semigroups for linear
evolution equations, vol. 194 of Graduate Texts in Mathematics,
Springer-Verlag, New York, 2000.
With contributions by S. Brendle, M. Campiti, T. Hahn, G. Metafune,
G. Nickel, D. Pallara, C. Perazzoli, A. Rhandi, S. Romanelli and R.

[44]
W. N. Everitt and L. Markus, Elliptic partial differential operators
and symplectic algebra, Mem. Amer. Math. Soc., 162 (2003), pp. x+111.

[45]
N. Filonov, On an inequality for the eigenvalues of the Dirichlet
and Neumann problems for the Laplace operator, Algebra i Analiz, 16
(2004), pp. 172–176.

[46]
R. L. Frank and A. Laptev, Inequalities between Dirichlet and
Neumann eigenvalues on the Heisenberg group, Int. Math. Res. Not. IMRN,
(2010), pp. 2889–2902.

[47]
L. Friedlander, Some inequalities between Dirichlet and Neumann
eigenvalues, Arch. Rational Mech. Anal., 116 (1991), pp. 153–160.

[48]
K. Furutani, Fredholm-Lagrangian-Grassmannian and the Maslov
index, J. Geom. Phys., 51 (2004), pp. 269–331.

[49]
P. R. Garabedian, Partial differential equations, Chelsea
Publishing Co., New York, second ed., 1986.

[50]
F. Gesztesy, Y. Latushkin, K. A. Makarov, F. Sukochev, and Y. Tomilov,
The index formula and the spectral shift function for relatively trace
class perturbations, Adv. Math., 227 (2011), pp. 319–420.

[51]
F. Gesztesy, Y. Latushkin, M. Mitrea, and M. Zinchenko, Nonselfadjoint operators, infinite determinants, and some applications,
Russ. J. Math. Phys., 12 (2005), pp. 443–471.

[52]
F. Gesztesy and M. Mitrea, Generalized Robin boundary conditions,
Robin-to-Dirichlet maps, and Krein-type resolvent formulas for
Schrödinger operators on bounded Lipschitz domains, in Perspectives
in partial differential equations, harmonic analysis and applications,
vol. 79 of Proc. Sympos. Pure Math., Amer. Math. Soc., Providence, RI, 2008,
pp. 105–173.

[53]
height 2pt depth -1.6pt width 23pt, Nonlocal Robin
Laplacians and some remarks on a paper by Filonov on eigenvalue
inequalities, J. Differential Equations, 247 (2009), pp. 2871–2896.

[54]
height 2pt depth -1.6pt width 23pt, A description of all
self-adjoint extensions of the Laplacian and Kreĭn-type resolvent
formulas on non-smooth domains, J. Anal. Math., 113 (2011), pp. 53–172.

[55]
G. R. Goldstein, Derivation and physical interpretation of general
boundary conditions, Adv. Differential Equations, 11 (2006), pp. 457–480.

[56]
V. I. Gorbachuk and M. L. Gorbachuk, Boundary value problems for
operator differential equations, vol. 48 of Mathematics and its Applications
(Soviet Series), Kluwer Academic Publishers Group, Dordrecht, 1991.
Translated and revised from the 1984 Russian original.

[57]
P. Grinfeld, Hadamard's formula inside and out, J. Optim. Theory
Appl., 146 (2010), pp. 654–690.

[58]
G. Grubb, A characterization of the non-local boundary value
problems associated with an elliptic operator, Ann. Scuola Norm. Sup. Pisa
Cl. Sci. (3), 22 (1968), pp. 425–513.

[59]
G. Grubb, Krein resolvent formulas for elliptic boundary problems in
nonsmooth domains, Rend. Semin. Mat. Univ. Politec. Torino, 66 (2008),
pp. 271–297.

[60]
J. Hadamard, Œuvres de Jacques Hadamard. Tomes I, II,
III, IV, Comité de publication des oeuvres de Jacques Hadamard: M.
Fréchet, P. Levy, S. Mandelbrojt, L. Schwartz, Éditions du Centre
National de la Recherche Scientifique, Paris, 1968.

[61]
A. Hassell, Ergodic billiards that are not quantum unique ergodic,
Ann. of Math. (2), 171 (2010), pp. 605–619.
With an appendix by the author and Luc Hillairet.

[62]
D. Henry, Perturbation of the boundary in boundary-value problems of
partial differential equations, vol. 318 of London Mathematical Society
Lecture Note Series, Cambridge University Press, Cambridge, 2005.
With editorial assistance from Jack Hale and Antônio Luiz

[63]
P. Howard and A. Sukhtayev, The Maslov and Morse indices for
Schrödinger operators on $[0,1]$, J. Differential Equations, 260
(2016), pp. 4499–4549.

[64]
P. Howard and A. Sukhtayev, The Maslov and Morse indices for
Sturm-Liouville systems on the half-line, Discrete Contin. Dyn. Syst.,
40 (2020), pp. 983–1012.

[65]
T. Kato, Perturbation theory for linear operators, Classics in
Mathematics, Springer-Verlag, Berlin, 1995.
Reprint of the 1980 edition.

[66]
J. B. Keller, Corrected Bohr-Sommerfeld quantum conditions for
nonseparable systems, Ann. Physics, 4 (1958), pp. 180–188.

[67]
S. Kesavan, Topics in functional analysis and applications, John
Wiley & Sons, Inc., New York, 1989.

[68]
M. Krein, The theory of self-adjoint extensions of semi-bounded
Hermitian transformations and its applications. I, Rec. Math. [Mat.
Sbornik] N.S., 20(62) (1947), pp. 431–495.

[69]
M. G. Kreĭn, The theory of self-adjoint extensions of
semi-bounded Hermitian transformations and its applications. II, Mat.
Sbornik N.S., 21(63) (1947), pp. 365–404.

[70]
Y. Latushkin and A. Pogan, The dichotomy theorem for evolution
bi-families, J. Differential Equations, 245 (2008), pp. 2267–2306.

[71]
Y. Latushkin and S. Sukhtaiev, The Maslov index and the spectra of
second order elliptic operators, Adv. Math., 329 (2018), pp. 422–486.

[72]
Y. Latushkin, S. Sukhtaiev, and A. Sukhtayev, The Morse and
Maslov indices for Schrödinger operators, J. Anal. Math., 135
(2018), pp. 345–387.

[73]
Y. Latushkin and A. Sukhtayev, Hadamard-type formulas via the
Maslov form, J. Evol. Equ., 17 (2017), pp. 443–472.

[74]
Y. Latushkin and Y. Tomilov, Fredholm differential operators with
unbounded coefficients, J. Differential Equations, 208 (2005), pp. 388–429.

[75]
V. È. Lyantse and O. G. Storozh, Metody teorii neogranichennykh
operatorov (Methods of the theory of unbounded operators), “Naukova
Dumka”, Kiev, 1983.

[76]
M. M. Malamud and V. I. Mogilevskii, Kreĭn type formula for
canonical resolvents of dual pairs of linear relations, Methods Funct. Anal.
Topology, 8 (2002), pp. 72–100.

[77]
V. Maslov, Theory of Perturbations and Asymptotic Methods, Izdat.
Moskov. Gos. Univ., Moscow 1965 (French translation Dunod, Paris, 1972).

[78]
W. McLean, Strongly elliptic systems and boundary integral
equations, Cambridge University Press, Cambridge, 2000.

[79]
S. Nakamura, A remark on the Dirichlet-Neumann decoupling and
the integrated density of states, J. Funct. Anal., 179 (2001), pp. 136–152.

[80]
K. Pankrashkin, Resolvents of self-adjoint extensions with mixed
boundary conditions, Rep. Math. Phys., 58 (2006), pp. 207–221.

[81]
A. Posilicano, Self-adjoint extensions of restrictions, Oper.
Matrices, 2 (2008), pp. 483–506.

[82]
A. Posilicano and L. Raimondi, Krein's resolvent formula for
self-adjoint extensions of symmetric second-order elliptic differential
operators, J. Phys. A, 42 (2009), pp. 015204, 11.

[83]
J. W. S. Rayleigh, Baron, The Theory of Sound, Dover
Publications, New York, N. Y., 1945.
2d ed.

[84]
M. Reed and B. Simon, Methods of modern mathematical physics. I.
Functional analysis, Academic Press, New York-London, 1972.

[85]
height 2pt depth -1.6pt width 23pt, Methods of modern
mathematical physics. IV. Analysis of operators, Academic Press
[Harcourt Brace Jovanovich Publishers], New York, 1978.

[86]
F. Rellich, Darstellung der Eigenwerte von $\Delta u+\lambda
u=0$ durch ein Randintegral, Math. Z., 46 (1940), pp. 635–636.

[87]
height 2pt depth -1.6pt width 23pt, Perturbation theory
of eigenvalue problems, Assisted by J. Berkowitz. With a preface by Jacob T.
Schwartz, Gordon and Breach Science Publishers, New York-London-Paris, 1969.

[88]
J. Robbin and D. Salamon, The spectral flow and the Maslov index,
Bull. London Math. Soc., 27 (1995), pp. 1–33.

[89]
J. Rohleder, Strict inequality of Robin eigenvalues for elliptic
differential operators on Lipschitz domains, J. Math. Anal. Appl., 418
(2014), pp. 978–984.

[90]
Y. Safarov, On the comparison of the Dirichlet and Neumann
counting functions, in Spectral theory of differential operators, vol. 225
of Amer. Math. Soc. Transl. Ser. 2, Amer. Math. Soc., Providence, RI, 2008,
pp. 191–204.

[91]
B. Sandstede and A. Scheel, On the structure of spectra of modulated
travelling waves, Math. Nachr., 232 (2001), pp. 39–93.

[92]
height 2pt depth -1.6pt width 23pt, Relative Morse
indices, Fredholm indices, and group velocities, Discrete Contin. Dyn.
Syst., 20 (2008), pp. 139–158.

[93]
K. Schmüdgen, Unbounded self-adjoint operators on Hilbert
space, vol. 265 of Graduate Texts in Mathematics, Springer, Dordrecht, 2012.

[94]
B. Simon, On the genericity of nonvanishing instability intervals in
Hill's equation, Ann. Inst. H. Poincaré Sect. A (N.S.), 24 (1976),
pp. 91–93.

[95]
M. E. Taylor, Partial differential equations I. Basic theory,
vol. 115 of Applied Mathematical Sciences, Springer, New York, second ed.,

[96]
G. Teschl, Mathematical methods in quantum mechanics, vol. 157 of
Graduate Studies in Mathematics, American Mathematical Society, Providence,
RI, second ed., 2014.
With applications to Schrödinger operators.

[97]
M. I. Višik, On boundary problems for systems of elliptic
differential equations and on the stability of their solutions, Doklady
Akad. Nauk SSSR (N.S.), 86 (1952), pp. 645–648.

In addition to Hypothesis <ref> we assume, for a given $t_0\in[0,1]$, that

\begin{equation}\label{2.2b}
\|(\cA_t-\bfi)^{-1}\|_{\cB(\cH, \cD)}=\cO(1) \text{ as $ t\rightarrow t_0$}.
\end{equation}

hyp1.3i3 Assume that  for a given $t_0\in[0,1]$ one has

\begin{equation}\lb{2.2bo}
\big\|(\cA_t-\bfi)^{-1}-(\cA_{t_0}-\bfi)^{-1}\big\|_{\cB(\cH, \cD)}=o(1),\ t\rightarrow t_0.
\end{equation}

Fix $t_0\in[0,1]$, $\zeta_0\not\in\Sp(H_{t_0})$ and define
\[\cU_\epsilon=\{(t,\zeta)\in[0,1]\times\bbC: |t-t_0|\le\epsilon, |\zeta-\zeta_0|\le\epsilon\} \text{ for $\epsilon>0$}.\]

* 18i1Assume Hypothesis <ref> and suppose that the mappings $t\mapsto \tr_t$, $t\mapsto V_t$, $t\mapsto Q_t$ are continuous at $t_0$. Then there exists an $\varepsilon>0$ such that if $(t,\zeta)\in\cU_\varepsilon$ then $\zeta\not\in\Sp(H_t)$ and
the operator valued function $t\mapsto R_t(\zeta)\in\cB(\cH)$ is continuous at $t_0$ uniformly for $|\zeta-\zeta_0|<\varepsilon$.

* 18i2Assume Hypothesis <ref> and suppose that the mappings $t\mapsto \tr_t$, $t\mapsto V_t$, $t\mapsto Q_t$ are continuous at $t_0$. Then there exists a constant $c>0$ such that for all $(t,\zeta)\in\cU_\varepsilon$ one has

\begin{align}\lb{new1.29n}
&\| R_{t}(\zeta)- R_{t_0}(\zeta) \|_{\cB(\cH)}\le c|t-t_0|.
\end{align}

* 18i3
Assume Hypothesis <ref> and suppose that the mappings $t\mapsto \tr_t$, $t\mapsto V_t$, $t\mapsto Q_t$ are differentiable at $t_0$.  Then for some $\varepsilon>0$ the following asymptotic expansion holds uniformly for $|\zeta-\zeta_0|<\varepsilon$

\begin{align}
\begin{split}\lb{new1.43nn}
R_t(\zeta)\underset{t\rightarrow t_0}{=}R_{t_0}(\zeta)&+
\big[R_{t_0}(\zeta)\dot V_{t_0}R_{t_0}(\zeta)+(\tr_{t_0}  R_{t_0}(\overline{\zeta}))^*\dot Q_{t_0}J \tr_{t_0}    R_{t_0}(\zeta)\\
&+(\tr_{t_0}  R_{t_0}(\overline{\zeta}))^*J\dot \tr_{t_0}  R_{t_0}(\zeta)\big](t-t_0)+o(t-t_0),\ \text{\ in\ }\cB(\cH).
\end{split}
\end{align}

In particular, the function $t\mapsto R_t(\zeta_0)$ is differentiable at $t=t_0$ and satisfies the following Riccati equation

\begin{align}
\begin{split}\lb{derR}
\dot{R}_{t_0}(\zeta_0)=R_{t_0}(\zeta_0)\dot V_{t_0}R_{t_0}(\zeta_0)&+(\tr_{t_0}  R_{t_0}(\overline{\zeta_0}))^*\dot Q_{t_0}J \tr_{t_0}    R_{t_0}(\zeta_0)\\
&+(\tr_{t_0}  R_{t_0}(\overline{\zeta_0}))^*J\dot \tr_{t_0}  R_{t_0}(\zeta_0).
\end{split}
\end{align}

§.§ The Krein formula via Theorem <ref>
In this section, we derive the well-known Krein formula, <cit.>,  from Theorem <ref>.  Following common convention we define one self-adjoint extension by

\begin{equation}\lb{1.35a}
\cA_1:=A^*\upharpoonright_{\ker(\Gamma_0)},
\end{equation}

and compare its resolvent with that of another arbitrary self-adjoint extension.  For $\zeta\in\bbC\setminus\bbR$ we define $\gamma(\zeta):=(\Gamma_0\upharpoonright_{\ker(A^*-\zeta)})^{-1}$ and notice that $\gamma(\zeta)\in\cB(\mathfrak H, \cH)$. This operator-valued function can be extended analytically to $\bbC\setminus\spec (\cA_1)$ giving rise to $M(\zeta):=\Gamma_1\gamma(\zeta)$, $\zeta\in \bbC\setminus\spec(\cA_1)$. Slightly abusing notation we do not distinguish between the operator $M(\zeta)$ and its graph, in particular,  for $\cG_2\in\Lambda(\bndra)$ we write $\cG_2-M(\zeta):=\cG_2-\gr(M(\zeta))$ and treat both terms in the right-hand side as linear relations, see, e.g., <cit.>. We write

\begin{equation}
\end{equation}

\[\ker(\cG_2-M(\zeta)):=\{x\in\mathfrak{H}: (x,0)\in(\cG_2-M(\zeta))\}=\{0\},\]

\[\ran(\cG_2-M(\zeta)):=\{y\in\mathfrak H:\exists x\in\mathfrak{H}\  {s.t.}\  (x,y)\in(\cG_2-M(\zeta)) \}=\mathfrak H.\]

In this case there exists an operator in $\cB(\mathfrak{H})$ whose graph is given by

\begin{equation}
\{(y,x) \in \bndra: (x,y)\in(\cG_2-M(\zeta)) \},
\end{equation}

such an operator is denoted by   $(\cG_2-M(\zeta))^{-1}$.
Assume (<ref>), and let $\cA_2$ be an arbitrary self-adjoint extension of $A$, then

\begin{equation}\lb{1.23}
\end{equation}

for $\zeta\not\in\Sp(\cA_1)\cup \Sp(\cA_2)$.

We denote $R_j(\zeta):=(\cA_j-\zeta)^{-1}$, $\cG_j:=\tr (\dom(\cA_j))$, $j=1,2$. The following assertions will be used throughout the proof

\begin{align}
&(\Gamma_0R_2(\zeta)f, \Gamma_1 R_2(\zeta)f)\in\cG_2, \text{\ by definition of $\cA_2$,}\lb{1.29}\\
&\Gamma_0 (I_{\cH}-\gamma(\zeta)\Gamma_0)=0,\text{\,that is,\,}\ran(I_{\cH}-\gamma(\zeta)\Gamma_0)\subset \dom(\cA_1),\lb{1.30}\\
&(A^*-\zeta)\gamma(\zeta)=0, \text{\ by definition of $\gamma(\zeta)$},\lb{1.31}\\
&\gamma(\zeta)^*=\Gamma_1 R_1(\overline{\zeta})\text{\ cf. \cite[Lemma 1]{DM}},\lb{1.31a}
\end{align}

Next, let us note the identity

\begin{equation}
\Gamma_1R_1(\zeta)f= \Gamma_1R_2(\zeta)f-M(\zeta) \Gamma_0 R_2(\zeta)f , f\in\cH,
\end{equation}

which can be derived by combining (<ref>) and (<ref>) as follows

\begin{align*}
\begin{split}
&\Gamma_1R_2(\zeta)f-M(\zeta) \Gamma_0 R_2(\zeta)f=\Gamma_1R_2(\zeta)f-\Gamma_1\gamma(\zeta) \Gamma_0 R_2f\\
&\quad=\Gamma_1( I_{\cH}-\gamma(\zeta)\Gamma_1)R_2(\zeta)f{=}\Gamma_1R_1(\zeta)(\cA_1-\zeta)( I_{\cH}-\gamma(\zeta)\Gamma_1)R_2(\zeta)f\\
&\quad=\Gamma_1R_1(\zeta)(A^*-\zeta)R_2(\zeta)f=\Gamma_1R_1(\zeta)f,\ f\in\cH.
\end{split}
\end{align*}

We also note that

\begin{align}
(\Gamma_0R_2(\zeta)f, \Gamma_1R_1(\zeta)f)\in (\cG_2-M(\zeta))\text{\ for all\ }f\in\cH,\lb{1.26}
\end{align}

Indeed, combining (<ref>) and (<ref>) we get

\begin{equation*}
(\Gamma_0R_2(\zeta)f, \Gamma_1R_1(\zeta)f)=(\Gamma_0R_2(\zeta)f,\Gamma_1R_2(\zeta)f-M(\zeta) \Gamma_0 R_2(\zeta)f)\in \cG_2-M(\zeta),
\end{equation*}

as asserted in (<ref>).

We are now ready to prove (<ref>) and consequently justify the existence and boundedness  of the second term in the right-hand side of (<ref>). First, we claim that

\begin{equation}\lb{1.34}
\ker(\cA_2-\zeta)=\gamma(\zeta) \ker (\cG_2-M(\zeta)),
\end{equation}

which yields  $\ker (\cG_2-M(\zeta))=\{0\}$ since $\zeta\in\bbC\setminus \Sp(\cA_2)$.
Let us notice that

\begin{equation*}
\ker (\cG_2-M(\zeta))=\{x\in \mathfrak{H}: (x, M(\zeta)x)\in \cG_2\}.
\end{equation*}

Then for $x\in \ker (\cG_2-M(\zeta))$ we have

\begin{equation*}
\tr \gamma(\zeta) x=(\Gamma_0\gamma(\zeta) x, \Gamma_1\gamma(\zeta) x)=(x, M(\zeta) x)\in\cG_2,
\end{equation*}

hence $\gamma(\zeta) x\in\dom(\cA_2)$ and

\[\gamma(\zeta) \ker (\cG_2-M(\zeta))\subset\ker(\cA_2-\zeta).\]

To show the opposite inclusion, fix $f\in \ker(\cA_2-\zeta)$ and notice that

\[(\Gamma_0f, M(\zeta)\Gamma_0f)= (\Gamma_0f, \Gamma_1f)\in\cG_2.\]

Hence $\Gamma_0f\in \ker(\cG_2-M(\zeta))$, which in turn implies
\[f=\gamma(\zeta)\Gamma_0f\in \gamma(\zeta) \ker (\cG_2-M(\zeta)),\]
as asserted.

Next, we claim that $\ran(\cG_2-M(\zeta))=\mathfrak{H}$. Indeed, by (<ref>) for any $x\in\mathfrak{H}$ there exists $g\in\dom(A^*)$ such that $\Gamma_1g=x$ and $\Gamma_0g=0$. Since $g\in\dom(\cA_1)$ there is some $f=f(x)$ such that $g=R_1(\zeta)f$. Then (<ref>) and (<ref>) yield

\begin{equation}\lb{1.33.5}
\end{equation}

as claimed. We also note that the vector $f=f(x)$ in (<ref>) is unique as if $f$ satisfies

\[0=\Gamma_1R_2(\zeta)f-M(\zeta)\Gamma_0R_2(\zeta)f\]

then $f\in\ker(\cG_2-M(\zeta))=\{0\}$.

To finish the proof of (<ref>) we now define the operator $(\cG_2-M(\zeta))^{-1}\in\cB(\mathfrak{H})$ by

\[(\cG_2-M(\zeta))^{-1}:x\mapsto z:=\Gamma_0R_2(\zeta)f\]

with $f$ from (<ref>). Then (<ref>) yields

\begin{equation*}\begin{split}
&\gr\big((\cG_2-M(\zeta))^{-1}\big)=\big\{(x,z): x\in\mathfrak{H}\big\}
\\&\quad=\big\{ \big(\Gamma_1R_2(\zeta)f-M(\zeta)\Gamma_0R_2(\zeta)f,\Gamma_0R_2(\zeta)f\big): x\in\mathfrak{H}\big\}=(\cG_2-M(\zeta))^{-1},
\end{split}\end{equation*}

thus finishing the proof of (<ref>).

We are ready to derive (<ref>). Combining (<ref>), $\Gamma_0\hatt R_1(\zeta)=0$, and (<ref>) we get

\begin{align}
\begin{split}
\lb{1.25new}
R_{1}(\zeta)- R_{2}(\zeta)&= (\Gamma_0\hatt R_1(\overline{\zeta}))^*\Gamma_1 R_2(\zeta)-(\Gamma_1\hatt R_1(\overline{\zeta}))^*\Gamma_0 R_2(\zeta)\\
&=-\gamma(\zeta)\Gamma_0 R_2(\zeta).
\end{split}
\end{align}

Assertions (<ref>) and (<ref>)  imply

\begin{equation}
\Gamma_0R_2(\zeta)f=(\cG_2-M(\zeta))^{-1}\Gamma_1R_1(\zeta)f=(\cG_2-M(\zeta))^{-1}\gamma^*(\overline{\zeta})f.
\end{equation}

Plugging this into (<ref>) yields (<ref>).

§.§ Old Lagrangian material

Selim, I put the old Lagrangian parametrization here as
I did not want just to delete it

It is well known that self-adjoint extensions of $A$ can be parameterized by Lagrangian planes.  Such parameterization depends on the choice of the trace operator $\tr$ and the “boundary" space $\mathfrak{H}$, see, e.g., <cit.> and <cit.>. The following theorem is yet another variant of the parameterization.

Let $\omega_1, \omega_2$ be  bounded, sesquilinear, skew-adjoint forms defined on Hilbert spaces $H_1, H_2$ correspondingly. Assume that $T:H_1\rightarrow H_2$ is a closed linear operator with dense domain,  dense range, and such that

\begin{equation}\lb{sympl}
\omega_1(u,v)=\omega_2(Tu, Tv), u,v\in\dom(T).
\end{equation}

Finally, fix a subspace $\cG\subset \dom(T)$. Then the subspace $\cG+\ker(T)$ is Lagrangian in $(\omega_1, \cH_1)$ if and only if $T(\cG)$ is Lagrangian in  $(\omega_2, \cH_2)$.

Assume that  $T(\cG)$ is Lagrangian in $(\omega_2, H_2)$. Then by (<ref>), the subspace $\cG+\ker(T)$ is isotropic in $(\omega_1, H_1)$. In order to prove that it is maximal isotopic, let us pick $v\in (\cG+\ker(T))^{\circ}$, i.e. $\omega_1(v, u)=0$ for all $u\in(\cG+\ker(T))$. Since $\dom(T)$ is dense in $\cH_1$ there exists a sequence $\{v_n\}_{n=1}^{\infty}\subset \dom(T)$ with $v_n\rightarrow v$ in $H_1$. For such a sequence we have

\begin{equation}\lb{omcon}
\omega_1(v_n, u)\rightarrow \omega_1(v, u)=0, u\in \cG.
\end{equation}

Let $P_2$ be the orthogonal projection onto  $T(\cG)$. Then, since $T(\cG)$ is Lagrangian, we have

\begin{equation}\omega_2(Tv_n, T u)=\omega_2((I_{H_2}-P_2)Tv_n, T u)=\langle(I_{H_2}-P_2)Tv_n, J_2T u\rangle_{H_2}\rightarrow 0, n\rightarrow\infty,   u\in \cG.
\end{equation}

Hence, $\{(I_{H_2}-P_2)Tv_n\}_{n=1}^{\infty}$ converges weakly to $0$ in the Hilbert space $\cJ_2T(\cG)$. Next, since the operator  $(I_{H_2}-P_2)T$ is closed (as a product of a bounded and closed ones) and since  $v_n\rightarrow v$, $(I_{H_2}-P_2)Tv_n\rightharpoonup 0$, we get $w\in\dom(T)$. Therefore (<ref>) yields $Tv\in (T(\cG))^{\circ}$. This gives $Tv\in T(\cG)$ and $(v-v_0)\in \cG$ for some $v_0\in \ker(T)$. Hence $v\in\cG+\ker(T)$, which shows the “if" part.

Assume that $\cG+\ker(T)$ is Lagrangian in $\cH_1$. Then by (<ref>) the subspace $T(\cG)\subset H_2$ is isotropic. In order to prove that it is maximal isotopic, let us pick $v\in (T(\cG))^{\circ}$, i.e. $\omega_2(v, Tu)=0$ for all $u\in(\cG+\ker(T))$. Since $\ran(T)$ is dense in $\cH_2$ there exists a sequence $\{u_n\}_{n=1}^{\infty}\subset \dom(T)$ with $Tu_n\rightarrow v$ in $H_2$. For such a sequence we have

\begin{equation}\lb{omcon2}
\omega_2(Tu_n, Tu)\rightarrow \omega_2(v, Tu)=0, u\in \cG.
\end{equation}

Let $P_1$ be the orthogonal projection onto  $\cG+\ker(T)$. Then, since $\cG+\ker(T)$ is Lagrangian, we have

\begin{equation}
\omega_1(u_n, u)=\omega_1((I_{H_1}-P_1)u_n,  u)=\langle(I_{H_1}-P_1)u_n, J_1 u\rangle_{H_1}\rightarrow 0, n\rightarrow\infty,   u\in \cG.
\end{equation}

Hence, $\{(I_{H_1}-P_1)u_n\}_{n=1}^{\infty}$ converges weakly to $0$ in the Hilbert space $\cJ_1(\cG+\ker(T))$. We observe that $(I_{H_1}-P_1)u_n \in \dom(T)$ and $T(I_{H_1}-P_1)u_n=(I_{H_1}-P_1)u_n$ since

thm1.4 Under Hypothesis <ref>, the self-adjoint extensions of $A$ whose domains are contained in $\cD=\dom(\tr)$ are in one-to-one correspondence with the Lagrangian planes in $\bndra$. Namely, the following two assertions hold.

* thm1.4part1 Let $\mathscr D\subset\cD$. If $\cA_{\mathscr D}$ is a self-adjoint extension of $A$ with $
\dom(\cA_{\mathscr D}):=\mathscr D$ then

\begin{equation}
\cG:=\overline{\tr(\mathscr D)}^{\bndra}\in\Lambda(\bndra).
\end{equation}

* thm1.4part2
Let $\cG\in\Lambda(\bndra)$ and $\tr^{-1}(\cG)=\{u\in\cD: \tr u\in\cG\}$ denote the preimage of $\cG$. Then the linear operator $\cA_{\tr^{-1}(\cG)}$ acting in $\cH$ and given by

\begin{equation}
\cA_{\tr^{-1}(\cG)}u:=A^* u \text{ for }  u\in \dom\left(\cA_{\tr^{-1}(\cG)}\right):={\tr^{-1}(\cG)}
\end{equation}

is essentially self-adjoint.

Part (<ref>). In order to show that $\cG$ is isotropic we employ the Green identity (<ref>): For arbitrary  $u,v \in\mathscr{D}$ one has

\begin{align}
\begin{split}
\omega\left(\tr u, \tr v \right) =\langle A^* u, v\rangle_{\cH}- \langle u,A^* v\rangle_{\cH}=0
\end{split}
\end{align}

where the latter equality follows since $\cA_{\mathscr D}$ is self-adjoint. Next, we show
$\cG ^{\circ}\subset\cG$.
To this end, we establish an intermediate inclusion

\begin{equation}\lb{b3}
\cG ^{\circ}\cap\tr (\cD )\subset \cG.
\end{equation}

For $u_0\in \cD$ such that $\tr u_0\in \cG ^{\circ}\cap\tr (\cD )$ one has
$\omega\left(\tr u_0, \tr v\right)=0$ for all $v\in\mathscr D$.
Then by Green's identity we have

\begin{equation}\lb{b8}
\langle A^* u_0, v\rangle_{\cH}- \langle u_0,\cA_{\mathscr D} v\rangle_{\cH}=\omega\left(\tr u_0,\tr v \right)=0 \text{ for all  $v\in\mathscr D$.}
\end{equation}

Hence, $\langle\cA^* u_0, v\rangle_{\cH}=\langle u_0,\cA_{\mathscr D} v\rangle_{\cH}$ for all $v\in\mathscr D.$ Thus $u_0\in \dom\left(\cA_{\mathscr D}^*\right)=\dom (\cA_{\mathscr D})$
since $\cA_{\mathscr D}$ is self-adjoint by the assumption. Consequently $\tr u_0\in\cG$ and (<ref>) holds.
Finally, employing (<ref>) and the standard properties of annihilator, we obtain
$\cG ^{\circ}\subset \overline{{\cG}+\left\{\tr(\cD )\right\}^{\circ}}^{\bndra}$.
This inclusion together with $\left\{\tr(\cD )\right\}^{\circ}=\{0\}$ and the fact that $\cG$ is isotropic yield $\cG^{\circ}=\overline{\cG}$ as asserted.

Part <ref>.
Following  <cit.> we introduce the space of abstract boundary values $\cH_{A}:= \dom(A^*)/\dom(A)$ equipped with the norm

\begin{equation}
\|[u]\|_{\cH_{A}}:=\inf\{\|u+w\|_{\cH_+}: w\in\dom(A) \},
\end{equation}

where $[u]$ denotes the equivalence class of $u\in\dom(A^*)$ and $\|\cdot\|_{\cH_+}$ is the graph norm of $A^*$. We define the symplectic form on $\cH_{A}$ by the formula

\begin{align}
\begin{split}
\wti \omega([u]&,[v]):=\langle A^*u, v \rangle_{\cH}-\langle u,A^*v
\rangle_{\cH},\ [u],[v]\in\dom(A^*)/\dom(A)\lb{ee2}.
\end{split}
\end{align}

Let $\mathscr D_{\cG}:=\tr^{-1}(\cG)$. By <cit.>,

Selim, is there a way to make a direct proof without going to 	<cit.>?

it suffices to show that the closure of the subspace
[{\mathscr D_{\cG}}]:=\{[u]: u\in{\mathscr D_{\cG}}\}
is Lagrangian in $\cH_{A}$ with respect to $\wti \omega$. Green's identity yields
[{\mathscr D_{\cG}}]\subset [{\mathscr D_{\cG}}]^{\circ}.
Next, we show that $[{\mathscr D_{\cG}}]^{\circ}\subset\overline{[{\mathscr D_{\cG}}]}$. First, we will obtain an auxiliary inclusion

\begin{equation}\lb{c3}
[{\mathscr D_{\cG}}]^{\circ}\cap [\cD ]\subset [{\mathscr D_{\cG}}].
\end{equation}

To that end, we notice that if  $[u_0]\in[{\mathscr D_{\cG}}]^{\circ}\cap [\cD ]$ then
$\wti\omega\left([u_0],[v]\right)=0 \text{\, for all  $[v]\in [{\mathscr D_{\cG}}]$}$,
that is,
$0=\langle A^* u_0, v \rangle_{\cH}-\langle u_0,A^* v \rangle_{\cH}$ for all  $[v]\in [{\mathscr D_{\cG}}]$.
Then by Green's identity one has

\begin{align}
\begin{split}
&\langle A^* u_0, v \rangle_{\cH}-\langle u_0,A^* v \rangle_{\cH}=\omega(\tr u_0, \tr v)=0 \text{ for all $[v]\in[\mathscr D_\cG]$}.
\end{split}
\end{align}

Therefore $\tr u_0\in(\cG\cap \tr(\cD))^{\circ}$.
Since $\cG$ is Lagrangian, we have

\begin{equation}\lb{c7}
\big(\cG\cap\tr\cD\big)^{\circ}=\overline{\cG+(\tr\cD)^{\circ}}=\overline{\cG+\{0\}}=\cG.
\end{equation}

Therefore $\tr u_0\in \cG$ and (<ref>) holds as asserted. Finally, using (<ref>) we infer

\begin{equation}\lb{c8}
\left[{\mathscr D_{\cG}}\right]^{\circ}\subset \overline{\overline{[{\mathscr D_{\cG}}]}+ \left[\cD\right]^{\circ}}=\overline{[{\mathscr D_{\cG}}]},
\end{equation}

where we used density of $\cD$.

We should reword next paragraph and include it in the Introduction

The spectral flow  of the family $\{B_x\}_{x=-\infty}^\infty$ of the self-adjoint operators $B_x$ is of interest as it is related to the spectrum of the Schröringer operator $-\Delta_{(x,y)}+V$ in $L^2(\Omega; \bbC^N)$. We conjecture that this relation could be established using spacial dynamics, cf. [70, 92, 91], via a connection to a first order differential operator, cf. [74] and [50].
We will just sketch possible arguments to describe the connection.
Indeed, the second order in $x$ operator $-\partial_x^2+B_x$ is Fredholm in $L^2(\Omega; \bbC^N)$ if and only if the first order in $x$ operator $-\partial_x+\widehat{B}_x$, with $\widehat{B}_x:=\begin{bmatrix}0&I\\ B_x&0\end{bmatrix}$, is Fredholm
in $L^2(\bbR; L^2({\mathbb B}_{r(x)}^n; \bbC^{2N}))$, and their Fredholm indices are equal, see, e.g., [92] and <cit.>. In turn, by a generalization of the celebrated Palmer Theorem,  see, e.g., [70, 92, 91], the Fredholm index of the operator $-\partial_x+\widehat{B}_x$ is equal to the spectral flow through the imaginary axis of the operator family
$\{\widehat{B}_x\}_{x=-\infty}^\infty$ or, equivalently,  to the so called relative Morse index $\dim(E_-^s\cap E_+^u)-\dim(E_-^u\cap E_+^s)$,
cf. [92] and <cit.>, <cit.>; here, we denote by $E_\pm^s$, respectively,
$E_\pm^u$ the spectral subspace of the asymptotic operator  $\widehat{B}_{\pm\infty}$ that corresponds to the part of its spectrum located in the left, respectively, right half-plane. Since $\spec(B_x)=(\spec(\widehat{B}_x))^2$, the spectral flow through the imaginary axis of the operator family
$\{\widehat{B}_x\}_{x=-\infty}^\infty$ is equal to the spectral flow through zero of the operator family
$\{{B}_x\}_{x=-\infty}^\infty$, as required.

The foundational for this theory finite dimensional results can be found, e.g., in Theorem II.5.4 and Theorem II.6.8 of [65] and sometimes called in a related Maslov index literature, cf. <cit.>, the Kato selection theorem as they allow one to properly choose the $m$ branches of the eigenvalue curves for $H_t$ and compute their slopes.

The fundamental relation between the spectral flow, a net number of eigenvalues crossing a fixed point, and the Maslov index, a topological invariant of paths of Lagrangian planes, relies on an explicit asymptotic expansion of eigenvalues near the crossing point and the fact that the first order coefficient in the expansion is precisely given by the Maslov crossing form. Yet another source of inspiration for this work is rooted in the classical Rayleigh–Hadamard–Rellich type formulas for the derivatives of eigenvalues curves for elliptic operators posed of a $t$-dependent family of multidimensional domains which are central in geometric spectral theory of elliptic operators on bounded domains, cf.

This paper can be viewed as a part of a wider effort to build full-scale perturbation theory for the case when one perturbs not an operator itself but also its domain. On one hand, this theory is rooted in the classical extension theory of self adjoint operators that goes back to Krein, Birman, Vishik, cf., e.g., [41, 56, 93], and has been an exceptionally active area of research
[1, 7, 13, 14, 23, 42, 56, 76, 80] culminating in the comprehensive monograph [12]. On the other hand it is related to the classical variational formulas for eigenvalues that go back to Rayleigh, Hadamard and Rellich [60, 83, 86], see also  <cit.> for a recent exposition. Thirdly, it is related to the progress in our understanding of the relation between the Malsov index and the spectral flow formulas, yet another classical subject originated in Atiayh-Padoti-Singer Theory and recently developed by a number of authors [19, 20, 21, 22, 26, 35, 36, 48, 71, 72, 73, 88] [add more].

One of the sources of inspiration for the current work is in the subject of variational formulas for partial differential equations on variable domains. The study of dependence on domains of  the eigenvalues, eigenfunctions, resolvent, etc., for partial differential operators has a distinguished history that goes back to the classical variational formulas for the derivative of the eigenvalues of the Laplacian restricted to a one-parametric family of domains, see [60, 83, 86] and <cit.>.
One can find further information and literature in the book [62], the review papers <cit.> and more recent sources <cit.> and consult <cit.> for the analytic theory for boundary perturbations. Of course, the literature on this classical subject is enormous and we do not make any attempts to review it. A typical example that we have in mind is furnished by certain second order elliptic operators, $\cL$, on a star-shaped domain $\Omega\subset\bbR^n$ with, say, Robin-type boundary conditions. Linear rescaling $x\mapsto tx$, $x\in\Omega$, $t\in(0,1]$, produces a family of shrinking
domains $\Omega_t\subseteq\Omega$. One restricts $\cL$ to $\Omega_t$ and then rescales the restricted operator back to $\Omega$ thus obtaining a family $\{\cL_t\}_{t\in(0,1]}$ of operators on $\Omega$. We emphasize that both the coefficients of $\cL_t$ and the boundary conditions will depend on $t$ as the rescaling forces us to multiply every derivative by $t^{-1}$. Thus, one needs to study the eigenvalues, eigenfunctions, and the resolvent operators of $\cL_t$ as functions of $t$, and the results of the current papers are designed to cover this class of problems. As an illustration, in Subsection <ref> we derive from the general results of the current paper the analogs of the classical Hadamard-Rayleigh-Rellich formula for the derivative of the eigenvalues of the Schrödinger operator with matrix-valued potential on the family of star-shaped domains. Also, we refer the reader to the papers [35, 36, 73] where these formulas appeared in connections with the Maslov index.

A ridiculously simple example of what we do in this paper is furnished by the one-dimensional Laplace operators $-\partial^2_D$ and $-\partial^2_N$ in $L^2([0,1])$ equipped with the Dirichlet and Neumann boundary conditions. The difference of the two operators on the intersection of their domains is zero; thus neither of them is an additive perturbation of another. At the same time, their resolvent operators
are: Indeed, Krein's type formula for the difference of the resolvents proved below gives ....

finish the calculation

We refer to [31, 32] for an illuminating discussion of this subject. Going further, we consider a family of the Laplacians, $-\partial^2_t$, $t\in[0,1]$, equipped with the boundary condition  $\big(t\gaD+(1-t)\gaN\big)u=0$ such that $-\partial^2_0=-\partial^2_D$ and $-\partial^2_1=-\partial^2_N$. As proved below, the resolvent operator of $-\partial^2_t$ satisfies the following Riccati-type differential equation ... while the derivatives of the eigenvalues $\lambda_j(t)$ are given by ....

The classical abstract Krein formula concerns the situation of a fixed self-adjoint extension $\cA_1:=A^*\upharpoonright_{\ker(\Gamma_0)}$ and arbitrary self-adjoint extesion $\cA_2$, defined by a Lagrangian plane $\cG_2$, and the difference of resolvents is computed as follows
\begin{equation}\lb{int3}
\text{ for $\zeta\not\in\Sp(\cA_1)\cup \Sp(\cA_2)$},
\end{equation}

where $\gamma(\zeta):= (\Gamma_0\upharpoonright_{\ker(A^*-\zeta)})^{-1}$ and the abstract Weyl function is defined by $M(\zeta):=\Gamma_1\gamma(\zeta)$, see <cit.>. In the PDE context, the Krein formulas has been derived for pairwise differences of Dirichlet (i.e. Friedrichs extensions), Neumann, and Robin realizations of elliptic operators, e.g., for the Robin Laplacian $-\Delta_{\Theta}$ acting in $L^2(\Omega)$ defined by the boundary condition $\wti\gaN u+\Theta \gaD u=0$ and Dirichlet Laplacian $-\Delta_D$, one has

\begin{equation}\lb{RtoD}
(-\Delta_{\Theta}-\zeta)^{-1}-(-\Delta_{D}-\zeta)^{-1}=[\wti \gaN(-\Delta_{\Theta}-\zeta)^{-1}]^*M_{\Theta, D}(\zeta)[\wti \gaN(-\Delta_{D}-\zeta)^{-1}]
\end{equation}

where $\wti \gaN$, $\gaD$ are versions of the Neumann and Dirichlet traces, and $M_{\Theta, D}$ is a Robin-to-Dirichlet map, see [52].


11institutetext: University of Western Macedonia, Kozani 50100, Greece
11email<EMAIL_ADDRESS>11email<EMAIL_ADDRESS>
# Integer programming models for the semi-obnoxious p-median problem

Nikolaos Ploskas 11 0000-0001-5876-9945 Kostas Stergiou 11 0000-0002-5702-9096

###### Abstract

The p-median problem concerns the location of facilities so that the sum of
distances between the demand points and their nearest facility is minimized.
We study a variant of this classic location problem where minimum distance
constraints exist both between the facilities and between the facilities and
the demand points. This specific type of problem can be used to model
situations where the facilities to be located are semi-obnoxious. But despite
its relevance to real life scenarios, it has received little attention within
the vast literature on location problems. We present twelve ILP models for
this problem, coupling three formulations of the p-median problem with four
formulations of the distance constraints. We utilize Gurobi Optimizer v9.0.3
in order to compare these ILP models on a large dataset of problems.
Experimental results demonstrate that the classic p-median model proposed by
ReVelle & Swain and the model proposed by Rosing et al. are the best
performers.

###### Keywords:

Facility location p-median problem Integer optimization Distance constraints.

##  1 Introduction

Facility location problems are among the most widely studied problems in OR,
AI, computational geometry, and other disciplines. Most such problems include
“pull” objectives, where clients wish to have the facility sites (e.g.,
pharmacies and stores) located close to them. For instance, the p-median
problem concerns the location of p facilities so that the sum of (weighted)
distances between the demand points and their nearest facility is minimized.

Another widely studied class includes problems with “push” objectives, where
we seek to place the facility sites away from demand points, and in some cases
from each other. Such facilities are known as obnoxious (or undesirable)
because they have hazardous properties (i.e., dump sites for chemicals,
prisons, etc.). In between, we have the class of semi-obnoxious (or semi-
desirable) location problems where the facilities have both desirable and
undesirable properties, meaning that both pull and push objectives are present
[16]. In this case, clients wish to have facilities located close to them, but
not too close. For example, the residents of a suburb wish to have restaurants
and bars situated near them, but not too close because of the noise and
traffic associated with such facilities.

One way to model the requirements that arise when trying to locate semi-
obnoxious facilities is by placing minimum allowed distance constraints
between facilities and/or between facilities and demand points [16]. For
instance, gas stations (which are typical semi-obnoxious facilities) need to
be at a minimum distance from one another and also from residencies, for
safety reasons. At the same time, they must not be placed too far away from
residencies for practical reasons. Such a scenario can be captured as a
p-median problem with distance constraints.

The p-median problem, which is $\mathcal{NP}$-hard on general networks for an
arbitrary p, was originally proposed by Hakimi [14, 15]. In 1970 ReVelle &
Swain [19] presented the first IP for the p-median problem. They utilized a
structure proposed by Balinski [1] in a plant location problem. Despite the
large number of variables and constraints of this problem, it has many
integer-friendly properties.

In 1979 Rosing et al. [20] proposed a formulation by using fewer Balinski
constraints and adding a constraint used by Efroymson & Ray [13] in a plant
location problem. This model reduces the number of constraints needed while
retaining many of the integer-friendly properties of the ReVelle & Swain
formulation. In 2003 Church [6] proposed a new formulation of the p-median
problem by building upon the Rosing et al. model and reducing variables based
upon a property of “equivalent assignment conditions”. There are also other
works that propose reformulations of the classic p-median problem [9, 10]. For
a thorough literature review on formulations of the p-median problem, see [6].

Location problems with distance constraints have received little attention
within the vast literature on location problems in OR and related areas. There
are some early works that considered maximum distance constraints between the
demand nodes and the facility locations [7, 8, 21]. Moon and Chaudhry were the
first to systematically study location problems with distance constraints
[17]. Among the problems studied was the Median$/l/l$, in their terminology,
where we seek to minimize the sum of the distances between the clients and
their closest facility (as in p-median) subject to constraints that impose
lower bounds (hence the $l$) in the distance between any two facilities and
also between any facility and any client. This is the problem that we revisit
in this paper, simply calling it p-median with distance constraints.

This specific problem has not received attention since then, but other
location problems with distance constraints have been studied [5, 18, 3, 11,
12]. In a work that is relevant to ours, Berman & Huang [3] presented various
formulations for the minimum weighted covering location problem aiming to
minimize the total demand covered subject to the condition that no two
facilities are allowed to be closer than a pre-specified distance.
Interestingly, they were the first to experimentally compare four,
increasingly more elaborate, models for the distance constraints in a location
problem.

In this work we make a detailed study of the p-median problem with distance
constraints for the first time. Given that the ILP modeling of this problem is
based on two decisions: how to model the p-median objective and how to model
the distance constraints, we develop and compare 12 different formulations by
combining the three main formulations of the p-median problem with the four
formulations of the distance constraints presented by Berman & Huang [3],
extending them to include minimum distance constraints between facilities and
clients.

To compare these ILP models, we utilize Gurobi Optimizer v9.0.3 and apply the
solver on a set of instances generated using as basis the standard p-median
benchmark dataset [2]. The experimental evaluation shows that the classic
p-median model proposed by ReVelle & Swain and the model proposed by Church
are the best performers. Regarding the distance constraint formulations, there
is no clear winner. An interesting conjencture that can be made based on our
results is that simple models that were considered inefficient in the past are
quite competitive on our studied problem due to the advancements in ILP solver
technology.

##  2 Problem definition

We assume that p homogeneous facilities are to be placed on a weighted network
$G=\\{V,E\\}$ embedded on a plane, where a set I of demand points (clients)
are already located. We also assume that the set J of nodes in $G$ where
facilities can potentially be located is known. Hence, we deal with a
discrete/network location problem. The weight $w_{ij}$ on an edge $(i,j)\in E$
denotes the symmetric service cost (typically the distance) between nodes $i$
and $j$.

Between each pair of facilities $i$ and $j$ there is a distance constraint
$Euc_{d}(i,j)>d_{1}$ specifying that the Euclidean distance between the nodes
where the facilities $i$ and $j$ are located must be greater than $d_{1}$,
where $d_{1}$ is a real constant. Also, between each facility $i$ and client
$k$ there is a distance constraint $Euc_{d}(i,k)>d_{2}$ specifying that the
Euclidean distance between the node where facility $i$ is located and node
where client $k$ is located must be greater than $d_{2}$, where $d_{2}$ is a
real constant. The constants $d_{1}$ and $d_{2}$ are common for all facilities
and clients since we consider the homogeneous case in this work.

At this point let us clarify that the Euclidean distance $Euc_{d}(i,j)$
between two nodes $i$ and $j$ should not be confused with the minimum distance
between the two nodes in the network, which is given by the sum of weights of
the edges that belong to the shortest path between $i$ and $j$. It is natural
to model the costs of service between $i$ and $j$ through the cost of the
paths between the two nodes, but we believe that any distance constraint
established due to, say, safety reasons, should consider a metric like the
Euclidean distance instead.

The goal is to minimize the sum of distances between each demand point and its
nearest facility subject to the satisfaction of all the distance constraints.
The existence of distance constraints between clients and facility locations
means that, in contrast to the standard p-median problem where all nodes in
the network are clients, in our case the set of potential facility sites is
discrete from the set of client nodes.

##  3 ILP models

In this section, we present various binary ILP models for solving the
homogeneous p-median problem with distance constraints so that the sum of
distances between the demand points and their nearest facility is minimized.
The following notation is used:

  * •

$I$: set of demand nodes.

  * •

$J$: set of candidate facility sites.

  * •

$p$: the number of facilities to be located.

  * •

$s_{ij}$: the shortest distance between any two nodes (demand nodes or
facility sites).

  * •

$d_{1}$: the minimum distance between each pair of facilities
$\left(i,j\right)$, where $i\in J$ and $j\in J$.

  * •

$d_{2}$: the minimum distance between each facility $i\in J$ and all demand
nodes.

  * •

$x_{j}=1$ if a facility is located at facility site $j\in J$ and 0 otherwise.

  * •

$y_{ij}=1$ if a demand node $i\in I$ is assigned to facility site $j\in J$ and
0 otherwise.

In subsection 3.1 we present the three p-median formulations used in this
work, while subsection 3.2 includes the four distance constraint models.

###  3.1 p-median models

We can formulate the p-median model proposed by ReVelle & Swain [19] as:

$\displaystyle\min$ $\displaystyle\sum\limits_{i\in I}\sum\limits_{j\in
J}s_{ij}y_{ij}$ (1) $\displaystyle\mathop{\rm s.t.}$
$\displaystyle\sum\limits_{j\in J}x_{j}=p$ (2) $\displaystyle\sum\limits_{j\in
J}y_{ij}=1$ $\displaystyle\forall i\in I$ (3) $\displaystyle y_{ij}\leq x_{j}$
$\displaystyle\forall i\in I,\forall j\in J$ (4) $\displaystyle
x_{j}\in\left\\{0,1\right\\}$ $\displaystyle\forall j\in J$ (5) $\displaystyle
y_{ij}\in\left\\{0,1\right\\}$ $\displaystyle\forall i\in I,\forall j\in J$
(6)

The objective function 1 aims at minimizing the sum of the shortest distances
between the demand points and their nearest open facility. Constraint 2
specifies that $p$ facilities are to be opened, while Constraint 3 ensures
that each demand node will be served by one facility site. Constraint 4, i.e.,
the Balinski constraints, guarantees that each demand node will be served by a
facility site that is opened.

This model has $\left|J\right|+\left|I\right|\times\left|J\right|$ variables
and $\left|I\right|+\left|I\right|\times\left|J\right|+1$ constraints.
Although it was almost impossible to solve problems of large size using this
model in the past, ILP solvers have become significantly faster and capable of
solving very large problems in the last two decades.

Rosing et al. [20] noticed that the Balinski constraints make the constraint
set too large and they could use the Efroymson & Ray constraints instead:

$\sum\limits_{i\in I}y_{ij}\leq nx_{j},\forall j\in J$ (7)

This constraint form adds only $\left|I\right|$ constraints to the model,
whereas the Balinski form adds $\left|I\right|\times\left|J\right|$
constraints. However, the Efroymson & Ray constraints are not integer-
friendly, requiring many time consuming branch and bound operations. Rosing et
al. formulated a model which balances between the use of Balinski constraints
and Efroymson & Ray constraints with the aim to generate a smaller integer-
friendly model. The proposed model is the following:

$\displaystyle\min$ $\displaystyle\sum\limits_{i\in I}\sum\limits_{j\in
J}s_{ij}y_{ij}$ (8) $\displaystyle\mathop{\rm s.t.}$
$\displaystyle\sum\limits_{j\in J}x_{j}=p$ (9) $\displaystyle\sum\limits_{j\in
F_{i}}y_{ij}=1$ $\displaystyle\forall i\in I$ (10) $\displaystyle y_{ij}\leq
x_{j}$ $\displaystyle\forall i\in I,\forall j\in K_{ir},\text{ where }i\neq j$
(11) $\displaystyle\sum\limits_{i\in I}y_{ij}\leq nx_{j}$
$\displaystyle\forall j\in J$ (12) $\displaystyle
x_{j}\in\left\\{0,1\right\\}$ $\displaystyle\forall j\in J$ (13)
$\displaystyle y_{ij}\in\left\\{0,1\right\\}$ $\displaystyle\forall i\in
I,\forall j\in J$ (14)

where $r$ is the Balinski constraint rank-cutoff level, $F_{i}$ is the set of
all facility sites except the $p-1$ furthest sites from client $i$, and
$K_{ir}$ is the set of the r-closest sites to client $i$.

This model is the same as the original formulation except that: (i) Efroymson
& Ray constraints 12 have been included for each facility site, (ii) Balinski
constraints 11 are included for only the r-closest facility sites to each
client, and (iii) $p-1$ assignments to each client are eliminated. This model
has $\left|J\right|+\left|I\right|\times\left(\left|J\right|-p+1\right)$
variables and $2\left|I\right|+r\left|I\right|+1$ constraints.

Church [6] built upon Rosing et al. model and introduced a new formulation by
reducing variables based upon a property of “equivalent assignment
conditions”. The proposed model is the following:

$\displaystyle\min$ $\displaystyle\sum\limits_{\left(i,j\right)\in
A}C_{ij}y_{ij}$ (15) $\displaystyle\mathop{\rm s.t.}$
$\displaystyle\sum\limits_{j\in J}x_{j}=p$ (16)
$\displaystyle\sum\limits_{k\leq\left|J\right|-p+1}\sum\limits_{\left(l,j\right)\in
M_{ik}}y_{lj}=1$ $\displaystyle\forall i\in I$ (17)
$\displaystyle\sum\limits_{i\in
A_{j}}y_{ij}\leq\left(\left|J\right|-p+1\right)x_{j}$ $\displaystyle\forall
j\in J$ (18) $\displaystyle y_{ij}\leq x_{j}$ $\displaystyle\forall i\text{
and }j\in\left\\{O_{ir}\right\\},\text{ where }i\neq j$ (19) $\displaystyle
x_{j}\in\left\\{0,1\right\\}$ $\displaystyle\forall j\in J$ (20)
$\displaystyle y_{ij}\in\left\\{0,1\right\\}$ $\displaystyle\forall i\in
I,\forall j\in J$ (21)

where:

  * •

$k$: index that refers to the closest facility site order.

  * •

$O_{ik}$: index of the $k\in J$ closest facility site to a client $i\in I$.

  * •

$\left\\{O_{ik}\right\\}=\left\\{O_{i1},O_{i2},\cdots,O_{ik}\right\\}$: the
set of facility sites that comprise the $k\in J$ closest facility sites to
client $i\in I$.

  * •

$a_{iklj}=\left\\{\begin{array}[]{ll}1,&\text{if variable }x_{lj}\text{
represents the }k\in J\text{ closest assignment variable }\\\ &\text{for
client }i\in I\\\ 0,&\text{ otherwise}.\end{array}\right.$

  * •

$M_{ik}=\left\\{\left(l,j\right)\mid x_{lj}\text{ represents the }k\in J\text{
closest facility assignment for client }i\in I\right\\}$.

  * •

$C_{lj}=\sum\limits_{i\in I}\sum\limits_{k\in J}d_{ik}a_{iklj}$.

  * •

$A=\left\\{\left(i,j\right)\mid\text{ variable }x_{ij}\text{ is included in
the model}\right\\}$, after variable substitution and eliminating any
assignment to the furthest $p-1$ facility sites for a given client $i\in I$.

  * •

$A_{j}=\left\\{i\mid\left(i,j\right)\in A\right\\}$.

Each constraint in this model mirrors the associated constraint in the Rosing
et al. model. The additional notation used to formulate this model is
associated with the concept of equivalent variables and their substitution.
For every equivalent assignment condition present in the model, a variable is
eliminated and a constraint is eliminated when that variable was present in a
Balinski constraint. This model is no larger than the Rosing et al. model.

###  3.2 Distance constraints

We now present four formulations for the minimum distance constraints between
each pair facilities as used by Berman & Huang [3] for the minimum weighted
covering location problem. We modify these formulations and use them to model
the p-median problem with distance constraints. Let us introduce two new sets:

  * •

$M=\left\\{(i,j)\mid i,j\in J,Euc_{d}(i,j)\leq d_{1}\right\\},\forall i\in
J,\forall j\in J$: the set of pairs $(i,j)$ s.t. facilities cannot be in
facility sites $i$ and $j$ because $i$ and $j$ are not in a safe distance
between each other.

  * •

$N=\left\\{(j)\mid j\in J,\exists k\in I,Euc_{d}(j,k)\leq
d_{2}\right\\},\forall j\in J$: a facility cannot be placed in facility site
$j$ because there exists a demand node $k$ that is not in safe distance from
$j$.

Using this notation the first, and simplest, distance constraint that can be
added to the p-median models is the following:

$x_{i}+x_{j}\leq 1,\forall\left(i,j\right)\in M$ (22)

Constraint 22 states that two facilities cannot be established if the distance
between them is less than $d_{1}$.

Another possible formulation for the distance constraint between two
facilities is the following:

$p\left(1-x_{j}\right)\geq\sum\limits_{k\in Q_{j}}x_{k},\forall j\in M$ (23)

where $Q_{j}=\left\\{k\in M\mid Euc_{d}(k,j)\leq d_{1},k\neq j,\forall j\in
M\right\\}$.

If a facility site is opened at $j\in M$, $x_{k}$ must take a value of $0$ for
all $k\in Q_{j}$. When $x_{j}=0$, at most $p$ facilities can be sited in
$Q_{j}$. This constraint was originally proposed by Moon & Chaudhry [17] to
formulate the anti-covering location problem. Constraint set 23 has fewer
constraints than constraint set 22.

The third formulation for this distance constraint is the following:

$n_{j}\left(1-x_{j}\right)\geq\sum\limits_{k\in Q_{j}}x_{k},\forall j\in M$
(24)

where $n_{j}$ is the minimum coefficient necessary to impose spatial
restrictions between members of the set $Q_{j}$. In order to calculate $n_{j}$
a complementary problem needs to be solved. Let us denote the graph
$G^{\prime}=\left(N^{\prime},L^{\prime}\right)$, where $N^{\prime}=M$ and
$L^{\prime}=\left\\{\left(j,k\right)\in M\mid Euc_{d}(j,k)\leq d_{1},j\neq
k\right\\}$. Let $G^{{}^{\prime}}_{j}$ be the subgraph of $G^{{}^{\prime}}$
induced by $Q_{j}\cup{j}$ and $q_{j}$ is the cardinality of the maximum
independent set of $G^{{}^{\prime}}_{j}$. Although the maximum independent set
problem is $\mathcal{NP}$-complete, $G^{{}^{\prime}}_{j}$ is usually a small
graph. In this work we calculate $q_{j}$ with the second version of the Bron &
Kerbosch algorithm [4]. Once $q_{j}$ is calculated, the value of $n_{j}$ can
be determined by $n_{j}=\min\left(p,q_{j}\right)$. Constraint 24 is tighter
than constraint 23.

The fourth formulation for the distance constraint between two facility sites
is the following:

$\displaystyle\sum\limits_{k\in H_{j}}x_{k}\leq 1,\forall H_{j}\in L_{1}$ (25)
$\displaystyle n^{{}^{\prime}}_{j}\left(1-x_{j}\right)\geq\sum\limits_{k\in
Q^{{}^{\prime}}_{j}}x_{k},\forall Q^{{}^{\prime}}_{j}\in L_{2}$ (26)

where $H_{j}$ is a maximal clique which includes node $j$ of $G^{\prime}$,
$Q^{{}^{\prime}}_{j}$ is the subset of $Q_{j}\cup j$ after removing all
elements in $H_{j}$. $L_{1}$ is the set of all $H_{j}\left(j\in M\right)$ and
$L_{2}$ the set of all nonempty $Q^{{}^{\prime}}_{j}\left(j\in M\right)$. In
order to obtain $n^{{}^{\prime}}_{j}$ the following ILP problem needs to be
solved:

$\displaystyle\min$ $\displaystyle\sum\limits_{v\in Q^{{}^{\prime}}_{j}}x_{v}$
(27) $\displaystyle\mathop{\rm s.t.}$ $\displaystyle x_{u}+x_{v}\leq 1$
$\displaystyle\forall u,v\in Q^{{}^{\prime}}_{j}\text{ and
}\left(u,v\right)\in G^{\prime}$ (28) $\displaystyle
x_{v}\in\left\\{0,1\right\\}$ $\displaystyle\forall v\in Q^{{}^{\prime}}_{j}$
(29)

Since $Q^{{}^{\prime}}_{j}$ is usually small, the solution of the ILP model
needs only a very small amount of time compared to the solution of the
p-median problem.

Finally, we also consider a distance constraint between each facility and each
client, where all clients should be at a safe distance ($d_{2}$) for all
opened facilities. This constraint can be formulated as:

$x_{j}=0,\forall j\in N$ (30)

This constraint set ensures that each demand node is on a safe distance from
each facility site if it is opened. Instead of setting variables $x_{j}$,
where $\left(j\right)\in N$, equal to $0$, these variables can be removed from
the model.

Coupling the three formulations of the p-median problem presented in
subsection 3.1 with the four formulations of the distance constraints
presented in this subsection, we get twelve different formulations of the
p-median problem with distance constraints.

##  4 Computational study

As this is the first time that a computational study of the p-median problem
with distance constraints is performed, there is a lack of benchmarks for this
problem, and location problems with distance constraints in general. Hence, we
experimented with random instances. We first describe the random generation
model, then we give results from the twelve formulations. Computations were
performed on an Intel i7 CPU 8700 with 16 GB of main memory, a clock of 3.2
GHz, an L1 cache of 348 KB, an L2 cache of 2 MB, and an L3 cache of 12 MB,
running under CentOS 8.4. The ILP models were solved using Gurobi v9.0.3.

###  4.1 Random generation model

The random generator takes as input the p-median benchmark dataset [2], which
consists of $40$ problems with a set $V$ of $100$–$900$ nodes and a value of
$p$ between $5$ and $200$. We present results from twelve selected problems
that were reasonably hard for the presented formulations. In the classical
p-median problem, all nodes are considered to be clients. However, in our
case, there is a minimum distance constraint between facilities and clients.
Therefore, we select $\left|J\right|<\left|V\right|$ nodes to be candidate
facilities, while the remaining nodes
$\left|I\right|=\left|V\right|-\left|J\right|$ are clients. We have considered
two cases:

  1. 1.

We start by randomly selecting $80$% of the nodes as candidate facility sites
and the remaining $20$% as clients. In case the resulting candidate facility
sites is less than or equal to $p$, then we progressively increase the number
of candidate facility sites until $\left|J\right|>p$. Also, in case the
generated instances are infeasible, we again increase the number of candidate
facilities (while decreasing the number of clients) until we get feasible
instances.

  2. 2.

We start by selecting $20$% of the nodes as candidate facility sites and the
remaining $80$% as clients. Similarly to the previous case, we progressively
increase the number of candidate facility sites until $\left|J\right|>p$ and
the generated instances are feasible.

In order to select the parameter $d_{1}$, we find the minimum and maximum
distance between all pairs of candidate facility sites and we set $d_{1}$
equal to a randomly generated number in the range
$\left[\min,\min+(\max-\min)/10\right]$. This interval generated feasible
problems in most cases. We perform a similar procedure to select the parameter
$d_{2}$. In this case, we calculate the minimum and maximum distance between
all candidate facility sites and all clients.

Finally, the p-median benchmark problems include only the distance for each
edge in the graph. We used Dijkstra’s algorithm to calculate the shortest path
$s_{ij},i\in V,j\in V$ between each pair of nodes. We set the Euclidean
distance between each pair of nodes to a random number in the interval
$\left[s_{ij}/2,s_{ij}\right]$.

###  4.2 Experimental results

Tables 2 \- 2 present the execution times of the twelve ILP formulations for
cases 1 and 2 described above. We use the following abbreviations: (i) RS1 -
RS4: the ReVelle & Swain formulation using the distance constraints in 22, 23,
24, and 25 \- 26, respectively, (ii) RRR1 - RRR4: the Rosing et al.
formulation using the constraints in 22, 23, 24, and 25 \- 26, respectively,
and (iii) CHU1 - CHU4: the Church formulation using the constraints in 22, 23,
24, and 25 \- 26, respectively. Results from the best method in each class are
given in bold.

All formulations were executed on ten instances for each of the problems.
Tables 2 \- 2 present the average execution times. As the Rosing et al. and
Church formulations require as input the parameter $r$, i.e., the Balinski
constraint rank-cutoff level, we ran these models for values of $r$ starting
from $10$ to $60$ with a step of $10$. As there is no value for $r$ that gives
the best results on all problem classes, Tables 2 \- 2 report the results for
$r=50$, which resulted in the best results on most classes. But in the two
rows giving the average execution times we also include, in parentheses, the
average execution times obtained by running the two models on each problem
class with the best $r$ value for the specific class.

The average results show that the ReVelle & Swain formulation is the best one
for problem of Case 1 (Table 2). This is contrary to the literature, since
computational studies on p-median formulations concluded that the ReVelle &
Swain cannot be used in practice due to its size [6]. However, current state-
of-the-art ILP solvers have become much faster and are able to solve very
large instances, using sophisticated preprocessing routines to greatly reduce
the size of the problems. The main reason for the superiority of the classic
p-median problems in the first case is that it performs much fewer simplex
iterations that the other p-median formulations.

When we customize the value of $r$ for each problem class, the Church
formulation is the best one on problems of Case 2 (Table 2). However, the
other two formulations are not far off. But given that the customization of
$r$ would require extensive experimentation, the results from the case where
$r$ is set to a good overall value (50 in this case) give a more accurate
picture.

Table 1: Comparison of ILP formulations for problems where
$\left|J\right|>\left|I\right|$.

Problem \ Model | RS1 | RS2 | RS3 | RS4 | RRR1 | RRR2 | RRR3 | RRR4 | CHU1 | CHU2 | CHU3 | CHU4
---|---|---|---|---|---|---|---|---|---|---|---|---
pmed21 | 9.11 | 9.52 | 10.64 | 10.32 | 10.59 | 9.92 | 11.03 | 8.51 | 9.98 | 9.03 | 11.10 | 8.61
pmed22 | 9.39 | 8.54 | 10.04 | 7.68 | 7.78 | 7.54 | 8.36 | 8.47 | 7.63 | 7.64 | 7.34 | 8.03
pmed23 | 0.90 | 1.89 | 2.17 | 1.04 | 0.61 | 1.00 | 1.07 | 0.67 | 0.68 | 1.17 | 1.03 | 0.53
pmed26 | 7.86 | 10.61 | 7.29 | 8.10 | 9.82 | 11.30 | 12.20 | 9.17 | 10.45 | 11.26 | 8.01 | 12.16
pmed27 | 6.86 | 7.92 | 6.27 | 6.71 | 5.35 | 8.18 | 6.98 | 4.52 | 6.29 | 5.72 | 5.95 | 5.60
pmed28 | 0.55 | 0.86 | 0.81 | 0.78 | 0.22 | 0.45 | 0.43 | 0.31 | 0.28 | 0.62 | 0.33 | 0.31
pmed31 | 55.19 | 52.44 | 50.56 | 53.23 | 60.94 | 63.10 | 64.65 | 80.31 | 60.14 | 56.92 | 55.25 | 54.36
pmed32 | 9.58 | 12.87 | 11.43 | 11.42 | 12.91 | 12.69 | 12.77 | 8.51 | 10.80 | 9.40 | 10.34 | 11.68
pmed33 | 1.25 | 2.20 | 2.88 | 1.95 | 0.63 | 1.47 | 1.13 | 0.78 | 0.75 | 1.37 | 1.41 | 0.81
pmed36 | 38.63 | 37.53 | 39.54 | 40.12 | 48.31 | 38.40 | 40.67 | 43.15 | 42.67 | 37.89 | 38.53 | 40.15
pmed37 | 2.99 | 5.61 | 5.58 | 3.64 | 1.85 | 3.76 | 4.63 | 2.85 | 2.21 | 4.44 | 3.62 | 2.70
pmed38 | 19.60 | 20.74 | 20.13 | 20.43 | 86.73 | 88.54 | 69.73 | 85.53 | 86.01 | 94.56 | 72.81 | 87.88
Average | 13.49 | 14.23 | 13.94 | 13.79 |  20.48 (18.50) |  20.53 (16.20) |  19.47 (17.27) |  21.07 (17.09) |  19.82 (15.38) |  20.00 (15.54) |  17.98 (15.94) |  19.40 (15.22)

Table 2: Comparison of ILP formulations for problems where
$\left|I\right|\geq\left|J\right|$.

Problem \ Model | RS1 | RS2 | RS3 | RS4 | RRR1 | RRR2 | RRR3 | RRR4 | CHU1 | CHU2 | CHU3 | CHU4
---|---|---|---|---|---|---|---|---|---|---|---|---
pmed21 | 4.56 | 5.93 | 4.58 | 4.29 | 4.36 | 6.25 | 6.24 | 6.20 | 4.75 | 4.28 | 4.35 | 5.88
pmed22 | 7.43 | 7.80 | 6.27 | 9.34 | 7.57 | 8.01 | 6.59 | 10.52 | 8.83 | 7.31 | 5.84 | 7.79
pmed23 | 0.38 | 0.44 | 0.38 | 0.36 | 0.38 | 0.36 | 0.36 | 0.39 | 0.38 | 0.34 | 0.31 | 0.32
pmed26 | 25.15 | 23.92 | 26.58 | 20.49 | 24.03 | 21.56 | 23.41 | 30.76 | 21.26 | 27.61 | 24.99 | 21.61
pmed27 | 17.74 | 20.47 | 23.43 | 17.33 | 25.29 | 21.49 | 24.59 | 24.48 | 18.36 | 24.41 | 24.44 | 28.77
pmed28 | 0.67 | 0.82 | 0.88 | 0.70 | 0.52 | 0.78 | 0.68 | 0.62 | 0.62 | 0.68 | 0.54 | 0.57
pmed31 | 55.19 | 52.44 | 50.56 | 53.23 | 60.94 | 63.10 | 64.65 | 80.31 | 60.14 | 56.92 | 55.25 | 54.36
pmed32 | 88.10 | 72.91 | 86.84 | 91.02 | 57.98 | 50.66 | 57.51 | 55.55 | 49.06 | 70.38 | 55.19 | 44.90
pmed33 | 1.00 | 1.19 | 1.13 | 0.99 | 0.75 | 1.09 | 1.14 | 0.73 | 0.68 | 1.00 | 0.89 | 0.73
pmed36 | 51.52 | 60.16 | 42.44 | 54.91 | 28.63 | 30.34 | 26.89 | 25.44 | 31.45 | 24.32 | 24.73 | 30.26
pmed37 | 1.69 | 1.91 | 2.34 | 1.63 | 1.53 | 1.95 | 1.62 | 1.21 | 1.21 | 1.39 | 1.38 | 1.26
pmed38 | 146.40 | 144.00 | 111.81 | 180.86 | 149.10 | 116.81 | 165.18 | 168.31 | 170.86 | 172.81 | 139.64 | 168.68
Average | 33.32 | 32.67 | 29.77 | 36.26 |  30.09 (24.72) |  26.87 (23.41) |  31.57 (25.94) |  33.71 (24.95) |  30.63 (23.32) |  32.62 (24.63) |  28.13 (23.13) |  30.43 (26.89)

Case 2 problem instances are harder on average and the Rosing et al.
formulation is the best performer in this case. Case 2 instances is harder
than Case 1 because they are much larger in size (both in terms of variables
and constraints), thus leading to more explored nodes in the branch-and-bound
tree and more simplex iterations.

Overall, if we consider only relatively hard problem instances that all
formulations need more than $10$ seconds to solve, then RS3 is $1.51$ \-
$1.90$ times faster than all Rosing et al. and Church formulations in Case 1
problems, whereas RRR2 is $1.07$ \- $1.37$ times faster than all ReVelle &
Swain and Church formulations in Case 2 ones. Note that some of the generated
instances are quite hard. The most difficult instance of Case 1 (resp. Case 2)
required $1,548$ (resp. $1,666$) seconds.

Regarding the distance constraints, there is no clear winner, indicating that
the more complex formulations do not always pay off. Most best execution times
are achieved by the distance constraints in 22, 24 and 25\- 26. The distance
constraint in 24 leads to better average execution times with the Rosing et
al. and Church formulations both in Case 1 problems and with the ReVelle &
Swain and Church formulations in Case 2 problems. Therefore, distance
constraint in 24 can be regarded as the best candidate for most cases on
average.

##  5 Conclusion

We have investigated for the first time a variant of the p-median problem
which concerns the location of semi-obnoxious facilities in the presence of
distance constraints between facilities and also between facilities and demand
points. We have presented twelve ILP models, coupling three formulations of
the p-median problem with four formulations of the distance constraints.
Experiments have demonstrated that, perhaps surprisingly, the classic p-median
formulation coupled with a simple model for the distance constraints is quite
competitive as it is the best one in problems with more candidate facility
sites than clients, and it is not far off the the best performer (i.e. the
Rosing et al.) in problems with more clients than candidate facility sites.

As future work we would like to consider specialized branch-and-bound
algorithms for this problem and extend our study to the case where the
facilities are not homogeneous. We can also extend/modify the proposed
techniques to related types of location problems, and of course we can
investigate the applicability of the methods to real-world problems.

## References

  * [1] Balinski, M.L.: Integer programming: methods, uses, computations. Management Science 12(3), 253–313 (1965)
  * [2] Beasley, J.E.: A note on solving large p-median problems. European Journal of Operational Research 21(2), 270–273 (1985)
  * [3] Berman, O., Huang, R.: The minimum weighted covering location problem with distance constraints. Computers and Operations Research 35(2), 356–372 (2008)
  * [4] Bron, C., Kerbosch, J.: Algorithm 457: finding all cliques of an undirected graph. Communications of the ACM 16(9), 575–577 (1973)
  * [5] Chaudhry, S., McCormick, T., Moon, I.D.: Locating independent facilities with maximum weight: Greedy heuristics. International Journal of Management Science 14(5), 383–389 (1986)
  * [6] Church, R.L.: Cobra: a new formulation of the classic p-median location problem. Annals of Operations Research 122(1), 103–120 (2003)
  * [7] Church, R.L., Meadows, M.E.: Results of a new approach to solving the p-median problem with maximum distance constraints. Geographical Analysis 9(4), 364–378 (1977)
  * [8] Chvatal, V.: A greedy heuristic for the set-covering problem. Mathematics of Operations Research 4(3), 233–235 (1979)
  * [9] Cornuejols, G., Nemhauser, G.L., Wolsey, L.A.: A canonical representation of simple plant location problems and its applications. SIAM Journal on Algebraic Discrete Methods 1(3), 261–272 (1980)
  * [10] Densham, P.J., Rushton, G.: Strategies for solving large location-allocation problems by heuristic methods. Environment and Planning A 24(2), 289–304 (1992)
  * [11] Drezner, T., Drezner, Z., Schöbel, A.: The weber obnoxious facility location model: a big arc small arc approach. Computers & Operations Research 98, 240–250 (2018)
  * [12] Drezner, Z., Kalczynski, P., Salhi, S.: The planar multiple obnoxious facilities location problem: A voronoi based heuristic. Omega 87, 105–116 (2019)
  * [13] Efroymson, M., Ray, T.: A branch-bound algorithm for plant location. Operations Research 14(3), 361–368 (1966)
  * [14] Hakimi, S.L.: Optimum locations of switching centers and the absolute centers and medians of a graph. Operations Research 12(3), 450–459 (1964)
  * [15] Hakimi, S.L.: Optimum distribution of switching centers in a communication network and some related graph theoretic problems. Operations Research 13(3), 462–475 (1965)
  * [16] Krarup, J., Pisinger, D., Plastria, F.: Discrete location problems with push–pull objectives. Discrete Applied Mathematics 123(1-3), 363–378 (2002)
  * [17] Moon, I.D., Chaudhry, S.: An analysis of network location problems with distance constraints. Management Science 30(3), 290–307 (1984)
  * [18] Moon, I.D., Papayanopoulos, L.: Minimax location of two facilities with minimum separation: Interactive graphical solutions. Journal of the Operations Research Society 42, 685–694 (1991)
  * [19] ReVelle, C.S., Swain, R.W.: Central facilities location. Geographical Analysis 2(1), 30–42 (1970)
  * [20] Rosing, K.E., ReVelle, C., Rosing-Vogelaar, H.: The p-median and its linear programming relaxation: An approach to large problems. Journal of the Operational Research Society 30(9), 815–823 (1979)
  * [21] Tansel, B.C., Francis, R.L., Lowe, T.J., Chen, M.: Duality and distance constraints for the nonlinear p-center problem and covering problem on a tree network. Operations Research 30(4), 725–744 (1982)