# Stability and optimality of multi-scale transportation networks with
distributed dynamic tolls

Rosario Maggistro1 and Giacomo Como1,2 *This reasearch was carried on within
the framework of the MIUR-funded Progetto di Eccellenza of the Dipartimento di
Scienze Matematiche G.L. Lagrange, CUP: E11G18000350001, and was partly
supported by the Compagnia di San Paolo and the Swedish Research
Council.1Department of Mathematical Sciences, Politecnico di Torino, Corso
Duca degli Abruzzi 24, 10129 Torino, Italy {rosario.maggistro,
<EMAIL_ADDRESS>Department of Automatic Control, Lund University,
BOX 118, SE-22100 Lund, Sweden<EMAIL_ADDRESS>
###### Abstract

We study transportation networks controlled by dynamical feedback tolls. We
consider a multiscale transportation network model whereby the dynamics of the
traffic flows are intertwined with those of the drivers’ route choices. The
latter are influenced by the congestion status on the whole network as well as
dynamic tolls set by the system operator. Our main result shows that a broad
class of decentralized congestion-dependent tolls globally stabilise the
transportation network around a Wardrop equilibrium. Moreover, using dynamic
marginal cost tolls, stability of the transportation network can be guaranteed
around the social optimum traffic assignment. This is particularly remarkable
as the considered decentralized feedback toll policies do not require any
global information about the network structure or the exogenous traffic load
on the network or state and can be computed in a fully local way. We also
evaluate the performance of these feedback toll policies both in the
asymptotic and during the transient regime, through numerical simulations.

###### Index Terms:

Transportation networks; traffic control; dynamic pricing; social-optimum
traffic assignment; Wardrop equilibrium; marginal cost pricing; dynamical flow
networks; robust distributed control.

##  I INTRODUCTION

In recent years. controlling the roadway congestion has become one of the main
target of the transportation research community. Proposed strategies include
imposing constraints on traffic flow through mechanisms such as variable speed
limits, ramp metering, or traffic signal control (see [1]–[4] and references
therein). However, such mechanisms do not consider neither the drivers’
perspective nor affect the total amount of vehicles. There has been also a
significant research effort to understand the drivers’ answer to external
communications from intelligent traveller information devices (see, e.g.,
[5]–[6]) and, in particular, studying the effect of such technologies on the
drivers’ route choice behaviour and on the dynamical properties of the
transportation network [7]. A traffic recommender which can announce
potentially misleading travel time information and a new class of latency
functions so as to influence the drivers’ behaviour was studied in [8] and
[9], respectively. Moreover, it is known that if individual drivers make their
own routing decisions to minimize their own experienced delays, overall
network congestion can be considerably higher than if a central planner had
the ability to explicitly direct traffic. Accordingly, to charge tolls for the
purpose of influencing drivers to make routing choices that result in globally
optimal routing was a central research focus (see [10]–[15]).

In this paper, we extend the model and results of [7] by introducing
decentralized congestion-dependent tolls in order to influence the driver’s
route choice behaviour. Specifically, we consider a multiscale dynamical model
of the transportation network whereby the traffic dynamics describing the real
time evolution of the local congestion level are coupled with those of the
drivers’ path preferences. We assume that the latter evolve following a
perturbed best response to global information about the congestion status of
the whole network and to decentralized flow-dependent tolls.

Our main result shows that by using non-decreasing decentralized flow-
dependent tolls and in the limit of a small update rate of the aggregate path
preferences, the transportation network globally stabilises around the Wardrop
equilibrium [16]. As in [7], we assume that the drivers’ path preferences
evolve at a slower time scale than the physical traffic flows and adopt a
singular perturbation approach [17] to the stability analysis of the ensuing
multiscale closed-loop traffic dynamics. Note that classic results of
evolutionary game theory and population dynamics [18]–[19] cannot be applied
to our framework since they suppose that the access to information take place
at a single temporal and spatial scale and that the traffic dynamics are
neglected by assuming that they are instantaneously equilibrated.

The introduction of tolls has long been studied as a way to influence the
rational and selfish behaviour of drivers so that the associated Wardrop
equilibrium can align with the system optimum network flow. A well-studied
taxation mechanism that guarantees this alignment is marginal-cost pricing
(see, e.g., [20] and [21]). Marginal-cost tolls do not require any global
information about the network structure, user demands or state and can be
computed in a fully local way. Using marginal-cost tolls we prove that our
transportation network stabilizes around the social optimum traffic
assignment. It is worth observing that our results go well beyond the
traditional setting [20] where only static frameworks are considered as well
as [21] where only path preference dynamics are consider, neglecting the
physical ones that are assumed equilibrated. In fact, our analysis is carried
over in a fully dynamical flow network setting. In this respect, the global
optimality guarantees that are obtained in this paper through decentralized
feedback toll policies should be compared with other recent results on global
performance and resilience results on robust distributed control of dynamical
flow networks [22]–[26].

In the last part of the paper through numerical simulations we compare the
performance both asymptotic and during the transient of the system by using
distributed marginal cost tolls and constant marginal cost ones. The latter,
know in the literature as “fixed” tolls (being the tolling function on each
edge a constant function of edge flow) have been well studied, and it is known
that they can be computed to enforce the social optimum equilibrium provided
that the system planner has a complete knowledge of the network topology, user
demand profile and delay functions. We show that not only is more convenient
take into account the marginal cost tolls at convergence speed level but also
they are strongly robust to variation of network topology, user demand and
traffic rate (see [27] and [28]).

The rest of this paper is organized as follows. In Section II, we describe the
model and observe the influence of distributed dynamics tolls on the network
dynamics. In Section III we state the main results of the paper. In Section IV
we provide a numerical study of the different time and asymptotic convergences
of the system. Section V draws conclusions and suggests future works. Due to
space limitations, we do not include any proofs of our results here and refer
the reader to a forthcoming journal publication [29].

###  I-A Notation

Let $\mathbb{R}$ and $\mathbb{R}_{+}:=\\{x\in\mathbb{R}:x\geq 0\\}$ be the set
of real and nonnegative real numbers, respectively. Let $\mathcal{A}$ and
$\mathcal{B}$ be finite sets. Then $|\mathcal{A}|$ denotes the cardinality of
$\mathcal{A}$, $\mathbb{R}^{\mathcal{A}}$ the space of real-valued vectors
whose components are indexed by elements of $\mathcal{A}$, and
$\mathbb{R}^{\mathcal{A}\times\mathcal{B}}$ the space of real-valued matrices
whose entries are indexed by pairs in $\mathcal{A}\times\mathcal{B}$. The
transpose of a matrix $Q\in\mathbb{R}^{\mathcal{A}\times\mathcal{B}}$ is
denoted by $Q^{\prime}\in\mathbb{R}^{\mathcal{B}\times\mathcal{A}}$, $I$ is an
identity matrix and $\mathbf{1}$ an all one vector whose size depends on the
context. We use the notation
$\Phi:=I-|\mathcal{A}|^{-1}\mathbf{11^{\prime}}\in\mathbb{R}^{\mathcal{A}\times\mathcal{A}}$
to denote the projection matrix of the space orthogonal to $\mathbf{1}$. The
simplex of a probability vector over $\mathcal{A}$ is denoted by
$S(\mathcal{A})=\\{x\in\mathbb{R}_{+}^{\mathcal{A}}:\mathbf{1}^{\prime}x=1\\}$.
Let $\|\cdot\|_{p}$ be the class of $p$-norms for $p\in[1,\infty]$, and by
default, let $\|\cdot\|:=\|\cdot\|_{2}$. Let now
$\text{sgn}:\mathbb{R}\to\\{-1,0,1\\}$ be the sign function, defined by
$\text{sgn}(x)=1$ if $x>0$, $\text{sgn}(x)=-1$ if $x<0$ and $\text{sgn}(x)=0$
if $x=0$. By convention, we will assume the identity $d|x|/dx=\text{sgn}(x)$
to be valid for every $x\in\mathbb{R}$, including $x=0$. Finally, given the
gradient $\nabla f$ of a function $f:D\to\mathbb{R}$ with
$D\subseteq\mathbb{R}^{\mathcal{A}}$, we denote with
$\tilde{\nabla}f=\Phi\nabla f$ the projected gradient on $S(\mathcal{A})$.

##  II MODEL DESCRIPTION

###  II-A Network characteristics

We describe the topology of the transportation network by a directed multi-
graph $\mathcal{G=(V,E)}$, where $\mathcal{V}$ is a finite set of nodes and
$\mathcal{E}$ is the set of links $e$, each directed from its tail node
$\theta_{e}$ to its head node $\kappa_{e}\neq\theta_{e}$. We shall allow for
parallel links, i.e. $\theta_{e}=\theta_{j}$ and $\kappa_{e}=\kappa_{j}$ with
$e\neq j$, but not for self loops, i.e., we shall assume that
$\theta_{e}\neq\kappa_{e}$ for every $e\in\mathcal{E}$. We shall denote by
$B\in\\{-1,0,1\\}^{\mathcal{V}\times\mathcal{E}}$ the node-link incidence
matrix of $\mathcal{G}$, whose entries are defined as $B_{ie}=1$ if
$i=\theta_{e}$, $B_{ie}=-1$ if $i=\kappa_{e}$, and $B_{ie}=0$ otherwise. For
two nodes $o\neq d$ in $\mathcal{V}$, an $o$-$d$ path is a length-$\ell$
string of links $p=(e_{1},e_{2},\ldots,e_{\ell})$ such that
$\theta_{e_{s+1}}=\kappa_{e_{s}}$ for $s=1,\ldots,\ell-1$, $\theta_{e_{1}}=o$,
$\kappa_{e_{\ell}}=d$, and no node is touched twice, i.e., $i_{r}\neq i_{s}$
for all $0\leq r<s\leq\ell$. The set of $o-d$ paths in $\mathcal{G}$ of any
length $\ell$ will be denoted by $\mathcal{P}$. Moreover, we shall denote the
corresponding link-path incidence matrix by
$A\in\\{0,1\\}^{\mathcal{E}\times\mathcal{P}}$ with entries

$A_{ep}=\begin{cases}1&\text{ if }e\text{ is along }p\\\
0&\text{otherwise}\end{cases}$

and assume that each link $e\in\mathcal{E}$ lies on at least one path from
node $o$ to node $d$. A path of length greater than or equal to $2$ from a
node to itself is referred to as a cycle. Observe that, in contrast to [7]
where the transportation network was assumed acyclic, we allow for the
presence of cycles in the network topology $\mathcal{G}$. For every link
$e\in\mathcal{E}$ and time instant $t\geq 0$ we denote the current traffic
density and flow by $x_{e}(t)$ and $f_{e}(t)$ respectively, and assume the
following functional dependence

$f_{e}=\mu_{e}(x_{e}),\quad e\in\mathcal{E},$ (1)

such that $\mu_{e}:\mathbb{R}_{+}\to\mathbb{R}_{+}$ is continuously
differentiable, strictly increasing, strictly concave and
$\mu_{e}(0)=0,\quad\mu_{e}^{\prime}(0)<\infty.$ Note that in classical
transportation theory the flow-density function are typically not strictly
increasing, but here our assumption is valid as long as we confine ourselves
to the free-flow region, as is done in [7]. Then, for every link
$e\in\mathcal{E}$, let $C_{e}:=\sup\\{\mu_{e}(x_{e}):x_{e}\geq 0\\}$ be its
maximum flow capacity and let $\mathcal{F}:=\prod_{e\in\mathcal{E}}[0,C_{e})$
be the set of feasible flow vectors. We shall use the delay functions

$T:\mathbb{R}_{+}^{\mathcal{E}}\to[0,+\infty]^{\mathcal{E}}\,,$
$T_{e}(f_{e}):=\begin{cases}\displaystyle+\infty&\text{if}\ f_{e}\geq
C_{e},\\\\[10.0pt] \displaystyle\frac{\mu_{e}^{-1}(f_{e})}{f_{e}}&\text{if}\
f_{e}\in(0,C_{e}),\\\\[10.0pt]
\displaystyle\frac{1}{\mu^{\prime}_{e}(0)}\quad&\text{if}\ f_{e}=0\end{cases}$
(2)

returning the delay incurred by drivers traversing link $e\in\mathcal{E}$,
when the current flow out of it is $f_{e}$. Note that, by the properties of
$\mu_{e}$, $T_{e}(f_{e})$ is continuous, strictly increasing, and such that
$T_{e}(0)>0$.

###  II-B Paths choice and traffic dynamics

We assume that the physical traffic flow consist of indistinguishable
homogeneous drivers which enter in the network through the origin node, travel
through it using the different paths and finally exit from the network through
the destination node. The relative appeal of the different paths to the
drivers is modelled by a time-varying probability vector over $\mathcal{P}$,
which will be referred as the current aggregate path preference and denoted by
$z(t)$. Assuming a constant unit in-flow in the origin node, we consider the
vector

$f^{z}:=Az$

of the flows associated to the path preference $z(t)$ and define

$Z:=\\{z\in\mathcal{S(\mathcal{P})}:f_{e}^{z}<C_{e}\ \forall
e\in\mathcal{E}\\}$

the set of feasible path preference. The vector $z(t)$ is updated as drivers
access global information about the current congestion status of the whole
network (that is embodied by the flow vector $f(t)$) and is influenced by a
vector of decentralized congestion-dependent tolls

$w:\mathbb{R}_{+}^{\mathcal{E}}\to[0,+\infty]^{\mathcal{E}},\quad
w_{e}(f_{e})\geq 0\quad\forall e\in\mathcal{E},$ (3)

that are charged to users traversing link $e$. In particular, we shall assume
that the tolls $w_{e}$ are continuous and non-decreasing functions of the
current flow for every link $e\in\mathcal{E}$.

We shall assume that the cost perceived by each user crossing a link
$e\in\mathcal{E}$ is given by the sum of the the delay $T_{e}(f_{e})$ and the
toll $w_{e}(f_{e})$. Moreover, as in [7], we shall assume that path
preferences are updated at some rate $\eta>0$ which is small with respect to
the time scale of the network flow dynamics. Then, from $f(t)$, the drivers
evaluate the vector $A^{\prime}(T(f(t))+w(f(t)))$, whose $p$th entry,
$\sum_{e}A_{ep}(T_{e}(f_{e}(t))+w_{e}(f_{e}(t)))$, coincides with the
perceived total cost that a driver expects to incur on path $p$ assuming that
the congestion levels on that path won’t change during the journey. Hence,
according to some feasible path preference $F^{h}(f(t))\in Z$, $z(t)$ evolves
as

$\dot{z}(t)=\eta(F^{h}(f(t))-z(t)),$ (4)

where $F^{h}:\mathcal{F}\to Z$ is a perturbed best response function,

$F^{h}(f):=\operatornamewithlimits{\arg\min}_{\alpha\in
Z_{h}}\\{\alpha^{\prime}A^{\prime}(T(f)+w(f))+h(\alpha)\\},\quad
f\in\mathcal{F},$ (5)

and $h:Z_{h}\to\mathbb{R}$ is an admissible perturbation such that
$Z_{h}\subseteq Z$ is a closed convex set, $h(\cdot)$ is strictly convex,
twice differentiable in $int(Z_{h})$, and is such that $\lim_{z\to\partial
Z_{h}}\lVert\tilde{\nabla}h(z)\rVert=\infty$. The definition of $F^{h}$ and
the conditions on $h$ imply that $F^{h}(f)\in int(Z_{h})$ and that $F^{h}(f)$
is differentiable on $\mathcal{F}$.
We now describe the local route decisions, characterizing the fraction of
drivers choosing each outgoing link when traversing a nondestination node.
Such a fraction is the function $G_{e}(z)$ defined as

$G_{e}(z)=\begin{cases}\displaystyle\frac{f_{e}^{z}}{\displaystyle\sum_{j\in\mathcal{E}:\theta_{j}=\theta_{e}}f_{j}^{z}}&\
\text{if}\ f_{e}^{z}>0,\\\
\displaystyle\frac{1}{|\\{j\in\mathcal{E}:\theta_{j}=\theta_{e}\\}|}&\
\text{if}\
\displaystyle\sum_{j\in\mathcal{E}:\theta_{j}=\theta_{e}}f_{j}^{z}=0,\end{cases}$
(6)

for every $e\in\mathcal{E}$. Note that $\sum_{k}G_{k}(z)=1$, where $k$ are the
outgoing links from the same node.
We refer to $G:Z\to\mathbb{R}^{\mathcal{E}}$ as the local decision function
that is continuously differentiable on $Z$.
Now, for every $e\in\mathcal{E}$ conservation of mass implies that

$\dot{x}_{e}(t)=H_{e}(f(t),z(t)),$ (7)

where for all $z\in Z$ and $f\in\mathcal{F}$,

$H_{e}(f,z):=G_{e}(z)\bigg{(}\delta_{\theta_{e}}^{(o)}+\sum_{j:\kappa_{j}=\theta_{e}}f_{j}\bigg{)}-f_{e}.$
(8)

We now consider the evolution of the coupled dynamics

$\begin{cases}\dot{z}(t)=\eta(F^{h}(f(t))-z(t)),\\\
\dot{x}(t)=H(f(t),z(t))\end{cases}$ (9)

where $F^{h}$ is defined in (5), $\eta>0$ is the rate at which $z(t)$ is
updated and $H(f,z)=\\{H_{e}(f,z):e\in\mathcal{E}\\}$.

##  III Main results

In this section we give the main results of the paper. We shall prove that for
small $\eta$ and $h$, the long-time behaviour of the system (9) is
approximately at Wardrop equilibrium [16] which, under proper distributed
dynamic tolls, coincides with the social optimum equilibrium.

###### Definition 1

(Social optimum equilibrium). A feasible flow vector $f^{*}\in\mathcal{F}$ is
a Social optimum equilibrium if and only if is the unique solution of the
following network flow optimization problem

$f^{*}=\operatornamewithlimits{\arg\min}_{\begin{subarray}{c}f\geq 0\\\
Bf=(\delta^{(o)}-\delta^{(d)})\end{subarray}}\sum_{e\in\mathcal{E}}f_{e}T_{e}(f_{e}).$
(10)

###### Definition 2

(Wardrop equilibrium). For a given vector $w\in\mathbb{R}_{+}^{\mathcal{E}}$
of decentralized link tolls, a feasible flow vector $f^{(w)}\in\mathcal{F}$ is
a Wardrop equilibrium if $f^{(w)}=f^{z}$ for some $z\in Z$ such that for all
$p\in\mathcal{P}$,

$z_{p}>0\quad\Longrightarrow\\\
\begin{array}[]{ll}\left(A^{\prime}\left(T(f^{z})+w(f^{z})\right)\right)_{p}\leq\\\
\left(A^{\prime}\left(T(f^{z})+w(f^{z})\right)\right)_{q}\quad\forall
q\in\mathcal{P}.\end{array}$ (11)

Existence and uniqueness of a Wardrop equilibrium are guaranteed considering
the direct multi-graph $\mathcal{G}$ and under the assumption on $\mu_{e}$ and
$w_{e}$. (See Theorem 2.4 and 2.5 in [30] for a complete proof).

###### Theorem 3

Let $\mathcal{G}$ be the direct multi-graph, $\mu$ be as in (1) and $w$ as in
(3). Then for every initial condition $(z(0),x(0))\in
Z\times[0,+\infty)^{\mathcal{E}}$ there exists a unique solution of (9).
Moreover, there exists a perturbed equilibrium flow $f^{(h)}\in\mathcal{F}$
such that for all $\eta>0$

$\limsup_{t\to\infty}\lVert f(t)-f^{(h)}\rVert\leq\delta(\eta),$ (12)

where $\delta(\eta)$ is a non negative real-valued, nondecreasing function
such that $\lim_{\eta\to 0}\delta(\eta)=0$. Moreover, for every sequence of
admissible perturbations $\\{h_{k}\\}$ such that $\lim_{k}\lVert
h_{k}\rVert=0$ and $\lim_{k}Z_{h_{k}}=\overline{Z}$ 111The convergence
$\lim_{k}Z_{h_{k}}=\overline{Z}$ holds with respect to the Hausdorff metric
and $\overline{Z}$ is the closure of $Z$., one has

$\lim_{k\to\infty}f^{(h_{k})}=f^{(w)}.$ (13)

Theorem 3 states that the system planner globally stabilises the
transportation network around the Wardrop equilibrium using increasing
decentralised congestion-dependent tolls.

###### Remark 4

Note that Theorem 3 is not a Corollary of Theorem 2.5 in [7], because,
although the functions $T$ and $w$ both depend on the flow $f$, it is not
possible consider an auxiliary function $\overline{T}=T+w$ and directly
applying the result from [7] due to the specific structure imposed on T in
(2).

Now, we choose as decentralized tolls the marginal cost ones, namely,

$w_{e}(f_{e})=f_{e}T^{\prime}_{e}(f_{e})\quad\forall e\in\mathcal{E}.$ (14)

Due the properties of the delay function $T_{e}(f_{e})$, the above tolls (14)
are increasing, then the Theorem 3 continue to hold. Moreover the following
holds

###### Corollary 5

Considering (14) one gets that the system (9) globally stabilises the
transportation network around the social optimum traffic assignment $f^{*}$
without knowing arrival rates or the network structure.

In order to prove the above we observe that considering proper costs on the
links, the vector $f^{(w)}$ is the solution of a network flow optimization
problem. Let

$D_{e}(f_{e}):=\int_{0}^{f_{e}}\Big{(}T_{e}(s)+sT^{\prime}_{e}(s)\Big{)}\,ds\quad
e\in\mathcal{E},$

be the integral of the perceived cost on link $e$ using (14). Then, the
network flow $f^{(w)}\in\mathbb{R}_{+}^{\mathcal{E}}$ is a Wardrop equilibrium
if and only if is the unique solution of the network flow optimization problem

$f^{(w)}=\operatornamewithlimits{\arg\min}_{\begin{subarray}{c}f\geq 0\\\
Bf=(\delta^{(o)}-\delta^{(d)})\end{subarray}}\sum_{e\in\mathcal{E}}D_{e}(f_{e}),$
(15)

where $Bf=(\delta^{(o)}-\delta^{(d)})$ is the mass conservation law. Moreover,
the Wardrop equilibrium coincides with the system optimum flow,

$f^{(w)}=f^{*}.$ (16)

The proof of such result is very simple and use the Lagrange techniques.

###### Remark 6

The tolls (14) differ by the well now decentralized constant marginal cost
tolls $w_{e}^{*}=f_{e}^{*}T^{\prime}_{e}(f_{e}^{*})\quad\forall
e\in\mathcal{E}$, since the latter, in order to be used, require the knowledge
of both of the social optimum flow and the inflow vector. Anyway taking into
account such $w_{e}^{*}$, condition (16) continue to hold.

##  IV Asymptotic and transient performances

In this section, through numerical simulations we will compare the different
performances both asymptotic and during the transient given by using the
marginal cost tolls (14) and the constant marginal cost ones (see the Remark
6). We performed several experiments with different graph topologies for
$\eta$ ranging from 0.1 to 50. In all these cases we found that the use of the
decentralized marginal cost tolls is more convenient than the constant
marginal ones. Indeed:
$-$ concerning the transient convergence, one shows that the time needed to
reach the perturbed equilibrium associated to the marginal cost tolls is lower
than the one to reach the equilibrium associated to the constant marginal
ones;
$-$ when the admissible perturbation goes to zero, the perturbed equilibrium
associated to marginal cost tolls, asymptotically converges to the social
optimum flow faster than the one associated to the constant marginal cost
ones.
We demonstrate these findings through the following example. The parameters
were selected as follows:

  * •

graph topology $\mathcal{G}$ as in Fig. 1;

  * •

the flow-density function is

$\mu_{e}(x_{e})=2(1-e^{-x_{e}})\quad\forall e\in\mathcal{E},$

and the corresponding delay function, according to (2) is given by

$T_{e}(f_{e})=\begin{cases}+\infty&\text{if}\ f_{e}\geq 2,\\\
\displaystyle\frac{1}{f_{e}}\log\left(\frac{2}{2-f_{e}}\right)&\text{if}\
f_{e}\in(0,2),\\\ 1/2&\text{if}\ f_{e}=0.\end{cases}$ (17)

  * •

$F^{h}$ as the logit function

$F_{p}^{h}(f)=\frac{\exp(-\beta(A^{\prime}(T(f)+w(f)))_{p})}{\sum_{q\in\mathcal{P}}\exp(-\beta(A^{\prime}(T(f)+w(f)))_{q})},\
p\in\mathcal{P},$ (18)

with $\beta>0$ the fixed noise parameter.

  * •

$\eta=0.1$, $G$ as in (6),

  * •

initial conditions: $z_{p_{1}}(0)=1/2$, $z_{p_{2}}(0)=1/6$,
$z_{p_{3}}(0)=1/3$, $x_{e_{1}}(0)=4$, $x_{e_{2}}(0)=2$, $x_{e_{3}}(0)=3$,
$x_{e_{4}}(0)=1$, $x_{e_{5}}(0)=5$.

By the implementations follows that for $t\in[0,350]$ and $\beta=1$, the first
time in which the system reaches the equilibrium associated to (14) is
$t=2.17\cdot 10^{2}$, while it is $t=2.5\cdot 10^{2}$ the one to approach the
equilibrium relative to $w_{e}^{*}$.
The 1-norm distance of $f^{\beta}$ (that is the perturbed equilibrium flow
corresponding to the system (9) using (18)) computed at final time $T=350$,
from the social optimum flow $f^{*}$ for $\beta$ ranging from 1 to 12 is
plotted in Fig. 2. This is done both considering (14) and $w_{e}^{*}$.

oabd$e_{1}$$e_{2}$$e_{3}$$e_{4}$$e_{5}$$1$$1$ Figure 1:  The graph topology
used for the simulations.

Note that the parameter $\beta$ should takes very large values in order to
completely vanish the norm of the difference between $f^{\beta}$ and $f^{*}$;
but, in our numerical example, we can see in Fig. 2 that already for
$\beta=12$ the previous norm is almost null and also the asymptotic
convergence of $f^{\beta}$ associated to (14) is slightly faster than the one
of $f^{\beta}$ associated to $w_{e}^{*}$.

Figure 2:  Plot of $\|f^{\beta}(T)-f^{*}\|_{1}$ for decentralised marginal and
constant marginal tolls .

###  IV-A Robustness

To investigate the robustness of the marginal cost tolls to variations of
network’s parameters, a system planner can study the effect of the variation
on the total latency computed in $f^{(w)}$, where the total latency is defined
as

$\mathcal{L}(f)=\sum_{e\in\mathcal{E}}f_{e}T_{e}(f_{e}).$

By corollary 5 follows that the efficiency guarantees provided by the marginal
cost tolls are robust to variation in network and demand structure. Indeed the
following hold:

###### Proposition 7

(See [20]) For homogeneous populations, the marginal cost tolls (14)
incentives optimal flows on all networks, i.e.,

$\mathcal{L}(f^{(w)})=\mathcal{L}(f^{*}).$ (19)

Hence, the marginal cost tolls are strongly robust to variations of network
topology, user demand structure and overall traffic rate. In the following we
will show (see Fig. 3), still using the graph topology in Fig. 1 and its
parameters, that

$\lim_{\beta\to+\infty}\mathcal{L}(f^{\beta})=\mathcal{L}(f^{*})$

and the asymptotic convergence using $f^{\beta}$ associated to (14) is lightly
faster than the one in which using $f^{\beta}$ associated to $w_{e}^{*}$.

Figure 3:  Plot of the difference
$\mathcal{L}(f^{\beta}(T))-\mathcal{L}(f^{*})$ as $\beta$ increases.

##  V CONCLUSIONS

In this paper, we studied stability of Wardrop equilibria of multi-scale
transportation networks with distributed dynamic tolls. We prove that if the
frequency of updates of path preferences is sufficiently small and considering
positive, non-decreasing decentralized flow-dependent tolls, then the state of
the network ultimately approaches a neighborhood of the Wardrop equilibrium.
Then, using a particular class of tolls, i.e., the marginal cost ones, we
observe that the stability is around the social optimum equilibrium and,
thanks to numerical experiments, the performances both asymptotic and during
the transient of the system is better than the one obtained considering the
constant marginal tolls. In future research, inspired by the numerical results
we will provide analytic estimates about the different convergence rates. We
also plan to define a more general class of tolls that does not require the
knowledge of the delay functions and at the same time guarantees the
convergence to the social optimum.

## References

  * [1] A. Hegyi, B. De Schutter, and H. Hellendoorn, “Model predictive control for optimal coordination of ramp metering and variable speed limits,” Transport Res C: Emer, vol. 13, no. 3, pp. 185-209, 2005.
  * [2] G. Gomes and R. Horowitz, “Optimal freeway ramp metering using the asymmetric cell transmission model,” Transport Res C: Emer, vol. 14, no. 4, pp. 244-262, 2006.
  * [3] G. Como, E. Lovisari, and K. Savla, “Convexity and robustness of dynamic traffic assignment and freeway network control,” Transp. Res. B: Methodol., vol. 91, pp. 446–465, 2016.
  * [4] P. Varaiya, “Max pressure control of a network of signalized intersections,” _Transport Res C: Emer._ , vol. 36, pp. 177–195, 2013.
  * [5] K. Srinivasan and H. Mahmassani,“ Modeling inertia and compliance mechanisms in route choice behavior under real-time information,” Transport. Res. Rec., no. 1725, pp. 45–53, 2000.
  * [6] A. Khattak, A. Polydoropoulou, and M. Ben-Akiva, “Modeling revealed and stated pretrip travel response to advanced traveler information systems,” Transport. Res. Rec., no. 1537, pp. 46–54, 1996.
  * [7] G. Como, K. Savla, D. Acemoglu, M.A. Dahleh, and E. Frazzoli, “Stability analysis of transportation networks with multiscale driver decisions,” SIAM J. Control Optim., vol. 51, no. 1, pp. 230–252, 2013.
  * [8] Y. Cheng and C. Langbort, “A model of informational nudging in transportation networks,” 55th IEEE Conference on Decision and Control, pp. 7598–7604, 2016.
  * [9] W. Krichene, J.D. Reilly, S. Amin, and A.M. Bayen, “Stackelberg Routing on Parallel Transportation Networks. In: Basar T., Zaccour G. (eds) Handbook of Dynamic Game Theory. Springer, Cham, 2017.
  * [10] M.J. Smith, “The marginal cost taxation of a transportation network,” Transp. Res. B: Methodol., vol. 13, no. 3, pp. 237–242, 1979.
  * [11] S. Morrison, “A survey of road pricing,” Transp. Res. A: Gen., vol. 20, no. 2, pp. 87–97, 1986.
  * [12] R.B. Dial, “Network-optimized road pricing: Part I: A parable and a model,” Oper. Res., vol. 47, pp. 54–64, 1999.
  * [13] R. Cole, Y. Dodis, and T. Roughgarden, “How much can taxes help selfish routing?,” J. Comput. Syst. Sci., vol. 72, pp. 444–467, 2006.
  * [14] L. Engelson and P. Lindberg, “ Congestion pricing of road networks with users having different time values,” Appl. Optim., vol. 101, pp. 81-104, 2006.
  * [15] G. Christodoulou, K. Mehlhorn, and E. Pyrga, “Improving the price of anarchy for selfish routing via coordination mechanisms ,” Algorithmica, vol. 69, no. 3, pp. 619–640, 2014.
  * [16] J.G. Wardrop, “Some theoretical aspects of road traffic research,” ICE Proc. Engrg. Divisions, vol. 1, no. 3, pp. 325–362, 1952.
  * [17] H.K. Khalil, Nonlinear Systems, 2nd ed., Prentice-Hall, Englewood Cliffs, NJ, 1996.
  * [18] J. Hofbauer and K. Sigmund, “Evolutionary game dynamics,” Bull. Amer. Math. Soc., vol. 40, pp. 479–519, 2003.
  * [19] W.H. Sandholm, Population Games and Evolutionary Dynamics, MIT Press, Cambridge, MA, 2011.
  * [20] M. Beckmann, C. McGuire, and C.B. Winsten, Studies in the Economics of Transportation, New Haven, CT: Yale University Press, 1956.
  * [21] W. Sandholm, “Evolutionary implementation and congestion pricing,” Rev. Econ. Stud., vol. 69, no. 3, 667-689, 2002.
  * [22] G. Como, K. Savla, D. Acemoglu, M.A. Dahleh, and E. Frazzoli. “Robust distributed routing in dynamical networks-Part I: Locally responsive policies and weak resilience,” IEEE Trans. Automat. Control, vol. 58, no. 2, pp. 317-332, 2013.
  * [23] G. Como, K. Savla, D. Acemoglu, M.A. Dahleh, , and E. Frazzoli, “Robust distributed routing in dynamical networks-Part II: Strong resilience, equilibrium selection and cascaded failures,” IEEE Trans. Automat. Control, vol. 58, no. 2, 333-348, 2013.
  * [24] G. Como, E. Lovisari, and K. Savla, “ Throughput optimality and overload behaviour of dynamical flow networks under monotone distributed routing,” IEEE Trans. Control Netw. Syst., vol. 2, no. 1, 57-67, 2015.
  * [25] G. Como, “On resilient control of dynamical flow networks”, _Annual Reviews in Control_ , vol. 43, pp. 70–80, 2017.
  * [26] A. Y. Yazicioglu, M. Roozbehani, and M. A. Dahleh, “Resilient Control of Transportation Networks by Using Variable Speed Limits,” IEEE Trans. Control Netw. Syst., DOI 10.1109/TCNS.2017.2782364, 2017.
  * [27] P.N. Brown and J.R. Marden, “Studies on robust social influence mechanisms: Incentives for efficient network routing in uncertain settings,” IEEE Control Systems, vol. 37, no. 1, pp. 98-115, 2017.
  * [28] P.N. Brown and J.R. Marden,“The robustness of marginal-cost taxes in affine congestion games ,” IEEE Trans. Autom. Control., vol. 62, no. 8, pp. 3999-4004, 2017.
  * [29] G. Como and R. Maggistro, “On Robust Distributed Dynamic Pricing in Multiscale Transportation Networks,” 2018.
  * [30] M. Patriksson, The Traffic Assignment Problem: Models and Methods, VSP International Science, Leiden, Netherlands, 1994.


# Optimal approximation of spherical squares by tensor product quadratic
Bézier patches

Aleš Vavpetič<EMAIL_ADDRESS>Emil Žagar<EMAIL_ADDRESS>lj.si Faculty of Mathematics and Physics, University of Ljubljana, Jadranska
19, Ljubljana, Slovenia  Institute of Mathematics, Physics and Mechanics,
Jadranska 19, Ljubljana, Slovenia

###### Abstract

In [1], the author considered the problem of the optimal approximation of
symmetric surfaces by biquadratic Bézier patches. Unfortunately, the results
therein are incorrect, which is shown in this paper by considering the optimal
approximation of spherical squares. A detailed analysis and a numerical
algorithm are given, providing the best approximant according to the
(simplified) radial error, which differs from the one obtained in [1]. The
sphere is then approximated by the continuous spline of two and six tensor
product quadratic Bézier patches. It is further shown that the $G^{1}$ smooth
spline of six patches approximating the sphere exists, but it is not a good
approximation. The problem of an approximation of spherical rectangles is also
addressed and numerical examples indicate that several optimal approximants
might exist in some cases, making the problem extremely difficult to handle.
Finally, numerical examples are provided that confirm theoretical results.

###### keywords:

Bézier patch , spherical square , optimal approximation , sphere approximation

###### MSC:

[2020] 65D05 , 65D07 , 65D15 , 65D17

††journal: Computer Aided Geometric Design

##  1 Introduction

One of the most important issues of computer-aided geometric design (CAGD) is
the approximation of curves and surfaces by simple polynomial-based objects.
It is well known that even fundamental geometric objects such as circular arcs
or spherical patches do not possess exact polynomial representations. It is
thus an interesting and practical issue to find their best polynomial
approximants. Much work has been done on the optimal polynomial approximation
of circular arcs (see, e.g., [4] and the references therein). However, much
less is known about the optimal approximation of spherical patches. The main
two classes of the polynomial (or spline) surfaces used in CAGD are triangular
patches and tensor product patches. The optimal approximation of equilateral
spherical triangles by triangular Bézier patches was studied recently in [5].
Approximation of rational tensor-product biquadratic Bézier surfaces (which
also includes the sphere) by polynomial tensor product patches was proposed in
[2]. An exact sphere representation by rational S-patches was considered
recently in [3]. In [1], the author studied the problem of the best
approximation of symmetric surfaces by biquadratic Bézier surfaces. However,
we show in this paper that the obtained results are incorrect since the
constructed approximants are not optimal even in the spherical case. Our
primary purpose is to characterize the optimal polynomial approximant and
provide an algorithm for its construction.

The paper is organized as follows. After the introduction in Section 1 some
basic preliminaries are explained in Section 2. The detailed analysis of the
construction of the best approximant together with a numerical algorithm are
given in Section 3. The next section considers $G^{1}$ smooth optimal
approximants of the sphere. In Section 5 some basic observations of the
approximation of spherical rectangles are provided. Numerical examples are
shown in Section 6, and some closing remarks are given in the last section.

##  2 Preliminaries

Let ${\mathcal{S}}$ be the unit sphere in $\mathbb{R}^{3}$ centered at
$\boldsymbol{0}=(0,0,0)$ and let ${\mathcal{R}}$ be the spherical square for
which the projection of its vertices along the vector $[0,0,1]^{T}$ are
vertices of the square centered at $(0,0)$ with its side equal to $2a$, where
$0<a\leq\tfrac{\sqrt{2}}{2}$. Let ${\mathcal{P}}$ be a tensor product
quadratic Bézier patch parameterized as

$\boldsymbol{p}(u,v)=\sum_{i=0}^{2}\sum_{j=0}^{2}B_{i}^{2}(u)B_{j}^{2}(v)\boldsymbol{b}_{ij},\quad
u,v\in[-1,1],$ (1)

where

$B_{i}^{2}(u)=\sum_{i=0}^{2}\binom{2}{i}\left(\frac{1+u}{2}\right)^{i}\left(\frac{1-u}{2}\right)^{2-i},\quad
i=0,1,2,$

are quadratic Bernstein polynomials parametrized on $[-1,1]$ and
$\boldsymbol{b}_{ij}\in\mathbb{R}^{3}$, $i,j=0,1,2$, are the corresponding
control points. Our goal is to find the optimal approximation of
${\mathcal{R}}$ by ${\mathcal{P}}$ according to the radial error

$\max_{(u,v)\in[-1,1]^{2}}\left|\|\boldsymbol{p}(u,v)\|_{2}-1\right|$

or according to the simplified radial error

$\max_{(u,v)\in[-1,1]^{2}}\left|f(u,v)\right|,\quad
f(u,v)=\|\boldsymbol{p}(u,v)\|_{2}^{2}-1.$ (2)

Note that the latter one is easier to handle since $f$ is a polynomial
function of the coefficients of $\boldsymbol{p}$. Due to this reason, we will
consider this error in the following. However, once the problem is solved for
the simplified radial error, an identical approach with $f$ replaced by
$g=\sqrt{f+1}-1$ can be used to find the best approximant according to the
radial error.

In order to get an appropriate approximation of ${\mathcal{R}}$ it is natural
to require that ${\mathcal{P}}$ interpolates ${\mathcal{R}}$ at least at its
four vertices and that the control points of the boundary curves of
${\mathcal{P}}$ lie in the planes passing through corresponding two endpoints
and the origin $\boldsymbol{0}$. Consequently, by symmetry we have

$\displaystyle\boldsymbol{b}_{00}$
$\displaystyle=(-a,-a,\sqrt{1-2a^{2}}),\quad\boldsymbol{b}_{20}=(a,-a,\sqrt{1-2a^{2}}),\quad\boldsymbol{b}_{02}=(-a,a,\sqrt{1-2a^{2}}),\quad\boldsymbol{b}_{22}=(a,a,\sqrt{1-2a^{2}}),$
$\displaystyle\boldsymbol{b}_{10}$
$\displaystyle=\alpha(\boldsymbol{b}_{00}+\boldsymbol{b}_{20}),\quad\boldsymbol{b}_{01}=\alpha(\boldsymbol{b}_{00}+\boldsymbol{b}_{02}),\quad\boldsymbol{b}_{21}=\alpha(\boldsymbol{b}_{20}+\boldsymbol{b}_{22}),\quad\boldsymbol{b}_{12}=\alpha(\boldsymbol{b}_{02}+\boldsymbol{b}_{22}),$
(3)

and $\boldsymbol{b}_{11}=\beta(0,0,1)$, for some unknown real parameters
$\alpha,\beta$. It is quite clear that $\alpha,\beta>0$. Thus
$\boldsymbol{p}=\boldsymbol{p}(\cdot,\cdot,\alpha,\beta)$ and consequently
$f=f(\cdot,\cdot,\alpha,\beta)$ which implies that we are looking for a
minimum of (2), i.e.,

$\min_{\alpha,\beta>0}\max_{(u,v)\in[-1,1]^{2}}\left|f(u,v,\alpha,\beta)\right|.$

Let $\angle=\\{(u,v)\in[-1,1]^{2};u=v\text{ or }v=-1\\}$ and define

$f_{s}(u,\alpha)=f(u,-1,\alpha,\beta),\quad
f_{d}(u,\alpha,\beta)=f(u,u,\alpha,\beta),\quad{\rm and\
}f_{s,d}=f|_{\angle}.$

We shall first find parameters $\alpha^{*}$ and $\beta^{*}$ which minimize the
maximum of $|f|$ over $\angle$ and finally show that they provide global
minimization of $|f|$ over $[-1,1]^{2}$, i.e.,

$(\alpha^{*},\beta^{*})={\rm{argmin}}_{\alpha,\beta>0}\max_{(u,v)\in[-1,1]^{2}}|f(u,v,\alpha,\beta)|.$

Note that $f_{s}(\cdot,\alpha)$ and $f_{d}(\cdot,\alpha,\beta)$ are even
functions on $[-1,1]$ and it is thus enough to analyse them only on $[0,1]$ if
it is more convenient.

The author in [1] claims that for the best approximant $\boldsymbol{p}$ the
error $f_{s}$ equioscillates. Although this looks reasonable, it is not true
(neither for $f$ nor for $g$) and we shall show in the following that one can
find a better approximant which does not possess this property (see Figure 1).

Figure 1: Graphs of $f_{s}$ (blue) and $f_{d}$ (orange) for “the optimal”
approximant from [1] (left) and for the optimal approximant constructed in
this paper (right). In all cases $a=1/\sqrt{3}$. It is clearly seen that the
one on the right has a smaller error on $\angle$.

##  3 Main result

To prove that the optimal approximant $f_{s}$ can not equioscillate, let us
start with the following important observation.

###### Lemma 3.1

For every $u\in[0,1)$ and $\alpha,\beta>0$ functions $f_{s}(u,\cdot)$,
$f_{d}(u,\cdot,\beta)$, and $f_{d}(u,\alpha,\cdot)$ are strictly increasing.

###### Proof 1

The result of the lemma follows directly from the observations

$\displaystyle\frac{\partial f_{s}}{\partial\alpha}(u,\alpha)$
$\displaystyle=\left(1-a^{2}\right)\left(1-u^{2}\right)\left(1+u^{2}+2\alpha\left(1-u^{2}\right)\right),$
$\displaystyle\frac{\partial f_{d}}{\partial\alpha}(u,\alpha,\beta)$
$\displaystyle=\frac{1}{2}\sqrt{1-2a^{2}}\beta\left(1-u^{4}\right)\left(1-u^{2}\right)^{2}+2\alpha\left(1-2a^{2}\right)\left(1-u^{4}\right)^{2}$
$\displaystyle+4a^{2}\alpha
u^{2}\left(1-u^{2}\right)^{2}+\frac{1}{2}\left(1-2a^{2}\right)\left(1-u^{8}\right)+u^{2}\left(1-u^{4}\right),$
$\displaystyle\frac{\partial f_{d}}{\partial\beta}(u,\alpha,\beta)$
$\displaystyle=\frac{1}{8}\left(1-u^{2}\right)^{2}\left(\beta\left(1-u^{2}\right)^{2}+\sqrt{1-2a^{2}}\left(1+u^{2}\right)^{2}+4\sqrt{1-2a^{2}}\alpha\left(1-u^{4}\right)\right).$

Quite clearly this lemma implies that $f_{s,d}(\cdot,\alpha^{*},\beta^{*})$
must equioscillate. Moreover, we will show that
$f_{d}(\cdot,\alpha^{*},\beta^{*})$ equioscillates and $f_{s}$ does not, which
contradict the results in [1]. Let us prove a particular relation between
$f_{d}$ and $f_{s}$ first.

###### Lemma 3.2

Let $\alpha\in\left[\tfrac{1}{2},\tfrac{1}{2-2a^{2}}\right]$ and let
$f_{d}(0,\alpha,\beta)\leq f_{s}(0,\alpha)$. Then
$f_{d}(u,\alpha,\beta)<f_{s}(u,\alpha)$ for all $u\in(0,1)$.

###### Proof 2

Let us assume first that $f_{d}(0,\alpha,\beta)=f_{s}(0,\alpha)$. Then

$\beta=\beta(\alpha)=2(1+2\alpha)\sqrt{1-a^{2}}-(1+4\alpha)\sqrt{1-2a^{2}}.$

If we define $\alpha(t)=\frac{a^{2}}{2\left(1-a^{2}\right)}t+\frac{1}{2}$,
then for $t\in[0,1]$ we have
$\alpha(t)\in\left[\tfrac{1}{2},\tfrac{1}{2-2a^{2}}\right]$ and it is eqnough
to show that $h(t)=f_{s}(u,\alpha(t))-f_{d}(u,\alpha(t),\beta(\alpha(t)))>0$
for all $t\in[0,1]$ and $u\in(0,1)$. But this follows directly from

$\displaystyle h(t)$
$\displaystyle=u^{2}\left(1-u^{2}\right)\left(2\left(1-\sqrt{1-2a^{2}}\sqrt{1-a^{2}}-\frac{3a^{2}}{2}\right)\left(1-u^{2}\right)\left(2-u^{2}\right)+a^{2}\right)(1-t)$
$\displaystyle+\frac{u^{2}\left(1-u^{2}\right)^{2}}{1-a^{2}}\left(\left(1-\sqrt{1-2a^{2}}\sqrt{1-a^{2}}-\frac{3a^{2}}{2}\right)\left(4-2u^{2}-a^{2}u^{2}\right)+\frac{a^{4}u^{2}}{2}\right)t(1-t)$
$\displaystyle+\frac{u^{2}\left(1-u^{2}\right)^{2}}{4\left(1-a^{2}\right)^{2}}\left(\left(1-\sqrt{1-2a^{2}}\sqrt{1-a^{2}}-\frac{3a^{2}}{2}-\frac{a^{4}}{8}\right)\left(8-4a^{2}\right)\left(2-a^{2}-u^{2}\right)+\frac{1}{2}a^{6}\left(a^{2}+3u^{2}\right)\right)t^{2}.$

If $f_{d}(0,\alpha,\beta)<f_{s}(0,\alpha)=f_{d}(0,\alpha,\beta(\alpha))$ then
by Lemma 3.1 $\beta<\beta(\alpha)$ and we conclude that
$f_{d}(u,\alpha,\beta)<f_{d}(u,\alpha,\beta(\alpha))\leq f_{s}(u,\alpha)$ for
all $u\in(0,1)$.

Note that

$f_{s}\left(u,\frac{1}{2}\right)=-a^{2}(1-u^{2})\leq 0\quad\text{and}\quad
f_{s}\left(u,\frac{1}{2(1-a^{2})}\right)=\frac{a^{4}(1-u^{2})^{2}}{4(1-a^{2})}\geq
0,$ (4)

thus the parameter $\alpha_{e}$ for which $f_{s}(\cdot,\alpha)$ equioscillates
must be on the interval $I=\left[\tfrac{1}{2},\tfrac{1}{2(1-a^{2})}\right]$.
But unfortunately, it may happen that $\alpha^{*}\not\in I$. Thus we define an
enlarged interval $I_{0}=\left[\tfrac{1}{2},\tfrac{3}{2}\right]\supset I$
which will be used in further analysis.

###### Corollary 3.3

If $f_{s}(\cdot,\alpha_{e})$ equioscillates, then for every $\beta>0$

$\max_{u\in[-1,1]}\left|f_{d}(u,\alpha_{e},\beta)\right|>\max_{u\in[-1,1]}\left|f_{s}(u,\alpha_{e})\right|.$

###### Proof 3

Let us suppose that

$\max_{u\in[-1,1]}\left|f_{d}(u,\alpha_{e},\beta)\right|\leq\max_{u\in[-1,1]}\left|f_{s}(u,\alpha_{e})\right|.$

Since $\max_{u\in[-1,1]}\left|f_{s}(u,\alpha_{e})\right|=f_{s}(0,\alpha_{e})$,
we have $f_{d}(0,\alpha_{e},\beta)\leq f_{s}(0,\alpha_{e})$ and by Lemma 3.2

$f_{d}(u,\alpha_{e},\beta)<f_{s}(u,\alpha_{e}),\quad u\in(0,1).$

Consequently
$\min_{u\in[-1,1]}f_{d}(u,\alpha_{e},\beta)<\min_{u\in[-1,1]}f_{s}(u,\alpha_{e})$
which is a contradiction.

###### Lemma 3.4

If $(\alpha^{*},\beta^{*})$ is an optimal pair of parameters then

$\min_{u\in[-1,1]}f_{s,d}(u,\alpha^{*},\beta^{*})=\min_{u\in[-1,1]}f_{d}(u,\alpha^{*},\beta^{*})\quad\text{and}\quad\max_{u\in[-1,1]}f_{s,d}(u,\alpha^{*},\beta^{*})=\max_{u\in[-1,1]}f_{d}(u,\alpha^{*},\beta^{*}).$

###### Proof 4

Since minimum and maximum of $f_{s,d}$ is attained either by $f_{s}$ or
$f_{d}$, four possibilities have to be considered.
If

$\min_{u\in[-1,1]}f_{s}(u,\alpha^{*})<\min_{u\in[-1,1]}f_{d}(u,\alpha^{*},\beta^{*})\quad{\rm
and}\quad\max_{u\in[-1,1]}f_{d}(u,\alpha^{*},\beta^{*})>\max_{u\in[-1,1]}f_{s}(u,\alpha^{*})$

or

$\min_{u\in[-1,1]}f_{d}(u,\alpha^{*},\beta^{*})<\min_{u\in[-1,1]}f_{s}(u,\alpha^{*})\quad{\rm
and}\quad\max_{u\in[-1,1]}f_{s}(u,\alpha^{*})>\max_{u\in[-1,1]}f_{d}(u,\alpha^{*},\beta^{*})$

then since $f_{s}$ does not depend on $\beta$ and $f_{d}$ is an increasing
function of $\beta$, we can find $\beta^{**}<\beta^{*}$ or
$\beta^{**}>\beta^{*}$ that $(\alpha^{*},\beta^{**})$ implies a better
approximant, respectively.
If
$\max_{u\in[-1,1]}f_{s}(u,\alpha^{*})>\max_{u\in[-1,1]}f_{d}(u,\alpha^{*},\beta^{*})$
then

$\max_{u\in[-1,1]}f_{s}(u,\alpha^{*})=f_{s}(0,\alpha^{*})>\max_{u\in[-1,1]}f_{d}(u,\alpha^{*},\beta^{*})$

and the case
$\min_{u\in[-1,1]}f_{s}(u,\alpha^{*})<\min_{u\in[-1,1]}f_{d}(u,\alpha^{*},\beta^{*})$
is not possible due to the Lemma 3.2. This confirms the result of the lemma.

The previous lemma reveals that we are looking for an optimal error $f_{d}$,
i.e., an equioscillating $f_{d}$. This is still a two-dimensional optimization
problem. But we shall prove that for the optimal pair of parameters
$\alpha^{*}$ and $\beta^{*}$ we must have
$f_{s}(0,\alpha^{*})=f_{d}(0,\alpha^{*},\beta^{*})$ which leads to one-
dimensional optimization. Namely, the quadratic function $\beta\mapsto
f_{d}(0,\alpha,\beta)-f_{s}(0,\alpha)$ has positive leading coefficient and
negative constant term, therefore it possesses the unique positive zero

$\beta=\beta_{0}(\alpha)=2(1+2\alpha)\sqrt{1-a^{2}}-(1+4\alpha)\sqrt{1-2a^{2}}.$

This leads to the following observation.

###### Lemma 3.5

If $f_{s}(0,\alpha)=f_{d}(0,\alpha,\beta)=0$, then $f_{s}(u,\alpha)\leq 0$ and
$f_{d}(u,\alpha,\beta)\leq 0$ for all $u\in[-1,1]$.

###### Proof 5

The solution of the system $f_{s}(0,\alpha)=0$ and $f_{d}(0,\alpha,\beta)=0$
is

$\alpha=\frac{1}{\sqrt{1-a^{2}}}-\frac{1}{2}\quad{\rm
and}\quad\beta=\beta_{0}\left(\frac{1}{\sqrt{1-a^{2}}}-\frac{1}{2}\right)=4+\sqrt{1-2a^{2}}\left(1-\frac{4}{\sqrt{1-a^{2}}}\right).$

Consequently

$f_{s}(u,\alpha)=-\left(1-\sqrt{1-a^{2}}\right)^{2}u^{2}(1-u^{2})\leq 0$

and

$f_{d}(u,\alpha,\beta)=-\frac{2u^{2}\left(1-u^{2}\right)}{1-a^{2}}\left(2\left(1-\sqrt{1-a^{2}}-\frac{a^{2}}{2}\right)u^{2}\left(u^{2}\left(1-a^{2}\right)+a^{2}\left(1-u^{2}\right)\right)\right.\\\
+\left.\left(1-u^{2}\right)^{2}\left(\sqrt{1-2a^{2}}-\sqrt{1-a^{2}}\right)^{2}+\left(1-\sqrt{1-2a^{2}}-a^{2}\right)u^{2}\left(1-u^{2}\right)\left(1-a^{2}\right)\right)\leq
0$

for $u\in[-1,1]$.

If the assumption $f_{s}(0,\alpha)=f_{d}(0,\alpha,\beta)=0$ is extended to
$f_{s}(0,\alpha)=f_{d}(0,\alpha,\beta)\geq 0$, the following result can be
confirmed.

###### Lemma 3.6

Let $f_{d}(0,\alpha,\beta)=f_{s}(0,\alpha)\geq 0$. Then

$\frac{df_{d}}{d\alpha}(0,\alpha,\beta_{0}(\alpha))\geq\frac{df_{d}}{d\alpha}(u,\alpha,\beta_{0}(\alpha))$

for all $u\in[-1,1]$.

###### Proof 6

Recall that $\alpha\geq\frac{1}{2}$ and observe that

$\frac{df_{d}}{d\alpha}(0,\alpha,\beta_{0}(\alpha))-\frac{df_{d}}{d\alpha}(u,\alpha,\beta_{0}(\alpha))=2u^{2}\left(k(\alpha-\tfrac{1}{2})+n\right),$

where

$\displaystyle k=$
$\displaystyle\left(\sqrt{1-2a^{2}}-\sqrt{1-a^{2}}\right)^{2}\left(1-u^{2}\right)\left(2\left(1-u^{2}\right)^{2}+u^{2}+u^{4}\right)$
$\displaystyle+2\left(\sqrt{1-2a^{2}}\sqrt{1-a^{2}}\right)u^{2}\left(1-u^{4}\right)+u^{6}\left(1-a^{2}\right)+a^{2}u^{2}\left(1-u^{2}\right)\left(3-u^{2}\right)\geq
0,$ $\displaystyle n=$ $\displaystyle
2\left(\sqrt{1-2a^{2}}-\sqrt{1-a^{2}}\right)^{2}\left(2-5u^{2}+5u^{4}-\frac{3u^{6}}{2}\right)$
$\displaystyle+\sqrt{1-2a^{2}}\sqrt{1-a^{2}}u^{2}+a^{2}\left(u^{2}-u^{4}+\frac{u^{6}}{2}\right)\geq
0.$

###### Corollary 3.7

If $f_{s}(0,\alpha)=f_{d}(0,\alpha,\beta)\geq 0$ then $\max
f_{s}(\cdot,\alpha)=\max f_{d}(\cdot,\alpha,\beta)=f_{s}(0,\alpha)$.

###### Proof 7

Since

$f_{s}\left(u,\frac{1}{\sqrt{1-a^{2}}}-\frac{1}{2}\right)=-\left(1-\sqrt{1-a^{2}}\right)^{2}u^{2}\left(1-u^{2}\right)\leq
0$

the inequality $f_{s}(0,\alpha)\geq 0$ implies
$\alpha\geq\alpha_{0}=\frac{1}{\sqrt{1-a^{2}}}-\frac{1}{2}$. Since
$\alpha_{0}\in[\frac{1}{2},\frac{1}{2-2a^{2}}]$, Lemma 3.2 implies
$f_{d}(u,\alpha_{0},\beta_{0}(\alpha_{0}))\leq 0$ for all $u\in[-1,1]$. The
result now follows from Lemma 3.6.

###### Lemma 3.8

Let $f_{d}(0,\alpha,\beta)>f_{s}(0,\alpha)\geq 0$. Let $\alpha^{\prime}$ and
$\beta^{\prime}$ be such that
$f_{s}(0,\alpha^{\prime})=f_{d}(0,\alpha,\beta)=f_{s}(0,\alpha^{\prime},\beta^{\prime})$.
Then $f_{d}(u,\alpha,\beta)\leq f_{d}(u,\alpha^{\prime},\beta^{\prime})\leq
f_{d}(0,\alpha,\beta)$ for all $u\in[-1,1]$.

###### Proof 8

For every $\mu\in[f_{s}(0,\alpha),f_{d}(0,\alpha,\beta)]$ let $\alpha(\mu)$ be
the only positive solution of $f_{s}(0,\alpha(\mu))=\mu$, let $\beta(\mu)$ be
the only positive solution of $f_{d}(0,\alpha(\mu),\beta(\mu))=\mu$, and
$\beta_{\alpha}(\mu)$ be the only positive solution of
$f_{d}(0,\alpha,\beta(\mu))=\mu$. Then

$\displaystyle\alpha(\mu)$
$\displaystyle=\frac{2\sqrt{1+\mu}\sqrt{1-a^{2}}-1+a^{2}}{2\left(1-a^{2}\right)},$
$\displaystyle\beta(\mu)$
$\displaystyle=\frac{\left(\sqrt{1-2a^{2}}+4\sqrt{1+\mu}\right)\left(1-a^{2}\right)-4\sqrt{1+\mu}\sqrt{1-2a^{2}}\sqrt{1-a^{2}}}{1-a^{2}},$
$\displaystyle\beta_{\alpha}(\mu)$
$\displaystyle=4\sqrt{1+\mu}-(1+4\alpha)\sqrt{1-2a^{2}}.$

For $\mu_{0}=f_{s}(0,\alpha)$, we have $\alpha=\alpha(\mu_{0})$, hence
$f_{d}(u,\alpha,\beta_{\alpha}(\mu_{0}))=f_{d}(u,\alpha(\mu_{0}),\beta(\mu_{0}))$
for all $u\in[-1,1]$. So it is enough to prove that $\tfrac{\partial
f_{d}}{\partial\mu}(u,\alpha,\beta(\mu))\leq\tfrac{\partial
f_{d}}{\partial\mu}(u,\alpha(\mu),\beta(\mu))$. Since

$\tfrac{\partial f_{d}}{\partial\mu}(u,\alpha(\mu),\beta(\mu))-\tfrac{\partial
f_{d}}{\partial\mu}(u,\alpha,\beta(\mu))=\frac{u^{2}\left(1-u^{2}\right)\left(A_{1}\left(\frac{3}{2}-\alpha\right)+A_{2}\left(\sqrt{1+\mu}-1\right)+A_{3}\right)}{\left(1-a^{2}\right)\sqrt{1+\mu}},$

where

$\displaystyle A_{1}$
$\displaystyle=2\sqrt{1-2a^{2}}\left(1-a^{2}\right)\left(1-u^{2}\right)^{2}\geq
0,$ $\displaystyle A_{2}$
$\displaystyle=2\left(1-u^{2}\right)\left(2\left(1-u^{2}\right)\sqrt{1-2a^{2}}\sqrt{1-a^{2}}+2u^{2}\left(1-2a^{2}\right)+a^{2}\right)\geq
0,$ $\displaystyle A_{3}$
$\displaystyle=2\sqrt{1-2a^{2}}\sqrt{1-a^{2}}2\left(1-u^{2}\right)^{2}\left(1-\sqrt{1-a^{2}}\right)$
$\displaystyle+2\left(\sqrt{1-a^{2}}\left(u^{4}\left(1-2a^{2}\right)+a^{2}u^{2}\right)+\left(1-u^{2}\right)\left(\left(1-2a^{2}\right)2u^{2}+a^{2}\right)\right)\geq
0,$

the lemma follows.

It is clear now that parameters $(\alpha^{*},\beta^{*})$ are such that $|f|$
restricted on $\angle$ has the smallest possible minimum. We will now prove
that $f(\cdot,\cdot,\alpha^{*},\beta^{*})$ has local extrema only on the
diagonals and on the sides of the square $[-1,1]^{2}$.

The function $f(\cdot,\cdot,\alpha,\beta)$ is a symmetric polynomial. Let us
reparameterize it by $e_{1}=u^{2}+v^{2}$ and $e_{2}=u^{2}v^{2}$. The map
$(u,v)\mapsto(u^{2}+v^{2},u^{2}v^{2})=(e_{1},e_{2})$ is a bijection between
$T=\\{(u,v)\in\mathbb{R}^{2};0<u<1,0<v<u\\}$ and
$S=\\{(e_{1},e_{2})\in\mathbb{R}^{2};0<e_{1}<2,e_{1}-1<e_{2}<\tfrac{1}{4}e_{1}^{2},e_{2}>0\\}$.
If we write $h(e_{1},e_{2},\alpha)=f(u,v,\alpha,\beta_{0}(\alpha))$, the map
$f(\cdot,\cdot,\alpha,\beta_{0}(\alpha))$ has local extrema on $T$ if and only
if $h(\cdot,\cdot,\alpha)$ has local extrema on $S$. But
$h(\cdot,\cdot,\alpha)$ as a function on $\mathbb{R}^{2}$ has only one local
extreme point $(x,y)$, namely

$\displaystyle x$
$\displaystyle=\frac{4(1+2\alpha)\left(\alpha(1-2\alpha)\left(2-4\alpha+a^{2}(3+2\alpha)\right)+\left(1+7\alpha+8\alpha^{2}-20\alpha^{3}\right)\left(\sqrt{1-2a^{2}}-\sqrt{1-a^{2}}\right)^{2}\right)}{(1-2\alpha)^{2}(6\alpha+1)\left((14\alpha+5)a^{2}+4\left(\sqrt{1-2a^{2}}\sqrt{1-a^{2}}-1\right)(2\alpha+1)\right)},$
$\displaystyle y$
$\displaystyle=\frac{(1+2\alpha)^{2}\left(a^{2}\left(1-4\alpha^{2}\right)+2\left(5+2\alpha-8\alpha^{2}\right)\left(\sqrt{1-2a^{2}}-\sqrt{1-a^{2}}\right)^{2}\right)}{(1-2\alpha)^{2}(1+6\alpha)\left(a^{2}(2\alpha-1)-2(1+2\alpha)\left(\sqrt{1-2a^{2}}-\sqrt{1-a^{2}}\right)^{2}\right)}.$

We shall see that for $\alpha=\alpha^{*}$ we have $y>1$, which implies that
$h(\cdot,\cdot,\alpha^{*})$ has no local extrema on $T$.

Let us first prove that
$\tfrac{1-\tfrac{1}{5}a^{4}}{2(1-a^{2})}=:\alpha_{l}<\alpha^{*}<\tfrac{1+\tfrac{1}{5}a^{4}}{2(1-a^{2})}=:\alpha_{r}$.
Since
$f_{s}(0,\alpha_{l})=\frac{a^{4}\left(5+10a^{2}+a^{4}\right)}{100\left(1-a^{2}\right)}\geq
0$, the functions $f_{s}(\cdot,\cdot,\alpha)$ and
$f_{d}(\cdot,\cdot,\alpha,\beta_{0}(\alpha))$ have the maximum at $u=0$ for
all $\alpha\in[\alpha_{l},\alpha_{r}]$. We can write

$\displaystyle
f_{d}\left(\frac{3}{4},\alpha_{l},\beta_{0}(\alpha_{l})\right)+f_{d}\left(0,\alpha_{l},\beta_{0}(\alpha_{l})\right)=-\frac{\eta}{6553600(1-a^{2})^{2}(\zeta\sqrt{1-a^{2}}\sqrt{1-2a^{2}}+\xi)}\leq
0,$

where

$\displaystyle\eta=a^{4}\sum_{i=0}^{8}\eta_{i}a^{2i}(1-2a^{2})^{8-i},\quad\zeta=\sum_{i=0}^{4}\zeta_{i}a^{2i}(1-2a^{2})^{4-i},\quad\xi=\sum_{i=0}^{5}\xi_{i}a^{2i}(1-2a^{2})^{5-i},$

such that $\eta_{i},\zeta_{i},\xi_{i}>0$ for all $i$. Hence the minimum of
$\left|f_{d}\left(\cdot,\cdot,\alpha_{l},\beta_{0}(\alpha_{l})\right)\right|$
is larger then the maximum, therefore $\alpha^{*}>\alpha_{l}$. On the other
hand, we can write

$\displaystyle
f_{d}\left(u,\alpha_{r},\beta_{0}(\alpha_{r})\right)+f_{d}\left(0,\alpha_{r},\beta_{0}(\alpha_{r})\right)=\frac{\zeta\sqrt{1-a^{2}}\sqrt{1-2a^{2}}+\xi}{100(1-a^{2})^{2}},$

where
$\zeta=4\left(10-5a^{2}+a^{4}\right)u^{2}\left(1-u^{2}\right)^{2}\left(5\left(2-a^{2}-u^{2}\right)+a^{4}\left(1-u^{2}\right)\right)\geq
0$. Since $\sqrt{1-a^{2}}\sqrt{1-2a^{2}}\geq 1-\tfrac{3}{2}a^{2}-a^{4}$ for
all $a\in[0,\tfrac{\sqrt{2}}{2}]$, we have

$\displaystyle
f_{d}\left(u,\alpha_{r},\beta_{0}(\alpha_{r})\right)+f_{d}\left(0,\alpha_{r},\beta_{0}(\alpha_{r})\right)=\frac{\zeta\sqrt{1-a^{2}}\sqrt{1-2a^{2}}+\xi}{100(1-a^{2})^{2}}\geq\frac{\zeta(1-\tfrac{3}{2}a^{2}-a^{4})+\xi}{100(1-a^{2})^{2}}.$

Since we can write

$\displaystyle\zeta(1-\tfrac{3}{2}a^{2}-a^{4})+\xi=\sum_{i=0}^{4}\left(\sum_{j=0}^{9-i}k_{i,j}u^{2j}(1-u^{2})^{9-i-j}\right)a^{12-2i}(1-2a^{2})^{i},$

where all coefficients $k_{i,j}>0$, we have
$f_{d}\left(u,\alpha_{r},\beta_{0}(\alpha_{r})\right)+f_{d}\left(0,\alpha_{r},\beta_{0}(\alpha_{r})\right)\geq
0$ for all $u\in[-1,1]$. Hence the maximum of
$\left|f_{d}\left(\cdot,\cdot,\alpha_{r},\beta_{0}(\alpha_{r})\right)\right|$
is larger then the minimum, therefore $\alpha^{*}<\alpha_{r}$.

Let us write $\alpha=(1-t)\alpha_{l}+t\alpha_{r}$ and prove that $y>1$ for all
$t\in[0,1]$, in particular, $y>1$ for $\alpha=\alpha^{*}$. Let us write
$y=\tfrac{y_{n}}{y_{d}}$. Let us first show that $y_{d}<0$. We can write

$-\frac{5(1-2a^{2})}{(1-2\alpha)^{2}(1+6\alpha)}y_{d}=\xi+\zeta\sqrt{1-2a^{2}}\sqrt{1-a^{2}},$

where
$\xi=40(1-2a^{2})^{3}+160a^{2}(1-2a^{2})^{2}+181a^{4}(1-2a^{2})+49a^{6}+8a^{4}(1-2a^{2})t+2a^{6}t$
and $\zeta=4(a^{4}(1-2t)-5(2-a^{2})$. Since $\xi>0$ for all $t\in[0,1]$ it is
enough to show that
$\xi^{2}-\left(\zeta\sqrt{1-2a^{2}}\sqrt{1-a^{2}}\right)^{2}>0$. This follows
from the equality

$\xi^{2}-\left(\zeta\sqrt{1-2a^{2}}\sqrt{1-a^{2}}\right)^{2}=17a^{12}(1-2t)^{2}+480a^{6}(1-a^{2})(1-t)+320ta^{6}+12a^{10}t+\\\
25a^{8}(1-2a^{2})t+105a^{8}(1-2a^{2})(1-t)+232a^{10}(1-t)+32ta^{10}(1-t)>0.$

Hence to prove $y>1$, we have to show that $y_{d}-y_{n}>0$. We can write

$\displaystyle y_{d}-y_{n}$
$\displaystyle=\xi\sqrt{1-a^{2}}\sqrt{1-2a^{2}}-\zeta=\frac{(\xi\sqrt{1-a^{2}}\sqrt{1-2a^{2}}+k)^{2}-(\zeta+k)^{2}}{\xi\sqrt{1-a^{2}}\sqrt{1-2a^{2}}+\zeta+2k}$
$\displaystyle\geq\frac{2k\xi\sqrt{1-a^{2}}\sqrt{1-2a^{2}}+\xi^{2}(1-a^{2})(1-2a^{2})-\zeta^{2}-2k\zeta}{\xi\sqrt{1-a^{2}}\sqrt{1-2a^{2}}+\zeta+2k}$
$\displaystyle\geq\frac{2k\xi(1-\frac{3}{2}a^{2}-\frac{1}{8}a^{4}-\frac{7}{4}a^{6})+\xi^{2}(1-a^{2})(1-2a^{2})-\zeta^{2}-2k\zeta}{\xi\sqrt{1-a^{2}}\sqrt{1-2a^{2}}+\zeta+2k},$

where $k=\tfrac{103901}{256}$ and

$\displaystyle\xi$
$\displaystyle=\frac{1}{625(1-a^{2})^{4}}\left(\sum_{i=0}^{15}\sum_{j=0}^{3}\xi_{i,j}a^{2i}(1-2a^{2})^{15-i}t^{j}(1-t)+\sum_{i=0}^{15}\xi_{i,4}a^{2i}(1-2a^{2})^{15-i}t^{4}\right),$
$\displaystyle\zeta$
$\displaystyle=\frac{1}{625(1-a^{2})^{4}}\left(\frac{1}{256}\left(\sum_{i=0}^{15}\sum_{j=0}^{3}\zeta_{i,j}a^{2i}(1-2a^{2})^{15-i}t^{j}(1-t)+\sum_{i=0}^{15}\zeta_{i,4}a^{2i}(1-2a^{2})^{15-i}t^{4}\right)+2k\right),$
$\displaystyle 2k\xi$
$\displaystyle(1-\frac{3}{2}a^{2}-\frac{1}{8}a^{4}-\frac{7}{4}a^{6})+\xi^{2}(1-a^{2})(1-2a^{2})-\zeta^{2}-2k\zeta$
$\displaystyle=\sum_{i=0}^{15}\sum_{j=0}^{8}\eta_{i,j}a^{2i+6}(1-2a^{2})^{15-i}t^{j}(1-t)+\sum_{i=0}^{15}\eta_{i,9}a^{2i+6}(1-2a^{2})^{15-i}t^{9},$

and all coefficients $\xi_{i,j}$, $\zeta_{i,j}$ and $\eta_{i,j}$ are positive.
Hence $y>1$ and the maximum and minimum of $f(\cdot,\alpha^{*},\beta^{*})$ are
on the diagonals or on the sides of the square $[-1,1]^{2}$. Thus we have
proved the following crucial theorem.

###### Theorem 3.9

The maximum and the minimum of $f(\cdot,\cdot,\alpha,\beta)$ are on the set
$\angle$.

It is quite clear now how to find the best approximant. The procedure is
identical for the simplified radial error and for the radial error by
replacing $f$ by $g$ so we consider only the first one.

  * a)

Solve the equation $f_{d}(0,\alpha,\beta)=f_{s}(0,\alpha)$ on $\beta$ to get
the admissible $\beta_{0}(\alpha):=\beta(\alpha)$.

  * b)

Solve the system of nonlinear equations

$\displaystyle\frac{\partial f_{d}}{\partial u}(u,\alpha,\beta_{0}(\alpha))$
$\displaystyle=0,$ $\displaystyle
f_{d}(0,\alpha,\beta_{0}(\alpha))+f_{d}(u,\alpha,\beta_{0}(\alpha))$
$\displaystyle=0$

on $u$ and $\alpha$ to get the admissible solution $u_{0}$ and $\alpha_{0}$.

  * c)

Return the best approximant $\boldsymbol{p}$ constructed by $\alpha_{0}$ and
$\beta_{0}(\alpha_{0})$.

Note that the above algorithm must be performed numerically since the equation
in b) can not be solved in a closed form in general. The Newton-Raphson method
might be used or any other appropriate algorithm for solving the system of
nonlinear equations.

##  4 Optimal $G^{1}$ approximation

A natural question arises whether the optimal approximant of the spherical
square constructed in the previous section can be used as a $G^{1}$ tensor
product quadratic Bézier spline approximant of the unit sphere. In this case,
either two or six patches can be put together. We shall focus only on the
spline of six patches since even in this case the obtained $G^{1}$ surface is
not a good approximation of the sphere.
If six patches are put together we have $a=\tfrac{\sqrt{3}}{3}$. The $G^{1}$
condition implies that the tangent plane of each of the three patches meeting
at the common corner point must coincide with the tangent plane of a sphere.
Hence $\tfrac{\partial\boldsymbol{p}}{\partial
u}(-1,-1)\times\tfrac{\partial\boldsymbol{p}}{\partial
v}(-1,-1)=[\tfrac{2}{3}(1-2\alpha)\alpha,\tfrac{2}{3}(1-2\alpha)\alpha,\tfrac{4}{3}(1-\alpha)\alpha]^{T}$
must be parallel to the unit normal of $\boldsymbol{p}$ at $(-1,-1)$, i.e.,
parallel to
$[-\tfrac{\sqrt{3}}{3},-\tfrac{\sqrt{3}}{3},\tfrac{\sqrt{3}}{3}]^{T}$. This
implies $\alpha_{G}=\tfrac{3}{4}$. For every $v\in[0,1]$ the tangent plane at
points $\boldsymbol{p}(1,v)$ and $R\boldsymbol{p}(0,v)$, where $R$ is the
rotation matrix around $y$-axis for the angle $\tfrac{\pi}{2}$, must coincide.
Due to the symmetry, the normal vector of the tangent plane must lie in the
plane defined by $\boldsymbol{b}_{20}$, $\boldsymbol{b}_{22}$, and the origin.
Therefore $\tfrac{\partial\boldsymbol{p}}{\partial
u}(-1,v)\times\tfrac{\partial\boldsymbol{p}}{\partial v}(-1,v)$ must be
perpendicular to $\boldsymbol{p}(1,-1)\times\boldsymbol{p}(1,1)$ which implies
$\beta_{G}=\tfrac{7\sqrt{3}}{6}$. So the $G^{1}$ spline approximant is unique
and its radial error is
$\left|g(0,0,\alpha_{G},\beta_{G})\right|=\tfrac{5\sqrt{3}-8}{8}\approx
0.0825$ which is much bigger than the radial error of the corresponding
$G^{0}$ spline of six optimal approximants constructed in the previous section
(see Figure 4 and Figure 5).

##  5 Rectangular case

Another interesting issue is the optimal approximation of a spherical
rectangle. Recall that in the case of the approximation of a spherical square,
the minimum of the simplified radial error was never attained at the boundary
of a patch. The situation is quite different in the case of the approximation
of a spherical rectangle. The explanation will be supported by numerical
evidence but no formal proof will be provided.
Let the projection of the vertices of a spherical rectangle along the vector
$[0,0,1]^{T}$ be vertices of a planar rectangle with its larger edge $2a$
fixed and let $2b<2a$ be it shorter edge. Similarly, as in the case of a
spherical square (1), we define a tensor product quadratic Bézier patch
$\boldsymbol{p}_{r}$ as

$\boldsymbol{p}_{r}(u,v,\alpha_{1},\alpha_{2},\beta)=\sum_{i=0}^{2}\sum_{j=0}^{2}B_{i}^{2}(u)B_{j}^{2}(v)\boldsymbol{c}_{ij},\quad
u,v\in[-1,1],$

where, analogously as in (3),

$\displaystyle\boldsymbol{c}_{00}=(-a,-b,\sqrt{1-a^{2}-b^{2}}),\quad$
$\displaystyle\boldsymbol{c}_{10}=\alpha_{1}(\boldsymbol{c}_{00}+\boldsymbol{c}_{20}),\quad$
$\displaystyle\boldsymbol{c}_{20}=(a,-b,\sqrt{1-a^{2}-b^{2}}),\quad$
$\displaystyle\boldsymbol{c}_{01}=\alpha_{2}(\boldsymbol{c}_{00}+\boldsymbol{c}_{02}),\quad$
$\displaystyle\boldsymbol{c}_{11}=\beta\,(0,0,1),\quad$
$\displaystyle\boldsymbol{c}_{21}=\alpha_{2}(\boldsymbol{c}_{20}+\boldsymbol{c}_{22}),\quad$
(5) $\displaystyle\boldsymbol{c}_{02}=(-a,b,\sqrt{1-a^{2}-b^{2}}),\quad$
$\displaystyle\boldsymbol{c}_{12}=\alpha_{1}(\boldsymbol{c}_{02}+\boldsymbol{c}_{22}),\quad$
$\displaystyle\boldsymbol{c}_{22}=(a,b,\sqrt{1-a^{2}-b^{2}}),$

with some $\alpha_{1},\alpha_{2},\beta>0$. As in (2), we define
$f_{r}(\cdot,\cdot,\alpha_{1},\alpha_{2},\beta)=\left\|\boldsymbol{p}_{r}(\cdot,\cdot,\alpha_{1},\alpha_{2},\beta)\right\|^{2}-1.$
We have seen in the analysis of the square case that the global minimum of $f$
can not appear on the boundary of the square $[-1,1]^{2}$. We shall see in the
following that this can actually happen in the rectangular case for some
particular ratio of sides $a$ and $b$. Once the global minimum is on the
boundary, numerical examples confirm that several optimal approximants of the
spherical rectangle might exist. This indicates that the rectangular case is a
much more challenging problem.
Let us now try to find the condition on $b$, which implies that the global
minimum and maximum of $f$ are on the boundary. If this is the case, we can
assume that the restriction of $f_{r}$ on $(u,-1)$, $u\in[-1,1]$,
equioscillate. This implies the parameter $\alpha_{1}$ and the point
$u_{0}\in(-1,0)$ such that $\tfrac{\partial f_{r}}{\partial
u}(u_{0},-1,\alpha_{1},\alpha_{2},\beta)=0$ and
$f_{r}(u_{0},-1,\alpha_{1},\alpha_{2},\beta)$ is the global minimum of $f_{r}$
for any $\alpha_{2}$ and $\beta$. A necessary condition for
$f_{r}(u_{0},-1,\alpha_{1},\alpha_{2},\beta)$ being the global minimum is
$\tfrac{\partial f_{r}}{\partial v}(u_{0},-1,\alpha_{1},\alpha_{2},\beta)\geq
0$. Suppose that the triple $(b_{a},\alpha_{a},\beta_{a})$ is the solution of
the system of nonlinear equations

$\displaystyle\tfrac{\partial f_{r}}{\partial
v}(u_{0},-1,\alpha_{1},\alpha_{2},\beta)=0,$ $\displaystyle
f_{r}(0,-1,\alpha_{1},\alpha_{2},\beta)=f_{r}(0,0,\alpha_{1},\alpha_{2},\beta),$
$\displaystyle
f(0,-1,\alpha_{1},\alpha_{2},\beta)=f(-1,0,\alpha_{1},\alpha_{2},\beta),$

with respect to $b$, $\alpha_{2}$ and $\beta$. For every $b=b_{a}$, the
candidate for the best approximant of the spherical rectangle related to $a$
and $b_{a}$ is determined by the parameters
$(\alpha_{1},\alpha_{a},\beta_{a})$. If $b<b_{a}$ all pairs of parameters
$(\alpha_{2},\beta)$ satisfying the inequalities

$\displaystyle\tfrac{\partial f}{\partial
v}(u_{0},-1,\alpha_{1},\alpha_{2},\beta)$ $\displaystyle\geq 0,$
$\displaystyle f(0,-1,\alpha_{1},\alpha_{2},\beta)$ $\displaystyle\geq
f(0,0,\alpha_{1},\alpha_{2},\beta),$ (6) $\displaystyle
f(0,-1,\alpha_{1},\alpha_{2},\beta)$ $\displaystyle\geq
f(-1,0,\alpha_{1},\alpha_{2},\beta),$

might induce the best approximant. Thus, in this case, we might have an
uncountable number of optimal polynomial approximants of the spherical
rectangle. Indeed, let $\boldsymbol{v}_{i}=(\alpha_{2,i},\beta_{i})$,
$i=1,2,3$, be the intersections of two pairs of curves in (6) which are
implicitly defined by taking equalities instead of inequalities. Numerical
computations indicate that all pairs $(\alpha_{2},\beta)$ from the triangle
defined by $\boldsymbol{v}_{i}$, $i=1,2,3$, imply an optimal approximant.
Moreover, numerical computations also indicate that several optimal
approximant might exist if a pair $(a,b)$ is taken from the grey region in
Figure 2. On Figure 3, radial errors of three optimal approximants of the
spherical rectangle with $a=0.75$ and $b=0.2$ are shown.

Figure 2: The domain (light grey) in $[0,1]^{2}$ for which pairs $(a,b)$ might
imply several optimal approximants of the spherical rectangle. The grey dashed
curve is the graph of the function $\min\\{a,\sqrt{1-a^{2}}\\}$ and together
with the $a$-axis determine the admissible domain of pairs $(a,b).$ Figure 3:
Graphs of $f$ for three optimal approximants of the spherical rectangle given
by parameters $a=0.75$ and $b=0.2$. Here $\alpha_{1}=1.0277$ and the left
graph corresponds to the pair $(\alpha_{2},\beta)=(0.5313,1.4456)$, the middle
one to $(\alpha_{2},\beta)=(0.5313,1.3881)$, and the right one to
$(\alpha_{2},\beta)=(0.5239,1.4550)$.

##  6 Numerical examples

Some numerical examples will be presented in this section, confirming proven
theoretical results. As the first example, let us consider the optimal
approximation of the unit sphere by the quadratic Bézier tensor product spline
approximant. It is easy to see that only the spline of two patches, each
approximating the hemisphere, and the spline of six patches, each
approximating the one-sixth of the unit sphere, are possible. Their plots and
graphs of radial errors are in Figure 4.

Figure 4: Optimal approximation of the unit sphere by the spline of two tensor
product quadratic Bézier patches (upper left) and by the spline of six tensor
product quadratic Bézier patches (upper right) together with the graphs of the
corresponding radial errors for the single patch (bottom). The colours of the
approximants indicate the distance from the sphere (red regions are out of the
sphere).

In Table 1, optimal parameters according to the radial error, radial distances
and the numerical convergence rates are collected for a set of chosen
parameters $a$. The numerical rate of convergence is estimated as follows. Let
$d_{1}$ and $d_{2}$ be two consecutive radial distances according to the
parameters $a_{1}$ and $a_{2}$. Assuming that the radial distance is of the
form $d={\rm const}\,a^{r}$, we can estimate
$r\approx\log\left(d_{1}/d_{2}\right)/\log\left(a_{1}/a_{2}\right)$. It is
clearly seen that the distance converges to zero as the square of the area of
the spherical square.

$a$ | $(\alpha_{r}^{*},\beta_{r}^{*})$ | $|g|$ | $r$
---|---|---|---
$a_{max}$ | $(1.0306,4.3393)$ | $8.2331\times 10^{-2}$ | $-$
$a_{max}/2$ | $(0.5698,1.1630)$ | $6.9966\times 10^{-4}$ | $6.87$
$a_{max}/4$ | $(0.5160,1.0333)$ | $3.7421\times 10^{-5}$ | $4.23$
$a_{max}/8$ | $(0.5039,1.0079)$ | $2.2596\times 10^{-6}$ | $4.05$
$a_{max}/16$ | $(0.5010,1.0020)$ | $1.4005\times 10^{-7}$ | $4.01$
$a_{max}/32$ | $(0.5002,1.0005)$ | $8.7349\times 10^{-9}$ | $4.00$
$a_{max}/64$ | $(0.5001,1.0001)$ | $5.4565\times 10^{-10}$ | $4.00$
Table 1: Optimal parameters $\alpha_{r}^{*}$ and $\beta_{r}^{*}$ according to
the radial error $g$, the corresponding radial distance $|g|$ and the
numerical rate of convergence $r$ for several parameters $a=a_{max}/2^{i}$,
$i=0,1,\dots,6$, with $a_{max}=1/\sqrt{2}$.

The approximation of the unit sphere by the $G^{1}$ spline of six tensor
product quadratic Bézier patches constructed in Section 4 is given in Figure
5. It is clearly seen that its radial error is much bigger than the radial
error of the optimal approximant in Figure 4. Moreover, the $G^{1}$
approximant is one-sided, i.e., the whole approximant is out of the sphere.
This suggests that omitting the interpolation conditions at the vertices of
the spherical square would imply better approximant by pulling it towards the
origin.

Figure 5: Approximation of the unit sphere by the $G^{1}$ spline of six tensor
product quadratic Bézier patches (left) together with the graph of the radial
error for the single patch (right). The colours of the spline approximant
indicate its distance from the sphere. The red indicates a bigger distance.

##  7 Conclusion

Finding the optimal polynomial approximant of a given surface is a challenging
nonlinear optimization problem. There are only a few references dealing with
this problem available. In this paper, we have shown that the results obtained
in [1] are not correct. As a counterexample, we found a better approximant of
the spherical square and provided an efficient algorithm for its construction.
It is natural to consider higher degree polynomial approximants of spherical
squares or at least approximants providing smoother polynomial spline patches
($G^{k}$ continuous tensor product spline patches). Both problems might be
considered future work, but they lead to much more complicated nonlinear
optimisation issues. On the other hand, the approximation of spherical
rectangles by tensor product quadratic patches might be of some interest.
Preliminary results reveal that, in some cases, it leads to several
(infinitely many) optimal solutions. Thus the square and the rectangular case
deeply differ in their nature.

Acknowledgments. The first author was supported by the Slovenian Research
Agency program P1-0292 and the grant J1-4031. The second author acknowledges
financial support from the Slovenian Research Agency program grant P1-0288 and
the grants N1-0137 and J1-3005.

## References

  * Eisele [1994] Eisele, E.F., 1994.  Best approximations of symmetric surfaces by biquadratic Bézier surfaces.  Comput. Aided Geom. Design 11, 331–343.
  * Floater [1997] Floater, M.S., 1997.  An $O(h^{2n})$ Hermite approximation for conic sections.  Comput. Aided Geom. Design 14, 135–151.
  * Grošelj and Šadl Praprotnik [2022] Grošelj, J., Šadl Praprotnik, A., 2022\.  Exact sphere representations over Platonic solids based on rational multisided Bézier patches.  Comput. Aided Geom. Design 98, Paper No. 102148, 17.
  * Vavpetič [2020] Vavpetič, A., 2020.  Optimal parametric interpolants of circular arcs.  Comput. Aided Geom. Design 80, Paper No. 101891, 9.
  * Vavpetič and Žagar [2022] Vavpetič, A., Žagar, E., 2022\.  Geometric approximation of the sphere by triangular polynomial spline patches.  Comput. Aided Geom. Design 92, Paper No. 102061, 12.