[EN]

[EN]

# Topology Optimization for Uniform Flow Distribution in Electrolysis Cells

Leon Baeck111Corresponding author: e-mail<EMAIL_ADDRESS>11
[DE]Fraunhofer Institute for Industrial Mathematics ITWM, Kaiserslautern,
Germany  Sebastian Blauth 11 Christian Leithäuser 11 René Pinnau 22 [DE]RPTU
Kaiserslautern-Landau, Technomathematics Group, Kaiserslautern, Germany  Kevin
Sturm 33 [AT]TU Wien, Institute of Analysis and Scientific Computing, Vienna,
Austria

###### Abstract

In this paper we consider the topology optimization for a bipolar plate of a
hydrogen electrolysis cell. We present a model for the bipolar plate using the
Stokes equation with an additional drag term, which models the influence of
fluid and solid regions. Furthermore, we derive a criterion for a uniform flow
distribution in the bipolar plate. To obtain shapes that are well-
manufacturable, we introduce a novel smoothing technique for the fluid
velocity. Finally, we present some numerical results and investigate the
influence of the smoothing on the obtained shapes.

##  1 Introduction

Hydrogen plays an important role for achieving a climate neutral industry as
well as climate neutral means of transportation. For that the hydrogen needs
to be produced, for example, by means of hydrogen electrolysis cells.
Typically, those cells split water into oxygen and the desired hydrogen by
using (green) electrical energy. We focus on a special kind of electrolysis
cell here, the so called proton exchange membrane (PEM) electrolysis cell [1],
which works as follows: First, water enters the cell on the anode side and is
distributed over the cell by the bipolar plate. The water wanders through the
so-called porous transport layer (PTL) to the membrane where the
electrochemical reaction takes place. The oxygen is then discharged from the
cell again by the bipolar plate. Whereas positive hydrogen ions can pass
through the membrane and form hydrogen at the cathode side. To guarantee the
efficiency of such a PEM electrolysis cell, the flow through the anode side
bipolar plate is essential. The water needs to be distributed uniformly over
the whole cell so that the entire membrane area is reached. For more details
on PEM electrolysis cells, we refer the reader, e.g., to [1].

Here, we only investigate the anode side of the bipolar plate. For simplicity,
we neglect the electrochemical reaction and phase changes and only consider
the distribution of a homogeneous fluid. To obtain a uniform flow through the
bipolar plate, we use techniques from topology optimization based on
topological sensitivity analysis [2].

Topology optimization considers the optimization of a domain by changing its
topological features by either adding or removing material. It was first
introduced in the context of solid mechanics, but has since been applied to a
wide variety of applications, for example, in elasticity [3, 4] and fluid
mechanics [5, 6]. We use the concept of the topological derivative, which was
first introduced in [3] and later mathematically justified in [2]. The
topological derivative measures the sensitivity of a shape functional under
infinitesimal topological changes. There are different approaches for solving
topology optimization problems numerically and here we utilize a level-set
approach which was introduced in [4].

This paper is structured as follows. In Section 2, we introduce a model for
the anode side bipolar plate using the Borvall-Petersson model from [5].
Additionally, we derive a criterion for a uniform flow distribution in the
bipolar plate. To obtain geometries suitable for manufacturing, we consider a
smoothing of the fluid velocity which is incorporated in the cost functional
for the optimization. In Section 3, we summarize some basic concepts of
topology optimization including the topological derivative as well as a
gradient-based algorithm from [4]. Further, we formally present the
topological derivative of the cost functional for our problem. Finally, in
Chapter 4, we numerically solve the topology optimization problem and present
our numerical results. Particularly, we focus on the influence of the
smoothing on the optimized geometries. Our results show that topology
optimization can be used to create novel bipolar plate designs.

##  2 Model problem for the anode side bipolar plate

Let us begin with presenting our mathematical model for the anode side of the
bipolar plate as well as the corresponding topology optimization problem for a
uniform flow distribution.

###  2.1 The Borvall-Petersson model

First, let the hold-all domain $D\subset\mathbb{R}^{d}$ be an open and bounded
set. In this paper, we focus on the case $d=2$, but an extension to $d=3$ is
also possible. Additionally, let $\Omega\subset D$ be measurable and open.
Then, we identify $\Omega$ as the fluid region and $D\setminus\bar{\Omega}$
consequently as the solid part of $D$. To model the fluid flow, we use the
Borvall Petersson model, which was first introduced in [5] for topology
optimization of fluids. We assume that the boundary $\Gamma=\partial D$ of the
hold-all domain $D$ is divided into three parts: The inlet
$\Gamma_{\text{in}}$, where the water enters the plate with velocity
$u_{\text{in}}$, the boundary $\Gamma_{\text{wall}}$, where we have a no-slip
condition and $\Gamma_{\text{out}}$, where a natural outflow condition is
applied. The model reads

$\begin{split}-\Delta u+\alpha u+\nabla p&=0&&\text{in }D,\\\
\text{div}(u)&=0&&\text{in }D,\\\ u&=u_{\text{in}}&&\text{on
}\Gamma_{\text{in}},\\\ u&=0&&\text{on }\Gamma_{\text{wall}},\\\
\partial_{n}u-pn&=0&&\text{on }\Gamma_{\text{out}},\end{split}$ (1)

where $u$ and $p$ denote the fluid velocity and pressure, respectively. Here,
$n$ denotes the outer unit normal vector on $\Gamma$. To distinguish between
the solid and fluid part of the domain the inverse permeability $\alpha$ is
used, which is given by

$\alpha(x)=\begin{cases}\alpha_{U}&\text{if $x\in D\setminus\bar{\Omega}$},\\\
\alpha_{L}&\text{if $x\in\Omega$},\end{cases}$

where $\alpha_{L}$ and $\alpha_{U}$ are positive constants. In particular,
$\alpha$ is chosen small inside the fluid region and large in the solid part.
For a detailed derivation of the model we refer to [5].

###  2.2 A criterion for uniform flow distribution

To increase the efficiency of the electrolysis cell, we want to assure that
the water gets distributed uniformly all over the cell. Without knowledge of
the shape of $\Omega$, it is not obvious how to characterize a flow as uniform
mathematically. Therefore, we take a different approach and introduce a
threshold velocity magnitude $u_{t}>0$. The goal is now to achieve a flow
whose magnitude reaches at least the threshold. This criterion ensures that
each part of the bipolar plate receives a sufficiently large flow and no dead
spots occur. Additionally, as our model also yields, albeit small, flow
velocities in the solid part, we have to restrict our criterion to the fluid
part $\Omega$. The optimization goal then reads

$\left\lVert u(x)\right\rVert\geq u_{t}\text{ for all }x\in\Omega,$ (2)

where $\left\lVert\cdot\right\rVert$ denotes the Euclidean norm on
$\mathbb{R}^{d}$.

There are two aspects to consider here. First, the target velocity constraint
might not be fulfillable close to the boundary $\partial\Omega$ of the fluid
domain due to small velocities in the solid part. Additionally, as we do not
have any restrictions on the solid part, large solid areas might appear in the
design structures. In fact, we observed this behavior in numerical tests.
Usually, we want to avoid solid structures above a certain size as the flow
inside the porous transport layer is hindered by such objects and this leads
to a degradation of the cell efficiency.

###  2.3 Smoothing of the velocity

To avoid large solid structures and diminish the effect of small flow
velocities in the solid parts, our goal is to extend the target velocity goal
$(\ref{optgoal1})$ to the hold-all domain $D$. We consider a smoothed flow
velocity which is defined on the entire hold-all domain D. For that we use the
smoothing characteristics of the heat equation. We consider the following
equation

$\begin{split}\partial_{t}u_{s}-\Delta u_{s}&=0&&\text{in }D,\\\ \nabla
u_{s}\cdot n&=0&&\text{on }\partial D,\\\ u_{s}(x,0)&=u(x)&&\text{in
}D.\end{split}$ (3)

Here, we have the velocity solution $u$ of the Stokes Darcy equation
$(\ref{stokesdarcy})$ as the initial condition. To reduce the numerical
effort, we discretize $(\ref{heat})$ with a single implicit Euler step with
step length $\Delta t>0$ and arrive at

$\begin{split}\frac{u_{s}-u}{\Delta t}-\Delta u_{s}&=0&&\text{in }D,\\\ \nabla
u_{s}\cdot n&=0&&\text{on }\partial D.\end{split}$ (4)

By modifying the time step length $\Delta t$ we can control the influence of
the smoothing, where the smoothing effect gets larger for increasing step
lengths. With the smoothed velocity $u_{s}$, we can now extend the criterion
$(\ref{optgoal1})$ to the entire hold-all domain by considering

$\left\lVert u_{s}(x)\right\rVert\geq u_{t}\text{ for all }x\in D.$ (5)

Finally, we apply the Moreau-Yosida regularization [7] to $(\ref{optgoal2})$
to arrive at our objective function

$J(\Omega,u)=\int_{D}\min\left(0,\left\lVert u_{s}\right\rVert-
u_{d}\right)^{2}\mathop{}\\!\mathrm{d}x.$ (6)

######  2.1

Using the smoothed velocity has several advantages for the topology
optimization. First, when we consider a small solid inclusion, we expect that
it does not severely effect the flow in the underlying porous transport layer.
In the smoothed velocity $u_{s}$ with a suitable time step $\Delta t$, the
effect of the small solid inclusion will not be resolved anymore, so that the
smoothed velocity $u_{s}$ is not influenced anymore by the small velocity
inside the solid region and the target velocity goal $(\ref{optgoal2})$ will
be fulfilled. On the other hand, if we consider a large solid inclusion, it
will still be resolved by the smoothed velocity $u_{s}$ and the criterion
$(\ref{optgoal2})$ will not be satisfied. Therefore, large solid structures
can be avoided with this approach.

###  2.4 Topology optimization problem

We summarize $(\ref{stokesdarcy})$, $(\ref{heat1step})$ and
$(\ref{objective})$ to state our topology optimization problem

$\begin{split}\min_{\Omega}J(\Omega,u)&=\int_{D}\min\left(0,\left\lVert
u_{s}\right\rVert-u_{t}\right)^{2}\mathop{}\\!\mathrm{d}x\\\ &\text{such that
}(\ref{stokesdarcy})\text{ as well as }(\ref{heat1step})\text{ are fulfilled
in the weak sense},\\\ &V_{L}\leq|\Omega|\leq V_{U}.\end{split}$ (7)

Here, we introduce an additional volume constraint for the fluid volume with
$V_{L}$ and $V_{U}$ being the lower and the upper border, respectively. The
goal is now to perform a topology optimization for problem
$(\ref{optproblem})$ using the topological derivative.

##  3 Topology optimization

We seek to apply a gradient-based topology optimization approach for
$(\ref{optproblem})$. For that we introduce the topological derivative and
state the corresponding topological derivative for the objective
$(\ref{objective})$. Furthermore, we briefly present the gradient-based
topology optimization algorithm we use for our numerical experiments.

###  3.1 Topological derivative

We consider an analogous setting to before, where $D$ is an open subset of
$\mathbb{R}^{d}$. We define $\mathcal{A}=\\{\Omega\subset D\text{ s.t.
}\Omega\text{ open}\\}$ as the set of all open subsets of $D$. To present the
topological derivative, we consider a general shape functional $S$, given by

$\Omega\in\mathcal{A}\mapsto S(\Omega)\in\mathbb{R}.$

Additionally, we assume to have a fixed shape $\Omega$ and a fixed point $z\in
D\setminus\partial\Omega$. The topological derivative then measures the
sensitivity of the shape functional $S$ with respect to infinitesimal
topological changes of the shape $\Omega$ around the point $z$. We
characterize the perturbed shape by

$\Omega_{z,\epsilon}=\begin{cases}\Omega\setminus\bar{\omega}_{z,\epsilon},&z\in\Omega\\\
\Omega\cup\omega_{z,\epsilon},&z\in D\setminus\bar{\Omega},\end{cases}$

where $\omega_{z,\epsilon}=z+\epsilon\omega$ with
$\omega\subset\mathbb{R}^{d}$ and $0\in\omega$ represents the shape of the
perturbation. Furthermore, we assume to have a positive function $l$ with
$\lim_{\epsilon\rightarrow 0}l(\epsilon)=0$. The topological derivative, see
for example [2], is then defined as

$D_{T}S(\Omega,\omega)(z)=\lim_{\epsilon\rightarrow
0}\frac{S(\Omega_{z,\epsilon})-S(\Omega)}{l(\epsilon)},$ (8)

if the limit exists. Typical choices are $w=B_{1}(0)$ with
$l(\epsilon)=|\omega_{z,\epsilon}|$, where $B_{1}(0)$ is the open unit ball
with center $0$.

###  3.2 Topological derivative for our model problem

We formally present the topological derivative of our model problem
$(\ref{optproblem})$, which can be derived, for example, with an averaged
adjoint approach [8]. The topological derivative for $(\ref{objective})$ is
given by

$D_{T}J(z)=-(\alpha_{U}-\alpha_{L})u(z)v(z)$ (9)

for all $z\in D\setminus\partial\Omega$. Here, $u$ is the weak solution of
$(\ref{stokesdarcy})$ and $v$ is the adjoint velocity which solves

$\begin{split}-\Delta v+\alpha v+\nabla q-\frac{1}{\Delta t}v_{s}&=0&&\text{in
}D,\\\ \text{div}(q)&=0&&\text{in }D,\\\ v&=0&&\text{on
}\Gamma_{\text{in}},\\\ v&=0&&\text{on }\Gamma_{\text{wall}},\\\
\partial_{n}v-qn&=0&&\text{on }\Gamma_{\text{out}},\end{split}$ (10)

and $v_{s}$ is the adjoint of the smoothed velocity which solves the equation

$\begin{split}\frac{1}{\Delta t}v_{s}-\Delta v_{s}&=2\frac{u_{s}}{\left\lVert
u_{s}\right\rVert}\min(0,\left\lVert u_{s}\right\rVert-u_{t})&&\text{in }D,\\\
\nabla v_{s}\cdot n&=0&&\text{on }\partial D.\end{split}$ (11)

The actual computation of the topological derivative is beyond the scope of
this paper and a topic of future research. For more details regarding
topological sensitivity analysis, we refer the reader to [2].

###  3.3 A gradient-based solution algorithm for topology optimization

We present an algorithm using the topological derivative as a search direction
in order to solve topology optimization problems. The method was introduced by
Amstutz and Andrä [4]. For the sake of brevity, we only give a short
presentation here and refer the reader to [8] for an overview over established
topology optimization algorithms as well as novel quasi-Newton methods for
topology optimization. First, we introduce a continuous level-set function
$\psi:D\rightarrow\mathbb{R}$ to characterize the fluid as well as the solid
part of the hold-all domain $D$ by

$\psi(z)\begin{cases}<0,&z\in\Omega,\\\ =0,&z\in\partial\Omega,\\\ >0.&z\in
D\setminus\bar{\Omega}.\end{cases}$

Additionally, we define the generalized topological derivative

$g(z)=\begin{cases}-D_{T}S(\Omega,\omega)(z),&z\in\Omega,\\\
+D_{T}S(\Omega,\omega)(z),&z\in D\setminus\bar{\Omega}.\end{cases}$

The algorithm then updates the level-set function consequently by a linear
combination of itself and the generalized topological derivative on the
$L^{2}$-sphere. The update formula reads

$\psi_{n+1}=\frac{1}{\sin(\theta_{n})}\left[\sin((1-\theta_{n})\kappa_{n})\psi_{n}+\sin(\theta_{n}\kappa_{n})\frac{g_{n}}{\left\lVert
g_{n}\right\rVert_{L^{2}(D)}}\right].$ (12)

Here, $\theta_{n}$ is the $L^{2}$-angle between the level-set function and the
generalized topological derivative

$\theta_{n}=\arccos\left[\frac{\langle g,\psi\rangle_{L^{2}(D)}}{\left\lVert
g\right\rVert_{L^{2}(D)}\left\lVert\psi\right\rVert_{L^{2}(D)}}\right].$

Additionally, $\kappa_{n}$ is chosen by a line search approach to guarantee
that the objective function value decreases each iteration. This procedure is
carried out until a local optimality condition is reached. Such an optimality
condition has been derived e.g. in [9]. Numerically, the algorithm is stopped,
when the angle $\theta_{n}$ becomes smaller than a certain numerical tolerance
$\epsilon_{\theta}$. This then implies a local optimality condition, again the
reader is referred to [4] or [9]. We want to conclude by stating that all
$\psi_{n}$ fulfill $\left\lVert\psi_{n}\right\rVert_{L^{2}(D)}=1$ if the
condition holds for the initial iteration.

##  4 Numerical results

We turn back to optimization problem $(\ref{optproblem})$ and describe our
numerical setting. First, the hold all domain is given by
$D=(0,1)\times(0,1)$. We define a quadratic inflow profile $u_{\text{in}}$
given by

$u_{\text{in}}=\begin{bmatrix}-\frac{400}{9}(y-0.35)(y-0.65)\\\
0\end{bmatrix}\;\text{for }x=0\text{ and }0.35\leq y\leq 0.65$

on the inflow boundary $\Gamma_{\text{in}}$. The situation is displayed in
Figure 4. The inverse permeability $\alpha$ for the fluid as well as the solid
part as is defined as

$\alpha_{L}=\frac{2.5}{125^{2}},\;\;\alpha_{U}=\frac{2.5}{0.008^{2}}.$

Furthermore, the stopping criterion for the algorithm is given by
$\epsilon_{\theta}=0.035$ which is equivalent to an angle of $2$ degrees. The
target velocity is chosen as $u_{t}=0.1$ and the lower and upper values for
the volume constraint are given by $V_{L}=0.5$ and $V_{U}=0.7$, respectively.

$D$$l=1.0$$l=1.0$$0.35$$0.5$

Schematic setup for the bipolar plate.

###  4.1 Numerical implementation

We give a short overview of the software used for the numerical
implementation. First, all underlying PDEs are solved using the finite element
package FEniCS [10]. We use a uniform triangular mesh consisting of 10.000
elements. The Stokes Darcy equation $(\ref{stokesdarcy})$ is discretized using
LBB-stable Taylor-Hood finite elements which consist of quadratic Lagrange
elements for the velocity and linear Lagrange elements for the pressure.
Additionally, for the smoothing equation $(\ref{heat1step})$, continuous
quadratic Lagrange elements are used.

Furthermore, we use the software package cashocs [11] for adjoint
computations. Cashocs is an open source software package with a python
interface which is based on FEniCS. The software can be used to solve
optimization problems constrained by partial differential equations in the
context of shape optimization and optimal control in an automated fashion.
Furthermore, cashocs implements a discretization of the continuous adjoint
approach which allows for adjoints to be computed automatically.

As we also introduced a volume constraint for $(\ref{optproblem})$, we give an
explanation on how the constraint is handled numerically. After each update of
the level-set function we perform a projection of the level-set function onto
the set of the admissible shapes. First, we assume to have a newly computed
update of the level-set function $\psi_{n+1}$, which was generated with
$(\ref{update})$, and the corresponding shape $\Omega_{n+1}$. Additionally, we
assume that the volume constraint is not fulfilled, otherwise nothing needs to
be done. Without loss of generality we assume that $|\Omega_{n+1}|>V_{U}$. To
perform the projection, we move the level-set function upwards by a strictly
positive constant until the upper border $V_{U}$ of the volume constraint is
reached. This means we search for a $c>0$, to move the level-set function
$\tilde{\psi}_{n+1}=\psi_{n+1}+c$, such that the corresponding shape
$\tilde{\Omega}_{n+1}$ fulfills $|\tilde{\Omega}_{n+1}|=V_{U}$. The case
$|\Omega_{n+1}|<V_{L}$ is handled analogously by taking a strictly negative
constant to alter the level-set function. Numerically, we compute the constant
$c$ by a bisection approach with numerical accuracy $\epsilon_{c}=10^{-4}$. As
we do not need to solve any additional PDEs, the computational cost for the
projection procedure is reasonable.

###  4.2 Numerical results

a)  b)  c)

Figure 1: Result of the topology optimization problem $(\ref{optproblem})$
with $\Delta t=10^{-3}$. In Figure a the final shape is displayed, in Figure b
the corresponding norm of the velocity field and in Figure c the norm of the
smoothed velocity field.

We start by choosing $\Delta t=10^{-3}$ as the step length for the smoothing
equation $(\ref{heat1step})$. The resulting optimized geometry and flow
velocity are displayed in Figure 1. The algorithm took 137 iterations to reach
the target accuracy. The objective function value is given by $8.6\cdot
10^{-8}$. From the obtained velocity fields and the previous investigations,
we observe that the optimization algorithm indeed produced a geometry with a
very uniform flow distribution. To do so, six fluid channels have been formed
and the flow is distributed rather uniformly between them. In more detail, the
target velocity goal $(\ref{optgoal2})$ for the smoothed velocity $u_{s}$ is
reached on $93.84\%$ of the whole domain. Furthermore, this extends to a
$91.79\%$ fulfillment of $(\ref{optgoal1})$ for $u$ on the fluid area.
Additionally, all obstacle sizes in fluid flow direction are about the same
order of magnitude.

We want to investigate the obstacle sizes for different values of the
smoothing parameter $\Delta t$. For that we choose two additional values for
the step length given by $\Delta t_{1}=5\cdot 10^{-3}$ as well as $\Delta
t_{2}=5\cdot 10^{-4}$. The results can be seen in Figure 2 for $\Delta t_{1}$
and in Figure 3 for $\Delta t_{2}$. It is apparent that the solid structures
for the higher time step length get larger which can be explained by the
higher influence of the smoothing. The larger the smoothing parameter is
chosen, the larger a solid inclusion has to be in order to be resolved in the
smoothed velocity $u_{s}$. In particular, we now observe that only 3 major
fluid channels can be found in the optimized geometry. Here, the target
velocity goal $(\ref{optgoal2})$ is fulfilled on $98.26\%$ of the domain $D$.
For the smaller time step length $\Delta t_{2}$ on the contrary, the solid
structures get smaller, which is again a consequence of the smaller influence
of the smoothing and we observe that the geometry branches out into several
channels which are then merged back together towards the outlet. We achieve a
$89.46\%$ realization of the target velocity goal for the smoothed velocity
$u_{s}$.

##  5 Conclusion

Summarizing, we presented a model for reaching a uniform flow distribution in
the anode side bipolar plate for hydrogen electrolysis cells. We showed in
numerical tests that a uniform flow distribution is indeed reached in the
optimized shapes. Additionally, we investigated the influence of the time step
length $\Delta t$ on the size of the solid structures, giving us a tool to
control the object sizes in the resulting shapes.

a)  b)  c)

Figure 2: Result of the topology optimization problem $(\ref{optproblem})$
with $\Delta t_{1}=5\cdot 10^{-3}$. In Figure a the final shape is displayed,
in Figure b the corresponding norm of the velocity field and in Figure c the
norm of the smoothed velocity field.

a)  b)  c)

Figure 3: Result of the topology optimization problem $(\ref{optproblem})$
with $\Delta t_{2}=5\cdot 10^{-4}$. In Figure a the final shape is displayed,
in Figure b the corresponding norm of the velocity field and in Figure c the
norm of the smoothed velocity field.

L. Baeck gratefully acknowledges financial support from the Fraunhofer
Institute of Industrial Mathematics ITWM.

## References

  * [1] S. Metz, T. Smolinka, C. I. Bernäcker, S. Loos, T. Rauscher, L. Röntzsch, M. Arnold, A. L. Görne, M. Jahn, M. Kusnezoff, G. Kolb, U. P. Apfel, and C. Doetsch, in: Producing hydrogen through electrolysis and other processes, edited by R. Neugebauer (Springer International Publishing, Cham, 2023), pp. 203–252.
  * [2] J. Sokolowski and A. Zochowski, SIAM Journal on Control and Optimization 37(4), 1251–1272 (1999).
  * [3] H. A. Eschenauer, V. V. Kobelev, and A. Schumacher, Structural optimization 8, 42–51 (1994).
  * [4] S. Amstutz and H. Andrä, J. Comput. Phys 216(2), 573–588 (2006).
  * [5] T. Borrval and J. Petersson, Internat. J. Numer. Methods Fluids 41(1), 77–107 (2003).
  * [6] L. N Sá, R. R Amigo, A. Novot, and E. N Silva, Structural and Multidisciplinary Optimization 54(2), 249–264 (2016).
  * [7] M. Hinze, R. Pinnau, M. Ulbrich, and S. Ulbrich, Optimization with PDE constraints, Mathematical Modelling: Theory and Applications, Vol. 23 (Springer, New York, 2009).
  * [8] S. Blauth and K. Sturm, Quasi-Newton Methods for Topology Optimization Using a Level-Set Method, 2023\.
  * [9] S. Amstutz, Optim. Methods Softw 26(4-5), 555–573 (2011).
  * [10] M. S. Alnæs, J. Blechta, J. Hake, A. Johansson, B. Kehlet, A. Logg, C. Richardson, J. Ring, M. E. Rognes, and G. N. Wells, Archive of Numerical Software 3(100) (2015).
  * [11] S. Blauth, Version 2.0 – cashocs: A computational, adjoint-based shape optimization and optimal control software, 2023.


Qingan<EMAIL_ADDRESS>Pan<EMAIL_ADDRESS>Nitin
<EMAIL_ADDRESS>Yuxin<EMAIL_ADDRESS>Yuan
<EMAIL_ADDRESS>Yi<EMAIL_ADDRESS>OPPO US Research
Center
InnoPeak Technology, Inc.
Palo Alto, California,US Wing LLC
Mountain View, California, US
(†Work done while at OPPO US Research Center.)  ∗Corresponding author
FisheyeDistill

# FisheyeDistill: Self-Supervised Monocular Depth Estimation with Ordinal
Distillation for Fisheye Cameras

###### Abstract

In this paper, we deal with the problem of monocular depth estimation for
fisheye cameras in a self-supervised manner. A known issue of self-supervised
depth estimation is that it suffers in low-light/over-exposure conditions and
in large homogeneous regions. To tackle this issue, we propose a novel ordinal
distillation loss that distills the ordinal information from a large teacher
model. Such a teacher model, since having been trained on a large amount of
diverse data, can capture the depth ordering information well, but lacks in
preserving accurate scene geometry. Combined with self-supervised losses, we
show that our model can not only generate reasonable depth maps in challenging
environments but also better recover the scene geometry. We further leverage
the fisheye cameras of an AR-Glasses device to collect an indoor dataset to
facilitate evaluation.

##  1 Introduction

Fisheye cameras have gradually gained popularity on Head Mounted Display
(HMD), as it offers a large field-of-view (FoV) that is well suited for
immersive AR/VR experience. To estimate the depth map from a distorted fisheye
image, a typical prerequisite is to perform image rectification. Such a
rectification step, however, undermines the benefits of using fisheye cameras
because it will inevitably reduce the FoV of the camera. Moreover, image
rectification itself is expensive on HMD as there are only very limited
computational resources available on the device. Therefore, it’s desirable to
directly estimate a depth map from the fisheye image, which is the main theme
of this paper.

The past few years have seen a growing interests in self-supervised depth
estimation [Garg et al.(2016)Garg, Bg, Carneiro, and Reid, Zhou et
al.(2017)Zhou, Brown, Snavely, and Lowe, Godard et al.(2019)Godard, Mac Aodha,
Firman, and Brostow] because it requires no groundtruth depth during training
and yet achieves comparable results with some of the well known supervised
methods. However, due to dearth of effective benchmarking datasets, depth
estimation for fisheye cameras remains under-explored, with only a few
exceptions [Kumar et al.(2020b)Kumar, Yogamani, Bach, Witt, Milz, and Mader,
Kumar et al.(2020a)Kumar, Hiremath, Bach, Milz, Witt, Pinard, Yogamani, and
Mäder]. Nonetheless, those fisheye depth estimation methods [Kumar et
al.(2020b)Kumar, Yogamani, Bach, Witt, Milz, and Mader, Kumar et
al.(2020a)Kumar, Hiremath, Bach, Milz, Witt, Pinard, Yogamani, and Mäder] are
only evaluated on well-conditioned outdoor self-driving datasets, including
the WoodScape dataset [Yogamani et al.(2019)Yogamani, Hughes, Horgan, Sistu,
Varley, O’Dea, Uricár, Milz, Simon, Amende, et al.], for which the groundtruth
depths are not yet publicly available. In this work, we follow this trend to
take a self-supervised approach to estimating a depth map from a distorted
image, but under a more challenging environment, which is captured by HMD.

Self-supervised depth estimation is not without challenges because the key
supervision signal of photo-consistency can become ineffective in low-
light/over-exposure conditions and in large areas of no/little textures as
commonly seen in indoor environments. On the other hand, there has been depth
models that are trained on a variety of datasets and can generalize well to
new unseen domains [Ranftl et al.(2019)Ranftl, Lasinger, Hafner, Schindler,
and Koltun, Ranftl et al.(2021)Ranftl, Bochkovskiy, and Koltun, Yin et
al.(2021)Yin, Zhang, Wang, Niklaus, Mai, Chen, and Shen]. Although such depth
models can generate visually good-looking depth maps, the 3D scene geometry in
those depth maps are often not well recovered [Yin et al.(2021)Yin, Zhang,
Wang, Niklaus, Mai, Chen, and Shen]. For example, when visualizing the depth
maps of an indoor room in 3D, the horizontal and vertical walls may not be
perpendicular to each other, showing distorted 3D scene geometry [Yin et
al.(2021)Yin, Zhang, Wang, Niklaus, Mai, Chen, and Shen].

In this work, we propose to combine the best of self-supervised depth
estimation and knowledge distillation [Hinton et al.(2015)Hinton, Vinyals, and
Dean] to build a robust depth model for fisheye cameras especially in
challenging indoor environments. In particular, we adapt the photometric loss
[Godard et al.(2019)Godard, Mac Aodha, Firman, and Brostow], which is usually
based on a pin-hole camera model, to a fisheye camera model. On top of the
self-supervised losses, we further devise a novel ordinal distillation loss
which aims to transfer the depth ordering information from a teacher model
into our target model. The intuition of ordinal distillation is that the
teacher model [Ranftl et al.(2019)Ranftl, Lasinger, Hafner, Schindler, and
Koltun, Ranftl et al.(2021)Ranftl, Bochkovskiy, and Koltun] is usually good at
predicting the depth ordering relationships of pixels (_e.g._ , a certain
pixel is closer or farther than other pixels), albeit not reflecting the exact
3D geometry. Our final depth model becomes more robust and accurate by (i)
learning the ordering relationships of all pixels via ordinal distillation and
(ii) respecting the 3D geometry enforced by the photometric loss.

To perform quantitative evaluation, we collect a new dataset which contains
stereo image sequences captured by a pair of stereo fisheye cameras mounted on
an AR-Glasses device, and generate the pseudo groundtruth depth maps using
stereo matching [Hirschmuller(2007)]. On this dataset, we show that our
proposed model leads to significant improvements over the baseline models.

##  2 Related Work

In this section, we briefly review a few related works on self-supervised and
supervised monocular depth estimation.

### Self-Supervised Depth Estimation

Most of self-supervised depth estimation methods assume that the camera is
calibrated and the images are rectified such that a photometric loss can be
constructed via a combination of unprojection and projection operations. Garg
_et al._ [Garg et al.(2016)Garg, Bg, Carneiro, and Reid] are the first to
leverage the photo-consistency between stereo images to build a self-
supervised loss for training a deep depth model. Zhou  _et al._ [Zhou et
al.(2017)Zhou, Brown, Snavely, and Lowe] extend to using temporal sequences to
train a depth model and a pose model with a temporal photometric loss. A lot
of follow-up methods are then proposed to improve [Garg et al.(2016)Garg, Bg,
Carneiro, and Reid, Zhou et al.(2017)Zhou, Brown, Snavely, and Lowe] by new
losses. Depth consistency losses are introduced to enforce the network to
predict consistent depth maps in stereo [Godard et al.(2017)Godard, Mac Aodha,
and Brostow] or temporal [Bian et al.(2019)Bian, Li, Wang, Zhan, Shen, Cheng,
and Reid, Luo et al.(2020)Luo, Huang, Szeliski, Matzen, and Kopf] images. Wang
_et al._ [Wang et al.(2018)Wang, Buenaposada, Zhu, and Lucey] observe a scale
diminishing issue in monocular training and propose a depth normalization
method to counter this issue. A few methods [Wang et al.(2019)Wang, Pizer, and
Frahm, Zou et al.(2020)Zou, Ji, Tran, Huang, and Chandraker] employ recurrent
neural networks to model long-term temporal dependencies for self-supervised
training. Some other methods [Chen et al.(2019)Chen, Schmid, and Sminchisescu,
Yin and Shi(2018), Zou et al.(2018)Zou, Luo, and Huang] introduce an
additional optical flow network to build cross-task consistencies in a
monocular training setup. Godard  _et al._ [Godard et al.(2019)Godard, Mac
Aodha, Firman, and Brostow] comprehensively analyze the challenges faced in
self-supervised depth learning, including occlusion, static pixels, and
texture-copying artifacts, and propose a set of novel techniques to handle
those challenges. Tiwari  _et al._ [Tiwari et al.(2020)Tiwari, Ji, Tran,
Zhuang, Anand, and Chandraker] propose to combine self-supervised depth and
geometric SLAM in a self-improving loop. Ji  _et al._ [Ji et al.(2021)Ji, Li,
Bhanu, and Xu] observe the challenges for indoor self-supervised depth
estimation and come up with a depth factorization module and a residual pose
estimation module to improve the performance under indoor environments. Watson
_et al._ [Watson et al.(2021)Watson, Mac Aodha, Prisacariu, Brostow, and
Firman] use multiple temporal images as input to the depth model and put forth
a teacher-student framework to deal with moving objects in the scene. Ji  _et
al._ [Ji et al.(2022)Ji, Yan, Ma, and Xu] propose an online system for
refining a depth model in a self-supervised manner, combining a robustified
SLAM system and a monocular depth model.

Figure 1: The system workflow of our FisheyeDistill. Our system consists of
two main parts, _i.e._ , a student model trained with self-supervised losses
and distillation losses, and a teacher model (MiDaS) [Ranftl et
al.(2019)Ranftl, Lasinger, Hafner, Schindler, and Koltun, Ranftl et
al.(2021)Ranftl, Bochkovskiy, and Koltun] that predicts depths to guide the
training of the student model. The SID stands for the scale-invariant
distillation and $SID_{M}$ is the corresponding function module described in
Sec. 3.3. ORD and $ORD_{M}$ are the ordinal distillation loss and the module.

### Supervised Depth Estimation

Ever since EigenDepth [Eigen et al.(2014)Eigen, Puhrsch, and Fergus], many
methods have been proposed to improve supervised monocular depth estimation on
a specific dataset, either by using better loss terms [Fu et al.(2018)Fu,
Gong, Wang, Batmanghelich, and Tao, Liu et al.(2015)Liu, Shen, Lin, and Reid]
or via multi-task learning [Ummenhofer et al.(2017)Ummenhofer, Zhou, Uhrig,
Mayer, Ilg, Dosovitskiy, and Brox, Teed and Deng(2018), Zhou et al.(2018)Zhou,
Ummenhofer, and Brox, Qi et al.(2018)Qi, Liao, Liu, Urtasun, and Jia, Guo et
al.(2018)Guo, Li, Yi, Ren, and Wang, Qiao et al.(2021)Qiao, Zhu, Adam, Yuille,
and Chen]. Instead of focusing on the depth prediction of one single dataset
(or domain), some other methods seek to train a depth model that can learn
across different domains. Along this line, Megadepth [Li and Snavely(2018)]
exploits internet photos to train a depth model with structure-from-motion
reconstructions [Schonberger and Frahm(2016)]. Li  _et al._ [Li et
al.(2019)Li, Dekel, Cole, Tucker, Snavely, Liu, and Freeman] leverage internet
videos with frozen people to predict sharp depth maps for images with people.
MiDaS [Ranftl et al.(2019)Ranftl, Lasinger, Hafner, Schindler, and Koltun]
proposes a scale-invariant loss to train the depth model on a large and
diverse training sets, which facilitates generality and cross-dataset learning
by providing a dataset-agnostic depth model. DPT [Ranftl et al.(2021)Ranftl,
Bochkovskiy, and Koltun] further improves MiDaS by using a vision transformer
network [Dosovitskiy et al.(2020)Dosovitskiy, Beyer, Kolesnikov, Weissenborn,
Zhai, Unterthiner, Dehghani, Minderer, Heigold, Gelly, et al.]. Although
visually pleasing depth maps can be generated by those models [Ranftl et
al.(2019)Ranftl, Lasinger, Hafner, Schindler, and Koltun, Ranftl et
al.(2021)Ranftl, Bochkovskiy, and Koltun], their quantitative performance on a
specific dataset is usually worse than those that are trained specifically on
the same dataset. This means that the accurate 3D geometry information is not
well preserved in those models [Ranftl et al.(2019)Ranftl, Lasinger, Hafner,
Schindler, and Koltun, Ranftl et al.(2021)Ranftl, Bochkovskiy, and Koltun].

Our method tries to reap the benefits of both self-supervised and supervised
models in the sense that it respects the 3D geometry via the use of self-
supervised losses and enforces ordering consistency between neighboring depth
pixels by distilling information from a diversely trained teacher model.

##  3 Method

In this work, we aim at learning a self-supervised monocular depth estimation
model for fisheye images (see Fig. 1). In all the previous works [Garg et
al.(2016)Garg, Bg, Carneiro, and Reid, Zhou et al.(2017)Zhou, Brown, Snavely,
and Lowe, Godard et al.(2019)Godard, Mac Aodha, Firman, and Brostow], view
synthesis is usually used in self-supervision by learning a depth and pose
relationship to synthesize virtual target images from neighboring views
$I_{t-1}$ and $I_{t+1}$. For this, it requires a projection function $\Phi$
that maps 3D points $P_{i}$ in 3D space to image coordinates
$p_{i}=\Phi(P_{i})$, and accordingly the corresponding unprojection function
$\Phi^{-1}$, which converts image pixels, based on the estimated depth map
$D$, into 3D space $P_{i}=\Phi^{-1}(p_{i},D)$ in order to acquire color
information from other views.

###  3.1 Fisheye Geometry Model

For fisheye cameras, given a 3D point $P_{i}=(X_{i},Y_{i},Z_{i})^{T}$ in
camera coordinates and with $(x_{i}=X_{i}/Z_{i},y_{i}=Y_{i}/Z_{i})$, the
projection function $p_{i}=\Phi(P_{i})$ from 3D point $P_{i}$ to distorted
image pixel $p_{i}=(u,v)^{T}$ can be obtained through the following mapping
equation:

$\Phi(P_{i})=\begin{bmatrix}u=f_{x}\cdot x_{d}+c_{x}\\\ v=f_{y}\cdot
y_{d}+c_{y}\end{bmatrix},$ (1)

where
$(x_{d}=\dfrac{x_{i}}{r}\phi_{\theta},y_{d}=\dfrac{y_{i}}{r}\phi_{\theta})$
with $r=\sqrt{{x_{i}}^{2}+{y_{i}}^{2}}$, $\theta=\arctan(r)$ is the angle of
incidence,
$\phi_{\theta}=\theta\cdot(1+k_{1}\cdot\theta^{2}+k_{2}\cdot\theta^{4}+k_{3}\cdot\theta^{6}+k_{4}\cdot\theta^{8})$
is the polynomial radial distortion model mapping the incident angle to the
image radius, $(f_{x},f_{y})$ and $(c_{x},c_{y})$ stand for the focal length
and the principal point derived from the intrinsic matrix $K$, and
$\\{k_{i}\\}$ denote the set of fisheye distortion coefficients.

For the unprojection part, due to the existence of fisheye distortion, we are
unable to directly transform a pixel $p_{i}$ into camera coordinates via the
pinhole model. In order to achieve that, we first generate an intermediate
rectified depth map $\hat{D}$ from the depth map $D$ estimated from our
network. This can be done by warping a pixel grid according to Eq. (1). Then
we leverage the rectified $\hat{D}$ to unproject the grid into 3D by simply
applying $\Phi^{-1}(p_{i},\hat{D})=\hat{D}K^{-1}\cdot p_{i}$. 111Here and
below, we omit a necessary step of converting to the homogeneous coordinates
for notation simplicity. The view synthesis process in our fisheye pipeline
can be summarized as: (i) Unproject a uniform pixel grid which has the same
resolution with input frames through $\Phi^{-1}(p_{i},\hat{D})$; (ii) Project
those 3D points by $\Phi(P_{i})$ and the associated pose information from pose
network module, to obtain distorted synthesis images. This process is
different from the FisheyeDistanceNet [Kumar et al.(2020a)Kumar, Hiremath,
Bach, Milz, Witt, Pinard, Yogamani, and Mäder] which needs to numerically
calculate $\theta$ out from the $4^{th}$ order polynomial.

###  3.2 Self-Supervised Losses

Following [Godard et al.(2019)Godard, Mac Aodha, Firman, and Brostow], we use
a photometric reprojection loss between two neighboring fisheye images and an
edge-aware depth smoothness loss as the self-supervised losses to train the
models. Specifically, given a target fisheye image $I_{t}$ and a source
fisheye image $I_{t}^{\prime}$, two models (_i.e._ , a depth model and a pose
model) are jointly trained to predict a dense depth map $D_{t}$ and a relative
transformation matrix $T_{tt^{\prime}}$. The per-pixel minimum photometric
reprojection loss [Godard et al.(2019)Godard, Mac Aodha, Firman, and Brostow]
can be computed as

$L^{ph}=\min\limits_{t^{\prime}}\rho(I_{t},I_{tt^{\prime}}),$ (2)

and

$I_{tt^{\prime}}=I_{t^{\prime}}\big{\langle}\Phi\big{(}T_{tt^{\prime}}\Phi^{-1}(p_{i},\hat{D})\big{)}\big{\rangle},$
(3)

where $t^{\prime}\in\\{t-1,t+1\\}$, $\rho(\cdot)$ denotes a weighted
combination of the L1 and Structured SIMilarity (SSIM) losses [Godard et
al.(2019)Godard, Mac Aodha, Firman, and Brostow], and
$\big{\langle}\cdot\big{\rangle}$ is the bilinear sampling operator. The edge-
aware smoothness loss is defined as

$L^{sm}=|\partial_{x}d^{\ast}_{t}|e^{-|\partial_{x}I_{t}|}+|\partial_{y}d^{\ast}_{t}|e^{-|\partial_{y}I_{t}|},$
(4)

where $d^{\ast}_{t}=d/\bar{d}_{t}$ is the mean-normalized inverse depth from
[Wang et al.(2018)Wang, Buenaposada, Zhu, and Lucey]. An auto-masking
mechanism [Godard et al.(2019)Godard, Mac Aodha, Firman, and Brostow] is also
adopted to mask out static pixels when computing the photometric loss.

###  3.3 Distillation Losses

Self-supervised losses are volatile when the brightness or RGB values of a
pixel become indistinguishable with its neighboring pixels. Previously, such
color ambiguity is implicitly handled via the smoothness loss and a multi-
scale strategy [Godard et al.(2019)Godard, Mac Aodha, Firman, and Brostow],
but they do not work well if the textureless regions are large. For example,
in indoor environments, textureless walls often occupy a large area in the
image; or in low-light/over-exposure conditions, most of the regions in an
image are near-textureless. To provide supervision in those scenarios, we
apply knowledge distillation [Hinton et al.(2015)Hinton, Vinyals, and Dean] to
distill depth information from a diversely trained teacher model [Ranftl et
al.(2019)Ranftl, Lasinger, Hafner, Schindler, and Koltun, Ranftl et
al.(2021)Ranftl, Bochkovskiy, and Koltun] into a smaller target model. In
particular, we propose to combine a novel depth ordinal distillation loss and
a scale-invariant distillation loss [Ranftl et al.(2019)Ranftl, Lasinger,
Hafner, Schindler, and Koltun], which we detail as below.

Pretrained models such as MiDaS [Ranftl et al.(2019)Ranftl, Lasinger, Hafner,
Schindler, and Koltun, Ranftl et al.(2021)Ranftl, Bochkovskiy, and Koltun] are
good at predicting relative depth relationships as they are trained on large
and diverse datasets. In light of this, we aim to distill from the teacher
model the ordinal information between neighboring pixels into our student
model. Inspired by [Chen et al.(2016)Chen, Fu, Yang, and Deng], we propose to
use a ranking loss to build an ordinal distillation loss between the depth
maps from the teacher model and the target student model.

Given a depth map $D^{MONO}$ from the student model and a depth map
$D^{MIDAS}$ from the teacher model, we compute a ranking loss from each pixel
$I(i,j)$ to its left neighbor $I(i-1,j)$ as follows,

$L_{l}^{od}=\begin{cases}\log\big{(}1+\exp(-D_{i,j}^{MONO}+D_{i-1,j}^{MONO})\big{)},&\text{if
$D_{i,j}^{MIDAS}>\alpha D_{i-1,j}^{MIDAS}$}\\\
\log\big{(}1+\exp(D_{i,j}^{MONO}-D_{i-1,j}^{MONO})\big{)},&\text{if
$D_{i,j}^{MIDAS}<\beta D_{i-1,j}^{MIDAS}$}\\\
|D_{i,j}^{MONO}-D_{i-1,j}^{MONO}|,&\text{otherwise,}\end{cases}$ (5)

where $\alpha$ and $\beta$ are hyper-parameters controlling the ranking gap
and are empirically set to $\alpha=1.1$ and $\beta=0.9$ in our experiments.
Similarly, we can compute a ranking loss from each pixel $I(i,j)$ to its
neighbor above it: $I(i,j-1)$ as $L_{t}^{od}$ . The final ordinal distillation
loss is then the sum of $L_{l}^{od}$ and $L_{t}^{od}$, _i.e._ ,

$L^{od}=L_{l}^{od}+L_{t}^{od}.$ (6)

Intuitively, if a pixel is predicted to be farther (or closer) than its
neighbors by the teacher model, by minimizing the ordinal distillation loss,
the student model is encouraged to predict similar depth ordering
relationships. If in the teacher model two neighboring pixels are predicted to
be close, a smoothness term is enforced for the student model.

We further add a scale-invariant distillation loss $L^{sd}$ between $D^{MONO}$
and $D^{MIDAS}$ to strengthen the supervision in textureless regions.
Following [Ranftl et al.(2019)Ranftl, Lasinger, Hafner, Schindler, and
Koltun], $D^{MONO}$ and $D^{MIDAS}$ are respectively normalized to be scale-
and shift-invariant (subtracted by the median scale and then divided by the
mean shift). The normalized versions are denoted as $\tilde{D}^{MONO}$ and
$\tilde{D}^{MIDAS}$. The scale-invariant loss is then defined as:

$L^{sd}=|\tilde{D}^{MONO}-\tilde{D}^{MIDAS}|.$ (7)

Our final training loss is a weighted combination of the photometric loss, the
smoothness loss and the distillation losses, _i.e._ ,

$L=L^{ph}+w^{sm}L^{sm}+w^{od}L^{od}+w^{sd}L^{sd}.$ (8)

##  4 AR-Glasses Fisheye Dataset

In order to evaluate the model, we contribute a large fisheye dataset for
training and evaluation.The dataset contains $321,300$ raw fisheye images
collected from a well calibrated AR-Glasses device in a stereo manner. All
images are acquired in an indoor environment with $8$ different scenes.
Besides various chairs, desks and decorations, there are also many lights,
texture-less walls and glasses therein. In addition, due to the computational
limit on the headset, all frames are captured under the resolution of
$640\times 400$ in gray-scale with a sub-optimal auto-exposure mechanism.
These make our dataset very challenging for the task of self-supervised depth
estimation as the photometric constraints would act poorly in such scenes. We
will have more discussion on it in Sec. 5.2.

More specifically, in our experiments, we use $238,186$ images for training
and $42,688$ and $20,213$ images as the validation and test sets respectively.
The $8$ scenarios are randomly shuffled and fed into the network while
training. Similar to [Godard et al.(2019)Godard, Mac Aodha, Firman, and
Brostow], we also apply the static frame filtering to clean our data and get a
final training and a validation set of $71,688$ and $6,265$ images. One major
limitation of other fisheye datasets [Matsuki et al.(2018)Matsuki, von
Stumberg, Usenko, Stückler, and Cremers] is that there is no easily accessible
ground truth for evaluation. Instead, we utilize the semi-global stereo
matching method [Hirschmuller(2007)] based on OpenCV’s implementation to
calculate pseudo ground truth for our dataset. In our experiments, we select
$300$ well calculated stereo depth maps from different viewpoints as ground
truth. We conduct all the experiments in the monocular mode.

##  5 Experiments

###  5.1 Implementation Details

We implement our system using Pytorch and employ Adam optimizer to minimize
the objective function in Eq. (8). We train the model for 20 epochs with a
learning rate of $10^{-4}$ on a 12GB GeForce RTX. Due to the limitation of
graphic memory, a small batch size of 6 is leveraged while training and batch
normalization [Ioffe and Szegedy(2015)] layers are also frozen to avoid
unexpected degeneration. We convert the sigmoid output $\sigma$ from the
network to depth via $D=1/(a\cdot\sigma+b)$, where $a$ and $b$ are chosen to
bound $D$ between $0.1$ and $100$ units. As the resolution of our AR-Glasses
images is $640\times 400$, we feed $640\times 384$ pixels as the network input
to approximate the original aspect ratio. The smoothness loss weight $w^{sm}$
is set to $10^{-3}$ and distillation loss weights $w^{od}=w^{sd}=1.0$. For
generating the rectified intermediate depth $\hat{D}$, we apply a sightly
small focal length than the one given in $K$. This would help the network to
better learn the content on distorted image boundaries. In all our
experiments, we use $0.8$ as the scale factor.

As previously discussed, while the MiDaS teacher model [Ranftl et
al.(2021)Ranftl, Bochkovskiy, and Koltun] provides a good guidance for depth
estimation even on challenging structures, it suffers from the limitation in
keeping true geometric relationship, due to the lack of cross-view photometric
constraints. Therefore, directly adopting a constant weight for the knowledge
distillation through all epochs would weaken the ability of learning better
geometry information. To circumvent this issue, we adopt a decay weighting
strategy for the distillation losses in training. After each $10k$ steps, we
decay the weights $w^{od}$ and $w^{sd}$ by applying a factor
$s=0.9^{2\lambda}$, where $\lambda=\dfrac{steps}{10000}$. Such a strategy
enables a major contribution from the distillation losses when the uncertainty
is high and facilitates the convergence of photometric constraints in the
later training stage. After five epochs, photometric losses begin to play a
more dominant role in model training, which ensures the acquisition of
geometric relationship between temporal frames. Similar to other self-
supervised works, we also set the length of the training snippets to 3 and
build the losses over 4 image scales.

###  5.2 Evaluation

We evaluate our model mainly on the proposed AR-Glasses fisheye dataset. We do
not leverage KITTI [Geiger et al.(2013)Geiger, Lenz, Stiller, and Urtasun] as
it is not a fisheye dataset. Woodscape [Yogamani et al.(2019)Yogamani, Hughes,
Horgan, Sistu, Varley, O’Dea, Uricár, Milz, Simon, Amende, et al.] dataset
could be a good candidate, but it is still not fully released. Currently only
around $8k$ images are available and without ground truth for evaluation.
While the image resolution of AR-Glasses dataset is smaller than KITTI
dataset, the ground truth is much denser. Moreover, since there is no explicit
sky region in the indoor environment, we do not perform any cropping for each
test image. We evaluate the models on all pixels whose depth values are not
zero in the computed groundtruth.

The quantitative results are reported in Table 1, which shows the performance
of our model on the AR-Glasses dataset. From the table, we can see that
without considering the fisheye distortion, Monodepth2 [Godard et
al.(2019)Godard, Mac Aodha, Firman, and Brostow] is unable to extract
plausible depth outputs from distorted images. MiDaS [Ranftl et
al.(2021)Ranftl, Bochkovskiy, and Koltun] shows a good generalization ability
on our dataset, despite being trained on rectified images. We are unable to
compare with the fisheye model [Kumar et al.(2020a)Kumar, Hiremath, Bach,
Milz, Witt, Pinard, Yogamani, and Mäder], as its source code is not publicly
available. As fisheye cameras are designed for near-field sensing and deployed
on AR-Glasses for indoor applications, the groundtruth depth values are capped
at $20m$ for evaluation. In addition, we also demonstrate how the decay factor
$s$ affects the performance in training. We report the comparison results by
employing $s=0.9^{2\lambda}$ and $\hat{s}=0.98^{2\lambda}$ respectively, where
$\hat{s}$ enforces more emphasis on the distillation part. It can be found
that letting the network learn more geometric information through the temporal
photometric loss helps to improve the results.

We also render the qualitative results to demonstrate our advantages in
challenging indoor environments, as shown in Fig. 4. Directly training upon
photometric constraints without knowledge distillation is prone to producing
halo and texture-copy artifacts. This is because that large texture-less
continuous surfaces, such as planar white walls and desks, and the low light
condition make photometric losses unable to well establish pixel
correspondences from the very beginning. In contrast, due to the good
generalization ability on different datasets, the MiDaS teacher model [Ranftl
et al.(2021)Ranftl, Bochkovskiy, and Koltun], provides a smooth and sharp
depth ordering hint to let the network approximate plausible depth even on
texture-less planes. However, the geometric constraint within MiDaS is
insufficient. This is notable from Fig. 3, where we visualize the 3D geometry
unprojected from corresponding depth maps of different methods. While the
depth map is visually-pleasant, MiDaS suffers from severe distortion in 3D
space. Our method overcomes the limitations existing in both [Godard et
al.(2019)Godard, Mac Aodha, Firman, and Brostow] and [Ranftl et
al.(2021)Ranftl, Bochkovskiy, and Koltun] and is able to give more consistent
results on the challenging indoor dataset.

Table 1: Quantitative results of different algorithms on AR-Glasses fisheye dataset. All the approaches are evaluated in the monocular mode and the estimated depths are scaled using median scaling. Method | Abs Rel | Sq Rel | RMSE | RMSE_log | $\delta<1.25$ | $\delta<1.25^{2}$ | $\delta<1.25^{3}$
---|---|---|---|---|---|---|---
Monodepth2 [Godard et al.(2019)Godard, Mac Aodha, Firman, and Brostow] | 0.348 | 0.108 | 0.198 | 0.398 | 0.449 | 0.749 | 0.895
MiDaS [Ranftl et al.(2021)Ranftl, Bochkovskiy, and Koltun] | 0.216 | 0.131 | 0.252 | 0.267 | 0.756 | 0.918 | 0.965
Ours with $s=0.98^{2\lambda}$ | 0.136 | 0.022 | 0.103 | 0.183 | 0.843 | 0.971 | 0.991
Ours with $s=0.9^{2\lambda}$ | 0.131 | 0.020 | 0.098 | 0.177 | 0.854 | 0.973 | 0.992

Figure 2: Qualitative results on the AR-Glasses dataset. From left to right, the images are respectively input fisheye images, MiDaS estimated disparity maps, disparity maps from our network without any distillation losses, and disparity maps with distillation losses. Figure 3: Visual comparison of the 3D point clouds generated from depth maps by MiDaS (middle) and our method (right). Compared to MiDaS, our results generate more geometrically correct point clouds. Table 2: Ablation study by using different variants on the proposed AR-Glasses dataset. Depths are all capped at $20m$. SID, ORD, DW respectively represent the scale-invariant distillation, ordinal distillation and decay weighting. Method | SID | ORD | DW | Abs Rel | Sq Rel | RMSE | RMSE_log | $\delta<1.25$ | $\delta<1.25^{2}$ | $\delta<1.25^{3}$
---|---|---|---|---|---|---|---|---|---|---
Ours | ✗ | ✗ | ✗ | 0.151 | 0.021 | 0.103 | 0.200 | 0.780 | 0.961 | 0.992
Ours | ✓ | ✓ | ✗ | 0.134 | 0.020 | 0.098 | 0.180 | 0.846 | 0.972 | 0.992
Ours | ✓ | ✗ | ✓ | 0.134 | 0.019 | 0.097 | 0.179 | 0.844 | 0.972 | 0.992
Ours | ✓ | ✓ | ✓ | 0.131 | 0.020 | 0.098 | 0.177 | 0.854 | 0.973 | 0.992

###  5.3 Ablation Study

We conduct an ablation study to illustrate the benefits of different
components used in our system and show the results in Table 2. (i) Remove both
distillation losses: The network is trained only with distortion-based self-
supervised photometric losses. As previously discussed, the proposed AR-
Glasses dataset is very challenging, due to under-exposure and textureless
structures. Only relying on photometric constraints is prone to producing halo
and texture-copy artifacts. In contrast, the distillation strategy offers a
huge quantitative improvement. Fig. 3 also provides a visual comparison for
the sake of illustration. (ii) Remove the ordinal distillation loss: Ordinal
distillation loss is good at the supervision in textureless regions. Removing
the ordinal distillation loss from the network would diminish the performance
gain over the baseline. (iii) Utilize constant distillation loss weights:
Apart from the comparison of using different $s$ as illustrated before, we
also launch an experiment with a constant distillation weighting strategy. For
the trade-off between keeping depth ordering and learning true geometry, we
set the loss weights to $w^{od}=w^{sd}=0.3$ without decaying here for testing.
It can be seen that using constant distillation loss weights decrease the
performance on our dataset, which verifies the effectiveness of our decay-
weighting strategy.

##  6 Conclusions

In this work, we have presented FisheyeDistill, a self-supervised monocular
depth estimation method for fisheye cameras. Our key insight is to combine
self-supervised fisheye photometric losses with a novel ordinal distillation
loss to train a robust depth model that can work well in challenging
environments for fisheye cameras.

## Appendix

We adopt the same network structures with [Godard et al.(2019)Godard, Mac
Aodha, Firman, and Brostow]. In all experiments we use a standard ResNet18 [He
et al.(2016)He, Zhang, Ren, and Sun] encoder for both depth and pose networks.
The pose encoder takes fisheye frames as input, as our network is trained to
learn the pose under the distortion circumstance. We also empirically found
that using corresponding rectified images as input to the pose network
decreases its performance. For the distillation part, we leverage the
pretrained DPT [Ranftl et al.(2021)Ranftl, Bochkovskiy, and Koltun] large
model as our teacher model. For further demonstration, we also upload a demo
video, which respectively shows a portion of our contributed AR-glass fisheye
dataset and the corresponding depth estimation results (with and without
knowledge distillation). It is evident that without the proposed distillation
strategies, the results suffer from significant halo artifacts, as is also
shown below.

Figure 4: Depth estimation results on the AR-Glasses dataset. From left to
right are respectively the fisheye frames in our dataset, and the color-coded
depth predictions from our network without and with using the proposed
distillation strategies.

## References

  * [Bian et al.(2019)Bian, Li, Wang, Zhan, Shen, Cheng, and Reid] Jia-Wang Bian, Zhichao Li, Naiyan Wang, Huangying Zhan, Chunhua Shen, Ming-Ming Cheng, and Ian Reid.  Unsupervised scale-consistent depth and ego-motion learning from monocular video.  _arXiv preprint arXiv:1908.10553_ , 2019.
  * [Chen et al.(2016)Chen, Fu, Yang, and Deng] Weifeng Chen, Zhao Fu, Dawei Yang, and Jia Deng.  Single-image depth perception in the wild.  _arXiv preprint arXiv:1604.03901_ , 2016.
  * [Chen et al.(2019)Chen, Schmid, and Sminchisescu] Yuhua Chen, Cordelia Schmid, and Cristian Sminchisescu.  Self-supervised learning with geometric constraints in monocular video: Connecting flow, depth, and camera.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 7063–7072, 2019.
  * [Dosovitskiy et al.(2020)Dosovitskiy, Beyer, Kolesnikov, Weissenborn, Zhai, Unterthiner, Dehghani, Minderer, Heigold, Gelly, et al.] Alexey Dosovitskiy, Lucas Beyer, Alexander Kolesnikov, Dirk Weissenborn, Xiaohua Zhai, Thomas Unterthiner, Mostafa Dehghani, Matthias Minderer, Georg Heigold, Sylvain Gelly, et al.  An image is worth 16x16 words: Transformers for image recognition at scale.  _arXiv preprint arXiv:2010.11929_ , 2020.
  * [Eigen et al.(2014)Eigen, Puhrsch, and Fergus] David Eigen, Christian Puhrsch, and Rob Fergus.  Depth map prediction from a single image using a multi-scale deep network.  _arXiv preprint arXiv:1406.2283_ , 2014.
  * [Fu et al.(2018)Fu, Gong, Wang, Batmanghelich, and Tao] Huan Fu, Mingming Gong, Chaohui Wang, Kayhan Batmanghelich, and Dacheng Tao.  Deep ordinal regression network for monocular depth estimation.  In _CVPR_ , pages 2002–2011, 2018.
  * [Garg et al.(2016)Garg, Bg, Carneiro, and Reid] Ravi Garg, Vijay Kumar Bg, Gustavo Carneiro, and Ian Reid.  Unsupervised cnn for single view depth estimation: Geometry to the rescue.  In _ECCV_ , pages 740–756, 2016.
  * [Geiger et al.(2013)Geiger, Lenz, Stiller, and Urtasun] Andreas Geiger, Philip Lenz, Christoph Stiller, and Raquel Urtasun.  Vision meets robotics: The kitti dataset.  _The International Journal of Robotics Research_ , 32(11):1231–1237, 2013.
  * [Godard et al.(2017)Godard, Mac Aodha, and Brostow] Clément Godard, Oisin Mac Aodha, and Gabriel J Brostow.  Unsupervised monocular depth estimation with left-right consistency.  In _CVPR_ , pages 270–279, 2017.
  * [Godard et al.(2019)Godard, Mac Aodha, Firman, and Brostow] Clément Godard, Oisin Mac Aodha, Michael Firman, and Gabriel J Brostow.  Digging into self-supervised monocular depth estimation.  In _ICCV_ , pages 3828–3838, 2019.
  * [Guo et al.(2018)Guo, Li, Yi, Ren, and Wang] Xiaoyang Guo, Hongsheng Li, Shuai Yi, Jimmy Ren, and Xiaogang Wang.  Learning monocular depth by distilling cross-domain stereo networks.  In _ECCV_ , pages 484–500, 2018.
  * [He et al.(2016)He, Zhang, Ren, and Sun] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 770–778, 2016.
  * [Hinton et al.(2015)Hinton, Vinyals, and Dean] Geoffrey Hinton, Oriol Vinyals, and Jeff Dean.  Distilling the knowledge in a neural network.  _arXiv preprint arXiv:1503.02531_ , 2015.
  * [Hirschmuller(2007)] Heiko Hirschmuller.  Stereo processing by semiglobal matching and mutual information.  _TPAMI_ , 30(2):328–341, 2007.
  * [Ioffe and Szegedy(2015)] Sergey Ioffe and Christian Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  In _ICML_ , pages 448–456. PMLR, 2015.
  * [Ji et al.(2021)Ji, Li, Bhanu, and Xu] Pan Ji, Runze Li, Bir Bhanu, and Yi Xu.  Monoindoor: Towards good practice of self-supervised monocular depth estimation for indoor environments.  In _ICCV_ , pages 12787–12796, 2021.
  * [Ji et al.(2022)Ji, Yan, Ma, and Xu] Pan Ji, Qingan Yan, Yuxin Ma, and Yi Xu.  Georefine: Self-supervised online depth refinement for accurate dense mapping.  _arXiv preprint arXiv:2205.01656_ , 2022.
  * [Kumar et al.(2020a)Kumar, Hiremath, Bach, Milz, Witt, Pinard, Yogamani, and Mäder] Varun Ravi Kumar, Sandesh Athni Hiremath, Markus Bach, Stefan Milz, Christian Witt, Clément Pinard, Senthil Yogamani, and Patrick Mäder.  Fisheyedistancenet: Self-supervised scale-aware distance estimation using monocular fisheye camera for autonomous driving.  In _ICRA_ , pages 574–581, 2020a.
  * [Kumar et al.(2020b)Kumar, Yogamani, Bach, Witt, Milz, and Mader] Varun Ravi Kumar, Senthil Yogamani, Markus Bach, Christian Witt, Stefan Milz, and Patrick Mader.  Unrectdepthnet: Self-supervised monocular depth estimation using a generic framework for handling common camera distortion models.  _arXiv preprint arXiv:2007.06676_ , 2020b.
  * [Li and Snavely(2018)] Zhengqi Li and Noah Snavely.  Megadepth: Learning single-view depth prediction from internet photos.  In _CVPR_ , pages 2041–2050, 2018.
  * [Li et al.(2019)Li, Dekel, Cole, Tucker, Snavely, Liu, and Freeman] Zhengqi Li, Tali Dekel, Forrester Cole, Richard Tucker, Noah Snavely, Ce Liu, and William T Freeman.  Learning the depths of moving people by watching frozen people.  In _CVPR_ , pages 4521–4530, 2019.
  * [Liu et al.(2015)Liu, Shen, Lin, and Reid] Fayao Liu, Chunhua Shen, Guosheng Lin, and Ian Reid.  Learning depth from single monocular images using deep convolutional neural fields.  _TPAMI_ , 38(10):2024–2039, 2015.
  * [Luo et al.(2020)Luo, Huang, Szeliski, Matzen, and Kopf] Xuan Luo, Jia-Bin Huang, Richard Szeliski, Kevin Matzen, and Johannes Kopf.  Consistent video depth estimation.  _TOG_ , 39(4):71–1, 2020.
  * [Matsuki et al.(2018)Matsuki, von Stumberg, Usenko, Stückler, and Cremers] Hidenobu Matsuki, Lukas von Stumberg, Vladyslav Usenko, Jörg Stückler, and Daniel Cremers.  Omnidirectional dso: Direct sparse odometry with fisheye cameras.  _IEEE Robotics and Automation Letters_ , 3(4):3693–3700, 2018.
  * [Qi et al.(2018)Qi, Liao, Liu, Urtasun, and Jia] Xiaojuan Qi, Renjie Liao, Zhengzhe Liu, Raquel Urtasun, and Jiaya Jia.  Geonet: Geometric neural network for joint depth and surface normal estimation.  In _CVPR_ , pages 283–291, 2018.
  * [Qiao et al.(2021)Qiao, Zhu, Adam, Yuille, and Chen] Siyuan Qiao, Yukun Zhu, Hartwig Adam, Alan Yuille, and Liang-Chieh Chen.  Vip-deeplab: Learning visual perception with depth-aware video panoptic segmentation.  In _CVPR_ , pages 3997–4008, 2021.
  * [Ranftl et al.(2019)Ranftl, Lasinger, Hafner, Schindler, and Koltun] René Ranftl, Katrin Lasinger, David Hafner, Konrad Schindler, and Vladlen Koltun.  Towards robust monocular depth estimation: Mixing datasets for zero-shot cross-dataset transfer.  _arXiv preprint arXiv:1907.01341_ , 2019.
  * [Ranftl et al.(2021)Ranftl, Bochkovskiy, and Koltun] René Ranftl, Alexey Bochkovskiy, and Vladlen Koltun.  Vision transformers for dense prediction.  _arXiv preprint arXiv:2103.13413_ , 2021.
  * [Schonberger and Frahm(2016)] Johannes L Schonberger and Jan-Michael Frahm.  Structure-from-motion revisited.  In _CVPR_ , pages 4104–4113, 2016.
  * [Teed and Deng(2018)] Zachary Teed and Jia Deng.  Deepv2d: Video to depth with differentiable structure from motion.  _arXiv preprint arXiv:1812.04605_ , 2018.
  * [Tiwari et al.(2020)Tiwari, Ji, Tran, Zhuang, Anand, and Chandraker] Lokender Tiwari, Pan Ji, Quoc-Huy Tran, Bingbing Zhuang, Saket Anand, and Manmohan Chandraker.  Pseudo rgb-d for self-improving monocular slam and depth prediction.  In _ECCV_ , pages 437–455, 2020.
  * [Ummenhofer et al.(2017)Ummenhofer, Zhou, Uhrig, Mayer, Ilg, Dosovitskiy, and Brox] Benjamin Ummenhofer, Huizhong Zhou, Jonas Uhrig, Nikolaus Mayer, Eddy Ilg, Alexey Dosovitskiy, and Thomas Brox.  Demon: Depth and motion network for learning monocular stereo.  In _CVPR_ , pages 5038–5047, 2017.
  * [Wang et al.(2018)Wang, Buenaposada, Zhu, and Lucey] Chaoyang Wang, José Miguel Buenaposada, Rui Zhu, and Simon Lucey.  Learning depth from monocular videos using direct methods.  In _CVPR_ , pages 2022–2030, 2018.
  * [Wang et al.(2019)Wang, Pizer, and Frahm] Rui Wang, Stephen M Pizer, and Jan-Michael Frahm.  Recurrent neural network for (un-) supervised learning of monocular video visual odometry and depth.  In _CVPR_ , pages 5555–5564, 2019.
  * [Watson et al.(2021)Watson, Mac Aodha, Prisacariu, Brostow, and Firman] Jamie Watson, Oisin Mac Aodha, Victor Prisacariu, Gabriel Brostow, and Michael Firman.  The temporal opportunist: Self-supervised multi-frame monocular depth.  In _CVPR_ , pages 1164–1174, 2021.
  * [Yin et al.(2021)Yin, Zhang, Wang, Niklaus, Mai, Chen, and Shen] Wei Yin, Jianming Zhang, Oliver Wang, Simon Niklaus, Long Mai, Simon Chen, and Chunhua Shen.  Learning to recover 3d scene shape from a single image.  In _CVPR_ , pages 204–213, 2021.
  * [Yin and Shi(2018)] Zhichao Yin and Jianping Shi.  Geonet: Unsupervised learning of dense depth, optical flow and camera pose.  In _CVPR_ , pages 1983–1992, 2018.
  * [Yogamani et al.(2019)Yogamani, Hughes, Horgan, Sistu, Varley, O’Dea, Uricár, Milz, Simon, Amende, et al.] Senthil Yogamani, Ciarán Hughes, Jonathan Horgan, Ganesh Sistu, Padraig Varley, Derek O’Dea, Michal Uricár, Stefan Milz, Martin Simon, Karl Amende, et al.  Woodscape: A multi-task, multi-camera fisheye dataset for autonomous driving.  In _ICCV_ , pages 9308–9318, 2019.
  * [Zhou et al.(2018)Zhou, Ummenhofer, and Brox] Huizhong Zhou, Benjamin Ummenhofer, and Thomas Brox.  Deeptam: Deep tracking and mapping.  In _ECCV_ , pages 822–838, 2018.
  * [Zhou et al.(2017)Zhou, Brown, Snavely, and Lowe] Tinghui Zhou, Matthew Brown, Noah Snavely, and David G Lowe.  Unsupervised learning of depth and ego-motion from video.  In _CVPR_ , pages 1851–1858, 2017.
  * [Zou et al.(2018)Zou, Luo, and Huang] Yuliang Zou, Zelun Luo, and Jia-Bin Huang.  DF-Net: Unsupervised joint learning of depth and flow using cross-task consistency.  In _ECCV_ , 2018.
  * [Zou et al.(2020)Zou, Ji, Tran, Huang, and Chandraker] Yuliang Zou, Pan Ji, Quoc-Huy Tran, Jia-Bin Huang, and Manmohan Chandraker.  Learning monocular visual odometry via self-supervised long-term modeling.  _arXiv preprint arXiv:2007.10983_ , 2020.


# Opinion Dynamic with agents immigration

Zhong-Lin Han  Department of Physics, University of Science and Technology of
China, Hefei 230026, China  Yu-Jian Li  Department of Modern Physics,
University of Science and Technology of China, Hefei 230026, China  Bing-Hong
Wang  Department of Modern Physics, University of Science and Technology of
China, Hefei 230026, China

###### Abstract

Abstract

###### pacs:

89.65.-s, 02.50.Le, 07.05.Tp, 87.23.Kg

##  I Introduction

Recent years, a large class of interdisciplinary problems has been
successfully studied with statistical physics methods. Statistical physics
establishes the bridge from microscopic characteristics to macroscopic
behaviors, for systems containing a large number of interacting components.
Using both analytical and numerical tools, it has contributed greatly to our
understanding of various complex systems. In this paper, we are motivated by
the statistical physics of a sociological problem, namely, opinion dynamics.

As one of the classical and traditional research areas in both social science
and theory physics, opinion dynamics has attracted much attention. A lot of
models concerning the process of opinion formation, such as voter model,
bounded confidence model, have been proposed previously. Meanwhile, some of
recent studies discussed and described the opinion dynamics on both common
conditions and various complex networks.

The issue of individual mobility has become increasingly fundamental due to
the Human migration and human dynamic. The issue is also important in other
contexts such as the emergence of Cooperation among individuals [20] and
species coexistence in cyclic competing games [21]. Recently, some empirical
data of human movements have been collected and analyzed [22,23]. From the
standpoint for dynamic of complex systems, when individuals (nodes, agents)
are mobile, the edges in the topological structures are no longer fixed,
yielding more different results on that than before.

In our paper, we try to propose a new model combining conventional opinion
dynamics with agents immigration according to information transmission and
evolution. In our simulation, we finally find a series of results reflecting
special and different features of opinion dynamics with immigration. By
introducing a parameter $\alpha$ to control the weight of influence of
individual opinions, according to a recent study considering weight influence,
we find that there also exist an optimal value of $\alpha$ leading to the
shortest consensus time for all individuals on a isotropic plane we concern.
After presenting the results of simulation in different situations, we also
analysis the results of our model in mathematical way, which leads us finding
out what are the exactly direct factors impacting the exponent of weight of
individual opinions.

Figure 1: (Color online) For $N=2000$ and $\langle k\rangle=4$, the density
$\rho_{c}$ as a function of $\alpha$ for different values of $r$ in the case
where all cooperators contribute the same cost c per game. Every cooperator
contributes a cost $c=1$ in every neighborhood that it plays.

In this paper we found up a new type of model for information dynamics with
immigration and we state related parameters and rules of our model. In order
to demonstrate the rationality of it, we presented both computational
simulation and mathematical analysis. Another most significant thing is that
we have designed a new mathematical method for model with linear
algebra.Compared with the previous method,we finally finish a complete model
for opinion and information dynamics with immigration.

##  II Model

As previous classical model focusing on the material process of spread and
formation of opinions, we spend more efforts on finding special results when
individuals carry opinions with immigration. To focus on a more efficient
situation, we just confine our discussion on an isotropic plane, without
special network effects. On the other hand, the individuals we concern are
just holding two kinds of opinion, the positive opinion $\psi_{+}>0$ and the
negative opinion $\psi_{-}<0$. According to one model on opinion dynamics
proposed before (.), we introduce the weight exponent to control the weight of
influence of each individual. We describe that all of the opinions of
individuals evolve simultaneously completely rely on its neighbors’ opinion
and neighbors’ weight. Here we describe the evolution process of the whole
individuals on the plane in mathematical way, which could be denoted as follow
(yang han xing)

$\displaystyle
p_{+}=\frac{\sum_{i}^{u}\omega_{i}^{\alpha}}{\sum_{i}^{u}\omega_{i}^{\alpha}+\sum_{j}^{v}\omega_{j}^{\alpha}},$
(1) $\displaystyle
p_{-}=\frac{\sum_{i}^{v}\omega_{i}^{\alpha}}{\sum_{i}^{u}\omega_{i}^{\alpha}+\sum_{j}^{v}\omega_{j}^{\alpha}}$
(2)

where $p_{+}$ and $p_{-}$ denote the probability of choosing positive opinion
and negative opinion,and the number of neighbors holding positive opinion is
$u$ while the number of negative ones is $v$. Here the model considers the
agents with weight impact $\omega_{i}$, which is controlled by weight exponent
$\alpha$. If the probability $p_{+}$ is lager than $p_{-}$, the agent we
concern will choose positive opinion at the next step. And it will be same as
choosing negative opinion.

Figure 2: (Color online) Cumulative payoff distribution for different values
of $\alpha$. The distribution is obtained after the cooperation density
becomes stable. The multiplication factor is set to be $r=1.6$. Solid curves
are theoretical predictions from Eq. (LABEL:eq:wealthdis).

In this model, the individual we concern evolves its opinion at $t+1$
according to its neighbors’ opinion in its view radius $r$, which is shown in
Figure.1. In Figure.1, the red agents hold positive opinions and black agents
hold negative opinions. There are $u+v$ neighbors in the view range of
individual we concern, while here are $u$ individuals hold positive opinion
and $v$ individuals hold negative opinion at $t$ step. After comparing the
weight of positive opinion and negative opinion, the individual we concern
evolves its opinion at $t+1$ step as this

$\psi_{i}^{(t+1)}=\sum_{r}\psi_{j}^{(t)}$ (3)

where $\psi_{j}^{(t)}$ denotes the opinion state of the $j$th neighbor of the
$i$th agent we concern at the $t$ step.

Figure 3: (Color online) Times series of cooperator density in hubs’
neighborhoods for (a) The multiplication factor is $r=1.2$ and each data point
is obtained by averaging over 50 runs.

After changing their opinion in the way above, all the individuals immigrate
on the plane. All the agents would be confined in the plane by periodic
boundary condition. The velocity and direction angle of each agent are
randomly distributed, which are kept by each agent all the time. After enough
period of time, number of the individuals holding positive opinions $N_{+}$
and the number of other individuals holding negative ones $N_{-}$ reach a
plateaus and dynamic equilibrium. At that certain point, we believe that the
process of opinion dynamics would be terminated. And we could find that if all
of the individuals enter the plateaus, the total number of individuals who
hold positive opinion at $t$ step $\psi_{+}^{t}$ would be approximately equal
to the number of individuals who also hold that at $t+1$ step
$\psi_{+}^{t+1}$. And we could carry on this description with mathematical
language,

$\sum_{i}^{N}\psi_{i}^{(t+1)}=\sum_{i}^{N}\psi_{i}^{(t)}$ (4)

The total time steps the system took could be defined as $T_{c}$ for
convergent time.

###  II.1 Results and analysis

In the following discussion and simulation, we confines our individuals on an
isotropic plane $(L\times L)$. The length of the boundary of this plane $L$ is
20, and the total number of individuals on the plane would be $N$. Here we
simulate the individuals have their initial velocity under Gauss distribution,
which would be more rational and close to facts. Each individual hold their
opinions (positive one or negative one) and their fixed weight of opinion with
random probability. The distribution of agents’ weight was established in a
random way at the beginning of evolution. The exponent $\alpha$ in equations
(1) and (2) controls the evolution process. And here we define $\rho$ and
$\Delta\rho$ to describe the changing process of the individuals holding
positive opinion. They are denoted in equations as follow

$\rho_{c}=\frac{N_{+}}{N}$ (5) $\Delta\rho=\frac{\Delta N_{+}}{N}$ (6)

In Figure 2,we show $\rho_{c}$ as function of evolution time $t$ for different
view radius $r$, both $r$=1.2 and $r$=1.5. The most interesting thing we could
find in this figure is that when evolution time $t$ is around 6500, the value
of $\Delta\rho$ plummets obviously, which finally reach the level under 0.1.
In fact, when changes of $\Delta\rho$ has lower amplitude of variation, it
also means that the individuals holding positive opinion enter the period of
dynamic equilibrium.

Figure 4: (Color online) For $r=1.6$, cooperator density $\rho_{c}$ as a
function of degree for different values of $\alpha$. Figure 5: (Color online)
For $r=1.6$, cooperator density $\rho_{c}$ as a function of degree for
different values of $\alpha$. Figure 6: (Color online) (a)[Initial
distribution of agents with two kinds of opinions,$N_{+}$:$N_{-}$=1.19:1.
(b)Final distribution of agents with two kinds of opinions,
$N_{+}$:$N_{-}$=2.91:1.

In this situation, we could find the consensus time in Figure 3. If $t<4000$
or $t>6500$, the $\Delta\rho$ changes in a very small range which would also
be shown in the figure. But the consensus time is not directly impacted by
only view radius $r$ of each individual. In Figure 4, here converge time
$T_{c}$ is taken as a function of average velocity $\alpha$ , while the view
radius $r$ is equal to 1.2. Interestingly, we show that would reach a minimum
value when $\alpha$ is around 2 under different values of total number of
individuals on the plane $N$. Here we present that consensus time $T_{c}$
changes with $\alpha$ in a ”smile curve”. And certainly the value of would be
higher if the $N$ is more. In fact, it is obvious to explain that when there
are more individuals holding different opinions, they would take more time to
reach consensus or dynamic equilibrium. In that we show that $N$ and weight
exponent $\alpha$ could both directly determine the consensus time $T_{c}$.
The more cogent demonstration is shown in Figure 2, which presents as a
function of $\alpha$. Here $\rho_{c}$ is the density of individuals holding
positive opinions at the consensus time.

In the Figure 2, we present that $\rho_{c}$ will reach a maximum when is
around 2. Meanwhile, the value of is greater if view radius is lager. To show
the result in a more intuitive way, we present the distribution map in Figure
6. In Figure 6, the red points are the ones holding positive opinions while
the black points present negative ones. In this figure, we present the
specific distribution.

To support former results, we mainly focus on the results that positive
opinions take dominant rate. In order to discuss the parameters reflecting
immigration of individuals, we present consensus time $T_{c}$ as a function of
average velocity of individuals with different $N$ in Figure 5. In this
figure, we find that consensus time $T_{c}$ would increase approximately in a
linear way when $v$ is less than 1.5. After that , it decreases in a certain
range without sharp changes.

To discuss the model in a more reliable way, we try to analysis the process by
founding up a series of equations for $m$ agents in total as follow. In
equations, we define that the $i$th individual we concern has a view radius
$r$, and at the t step there are $s_{m}$ individuals in its view range as its
neighbors,and here we define that $s_{ij}$ as the $i$th neighbor of the $j$th
agent we concern. In that, it is $j$th opinion state of agent’s neighbor at
$t$ step that determine the opinion updating of this individuals at next time
step $t+1$. If $\psi$ is positive, individuals who holding positive opinions
would have greater weight than those who hold negative opinions. As a
consequence, the equations could be formed as follow:

$\displaystyle\left\\{\begin{array}[]{c}\psi_{1}^{(t+1)}=\psi^{(t)}(r,s_{11})\cdot\omega_{s_{11}}^{\alpha}+\ldots+\psi^{(t)}(r,s_{m-1,1})\cdot\omega_{s_{m-1,1}}^{\alpha}+\psi^{(t)}(r,s_{m,1})\cdot\omega_{s_{m,1}}^{\alpha}\\\
\vdots\\\
\psi_{m-1}^{(t+1)}=\psi^{(t)}(r,s_{1,m-1})\cdot\omega_{s_{11}}^{\alpha}+\ldots+\psi^{(t)}(r,s_{m-1,m-1})\cdot\omega_{s_{m-1,m-1}}^{\alpha}+\psi^{(t)}(r,s_{m,m-1})\cdot\omega_{s_{m,m-1}}^{\alpha}\\\
\psi_{1}^{(t+1)}=\psi^{(t)}(r,s_{1,m})\cdot\omega_{s_{1,m}}^{\alpha}+\ldots+\psi^{(t)}(r,s_{m-1,m})\cdot\omega_{s_{m-1,m}}^{\alpha}+\psi^{(t)}(r,s_{m,m})\cdot\omega_{s_{m,m}}^{\alpha}\\\
\end{array}\right.$ (11)

To describe the model in a simpler way, we try to apply linear algebra instead
of these traditional equations. In order to write in that way, we also
introduce a new parameter $n_{ij}^{(t)}$ into this matrix description. Here
$n_{ij}^{(t)}$ reflects that the times of opinion exchanging or sharing of
$j$th individual we concern at $t$ step. In other word, $n_{ij}^{(t)}$ is a
standard that concerns how many times the $i$th individual impacts others
opinion updating choice of next time step at $t$ step. When the whole
individuals get into the plateaus of dynamic equilibrium, we discussed in
Sec.2, the opinions individuals holding would be described as where
$\omega(i)$ is the weight of the ith agent, which is

$\left(\begin{array}[]{c}\psi_{1}^{(t+1)}\\\ \psi_{2}^{(t+1)}\\\ \vdots\\\
\psi_{m-1}^{(t+1)}\\\ \psi_{m}^{(t+1)}\\\
\end{array}\right)=\left(\begin{array}[]{ccccc}\psi_{1}^{(t)}&\psi_{2}^{(t)}&\ldots&\psi_{m-1}^{(t)}&\psi_{m}^{(t)}\\\
\psi_{1}^{(t)}&\vdots&\vdots&\vdots&\vdots\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
\vdots&\vdots&\vdots&\vdots&\psi_{m}^{(t)}\\\
\psi_{1}^{(t)}&\psi_{2}^{(t)}&\ldots&\ldots&\psi_{m}^{(t)}\end{array}\right)\times\left(\begin{array}[]{cccc}n_{11}^{(t)}\cdot\omega_{1}^{\alpha}&n_{12}^{(t)}\cdot\omega_{1}^{\alpha}&\ldots&n_{1m}^{(t)}\cdot\omega_{1}^{\alpha}\\\
n_{21}^{(t)}\cdot\omega_{2}^{\alpha}&n_{22}^{(t)}\cdot\omega_{2}^{\alpha}&\ldots&n_{2m}^{(t)}\cdot\omega_{2}^{\alpha}\\\
\vdots&\vdots&\ldots&\vdots\\\
n_{m-1,1}^{(t)}\cdot\omega_{m-1}^{\alpha}&n_{m-1,2}^{(t)}\cdot\omega_{m-1}^{\alpha}&\ldots&n_{m-1,m}^{(t)}\cdot\omega_{m-1}^{\alpha}\\\
n_{m,1}^{(t)}\cdot\omega_{m}^{\alpha}&n_{m,2}^{(t)}\cdot\omega_{m}^{\alpha}&\ldots&n_{m,m}^{(t)}\cdot\omega_{m}^{\alpha}\end{array}\right)$
(12)

Here we denote that $n_{ij}^{(t)}$ as the times for opinion exchanges between
the $i$th and $j$th agents at the $t$ time step. At that certain situation, we
would find that When the whole system has entered the final homeostasis, the
whole opinions of agents we concern would be invariable, which means that .
And we could finally find that the time of opinion sharing or exchanging is
related to the weight of the agent and exponent in the equation. After
simplified such matrix equation, we finally get a direct function for
$n_{ij}^{(t)}$ and $\alpha$

$\sum_{j=1}^{m}n_{ij}^{(t)}\omega_{i}^{\alpha}=1$ (13)

where $\omega_{i}$ is weight of the $i$th agent we concern.
By choosing five different groups $n_{ij}^{(t)}$ when we fixed $\alpha=2$, we
finally calculate the $\alpha$ with the function for support that, from which
we get $\alpha=1.98,1.93,1.82,1.88,1.91$.By the calculation, the function
reflects a very important and simple relationship between $n_{ij}^{(t)}$ and
weight $\omega_{i}$

##  III Conclusion and Discussions

The new mode of opinion evolution with immigration is different from the
conventional opinion dynamics. It presents that density of positive opinion
agents would be maximum when the weight exponent $\alpha$ is around 2. In
summary, we found up a new model for opinion exchange and communication among
agents with immigration. The state matrix we present for analysis and
quantitative simulation could also be widely used for more complex situation.
The opinion carried by agents represent a kind of state or parameter of agents
in motion.So more application and analysis could be carried on with this model
and method in future. Discussion we present above is not only to demonstrate
our model, but also open up a new combination between opinion communication
and agent-based motion. State consensus time is also a very important
parameter to describe a system or a group of agents, which could also be one
certain standard for different situations.

###### Acknowledgements.

This work is funded by the National Basic Research Program of China (973
Program No.2006CB705500), the National Natural Science Foundation of China
(Grant Nos. 60744003, 10635040, 10532060) and by the Special Research Funds
for Theoretical Physics Frontier Problems (NSFC No.10547004 and A0524701). WXW
and YCL are supported by AFOSR under Grant No. FA9550-07-1-0045.

## References

  * (1) A. M. Colman, Game Theory and Its Applications in the Socia and Biological Sciences (Butterworth-Heinemann, Oxford, 1995).
  * (2) J. M. Smith, Evolution and the Theory of Games (Cambridge University Press, Cambridge, England, 1982).
  * (3) H. Gintis, Game Theory Evolving (Princeton University Press, Princeton, NJ, 2000).
  * (4) R. Axelrod and W. D. Hamilton, Science 211, 1390 (1981).
  * (5) C. Hauert and M. Doebeli, Nature (London) 428, 643 (2004); M. Sysi-Aho, J. Saramäki, J. Kertész, and K. Kaski, Eur. Phys. J. B 44, 129 (2005); L.-X. Zhong, D.-F. Zheng, B. Zheng, C. Xu, and P. M. Hui, Europhys. Lett. 76, 724 (2006); A. Szolnoki and G. Szabó, ibid. 77, 30004 (2007); G. Abramson and M. Kuperman, Phys. Rev. E 63, 030901(R) (2001); B. J. Kim, A. Trusina, P. Holme, P. Minnhagen, J. S. Chung, and M. Y. Choi, ibid. 66, 021907 (2002); H. Ebel and S. Bornholdt, ibid. 66, 056118 (2002); Z.-X. Wu, X.-J. Xu, Y. Chen, and Y.-H. Wang, ibid. 71, 037103 (2005); W.-X. Wang, J. Ren, G. Chen, and B.-H. Wang, ibid. 74, 056113 (2006); C.-L. Tang, W.-X. Wang, X. Wu, and B.-H. Wang, Eur. Phys. J. B 53, 411 (2006); J. Ren, W.-X. Wang, and F. Qi, Phys. Rev. E 75, 045101(R) (2007); Z. Rong, X. Li, and X. Wang, ibid. 76, 027101 (2007); X. Chen and L. Wang, ibid. 77, 017103 (2008); F. Fu and L. Wang, ibid. 78, 016104 (2008); F. Fu, C. Hauert, M. A. Nowak, and L. Wang, ibid. 78, 026117 (2008); W.-X. Wang, J. Lü,2 G. Chen, and P. M. Hui, ibid. 77, 046109 (2008).
  * (6) J. H. Kagel and A. E. Roth, $The$ $Handbook$ $of$ $Experimental$ $Economics$, (Princeton University Press, Princeton, NJ, 1997).
  * (7) C. Hauert, S. De Monte, J. Hofbauer, and K. Sigmund, Science 296, 1129 (2002).
  * (8) G. Szabó and C. Hauert, Phys. Rev. Lett. 89, 118101 (2002).
  * (9) C. Hauert, S. De Monte, J. Hofbauer, and K. Sigmund, J. Theor. Biol. 218, 187 (2002).
  * (10) D. Semmann, H. J. Krambeck, and M. Milinski, Nature (London) 425, 390 (2003).
  * (11) H. Brandt, C. Hauert, and K. Sigmund, Proc. Natl. Acad. Sci. 103, 495 (2006).
  * (12) J. Y. Guan, Z. X. Wu, and Y. H. Wang, Phys. Rev. E 76, 056101 (2007).
  * (13) F. C. Santos, M. D. Santos, and J. M. Pacheco, Nature (London) 454, 213 (2008).
  * (14) See, for example, K. Wiesenfeld and F. Moss, Nature (London) 373, 33 (1995); L. Gammaitoni, P. Hänggi, P. Jung, and F. Marchesoni, Rev. Mod. Phys. 70, 223 (1998).
  * (15) G. Szabó and C. Töke, Phys. Rev. E 58, 69 (1998); G. Szab and C. Hauert, Phys. Rev. Lett. 89, 118101 (2002); G. SzabÓ and J. Vukov, Phys. Rev. E 69, 036107 (2004).
  * (16) M. Perc and A. Szolnoki, Phys. Rev. E 77, 011904 (2008).
  * (17) Z. X. Wu, X. J. Xu, Z. G. Huang, S. J. Wang, and Y. H. Wang, Phys. Rev. E 74, 021107 (2006).
  * (18) J. Y. Guan, Z. X. Wu, Z. G. Huang, X. J. Xu and Y. H. Wang, Europhys. Lett. 76, 1214 (2006).
  * (19) J. Ren, W. X. Wang, G. Yan, and B. H. Wang, arXiv: physics/0603007v1.
  * (20) G. Szabó and J. Vukov, Phys. Rev. E 69, 036107 (2004).
  * (21) G. Szabó, J. Vukov, and A. Szolnoki, Phys. Rev. E 72, 047107 (2005).
  * (22) M. Perc, Phys. Rev. E 75, 022101 (2007).
  * (23) Z. X. Wu and Y. H. Wang, Phys. Rev. E 75, 041114 (2007).
  * (24) A. L. Barabási and R. Albert, Science 286, 509 (1999).
  * (25) G. Szabó and G. Fáth, Phys. Rep. 446, 97 (2007).
  * (26) C. Hauert and G. Szabó, Complexity 8, 31 (2003).
  * (27) C. Hauert, Adv. Complex Syst. 9, 315 (2006).
  * (28) F. C. Santos and J. M. Pacheco, Phys. Rev. Lett. 95, 098104 (2005).
  * (29) F. C. Santos, J. M. Pacheco, and T. Lenaerts, Proc. Natl Acad. Sci. 103, 3490 (2006).
  * (30) V. Pareto, Le Cours d’ Économie Politique (Macmillan, Lausanne, Paris, 1987).
  * (31) M.-B. Hu, W.-X.Wang, R. Jiang, Q.-S. Wu, B.-H. Wang, and Y.-H. Wu, Eur. Phys. J. B 53, 273 (2006).
  * (32) We have also studied the prisoner’s dilemma game on scale-free networks and found the phenomenon of diversity-optimized cooperation.